Bi-directional evolutionary optimization for photonic band gap structures

Journal of Computational Physics 302 (2015) 393–404

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Bi-directional evolutionary optimization for photonic band gap structures Fei Meng a,d , Xiaodong Huang a,b,∗ , Baohua Jia c a

Centre for Innovative Structures and Materials, School of Civil, Environmental and Chemical Engineering, RMIT University, GPO Box 2476, Melbourne, VIC 3001, Australia b Key Laboratory of Advanced Technology for Vehicle Body Design & Manufacture, Hunan University, Changsha, 410082, China c Centre for Micro-Photonics, Faculty of Engineering & Industrial Science, Swinburne University of Technology, PO Box 218, Hawthorn, VIC 3122, Australia d School of Civil Engineering, Central South University, Changsha 410075, China

a r t i c l e

i n f o

Article history: Received 19 January 2015 Received in revised form 6 September 2015 Accepted 7 September 2015 Available online 14 September 2015 Keywords: Topology optimization Photonic band gap Bi-directional evolutionary structural optimization (BESO) Periodic unit cell

a b s t r a c t Toward an efficient and easy-implement optimization for photonic band gap structures, this paper extends the bi-directional evolutionary structural optimization (BESO) method for maximizing photonic band gaps. Photonic crystals are assumed to be periodically composed of two dielectric materials with the different permittivity. Based on the finite element analysis and sensitivity analysis, BESO starts from a simple initial design without any band gap and gradually re-distributes dielectric materials within the unit cell so that the resulting photonic crystal possesses a maximum band gap between two specified adjacent bands. Numerical examples demonstrated the proposed optimization algorithm can successfully obtain the band gaps from the first to the tenth band for both transverse magnetic and electric polarizations. Some optimized photonic crystals exhibit novel patterns markedly different from traditional designs of photonic crystals. © 2015 Elsevier Inc. All rights reserved.

1. Introduction Photonic crystals are optical structures consist of dielectric materials with different refractive indexes. They have lattice constants around the wavelength of light and periodicity in one, two or three dimensions. Photonic crystals are also called photonic band gap (PBG) structures because of their ability of prohibiting the propagation of electromagnetic waves within certain frequency ranges [1,2]. In practice, a broader band gap means broader available bandwidth of signals and applications, so it is of great significance to design photonic crystals with large band gaps. The optical properties of photonic crystals depend not only on the properties of dielectric materials but also their spatial distributions [3]. Therefore, for given materials, the design of photonic crystals becomes a typical topology optimization question: how to periodically distribute the materials to maximize the specific band gap. Due to the polarization of electromagnetic waves, transverse magnetic polarization (TM mode) and transverse electric polarization (TE mode) can be considered separately. The traditional design approach of photonic band gap (PBG) structures is a trial-and-error process based on physical intuitions and parametric studies. This process would be inefficient and time-consuming [4–8], and the resulting design may also be away from the optimum. The systematic way to find the optimal design of photonic crystals is formulating the

*

Corresponding author at: Centre for Innovative Structures and Materials, School of Civil, Environmental and Chemical Engineering, RMIT University, GPO Box 2476, Melbourne, VIC 3001, Australia. Tel.: +61 3 99253320; fax: +61 3 96390138. E-mail address: [email protected] (X. Huang). http://dx.doi.org/10.1016/j.jcp.2015.09.010 0021-9991/© 2015 Elsevier Inc. All rights reserved.

394

F. Meng et al. / Journal of Computational Physics 302 (2015) 393–404

