Multi-objective, robust constraints enforced global topology optimizer for optical devices

Abstract

A method for optimization of photonic devices is disclosed. The method includes receiving a set of unconstrained latent variables; mapping the set of unconstrained latent variables to a constrained space to generate a constrained device; calculating the permittivity across each element of the constrained device; determining a permittivity-constrained width gradient based at least partially on the permittivity across each element; and optimizing the set of unconstrained latent variables by at least partially using the permittivity-constrained width gradient.

Description

FIELD OF THE INVENTION

This disclosure relates to systems and methods for optimizing fabrication of optical devices.

BACKGROUND

Optical devices such as metasurfaces may benefit from optimization. Metasurfaces are wavefront shaping devices with subwavelength-scale structures which have been demonstrated for a broad range of applications in sensing, imaging, and information processing. Inverse design using topology optimization has proven to be highly effective for designing metasurfaces with optimized performance, which uses adjoint based gradient computations to iteratively update the dielectric distribution of the device to achieve better performance. A challenge in the inverse design of metasurfaces is to not only optimize device performance but also guarantee device fabricability, for example, enforcing control on the minimum feature size of the pattern so that the fabrication requirement will not go beyond the ability of lithography and etching. Although existing strategies to restrict the minimum feature size such as applying filters and employing analytical penalty functions are able to tackle the problem to some extent, they all lack the ability to strictly avoid the existence of tiny features.

SUMMARY OF THE INVENTION

Various embodiments are directed to a method for optimization of photonic devices. The method receives a set of unconstrained latent variables, maps the set of unconstrained latent variables to a constrained space to generate a constrained device, and calculates the permittivity across each element of the constrained device. The method further determines a permittivity-constrained variable gradient based at least partially on the permittivity across each element and optimizes the set of unconstrained latent variables by at least partially using the permittivity-constrained variable gradient.

In various other embodiments, mapping the unconstrained space to a constrained space includes reparameterization of the unconstrained latent variables.

In still various other embodiments, the reparameterization of the unconstrained latent variables includes solving for a set of constrained width space {w_i} through the following equations:

$v_i=\mathrm{Sigmoid}\left(u_i\right),v_M=1;k_i=v_i\sqrt{1-{\sum }_{j=1}^{i-1}{k}_j^2},k_1=v_1;\mathrm{and}{w}_i=\frac{k_i}{{\sum }_{j=1}^M{k}_j}\left(L-M{w}_{\min }\right)+w_{\min },$ where the set of unconstrained latent variables are {u_i} and wherein w_min is a minimum width of each feature of the constrained width space, M is the total number of features, and Σ_i=1^M w_i=L.

In still various other embodiments, calculating the permittivity across each element of the constrained device includes performing a gray-space relaxation of the set of constrained width space to determine a function of refractive index across each element in the set of constrained width space.

In still various other embodiments, the gray-space relaxation comprises solving for

$\begin{array}{cc}{n}_i\left(x\right)={\left[\exp \left(\pm \frac{x^2-{\left(\frac{w_i}{2}\right)}^2}{T}\right)+1\right]}^{-1},& \end{array}$ where x is the position across the width of a certain element, where x=0 is at the middle of the element, and T controls the binarization of the pattern.

In still various other embodiments, the method further includes determining an efficiency-constrained variable gradient of the constrained device based on the permittivity-constrained variable gradient.

In still various other embodiments, the efficiency-constrained variable gradient is further based on an efficiency-permittivity gradient.

In still various other embodiments, the efficiency-constrained variable gradient is calculated using the following equation:

$\frac{\partial \mathrm{Eff}}{\partial w}=\frac{\partial \mathrm{Eff}}{\partial \epsilon }\times \frac{\partial \epsilon }{\partial w}$ where

$\frac{\partial \epsilon }{\partial w}$ is the permittivity-constrained width gradient and

$\frac{\partial \mathrm{Eff}}{\partial \epsilon }$ is the efficiency-permittivity gradient.

In still various other embodiments, an efficiency-latent variable gradient is calculated at least partially using a latent variable gradient of the latent unconstrained variables in relation to constrained device calculated using the following equations:

$v_i=\mathrm{Sigmoid}\left(u_i\right),v_M=1;k_i=v_i\sqrt{1-{\sum }_{j=1}^{i-1}{k}_j^2},k_1=v_1;\mathrm{and}{w}_i=\frac{k_i}{{\sum }_{j=1}^M{k}_j}\left(L-M{w}_{\min }\right)+w_{\min },$ where the set of unconstrained latent variables are {u_i} and wherein w_min is a minimum width of each feature of the constrained width space, M is the total number of features, and Σ_i=1^M w_i=L.

In still various other embodiments, the efficiency grating in relation to the unconstrained latent variables is calculated at least partially using the following equation:

$\frac{\partial \mathrm{Eff}}{\partial u}=\frac{\partial \mathrm{Eff}}{\partial w}\times \frac{\partial w}{\partial u},$ where

$\frac{\partial w}{\partial u}$ is the constrained width-latent variables gradient and

$\frac{\partial \mathrm{Eff}}{\partial w}$ is the efficiency-constrained width gradient.

In still various other embodiments, the efficiency-permittivity gradient is calculated using an adjoint variable method.

In still various other embodiments, the adjoint variable method includes forward and adjoint simulation.

In still various other embodiments, the method of further includes using a generator to generate the set of unconstrained latent variables; and iteratively updating the generator at least partially based on the permittivity-constrained variable gradient.

In still various other embodiments, the unconstrained latent variables comprise unconstrained latent widths and mapping the set of unconstrained latent variables comprises mapping the unconstrained latent widths to the constrained space to generate constrained widths. The method further includes receiving a set of unconstrained latent permittivities and mapping the set of unconstrained latent permittivities to a constrained space to generate constrained average permittivities of the constrained device. Calculating the permittivity profile across each element of the constrained device is further based on the constrained widths and the constrained average permittivities of the constrained device.

In still various other embodiments, the method further includes determining a permittivity-constrained average permittivity gradient of the constrained device and optimizing the set of unconstrained latent permittivities at least partially based on the permittivity-constrained average permittivity gradient of the constrained device.

In still various other embodiments, the method further includes binarizing the permittivity across each element of the constrained device, where determining a permittivity-constrained variable gradient is based on the binarized permittivity across each element.

Various embodiments are directed to a method for optimization of photonic devices, the method including receiving a set of variables including a set of refractive indices including a set of corresponding thicknesses, calculating a reflection gradient with respect to the set of refractive indices and their set of corresponding thicknesses, and optimizing the set of variables at least partially using the reflection gradient to obtain an optimized device including a set of optimized refractive indices with corresponding optimized thicknesses.

In various other embodiments, the method further includes calculating the reflection spectrum of the set of variables, where calculating the reflection gradient includes calculating a partial differential of the reflection spectrum.

In still various other embodiments, calculating the reflection spectrum comprises using a transfer matrix method.

In still various other embodiments, the method further includes: using a generator to generate the set of variables; using the set of variables to calculate a corresponding reflection gradient; and iteratively updating the generator at least partially using the reflection gradient.

BRIEF DESCRIPTION OF THE DRAWINGS

The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawing(s) will be provided by the Office upon request and payment of the necessary fee.

The description will be more fully understood with reference to the following figures and data graphs, which are presented as various embodiment of the disclosure and should not be construed as a complete recitation of the scope of the disclosure.

FIG. 1 illustrates an overview of reparameterization for an example metagrating design.

FIG. 2A illustrates an example schematic of metagratings including silicon nanoridges in accordance with an embodiment of the invention.

FIG. 2B illustrates the permittivity profile ε(x) of the metagratings disclosed in FIG. 2A.

FIG. 2C illustrates a cross sectional view of metagratings described in connection with FIG. 2A.

FIG. 3 illustrates a schematic of an inverse design method with reparameterization in accordance with an embodiment of the invention.

FIG. 4 is flowchart of the reparameterized gradient-based local optimizer in accordance with an embodiment of the invention.

FIG. 5 illustrates an example pseudo-code describing the algorithm of FIGS. 3 and 4.

FIG. 6A illustrates a graphical representation of the effect of different levels of T.

FIG. 6B is various graphs illustrating a study between the efficiency of binarized structure and gray-scale structure of devices with similar width for different T values.

FIG. 7A is a plot illustrating an optimization process using the method described in connection with FIGS. 3 and 4 with constraints on the minimum feature size of 40 nm.

FIG. 7B is a plot illustrating the efficiency of the best optimized devices with constraints on different minimum feature sizes and their corresponding device patterns.

FIG. 8A illustrates histograms of device efficiency and minimum feature size for 100 locally optimized metagratings.

FIG. 8B illustrates charts describing efficiency, binarization, and T as a function of iteration number for a representative device optimization, together with device layouts and gradients for specific iterations.

FIG. 9A illustrates an example flowchart of reparameterized GLOnets in accordance with an embodiment of the invention.

FIG. 9B illustrates pseudocode for the reparameterization technique of FIG. 9A.

FIG. 10 illustrates efficiency, diversity, and binarization as a function of iteration number for a representative GLOnets training run.

FIG. 11 illustrates histograms of generated device efficiencies and images of representative physical devices as a function of iteration number.

FIG. 12A illustrates a benchmark comparison of devices designed using reparameterized GLOnets and the reparameterized local optimizer, for devices comprising three silicon ridges and parameters matching those in FIG. 7A.

FIG. 12B illustrates a plot of GLOnets-designed device efficiencies for differing device topologies. The topology of the globally optimized device strongly depends on the MFS constraint.

FIG. 13 illustrates a schematic for an inverse design method with reparameterization which receives input from a generator in accordance with an embodiment of the invention.

