1. Introduction

The need for miniaturization is a current challenge for optical engineering. Optical metasurfaces [1,2] are precision-fabricated nanostructured thin-films which exploit the rich optical physics of nanoscale resonators. Metasurfaces could potentially enable drastic miniaturization of optical devices and enable novel applications [1,3–5]. In particular, the possibility to create arrayed spectral filters is highly attractive for optical imaging and displays [6] where metasurfaces enjoy many advantages compared to conventional technologies [7]. They are more resistant to age and environmental degradation and can better withstand UV radiation and higher temperatures. Initial investigations into metasurface spectral filters used plasmonic metals [8–13], but there has been a gradual shift towards all-dielectric metasurfaces (ADM) [14–17] owing to their comparatively lesser optical absorption.

Metasurface spectral filters can be designed to operate in reflection-mode or transmission-mode. Operating in transmission-mode, they can potentially replace conventional color filter arrays (CFA) for realizing submicron-size pixels in backside-illuminated CMOS image sensors [18–21] and enable novel multispectral image sensors. Despite their potential, the study of ADM CFA for transmission-mode has received considerably less attention compared to reflection-mode structural color filters. The noticeable absorption of silicon (especially at lower wavelengths) is an impediment in realizing transmission-mode metasurface devices for visible wavelength operation. Although, good transmission efficiencies have been reported in [20–24], the colors were subtractive in nature (CMYK colors are more appropriate for reflected light applications). Kanamori and coworkers [25] report good transmission efficiency for RGB filters using silicon gratings, but the polarisation sensitivity of this structure is a limiting factor. Additionally, they have not considered the absorptive silicon back-plane. To the best of our knowledge, the reports by Horie and coworkers [18] and by Berzins and coworkers [19] are the only ones that focus on transmission-mode dielectric metasurfaces for primary color image sensor CFAs. The cylindrical dielectric particle metasurfaces reported in [19] gave filters which lacked purity and were found near the white point of the CIE 1931 color space. Although the void metasurfaces reported in [18] tend to have an improved color gamut as compared to [19], the purity of blue color is extremely poor because of reduced transmission efficiency. In order to compete with conventional technologies, the performance parameters of silicon transmission-mode color filters have to further improve in order to meet or surpass those of conventional filters in terms of color purity, color gamut and intra-pixel uniformity.

The meta-atom shape provides additional design degrees of freedom which could potentially enhance performance for a given material platform. Although, the space of potential shapes within the reach of modern nanofabrication is very large, an overwhelming amount of experimental reports use high-symmetry shapes like circles, ellipses and rectangles; although simple, these geometries may not be optimal. To construct all-dielectric metasurfaces, there are several choices for medium to high-index CMOS-compatible optical dielectrics [20]. Computational inverse design methods [26] are neccessary to explore the vast design space of material-structure combinations. A plethora of inverse design techniques have been reported for metasurface design [27] including data-driven methods like Deep Learning [28–34]. In particular, topology optimization [35–37] is a popular local optimization technique that can search over a broader set of shapes. However, local optimizers may get stuck at suboptimal designs which is especially problematic in non-convex optimization scenarios like metasurface design [38]. Evolutionary techniques [39,40] are global optimizers and have been considered for metasurface inverse design [27,38]. Global optimizers can better harness heterogenous computing architectures because they are inherently parallel.

The nanofabrication process introduces stochastic perturbations that may degrade the performance of metasurfaces due to the inherent shape-dependent properties of constituent meta-atoms. Thus, the inverse design process must ensure that the resulting designs exhibit some robustness to fabrication errors. Robust design has been considered in topology optimizers; but, in addition to their tendency to get stuck at local optima, they require special care to restrict the search space to feasible designs that meet fabrication constraints [41,42]. Also, as the resulting shapes are non-trivial, there is no easy way to assess the fidelity of nanofabrication of a particular design which makes it difficult to guide the optimizer towards robust designs [43]. On the other hand, the disadvantage of evolutionary approaches to robust design is the vastly increased convergence time for structures with high dimensionality [44].

We report a robust inverse design method for meta-atom design considering the example of polygon-shaped meta-atoms [40,45]. This work builds on our previous contribution [40] which did not feature the robust design methodology and only considered idealized polygon shapes without smooth boundaries. To the best of our knowledge, this is the first paper that addresses the issue of intrapixel color nonuniformity, an important performance metric for CFAs. We propose a simple algorithm based on vertex dithering and spline interpolation based polygon smoothing that can recreate the fabrication imperfections of nanofabricated metasurfaces. Specifically, the statistical model is seen to recreate the shape perturbations associated with electron beam fabrication of silicon metasurfaces. The modified fitness evaluation step is designed to penalize designs whose spectra are highly sensitive to fabrication imperfections. By using the concept of orthogonal arrays, the computational burden of fitness evaluation is reduced by resorting to a partial fraction study instead of a full fraction study including all corner case designs.

