Constraining Continuous Topology Optimizations to Discrete Solutions for Photonic Applications

Photonic topology optimization is a technique used to find the permittivity distribution of a device that optimizes an electromagnetic figure-of-merit. Two common versions are used: continuous density-based optimizations that optimize a gray scale permittivity defined over a grid, and discrete level-set optimizations that optimize the shape of the material boundary of a device. In this work we present a method for constraining a continuous optimization such that it is guaranteed to converge to a discrete solution. This is done by inserting a constrained suboptimization with low computational overhead cost at each iteration of an overall gradient-based optimization. The technique adds only one hyperparameter with straightforward behavior to control the aggressiveness of binarization. Computational examples are provided to analyze the hyperparameter behavior, show this technique can be used in conjunction with projection filters, show the benefits of using this technique to provide a nearly discrete starting point for subsequent level-set optimization, and show that an additional hyperparameter can be introduced to control the overall material/void fraction. This method excels for problems where the electromagnetic figure-of-merit is majorly affected by the binarization requirement and situations where identifying suitable hyperparameter values becomes challenging with existing methods.


Derived solution to Lagrange dual problem from main text
To obtain the dual problem, we start with the Lagrangian. This is defined as: The dual function is defined as the infimum of the Lagrangian over ⃗ x: Since (S1) is linear, the infimum is −∞ if any slope is non-zero.
The optimization we wish to solve is called the Lagrange dual problem: Incorporating (S3) into Eq. S4, we get: The fourth constraint of (S5) can be used to write the FoM and the second constraint in terms of only ν and λ 1 .
For any ν, the optimal elements of ⃗ λ 1 are as small as possible since ⃗ u ≥ 0 and ⃗ l ≤ 0. The solution to ⃗ λ 1 is obtained by satisfying the constraints while minimizing ⃗ λ 1 . The solution for the i-th element of ⃗ λ 1 is thus: where H(x) is the Heaviside step function. ⃗ λ 2 is solved using the fourth constraint of (S5).
Here, H(⃗ x) represents the Heaviside function applied element-wise to ⃗ x and ⊙ represents element-wise multiplication. Since the optimal ⃗ λ 1 and ⃗ λ 2 are dependent on ν, the dual problem is a function of only ν ∈ R. The optimization can be reduced to a single dimension as follows.
Note that we have converted this from a maximization problem to a minimization problem, and omitted a −c · l term since it is constant. This problem can be solved numerically using open-source functions. We denote the optimal dual variables as (ν ⋆ , ⃗ λ ⋆ 1 , ⃗ λ ⋆ 2 ). The optimal solutions of the dual problem can be mapped back to original variable x using the Karush-Kuhn-Tucker (KKT) conditions, which are a set of conditions that hold under certain regularity conditions. Fortunately, linearity of the constraints is a sufficient condition for the KKT conditions to hold. We use the KKT condition of complementary slackness The solution is

Computational complexity
The Matlab code in Section 4 can be used to evaluate the computational complexity of dual-simplex, interiorpoint, and the method used in this manuscript. The analysis was done on a computer with the following specifications: • PC: AMD Ryzen Threadripper 3990X 64-Core Processor 2.90 GHz • RAM: 128 GB (4x32 GB, 3200 MHz DDR4) • GPU: GeForce RTX 3070 • OS: Windows 10 (21H2) • MATLAB version: R2020B Update 7 Figure S1 shows the computation time for each algorithm as a function of the input dimension size. Similar plots were obtained with Python's SciPy library, and code for this is available on request. Figure S1: Computational complexity of various algorithms that can be used to solve Eq. 3 in the main manuscript. The red, green, and blue lines plot the computation time with respect to the number of dimensions when using dual-simplex, interior-point, and the method described in Section 1. The grey dashed lines are representative of N , N 2 , and N 3 scaling.

Convergence tests
Optimizations performed in the main manuscript used a λ 0 /(10 × n max ) FDTD grid resolution, where λ 0 is the center wavelength of the optimized band and n max is the maximum refractive index that the device voxels can assume (which varies from 1.0 to 3.0 in the main manuscript). To confirm that the final results are realistic, we simulate the final structure on a finer grid of λ 0 /(20 × n max ). The transmissions through the three different apertures of the focal plane are drawn in Fig. S2, similar to how the transmissions were plotted in Fig. 1 of the main manuscript.