# Inverse design of bright, dielectric metasurfaces color filters based on back-propagation and multi-valued artificial neural networks

Arthur Clini de Souza, Stéphane Lanteri, Hugo Enrique Hernandez-Figueroa, Marco Abbarchi, David Grosso, Badre Kerzabi, Mahmoud Elsawy

► To cite this version:

Arthur Clini de Souza, Stéphane Lanteri, Hugo Enrique Hernandez-Figueroa, Marco Abbarchi, David Grosso, et al.. Inverse design of bright, dielectric metasurfaces color filters based on back-propagation and multi-valued artificial neural networks. SPIE Photonics Europe - Machine Learning in Photonics, Apr 2024, Strasbourg, France. pp.24, 10.1117/12.3016617. hal-04884765

HAL Id: hal-04884765

https://inria.hal.science/hal-04884765v1

Submitted on 13 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Inverse design of bright, dielectric metasurfaces color filters
based on back-propagation
and multi-valued artificial neural networks

Arthur Clini de Souza<sup>1,2,3</sup>, Stéphane Lanteri<sup>1</sup>, Hugo Enrique Hernandez-Figueroa<sup>2</sup>, Marco Abbarchi<sup>3,4</sup>, David Grosso<sup>3,4</sup>, Badre Kerzabi<sup>3</sup>, and Mahmoud Elsawy*<sup>1</sup>

<sup>1</sup>Université Côte d'Azur, Inria, CNRS, LJAD, 06902 Sophia Antipolis Cedex, France

<sup>2</sup>Laboratory of Applied and Computational Electromagnetism (LEMAC), School of Electrical and Computer Engineering (FEEC), University of Campinas (UNICAMP) Campinas, São Paulo, Brazil

<sup>3</sup>Solnil, 95 Rue de la République, 13002 Marseille, France

<sup>4</sup>Université Aix Marseille, CNRS, Université de Toulon, IM2NP, UMR 7334, F-13397 Marseille, France

ABSTRACT

The present work showcases an innovative optimization methodology based on deep learning that combines Multi-Valued Artificial Neural Networks and back-propagation optimization. The methodology addresses the inherent limitations of conventional approaches when employed in isolation. We applied the proposed methodology to design structural color filters that surpasses the sRGB gamut while preserving fabrication constraints.

**Keywords:** Metamaterials, Deep Learning, Photonics

1. INTRODUCTION

In this paper, we build upon our previous work published in Ref. 1. Structural color filters are devices whose reflection or transmission spectra are engineered to achieve a desired perceived color.¹ For this purpose, dielectric photonic metasurfaces have recently emerged as a promising platform for high quality filters.¹⁻⁵ By changing the shape of the resonators, it is possible to control its optical response, such as reflection, transmission and phase.¹,⁶⁻¹⁰ However, the behavior of resonances makes the inverse design process non-intuitive, requiring advanced optimization algorithms such as those based on Deep Learning (DL).¹,⁹⁻¹⁸ Classical optimization approaches rely on a large number of simulations for each color target¹⁹⁻²¹ and require a complete re-initialization of the process for each case.

Here we present a modification of a recently introduced data-driven methodology for efficiently designing vivid and narrow color filter based on photonic metasurfaces.¹ Our approach combines the ability to find suitable designs of Multi-Valued Artificial Neural Networks (MVANNs) with a fully differentiable surrogate model. By propagating back the gradients, the surrogate model is used to enhance the training of the MVANN, preventing poor performance associated with extrapolation while increasing its accuracy when prompted with one of the targets.¹², ¹³, ²² Furthermore, our methodology is capable of finding solutions while considering fabrication constraints (e.g. minimal edge-to-edge-distance, maximal aspect ratio, constant height for each structure, refractive index of the underlying material). Relying solely on 585 simulations we optimize 200 different structures spanning a spectrum of color targets thus producing highly vivid color filters. The optimized geometries exhibit a single sharp resonance response representing all the colors across the visible spectrum, standing as a strong competitor in the field of metasurface color filters.¹, ¹⁸, ²³⁻²⁷

Mahmoud Elsawy: E-mail: mahmoud.elsawy@inria.fr
## 2. GEOMETRY AND SURROGATE MODEL

