Open Access
Add to BibSonomy
Abstract
The complexity of applications addressed with photonic integrated circuits is steadily rising and poses increasingly challenging demands on individual component functionality, performance and footprint. Inverse design methods have recently shown great promise to address these demands using fully automated design procedures that enable access to non-intuitive device layouts beyond conventional nanophotonic design concepts. Here we present a dynamic binarization method for the objective-first algorithm that lies at the core of the currently most successful inverse design algorithms. Our results demonstrate significant performance advantages over previous implementations of objective first algorithms, which we show for a fundamental TE00 to TE20 waveguide mode converter both in simulation and in experiments with fabricated devices.
© 2023 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
The design process of photonic integrated circuit components typically relies on brute force optimization of critical parameters in an initial layout that was based on some physical intuition deriving from established concepts such as interference, diffraction, coupled mode theory or similar. While a wide range of devices like directional couplers [1], resonators [2], photonic crystals [3], adiabatic tapers [4] or coupling structures from optical fibers to on-chip integrated photonic waveguides have been created following this paradigm, they all suffer from large device footprints and a limited number of tunable design parameters. The few degrees of freedom available to the designer (typically less than ten) usually only provide access to a tiny fraction of the solution space, which may not contain layouts combining optimal efficiency and small footprint. Nanophotonic circuits offer solutions to problems of ever-increasing complexity [5]. In order to guarantee future scalability, fully automated design processes for more efficient and compact devices that exhibit novel functionalities based on complex interference phenomena not accessible by human intuition are needed.
Recently, various techniques to design nanophotonic devices following a top-down approach, such as topology optimization [6–8], genetic algorithms [9,10], particle swarm optimization [11] or the objective-first approach [12], have emerged. The objective-first algorithm has proven particularly successful for designing high-efficiency and low-footprint nanophotonic devices [12–15] and currently constitutes one of the most promising nanophotonic inverse design methods.
A problem of the objective-first approach however is that the resulting permittivity distribution consists of continuously varying values, while modern nanophotonic fabrication methods only allow for permittivity distributions containing two (few) discrete values, as each position in the design can only be occupied by either one material, e.g. the dielectric, or another, e.g. air or silicon dioxide. In the originally proposed version of the objective first algorithm, reconciliation with a discrete permittivity distribution is accomplished by thresholding the continuous permittivity map. This constitutes a dramatic intervention in the optimization process because thresholding is not being guided by any gradient information. Even subsequently applied adjoint sensitivity analysis [6,13] is inadequate to compensate for the loss of the objective-first-specific "global" optimization characteristic and will only yield solutions that are far from any optimum the objective-first algorithm had initially converged to.
Static binarization penalty terms to the objective-first method have been proposed to overcome this problem by steadily increasing the binarity (dielectric vs. air) of the structure with each iteration [16]. While this modification has shown improvements, it tends to limit the ability of the algorithm to develop freely during optimization. The corresponding penalty term uses a thresholded version of the device structure in a previous iteration as a reference for calculating the update of the binarization level, which is why we refer to this technique as the "static" binarization method. We here remove the aforementioned structural bias by calculating the binarization level update in real time in each iteration of the convex optimization step. To do so, we exploit so called mixed-integer optimization capabilities of state-of-the-art convex optimization software. Our "dynamic" binarization method thus removes critical limitations of "static" techniques on the iterative evolution of the device structure towards an efficient solution.
We here demonstrate the superior performance of our approach for the example of a spatial waveguide mode converter, transforming fundamental transverse electric modes, TE$_{00}$, into TE$_{20}$ modes, which is challenging to achieve with other design methods. Such mode converters are fundamental components as they are essential for various tasks such as optical fiber coupling [17,18], phase matching [19] or coupling between different waveguide geometries [12]. However, our approach straightforwardly applies to a wide range of other nanophotonic circuit components.
We obtain a highly efficient and compact device optimized in full 3D exhibiting a conversion efficiency of 81.1 % and a footprint of just 10 µm$^2$. In direct comparison to a device optimized using the static binarization method we observe an increase of $9.8$ % in conversion efficiency. We fabricate the device on a 200 nm thin Si$_3$N$_4$-on-insulator platform using electron beam lithography and verify its functionality using mode imaging techniques.
2. Inverse design algorithm
While manual design methods usually follow a bottom-up approach relying predominantly on phenomena accessible by human intuition and brute force parameter optimization, inverse design algorithms follow a top-down approach by specifying given input and desired output waveforms. The design region is considered as a black box, spatially discretized into pixels. The algorithm then tries to find a permittivity distribution, $\epsilon$, within the design region that fulfills the design objectives. The structure of our algorithm is shown in Fig. 1. In general, the algorithm is split into two phases, where the first phase follows a global optimization paradigm under allowance of continuous permittivity values. By splitting the global nanophotonic black-box optimization problem where both, the field distribution and the structure are unknown, into two sub-problems, namely the composition of the target field (field sub-problem) and solving for the corresponding structure (structural sub-problem), we can apply the alternating direction method of multipliers (ADMM) [20]. Following the "objective-first" approach, we only allow for small changes of the target field within this process. We derive the two sub-problems starting from a transformation of the time-harmonic Maxwell equations into a linear algebra form shown in Eq. (1) and (2).
where $\mu _0$ is the permeability, $z$ is representative of the permittivity $\epsilon$, $x$ the electric field $E$, and $b$ depends on the electric current density $J$. The field sub-problem can be expressed asTo obtain Eq. (4) we utilize the bi-linearity of Eq. (2), and find that $A(z)x - b(z) = B(x)z - d(x)$. It is this relation, that qualifies the problem to be solved using ADMM. Here, the structure parameter $z$ is constrained to fulfill the structural objectives, namely to reside between $z_{\mathrm {min}}=0$ and $z_{\mathrm {max}}=1$. $z$ is linearly mapped to the material-dependent minimum and maximum permittivity $\epsilon _{\mathrm {min}}$ and $\epsilon _{\mathrm {max}}$, respectively. Enforcing the specified design objectives $f(x)$ in Eq. (3) leads to Eq. (2) not being fulfilled. The residual of this function ("physics residual") is the figure of merit of the optimization procedure. This method can be extended to support multiple objectives by transforming Eq. (2) into a sum over all desired objectives involved.
After phase 1 has converged, we leave the continuous variable regime and transition to phase 2, where we convert the structure to a boundary parameterization using a level-set representation and fine tune the device following a first order gradient in the direction of steepest descent. The gradient field with respect to the objective function is calculated by performing an adjoint sensitivity analysis involving one forward and one adjoint simulation per iteration [21,22].
Fig. 1. Our optimization algorithm can be split into two separate stages. In the first step we optimize the structure with continuous permittivity values using an Alternating Direction Method of Multipliers-implementation of the "objective-first" algorithm. We alternate between solving the structural sub-problem and solving the field sub-problem until a termination criterion is met, which is either tied to the physics residual or the maximum number of iterations. Afterwards, we convert the resulting permittivity distribution into a level-set representation to fine tune the structure using the adjoint sensitivity analysis and following the gradient’s steepest descent. The second stage is terminated after a fixed number of iterations.
3. Binarization of the continuous structure
The objective-first optimization strategy suffers from a significant drawback. In its original form [12] the resulting structure features continuous permittivity values ranging from user-specified minimum to maximum permittivities. As typically only distributions containing two discrete permittivity values can be realized in most modern nanofabrication processes, such as electron beam lithography, the resulting structure has to be converted from continuous to binary values.
The initial solution [12] was to threshold the continuous distribution, i.e. to introduce a critical permittivity value between $\epsilon _{\mathrm {min}}$ and $\epsilon _{\mathrm {max}}$. Each $\epsilon$ that lies below the critical value is being set to $\epsilon _{\mathrm {min}}$ and each value that lies above it is being set to $\epsilon _{\mathrm {max}}$. As this is a critical intervention in the optimization flow, which is not based on any gradient information, the resulting structure can be expected to differ significantly from the result of the first optimization step in both permittivity distribution and performance.
Efforts to overcome this issue have relied on introducing a static binarization penalty term to Eq. (4), such that the minimization task can be expressed as [16,23]
To remove this limitation, we calculate $z_\mathrm {bin}$ dynamically while solving the structural sub-problem. We configure the variable to be constrained to binary values and make it subject to change while solving the problem with a convex optimization toolkit. We therefore refer to this method as dynamic bynarization. As the problem now involves a non-continuous variable, it has become a mixed-integer-optimization problem, which requires specific solver capabilities and care to be taken while formulating the minimization problem. The hereby obtained flexibility in $z_{\mathrm {bin}}$ enables stricter formulations of the binarization function such as shown in Eq. (6). This problem statement involving a hard constraint enforces a certain level of binarity in each iteration.
Here $\lambda$ is being decreased during the optimization process with $\lambda = 1$ allowing an unrestricted continuous structure and $\lambda =0$ enforcing completely binary structure variables. To technically implement the hard constraint we introduced the normalization term $N=||z_\mathrm {even}||$ with $z_\mathrm {even}$ representing evenly distributed structure variables $z_i=0.5$.
4. Simulation procedure and computational methods
To calculate the electromagnetic target field within the device that we aim to find a correspondent permittivity distribution for, we need to perform two electromagnetic simulations for each input mode involved. One forward simulation using the specified input mode as a source and an adjoint simulation where the desired output mode provides the source. We simulate these fields using the finite-difference frequency-domain (FDFD) framework OpenCLFDFD [24], which enables hardware accelerated parallel computation. We run the FDFD computations on High Performance Computing resources (Nvidia Tesla A100 GPUs on the PALMA II cluster) to reduce the computation time for each simulation.
To solve the minimization problems stated in Eqs. (5) and (6), we use the Matlab package CVX [25]. To perform convex optimization involving binary variables we employ the mixed-integer optimization capabilities of the MOSEK solver [26] supported by CVX. We configure the solver to use sophisticated heuristics [27] to ensure that a good initial solution to the mixed integer problem can be found quickly.
5. Design of a TE$_{00}$ to TE$_{20}$ waveguide mode converter
To demonstrate the benefits of using dynamic over static binarization functions, we design a TE$_{00}$ to TE$_{20}$ waveguide mode converter for a wavelength of $775$ nm in full 3D. We chose 200 nm thin silicon nitride (Si$_3$N$_4$) on 2 µm thick silicon dioxide (SiO$_2$) as a material platform. The pixel size is set to 40 nm and the initial structure is composed of uniformly distributed permittivity values, i.e. $z_i = 0.5$. We specify the device width and length as 3.64 µm and 2.84 µm, respectively. The objective function $f(x)$ is configured such that the algorithm will maximize the power present in the desired TE$_{20}$ mode and suppress other transverse electric modes:
Here, $\eta _m(x)$ is the overlap integral of TE-mode number $m$ with the electric field $x$ at the output waveguide of the device. For a negligible reflection from the output port, which we can assume due to the usage of surrounding PMLs, we calculate the mode overlap as follows:
with $P_m$ and $P_{src}$ being the power in mode $m$ and the power of the source mode, respectively. $H_m$ is the magnetic field of the target mode $m$.We conducted 200 iterations of the objective-first step followed by 200 iterations using the steepest-descent method. The resulting performance of devices simulated simulated with static and dynamic binarization are shown in Fig. 2. The data shows that the conversion efficiency into the TE$_{20}$ mode of the dynamically binarized device is higher over the entire simulated range from $475$ nm to $1025$ nm wavelength. While the device was only optimized for the design wavelength of $775$ nm, its $3$ dB bandwidth ranges from $575$ nm to $925$ nm. We further observe that the dynamic binarization method is capable of reducing the physics residual more efficiently than static binarization approaches. Figure 3 compares the evolution of the physics residual for a soft constraint implementation (following Eq. (5)) of a dynamic and a static binarization term. We find that, especially in later stages of the optimization, the dynamic binarization term is able to lower the physics residual further than the static version while even producing a slightly higher level of binarity. Note that small differences in the physics residual can have a significant impact in the final device performance because the result of this optimization step is used as an input to fine tune the structure using a steepest descent method.
Fig. 2. Conversion efficiencies of the fundamental TE$_{00}$ input-to lowest-order output modes, simulated for a range of wavelengths centered around the design wavelength of $775$ nm. Results for the device optimized using the dynamic and static binarization methods are depicted in dots and crosses, respectively.
Fig. 3. Evolution of the physics residual during the optimization for both the dynamic and static binarization methods, implemented as a soft constraint with the corresponding binarity of the structure. The oscillations in the residuals are related to structure resets because the algorithm did not succeed in lowering the residual in a certain direction of the optimization.
6. Fabrication and experimental verification
To verify the functionality of the device, we fabricate the calculated binary permittivity distributions on a $200$ nm silicon nitride on $2$ µm silicon dioxide on silicon platform using high-resolution electron beam lithography (EBPG 5150, Raith) and reactive ion etching. The left side of Fig. 4 (a) shows a scanning electron micro-graph of the fabricated device, which has been optimized using the dynamic binarization method. The output waveguide is expanded adiabatically until it reaches a width of $10$ µm to enlarge the mode profile while maintaining the mode composition. The enlarged profile emerging from the cleaved waveguide facet, as depicted on the right side of Fig. 4 (a), can then be examined using a microscope objective. We illuminate the device with a Ti:Sapphire laser operating at $\lambda =775$ nm wavelength, connected to a polarization-controlled single-mode optical fiber array, which interfaces with the nanophpotonic chip via grating couplers, as shown in Fig. 4 (b). The filling factor and period of the grating coupling structures are optimized for TE$_{00}$ modes at $\lambda =775$ nm. A 100x microscope objective is used to further enlarge the mode profile emerging from the chip facet, which is then imaged onto a CCD-camera using a second 7x magnifying lens configuration. A polarization filter between the aforementioned magnification components allows us to filter out either the horizontal $H$- or vertical $V$-components of the electric field vector.
Fig. 4. (a) Left: Scanning electron micrograph of fabricated SiN-on-insulator mode converters, optimized using a dynamic binarization function. The simulated electric field distribution ($E_y$-component, perpendicular to direction of propagation). Right: Cleaved edge of the chip with waveguide facet after adiabatic tapering. (b) Measurement setup to determine the relative mode composition using a CCD-camera and magnification lenses. (c) Setup to determine the insertion loss per device. (d) Mode profile with intensity distribution measured along the red lines. The fit used to numerically determine the mode composition is visualized using a dashed line overlay. (e) Measured absolute transmission of the TE$_{20}$ mode for the dynamically and statically binarized devices.
To determine the insertion loss of the device, we employ the measurement setup shown in Fig. 4 (c). In this case a polarization-controlled and single-mode fiber-coupled white light source is used to illuminate the input grating coupler through the fiber array, while the output port is analyzed with a spectrometer. The on-chip nanophotonic circuitry consists of two identical mode conversion devices which are connected through their output-waveguides. The second converter reverses the mode conversion process such that the insertion losses of the two devices accumulate. The total loss of an individual mode converter device can then be inferred after measuring the insertion losses of the grating couplers in similar reference devices. To ensure compatibility between the total transmission of the reference devices and the relative mode composition of the measured device, we determine the mean total transmission and the standard deviation of mode converters and grating couplers for multiple fabricated circuits. The mean insertion loss and the relative mode composition determined by measuring the intensity profile with a CCD-camera, as visualized in Fig. 4 (d), allows us to calculate the absolute conversion efficiency into the TE$_{20}$ mode.
We utilize a least square fitting algorithm to match simulated intensity profiles to the measured data along the cross-section indicated by red horizontal lines in Fig. 4 (d). The amplitudes and relative phases of the individually calculated waveguide modes are free parameters determined by the fitting algorithm, allowing for a quantitative assessment of the mode purity of the output from fabricated devices.
Though limited by the narrow effective wavelength range of the employed grating couplers, we find performance differences between the statically and dynamically binarized device layouts, as depicted in Fig. 4 (e), which show slightly lower overall but otherwise comparable relative efficiency differences to our simulation results (see Fig. 2). Differences between experimental and simulation results are likely due to fabrication imperfections. The uncertainties in wavelength for each measurement shown in Fig. 4 (e) result from the FWHM of the Ti:Sapphire laser, while the uncertainty of the Transmission into the TE$_{20}$ mode results from the standard deviation of the fitting parameters, as well as the standard deviations of the mean total transmission of mode converter and grating coupler devices. Our results provide experimental evidence for significantly enhanced conversion efficiencies from the TE$_{00}$ into the TE$_{20}$ mode for dynamically binarized devices compared to their statically binarized counterparts, as expected from our simulation results.
7. Conclusion
We have introduced a novel dynamic binarization bias for gradient-based continuous inverse design methods in nanophotonics and demonstrated significant performance improvements over previous static implementations. Our approach leverages state-of-the-art mixed-integer optimization techniques and is integrated into established objective-first inverse design frameworks with only minor modifications, while not imposing any new application constraints. As a specific example of a non-trivial photonic integrated circuit component, we show the optimization of a TE$_{00}$ $\rightarrow$ TE$_{20}$ waveguide mode converter in full 3D. Our dynamic binarization method achieves a $9.8$ % increase in conversion efficiency compared to the currently used static binarization approaches.
In future implementations, we anticipate the first phase of our algorithm to benefit from the recently demonstrated method of including an energy constraint in the objective function [28], which at the same time may reduce losses and realize fabrication constraints. Fabrication constraints may also be accounted for by modifying the level-set function like demonstrated in [29] or by level-set related methods which show characteristics similar to mathematical morphology operations shown in [30]. The combination of dynamic binarization functions and neighborhood biases demonstrated in [23] seems promising as a reduction of scattering centers and non-fabricable features might emerge similarly to the method employing an energy constraint mentioned above. Further, we anticipate benefits in the performance of the second phase by employing the L-BFGS-B algorithm [31].
We fabricated and experimentally tested the functionality of the device through direct imaging of the waveguide mode after adiabatic tapering and find qualitative agreement with the simulated conversion efficiency, thus confirming the superiority of the dynamically binarized device. Our work combines progress in two rapidly evolving fields of research, namely mixed-integer optimization and nanophotonic inverse design, thus highlighting the potential for improving photonic integrated circuit performance. We anticipate that our work benefits a wide range of applications in nanophotonics and integrated quantum technology [32].
Funding
Deutsche Forschungsgemeinschaft (CRC-1459 (C05)); Ministerium für Kultur und Wissenschaft des Landes Nordrhein-Westfalen (421-8.03.03.02–130428).
Acknowledgments
We would like to thank the Münster Nanofabrication Facility (MNF) for their support in nanofabrication related matters. C.S. acknowledges support from the Ministry for Culture and Science of North Rhine-Westphalia (421-8.03.03.02–130428). The authors acknowledge support by the German Research Foundation (DFG, CRC 1459).
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. P. D. Trinh, S. Yegnanarayanan, and B. Jalali, “Integrated optical directional couplers insilicon-on-insulator,” Electron. Lett. 31(24), 2097–2098 (1995). [CrossRef]
2. W. Bogaerts, P. De Heyn, T. Van Vaerenbergh, K. De Vos, S. Kumar Selvaraja, T. Claes, P. Dumon, P. Bienstman, D. Van Thourhout, and R. Baets, “Silicon microring resonators,” Laser Photonics Rev. 6(1), 47–73 (2012). [CrossRef]
3. J. Olthaus, P. P. Schrinner, D. E. Reiter, and C. Schuck, “Optimal photonic crystal cavities for coupling nanoemitters to photonic integrated circuits,” Adv. Quantum Technol. 3(2), 1900084 (2020). [CrossRef]
4. J. J. Fijol, E. E. Fike, P. B. Keating, D. Gilbody, J. J. LeBlanc, S. A. Jacobson, W. J. Kessler, and M. B. Frish, “Fabrication of silicon-on-insulator adiabatic tapers for low-loss optical interconnection of photonic devices,” Photonics Packaging and Integration III 4997, 1 (2003). [CrossRef]
5. G. R. Steinbrecher, J. P. Olson, D. Englund, and J. Carolan, “Quantum optical neural networks,” npj Quantum Inf 5(1), 60 (2019). [CrossRef]
6. J. S. Jensen and O. Sigmund, “Topology optimization of photonic crystal structures: a high-bandwidth low-loss t-junction waveguide,” J. Opt. Soc. Am. B 22(6), 1191–1198 (2005). [CrossRef]
7. L. F. Frellsen, Y. Ding, O. Sigmund, and L. H. Frandsen, “Topology optimized mode multiplexing in silicon-on-insulator photonic wire waveguides,” Opt. Express 24(15), 16866–16873 (2016). [CrossRef]
8. E. Hassan and A. C. Lesina, “Topology optimization of dispersive plasmonic nanostructures in the time-domain,” Opt. Express 30(11), 19557–19572 (2022). [CrossRef]
9. L. Sanchis, A. Håkansson, D. López-Zanón, J. Bravo-Abad, and J. Sánchez-Dehesa, “Integrated optical devices design by genetic algorithm,” Appl. Phys. Lett. 84(22), 4460–4462 (2004). [CrossRef]
10. J. Jiang, J. Cai, G. P. Nordin, and L. Li, “Parallel microgenetic algorithm design for photonic crystal and waveguide structures,” Opt. Lett. 28(23), 2381–2383 (2003). [CrossRef]
11. J. C. Mak, C. Sideris, J. Jeong, A. Hajimiri, and J. K. Poon, “Binary particle swarm optimized 2× 2 power splitters in a standard foundry silicon photonic platform,” Opt. Lett. 41(16), 3868–3871 (2016). [CrossRef]
12. J. Lu and J. Vučković, “Objective-first design of high-efficiency, small-footprint couplers between arbitrary nanophotonic waveguide modes,” Opt. Express 20(7), 7221–7236 (2012). [CrossRef]
13. J. Lu and J. Vučković, “Nanophotonic computational design,” Opt. Express 21(11), 13351–13367 (2013). [CrossRef]
14. M. Tutgun, A. Alpkilic, Y. Yilmaz, A. Yeltik, D. Yilmaz, and H. Kurt, “On-chip photonic filter design via objective-first algorithm,” in Optical Modeling and Performance Predictions X, vol. 10743 (SPIE, 2018), pp. 154–161.
15. Y. A. Yılmaz, A. M. Alpkılıç, A. Yeltik, D. Yılmaz, M. Tutgun, and H. Kurt, “Parametric study of wavelength demultiplexers designed via objective-first algorithm,” in Nanophotonic Materials XV, vol. 10720 (SPIE, 2018), pp. 8–15.
16. F. Callewaert, S. Butun, Z. Li, and K. Aydin, “Inverse design of an ultra-compact broadband optical diode based on asymmetric spatial mode conversion,” Sci. Rep. 6(1), 32577 (2016). [CrossRef]
17. M. B. Frish, J. Fijol, E. E. Fike, S. A. Jacobson, P. B. Keating, W. J. Kessler, J. LeBlanc, C. Bozler, M. Fritze, C. Keast, J. Knecht, R. Williamson, and C. Manolatou, “Coupling of single mode fibers to planar si waveguides using vertically tapered mode converters,” in Integrated Photonics Research, (Optica Publishing Group, 2002), p. IFB2.
18. Q. Fang, T.-Y. Liow, J. F. Song, C. W. Tan, M. B. Yu, G. Q. Lo, and D.-L. Kwong, “Suspended optical fiber-to-waveguide mode size converter for silicon photonics,” Opt. Express 18(8), 7763–7769 (2010). [CrossRef]
19. T. L. Penner, H. R. Motschmann, N. J. Armstrong, M. C. Ezenyilimba, and D. J. Williams, “Efficient phase-matched second-harmonic generation of blue light in an organic waveguide,” Nature 367(6458), 49–51 (1994). [CrossRef]
20. Z. Wen, D. Goldfarb, and W. Yin, “Alternating direction augmented lagrangian methods for semidefinite programming,” Math. Prog. Comp. 2(3-4), 203–230 (2010). [CrossRef]
21. M. B. Giles and N. A. Pierce, “An introduction to the adjoint approach to design,” Flow, Turbul. Combust. 65(3/4), 393–415 (2000). [CrossRef]
22. C. M. Lalau-Keraly, S. Bhargava, O. D. Miller, and E. Yablonovitch, “Adjoint shape optimization applied to electromagnetic design,” Opt. Express 21(18), 21693–21701 (2013). [CrossRef]
23. L. Su, A. Y. Piggott, N. V. Sapra, J. Petykiewicz, and J. Vuckovic, “Inverse design and demonstration of a compact on-chip narrowband three-channel wavelength demultiplexer,” ACS Photonics 5(2), 301–305 (2018). [CrossRef]
24. J. Petykiewicz, “Jan Petykiewicz - Software,” https://janp.me/software.html.
25. M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” http://cvxr.com/cvx.
26. M. ApS, “Mosek,” https://www.mosek.com/.
27. P. Bonami, G. Cornuéjols, A. Lodi, and F. Margot, “A feasibility pump for mixed integer nonlinear programs,” Math. Program 119(2), 331–352 (2009). [CrossRef]
28. G. Zhang, D.-X. Xu, Y. Grinberg, and O. Liboiron-Ladouceur, “Experimental demonstration of robust nanophotonic devices optimized by topological inverse design with energy constraint,” Photonics Res. 10(7), 1787–1802 (2022). [CrossRef]
29. M. Zhou, B. S. Lazarov, F. Wang, and O. Sigmund, “Minimum length scale in topology optimization by geometric constraints,” Comput. Methods Appl. Mech. Eng. 293, 266–282 (2015). [CrossRef]
30. F. Wang, B. S. Lazarov, and O. Sigmund, “On projection methods, convergence and robust formulations in topology optimization,” Struct. Multidisciplinary Optimization 43(6), 767–784 (2011). [CrossRef]
31. P.-I. Schneider, X. Garcia Santiago, V. Soltwisch, M. Hammerschmidt, S. Burger, and C. Rockstuhl, “Benchmarking five global optimization approaches for nano-optical shape optimization and parameter reconstruction,” ACS Photonics 6(11), 2726–2733 (2019). [CrossRef]
32. G. Moody, V. J. Sorger, D. J. Blumenthal, et al., “2022 roadmap on integrated quantum photonics,” J. Physics: Photonics 4, 012501 (2022). [CrossRef]