problem with appropriate objective functions and/or constraints, and then solving it by topology optimization methods. The topology optimization methods were originally used in elastic field to seek the best materials layout within a given design space so that the resulting structures had the maximum stiffness or minimum compliance. So far, the commonly used topology optimization methods include Solid Isotropic Material with Penalization (SIMP) method [9–12], level set method [13–17], Evolutionary Structural Optimization (ESO) [18,19] and its later version bidirectional ESO (BESO) [20–22]. In order to obtain photonic or phononic crystals with large band gaps, many researchers have investigated into the design of PBG structures using various topology optimization methods, e.g. genetic algorithms [23–28], level set method [29–31], SIMP [32–35]. These methods have proven to be useful and some interesting designs have been obtained. However, because of the complicate nature of the solution space, a systemic topology optimization for PBG structures with any band gap is still a challenge. Sigmund and Jensen [34,35] proposed two-stage optimization where the first stage using coarse grids generated five best topologically different candidates for each band. Starting from the candidate topologies with the opening of the band gap, the SIMP method was then employed at the second stage to find the optimal PBG structures. Men et al. [36,37] systematically reported the design of PBG structures using semidefinite programming and subspace method. The optimization starts from a series of randomly generated initial designs, the photonic band gap could be obtained but the occurrence of the band gap could not be guaranteed. While these methods are attractive and have successfully obtained the optimized solutions of photonic band gap crystals, the computational cost for the optimization based on finite element analysis (FEA) is still expensive because hundreds, even thousands of iterations are needed. More importantly, the previous research also indicated that the obtained solutions depend highly on the used optimization parameters and algorithm [34–37]. This is because of the fact that a number of different topologies could possess the same band gap property and there is no unique solution for the design of PBG structures. Therefore, it is important to attempt new and different optimization algorithms, such as BESO, in order to find a much wider range of possible solutions to the PBG structures. Toward an efficient but easy-implement optimization of PBG structures, a new approach based on the BESO method [20–22] is proposed in this paper. Based on FEA, BESO gradually removes elements from or adds elements to the design domain so that the resulting topology evolves towards an optimum. It has been demonstrated that the current BESO method is capable of generating reliable and practical topologies for optimization problems with various constraints such as natural frequency [38], compliant mechanisms [39] and nonlinear structures [40,41]. In recent years, the BESO algorithm has also been extended to the design of materials with mechanical or electromagnetic properties [42,43]. In this paper, we start from finite element discretization of photonic crystals for the propagation of electromagnetic waves. Based on FEA and the sensitivity analysis for the single or multiple eigenfrequencies, a BESO method is established by using discrete design variables. Starting from a simple initial design without any band gap, BESO evolves the topology of the unit cell to the optimized structures with the desired band gap. Finally, numerical examples are given to demonstrate the effectiveness and efficiency of the proposed optimization algorithm for the design of photonic crystals. 2. Finite element discretization Photonic crystals are usually composed of two or more homogeneous materials distributed periodically as shown in Fig. 1(a). The propagation of electromagnetic waves in photonic crystals is generally governed by Maxwell’s equations [1]. However, there are two possible polarizations of the magnetic and electric fields for 2D cases, namely TM (transverse magnetic) and TE (transverse electric) modes. In the TM mode, the magnetic field is confined to the plane of wave propagation and the electric field E = (0, 0, E (k, r)) is perpendicular to this plane, where the vector r = (x, y ) denotes the coordinates in the plane and k is the wave vector. In contrast, the electric field of the TE mode is confined to the plane of wave propagation and the magnetic field H = (0, 0, H (k, r)) is perpendicular to this plane. Since there are no point sources or sinks of electric displacement and magnetic fields in photonic crystal, the time-harmonic Maxwell equations can be reduced to two decoupled equations as

 2   ω −∇ · ∇ E (k, r) = ε (r) E (k, r) for TM mode  −∇ ·

c

1

ε(r)

  2 ω H (k, r) for TE mode ∇ H (k, r) = c

(1a) (1b)

where c is the speed of light, ω is the angular frequency of the electromagnetic wave and ε (r) is the dielectric function. Due to the periodicity of the crystal, the dielectric function satisfies ε (r) = ε (r + R) and R is the lattice translation vector. According to the Bloch–Floquet theory [44], the magnetic and electric fields can be represented as the product of a periodic function and an exponential factor as

H (k, r) = H (r) exp(ik · r) E (k, r) = E (r) exp(ik · r)

for TM mode for TE mode

(2a) (2b)

Thus, Maxwell equations can be further converted to the eigenvalue problems within a unit cell.

 2   ω −(∇ + ik) · (∇ + ik) E (r) = ε E (r) for TM mode c

(3a)

F. Meng et al. / Journal of Computational Physics 302 (2015) 393–404

395

Fig. 1. (a) Photonic crystals with 3 × 3 unit cells. (b) Irreducible Brillouin zone ( -X-M- ). (c) Band diagram.

 −(∇ + ik) ·

1

ε

  2 ω (∇ + ik) H (r) = H (r) for TE mode c

(3b)

Due to the periodicity of the unit cell, only the wave vectors on the boundary of irreducible Brillouin zone are considered for the calculation of the band diagram. Fig. 1(b) shows the irreducible Brillouin zone and its boundary ( -X-M- ) for the photonic crystal with a square lattice vector R = (a, a). However, it is impossible to sweep all points along the boundary of Brillouin zone [44]. Numerical examples show that dividing each side into 10 equal segments with 11 points is enough for the accurate band diagram. The wave vector ki = (kx , k y ) (i = 1, · · · , 31) starts from the point (0, 0), to X(π /a, 0), then M(π /a, π /a) and finally comes back to the point  . For the given wave vector ki = (kx , k y ), Eqs. (3a) and (3b) become the typical eigenvalue problems and can be easily solved through finite element analysis (FEA). To conduct FEA, the unit cell is discretized into finite elements. The components E (r) and H (r) are expressed as [45] n 

E (r) =

N i (r) E ei

