Establishing exhaustive metasurface robustness against fabrication uncertainties through deep learning

2.1 Metasurface parameterization

Here, we step away from TO in order to address the question of how to make an exhaustive guarantee of metasurface robustness. Because our designs are not parameterized per pixel, we exchange the TO limitation in exhaustive robustness for a limitation in structural flexibility.

Whereas a traditional TO supercell has the potential to produce entirely contiguous structures, we chose to divide our supercell into discrete unit cells (although our method is flexible enough to handle contiguous supercells). This partitioning helps to preserve the analogy with the more traditional library-sourced supercells in the literature, as developing a library of individual unit cells is necessary to scale metasurfaces up for applications such as metalenses [42], [43], [44]. Traditional metasurface designs select an ordering of unit cells with different phase responses to tile into the metasurface. This causes a transmitted wave to be deflected anomalously into a chosen diffraction order m. A metasurface can be made to focus by engineering a phase gradient that changes linearly from its center to the exterior. The core of the focusing metasurface is the development of supercells. Each supercell is composed of several unit cells, with the combination of unit cells providing the phase gradient needed. The diffraction angle is chosen with Equation (1):

where d is the supercell width. This applies to all diffraction orders m where θ ∈ ℝ.

In [45], the authors demonstrated how supercell construction is complicated by coupling effects between adjacent unit cells which interfere with the supercell’s diffraction efficiency. Indeed, constructing a supercell with maximized efficiency from individually optimized unit cells is often nontrivial. In the case of this study, such an approach could also serve to obscure the source of losses between fabrication error and coupling effects. As our goal is to clearly identify losses in efficiency that are directly due to fabrication error, we simulate and optimize full supercells rather than individual unit cells.

For this work, we chose a freeform geometrical parameterization for our unit cells which has been used successfully in [46]. As shown in Figure 1, unit cell masks are created using a fourfold symmetric spline surface function defined in terms of several control points.

Figure 1:

Process for creating unit masks using thresholding of spline surfaces, buffering, and MFS filtering. Each unit cell is assumed to be fourfold symmetric. After constructing these masks, the unit cells are combined together in sets of four to create supercells.

The spline surface function interpolates between the control points, producing a smooth freeform surface that can be manipulated by varying the control point heights. In this work, each unit cell is defined by a 3 × 3 uniformly spaced grid of control points in one quadrant which is then translated to the other quadrants to enforce fourfold symmetry (although twofold symmetry and completely arbitrary geometries are also trivial to generate with this approach). To produce a mask, the function is rasterized and then a thresholding value is used to create the final binary mask. Each unit cell is given a buffer of pixels around its exterior to enforce discontinuity between the unit cells. To ensure that unit cell structures have a realistic MFS, each cell is further refined using a series of image processing techniques that enforce the desired MFS. This process is performed for each unit cell before they are combined together to form the final four-element supercell. The spline-based structural parameterization used in this study serves as a proof-of-concept, but the techniques that follow can be extended to individual unit cells, contiguous supercells, and metadevices based on true 3D geometries [47].

For a chosen operational wavelength of λ = 1.55 µm, the supercell width was selected to be 4λ/3 (2067 nm). Unit cells had equal side lengths of λ/3 (516.7 nm). Silicon was chosen for the patterned layer and it is assumed to be on top of an infinite SiO₂ substrate. The Si pattern height was chosen to be λ/2 (775 nm). Transmission through the supercell was assumed in the +z direction from the substrate into the air. This arrangement produces three possible diffraction orders, −1, 0, and +1 at angles of −48.6°, 0°, and 48.6°, respectively. Each supercell is evaluated using a rigorous coupled-wave analysis (RCWA) code, and diffraction orders were calculated based on the fields exiting the structure [48]. The structures were rasterized at a resolution of 256 px by 1024 px, leading to a pixel side length of approximately 2 nm.

2.2 Tolerance analysis

There are many ways to model over/under dosing/etching numerically, ranging from highly physics informed to simpler image processing methods [9, 10, 49, 50]. Nevertheless, they all produce similar geometrical erosion or dilation of a structure’s mask, also known as edge deviation.