The simulation setup and dataset are the same as presented in Ref. 1. Fig. 1 shows the metasurface geometry composed of slanted titanium dioxide (TiO$_2$) ridges sitting atop of a silicon dioxide (SiO$_2$) substrate with refractive index $n_s = 1.45$. The refractive index of TiO$_2$ is determined through ellipsometry (available in the supplementary information section and is about $n_{TiO_2} = 2.4$). By slanting the ridges, the symmetry in the propagation direction is broken, resulting in high-quality resonant modes.$^{28,29}$ In this study, four design parameters are considered: the metasurface period (P), resonator height (H), base width (W), and the slant (S), representing how much smaller is the top compared to the base width of the resonator.

<figure><img src="image_2.png" /><figcaption>Figure 1. Schematic representation of the metasurface. Incident light is injected from top with normal incidence with electric field polarized on the x-axis. The simulation comprises of a single resonator with periodic boundary conditions.¹</figcaption></figure>

The details of the dataset generation are shown in detail in Ref. 1. The software used to simulate is the open source python library MEEP.³⁰ In total, 585 different geometries were simulated in 33 minutes and 40 seconds.

Then, we train a feedforward Artificial Neural Network (ANN) surrogate model to predict the spectral responses. Its purpose is twofold. Firstly as a faster estimator of the spectral responses and secondly as a differentiable solver, enabling gradient descent optimization. Therefore, the performance of the surrogate model is of paramount importance as it serves as the base of the inverse design methodology. The details of the surrogate model can be found in Ref. 1.

Figure 2 a) shows how the surrogate model generalizes by comparing its performance on the training data and on the test data. Integrating the wavelength as an input variable has proven to be an efficient way to reduce the dataset size while maintaining acceptable generalization performance.¹

## 3. INVERSE DESIGN

Applying gradient descent over the surrogate model is not a robust inverse design methodology, as the result strongly depends on the initial guess. On the other hand, an ANN could be used to estimate. A regular ANN would be ineffective as the inverse problem presents multiple solutions for a single target. There are numerous solutions in the literature, and the one chosen for this problem is MVANN²² due to its relative simplicity.

MVANN are simply neural networks that branches into a fixed number of outputs. In the scenario of this study, once prompted a target spectrum, each branch will output a different design. The training loss is equivalent as the analogous resistance of resistors arranged in parallel. If one output's loss is low, it will drive the entire model loss to an equivalent value, implying that each output will focus on different sets of solutions.²²

The first step is choosing our target line shape. We utilized two different lineshapes, a Lorentzian function shape resembling a Fano resonance-like shape and a broad resonance. The Lorentzian is defined as follows:

$$L(f) = \frac{\omega^2}{\omega^2 + 4(f - f_0)^2}, \qquad (1)$$
<figure><img src="image_3.png" /><figcaption>Figure 2. a) Comparison between full-wave simulation and the prediction from the surrogate model. Figures on the top row indicate the reflection spectra of three randomly selected samples from the training dataset, while the bottom, from the test dataset. Notably, the surrogate model demonstrates robust generalization capabilities despite the sparsity of the dataset. b) Training history of the Surrogate model.</figcaption></figure>

where $\omega$ represents the Full Width at Half Maximum (FWHM) while $f_0$ corresponds to the central frequency. The units employed are in accordance with MEEP's configuration, where both frequency and FWHM are measured in $\mu m^{-1}$. Additionally, it is worth noting that throughout the optimization process, the FWHM in frequency remains unchanged for all targets. This means that the Quality Factor (QF) of the desired spectra undergoes variations as a function of $f_0$, in accordance with Equation 2.

$$QF := \frac{f_0}{\omega} \qquad (2)$$

The secondary broad target is defined as the sum of two Gaussians

$$B(f) = 0.36 * \exp\left(-\left(\frac{f - f_c - 0.06}{0.045}\right)^2\right) + \exp\left(-\left(\frac{f - f_c + 0.055}{0.077}\right)^2\right) \qquad (3)$$

where $f$ is the frequency and $f_c$ is the center frequency. This target was chosen by observing and comparing different parameters in Eq. 3 with some spectra present in the dataset.

In the original study, one of the largest difficulties was training the MVANN, as the prompted target spectrum was completely outside of the training points forcing us to rely on other techniques to reduce overfitting to increase the network accuracy outside the training points bubble.¹ Here we addressed this problem by coupling the surrogate model on the output of the MVANN. The model could be trained to lower losses, with a lower number of parameters and no regularization. Rendering the MVANN so accurate that the back-propagation step on the inverse design loop was no longer necessary to refine the solutions.

