Analytical level set fabrication constraints for inverse design

Inverse design methods produce nanophotonic devices with arbitrary geometries that show high efficiencies as well as novel functionalities. Ensuring fabricability during optimization of these unrestricted device geometries is a major challenge for these design methods. In this work, we construct a fabrication constraint penalty function for level set geometry representations of these devices. This analytical penalty function limits both the gap size and boundary curvature of a device. We incorporate this penalty in a fully automated optical design flow using a quasi-Newton optimization method. The performance of our design method is evaluated by designing a series of waveguide demultiplexers (WDM) and mode converters with various footprints and minimum feature sizes. Finally, we design and experimentally characterize three WDMs with a 80 nm, 120 nm and 160 nm feature size.


Analytical level set fabrication constraints for inverse design
Dries Vercruysse 1,2 , Neil V. sapra 1 , Logan su 1 , Rahul trivedi 1 & Jelena Vučković 1 Inverse design methods produce nanophotonic devices with arbitrary geometries that show high efficiencies as well as novel functionalities. Ensuring fabricability during optimization of these unrestricted device geometries is a major challenge for these design methods. In this work, we construct a fabrication constraint penalty function for level set geometry representations of these devices. This analytical penalty function limits both the gap size and boundary curvature of a device. We incorporate this penalty in a fully automated optical design flow using a quasi-Newton optimization method. The performance of our design method is evaluated by designing a series of waveguide demultiplexers (WDM) and mode converters with various footprints and minimum feature sizes. Finally, we design and experimentally characterize three WDMs with a 80 nm, 120 nm and 160 nm feature size.
Photonic design is becoming more complex and demanding as a growing number of applications are relying on nanophotonic devices. To meet this challenge, designers increasingly rely on advanced optimization techniques rather than classical photonic design approaches [1][2][3][4] . Instead of tuning relatively simple known geometries with a limited set of parameters, as is traditionally done, these new methods explore devices with completely arbitrary geometries. Capitalizing on the increased degrees of freedom, devices have been designed with extremely small footprints, high efficiencies, and novel functionalities that can not be achieved with classical methods [5][6][7][8][9][10][11][12][13][14] .
A robust design method not only needs to maximize performance, but also needs to ensure that the final device is fabricable. Various fabrication processing limits, e.g., lithography resolution or etch aspect ratio, restrict the minimum feature size that can be achieved practically. Effectively enforcing a minimum feature size on an arbitrary design geometry is challenging yet essential to produce a fully automated design algorithm. Often, optimization problems maintain fabricability by parametrizing simple geometric shapes that cannot violate the minimum feature size, such as hole arrays or pixelated structures 15,16 . This, however, restricts the geometry and thus limits the design space. For a fully arbitrary structure, feature size can be controlled by systematically projecting intermediate results on to the space of fabricable designs 3,17,18 . Such projection steps impede the use of higher order optimization methods, which decrease the optimization time and computational requirements. Another approach is to evaluate a dilated and eroded version of the device during the optimization process. Extensions of this method can allow for fabrication robustness 9,19,20 , or a length scale constraint 21,22 .
In this work, we introduce an analytical fabrication constraint for arbitrary geometries defined by level set functions. This constraint has the advantage that it can be added to the objective as a penalty term, allowing for simultaneous optimization of performance and fabricability with quasi-Newton methods. This is in contrast to previous work where fabrication constraints were enforced through systematic corrections to the device geometry 18 . Expanding on the design methodology established in our group 2,6,11,18 , we present a completely automated optimization process which implements fabricabilty through this penalty term. The performance of this optimization process is evaluated for wavelength demultiplexers (WDMs) and mode converters with varied footprints and minimum feature sizes. Finally, we employ our method to design three WDMs with different minimum feature sizes, which are fabricated and experimentally characterized.