Certain types of designs are less prone to sensitivity than others. However, as will be shown, even designs with relatively relaxed goals can be overly sensitive despite being optimized for robustness in a nonexhaustive manner. Capturing edge deviation with a granularity relevant to a thorough sensitivity analysis produces a demanding resolution requirement. Whereas previous work has carried out inverse-design using only nominal and extreme erosion/dilations, we found that relevant fluctuations in the performance can occur due to very small edge deviation increments between the extrema. Even at the maximum resolution feasible due to our hardware constraints, pixel-by-pixel changes (±0.0013λ or ±2 nm for our problem) could still produce noticeable changes in a design’s performance. Combining this level of detail with a wide range of target edge deviations can cause a dramatic increase in the number of simulations required to test a given design. Beyond the nominal performance and performances at the extreme positive or negative edge deviations, an exhaustive approach must evaluate a significant number of designs in-between.

For a proof-of-concept demonstration, we implemented binary mask erosion and dilation using standard image processing algorithms. For n pixels of erosion or dilation, the masks were iteratively modified using a single-pixel square structuring element n times. Since our goal was to measure changes on a per-pixel level, this method made it easy to track edge deviation length.

2.3 Deep learning

To properly stress test the proposed process, we resolved to test a very wide range of possible edge deviations, out to ±20 nm (±0.013λ). At our resolution, this corresponds to a total of 21 variations of the base design which must be simulated per supercell. Despite using a fast RCWA solver to evaluate designs and their variants, the resolution and variation count requirements of this problem present a significant computational challenge, especially when contextualized inside an optimization. A consequence of more complex shape parameterizations than canonical shapes such as rectangles and disks is that, in addition to the improved performance opportunity [46], they are simultaneously more challenging to optimize due to the increased degrees of freedom. This supercell configuration specifically has 36 degrees of freedom (nine control points per unit cell) which dictates the supercell mask geometry. We found that the optimizer NSGA-II converged after more than 150,000 samples, a substantial number. Indeed, this may be the lower end of how many function evaluations may actually be needed to exhaustively optimize a problem in 36 dimensions such as the aforementioned.

On the one hand, the scale of a supercell optimization where each sample requires 21 different simulations anticipates a total number of simulations in excess of 3,000,000; a number that would be utterly intractable with commercial finite element (FEM) or finite difference time domain (FDTD) solvers. However, the total number of samples needed to train a neural network to approximate a supercell simulation, even with edge deviations, is only a fraction of that amount; on the order of ∼65,000 simulations. This is the core motivation for introducing DL into this problem.

As has already been discussed, there are many different potential DL models available for estimating electromagnetic phenomena, and even metasurface behavior specifically. This work employs a model similar to others found in the literature and relies specifically on a combination of the U-Net and convolutional neural network (CNN) topologies [51, 52]. Figure 2 shows the full structure and details of the model, with subnetworks labeled.

Figure 2:

Diagram of network topology, composed of two subnetworks. The first subnetwork (top) converts a binary supercell mask to electric fields via a U-Net. Crosslinks from the U-net help ensure that the geometry of the structure is preserved in the E-fields. The second subnetwork (bottom) converts electric fields to diffraction coefficients using a standard multi-layered CNN. Both subnetworks exploit mirror/anti-mirror y-symmetry of the mask and fields to cut the size of the network by half. Each layer or set of layers is labeled to show its shape: y-pixels by x-pixels by # of features.

As shown in the figure, the first half of the network (Si to E-field) converts a supercell binary mask to E-fields within the Si pattern and is based on the U-Net architecture. The U-Net architecture shares some structural similarities with an autoencoder in that it compresses an input image down to a latent space (or set of features) before decompressing it back to an image. However, a key innovation of the U-Net is the introduction of crosslinks which copy values from the compression half of the U-Net to the decompression half. Developed for image segmentation and tagging purposes, the U-Net is useful for preserving sharp geometrical features despite the compression process, which makes it ideal for converting a binary mask to E-fields.

The second half of the network (E-fields to Diffraction Coefficients) takes the interior E-fields and predicts diffraction coefficients for the transmitted fields. This subnetwork is a standard CNN that uses a series of convolution layers to extract features from the image before performing a regression on the resulting feature set.

A unique attribute of this work is the use of an intermediate E-field representation that precedes the final prediction of diffraction coefficients. This structure was selected after experimentation with a wide variety of alternate structures due to its comparatively low error. It is worth noting that a very low error in the diffraction coefficient is necessary as computing the efficiency will tend to exaggerate any discrepancies due to its squaring of the coefficient’s magnitude.