A standard multi-island Differential Evolution (DE) is the evolutionary optimizer used here with a modified fitness evaluation step that uses the statistical model discussed above. Differential Evolution is chosen as it is simple to implement real valued vector encoding compared to other approaches. Additionally, reports in literature have suggested that DE has the following advantages [46]: (1) it can find the true global minimum of a multi modal search space regardless of the initial parameter values; (2) it entails fast convergence; and (3) it requires a small number of control parameters. The approach proposed here can also be implemented with several other evolutionary algorithms like genetic algorithms and particle swarm optimization. For the design of primary RGB color filters, our inverse design procedure predicts performance improvement over previously published reports. The paper is organized as follows. In section 2, the geometry and the inverse design algorithm is described. In section 3, the results are presented and discussed before concluding the paper in section 4.

2. Robust inverse design methodology

2.1 Description of geometry

As seen in Fig. 1, we consider a simplified layered model [18,19] of a backside-illuminated CMOS image sensor integrated with ADM color filters with the following assumptions: (1) pixel-size and inter-pixel effects are neglected (as reported in [19], this is a good approximation for pixels larger than 0.5 µm); (2) the photodetector layer is replaced by a semi-infinite crystaline silicon backlayer. This considers reflections at the first interface between the silica spacer and the underlying silicon layer but neglects role of interfaces further below. Both particle resonators and the complementary void resonators (shown in Fig. 1(b)) are considered for the design process.

 

Fig. 1. Geometry of the proposed structure. (a) Construction of an exemplar polygonal meta-atom of 8-fold symmetry. Insets show the complete polygon and its edge-refined version. (b) Perspective view of the exemplar particle and void polygon-shaped resonators with layer descriptions. (c) illustrates that the set generated by the construction procedure contains many common shapes.

Download Full Size | PDF

Figure 1(a) depicts the meta-atom construction. The set of polygon shapes chosen for the meta-atoms is bigger than the set containing high-symmetry shapes but smaller than that considered in topology optimization procedures. The meta-atoms are described explicitly by the coordinates of the polygon vertices which is simpler than the implicit shape definitions used in topology optimizations; furthermore, as the polygons remain as individual resonators, their optical response is easier to interpret. The polygon is specified by its vertices in a curvilinear coordinate system in terms of $r$, $\alpha$ and $\theta$, where $r$ represents the distance from the centre of the unit-cell. $\alpha$ is the angle of the last point of the octant (yellow point in Fig. 1(a)) and the angles $\theta$ depict a point’s angle with respect to it’s previous point. For $n$ points per octant, the values of $r$, $\alpha$ and $\theta$ are chosen as follows: $\displaystyle r_1,r_2,\ldots .., r_n < p/2, (p =$ periodicity), $\displaystyle \theta _1, \theta _2,\ldots .., \theta _n < \alpha$ with $\displaystyle \sum _{i=1}^{n} \theta _i = \alpha$ and $\displaystyle \alpha < 45^{\circ }$. The angles $\theta$ sum up to $\alpha$, they are chosen according to $\displaystyle \theta _i = \frac {\phi _i}{\sum _{i = 1}^{n}\phi _i} \hspace {0.1cm}\alpha$ where the weights $\phi _i$ are chosen randomly. Finally, the polygons’ vertices are connected in a 2-D spline curve in order to generate the smooth polygons. The insets of Fig. 1(a) shows the full polygon along with its edge-smoothed version. The shapes considered throughout this article are limited to an eight-fold symmetry to ensure polarisation-independent optical response. The thickness of the metasurface layer has been kept constant to 80 nm as a trade-off between purity of red and blue colors as discussed in [40].

Although this method of defining polygon shapes does not generate shapes like H and T metasurfaces, they still cover a number of general shapes (see Fig. 1(c)) like triangles, rectangles, squares, crosses and stars. Embedding the polygon design with a spline joining of vertices, the set of shapes can be broadened to include ellipses, Steiner ellipses, circles and other round-edged shapes. The proposed polygon shapes describe a much broader set of polygons than the approach in [45] and our previous work [40]. Inampudi and coworkers [45] restrict the shapes by evenly spacing the angles $\theta$.