At each training epoch of the MVANN, we compute the regular MVANN loss described in Ref. 1 and update its weights. The next step is prompting the desired targets, computing the solutions with the surrogate model, applying the MVANN loss on the computed spectra and updating its weights again. This extra step moved the back-propagation optimization from the inverse design loop to the training of the MVANN.

The targets added to the training of the MVANN were 100 Lorentzian and 100 broad targets engineered to mimic a resonance found in the dataset.

In this study, rather than considering an uneven height across different structures, we changed the methodology for optimizing the designs with a constrained height (that is more reasonable from a realistic fabrication
point of view). Rather than manually choosing the designs and optimizing them using back-propagation, we created an ANN that condenses all the solutions from the MVANN for all targets and computes the best design for each target, enabling the simultaneous optimization of 200 targets while preserving fabrication constraints such as uniform height and dimensions constraints.

The optimization methodology starts by computing the initial weights matrix of the condenser, which for each target should be a combination of the best designs proposed by the MVANN. To achieve that, we prompt the targets to the MVANN and compute the Mean Square Error (MSE) loss of the spectrum predicted by the surrogate model and the target. We normalize the loss using

$$E_{norm} = 2 \times \frac{E - \mu}{\sqrt{\sigma^2 + \epsilon}} \qquad (4)$$

where $E_{norm}$ is the normalized error, $\mu$ is the average loss for a specific target, $\sigma$ is the variance of the loss for the same target and $\epsilon = 10^{-5}$ a constant to avoid division by zero, meaning that for each target, different $\mu$'s and $\sigma$'s are computed to normalize the error for each target. Then we applied the Softmax function target-wise. Thus, for each target we obtain the percentage of how much each solution of the MVANN should contribute on the final solution based on its accuracy. Therefore, we created an ANN and initialize its weights such that the given solution for each target is the combination of each solution weighted by its accuracy to the target. Then, the solutions are passed to a secondary network, hereby called post-process network, the heights of all solutions are averaged using a linear layer and impose hard limits on the parameters using activation functions, forcing $S < 0.8$ and $W < P - H/2$, on this case, using the Softplus ($\frac{1}{\beta} \log(1 + \exp(\beta x))$) activation function with $\beta = 20$. Softplus behavior is similar to Rectified Linear Unit (ReLU) except it has continuous derivative. Then, the ANN is trained to reduce the loss between the target and the predicted spectrum by the surrogate model.

We present two different scenarios, one in which only the height is constrained and the other in which both the height and the slant angle are constrained. The schematic of the condenser network is shown in Fig. 3.

<figure><img src="image_4.png" /><figcaption>Figure 3. Schematic of the condenser network. We engineer the initial weights of the network such that the solution is a weighted combination of the best MVANN predictions. To do so, each design is fed to an ANN with 4 inputs and 4 outputs and its weight initialization is the identity matrix multiplied by its corresponding contribution. For each target, each solution is post-processed, assuring the parameters stay within specifications and also uniforming variables in all optimizations, such as the height or the slant angle. The proposed results are then processed by the surrogate model and the loss is the MSE between the predicted spectrum and the target. The network is trained with 1000 epochs with Adam optimizer, initial learning rate of $3 \times 10^{-4}$ that is multiplied by 0.995 by each epoch. A patience of 30 was implemented.</figcaption></figure>
Figure 4 represents the results of the optimization with uniform height. And as could be seen, the methodology was able to optimize all the designs efficiently and achieve results as bright as the original results without that constrain.¹ The color of the sharp Lorentzian and the broad target are similar, however, a small difference can be highlighted. The Lorentzian is able to provide brighter colors from blue to the green, while the broader resonance is able to provide better results from yellow to red.

<figure><img src="image_5.png" /><figcaption>Figure 4. (a) Chromaticity diagram depicting the outcomes of the optimization with uniform height. The arrows in the diagram indicate the corresponding designs showcased in (b). The gray lines represents the calculated perceived color corresponding to the target spectrum of the Lorentzian target and the broad target. (b) detailed examination of 6 different targets, displaying their dimensions (in nm) of each geometrical parameter and the corresponding sRGB color representation of the simulated spectrum. Condenser training time: 8 min 33s. Simulation time: 2 min 49s.</figcaption></figure>