3.1 Data and training

Successfully training a DNN requires the selection of the right terms, which dictate the learning process, known as hyperparameters. While there are some general guidelines available for how to choose these (i.e., learning rate, minibatch size, dataset size, etc.), ultimately the right choice can vary from problem to problem.

For image-based computational electromagnetics DL models, recent work has typically used training sets on the order of 20,000 to 50,000 samples [18], [19], [20, 22, 23, 25]. For our problem, we used RCWA to sample a total of 81,408 random supercell designs. 65,024 of these samples were used for training while 16,384 were reserved for validation; an 80:20 training to validation set ratio.

Beyond the 36 variables defining the supercell mask control points, an additional variable representing the erosion/dilation was randomly sampled as well. Each data point was a set of three terms: a mask, E-fields (real and imaginary), and diffraction coefficients (real and imaginary). Because the network is constructed as a pair, both subnetworks needed to be trained separately. The hyperparameters governing the learning can be seen in Table 1. As is typical with DL models, their applicability from one problem to the next is highly sensitive to the proper selection of hyperparameters. These were selected through a process of trial-and-error, a reality of DL development that is itself a point of much research.

Some resources provided critical guidance in the selection of our hyperparameters. For example, for some time now it has been known that DL models can generalize better when trained with smaller minibatch sizes [53]. Owing to the high-resolution requirements of fine-grained tolerancing, the pair of models took up substantial memory during training. Thus, we used a relatively small minibatch size of 16 designs for both the Si to E-field and E-field to Diffraction Coefficient networks. Another way GPU memory overuse was ameliorated was to exploit mirror symmetry in the y-plane of the supercell. This allowed us to cut the data size in half and was achieved by introducing special boundary conditions in the convolution layers which mixed circular and reflection padding.

Figure 3A and B show the training progression of the two networks for both the training and validation datasets. In both cases, the learning rate was initialized to 1e−4 and then made to decay with a scalar of 0.9925 at an interval of 200 iterations. This decay helped both the training and validation mean-squared errors (MSEs) to converge to smaller values at the conclusion of training. Before training, all data were normalized to have a mean of 0 and a variance of 1. The model was developed and trained using the PyTorch library on a machine with four Nvidia RTX 2080 ti graphics cards, each providing 11 GB of RAM for a total of 44 GB of GPU RAM.

Figure 3:

MSE curves show the training of both halves of the network.

(A) Si to E-field U-Net training convergence curves. (B) E-field to diffraction efficiency CNN training convergence curves. MSE for training and validation sets are shown in both figures. Training curves have been smoothed across multiple minibatches to emphasize convergence trends. (C) Normalized cross-correlation of the maximum E-field predicted by the Si to E-field network for the validation dataset, with a median value of 0.993 shown. (D) Percent errors for both the Si to E-field network and E-field to diffraction efficiency network tested with the validation set. For the Si to E-field network, the relative error is reported between prediction and ground truth max E-field yielding a median value of 6.08%. Diffraction efficiency errors are also reported but as an absolute error. The E-field to diffraction efficiency network alone has a median absolute percent error of 0.76%, while the full network chain has a median absolute percent error of 2.45%.

To evaluate the accuracy of the networks, we employed several error metrics. Additionally, the networks had to be tested individually and in series. For the first network which outputs E-fields, we evaluated two error metrics for each design in the validation set: normalized cross-correlation (NCC) and relative percent error. Both metrics were evaluated by comparing the max E-field within the structures. Figure 3C shows the NCC of the validation set E-fields. Since an ideal value of NCC is 1 and any value over 0.8 is considered good, a median NCC of 0.993 indicates that the network is a good predictor. Because the interior E-fields are unbounded, a relative error metric is more helpful in understanding the network’s accuracy as opposed to an absolute error metric. Figure 3D shows the validation set’s median relative percent error of 6.08%, which is in good agreement with similar networks reported elsewhere in the literature [26]. For the second network, which outputs diffraction efficiencies, we compute an absolute percent error for each design. Figure 3D shows this metric for both the second network by itself and when put in series with the first network. These arrangements both indicate very good agreement, achieving median absolute errors of 0.76% and 2.45%, respectively.

