# Automated Discovery and Optimization of 3D Topological Photonic Crystals

Samuel Kim<EMAIL_ADDRESS>Department of Electrical Engineering and Computer
Science, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA
Thomas Christensen<EMAIL_ADDRESS>Department of Physics, Massachusetts
Institute of Technology, Cambridge, Massachusetts, USA  Steven G. Johnson
Department of Physics, Massachusetts Institute of Technology, Cambridge,
Massachusetts, USA  Department of Mathematics, Massachusetts Institute of
Technology, Cambridge, Massachusetts, USA  Marin Soljačić  Department of
Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts, USA

###### Abstract

Topological photonic crystals have received considerable attention for their
ability to manipulate and guide light in unique ways. They are typically
designed by hand based on careful analysis of their bands and mode profiles,
but recent theoretical advances have revealed new and powerful insights into
the connection between band symmetry, connectivity, and topology. Here we
propose a combined global and local optimization framework that integrates a
flexible symmetry-constrained level-set parameterization with standard
gradient-free optimization algorithms to optimize topological photonic
crystals, a problem setting where the objective function may be highly non-
convex and non-continuous. Our framework can be applied to any symmetry-
identifiable band topology, and we demonstrate its applicability to several
prominent kinds of three-dimensional band topology, namely $\Gamma$-enforced
nodal lines, Weyl points, and Chern insulators. Requiring no prior examples of
topological photonic crystals or prior knowledge on the connection between
structure and band topology, our approach indicates a path towards the
automated discovery of novel topological photonic crystal designs.

##  I Introduction

The past few decades have seen tremendous advances in optimization and inverse
design techniques for nanophotonic components and, in particular, photonic
crystals (PhCs) [1, 2, 3]. These advances have been spurred partly by the
increase in computing capacity and methodological innovations, and partly by
emerging micro- and nanofabrication capabilities that have significantly
expanded the overall design space for device structures. A central and
persisting effort has focused on the optimization of PhCs with large photonic
band gaps, with approaches including gradient-based [4, 5, 6], semi-definite
programming (SDP) [7, 8], and gradient-free [9, 10] techniques. In more recent
years, topological PhCs have emerged as a particularly exciting research
direction, and with the more recent introduction of topological band theory to
photonics [11, 12, 13], the PhC’s topology has also emerged as a quantity of
interest for design. Here, we revisit the question of optimizing photonic band
gaps in this new context, seeking the automated discovery and optimization of
large, topologically nontrivial band gaps and well-isolated, robust
topological degeneracies.

Topological PhCs are attractive due to their promised robustness for photonic
devices and the plethora of associated exotic optical phenomena [11, 12, 13].
In contrast to their electronic counterparts, topological PhCs offer the
exciting prospect of a flexible and tunable design space, with structural
morphology subject only to the limits of fabrication constraints. This design
freedom has not only endowed photonics with the underlying technological
promise of topological protection, but also positioned it as a promising
platform for the fundamental exploration of spinless topological physics.
Exemplifying this, the first experimental realization of a Chern insulator
without Landau levels— the quantum anomalous Hall effect [14, 15]—was achieved
in a gyromagnetic 2D PhC [16]. Similarly, the first experimental observation
of a Weyl point was achieved in a 3D high-index PhC [17], simultaneously with
its observation in TaAs [18, 19]. Gapless photonic topology has since grown to
encompass a wide variety of phases, including unconventional Weyl points [20,
21, 22], Dirac points [23, 24, 25], and nodal lines [26, 27, 28, 29, 30, 31,
32, 33].

These advances in photonic topology have been achieved largely by using
analogy, physical intuition, and symmetry requirements. Similarly,
optimization of associated designs have depended mainly on parameter sweeps
[34] and trial-and-error. Although recent studies have applied topology
optimization to topological PhCs [35, 36, 37], these efforts have all relied
on optimization of a continuous proxy objective to implicitly encourage
topology-associated behavior (e.g., power flow [35], directionality of Purcell
enhancement [37], or density of states [36]) in lieu of explicitly enforcing a
discrete topological constraint. In addition to requiring a different
formulation of this proxy for each topological effect, a significant downside
of this approach is the frequent need for an initial candidate for
optimization which is already topological [35, 37] to avoid designs that
maximize the proxy objective but lack the intended band topology. There are
two significant challenges to explicitly enforcing band topology. First, they
are inherently discrete and discontinuous, rendering commonly-used
optimization algorithms that assume differentiability nominally inapplicable.
Second, the conventional evaluation of topological invariants from band
holonomy [38] (i.e., parallel transport of wave functions, involving notions
of band connections and curvatures) is computationally costly since they
require simulating the full Brillouin zone (BZ) (e.g., 2D Chern number) or a
dense subset of it (e.g., Berry phase or 3D Chern vector).

In this work, we propose to use a combination of gradient-free optimization
algorithms, a symmetry-constrained level-set parameterization of the
structure, and a computationally efficient symmetry-based evaluation of band
topology to automate the discovery and optimization of novel topological PhCs.
An overview of the optimization process’ elements is shown in Fig. 1, the flow
of which is summarized in Fig. 1(a). By using a level-set function expressed
as a symmetry-constrained Fourier sum, we can parameterize even highly complex
geometries using only a low-dimensional parameter space [Fig. 1(b)], which in
turn makes both global and local gradient-free optimization algorithms
feasible. To efficiently evaluate the band topological constraints, we
incorporate the recently introduced frameworks of topological quantum
chemistry [39] and symmetry indicators [40, 41], which have also been applied
to PhCs [42, 43]. Furthermore, we use a global optimization algorithm in
conjunction with a stochastic local optimization stage to escape local optima
during optimization while also handling the non-continuous objective function
stemming from the topology of the PhC. We focus on 3D topological PhCs and
present novel examples in three settings:  (i) $\Gamma$-enforced topological
nodal lines,  (ii) ideal (i.e., frequency-isolated) Weyl points in the
interior of the BZ, and  (iii) photonic Chern insulators with chiral surface
states.  In the last setting, we find a photonic Chern insulator with the
largest known complete bandgap, and achieve this without explicitly relying on
the supercell modulation technique that underlies earlier designs. While we
focus on a handful of concrete topological properties, our method can in
principle be applied to any symmetry-identifiable band topology.

Figure 1: Overview of topological PhC optimization. (a) Flowchart of the
optimization process in each iteration. Optimization continues until a
convergence criterion or reaching a maximum number of iterations. (b) The
structure of the PhC is parameterized by a continuous vector consisting of the
geometry coefficients (the Fourier sum coefficients), the filling fraction,
and the unit cell lattice parameters. (c) Symmetry-based tools can be used to
calculate the band connectivity and topology from the high-symmetry (HS)
$\mathbf{k}$-points, providing a computationally efficient evaluation of the
topology constraint (here, exemplified for a hypothetical PhC in SG 148). (e)
If the PhC is found to have the desired topological indicator as computed from
the band symmetries, then the bandstructure is calculated along the HS
$\mathbf{k}$-lines (shown in purple in the inset). This example shows a
complete bandgap highlighted in yellow. (d) Optimization performance of
complete bandgaps in SG 27 over multiple trials, as measured by the best value
of the objective found so far as a function of iteration. Optimization
consists of a 2-step process: global optimization with multiple trials,
followed by local optimization using the best candidates from global
optimization. Colors represent different trials from random initializations.

##  II Methods

###  II.1 Photonic Crystal Parameterization

A key choice in the optimization of any system is the parameterization, i.e.,
the degrees of freedom. In the context of PhCs, this entails a
parameterization of the dielectric function $\varepsilon(\mathbf{r})$ across
the unit cell. Historically, a widely employed approach has involved
restricting the PhC structure to simple “basic” geometric shapes (e.g.,
circles, squares, spheres, or cubes of dielectric material) requiring very few
geometric parameters (e.g., radius or width), which in turn enables
optimization using gradient-free techniques or even plain parameter sweeps
[44]. Such approaches have been surprisingly effective, achieving for example
(trivial) complete three-dimensional bandgaps of nearly $35\%$ with an index
contrast of 4 in a diamond lattice of spherical holes [45]. At the opposite
end of the spectrum, topology optimization (referring to the optimization of
the geometric and topological shape of a _structure_ ; not to be confused with
topology of _bands_) is an extremely flexible approach that divides the design
region into a grid of pixels or voxels, each representing a continuous design
parameter. The associated dimensionality of the design space can be very high,
ranging from thousands to millions. To make this optimization tractable,
gradient-based algorithms [4, 5] or subspace methods combined with SDP [7, 8]
have been used, enabling optimization of PhC bandgaps in both 2D and 3D.

To represent the permittivity $\varepsilon(\mathbf{r})$ over the coordinates
$\mathbf{r}$ of the unit cell, we adopt a level-set parameterization. Level-
set functions are significantly more flexible than combinations of simple
shapes while making the problem more tractable than voxel-based
discretizations for gradient-free optimization algorithms. Concretely, we
introduce a level-set function $\phi(\mathbf{r})$ whose intersection with a
level-set offset $\Delta$ determines the boundary between regions of
permittivity $\varepsilon_{1}$ and $\varepsilon_{2}$:

$\varepsilon(\mathbf{r})=\begin{cases}\varepsilon_{2}&\operatorname{Re}\phi(\mathbf{r})>\Delta\\\
\varepsilon_{1}&\operatorname{Re}\phi(\mathbf{r})\leq\Delta\end{cases},$ (1)

Going forward, we will take $\varepsilon_{1}=1$, corresponding to vacuum or
air, and denote $\varepsilon_{2}$ simply by $\varepsilon$. The advantages of a
level-set function is that it allows a relatively low-dimensional
parameterization while retaining a versatile geometric design space. The
technique has been widely employed in photonic optimization, with level-set
functions parameterized by spherical harmonics [46], Hamilton-Jacobi
formulations [6], and eigenfunctions of a correlation function [10].

Here, exploiting the periodic nature of PhCs, we parameterize the level-set
function $\phi(\mathbf{r})$ by a finite Fourier summation of plane waves with
spatial frequencies at the reciprocal lattice vectors
$\\{\mathbf{G}\\}=\\{n_{1}\mathbf{G}_{1}+n_{2}\mathbf{G}_{2}+n_{3}\mathbf{G}_{3}\mid
n_{i}\subset\mathbb{Z}\\}$, i.e., by:

$\phi(\mathbf{r};\mathbf{c_{\mathbf{G}}})=\sum_{\\{\mathbf{G}\\}}c_{\mathbf{G}}\mathrm{e}^{\mathrm{i}\mathbf{G}\cdot\mathbf{r}},$
(2)

with complex expansion coefficients $\mathbf{c}=\\{c_{\mathbf{G}}\\}$. The
Fourier parameterization has a number of advantages:  (i) it automatically
incorporates lattice periodicity,  (ii) imposing an upper cutoff
$G_{\mathrm{max}}$ on the norm of included reciprocal lattice vectors
translates to limiting the spatial variation of the permittivity profile
(i.e., small feature sizes) and approximates feature-size constraints, and
(iii) symmetry constraints can be straightforwardly incorporated which
additionally translates to a reduction in the number of free expansion
coefficients  .

The third point—symmetry-constraints—is crucial in our context, given the
close connection between symmetry and band topology that we aim to exploit. We
briefly summarize how symmetry constraints can be imposed on the Fourier
coefficients, following Ref. 43. Concretely, the permittivity
$\varepsilon(\mathbf{r})$ must be invariant under every symmetry operation $g$
in the space group (SG) $\mathcal{G}$, such that
$g\varepsilon(\mathbf{r})=\varepsilon(g^{-1}\mathbf{r})=\varepsilon(\mathbf{r})$—or
equivalently $\phi(g^{-1}\mathbf{r})=\phi(\mathbf{r})$ for every
$g\in\mathcal{G}$. This imposes the coefficient constraints
$c_{\mathbf{G}}\mathrm{e}^{-\mathrm{i}g\mathbf{G}\cdot\mathbf{w}}=c_{g\mathbf{G}}$,
with $g$ expressed in Seitz notation $g=\\{\mathbf{W}|\mathbf{w}\\}$ with
rotation part $\mathbf{W}$ and translation part $\mathbf{w}$. These
constraints impose a set of interrelations among the Fourier coefficients,
linking $c_{\mathbf{G}}$ of symmetry-related reciprocal lattice vectors
$\mathop{\mathrm{star}}(\mathbf{G})=\\{g\mathbf{G}\mid g\in\mathcal{G}\\}$,
i.e., linking coefficients associated with distinct stars of $\mathbf{G}$.
Overall, this reduces the unconstrained Fourier summation in Eq. 2 to a
symmetry-constrained sum:

$\phi(\mathbf{r};\mathbf{c_{\text{star}(\mathbf{G})}})=\sum_{\\{\text{star}(\mathbf{G})\\}}c_{\text{star}(\mathbf{G})}\sum_{\mathbf{G}\in\text{star}(\mathbf{G})}c_{\mathbf{G}}^{\text{star}(\mathbf{G})}\mathrm{e}^{\mathrm{i}\mathbf{G}\cdot\mathbf{r}}.$
(3)

Here, $c_{\mathbf{G}}^{\mathop{\mathrm{star}}(\mathbf{G})}$ denotes fixed,
symmetry-determined coefficients among reciprocal lattice vectors from the
same star while $c_{\mathop{\mathrm{star}}(\mathbf{G})}$ denotes the remaining
free coefficients, each associated with a specific star. We only need to
optimize the free parameters $\mathbf{c_{\text{star}(\mathbf{G})}}$. Note that
the number of parameters thus depends on the SG: roughly, more symmetries
(larger group order) translate to fewer free parameters. Generally, these
parameters are complex; in the presence of inversion symmetry, however, they
can be restricted to be real parameters.

In practice, we limit the included stars of $\mathbf{G}$ to those that have a
representative element with components $n_{i}\in\\{0,\pm 1,\pm 2\\}$. In
addition, we parameterize the level-set boundary using the filling fraction
$f$ rather than the level-set offset $\Delta$, since the mapping between
$\Delta$ and $f$ is monotonic and the filling fraction is more intuitive to
interpret. This choice is also motivated by the empirical observation that a
filling fraction parameterization leads to better optimization performance
(likely because uniform sampling of $\Delta$ and
$c_{\mathop{\mathrm{star}}(\mathbf{G})}$ produces structures with normally
distributed filling fractions of small variance and mean close to 0.5, which
is an undesirable bias). We bound the values of $f$ and
$c_{\mathop{\mathrm{star}}(\mathbf{G})}$ to $f\in[0,1]$ and
$c_{\mathop{\mathrm{star}}(\mathbf{G})}\in[-1,1]$. Altogether, the symmetry-
constrained level-set Fourier parameterization of Eq. 3 allows a great deal of
shape design freedom, as illustrated e.g., in Fig. 1(b), while requiring only
very few parameters (in the range 11–69 for the problems considered here).

We also include the unit cell’s lattice vectors
$\\{\mathbf{R}_{1},\mathbf{R}_{2},\mathbf{R}_{3}\\}$ as part of the overall
structure specification, which can be parameterized (up to an immaterial
overall rotation) by their lengths $\\{a,b,c\\}$ and mutual angles
$\\{\alpha=\angle(\mathbf{R}_{2},\mathbf{R}_{3}),\beta=\angle(\mathbf{R}_{3},\mathbf{R}_{1}),\gamma=\angle(\mathbf{R}_{1},\mathbf{R}_{2})\\}$.
We exploit the scale-invariance of the dispersion-free Maxwell’s equations to
eliminate one of these parameters, setting $a=1$ without loss of generality.
The remaining 5 parameters are typically constrained further by the choice of
SG: e.g., SGs in the cubic crystal system (195–230) have no free lattice
parameters whereas SGs in the monoclinic crystal system (3–15) have 3 free
parameters ($b$, $c$, and $\beta$) [47]. As a practical matter, we restrict
the free lattice vector lengths and angles to $b,c\in[0.75,1.25]$ and
$\alpha,\beta,\gamma\in[90^{\circ},150^{\circ}]$ during optimization.
Intuitively, we expect that a nearly isotropic BZ will tend to host larger HS
bandgaps; empirically, we observe that restricting $b,c$ to near unity indeed
improves optimization convergence.

Unless otherwise specified, we adopt a scalar, fixed permittivity of
$\varepsilon=16$, roughly representative e.g., of silicon at visible
frequencies or certain ceramic-filled plastics at microwave frequencies [17].
The permittivity could in principle be included as an additional optimization
parameter; in practice, however, bandgap optimization tends to automatically
favor increased permittivity, rendering its optimization trivial.

###  II.2 Objective

In the topological settings explored in this work, we wish to maximize the
bandgap along high-symmetry (HS) $\mathbf{k}$-lines in the BZ. However,
bandgap optimization is a notably difficult problem as the bandgap objective
is not fully differentiable (especially at points of band crossing) and is a
highly non-convex problem with numerous local optima [5, 7]. Earlier works
were only able to optimize the bandgap if the structure used as an initial
point for optimization already had at least an incomplete gap [4, 5]. To
address this limitation, Cox and Dobson [5] combine gradients with
evolutionary algorithms to better search the entire parameter space, while Men
_et al._ [7] transform the objective and constraints to improve the
differentiability of the problem. However, these methods are still prone to
getting trapped in local maxima. We note that using the SDP and subspace
method developed by Men  _et al._ [8] to maximize the bandgap in SG 13, we
were not able to find a single example with a complete HS bandgap over dozens
of trials. Furthermore, in optimizing topological PhCs we introduce additional
constraints in the form of the topological indicator of the bands of interest.
Topology is inherently a discrete quantity,and so the constraints are non-
continuous, making this problem difficult to address using existing
approaches.

To address these challenges, we propose an optimization framework shown in
Fig. 1(a) which uses gradient-free optimization algorithms. This optimization
framework allows us to (1) simultaneously optimize the lattice parameters as
we are interested in SGs that do not have cubic lattices, (2) easily search
the global space which is highly non-convex and has numerous local maxima, (3)
handle non-differentiable objectives, and (4) incorporate non-continuous
constraints corresponding to the topological indicators of the bands.

In particular, we choose the objective function
$L(\omega_{n\mathbf{k}})$—i.e., the quantity we wish to maximize—as the
relative bandgap:

$\Delta\omega_{n}=2\frac{\min\limits_{\mathbf{k}}\omega_{n+1}(\mathbf{k})-\max\limits_{\mathbf{k}}\omega_{n}(\mathbf{k})}{\min\limits_{\mathbf{k}}\omega_{n+1}(\mathbf{k})+\max\limits_{\mathbf{k}}\omega_{n}(\mathbf{k})}$
(4)

where $\omega_{n}(\mathbf{k})$ is the $n$th band from the bottom.
$L(\omega_{n\mathbf{k}})$ can be in general a function of $\mathbf{k}$ across
the entire BZ. However, to make optimization of 3D PhCs computationally
efficient, we restrict the domain of $\mathbf{k}$ to the HS $\mathbf{k}$-lines
as shown in Fig. 1(e), discretized into approximately 300 equally spaced
$\mathbf{k}$-points.

For some applications, such as $\Gamma$-enforced topology, $n$ is fixed.
However, for other applications where we are agnostic to the particular band
and simply wish to maximize _a_ gap, we search over multiple bands by using
the objective function:

$\Delta\omega=\max_{n}\Delta\omega_{n}$ (5)

where $n$ can vary from 1 up to an upper limit we set based on computational
constraints.

In the context of topological PhCs, we want to constrain bands 1 through $n$
(where $n$ may be fixed or variable) to have a particular topological index.
To enforce the topology of the bands, we check the topological index in each
iteration as shown in Fig. 1(c), as symmetry indicators allow us to quickly
calculate the index from only the HS $\mathbf{k}$-points (the details of which
are explained in the section on “PhC Simulation”). If the bands are indeed
topological, we then calculate the entire band structure along the HS
$\mathbf{k}$-lines from which we set $L(\omega_{nk})=\Delta\omega_{n}$ or
$L(\omega_{nk})=\Delta\omega$, as appropriate. Otherwise, we assign a penalty
term of $-2$ to the objective, as this is the minimum possible value of the
relative bandgap. Incorporating the constraint into the objective function
allows us to use a variety of optimization algorithms that are not necessarily
designed for nonlinear constraints. The overall process in each iteration is
summarized in Fig. 1(a)

###  II.3 Optimization

As shown in Fig. 1(d), the optimization framework consists of 2 stages: a
global optimization stage to search the entire parameter space, and a local
optimization stage that takes the best candidates from the first stage and
fine-tunes them within a small region in the parameter space.