and

H (r) =

i =1

e

n  e

N i (r) H ie

(4)

i =1

where e denotes element and n is the total number of nodes for each element. N = { N 1 N 2 N 3 N 4 } is the shape function for a 4-node square element. Following the standard formulation of finite element method, Eqs. (3a) and (3b) can be written in matrix form as





e

E=

K

e

 e

1

εe

 2  

ω c

 e

K

H=



e

εe M E for TM mode

e

 2 

ω c

(5a)

 e

M

H

for TE mode

(5b)

e

where H and E are the eigenvectors for the magnetic and electric fields, respectively. Ke can be expressed by

Ke = K1 + K2 + K3 + K4 where,

 

K1 = A

 ∂ NT ∂ N ∂ NT ∂ N + dA ∂x ∂x ∂y ∂y

 

K2 = ikx A

 ∂N ∂ NT dA N − NT ∂x ∂x

 

K3 = ik y A

K4 =



k2x

 ∂N ∂ NT dA N − NT ∂y ∂y   2 T

+ ky

N Nd A

A

(6)

396

F. Meng et al. / Journal of Computational Physics 302 (2015) 393–404

where A denotes the total area of an element. Me is expressed by



Me =

NT Nd A

(7)

A

By sweeping the wave vector along the boundary of Brillouin zone, the band diagram such as Fig. 1(c) can be obtained. 3. Optimization formulation 3.1. Objective function Our aim for designing a photonic crystal is to maximize the band gap at the desired frequency band. Due to the lack of fundamental length scale in Maxwell’s equation, the band gap–midgap ratio, which is independent of the lattice constant of the photonic crystal, is more meaningful than the absolute value of the band gap. Therefore, the objective function for designing photonic crystals can be changed to maximize the band gap–midgap ratio between two adjacent bands (referred to as band m and band m + 1) as

max f (xe ) = 2

min ωm+1 (k) − max ωm (k)

(8)

min ωm+1 (k) + max ωm (k)

where xe is the design variable which will be explained later. It should be noted that the frequency is the function of the wave vector. The above optimization problem is hardly solved directly because the locations of max ωm and min ωm+1 are changing with the change of the topology during the optimization process. To overcome this difficulty, it is assumed that any point from the (m + 1)th to the highest band whose frequency is less than c 2 = 1.1 × min ωm+1 (k) as line shown in Fig. 1(c) may become the location of min ωm+1 in the next iteration. Similarly, any point from the mth to the lowest band whose frequency larger than c 1 = 0.9 × max ωm (k) may become the location of max ωm in the next iteration. Thus, a new objective function is established as

31 m 31 ∞ ˆ i j ω j (ki ) j =m+1 w i j ω j (ki ) − i =1 i =1 j =1 w ¯ max f (xe ) = 2 31 ∞ 31 m ˆ i j ω j (ki ) j =m+1 w i j ω j (ki ) + i =1 i =1 j =1 w

(9)

ˆ i j are weight factors. Firstly, we define the following two parameters where w i j and w



c 2 − ω j (ki ) 0 ˆ i j = ω j (ki ) − c 1 A 0 Ai j =

when ω j (ki ) < c 2 otherwise when ω j (ki ) > c 1 otherwise

(i = 1, · · · , 31; j = m + 1, · · · , ∞)

(10a)

(i = 1, · · · , 31; j = 1, · · · , m)

(10b)

Then, the weight factors are given as

Ai j w i j = 31 ∞ i =1

j =m+1

ˆ ij A ˆ i j = 31 m w i =1

As a result,

j =1

31 ∞ i =1

(11a)

Ai j

(11b)

ˆ ij A

j =m+1

w i j = 1 and

31 m i =1

j =1

ˆ i j = 1. w

3.2. Design variable The design variable xe is the artificial variable of an element which will be used to represent the topology of a unit cell. It is assumed that a photonic crystal is composed of two materials with permittivity ε1 and ε2 (where ε1 < ε2 ). The elemental design variable xe can be constructed by following that the element is composed of material 1 with permittivity ε1 when xe = 0 and material 2 with permittivity ε2 when xe = 1. The traditional BESO method [20–22] used the discrete design variable xe = 0 or 1 only. However, numerical experience shows that the design of photonic crystals is very sensitive to the change of the design variable so that a very fine mesh should be used. Instead of using an extreme fine mesh, the design variable of the BESO method in this paper can be assigned with discrete intermediate design values. In other words, the variation of the design variable in each iteration is limited to be xe = 0.1. For the compliance minimization of structures, the well-known SIMP model [12] makes elements with intermediate design variables density uneconomical in the optimization process and thus the solution naturally tends to be 0/1. In the design of photonic crystals, a 0/1 solution results in a larger band gap since intermediate design variables reduce the contrast between the materials. Therefore, the linear material interpolation law has been applied successfully to the design of photonic crystals for TM modes [23,34] as