As will be discussed in the next section on optimization, training our network to have a low absolute DE error was crucial to making this solution work. Even under simpler conditions, applying a state-of-the-art evolutionary optimizer to a DNN can often be met with challenges. However, our task was made significantly more difficult than is typical in nanophotonic DL because the high-resolution requirement (necessitated by tolerance analysis) scaled up the network size dramatically. This made training slower and more difficult. Nevertheless, our training regimen was successful in producing a suitable network for use in optimization.

3.2 Optimization and filtering

With the networks trained, we approached the inverse design of exhaustively robust supercells using a multiobjective evolutionary algorithm (MOEA). MOEAs are a category of optimizers designed to characterize the tradeoff between two or more objective functions at the same time [54]. Rather than identifying a single solution as is typical for a single objective optimizer, the goal of a multiobjective optimizer is to produce a family of solutions, called a Pareto set, that represent the best tradeoff between the objectives. Like their single objective counterparts, MOEAs are typically global optimizers. This study relies on this global characteristic to overcome the inherent multimodality which arises when using more complex supercell parameterizations and which challenges TO-based inverse-design. For this work, the NSGA-II optimizer was selected due to its popularity and proven applicability to a wide range of problems [54].

For this study, two cost functions were used to explore the space of possible supercells, summarized in Equation (2). First, the nominal diffraction efficiency DE₀ is used to search for designs that are performant in the base case. Second, a worst-case change is computed by taking the difference of the diffraction efficiencies DE _i of all possible variations of the structure with its nominal diffraction efficiency. By optimizing against these two costs, we can clearly show a tradeoff between nominal performance and guaranteed performance under edge deviation.

However, combining optimization with DL is not without its challenges. Because DL models can only be used to produce good approximations of an objective function, there are often multiple subdomains of inputs to the neural network which will produce artificially attractive cost values. Therefore, pairing neural networks with optimization can regularly lead to situations where inaccuracies of the network are exploited by the optimizer. Fortunately, there are different methods to overcome this problem. First, neural networks can be integrated into MOEAs to form a surrogate-assisted optimization. This has the advantage of adapting the neural network to better fit the objective function where it matters. Another option—the one chosen for this work—is to filter the cost values as they are computed. By comparing the network’s prediction with the full-wave solution at just the nominal case, designs which have inherently high error can be filtered from the set. Because there are so many more variations to be tested than just the nominal performance for any given supercell, there is still a substantial speedup to be had from this approach. If a design’s nominal performance exceeds some maximum error value (2% in our case), then that cost can be modified according to Equation (3):

This method exposes and disincentivizes those input subdomains of the objective function which have large MSE in a way that will guide the optimizer away from them and toward well-fit solutions. It is of note that this issue can arise even with networks that have relatively high accuracy (e.g., <5% error in absolute DE). Furthermore, by testing nominal designs and culling those with high error, it is possible to short-circuit the exhaustive test with the DNN on such samples, which can afford further time savings.

3.3 Edge deviation study

Multiobjective optimization of the spline-cut supercells for nominal diffraction efficiency and stability reveals very important robustness characteristics of optimized supercells. Figure 4 summarizes the differences that can result from selecting supercell designs with or without an exhaustive performance guarantee. Within the Pareto front, there is a clear trend where the highest efficiency designs are unstable against perturbations, but lower efficiency designs can be more stable. While this difference would be expected to a certain extent, we can now explicitly quantify this tradeoff. Indeed, as shown in Figure 4E, the improvement in guaranteed performance can be significant, on the order of >35% absolute diffraction efficiency. The degree of the improvement depends on the extent of perturbation considered.

Figure 4:

Robustness tradeoff analysis.

(A) Multiobjective study of supercell nominal performance versus loss under ±20 nm of perturbation. The nominal performance of each design is the y-axis, and the guaranteed performance is the x-axis. The Pareto front is marked with a dashed green line—an ideal design would lie in the upper-right corner of the plot. Three designs marked a, b, and c show the best selection against different metrics: Best nominal, best weighted, and best guarantee + nominal, respectively. (B) Replotting the samples from (A) with a new y-axis now showing the weighted metric. (C) Comparison of edge deviation curves for options a, b, and c, including prediction made by DNN and validation from RCWA solver. Paying a loss of 4.5% in absolute diffraction efficiency for the nominal case can lead to an improvement of 35% absolute diffraction efficiency in the guarantee. (D) Supercell masks for the three designs shown. (E) Table comparing the performance difference between the three designs in terms of nominal +1 DE and guaranteed +1 DE. Design found using the neural-network-predicted guarantee (c) achieves a substantial boost in guaranteed performance over a traditionally optimized design (a).

