二維光子晶體瑕疵共振腔內模態之分析

全文

(1)國立交通大學光電工程研究所碩士論文. 二維光子晶體瑕疵共振腔內模態之分析 Analysis of Defect Modes of Two-Dimensional Photonic Crystal Nano-Cavities. 研究生 : 趙竑鈞指導教授 : 李柏璁. 教授. 中華民國九十四年七月.

(2) Abstract in Chinese 在這篇論文裡研究的是二維光子晶體瑕疵共振腔內的模態。首先, 我們介紹時域有限差分法和相關的邊界值條件。其次, 一些對稱性的分析和共振腔的設計原則將被討論。最後, 我們將展示一些模擬的技巧和光子晶體模擬的結果。. i.

(3) Abstract in English In this thesis, defect modes of two-dimensional photonic crystal nano-cavity are analyzed. The defect is formed by removing a single air hole from array of air holes of hexagonal lattice arrangement on a dielectric slab. First, the finite difference time domain method with various boundary conditions is introduced. Second, symmetry analysis of defect modes and design rules for high quality factor cavities are presented. Finally, techniques of the simulation and the results of the simulation, like photonic band structures, resonant frequencies, mode profiles, and quality factors are exhibited.. ii.

(4) Acknowledgements 首先要衷心感謝我的指導老師李柏璁教授。這兩年來, 有幸跟隨老師學習與研究, 不但令我於光電此一領域有著更深入的體驗之外, 也習得了更嚴謹的做研究態度。而這一段日子以來, 自己的學習總是不夠積極, 但老師仍然不厭其煩地指導我, 並包容我的怠惰, 在此要深深地對老師致上最深的謝意。此外, 亦要感謝行政院國家科學委員會(計畫編號:NSC-92-2218-E-009-023, NSC-93-2215-E009-059), 以及經濟部學界科專計畫 (計畫編號:93-EC-17-A-07-SI-0011), 讓我在研究的過程中得到足夠的研究資源及經費。而我也要向實驗室中諸位夥伴們致謝。兩年來, 實驗室的氣氛總是如此愉悅融洽, 大夥兒一同埋頭研究的感覺如此美好。每當遭遇難關之際, 大家總是並肩苦思、齊心協力解決難題, 而在突破瓶頸的那一剎那, 大夥兒亦是一塊兒分享那成功的喜悅。我之所以能夠完成這篇論文, 同儕之間的鼓勵與打氣亦是居功厥偉。最後, 也要感謝我的家人與關心我的朋友們, 非常感謝大家的支持與鼓勵。感謝大家。. iii.

(5) Contents. Abstract in Chinese. i. Abstract in English. ii. Acknowledgements. iii. Contents. iv. List of Tables. vi. List of Figures. viii. 1 Introduction. 1. 1.1. Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 1. 1.2. Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 2. 1.3. Thesis Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 3. 2 Simulation Principles for the FDTD Method. 4. 2.1. Maxwell’s Equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 4. 2.2. Yee’s Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 6. 2.3. Heaviside-Lorentz Unit System . . . . . . . . . . . . . . . . . . . . . . . .. 7. 2.4. Sullivan’s Implement for Perfectly Matched Layer Absorbing Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 8. 2.5. My Implement for PML ABC . . . . . . . . . . . . . . . . . . . . . . . . . 10. 2.6. Bloch’s Boundary Conditions for Periodic Structures . . . . . . . . . . . . 11. 2.7. Symmetric Boundary Conditions . . . . . . . . . . . . . . . . . . . . . . . 11. iv.

(6) 2.8. Poynting Vector and EM-Energy in Yee’s Lattice . . . . . . . . . . . . . . 13. 2.9. Time Step Issues . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14. 3 Symmetry Analysis 3.1. 3.2. 3.3. 16. Some Terminologies and Useful Theorems . . . . . . . . . . . . . . . . . . 16 3.1.1. Basic Concepts of Group . . . . . . . . . . . . . . . . . . . . . . . . 16. 3.1.2. Representation of a Group . . . . . . . . . . . . . . . . . . . . . . . 18. 3.1.3. Basis Functions and Projection Operators . . . . . . . . . . . . . . 20. Classification of Defect Modes in Hexagonal Lattice . . . . . . . . . . . . . 22 3.2.1. Symmetry Group of Photonic Crystals . . . . . . . . . . . . . . . . 22. 3.2.2. Symmetry of 2D Photonic Crystals of Finite Thickness . . . . . . . 23. 3.2.3. Symmetry of Defect Modes. . . . . . . . . . . . . . . . . . . . . . . 24. Some Design Rules . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25. 4 Simulation Results 4.1. 4.2. 27. Photonic Band Structures . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 4.1.1. 2D PWE Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27. 4.1.2. 3D PWE Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28. 4.1.3. 3D FDTD Method . . . . . . . . . . . . . . . . . . . . . . . . . . . 29. Defect Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 4.2.1. Resonant Frequencies of Defect Modes . . . . . . . . . . . . . . . . 31. 4.2.2. Mode Profiles . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35. 4.2.3. Quality Factors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39. 5 Summary. 40. References. 41. v.

(7) List of Tables 2.1. Comparison of Maxwell’s equations in different unit systems. . . . . . . . .. 8. 3.1. Character table for the C6v . . . . . . . . . . . . . . . . . . . . . . . . . . . 19. 3.2. Character table for the C1h . . . . . . . . . . . . . . . . . . . . . . . . . . . 21. 4.1. Comparison of photonic band gaps using different methods. The frequencies are taken to be normalized (a/λ). . . . . . . . . . . . . . . . . . . . . . 31. 4.2. Parameter table. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35. 4.3. Normalized peak frequencies. . . . . . . . . . . . . . . . . . . . . . . . . . . 36. 4.4. The quality factors. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39. vi.

(8) List of Figures 1.1. Painter’s structure. (Adopted from reference [3]) . . . . . . . . . . . . . . .. 2. 2.1. Yee’s unit cell. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. 6. 2.2. Bloch’s boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . . . 11. 2.3. Symmetric boundary conditions. . . . . . . . . . . . . . . . . . . . . . . . . 12. 3.1. Symmetry operations for two-dimensional hexagonal lattice. . . . . . . . . 17. 3.2. First Brillouin zone of a hexagonal lattice. K stands for the set of all Ki and M stands for the set of all Mi for i = 1 · · · 6. . . . . . . . . . . . . . . . 18. 3.3. light cone . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26. 4.1. The TE band structure by 2D PWE method.. 4.2. Hz mode profiles at band edges. . . . . . . . . . . . . . . . . . . . . . . . . 28. 4.3. The TE-like band structure by 3D PWE. . . . . . . . . . . . . . . . . . . . 29. 4.4. Double unit cell and the locations of sources and monitors. . . . . . . . . . 29. 4.5. The band structure by the FDTD method. (r/a = 0.3, d/a = 0.4, n = 3.4). 4.6. Various patterns by changing hole radii. The radius of holes of most inner. . . . . . . . . . . . . . . . . 28. 30. layer is 0.3a. (a) r/a gradually decrease by 0.02 per layer outward. (b) r/a gradually decrease by 0.01 per layer outward. (c) typical pattern without changing the hole radii. (d) r/a gradually increase by 0.01 per layer outward. (e) r/a gradually increase by 0.02 per layer outward. . . . . . . . 32 4.7. The band structures of defect system Type C( ar is fixed to be 0.3) by the 2D and 3D PWE methods. . . . . . . . . . . . . . . . . . . . . . . . . . . . 33. 4.8. Excitation Signal. (α = 1024 and f0 = 0.015) . . . . . . . . . . . . . . . . . 34. 4.9. The locations of source and monitors. . . . . . . . . . . . . . . . . . . . . . 34. 4.10 The corresponding spectrums. . . . . . . . . . . . . . . . . . . . . . . . . . 36. vii.

(9) 4.11 Definition of domain of supercell and profiles of doubly degenerate defect modes by 2D PWE method. . . . . . . . . . . . . . . . . . . . . . . . . . . 37 4.12 Dipole modes by 3D FDTD method. Light cone is shown in red circle. . . 38. viii.

(10) Chapter 1 Introduction 1.1. Background. In 1987, Yablonovitch and John generalized the idea of distributed Bragg reflection(DBR), which is the reflection of electromagnetic waves in periodic arrangement of multilayer films of different dielectric constants, from one dimension to multi-dimension[1]. The generalization inspired the name “photonic crystal”. DBR is an example of onedimensional photonic crystals. If the materials are periodic in two directions, they are called two-dimensional photonic crystals. Three-dimensional photonic crystals imply that the periodicity of material is in all directions[2]. Photonic crystals have many interesting phenomenons like photonic band gap, defect modes, negative refraction, superprism effect, etc. Photonic band gap prohibits light with certain frequencies to propagate, which is the analogue to band gap of electron states in crystalline materials. Here we term the “air band” corresponding to the conduction band and “dielectric band” to the valence band. Defect modes result from locally changing of dielectric constant as a perturbation of the periodic dielectric structure, which is the analogue to defect states of electron states in crystalline materials. Applications of photonic crystals to the development of semiconductor lasers have been addressed in recent years. In 1999, Painter proposed a new structure of photonic crystal laser by the use of defects in two-dimensional(2D) photonic crystal slabs[3]. The optical cavity is formed by removing a single air hole from array of air holes of hexagonal lattice arrangement on a dielectric slab. And the optical gain is provided by the layers of quantum wells within the slab pumped by light. Defect cavities in 2D photonic crystal. 1.