FIG. 14A shows an example generator training process for a 65° deflector at wavelength 850 nm with a minimum feature size of 30 nm.

FIG. 14B illustrates globally optimized device efficiencies as a function of minimum feature size.

FIG. 14C illustrates a plan view of an optimized metagrating with minimum feature size of 20 nm with the highest absolute efficiency.

FIG. 14D illustrates a plan view of an optimized metagrating with minimum feature size of 50 nm with the highest absolute efficiency.

FIG. 15A illustrates various schematics of metagratings including fabrication imperfections modeled by the erosion and dilation of device layouts, which is parametrized by δw.

FIG. 15B illustrates various plots comparing performance of devices globally optimized with and without robustness.

FIG. 16A is a plot illustrating efficiency of a globally optimized device generated with reparameterization and accounting for fabrication imperfections at a minimum feature size of 30 nm.

FIG. 16B illustrates a plan view of the optimized device without fabrication robustness.

FIG. 16C illustrates a plan view of the optimized device with fabrication robustness and a minimum feature size of 30 nm.

FIG. 17 illustrates an example schematic of a photonic thin film stack.

FIG. 18 illustrates a flow chart of an example global optimization algorithm using a TMM in accordance with an embodiment of the invention.

FIG. 19 illustrates a generative neural network used to produce a distribution of thin film stack configurations.

FIG. 20 illustrates a cross sectional view of the resultant structure from the optimization.

FIG. 21 illustrates various plots of the results for various wavelengths across angles between 0° and 60°.

FIG. 22 illustrates the results from another example optimization performed to an incandescent light bulb filter that transmits visible light and reflects infrared light.

FIGS. 23 and 24 illustrate the broadband reflection characteristics of a 45-layer generative neural network optimized device.

FIG. 25 illustrates the emission intensity spectrum of the light bulb with and without the thin film filter.

FIG. 26 is a chart comparing the results of a comparison of thin film filters generated using a residual global neural network, a deep global neural network, and a generic algorithm.

FIG. 27A illustrates a flowchart of a reparameterization process for a variable-topology global optimizer.

FIG. 27B illustrates schematics of the piecewise transformation used for device reparameterization.

FIG. 27C illustrates various device layouts over the course of GLOnets training for a representative optimization run.

FIG. 27D illustrates a plot of device efficiencies as a function of the minimum feature size (MFS) constraint, designed using fixed-topology and variable-topology reparameterized GLOnets.

FIG. 27E illustrates scatter plots of device efficiency and MFS for devices designed using variable-topology reparameterized GLOnets, original GLOnets, and adjoint-based local optimization, for three different MFS constraints.

FIG. 28 illustrates a piecewise transformation used in variable-topology reparameterized GLOnets.

FIG. 29 illustrates an example algorithm for newly configured GLOnets where physical devices in this form are generated using two separate generative networks, one for section widths and the other for refractive indices being optimized.

FIG. 30 illustrates a schematic for an optimization process including a binarization step in accordance with an embodiment of the invention.

FIGS. 31A and 31B illustrate a comparison between the reparameterization processes described in connection with FIGS. 3-5 and the reparameterization process described in connection with FIG. 30.

FIG. 32 illustrates an example of optimization system that optimizes photonic devices while enforcing fabrication constraints in accordance with an embodiment of the invention.

FIG. 33 illustrates an example of an optimizer element that executes instructions to perform processes that optimizes photonic devices while enforcing fabrication constraints in accordance with an embodiment of the invention.

FIG. 34 illustrates an example optimization method in accordance with embodiments of the invention.

DETAILED DESCRIPTION

Nanophotonic devices may be capable of manipulating and guiding electromagnetic waves propagating in free space and on chip, and they have a broad range of applications in imaging, sensing, and optical communications. Among the most effective methods to design these devices is gradient-based topology optimization, which has been used to realize metagratings, metalenses, and on-chip photonic devices that utilize complex electromagnetic wave dynamics to achieve exceptional performance. Topology optimization is performed by discretizing the device structure into small voxels, initializing each voxel with grayscale dielectric permittivity values c, and then iteratively modifying these permittivity values in a manner that improves a figure of merit (FoM). These modifications are based on gradient terms

$\frac{\partial \left(\mathrm{FoM}\right)}{\partial \epsilon }$ that can be computed for each pixel using the adjoint variables method or autodifferentiation. Topology optimization can be performed in the context of local optimization, in which local gradients are directly used to perform gradient descent on grayscale device structures, or global optimization with global topology optimization networks (GLOnets), in which local gradient calculations are combined with generative neural network training to perform global population-based optimization.

A critical concern in the inverse deign methods for optical devices may be guaranteeing device fabricability, making it desirable to incorporate practical experimental constraints such as the minimum feature size or robustness to fabrication imperfections. Existing techniques such as applying threshold filters and adding penalties to the figure of merit may be able to push the output devices towards the desired constrained design space. However, they may lack the ability to enforce hard constraints such as avoiding the appearance of tiny features smaller than the minimum specified feature size.

Further, a critical issue concerning topology-optimized devices is the practical implementation of hard geometric constraints imposed by experimental considerations. A particularly important geometric constraint is minimum feature size (MFS), which arises due to limitations in lithography patterning resolution and etching aspect ratio. The imposition of MFS constraints also enables proper base patterns to be defined and used in algorithms that enforce fabrication robustness, in which base patterns are codesigned with their geometrically eroded and dilated forms. Without these constraints, topology-optimized devices often possess complex geometric shapes with very small feature sizes, making them difficult if not impossible to experimentally fabricate in a reliable manner.

Current methods to impose MFS constraints fall in one of three classes. The first is to set the device voxel dimensions or spacing between features to match the desired MFS. While this method works, it adds significant granularity to the device design space, limiting the final device performance. The second method is to use regularization terms that penalize the FoM when the MFS constraint is violated. While this technique will generally push devices toward regions of the design space that satisfy the desired constraints, it does not guarantee their enforcement. A third method is to optimize the device in an unconstrained manner and then use threshold filters to incorporate constraints. This method can be applied during the iterative optimization process or after unconstrained optimization is performed, and while it has the potential to work well, it only works when optimized devices in the unconstrained space locally map onto high performance optima in the constrained space.

Various embodiments of this disclosure relate enforcing strict geometric constraints within an inverse optimizer, while maintaining the fine pixel-level granularity of the physical design space, by reparameterizing the physical design space itself. FIG. 1 is an overview of reparameterization for an example metagrating design. The freeform metagratings including silicon nanoridges are described by 1×N vectors. To enforce hard physical constraints, mathematical transformations are used to map devices represented in an unconstrained latent space to a constrained device space. The idea, outlined in FIG. 1, is to initially define devices in a latent space free of constraints and then transform these devices to a physical device space with hard constraints using analytic and differentiable mathematical transformations. These devices are evaluated and optimized within this physical device space, and modifications to these devices are mapped back to their latent representations. The implementation of reparameterization for the design of high efficiency metagratings using local optimization are demonstrated based on the adjoint variables method, and GLOnets-based global optimization. While much of the disclosure relates to the application of reparameterization to the enforcement of feature size constraints, the concept provides a general formalism to incorporate physical information into black box inverse design tools.

Further, various embodiments of this disclosure relate to a method that accurately imposes constraints through reparametrization. In some embodiments, the method may impose hard constraints in the optimization process through reparametrizing the problem. The pixelated device pattern is reparametrized by a new set of parameters where the constraints in the old design space are naturally eliminated, and optimization is performed in this new design space instead. This method may be used to design metasurfaces such as metagratings including silicon nanoridges with strict restrictions on minimum feature size.

In some embodiments, via a sequence of mathematical transformations, the method uses variables in an unconstrained latent space to reparametrize the original constrained device space such that the constraints are naturally imposed through the reparameterization process. Adjoint based optimization may then be performed with this newly introduced unconstrained variable space where performance gradients are first calculated regarding the original constrained device space and then backpropagated to the unconstrained latent space.

Example Photonic Devices

Photonic devices including metasurfaces such as high efficiency metagratings including silicon nanoridges may be fabricated with accurate control on the minimum feature size are designed. Silicon nanoridges may deflect normally incident light to the +1 diffraction order. FIG. 2A illustrates an example schematic of metagratings including silicon nanoridges in accordance with an embodiment of the invention. Deflection efficiency of the metagrating may be the figure of merit here which may be defined as the electromagnetic power deflected into the first diffraction order normalized with the total input power with normal incidence. In the original device space, the metagrating may be defined by the permittivity profile c(x) within each grating period which may denotes the material distribution at a pixel level, where c(x) is normalized so that 0 represents air and 1 represents silicon. FIG. 2B illustrates the permittivity profile ε(x) of the metagratings disclosed in FIG. 2A. In adjoint-based topology optimization, ε(x) may be iteratively adjusted to increase deflection efficiency according to the gradient calculated from a forward and an adjoint simulation.

At all stages of optimization, an individual metagrating period may be defined by a permittivity profile ε=ε(X), which denotes the material distribution at each voxel located at x within an individual grating period. ε(x) is normalized so that 0 represents air and 1 represents silicon. The metagrating period has a width L that is subdivided into N=256 voxels, so that both x and E are vectors with dimension N. Our FoM is deflection efficiency, which is defined as the electromagnetic power deflected into the +1 diffraction order normalized to the total incident beam power. The deflection efficiency is denoted as Eff=Eff(ε) and calculated with a rigorous coupled-wave analysis (RCWA) electromagnetic solver such as Reticolo.

The objective is to find a permittivity profile that maximizes the deflection efficiency, which can be cast as the following optimization:

• • maximize Eff(ε)
  • subject to εϵ{0,1}^N

To define feature size constraints on the pixelized pattern ε(x), the pixelized pattern may be transformed to a vector of discrete width values, w, with which MFS constraints can be readily defined and incorporated. For devices of fixed topology, each element in w represents ridge widths and air gap values which may generalize to devices of varying topology. By reframing the pixelized patterns to a vector of geometric values, the original optimization problem may be reformulated for devices of fixed topology to be:

$\mathrm{maximize}\mathrm{Eff}\left(w\right)\mathrm{subject}\mathrm{to}{w}_i\ge w_{\min },i=1,2,\dots ,M\sum _{i=1}^M{w}_i=L$ w_min denotes the MFS, w_i denotes the width of each structural feature, which can be a silicon ridge or air gap, and M denotes the total number of silicon and air features. The total device width should be equal to the grating period L.

By first transforming the pixelated pattern into a pattern represented by the ridge width and ridge separation values, the minimum feature size can be described as w_i≥w_min with i=1, 2, . . . , M where w_min denotes the minimum feature size, w_i denotes the width of each bar (silicon ridge or air gap), and M denotes the total number of bars. FIG. 2C illustrates a cross sectional view of metagratings described in connection with FIG. 2A. Another constraint to this width space is Σ_i=1^N w_i=L where L denotes the periodicity of the metagrating.

Processes Including Reparameterization

The reparametrized inverse design process may be used to design high efficiency metagratings with accurate control on the minimum feature size. FIG. 3 illustrates a schematic of an inverse design method with reparameterization in accordance with an embodiment of the invention. FIG. 4 is flowchart of the reparameterized gradient-based local optimizer in accordance with an embodiment of the invention. FIG. 5 illustrates an example pseudo-code describing the algorithm of FIGS. 3 and 4. Latent vectors are transformed to physical devices, and gradients calculated for physical devices are backpropagated to update the latent vectors. Device patterns may be reparametrized by the base parameters (u₁, u₂, . . . , u_N-1) without constraints, and constraints may be naturally imposed through mapping from base parameters to device patterns. Shown in the flow chart of the upper panel, each iteration includes a forward simulation from parameters to objective and a backward propagation of gradients to update the parameters. The method may start with a base parameter set {u_i} that may have no constraints and map it into a device pattern described by a pixilated refractive index profile n(x) with restrictions on various parameters such as the minimum feature size through a series of transformations, also known as the reparameterization, where the constraints may be naturally imposed. The refractive index profile n(x) may be normalized so that 0 represents air and 1 represents silicon. The constraint-enforced device pattern may be inputted into an electromagnetic simulator such as a rigorous coupled-wave analysis solver to evaluate various objectives such as the deflector efficiency in the +1 diffraction order. As for the backward calculation, the gradient of the refractive index profile may be first attained from the adjoint variable method and then backpropagated to the base parameters. Therefore, gradient descent may be used to update those base parameters instead of the device pattern.

The reparameterization process of mapping {u_i} into n(x) is elaborated in the lower panel of FIG. 3. The base parameter set {u_i} may be mapped through a set of differentiable functions into {w_i} which describes the width of the air/silicon bars which may describe the refractive index profile of the metagrating. The reparameterization technique may incorporate experimental constraints and may include but is not limited to two steps:

• • the first step may be to map from an unconstrained variable space {u_i} to a constrained width space {w_i} where the constraint is naturally imposed;
  • the second step may be to map from the width space {w_i} to the original device space ε(x) where a binarized pattern is relaxed to gray-scale such that the gradient is able to propagate.

As illustrated in the leftmost box on the lower panel of FIG. 3, the reparameterization process may start with an unconstrained variable space {u_i} with i=1, 2, . . . , M−1. A set of differentiable transformation functions may be performed:

$\begin{array}{cc}{v}_i=\mathrm{Sigmoid}\left(u_i\right)=\frac{e^{u_i}}{e^{u_i}+1},v_M=1,v_M=1,i=1,\dots ,M-1& \left(1\right)\end{array}$ $\begin{array}{cc}{k}_i=v_i\sqrt{1-{\sum }_{j=1}^{i-1}{k}_j^2},k_1=v_1,i=2,\dots ,M& \left(2\right)\end{array}$ $\begin{array}{cc}{w}_i=\frac{k_i}{{\sum }_{j=1}^M{k}_j}\left(L-M{w}_{\min }\right)+w_{\min },i=1,\dots ,M& \left(3\right)\end{array}$

The transformation of u to w may be defined as w=ƒ(u, w_min). In this manner, MFS constraints on the metagrating pattern, have been imposed onto the width vector w.

Through the transformations described in the leftmost box on the lower panel of FIG. 3, constraints on the metagrating pattern may be imposed as w_min≤w_i and Σ_i=1^N w_i=L, where N is the total number of air/silicon bars, and L is the period of the grating. Next, a gray-scale relaxation may be performed using equation 4 to map from {w_i} to ε=n(x)=f(x,w,T):

$\begin{array}{cc}{n}_i\left(x\right)={\left[\exp \left(\pm \frac{x^2-{\left(\frac{w_i}{2}\right)}^2}{T}\right)+1\right]}^{-1}& \left(4\right)\end{array}$

FIG. 6A illustrates a graphical representation of the effect of different levels of T. Equation 4 characterizes the gray-scale pattern, where the ± sign indicates + for Si and − for air, and T controls the binarization of the pattern and decreases during the course of the iterative optimization so that the refractive index is binarized to either air or Si. T is a hyperparameter, analogous to a temperature term, which controls the degree of binarization of the pattern. {tilde over (x)}⁽ⁱ⁾ represents the position of pixels corresponding to the ith section, and it is defined such that the center of w_i falls onto the position of {tilde over (x)}⁽ⁱ⁾=0. Optimizing in a binarized space can be extremely hard as it is hard to propagate the efficiency gradient from the original device space ε(x) to the width space {w_i} in the iterative optimization process. Thus, equation 4 may be used to relax the binarized structure into grey scale. Although the more gray-scale the easier the gradient can be propagated, the landscape of the design space can be changed dramatically when the structure is relaxed from binarized to gray-scale. Thus, there is a tradeoff between propagating the gradient and how accurate the design space may be described with a gray-scale relaxation. The top images of FIG. 6A illustrate the reparameterized pattern after gray-scale relaxation. The bottom plots of FIG. 6A illustrate various results after gray scale relaxation for different T values. FIG. 6B is various graphs illustrating a study between the efficiency of binarized structure and gray-scale structure of devices with similar width for different T values. As illustrated, for a population of 60,000 different devices and for T≤1000, higher than 0.9 correlation may be achieved, indicating a high-level of similarity between the gray-scale space and the binarized space.

With a physical device in hand,

$\frac{\partial \mathrm{Eff}}{\partial \epsilon }$ may be evaluated using the adjoint variables method, in which forward and adjoint simulations are performed using the RCWA simulator. This gradient term may be then used to calculate

$\frac{\partial \mathrm{Eff}}{\partial u}$ using back-propagation, which may be based on the chain rule. Backpropagation can be used in the algorithm because the entire reparameterization process can be described as a continuous and differentiable computation graph. In practice, backpropagation may be performed by programming the algorithm in PyTorch and using built-in autodifferentiation packages. Finally, u may be updated using the adaptive moment estimation (Adam) algorithm, which may be a variant of gradient descent, and the entire process is iteratively repeated.

A metagrating may be designed that deflects normally incident TM-polarized light at a wavelength of 850 nm to 65 degrees. FIG. 7A is a plot illustrating an optimization process using the method described in connection with FIGS. 3 and 4 with constraints on the minimum feature size of 40 nm. The plot on top illustrates deflection efficiency and binarization over the course of the optimization process, and corresponding device patterns. The topology of the metagratings may be fixed at 6 air/Si bars and the pattern may gradually become binary over the course of optimization. The deflection efficiency of the final design reaches 89.4%.

FIG. 7B is a plot illustrating the efficiency of the best optimized devices with constraints on different minimum feature sizes and their corresponding device patterns. This plot illustrates the design of a series of metagratings with different minimum feature sizes. As illustrated, when increasing the minimum feature size constraint from none to 100 nm, the best optimized devices experience an efficiency drop from 97.3% to 71.1% while keeping the similar geometry features. The best device achieved with no minimum feature size constraint may performance at a 97.3% efficiency.

In some embodiments, the method may be used to enforce fabrication constraints through reparameterization. In the new reparameterized design space, those constraints may be naturally eliminated, and fabrication optimization may be achieved. The device designed by the method may be able to achieve extremely high performance and strictly satisfy fabrication constraints. Implementations of this design scheme may include variations in the topology, enrichment in constraints, and diversification of structural geometries.

Reparameterized local optimization may be used to design metagratings that deflect normally incident transversemagnetic-polarized light at a wavelength of 850 nm to 65 degrees. In one example, the MFS constraint may be set to 60 nm, the topology is fixed to contain three silicon ridges (M=6), and the device thickness is 335 nm. The refractive index of silicon may be used and only the real part of the index may be used to simplify the design problem. Complex refractive indices can be readily incorporated into the optimization procedure without loss of generality. For the disclosed system, the role of the imaginary part of the index may be small, and its inclusion results in a less than 2% drop in the total deflection efficiency for all of the optimized devices.

FIG. 8A illustrates histograms of device efficiency and minimum feature size for 100 locally optimized metagratings. The inset shows the device operating parameters. FIG. 8A shows a wide distribution of efficiency values, indicating that the design space is highly nonconvex with many local optima. The MFS distribution of these devices, also plotted in FIG. 8A, indicates that all devices have feature sizes larger than 60 nm, demonstrating that reparameterization enforces hard geometric constraints.