Level Set Fabrication Constraint
To facilitate optimization, a device geometry is often represented by binary pixels on a grid that matches the simulation mesh or the minimum fabricable feature size 7,16,23 . This Manhattan representation, however, limits the design space. Preferably, the design border should be able to move continuously. This can be achieved by directly parameterizing the device boundary as a polygon 24 , or indirectly by using a level set function, as is done in this work. A level set function is a continuous function that defines material where the function is positive, and etch (2019) 9:8999 | https://doi.org/10.1038/s41598-019-45026-0 www.nature.com/scientificreports www.nature.com/scientificreports/ where the function is negative, thereby setting the boundary to be the zero-crossing 25 . An example of a 2D level set function can be seen in Fig. 1(c).
For the device to be considered fabricable, the geometry needs to adhere to a target minimum feature size, d, which the fabrication process can resolve. This requires, firstly, that there are no gaps smaller than the minimum feature size, and secondly, that the radius of curvature is larger than half the minimum feature size ( Fig. 1(b)). Without enforcing these constraints during the optimization, the final designs typically have small features that can be difficult to fabricate, e.g., features smaller than 80nm (Suppl. Information). To ensure fabrication requirements, we introduce two level set specific constraints to the optimization problem: a minimal gap and radius of curvature constraint.
The minimal gap size is enforced by the constraint: where p is a vector describing the parametrization; φ x y p ( , ; ) is the level set function; φ x y p ( , ; ) v and x y p ( , ; ) vv φ are the first and second derivatives, respectively, in the gradient direction 26 ; and β > 0 is a constant. This constraint can be understood intuitively when ignoring the second term. Since the second derivative indicates how quickly the level set function bends back towards the zero-plane, the feature size is controlled by limiting the second derivative based on the value of the function at that point. For a 1D level set function, this constraint is tight for a sinusoidal function with a periodicity 2d, which corresponds to a grating with feature size d. The additional term, β φ ⋅ x y p ( , ; ) v , is added to relax the constraint near the zero-plane, so that the second derivative would not have to be exactly zero (for numerical reasons). β is typically set to 1 3 . An example of where this constraint is violated can be seen in Fig. 2(a).
The minimum radius of curvature, r, is enforced by the constraint: Although the curvature is only relevant at the device boundary (φ = 0), the constraint is evaluated over the entire design region. The constraint is formulated this way because only penalizing curvature at grid points near the boundary results in a highly non-differentiable penalty function, which hinders the optimization process. At points where φ = x y p ( , ; ) 0, the arctan-term in Equation (2) will be π 2 , in which case the constraint becomes r p ( ) . At the extrema, where φ = x y p ( , ; ) 0 v and the radius of curvature is 0, the arctan-term will be 0 effectively removing the constraint. As such Equation (2) constrains the radius of curvature to d 2 at the device border (i.e., at the level set zero-plane) and relaxes it away from the device boundary. Figure 2(b) shows where this constraint is violated in an example structure.
Both constraints can be combined in a penalty function: ( , 0). The penalty weight, ζ, sets the relative importance between the gap and curvature penalty terms, which was set to two in the present implementation. If both constraints 1 and 2 are met, the penalty function will be zero. We observe that optimized devices that meet the constraints do not meet the target minimum feature size, d, set in these equations, but tend to be consistently smaller. Therefore, during optimization, we set d to be 15% higher than the target feature size.

Inverse Design
An optical design problem with the fabrication penalty in Equation (3) can be formulated as: where f EM is a figure of merit for the optical performance as a function of the fields E i and the parametrization vector, p. This optimization problem is subject to Maxwell equations, where ω i and J i are the angular frequency and current sources of mode i, respectively, and p ( ) ε is permittivity as a function of the parametrization. Our inverse design method optimizes the optical problem in Equation (4) in three stages: a continuous optimization stage, where the permittivity of the design region can take any value in between the waveguide and the cladding; a discretization stage, where the continuous result is discretized; and a discrete optimization stage, where the level set representation of the device is optimized with fabrication constraints (Fig. 3a-d) 2,6,11 .
In the first stage (continuous) the parametrization is not a level set function, and as such we omit the penalty term from Equation (4). Yet, the feature size in this stage still needs to be controlled, since small features at the end of the continuous optimization will result in a poor starting condition for the discrete stage. We, therefore, use a coarser grid for the parameterization and interpolate onto the finer simulation grid using cubic interpolation (Suppl. Information). This coarse grid has a pitch of 1.75 times the required minimum feature size. In addition, we apply a sigmoid function over the interpolated result, in order to make the continuous structure more discrete (Suppl. Information). An example of the permittivity distribution for a WDM at the end of a continuous stage can be seen in Fig. 3(b).
The non-fabricable structure is subsequently discretized as illustrated in Fig. 3(b,c) to form the initial guess for the second optimization stage (discrete). Here, the device is strictly composed of regions with the waveguide permittivity and regions with the cladding permittivity, i.e., etched regions. To discretize, we fit a level set parametrization to the continuous optimization result, taking f fab (p) into account to assure fabricability. The fitted result is subsequently used as a starting condition for the discrete optimization stage. In the final discrete stage, we solve the optimization problem shown in Equation (4) (Fig. 3(c,d)). The weight factor τ is increased ten times over the optimization process, and each sub-optimization is solved using the quasi-Newton method L-BFGS-B in order to yield fast convergence.
The penalty term in Equation (4) introduces a challenge to solve the problem directly on a fine simulation grid. Rather than directly parameterizing on the simulation grid, we again use a coarse grid and interpolate on the finer simulation grid to smoothen the optimization landscape (Suppl. Information).