2.2 Inverse design algorithm

A multi-island Differential Evolution (DE) (see Fig. 2(a)) optimization is employed in this work for the robust inverse design of the ADM color filters. DE is a stochastic optimization technique that takes a randomly initialized set of solutions (called the "population") and performs nature inspired techniques like mutation and crossover to produce a new set of solutions [47,48]. A selection process ensues where parent solutions are compared with the offspring and fittest solutions among them are promoted to the next generation. These steps are repeated until a predetermined termination criterion is fulfilled. In the multiple island version, the total population is divided evenly among several islands and DE is performed in each islands for multiple number of epochs (each epoch contains multiple number of iterations). At the end of each epoch, the best solution are exchanged among the islands in a round-robin fashion. The multi-island approach can harness heterogeneous computing hardware and prevents the optimizer from getting stuck at local optima.

 

Fig. 2. (a) Schematic of the Multi-island Differential Evolution algorithm. (b) Block diagram of fitness evaluation procedure for robust inverse design. (c) Typical progression of a 3-island DE optimization (population size = 35) on the set of polygons with 24 edges and 8-fold symmetry. Insets show the designs with the best fitness at the end of each epoch in each of the 3 islands.

Download Full Size | PDF

In our case, an individual in the population is represented by a solution vector $\vec {x}$ that contains the radii ($r$), angle weights ($\phi$) and maximum angles ($\alpha$):

However, in a robust inverse design process, the fitness of a solution vector should consider not only the design performance, but also its sensitivity to fabrication imperfections. In other words, the fitness estimation should be carried out on a set of polygons rather than the single polygon because a solution vector can result in different polygons due to stochasticity of the fabrication process. The set constituted by the base polygon $P_0$ and its perturbed variants $P_1, P_2 \cdots$ is infinitely large. For practical purposes, we consider a representative subset $P$ of size $v$ given by:

The set of transmittance spectra $T$ of the polygons in $P$, $\displaystyle \{T_0(\lambda ), T_1(\lambda ) \hspace {0.1cm} \cdots \hspace {0.1cm}. T_{v-1}(\lambda )\}$, can be used to derive the following mean and variance estimates over $P$:

For constructing the set $P$, we need to determine the range of variation of each of the parameters of the solution vector $\vec {x}$ which is dictated by the lithographic process. Within the entire swing of a parameter, the possible errors can be binned into predetermined number of error levels. In a problem with $a$ parameters and $b$ number of error levels each, the total number of corner cases is thus $a^b$. A full factorial design covering all corner cases is prohibitively expensive. On the other hand, a partial fraction study (e.g. random Gaussian distribution of all the errors) may result in the loss of valuable information. Following the proposal of Genici Taguchi [50], we use a partial fraction study, where full information of the control variables is acquired with a minimum number of case studies using a special set of arrays called the orthogonal arrays. A column in the orthogonal array table represents the different error levels taken by a variable, while a row represents a single case with a set of one variant from each variable. All the columns in the table are mutually orthogonal. The table is designed in such a way that each error levels in a single variable appears equal number of times, obeying the balancing property of the arrays. While there are many standard orthogonal arrays available, each of the array is meant for a specific number of independent design variables and levels. The minimum number of case studies that needs to be carried out is determined by the formula:

Figure 2(b) shows the block diagram for constructing the set $P$ and evaluating the fitness for a given solution vector $\vec {x}$. In this work, the radius error bound was set to ±15 nm and we binned the errors into three levels (0, +15, -15). The alterations in angles are not considered in the modelling of perturbations, because, in the presence of smoothing, it was found that the shape variations caused by small angle perturbations were not noticeable. Using the formula in Eq. (8), it can be determined that for 13 vertices with 3 levels and 3 vertices with 2 levels of errors, we need 30 sets of perturbed variants. The standard Taguchi’s orthogonal array table nearest to 30 that can accommodate this variation is the L36 table. Similarly, for 24 vertices with 3 levels of errors, the L54 table was suitable to construct the orthogonal array. In section 3.3, the efficacy of the Taguchi’s orthogonal arrays is demonstrated using the example of a cross-shaped meta-atom.

3. Results and discussion