F. Meng et al. / Journal of Computational Physics 302 (2015) 393–404

397

Fig. 2. Material interpolation schemes for TM and TE modes.

ε(xe ) = ε1 (1 − xe ) + ε2 xe for TM mode

(12)

In the case of TE modes, the above linear material interpolation has also been applied in the literature [34]. However, our numerical experience indicated that this material interpolation scheme worked only for maximizing some TE band gaps. Here, we adopt the following inverse linear material interpolation for TE modes

1

ε(xe )

=

1

ε1

(1 − xe ) +

1

ε2

xe

for TE mode

(13)

The above material interpolation scheme has also been used in Refs. [24,35]. It is noted that this inverse linear material interpolation implies the penalization of permittivity for intermediate elements as shown in Fig. 2. 3.3. Sensitivity analysis For the given wave vector, ki , the corresponding frequencies, ω j (ki ) and eigenvectors, u j (ki ) ( j = 1, 2, · · ·) can be extracted by solving Eq. (1a) or Eq. (1b). The derivation of ω j (ki ) with regard to xe can be easily formulated as

  ∂ ω j (ki ) ∂M 1 ∂K = uTj c 2 − ω2j (ki ) uj ∂ xe 2ω j (ki ) ∂ xe ∂ xe

(14)

Using the material interpolation scheme, Eq. (12), for TM modes, the derivations of matrix K and M are

∂K =0 ∂ xe ∂M = (ε2 − ε1 )Me ∂ xe

(15a) (15b)

where the dimensions of Me are the same as those of M but all components in Me which are not related to element e are zero. Similarly, using the material interpolation scheme, Eq. (13), for TE modes, the derivations of matrix K and M are

  ∂K 1 1 = − Ke ∂ xe ε2 ε1 ∂M =0 ∂ xe

(16a) (16b)

Substituting the above equations into Eq. (14), the sensitivity of all eigenfrequencies can be easily obtained since their corresponding eigenvectors can be outputted from FEA. Therefore, no additional computational cost is needed for sensitivity analysis, excluding storage consideration. It should be noted that the above sensitivity analysis is only suitable for the case of single eigenvalue. Multiple eigenvalues are widely existed in the design of photonic crystals. Investigation of sensitivity analysis of multiple eigenvalues is available in many papers [46–48]. Take the case of a double eigenvalue with two corresponding eigenvectors (ω, u1 , u2 ) as an example, but the extension of a higher number of multiplicities is straightforward. The multiplicity of the eigenvalue implies that any linear combination of the eigenvectors u1 and u2 (e.g. c 1 u1 + c 2 u2 ) corresponding to ω , will also satisfy the original eigenvalue problem. The sensitivities of double eigenvalues are eigenvalues of a two-dimensional algebraic subeigenvalue problem [46] as

fsk c = 0

(17)

where c = {c 1 , c 2 } and fsk (s, k = 1, 2) denote generalized gradient vectors as T

398

F. Meng et al. / Journal of Computational Physics 302 (2015) 393–404

fsk =



1 2ω

uTs c 2



∂K ∂M − ω2 uk ∂ xe ∂ xe

(18)

In the practical implementation in each iteration step of the optimization we need to check if there are multiple eigenvalues for any given wave vectors and in this case calculate the sensitivities according to Eq. (17). With the sensitivities for multiple eigenvalues, it is natural to assign the lowest sensitivity to the lowest order eigenvalue and the highest sensitivity to the highest order eigenvalue. Thus, the sensitivity number defined with the derivation of the modified objective function can be calculated accordingly

αe =

∂ ¯f (xe ) ∂ xe

(19)

4. Numerical implementation for BESO In the topology optimization of mechanical structures, the filtering scheme can effectively alleviate the numerical instabilities of the checkerboard pattern and mesh-dependency [20–22]. For the design of photonic crystals, the usage of the filter scheme is arguable [35]. Our numerical experience indicated that the checkerboard pattern does not exist in the design photonic crystals, but a different mesh may results in a different optimized design. This mesh-dependency problem is different from that of mechanical structures, where a fine mesh results in a more detail design of the structures. To avoid this mesh-dependency problem, the filter scheme [20–22] is still used for the design of photonic crystals in this paper. The filter scheme modifies the sensitivity number as

n

αˆ e = i=n1

w i αe

i =1

(20)

wi

where the weight factor w i is defined by



wi =

rmin − r ie , 0,

if rmin > r ie if rmin ≤ r ie

(21)

where r ie denotes the distance between the center of element e and i. rmin is the radius of the filter, defined to identify the

√