In order to obtain a uniform slant angle, we created a layer on the post-process network that computes the cosine of the slant angle given W, S and H

$$ \cos(\theta) = \frac{WS}{\sqrt{4H^2 + W^2 S^2}} \qquad (5) $$

We then average the values with a linear layer, followed by an update of the value of S as

$$ S = \frac{2H}{W \sqrt{1 / \cos(\theta)^2 - 1 + \epsilon}} \qquad (6) $$

where $\epsilon = 10^{-7}$ is a value to prevent division by zero.

Figure 5 represents the optimization results for constraining the height and the slant angle. With the increase of constraints, it is observed a slight decrease in the quality of the results, specially in the red. It is also observed a clustering of designs located in the blue region. Those designs represent solutions not found by the network, rendering solutions without a strong resonance, rendering dim reflections and non-ideal colors.
<figure><img src="image_6.png" /><figcaption>Figure 5. Optimization results for fixing the height and the angle. "A" depicts the base angle of the trapezoid. And as could be seen, the optimization was able to find multiple high efficiency design while sharing the height and the slant angle. Condenser training time: 8 min 6s. Simulation time: 3 min 3s.</figcaption></figure>

## 4. CONCLUSION

In conclusion, we improved our previous methodology,¹ achieving vivid colors through the visible spectrum while preserving strict fabrication constraints. By relying solely on 585 simulations, our approach demonstrates exceptional efficiency in the optimization process, enabling the conception of numerous optimized and cost-effective designs while not relying on providing further high fidelity simulation data. By modeling a resonant metasurface with a broken symmetry in the z-axis and coupling with an efficient optimization method, we were able to conceive structural color filters that ventures beyond the sRGB colors. Furthermore, our DL optimization method was built to provide designs within tight fabrication constraints. Additionally, the optimization methodology is not constrained on the optimization of metasurface color filters and can be used in any data-driven optimization approach, opening new possibilities for the development of various devices.

## APPENDIX A. TiO₂ REFRACTIVE INDEX

The characterization of the two-dimensional flat layers made of TiO₂ was performed in order to determine their optical constants. The measurements were carried out with a Wollam M200V ellipsometer (350-1000 nm) on a layer obtained via sol-gel chemistry and dip-coating on a bulk Si substrate.³¹ The ellipsometric data (Δ and Ψ) are fitted by a Cauchy model to extract refractive index n and extinction coefficient k (being the latter equal to zero in the considered wavelength range). These layers are of optical quality being their surface roughness below 0.1 nm.³²

Measuring the refractive index evolution of a TiO₂ thin layer kept in a ellipsometry-porosimetry chamber with controlled humidity (atmosphere where air and water were progressively mixed) we assessed the adsorption/desorption of water in order to measure its porosity. From the corresponding isotherms we obtain the mesoporosity of the layers³³ (pore volume, size distribution, and interconnection through adapted Kelvin models). Measurements on dense TiO₂ (n > 2.4) account for a porosity well-below 10%.

The dielectric constant of the TiO₂ was fitted by a Lorentzian susceptibility function to be utilized by MEEP. The difference between the ellipsometry data and the fitted curve can be seen in Fig. 6.¹
<figure><img src="image_7.png" /><figcaption>Figure 6. Comparative analysis of the measured refractive index of TiO₂ and the approximation generated through modeling using MEEP. The experimental data capturing the dispersion characteristics were obtained through Ellipsometry spectroscopy detailed below. Subsequently, the collected data were fitted to a Lorentzian susceptibility function. Notably, the approximation achieved through the modeling process closely aligns with the measured dispersion, demonstrating the efficacy of the approach in accurately representing the dispersion behavior of TiO₂.</figcaption></figure>

## APPENDIX B. SURROGATE MODEL TRAINING

The Surrogate model was programmed to train in 1000 epochs, however with patience of 20 and minimum delta of $10^{-6}$, meaning if the validation loss was not reduced by $10^{-6}$ during 20 epochs of training, the training stops. When the validation loss is lower than the lowest loss minus $10^{-6}$, the model parameters were saved and we load the last model parameters saved when the training finishes which happened at epoch 85. Figure 2 b) show the loss history for the training of the surrogate model. In the present study, the difference with our original model, is that the model was trained on pytorch³⁴ instead of tensorflow³⁵ and on a NVIDIA RTX A1000 6GB Laptop GPU instead of google colab NVIDIA T4 GPU. Achieving the same performance after 1 minute of training rather than 10 minutes. The dataset (585 simulations) was divided into training (527), validation (28) and test (30). The Final training loss was $6.7 \times 10^{-4}$, validation loss was $6.5 \times 10^{-4}$, test loss was $1.2 \times 10^{-3}$. The total training time was 29 seconds. Total number of epochs: 69. The optimizer is Adam with initial learning rate of $10^{-3}$. A scheduler was used that reduced the learning rate linearly to $2 \times 10^{-4}$ at epoch 120.

