Optics Communications 284 (2011) 226–230
Contents lists available at ScienceDirect
Optics Communications j o u r n a l h o m e p a g e : w w w. e l s ev i e r. c o m / l o c a t e / o p t c o m
The optimization of large gap–midgap ratio photonic crystal with improved Bisection-Particle Swarm Optimization Bin Jiang, Anjin Liu, Wei Chen, Mingxin Xing, Wenjun Zhou, Wanhua Zheng ⁎ Nano-optoeletronics, Institute of Semiconductors, Chinese Academic of Science, Beijing 100083, China
a r t i c l e
i n f o
Article history: Received 22 January 2010 Received in revised form 6 July 2010 Accepted 23 August 2010 Keywords: Particle Swarm Optimization Bisection-Particle Swarm Optimization Gap–midgap ratio
a b s t r a c t Photonic crystals possessing large gap–midgap ratios are designed using the improved Bisection-Particle Swarm Optimization (BPSO) algorithm and plane-wave expansion method. The optimal gap–midgap ratios of three kinds of photonic crystals, such as ring rods, square rods connected by veins, and rotating square rods, are calculated, and the corresponding geometry parameters are also represented. Compared with the traditional dimension-by-dimension scanning method, the improved BPSO algorithm greatly reduces the computation time, and its accuracy is controllable. © 2010 Elsevier B.V. All rights reserved.
1. Introduction Photonic crystal is a dielectric material with periodic or quasiperiodic dielectric distribution, and possesses some unique optical characteristics, such as negative refraction and band gap [1,2]. Especially, photonic crystal with complete (transverse electric and transverse magnetic, TE and TM) band gap can realize the modulation of photon, and has wild application perspectives in optical communication, such as photonic crystal filters, lasers, and waveguides [3–6]. In these applications, the size of the band gap plays an important role. The larger the band gap is, the wider applicable frequency it has. One of the main factors that will restrict the size of the band gap is the degeneracy at the high symmetry points in the Brillouin zone. The structure with low symmetry, for example, rectangular lattice [7], the combination of two circular scatterers [8], and no-circular scatterers, such as square and hexagonal scatterers [9], can destroy the degeneracy and enlarge the size of the band gap. From the viewpoint of the materials, the symmetry of the scatterer can be further reduced. For example, the anisotropic materials can increase the size of the band gap [10]. However, this method requires some special materials, thus its application is constrained. Moreover, the rotation of the scatterer can reduce the symmetry, and then increase the size of the band gap [11]. In order to design the photonic crystal structures with large gap– midgap ratios, we have to design the suitable scatterer structures, and optimize the geometry parameters. In the condition of the multiparameter optimization, the traditional dimension-by-dimension
⁎ Corresponding author. E-mail address:
[email protected] (W. Zheng). 0030-4018/$ – see front matter © 2010 Elsevier B.V. All rights reserved. doi:10.1016/j.optcom.2010.08.048
scanning method requires a huge amount of time. Therefore, a high efficient and accurate optimization method is needed. In recent years, in order to obtain the optimal structure parameters, a lot of studies on the computation method have been done and artificial intelligence has been introduced to the optimization of photonic crystal band gap. For example, Shen et al. [12] used the genetic algorithm to design the photonic crystal structure with a large gap–midgap ratio. But the structure obtained by genetic algorithm sometimes has micro structures, and it is difficult for fabrication. Furthermore, the gap– midgap ratio is quite sensitive to the lattice constant and air dutyratio. When the gap–midgap ratio is chosen as the fitness, the variance will be large, and many individuals with zero fitness will exist in the population, thus excellent genes will be aggregated only in some certain individuals. Although new excellent individuals can be generated through hybridization and mutation, it is still possible to converge to the local optimal solution prematurely. So some work will be needed to transform the fitness. In this paper, We propose the improved Bisection-Particle Swarm Optimization (BPSO) algorithm according to the characteristics of the photonic crystal band gap, and use it together with plane-wave expansion (PWE) method to design the photonic crystal structures with large gap–midgap ratios, such as, ring rods, square rods connected by veins, and rotating square rods. 2. Theory method and model 2.1. Bisection-Particle Swarm Optimization algorithm Particle Swarm Optimization (PSO) algorithm originated from the observation of the flock feeding behavior was first proposed by Dr. Eberhart and Dr. Kennedy in 1995 [13,14]. It belongs to a branch
B. Jiang et al. / Optics Communications 284 (2011) 226–230
of evolutionary algorithm. In the PSO algorithm, the particles can dynamically adjust their flight speed and direction according to their companions' flying experience and their own flying experience, then constantly move closer to the optimal particle, and eventually reach the optimal solution to the problem. In this paper, with the full study of the characteristics of the photonic crystal band gap, we propose the Bisection-Particle Swarm Optimization (BPSO) algorithm, which can dynamically adjust the speed of the particle swarm, and fully consider the impact between the particles. In addition, after each halving interval N, the algorithm halves the span of each search dimension with the global optimal particle as the center. In this way, the algorithm can achieve the higher accuracy and faster search speed. 2.2. BPSO algorithm proceeding
2.
3.
In Fig. 1, we provide the flow chart of BPSO algorithm, the detailed description of the algorithm is represented as follows. 1. Initialize the global parameters. Global parameters, such as particle swarm size, evolution generation, halving interval N, termination conditions, position and speed limits, initial position and velocity of 4. 5. 6.
7.
8.
227
each particle, optimization parameters (For the two dimensional (2D) triangular lattice of air hole rods, the only optimization parameter is the radius r of the air hole rod (0 b r b 0.5a, a is the lattice constant)) are initialized. Initialize the optimal value of each individual in the particle swarm and the optimal value of the particle swarm. The gap–midgap ratio of the each individual is calculated and then defined as the fitness of each individual (the higher the fitness is, the better the individual is). The fitness of each particle i is assigned to the optimal fitness Fif of each particle, and the initial location of each particle i is recorded as the optimal position Pif. Furthermore, the largest fitness of the particle swarm is assigned to the global optimal fitness Fgf, and the position of the corresponding particle is recorded as the global optimal position Pgf. Determine whether to halve the search space. If the evolution generation is the multiple of the halving interval N, the search space will be halved as shown in Fig. 2. We make the position of the optimal individual as the center, halve each dimension of the search space, and also halve the speed limit. After that, new swarm is randomly generated within the newly generated search space, and the fitness and the position of the global optimum are assigned to the first particle of the newly generated swarm. Update the velocity and position of each particle in accordance with the position and velocity renewal equations. Calculate the gap–midgap ratio of each particle. Compare the fitness of each particle i. If the current fitness of each particle is larger than its ever experienced optimal fitness Fif, the current position is recorded as Pif. Also, if the current optimal fitness of the particle swarm is larger than the global optimal fitness Fgf, the global optimal position Pgf is updated with the position Pif of the individual with the current global optimal fitness. Determine whether the termination conditions (For example, the expected fitness is achieved, or the predetermined maximal evolution generation is reached) are achieved. If not, return to Step 3. Otherwise go to step 8. Output the result.
In the PSO algorithm, the velocity of particle is essential to the convergence speed and the accuracy of the algorithm. In this paper, we adjust the renewal equations on the basis of the characteristics of photonic crystal band gap. Particle's velocity is composed of the following three components: 1) The inertia part of the velocity. It embodies the trust of the particle's current state, and the particle will move inertially according to its former speed. 2) The cognitive part of the velocity. It embodies the particle's own thinking. The particles will move referring to their prior moving experience, and enable particles to have a sufficiently strong global search ability to avoid the local minimum. 3) The social part of the velocity. It embodies the information sharing and mutual cooperation between various particles, that is, the movement of particles will be affected by the other particles in the particle swarm. Each particle can imitate the behavior of its better companions, and reach the optimal value faster.
Fig. 1. The flow chart of Bisection-Particle Swarm Optimization (BPSO) algorithm.
Since the band structure is sensitive to the change of the geometry parameters, the use of gap–midgap ratio as the fitness will make the variance of the population fitness large. As a result, the excellent information will be concentrated only in a small number of particles, which makes the optimal solution prematurely converge to the local optimal solution. Thus, in the early period, the cognitive part should be emphasized, which helps the particle swarm extend to a more broad scope, and find the global optimal solution. So the cognitive part is multiplied by a gradually decreasing factor C1. On the other hand, as the particle swarm evolves after some generations and gradually converges to the global optimal solution, in order to improve the accuracy of the
228
B. Jiang et al. / Optics Communications 284 (2011) 226–230
Fig. 2. The schematic diagram of search space halving in the 2D search using BPSO algorithm. (a) The optimal individual is selected at the end of each halving interval (N generations) and (b) the position of the optimal individual is set as the center, and each dimension of all search span is halved.
algorithm, social part shall be weighted more. Therefore, the social part is multiplied by a gradually increasing factor C2. To sum up, we obtain the improved speed and displacement renewal equations: Particle's velocity: vid: = w ⁎ vid + C1 ⁎ rand () ⁎ (Pid − xid) + C2 ⁎ rand () ⁎ (Pgd − xid). Particle displacement: xid: = xid + vid. Among them, acceleration factors (C1 is larger than 1, and is regressive; C2 is smaller than 1, and is incremental) are formulated as follows: C1 =1 + e− i/nmax C2 = 1−e− i/nmax (i is the evolution generation, and nmax is the predetermined maximal evolution generation). In the BPSO algorithm, after every halving interval N, the search space will be halved in each dimension. As shown in Fig. 3, if the size
of swarm remains the same in each generation, the search density will be increased exponentially. In the case of two-parameter optimization, when the evolution generation reaches 4 N, the search density is as large as 256 times of that of the initial generation. In the case of the three-parameter optimization, when the evolution generation equals to 4 N, the search density is as large as 4096 times of that of the initial generation. So the improved BPSO algorithm can achieve high search precision. 2.3. Plane-wave expansion method In this paper, plane-wave expansion method is used to calculate the photonic crystal gap–midgap ratio. Using Bloch theorem, we can get the following master equations from Maxwell's equations [15,16]: 1 ω2 ∇ × ½∇ × EðrÞ = 2 Eðr Þ εðr Þ c ∇×
1 ω2 ∇ × H ðr Þ = 2 H ðr Þ: εðr Þ c
Particularly, for 2D photonic crystal, the master equations for the E-polarization and H-polarization are as follows:
ω2 EðGÞ c2
ω2 HðGÞ c2
E polarization ðTM modeÞ: ∑ jk+ Gj jk +G′ jε−1 G−G′ E G′ = ′
G
H polarization ðTE modeÞ: ∑ðk+ GÞ⋅ k+ G′ ε−1 G−G′ H G′ = ′
G
−1
ε (G–G′) is the inverse of the Fourier transform matrix of the distribution function of the dielectric constant index. E(G) and H(G) are the Fourier transform matrixes of the distribution functions of the electric field and magnetic field, respectively. Based on the above equations, the dispersion relation of Epolarization and H-polarization can be calculated. Here, we take the band gap of triangular lattice of circular dielectric rods as an example to verify the correctness of our PWE calculation program. In the calculation, the number of plane-waves is 625, and the obtained result is compared with the result calculated by the RSoft software. The calculated results are shown in Table 1. The obtained results are consistently well with that calculated by the RSoft software, and the calculation error is less than 0.1%. Therefore, in
Table 1 Taking band gap of the triangular lattice of dielectric rods for example, the result obtained by our program is compared with the result calculated by the RSoft software. R/a = 0.17601 Lower band Upper band Gap size ε = 11.4 edge of the gap edge of the gap Fig. 3. As the number of halving increases, the search density increases. (a) Twoparameter and (b) three-parameter.
This paper RSoft
0.301793 0.301871
0.490477 0.490436
Midgap
Gap–midgap ratio
0.18868 0.39613 0.476311 0.18856 0.39615 0.475988
B. Jiang et al. / Optics Communications 284 (2011) 226–230 Table 2 The scatterer structure is a circular rod. The largest gap–midgap ratio of TE, TM, and TE and TM modes and corresponding parameters obtained by BPSO algorithm.
Square lattice of circular dielectric rods Square lattice of circular air hole rods Triangular lattice of circular dielectric rods Triangular lattice of circular air hole rod
R(/a) (TE)
Gap–midgap ratio
R(/a) (TM)
Gap–midgap ratio
0.19369
0.03457
0.19464
0.37979
0.44010
0.13865
0.49157
0.25445
0.17574
0.09453
0.17602
0.47631
0.42494
0.49419
0.49555
0.18960
order to save computing time, 625 plane waves are used while searching for the optimal parameters by BPSO algorithm, and later 1681 plane waves are used to obtain more precise results. 3. Analysis of simulation results In the calculation, different photonic crystal scatterer structures (triangular lattice and square lattice) are optimized using BPSO algorithm. The dielectric material GaAs (ε = 11.4) is used, and the lattice constant a equals to 1 μm. Three kinds of photonic crystals, such as ring rods, square rods connected by veins, and rotating square rods, are calculated. In the case of the circular rod, the only optimized parameter is the radius R. In the BPSO algorithm, R is provided as 0 b R b 0.5a, the size of swarm is set to 25, the evolution generation is set to 10, and the halving interval N is set to 7. According to the rule of thumb [17], TE band gap is easier formed in the structure with the continuous distribution of the low refractive index materials; TM band gap is easier generated in the structure with the continuous distribution of the high refractive index materials. From Table 2, we can see that the optimized dielectric rod has a larger TM band gap and a smaller TE band gap, while the optimized air hole
229
Table 3 The scatterer structure is a ring rod. The optimal gap–midgap ratio of the TE mode and the corresponding parameters obtained by BPSO algorithm. R1 is the radius of the outer circular dielectric rod and R2 is the radius of the inner circular air holes rod.
Triangular lattice (TE) Square lattice(TE)
R1(/a)
R2(/a)
Gap–midgap ratio
0.45311 0.49745
0.29375 0.40126
0.27764 0.35732
rod has a larger TE band gap and a smaller TM band gap. The TE gap– midgap ratios of square lattice of circular dielectric rods and triangular lattice of circular dielectric rods are quite small, and both of them are less than 0.1. Then it means that there leaves a large room for improvement. In Fig. 4(a) and (b), a circular air hole rod is introduced in the center of the circular dielectric rod to form a ring rod, which changes the dispersion relation of the original circular dielectric rod and can obtain a larger gap–midgap ratio. In the case of ring rod, the optimized parameters are the radius R1 of the outer circular dielectric rod, and the radius R2 of the inner circular air hole rod. In the BPSO algorithm, R1 and R2 are provided as 0 b R2 b R1 b 0.5a, the size of swarm is set to 25, the evolution generation is set to 20, and the halving interval N is set to 8. As can be seen from Table 3, the TE gap–midgap ratio of square lattice and triangular lattice of ring dielectric rods are significantly improved. The TE gap–midgap ratio of square lattice ring dielectric rod is increased by 8 times, from 3.45% to 27.76%, and the TE gap–midgap ratio of triangular lattice of ring dielectric rod is increased by 4 times, from 9.45% to 35.7%. The increase of the gap–midgap ratio is because the ring structure can be considered as the combination of the dielectric array and the air hole array. The introduction of air hole array makes contribution to the generation of the larger TE gap–midgap ratio [17]. For the square lattice of square dielectric rods shown in Fig. 4(c), due to the isolated distribution of the high refractive index material, it has only a small TE band gap. Therefore, it only has a very small complete (TE and TM) gap–midgap ratio. In Fig. 4(d), because of the continuous distribution of the high refractive index material, it is easy
Fig. 4. Different photonic crystal scatterer structures optimized by BPSO algorithm in this paper. (a) Square lattice of ring dielectric rods (b) triangular lattice of ring dielectric rods (c) square lattice of square dielectric rods (d) square lattice of cross dielectric rods (e) square lattice of square dielectric rods connected by veins and (f) square lattice of rotating square dielectric rods.
230
B. Jiang et al. / Optics Communications 284 (2011) 226–230
Table 4 The scatterer is a square rod. The optimal gap–midgap ratios of TE, TM, and TE and TM modes and the corresponding parameters obtained by BPSO algorithm.
Square dielectric rod Square air hole rod
Length l(/a)(TE)
Gap–midgap ratio
Length l(/a)(TM)
Gap–midgap ratio
Length l(/a)(TE and TM)
Gap–midgap ratio
0.636861 0.815171
0.095681 0.271726
0.342771 0.934202
0.372135 0.103128
0.640027 0.755437
0.051707 0.003029
to generate a large TE gap–midgap ratio. In Fig. 4(e), combining the square rods with cross rods can make up the defect of the low TE gap– midgap ratio of square dielectric rods, and then enable the new structure to possess a larger complete (TE and TM) gap–midgap ratio. In the BPSO algorithm, the optimized parameters are the side length l (0 b l b a) of the square rod, and the width y (0 b y b a) of the cross rod. The size of swarm is set to 25, the evolution generation is set to 20, and the halving interval is set to 8. The obtained optimized results show that when l equals to 0.54503a, y equals to 0.06610a, the complete (TE and TM) band gap size equals to 0.0567 (ωa/2πc), which is slightly better than the result (0.0548 (ωa/2πc)) in the literature [8] using two circular rods in the square lattice, and the gap–midgap ratio equals to 12.94%. In addition to the last two methods (adding a circular airhole rod to the circular dielectric rod or adding a cross rod to the square dielectric rod), the rotation of the scatterers with low symmetry can also increase the gap–midgap ratio. Table 4 shows the optimized TE, TM, and TE and TM gap–midgap ratios and the corresponding side lengths for square dielectric rod and square air hole rod. The square dielectric rod obtains a large TM gap–midgap ratio and a small TE gap–midgap ratio, while the square air hole rod obtains a large TE gap–midgap ratio and a small TM gap–midgap ratio. Both of these two types of scatterers have extremely low TE and TM gap–midgap ratios, and the gap–midgap ratio of square air hole rod is almost 0. In Fig. 4(f), as the scatterer is a square rod with four-fold symmetry, the rotation of the scatterer can change the distribution of the refractive index material, thus can p modify the dispersion ffiffiffi relation of the scatterer. The length l 0blb 2 = 2a of the square scatterer and the rotation angle θ (0 b θ b 2π) are optimized by BPSO algorithm. In the BPSO algorithm, the size of swarm is set to 25, the evolution generation is set to 20, and the halving interval N is set to 8. As shown in Table 5, after the rotation of the square scatterer, the complete (TE and TM) gap–midgap ratio of square lattice of square dielectric rods is increased from 5.17% to 8.89%, while the complete (TE and TM) gap–midgap ratio of square lattice of square air hole rods is increased from 0.303% to 6.79%. These results show that, the rotation of the scatterers with low symmetry can improve the gap– midgap ratios to some extent. For multi-parameter optimization problem, the improved BPSO algorithm obtains more obvious advantages. Taking the rectangular dielectric rod for example, there exist three optimized geometry parameters: the length of rectangular dielectric rod la (0 b la b a), the width of rectangular dielectric rod lb (0 b lb b a), and the rotation angle θ (0 b θ b π/2). If the traditional dimension-by-dimension scanning method is used and the scanning step lengths are set to 0.01a, 0.01a, and 0.01, respectively, 1.57 million times of plane-wave expansion calculation are needed, which will cost 3052 hours by using our computing system. However, using BPSO algorithm, when the size of swarm is set to 100, and the evolution generation is set to 100, only 10,000 times of the plane-wave expansion calculation are needed,
Table 5 The scatterer is a rotating square rod. The optimal gap–midgap ratios and the corresponding parameters obtained by BPSO algorithm.
Square dielectric rod Square air hole rod
TE and TM TE and TM
Length l(/a)
θ (rad)
Gap–midgap ratio
0.561179 0.707107
2.298560 5.495716
0.088866 0.067876
which will cost just 20 hours, and is only 1/157 of that used by the traditional dimension-by-dimension scanning method. Obviously, the improved BPSO algorithm saves a lot of time and computing resources. Therefore, using BPSO algorithm to design photonic crystal with a large gap–midgap ratio can have much significance. 4. Conclusions In this paper, the improved BPSO algorithm and plane-wave expansion method are used to optimize the geometry parameters and spatial location of the scatterer in the photonic crystal to obtain a large gap–midgap ratio. The improved BPSO algorithm can dynamically adjust the search speed according to the characteristics of the photonic band gap, for example, the sensitivity of the structure parameters. Also, the introduction of halving interval N can continuously condense the search space and can improve search accuracy. Compared with traditional dimension-by-dimension scanning method, the improved BPSO algorithm can save a lot of time and computing resources, and has a higher search accuracy that can be also dynamically adjusted. Then it offers useful guidance for the further optimization of the structure of photonic crystal with a large gap–midgap ratio. In addition, the improved BPSO algorithm can also be used to design small mode volume and high-Q photonic crystal cavities, high transfer rate photonic crystal waveguide devices, etc. So it is a powerful tool for designing photonic crystal devices. Acknowledgements This work is supported by the National Basic Research Program of China (Grant No. 2006CB921705), the National Natural Science Foundation of China (Grant Nos. 10634080, 60677046 and 60838003), and the National High Technology Research and Development Program of China (Grant Nos. 2007AA03Z410 and 2007AA03Z408). References [1] S. John, Phys. Rev. Lett. 58 (1987) 2486. [2] E. Yablonvitch, Phys. Rev. Lett. 58 (1987) 2059. [3] K. Sakoda, Optical Properties of Photonic Crystals, Spriger-Verlag, Berlin Heidelberg, 2005, p. 187. [4] S. Noda, A. Chutinan, M. Imada, Nature 407 (2000) 608. [5] O. Painter, R.K. Lee, A. Scherer, A. Yariv, J.D. O' Brien, P.D. Dapkus, I. Kim, Science 284 (1999) 1819. [6] A. Chutinan, S. Noda, Phys. Rev. B 62 (2000) 4488. [7] M. Qiu, S.L. He, Phys. Rev. B 60 (1999) 10610. [8] C.M. Anderson, K.P. Giapis, Phys. Rev. Lett. 77 (1996) 2949. [9] Pierre R. Villeneuve, Michel Piche, Phys. Rev. B 46 (1992) 4969. [10] Z.Y. Li, B.Y. Gu, G.Z. Yang, Phys. Rev. Lett. 81 (1998) 2574. [11] X.H. Wang, B.Y. Gu, Z.Y. Li, G.Z. Yang, Phys. Rev. B 60 (1999) 11417. [12] L.F. Shen, Y. Zhuo, S.L. He, Phys. Rev. B 68 (2003) 035109. [13] J. Kennedy, R. Eberhart, Particle swarm optimization [C], Proc. of the IEEE Int. Conf. on Neural Networks, Piscataway, NJ, 1995, p. 1942. [14] R. Eberhart, J. Kennedy, A new optimizer using particle swarm theory [C], Proc of the Sixth International Symposium on Micro Machine and Human Science, Nagoya,Japan, 1995, p. 39. [15] K.M. Ho, C.T. Chan, C.M. Soukoulis, Phys. Rev. Lett. 65 (1990) 3152. [16] S.G. Johnson, J.D. Joannopoulos, Opt. Express 8 (2001) 173. [17] J.D. Joannopoulos, S.G. Johnson, J.N. Winn, R.D. Meade, Photonic Crystals: Molding the Flow of Light, Princeton University Press, 2008, p. 66.