neighboring elements that affect the sensitivity number of element i. rmin = 2a/50 is used for all numerical examples in this paper. In order to improve the stability and convergence of solution, elemental sensitivity numbers can be further averaged with their corresponding values in the previous iteration as [20–22]

α˜ ek =

1 2

α˜ ek−1 + αˆ ek



(22)

where k is the current iteration number. Different from most of topology optimization problems, there is no volume constraint for the optimal design of photonic crystals. It is understandable as there is no band gap when the unit cell is full of either material 1 or material 2. Therefore, the optimal volume fraction of material 2 must be between 0% and 100%. For simplicity, the volume fraction, V f , is for material 2 in the thereafter sections. BESO starts from an initial design which is nearly full of material 2 and there is no band gap for the initial design. The optimization for the variation of the volume fraction has three stages. At the first stage, the volume fraction gradually decreases when it is larger than a prescribed volume, V ∗f , as

V kf +1 = V kf (1 − ER) when V kf +1 > V ∗f

(23)

where ER is the evolution rate. ER = 2% and V ∗f = 40% are used throughout this paper. Due to the complexity of the optimization of photonic band gap crystals, many local optimal solutions may exist. The setting of V ∗f means we want to find the optimal solution with the volume fraction around V ∗f . The volume fraction at the second stage is determined according to the variation history of the band gap as



V kf +1 = V kf 1 +

( f (xke ) − f (xke−1 )) · ( V kf − V kf −1 ) | f (xke ) − f (xke−1 )| · | V kf − V kf −1 |

 ER

(24)

The second stage is ended until the volume fraction is vibrantly convergent to a constant value, V ∗∗ . Then the volume f fraction keeps a constant at the third stage until the topology of the unit cell and objective function are stably convergent. Similar to the conventional BESO method [20–22], the design variables are updated according to the relative values of ˜ ek , a threshold sensitivity numbers and volume fraction. Based on the relative ranking of the calculated sensitivity numbers α k of the sensitivity number, αth , is determined by using the bi-section method so that the volume fraction in the next iteration is equal to V kf +1 as described in the above section. The design variable of each element is modified by comparing its sensitivity number with the threshold as

F. Meng et al. / Journal of Computational Physics 302 (2015) 393–404

399

Fig. 3. Flow chart of the BESO procedure.

xke+1 =



min(xke + x, 1), max(xke

k ˜ ek > αth if α

k − x, 0), if α˜ ek < αth

(25)

where x = 0.1 throughout the paper, which means BESO uses discrete design variables. Although BESO may take some discrete intermediate design variables during the optimization process, the final design is naturally convergent to an almost 0/1 design due to the adoption of the material interpolation schemes. Fig. 3 depicts the proposed BESO procedure for the design of photonic crystals with a desired band gap. Firstly, the FE analysis for the unit cell of a photonic crystal is conducted by COMSOL Multiphysics, to obtain the eigenfrequencies, corresponding eigenvectors, and the band diagram. Then, the sensitivity numbers are calculated and the topology of the unit cell is updated accordingly. The thereafter iterative process evolves the topology of the unit cell towards its optimum until both the topology and objective function are convergent. 5. Numerical results and discussion In this section, numerical examples will be presented to illustrate the effectiveness of the proposed optimization method. Here, we consider the design of photonic crystals with a square lattice a = 1, however, the proposed method can be equally applied for photonic crystals with other lattices. It is assumed that the photonic crystals are composed of two materials: air (ε1 = 1) and gallium arsenide (GaAs) (ε2 = 11.4). The unit cell represented the structure of photonic crystals is meshed √ with 64 × 64 bilinear square elements. The BESO parameters are ER = 2% and rmin = 2a/50. BESO starts from the initial design which is full of material 2 (GaAs) except for four elements in the center of the unit cell with material 1 (air). To illustrate microstructures in figures, white and black colors denote air and GaAs, respectively. Fig. 4 gives the optimized 3 × 3 unit cells of photonic crystals for TM modes with maximum band gaps from the first to tenth band (where the dashed box represents the unit cell). Their corresponding band diagrams are also given in

400

F. Meng et al. / Journal of Computational Physics 302 (2015) 393–404

Fig. 4. Optimized 3 × 3 unit cells and their band diagrams for TM band gaps. (a) The first band. (b) The second band. (c) The third band. (d) The fourth band. (e) The fifth band. (f) The sixth band. (g) The seventh band. (h) The eighth band. (i) The ninth band. (j) The tenth band.

F. Meng et al. / Journal of Computational Physics 302 (2015) 393–404

401

Fig. 5. Optimized 3 × 3 unit cells and their band diagrams for TE band gaps. (a) The first band. (b) The second band. (c) The third band. (d) The fourth band. (e) The fifth band. (f) The sixth band. (g) The seventh band. (h) The eighth band. (i) The ninth band. (j) The tenth band.