## APPENDIX C. MVANN MODEL TRAINING

The MVANN training was programmed similarly to the surrogate model, the only difference was that it was programmed to stop at epoch 300 and it stopped at epoch 191. To train the MVANN, we created a dataset with 7839 predictions from the surrogate model (Prediction time: 1.6s). The dataset was divided into train (7055) and validation (784). The architecture of the MVANN is shown in Fig. 7.

a)

<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

b)

<figure><img src="image_10.png" /><figcaption>Figure 7. (a) Architecture of the MVANN. (b) Training history of the MVANN. The total training time was 7 minutes and 31 seconds. The surrogate loss represents the loss after prompting the 200 targets and applying the MVANN loss to the computed spectrum.</figcaption></figure>
The optimizer is Adam with initial learning rate of $2 \times 10^{-3}$ and an exponential scheduler that multiplies the learning rate by 0.99 at every epoch.

## APPENDIX D. ANGLE SENSITIVITY

In this section, we analyse the sensitivity for different angles of incidence. The computations were performed with COMSOL Multiphysics and the results can be seen in Fig. 8.

<figure><img src="image_11.png" /><figcaption>Figure 8. (a) Chromaticity diagram depicting the evolution of a color for a given incident angle. (b) Evolution of the spectra for different angles of incidence. Total simulation time: 30 min 11 s.</figcaption></figure>

For the showcase in Fig 8, we chose the highlighted designs from Fig. 5. And as can be seen, it is very sensitive, considering maximum 5 degrees before it losing completely its original color. In all cases it is observed the appearance of a higher frequency resonance that shifts the color to the blue. It is observed that the sharp Lorentzian resonance is more sensitive to the incident angle compared to the broad resonance.

## APPENDIX E. TIME TABLE

<table><caption>Table 1. The original dataset simulation was conducted on 4 threads of an Intel® Xeon® W-2125 Processor operating at 4.0 GHz. While the simulations performed for this study to validate the results were performed on 10 threads of a 13th Gen Intel(R) Core(TM) i7-13700H. The ANN were trained using pytorch on a NVIDIA RTX A1000 6GB Laptop GPU.</caption><thead><tr><th colspan="2">Simulations</th><th colspan="2">Models Training</th><th colspan="2">Constrained height</th><th colspan="2">Constrained height, angle</th></tr></thead><tbody><tr><td>Total</td><td>33 min 40 s</td><td>Surrogate</td><td>29 s</td><td>Simulations</td><td>2 min 49 s</td><td>Simulations</td><td>3 min 3 s</td></tr><tr><td>Average</td><td>3.4 s</td><td>Generation</td><td>1.6 s</td><td>Training</td><td>8 min 33 s</td><td>Training</td><td>8 min 6 s</td></tr><tr><td></td><td></td><td>MVANN</td><td>8 min 22 s</td><td>Other</td><td>1 s</td><td>Other</td><td>1 s</td></tr><tr><td></td><td></td><td>Total</td><td>8 min 53 s</td><td>Total</td><td>11 min 22 s</td><td>Total</td><td>11 min 10 s</td></tr></tbody></table>

## ACKNOWLEDGMENTS

Hugo EHF would like to thank the support of Brazilian government agencies CAPES and CNPq (through project: 312714/2019-0), and the São Paulo Research Foundation (FAPESP) agency through the projects: 2021/11380-5 (CCD/CPTEn) and 2021/06506-0.
REFERENCES

[1] Clini de Souza, A., Lanteri, S., Hernández-Figueroa, H. E., Abbarchi, M., Grosso, D., Kerzabi, B., and Elsawy, M., "Back-propagation optimization and multi-valued artificial neural networks for highly vivid structural color filter metasurfaces," *Scientific Reports* **13**(1), 21352 (2023).