The first stage of optimization uses a randomized variant of the DIRECT-L
algorithm (which we will refer to as DIRECT-L-RAND) as implemented by NLopt
[48, 49]. The randomized variant allows us to run the algorithm over multiple
trials to take advantage of high-performance computing cluster resources. We
also use a random initial point to further introduce diversity over different
trials. In principle, the non-randomized variants, DIRECT or DIRECT-L, can be
used in place of DIRECT-L-RAND, although they are insensitive to the initial
point and different trials will generally converge to the same solution.

As shown in Fig. 3(a), the global optimization stage not only is able to
quickly find candidates with the desired topological index, but also produces
a variety of possible candidates from which we can start the local search. The
PhC geometries as well as their band structures are quite diverse,
demonstrating the effectiveness of stochastic global exploration.

Figure 2: Flowchart of the local optimization algorithm. Due to the non-
differentiable objective function, Sbplx is subject to premature convergence.
To ameliorate this, we periodically switch to ISRES to escape local maxima and
continue optimization.

The second stage is a local optimization to refine the best candidates found
from global optimization. We primarily use the Sbplx algorithm, a variant of
Nelder-Mead, as implemented by the NLopt library [48, 50]. The challenge with
applying popular local optimization algorithms to the problem of topological
PhCs is that they often model the objective and constraints inside trust
regions using linear approximations and thus require differentiability of both
the objective and constraints. However, the discrete nature of topology means
that the constraint (or objective) is non-differentiable, which leads to the
failure modes in the optimization algorithm, including getting stuck in non-
feasible regions of the parameter space, decreasing the step size to the point
where it no longer improves the objective, or even diverging and jumping to a
point far away from the initial point. Thus, to avoid the local optimization
prematurely converging, we periodically restart the search algorithm starting
from the best point found so far to reset the trust region and step size. We
also occasionally switch to ISRES [51, 48], a stochastic evolutionary-based
global optimization algorithm, to escape local maxima. The bounds of ISRES are
set to be a small hypercube around the best point found so far so that it
behaves like a stochastic local optimization. This process is shown in Fig. 2.

As seen in Fig. 3(b), the local optimization stage is effective in quickly and
significantly increasing the bandgap. The vast difference in performance
between global and local optimization stages indicates the difficulty of the
bandgap problem and suggests that the set of feasible parameters that support
complete HS bandgaps is small relative to the parameter space. The diversity
of the results from global optimization and the fact that the different trials
do not converge to the same structure during local optimization also
demonstrates that the bandgap problem likely has many local optima, many of
which do not have complete HS bandgaps. Thus, multiple trials with different
random initializations are necessary to find the globally optimal structure.

Our framework is similar to multi-start local optimization algorithms which
apply local optimization algorithms multiple times from different points to
find the global maxima in non-convex problems. The difference is that we apply
multiple trials of stochastic global optimization to pick the starting points
for local optimization rather than randomly sampling the starting points. This
is necessary due to (1) the extremely large number of local maxima in the
parameter space and (2) the topology constraints. Only a small fraction of
randomly sampled points obey the desired topology constraints (typically
$<10\%$ and often $<1\%$), so the global search is necessary to find valid
starting points with the desired topology. Thus, unlike prior work that
optimizes topological PhCs from a known topological structure [34, 37], our
method is able to automatically discover new topological PhCs.

Figure 3: Global and local optimization of PhCs in SG 13 with
$\Gamma$-enforced topology. (a–b) Evolution of the the complete HS gap, i.e.,
the optimization objective, during the (a) global and (b) local optimization
stages. Different curves correspond to different trials from random
initializations. Of the 10 trials shown in (a), the 5 best trials are used as
initialization for the local optimization stage. The three trials with best
eventual local optimization result are color-highlighted (red, purple, and
green). The unit cell evolution during local optimization for the green trial
is shown at select iterations (stars). (c–d) Band structures and unit cells
for the color-highlighted (red, purple, and green) optimization curves at the
final iteration of global (c) and local (d) optimization.

###  II.4 PhC Simulation

We use the MIT Photonics Bands (MPB) software to solve for the eigenmodes of
the PhC unit cell and calculate the band structure and associated symmetry
eigenvalues [52]. Specifically, we use a $16\times 16\times 16$ spatial
resolution during optimization and a $32\times 32\times 32$ resolution when
calculating the full BZ. For the symmetry-based evaluation of band topology,
we use the Julia packages Crystalline.jl [43], which implements the
methodology of topological quantum chemistry [39] and symmetry indicators [41,
53]; and MPBUtils.jl, which implements additional MPB-specific functionality
[54]. The choice of HS $\mathbf{k}$-paths is not unique for a particular SG;
we use the SeeK paths as defined by Hinuma  _et al._ [55] and as implemented
by Brillouin.jl [56].

We briefly summarize the salient ideas of the symmetry-based evaluation of
band topology (see, e.g., Refs. 57 and 40 for existing reviews of the general
methodology and Ref. 43 for its adaptation to PhCs). Given a selection of
bands $\\{n\\}$, we first compute the band symmetry eigenvalues
$\langle\mathbf{E}_{n\mathbf{k}}|g\mathbf{D}_{n\mathbf{k}}\rangle$ at the
special $\mathbf{k}$-points of the BZ, using the Bloch eigenmodes of the
$\mathbf{E}$\- and $\mathbf{D}$-field, over the symmetry operations $g$ of
each $\mathbf{k}$-point’s little group (i.e., the operations that leaves the
$\mathbf{k}$-point invariant, modulo reciprocal lattice vectors). From the set
of symmetry eigenvalues, we next compute the combination of irreducible
representations (irreps) that the bands transform as at each
$\mathbf{k}$-point. The overall information is aggregated into an integer-
valued symmetry vector:

$\mathbf{n}=[n_{\mathbf{k}_{1}}^{1},n_{\mathbf{k}_{1}}^{2},\ldots,n_{\mathbf{k}_{N}}^{1},n_{\mathbf{k}_{N}}^{2},\ldots,\mu],$
(6)

where $n_{\mathbf{k}_{i}}^{\alpha}$ denotes the number of bands in the
selection that transform as the $\alpha$th irrep $D_{\mathbf{k}_{i}}^{\alpha}$
of the little group at $\mathbf{k}_{i}$, and $\mu$ denotes the total number of
bands in the selection.

For a given PhC, we build up a set of symmetry vectors iteratively, starting
from the lowest bands. For each candidate symmetry vector, we test its
consistency against the set of compatibility relations imposed by the SG—i.e.,
constraints on how band symmetries can connect across the Brillouin zone [58,
59]—and if incompatible, additional bands are included until a compatible band
grouping is found. Each such compatible symmetry vector signals a set of bands
which are (pointwise) gapped from all other bands along all HS
$\mathbf{k}$-lines. The associated assignment of topology to each such
symmetry vector can be achieved using the formalism of symmetry indicators
[40, 41] and topological quantum chemistry [39]. In brief, for each SG, a
choice of a set of trivial generators $\\{\mathbf{a}_{i}\\}$ and a set of
nontrivial generators $\\{\mathbf{b}_{j}\\}_{j=1}^{d^{\text{BS}}}$ exists such
that any compatible vector can be expressed as their integer combination
according to [60]:

$\mathbf{n}=\sum_{i}c_{i}\mathbf{a}_{i}+\sum_{j}\nu_{j}\mathbf{b}_{j},$ (7)

where $c_{i}\in\mathbb{Z}$ are integer coefficients and
$\nu_{j}\in\mathbb{Z}_{\lambda_{j}}$ are coefficients from the ring of
integers modulo $\lambda_{j}$, i.e., from
$\mathbb{Z}_{\lambda_{j}}=\\{0,1,\ldots,\lambda_{j}-1\\}$. Here, $\lambda_{j}$
are the minimal integers that allow an integer-coefficient expansion of
$\lambda_{j}\mathbf{b}_{j}$ in $\\{\mathbf{a}_{i}\\}$ (i.e.,
$\lambda_{j}\mathbf{b}_{j}$ is trivial from the perspective of symmetry). The
symmetry indicator, which characterizes the symmetry-identifiable topology, is
precisely the expansion coefficients $\nu_{i}$: that is, we associate with
every $\mathbf{n}$ the symmetry indicator
$\bm{\nu}=(\nu_{1},\ldots,\nu_{d^{\text{BS}}})$ which is an element of the
indicator group
$\mathbb{Z}_{\lambda_{1}}\times\ldots\times\mathbb{Z}_{d^{\text{BS}}}$.
Nontrivial band topology corresponds to nonzero $\bm{\nu}$ and symmetry-
identifiable topology consequently exists only in SGs with indicator group
different from $\mathbb{Z}_{1}$. In practice, we compute the compatibility
relation tests as well as the symmetry indicators using a procedure described
elsewhere [61, 43]. A sketch of the application of the framework to a
hypothetical PhC in SG 148 (R$\overline{\text{3}}$) with indicator group
$\mathbb{Z}_{2}\times\mathbb{Z}_{4}$ is shown in Fig. 1(c).

We emphasize that the computational efficiency of this symmetry-based approach
is crucial to our optimization technique. With it, we can evaluate both the
band connectivity and topology using just a handful of special
$\mathbf{k}$-points (6–8 for the SGs considered here); without it, this would
require simulating a densely discretized subset of the BZ, rendering the
evaluation of the penalty term computationally prohibitive.

##  III Results and discussion

This section presents a number of optimized PhCs with nontrivial topology,
including $\Gamma$-enforced nodal lines, Weyl points, and Chern insulating
gaps.

###  III.1 $\Gamma$-enforced topological nodal lines

We start by considering the optimization of PhCs exhibiting so-called
$\Gamma$-enforced topological nodal lines, a recently introduced type of
photonic band topology that can be inferred solely from the number of bands
below the first gap along the HS $\mathbf{k}$-lines of the BZ [43]. In the
presence of time-reversal symmetry, this band topology can exist in six SGs
(13, 48, 49, 50, 68, and 86) and arises as a direct result of a uniquely
photonic zero-frequency polarization singularity at the $\Gamma$-point. It is
associated with an apparent bandgap between bands 2 and 3 along HS lines of
the BZ which, however, vanishes along nodal lines in the interior of the BZ.
Conceptually, it is analogous to filling-enforced topology [62], i.e.,
topology that can be inferred from band connectivity alone. Here, we focus on
optimizing the gap along the HS $\mathbf{k}$-lines (which we denote as the “HS
bandgap”) with the aim of obtaining “ideal” frequency-isolated nodal lines. We
explore each of the six SGs that support $\Gamma$-enforced nodal lines. In all
cases, $\Gamma$-enforced topology occurs only between bands 2 and 3, so we use
the objective in Eq. 4 with $n=2$.