In Figure 4A, all the samples from optimization are plotted with respect to their Nominal and Guaranteed +1 DE, as predicted by the DNN. The Pareto front, which is the optimal tradeoff between these two objectives, is marked with a dashed green line. The extrema for each objective, marked a and c, demonstrate how much difference there is between a design selected for performance-only and robustness, respectively.

Also included is point b, selected as the maximum of a different objective. Figure 4B recontextualizes the samples from the optimization on a new y-axis which shows a weighted sum of the nominal +1 DE with the extreme perturbation (±20 nm) +1 DEs. While design b does ultimately achieve a higher guaranteed +1 DE than design a, it is clear that there are even more robust solutions like design c.

Figure 4C is an edge deviation plot that shows how these three designs compare to each other, as predicted by the DNN and validated with RCWA. For example, with just ±5 nm of perturbation, design a’s performance degrades to such an extent as to make it interchangeable with designs b or c in terms of guaranteed performance. Indeed, the full curve showing the falloff trends for the two designs indicates just how much better design c does across the board, culminating in a >35% absolute performance gain over at 10 nm of perturbation. All that for <4.5% loss in nominal performance makes design c a potential preference in terms of overall performance.

This study affirms what has been shown in other studies, that a naïve single objective optimization of nominal performance will not only yield unstable results but significantly worse guaranteed performance as well. However, this work expands upon this result by showing, for the first time, that establishing exhaustive metasurface robustness requires fine-grained sampling of variations. Considering a weighted metric of the nominal case and extreme perturbations provides some advantage over purely nominal-performance focused optimization. However, failures in the edge deviation curve can occur at scales missed by this kind of analysis. Therefore, designers who want to guarantee performance over a range of edge deviations can take advantage of the proposed DL-augmented technique to affect optimizations that achieve exhaustive tolerancing.

3.4 Speedup performance summary

A thorough analysis of the performance of this approach is presented in Table 2. The expected time that would be taken by a purely RCWA version of our study is listed in one column, with timing for our DL augmented study shown next to it. A speedup column further clarifies the improvement we get over just RCWA.

In this analysis, we include the initial time investment required to collect data and train the networks as part of the total time attributed to the DL augmented approach. This extra time is often not included nor discussed in other DL papers involving optimization in nanophotonics, but it makes a significant difference in realistically understanding the overall utility of the technique. For example, the speedup for a single optimization without accounting for the initial time investment is 14.8 times faster than RCWA alone. However, when the real-time investment is included, this amount is reduced to a 4.37 times speedup.

In reality, oftentimes a single optimization is not enough to solve a problem. In the course of collecting results for the studies in this paper, for example, we ran no fewer than 6 DL augmented optimizations. From this point of view, our full study’s actual speedup was 10.6 times over a purely RCWA-based study. We envision this technique providing a more powerful platform for studying robustness in general than a purely full-wave approach.

Additionally, the speedup reported is the lower bound of what is achievable using our technique, even including the startup time cost. Our choice of RCWA represents a best-case choice for full-wave solver speed. If this solver were replaced with FEM or FDTD techniques, the speedup would significantly increase.

This DL augmented approach to metasurface design opens the doors to more complex studies in the future. As shown in Figure 5, further increases to optimization count or a number of variation tests per supercell make the problem even harder to solve, which increasingly favors our DL method. With finer resolved structures and/or additional tolerance considerations (e.g., edge roughness and sidewall angle), the number of variation tests required per unit cell would increase combinatorically, offering speedups on the order of 20–30 times over a purely full-wave approach.

Figure 5:

Speedup curves using DL augmentation versus pure full-wave solvers.

(A) Speedup when increasing the number of variation tests required by each supercell. (B) Speedup when increasing the number of total evaluations performed in an optimization.

4.1 Future work

Many possible next steps exist to extend this work further. To overcome the necessity to filter designs during optimization, the learning process could be integrated directly into the optimizer. As a now surrogate-assisted optimizer, the model may converge to even higher performances than previously discovered. Other improvements might include developing a smaller DNN model, applying a more physically realistic erosion/dilation model to the data generation, extending the complexity of the supercell parameterization to include more control points per unit cell, and training the model to generalize across more types of fabrication errors (edge roughness, sidewall angle, etc.).