The design, simulation, and optimization were run in a system with the following specifications. Processor: Intel Core i7-5960X CPU at 3.00 GHz; number of cores: 8; and RAM: 128 GB. The computer code used in this work have been made publicly available at [52]. The time taken for getting the transmittance spectrum of a single shape with 30 discrete samples of wavelength between 400 nm and 700 nm was 30 seconds. 37 number of variations of each shape (one original and 36 generated from L36 table) were run; thus making the time required for one simulation to be $\approx$ 18.5 minutes. The DE schemes implemented here are DE/rand/1/bin (mutant is derived through one difference operation between two random solutions and a binary crossover takes place between parent and mutant solution vector to generate the offspring) for exploration and DE/best/1/bin (instead of random solution vector, the best solution vector participates in deriving mutant vector) for exploitation. The optimization was operated on three islands with population size of 35 each. The progression of a typical optimization run is seen in Fig. 2(c). A good rate of convergence was observed with mutation rate $(m_r) = 0.4$ and crossover rate $(c_r) = 0.5$, with the convergence criteria being satisfied within eight generations (20 iterations each). The optimizations were run for 10 generations, taking the time frame for a single optimization to be $\approx$ 2.5 days. Best results for Red and Green filters were obtained using 2 points per octant and for blue filters, using 3-points per octant. No significant improvement in performance was observed when the number of points per octant was increased beyond 3.

The fitness evaluations during optimization were performed using the Stanford Stratified Structure Solver ($S^4$) which uses the Rigorous Coupled Wave Analysis (RCWA) technique along with S-matrix algorithm to solve Maxwell’s equations in layered periodic structures. The number of basis functions were set to 50 in all cases. The optical responses of the final designs were verified with full-wave simulations in CST Microwave Studio (MWS). The full wave verification simulations used the frequency domain finite element method (FEM) solver of CST MWS with unit-cell boundary conditions and Floquet port excitation. The FEM solver used a tetrahedral meshing scheme, started out with a step size of 4 per smallest simulation wavelength and a minimum of 10 for any given edge. The solver then iteratively refined the mesh before converging to a final mesh geometry. A good agreement was achieved between the optical responses given by both tools when compared with published results of Horie et. al [18].

3.1 Optical response of the proposed structure

The optical response of the proposed structure is influenced by several factors including the dispersion of silicon, the modes of the particle/void nanoresonators, the periodic interaction among resonators, and, the multiple reflections at the various interfaces. The spectral response for a cross-shaped void meta-atom is considered in Fig. 3(a,b). The high-index dielectric material considered for this work is poly-crystalline silicon (poly-Si) on top of a silica spacer. Silicon has high optical absorption in the lower visible wavelength and nearly flat dispersion with negligible absorption in the the near IR range (see inset in Fig. 3(a)). Thus the low-wavelength response is dominated by absorption leading to a low transmittance for most polygon shapes. It was observed that shapes with low values of the fill fraction (areal fraction of silicon in a unit cell) exhibited relatively reduced absorption in the low-wavelength region. For periodicity values less than 350 nm, the transmittance spectra at higher wavelengths (beyond 900 nm) were smooth and did not exhibit sharp features regardless of the polygon shape. In the mid and higher visible wavelengths, the spectra were highly shape dependent and exhibit several sharp spectral features. This was particularly evident for the case when the backlayer is absent with sharp transmission dips accompanied by high absorption (see markings $P_1$ to $P_4$ in Fig. 3(a)).

 

Fig. 3. Numerically simulated optical response (absorption (a) and transmittance (b) spectra) of a structure using an exemplar cross-shaped void meta-atom generated using the process described in Fig. 1. Structures with and without the backlayer are considered for comparison as seen in the insets of subfigure (b). Inset in (a) shows the $n-\kappa$ data of the poly-Si material. Periodicity, length and width of the cross-shaped arms are 333 nm, 200 nm and 40 nm respectively. Silica spacer is 115 nm thick, metasurface is 80 nm thick.

Download Full Size | PDF

The sharp spectral features seen here exhibit much smaller linewidths than the Mie resonances in individual particle [53] and void resonators [54]. In the case of structures without a backlayer, these resonances can be interpreted as the vertical Fabry-Perot resonances of the Bloch modes of the periodic structure [55]. In the structure with a backlayer, we have three participating interfaces: (1) air-metasurface interface, (2) metasurface-spacer interface and (3) spacer-backlayer interface. In the structure without a backlayer, we only have two participating interfaces: (1) top air-metasurface interface and (2) bottom metasurface-air interface. For structures with a backlayer, reflections arising from the third interface alter the resonance conditions leading to spectral shifts and suppression of the sharp peaks. Since there is no simple relationship between the polygon shape and its spectral response, numerical optimization is needed to determine the optimal shape.