Running our optimization algorithm for each of these six SGs, across 10–20
global trials and 5000 subsequent local iterations for the top 3–5 global
candidates, we find $\Gamma$-enforced topological designs in every SG.
However, only in SG 13 (P2/c; generated by inversion and a glide operation
$\\{2_{010}|0,0,\tfrac{1}{2}\\}$ in the monoclinic crystal system) are we able
to find full band gaps across the HS $\mathbf{k}$-path. Notably, SG 13 is a
(maximal) subgroup of the remaining SGs: the inability to obtain HS bandgaps
in these supergroups suggests that their additional symmetry constraints act
counter to gap formation. Additionally, as shown in Fig. 3, the optimization
framework is able to find full HS bandgaps in nearly half of the trials.

Figure 4: Nodal line location and dispersion in optimized candidates of SG 13.
(a) Nodal line frequency dispersion compared to the HS bands. The right-most
candidate exhibits a maximally-isolated nodal line. (b) The nodal lines
(color-coded by frequency) in the BZ (outlined in black) run along the $k_{y}$
direction. Purple lines represent the HS $\mathbf{k}$-lines over which the
bandgap is maximized

The best candidate we have found in SG 13 is presented in the bottom-most
panel of Fig. 3(b), which contains a $12.7\%$ HS gap between bands 2 and 3. We
can see that the structure converges to what resembles a double gyroid
structure, a well-known design for complete bandgaps [26]. The symmetry
indicator of the bottom 2 bands is $(\nu_{1},\nu_{2})=(1,1)$, indicating the
presence of a pair of nodal lines, each protected by a $\pi$-Berry phase [60].
Performing a full BZ computation over a dense $\mathbf{k}$-grid, we see a pair
of nodal lines running along the $k_{y}$ axis in Fig. 4(b), consistent with
the predictions by Song  _et al._ [60]. Unsurprisingly, the nodal lines are
far away in the BZ from the HS $\mathbf{k}$-lines along which we maximized the
bandgap, since the band crossings that form nodal lines are antithetical to
the optimization objective. In other words, optimizing the HS bandgap pushes
the nodal line away from the HS $\mathbf{k}$-lines.

The nodal line is frequency dispersive, exhibiting a relative $4.2\%$
variation across its span. While there exist isofrequency contours around the
nodal line at the nodal line frequencies due to its dispersion, there do not
exist any nodal line frequency pockets elsewhere in the BZ, making the nodal
line well isolated in frequency. Ideally, we would like to also minimize the
frequency dispersion of the nodal line [27], as it is necessary to realize the
characteristic linear DOS dispersion of nodal lines. However, this would
require formulating an objective that is a function of the entire BZ rather
than just the HS $\mathbf{k}$-lines, which is computationally prohibitive; we
leave such efforts to future work.

Because $\Gamma$-enforced topology is diagnosable from the connectivity alone
(i.e., if the appropriate bands are connected, they must be topological),
optimization methods such as the SDP solver [8] that were not formulated
specifically for topological PhCs can be applied for this setting without
requiring additional constraints. However, starting from random
initializations, the SDP solver was not able to find any complete HS gaps in
SGs 13, 48, and 68 over multiple trials and over a range of lattice bases. The
difficulty in finding bandgaps points towards the nonconvexity of the problem
and the necessity of global optimization methods.

Figure 5: Permittivity dependence of nodal line properties. (a) We compare the
HS gap attainable by simply scaling the permittivity of an original design
(brown lines) versus the HS gap attainable by re-optimizing the structure at
the new target permittivity (green lines), using the original structure as
initialization point. The starting design is that associated with Fig. 4(b),
originally optimized at $\varepsilon=16$. Substantial improvements are
attainable by re-optimizing. (b) Nodal line frequency dispersion as a function
of the HS bandgap for the re-optimized structures, showing a monotonic
decrease of nodal line dispersion with HS gap size.

#### Optimization under variable permittivity

To study the minimum index required to support an ideal nodal line, we take
the best candidate (which has a $12.7\%$ HS gap) and calculate the band
structure as a function of index for the optimized structure in Fig. 5. As
expected, the HS bandgap increases monotonically with the index contrast. The
same structure has a HS bandgap of 0 at a relative permittivity of
$\varepsilon\approx 10.6$. Alternatively, if we take the same structure and
re-optimize it using our local optimization algorithm in Fig. 2, we can
further improve the bandgap at a fixed permittivity, decreasing the minimum
permittivity required for a complete HS bandgap to $\varepsilon\approx 8.3$.
Looking at the nodal line for the optimized structures at different
permittivities in Fig. 5(b), the nodal line frequency dispersion scales
inversely with the attainable HS bandgap. In other words, the nodal line spans
a greater fraction of the bandgap in frequency as permittivity decreases,
remaining isolated inside the bandgap down to $\varepsilon=12$.

Many previous proposals for PhCs containing nodal lines either require
metallic materials or very high index-contrast dielectrics [26, 31, 33]. Our
work paves the way for demonstrating nodal lines in PhCs with lower index,
potentially enabling experimental realization in a wider range of frequencies.

###  III.2 Weyl points

Figure 6: Weyl point optimization in noncentrosymmetric SGs. Selected results
for optimized PhCs designs with Weyl points in SGs 27, 37, 81, and 82,
visualized by their dispersions along HS $\mathbf{k}$-lines (red lines,
nontrivial valence bands; blue lines, conduction bands; yellow shading, HS
gap; dashed black line, Weyl point frequency), density of states (DOS), and
Weyl point locations in the BZ. (a–b) Optimized PhC Weyl point designs in SGs
81 and 82: the unit cells of each design are shown (for SG 82, we show the
conventional unit cell, not the primitive). The BZs feature 4 ideal Weyl
points in the $k_{z}=\pi/c$ plane (with equivalent copies in the
$k_{z}=-\pi/c$ plane) in SG 81 and in the $k_{z}=0$ plane in SG 82, without
intersecting Fermi pockets. The density of states (DOS) exhibits a parabolic
frequency dependence at the Weyl point frequency. (c,d) Optimized results in
SG 82 with HS bandgaps, but non-ideal Weyl points. (c) The Weyl point
frequency lies on the valence HS band edge, preventing frequency-isolation.
(d) Although the Weyl point frequency lies in the center of the HS bandgap,
the DOS is asymmetric and non-parabolic, due to conduction band minima right
above the Weyl point frequency that do not intersect the considered HS
$\mathbf{k}$-lines. (e-f) Optimized PhC designs in SGs 27 and 37. In both
cases, the Weyl points overlap spectrally with large Fermi pockets that extend
over the interior of the BZ without intersecting any HS $\mathbf{k}$-line.

Next, we pursue designs of PhCs with frequency-isolated Weyl points, also
known as ideal Weyl points. Weyl points can arise in three different regions
of the BZ:  (i) at HS $\mathbf{k}$-points, stabilized by spatial or time-
reversal symmetries [63, 64],  (ii) along HS $\mathbf{k}$-lines, corresponding
to the intersection of bands transforming as different irreps [65, 66], and
(iii) at generic $\mathbf{k}$-points in the interior of BZ, stabilized by a
nontrivial symmetry indicator [60].  To date, proposed designs of frequency-
isolated ideal photonic Weyl points fall in the first two classes. Here, we
seek designs in the latter class, i.e., Weyl points at generic momenta. As
before, we use the bandgap metric along HS $\mathbf{k}$-lines as a proxy for
isolating the Weyl points from other bulk bands. Additionally, because we are
agnostic to the number of bands below the gap, we allow the objective
$L(\omega_{kn})$ to look over multiple bands, as in Eq. 5.

#### Space group constraints

Song  _et al._ [60] showed that for spinless, time-reversal invariant
particles in non-centrosymmetric crystals (i.e., lacking inversion centers),
all symmetry-indicated nontrivial band topology is associated with Weyl points
in the interior of the BZ. We therefore initially restrict our attention to
the 12 non-centrosymmetric SGs with nontrivial indicator group. However, only
a subset of these groups are suitable for finding well-isolated Weyl points
since the Nielsen–Nimomiya theorem requires Weyl points to occur in pairs of
opposite chirality, which in general do not need to have the same frequency
[67]. However, the existence of additional symmetries, e.g., mirrors or
rotations, can map opposite-chirality Weyl points onto each other, thereby
pinning them to the same frequency and allowing ideal single-frequency Weyl
points. Consulting the results of Ref. 60 for each SG individually, we find
that only six SGs (27, 37, 81, 82, 103, and 184) can host ideal, symmetry-
indicated Weyl points. We perform global and local optimizations for PhCs in
each of these SGs, and find designs with HS bandgaps in SGs 27 (Pcc2), 37
(Ccc2), 81 (P$\overline{\text{4}}$), and 82 (I$\overline{\text{4}}$). SGs 27
and 37 are centering-variants, both generated by a 2-fold rotation and an
orthogonal glide plane; similarly, 81 and 82 are centering-variants, both
generated by 4-fold rotoinversion.

#### Ideal Weyl point designs

We show results for a range of optimized designs from these SGs in Fig. 6. In
SGs 81 and 82 we find several designs with ideal Weyl points and without any
intersecting trivial bulk bands, i.e., without any “Fermi pockets”. Figure
6(a) highlights one such design from SG 81. The band structure has a $6.3\%$
HS gap between bands 4 and 5, with Weyl points closing the gap in the interior
of the BZ. The ideality of the Weyl points is also clearly exhibited in the
density of states (DOS), calculated using the tetrahedron method
implementation from Ref. 68, which displays the characteristic parabolic
dependence $\mathop{\mathrm{DOS}}(\omega)\propto(\omega-\omega_{0})^{2}$
associated with ideal Weyl points at frequency $\omega_{0}$. To determine the
Weyl point locations, we perform a full BZ dispersion calculation and identify
the Weyl points with the $\mathbf{k}$-points where the conduction–valence
frequency difference vanishes. We observe 4 Weyl points in the BZ at
$k_{z}=\pi/c$ (8 are shown; note, however, that the $k_{z}=\pm\pi/c$ planes
are equivalent), consistent with expectations for a symmetry indicator
$(\nu_{1},\nu_{2})=(0,1)$ [60]. Similar to the nodal line optimization, the
Weyl points are pushed away from the HS $\mathbf{k}$-lines in the $k_{z}=\pi$
plane, since the Weyl points represent a band-closing point which, if near a
HS $\mathbf{k}$-line, would act counter to the HS gap metric.