(11) slab confine light through two mechanisms. One is total internal reflection in the vertical direction, and the other is DBR in plane. Schematic depiction of Painter’s work is shown is Figure 1.1.. Figure 1.1: Painter’s structure. (Adopted from reference [3]). Several improvements of the photonic crystal lasers have been fabricated and investigated by modifying the geometry of the photonic crystal patterns. For examples, Painter et al. split the doubly degenerate dipole modes by enlarging or shifting the air holes next to the single defect[3], [4]. Hong-Gyu Park et al. modified six nearest neighbor holes around the defect to observe the monopole defect mode[5], [6]. Vˇ ucković et al. used a type of defect by elongating holes along the symmetry axes to improve the vertical loss[7], [8]. Po-Tsung Lee et al. made low threshold photonic crystal lasers for room temperature operation by lithographic tuning of cavities of 19 missing holes[9], [10]. Future efforts will be focused on electrically driven structures with a good heat sink to solve heat problems.. 1.2. Motivation. Purcell made a prediction that the rate of spontaneous emission may be enhanced in a resonant cavity[11]. The Purcell factor Fp is defined as Fp ≡. τ0 Γ ≡ , Γ0 τ. (1.1). where Γ is the spontaneous emission rate and τ is the lifetime. In optimal quantum well structures it can be derived explicitly by[12] Fp =. 3Qg λ 3 , 2πVef f 2n 2. (1.2).

(12) where Q is the quality factor of the resonant cavity, g is the mode degeneracy, λ is the resonant wavelength, n is the refractive index of the materials, and Vef f is the effective mode volume defined by R. Vef f =. ε(r)E 2 (r)d3 r , max(ε(r)E 2 (r)). (1.3). where E is the electric field, and ε is the permittivity of the material. Therefore the design of an optical cavity for higher quality factor and smaller mode volume is the key to refine photonic crystal lasers.. 1.3. Thesis Overview. In this chapter, we have a brief introduction of photohic crystals and their application to semiconductor lasers. We begin in chapter 2 by introducing the method of finite difference time domain(FDTD) developed by Yee, and some techniques used in my simulations like Heaviside Lorentz unit system. The implement of various boundary conditions and expressions of Poynting vector and electromagnetic energy density in Yee’s lattice will be addressed. In chapter 3 a study of symmetry analysis and design rules for high quality factor photonic crystal cavities will be discussed. Simulation procedure and results by both the methods of plane wave expansion(PWE) and FDTD are demonstrated in chapter 4. In the end, the final conclusion will be presented in chapter 5.. 3.

(13) Chapter 2 Simulation Principles for the FDTD Method 2.1. Maxwell’s Equations. Consider a region of space which is source-free, Maxwell’s curl equations are given by 1 ρ0 ∂H = − ∇ × E − H, ∂t µ µ 1 σ ∂E = ∇ × H − E. ∂t ε ε. (2.1). in MKS unit system where E is the electric field in volts/meter; H is the magnetic field in amperes/meter, ε is the electrical permittivity in farads/meter, σ is the electrical conductivity in ohms/meter(siemens/meter), µ is the magnetic permeability in henrys/meter, and ρ0 is an equivalent magnetic resistivity in ohms/meter. The magnetic resistivity term is provided to yield symmetric curl equations and allow for the possibility of a magnetic field loss mechanism. Because of small variation of magnetic permeability in materials, the effect of magnetization is ignored. So the magnetic permeability of materials µ is assumed to be a constant as the magnetic permeability in vacuum, say µ0 . The material is taken to be lossless for simplicity, which implies that σ and ρ0 are 0. Moreover, ε is assumed to be nondispersive and isotropic. The following scalar equations is equivalent to Maxwell’s curl equations in the rectangular coordinate system (x, y, z). 1 ∂Ey ∂Ez ∂Hx = ( − ), ∂t µ0 ∂z ∂y. 4.