402

F. Meng et al. / Journal of Computational Physics 302 (2015) 393–404

Fig. 6. Evolution histories of objective function and volume fraction for maximizing the seventh band gap of TE mode.

Fig. 4. Generally, GaAs with the high permittivity is isolated by air with low permittivity which is consistent with the previous reports [34,36]. The distribution shapes of GaAs can be solid circles or circle rings. The obtained maximum band gap–midgap ratio for TM mode is 42.89% which is the gap between the sixth and seventh bands. Although the lattice constant a = 1 is used for all examples, it is interesting to find that some optimized topologies are similar but alike with a half of the lattice constant, e.g. the optimized topologies between the first and fourth bands or between the second and ninth bands. It should be noted that the resulting topology is not a unique solution for the design of photonic crystals. For instance, the equivalent solution can also be found by shifting the optimized unit cell with half the lattice along x and y directions simultaneously. It indicates well that the optimized design of photonic crystals highly depends on the initial design. Fig. 5 shows the optimized 3 × 3 unit cells of photonic crystals for TE modes with maximum band gaps from the first to tenth band and their corresponding band diagrams. It should be pointed out that we fail to obtain the third band gap using the given initial design. Instead, a new initial design which is full of material 2 (GaAs) except for two elements at the middle of four edges with material 1 (air) is used. Opposite to the optimized designs for TM modes, GaAs tends to be connected together for the optimized designs for TE modes except for those of the second, fourth and eighth band gaps. The obtained maximum band gap–midgap ratio for TE mode is 44.27% which is the gap between the seventh and eighth bands. It seems that maximizing band gaps at high frequency level usually generate more complex topologies which are hardly figured out by scientists. In order to illustrate the whole optimization process, Fig. 6 shows the evolution histories of objective function and volume fraction for maximizing the seventh band gap of TE mode. It can be seen that the band gap–midgap ratio for the initial design is −10.47%, which means no band gap. When the volume fraction of GaAs gradually decreases, the band gap–midgap ratio gradually increases from a negative value to a positive one which means the band gap turns up. This band gap continues to increase and then converges to the maximum value, 44.27%, meanwhile the volume fraction is also convergent to a stable value, 39.82%. The whole optimization process only needs 60 iterations. Fig. 7 gives the evolution history of the unit cell. It can be seen that the optimized topology is very close to a 0/1 design although there are a lot of gray areas during the optimization process. It can be seen that all optimized solutions in Figs. 4 and 5 have clear topologies without any blur areas although some discrete intermediate design variables are used. The obtained band gaps demonstrate the effectiveness of the proposed optimization method. Some optimized topologies are similar to those in the literature [23,29,34,36], but some topologies such as Figs. 4(f) and (h) for TM mode and Figs. 5(d), (e), (h), (i), (j) for TE modes are firstly reported. It should be noted that the computational burden for the optimization of photonic crystals mainly comes from a large number of FE analysis. Compared with other optimization methods with hundreds, even thousands of iterations, all of the present optimized designs are obtained with around 100 iterations and the computational efficiency of the proposed optimization method is therefore obvious. It should be noted that the optimization problem for PBG structures may have many local optima and the optimized solution may highly depend on the initial guess as mentioned before. To successfully obtain the desirable band gap, the optimization algorithms developed by Sigmund and Jensen [34,35] were only applied to the initial design with band gaps. Men et al. [36,37] examined two different types of initial configurations: photonic crystals exhibiting band gaps and random distribution. Starting from an initial random design without any band gap, the optimization algorithm often failed to find a proper design with band gaps, for example, the successful rates for the ninth TE and TM band gaps were only 13.3% and 16.7% respectively. In this paper, a simple initial design without any band gap is used and BESO can successfully find all band gaps except for the third TE band gap. By slightly changing the initial design, the TE third band gap can still be successfully obtained. Certainly, other solutions could also be found if BESO starts from another initial design. However, further investigation into seeking the global optima of photonic band gap crystals is still needed.

F. Meng et al. / Journal of Computational Physics 302 (2015) 393–404

403

Fig. 7. Evolution history of the unit cell for maximizing the seventh band gap of TE mode. (a) Initial design. (b) Iteration 10. (c) Iteration 20. (d) Iteration 30. (e) Iteration 40. (f) Final optimized design.