We have also found an example of a well-isolated, ideal Weyl point in SG 82,
as shown in Fig. 6(b). The bandstructure displays a gap of $5.0\%$ between
bands 6 and 7, and the DOS drops to 0 at the Weyl point frequency as expected.
In this SG, the Weyl points are fixed to the $k_{z}=0$ plane in the BZ. Again,
the Weyl points are pushed away from the HS $\mathbf{k}$-lines by the
optimization process.

#### Limitations of the high-symmetry bandgap metric

While our optimization procedure has discovered several designs for ideal Weyl
points, the HS bandgap metric is not infallible. Concretely, even when large
complete HS bandgaps are achieved, the corresponding Weyl points may not
necessarily be truly frequency-isolated. We identify two related failure
modes. First, the success of the HS bandgap metric hinges on an assumption of
minimally varying dispersions between HS lines, i.e., on the absence of non-
Weyl-related local optima of the dispersion at generic $\mathbf{k}$-points.
Although this is the tacit assumption underlying all the visualizations of the
band structure along HS lines, it is well-known that this assumption can be
violated [69, 70, 71]. Second, and relatedly, the frequency-location of the
Weyl points relative to the HS bandgap is not pinned; specifically, although
the _existence_ of Weyl points is guaranteed by topology, their frequency
_location_ is not guaranteed to fall near the center of the HS bandgap.

Figure 6(c-g) summarizes a few examples of these failure modes, where despite
obtaining a large HS gap, the optimization procedure does not yield truly
ideal Weyl points. Figure 6(c,d) shows the dispersions and DOS for two designs
in SG 82 with bandgaps of $2.0\%$ and $3.3\%$, respectively. In Fig. 6(c) the
Weyl point lies on the edge of the HS gap, illustrating the second failure
mode mentioned above. Despite this, the DOS still exhibits an approximately
parabolic dependence near the Weyl point and vanishes exactly at the Weyl
point frequency. Even when the Weyl point occurs near the center of the HS
gap, it may still be non-ideal, as illustrated by the design candidate
considered in Fig. 6(d): the DOS is both asymmetrical and distorted from the
expected parabolic shape due to the conduction band dropping close to, but not
intersecting with, the Weyl point frequency elsewhere in the BZ, in
illustration of the first failure mode mentioned earlier.

A more pronounced type of this failure mode, i.e., spectral intersection with
Fermi pockets despite an otherwise “centered” Weyl point, was widespread in
the optimized designs found in SG 27 and 37. Exemplifying this, Fig. 6(e)
shows an optimized PhC in SG 27 with a HS gap of $6.3\%$, and Fig. 6(f) shows
a PhC in SG 37 with a HS gap of $6.2\%$. In both settings, the Weyl points are
pinned to the $k_{z}=\pi$ plane (and, equivalently, $k_{z}=-\pi$ plane), and
the Weyl point frequency lies in the middle of the HS gap. Despite this,
neither case exhibits ideal Weyl points; instead, the Weyl points are occluded
by Fermi pockets that reside in the interior of the BZ, i.e., there exist
local optima of the bands inside the BZ. We visualize these Fermi pockets in
Fig. 6(e,f) by plotting the isofrequency contours at the Weyl point frequency
$\omega_{0}$.

Overall, the HS bandgap metric serves as a meaningful and largely effective
proxy for optimizing Weyl points, as evidenced by our findings of ideal Weyl
point designs in SGs 81 and 82. Overcoming the limitations of the HS gap
metric is in principle possible, e.g., by replacing it by a full BZ gap metric
(at considerable computational cost) or by supplementing the HS paths with
additional paths, either from the outset or adaptively. The DOS, which can be
computed efficiently using the tetrahedron or Gilat–Raubenheimer methods [68],
could potentially also serve as an effective proxy.

###  III.3 Chern insulators

Figure 7: Optimization of photonic Chern insulators. BZs and HS
$\mathbf{k}$-lines for (a) SG 75 and (b) SG 168, respectively. PhC structure
and band structure (complete gap highlighted in yellow) for (c) SG 75, $C=1$,
(d) SG 75, $C=2$, and (e) SG 168, $C=1$. (f) Supercell calculation for the SG
168 candidate in (e) to calculate surface states. The supercell consists of
$1\times 16\times 1$ unit cells. At the interface, the sign of the applied
applied magnetic field $\mathbf{B}$ is flipped: i.e., the orange units cells
have $\mathbf{B}=+B\hat{\mathbf{z}}$ ($C=+1$) and the blue unit cells have
$\mathbf{B}=-B\hat{\mathbf{z}}$ ($C=-1$). (g) Surface state dispersion at the
interface, extracted from supercell calculation (projected bulk states
indicated in gray).

As a final showcase, we pursue designs of time-reversal broken 3D PhCs with
nonzero Chern numbers, i.e., photonic 3D Chern insulators. There have been
numerous proposals [15, 72, 73, 74, 75] and associated experiments [16, 76,
77, 78, 79] of photonic 2D Chern insulators. The 3D generalization, however,
has seen substantially less development, with only two designs known to us
[80, 81]. Both of these designs require supercell modulations, which are
undesirable from a fabrication perspective, and feature relatively small
(${\lesssim}\,2\%$) bandgaps. 3D Chern insulators are characterized by a Chern
vector $\mathbf{C}=(C_{1},C_{2},C_{3})$ rather than a single Chern number,
with components $C_{i}$ giving the Chern number of a $\mathbf{k}$-slice
orthogonal to $\mathbf{G}_{i}$ [67]. Since the associated bands are assumed to
be fully gapped, the slice-Chern numbers $C_{i}$ are invariant with respect to
the slice position $k_{i}$.

Nontrivial Chern phases require breaking time-reversal symmetry. To achieve
this in the PhC context, we consider a gyroelectric material under a
$\hat{\mathbf{z}}$-oriented external magnetic field with a permittivity tensor
of the form [26]:

$\bm{\varepsilon}(g)=\begin{bmatrix}\varepsilon_{\perp}&\mathrm{i}g&0\\\
-\mathrm{i}g&\varepsilon_{\perp}&0\\\
0&0&\varepsilon_{\parallel}\end{bmatrix},$ (8)

where $\varepsilon_{\perp}=(\varepsilon_{\parallel}^{2}+g^{2})^{1/2}$. We take
$\varepsilon_{\parallel}=16$ and $g=12$ corresponding to a dimensionless
effective magnetic field intensity of
$|\mathbf{B}|={g}/{\varepsilon_{\parallel}}=0.75$. Note that the determinant
of $\bm{\varepsilon}(g)$ is independent of $g$ and $|\mathbf{B}|$, ensuring
that the overall band structure is not shifted as a whole by the applied
magnetic field [26].

Although the slice-Chern numbers can be computed by discretization [82], the
computational cost is prohibitive for optimization purposes. Instead, we focus
our attention to PhCs in SGs 75 (P4) and 168 (P6), where the slice-Chern
number $C_{z}$ can be computed from the bands’ $n$-fold rotation-eigenvalues
modulo $n$ [83], or equivalently, from the symmetry indicators
$\nu_{1}^{(75)}\in\mathbb{Z}_{4}$ and $\nu_{1}^{(168)}\in\mathbb{Z}_{6}$.
Furthermore, under suitable convention choices, the symmetry indicators map
directly to $C_{z}$, i.e., $C_{z}^{(75)}=\nu_{1}^{(75)}\pmod{4}$ and
$C_{z}^{(168)}=\nu_{1}^{(168)}\pmod{6}$ 111This mapping involves a fixed
convention for the choice of nontrivial generators in Eq. 7 [60].
Specifically, we take
$\mathbf{b}^{(75)}_{1}=-\hat{\mathbf{e}}_{\Gamma_{2}}+\hat{\mathbf{e}}_{\Gamma_{4}}-\hat{\mathbf{e}}_{Z_{2}}+\hat{\mathbf{e}}_{Z_{4}}$
and
$\mathbf{b}^{(168)}_{1}=+\hat{\mathbf{e}}_{H_{2}}-\hat{\mathbf{e}}_{H_{3}}+\hat{\mathbf{e}}_{K_{2}}-\hat{\mathbf{e}}_{K_{3}}-\hat{\mathbf{e}}_{L_{1}}+\hat{\mathbf{e}}_{L_{2}}-\hat{\mathbf{e}}_{M_{1}}+\hat{\mathbf{e}}_{M_{2}}$.
Equivalently, the symmetry indicators are given by
$\nu_{1}^{(75)}=2n_{M_{2}}+n_{M_{3}}-n_{M_{4}}-2n_{\Gamma_{1}}-n_{\Gamma_{3}}+n_{\Gamma_{4}}-2n_{X_{1}}\bmod
4$ and
$\nu_{1}^{(168)}=n_{\Gamma_{1}}-2n_{\Gamma_{2}}+5n_{\Gamma_{3}}+2n_{\Gamma_{4}}+3n_{\Gamma_{5}}+2n_{H_{1}}-2n_{H_{2}}-3n_{L_{1}}\bmod
6$, with irrep labels given in the CDML notation [98]. . Lacking a
corresponding $\mathbf{k}$-slice with nontrivial rotation symmetry, the
associated $C_{x,y}$ numbers are indeterminable from symmetry.

Figure 7 summarizes the results of applying our technique to the optimization
of 3D Chern insulating phases with target Chern vectors
$\mathbf{C}=(0,0,C_{z})$ for $C_{z}=1$ and $C_{z}=2$. For SG 75, we find PhC
designs with $C_{z}=1$ and a HS gap of $1.5\%$ between bands 4 and 5 [Fig.
7(c)]. Interestingly, the optimization finds a design with larger HS gap of
$2.1\%$ between bands 6 and 7 for a target Chern number of $C_{z}=2$ [Fig.
7(d)]. Conversely, in SG 168, we are unable to find a $C_{z}=2$ design with a
complete HS gap, but find a PhC design with a large $5.1\%$ $C_{z}=1$ HS gap
between bands 7 and 8 [Fig. 7(e)].