[2] Yang, Y., Kravchenko, I. I., Briggs, D. P., and Valentine, J., "All-dielectric metasurface analogue of electromagnetically induced transparency," *Nature communications* **5**(1), 5753 (2014).

[3] Arbabi, A., Horie, Y., Bagheri, M., and Faraon, A., "Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission," *Nature nanotechnology* **10**(11), 937-943 (2015).

[4] Shalaev, M. I., Sun, J., Tsukernik, A., Pandey, A., Nikolskiy, K., and Litchinitser, N. M., "High-efficiency all-dielectric metasurfaces for ultracompact beam manipulation in transmission mode," *Nano letters* **15**(9), 6261-6266 (2015).

[5] Lin, D., Fan, P., Hasman, E., and Brongersma, M. L., "Dielectric gradient metasurface optical elements," *science* **345**(6194), 298-302 (2014).

[6] Decker, M., Staude, I., Falkner, M., Dominguez, J., Neshev, D. N., Brener, I., Pertsch, T., and Kivshar, Y. S., "High-efficiency dielectric huygens' surfaces," *Advanced Optical Materials* **3**(6), 813-820 (2015).

[7] Yu, Y. F., Zhu, A. Y., Paniagua-Domínguez, R., Fu, Y. H., Luk'yanchuk, B., and Kuznetsov, A. I., "High-transmission dielectric metasurface with $2\pi$ phase control at visible wavelengths," *Laser & Photonics Reviews* **9**(4), 412-418 (2015).

[8] Um, H.-D., Choi, D., Solanki, A., Huang, E., Seo, K., and Habbal, F., "Dynamic selection of visible wavelengths using resonant tio2 nanostructures," *Nanophotonics* (0) (2023).

[9] Lin, R., Valuckas, V., Do, T. T. H., Nemati, A., Kuznetsov, A. I., Teng, J., and Ha, S. T., "An efficient neural optimizer for resonant nanostructures: demonstration of highly-saturated red silicon structural color," *arXiv preprint arXiv:2304.13516* (2023).

[10] Shang, X., Niu, J., Li, H., Li, L., Hu, H., Lu, C., and Shi, L., "Polarization-sensitive structural colors based on anisotropic silicon metasurfaces," in *Photonics*, **10**(4), 448, MDPI (2023).

[11] Ma, W., Cheng, F., Xu, Y., Wen, Q., and Liu, Y., "Probabilistic representation and inverse design of metamaterials based on a deep generative model with semi-supervised learning strategy," *Advanced Materials* **31**(35), 1901111 (2019).

[12] Peurifoy, J., Shen, Y., Jing, L., Yang, Y., Cano-Renteria, F., DeLacy, B. G., Joannopoulos, J. D., Tegmark, M., and Soljačić, M., "Nanophotonic particle simulation and inverse design using artificial neural networks," *Science advances* **4**(6), eaar4206 (2018).

[13] Asano, T. and Noda, S., "Optimization of photonic crystal nanocavities based on deep learning," *Optics express* **26**(25), 32704-32717 (2018).

[14] Jiang, J. and Fan, J. A., "Global optimization of dielectric metasurfaces using a physics-driven neural network," *Nano letters* **19**(8), 5366-5372 (2019).

[15] Wen, F., Jiang, J., and Fan, J. A., "Robust freeform metasurface design based on progressively growing generative networks," *ACS Photonics* **7**(8), 2098-2104 (2020).

[16] Wiecha, P. R. and Muskens, O. L., "Deep learning meets nanophotonics: a generalized accurate predictor for near fields and far fields of arbitrary 3d nanostructures," *Nano letters* **20**(1), 329-338 (2019).

[17] Wiecha, P. R., Arbouet, A., Girard, C., and Muskens, O. L., "Deep learning in nano-photonics: inverse design and beyond," *Photonics Research* **9**(5), B182-B200 (2021).

[18] Hemmatyar, O., Abdollahramezani, S., Kiarashinejad, Y., Zandehshahvar, M., and Adibi, A., "Full color generation with fano-type resonant hfo 2 nanopillars designed by a deep-learning approach," *Nanoscale* **11**(44), 21266-21274 (2019).

[19] Elsawy, M. M., Lanteri, S., Duvigneau, R., Brière, G., Mohamed, M. S., and Genevet, P., "Global optimization of metasurface designs using statistical learning methods," *Scientific reports* **9**(1), 1-15 (2019).