FIG. 8B illustrates charts describing efficiency, binarization, and T as a function of iteration number for a representative device optimization, together with device layouts and gradients for specific iterations. Over the course of device optimization, the device efficiency increases, the pattern becomes more binary, and gradient contributions become more localized to the air-silicon boundaries. The initial device pattern is random and set to have a large T. In the high T regime, the device possesses substantial regions of grayscale dielectric material, and the degree of binarization, defined as Σ_j=1^N|2ε(x_j)−1|/N, is low (as the device becomes more binary, this term gets closer to 1). In addition,

$\frac{\partial \epsilon }{\partial w}$ spans many voxels, which translate to large

$\frac{\partial \mathrm{Eff}}{\partial w}=\frac{\partial \mathrm{Eff}}{\partial \epsilon }\times \frac{\partial \epsilon }{\partial w}$ upon backpropagation. These gradients allow relatively large changes to the device layout to be made each iteration in the early stages of optimization, while the algorithm is broadly searching for a local optimum. It is noted that the initial magnitude of T cannot be too large and must be judiciously chosen because the grayscale and binary design spaces are different, and these differences become more substantial as T increases. Sufficient correlation between these two design spaces is required for gradients within the grayscale design space to reliably improve the device over the course of binarization. As the local optimizer evolves, T is gradually reduced and the device becomes more binarized. At this stage,

$\frac{\partial \varepsilon }{\partial w}$ becomes more localized to the air-silicon interfaces and the gradients provide fine-tuning of the device layout, in a manner akin to conventional boundary optimization. Upon the completion of optimization, T is decreased to zero, and the device is fully binarized and contains only silicon and air material.

Reparameterization not only is able to enforce robust physical constraints, it also reduces the computational cost for the optimization process because it serves to reduce the dimensionality of the physical design space itself. Through reparameterization, N dimensional (N=256) pixelized device patterns are reduced to a M−1-dimensional latent vector representations, leading to over an order of magnitude reduction in the dimensionality of the design parameters. This results in a physical design space that is dramatically smaller with fewer local optima, yielding a simplified optimization problem that requires less computational cost for both local and global optimization algorithms.

Embodiments Including Generators

In some embodiments, reparameterization may be applied to GLOnets to enable population-based global optimization of fixed-topology devices with hard geometric constraints. FIG. 9A illustrates an example flowchart of reparameterized GLOnets in accordance with an embodiment of the invention. A distribution of latent vectors is produced by a generative neural network and reparameterized to constrained physical devices. Network weights are iteratively updated to minimize the loss function using backpropagation. FIG. 9B illustrates pseudocode for the reparameterization technique of FIG. 9A. Reparameterized GLOnets involve the iterative generation of devices from the network, evaluation of a loss function, and update of network weights that minimize the loss function. Unlike the original GLOnets architecture, in which a generative neural network produces a distribution of physical devices, the generator in the reparameterized GLOnets generates a distribution of latent vectors u from uniformly distributed noise vectors z∈U^M-1(−1,1). This nonlinear mapping from z to u is performed using a series of fully connected layers: u=G_ϕ(z), where ϕ are the network weights. From a statistical perspective, the generator maps a uniform distribution of noise vectors to a probability distribution:G_ϕ:U^M-1(−1,1)→P_ϕ(u)  (5)where P_ϕ(u) denotes the probability of generating u in the latent device space. The reparameterization process then follows, which maps u to the physical device pattern ε using the same formalism specified for the reparameterized local optimizer.

The next step is the evaluation of the loss function using metrics calculated from the generated devices. The loss function is engineered so that minimizing the loss function maximizes the probability that the neural network generates the optimal latent vector u*, which maps onto the globally optimized device ε*. For reparameterized GLOnets, the loss function is defined to be:

$\begin{array}{cc}L\left(\epsilon ,g_{\epsilon },\mathrm{Eff}\right)=-\frac{1}{K}{\sum }_{k=1}^K\frac{1}{\sigma }\exp \left(\frac{{\mathrm{Eff}}^{\left(k\right)}}{\sigma }\right){\epsilon }^{\left(k\right)}·g_{\epsilon }^{\left(k\right)}& \left(6\right)\end{array}$

The efficiencies, Eff, and efficiency gradients,

$g_{\epsilon }^{\left(k\right)}=\frac{\partial \mathrm{Eff}}{\partial \epsilon },$ of the physical devices are calculated with forward and adjoint electromagnetic simulations. K is the batch size, and σ is a hyperparameter that biases network training toward devices possessing relatively high efficiencies and large gradients. To minimize the loss function, backpropagation is performed to modify network weights ϕ in the generator. As the mapping functions that link the noise vector, latent vector, and physical device pattern profile are continuous and differentiable, backpropagation is performed via the chain rule in a straightforward manner. Upon the completion of network training, the network generates latent vectors that map to physical devices clustered around the global optimum.

Reparameterized GLOnets are used to optimize metagratings with the same specifications as those designed by local optimization which are discussed in connection with FIGS. 7A and 7B. FIGS. 10 and 11 are summaries of results for a single network trained over the course of 100 iterations with a batch size of 100 devices. FIG. 10 illustrates efficiency, diversity, and binarization as a function of iteration number for a representative GLOnets training run. FIG. 11 illustrates histograms of generated device efficiencies and images of representative physical devices as a function of iteration number. Over the course of optimization, the devices become more binary, and the distribution of generated devices narrows and shifts to high efficiency values. At the start of network training, the reparameterized generator has no knowledge of good metagrating designs and outputs a diversity of grayscale devices with modest efficiencies. Device diversity for a given batch is computed as

$\frac{1}{M}{\sum }_{i=1}^M\sqrt{\frac{{\sum }_{k=1}^K{\left(w_i^{\left(k\right)}-{\overline{w}}_i\right)}^2}{K-1}}.$ Over the course of network training, the best generated device and the overall device distribution shift toward higher efficiency regimes (top panel, FIG. 10), and the devices become more binarized and less diverse (bottom panel, FIG. 10). Upon the completion of network training, the best sampled device is simulated to have an efficiency of 90%, which is higher than the best locally optimized device in FIG. 7A.

The evolution of the device distribution over the course of network training is examined in more detail in FIG. 11, which shows representative device layouts and performance histograms at different moments of optimization. In early iterations, the reparameterized generative network outputs a diverse distribution of grayscale devices, indicating that the network is broadly sampling the device design space in its search for good optima. The performance histograms show mostly modest to low efficiency devices. At later iterations, the network has identified and converged to more narrow regions of the design space and outputs less diverse distributions of devices. The final network produces binary devices that are clustered around the global optimum.

FIG. 12A illustrates a benchmark comparison of devices designed using reparameterized GLOnets and the reparameterized local optimizer, for devices comprising three silicon ridges and parameters matching those in FIG. 7A. To further map out the metagrating design space, sets of devices with differing MFSs and topologies are optimized. Globally and locally optimized devices containing three silicon ridges, and the results are summarized in FIG. 12A. Each data point in the reparameterized GLOnets curve represents the best device from a single trained network, while each data point in the reparameterized local optimizer curve represents the best of 100 individual optimizations. All data points are calculated using the same number of RCWA simulations. A comparison between these curves indicates that GLOnets consistently outperform the local optimizer for nearly all MFSs, demonstrating its ability to effectively perform global optimization. Both curves show a decrease in efficiency as the MFS constraint is increased, indicating that small features are advantageous in enhancing light diffraction efficiency. The presence of small features in high performing devices is visualized in images of the device layouts, which show that the best globally optimized devices possess at least one feature near or at the MFS length scale. For MFS constraints smaller than 20 nm, the reparameterized GLOnets generate the same optimal device with a deflection efficiency of 98%. This plateauing follows because the unconstrained globally optimal device has a MFS of 26 nm.

FIG. 12B illustrates a plot of GLOnets-designed device efficiencies for differing device topologies. The topology of the globally optimized device strongly depends on the MFS constraint. Globally optimized efficiencies for devices with differing topologies are summarized in FIG. 12B, together with layouts of the best devices for a given MFS constraint. These data indicate that the topology and layout of the best device is a strong function of MFS. For no constraints up to a 30 nm MFS, devices containing three silicon ridges produce the highest overall efficiencies. Devices with four and five ridges also have efficiencies near the global optimum and have the benefit of larger feature sizes. Above a MFS of 30 nm, the efficiencies of three ridge devices drop dramatically and the four and five ridge devices are optimal. Above an MFS of 90 nm, two ridge devices have the highest overall efficiency.

FIG. 13 illustrates a schematic for an inverse design method with reparameterization which receives input from a generator in accordance with an embodiment of the invention. This inverse design method is similar to the method described in connection with FIG. 3 and the overlapping steps will not be repeated. A generator may be used to generate base parameters (u₁, u₂, . . . , u_N-1) without constraints, and constraints may be naturally imposed through mapping from base parameters to device patterns. Instead of directly generating a distribution of pixelized patterns, the generator uses generative networks to generate a distribution of devices in the latent variable space {u_i}. The generator iteratively trains with gradients backpropagated from the latent variables. Since the dimensionality of the latent variable space is reduced by orders of magnitude when compared to the original device space described with pixels, the generative network may be trained with a reduced batch size while still maintaining a thorough search of the design space. The differentiable nature of all the reparameterization functions improves the feasibility of gradient backpropagation in the proposed scheme.