(14) ∂Hy ∂t ∂Hz ∂t ∂Ex ∂t ∂Ey ∂t ∂Ez ∂t. = = = = =. 1 ∂Ez ∂Ex ( − ), µ0 ∂x ∂z 1 ∂Ex ∂Ey − ), ( µ0 ∂y ∂x 1 ∂Hz ∂Hy ( − ), ε ∂y ∂z 1 ∂Hx ∂Hz ( − ), ε ∂z ∂x 1 ∂Hy ∂Hx ( − ). ε ∂x ∂y. (2.2). The six coupled partial differential equations forms the basis of the FDTD algorithm for electromagnetic field interactions with general three-dimensional objects. Before proceeding with the details of the algorithm, it is informative to consider one important simplification of the full three-dimensional case. Assuming that neither the incident plane wave excitation nor the modeled geometry has any variation in the z-direction (i.e., all partial derivatives with respect to z equal zero), Maxwell’s curl equations reduce to two decoupled sets of scalar equations. These decoupled sets which are termed the transverse magnetic (TM) mode and the transverse electric (TE) mode describe two-dimensional wave interactions with objects. Another viewpoint of classification of electromagnetic fields which gets the same result will be addressed in Section 3.2. The relevant equations for each case are as follow: -TM case (Ez, Hx, and Hy field components only) 1 ∂Ez ∂Hx =− , ∂t µ0 ∂y 1 ∂Ez ∂Hy = , ∂t µ0 ∂x 1 ∂Hy ∂Hx ∂Ez = ( − ). ∂t ε ∂x ∂y. (2.3). -TE case (Hz, Ex, and Ey field components only) 1 ∂Hz ∂Ex = , ∂t ε ∂y ∂Ey 1 ∂Hz =− , ∂t ε ∂x 1 ∂Ex ∂Ey ∂Hz = ( − ). ∂t µ0 ∂y ∂x. 5. (2.4).

(15) 2.2. Yee’s Algorithm. In 1966, Yee introduced a set of finite-difference equations for equations (2.2). Following Yee’s notation, we denote a space point in a rectangular lattice as (i, j, k) = (i∆x, j∆y, k∆z),. (2.5). F n (i, j, k) = F (i∆x, j∆y, k∆z, n∆t),. (2.6). where ∆x, ∆y, and ∆z are the lattice space increments in the x-, y-, z-coordinate directions, ∆t is the time increment, and i, j, k, and n are integers. Yee used centered finite-difference expressions for the space and time derivatives that are both simply programmed and second-order accurate in the space and time increments respectively: F n (i + 12 , j, k) − F n (i − 12 , j, k) ∂F n (i, j, k) = + O(∆x2 ), ∂x ∆x 1 n+ 21 n ∂F (i, j, k) (i, j, k) − F n− 2 (i, j, k) F = + O(∆t2 ). ∂t ∆t. (2.7) (2.8). To achieve the accuracy of (2.7) and to realize all of the required space derivatives of equations (2.2), Yee positioned the components of E and H in an unit cell of the lattice as shown in Figure 2.1. To achieve the accuracy of (2.8), Yee used a leap-frog algorithm which evaluates E and H at alternate half time steps. The followings are sample finite-difference time-stepping expressions for a magnetic and an electric field component resulting from. /2. ,k. +1. /2. ). Hy(i-1/2,j,k+1/2). ) ,k ,j /2 i+ 1. Hz(i-1/2,j+1/2,k). Ey(i,j+1/2,k). Hx (. i,j. +1. /2. ,k. -1 /. 2). Ex (. Hz(i+1/2,j-1/2,k). Hy(i+1/2,j,k+1/2). Hz(i+1/2,j+1/2,k). Hx (. i,j. Hx. +1. (i, j. -1. /2. ,k. +1. /2. ). Z. Ez(i,j,k+1/2). these assumptions[13]:. Hy(i+1/2,j,k-1/2). Y X. Figure 2.1: Yee’s unit cell.. 6.

(16) n+ 21. Hx. 1 1 ∆t 1 1 n− 1 (i, j + , k + ) = Hx 2 (i, j + , k + ) + 2 2 2 2 µ0 1 1 1 ×{ [Eyn (i, j + , k + 1) − Eyn (i, j + , k)] ∆z 2 2 1 1 1 [E n (i, j, k + ) − Ezn (i, j, k + )]}, + ∆y z 2 2. (2.9). ∆t 1 1 Ezn+1 (i, j, k + ) = Ezn (i, j, k + ) + 2 2 ε(i, j, k + 12 ) 1 1 1 1 1 n+ 1 n− 1 ×{ [Hy 2 (i + , j, k + ) − Hy 2 (i − , j, k + )] ∆x 2 2 2 2 1 1 1 1 1 n+ 21 n+ 12 [Hx (i, j − , k + ) − Hx (i, j + , k + )]}. + ∆y 2 2 2 2 (2.10) With the system of finite-difference equations, the new value of a field vector component at any lattice point depends only on its previous value and on the previous values of the components of the other field vector at adjacent points.. 2.3. Heaviside-Lorentz Unit System. Under the consideration for numerical precision and the efficiency of calculation, we will use Maxwell’s equations in the Heaviside-Lorentz unit system and set the velocity of light in vacuum c0 equal to 1 in practice. Table 2.1 shows the definitions of ε0 , µ0 , D, H, macroscopic Maxwell’s equations, and Lorentz force per unit charge in various systems of units. Now the sample update expressions (2.9) and (2.10) in Heaviside-Lorentz unit system reduce to: 1 1 1 1 n− 1 (i, j + , k + ) = Hx 2 (i, j + , k + ) + ∆t 2 2 2 2 1 1 1 ×{ [Eyn (i, j + , k + 1) − Eyn (i, j + , k)] ∆z 2 2 1 1 1 + [E n (i, j, k + ) − Ezn (i, j, k + )]}, (2.11) ∆y z 2 2 ∆t 1 1 Ezn+1 (i, j, k + ) = Ezn (i, j, k + ) + 2 2 εr (i, j, k + 21 ) 1 1 1 1 1 n+ 1 n− 1 ×{ [Hy 2 (i + , j, k + ) − Hy 2 (i − , j, k + )] ∆x 2 2 2 2 1 1 1 1 1 n+ 21 n+ 21 + [Hx (i, j − , k + ) − Hx (i, j + , k + )]}, (2.12) ∆y 2 2 2 2. n+ 12. Hx. where εr is dielectric constant of the material.. 7.

(17) Table 2.1: Comparison of Maxwell’s equations in different unit systems. System. MKS. Gaussian. Heaviside-Lorentz. 107. ε0. 4πc2. 1. 1. µ0. 4π × 10−7. 1. 1. D. D = ε0 E + P. D = E + 4πP. D=E+P. H. H=. 1 B µ0. H = B − 4πM. H=B−M. ∇ · D = 4πρ. ∇·D=ρ. −M. ∇·D=ρ Maxwell’s. ∇×H=J+. Equations. ∇×E+. Lorentz Force. 2.4. ∂B ∂t. ∂D ∂t. =0. ∇×H=. 4π J c. ∇×E+. +. 1 ∂B c ∂t. ∇·B=0. ∇·B=0. E+v×B. E+. v c. 1 ∂D c ∂t. =0. ∇ × H = 1c (J + ∇×E+. ×B. 1 ∂B c ∂t. ∂D ) ∂t. =0. ∇·B=0 E+. v c. ×B. Sullivan’s Implement for Perfectly Matched Layer Absorbing Boundary Conditions. The size of domain that can be simulated using the FDTD method is limited by the computer resources. And the vector field components at the lattice truncation planes cannot be computed using the centred-differencing approach discussed earlier because of the absence of known field data at points outside of the lattice truncation. In our problem the boundary condition must be consistent with Maxwell’s equations in that an outgoing scattered wave must exit the simulated domain without nonphysical reflection, to simulate extending the finite size domain to infinity, just as if the lattice truncation is invisible. Therefore a suitable absorbing boundary condition(ABC) of the computation domain must be employed. One of the most efficient ABCs is the perfectly matched layer(PML) developed by Berenger[14], who employed a fictitious, directionally dependent pair of electric and magnetic conductivities for the purpose of absorbing outgoing waves and minimizing the reflection back into the problem domain. In 1996, Sullivan proposed an easier, and more efficient method to implement the PML ABC[15]. The initial implementation of PML required the E and H fields to be split. By the introduction and delicate arrangement of the g’s and f ’s coefficients in his paper, Sullivan used an unsplit PML scheme in the three-dimensional(3D) FDTD method.. 8.

(18) In addition, he used displacement fields instead of electric fields in the traditional FDTD iteration. This required one more step during the iteration to calculate electric fields from the displacement fields, thus more time of calculation. The following shows the update expressions of Dx and Hx in x-direction for example: -Dx case: h 1 1 1 1 n− 1 n− 1 curl h = Hz 2 (i + , j + , k) − Hz 2 (i + , j − , k) 2 2 2 2 1 1 1 1 i n− 12 n− 12 − Hy (i + , j, k + ) + Hy (i + , j, k − ) , 2 2 2 2 1 1 1 n−1 n (i + , j, k) + 0.5 × gi1(i + ) × curl h, ID (i + , j, k) =ID x x 2 2 2 1 1 1 n Dxn (i + , j, k) =Dxn−1 (i + , j, k) + 0.5 × (curl h + ID (i + , j, k)) x 2 2 2. (2.13). -Hx case: h 1 1 curl e = Eyn (i, j + , k + 1) − Eyn (i, j + , k) 2 2 1 i 1 − Ezn (i, j + 1, k + ) + Ezn (i, j, k + ) , 2 2 1 1 1 1 n+ 21 n− 12 IHx (i, j + , k + ) =IHx (i, j + , k + ) + 0.5 × f i1(i) × curl e, 2 2 2 2 1 1 1 1 n− 21 n+ 21 Hx (i, j + , k + ) =Hx (i, j + , k + ) 2 2 2 2 1 1 n+ 1 + 0.5 × (curl e + Ix 2 (i, j + , k + )) 2 2. (2.14). For good absorption, Sullivan’s g’s and f ’s coefficients in the PML region are verified experimentally to be[16]: f i1(i) = xn(i), 1 , gi2(i) = 1 + xn(i) 1 − xn(i) , gi3(i) = 1 + xn(i) where xn(i) = 13 ×. 3 i length pml. (2.15) (2.16) (2.17). for i = 1, 2, . . . , length pml and length pml is the number. of layers of PML as it goes from the PML region in x-direction into main problem domain. Throughout the main problem domain, f i1 is zero, gi2 and gi3 are 1, thus the update expressions reduce to the typical update expressions (2.9) and (2.10) in Yee’s algorithm. gi1, f i2, and f i3 are similar but with some slight differences, only in that they are computed at the half intervals (i + 21 ).. 9.

(19) 2.5. My Implement for PML ABC. Sullivan used equal grid sizes in x-, y-, and z-directions and MKS unit system in his paper. I’ve made some modifications in his code to adapt the conditions of unequal grid sizes and Heaviside-Lorentz unit system. In addition, I’ve made a simplification in my simulation. The electric fields are directly calculated without calculating displacement fields first. The followings are sample code for implements for Ex and Hx for all directions. -Ex case: n− 12 n− 1 Hz (i + 21 , j + 12 , k) − Hz 2 (i + 21 , j − 12 , k) , dif hz = ∆y n− 1 n− 1 − Hy 2 (i + 21 , j, k + 12 ) + Hy 2 (i + 12 , j, k − 21 ) dif hy = , ∆z 1 1 1 n−1 n IEx,z (i + , j, k) + dt × gi1(i + ) × dif hz, (i + , j, k) =IEx,z 2 2 2 1 1 1 n−1 n IEx,y (i + , j, k) =IEx,y (i + , j, k) + dt × gi1(i + ) × dif hy, 2 2 2 1 1 n n−1 Ex (i + , j, k) =gj3(j) × gk3(k) × Ex (i + , j, k) 2 2 dt +gj2(j) × gk2(k) × εrx (i + 21 , j, k) i h 1 1 n n × dif hz + dif hy + IEx,z (i + , j, k) + IEx,y (i + , j, k) . 2 2. (2.18). -Hx case: dif dif 1 n+ 1 IHx,y2 (i, j + , k + 2 1 n+ 21 IHx,z (i, j + , k + 2 1 n+ 12 Hx (i, j + , k + 2. n Ey (i, j + 12 , k + 1) − Eyn (i, j + 12 , k) ey = , dz − Ezn (i, j + 1, k + 12 ) + Ezn (i, j, k + 21 ) , ez = dy 1 1 1 n− 1 ) =IHx,y2 (i, j + , k + ) + dt × f i1(i) × dif ey, 2 2 2 1 1 1 n− 21 ) =IHx,z (i, j + , k + ) + dt × f i1(i) × dif ez, 2 2 2 1 1 1 n− 1 ) =f j3(j) × f k3(k) × Hx 2 (i, j + , k + ) 2 2 2 h +f j2(j) × f k2(k) × dt × dif ey + dif ez. 1 1 1 1 n+ 1 n+ 1 +IHx,y2 (i, j + , k + ) + IHx,z2 (i, j + , k + )). 2 2 2 2. 10. (2.19).

(20) 2.6. Bloch’s Boundary Conditions for Periodic Structures. Bloch’s theorem tells us that if the electromagnetic fields are in a dielectric structure which is infinitely extended periodic in some direction, we can express the fields as the compositions of plane waves propagating along that direction and modulated by a periodic function in that direction. Figure 2.2 shows an illustration of the implement of Bloch’s boundary conditions in the FDTD method. The structure is infinitely extended periodic in x- and y-directions with unit cell of length Lx and width Ly . The specified wave vector is k = (kx , ky ). For +x direction, the field leaves the right boundary and re-enter from the left boundary with a phase shift of e−ikx Lx . The phase shift is e+ikx Lx in the opposite direction. The argument for y-direction is similar. My fyq).jlyMy*. fyq).jlzMz*. fyq)jlzMz*. Mz. fyq)jlyMy*. Figure 2.2: Bloch’s boundary conditions.. Bloch’s boundary conditions are usually used in calculating the band structures of photonic crystals. More details will be discussed in Section 4.1.3.. 2.7. Symmetric Boundary Conditions. Symmetric boundary conditions are used when one is interested in a problem that exhibits one or more planes of symmetry. For the case of symmetry, the reflected tan-. 11.