[20] Xu, D., Luo, Y., Luo, J., Pu, M., Zhang, Y., Ha, Y., and Luo, X., "Efficient design of a dielectric metasurface with transfer learning and genetic algorithm," *Optical Materials Express* **11**(7), 1852-1862 (2021).

[21] Jafar-Zanjani, S., Inampudi, S., and Mosallaei, H., "Adaptive genetic algorithm for optical metasurfaces design," *Scientific reports* **8**(1), 11040 (2018).
[22] Zhang, C., Jin, J., Na, W., Zhang, Q.-J., and Yu, M., "Multivalued neural network inverse modeling and applications to microwave filters," *IEEE Transactions on Microwave Theory and Techniques* **66**(8), 3781-3797 (2018).

[23] Liu, H., Peng, H., Li, K., Lu, L., Deng, J., Liu, Y., Qiu, C., Li, G., and Cheng, X., "Transfer printing of solution-processed 3d zno nanostructures with ultra-high yield for flexible metasurface color filter," *Advanced Materials Interfaces* **9**(21), 2101963 (2022).

[24] Panda, S. S. and Hegde, R. S., "Transmission-mode all-dielectric metasurface color filter arrays designed by evolutionary search," *Journal of Nanophotonics* **14**(1), 016014-016014 (2020).

[25] Wang, Y., Huang, W., Lin, Y.-S., and Yang, B.-R., "A tunable color filter using a hybrid metasurface composed of zno nanopillars and ag nanoholes," *Nanoscale Advances* **4**(17), 3624-3633 (2022).

[26] Khaidarov, E., Eschimese, D., Lai, K. H., Huang, A., Fu, Y. H., Lin, Q., Paniagua-Dominguez, R., and Kuznetsov, A. I., "Large-scale vivid metasurface color printing using advanced 12-in. immersion photolithography," *Scientific Reports* **12**(1), 14044 (2022).

[27] Yang, W., Xiao, S., Song, Q., Liu, Y., Wu, Y., Wang, S., Yu, J., Han, J., and Tsai, D.-P., "All-dielectric metasurface for high-performance structural color," *Nature communications* **11**(1), 1864 (2020).

[28] Barton III, D., Lawrence, M., and Dionne, J., "Wavefront shaping and modulation with resonant electro-optic phase gradient metasurfaces," *Applied Physics Letters* **118**(7), 071104 (2021).

[29] Lawrence, M., Barton III, D. R., Dixon, J., Song, J.-H., van de Groep, J., Brongersma, M. L., and Dionne, J. A., "High quality factor phase gradient metasurfaces," *Nature Nanotechnology* **15**(11), 956-961 (2020).

[30] Oskooi, A. F., Roundy, D., Ibanescu, M., Bermel, P., Joannopoulos, J. D., and Johnson, S. G., "Meep: A flexible free-software package for electromagnetic simulations by the fdtd method," *Computer Physics Communications* **181**(3), 687-702 (2010).

[31] Bindini, E., Naudin, G., Faustini, M., Grosso, D., and Boissière, C., "Critical role of the atmosphere in dip-coating process," *The Journal of Physical Chemistry C* **121**(27), 14572-14580 (2017).

[32] Chehadi, Z., Bouabdellaoui, M., Modaresialam, M., Bottein, T., Salvalaglio, M., Bollani, M., Grosso, D., and Abbarchi, M., "Scalable disordered hyperuniform architectures via nanoimprint lithography of metal oxides," *ACS Applied Materials & Interfaces* **13**(31), 37761-37774 (2021).

[33] Loizillon, J., Baumgartner, B., Sinturel, C., Abbarchi, M., Lendl, B., and Grosso, D., "In-depth study of coating multimodal porosity using ellipsometry porosimetry in desorption scanning mode," *The Journal of Physical Chemistry C* **123**(38), 23464-23479 (2019).

[34] Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., DeVito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S., "Pytorch: An imperative style, high-performance deep learning library," in [*Advances in Neural Information Processing Systems* 32], 8024-8035, Curran Associates, Inc. (2019).

[35] Abadi, M., Barham, P., Chen, J., Chen, Z., Davis, A., Dean, J., Devin, M., Ghemawat, S., Irving, G., Isard, M., et al., "Tensorflow: a system for large-scale machine learning.," in [*Osdi*], **16**(2016), 265-283, Savannah, GA, USA (2016).