FIG. 14A shows an example generator training process for a 65° deflector at wavelength 850 nm with a minimum feature size of 30 nm. The purple curve shows the efficiency of the best performed device out of a batch of 100 devices over the course of training. By the end of a 100-itertation-training, the maximal efficiency goes to above 0.95. Meanwhile, the binarization of the device shown as the green curve ends at 1, indicating a fully binarized pattern.

FIGS. 14B-14D illustrates various results of reparametrized with various minimum feature sizes. FIG. 14B illustrates globally optimized device efficiencies as a function of minimum feature size. Different colored curves indicate different topologies of the metagrating where M=the total number of bars of air and silicon. As illustrated, with a fixed topology, the efficiency of the globally optimized devices may decrease as the minimum feature size constraint increases. With different minimum feature size constraint, the globally optimized devices can have different topology, and certain topology may be favored at certain minimum feature sizes. For example, the globally optimized devices has three silicon bars or ridges (M=6) at 20 nm minimum feature size and four silicon bars or ridges (M=8) at 50 nm minimum feature size.

FIG. 14C illustrates a plan view of an optimized metagrating with minimum feature size of 20 nm with the highest absolute efficiency. This metagrating has three silicon bars 602 and three air bars 604 and has an absolute efficiency of 98.1%.

FIG. 14D illustrates a plan view of an optimized metagrating with minimum feature size of 50 nm with the highest absolute efficiency. This metagrating has four silicon bars 606 and four air bars 608 and has an absolute efficiency of 96.5%.

Fabrication Imperfections Enforced with Reparameterization

Devices that satisfy MFS constraints still may not guarantee good experimental performance as the device could be sensitive to other types of fabrication imperfections. Robustness to fabrication imperfections can be readily incorporated into reparameterized GLOnets. In some embodiments robustness may be incorporated in fixed-topology reparameterized GLOnets. In some embodiments, robustness may be incorporated in variable-topology reparameterized GLOnets. Robustness criteria are practically important due to geometric imperfections arising from all experimental fabrication processes. In some cases, geometrically eroded and dilated versions of the devices in the FoM may be used to account for fabrication imperfections.

FIG. 15A illustrates various schematics of metagratings including fabrication imperfections modeled by the erosion and dilation of device layouts, which is parametrized by δw. The pixelized patterns ε(x) generated from the width vector w, eroded devices may be realized by decreasing the silicon ridge widths by a width variation εw and increasing the air gap widths by Ow, and they may be represented by width vector w_e. Dilated devices may be realized by performing the opposite geometric transformations and are represented by width vector w_d. The new FoM used in the GLOnets loss function is the weighted average deflection efficiency and is defined to be:Eff_avg=0.5×Eff_original+0.25×Eff_eroded+0.25×Eff_dilated  (7)Eff_avg=0.5×Eff(w)+0.25×Eff(w_e)+0.25×Eff(w_d)  (11)

Fixed-topology reparameterized GLOnets with robustness may be used to globally optimize a metagrating with a MFS constraint of 30 nm and δw=10 nm. The device may operate with a wavelength of 850 nm and deflection angle of 65 deg, and the topology may be fixed to contain three silicon ridges.

FIG. 15B illustrates various plots comparing performance of devices globally optimized with and without robustness. The device has a minimum feature size constraint of 30 nm, wavelength of 850 nm, and a deflection angle of 65 deg, and the robust devices are designed to accommodate δw=10 nm. The deflection efficiency of the globally optimized device as a function of width variation is plotted in FIG. 15B and shows that the device operates with a relatively high efficiency above 90% for width variations within ±10 nm. In comparison, the globally optimized device without robustness has a higher overall peak efficiency when no width variation is present, but the efficiency rapidly decreases for negative width variations. Variable-topology reparameterized GLOnets with robustness may be as effective as fixed-topology reparameterized GLOnets with robustness at generating robust devices.

Given that the final pixelized pattern ε(x) may be generated from the width space {w_i}, fabrication imperfections may be defined in terms of pattern erosion (width of Si ridge decreases and width of air gap increases) and dilation (width of Si ridge increases and width of air gap decreases). In some embodiments, an averaged deflection efficiency may be defined as the figure of merit to reflect fabrication robustness where during each iteration the efficiency of multiple devices (e.g. three devices) may be averaged. The multiple devices may include:

• • The device with the original pattern,
  • The device with width of Si ridges increasing 10 nm each, and
  • The device with width of air gaps increasing 10 nm each.

In this way, robustness to fabrication imperfections may be introduced to the design space constrained with minimum feature size. FIG. 16A is a plot illustrating efficiency of a globally optimized device generated with reparameterization and accounting for fabrication imperfections at a minimum feature size of 30 nm. FIG. 16B illustrates a plan view of the optimized device without fabrication robustness. The device includes three silicon bars 702 which are interleaved with three air bars 704. This optimized device includes an absolute efficiency of 97.7%.

FIG. 16C illustrates a plan view of the optimized device with fabrication robustness and a minimum feature size of 30 nm. The device includes three silicon bars 706 which are interleaved with three air bars 708. Comparing with the non-robust device, the robust device shares a very similar pattern and maintains a high-level performance of 0.928 average efficiency within ±10 nm of pattern variation.

Optimization of Photonics Thin Film Stacks

In some embodiments, the optimization techniques may be used to optimize optical devices such as photonic thin film stacks. Thin film stacks have been used in many optical systems. FIG. 17 illustrates an example schematic of a photonic thin film stack. The materials and thickness of thin film layers may be optimized to achieve an advantageous reflection spectrum. The number of thin film layers may be fixed to be N. The refractive indices of N-layer thin film stack are denoted as a vector n(λ)=(n₁(λ), n₂(λ), . . . , n_N(λ))^T, while the thicknesses are t=(t₁, t₂, . . . , t_N)^T. Each layer may be one of M types of materials n_i(λ)∈{m₁(λ), m₂(λ), . . . , m_M(λ)}. Optimization may be performed to determine the proper n and t to achieve advantageous reflection over the wavelength bandwidth, incident angle range and specified polarization based on Equation 5:{n*,t*}=arg min_{n,t}Σ_λ,θ,pol((n,t|λ,θ,pol)−*(λ,θ,pol))²  (5)

The desired reflection spectrum is denoted as *(λ,θ,pol). The objective may be described to determine Equation 6:O(n,t)=Σ_λ,θ,pol((n,t|λ,θ,pol)−*(λ,θ,pol))²  (6)

To determine the optimization, a global optimization algorithm may be used. In some embodiments the optimization algorithm may be a transfer matrix method (TMM). FIG. 18 illustrates a flow chart of an example global optimization algorithm using a TMM in accordance with an embodiment of the invention. The first step may be to figure out the calculation of the reflection gradient with respect to refractive indices 902 and the thicknesses 904. Since the reflection spectrum of the thin film stack can be calculated using TMM, the gradient can be directly obtained by the chain rule. Therefore, a TMM solver 906 may be programmed using an auto-differentiation package such as Pytorch.

In some embodiments, a generative neural network may be used to determine a distribution of thin film stack configurations. FIG. 19 illustrates a generative neural network used to produce a distribution of thin film stack configurations. An input 1002 may be a uniformly distributed random vector z. Different z^(k) may map onto different thin film stack configurations specified as refractive indices and thicknesses, denoted as {n^(k), t^(k)}. And their reflection spectra (n^(k), t^(k)|λ, θ, pol) may then be calculated by the TMM solver 906. The global optimization problem may be reframed as optimizing the distribution of outputted configurations. The optimization objective or the loss function may be defined by Equation 7:

$\begin{array}{cc}L=𝔼\left[\exp \left(-\frac{o\left(n,t\right)}{\sigma }\right)\right]\approx \frac{1}{K}{\sum }_{k=1}^K\exp \left(-\frac{o\left(n^{\left(k\right)},t^{\left(k\right)}\right)}{\sigma }\right)& \left(7\right)\end{array}$

The gradient of the loss function with respect to the neuron weights may be calculated by propagation, and the generator may then be optimized by gradient descent. After optimization, the distribution of thin film stack configurations may converge to a delta function located at the global optimum.

In practice, the refractive indices of thin film layers may not take arbitrary values and may be chosen from a material database and thus n may be a categorical variable. However, the gradients may not be calculated for categorical variables. To overcome this issue, a new scheme to generate refractive indices is illustrated in FIG. 18. As illustrated, the network first maps the random vector z onto an N-by-M matrix A 908 and then applies a Softmax function 910 to each row of matrix A 908 to generate a probability matrix P 912. An example Softmax function 910 is represented as Equation 8:

$\begin{array}{cc}{P}_{ij}=\frac{\exp \left(\alpha ·A_{ij}\right)}{{\sum }_{j=1}^M\exp \left(\alpha ·A_{ij}\right)}& \left(8\right)\end{array}$ α is a hyperparameter to tune the sharpness of the Softmax function 910. The i^th row of probability matrix P 912 is a 1×M vector, which represents the probability distribution over the material database for the i^th layer. Then the expected refractive index of each layer given by this probability distribution is calculated by n_i(λ)=Σ_j=1^M m_j(λ)·P_ij, which is a continuous variable with the valid gradient. As such, the gradient with respect to the refractive index n_i may be backpropagated to the probability matrix P 912 and then the neurons. In addition, to ensure the material probability distribution of each layer converges to only one material, the hyperparameter a may increase over the course of optimization.

Since the desired reflection spectra usually cover a wide range of wavelengths and incident angles, finding the proper thin film configuration may be complex and challenging. A deep neural network may be beneficial to produce complex distributions in the design space in the early stage of optimization, while the ability to express a relatively simple distribution is also important for smooth convergence at the end. FIG. 19 illustrates an example deep residual network for the generator architecture. The generator may include many residual blocks 1004. The input and the output of each residual block 1004 may have the same dimension and a skip connection. The skip connection may be represented as x_out=ƒ(x_in)+x_in, where ƒ(x_in) represents the fully connected layer, batch normalization layer, and a non-linear activation layer. In some embodiments, the non-linear activation layer may be a Leaky ReLU. Skip connections may allow a deep neural network to degrade to a shallow neural network by skipping many network layers.