6. Conclusions This paper systematically investigates the topology optimization of 2D photonic crystals for both TM and TE modes. A new optimization algorithm based on BESO is proposed for the design of PBG structures. Through iterative finite element analysis, the optimization algorithm evolves the initial structure without any band gap to the optimized PBG structures with desired band gap. Numerical results indicate that the algorithm proposed in this paper is effective and applicable to any frequency level. Moreover, the proposed algorithm is computationally efficient as the solution usually converges with about 100 iterations. It can be easily implemented as a post-processor to standard commercial finite element analysis software packages, e.g. COMSOL Multiphysics, which has been used in this study. Further work is recommended for the manufacture and testing of the optimized PBG structures and topology optimization of 3D PBG structures. Acknowledgement The authors wish to acknowledge the financial support from the Australian Research Council (FT130101094 and DE120100291) and the China Scholarship Council. References [1] J.D. Joannopoulos, S.G. Johnson, J.N. Winn, R.D. Meade, Photonic Crystals: Molding the Flow of Light, Princeton University Press, 2008. [2] J.G. Fleming, S.-Y. Lin, I. El-Kady, R. Biswas, K.M. Ho, All-metallic three-dimensional photonic crystals with a large infrared bandgap, Nature 417 (2002) 52–55. [3] E.E. Hart, A. Sobester, K. Djidjeli, M. Molinari, K.S. Thomas, S.J. Cox, A geometry optimization framework for photonic crystal design, Photonics Nanostruct. Fundam. Appl. 10 (2012) 25–35. [4] M. Doosje, B.J. Hoenders, J. Knoester, Photonic bandgap optimization in inverted fcc photonic crystals, J. Opt. Soc. Am. B 17 (2000) 600–606. [5] J.H. Wu, A.Q. Liu, L.K. Ang, Band gap optimization of finite photonic structures using apodization method, J. Appl. Phys. 100 (2006) 084309. [6] M.M. Hossain, G. Chen, B. Jia, X.H. Wang, M. Gu, Optimization of enhanced absorption in 3D-woodpile metallic photonic crystals, Opt. Express 18 (2010) 9048–9054. [7] D. Wang, Z. Yu, Y. Liu, P. Lu, L. Han, H. Feng, X. Guo, H. Ye, The optimal structure of two dimensional photonic crystals with the large absolute band gap, Opt. Express 19 (2011) 19346–19353. [8] P. Shi, K. Huang, Y. Li, Photonic crystal with complex unit cell for large complete band gap, Opt. Commun. 285 (2012) 3128–3132. [9] M.P. Bendsøe, Optimal shape design as a material distribution problem, Struct. Optim. 1 (1989) 193–202. [10] G.I.N. Rozvany, M. Zhou, T. Birker, Generalized shape optimization without homogenization, Struct. Optim. 4 (1992) 250–254. [11] M. Zhou, G.I.N. Rozvany, The COC algorithm, part II: topological, geometry and generalized shape optimization, Comput. Methods Appl. Mech. Eng. 89 (1991) 197–224. [12] M.P. Bendsøe, O. Sigmund, Topology Optimization: Theory, Methods and Applications, Springer-Verlag, Berlin, 2003. [13] M.Y. Wang, X. Wang, D. Guo, A level set method for structural topology optimization, Comput. Methods Appl. Mech. Eng. 192 (2003) 227–246.

404

F. Meng et al. / Journal of Computational Physics 302 (2015) 393–404

