Band Gap Optimization of Two-Dimensional Photonic Crystals Using Semidefinite Programming and Subspace Methods
Han Men, Ngoc-Cuong Nguyen, Robert M. Freund, Pablo A. Parrilo, Jaume Peraire
BBand Gap Optimization of Two-Dimensional Photonic CrystalsUsing Semidefinite Programming and Subspace Methods ∗ H. Men † , N. C. Nguyen ‡ , R. M. Freund § , P. A. Parrilo ¶ , and J. Peraire (cid:107) October 24, 2018
Abstract
In this paper, we consider the optimal design of photonic crystal band structures for two-dimensional square lattices. The mathematical formulation of the band gap optimization prob-lem leads to an infinite-dimensional Hermitian eigenvalue optimization problem parametrized bythe dielectric material and the wave vector. To make the problem tractable, the original eigen-value problem is discretized using the finite element method into a series of finite-dimensionaleigenvalue problems for multiple values of the wave vector parameter. The resulting optimizationproblem is large-scale and non-convex, with low regularity and non-differentiable objective. Byrestricting to appropriate eigenspaces, we reduce the large-scale non-convex optimization prob-lem via reparametrization to a sequence of small-scale convex semidefinite programs (SDPs)for which modern SDP solvers can be efficiently applied. Numerical results are presented forboth transverse magnetic (TM) and transverse electric (TE) polarizations at several frequencybands. The optimized structures exhibit patterns which go far beyond typical physical intuitionon periodic media design.
The propagation of waves in periodic media has attracted considerable interest in recent years.This interest stems from the possibility of creating periodic structures that exhibit band gapsin their spectrum, i.e., frequency regions in which the wave propagation is prohibited. Bandgaps occur in many wave propagation phenomena including electromagnetic, acoustic and elasticwaves. Periodic structures exhibiting electromagnetic wave band gaps, or photonic crystals, haveproven very important as device components for integrated optics including frequency filters [11],waveguides [10], switches [20], and optical buffers [27].The optimal conditions for the appearance of gaps were first studied for one-dimensional crys-tals by Lord Rayleigh in 1887 [18]. In a one-dimensional periodic structure, one can widen theband gap by increasing the contrast in the refractive index and difference in width between thematerials. Furthermore, it is possible to create band gaps for any particular frequency by changing ∗ This research has been supported through AFOSR grant FA9550-08-1-0350 and the Singapore-MIT Alliance. † National University of Singapore, Center for Singapore-MIT Alliance, Singapore 117576, email: [email protected] ‡ MIT Department of Aeronautics and Astronautics, 77 Massachusetts Ave., Cambridge, MA 02139, USA, email: [email protected] § MIT Sloan School of Management, 50 Memorial Drive, Cambridge, MA 02142, USA, email: [email protected] ¶ MIT Department of Electrical Engineering and Computer Science, 77 Massachusetts Ave., Cambridge, MA 02139,USA, email: [email protected] (cid:107)
MIT Department of Aeronautics and Astronautics, 77 Massachusetts Ave., Cambridge, MA 02139, USA, email: [email protected] a r X i v : . [ m a t h . O C ] J u l he periodicity length of the crystal. Unfortunately, however, in two or three dimensions one canonly suggest rules of thumb for the existence of a band gap in a periodic structure, since no rigorouscriteria have yet been determined. This made the design of two- or three-dimensional crystals a trialand error process, being far from optimal. Indeed, the possibility of two- and three-dimensionallyperiodic crystals with corresponding two- and three-dimensional band gaps was not suggested until100 years after Rayleigh’s discovery of photonic band gap in one dimension, by Yablonovitch [25]and John [14] in 1987.From a mathematical viewpoint, the calculation of the band gap reduces to the solution of aninfinite-dimensional Hermitian eigenvalue problem which is parametrized by the dielectric functionand the wave vector. In the design setting, however, one wishes to know the answer to the question:which periodic structures, composed of arbitrary arrangements of two or more different materials,produce the largest band gaps around a certain frequency? This question can be rigorously ad-dressed by formulating an optimization problem for the parameters that represent the materialproperties and geometry of the periodic structure. The resulting problem is infinite-dimensionalwith an infinite number of constraints. After appropriate discretization in space and considerationof a finite set of wave vectors, one obtains a large-scale finite-dimensional eigenvalue problem whichis non-convex and is known to be non-differentiable when eigenvalue multiplicities exist. The cur-rent state-of-the-art work done on this problem falls into two broad categories. The first kind triesto find the “optimal” band structure by parameter studies – based on prescribed inclusion shapes(e.g., circular or hexagonal inclusions) [9] or fixed topology [26]. The second kind attempts to useformal topology optimization techniques [19, 7, 4, 15]. Both approaches typically use gradient-based optimization methods. While these methods are attractive and have been quite successfulin practice, the optimization processes employed explicitly compute the sensitivities of eigenval-ues with respect to the dielectric function, which are local subgradients for such non-differentiableproblem. As a result, gradient-based solution methods often suffer from the lack of regularity ofthe underlying problem when eigenvalue multiplicities are present, as they typically are at or nearthe solution.In this paper we propose a new approach based on semidefinite programming (SDP) and sub-space methods for the optimal design of photonic band structure. In the last two decades, SDP hasemerged as the most important class of models in convex optimization; see [1, 2, 16, 22, 24]. SDPencompasses a huge array of convex problems as special cases, and is computationally tractable(usually comparable to least-square problems of comparable dimensions). There are three distinctproperties that make SDP very suitable for the band gap optimization problem. First, the un-derlying differential operator is Hermitian and positive semidefinite. Second, the objective andassociated constraints involve bounds on eigenvalues of matrices. And third, as explained below,we can approximate the original non-convex optimization problem by a semidefinite program forwhich SDP can be well applied, thanks to its efficiency and robustness of handling this type ofspectral objective and constraints.In our approach, we first reformulate the original problem of maximizing the band gap betweentwo consecutive eigenvalues as an optimization problem in which we optimize the gap in eigenvaluesbetween two orthogonal subspaces. The first eigenspace consists of eigenfunctions correspondingto eigenvalues below the band gap, whereas the second eigenspace consists of eigenfunctions whoseeigenvalues are above the band gap. In this way, the eigenvalues are no longer present in ourformulation; however, like the original problem, the exactly reformulated optimization problemis large-scale. To reduce the problem size, we truncate the high-dimensional subspaces to onlya few eigenfunctions below and above the band gap [5, 17], thereby obtaining a new small-scaleyet non-convex optimization problem. Finally, we keep the subspaces fixed at a given decisionparameter vector and use a reparametrization of the decision variables to obtain a convex semidef-2nite optimization problem for which SDP solution methods can be effectively applied. We applythis approach to optimize band gaps in two-dimensional photonic crystals for either the transversemagnetic (TM) or the transverse electric (TE) polarizations.The rest of the paper is organized as follows. In Section 2 we introduce the governing differentialequations and the mathematical formulation of the band gap optimization problem. We then discussthe discretization process and present the subspace restriction approach. In Section 3 we introducethe semidefinite programming formulation of the band structure optimization, and lay out theoptimization steps involved in solving the problem. Numerical results are presented in Section 4for both the TE and TM polarizations in square lattices. Finally, in Section 5 we conclude withseveral remarks on anticipated future research directions. Our primary concern is the propagation of electromagnetic linear waves in periodic media, andthe design of such periodic structures, or photonic crystals, to create optimal band gaps in theirspectrum. The propagation of electromagnetic waves in photonic crystals is governed by Maxwell’sequations. The solutions to these equations are in general very complex functions of space andtime. Due to linearity however, it is possible to separate the time dependence from the spatialdependence by expanding the solution in terms of harmonic modes – any time-varying solution canalways be reconstructed by a linear combination of these harmonic modes using Fourier analysis.By considering only harmonic solutions, the problem is considerably simplified since it reduces toa series of eigenvalue problems for the spatially varying part of the solutions (eigenfunctions) andthe corresponding frequencies (eigenvalues).In the absence of sources and assuming a monochromatic wave, i.e., with magnetic field H ( r , t ) = H ( r ) e − iωt , and electric field E ( r , t ) = E ( r ) e − iωt , Maxwell’s equations can be written in the fol-lowing form: ∇ × (cid:18) ε ( r ) ∇ × H ( r ) (cid:19) = (cid:16) ωc (cid:17) H ( r ) , in R , ε ( r ) ∇ × ( ∇ × E ( r )) = (cid:16) ωc (cid:17) E ( r ) , in R , where c is the speed of light, and ε ( r ) is the dielectric function. In two dimensions, there are twopossible polarizations of the magnetic and electric fields. In TE (transverse electric) polarization,the electric field is confined to the plane of wave propagation and the magnetic field H = (0 , , H )is perpendicular to this plane. In contrast, in TM (transverse magnetic) polarization, the magneticfield is confined to the plane of wave propagation and the electric field E = (0 , , E ) is perpendicularto this plane. In such cases, the Maxwell’s equations can be reduced to scalar eigenvalue problemsTE : − ∇ · (cid:18) ε ( r ) ∇ H ( r ) (cid:19) = (cid:16) ωc (cid:17) H ( r ) , in R , (1)TM : − ∇ · ( ∇ E ( r )) = (cid:16) ωc (cid:17) ε ( r ) E ( r ) , in R . (2)Note that the reciprocal of the dielectric function is present in the differential operator for the TEcase, whereas the dielectric function is present in the right-hand side for the TM case.3igure 1: Left: A photonic crystal on a square lattice. The dashed box represents the primitiveunit cell (Ω), where a is the periodicity length of the lattice. Right: The reciprocal lattice, andthe dashed box represents the first Brillouin zone ( B ). The irreducible zone is the green triangularwedge, and its boundary is denoted by ∂ B .For two-dimensional square lattices the dielectric function satisfies ε ( r ) = ε ( r + R ), where R are the crystal lattice vectors . By applying the Bloch-Floquet theory [3, 12] for periodic eigenvalueproblems we obtain that H ( r ) = e i k · r H k ( r ) , and E ( r ) = e i k · r E k ( r ) , where H k ( r ) and E k ( r ) satisfyTE : ( ∇ + i k ) · (cid:18) ε ( r ) ( ∇ + i k ) H k ( r ) (cid:19) = (cid:16) ωc (cid:17) H k ( r ) , in Ω , (3)TM : ( ∇ + i k ) · (( ∇ + i k ) E k ( r )) = (cid:16) ωc (cid:17) ε ( r ) E k ( r ) , in Ω , (4)respectively. Thus, the effect of considering periodicity is reduced to replacing the indefinite periodicdomain by the unit cell Ω and ∇ by ∇ + i k in the original equation, where k is a wave vectorin the first Brillouin zone B . Note that the unit cell Ω and the Brillouin zone B depend on thelattice type (e.g., square or triangular lattices) as well as the crystal lattice vectors R . If we furthertake into consideration the symmetry group of the square lattice [23], we only need to consider allpossible wavevectors k on the irreducible Brillouin zone, or (under certain conditions) its boundary[13]. Figure 1 shows an example of the unit cell and the Brillouin zone for a square lattice.For notational convenience, we write the above equations in the following operator form A u = λ M u, in Ω , (5)where, for the TE case, u ≡ H k ( r ), λ ≡ ω /c , and A ( ε, k ) ≡ − ( ∇ + i k ) · (cid:18) ε ( r ) ( ∇ + i k ) (cid:19) , M ≡ I ; (6)whereas, for the TM case, u ≡ E k ( r ), λ ≡ ω /c , and A ( k ) ≡ − ( ∇ + i k ) · ( ∇ + i k ) , M ( ε ) ≡ ε ( r ) I. (7)Here I denotes the identity operator. We denote by ( u m , λ m ) the m -th pair of eigenfunction andeigenvalue of (5) and assume that these eigenpairs are numbered in ascending order: 0 < λ ≤ λ ≤ · · · ≤ λ ∞ . For a square lattice, R denotes the vectors spanned by { a e , a e } , where e and e are the unit basis vectorsand a is the periodicity length of the crystal [13]. .2 The Optimization Problem The objective in photonic crystal design is to maximize the band gap between two consecutivefrequency modes. Due to the lack of fundamental length scale in Maxwell’s equations, it can beshown that the magnitude of the band gap scales by a factor of s when the crystal is expanded bya factor of 1 /s . Therefore, it is more meaningful to maximize the gap-midgap ratio instead of theabsolute band gap [13]. The gap-midgap ratio between λ m and λ m +1 is defined as J ( ε ( r )) = inf k ∈ ∂ B λ m +1 ( ε ( r ) , k ) − sup k ∈ ∂ B λ m ( ε ( r ) , k )inf k ∈ ∂ B λ m +1 ( ε ( r ) , k ) + sup k ∈ ∂ B λ m ( ε ( r ) , k ) , where ∂ B represents the irreducible Brillouin zone boundary; see Figure 1 for example.A typical characterization of the dielectric function ε ( r ) is the distribution of two differentmaterials. Suppose that we are given two distinct materials with dielectric constants ε min and ε max where ε min < ε max . We wish to find arrangements of the materials within the unit cell Ω whichresult in maximal gap-midgap ratio. To this end, we decompose the unit cell Ω into N ε disjointsubcells K i , ≤ i ≤ N ε , such that Ω = ∪ N ε i =1 K i and K i ∩ K j = ∅ for i (cid:54) = j . Here we take thissubcell grid to be the same as the finite element triangulation of the unit cell as we are going todiscretize the continuous eigenvalue problem by the finite element method. Our dielectric function ε ( r ) takes a unique value between ε min and ε max on each subcell, namely, ε ( r ) = ε i ∈ R on K i and ε min ≤ ε i ≤ ε max . However, due to the symmetry of square lattice, we only need to definethe dielectric function ε ( r ) over part of the unit cell (1 / ε ( r ) is discretized into a finite dimensional vector ε = ( ε , . . . , ε n ε ) ∈ R n ε (with n (cid:15) ≤ N (cid:15) ) which resides in the following admissible region: Q ad ≡ { ε = ( ε , . . . , ε n ε ) ∈ R n ε : ε min ≤ ε i ≤ ε max , ≤ i ≤ n ε } . This region consists of piecewise-constant functions whose value on every subcell varies between ε min and ε max . Moreover, to render this problem computationally tractable, we replace the irreducibleBrillouin zone boundary ∂ B by a finite subset S n k = { k t ∈ ∂ B , ≤ t ≤ n k } , where k t , ≤ t ≤ n k , are wave vectors chosen along the irreducible Brillouin zone boundary. As aresult, the band gap optimization problem that maximizes the gap-midgap ratio between λ m and λ m +1 can be stated as follows:max ε J ∗ ( ε ) = min k ∈S nk λ m +1 ( ε , k ) − max k ∈S nk λ m ( ε , k )min k ∈S nk λ m +1 ( ε , k ) + max k ∈S nk λ m ( ε , k )s.t. A ( ε , k ) u j = λ j M ( ε ) u j , j = m, m + 1 , k ∈ S n k ,ε min ≤ ε i ≤ ε max , ≤ i ≤ n ε . (8)In this problem a subtle difference between TE and TM polarizations lies in the operators of theeigenvalue problem: A and M take the form of either (6) for the TE case or (7) for the TM case.In either case, note that the eigenvalue problems embedded in (8) must be addressed as part of anycomputational strategy for the overall solution of (8). We consider here the finite element method to discretize the continuous eigenvalue problem (5).This produces the following discrete eigenvalue problem A h ( ε , k ) u jh = λ jh M h ( ε ) u jh , j = 1 , . . . , N , k ∈ S n k , (9)5here A h ( ε , k ) ∈ C N ×N is a Hermitian stiffness matrix and M h ( ε ) ∈ R N ×N is a symmetric positivedefinite mass matrix. These matrices are sparse and typically very large (
N (cid:29) λ h ≤ λ h ≤ · · · ≤ λ N h .It is important to note that the dependence of the above matrices on the design parametervector ε is different for the TE and TM polarizations. In the TE case, A TE h depends on ε and M TE h does not, whereas in the TM case M TM h depends on ε and A TM h does not. More specifically, since ε ( r ) is a piecewise-constant function on Ω, the ε -dependent matrices can be expressed as A TE h ( ε , k ) = n ε (cid:88) i =1 ε i A TE h,i ( k ) , M TM h ( ε ) = n ε (cid:88) i =1 ε i M TM h,i , (10)where the matrices A TE h,i ( k ) and M TM h,i , ≤ i ≤ n ε are independent of ε . We note that A TE h ( ε , k )is linear with respect to 1 /ε i , 1 ≤ i ≤ n ε , while M TM h ( ε ) is linear with respect to ε i , 1 ≤ i ≤ n ε . The affine expansion (10) is a direct consequence of the fact that we use piecewise-constantapproximation for the dielectric function ε ( r ). (In the TE case, we will shortly change our decisionvariables to y i = 1 /ε i , 1 ≤ i ≤ n ε , so as to render A TE h affine in the variables y , . . . , y n ε .)After discretizing the eigenvalue problem (5) by the finite element method, we obtain thefollowing band gap optimization problem:max ε J h ( ε ) = min k ∈S nk λ m +1 h ( ε , k ) − max k ∈S nk λ mh ( ε , k )max k ∈S nk λ m +1 h ( ε , k ) + max k ∈S nk λ mh ( ε , k )s.t. A h ( ε , k ) u jh = λ jh M h ( ε ) u jh , j = m, m + 1 , k ∈ S n k ,ε min ≤ ε i ≤ ε max , ≤ i ≤ n ε . (11)Unfortunately, this optimization problem is non-convex; furthermore it suffers from lack of regu-larity at the optimum. The reason for this is that the eigenvalues λ mh and λ m +1 h are typically notsmooth functions of ε at points of multiplicity, and multiple eigenvalues at the optimum are typicalof structures with symmetry. As a consequence, the gradient of the objective function J ( ε ) withrespect to ε is not well-defined at points of eigenvalue multiplicity, and thus gradient-based descentmethods often run into serious numerical difficulties and convergence problems. In this section we describe our approach to solve the band gap optimization problem based on asubspace method and semidefinite programming (SDP). In our approach, we first reformulate theoriginal problem as an optimization problem in which we aim to maximize the band gap obtained byrestriction of the operator to two orthogonal subspaces. The first subspace consists of eigenfunctionsassociated to eigenvalues below the band gap, and the second subspace consists of eigenfunctionswhose eigenvalues are above the band gap. In this way, the eigenvalues are no longer explicitlypresent in the formulation, and eigenvalue multiplicity no longer leads to lack of regularity. Thereformulated optimization problem is exact but non-convex and large-scale. To reduce the problemsize, we truncate the high-dimensional subspaces to only a few eigenfunctions below and above theband gap [5, 17], thereby obtaining a new small-scale yet non-convex optimization problem. Finally,we keep the subspaces fixed at a given decision parameter vector to obtain a convex semidefiniteoptimization problem for which SDP solution methods can be efficiently applied.6 .1 Reformulation of the Band Gap Optimization Problem using Subspaces
We first define two additional decision variables: λ uh := min k ∈S nk λ m +1 h ( ε , k ) , λ (cid:96)h := max k ∈S nk λ mh ( ε , k ) , and then rewrite the original problem (11) as P : max ε ,λ uh ,λ (cid:96)h λ uh − λ (cid:96)h λ uh + λ (cid:96)h s.t. λ mh ( ε , k ) ≤ λ (cid:96)h , λ uh ≤ λ m +1 h ( ε , k ) , ∀ k ∈ S n k ,A h ( ε , k ) u mh = λ mh M h ( ε ) u mh , ∀ k ∈ S n k ,A h ( ε , k ) u m +1 h = λ m +1 h M h ( ε ) u m +1 h , ∀ k ∈ S n k ,ε min ≤ ε i ≤ ε max , i = 1 , . . . , n ε ,λ uh , λ (cid:96)h > . (12)Next, we introduce the following matrices:Φ ε ( k ) := [Φ ε (cid:96) ( k ) | Φ ε u ( k )] := [ u h ( ε , k ) . . . u mh ( ε , k ) | u m +1 h ( ε , k ) . . . u N h ( ε , k )] , where Φ ε (cid:96) ( k ) and Φ ε u ( k ) consist of the first m eigenvectors and the remaining N − m eigenvectors,respectively, of the eigenvalue problem: A h ( ε , k ) u jh = λ jh M h ( ε ) u jh , ≤ j ≤ N . We will also denote the subspaces spanned by the eigenvectors of Φ ε (cid:96) ( k ) and Φ ε u ( k ) as sp (Φ ε (cid:96) ( k ))and sp (Φ ε u ( k )), respectively.The first three sets of constraints in (12) can be represented exactly asΦ ε ∗ (cid:96) ( k )[ A h ( ε , k ) − λ (cid:96)h M h ( ε )]Φ ε (cid:96) ( k ) (cid:22) , ∀ k ∈ S n k Φ ε ∗ u ( k )[ A h ( ε , k ) − λ uh M h ( ε )]Φ ε u ( k ) (cid:23) , ∀ k ∈ S n k , where “ (cid:23) ” is the L¨owner partial ordering on symmetric matrices, i.e., A (cid:23) B if and only if A − B is positive semidefinite. We therefore obtain the following equivalent optimization problem: P : max ε ,λ uh ,λ (cid:96)h λ uh − λ (cid:96)h λ uh + λ (cid:96)h s.t. Φ ε ∗ (cid:96) ( k )[ A h ( ε , k ) − λ (cid:96)h M h ( ε )]Φ ε (cid:96) ( k ) (cid:22) , ∀ k ∈ S n k , Φ ε ∗ u ( k )[ A h ( ε , k ) − λ uh M h ( ε )]Φ ε u ( k ) (cid:23) , ∀ k ∈ S n k ,ε min ≤ ε i ≤ ε max , i = 1 , . . . , n ε ,λ uh , λ (cid:96)h > . (13)Although the reformulation P is exact, there is however a subtle difference in the interpretationof P and P : P can be viewed as maximizing the gap-midgap ratio between the two eigenvalues λ mh and λ m +1 h ; whereas P can be viewed as maximizing the gap-midgap ratio between the twosubspaces sp (Φ ε (cid:96) ( k )) and sp (Φ ε u ( k )). The latter viewpoint allows us to develop an efficient subspaceapproximation method for solving the band gap optimization problem as discussed below.7 .2 Subspace Approximation and Reduction Let us assume that we are given a parameter vector ˆ ε . We then introduce the associated matricesΦ ˆ ε ( k ) := [Φ ˆ ε (cid:96) ( k ) | Φ ˆ ε u ( k )] = [ u h (ˆ ε , k ) . . . u mh (ˆ ε , k ) | u m +1 h (ˆ ε , k ) . . . u N h (ˆ ε , k )] , where Φ ˆ ε (cid:96) ( k ) and Φ ˆ ε u ( k ) consist of the first m eigenvectors and the remaining N − m eigenvectors,respectively, of the eigenvalue problem A h (ˆ ε , k ) u jh = λ jh M h (ˆ ε ) u jh , ≤ j ≤ N . Under the presumption that sp (Φ ˆ ε (cid:96) ( k )) and sp (Φ ˆ ε u ( k )) are reasonable approximations of sp (Φ ε (cid:96) ( k ))and sp (Φ ε u ( k )) for ε near ˆ ε , we replace Φ ε (cid:96) ( k ) with Φ ˆ ε (cid:96) ( k ) and Φ ε u ( k ) with Φ ˆ ε u ( k ) to obtain P ˆ ε : max ε ,λ uh ,λ (cid:96)h λ uh − λ (cid:96)h λ uh + λ (cid:96)h s.t. Φ ˆ ε ∗ (cid:96) ( k )[ A h ( ε , k ) − λ (cid:96)h M h ( ε )]Φ ˆ ε (cid:96) ( k ) (cid:22) , ∀ k ∈ S n k , Φ ˆ ε ∗ u ( k )[ A h ( ε , k ) − λ uh M h ( ε )]Φ ˆ ε u ( k ) (cid:23) , ∀ k ∈ S n k ,ε min ≤ ε i ≤ ε max , i = 1 , . . . , n ε ,λ uh , λ (cid:96)h > . (14)Note in P ˆ ε that the subspaces sp (Φ ˆ ε (cid:96) ( k )) and sp (Φ ˆ ε u ( k )) are approximations of the subspaces sp (Φ ε (cid:96) ( k )) and sp (Φ ε u ( k )) and are no longer functions of the decision variable vector ε .Note also that the semidefinite inclusions in P ˆ ε are large-scale, i.e., the rank of either the firstor second inclusion is at least N /
2, for each k ∈ S n k , and N will typically be quite large. In orderto reduce the size of the inclusions, we reduce the dimensions of the subspaces by considering onlythe “important” eigenvectors among u h ( ε , k ) . . . u mh ( ε , k ) , u m +1 h ( ε , k ) . . . u N h ( ε , k ), namely those a k eigenvectors whose eigenvalues lie below but nearest to λ mh ( ε , k ) and those b k eigenvectorswhose eigenvalues lie above but nearest to λ m +1 h ( ε , k ), for small values of a k , b k , typically chosenin the range between 2 and 5, for each k ∈ S n k . This yields reduced matricesΦ ˆ ε a k + b k ( k ) := [Φ ˆ ε a k ( k ) | Φ ˆ ε b k ( k )] = [ u m − a k +1 h (ˆ ε , k ) . . . u mh (ˆ ε , k ) | u m +1 h (ˆ ε , k ) . . . u m + b k h (ˆ ε , k )] . Substituting Φ ˆ ε a k ( k ) in place of Φ ˆ ε ∗ (cid:96) ( k ) and Φ ˆ ε b k ( k ) in place of Φ ˆ ε ∗ u ( k ) in the formulation P ˆ ε yieldsthe following reduced optimization formulation: P ˆ ε : max ε ,λ uh ,λ (cid:96)h λ uh − λ (cid:96)h λ uh + λ (cid:96)h s.t. Φ ˆ ε ∗ a k ( k )[ A h ( ε , k ) − λ (cid:96)h M h ( ε )]Φ ˆ ε a k ( k ) (cid:22) , ∀ k ∈ S n k , Φ ˆ ε ∗ b k ( k )[ A h ( ε , k ) − λ uh M h ( ε )]Φ ˆ ε b k ( k ) (cid:23) , ∀ k ∈ S n k ,ε min ≤ ε i ≤ ε max , i = 1 , . . . , n ε ,λ uh , λ (cid:96)h > . (15)In this way the formulation P ˆ ε seeks to model only the anticipated “active” eigenvalue con-straints, in exact extension of active-set methods in nonlinear optimization. The integers a k , b k are8etermined indirectly through user-defined parameters r l >
0, and r u >
0, where we retain onlythose eigenvectors whose eigenvalues are within 100 r l % beneath λ mh (ˆ ε , k ) or whose eigenvalues arewithin 100 r u % above λ m +1 h (ˆ ε , k ). This translates to choosing a k , b k ∈ N + as the smallest integersthat satisfy λ mh (ˆ ε , k ) − λ m − a k +1 h (ˆ ε , k ) λ mh (ˆ ε , k ) ≤ r l ≤ λ mh (ˆ ε , k ) − λ m − a k h (ˆ ε , k ) λ mh (ˆ ε , k ) ,λ m + b k h (ˆ ε , k ) − λ m +1 h (ˆ ε , k ) λ m +1 h (ˆ ε , k ) ≤ r u ≤ λ m + b k +1 h (ˆ ε , k ) − λ m +1 h (ˆ ε , k ) λ m +1 h (ˆ ε , k ) . The dimensions of the resulting subspaces sp (Φ ˆ y a k ( k )) and sp (Φ ˆ y b k ( k )) are typically very small( a k , b k ∼ , . . . , P ˆ ε has significantly smaller semidefinite inclusions than if the full subspaceswere used. Also, the subspaces are kept fixed at ˆ ε in order to reduce the nonlinearity of theunderlying problem. Furthermore, we show below that for the TE and TM polarizations that P ˆ ε can be easily re-formulated as a linear fractional semidefinite program, and hence is solvable usingmodern interior-point methods. We now show that by a simple change of variables for each of the TE and TM polarizations, problem P ˆ ε can be converted to a linear fractional semidefinite program and hence can be further convertedto a linear semidefinite program. TE polarization.
We introduce the following new decision variable notation for convenience: y := ( y , y , . . . , y n y ) := (1 /ε , . . . , /ε n ε , λ (cid:96)h , λ uh ) , and set y min = 1 /ε max and y max = 1 /ε min . We also amend our notation to write various functionaldependencies on y instead of ε such as Φ ˆ y (cid:96) ( k ), etc. Utilizing (10), we re-write P ˆ ε for the TEpolarization as P ˆ y TE : max y y n y − y n y − y n y + y n y − s.t. Φ ˆ y ∗ a k ( k ) (cid:104)(cid:80) n y − i =1 y i A TE h,i ( k ) − y n y − M TE h (cid:105) Φ ˆ y a k ( k ) (cid:22) , ∀ k ∈ S n k , Φ ˆ y ∗ b k ( k ) (cid:104)(cid:80) n y − i =1 y i A TE h,i ( k ) − y n y M TE h (cid:105) Φ ˆ y b k ( k ) (cid:23) , ∀ k ∈ S n k ,y min ≤ y i ≤ y max , i = 1 , . . . , n y − ,y n y − , y n y > . (16)We note that the objective function is a linear fractional expression and the constraint functionsare linear functions of the variables y . Therefore P ˆ y TE is a linear fractional SDP. Using a standardhomogenization [6, 8], a linear fractional SDP can be converted to a linear SDP. Indeed, for notational simplicity consider a linear fractional optimization problem of the form max x c T xd T x subjectto b − Ax ∈ K , x ∈ K , where d T x > x and K , K are convex cones. Then this problem isequivalent to the problem max w,θ c T w subject to bθ − Aw ∈ K , w ∈ K , d T w = 1, θ ≥
0, under the elementarytransformations x ← ( w/θ ) and ( w, θ ) ← ( x/d T x, /d T x ), see [6, 8]. M polarization.
We introduce slightly different decision variable notation for convenience: z := ( z , z , . . . , z n z ) := ( ε , . . . , ε n ε , /λ (cid:96)h , /λ uh ) , and set z min = ε min and z max = ε max . Similar to the TE case, we amend our notation to writevarious functional dependencies on z instead of ε such as Φ ˆ z (cid:96) ( k ), etc. Noting that λ uh − λ (cid:96)h λ uh + λ (cid:96)h = z n z − − z n z z n z − + z n z , utilizing (10), and multiplying the semidefinite inclusions of (15) by z n z − and z n z , respectively, were-write P ˆ ε for the TM polarization as P ˆ z TM : max z z n z − − z n z z n z − + z n z s.t. Φ ˆ z ∗ a k ( k ) (cid:104) z n z − A TM h ( k ) − (cid:80) n z − i =1 z i M TM h,i (cid:105) Φ ˆ z a k ( k ) (cid:22) , ∀ k ∈ S n k , Φ ˆ z ∗ b k ( k ) (cid:104) z n z A TM h ( k ) − (cid:80) n z − i =1 z i M TM h,i (cid:105) Φ ˆ z b k ( k ) (cid:23) , ∀ k ∈ S n k ,z min ≤ z i ≤ z max , i = 1 , . . . , n z − ,z n z − , z n z > . (17)Here again the objective function is a linear fractional form and the constraint functions are linearfunctions of the variables z . Therefore P ˆ z TM is a linear fractional SDP with format similar to thatof P ˆ y TE .Since both P ˆ y TE and P ˆ z TM are linear fractional semidefinite programs, they can be solved veryefficiently by using modern interior point methods. Here we use the SDPT3 software [21] for thistask. We summarize our numerical approach for solving the band gap optimization problem of the TEpolarization in the following table. Essentially the same algorithm (with the modifications describedin the previous section) is used to solve the band gap optimization problem of the TM polarization.
We consider a two-dimensional photonic crystal confined in the computational domain of a unit cellof the square lattice, and with square domain Ω ≡ [ − , × [ − , ×
64, which yields a mesh size of h = 1 /
32 and 4096 linear square elements.The dielectric function ε is composed of two materials with dielectric constants ε min = 1 (air)and ε max = 11 . / ε i , i = 1 , , . . . , n ε ) isthus reduced to n ε = (1 + 32) × / mplementation StepsStep 1. Start with an initial guess y and an error tolerance (cid:15) tol , and set ˆ y := y . Step 2.
For each wave vector k ∈ S n k , do:Determine the subspace dimensions a k and b k .Compute the matrices Φ ˆ y a k ( k ) and Φ ˆ y b k ( k ). Step 3.
Form the semidefinite program P ˆ y TE . Step 4.
Solve P ˆ y TE for an optimal solution y ∗ . Step 5. If (cid:107) y ∗ − ˆ y (cid:107) ≤ (cid:15) tol , stop and return the optimal solution y ∗ .Else update ˆ y ← y ∗ and go to Step 2.
Table 1: Main algorithm for solving the band gap optimization problem.(16 ×
16) and dielectric function for the square lattice to aid visualization; note that the actualcomputational mesh (64 ×
64) is finer than this one. The shaded cells represent those modeled by ε , and the rest are obtained through symmetry. Furthermore, in this case, the irreducible Brillouinzone B is the triangle shown in Figure 1, with n k = 12 k -points taken along the boundary of thisregion ( ∂ B ). Band diagrams plotted in the figures below show the eigenvalues moving along theboundary of B , from Γ to X to M and back to Γ.Figure 2: An illustration of a coarse mesh (16 ×
16) and dielectric function for the square lattice.The shaded cells indicate the decision variables relating the dielectric material ( ε i , i = 1 , , . . . , n ε ).Note that the actual computational mesh (64 ×
64) is finer than this one.
Because the underlying optimization problem may have many local optima, the performance ofour method can be sensitive to the choice of the initial values of the decision variables y , whichin turn depend on the initial configuration ε . Indeed, different initial configurations do lead todifferent local optima as shown in Figure 3 for the second TE band gap and in Figure 4 for the11ourth TM band gap. Therefore, the choice of the initial configuration is important. We examinehere two different types of initial configurations: photonic crystals exhibiting band gaps at the lowfrequency spectrum and random distribution.The well-known photonic crystals (e.g., dielectric rods in air – Figure 4(a), air holes in dielectricmaterial, orthogonal dielectric veins – Figure 3(d)) exhibit band-gap structures at the low frequencyspectrum. Such a distribution seems to be a sensible choice for the initial configuration as itresembles various known optimal structures [4]. When these well-known photonic crystals areused as the initial configuration, our method easily produces the band-gap structures at the lowfrequency mode (typically, the first three TE and TM modes). On the other hand, maximizing theband gap at the high frequency mode (typically, above the first three TE and TM modes) tendsto produce more complicated structures which are very different from the known photonic crystalsmentioned above. As a result, when these photonic crystals are used as the initial configurationsfor maximizing the band gap at the high frequency mode, the obtained results are less satisfactory.Random initial configurations such as Figures 3(a) and 4(d)) have very high spatial variationand may thus be suitable for maximizing the band gap at the high frequency mode. Indeed, weobserve that random distributions often yield larger band gaps (better results) than the knownphotonic crystals for the high frequency modes. Of course, the random initialization does noteliminate the possibility of multiple local optima intrinsic to the physical problem. In view of thiseffect, we use multiple random distributions to initialize our method. In particular, we start ourmain algorithm with a number of uniformly random distributions as initial configurations to obtainthe optimal structures in our numerical results discussed below. (a) Initial crystal configuration λ % = 27.45% Γ X M Γ F r equen cy λ = ( ω / c ) (c) Optimized band structure λ % = 31.68% Γ X M Γ F r equen cy λ = ( ω / c ) (f) Optimized band structure Figure 3: Two locally optimal band gaps between λ and λ in the square lattice12 a) Initial crystal configuration λ % = 32.77% Γ X M Γ F r equen cy λ = ( ω / c ) (c) Optimized band structure λ % = 36.61% Γ X M Γ F r equen cy λ = ( ω / c ) (f) Optimized band structure Figure 4: Two locally optimal band gaps between λ and λ in the square lattice The dimensions of the subspaces sp (Φ ˆ y a k ( k )) and sp (Φ ˆ y b k ( k )) are determined indirectly by theparameters r l and r u . A good choice of r l (and r u ) is one that returns a k (cid:28) N (and b k (cid:28) N ), andat the same time includes the “important” eigenvectors to enhance convergence to an optimum. Inour numerical experiments, we choose r u = r l = 0 . a k and b k in the range of [2 , r u and r l (e.g., r l = r u = 0 . a k and b k and hence increases computational cost,does not yield fewer iterations than choosing r u = r l = 0 . With all the programs implemented in MATLAB and the computation performed on a Linux PCwith Dual Core AMD Opteron 270, 1 . λ TE , ∆ λ TE , ∆ λ TE , ∆ λ TE , ∆ λ TE , ∆ λ TE , ∆ λ TE , ∆ λ TE , ∆ λ TE , ∆ λ TE , Execution time (min) 5 . . . . . . . . . . λ TM , ∆ λ TM , ∆ λ TM , ∆ λ TM , ∆ λ TM , ∆ λ TM , ∆ λ TM , ∆ λ TM , ∆ λ TM , ∆ λ TM , Execution time (min) 1 . . . . . . . . . . Table 2: Example of computation time and the number of outer iterations of a successful runfor optimizing various band gaps, for both TE and TM polarization. Here ∆ λ T Ei,i +1 denotes thegap-midgap ratio between the i th and ( i + 1) st eigenvalue for the TE polarization. (a) Γ X M Γ F r equen cy λ = ( ω / c ) TM TE (b)
Figure 5: (a) Random starting structure with translation, rotation, and reflection symmetry, 3 × eigs function in the current implementation). Another promisingapproach is to explore mesh adaptivity and incorporate non-uniform meshing for the representationof the dielectric function, as well as the eigenvalue calculation. As further discussed in Section 5,it should be possible to significantly reduce the number of decision variables and computation costwith this approach. We start the optimization procedure with a random distribution of the dielectric, such as the oneshown in Figure 5(a). The corresponding band structures of the TE and TM fields are shownin Figure 5(b). In Figure 6, we present an example of the evolution of the crystal structure asthe optimization process progresses. (The light color indicates the low dielectric constant and thedark color denotes the high dielectric constant.) As illustrated in Figure 7, the gap-midgap ratiostarts from a negative value ( − . .
90% corresponding to the optimal configuration (Figure 7(f)) at which timethe optimization process terminates successfully. Another example of the optimization evolutionfor TE polarization is shown in Figure 8 and Figure 9, in which the gap-midgap ratio increasesfrom − .
21% to +29 . a) (b) (c)(d) (e) (f) Figure 6: The evolution of the square lattice crystal structure for optimizing the gap-midgap ratiobetween λ and λ . Γ X M Γ F r equen cy λ = ( ω / c ) (a) λ % = −2.15% Γ X M Γ F r equen cy λ = ( ω / c ) (b) λ % = −0.9856% Γ X M Γ F r equen cy λ = ( ω / c ) (c) λ % = 14.58% Γ X M Γ F r equen cy λ = ( ω / c ) (d) λ % = 32.67% Γ X M Γ F r equen cy λ = ( ω / c ) (e) λ % = 44.12% Γ X M Γ F r equen cy λ = ( ω / c ) (f) Figure 7: The corresponding band structure (of Figure 6)) and the gap-midgap ratio between λ and λ in the square lattice. 15 a) (b) (c)(d) (e) (f) Figure 8: The evolution of the square lattice crystal structure for optimizing the gap-midgap ratiobetween λ and λ . Γ X M Γ F r equen cy λ = ( ω / c ) (a) λ % = −39.63% Γ X M Γ F r equen cy λ = ( ω / c ) (b) λ % = −19.24% Γ X M Γ F r equen cy λ = ( ω / c ) (c) λ % = 9.24% Γ X M Γ F r equen cy λ = ( ω / c ) (d) λ % = 28.71% Γ X M Γ F r equen cy λ = ( ω / c ) (e) λ % = 29.23% Γ X M Γ F r equen cy λ = ( ω / c ) (f) Figure 9: The corresponding band structure (of Figure 8)) and the gap-midgap ratio between λ and λ in the square lattice. 16 a) Optimal crystal structure λ % = 34.83% Γ X M Γ F r equen cy λ = ( ω / c ) (b) Optimal band structure Figure 10: Optimization of band gap between λ and λ in the square lattice. (a) Optimal crystal structure λ % = 43.82% Γ X M Γ F r equen cy λ = ( ω / c ) (b) Optimal band structure Figure 11: Optimization of band gap between λ and λ in the square lattice.In Figures 10 through 19, we present only plots of the final optimized crystal structures andthe corresponding band structures for the 6 th through 10 th optimized band gaps for TE and TMpolarizations. We see that the optimized TM band gaps are exhibited in isolated high- ε structures,while the optimized TE band gaps appear in connected high- ε structures. This observation has alsobeen pointed out in [13] (p75) “the TM band gaps are favored in a lattice of isolated high- ε regions,and TE band gaps are favored in a connected lattice ”, and observed in [15] previously. For both TEand TM polarizations, the crystal structures become more and more complicated as we progressto higher bands. It would be very difficult to create such structures using physical intuition alone.The largest gap-midgap ratio for the TM case is 43 .
9% between the seventh and eighth frequencybands, while the largest ratio for the TE case is 44 . a) Optimal crystal structure λ % = 39.75% Γ X M Γ F r equen cy λ = ( ω / c ) (b) Optimal band structure Figure 12: Optimization of band gap between λ and λ in the square lattice. (a) Optimal crystal structure λ % = 37.93% Γ X M Γ F r equen cy λ = ( ω / c ) (b) Optimal band structure Figure 13: Optimization of band gap between λ and λ in the square lattice. (a) Optimal crystal structure λ % = 38.79% Γ X M Γ F r equen cy λ = ( ω / c ) (b) Optimal band structure Figure 14: Optimization of band gap between λ and λ in the square lattice.18 a) Optimal crystal structure λ % = 28.01% Γ X M Γ F r equen cy λ = ( ω / c ) (b) Optimal band structure Figure 15: Optimization of band gap between λ and λ in the square lattice. (a) Optimal crystal structure λ % = 44.12% Γ X M Γ F r equen cy λ = ( ω / c ) (b) Optimal band structure Figure 16: Optimization of band gap between λ and λ in the square lattice. (a) Optimal crystal structure λ % = 33.9% Γ X M Γ F r equen cy λ = ( ω / c ) (b) Optimal band structure Figure 17: Optimization of band gap between λ and λ in the square lattice.19 a) Optimal crystal structure λ % = 34.9% Γ X M Γ F r equen cy λ = ( ω / c ) (b) Optimal band structure Figure 18: Optimization of band gap between λ and λ in the square lattice. (a) Optimal crystal structure λ % = 29.4% Γ X M Γ F r equen cy λ = ( ω / c ) (b) Optimal band structure Figure 19: Optimization of band gap between λ and λ in the square lattice.20 Conclusions and Future Work
We have introduced a novel approach, based on reduced eigenspaces and semidefinite programming,for the optimization of band gaps of two-dimensional photonic crystals on square lattices. Ournumerical results convincingly show that the proposed method is very effective in producing avariety of structures with large band gaps at various frequency levels in the spectrum.Since our computational techniques make essential use of the finite element method, we antic-ipate that notions of mesh adaptivity can be easily incorporated into our approach, and thus itscomputational efficiency will be improved even further. For example, one can start with a relativelycoarse mesh and converge to a near-optimal solution, and then judiciously refine the finite elementmesh (e.g., refining elements at the interface of dielectric materials) using the current optimal solu-tion at the coarser mesh as the new initial configuration. We intend to explore this approach andreport the details and results in a forthcoming paper.The main strengths of our proposed approach to solve eigenvalue gap optimization problemis the fact that SDP-based methods do not require explicit computation of (sub-)gradients of theobjective function (which are ill-defined in the case of eigenvalue multiplicities), hence maintainingthe regularity of the formulation. The approach proposed in this paper can also be readily extendedto deal with more general problems, such as the optimization of photonic crystals in combined TEand TM fields, optimizing multiple band gaps, dealing with other types of lattices (e.g. triangular),as well as modeling and optimizing the design of three-dimensional photonic crystals.
Acknowledgments.
We are grateful to Professor Steven Johnson of the Mathematics Departmentof MIT for numerous discussions on this research. We thank Professor Lim Kian Meng of NationalUniversity of Singapore for advising and supporting of H. Men.
References [1] F. Alizadeh. Interior point methods in semidefinite programming with applications to combi-natorial optimization.
SIAM Journal on Optimization , 5(1):13–51, 1995.[2] F. Alizadeh, J. P. A. Haeberly, and M. L. Overton. Primal-dual interior-point methods forsemidefinite programming: convergence rates, stability and numerical results.
SIAM Journalon Optimization , 8(3):746–768, 1998.[3] F. Bloch. ¨Uber die quantenmechanik der elektronen in kristallgittern.
Zeitschrift f¨ur PhysikA Hadrons and Nuclei , 52(7):555–600, 1929.[4] M. Burger, S. J. Osher, and E. Yablonovitch. Inverse problem techniques for the design ofphotonic crystals.
IEICE Trans. Electron. E , 87:258–265.[5] E. Cances, C. LeBris, N. C. Nguyen, Y. Maday, A. T. Patera, and G. S. H. Pau. Feasibilityand competitiveness of a reduced basis approach for rapid electronic structure calculations inquantum chemistry. In
Proceedings of the Workshop for High-dimensional Partial DifferentialEquations in Science and Engineering (Montreal) , volume 41, pages 15–57, 2007.[6] A. Charnes and W. W. Cooper. Programming with linear functionals.
Naval Research LogisticsQuarterly , 9, 1962.[7] S. J. Cox and D. C. Dobson. Band structure optimization of two-dimensional photonic crystalsin H-polarization.
Journal of Computational Physics , 158(2):214–224, 2000.218] B. D. Craven and B. Mond. The dual of a fractional linear program.
Journal of MathematicalAnalysis and Applications , 42(3):507–512, 1973.[9] M. Doosje, B. J. Hoenders, and J. Knoester. Photonic bandgap optimization in inverted fccphotonic crystals.
Journal of the Optical Society of America B , 17(4):600–606, 2000.[10] S. Fan, J. D. Joannopoulos, J. N. Winn, A. Devenyi, J. C. Chen, and R. D. Meade. Guidedand defect modes in periodic dielectric waveguides.
Journal of the Optical Society of AmericaB , 12(7):1267–1272, 1995.[11] S. Fan, P. Villeneuve, J. Joannopoulos, and H. Haus. Channel drop filters in photonic crystals.
Optics Express , 3(1):4–11, 1998.[12] G. Floquet. Sur les equations differentielles lineaires a coefficients periodiques.
Ann. EcoleNorm. Ser , 2(12):47–89, 1883.[13] J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade.
Photonic crystals: moldingthe flow of light . Princeton university press, 2008.[14] S. John. Strong localization of photons in certain disordered dielectric superlattices.
PhysicalReview Letters , 58(23):2486–2489, 1987.[15] C. Y. Kao, S. Osher, and E. Yablonovitch. Maximizing band gaps in two-dimensional photoniccrystals by using level set methods.
Applied Physics B: Lasers and Optics , 81(2):235–244, 2005.[16] Y. Nesterov and A. Nemirovskii. Interior-point polynomial algorithms in convex programming.
SIAM studies in applied mathematics , 13, 1994.[17] G. S. H. Pau. Reduced-basis method for band structure calculations.
Physical Review E(Statistical, Nonlinear, and Soft Matter Physics) , 76(4):046704, 2007.[18] L. Rayleigh. On the maintenance of vibrations by forces of double frequency, and on the propa-gation of waves through a medium endowed with a periodic structure.
Phil. Mag , 24(147):145–159, 1887.[19] O. Sigmund and J. S. Jensen. Systematic design of phononic band-gap materials and structuresby topology optimization.
Philosophical Transactions: Mathematical, Physical and Engineer-ing Sciences , pages 1001–1019, 2003.[20] M. Soljacic, S. G. Johnson, M. Ibanescu, Y. Fink, and J. D. Joannopoulos. Optimal bistableswitching in nonlinear photonic crystals.
Physical Review E , 66:055601, 2002.[21] R. H. T¨ut¨unc¨u, K. C. Toh, and M. J. Todd. Solving semidefinite-quadratic-linear programsusing SDPT3.
Mathematical Programming , 95(2):189–217, 2003.[22] L. Vandenberghe and S. Boyd. Semidefinite programming.
SIAM review , 38(1):49–95, 1996.[23] H. Weyl.
Symmetry . Princeton Univ. Press, 1952.[24] H. Wolkowicz, R. Saigal, and L. Vandenberghe.
Handbook of semidefinite programming: theory,algorithms, and applications . Kluwer Academic Publishers, 2000.[25] E. Yablonovitch. Inhibited spontaneous emission in solid-state physics and electronics.
PhysicalReview Letters , 58(20):2059–2062, 1987. 2226] X. L. Yang, L. Z. Cai, Y. R. Wang, C. S. Feng, G. Y. Dong, X. X. Shen, X. F. Meng, andY. Hu. Optimization of band gap of photonic crystals fabricated by holographic lithography.
EPL-Europhysics Letters , 81(1):14001–14001, 2008.[27] M. F. Yanik and S. Fan. Stopping and storing light coherently.