(21) gential components of electric fields with respect to the plane of symmetry keep the same and the reflected normal components of electric fields are reversed. Whereas the reflected tangential components of magnetic fields with respect to the plane of symmetry are reversed and the reflected normal components keep the same. In the case of asymmetric boundary conditions, the situations of reflected electric fields and magnetic fields are just exchanged. As a conclusion, symmetric boundary behaves like a mirror for the electric field, and anti-mirror for the magnetic field. On the contrary, an asymmetric boundary behaves like a mirror for the magnetic field, and anti-mirror for the electric field. A visual explanation of a symmetric boundary condition is shown in Figure 2.3 below. Careful consideration must be given to whether symmetric or asymmetric boundary conditions are requested, given the vector symmetry of the desired solution. More symmetry analyses will be addressed in the next chapter.. symmetric boundary. asymmetric boundary. Magnetic Field Components. Figure 2.3: Symmetric boundary conditions.. In the FDTD method, applying one symmetric/asymmetric boundary condition reduces the size of computation domain to half of the origin, thus half the time of simulating the whole domain is saved. In Yee’s algorithm, at most three symmetry/asymmetry boundary conditions can be applied because of the orthogonal lattice. Therefore we may improve the speed of calculation at most eight times faster in some situation. To calculate the tangential components of electric fields in the plane of symmetry, the field data at points just outside of the plane of symmetry are needed for the finite. 12.