Examples Optimization of Photonic Thin-Film Stack

A 3-layer thin film stack forming a three-layer AR coating for Si was designed to minimize average reflection from silicon over the incident angle range between 0° and 60° and the wavelength range between 400 nm and 1100 nm. The refractive index was between 1.09 and 2.60 and the thicknesses were between 0 nm and 200 nm of each layer. FIG. 20 illustrates a cross sectional view of the resultant structure from the optimization. FIG. 21 illustrates various plots of the results for various wavelengths across angles between 0° and 60°. Other algorithms have been observed to take over 19 days of CPU computation for the optimization. A multi-start gradient optimizer may take 15 minutes to find the global optimum. The disclosed generative neural network solved for the global optimum in less than one minute with a single GPU, demonstrating its efficiency and efficacy.

FIG. 22 illustrates the results from another example optimization performed to an incandescent light bulb filter that transmits visible light and reflects infrared light. The reflected infrared light can be recycled to heat up the thermal emitter and thus the emitter is able to achieve a higher temperature while consuming a lower input power. The emission efficiency in the visible range can be enhanced in this way. The material database includes seven different dielectric material types: magnesium fluoride, silicon dioxide, silicon carbide, silicon mononitride, aluminum oxide, hafnium dioxide, and titanium dioxide. FIGS. 23 and 24 illustrate the broadband reflection characteristics of a 45-layer generative neural network optimized device. As illustrated, the device operates with nearly ideal transmission at wavelengths between 500 nm and 700 nm and nearly ideal reflection at near-infrared wavelengths, for both normal incidence and for incidence angles averaged over all possible solid angles. FIG. 25 illustrates the emission intensity spectrum of the light bulb with and without the thin film filter. To measure the enhancement of visible light emission due to the thin film filter, the emissivity enhancement factor χ is computed:

$\chi =\frac{{\int }_0^{\infty }{E}_{\mathrm{emitter}+\mathrm{stack}}\left(P_0,\lambda \right)V\left(\lambda \right)d\lambda }{{\int }_0^{\infty }{E}_{\mathrm{emitter}}\left(P_0,\lambda \right)V\left(\lambda \right)d\lambda }$

P₀ is the input power and V(λ) denotes the eye's sensitivity spectrum. The generative neural network optimized device achieved χ of 18.0 compared with a genetic algorithm optimized device which achieved χ of 14.8. The generative neural network improved the performance by 22%.

FIG. 26 is a chart comparing the results of a comparison of thin film filters generated using a residual global neural network, a deep global neural network, and a generic algorithm. The application of residual neural network to thin film filters including different number of layers and a comparison with a fully connected deep neural networks shows that residual neural network perform better in terms of finding optimized devices for non-convex problems. And for optimized devices with 45 layers, the neural network based algorithms clearly outperformed the genetic algorithm in optimization of the thin film filters. Furthermore, the neural network based algorithms optimizes devices with similar performance as devices optimized with the genetic algorithm with approximately two thirds the number of layers, which may allow for high levels of optimization with less layers thus decreasing manufacturing complexity.

Embodiments Including Optimization with Variable Topology

A global optimizer may benefit from the ability to search for the proper device topology as well as the detailed layout for that topology. Optimizers for fixed-topology devices can perform this task by parametrically sweeping across a wide range of topologies, as performed in FIG. 12B. But this route can be computationally intractable when scaling to complex systems. In some embodiments, reparameterized global optimization can combine concepts in boundary and topology optimization to globally optimize geometrically constrained devices with variable topology.

FIG. 27A illustrates a flowchart of a reparameterization process for a variable-topology global optimizer. Two separate generators 2702a, 2702b are used to produce distributions of width vectors 2704a and permittivity vectors 2704b, which are combined to produce constrained physical devices. As before, with fixed-topology optimization, the device permittivity profile is subdivided into {circumflex over (M)} sections. However, instead of assuming that a given section is either silicon or air, each section was defined to have the parameters (w_i, ), where w_i denotes the section width 2706a and {circumflex over (ε)}_l denotes the permittivity 2706b at the center of the section. A range of topologies can now arise because adjacent sections in the final binarized devices can have either the same or different permittivity value. The generalized permittivity profile 2708 of a device with the new scheme is calculated from the mapping function ε=ĥ(x, w, {circumflex over (ε)}, T), where ĥ is a piece-wise transformation.

FIG. 27B illustrates schematics of the piecewise transformation used for device reparameterization. The permittivity profile of the entire structure may be calculated by smoothly connecting the permittivity at the center of each section using a function similar to in Equation 4. It is advantageous to have the mapping function ε=ĥ(x, w, {circumflex over (ε)}, T) work in a way such that the permittivity profile is a continuous and differentiable function of x that connects the permittivity at the center of each section, as shown in FIG. 27B. The permittivity profile may be used to determine a device efficiency as described above.

FIG. 28 illustrates a piecewise transformation used in variable-topology reparameterized GLOnets. The piecewise transformation ĥ deals with each section piece (w_i, ) by separately determining the left half and the right half of the permittivity profile of the corresponding section. The left half lies from the interface between the i−1^th and the i^th section to the center of the i^th section, e.g.

$\begin{array}{cc}{\tilde{x}}^{\left(i\right)}\in \left[-\frac{w_i}{2},0\right],& \end{array}$ with its function defined as:

$$\begin{gathered} \begin{array}{cc}\hat{h}\left(x,w,\hat{\epsilon },T\right)={\hat{\epsilon }}_{i-1}+\left({\hat{\epsilon }}_i-{\hat{\epsilon }}_{i-1}\right)·\frac{h_{+}\left({\tilde{x}}^{\left(i\right)},w_i,T\right)+h_{+}\left(0,w_i,T\right)-1}{2·h_{+}\left(0,w_i,T\right)-1},\mathrm{if}{\hat{\epsilon }}_{i-1}\le {\hat{\epsilon }}_i\text{ \\ h^(x,w,εˆ,T)=εˆi+(εˆi-1-εˆi)·h-(x˜(i),wi,T)2·-(0,wi,T),if⁢εˆi-1>εˆi⁢where⁢ \\ h±(x˜(i),wi,T)=1exp(±(x∼(i))2-(wi2)2T·(Lwi⁢M)2)+1}& \end{array} \end{gathered}$$

The right half is defined in a similar way on {tilde over (x)}⁽ⁱ⁾∈

$\left[0,\frac{w_i}{2}\right],$ from the center location to the interface between the i^th and the i+1th section:

$$\begin{gathered} \begin{array}{cc}\hat{h}\left(x,w,\hat{\epsilon },T\right)={\hat{\epsilon }}_i+\left({\hat{\epsilon }}_{i+1}-{\hat{\epsilon }}_i\right)·\frac{h_{-}\left({\tilde{x}}^{\left(i\right),}{w}_i,T\right)}{2·h_{-}\left(0,w_i,T\right)},\mathrm{if}{\hat{\epsilon }}_{i-1}\le {\hat{\epsilon }}_i\text{ \\ hˆ(x,w,εˆ,T)=εˆi+(εˆi-εˆi+1)·h+(x˜(i),wi,T)+h+(0,wi,T)-12·h+(0,wi,T)-1,if⁢εˆi-1>εˆi}& \left(i\right)\end{array} \end{gathered}$$ Thus, a continuous and differentiable mapping function ĥ is defined to transform from the piecewise vectors (w,{circumflex over (ε)}) to the device's permittivity profile ε, such that the gradient can be backpropagated during the optimization process. In some embodiments, a gradient based on the width vector (w_i) and the center permittivity vector ({circumflex over (ε)}_i) may be solved for. A process such as gradient descent may be used to optimize each of these vectors in order to optimize the device. In some embodiments, these width and permittivity may be solved for simultaneously. In some embodiments, the optimum width of the elements may be solved for followed by the optimum permittivity or vice versa.

FIG. 29 illustrates an example algorithm for newly configured GLOnets where physical devices in this form are generated using two separate generative networks, one for section widths and the other for refractive indices being optimized. The widths w may be produced the same way as with fixed-topology reparameterized GLOnets, where latent vectors u generated from G_1ϕ(z) are reparametrized to w with MFS constraints. The permittivity vector {circumflex over (ε)}∈^{{circumflex over (M)}}, which determines the permittivity at the center of each section, results from latent vectors v∈^{{circumflex over (M)}} generated from G_2ψ(z), followed by the transformation {circumflex over (ε)}=q(v,γ):

$\begin{array}{cc}{\hat{\epsilon }}_i=\mathrm{sigmoid}\left(\gamma v_i\right)=\frac{e^{\gamma v_i}}{e^{\gamma v_i}+1},i=1,\dots ,\hat{M}& \left(9\right)\end{array}$ γ is a tunable hyperparameter that controls the binarization of the permittivity vector {circumflex over (ε)} and is analogous to the hyperparameter T used to control the binarization of the device widths. As with T, γ is initially set to produce grayscale refractive index values and is manually increased in a gradual manner so that the final device possesses binary refractive index values of silicon or air. {circumflex over (M)} specifies the total number of sections and sets the upper limit in the topological complexity of generated devices. For example, if {circumflex over (M)}=10, physical devices with up to five nanoridges can be generated. In some embodiments, the physical device may include less than