3.2 Experimental verification of the shape perturbation model

In this subsection, we demonstrate the effectiveness of the algorithm for generating the shape perturbations against SEM micrographs of two metasurface samples fabricated using electron-beam lithography: (1) a cylindrical geometry and (2) a cross-shaped geometry reported in [17]. The metasurface with the cylindrical geometry [51] had a footprint of 50 µm×50 µm and contained single and arrayed reflection-mode filters of varying pixel widths. The metasurface was fabricated on an ITO coated glass substrate. The substrate was cleaned using acetone followed by isopropyl alcohol (IPA). The fabrication process begins with the thin film deposition of a 220 nm thick layer of amorphous silicon on the substrate by plasma enhanced chemical vapor deposition (PECVD) technique. Followed by this, a negative electron beam resist hydrogen silsesquioxane (HSQ) was deposited by using spin coating. The e-beam writing process was then performed on the resist using EBL tool (Raith – eLINE). After developing the resist, the nanoscale patterns were transferred using inductively coupled plasma based reactive ion etching (PlasmaLab system 100 ICP 180, Oxford instruments). The material characterization of the fabricated device was done by using scanning electron microscope (ULTRA 55, Carl Zeiss AG). A visible light microscope (Leica DM2500) was used to obtain the color image of the light reflected from the metasurface. Figure 4(a) and (b) depicts one 9-element filter array fabricated using the above process. The optical micrographs indicate noticeable nonuniformity in the pixels, particularly, at the pixel edges. In the SEM micrographs, it was observed that the side-walls were nearly vertical and top surfaces nearly flat; however the lateral shapes exhibit noticeable shape variations. The intrapixel nonuniformity can result from the shape variance as well as interpixel crosstalk at the boundaries.

 

Fig. 4. Verification of the efficacy of the algorithm for generating perturbed variants of a given shape in an experimental setting. (a) SEM and optical micrographs of reflection mode color filters with cylindrical meta-atoms reported in [51]. (b) Enhanced view of one pixel in (a) showing the irregularities in the shape and its impact on intra-pixel color uniformity. (c) (Left panel) Traced circumferences of selected resonators in (b) (Right panel) Shapes generated by the algorithm. (d, e) Spread of chromaticity coordinates for the shape sets seen in (c). (f) Verification results for cross shaped resonators reported in [17]. SEM micrographs in (f) are reproduced with permission from [17] ©2017, American Chemical Society.

Download Full Size | PDF

To assess the algorithm for generating shape variations, the outlines of the selected nanoresonators in the SEM micrograph of Fig. 4(b) are compared with those generated by the algorithm in Fig. 4(c) for the cylindrical geometry. For the cross shape, the comparison is shown in Fig. 4(f). The outlines were digitally traced using WebPlotDigitizer [56] and the highest radial error was found to be 15 nm. A good resemblance can be observed between the experimental and digital models in terms of shape for both geometries. Next, numerical simulations of the reflectance spectrum of the digitally traced outlines were compared with those of outlines generated by the algorithm for the cylindrical geometry. The chromaticity coordinates of the SEM outlines is presented in Fig. 4(d) and for the algorithm generated outlines in Fig. 4(e), respectively. Both these plots shows similar spreads and are in agreement with the color variations seen in the optical micrograph of Fig. 4(b).

3.3 Effectiveness of the orthogonal array approach

In this subsection, we take an example geometry and show that the orthogonal array can efficiently capture the variations in the optical response of a given shape subject to fabrication imperfections. The transmittance spectra (dotted lines) of twelve different perturbed variants (generation process illustrated in Fig. 5(a)) are compared to that of the base polygon (solid lines) in Fig. 5(b) along with their chromaticity coordinates. The chromaticity coordinates of all perturbed variants of a given shape considering a maximum error limit define a region in the chromaticity diagram whose extent determines its sensitivity to fabrication perturbations. With 16-sided polygons, considering 2 level of errors, the number of possible corner cases is $2^{16}$, which is impractical to calculate. In Fig. 5(c), the chromaticity spread is evaluated for three sets of shapes: (1) 500 randomly selected corner cases; (2) 50 shapes considering normally distributed error (zero mean and standard deviation of 10 nm); and (3) 16 shapes generated from the L16 Taguchi’s orthogonal array table. It can be observed that the 16 samples generated by the L16 table can estimate the spread as well as set (1) with a considerable reduction in computational burden.

 