Results and Discussion
WDM with fabrication constraints. The optimization results for a 1300 nm/1550 nm WDM is shown in Fig. 3(e-i). The device was optimized in 2D with a 2.5 μm × 2.5 μm design area, the refractive index of the waveguide and surrounding cladding is 3.48 and 1.44, respectively, and the target minimum feature size, d, was set to 120 nm. Figure 3(e) shows the electromagnetic (EM) objective and penalty function during the optimization process. The continuous optimization stage takes 48 iterations, during which the sigmoid function slope, k, is www.nature.com/scientificreports www.nature.com/scientificreports/ changed twice. This change results in a large increase in the EM-objective at the 16 th iteration and a small increase at the 32 nd iteration. After the continuous stage, the structures is discretized, which again causes a performance decrease. At this point, the optimization relies on a level set parametrization and the problem in Equation (4) is solved. The penalty function for the gap and the curvature in Equation (3) is shown in Fig. 3(a) by the orange curve. Over the discrete optimization stage from iteration 49 to 153, the EM-objective is reduced to 1.5 × 10 −3 , and the fabrication penalty function reaches 0. The final structure (Fig. 3(i)) has an efficiency of 93% and 92% for 1300 nm and 1550 nm, respectively (Fig. 3(f)). The field profiles for both wavelengths can be seen in Fig. 3(g,h). With the penalty function reaching 0, both fabrication constraints are completely met. The device has a minimum gap size of 121.2 nm and a minimum radius of curvature of 63.1 nm.
Minimum feature size vs device footprint. The inverse design method was evaluated by optimizing a series of optical devices with different device parameters. WDMs for 1300 nm/1550 nm were optimized for footprints ranging from 1.5 μm × 1.5 μm to 3 μm × 3 μm and target minimum feature sizes of 80 nm, 120 nm and 160 nm. For each device configuration, we optimized 50 devices starting with different random initial conditions. Of the two optimized wavelengths, the lowest efficiency value is given as the device efficiency in Fig. 4(a-d). The x-axis of these figures shows the minimum feature size of the optimized device, i.e., the minimum of the gap sizes and curvature diameters found in the device geometry (Suppl. Information). In all four figures, the device minimum feature size lies around the target minimum feature size. A small number of the designs still violate the requirement, yet this could be improved by setting the constraint slightly higher, or by increasing the pitch of the coarse grid. Alternatively, when increasing the penalty factor, τ, one could also perform an optimization on the penalty term only to obtain a fabricable structure as a starting condition. Several trends regarding the minimum feature size and the device footprint are visible. For each footprint, the maximum efficiency decreases as the minimum feature size increases. For example, 1.5 μm × 1.5 μm WDMs with an 80 nm feature size can reach efficiencies of 96.6%, but the highest efficiency for the 160 nm feature size is 77.1%. In addition, the spread on the efficiency also becomes larger as the target feature size increases. While all 1.5 μm × 1.5 μm WDMs with a minimum feature size of 80 nm reaches an efficiency between 91.8% and 96.6%, the efficiencies of devices with a 160 nm feature size spans from 11.2% to 77.1%. These trends appear consistently over all the device footprints. Between different device footprints, maximum efficiency increases as the footprint increases. For the 160 nm feature size, the maximum efficiency for a 1.5 μm × 1.5 μm device is 77.1%, while a 3 μm × 3 μm WDM can reach almost 96.7%. Additional optimization metrics, such as average iteration counts and computation time, can be found in the Supplementary Information.
A similar parameter sweep was done for a TE0-to-TE1 mode converter. The mode conversion efficiencies for 50 optimizations with different footprints and minimum feature sizes are shown in Fig. 4(e-h). The same trends as for the WDM can be observed. The effect of the device footprint on the efficiency spread is nevertheless more pronounced. For example, the 160 nm minimum feature size has an efficiency range from 33.6-89.3% for the 1 μm × 1 μm device, yet for the 3 μm × 3 μm device range is already reduced to 88.4-97.8%. The optimization www.nature.com/scientificreports www.nature.com/scientificreports/ landscape of the mode converter presumably has higher performing local minima to which the device can converge, making it more likely to find high-efficiency devices as compared to the WDM problem. 3D WDM designs. Three-dimensional WDM's with a footprint of 3 μm × 3 μm were developed in 220 nm thick SOI with top oxide cladding for 1300 nm/1550 nm wavelengths (Fig. 5). We fabricated a set of devices whose minimum feature size, d, was set to 80nm, 120nm and 160nm. Table 1 summarizes the efficiencies and the fabrication feature sizes of all three devices. These devices adhere to both fabrication constraints, and as expected, devices with a smaller minimum feature size have a higher efficiency (Full spectra can be found in the Suppl. Information). The experimental measurement follows the trend seen in simulations, yet show an overall lower efficiency which can be attributed to fabrication imperfections.

Conclusion
We have developed analytical constraints that limit the gap size and curvature of a structure defined by a level set function. These constraints were incorporated in a fully automated inverse design method in the form of a penalty function. Using this new penalty function, we optimized a series of 2D WDMs and mode converters with varied footprints and minimum feature sizes. Analysis of the feature size of a series of devices that are  www.nature.com/scientificreports www.nature.com/scientificreports/ optimized with certain target feature size, shows a distribution with a lower bound around this target. In addition, the optimization series demonstrate that the spread in the final efficiency of optimized devices increases as the fabrication constraint increases. To improve optimization with large minimum feature sizes, future work might focus on improving initial conditions and on early termination to reduce the computational load of sweeping initial conditions. The validity of this method was demonstrated through simulation in both 2D and 3D, as well as experimental confirmation of device performance. The robustness and flexibility afforded by the presented design approach illustrates the maturity of arbitrary geometry nanophotonic design for widespread use in the integrated photonics community. www.nature.com/scientificreports www.nature.com/scientificreports/