[14] X. Wang, M.Y. Wang, D. Guo, Structural shape and topology optimization in a level-set-based framework of region representation, Struct. Multidiscip. Optim. 27 (2004) 1–19. [15] J.A. Sethian, A. Wiegmann, Structural boundary design via level set and immersed interface methods, J. Comput. Phys. 163 (2) (2000) 489–528. [16] S. Zhou, W. Li, G. Sun, Q. Li, A level-set procedure for the design of electromagnetic metamaterials, Opt. Express 18 (2010) 6693–6702. [17] S. Zhou, W. Li, Y. Chen, G. Sun, Q. Li, Topology optimization for negative permeability metamaterials using level-set algorithm, Acta Mater. 59 (2011) 2624–2636. [18] Y.M. Xie, G.P. Steven, A simple evolutionary procedure for structural optimization, Comput. Struct. 49 (1993) 885–896. [19] Y.M. Xie, G.P. Steven, Evolutionary Structural Optimization, Springer, London, 1997. [20] X. Huang, Y.M. Xie, Evolutionary Topology Optimization of Continuum Structures: Methods and Applications, John Wiley & Sons, Chichester, 2010. [21] X. Huang, Y.M. Xie, Convergent and mesh-independent solutions for the bi-directional evolutionary structural optimization method, Finite Elem. Anal. Des. 43 (14) (2007) 1039–1049. [22] X. Huang, Y.M. Xie, Bi-directional evolutionary topology optimization of continuum structures with one or multiple materials, Comput. Mech. 43 (3) (2009) 393–401. [23] D.C. Dobson, S.J. Cox, Maximizing band gaps in two-dimensional photonic crystals, SIAM J. Appl. Math. 59 (1999) 2108–2120. [24] S.J. Cox, D.C. Dobson, Band structure optimization of two-dimensional photonic crystals in H-polarization, J. Comput. Phys. 158 (2000) 214–224. [25] L. Shen, Z. Ye, S. He, Design of two-dimensional photonic crystals with large absolute band gaps using a genetic algorithm, Phys. Rev. B 68 (2003) 035109. [26] S. Preble, M. Lipson, Two-dimensional photonic crystals designed by evolutionary algorithms, Appl. Phys. Lett. 86 (2005) 061111. [27] H. Li, L. Jiang, W. Jia, H. Qiang, X. Li, Genetic optimization of two-dimensional photonic crystals for large absolute band-gap, Opt. Commun. 282 (2009) 3012–3017. [28] H.W. Dong, X.X. Su, Y.S. Wang, C. Zhang, Topological optimization of two-dimensional phononic crystals based on the finite element method and genetic algorithm, Struct. Multidiscip. Optim. 50 (2014) 593–604. [29] C.Y. Kao, S. Osher, E. Yablonovitch, Maximizing band gaps in two-dimensional photonic crystals by using level set methods, Appl. Phys. B 81 (2005) 235–244. [30] L. He, C.-Y. Kao, S. Osher, Incorporating topological derivatives into shape derivatives based level set methods, J. Comput. Phys. 225 (2007) 891–909. [31] A. Takezawa, M. Kitamura, Phase field method to optimize dielectric devices for electromagnetic wave propagation, J. Comput. Phys. 257 (2014) 216–240. [32] O. Sigmund, J.S. Jensen, Systematic design of phononic band-gap materials and structures by topology optimization, Philos. Trans. R. Soc., Math. Phys. Eng. Sci. 361 (2003) 1001–1019. [33] P.I. Borel, A. Harpøth, L.H. Frandsen, M. Kristensen, P. Shi, J.S. Jensen, O. Sigmund, Topology optimization and fabrication of photonic crystal structures, Opt. Express 12 (2004) 1996–2001. [34] O. Sigmund, K. Hougaard, Geometric properties of optimal photonic crystals, Phys. Rev. Lett. 100 (2008) 153904. [35] J.S. Jensen, O. Sigmund, Topology optimization for nano-photonics, Laser Photonics Rev. 5 (2) (2011) 308–321. [36] H. Men, N.C. Nguyen, R.M. Freund, P.A. Parrilo, J. Peraire, Bandgap optimization of two-dimensional photonic crystals using semidefinite programming and subspace methods, J. Comput. Phys. 229 (2010) 3706–3725. [37] H. Men, N.C. Nguyen, R.M. Freund, K.M. Lim, P.A. Parrilo, J. Peraire, Design of photonic crystals with multiple and combined band gaps, Phys. Rev. E 83 (2011) 046703. [38] X. Huang, Z.H. Zuo, Y.M. Xie, Evolutionary topology optimization of vibrating continuum structures for natural frequencies, Comput. Struct. 88 (2010) 357–364. [39] X. Huang, Y. Li, S.W. Zhou, Y.M. Xie, Topology optimization of compliant mechanisms with desired structural stiffness, Eng. Struct. 79 (2014) 13–21. [40] X. Huang, Y.M. Xie, Evolutionary topology optimization of geometrically and materially nonlinear structures under prescribed design load, Struct. Eng. Mech. 34 (5) (2010) 581–595. [41] X. Huang, Y.M. Xie, G. Lu, Topology optimization of energy absorption structures, Int. J. Crashworthiness 12 (6) (2007) 663–675. [42] X. Huang, A. Radman, Y.M. Xie, Topological design of microstructures of cellular materials for maximum bulk or shear modulus, Comput. Mater. Sci. 50 (2011) 1861–1870. [43] X. Huang, Y.M. Xie, B. Jia, Q. Li, S.W. Zhou, Evolutionary topology optimization of periodic composites for extremal magnetic permeability and electrical permittivity, Struct. Multidiscip. Optim. 46 (2012) 385–398. [44] C. Kittel, Introduction to Solid State Physics, Wiley, New York, 1996. [45] J. Jin, The Finite Element Method in Electromagnetics, 2nd ed., John Wiley & Sons, New York, 2002. [46] A.P. Seyranian, E. Lund, N. Olhoff, Multiple eigenvalues in structural optimization problems, Struct. Optim. 8 (4) (1994) 207–227. [47] J.S. Jensen, N.L. Pedersen, On maximal eigenfrequency separation in two-materials structures: the 1D and 2D scalar cases, J. Sound Vib. 289 (2006) 967–986. [48] J. Du, N. Olhoff, Topological design of freely vibrating continuum structures for maximum values of simple and multiple eigenfrequencies and frequency gaps, Struct. Multidiscip. Optim. 34 (2007) 91–110.