In more detail, we consider a $(010)$ interface with the applied magnetic
field’s sign flipped at the interface [Fig. 7(f)]. This flip preserves the
band structure and, more specifically, the gap; but negates the associated
Chern number, yielding a gap Chern number of $\Delta C_{z}=1-(-1)=2$. By the
bulk–boundary correspondence principle, two chiral surface states must
therefore traverse the gap. To verify this, we computed the surface band
structure across the remaining good wave numbers $(k_{x},k_{z})$ using a
supercell of $8+8$ unit cells along $y$, as shown in Fig. 7(g). A pair of
chiral surface states traverse the gap, consistent with expectations, and
effectively behave as a pair of 1D edge states at each $k_{x}$ value.
Additionally, we observe that the complete HS gap metric is highly effective
in this example, in fact equaling the full BZ gap. For computational reasons,
our $y$-supercell was implemented as periodic; as a result, a set of redundant
oppositely-traversing dispersion copies, residing at a second interface, were
manually removed from the band structure visualization for clarity.

##  IV Conclusion

We have presented a combined global and local optimization framework which
exploits several off-the-shelf optimization algorithms to effectively handle
the non-continuous and non-differentiable objective functions that arise in
the bandgap optimization of topological PhCs. We incorporate a level-set
function expressed as a symmetry-constrained Fourier sum with a relatively low
parameter-space dimensionality that makes optimization by gradient-free
methods tractable, while simultaneously instating a symmetry-constrained
search space. We have applied our method in several distinct symmetry
settings, hosting several distinct types of nontrivial band topology. In doing
so, we have discovered several novel topological photonic designs featuring
nodal lines, ideal Weyl points, and nontrivial 3D Chern insulators. To our
knowledge, these are the first proposed structures with nodal lines and Weyl
points in the interior of the BZ, a design task that would be computationally
intractable without the use of the symmetry-based topological band analysis
tools that we use here. As an example of the applicability and potential of
our framework, our best Chern insulator design has the largest known bandgap
of any 3D photonic Chern insulator, and achieves this without explicitly
incorporating the supermodulation techniques employed in previous designs.

As we also discuss in the Results section, the HS bandgap maximization
objective serves as an effective proxy for frequency-isolation of nodal lines
and Weyl points. Nevertheless, this proxy can lead to non-ideal degeneracies
that that overlap spectrally with other “trivial” bulk band features. Future
work could focus on further optimization of the frequency isolation,
especially in the local optimization stage given the computational cost of
considering the full BZ. For example, an objective function for Weyl points
could optimize the shape of the DOS directly such that it realizes the
characteristic parabolic frequency-dependence over a desired range to exclude
the possibility of Fermi pockets at or near the Weyl point frequency.
Analogously, a modified objective function for nodal lines could incorporate
terms that minimize frequency dispersion across the line while simultaneously
maximizing the bandgap to nearby Fermi pockets.

While we explore just three different kinds of band topology, the method we
present is general and can be applied to any other topological feature that
can be analyzed from band symmetry. For example, the optimization scheme is
not restricted to nodal lines and Weyl points in the interior of the BZ, but
could be easily extended to nodal lines crossing or Weyl points lying on HS
$\mathbf{k}$-lines [85], simply by adjusting the objective function or
symmetry setting accordingly. More generally, other topological features of
potential interest include the quantum spin Hall effect, higher-order
topological insulators (e.g., hinge states and corner states) [86, 87], and
symmetry-protected degeneracies more broadly [64]. Furthermore, our
exploration of optimization techniques can provide insight into bandgap
optimization of non-topological (i.e., trivial) 3D and 2D PhCs, as well as of
other quasiparticles, such as in engineered phononic crystals.

The Fourier level-set parameterization and its incorporation of a maximal
(spatial) frequency cutoff has the advantage that it generally avoids rapid
structural variations and fine details. Nevertheless, such constraints could
potentially be incorporated with higher fidelity and stricter tolerances by
exploiting standard techniques in topology optimization [88, 89]. In the
context of PhCs, robust optimization techniques using min-max formulations for
the objective function have been used to ensure satisfactory performance in
the case of over- or underetching during fabrication [90, 8]. In photonics
design more broadly, density filters and projection steps can be applied
periodically throughout optimization to remove small features and gaps [91]. A
separate but related concern of particular importance for 3D structures is the
incorporation of connectivity constraints to avoid floating structures or
holes. We note that several of our proposed Weyl point and Chern insulator
structures exhibit floating features, reflecting the absence of a connectivity
enforcement in the Fourier level-set parameterization. While enforcement of
connectivity constraints remains relatively underexplored in the 3D PhC
context—-perhaps in reflection of earlier topology optimizations of trivial
bandgaps having automatically tended towards connected designs [8]—such
constraints can in principle be explicitly incorporated [92]. On the other
hand, even PhCs with unconnected features may be realizable with emerging
fabrication techniques such as implosion fabrication, wherein dielectric or
metallic structures are embedded in a low-index hydrogel scaffold [93, 94].

In principle, any global or local optimization algorithm can be used in our
framework. Although we have explored several choices for each stage, we have
not exhaustively compared different algorithms, and further improvements may
well be achievable by simply using more optimal algorithms, e.g., algorithms
that explicitly handle non-continuous constraints. Similarly, gradient-based
optimization algorithms could be used in the local optimization stage to speed
up convergence. As a particularly exciting outlook, we note opportunities to
exploit Bayesian optimization, a model-based algorithm typically used for
global optimization of expensive non-convex functions. Recently, Bayesian
optimization using Bayesian neural networks was demonstrated to outperform
several gradient-free global optimization algorithms in the design of PhCs
[95]. Bayesian optimization can also be extended to take advantage of gradient
information [96] and non-continuous constraints [97] to improve convergence.

## Acknowledgements

The authors thank Charlotte Loh for helpful discussions. This research is
based upon work supported in part by the National Defense Science &
Engineering Graduate (NDSEG) Fellowship Program, the Air Force Office of
Scientific Research under the award number FA9550-20-1-0115 and
FA9550-21-1-0299, the US Office of Naval Research (ONR) Multidisciplinary
University Research Initiative (MURI) grant N00014-20-1-2325 on Robust
Photonic Materials with High-Order Topological Protection and the U.S. Army
Research Office through the Institute for Soldier Nanotechnologies at MIT
under Collaborative Agreement Number W911NF-18-2-0048. MIT SuperCloud and
Lincoln Laboratory Supercomputing Center provided HPC resources that have
contributed to the results reported in this paper.

## Code availability

The code underlying the results presented here will be made publicly available
at https://github.com/samuelkim314/topo-phc-opt.

