Robust topology optimization of three-dimensional photonic-crystal band-gap structures

We perform full 3D topology optimization (in which"every voxel"of the unit cell is a degree of freedom) of photonic-crystal structures in order to find optimal omnidirectional band gaps for various symmetry groups, including fcc (including diamond), bcc, and simple-cubic lattices. Even without imposing the constraints of any fabrication process, the resulting optimal gaps are only slightly larger than previous hand designs, suggesting that current photonic crystals are nearly optimal in this respect. However, optimization can discover new structures, e.g. a new fcc structure with the same symmetry but slightly larger gap than the well known inverse opal, which may offer new degrees of freedom to future fabrication technologies. Furthermore, our band-gap optimization is an illustration of a computational approach to 3D dispersion engineering which is applicable to many other problems in optics, based on a novel semidefinite-program formulation for nonconvex eigenvalue optimization combined with other techniques such as a simple approach to impose symmetry constraints. We also demonstrate a technique for \emph{robust} topology optimization, in which some uncertainty is included in each voxel and we optimize the worst-case gap, and we show that the resulting band gaps have increased robustness to systematic fabrication errors.


Introduction
In this paper, we present the first fully three-dimensional (3D) robust topology optimization (in which every voxel is a degree of freedom) of complete photonic band gaps in 3D photonic crystals, in contrast to earlier band-gap topology-optimization work [1][2][3][4][5][6][7] that was limited to two dimensions (2D) and did not address robustness to manufacturing defects. Our results in §3.1 confirm that, for diamond symmetry, known "hand-designed" 3D crystal structures [8], appear to be close to optimal with respect to the fractional band gap. However, the optimization also appears to discover some previously unknown structures for other symmetry groups, including a new fcc-symmetry structure that has a larger gap for the same bands than the well-known "inverse-opal" design [9] (although its gap is still smaller than diamond). On the one hand, cases where the optimization yields structures that are reminiscent of previous hand designs with only slightly larger gaps, despite searching a large space of continuously varying structures without regard for ease of fabrication, suggest limitations on the prospects for future improvements in band-gap sizes with conventional dielectric materials, in which we find gaps for index contrasts ≥ 1.9 : 1 (in §3.3) similar to previous authors [10]. On the other hand, the ability of optimization to discover new designs for a given symmetry group may inspire exploration of new fabrication technologies that are better suited to those topologies (or a distorted version thereof) than to previous structures. Moreover, our optimization approach (in §2) builds on and illustrates our recent subspace and semidefinite-program (SDP) formulation [7] that is applicable to many other dispersion-relation design problems. Our approach also demonstrates our new "fabrication adaptivity" (FA) algorithm [11] for robust topology optimization (see §3.4)-optimization of the "worst case" subject to manufacturing and other uncertainties-which differs from previous work on robust optimization in electromagnetism [12][13][14][15][16][17][18][19][20] by taking into account the material constraints that arise in topology optimization. We also describe a new technique for imposing complex symmetry constraints (e.g. diamond symmetry) which retains the simplicity of an orthogonal grid of degrees of freedom (in §2), and allows us to explore a much wider variety of symmetry groups than were considered in most previous topology-optimization work.
Photonic band gaps are frequency ranges ∆ω in which there are no propagating electro-magnetic waves, and have many potential applications, e.g, waveguides, filters, resonant cavities, etc. Band gaps can be achieved in a variety of periodic dielectric structures (photonic crystals) in 3D. The earliest 3D structures were designed by hand or by optimizing over a few parameters [8,10,[21][22][23][24][25][26][27]. In contrast to few-parameter optimization, topology optimization [28] (in which the geometry is completely arbitrary, constrained only by the spatial resolution) involves qualitatively different computational methods and geometric parameterizations.  [1][2][3][4]32] or fractional gaps [5,6]. Although most previous works maximize the absolute gap ∆ω, the fractional gap size ∆ω/ω (whereω is the mid-gap frequency) is typically the preferred metric [10] due to its scale invariance, and optimizing ∆ω generally produces a sub-optimal fractional gap. One major difficulty with frequency-gap optimization, and with eigenvalue optimization in general, is that eigenvalues are not generally differentiable functions of the design variables when eigenvalues are degenerate, and this can cause gradient-based optimization methods to break down. Moreover, we found that these breakdowns are particularly problematic in 3D band structures where many accidental degeneracies can arise. Methods like the generalized gradient have been proposed to deal with the sensitivity of the repeated eigenvalues [42,43]. Our approach instead builds on a subspace projection to reformulate the optimization problem as a convex semidefinite program, in which eigenvalue sensitivities are not explicitly required. Schematically, the process is depicted in Fig. 1: we describe the unit cell by discretized ε i voxels (symmetrized in §2), solve it to find the band structure λ (k, ε) (λ = (ω/c) 2 are the Maxwell eigenvalues [10]), and then optimize the fractional gap f (ε) := ∆λ /λ by a sequence of SDPs. This nominal optimization problem is denoted by P N . (We show in §2 that optimizing ∆λ /λ is equivalent to optimizing ∆ω/ω, but we find the former more convenient. ) The optimization problem proposed above aims to find a nominal optimal solution ε * = arg max ε f (ε) (where ε ∈ R n parameterizes the geometry), in which the process assumes the fractional gap f is exact and deterministic, and the optimal solution ε * can be fabricated precisely. This clearly poses limitations in realistic settings where f may contain uncertainties or ε * cannot be exactly built due to fabrication errors. A methodology to address this problem is robust optimization, which is a broad category of optimization approaches that take uncertainties into account [15,[44][45][46]. Here, we adopt a maximin version of the robust optimization formulation, see P R in Fig. 1, in which the worst case is optimized, i.e.,ε * = arg max ε min ε f (ε, ε ), and the variable ε contains the uncertainties. Previous works on robust optimization either focused on convex robustification where the robust formulation retains the convex structure of the nominal problem [15,45,46], or on nonconvex robust-optimization problems with the worst-case unknowns ε residing in a space independent of ε [13,14,[16][17][18][19][20]. In the framework of topology optimization, in which each voxel is an independent degree of freedom, one important type of uncertainty is in the value of each voxel, i.e., ε, ε are in the same space S. The robust formulation hence becomesε * δ = arg max ε∈S min ε ∈S, ε −ε ≤δ f (ε, ε ). The modeling and computational issues of this formulation were addressed in 2D photonic crystals by our previous work [11]. is often known as the forward problem, in which given the photonic crystal dielectric function (or in a discrete representation ε i ), one computes the band structure by solving the Maxwell equations. The optimal design problem (P N ), or the inverse problem, seeks to compute the optimal dielectric function ε * that maximizes the frequency band gap. The robust optimal design problem (P R ) seeks to compute a more robust optimal solutionε * δ by solving a maxmin optimization problem, in case the nominal optimum ε * is not easily fabricable.

Problem Formulations and Solution Methods
The time-harmonic Bloch-wave magnetic-fields H(r)e ik·r−ωt of the Maxwell equations solve the following eigenproblem [10], where c is the vacuum speed of light, and ∇ k := ∇+ik for each Bloch wavevector k in Brillouin zone B. In the nominal design problem P N , we seek to maximize the frequency-gap ratio g(u) between two bands ω m (ε(u), k) and ω m+1 (ε(u), k), . ( The simplest parametrization of ε, given a computational solver for (1) on a discrete grid, would simply be ε i ∈ [ε min , ε max ] on each grid point. Our implementation is similar to this approach in spirit, but uses a slightly different parametrization ε(u) for variables u ∈ D := [0, 1] N u that is described in detail below.
It turns out that the fractional gap in eigenvalues λ = (ω/c) 2 , is a monotonic function of g(u); hence it is equivalent to optimize f (u) or g(u), and we found f (u) easier to work with. It is a straightforward exercise to prove this monotonicity, e.g., by ≥ 0 for all u 1 , u 2 ∈ D from the semi-definiteness [10] of the eigenproblem. Conversely, neither g nor f are monotonic in the absolute gap ∆ω, so our optimization problem is not equivalent to optimizing ∆ω.
In principle, the band gaps must be determined from the band extrema over the entire Brillouin zone B. However, the problem simplifies for the high-symmetry structures to which we constrain ourselves below. First, B can be reduced by the rotational/mirror symmetries to an irreducible Brillouin zone (IBZ) [10]. Furthermore, it can easily be shown that the vertices and edges of the IBZ (denoted as ∂ B) are local extrema of the bands, due to the rotational symmetries. Empirically, it has been observed that the global band extrema almost always fall on ∂ B, with rare exceptions [47]. As a practical matter, one typically designs photonic band gaps by considering only ∂ B, and then the interior of B is checked a posteriori [10]. We adopt that procedure here, as well: we only consider ∂ B when maximizing f (u), and we verify afterwards that none of the optimized structures have gap edges elsewhere in B.
The optimal solution of max u f (u) often cannot be fabricated directly, for example, due to fabrication errors, technological limitations on the resolution of the fine features, or unavailable materials as a result of a non-binary solution u, i.e., 0 < u i∈I < 1. This nonfabricability of the nominal optimal design occurs in many optimization algorithms that do not explicitly incorporate in the solution robustness, and is an especially common issue in topology optimization, because the large number of degrees of freedom may make it easier to find non-robust solutions. However, large band gaps tend to be inherently robust to non-systematic errors, e.g., surface roughness [48], and so we found that it was not necessary to design that type of robustness explicitly. Instead, we modified our algorithm to ensure robustness to systematic (i.e., periodic and symmetric) errors, by solving the following maximin fabrication-adaptivity (FA) problem [11]:ũ * δ = arg max where δ is the robustness parameter and · 1 denotes the L 1 norm, and signifies the percentage range of allowable perturbation. The basic optimization algorithm used here for solving (3) is a subspace-based semidefiniteprogram (SDP) formulation previously proposed in Ref. 7. The FA optimization problem in (4) is solved using a linear-program (LP) based iterative algorithm proposed in §3 of Ref. 11. Numerical solutions to the Maxwell eigenvalue equation (1) required by both optimization problems P N and P R are computed using an efficient preconditioned block-iterative eigensolver in a planewave basis (from the MPB package [49]).
Successful optimization of the 3D structures also relies on some efficient geometry modeling. The periodic unit cell in 3D photonic crystal is a parallelepiped, but in order to restrict ourselves to high-symmetry structures (e.g. the diamond symmetry in Fig. 2), we should, in principle, have design variables ε i only in the "asymmetric unit" [50] of the cell: this is a typically wedgeshaped polyhedron whose symmetry operations (e.g., rotations, reflections, etc.) yield the entire unit cell. However, it is computationally inconvenient to define a structured mesh of design variables over such a polyhedron. Instead, we employ the following transformation. We define a 3D grid of design variables u i ∈ [0, 1] over the smallest parallelepiped U that encloses the asymmetric unit, as shown at the upper left of Fig. 2. To obtain the material ε at any point r, we first perform all symmetry operations on U, interpolate u(r) from the grid u i for each of these transformations (e.g., 48 for diamond symmetry), compute the averageū of all of these u(r), and finally obtain ε(r) = ε min +ū(ε max − ε min ). This projection procedure allows us to easily impose any desired symmetry group while maintaining the simplicity of a Cartesian grid of unknowns. Besides the geometry modeling discussed above, we also explore efficient parallel computation whenever possible. The eigensolver MPB [49] includes a distributed-memory parallel version. The optimization solver we use for the SDP and LP problems, MOSEK [51], supports multithreaded computation. In addition, the FA problem (4) is also conveniently parallelized, thanks to the numerous and independent LPs to be solved. The resulting software allow us to solve a typical 3D problem P N in §3.1 in less than an hour with 8 CPUs, and a typical 3D problem P R in §3.4 in less than a day with 8 CPUs.

Results and Discussion
For most of the structures shown below (except for §3.3), we consider two contrasting dielectric materials, where the high refractive-index material has n hi = 3.6 (similar to GaAs in the near infrared), and the low refractive-index material has n lo = 1 (air). We also consider various prescribed symmetries, for example, simple cubic (space group no. 221), face-centered cubic (no. 225), body-centered cubic (no. 214) and diamond (no. 227), as they correspond to some known structures with sizable band gaps [8,23,27].
3.1. Optimal structures: SC5, Diamond2, Diamond8, BCC2, and FCC8 Starting from any random distribution of material configuration (which generically has a negative gap f ), the optimization algorithm is always able to increase the gap size. However, we are only able to obtain gaps between certain pairs of bands (m and m + 1), similar to previous authors (and unlike 2D [4, 6,36]). Furthermore, the non-convexity means that some starting points lead to suboptimal local optima (often small or negative gaps) even for the known separable bands. Nevertheless, by repeating the optimization for many random starting points, we were able to obtain several structures with large gaps. For the diamond structures, we only found a handful of local optima, so that the large-gap structure was found for about 20% of starting points. For the simple-cubic structures, the large-gap structure was found for only about 5% of the starting points. For the face-centered cubic structures, we find at least two optima, one inverse-opal-like structure about 5% of the time and another (larger-gap) structure about 15% of the time. As explained in §2, we optimize one symmetry group at a time, where constraining the symmetry group has the essential benefit of allowing us to evaluate only the edges of the irreducible Brillouin zone.
A structure constrained to have symmetry Pm3m (space group no. 221) in simple cubic lattice, denoted "SC5", is shown in Fig. 3(a) and resembles a cubic lattice of hollow spheres connected by cylindrical "bonds." SC5 has a frequency gap of 16.26% between the 5th and the 6th bands. Structures very similar to SC5, which also had a gap of ∼ 13% (for 3.6 : 1 index contrast) between the same bands, were found by few-parameter optimization in previous work [23, 24, 27], although the previous work did not identify the possibility of improving the gap by introducing air holes in the center of each dielectric sphere. SC5 was the only case of a substantial gap (> 10%) in a cubic lattice that we found after examining many pairs of bands, and is probably quite challenging to fabricate at infrared scales.
A structure constrained to have symmetry Fd3m (space group no. 227) in face-centered cubic lattice (or commonly known as diamond symmetry), denoted "Diamond2", is shown in Fig. 3(b), and has a 30.15% gap between the 2nd and 3rd bands. The topology of this design is very reminiscent of a large number of diamond-like photonic-crystal hand designs in previous works [8], all of which had 20-30% gaps between the second and third bands. (The connectivity resembles the bond pattern in atomic diamond structures [10].) The many published variations on this structure were intended to adapt the structure to various fabrication technologies. Our results show that, even if one removes the constraint of easy fabrication, only a slight improvement on the gap size can be obtained for this symmetry and this pair of bands.
We also examined many other pairs of bands in diamond-symmetry structures, and only found one other possibility for a large gap, which is denoted "Diamond8" and is shown in Fig. 3(c). This structure has a 15.32% complete gap between the 8th and 9th bands, and does not bear any obvious resemblance to previous designs. Since the gap is smaller than in Dia-mond2, however, there seems little reason to seek a variant of Diamond8 that might be easier to fabricate.
A structure constrained to have symmetry I4 1 32 (space group no. 214) in body-centered lattice, denoted "BCC2", is shown in Fig. 3(d), and has a gap of 27.39% between the 2nd and 3rd bands. The existence of this gap has been previously reported in [24,52], but with smaller sizes than the optimum we have here. The crystal structures correspond to a family of bicontinuous cubic structures, specifically the single gyroid, which have been explored for self-assembly of large-scale photonic materials.
A structure constrained to have symmetry Fm3m (space group no. 225) in face-centered cubic lattice, denoted "FCC8", is shown in Fig. 3(e), and has a gap of 17.42% between the 8th and 9th bands. As we go to higher bands and lower-symmetry structures, such as FCC8, we find that there are many more local optima in the gap-optimization problem. For example in Fig. 4, two locally optimal structures are obtained starting from different initial solutions. With a random initial structure (initial u) like the one shown on the top left of Fig. 4(a), which has a negative gap, we often obtain the optimal crystal structure shown on the top right. About 5% of the time, however, or alternatively if we had started the optimization with the "inverse opal" structure (top left of Fig. 4(b)) [9] (a close-packed fcc lattice of air holes in a dielectric matrix), we would obtain an inverse-opal-like local optimum shown at the top right of Fig. 4(b) with a gap of 15.35%. (Note that optimization also reproduces the later discovery that it is advantageous to introduce additional air voids and pores [53] which can be modeled by overlapping spherical shells [10].) The FCC8 structure seems to be topologically distinct from the inverse opal: if we view the hollow dielectric blobs at the faces and corners of the cubes as "atoms," then the atoms in the inverse opal are connected by 8 rod-like "bonds" per atom, whereas the FCC8 structure has 12 bonds per atom. Essentially, FCC8 is an fcc lattice of small hollow spherical shells, each of which is connected to all 12 nearest neighbors via dielectric rods (and we explicitly reparameterize the structure in this fashion below).
3.2. Gap vs. geometry parametrization: SC5(r 1 , r 2 , r 3 ), Diamond2(r) and FCC8(r 1 , r 2 , r 3 ) Although we used topology optimization to exploit the maximum number of degrees of geometric freedom, the resulting SC5 and Diamond2 structures are relatively simple and can ultimately be described by a much smaller number of parameters, and the gaps are relatively insensitive to small geometric variations (such as the exact shape of the diamond "bonds"). For example, we can reproduce a structure very similar to SC5 with a sphere-and-cylinder structure characterized by three parameters (r 1 , r 2 , r 3 ), shown in Fig. 5(a). Here r 1 denotes the radius of the air sphere, r 2 denotes the radius of the dielectric sphere, and r 3 denotes the radius of the dielectric cylinder of length a − 2r 2 , where a is the lattice constant. Optimizing over those three parameters, with values given in Fig. 5(a), yields a gap that is actually slightly larger (∼ 17%) than the one discovered by topology optimization (∼ 15%), but this appears to be an artifact of the finite spatial resolution-even at the same spatial resolution, the MPB software uses a more accurate subpixel-averaging technique [49] in the case of simple geometric shapes (sphere and cylinders) than for the grid-based parameterization of topology optimization. Similarly, Diamond2 can be roughly reproduced by a diamond lattice of dielectric cylinder "bonds" parameterized by their radius r, shown in Fig. 6(a). For a radius r * = 0.1a, this reconstructed structure has a 31.56% gap, versus the 30.15% of the optimal solution, and again the difference seems to be an artifact of the finite spatial resolution and different subpixel averaging techniques in the two simulations. FCC8 can also be reconstructed by a similar sphere-and-cylinder structure parameterized by three parameters (r 1 , r 2 , r 3 ), shown in Fig. 7(a). The main advantage of this kind of re-parameterization of the topology-optimized structure is that it makes the results easier to communicate: anyone can use the reconstructed parameters to reproduce our results, whereas reproducing the topology-optimized structure would be difficult without access to the electronic data files.   Photonic band-gap sizes are known to be dependent on the refractive-index contrast of the two constituent dielectric materials. For similar structures, one generally expects the optimal gap to reduce monotonically as the index contrast decreases (e.g. see Appendix C in Ref. 10), although the optimal parameters will also change with index contrast. An important question is the minimum index contrast for which a band gap is possible. To answer this and similar questions, we plot the optimal gap as a function of index contrast in Fig. 8 for different symmetries and pairs of bands. The smallest index contrast for which we found a gap was n hi /n lo = 1.9, for the Diamond2 structure, very similar to previous results for hand design of this type of structure [8,10,54,55]. The optimal structures at the highest and lowest index contrasts are shown as insets of Fig. 8, and the general trend is that at higher index contrasts, the dielectric veins become thinner (a well known phenomenon dating back to quarter-wave stacks in 1D [10]) and more regular (cylindrical or spherical) in shape.
3.4. Gap vs. FA parameter: Diamond2(δ ) As discussed in Sec. 2, we also considered a robust/fabrication-adaptive (FA) version of our gap-optimization problem, which maximizes the worst-case gap with respect to uncertainty in the parameters u. In this section, we present results from FA optimization of the Diamond2 structure as a function of the amount of uncertainty δ (the mean absolute error in u(r) at each point). The δ = 0 solution is equivalent to the original ("nominal") optimization problem, and with increasing δ one expect the optimum worst-case to have a smaller gap. However, if one simply added uncertainty to the nominal design, one would expect its gap to decrease faster with δ than the FA designs (which are redesigned for each δ ). As discussed in Sec. 2, increasing δ represents increasing systematic errors in the structure (i.e. errors that are periodic and symmetrical). Precisely this behavior is shown in Fig. 9. The top-most curve is the gap size of the FAoptimized structure as a function of δ . In particular, we letũ * δ denote the geometry parameters found in solution to the FA optimization problem (4) for a given δ (so thatũ * 0 is the nominal optimum). The top curve in Fig. 9 is g(ũ * δ ), the fractional gap from the theũ * δ , and as expected the gap decreases with δ , decreasing especially rapidly for δ 2.5%. This represents a tradeoff between robustness to increasing uncertainties versus performance (gap size). The insets show the corresponding structures: similar to our results in two dimensions [11], the optimal dielectric veins are thinner in the presence of systematic uncertainties.
After the FA-optimized structure is found, we verified that the structureũ * δ was much more robust to systematic errors than the non-robust structureũ * 0 . This is shown in the bottom two curves in Fig. 9. The middle curve shows min d g(ũ * δ + d): the worst-case gap in response to systematic errors d (of mean size δ ) in the robust structureũ * δ . Naturally this is worse than the performance ofũ * δ without errors, and the gap vanishes for δ > 3.5%. The bottom curve shows a similar but much larger degradation of the nominally optimal structureũ * 0 in response to similar perturbations, and in this case the gap typically vanishes for δ > 1.5%. For the nominal structure, for computational simplicity we only computed the mean (expected) gap size in response to randomly chosen perturbations d (in 30 trials); since the mean gap is obviously an upper bound on the worst-case gap, the degradation of the mean gap is sufficient to show that this structure is much worse than the FA structure.

Conclusion
Our results show that full 3D topology optimization of photonic band gaps is feasible, but suggest that little room remains for improving upon existing band-gap sizes even when fabrication constraints are removed. Of course, it is impossible to completely rule out the possibility that a much larger gap, or a much lower index contrast, is possible, for two reasons. First, the nonconvexity of the optimization problem makes it difficult to prove that one has obtained the global optimum in such a large parameter space, although the fact that we obtain only a limited number of local optima from a large number of random starting points suggests that a global optimum has been found. Second, our optimization only searches one symmetry group and one pair of bands at a time, and future work could continue searching more groups and band pairs-our method of imposing a symmetry group, via projections, makes this particularly easy (as illustrated by the variety of space groups we were able to explore in this paper, in contrast to previous topology-optimization work). In principle, one could extend our approach in order to avoid imposing the symmetry group, by making the lattice vectors degrees of freedom and omitting the symmetry projection of the grid. This would require a much more expensive calculation because it would entail searching the entire Brillouin zone (or equivalently the entire unit cell in k space), but this may be feasible with a supercomputer-scale calculation or by a more clever method that performs an inner optimization over k. However, we conjecture that such a tour de force would merely confirm that the high-symmetry structures are optimal, because the requirement of an omnidirectional gap tends to favor high-symmetry structures in order to have the same gap in multiple directions. (Nevertheless, the existence of large gaps in structures such as the single gyroid, which have only moderately large symmetry groups, lends  Fig. 9. Robust optimal designs for Diamond2, whereũ * δ solves the FA problem (4) with uncertainty δ in u, andũ * 0 is the nominal optimum. Top curve: gap size g(ũ * δ ) of FA structure versus δ , showing tradeoff between gap size and robustness to greater uncertainties. Middle curve: worst-case gap min d g(ũ * δ + d) of the FA structure plus errors, which degrades the gap size. Bottom curve: much greater degradation of gap size is found if we add uncertainty to the nominal (non-robust) structureũ * 0 . (In this case, we show the mean gap size for random perturbations d, with the standard deviation shown as error bars; this is an upper bound on the worst-case degradation.) some encouragement to a more thorough search.) Regardless of the need for more band-gap designs, which is likely to be mainly driven by the discovery of new fabrication methods, the feasibility of 3D gap optimization offers the prospect of 3D topology optimization for many other dispersion-engineering problems. For example, our methods could be easily adapted to optimize superprism effects [56][57][58][59][60][61], supercollimation [62][63][64][65][66][67][68], dispersion compensation [69][70][71][72], phase matching for nonlinear optics [73][74][75], negative refraction for imaging [76][77][78] or other dispersion constraints for various mode-coupling and mode-conversion problems [79]. The SDP approach is essentially the same regardless of whether one is maximizing a band gap or minimizing the difference between ω(k) and some desired dispersion shape, and the need for robustness to fabrication variation arises in many such applications.