1. 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–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–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 electromagnetic 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–28]. In contrast to few-parameter optimization, topology optimization [29] (in which the geometry is completely arbitrary, constrained only by the spatial resolution) involves qualitatively different computational methods and geometric parameterizations. Previous topology parameterizations include level-set descriptions [4, 30–32] and continuous relaxations in which ε (the dielectric permittivity) of every voxel can vary continuously in [ε_min, ε_max] [1–3,5,6,33–36]. In principle, the latter approach can yield structures with unphysical intermediate materials, but we find in practice that this does not occur in our 3D band-gap optimization here (we encountered rare counter-examples elsewhere, e.g., §4.3 in [37]), and in any case there are various regularization approaches to eliminate these artifacts if needed (e.g., the homogenization method in [38], the SIMP method introduced in [39] and [40], an alternative interpolation method proposed in [41], and the SINH method in [42]). A variety of optimization methods have been proposed for topology optimization of gaps [1–4, 33] or fractional gaps [5, 6]. Although many previous topology-optimization 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 [43, 44]. 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)² 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.)

 

Fig. 1 Schematic of optimization process of photonic-crystal band-gap structures. Computation from (a) to (b) 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 ${\tilde{\epsilon }}_{\delta }^{*}$ by solving a maxmin optimization problem, in case the nominal optimum ε^* is not easily fabricable.

Download Full Size | PDF

The optimization problem proposed above aims to find a nominal optimal solution ε^* = argmax_ε f (ε) (where ε ∈ ℝⁿ 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,45–47]. 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., ε̃^* = argmax_ε 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, 46, 47], or on nonconvex robust-optimization problems with the worst-case unknowns ε′ residing in a space independent of ε [13, 14, 16–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 ${\tilde{\epsilon }}_{\delta }^{*}={\text{argmax}}_{\epsilon \in S}{\text{min}}_{{\epsilon }^{\prime }\in S,\Vert {\epsilon }^{\prime }-\epsilon \Vert \le \delta }f(\epsilon ,{\epsilon }^{\prime })$. The modeling and computational issues of this formulation were addressed in 2D photonic crystals by our previous work [11].

2. 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],

It turns out that the fractional gap in eigenvalues λ = (ω/c)²,

In principle, the band gaps must be determined from the band extrema over the entire Brillouin zone ℬ. However, the problem simplifies for the high-symmetry structures to which we constrain ourselves below. First, ℬ 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 ∂ℬ) 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 ∂ℬ, with rare exceptions [48]. As a practical matter, one typically designs photonic band gaps by considering only ∂ℬ, and then the interior of ℬ is checked a posteriori [10]. We adopt that procedure here, as well: we only consider ∂ℬ when maximizing f (u), and we verify afterwards that none of the optimized structures have gap edges elsewhere in ℬ.

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∈ℐ < 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 [49], 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]:

The basic optimization algorithm used here for solving Eq. (3) is a subspace-based semidefinite-program (SDP) formulation previously proposed in 7. The FA optimization problem in Eq. (4) is solved using a linear-program (LP) based iterative algorithm proposed in §3 of [11]. Numerical solutions to the Maxwell eigenvalue equation Eq. (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 [50]).

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” [51] of the cell: this is a typically wedge-shaped 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 that encloses the asymmetric unit, denoted as U in 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.

 

Fig. 2 Material structure and the optimization design region in a diamond lattice. The optimization design region is limited to the rectangular parallelepiped, denoted by U in (1), and is only 1/64 of the cube. Through necessary symmetric operations, e.g., (2), ..., (48), the material ε(u) (top) is reconstructed.

Download Full Size | PDF

Besides the geometry modeling discussed above, we also explore efficient parallel computation whenever possible. The eigensolver MPB [50] includes a distributed-memory parallel version. The optimization solver we use for the SDP and LP problems, MOSEK [52], supports multithreaded computation. In addition, the FA problem Eq. (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.

3. 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, 28].

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, 37]). 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 Pm3̄m (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.3% 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, 28], 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.

 

Fig. 3 Optimized photonic crystals of prescribed symmetries (space groups no. 214, 221, 225, 227) with complete gaps between consecutive bands.

Download Full Size | PDF

A structure constrained to have symmetry Fd3̄m (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.2% 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.3% 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 Diamond2, however, there seems little reason to seek a variant of Diamond8 that might be easier to fabricate.

A structure constrained to have symmetry I4₁32 (space group no. 214) in body-centered lattice, denoted “BCC2”, is shown in Fig. 3(d), and has a gap of 27.4% between the 2nd and 3rd bands. The existence of this gap has been previously reported in [24,53], 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 Fm3̄m (space group no. 225) in face-centered cubic lattice, denoted “FCC8”, is shown in Fig. 3(e), and has a gap of 17.4% 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, shown on the top left of Fig. 4(b) as a close-packed fcc lattice of air holes in a dielectric matrix [9], we would obtain an inverse-opal–like local optimum shown at the top right of Fig. 4(b) with a gap of 15.4%. (Note that optimization also reproduces the later discovery that it is advantageous to introduce additional air voids and pores [54] 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).

 

Fig. 4 Band structures of initial and optimized photonic crystals of a face-centered cubic (fcc) lattice (no. 225). The two initial solutions and their band structures are shown on the left, and their corresponding optimal solutions and the band structures are shown on the right.

Download Full Size | PDF

3.2. Gap vs. geometry parametrization: SC5(r₁, r₂, r₃), Diamond2(r) and FCC8(r₁, r₂, r₃)

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₁, r₂, r₃), shown in Fig. 5(a). Here r₁ denotes the radius of the air sphere, r₂ denotes the radius of the dielectric sphere, and r₃ denotes the radius of the dielectric cylinder of length a − 2r₂, where a is the lattice constant. Optimizing over those three parameters, with values given in Fig. 5(a), yields a gap that is 16.7%, comparable to the one discovered by topology optimization (16.3%), with the 0.4% difference (relative to ω̄) due to differences in discretization (pixel boundaries in the projected degrees of freedom versus the parametrized objects) at the 32 × 32 × 32 resolution used here. 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 30.7% gap, versus the 30.2% of the topology optimal solution. FCC8 can also be reconstructed by a similar sphere-and-cylinder structure parameterized by three parameters (r₁, r₂, r₃), 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.

 