## References

  * Jensen and Sigmund [2010] J. Jensen and O. Sigmund, Topology optimization for nano-photonics, Laser Photonics Rev. 5, 308 (2010).
  * Molesky  _et al._ [2018] S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, Inverse design in nanophotonics, Nat. Photonics 12, 659 (2018).
  * Li  _et al._ [2019] W. Li, F. Meng, Y. Chen, Y. fan Li, and X. Huang, Topology optimization of photonic and phononic crystals and metamaterials: A review, Adv. Theory Simul. 2, 1900017 (2019).
  * Dobson and Cox [1999] D. C. Dobson and S. J. Cox, Maximizing band gaps in two-dimensional photonic crystals, SIAM J. Appl. Math. 59, 2108 (1999).
  * Cox and Dobson [2000] S. J. Cox and D. C. Dobson, Band structure optimization of two-dimensional photonic crystals in $H$-polarization, J. Comput. Phys. 158, 214 (2000).
  * Kao  _et al._ [2005] C. Y. Kao, S. Osher, and E. Yablonovitch, Maximizing band gaps in two-dimensional photonic crystals by using level set methods, Appl. Phys. B: Lasers Opt. 81, 235 (2005).
  * Men  _et al._ [2010] H. Men, N. C. Nguyen, R. M. Freund, P. A. Parrilo, and J. Peraire, Bandgap optimization of two-dimensional photonic crystals using semidefinite programming and subspace methods, J. Comput. Phys. 229, 3706 (2010).
  * Men  _et al._ [2014a] H. Men, K. Y. Lee, R. M. Freund, J. Peraire, and S. G. Johnson, Robust topology optimization of three-dimensional photonic-crystal band-gap structures, Opt. Express 22, 22632 (2014a).
  * Goh  _et al._ [2007] J. Goh, I. Fushman, D. Englund, and J. Vučković, Genetic optimization of photonic bandgap structures, Opt. Express 15, 8218 (2007).
  * Yan  _et al._ [2021] Y. Yan, P. Liu, X. Zhang, and Y. Luo, Photonic crystal topological design for polarized and polarization-independent band gaps by gradient-free topology optimization, Optics Express 29, 24861 (2021).
  * Lu  _et al._ [2014] L. Lu, J. D. Joannopoulos, and M. Soljačić, Topological photonics, Nat. Photonics 8, 821 (2014).
  * Ozawa  _et al._ [2019] T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, and I. Carusotto, Topological photonics, Rev. Mod. Phys. 91, 10.1103/revmodphys.91.015006 (2019).
  * Tang  _et al._ [2022] G.-J. Tang, X.-T. He, F.-L. Shi, J.-W. Liu, X.-D. Chen, and J.-W. Dong, Topological photonic crystals: Physics, designs, and applications, Laser Photonics Rev. 16, 2100300 (2022).
  * Haldane [1988] F. D. M. Haldane, Model for a quantum Hall effect without Landau levels: Condensed-matter realization of the “parity anomaly”, Phys. Rev. Lett. 61, 2015 (1988).
  * Raghu and Haldane [2008] S. Raghu and F. D. M. Haldane, Analogs of quantum-Hall-effect edge states in photonic crystals, Phys. Rev. A 78, 033834 (2008).
  * Wang  _et al._ [2009] Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, Observation of unidirectional backscattering-immune topological electromagnetic states, Nature 461, 772 (2009).
  * Lu  _et al._ [2015] L. Lu, Z. Wang, D. Ye, L. Ran, L. Fu, J. D. Joannopoulos, and M. Soljačić, Experimental observation of Weyl points, Science 349, 622 (2015).
  * Xu  _et al._ [2015] S.-Y. Xu, I. Belopolski, N. Alidoust, M. Neupane, G. Bian, C. Zhang, R. Sankar, G. Chang, Z. Yuan, C.-C. Lee, S.-M. Huang, H. Zheng, J. Ma, D. S. Sanchez, B. Wang, A. Bansil, F. Chou, P. P. Shibayev, H. Lin, S. Jia, and M. Z. Hasan, Discovery of a Weyl fermion semimetal and topological Fermi arcs, Science 349, 613 (2015).
  * Lv  _et al._ [2015] B. Q. Lv, H. M. Weng, B. B. Fu, X. P. Wang, H. Miao, J. Ma, P. Richard, X. C. Huang, L. X. Zhao, G. F. Chen, Z. Fang, X. Dai, T. Qian, and H. Ding, Experimental discovery of Weyl semimetal TaAs, Phys. Rev. X 5, 031013 (2015).
  * Jörg  _et al._ [2022] C. Jörg, S. Vaidya, J. Noh, A. Cerjan, S. Augustine, G. von Freymann, and M. C. Rechtsman, Observation of quadratic (charge-2) Weyl point splitting in near-infrared photonic crystals, Laser Photonics Rev. 16, 2100452 (2022).
  * Shi and Yang [2020] X. Shi and J. Yang, Spin-1 Weyl point and surface arc state in a chiral phononic crystal, Phys. Rev. B 101, 10.1103/physrevb.101.214309 (2020).
  * Zhang  _et al._ [2018] T. Zhang, Z. Song, A. Alexandradinata, H. Weng, C. Fang, L. Lu, and Z. Fang, Double-Weyl phonons in transition-metal monosilicides, Phys. Rev. Lett. 120, 016401 (2018).
  * Wang  _et al._ [2016] H. Wang, L. Xu, H. Chen, and J.-H. Jiang, Three-dimensional photonic Dirac points stabilized by point group symmetry, Phys. Rev. B 93, 235155 (2016).
  * Guo  _et al._ [2019] Q. Guo, O. You, B. Yang, J. B. Sellman, E. Blythe, H. Liu, Y. Xiang, J. Li, D. Fan, J. Chen, _et al._ , Observation of three-dimensional photonic Dirac points and spin-polarized surface arcs, Phys. Rev. Lett. 122, 203903 (2019).
  * Wang  _et al._ [2017] H.-X. Wang, Y. Chen, Z. H. Hang, H.-Y. Kee, and J.-H. Jiang, Type-II Dirac photons, npj Quant. Mater. 2, 54 (2017).
  * Lu  _et al._ [2013] L. Lu, L. Fu, J. D. Joannopoulos, and M. Soljačić, Weyl points and line nodes in gyroid photonic crystals, Nat. Photonics 7, 294 (2013).
  * Yan  _et al._ [2018] Q. Yan, R. Liu, Z. Yan, B. Liu, H. Chen, Z. Wang, and L. Lu, Experimental discovery of nodal chains, Nat. Phys. 14, 461 (2018).
  * Gao  _et al._ [2018] W. Gao, B. Yang, B. Tremain, H. Liu, Q. Guo, L. Xia, A. P. Hibbins, and S. Zhang, Experimental observation of photonic nodal line degeneracies in metacrystals, Nat. Commun. 9, 1 (2018).
  * Xia  _et al._ [2019] L. Xia, Q. Guo, B. Yang, J. Han, C.-X. Liu, W. Zhang, and S. Zhang, Observation of hourglass nodal lines in photonics, Phys. Rev. Lett. 122, 103903 (2019).
  * Yang  _et al._ [2020] E. Yang, B. Yang, O. You, H.-C. Chan, P. Mao, Q. Guo, S. Ma, L. Xia, D. Fan, Y. Xiang, _et al._ , Observation of non-Abelian nodal links in photonics, Phys. Rev. Lett. 125, 033901 (2020).
  * Park  _et al._ [2021] H. Park, S. Wong, X. Zhang, and S. S. Oh, Non-Abelian charged nodal links in a dielectric photonic crystal, ACS Photonics 8, 2746 (2021).
  * Wang  _et al._ [2021] D. Wang, B. Yang, Q. Guo, R.-Y. Zhang, L. Xia, X. Su, W.-J. Chen, J. Han, S. Zhang, and C. T. Chan, Intrinsic in-plane nodal chain and generalized quaternion charge protected nodal link in photonics, Light: Sci. Appl. 10, 1 (2021).
  * Park  _et al._ [2022] H. Park, W. Gao, X. Zhang, and S. S. Oh, Nodal lines in momentum space: topological invariants and recent realizations in photonic and other systems, Nanophotonics 11, 2779 (2022).
  * Xie  _et al._ [2021] X. Xie, J. Dang, S. Yan, W. Zhang, H. Hao, S. Xiao, S. Shi, Z. Zuo, H. Ni, Z. Niu, _et al._ , Optimization and robustness of the topological corner state in second-order topological photonic crystals, Opt. Express 29, 30735 (2021).
  * Christiansen  _et al._ [2019] R. E. Christiansen, F. Wang, O. Sigmund, and S. Stobbe, Designing photonic topological insulators with quantum-spin-Hall edge states using topology optimization, Nanophotonics 8, 1363 (2019).
  * Lin  _et al._ [2018] Z. Lin, L. Christakis, Y. Li, E. Mazur, A. W. Rodriguez, and M. Lončar, Topology-optimized dual-polarization Dirac cones, Phys. Rev. B 97, 081408 (2018).
  * Nussbaum  _et al._ [2022] E. Nussbaum, N. Rotenberg, and S. Hughes, Optimizing the chiral Purcell factor for unidirectional single-photon emitters in topological photonic crystal waveguides using inverse design, Phys. Rev. A 106, 033514 (2022).
  * Alexandradinata [2015] A. Alexandradinata,  _Spatially-protected Topology and Group Cohomology in Band Insulators_ , Ph.D. thesis, Princeton University (2015).
  * Bradlyn  _et al._ [2017] B. Bradlyn, L. Elcoro, J. Cano, M. G. Vergniory, Z. Wang, C. Felser, M. I. Aroyo, and B. A. Bernevig, Topological quantum chemistry, Nature 547, 298 (2017).
  * Po  _et al._ [2017] H. C. Po, A. Vishwanath, and H. Watanabe, Symmetry-based indicators of band topology in the 230 space groups, Nat. Commun. 8, 50 (2017).
  * Kruthoff  _et al._ [2017] J. Kruthoff, J. de Boer, J. van Wezel, C. L. Kane, and R.-J. Slager, Topological classification of crystalline insulators through band structure combinatorics, Phys. Rev. X 7, 041069 (2017).
  * de Paz  _et al._ [2019] M. B. de Paz, M. G. Vergniory, D. Bercioux, A. García-Etxarri, and B. Bradlyn, Engineering fragile topology in photonic crystals: Topological quantum chemistry of light, Phys. Rev. Research 1, 032005 (2019).
  * Christensen  _et al._ [2022] T. Christensen, H. C. Po, J. D. Joannopoulos, and M. Soljačić, Location and topology of the fundamental gap in photonic crystals, Phys. Rev. X 12, 021066 (2022).
  * Joannopoulos  _et al._ [2008] J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade,  _Photonic Crystals: Molding the Flow of Light_ , 2nd ed. (Princeton University Press, 2008).
  * Ho  _et al._ [1990] K. Ho, C. T. Chan, and C. M. Soukoulis, Existence of a photonic gap in periodic dielectric structures, Phys. Rev. Lett. 65, 3152 (1990).
  * Miller  _et al._ [2014] O. D. Miller, C. W. Hsu, M. T. H. Reid, W. Qiu, B. G. DeLacy, J. D. Joannopoulos, M. Soljačić, and S. G. Johnson, Fundamental limits to extinction by metallic nanoparticles, Phys. Rev. Lett. 112, 123903 (2014).
  * Aroyo [2016] M. I. Aroyo, ed.,  _International Tables for Crystallography Volume A: Space-group Symmetry_, 6th ed. (Wiley, 2016).
  * Johnson [2014] S. G. Johnson, The NLopt nonlinear-optimization package (2014).
  * Gablonsky and Kelley [2001] J. M. Gablonsky and C. T. Kelley, A locally-biased form of the DIRECT algorithm, J. Glob. Optim. 21, 27 (2001).
  * Rowan [1990] T. H. Rowan,  _Functional stability analysis of numerical algorithms_ , Ph.D. thesis, The University of Texas at Austin (1990).
  * Runarsson and Yao [2005] T. P. Runarsson and X. Yao, Search biases in constrained evolutionary optimization, IEEE Trans. Syst. Man Cybern. 35, 233 (2005).
  * Johnson and Joannopoulos [2001] S. G. Johnson and J. D. Joannopoulos, Block-iterative frequency-domain methods for maxwell’s equations in a planewave basis, Opt. Express 8, 173 (2001).
  * Po [2020] H. C. Po, Symmetry indicators of band topology, J. Phys.: Condens. Matter 32, 263001 (2020).
  * [54] MPBUtils.jl (v0.1.11), https://github.com/thchr/MPBUtils.jl.
  * Hinuma  _et al._ [2017] Y. Hinuma, G. Pizzi, Y. Kumagai, F. Oba, and I. Tanaka, Band structure diagram paths based on crystallography, Comput. Mater. Sci. 128, 140 (2017).
  * [56] Brillouin.jl (v0.5.10), https://github.com/thchr/Brillouin.jl.
  * Cano and Bradlyn [2021] J. Cano and B. Bradlyn, Band representations and topological quantum chemistry, Annu. Rev. Condens. Matter Phys. 12, 225 (2021).
  * Bouckaert  _et al._ [1936] L. P. Bouckaert, R. Smoluchowski, and E. Wigner, Theory of Brillouin zones and symmetry properties of wave functions in crystals, Phys. Rev. 50, 58 (1936).
  * Kaxirax and Joannopoulos [2019] E. Kaxirax and J. D. Joannopoulos,  _Quantum theory of materials_ (Cambridge University Press, 2019).
  * Song  _et al._ [2018] Z. Song, T. Zhang, and C. Fang, Diagnosis for nonmagnetic topological semimetals in the absence of spin-orbital coupling, Phys. Rev. X 8, 031069 (2018).
  * Elcoro  _et al._ [2020] L. Elcoro, Z. Song, and B. A. Bernevig, Application of induction procedure and Smith decomposition in calculation and topological classification of electronic band structures in the 230 space groups, Phys. Rev. B 102, 035110 (2020).
  * Po  _et al._ [2016] H. C. Po, H. Watanabe, M. P. Zaletel, and A. Vishwanath, Filling-enforced quantum band insulators in spin-orbit coupled crystals, Science Adv. 2, e1501782 (2016).
  * Mañes [2012] J. L. Mañes, Existence of bulk chiral fermions and crystal symmetry, Phys. Rev. B 85, 155118 (2012).
  * Yu  _et al._ [2022] Z.-M. Yu, Z. Zhang, G.-B. Liu, W. Wu, X.-P. Li, R.-W. Zhang, S. A. Yang, and Y. Yao, Encyclopedia of emergent particles in three-dimensional crystals, Sci. Bull. 67, 375 (2022).
  * Watanabe and Lu [2018] H. Watanabe and L. Lu, Space group theory of photonic bands, Phys. Rev. Lett. 121, 263903 (2018).
  * Zhang  _et al._ [2020a] T. Zhang, L. Lu, S. Murakami, Z. Fang, H. Weng, and C. Fang, Diagnosis scheme for topological degeneracies crossing high-symmetry lines, Phys. Rev. Research 2, 022066 (2020a).
  * Vanderbilt [2018] D. Vanderbilt,  _Berry Phases in Electronic Structure Theory_ (Cambridge University Press, 2018).
  * Liu  _et al._ [2018] B. Liu, J. D. Joannopoulos, S. G. Johnson, and L. Lu, Generalized Gilat–Raubenheimer method for density-of-states calculation in photonic crystals, J. Opt. 20, 044005 (2018).
  * Harrison  _et al._ [2007] J. M. Harrison, P. Kuchment, A. Sobolev, and B. Winn, On occurrence of spectral edges for periodic operators inside the Brillouin zone, J. Phys. A: Math. Theor. 40, 7597 (2007).
  * Craster  _et al._ [2012] R. V. Craster, T. Antonakakis, M. Makwana, and S. Guenneau, Dangers of using the edges of the Brillouin zone, Phys. Rev. B 86, 115130 (2012).
  * Maurin  _et al._ [2018] F. Maurin, C. Claeys, E. Deckers, and W. Desmet, Probability that a band-gap extremum is located on the irreducible Brillouin-zone contour for the 17 different plane crystallographic lattices, Int. J. Solids Struct. 135, 26 (2018).
  * Wang  _et al._ [2008] Z. Wang, Y. Chong, J. D. Joannopoulos, and M. Soljačić, Reflection-free one-way edge modes in a gyromagnetic photonic crystal, Phys. Rev. Lett. 100, 013905 (2008).
  * Ao  _et al._ [2009] X. Ao, Z. Lin, and C. T. Chan, One-way edge mode in a magneto-optical honeycomb photonic crystal, Phys. Rev. B 80, 033105 (2009).
  * Skirlo  _et al._ [2014] S. A. Skirlo, L. Lu, and M. Soljačić, Multimode one-way waveguides of large Chern numbers, Phys. Rev. Lett. 113, 113904 (2014).
  * Liu  _et al._ [2012] K. Liu, L. Shen, and S. He, One-way edge mode in a gyromagnetic photonic crystal slab, Opt. Lett. 37, 4110 (2012).
  * Skirlo  _et al._ [2015] S. A. Skirlo, L. Lu, Y. Igarashi, Q. Yan, J. Joannopoulos, and M. Soljačić, Experimental observation of large Chern numbers in photonic crystals, Phys. Rev. Lett. 115, 253901 (2015).
  * Poo  _et al._ [2011a] Y. Poo, R.-x. Wu, Z. Lin, Y. Yang, and C. T. Chan, Experimental realization of self-guiding unidirectional electromagnetic edge states, Phys. Rev. Lett. 106, 093903 (2011a).
  * Yang  _et al._ [2013] Y. Yang, Y. Poo, R. xin Wu, Y. Gu, and P. Chen, Experimental demonstration of one-way slow wave in waveguide involving gyromagnetic photonic crystals, Appl. Phys. Lett. 102, 231113 (2013).
  * Poo  _et al._ [2011b] Y. Poo, R.-x. Wu, Z. Lin, Y. Yang, and C. T. Chan, Experimental realization of self-guiding unidirectional electromagnetic edge states, Phys. Rev. Lett. 106, 093903 (2011b).
  * Lu  _et al._ [2018] L. Lu, H. Gao, and Z. Wang, Topological one-way fiber of second Chern number, Nat. Commun. 9, 1 (2018).
  * Devescovi  _et al._ [2021] C. Devescovi, M. García-Díez, I. Robredo, M. Blanco de Paz, J. Lasa-Alonso, B. Bradlyn, J. L. Mañes, M. G Vergniory, and A. García-Etxarri, Cubic 3D Chern photonic insulators with orientable large Chern vectors, Nat. Commun. 12, 1 (2021).
  * Fukui  _et al._ [2005] T. Fukui, Y. Hatsugai, and H. Suzuki, Chern numbers in discretized Brillouin zone: Efficient method of computing (spin) Hall conductances, J. Phys. Soc. Jpn. 74, 1674 (2005).
  * Fang  _et al._ [2012] C. Fang, M. J. Gilbert, and B. A. Bernevig, Bulk topological invariants in noninteracting point group symmetric insulators, Phys. Rev. B 86, 115112 (2012).
  * Note [1] This mapping involves a fixed convention for the choice of nontrivial generators in Eq. 7 [60]. Specifically, we take $\mathbf{b}^{(75)}_{1}=-\hat{\mathbf{e}}_{\Gamma_{2}}+\hat{\mathbf{e}}_{\Gamma_{4}}-\hat{\mathbf{e}}_{Z_{2}}+\hat{\mathbf{e}}_{Z_{4}}$ and $\mathbf{b}^{(168)}_{1}=+\hat{\mathbf{e}}_{H_{2}}-\hat{\mathbf{e}}_{H_{3}}+\hat{\mathbf{e}}_{K_{2}}-\hat{\mathbf{e}}_{K_{3}}-\hat{\mathbf{e}}_{L_{1}}+\hat{\mathbf{e}}_{L_{2}}-\hat{\mathbf{e}}_{M_{1}}+\hat{\mathbf{e}}_{M_{2}}$. Equivalently, the symmetry indicators are given by $\nu_{1}^{(75)}=2n_{M_{2}}+n_{M_{3}}-n_{M_{4}}-2n_{\Gamma_{1}}-n_{\Gamma_{3}}+n_{\Gamma_{4}}-2n_{X_{1}}\bmod 4$ and $\nu_{1}^{(168)}=n_{\Gamma_{1}}-2n_{\Gamma_{2}}+5n_{\Gamma_{3}}+2n_{\Gamma_{4}}+3n_{\Gamma_{5}}+2n_{H_{1}}-2n_{H_{2}}-3n_{L_{1}}\bmod 6$, with irrep labels given in the CDML notation [98].
  * Zhang  _et al._ [2020b] T. Zhang, L. Lu, S. Murakami, Z. Fang, H. Weng, and C. Fang, Diagnosis scheme for topological degeneracies crossing high-symmetry lines, Phys. Rev. Research 2, 022066 (2020b).
  * Benalcazar  _et al._ [2017] W. A. Benalcazar, B. A. Bernevig, and T. L. Hughes, Electric multipole moments, topological multipole moment pumping, and chiral hinge states in crystalline insulators, Phys. Rev. B 96, 245115 (2017).
  * Ota  _et al._ [2019] Y. Ota, F. Liu, R. Katsumi, K. Watanabe, K. Wakabayashi, Y. Arakawa, and S. Iwamoto, Photonic crystal nanocavity based on a topological corner state, Optica 6, 786 (2019).
  * Sigmund [2007] O. Sigmund, Morphology-based black and white filters for topology optimization, Struct. Multidiscip. Optim. 33, 401 (2007).
  * Lazarov  _et al._ [2016] B. S. Lazarov, F. Wang, and O. Sigmund, Length scale and manufacturability in density-based topology optimization, Arch. Appl. Mech. 86, 189 (2016).
  * Men  _et al._ [2014b] H. Men, R. M. Freund, N. C. Nguyen, J. Saa-Seoane, and J. Peraire, Fabrication-adaptive optimization with an application to photonic crystal design, Oper. Res. 62, 418 (2014b).
  * Hammond  _et al._ [2022] A. M. Hammond, A. Oskooi, M. Chen, Z. Lin, S. G. Johnson, and S. E. Ralph, High-performance hybrid time/frequency-domain topology optimization for large-scale photonics inverse design, Opt. Express 30, 4467 (2022).
  * Li  _et al._ [2016] Q. Li, W. Chen, S. Liu, and L. Tong, Structural topology optimization considering connectivity constraint, Structural and Multidisciplinary Optimization 54, 971 (2016).
  * Oran  _et al._ [2018] D. Oran, S. G. Rodriques, R. Gao, S. Asano, M. A. Skylar-Scott, F. Chen, P. W. Tillberg, A. H. Marblestone, and E. S. Boyden, 3d nanofabrication by volumetric deposition and controlled shrinkage of patterned scaffolds, Science 362, 1281 (2018).
  * Mills  _et al._ [2021] B. Mills, Y. Salamin, G. Yang, D. Oran, Y. Sun, S. Maayani, S. E. Kooi, A. A. E. Elawad, J. J. Lopez, C. Swain, _et al._ , Implosion fabrication as a platform for three-dimensional nanophotonics, in  _CLEO: Science and Innovations_ (Optical Society of America, 2021) pp. STh4J–1.
  * Kim  _et al._ [2022] S. Kim, P. Y. Lu, C. Loh, J. Smith, J. Snoek, and M. Soljacic, Deep learning for Bayesian optimization of scientific problems with high-dimensional structure, Transactions on Machine Learning Research (2022).
  * Ament and Gomes [2022] S. E. Ament and C. P. Gomes, Scalable first-order Bayesian optimization via structured automatic differentiation, in  _International Conference on Machine Learning_ (PMLR, 2022) pp. 500–516.
  * Gelbart [2015] M. A. Gelbart,  _Constrained Bayesian optimization and applications_ , Ph.D. thesis, Harvard University (2015).
  * Cracknell  _et al._ [1979] A. P. Cracknell, B. L. Davies, S. C. Miller, and W. F. Love,  _Kronecker Product Tables. General Introduction and Tables of Irreducible Representations of Space Groups._ , Vol. 1 (New York: IFI/Plenum, 1979).