Fig. 5. Demonstration of the efficacy of Taguchi orthogonal array method for the example of a cross-shaped polygon. (a) Illustration of the process for generating a perturbed variant of a base polygon (i) including vertex dithering (ii) and edge-smoothing (iii). (b) Comparison of the transmittance spectra and chromaticity coordinates of the base polygon (cyan cross) and 12 different perturbed variants (white circles). The color spread induced by shape perturbations is clearly noticeable. (c) Comparison of performance ambit for three different shape sets: (1) 500 randomly chosen corner cases, (2) 50 randomly generated shapes with a Gaussian variation and (3) 16 shapes generated by Taguchi’s L16 table.

Download Full Size | PDF

3.4 Performance of designed RGB filters

In Fig. 6, we summarize the performance results for the filters obtained through the optimization process described above. The results are compared with the results for transmission-mode RGB filters reported by other workers in both the plasmonic [12] and silicon [18,40] material platforms. The geometrical shapes of the optimal filters for blue, green and red colors are seen in the insets of Fig. 6(c) and geometrical parameters are listed in Table 1. It is seen that the color purity and color gamut is enhanced noticeably. Particularly, improved blue and green filters are observed. The performance of blue and green filters is likely to be problematic for the silicon material platform.

 

Fig. 6. Performance of the designed RGB filters for color gamut (a) and perturbation ambit (b) compared with Horie et. al. [18] and Panda et. al [40]. In (b), the ellipses represent the perturbation ambit of the filters from this work (solid) and filters in [18] (dotted). (c) Transmission spectra of the robust color filters compared with [18] (spectral variations due to perturbations is also seen along with the scaled CMFs shown for reference). (d) Performance of index-matched single and 2-layered metasurface color filters. (e) Pictorial comparison along several metrics with previously reported plasmonic [12] and ADM [18,19] color filters.

Download Full Size | PDF

Table 1. Geometry specifications of optimized design for all three configurations seen in Fig. 6(d).

View Table

In Fig. 6(b), perturbation ambit (generated using the Taguchi orthogonal array) for chromaticity in each filters are shown by drawing ellipses in the CIE 1931 color space chromaticity diagram of Fig. 6(b) and by shaded spectra in Fig. 6(c). A similar procedure is adopted for the cylindrical shapes reported in [18] and the corresponding perturbation ambits are shown for comparison. Optimizing for robustness may degrade the color gamut to some extent. But, it is seen that designs which do not explicitly optimize for robustness will lead to a high degree of intra-pixel color nonuniformity.

An extensive comparison between our work, plasmonic transmission-mode CFAs [12] and all dielectric transmission-mode CFAs [18,19] was carried out in terms of transmission efficiency, degree of freedom, RGB gamut and uniformity with respect to shape perturbation whose results are summarized in Fig. 6(e). It can be observed that the structures obtained in this work provide a good transmission efficiency along with an enhanced color gamut with respect to other dielectric color filters reported in [18,19]. However, the color gamut observed are not as good as the gamut of a plasmonic color filter arrays which could indicate the ultimate limitations imposed by the dielectric permittivity of silicon.

In [57], multiple layers of same filters in a cascaded configuration are used to widen the gamut; however, this technique will decrease the transmission efficiency and complicate the fabrication process. Instead, we considered single and two-layered metasurfaces (with a silica spacer of 115 nm in-between the layers) enclosed with a polymer index matching layers as suggested in [23]. After the optimization was run, a significantly enhanced gamut was observed (as plotted in Fig. 6(d)) for both single and double layered configurations in comparison to a single layered shape without the capping. Table 1 lists the optimal geometrical parameters for these configurations.

4. Conclusions

We present the robust inverse design of transmission-mode metasurface color filters and predict significant performance improvements in comparison to earlier reports which used simple shapes. Further extensions of the design process can be considered for minimizing the inter-pixel cross-talk for a predetermined pixel size [19] and for designing multispectral filters [58,59]. This technique can be easily extended to the design of high performance metasurface structural colors and, broadly, for the design of other pixellated metasurface design problems, notably, the design of extended unit cell meta-gratings [60,61]. The design approach is complementary to ongoing progress on improving the optical properties of silicon for metasurface applications [20] (i. e. hydrogenated silicon, silicon nitride (Si$_3$N$_4$), silicon rich silicon-nitride (SRN)). The robust design process involves more computation in comparison to standard evolutionary algorithms. Surrogate models, including Deep Neural Network surrogate models [62,63] could be used to further accelerate the robust inverse design algorithm.