$\frac{\hat{M}}{2}$ nanoridges. For example, it {circumflex over (M)}=10, the physical device may have 4 nanoridges.

The ability for variable-topology reparameterized GLOnets to search across different device topologies is demonstrated for metagratings with the same operating parameters specified in FIGS. 7A and 10. FIG. 27C illustrates various device layouts over the course of GLOnets training for a representative optimization run. The devices are initially grayscale with varying topology and they converge to a narrow distribution of binary devices. The network is trained using a total of 100 iterations with a batch size of 100, and device layouts at different stages of optimization are shown in FIG. 27C. Near the beginning of the optimization at iteration no. 5, GLOnets generate a wide range of random devices with different topologies. As the training process evolves, the device distribution gradually converges to a narrow range of topologies, and by iteration no. 50, G_2Ψ(z) has collapsed onto a single device topology. GLOnets subsequently focus only on boundary optimization to refine the four-ridge device layout. At the end of training, the generated devices are clustered around the same global optimum as that achieved with Algorithm 2 of FIG. 9 for M=8.

FIG. 27D illustrates a plot of device efficiencies as a function of the minimum feature size (MFS) constraint, designed using fixed-topology and variable-topology reparameterized GLOnets. To perform a more systematic benchmark analysis of the variable-topology reparameterized GLOnets, metagratings with parameters matching those in FIG. 12B were designed for a wide range of MFS constraints. The parameters match those in FIG. 7A. The results are plotted in FIG. 27D together with the efficiencies of the best overall devices from FIG. 12B, and they show that variable-topology reparameterized GLOnets are able to search for the correct device topology for a given MFS constraint. Furthermore, the best variable-topology devices (solid black line) largely follow the performance of the best fixed-topology devices (dashed color lines), indicating that this new variant of GLOnets can effectively perform reparameterized boundary optimization.

Variable-topology reparameterized GLOnets were benchmarked with two other methods. The first is local adjoint based optimization, where 200 unconstrained topology optimization iterations are first performed to produce binary devices from random grayscale devices, followed by 50 boundary optimization iterations to refine the binary device layouts and incorporate MFS constraints. The second is the original GLOnets, where a Gaussian filter is used at the network output in an attempt to impose a MFS constraint. The device designs are benchmarked by two criteria, deflection efficiency and the MFS in the physical device.

FIG. 27E illustrates scatter plots of device efficiency and MFS for devices designed using variable-topology reparameterized GLOnets, original GLOnets, and adjoint-based local optimization, for three different MFS constraints. The scatter plots illustrate the efficiencies and MFSs of devices designed using the three aforementioned techniques, for three MFS constraints. Variable-topology reparameterized GLOnets generate devices that are clustered around high efficiency values and that strictly satisfy the MFS constraint. The original GLOnets are also able to find clusters of devices with high deflection efficiency. However, these devices do not consistently satisfy the hard MFS constraints because the Gaussian filter is a soft constraint. The combination of local topology and boundary optimization may be able to produce devices that strictly satisfy the MFS constraint, but may not guarantee high performance. For small MFS constraints, the efficiencies of the best locally optimized devices are comparable with those of the best reparameterized GLOnets devices, as the globally optimal unconstrained and constrained devices are in similar design space regions. As such, enough instances of unconstrained local topology optimization may produce some devices near the globally optimal constrained device. However, for large MFS constraints, the local optimization approach may have difficulty in reliably finding the global optimum. The reason is because the topology of devices with large MFS constraints can be very different from the unconstrained devices (see FIG. 12B), so that performing unconstrained topology optimization followed by local boundary optimization may not be effective.

Embodiments Including all Binary Devices

In the version of reparameterization discussed above in connection with FIGS. 3-5, the pattern obtained from the latent vector may gradually binarize over the iterative optimization process depending on the hyperparameter T which controls the degree of binarization (see Equation 4). In this manner, only the device from the very last step (where T=1) after completion of the optimization process may be completely binary, while devices from all the intermediate steps are represented in grayscale. This gradual tuning of the grayscale level can post several challenges to the inverse design process: 1) the hyperparameter T may benefit from being tuned in a judicious manner such that it is neither too fast nor too slow to guarantee a stable and smooth transition between grayscale and binary structures. 2) Since all the intermediate devices are not binary, only the device from the very last step may be used based on the design constraints on material choices due to fabrication concerns (for example, in the case when the device is constrained between air and silicon).

In some embodiment, the binary devices may be obtained and evaluated at each step during the iterative optimization. Unlike the reparameterization process described in connection with FIGS. 3-5 where the grayscale pattern is used to directly solve for the efficiency (Eff) and efficiency-permittivity gradient

$\left(\frac{\partial \mathrm{Eff}}{\partial \epsilon }\right),$ the grayscale pattern may be binarized. FIG. 30 illustrates an example process including a binarization step in accordance with an embodiment of the invention. Similar to the reparameterization process of FIGS. 3-5, first a set of latent variables 3002 is reparametrized into a grayscale pattern (ε) 3004. The grayscale pattern (ε) 3004 may be binarized into a binary pattern (ε_b) 3006 by thresholding the dielectric constant. Then, an efficiency (Eff_b) 3008 of the binary pattern and gradient

$\left(\frac{\partial {\mathrm{Eff}}_b}{\partial {\epsilon }_b}\right)$ of the binary pattern may be simulated with an EM solver. Lastly, the gradient evaluated using the binary device may be directly applied onto the grayscale pattern to update the grayscale pattern, and eventually backpropagate onto the latent vector. The gradient may be used to optimize the latent vector through processes such as gradient descent. In some embodiments, the gradient on the binary device

$\left(\frac{\partial {\mathrm{Eff}}_b}{\partial {\epsilon }_b}\right)$ may be directly applied to the grayscale device because the grayscale device is close enough to the binary device under a small T. In this manner, the EM simulation at each step may be performed for the binary device. The hyperparameter T can be set as a constant value throughout the optimization as long as the binary device is similar enough to its corresponding grayscale version.

FIGS. 31A and 31B illustrate a comparison between the reparameterization process described in connection with FIGS. 3-5 and the reparameterization process utilizing all binary devices of FIG. 30. Here a 65 degree silicon metagrating that works at 850 nm wavelength with a minimum feature size constraint of 30 nm was optimized using GLOnets. FIG. 31A illustrates results from the reparameterization process described in connection with FIGS. 3-5 whereas FIG. 31B illustrates results from the reparameterization process including the binarization step described in connection with FIG. 30. These methods were able to obtain a similar performing output device efficiency of 96% after a similar number of iterations. Unlike the original version where T gradually decreases, the fully binary version has a constant T, which can be seen from the flat binarization (green) curve on the right panel. In this example, the binarization hyperparameter was set at 0.95. Moreover, since devices from each iteration are all binary, the optimization routine may be paused at any step with a fully binary device obtained at any step, for example as earlier as the 50 iteration step where the best device curve (purple) becomes relatively flat, the iteration may be paused.

Example Network System

FIG. 32 illustrates an example of optimization system that optimizes photonic devices while enforcing fabrication constraints in accordance with an embodiment of the invention. Network 3200 includes a communications network 3260. The communications network 3260 is a network such as the Internet that allows devices connected to the network 3260 to communicate with other connected devices. Server systems 3210, 3240, and 3270 are connected to the network 3260. Each of the server systems 3210, 3240, and 3270 is a group of one or more servers communicatively connected to one another via internal networks that execute processes that provide cloud services to users over the network 3260. One skilled in the art will recognize that the optimization system may exclude certain components and/or include other components that are omitted for brevity without departing from this invention.

For purposes of this discussion, cloud services are one or more applications that are executed by one or more server systems to provide data and/or executable applications to devices over a network. The server systems 3210, 3240, and 3270 are shown each having three servers in the internal network. However, the server systems 3210, 3240 and 3270 may include any number of servers and any additional number of server systems may be connected to the network 3260 to provide cloud services. In accordance with various embodiments of this invention, the optimization system that uses systems and methods that optimizes photonic devices while enforcing fabrication constraints in accordance with an embodiment of the invention may be provided by a process being executed on a single server system and/or a group of server systems communicating over network 3260.

Users may use personal devices 3280 and 3220 that connect to the network 3260 to perform processes that optimizes photonic devices while enforcing fabrication constraints in accordance with various embodiments of the invention. In the shown embodiment, the personal devices 3280 are shown as desktop computers that are connected via a conventional “wired” connection to the network 3260. However, the personal device 3280 may be a desktop computer, a laptop computer, a smart television, an entertainment gaming console, or any other device that connects to the network 3260 via a “wired” connection. The mobile device 3220 connects to network 3260 using a wireless connection. A wireless connection is a connection that uses Radio Frequency (RF) signals, Infrared signals, or any other form of wireless signaling to connect to the network 3260. In the example of this figure, the mobile device 3220 is a mobile telephone. However, mobile device 3220 may be a mobile phone, Personal Digital Assistant (PDA), a tablet, a smartphone, or any other type of device that connects to network 3260 via wireless connection without departing from this invention.

As can readily be appreciated the specific computing system used to optimizes photonic devices while enforcing fabrication constraints is largely dependent upon the requirements of a given application and should not be considered as limited to any specific computing system(s) implementation.

Example Optimizer

FIG. 33 illustrates an example of an optimizer element that executes instructions to perform processes that optimizes photonic devices while enforcing fabrication constraints in accordance with an embodiment of the invention. Optimizer elements in accordance with many embodiments of the invention can include (but are not limited to) one or more of mobile devices, cameras, and/or computers. Optimizer element 3300 includes processor 3305, peripherals 3310, network interface 3315, and memory 3320. One skilled in the art will recognize that the optimizer element may exclude certain components and/or include other components that are omitted for brevity without departing from this invention.