Fig. 5 (a) Geometry characterization of SC5 with three parameters (r₁, r₂, r₃). (b) Reconstructed photonic crystal with $\left(r_1^{*},r_2^{*},r_3^{*}\right)=\left(0.14a,0.36a,0.105a\right)$, and a frequency gap of 16.7%.

Download Full Size | PDF

 

Fig. 6 (a) Geometry characterization of Diamond2 with one parameters r. (b) Reconstructued photonic crystal with r^* = 0.1a, and a frequency gap of 30.7%.

Download Full Size | PDF

 

Fig. 7 (a) Geometry characterization of FCC8 with three parameters (r₁, r₂, r₃). (b) Reconstructued photonic crystal with $\left(r_1^{*},r_2^{*},r_3^{*}\right)=\left(0.12a,0.19a,0.08a\right)$, and a frequency gap of 17.8%.

Download Full Size | PDF

3.3. Gap vs. refractive index contrast n_hi/n_lo: SC5, Diamond2, and FCC8

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 [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,55,56]. 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.

 

Fig. 8 Gap vs. refractive-index contrast n_hi/n_lo

Download Full Size | PDF

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 expects 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 FA-optimized structure as a function of δ. In particular, we let ${\tilde{u}}_{\delta }^{*}$ denote the geometry parameters found in the solution to the FA optimization problem (4) for a given δ (so that ${\tilde{u}}_0^{*}$ is the nominal optimum). The top curve in Fig. 9 is $g\left({\tilde{u}}_{\delta }^{*}\right)$, the fractional gap from the the ${\tilde{u}}_{\delta }^{*}$, 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. We postulate that these robust structures could have been captured by formulations based on uniform or random boundary variations with controlled correlation lengths.

 

Fig. 9 Robust optimal designs for Diamond2, where ${\tilde{u}}_{\delta }^{*}$ solves the FA problem (4) with uncertainty δ in u, and ${\tilde{u}}_0^{*}$ is the nominal optimum. Top curve: gap size $g\left({\tilde{u}}_{\delta }^{*}\right)$ of FA structure versus δ, showing tradeoff between gap size and robustness to greater uncertainties. Middle curve: worst-case gap ${\text{min}}_{d}g\left({\tilde{u}}_{\delta }^{*}+d\right)$ 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 ${\tilde{u}}_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.)

Download Full Size | PDF

After the FA-optimized structure is found, we verified that the structure ${\tilde{u}}_{\delta }^{*}$ was much more robust to systematic errors than the non-robust structure ${\tilde{u}}_0^{*}$. This is shown in the bottom two curves in Fig. 9. The middle curve shows ${\text{min}}_{d}g\left({\tilde{u}}_{\delta }^{*}+d\right)$: the worst-case gap in response to systematic errors d (of mean size δ) in the robust structure ${\tilde{u}}_{\delta }^{*}$. Naturally this is worse than the performance of ${\tilde{u}}_{\delta }^{*}$ without errors, and the gap vanishes for δ > 3.5%. The insets show the corresponding worst-case structures. These structures contain many features of the robust solutions, The insets demonstrate the corresponding worst-case structures, which contain many features of the robust solutions, but with additional spatially-varying uncertainties that degrade the bandgap in the worst possible way.

The bottom curve shows a similar but much larger degradation of the nominally optimal structure u^* 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.

As a technical computational matter, one question is whether a straightforward iterative optimization of the robust problem converges: that is, whether d necessarily converges as ũ approaches an optimum. In the structures considered in this paper for Fig. 9, we found that d converges along with ũ without any special handling. However, from previous work [11] we know that this is problem-dependent: in some instances (outside of the scope of this paper) we have observed that d does not necessarily converge along with ũ, in which case we add a trust-region constraint on the values of d⁽ⁱ⁺¹⁾ at iteration (i + 1), of the form ‖d⁽ⁱ⁺¹⁾ − d⁽ⁱ⁾‖₁ ≤ ν_i, where ν_i is a sequence of descending tolerances. (This forces d to asymptotically remain in a single local worst-case optimum as we converge to the optimum ũ^*.)

4. 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 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 [57–62], supercollimation [63–69], dispersion compensation [70–73], phase matching for nonlinear optics [74–76], negative refraction for imaging [77–79] or other dispersion constraints for various mode-coupling and mode-conversion problems [80]. 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.