(22) difference scheme. The data can be derived from the symmetry property of the data at points just inside of the plane of symmetry. In the plane of asymmetry, the tangential components of electric fields are set to zero directly, just as a perfect electric conductor boundary. For example, if we take the plane of symmetry to be the z-plane in the middle of the whole problem domain, the computation domain is restricted to be the lower part of the domain. For the symmetric boundary, the first term in right hand side of the n− 21. second procedure in the sample code of Ex in previous section, expressed as −Hy 1 , j, nk 2 2. (i +. + 21 ) where nk is the maximum index at the symmetry boundary in +z direction,. is located outside the symmetric boundary. According to the principles of symmetry, it n− 21. should be transferred to Hy n− 1 2. 2Hy. − 12 ) (i+ 21 ,j, nk 2 , ∆z. 2.8. (i + 21 , j, nk − 12 ). Finally the modified dif hy term becomes 2. depending on only the data points inside the computation domain.. Poynting Vector and EM-Energy in Yee’s Lattice. Sometimes it is necessary to study the energy flows and electromagnetic energies stored in the structure. But the vector field components in Yee’s lattice aren’t positioned at same locations so it needs some tricks to implement the electromagnetic energy and Poynting vector in Yee’s lattice. First of all we define two operators. Let A be any component of electric fields and magnetic fields, the differencing operator dˆξ and averaging operator m ˆ ξ are defined as follows: A(. . . , ξ + 21 ∆ξ, . . .) − A(. . . , ξ − 12 ∆ξ, . . .) , dˆξ [A(. . . , ξ, . . .)] = ∆ξ A(. . . , ξ + 12 ∆ξ, . . .) + A(. . . , ξ − 12 ∆ξ, . . .) , m ˆ ξ [A(· · · , ξ, · · · )] = 2. (2.20) (2.21). with ξ is defined for x, y, z, t. The discrete Poynting vector now is written as [17]:  m ˆ [m ˆ (E )m ˆ (H )] − m ˆ z [m ˆ t (Hy )m ˆ x (Ez )]  y x y t z  S =E×H = m ˆ [m ˆ (E )m ˆ (H )] − m ˆ x [m ˆ t (Hz )m ˆ y (Ex )]  z y z t x m ˆ x [m ˆ z (Ex )m ˆ t (Hy )] − m ˆ y [m ˆ t (Hx )m ˆ z (Ey )]. 13.     . (2.22).

(23) Note that the spatial and temporal allocations of the individual field components of Poynting vectors are the same as those of electric fields in Yee’s algorithm. The electric field and the magnetic field energy density are expressed as: ε|E|2 = m ˆ x (εExn 2 ) + m ˆ y (εEyn 2 ) + m ˆ z (εEzn 2 ) n− 12. µ|H|2 = m ˆ ym ˆ z (µHx. n+ 12. Hx. n− 12. )+m ˆ xm ˆ z (µHy. n+ 12. Hy. (2.23) n− 12. )+m ˆ xm ˆ y (µHz. n+ 21. Hz. ).. (2.24). In continuous world, the power radiated from a close volume V can be expressed as the R H integral of Poynting vector on the boundaries of V, said S: Prad = V (∇ · S)dτ = S S · da. The radiated power from a box in Yee’s lattice, in analogy to real world, is expressed as a summation: Prad. ∆x∆y∆z = 2. X j1 k1 i1 X X. ˆ D·S ,. (2.25). i=i0 j=j0 k=k0. ˆ = (dˆx , dˆy , dˆz ) is the discretized divergence operator, i0 and i1 are the indices where D of lower and upper bounds of the box in x-direction, respectively, j0 , j1 , k0 , and k1 are similar for y− and z−directions. Simplifying this we get a much easier and faster way for computing energy flows in to or out of this box. Prad. j1 k1 ∆y∆z X X i=i1 + 1 Sx |i=i − 21 = 0 2 2 j=j k=k. +. 0. 0. k1 X. i1 X. 0. i=i0. ∆z∆x 2 k=k. j=j1 + 12. Sy |j=j. 1 0− 2. j1 i1 X ∆x∆y X k=k1 + 1 Sz |k=k − 21 + 0 2 2 i=i j=j 0. 0. 2.9. (2.26). Time Step Issues. There are some issues of the time step ∆t. First, ∆t must satisfy the inequality of Courant’s stability condition[18]: cmax ∆t ≤ q. 1 1 ∆x2. +. 1 ∆y 2. +. 1 ∆z 2. ,. (2.27). where cmax is the maximum electromagnetic wave phase velocity within the media being modeled. Here cmax is 1 in Heaviside-Lorentz unit system. For the purely TE and TM cases, it can be shown that the modified time-step limit for numerical stability is obtained simply by setting ∆z = ∞. Second, ∆t is relevant to the Nyquist frequency. The sampling. 14.

(24) time for digital filter process is equal to the time step ∆t in FDTD algorithm, therefore the corresponding sampling frequency is. c0 ∆t. =. 1 ∆t. Hz, since c0 is assumed to be 1 in Heaviside-. Lorentz unit system. The Nyquist frequency is then. 1 Hz 2∆t. by the sampling theorem,. which represents the highest frequency that the data can produce without the effect of aliasing. So if you want to get the information at higher frequencies, you will needed the smaller time step interval. Finally, the more time steps going, the finer resolution in the frequency domain. This can be derived directly from the basic theorem of discrete Fourier transform.. 15.

(25) Chapter 3 Symmetry Analysis Physical systems exhibit intrinsic symmetries which can be used to simplify the solution of the equations governing these systems. Group theory is a mathematical tool to investigate the classification of the solutions within the symmetry. In this chapter we will focus on the C6v group as an example which is the symmetry group of 2D photonic crystals of hexagonal lattice.. 3.1. Some Terminologies and Useful Theorems. In this section I will give some basic knowledge of group theory and some related useful theorems that are needed to analyze the symmetry property of photonic crystals. Instead of rigorous definitions and theorems, I shall give some loose definitions and theorems of group theory with examples to help us understanding the language of group theory. People who are interested in the details of group theory may consult books of group theory[19].. 3.1.1. Basic Concepts of Group. A set G of elements A, B, C, . . . is called a group if there exists an operation which associates any ordered pair in G with a third element in G. This operation, often called multiplication or product, and marked as “·”, satisfies the following requirements: 1. The associative law: A · (B · C) = (A · B) · C. 2. There exists an identity element E so that the product of E and any element in G, said A, is A itself.. 16.

(26) 3. For every element A in G, there exists an inverse element A−1 so that the product of A and A−1 is identity. An element B in G is said to be conjugate to A if there exists a group element R such that B = RAR−1 . A class is a collection of mutually conjugate elements of a group. The number of elements in G is the order of the G. A group is a finite group if this group has a finite number of elements. A physically important example of a finite group is the set of symmetry operations which are the covering operations of a symmetrical object. By a covering operation, we mean a transformation which would bring the object into a form indistinguishable from the origin one. In a crystalline structure, the covering operations are consisted of translations of multiples of integer of lattice vector(which is called primitive translations), rotations, reflections, and compositions of them. The complete set of covering operations is called the space group of the crystal. The group of operations which is obtained by setting all translations in the space group elements equal to zero is called the point group of the crystal. It has been shown that there are only 32 point groups in a crystal system[19]. The symmetry group of hexagonal lattice is the point group C6v = {E, C6 , C6−1 , C3 , C3−1 , C2 , σx , σx0 , σx00 , σy , σy0 , σy00 }. The symmetry operations are shown in Figure 3.1. The classes of C6v = {(E), (C6 , C6−1 ), (C3 , C3−1 ), (σx , σx0 , σx00 ), (C2 ), (σy , σy0 , σy00 )}, or written by C6v = {E, 2C6 , 2C3 , C2 , 3σx , 3σy }. sx. sy'. sy'' sx'. C2. C3 C6. sy. sx''. Figure 3.1: Symmetry operations for two-dimensional hexagonal lattice.. In momentum space the first Brillouin zone of a hexagonal lattice, as shown in Figure 3.2, is a hexagon. Therefore the point group in momentum space, denoted by M, is again C6v . The star of k (?k) is the set of wave vectors generated by applying all the operators. 17.

(27) of M to k. When k touches the zone boundary, two wave vectors of ?k may be equivalent if there is a displacement of a reciprocal lattice vector G between them. They must be treated as identical. For example, kK1 , kK3 and kK5 are treated as identical one because there is a difference of a reciprocal lattice between them. And the set of wave vectors kK2 , kK4 , and kK6 are similar. So ?kK is {kK1 , kK2 }. ky K6. M1. K1. M6. M2. K5. K2. G. M5. kx. M3 K4. M4. K3. Figure 3.2: First Brillouin zone of a hexagonal lattice. K stands for the set of all Ki and M stands for the set of all Mi for i = 1 · · · 6.. The point group of a given wave vector k, denoted by Mk , is the subgroup of the point group of M that consists of all the rotations of M that rotate k into itself or its equivalent vectors. We can easily verify that MkM is C2v and MkK is C3v .. 3.1.2. Representation of a Group. Two groups are isomorphic when an unique, one-to-one correspondence exists between their elements in such a way that products correspond to products. A (faithful) representation of a group is a matrix group to which the group to be represented is isomorphic. Thus, it consists of the assignment of a matrix Γ to each group element in such a way that Γ(R1 )Γ(R2 ) = Γ(R1 R2 ) holds true for all matrices Γ. The number of rows and columns in a representation matrix is called the dimension of the representation. Similarity transformation of Γ is to multiply a non-singular matrix S in front of Γ and S −1 after Γ. Γ1 and Γ2 are two representations   of G with dimensions l1 and l2 . A new representation can be generated. Γ = . Γ1 0. 0. 2. Γ.  The representation Γ is called the direct sum of the. 18.

(28) representations Γ1 and Γ2 , expressed as Γ = Γ1 ⊕ Γ2 . Representations which can be transformed into direct sums by similarity transformations are called reducible. If such a transformation is not possible, the representation is called irreducible. Each irreducible representation has its own spatial symmetry which is expressed by its character. We define the character of the representation Γ of group element R as the trace of the matrix Γ(R), indicated by χ(R). It is interesting to know that the character of the representation will be invariant during a similarity transformation. It is convenient to display the characters of the various representations in a character table for any given group. The columns are labeled by the various classes, and the number of elements in the class. The rows are labeled by the irreducible representations. The entries in the table are the characters of corresponding representations and classes. The character table of C6v is shown in Table 3.1. Table 3.1: Character table for the C6v . C6v. E. 2C6. 2C3. C2. 3σy. 3σx. A1. 1. 1. 1. 1. 1. 1. A2. 1. 1. 1. 1. -1. -1. B1. 1. -1. 1. -1. 1. -1. B2. 1. -1. 1. -1. -1. 1. E1. 2. 1. -1. -2. 0. 0. E2. 2. -1. -1. 2. 0. 0. Here are some important theorems for the characters of representations: • The number of inequivalent irreducible representations of a group G is equal to the number of classes of G. • The characters of group elements in the same class are equal. • The dimension of the representations is equal to the character of the identity element. • The sum of the squares of the dimensions of all irreducible representations is equal to the order of the group. • Two representations are equivalent if their character systems are equivalent.. 19.

(29) P. • A representation is irreducible if. |χ(R)|2 is equal to the order of the group G.. R∈G. • Two orthogonality relations: X. χi (R)∗ χj (R) = gδij ,. (3.1). gδij , Ni. (3.2). R∈G. X. χk (Ti )∗ χk (Tj ) =. k. where χi stands for the character of Γi , g is the order of the group, and Ni is the number of elements conjugate to Ti . For any reducible representation there is an unique way to reduce it into direct sum of irreducible representations. So if you can find a combination that works it is right. Here is a simple way to check the multiplicity mi (Γ) of the irreducible representation Γi contained in a reducible representation Γ. mi (Γ) =. 1X χ(R)χi (R)∗ , g R∈G. (3.3). where χ stands for the character of Γ, χi for the character of Γi , and g for the order of the group.. 3.1.3. Basis Functions and Projection Operators. Now we define the basis functions of the ith irreducible representation Γi for an arbitrary scalar function φ(r)(φ for simplicity). Two labels are needed for a basis function. One is for the irreducible representation and the other is for the row within the representation. Let a basis function belonging to the nth row of Γi be denoted by φin . The other functions φim required to complete the basis for the representation are called the partners of the given function. The acting of a symmetry operation, said R, on φ can be expressed as PˆR φ(r) = φ(R−1 r). When R acts on a vector field Φ(r), both the field vector and the argument are altered accorading to: PˆR Φ(r) = RΦ(R−1 r). By this definition the result of operating with an element R of the group on φin can be expressed as a linear combination of φin and its partners as follows: PˆR φin =. li X. φim Γi (R)mn ,. m=1. where li is the dimension of the representation.. 20. (3.4).

(30) φ may be expanded as a sum of basis functions: φ =. li PP. φin . We should note that. i n=1. this is not an expansion like the Fourier series. We don’t expand a function with respect to an orthonormal complete set of functions here. The set of functions φim depends on φ itself. They can be found by means of the so-called projection operators. i Pmn =. li X i Γ (R)∗ PˆR g R∈G mn. i Pnn φ = φin. (3.5) (3.6). i where g is the order of group G. The partners of φin are then obtained by φim = Pmn φin . i Also note that the operator Pmn requires the information of the full representation. matrices. There is an easier projection operator, denoted by P i , which brings less-detailed results by using the information of only the characters of the representation. P i carries φ into the part of the ith representation. Pi =. X. i Pnn =. n. li X i χ (R)∗ PˆR g R∈G. P i φ = φi. (3.7) (3.8). As a trivial example, consider the group C1h , consisting of the identity and the reflection operator σ which takes x into −x. This group has two classes and has two one-dimensional irreducible representations. The character table is Table 3.2. Table 3.2: Character table for the C1h . C1h. E. σ. A. 1. 1. B. 1. -1. Hence, the projection operator of Γ1 is P 1 = 21 (Eˆ + σ ˆ ), and that of Γ2 is P 2 = 12 (Eˆ − σ ˆ ). Operating on an arbitrary function φ(x) yields the basis functions P 1 φ(x) = 21 (φ(x) + φ(−x)) and P 2 φ(x) =. 1 (φ(x) 2. − φ(−x)). It reflects the fact that any function can be. expressed as the sum of odd and even functions constructed as above and that any odd function is orthogonal to any even function.. 21.

(31) 3.2. Classification of Defect Modes in Hexagonal Lattice. 3.2.1. Symmetry Group of Photonic Crystals. The eigenvalue equations of Maxwell’s equations in photonic crystals are given by[20]: 1 ω2 ∇ × H(r) = 2 H(r), εr (r) c ω2 1 ∇ × ∇ × E(r) = 2 E(r), LE E(r) ≡ εr (r) c. LH H(r) ≡ ∇ ×. (3.9) (3.10). Analogue to the case in quantum mechanics, the group of these equations is formed by the space group of the photonic crystal which transforms εr (r) into itself when the group of Schrödinger’s equation is determined by the potential function V (r). One can prove that any symmetry operator belonging to the point group of 3D photonic crystals commutes with LE and LH . This exhibits the ability to classify the modes of equations (3.10) Ekn (r) and Hkn (r) according to the irreducible representations of the group Mk in momentum space, just as the treatment to Schrödinger’s equation in quantum mechanics. However, there is something important to be noted. The symmetry of the magnetic field and that of the electric field are generally different from each other, as the former is an axial vector field whereas the latter is a polar vector field. The character for magnetic field differs from that of electric field by a multiple of the determinant of the operation. The determinant is +1 when the operation is proper transformation, wheras a determinant of −1 is derived for improper transformations. That means the characters of improper transformations for magnetic fields need to change sign from the origin. Electromagnetic fields can be classified into pure TE or TM modes when the structure has mirror symmetry and continuous translational symmetry in z-direction. In other words, the dielectric constant of the structure remains constant and extends to infinity in z direction. Now the eigenvalue problems for TE and TM modes in 2D photonic crystals are given by[20]:. 22.

(32) ∂ ω2 ∂ 1 ∂ 1 ∂ ≡− Hz (rk ) = 2 Hz (rk ), + ∂x εr (rk ) ∂x ∂y εr (rk ) ∂y c 2 2 2 ∂ ω 1 ∂ (2) LE Ez (rk ) ≡ − + 2 Ez (rk ) = 2 Ez (rk ), 2 εr (rk ) ∂x ∂y c. (2) LH Hz (rk ). (3.11) (3.12). respectively, where rk denotes the in-plane position vector (x, y). (2). (2). Like the 3D case of photonic crystals, we should note that LH and LE commute with the 2D symmetry operations belonging to the point group of the 2D photonic crystal. Therefore the electromagnetic fields can be classified according the irreducible representations of the k-group Mk . Since no 2D symmetry operation changes Ez or Hz , we can observe the symmetry properties of the modes just from the scalar field Ez for the TM polarization and the scalar field Hz for the TE polarization for 2D photonic crystals for convenience.. 3.2.2. Symmetry of 2D Photonic Crystals of Finite Thickness. In actual case the 2D photonic crystal slabs are surrounded by air in the z direction. The symmetry group of this structure is D6v , which is a direct product of C6c and C1h point groups: D6h = C6v × C1h where C1h consists of the identity operation and the mirror reflection by the middle plane of the slab, σ ˆz . The character table of C1h is shown in Table 3.2. Thus the middle plane of the slab can be taken as a symmetry boundary for σz = +1 to save the computation time. On the contrary, it can be taken as an asymmetry boundary for σz = −1. The modes for σz = +1 are the so-called TE-like modes. And those for σz = −1 are TM-like modes. The dominant components of electromagnetic fields of TE-like modes are Ex , Ey , and Hz while other components are small and close to zero for a slab thickness around λ/2. Here we take the λ to be around the mid-gap. The situation of TM-like modes is similar, the dominant components are Hx , Hy , and Ez . From now on we will focus on the discussion of symmetry to TE-like modes, which is pure TE mode and has only magnetic component Hz in the middle plane of the slab.. 23.

(33) 3.2.3. Symmetry of Defect Modes. A small variation of dielectric constant in a small region within the photonic crystals can be seen as a perturbation of modes of photonic crystals without defects. Enlarging a hole creates acceptor modes in photonic crystal band gap, which is the analogue to the acceptors in electronic states in a crystal. Similarly, reducing the radius of hole creates donor modes in analogy to donors in electronic states in a crystal[21]. The easiest way to create donor modes is to remove single air hole from array of holes of hexagonal lattice arrangement. For the coordinate system in which origin is at the defect center, the space group and the point group are still C6v as before. The defect modes are attributed to the irreducible representations of the C6v point group. There are four one-dimensional representations and two two-dimensional representations for the C6v point group. The modes at photonic band edges are used as a symmetry basis to generate approximated forms of defect modes[4]. The minimum in air band occurs at M -points. Therefore donor modes can be approximated by the unperturbed Bloch modes in air band edges. On the contrary, acceptor mode can be approximated by the unperturbed Bloch modes in dielectric band edges. Let us find the symmetry basis of unperturbed Bloch modes at the band edges first. A symmetry basis for the modes of the photonic crystal at the M -point can be found by projecting the seed Bloch modes B?kM to the representations of the group of the wave vector MkM (=C2v ). The BAM2 corresponds to the dielectric band mode and BBM1 to the air band mode. A set of basis functions for modes in air band edges can be found:   sin(kM1 · rk )     CBBM1 = zˆ  sin(kM2 · rk )  ,   sin(kM3 · rk ). (3.13). where the superscript stands for the wave vectors and the subscript stands for the representations of the corresponding group of those wave vectors. Similarly, the basis functions for modes in the dielectric band edges should be:   −ikK1 ·rk −ikK3 ·rk −ikK5 ·rk e + e + e  V BAK02 = zˆ  −ikK2 ·rk −ikK4 ·rk −ikK6 ·rk e +e +e. 24. (3.14).

(34) M basis under the operations of C6v can be written as The representation of the CBB1. the following six  1   ΓE1 (E) =  0  0. matrices:   0 0 −1 0 0     1 0  , ΓE1 (C2 ) =  0 −1 0   0 1 0 0 −1. . . . . . 0. 0. 1. .        , ΓE1 (C3 ) =  −1 0 0  ,    0 −1 0 . . −1 0 0 1 0 0 0 1 0           ΓE1 (C6 ) =  0 0 1  , ΓE1 (σy ) =  0 0 −1  , ΓE1 (σx ) =  0 0 1      0 1 0 0 −1 0 −1 0 0.    . . This representation equals to E1 ⊕ B100 . Using the projection operators on CBBM1 , a set of basis functions for E1 and B100 appears as follows: BBd1100 = zˆ(sin(kM1 · rk ) − sin(kM2 · rk ) + sin(kM3 · rk )),. (3.15). BEd11 ,1 = zˆ(2 sin(kM1 · rk ) + sin(kM2 · rk ) − sin(kM3 · rk )),. (3.16). BEd11 ,2 = zˆ(sin(kM2 · rk ) + sin(kM3 · rk )).. (3.17). The BBd100 mode transforms like a hexapole, whereas the doubly degenerate modes BEd1 1. transform as an (x, y)-dipole pair.. 3.3. Some Design Rules. In real space the design of cavities of high quality factors lies in maximalizing the overlap of the defect mode and the dielectric region to obtain the maximal optical gain. In actual case, the quality factor is limited by radiation loss in the vertical direction while the horizontal confinement can be greatly enhanced as the number of surrounding layers of air holes increases. In momentum space, one way to reduce the vertical radiation loss is to minimize the power of the electromagnetic fields within the light cone of the slab. The light cone is defined as ω ≥ c|kk |, as shown in Figure 3.3. The fields inside the light cone are radiative, whereas those outside the light cone are guided within the slab. If the slab thickness is too large, the structure has multi-mode in the vertical direction which is undesirable. If the slab is too thin, the mode is not confined well within it vertically. When the slab thickness increases, we would expect the increase of Q⊥ and the. 25.

(35) ω. light cone. guided mode region. radiation mode region. ω = ω0 ky. kx. Figure 3.3: light cone. decrease of Qk . The increase of Q⊥ is due to the lowering of frequencies of defect modes in a thicker slab causing the shrinking of the size of the light cone. The decrease of Qk is predicted because the in-plane bandgap shrinks in size as the slab thickness increases. Thus the control of slab thickness is important. A typical value of slab thickness of around λ/2 is usually used, where λ is the wavelength around the mid-gap in the material. There is a tradeoff between the mode volume and the quality factor. The localization of the mode in real space causes a spread of the mode in momentum space and a portion of the mode intersecting the light cone causes radiation loss vertically. We can improve this by utilizing larger cavities like seven missing holes or nineteen missing holes. But doing this may complicate the modes in the cavity.. 26.

(36) Chapter 4 Simulation Results 4.1. Photonic Band Structures. Band structures of photonic crystals are dispersive relations of photonic crystals which are often simulated by the PWE method because of their periodic dielectric structures. Besides, the FDTD method can be used to calculate the photonic band structures, which is widely used to investigate the scattering of waves in specific structures. In this section we exhibit the results of simulation of photonic band structures using both PWE and FDTD methods. A brief introduction of the PWE method can be found in the excellent book written by Sakoda[20].. 4.1.1. 2D PWE Method. 2D PWE method is used to simulate those structures which are periodic in two directions. Band gaps of TE modes are preferred by the virtue of etching periodic air holes on materials[2], so let us focus on the 2D PWE simulations of TE modes. Figure 4.1 shows a typical band structure of hexagonal lattice. The structure parameters are as follows: the reflective index of the material is 3.4, the hole radius is 0.3a, and slab thickness is 0.4a, where a is the lattice constant. The effective index of fundamental TE mode of this slab is evaluated to be 2.65 by the method of finite difference scheme[22]. Mode profiles can also be derived from the PWE method. Figure 4.2 shows the mode profiles at first band edge at K-point and M-point.. 27.

(37) 0.8. Frequency (ωa/2πc=a/λ). 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0.0. K. Γ. M. Γ. Figure 4.1: The TE band structure by 2D PWE method.. 5. 1.6. 5. 10. 1.4. 10. 15. 25. 0.8. 30. 0.6. 35. 40. 0.4. 40. 10. 15. 20. −0.5 −1. 45. 0.2 5. 0. 30. 35. 45. 0.5. 20. 1. 25. 1. 15. 1.2. 20. 1.5. 25. −1.5 5. (a) K-point. 10. 15. 20. 25. (b) M-point. Figure 4.2: Hz mode profiles at band edges.. 4.1.2. 3D PWE Method. 3D PWE method is typically used for the analysis of photonic band structures of 3D photonic crystals. However, it can also be used to calculate 2D photonic crystal slabs via some techniques. When defining the dimension of supercell in the direction perpendicular to the slab, a margin of air of at least twice broader than the slab thickness is needed above and below the slab. It benefits from that the coupling between the modes of adjacent slabs is negligible. Figure 4.3 shows the band structure of the identical structure as in Section 4.1.1. Note that we have chosen the symmetry plane to be at the middle of the slab. Thus it is a band structure of TE-like modes of the slab.. 28.

(38) 0.7. Frequency (ωa/2πc=a/λ). 0.6 0.5 0.4 0.3 0.2 0.1 0.0. Γ. K. M. Γ. Figure 4.3: The TE-like band structure by 3D PWE.. 4.1.3. 3D FDTD Method. In Yee’s lattice, it is difficult to form the unit cell of hexagonal lattice which is a rhombus in orthogonal grids. Thus we use a rectangular domain containing two unit cells shown in Figure 4.4: We assume the in-plane periodicity extends to infinity for the convenience to apply Bloch’s boundary conditions with a specific wave vector at all four sides of the domain. The procedure of Bloch’s boundary conditions is discussed in Section 2.6. Open boundary conditions such as PML ABCs are applied to the upper and lower boundaries in the. pha. se dual shi ft. directions outward and inward to this paper.. source monitor. Figure 4.4: Double unit cell and the locations of sources and monitors.. 29.

(39) A short impulse is used as the excitation source. One should note that the source should have a dual source at the other unit cell in the domain. And there is a phase shift of e−ik·a between them, where a is the lattice vector of unit cell. To avoid neglecting all possible modes, time monitors are randomly chosen at three positions to record the time signals. Then we collect the three time signals and add them as a new signal. Normalize this signal with maximal amplitude. Apodize them using Gaussian function to exclude the undesired frequencies. Chirp z-transform the signal with a range of frequency of interest to observe the frequency response. To obtain the band structure we have to repeat the procedure of Bloch’s boundary conditions with wave vectors along the edges of irreducible Brillouin zone in k-space. As an example, the band structure of a hexagonal lattice is shown in Figure 4.5. The TE-like photonic band gap is ranged from a/λ = 0.288 to a/λ = 0.364. Note that the missing points may result from the abandoning of small signals of 1/1000 smaller than the maximum peak. 0.5 0.45 0.4. 0.364 Photonic Band Gap. 0.3. 0.288. 0.25. ligh. 0.2. ne. t co. ne. t co. ligh. Normalized frequency (a/λ). 0.35. 0.15 0.1 0.05 0. 0 Γ. 19 M. 38 K. Γ. Figure 4.5: The band structure by the FDTD method. (r/a = 0.3, d/a = 0.4, n = 3.4). 30.

(40) 4.2. Defect Systems. In this section, the refractive index of the material is fixed to be 3.4. This implies that here we are not concerned with the dispersive media and lossy materials so that n is taken to be a frequency-independent real constant. The ratio between the slab thickness and the lattice constant d/a is chosen to be 0.4. The defect is formed by removing a single air hole from array of holes of hexagonal lattice arrangement on a dielectric slab. The number of layers surrounding the defect is five. We have also designed some new structures by gradually increasing or decreasing the hole radii as the surrounding layers go outward the defect. Figure 4.6 shows these designs. The ranges of photonic band gaps of photonic crystal slabs of hexagonal lattice without defects vs different r/a are listed in table 4.1 for references of our structures. The calculation procedure has been described in previous sections. Table 4.1: Comparison of photonic band gaps using different methods. The frequencies are taken to be normalized (a/λ). method. r/a. 2D PWE. 0.28. 0.30. air band edge. 0.303 0.318. 0.336. 0.356 0.378. (nef f = 2.65). dielectric band edge. 0.262 0.265. 0.269. 0.275 0.282. 3D PWE. air band edge. 0.330 0.335. 0.337. 0.340 0.343. dielectric band edge. 0.283 0.286. 0.291. 0.297 0.304. air band edge. 0.330 0.346. 0.364. 0.385 0.408. dielectric band edge. 0.279 0.282. 0.288. 0.293 0.305. 3D FDTD. 4.2.1. 0.26. 0.32. 0.34. Resonant Frequencies of Defect Modes. The PWE method is still useful to determine the resonant frequencies of the defect cavity by the use of so-called supercell scheme, which defines a larger size of unit cell containing the defect region. The spirit of supercell scheme lies in the periodic arrangement of supercell and the enough space apart from each other for small coupling between defect modes. Figure 4.7 shows an example of the band structure of the defect system by 2D and 3D PWE methods. A supercell containing 7 by 7 unit cells and a single defect is used. In 3D PWE method a margin of air of twice broader than the slab thickness is. 31.

(41) (a) Type A. (b) Type B. (c) Type C. (d) Type D. (e) Type E. Figure 4.6: Various patterns by changing hole radii. The radius of holes of most inner layer is 0.3a. (a) r/a gradually decrease by 0.02 per layer outward. (b) r/a gradually decrease by 0.01 per layer outward. (c) typical pattern without changing the hole radii. (d) r/a gradually increase by 0.01 per layer outward. (e) r/a gradually increase by 0.02 per layer outward.. 32.

(42) F re quency ( ωa /2 πc =a /λ). 0.4. 0.3. 0.2. 0.1. 0.0. Γ. K. M. Γ. (a) 2D PWE (nef f = 2.65). F re quency ( ωa /2 πc =a /λ). 0.4. 0.3. 0.2. 0.1. 0.0. Γ. K. M. Γ. (b) 3D PWE. Figure 4.7: The band structures of defect system Type C( ar is fixed to be 0.3) by the 2D and 3D PWE methods. used in the direction perpendicular to the slab. In 2D case, Two degenerate defect modes of normalized frequency around 0.31 in the photonic band gap are observed. In 3D case this frequency is around 0.323. Now adding a defect provides the existence of the defect modes inside the gap though it breaks the translational symmetry. So Bloch’s boundary conditions can no more be applied to obtain the band structure by the FDTD method. Although the band diagrams of defect structures cannot be simulated by the FDTD method with Bloch’s boundary. 33.

(43) 30 0.5 0.4. 25. 0.3 20 Intensity (a.u.). Amplitude. 0.2 0.1 0 −0.1 −0.2. 15. 10. −0.3 −0.4. 5. −0.5 100. 200. 300. 400. 500 600 Time steps. 700. 800. 900. 0. 1000. (a) impulse signal. 0. 0.5. 1. 1.5 2 2.5 Normalized frequency (a/λ). 3. 3.5. 4. (b) spectrum of signal. Figure 4.8: Excitation Signal. (α = 1024 and f0 = 0.015) conditions. There are still some tricks to locate the frequencies of the defect modes inside the photonic band gap. We use an excitation of a magnetic point source located at a low symmetric point in the cavity in order to excite all defect modes in the cavity. The polarization of the source is set to be TE polarization in which we are interested. The temporal distribution is a short Gaussian pulse covering the spectrum of the whole band gap and with a central frequency located at the mid-gap frequency. The source signal can be expressed as Hz (t) = sin(2πf0 t) · e−α(t−t0 ) , where f0 is taken around the mid-gap frequency, t and t0 is in unit of ∆t, and α is the coefficient for apodizing the original sine function exponentially. Figure 4.8(a) shows the signal as a function of time steps and Figure 4.8(b) shows the corresponding spectrum. Parameters used in the calculation are listed in table 4.2. Some data processes are performed to analyze the results. First of all, we record the time evolution of the Hz field at three chosen low symmetry positions inside the cavity for 33792 time steps, to detect all possible defect modes. The positions of the source and. monitor 1. monitor 3. source monitor 2. monitors. Figure 4.9: The locations of source and monitors.. 34.

(44) Table 4.2: Parameter table. parameter value ∆x. 1 √. ∆y. 3/2. ∆z. 1. ∆t. 0.54. size of domain (number of grids). 260 × 260 × 88. lattice constant a. 20. number of PML layers. 8. mesh refinement†. 1000. source shift*. (7, 3, 0). monitor 1 shift*. (5, 7, 0). monitor 2 shift*. (3, −8, 0). monitor 3 shift*. (−1, 1, 0). †Number of small grids for smoothing the distribution of dielectric constant in each Yee’s unit cell. *The shifts are from the center of the structure in unit of the corresponding grid sizes.. the monitors are shown in Figure 4.9. Then we do fast Fourier transformation to the collected time domain signals to get the frequency responses. The results are shown in Table 4.3 and Figure 4.10. Two interesting phenomenons are observed. There exists a shift of frequencies as the hole radii changes. And the donor mode seems to be enhanced as the hole radii increasing per layer outward, whereas the acceptor mode seems to be suppressed. Moreover, the changing of the radii of holes reduces the first peak greatly. To describe the effect explicitly, we have to use the quantity “Q” to measure the “quality” of the modes in the cavity[23]. We’ll discuss this later.. 4.2.2. Mode Profiles. In 2D PWE method, the eigenvectors of the problem correspond to the mode profiles. Some mode profiles at high symmetry points in momentum space are shown in Figure 4.11. The mode profiles can also be generated by the 3D FDTD method. Typically it is convenient to choose a random field as an excitation. The fields not belonging to resonant frequencies will disappear in a short time, and any fields which can oscillate in. 35.

(45) 14. Table 4.3: Normalized peak frequencies. Type. 12. first peak. second peak. Type A. 0.30857. 0.35604. Type B. 0.30970. 0.36169. Type C. 0.31083. 0.35717. 4. Type D. 0.31196. 0.36056. 2. Type E. 0.31309. 0.36282. monitor 1 monitor 2 monitor 3. Intensity (a.u.). 10 8 6. 0 0.2. 0.25. 0.3 Normalized frequency (a/λ). 0.35. 0.4. (a) Type A 25. 3500. monitor 1 monitor 2 monitor 3. 3000. 15. Intensity (a.u.). Intensity (a.u.). 20. 4000. monitor 1 monitor 2 monitor 3. 10. 2500 2000 1500 1000. 5. 500 0 0.2. 0.25. 0.3 Normalized frequency (a/λ). 0.35. 0 0.2. 0.4. 0.25. (b) Type B 70 60. 0.3 Normalized frequency (a/λ). 0.35. 0.4. 0.35. 0.4. (c) Type C 90. monitor 1 monitor 2 monitor 3. 80. monitor 1 monitor 2 monitor 3. 70 60 Intensity (a.u.). Intensity (a.u.). 50 40 30. 50 40 30. 20. 20 10 0 0.2. 10 0.25. 0.3 Normalized frequency (a/λ). 0.35. 0 0.2. 0.4. (d) Type D. 0.25. 0.3 Normalized frequency (a/λ). (e) Type E. Figure 4.10: The corresponding spectrums. the resonator become resonant modes. But it is difficult to remove the undesired modes completely. And the ratio of components of composite modes depends on the choice of initial field. Since Painter et al.[4] showed that the defect modes behave like E1 mode under the symmetry group of C6v . The cavity mode is excited by TE-polarized magnetic point sources(Hz ) narrowly peaked around the desired mode’s frequency and placed in the. 36.

(46) (a) supercell. (b) Γ-point(1). (c) Γ-point(2). (d) K-point(1), (intensity). (e) K-point(2), (intensity). (f) M-point(1). (g) M-point(2). Figure 4.11: Definition of domain of supercell and profiles of doubly degenerate defect modes by 2D PWE method. same symmetry as the mode of interest to avoid exciting the undesired modes. We can further classify these modes according to (σx , σy ) = (−1, +1), and (σx , σy ) = (+1, −1), which can be named as x-dipole mode and y-dipole mode respectively. Thus the point sources of x-dipole mode are located at a shift of magnitude +3∆y and −3∆y. 37.

(47) (a) x-dipole mode profile (Hz ). (b) y-dipole mode profile (Hz ). 0.2. 0.2 8. 7. 7. 6. 0.1. 0.1 6. 5. ky/2π. ky/2π. 5 4. 0. 0. 4. 3. 3. 2. −0.1. −0.1. 2. 1. −0.2 −0.2. −0.1. 0. kx/2π. 0.1. 1. −0.2 −0.2. 0.2. e z |) (c) x-dipole mode in momentum space (|H. −0.1. 0. kx/2π. 0.1. 0.2. e z |) (d) y-dipole mode in momentum space (|H. Figure 4.12: Dipole modes by 3D FDTD method. Light cone is shown in red circle. in y-direction from the center with a phase difference of. π 2. between them to match the. symmetry conditions. For y-dipole mode the point sources are shifted from center by +3∆x and −3∆x in x-direction. And a phase difference of. π 2. is still needed.. Figure 4.12 shows the doubly degenerate defect modes labelled by E1 symmetry and the corresponding Hz fields in momentum space. One should note that the dipole mode profiles in momentum space are consistent with the criterion of symmetry analysis of donor modes equations (3.16) and (3.17). The dominant Fourier component of x-dipole mode is ±kM1 . The second dominant Fourier components of x-dipole mode are ±kM2 and ±kM3 . While the dominant Fourier components of y-dipole mode are just ±kM2 and ±kM3 .. 38.

(48) 4.2.3. Quality Factors. The quality factor(Q) of a defect mode is defined as follows: Q = ω0 U/P , where ω0 is the resonant frequency of the mode, U is the stored electromagnetic energy, and P is the radiated power, both U and P are time-averaged. Q can be further divided into in-plane part Qk and vertical part Q⊥ : Pk 1 P P⊥ 1 1 = = + ≡ , + Q ω0 U ω0 U ω0 U Q⊥ Qk where Qk stands for the degree of confinement in-plane and Q⊥ stands for that in vertical directions. This definition of Q is equivalent to that found by fitting parameters from ω0. U = U0 e− Q t . The quality factor can be derived from the FDTD simulations. The excitation is the same as that described in the previous section. The defect mode is excited using point sources narrowly peaked at the resonant frequency of the mode and located according to the symmetry of the mode to excite the mode of interest. Then we record the electromagnetic energy stored in the cavity and the integral of the Poynting vectors over the surface of the domain as the radiated power. P⊥ is measured at ±2a above and below the slab, where a is the lattice constant. And Pk is measured on the other four sides of the domain. The simulation results are shown in Table 4.4. Table 4.4: The quality factors. Type. Q⊥. Qk. Qtot. Q0 *. Type A. 282. 443. 172. 170. Type B. 276. 1032. 218. 213. Type C. 287. 1325. 236. 231. Type D 303. 1524. 253. 247. Type E. 1671. 269. 264. 321. 0. *Q is derived from fitting parameters. We can observe that Qk increases as the radii increase per layer outward. It can be interpreted that this effect is originated from the increase of averaged. 1 εr (r). as r goes. away from the defect center, which acts as the character of potential function in quantum mechanics, providing the effect in analogy to quantum well structures.. 39.

(49) Chapter 5 Summary First the techniques for the FDTD method including Yee’s algorithm and various boundary conditions used for the simulation for photonic crystals are introduced. Then a brief introduction of symmetry analysis by group theory and design rules of defect modes of photonic crystal nano-cavities are investigated. The modes in single defect in 2D photonic crystal slab are doubly degenerate dipole modes. Numerical simulations of the band structures, mode profiles, and resonant frequencies are obtained by both the PWE method and the FDTD method with various boundary conditions. Some techniques like supercell scheme or domain containing two unit cells are employed. For the defect modes, quality factors are also derived by the FDTD method. Proper excitation is needed. The excitation of point sources narrowly peaked around the frequency of the desired mode are placed according to the symmetry property of the mode. The simulated mode profiles are consistent with the criterion of symmetry analysis. For the design of the cavity for high quality factors, we also try to compare slightly different structures by gradually increasing or decreasing the hole radii per layer outward. Numerical experiments showed that a gradually increase of radii of holes per layer outward may increase the in-plane quality factor of donor modes, just as the case of the enhancement of confining electrons by quantum wells.. 40.

(50) References [1] E. Yablonovitch, “Inhibited spontaneous emission in solid-state physics and electronics,” Phys. Rev. Lett., vol. 58, pp. 2059–2062, 1987. [2] J. D. Joannopoulos, R. D. Meade, and J. N. Winn, Photonic Crystals. Princeton, New York: Princeton Univ. Press, 1995. [3] O. Painter, R. K. Lee, A. Scherer, A. Yariv, J. D. O’Brien, P. D. Dapkus, and I. Kim., “Two-dimensional photonic band-gap defect mode laser,” Science, vol. 284, pp. 1819–1821, 1999. [4] O. Painter, K. Srinivasan, J. D. O’Brien. A. Scherer and P. D. Dapkus, “Tailoring of the resonant mode properties of optical nanocavities in two-dimensional photohic crystal slab waveguides,” J. Opt. A: Pure Appl. Opt., vol. 3, pp. 161–170, 2001. [5] H. G. Park, J. K. Hwang, J. Huh, H. Y. Ryu, Y. H. Lee, and J. S. Kim, “Nondegenerate monopole-mode two-dimensional photonic band gap laser,” Appl. Phys. Lett. vol. 79, pp. 3032–3034, 2001. [6] H. G. Park, J. K. Hwang, J. Huh, H. Y. Ryu, S. H. Kim, J. S. Kim, and Y. H. Lee, “Characteristics of modified single-defect two-dimensional photonic crystal lasers,” IEEE J. Quantum Electron., vol. 38, pp. 1353–1365, 2002. [7] Vˇ ucković, Marko Lonćar, Hideo Mabuchi, and Axel Scherer, “Design of photonic crystal microcavities for cavity QED,” Phys. Rev. E, vol. 65, 016608, 2001. [8] Vˇ ucković, Marko Lonćar, Hideo Mabuchi, and Axel Scherer, “Optimization of the Q factor in photonic crystal microcavities,” IEEE J. Quantum Electron., vol. 38, pp. 850–856, 2002.. 41.

(51) [9] P. T. Lee, J. R. Cao, S. J. Choi, Z. J. Wei, J. D. O’Brien, and P. D. Dapkus, “Operation of photonic crystal membrane lasers above room temperature,” Appl. Phys. Lett., vol. 81, pp. 3311–3313, 2002. [10] P. T. Lee, J. R. Cao, S. J. Choi, Z. J. Wei, J. D. O’Brien, and P. D. Dapkus, “Roomtemperature operation of VCSEL-pumped photonic crystal lasers,” IEEE Photon. Technol. Lett., vol. 14, pp. 435–437, 2002. [11] E. M. Purcell, “Spontaneous emission probabilities at radio frequencies,” Phys. Rev., vol. 69, p. 681, 1946. [12] M. Boroditsky, R. Vrijen, T. F. Krauss, R. Coccioli, R. Bhat, E. Yablonovitch, “Spontaneous emission extraction and Purcell enhancement from Thin-Film 2-D Photonic Crystals,” J. Lightwave Technol., vol. 17, pp. 2096–2112, 1999. [13] K. S. Yee, “Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media,” IEEE Trans. Antennas Propag., vol. 14, pp. 302–307, 1966. [14] J. P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys, vol. 114, pp. 185–200, 1994. [15] Sullivan, “An unsplit step 3-D PML for use with the FDTD method,” IEEE Microwave Guided Wave Lett., vol. 7, pp. 184–186, 1997. [16] Z. S. Sacks, D. M. Kingsland, R. Lee, and J. F. Lee, “A perfectly matched anisotropic absorber for use as an absorbing boundary condition,” IEEE Trans. Antennas Propagat., vol. 43, pp. 1460–1463, 1995. [17] J. D. Moerloose and D. D. Zuttler, “Poynting’s theorem for the finite-difference-timedomain method,” Micro. Opt. Technol. Lett. Vol. 8, pp. 257–260, 1995. [18] A. Taflove and S. C. Hagness, Computational Electrodynamics : the finite-difference time-domain method, 2nd ed., Boston: Artech House, 2000. [19] M. Tinkham, Group Theory and Quantum Mechanics, New York: McGraw-Hill, 1964. [20] K. Sakoda, Optical Properties of Photonic Crystals, Berlin: Springer, 2001.. 42.

(52) [21] E. Yablonovitch, T. J. Gmitter, R. D. Meade, A. M. Rappe, K. D. Brommer, and J. D. Joannopoulos, “Donor and acceptor modes in photonic band structure,” Phys. Rev. Lett., vol. 67, pp. 3380–3383, 1991. [22] L. A. Coldren and S. W. Corzine, Diode Lasers and Photonic Integrated Circuits, New York: John Wiley & Sons, 1995. [23] J. D. Jackson, Classical Electrodynamics, 3rd ed., New York: Wiley, 1999. [24] W. H. Press, et al., Numerical Recipes in C: The Art of Scientific Computing, 2nd ed., New York: Cambridge, 1992.. 43.

(53)