The processor 3305 can include (but is not limited to) a processor, microprocessor, controller, or a combination of processors, microprocessor, and/or controllers that performs instructions stored in the memory 3320 to manipulate data stored in the memory. Processor instructions can configure the processor 3305 to perform processes in accordance with certain embodiments of the invention.

Peripherals 3310 can include any of a variety of components for capturing data, such as (but not limited to) cameras, displays, and/or sensors. In a variety of embodiments, peripherals can be used to gather inputs and/or provide outputs. The optimizer element 3300 can utilize network interface 3315 to transmit and receive data over a network based upon the instructions performed by processor 3305. Peripherals and/or network interfaces in accordance with many embodiments of the invention can be used to gather inputs that can be used to optimize photonic devices while enforcing fabrication constraints.

Memory 3320 includes an optimization application 3325 and generator 3330. The optimization application 3325 and the generator 3330 in accordance with several embodiments of the invention can be used to optimize photonic devices while enforcing fabrication constraints.

The optimization application 3325 and generator 3330 may be used to perform the methods described above and the methods described below in FIG. 33. For example, the optimization application 3325 may receive a set of unconstrained latent variables; map the set of unconstrained latent variables to a constrained space to generate a constrained device; calculate the permittivity across each element of the constrained device; determine a permittivity-constrained width gradient based at least partially on the permittivity across each element; and optimize the set of unconstrained latent variables by at least partially using the permittivity-constrained width gradient.

In some embodiments, mapping the unconstrained space to a constrained space may include reparameterization of the unconstrained latent variables.

In some embodiments, calculating the permittivity across each element of the constrained device may include performing a gray-space relaxation of the set of constrained width space to determine a function of refractive index across each element in the set of constrained width space.

In some embodiments, the optimization application 3325 may further determine an efficiency-constrained width gradient of the constrained device based on the permittivity-constrained width gradient. The efficiency-constrained width gradient is further based on an efficiency-permittivity efficiency gradient.

In some embodiments, the generator 3330 may generate the set of unconstrained latent variables; and the generator 3330 may be iteratively updated at least partially based on the permittivity-constrained width gradient.

In various embodiments, the optimization application 3325 may receive a set of variables including a set of refractive indices including a set of corresponding thicknesses; calculate a reflection gradient with respect to the set of refractive indices and their set of corresponding thicknesses; optimize the set of variables at least partially using the reflection gradient to obtain an optimized device including a set of optimized refractive indices with corresponding optimized thicknesses.

In some embodiments, the optimization application 3325 may further calculate the reflection spectrum of the set of variables, where calculating the reflection gradient comprises calculating a partial differential of the reflection spectrum. Calculating the reflection spectrum may include using a transfer matrix method.

In some embodiments, the generator 3330 may be used to generate the set of variables; the set of variables may be used to calculate a corresponding reflection gradient; and the generator 3330 may be iteratively update at least partially using the reflection gradient.

Although a specific example of element 3300 is illustrated in this figure, any of a variety of optimizer elements can be utilized to perform processes optimize photonic devices while enforcing fabrication constraints similar to those described herein as appropriate to the requirements of specific applications in accordance with embodiments of the invention.

Example Optimization Method

FIG. 34 illustrates an example optimization method 3400 in accordance with embodiments of the invention. The optimization method 3400 includes receiving (3402) a set of unconstrained latent variables. A generator may be used to generate the set of unconstrained latent variables. The generator may be updated at least partially based on a permittivity-constrained width gradient. The set of unconstrained latent variables may be initially randomly generated.

The optimization method 3400 may further include mapping (3404) the set of unconstrained latent variables to a constrained space to generate a constrained device. The mapping may include reparameterization of the unconstrained latent variables. the reparameterization of the unconstrained latent variables may include solving for a set of constrained width space {w_i} through the following equations:

$v_i=\mathrm{Sigmoid}\left(u_i\right),v_M=1;k_i=v_i\sqrt{1-{\sum }_{j=1}^{i-1}{k}_j^2},k_1=v_1;\mathrm{and}{w}_i=\frac{k_i}{{\sum }_{j=1}^M{k}_j}\left(L-{\mathrm{Mw}}_{\min }\right)+w_{\min },$ where the set of unconstrained latent variables are {u_i} and wherein w_min is a minimum width of each feature of the constrained width space, M is the total number of features, and Σ_i=1^M w_i=L.

The optimization method 3400 may further include calculating (3406) the permittivity across each element of the constrained device. Calculating the permittivity across each element of the constrained device may include performing a gray-space relaxation of the set of constrained width space to determine a function of refractive index across each element in the set of constrained width space. The gray-space relaxation may include solving for:

$n_i\left(x\right)={\left[\exp \left(\pm \frac{x^2-{\left(\frac{w_i}{2}\right)}^2}{T}\right)+1\right]}^{-1},$ where x is the position across the width of a certain element, where x=0 is at the middle of the element, and T controls the binarization of the pattern.

The optimization method 3400 may further include determining (3408) a permittivity-constrained width gradient based at least partially on the permittivity across each element. The optimization method 3400 may further include optimizing (3410) the set of unconstrained latent variables at least partially using the permittivity-constrained width gradient. The permittivity-constrained width gradient may be at least partially used to calculate the efficiency-constrained width gradient.

In embodiments where a generator is used to generate the set of unconstrained latent variables, the generator may be iteratively updated at least partially based on the permittivity-constrained width gradient. The generator may be used to generate another set of unconstrained latent variables which may be used to update the constrained device.

DOCTRINE OF EQUIVALENTS

While the above description contains many specific embodiments of the invention, these should not be construed as limitations on the scope of the invention, but rather as an example of one embodiment thereof. It is therefore to be understood that the present invention may be practiced in ways other than specifically described, without departing from the scope and spirit of the present invention. Thus, embodiments of the present invention should be considered in all respects as illustrative and not restrictive. Accordingly, the scope of the invention should be determined not by the embodiments illustrated, but by the appended claims and their equivalents.

Claims

1. A method of manufacturing photonic devices, the method comprising: using a generator to generate a set of unconstrained latent variables;mapping the set of unconstrained latent variables to a constrained space to generate a constrained photonic device;calculating the permittivity across each element of the constrained photonic device;determining a permittivity-constrained variable gradient based at least partially on the permittivity across each element of the constrained photonic device;iteratively updating the generator at least partially based on the permittivity-constrained variable gradient;generating using the generator an optimized set of unconstrained latent variables; andmanufacturing the generated constrained photonic device.

2. The method of claim 1, wherein mapping the unconstrained space to a constrained space comprises reparameterization of the unconstrained latent variables.

3. The method of claim 2, wherein the reparameterization of the unconstrained latent variables comprises solving for a set of constrained width space {wi} through the following equations:

4. The method of claim 3 wherein calculating the permittivity across each element of the constrained device comprises performing a gray-space relaxation of the set of constrained width space to determine a function of refractive index across each element in the set of constrained width space.

5. The method of claim 4, wherein the gray-space relaxation comprises solving for

6. The method of claim 4, further comprising determining an efficiency-constrained variable gradient of the constrained device based on the permittivity-constrained variable gradient.

7. The method of claim 6, wherein the efficiency-constrained variable gradient is further based on an efficiency-permittivity gradient.

8. The method of claim 7, wherein the efficiency-constrained variable gradient is calculated using the following equation:

9. The method of claim 8, wherein an efficiency-latent variable gradient is calculated at least partially using a latent variable gradient of the unconstrained latent variables in relation to constrained device calculated using the following equations:

10. The method of claim 9, wherein the efficiency grating in relation to the unconstrained latent variables is calculated at least partially using the following equation:

11. The method of claim 10, wherein the efficiency-permittivity gradient is calculated using an adjoint variable method.

12. The method of claim 11, wherein the adjoint variable method comprises forward and adjoint simulation.

13. The method of claim 1, wherein the set of unconstrained latent variables comprise unconstrained latent widths, and wherein mapping the set of unconstrained latent variables comprises mapping the unconstrained latent widths to the constrained space to generate constrained widths, and the method further comprises: receiving a set of unconstrained latent permittivities; andmapping the set of unconstrained latent permittivities to a constrained space to generate constrained average permittivities of the constrained device,wherein calculating the permittivity profile across each element of the constrained device is further based on the constrained widths and the constrained average permittivities of the constrained device.

14. The method of claim 13, further comprising: determining a permittivity-constrained average permittivity gradient of the constrained device, andoptimizing the set of unconstrained latent permittivities at least partially based on the permittivity-constrained average permittivity gradient of the constrained device.

15. The method of claim 1, further comprising binarizing the permittivity across each element of the constrained device, wherein determining a permittivity-constrained variable gradient is based on the binarized permittivity across each element.

16. A method for manufacturing photonic devices, the method comprising: using a generator to generate a set of variables including a set of refractive indices and a set of thicknesses, wherein each of the refractive indices of the set of refractive indices corresponds to a thickness of the set of thicknesses;calculating a reflection gradient with respect to the set of refractive indices and the set of thicknesses;optimizing the set of variables at least partially using the reflection gradient to obtain an optimized device including a set of optimized refractive indices with corresponding optimized thicknesses;iteratively updating the generator at least partially using the reflection gradient; andmanufacturing the optimized device including the set of optimized refractive indices with corresponding optimized thicknesses.

17. The method of claim 16, further comprising calculating the reflection spectrum of the set of variables, wherein calculating the reflection gradient comprises calculating a partial differential of the reflection spectrum.

18. The method of claim 17, wherein calculating the reflection spectrum comprises using a transfer matrix method.