# Inverse design of arbitrary large-scale metalenses with a physically-tailored deep learning model

Arthur Clini de Souza, Paloma Elias da Silva Pellegrini, Silvia Vaz Guerra Nista, Stéphane Lanteri, Hugo Enrique Hernandez-Figueroa, Marco Abbarchi, Badre Kerzabi, Mahmoud Elsawy

► To cite this version:

Arthur Clini de Souza, Paloma Elias da Silva Pellegrini, Silvia Vaz Guerra Nista, Stéphane Lanteri, Hugo Enrique Hernandez-Figueroa, et al.. Inverse design of arbitrary large-scale metalenses with a physically-tailored deep learning model. Laser and Photonics Reviews, In press. hal-05333662

HAL Id: hal-05333662

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

Submitted on 27 Oct 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.

Public Domain
ORIGINAL ARTICLE

OPEN ACCESS

# Inverse design of arbitrary large-scale metalenses with a physically-tailored deep learning model

Arthur Clini de Souza<fnref n="1" />,² Ⓒ | Paloma Elias da Silva Pellegrini³ | Silvia Vaz Guerra Nista³ | Stéphane Lanteri<fnref n="1" /> | Hugo Enrique Hernandez-Figueroa⁴ | Marco Abbarchi²,⁵ | Badre Kerzabi² | Mahmoud Elsawy*<fnref n="1" /> Ⓒ

**Correspondence:** Mahmoud Elsawy (mahmoud.elsawy@inria.fr)

**Received:** 20 January 2023 | **Revised:** 21 January 2023 | **Accepted:** 22 January 2023

**Funding:**

**Keywords:** Large-scale metalens | Optimization | Inverse design | Deep Neural Network | Multi-Valued Artificial Neural Network | Metasurface | Augmented Reality | Virtual Reality

## ABSTRACT

Efficient large-scale metalenses with oblique incidence have drawn a lot of attention for designing future augmented and virtual reality devices. Classical methods based on traditional lookup tables for designing such metalenses lacked efficiency and presented several engineering challenges to achieve high-quality images. This study introduces a novel inverse design methodology to design large-scale metalenses that combines high performance, scalability, and exceptional design flexibility. Our innovative approach assembles multiple optimized beam deflectors that redirect incident light to the focal point. Leveraging deep learning models and incorporating a physics-aware assembling strategy into the optimization process, the high efficiency of the meta-deflectors is translated into high focusing efficiency of the final metalens design. Once the base models are trained, the inverse design process can generate large-scale metalenses with varying sizes, numerical apertures, and focal positions, without the need to re-initialize the entire optimization process, enabling rapid prototyping and versatility in designing optical devices.

## 1 | Introduction

A key technological challenge in Virtual and Augmented Reality (respectively, VR and AR) is the conception of lightweight and efficient smart glasses capable of displaying bright and high-resolution images [1]. However, current designs are based on bulky refractive elements or inefficient diffractive elements [2]. Disruptive advances in the latter design approaches will serve as a pivotal factor in the miniaturization of such devices. Nowadays, metasurfaces are the main focus of such advances, owing to their peculiar properties in controlling the incoming light from the micro-display towards the eyes [3].

Metasurfaces and in particular metalenses are planar photonic devices composed of subwavelength optical elements called meta-atoms that are shaped and organized in a wise way [4]. The applications of metalenses are extensive and span a wide range

of photonic devices, such as imaging, structural color filtering, holography, cloaking, AR, and VR [5, 6, 7, 8, 9, 10]. However, one of the main challenges in metalens engineering is the design of large-scale structures that maintain high efficiency and low aberrations without violating the material and fabrication constraints.

In general, two strategies can be used to design a metalens, namely Local Periodic Approximation (LPA) and inverse design such as topological optimization with gradient-based algorithms, and geometrical optimization using, e.g., evolutionary algorithms or Bayesian optimization methods [11, 12, 13, 14]. Nevertheless, each method has its own limitation and is suitable for different applications. On one hand, in LPA setting, the metalens is designed by manually assembling elements from a lookup table according to the required response. Each element of this lookup table

<footnote n="1">¹Université Côte d'Azur, Inria, CNRS, LJAD, 06902 Sophia Antipolis Cedex, France | ²SOLNIL, Pépinière d'entreprises Luminy Biotech, 163 avenue de Luminy, France | ³Center for Semiconductor Components and Nanotechnology, University of Campinas (UNICAMP), Campinas, São Paulo, Brazil | ⁴Laboratory of Applied and Computational Electromagnetism (LEMAC), School of Electrical and Computer Engineering (FEEC), University of Campinas (UNICAMP), Campinas, São Paulo, Brazil ⁵Université Aix Marseille, CNRS, Université de Toulon, IM2NP, UMR 7334, F-13397 Marseille, France |</footnote>

<footnote>This is an open access article under the terms of the Creative Commons Attribution-Non Commercial License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited and is not used for commercial purposes.</footnote>

<footnote>© 2024 The Author(s) Laser and Photonics Reviews published by Wiley Periodicals LLC on behalf of American Institute of Chemical Engineers.</footnote>
is a periodic meta-atom with its associated phase shift and reflection or transmission amplitude, which are computed with a fullwave solver such as Rigorous Coupled-Wave Analysis (RCWA) or Finite Difference Time-Domain (FDTD) method. The main limitation is related to the reduced accuracy when neighboring meta-atoms interact strongly and when the phase changes rapidly. This is typically the case for nonlocal resonant metasurfaces and for large Numerical Aperture (NA) metalenses associated with oblique incidence [15].

On the other hand, the common element between the different inverse design strategies is the need to perform fullwave electromagnetic simulations for the entire structure. This can be computationally expensive, especially when designing large-scale, e.g., centimeter-scale, metalens configurations. Topology Optimization (TO) is more suitable for high-dimensional, continuously defined design spaces because it can take advantage of gradient information to optimize efficiently. Nevertheless, TO is not a global optimization method, being vulnerable to bad local minima. Furthermore, TO generally relies on the adjoint state formalism, which requires solving both the forward and adjoint electromagnetic fields at each optimization step, therefore TO can become memory- and computation-intensive. Inverse design strategies that leverage evolutionary algorithms do not require access to the gradient of the objective function. This allows them to handle discrete parameters and avoid local minima. However, they scale poorly as the number of design variables grows, as the number of simulations needed rises exponentially. Bayesian Optimization (BO) is attractive in low-dimensional situations because, by building surrogate models, it can, with a reasonable budget of fullwave simulations, reach satisfying solutions. But, as the parameter space increases in size, BO faces issues in updating and storing the surrogate model and thus becomes impractical in large-scale situations.

Recently, an alternative design strategy based on the individual design of multiple metad deflectors that direct incident light to the focus has been proposed [16, 14]. One major drawback of this approach is related to the necessity of restarting the optimization process for each individual metagrating. Thus, this approach becomes computationally expensive for the design of large-scale metalenses [17]. Moreover, the final metalens design often shows rapid spatial variations, which can violate the periodic assumption used in simulating individual meta-atoms, reducing overall efficiency.

In order to overcome the previous limitations, in this study, we build on the latter approach, and we demonstrate an innovative and efficient deep learning inverse design methodology to construct highly efficient metalenses working in reflection with oblique incidence adapted to AR/VR applications. More precisely, we develop an advanced architecture to deliver highly optimized metagratings that deflect light at arbitrary different angles. Furthermore, our innovative approach allows us to intelligently select adjacent metagratings, providing a smooth spatial variation of the parameters while controlling the phase to achieve as close as approximation as possible to a guassian beam [18, 19, 7, 20]. Furthermore, our method optimizes the entire metalens directly, accounting for the strong near-field coupling between individual metad deflectors and nonlocal effects. In addition, our design strategy offers the flexibility to create structures of varying sizes and focal lengths without requiring a new optimization for each case.

<figure><img src="image_1.png" /><figcaption>**FIGURE 1** | (a) A schematic representation of the considered metalens composed of several beam deflectors, where each metadeflector is deflecting light to a given angle, illustrated by the red arrows. (b) An illustration of super-cell grating (beam deflector) where H denotes the height of the TiO₂ pillars, T is the thickness of the SiO₂ layer separating the pillars from the metal, S is the distance between each pillar and W is the width of each pillar. The parameters ranges considered in the simulations are: 50 nm ≤ H ≤ 200 nm, 50 nm ≤ T ≤ 150 nm and 50 nm ≤ S, W ≤ 1150 nm. A detailed discussion related to the choice of those ranges can be found in the supplementary information section A. The 2D simulations are conducted using RCWA S4 library [24].</figcaption></figure>

Our design methodology process consists of three steps. Firstly, we train a purely data-driven deep surrogate model, serving as a fast differentiable solver that is the basis for the optimization. Then, to address the multi-valued nature of inverse design, we employ a Multi-Valued Artificial Neural Network (MVANN), which is a generative artificial intelligence architecture consisting of a multi-layer perceptron with multiple outputs [7]. In the final stage, we condense all the predictions of the MVANN into a single prediction and use a tandem neural network to refine the solution, ensuring that the final structures respect both the underlying physics and fabrication constraints [20, 21, 22]. A preliminary fabrication validation of the design is conducted by electron beam lithography (EBL) and is presented in the supplementary information.

## 2 | Methodology

### 2.1 | Geometry and surrogate model

We consider a metalens working in reflection composed of several beam deflectors to efficiently direct the incident light towards the focus as it is shown in **Figure 1**(a). Each metagrating consists of two ridges with high refractive index TiO₂. The ridges are separated from the substrate by two layers made of SiO₂ (thickness is one of the optimized parameters) and Ag (100 nm), from top to bottom, respectively, see **Figure 1**(b). The light is injected at a 45° angle from the top with a P-polarization (i.e. the magnetic field is parallel to the y-axis) with a wavelength $λ = 633$ nm. The motivation for such a geometry is based on its high efficiency of phase and amplitude manipulation based on an asymmetric Gires–Tournois cavity [23].

The initial step in our methodology involves the creation of a dataset ($10^5$ simulations) to train a Numerical Network (NN) surrogate model. The fullwave simulations are divided into two sets. The first set comprises 25 000 simulations with $S_1 = S_2$ and $W_1 = W_2$, considering the behavior of a single resonator and targeting high diffraction angles. The second set comprises the
remaining 75 000 simulations that are conducted without the previously-mentioned constraints to target smaller deflection angles. The geometrical parameters are sampled from the Halton sequence with Owen scrambling to ensure a uniform distribution of the design space [25, 26]. The sampled designs are simulated and the diffraction efficiency for each reflected diffraction order is recorded.

The training has been performed by considering either the discrete representation of each mode or using a continuous angular format. In the latter case (our choice), the discrete modes 1 are post-processed to a continuous angular representation where we represent each mode as a Gaussian with non-zero angular Full Width at Half Maximum (FWHM). Once the dataset is defined, the next step is the training of a fully differentiable data-driven surrogate model as illustrated in Figure 2 where the NN maps the values of the geometric parameters and the angle reflected to the reflection efficiency. A typical example is provided on the far right of Figure 2. Being a fast differentiable solver, the surrogate model serves as the foundation of our inverse design strategy. In particular, it is used to aid in the training of subsequent NNs. Therefore, its performance dictates the quality of the results obtained in the inverse design. We refer to the supplementary information section B for additional details on the training and performance of this surrogate model.

## 2.2 | Inverse design strategy

The next important goal is to develop an approach that allows to design efficient large-scale metalenses by gathering different optimized metagratings that deflect the light at different angles. Therefore, as a first step, in this section, we develop an advanced inverse design methodology in which, given a desired deflection angle, we are able to obtain the corresponding geometrical parameters of candidate metagratings that maximize the deflection efficiency in that direction.

Yet, owing to the one-to-many nature of the inverse design, a specific NN needs to be applied. In our case, we employ a MVANN, a NN with multiple solutions designed to provide a robust solution to the one-to-many problem. The final fundamental step is to rely on the physical knowledge of the problem in order to organize the optimized metagratings while maximizing the overall focusing performance. In total, to achieve this goal, we deploy two different NNs, all trained with the aid of the surrogate model presented in Figure 2 and in the Supplementary Information section B.

### 2.2.1 | Multi-Valued Artificial Neural Network

The MVANN is a NN with multiple outputs to comply with the multi-valued nature of the inverse design process. Its input (x) consists of the angular reflection profile represented by a 150-dimensional vector, with each dimension representing the reflection at predefined angles from -90 deg to 90 deg, see the far most right of Figure 2 for an example. In other words, we decompose the angular response into 150 different angles and input them simultaneously to the MVANN. In our case, we chose the output of the MVANN as $N_{outs} = 25$. Each output i of the

MVANN (denoted as MVANN$_i$ (x)) gives the geometrical parameters of the corresponding optimized ith metadeflector, see the Supplementary Information section C for more details. The most important ingredient of the MVANN is its loss function ($L_{MVANN}$), which is here defined as

$$L_{MVANN}(\mathbf{x}, \mathbf{y}) = \frac{1}{\sum_{i=1}^{N_{outs}} \frac{1}{\text{MSE}_i(\mathbf{x}, \mathbf{y})}} \quad (1)$$

In general, the loss function is evaluated based on the Mean Squared Error (MSE) loss between each output $y_i \equiv MVANN_i(\mathbf{x})$ and any given defined response $\mathbf{y}$.

Two different training steps are employed in order to improve the performance and to ensure the convergence. The initial step involves training with a dataset generated by the surrogate model, as indicated in Figure 3 (a). At each epoch, in the first step, we employ mini-batch training, meaning that the weights of the MVANN are updated multiple times before passing to the second step to provide a better guided solution to increase the training efficiency. Therefore, at the end of the first step, the model is able to provide reasonable predictions. The second step given in Figure 3 (b) takes the targets (which will be used to build the final metalens) to the MVANN given the updated weights from the first step. The output of the second step is then passed to the surrogate, which predicts the angular responses of the 25 designs proposed by the MVANN. Then, the MVANN weights are updated to minimize the difference between the angular response predicted by the surrogate model and the target angular response. The second training step fine-tunes the performance to guaranty high-performance design responses, especially for the solutions that are not present in the training dataset [28]. It is worth mentioning that the same NN model is used in both steps. When the first step finishes, the weights are updated, and the new model with updated weights is used in the second training step. Moreover, both steps describe a single training epoch, and they are repeated until convergence, as illustrated in Figure 3. We refer to the Supplementary Information section C for additional information on the architecture and training of the MVANN. After training the MVANN, the required optimized metagrating geometries can be directly retrieved for any desired deflection angle.

## 2.3 | Physics-aware design of arbitrary large-scale metalenses

We remind the reader that our goal is, based on the outputs of the MVANN, to gather and organize the distribution of the optimized metagratings ot achieve highly efficient metalenses with different sizes and arbitrary focusing properties. Although MVANN is a solid NN-based inverse design methodology, four key challenges arise that hinder high-performance design when constructing the metalens from our optimized metagratings.

The first challenge is to ensure that our overall NN model provides the correct metagrating super-cell period required to direct light to the desired angle across the metalens. In our case, the period of the super-cell is defined by the sum of the geometrical optimization parameters $S_1, S_2, W_1$, and $W_2$. Thus, the MVANN might suggest metagratings with slight deviations from the ideal period required to deflect light at the correct angle. The accumulation of these small angular discrepancies is sufficient to cause a

<footnote n="1">¹ In a periodic array, the deflection angles are discrete</footnote>
<figure><img src="image_2.png" /><figcaption>**FIGURE 2** | Schematic representation of the surrogate model. The considered activation function is sine for the hidden layers, while the last layer employs sigmoid [27]. Furthermore, a comparative analysis between the surrogate model and simulation results are presented in the far most right figure. The given comparison is obtained from the test dataset as follows: H = 115 nm, T = 55 nm, S₁ = 500 nm, W₁ = 577 nm, S₂ = 570 nm, W₂ = 244 nm.</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**FIGURE 3** | Schematic of the training of the MVANN, where (a) and (b) represent the first and second training steps at each epoch, respectively. At each epoch, in the first step, we employ mini-batch training, meaning that the weights of the MVANN are updated multiple times before passing to the second step to provide a better guided solution to increase the training efficiency. In the second step, we pass the target objectives as input and obtain the estimated parameters. These parameters are then passed to the surrogate model to estimate the losses and update the weights of the MVANN. We refer to the text and the **section C** in the Supplementary Information file for more details.</figcaption></figure>

significant spatial error at the focal point, leading to a significant drop in performance.

The second challenge arises from the fact that the fabrication of different film thicknesses on a single wafer is impractical. In our model, parameters such as H and T are also optimized, and therefore the optimal designs provided by the MVANN might yield gratings with different H and T. However, in practice they must remain uniform across the whole metalens. A commonly used strategy is to fix H and T in advance based on basic physical intuition [17, 14]. However, the best scenario would be to optimize all parameters (including H and T) simultaneously, obtaining the optimal uniform film thicknesses together with the widths of the resonators to maximize the efficiency of the metalens. Nonetheless, imposing uniformity of H and T within different evaluations of the MVANN is not straightforward.

The third challenge is related to the discontinuity of the parameters between neighboring metagratings. The simulations used to optimize the metadeflectors were performed with periodic boundary conditions applied to at the super-cell level, implying an infinite repeating pattern. However, the metalens has a finite size, therefore this concept is not totally valid; moreover, in our design approach we select one metagrating super-cell for each desired angle. Consequently, to maximize the deflection efficiency and respect our assumption, it is essential to ensure a smooth parameter variation between neighboring deflectors.

The last challenge is related to the phase matching of the metadeflectors. Our ideal is to utilize different order modes to deflect incident light to the focus. The -1 order mode is more suitable for the deflection at positive angles, while -2 order mode is appropriate for deflection at negative angles. This creates a discontinuity
in the size of the metadeflectors and, more importantly, a discontinuity in the phase shift that needs to be taken into consideration. In order to solve the above issues and to design highly efficient large-scale metalenses, we perform a post-processing of the output solutions of the MVANN.

The first step is to define the characteristics of the desired metalens such as diameter and focal position. Then, we start our post-processing by initializing the position and period of each metadeflector that will be used to direct the light to different angles. It is worth mentioning that defining the exact position of the metadeflectors is not straightforward. Therefore, an iterative algorithm has been developed to precisely define positions and the periods (see the Supplementary Information section D). The position and the period are considered as hard constraints during the construction of the metalens.

The second step in our post-processing involves condensing the predictions of the MVANN into a single solution to apply the constraints as uniform H, T, and the correct required period, to get a good initial starting point. Thereafter, a tandem NN is deployed to refine the solutions (initial guess from the previous step) and apply a smooth transition of geometrical parameters across the metalens to obtain robust optimal designs. Finally, an innovative framework has been derived to compensate for the phase shift of the deflectors between the -1 and -2 zones to maximize focusing performance.

### 2.3.1 | Condensing the MVANN solutions

The main objective of condensing the predictions of the MVANN is to obtain a single optimized solution per target, considering the hard constraints defined in the previous section. To condense the solutions, we developed a function whose input is all the solutions provided by MVANN for all targets (25 per target), while its output is a single set of parameters per target. The output of the condenser function will then be used as a starting point to obtain the optimized metagrating distributions for the given metasurface structure. The condenser function is configured to combine the best solutions provided by MVANN while filtering the worst ones. For each target, a given weight is assigned to each output solution of MVANN. This weight is calculated from the MSE between the target response and the predicted response from the MVANN. We refer to the Supplementary Information section E for a detailed explanation of this process. This procedure ensures that only promising designs are retained.

The condenser enforces uniform values of H and T across all targets by averaging the respective values for all targets.

Furthermore, hard constraints are imposed on the geometrical parameters $S_1, W_1, S_2$, and $W_2$ for the period $P = S_1+W_1+S_2+W_2$ to be the exact one to direct light to the desired target angle.

The solutions outputted by the condenser will serve as an initial guess for the following step.

### 2.3.2 | Tandem NN with pre-training

Condensing the solutions provided by the MVANN leads to a discontinuous parameter distribution across the metalens. The final step in the methodology is to provide a gradual parameter variation across the metalens, and this is achieved by using a tandem NN, enforcing a continuous representation of the geometrical parameters.

The motivation behind such requirement stems from the non-local nature of the deflectors, and as simulations are performed assuming an infinite array of adjacent structures.

Therefore, to design a metalens it is essential to ensure the periodicity approximation by avoiding abrupt changes in the parameter space.

To address this issue, we introduce a tandem NN for inverse design. In this framework, the input is a single variable representing the position $x$ of the center of each meta-deflector in the metalens, while the outputs are the corresponding geometrical parameters [29]. This setup enables the tandem network to learn a continuous mapping between position and geometry, effectively regulating how quickly the parameters can vary across the lens. As a result, neighboring metagratings exhibit smooth parameter transitions, which preserves the validity of the periodic approximation used during training and leads to improved overall performance.

We use the suggested parameters provided by condensing the MVANN solutions as an initial guess for the tandem NN.

In other words, before training the tandem NN with the surrogate model<fnref n="2" />, we pre-train it with the output of the condenser, which provides a good first guess and speeds up convergence. This strategy allows us to initialize the network near favorable local minima, resulting in consistently high-performance devices. Additional details on the training and architecture of the tandem NN can be found in the Supplementary Information section H.

The proposed tandem NN presented in Figure 4 exploits hard constraints to obtain a smooth parameter variation between neighboring deflectors. In other words, we deploy a trainable NN with one input indicating the metadeflector center $x$ and three outputs. These output parameters, together with the desired metagrating period, are used to determine the geometrical parameters $S_1, W_1, S_2$ and $W_2$. The first output parameter indicates the occupation ratio of the leftmost deflector ($FF_{main}$), while the other two parameters represent the filling fractions of the left ($FF_1$) and right deflectors ($FF_2$) as shown in the right part of Figure 4 (a). The relation between the geometrical parameters and the three-dimensional output vector is given in Figure 4 (a).

Afterward, to enforce uniform values of H and T across the entire metalens, we define them as global learnable parameters, that are independent of the input position. In practice, once the network is trained, the same values of H and T are applied consistently to all input positions $x$.

The training process is illustrated in Figure 4 (b) where, as an example, we demonstrate the optimization of a metalens with 20 µm of diameter, with $x$-axis representing the position of the metadeflectors in the metalens and $y$ the value of the parameter $S_1$ across the metalens. We considered metagratings that have been optimized to maximize deflection in the -1 order mode for the negative angles (left part of the lens) and optimized metagratings that deflect in the -2 order mode for the positive angles (right part).

Figure 4b(i) exemplifies the training of the condenser, showing the parameter $S_1$ obtained after condensing (blue dots) and the solution provided by a randomly initialized tandem NN.

Figure 4b(ii) displays the next step, where we pre-train the tandem NN to align with the solutions provided by the condensed MVANN outputs.

<footnote n="2">²A tandem NN is usually trained by having its outputs processed by a surrogate model, and optimizing directly the objective function.</footnote>
<figure><img src="image_4.png" /><figcaption>FIGURE 4 | Schematic of the tandem NN, depicting its (a) schematic and (b) training. It is a NN that has the position x in µm as input and the parameters (H, T, S₁, W₁, S₂ and W₂) as output. The training process illustrated in (b). See the text for more details.</figcaption></figure>

Finally, Figure 4b(iii), shows the final solution provided by the tandem NN after the training, leveraging the surrogate model to maximize the deflection efficiency for each target angle.

As can be seen by Figure 4b(iii), the final parameters differs from the pre-trained parameters, which indicates that the condensed solution is not optimal, and our tandem NN refines the solution, converging to a more efficient design.

A parameter discontinuity is observed near the origin, which can be explained by the abrupt variation of the super-cell period. As, this region is the transition between the -1 and -2 diffraction order.

### 2.3.3 | Phase matching across metad deflectors

The metalens is sectioned into the most efficient order modes to deflect to the focus.

In our case, it is divided into 2 sections, -1 order modes for lower deflection angles and -2 order modes for higher.

The decision of where this transition is located is made based on a study detailed in the Supplementary Information, Section F.

However, the phase of each region must be matched to properly calculate the efficiency of the metalens.

Nevertheless, we only optimize for efficiency, without controlling the phase.

To ensure phase matching between both sections, we introduce a separation between them.

Equation 2 represents the phase at the focus for the two metad deflectors located at the boundary between the -1 and -2 order regions.

$$ \phi_1(x) = k\sqrt{x_{cl}^2 + f^2} = \phi_2(x) = k\sqrt{x_{cr}^2 + f^2} + \phi + k_x d + 2n\pi \quad (2) $$

Where $\phi$ represents the phase shift between both regions, which is calculated by simulating the output phase of each corresponding mode and averaging their complex amplitudes.

Then, the absolute phase difference is given by the complex angle between both averaged responses.

$n$ is the diffraction order, $d = \frac{P_l + P_r}{2} + S$ is the distance between the centers of the deflectors at the interface, $P_l$ and $P_r$ represents the periods of the left and right deflectors, respectively.

$S$ is the separation between the deflectors, $f$ is the focal distance and $x_{c,l}$ and $x_{c,r}$ correspond to the center positions of the leftmost and rightmost deflectors, respectively.

Figure G8 in the Supplementary Information section G represents a schematic representation of the path taken by incident light.

Finally, we want to find the separation $S$ such as $\phi_1(x) = \phi_2(x)$. To do so, we solve the quadratic equation represented by Equation 3 for $d$.

A detailed examination of Equation 3 can be found in the Supplementary Information section G.

$$ d^2(k^2 - k_x^2) + d(2k^2 x_{cl} + 2kk_x \sqrt{x_{cl}^2 + f^2} - 2k_x \phi + 4k_x n\pi) + \\ + 2(\phi - 2n\pi)k\sqrt{x_{cl}^2 + f^2 - \phi^2 + 4n\pi\phi - 4n^2\pi^2} = 0 \quad (3) $$

## 3 | Numerical results

In this section, we numerically demonstrate our optimization methodology by comparing it to state-of-the art methods. The first demonstrated example consists of a 50 µm lens with a NA of 0.3. Figure 5 (a) demonstrates that the lens shows focusing at the intended focal distance, while Figure 5 (b) shows the power flow profile at the focal plane for the ideal lens and the designed lens. The efficiency of the metalens is calculated based on the fraction of the input power and the power at the focal plane, in this case, it is observed an efficiency higher than 60%, with the remaining power being either scattered due to diffraction or absorbed by the silver layer.

To obtain the power flow distributions, we simulated the near-field responses with Finite Element Method (FEM) (COMSOL Multiphysics) and propagated them using the Rayleigh-Sommerfeld diffraction formula [30]. The details of the derivation
<figure><img src="image_5.png" /></figure>

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

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

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

<figure><img src="image_9.png" /><figcaption>**FIGURE 5** | Numerical demonstration of a 50 µm diameter metalens operating in reflection with NA = 0.3. (a) and (b) show the normalized z component of the power flow: (a) corresponds to a centralized focus, and the dashed line indicates the plane where the cross-section in (b) is displayed. The profiles are normalized to the peak flux of the ideal focal spot, which is obtained by propagating the ideal parabolic phase profile. (c) and (d) illustrate the capability of our method to focus at different focal positions, achieved by appropriately adjusting the target angles, grating pitches, and separation to account for the new lens center.</figcaption></figure>

of the 2D derivation for the Rayleigh-Sommerfeld diffraction formula are provided in the Supplementary Information section K. Figure 5(c–d) highlights the flexibility of our methodology by showing that it is capable of designing lenses with different focal positions in the x-axis.

One advantage that our methodology offers is the adaptability to different scenarios, as classical optimization algorithms would require a complete re-initialization of the optimization for each focal position, while our method requires only the training of the tandem NN which lasts seconds.

The focusing efficiencies achieved are 39% for (c) and 44% for (d). **Figure 6** illustrates the efficiencies of metalenses with different NA and compared with state-of-the-art methods. The observed results for our method (blue curves) shows invariance across different metalens sizes.

We compared our method with LPA (green curves; see Supplementary Information Section I for details).

It is observed that our design consistently outperforms LPA across all NAs and diameters.

We also tested a classical tandem NN (without the pre-training phase), shows by the orange curves in Figure 6.

Our methodology (blue curves) provides consistently the best possible results, shown by the reduced size of the error bar,

highlighting the importance of the pre-training of the tandem NN.

One important remark is that our method is not limited by size, with computational cost being invariant with metalens size.

In the supplementary information section M we demonstrate the fabrication and a preliminary attempt of characterization of a 700 µm lens.

## 4 | Discussion and conclusion

We developed an innovative inverse design methodology to optimize arbitrary large-scale metalenses. We incorporated physical knowledge of gratings and metasurfaces into a deep learning methodology that enables the inverse design of highly efficient large-scale metalenses based on an ensemble of metad deflectors. Traditional optimization approaches, such as topology optimized designs, require the simulation of the entire metalens. Such method is able to provide the best efficiencies, however, becomes computationally prohibitive for large, centimeter-scale metalenses. Our method overcomes this limitation by incorporating deep learning with physics-aware design principles. This
<figure><img src="image_10.png" /><figcaption>**FIGURE 6** | Comparison between our method and state-of-the-art methods (LPA and tandem NN without pre-training). Each graph represent the average and standard variation of 10 independent optimizations. For each NA, the optimal proportion of -1 and -2 order elements was selected, as detailed in Supplementary Information (Figure F7).</figcaption></figure>

approaches uses the concepts of coupled systems while ensuring that the approximation of infinite repeating systems remains valid.

Recently, novel approaches have been developed where each diffractive element is optimized individually and then assembled manually to construct the metalens [17, 14]. Although, these approaches achieve high diffraction efficiencies at the level of individual diffractive elements, they do not inherently account for coupling between adjacent elements when constructing the metalens. Moreover, the manual assembly process can be computationally expensive, particularly for large-scale applications. In contrast, our methodology naturally incorporates element coupling and phase matching, enabling the rapid design of high-performance metalenses with varied specifications without the need to re-optimize the entire structure.

Additionally, we performed a preliminary experimental characterization of a 700 µm lens present in the Supplementary Information. We demonstrated that the fabricated metalens focal distance matches the designed focal distance. In the follow-up study, the characterization process will be further developed to obtain a more robust characterization of the metalenses.

## Author Contributions

A.C.D.S. developed the optimization framework and performed the numerical simulations. M.E. proposed the metalens structure. A.C.D.S. and M.E. wrote the manuscript and discussed the physical results. P.E.d.S.P and S.V.G.N fabricated and characterized the structure and discussed the fabrication with M.A. A.C.D.S, S.L., H.E.H.F. and M.E. discussed the numerical simulations. All the authors reviewed the manuscript.

## Acknowledgments

The authors acknowledge the fund from the the Association Nationale de la Recherche et de la Technologie under the CIFRE contract 2022/1347. This research used facilities of the Brazilian Nanotechnology National Laboratory (LNNano), part of the

Brazilian Centre for Research in Energy and Materials (CN-PEM), a private non-profit organization under the supervision of the Brazilian Ministry for Science, Technology, and Innovations (MCTI). The Nanofab staff is acknowledged for the assistance during the experiments MNF-20250217. This research was funded by the National Council for Scientific and Technological Development (CNPq), project 314539/2023-9 (HEHF), and by the São Paulo Research Foundation (FAPESP), projects 2021/11380-5 (CPTEn), 2021/00199-8 (SMARTNESS), 2022/11596-0 (EMU).

## Financial Disclosure

None reported.

## Conflicts of Interest

The authors declare no conflicts of interest.

## References

1. Steven Vi, Tiago Silva da Silva, and Frank Maurer. 2019. "User experience guidelines for designing HMD extended reality applications." In *Human-Computer Interaction-INTERACT 2019: 17th IFIP TC 13 International Conference, Paphos, Cyprus, September 2-6, 2019, Proceedings, Part IV* 17, 319-341. Springer.

2. Younghwan Yang, Junhwa Seong, Minseok Choi, Junkyeong Park, Gyeongtae Kim, Hongyoon Kim, Junhyeon Jeong, Chunghwan Jung, Joohoon Kim, Gyoseon Jeon, et al. "Integrated metasurfaces for re-envisioning a near-future disruptive optical platform." *Light: Science & Applications* 12, no. 1 (2023): 152.

3. Yan Li, Xiaojin Huang, Shuxin Liu, Haowen Liang, Yuye Ling, and Yikai Su. "Metasurfaces for near-eye display applications." *Opto-Electronic Science* 2, no. 8 (2023): 230025-1.

4. Jie Hu, Sankhyabrata Bandyopadhyay, Yu-hui Liu, and Li-yang Shao. "A review on metasurface: from principle to smart metadevices." *Frontiers in Physics* 8 (2021): 586087.

5. Yingli Ha, Yu Luo, Mingbo Pu, Fei Zhang, Qiong He, Jinjin Jin, Mingfeng Xu, Yinghui Guo, Xiaogang Li, Xiong Li, et al. "Physics-data-driven intelligent optimization for large-aperture metalenses." *Opto-Electronic Advances* 6, no. 11 (2023): 230133-1.
6. Zhi-Bin Fan, Yun-Fan Cheng, Ze-Ming Chen, Xia Liu, Wen-
Long Lu, Shi-Hao Li, Shao-Ji Jiang, Zong Qin,, and Jian-Wen
Dong. "Integral imaging near-eye 3D display using a nanoim-
print metalens array." *eLight* 4, no. 1 (2024): 3.

7. Arthur Clini de Souza, Stéphane Lanteri, Hugo Enrique Hernandez-Figueroa, Marco Abbarchi, David Grosso, Badre Kerzabi,, and Mahmoud Elsawy. "Back-propagation optimization and multi-valued artificial neural networks for highly vivid structural color filter metasurfaces." *Scientific Reports* 13, no. 1 (2023): 21352.

8. Egor Khaidarov, Damien Eschimese, Keng Heng Lai, Ai-hong Huang, Yuan Hsing Fu, Qunying Lin, Ramon Paniagua-Dominguez,, and Arseniy I Kuznetsov. "Large-scale vivid metasurface color printing using advanced 12-in. immersion photolithography." *Scientific Reports* 12, no. 1 (2022): 14044.

9. Jing Cheng Zhang, Yubin Fan, Jin Yao, Mu Ku Chen, Shirong Lin, Yao Liang, Borui Leng,, and Din Ping Tsai. "Programmable optical meta-holograms." *Nanophotonics* 13, no. 8 (2024): 1201-1217.

10. He-Xiu Xu, Guangwei Hu, Yanzhao Wang, Chaohui Wang, Mingzhao Wang, Shaojie Wang, Yongjun Huang, Patrice Genevet, Wei Huang,, and Cheng-Wei Qiu. "Polarization-insensitive 3D conformal-skin metasurface cloak." *Light: science & applications* 10, no. 1 (2021): 75.

11. Mahmoud MR Elsawy, Stéphane Lanteri, Régis Duvigneau, Jonathan A Fan,, and Patrice Genevet. "Numerical optimization methods for metasurfaces." *Laser & Photonics Reviews* 14, no. 10 (2020): 1900445.

12. Lingcong Luo, Zhengyang Wang, Juntao Li,, and Haowen Liang. 2023. "Wide-field-of-view trans-reflective RGB-achromatic met-als for augmented reality." In *Photonics*, Vol 10, 590. MDPI.

13. Zhicheng Wu, Ming Zhou, Erfan Khoram, Boyuan Liu,, and Zongfu Yu. "Neuromorphic metasurface." *Photonics Research* 8, no. 1 (2019): 46-50.

14. P Dainese, L Marra, D Cassara, A Portes, J Oh, J Yang, A Palmieri, JR Rodrigues, AH Dorrah,, and F Capasso. "Shape optimization for high efficiency metasurfaces: theory and implementation.". *arXiv preprint arXiv:2405.03930*.

15. Villads Egede Johansen, Uğur Meriç Gür, Jade Martínez-Llinás, Jesper Fly Hansen, Akbar Samadi, Maria Skak Vestergaard Larsen, Theodor Nielsen, Fredrik Mattinson, Moritz Schmidlin, N Asger Mortensen, et al. "Nanoscale precision brings experimental metalens efficiencies on par with theoretical promises." *Communications Physics* 7, no. 1 (2024): 123.

16. Mahdad Mansouree, Andrew McClung, Sarath Samudrala,, and Amir Arbabi. "Large-scale parametrized metasurface design using adjoint optimization." *Acs Photonics* 8, no. 2 (2021): 455-463.

17. Brandon Born, Sung-Hoon Lee, Jung-Hwan Song, Jeong Yub Lee, Woong Ko,, and Mark L Brongersma. "Off-axis metasurfaces for folded flat optics." *Nature Communications* 14, no. 1 (2023): 5602.

18. Zhaoyi Li, Raphaël Pestourie, Joon-Suh Park, Yao-Wei Huang, Steven G Johnson,, and Federico Capasso. "Inverse design enables large-scale high-performance meta-optics reshaping virtual reality." *Nature communications* 13, no. 1 (2022): 2409.

19. Peter R Wiecha, Arnaud Arbouet, Christian Girard,, and Otto L Muskens. "Deep learning in nano-photonics: inverse design and beyond." *Photonics Research* 9, no. 5 (2021): B182-B200.

20. Arthur Clini de Souza, Stéphane Lanteri, Hugo Enrique Hernandez-Figueroa, Marco Abbarchi, David Grosso, Badre Kerzabi,, and Mahmoud Elsawy. 2024. "Inverse design of bright, dielectric metasurfaces color filters based on back-propagation and multi-valued artificial neural networks." In *Machine Learning in Photonics*, Vol 13017, 62-71. SPIE.

21. Arthur Clini de Souza, Paloma Elias da Silva Pellegrini, Silvia Vaz Guerra Nista, Stéphane Lanteri, Hugo Enrique Hernandez-Figueroa, Marco Abbarchi, Badre Kerzabi,, and Mahmoud Elsawy. 2024. "Cost-Efficient Deep Learning Method for Mass

Inverse Design of Photonic Devices.” In 2024 SBFoton International Optics and Photonics Conference (SBFoton IOPC), 1–3. IEEE.

22. Menglong Luo, and Sang-Shin Lee. "Tandem neural network-assisted inverse design of highly efficient diffractive slanted waveguide grating." *Optics Express* 32, no. 7 (2024): 12587-12600.

23. Mahmoud Elsawy, Christina Kyrou, Elena Mikheeva, Rémi Colom, Jean-yves Duboz, Khosro Zangeneh Kamali, Stéphane Lanteri, Dragomir Neshev,, and Patrice Genevet. "Universal active metasurfaces for ultimate wavefront molding by manipulating the reflection singularities." *Laser & Photonics Reviews* 17, no. 7 (2023): 2200880.

24. Victor Liu, and Shanhui Fan. "$S^4$: A free electromagnetic solver for layered periodic structures." *Computer Physics Communications* 183, no. 10 (2012): 2233-2244.

25. John H Halton. "On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals." *Numerische Mathematik* 2 (1960): 84-90.

26. Art B Owen. "A randomized Halton algorithm in R.". *arXiv preprint arXiv:1706.02808*.

27. Vincent Sitzmann, Julien Martel, Alexander Bergman, David Lindell,, and Gordon Wetzstein. "Implicit neural representations with periodic activation functions." *Advances in neural information processing systems* 33 (2020): 7462-7473.

28. Vignesh Gopakumar, Stanislas Pamela,, and Debasmita Samad- dar. "Loss landscape engineering via data regulation on PINNs." *Machine Learning with Applications* 12 (2023): 100464.

29. Christopher Yeung, Ju-Ming Tsai, Brian King, Benjamin Pham, Julia Liang, David Ho, Mark W Knight,, and Aaswath P Raman. "Designing multiplexed supercell metasurfaces with tandem neural networks." *Nanophotonics* 10 (2021): 1133-1143.

30. "COMSOL Multiphysics Reference Manual, version 6.2.0.415", COMSOL, Inc, www.comsol.com , "COMSOL Multiphysics Reference Manual, version 6.2.0.415", COMSOL, Inc, www.comsol.com.

31. Ilya Loshchilov, and Frank Hutter. "Decoupled weight decay regularization." . *arXiv preprint arXiv:1711.05101*.

32. Kaiming He, Xiangyu Zhang, Shaoqing Ren,, and Jian Sun. 2015. "Delving deep into rectifiers: Surpassing human-level performance on imagenet classification." In *Proceedings of the IEEE international conference on computer vision*, 1026-1034.

33. Xavier Glorot, and Yoshua Bengio. 2010. "Understanding the difficulty of training deep feedforward neural networks." In *Proceedings of the thirteenth international conference on artificial intelligence and statistics*, 249-256. JMLR Workshop and Conference Proceedings.

34. Pan Zhou, Xingyu Xie, Zhouchen Lin,, and Shuicheng Yan. "Towards understanding convergence and generalization of AdamW". *IEEE transactions on pattern analysis and machine intelligence*.

35. Kingma Diederik. "Adam: A method for stochastic optimization.". (No Title).

36. Joseph W Goodman 2005. *Introduction to Fourier optics*, Roberts and Company publishers.

37. Douglas B Murphy, and Michael W Davidson. 2012. *Fundamentals of light microscopy and electronic imaging*, John Wiley & Sons.

38. Ajay Shrestha, and Ausif Mahmood. "Review of deep learning algorithms and architectures." *IEEE access* 7 (2019): 53040-53065.

39. Dan Hendrycks, and Kevin Gimpel. "Gaussian error linear units (gelus)." . *arXiv preprint arXiv:1606.08415*.

40. Laurent Dinh, Razvan Pascanu, Samy Bengio,, and Yoshua Ben-
gio. 2017. "Sharp minima can generalize for deep nets." In *International Conference on Machine Learning*, 1019-1028. PMLR.

41. Yuhan Zheng, Mingfeng Xu, Mingbo Pu, Fei Zhang, Di Sang, Yinghui Guo, Xiong Li, Xiaoliang Ma,, and Xiangang Luo. "Designing high-efficiency extended depth-of-focus metalens via topology-shape optimization." *Nanophotonics* 11, no. 12 (2022): 2967-2975.
42. MG Moharam, and Thomas K Gaylord. "Rigorous coupled-wave analysis of planar-grating diffraction." *Journal of the Optical Society of America* 71, no. 7 (1981): 811-818.

43. Jie Wu, Guangxin Liu, Lingyan Li, Haitao Xu, Yuwei Lu, Jie Tan,, and Jingfeng Liu. "Tandem neural networks for designing nanoplasmonic structures to manipulate quantum dynamics." *Optics Express* 33, no. 5 (2025): 9790-9803.

44. Xumin Ding, Chunsheng Guan, Zhuochao Wang, Shengying Liu, Kuang Zhang, Xuemai Gu, Qun Wu,, and Ming Jin. "Meta-surface for bending the reflected wave under oblique incidence." *IEEE Transactions on Magnetics* 55, no. 11 (2019): 1-4.

45. Mohammad Yazdi, and Mohammad Albooyeh. "Analysis of metasurfaces at oblique incidence." *IEEE Transactions on Antennas and Propagation* 65, no. 5 (2017): 2397-2404.

46. Mohammad Albooyeh, Sergey Kruk, Christoph Menzel, Christian Helgert, Matthias Kroll, Adam Krysinski, Manuel Decker, Dragomir N Neshev, Thomas Pertsch, Christoph Etrich, et al. "Resonant metasurfaces at oblique incidence: interplay of order and disorder." *Scientific reports* 4, no. 1 (2014): 4484.

47. Shuming Wang, Pin Chieh Wu, Vin-Cent Su, Yi-Chieh Lai, Mu-Ku Chen, Hsin Yu Kuo, Bo Han Chen, Yu Han Chen, Tzu-Ting Huang, Jung-Hsi Wang, et al. "A broadband achromatic metals in the visible." *Nature nanotechnology* 13, no. 3 (2018): 227-232.

48. Shijie Xiong, and Xianguang Yang. "Optical color routing enabled by deep learning." *Nanoscale* 16, no. 19 (2024): 9284-9294.

49. Taewon Choi, Chulsoo Choi, Junseo Bang, Youngjin Kim, Hyun-woo Son, Changhyun Kim, Junhyeok Jang, Yoonchan Jeong,, and Byoungho Lee. "Multiwavelength achromatic deflector in the visible using a single-layer freeform metasurface." *Nano Letters* 24, no. 35 (2024): 10980-10986.

50. Raphaël Pestourie, Youssef Mroueh, Chris Rackauckas, Payel Das,, and Steven G Johnson. "Physics-enhanced deep surrogates for partial differential equations." *Nature Machine Intelligence* 5, no. 12 (2023): 1458-1465.

51. Reza Marzban, Hamed Abiri, Raphael Pestourie,, and Ali Adibi. "HiLAB: A Hybrid Inverse-Design Framework." *arXiv preprint arXiv:2505.17491*.

52. Adam Overvig, and Andrea Alù. "Diffractive nonlocal metasurfaces." *Laser & Photonics Reviews* 16, no. 8 (2022): 2100633.

53. Kunal Shastri, and Francesco Monticone. "Nonlocal flat optics." *Nature Photonics* 17, no. 1 (2023): 36-47.

54. Siyang Xiao, Fen Zhao, Dongying Wang, Junjie Weng, Yan Wang, Xin He, Huan Chen, Zhaojian Zhang, Yang Yu, Zhenfu Zhang, et al. "Inverse design of a near-infrared metalens with an extended depth of focus based on double-process genetic algorithm optimization." *Optics Express* 31, no. 5 (2023): 8668-8681.

## Supporting Information

Additional supporting information can be found online in the Supporting Information section.

## APPENDIX

### A | Deflector ranges choice

The parameter space can be split into two parts for the thickness and width of the resonators respectively. For the thickness T of the SiO₂ layer, we chose to span from 50 nm to 150 nm. The minimum value was chosen to provide a minimum separation between the Ag layer and the TiO₂ resonators. The reason behind this choice is to prevent interference between some modes with high field concentration on the edges of the resonators, and the metal layer. For the choice of the maximum thickness, we considered 150 nm to avoid Gires-Tournois excitation that can limit

<figure><img src="image_11.png" /><figcaption>**FIGURE B1 |** Training evolution of the surrogate model.</figcaption></figure>

the control of the coupling. The height H for the TiO₂ resonators ranges from 50 nm to 200 nm. The choice for the upper bound was made based on limitations on nano-imprint lithography, a cost-effective nanofabrication method with the tightest limitation on structures heights. We also confirmed that the chosen maximum height is enough to provide 2π phase manipulation using LPA. Considering the resonators width and spacing, we considered the minimum to be 50 nm due to fabrication limitations while the maximum to be 1150 nm, to maximize the deflection angle. By having a period of 4.6 µm of length, we can deflect light from the -1 order mode to 57 deg. We have observed larger values of period for high NA metalenses however, the tandem NN proved to generalize well despite extrapolating the dataset used for training the surrogate model.

### B | Surrogate model training

**Figure B1** shows the loss for each epoch of the surrogate model, indicating that the validation loss matches the train loss, showing good generalization performance. The optimizer used here is AdamW with batch size of 2880 [31]. We implemented a linear scheduler to reduce the step size linearly at every epoch. The step size starts at $5 \times 10^{-3}$ and ends at $2.5 \times 10^{-7}$ at epoch 150. We observe that by halting the step reduction, we observe a stagnation in the loss evolution. Additionally, we consider patience of 20 and a minimum delta of $10^{-6}$. That is, if the validation loss did not improve more than $10^{-6}$ over 20 epochs, the training will stop. The weight initialization utilized was Kaiming uniform, considering ReLu nonlinearity and fan in mode [32, 27]. After the training of the surrogate model, we demonstrate its accuracy in **Figure B2**, showing train and test performance.

### C | MVANN architecture and training details

The MVANN considered architecture is shown in **Figure C3**. Before the training starts (initial phase before Figure 3 (a)) the weights were initialized with Xavier uniform [33]. The optimizer of choice is AdamW with initial learning rate of $5 \times 10^{-3}$ with exponential decay of learning rate, reducing it 1% per batch [34]. We considered mini-batch training for the first phase (Figure 3 (a)). The batch size was set to 200 for both the dataset comprised of $5.6 \times 10^4$ designs generated using the surrogate model and the number of targets used in the second phase.
The dataset to train the MVANN in the first step was generated using the surrogate model to maintain consistency in the data used to train it in the second step (which uses the surrogate model).
<figure><img src="image_12.png" /><figcaption>**FIGURE B2** | Surrogate model performance after training. Here we map the discrete deflection efficiencies into a continuous pattern. In red, the solution given by the RCWA simulation, in blue, the surrogate model prediction, and in dashed green, the surrogate model prediction after interchanging $S_1$ with $S_2$ and $W_1$ with $W_2$, hereby called the rotated design. The response of the design and its rotated version should be identical, however this behavior is not set as a hard constraint for the surrogate model, and is instead learned.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**FIGURE C3** | Schematic representation of the architecture of the MVANN. The number in parentheses indicating the dimension of each layer.</figcaption></figure>

To determine the number of outputs, a test was performed with different numbers of outputs, shown in **Figure C4**. We considered two types of losses: firstly, the regular MVANN loss; secondly, the minimum MSE loss across all solutions. The regular MVANN loss is differentiable while the minimum MSE is not, meaning that it is not possible to use the latter to train the MVANN. The problem lies in the fact that the regular MVANN loss cannot be used to determine the number of outputs, as this loss naturally decreases by increasing the number of solutions. This effect is caused because we add another number to the denominator of the loss, therefore, reducing the MVANN loss without necessarily finding solutions with lesser MSE loss. This implies that, using the regular MVANN loss to determine the number of outputs render an inconclusive result as seen by Figure C4 represented by the blue and orange curves, where we observe a monotonic descending curve with the increase of $N_{outs}$. Therefore, we proposed an alternative loss to evaluate the number of outputs. By utilizing the minimum MSE loss across all solutions, it becomes possible to evaluate if the solutions themselves are improving, as this measure does not depend on the number of

solutions. Finally, the most important loss to consider is the minimum MSE loss of the targets, as the main purpose of the MVANN in this study is to predict the parameters given the targets, and it can be observed a stagnation with $N_{outs} > 30$. The final choice of 25 outputs is due to its performance being similar to 30 outputs while having less computational burden.

## D | Algorithm for computing the position of the defectors

A measure for normalizing the pitch of the supercell was required due to the nature of the data utilized to train the surrogate model. By utilizing a continuous angular representation for the deflection efficiencies, designs with similar pitches can yield low MSE in their profiles in comparison to the target profiles; however, the deflection angles are not identical, resulting in catastrophic performance for a metalens. **Figure D5** demonstrates the difference between controlling and not controlling the pitch of each element.

Therefore, it is essential to impose hard constraints on certain parameters to assert the best performance, implying that the primary step to conceive the metalens lies in the definition of the position and pitch of each metadeflector to subsequently, impose those values to the tandem NN as hard constraints. **Equation D1** shows the relation between the central position of the resonator with its corresponding pitch

$$P(x^c) = \frac{-m\lambda\sqrt{(x^c - x_f)^2 + f^2}}{\sin\theta_{inc}\sqrt{(x^c - x_f)^2 + f^2 + (x^c - x_f)^2}}, \quad (D1)$$

where $P(x^c)$ is the pitch, $f$ is the focal distance, $x_f$ x-coordinate of the focus, $m$ is the mode order, $\theta_{inc}$ is the incident angle and $x^c$ is the position of the center of the deflector, with $c$ being a notation rather than an index or number. It is important to note
<figure><img src="image_14.png" /><figcaption>**FIGURE C4** | Losses as a function of the number of outputs given by the MVANN. It is firstly compared the regular MVANN loss of the validation dataset with the targets' dataset (second step of Figure 3). Afterward, we considered an alternative loss, where it is calculated the minimum MSE loss from all solutions.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**FIGURE D5** | Example of a solution with and without normalization of the supercell pitch where we predefined the position of the deflectors, calculated the desired angles and performed the inverse design. The difference of the unnormalized case lies on the absence of the normalization of the parameters at the last layer, leaving the choice of the pitch to the NN. Evidently, if the hard constraint of the supercell pitch is not imposed, the deflectors will not point exactly at the focus, resulting in a wider focal spot and probably a deformed image.</figcaption></figure>

that as a meta-atom has a non-zero length, there are two methods to model its position, by using its leftmost boundary ($x^b$, with b also representing a notation), or its center, $x^c$. Additionally, we model the metadeflectors considering a different parametrization, where we translate the resonators by $S_1/2$, centralizing them in the supercell.

The following step involves centralizing the metalens, for which we have developed an iterative method to compute the position and pitch of each metadeflector. This iterative method starts by guessing the position $x_i^b$ of the leftmost boundary of the leftmost resonator.

In our case, we guessed $x_0^b = -0.4N$, $N$ being the number of unit cells used to conceive the metalens. Then, employing Equation D1, the pitch is calculated for the metadeflector at that

position. However, Equation D1 is calibrated to receive the position $x^c$ of the center rather than the left boundary of the supercell $x^b$. Therefore, to find its center, we perform another iterative method (using index $j$), in which we compute the length of the deflector, and from that, estimate the center $x_{i,j}^c$ of the cell using $x_{i,j+1}^c = x_i^b + P(x_{i,j}^c)/2$. This additional iterative method is repeated until $j = n$, where $n$ is the number of iterations. In our case, using $n = 4$, both the center $x_{i,j=n}^c$ and the relative length $P_i = P(x_{i,j=n}^c)$ converged. Then, to know the leftmost position of the following metadeflector $x_{i+1}^b$, we add the position of the current resonator $x_i^b$, its calculated pitch $P(x_{i,j=n}^c)$ and the separation between the resonators $s_i$. This process is repeated for each metadeflector in the metalens. Finally, with the knowledge of the total length for each unit cell, it is possible to assemble them to realize a metalens with focal distance $f$. Nevertheless, the focus
is in general not centered as a result of the guessed position of the leftmost resonator not being correct. To centralize the metalens, we first compute its total length $L = \sum_{i=0}^{N-1} P_i$, to then adjust the position of the leftmost resonator and position it at $-L/2$. We iterate until the center of the metalens is at the origin of the system (within machine precision) or the number of iterations reaches 50. Algorithm 1 summarizes the calculation process for defining the position of the metadeflectors along the metalens. Finally, by knowing the ideal period of each supercell, we use those values to normalize the condenser and the tandem NN.

<caption>**Algorithm 1** Computation of each metadeflector *i* position ($x_i^b$) and pitch ($P_i$).</caption>

```
x_{i=0}^{b} \leftarrow 0.4N \quad \triangleright \text{ The initial guess of the position of the left boundary of the resonator, in this case, } 0.4N.
while Lens is not centered or until 50 iterations do
    for Deflector $i = 0$ until $i = N - 1$ do
        $x_{i,j=0}^{c} \leftarrow x_i^b \quad \triangleright \text{ Initial guess for the position of the center of the deflector}$
        for Iteration $j = 0$ until $j = n$ do \quad \triangleright \text{ In our case, we considered } n = 4.
            $x_{i,j+1}^{c} \leftarrow x_i^b + P(x_{i,j})/2 \quad \triangleright \text{ Using Equation D1}$
            end for
            $x_{i+1}^{b} \leftarrow x_i^b + P(x_{i,j=n}^{c}) + s_i \quad \triangleright \text{ Where } s_i \text{ is the separation between deflector } i \text{ and } i+1.$
            end for
            $L \leftarrow \sum_{i=0}^{N-1} P_i$
            $x_{i=0}^{b} \leftarrow -L/2$
end while
```

## E | Condenser initialization

Equation E2 shows the equation used to compute the corresponding weight of a solution with MSE E.

$$W_{i,j} = \text{Softmax} \left( \frac{E_i - \mu_j}{\Theta \sqrt{\sigma_j^2 + \epsilon}} \right) \quad (E2)$$

where $W_{i,j}$ denotes the weight associated with solution *i* for target *j*, $\mu_j$ and $\sigma_j$ represent the mean and variance of the MSE for target *j*, and the constant $\epsilon = 10^{-5}$ is introduced to prevent division by zero. $\Theta$ is a hyperparameter that controls the smoothness of the weights' distribution; specifically, a lower value of $\Theta$ results in a sharper distribution of weights, emphasizing the best solutions, while higher values of $\Theta$ lead to a uniform distribution of weights. The quality of the solutions found by the condenser are highly sensitive to the selection of this hyperparameter. Here, we determined $\Theta = 2$ after evaluating different values. Figure E6 illustrates the schematic representation for the initialization of weights of the condenser.

## F | Definition of the optimal proportion of elements using the -1 and -2 order mode

While utilizing the strategy of exploiting the best performing modes to deflect light to the focus, it is necessary to define the position of the boundary between elements that exploits the -1 order mode and -2 order mode to deflect light to the focus. This position is chosen based on the proportion of elements in each section rather than the position $x$ on the metasurface, as the latter demonstrated convergence issues of Algorithm 1. The algorithm will diverge in some scenarios as the deflectors in the boundary between using the -1 and -2 order mode would change the sections they belong across iterations, never settling in a side. In our study, we considered the ratio of elements using the -1 order mode. We performed a sweep over the -1 order mode ratio and NAs to calculate their relative efficiencies and determine the optimal ratio for each NA as displayed on Figure F7. It is worth mentioning that we certify that the diameter of the lenses do not vary with the ratio of the -1 order mode, and therefore, being a fair comparison. It is important to note that this study purpose is to maximize the efficiency of the resonators. Nevertheless, such study is computationally expensive and was performed for scientific purposes, and can be substituted by simpler studies without much loss of efficiency, for example, considering a fixed angle of deflection as the limit between orders (for example, positive deflection angles uses -1 and negative, -2, using 0 degrees as the limit), a fixed ratio for the -1 order mode or even, choosing a ratio based on Figure F7. In future studies, we will analyze computationally efficient methods to define where to position the frontier between different modes. Moreover, to inverse design with intermediate values of numerical apertures, an interpolation of the provides a good approximation. However, changing the geometry and materials implies the necessity to perform this study from scratch. For this reason, plan to study methods to analytically calculate the optimal ratio.

## G | Validation of Equation 3

Figure G8 displays the schematic representation of the path taken by the incident light to the focus, leading to the calculation of $\phi_1$ and $\phi_2$ and the derivation of Equation 2.

To validate Equation 3 we performed a separation sweep across different lenses with different NAs. The result is displayed at Figure G9.

We can observe that by choosing the ratio that results in the best efficiencies, it is observed that our method returns the best point on the graph while using $n = 2$. This holds true mostly due to the bias of the choice, as our method accounts for a single effect (zero order Gaussian beam). This ensures that light from the metalens meet at the focus in phase. Other effects might happen that stray out of the classical Gaussian beam, such as by increasing the order or transforming the beam into a Bessel beam, or a combination of different effects. Those effects could be modeled and studied, however, they would unnecessarily increase the complexity of Equation 3 that is precise enough to provide us a highly efficient metalens.

## H | Tandem NN architecture and training

For the tandem NN, we employ an architecture comprising 4 layers, 64 neurons each, with sine activation function on hidden layers and sigmoid for the last layer. Prior to pre-training, the weights are initialized with Pytorch (2.2.2+cu121) default
<figure><img src="image_16.png" /><figcaption>**FIGURE E6** | Diagram of the condenser weight initialization. The condenser is a neural network that condenses all solutions from all targets into a single solution for each target. The process starts with the definition of the targets. Afterward, the solutions generated by the MVANN (25 for each target) are presented to the surrogate model, subsequently calculating the MSE between the target and each solution. Finally, we utilize Equation E2 to find the initial weights of the condenser. Each column represent a different target and each row represent a different output from the MVANN. In other words, horizontally we vary the target and vertically, the output of the MVANN. Indicating that the sum of each column equals 1.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**FIGURE F7** | Results obtained after performing a sweep over the element ratio utilizing the -1 order mode. The considered lens have a 50 µm diameter to reduce computational costs. For each point, 10 solutions were obtained with the method to compute the mean and standard variation of the points to build the error bars. The best chosen percentages are: 67.5 % for 0.2, 57.5 % for 0.3 and 0.4, 47.5 % for 0.5 and 40 % for 0.6. Total time: 07:08:10</figcaption></figure>

weight initialization. The pre-training consists of full batch and 100 epochs of training. In the following step (i.e. the training with the surrogate model), we considered a batch size of 512. For both parts, we used Adam optimizer with learning rate of 3 × 10⁻⁴ and exponential decay with γ = 0.998 [35] (Figure 4). The choice of those parameters was derived from a process of trial and error, since an excessively flexible tandem network defeats the purpose of having a continuous condenser, as it causes an overly rapid spatial variation in parameters, diminishing the deflection efficiencies on the final assembled lens. This principle is applied to a less flexible model, where the parameters do not vary fast enough to match the demands of high-performance deflectors.

## I | LPA simulations

LPA method is based on optimizing subwavelength phase shifters to match a desired phase shift profile. In our case, we considered a super-cell composed of one pillar with fixed film thicknesses (H and T) and varied the width of the pillar to be able to provide 2π phase shift. Therefore, relying on few simulations, it is possible to build a model that maps the width of the ridge to its respective phase shift. By knowing the ideal phase map, it is a matter of asking our model with the desired phase to know the width of the ridges.

For the LPA metalenses, the simulations of the single cell were performed using S4 [24]. The supercell is composed of a single TiO₂ pillar over a layer of SiO₂ followed by a layer of Ag with

<figure><img src="image_18.png" /><figcaption>**FIGURE G8** | Light propagating and interacting with the meta-deflectors at the boundaries between regions that leverages different modes to deflect the incident light to the focus. The red color represents the wavefront of the incident light. The green color represents the path taken by the light that interacts with the deflector that utilizes the -1 order mode, while in yellow, the -2 order mode. We consider a relative average phase shift between both regions to be φ, represented in yellow.</figcaption></figure>

Si substrate. We computed the reflection efficiencies by dividing the 0th order mode power efficiency by the incident power. The phases are computed using the relative argument to the solution with the least filling fraction, as shown by Equation 13.

$$ \phi = \arg \left( \frac{-S_{11}}{S_{11(i)}} \right) \qquad (13) $$

Where $S_{11}$ is the complex mode of the 0th order mode and $S_{11(i)}$ refers to the first initial filling fraction, in our case 30%. By trial and error, we varied the parameters, height (H) of the pillar, pitch (P) of the supercell and thickness (T) of the SiO₂ layer to find a combination that results in 2π phase shift while boasting reasonable parameter values. In our case, we used H = 150 nm, P = 350 nm, T = 100 nm. Obtaining Figure 110.

Then, to build a lens, we employ the parabolic phase profile accounting for the off-axis illumination represented by Equation 14.

$$ -k(\sqrt{x^2 + f^2} - \sqrt{f^2 + (D/2)^2} + x \sin \theta) \qquad (14) $$

Where $k = 2\pi/\lambda$ is the wave vector, $x$ is the element position, $f$, the focal length, $D$, the lens diameter and $\theta$ the incident angle.
<figure><img src="image_19.png" /><figcaption>**FIGURE G9** | Study of different separations between the regions using the -1 and -2 order mode over different NAs. We utilized a 50 µm metalens and considered the highest performing ratios for each NA. Indicating that using n = 2, our choice proved to find the optimal point of efficiency. An interesting exception can be found at the NA 0.6, where we theorize that other effects might be happening, hampering the efficiency of a pure Gaussian mode efficiency.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**FIGURE I10** | Deflection efficiencies to the 0-order mode and phase shifts of a single resonator in a supercell with 350 nm of pitch, 150 nm of resonator height, 100 nm of SiO₂ thickness. Indicating high efficiencies and a 2π phase modulation.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**FIGURE J11** | Output phase of the designated mode for the 50 µm calculated with RCWA.</figcaption></figure>

## J | Phase studies of the meta-deflectors on the final lens

One important aspect that requires attention in this study is the phase matching between deflectors within a region (i.e. utilizing -1 mode or -2 mode). **Figure J11** demonstrates the phase of each section utilized to redirect the incoming light to the focus.

It can be observed that during a continuous section of the metalens, the solution found by our method proved to maintain continuity on the phase, alleviating the necessity to add a separation between every deflector.

**Figure J12** shows how the spatial distribution of the ideal phase map and the one obtained by our design: by manipulating the spatial change in phase rather than the phase itself, we can obtain a much clearer phase distribution compared to LPA, resulting in much higher efficiencies. Due to the abrupt change in the parameter space, we can observe a small discontinuity in the phase near the origin, where the transition from the -1 mode to the -2 order mode occurs.

## K | Derivation of the Rayleigh-Sommerfeld diffraction formula for 2-D simulations

The deduction for the Rayleigh-Sommerfeld diffraction formula is as follows:

We consider the general Rayleigh-Sommerfeld diffraction formula in 3-D considering COMSOL's convention of $E = E_0e^{-i(kx-\omega t)}$ shown by **Equation K5** [30, 36].

$$U_1(\mathbf{r}) = \frac{1}{2\pi} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} U_0(\mathbf{r}') \frac{\partial}{\partial \mathbf{n}} \left( \frac{e^{-ik|\mathbf{r}-\mathbf{r}'|}}{|\mathbf{r}-\mathbf{r}'|} \right) dx'dy' \quad (K5)$$

Where $\mathbf{r} = (x, 0, z)$ is the vector linking a point anywhere in space where we know the fields (in our case, a point above the metalens) with the point where we want to compute the field (e.g. a point in the focal plane), $\lambda = 633$ nm is the wavelength of the light, $d$ is the distance separating the metalens and the point where we want to compute the field and $k$ is the wave vector $2\pi/\lambda$.

Considering $U_0(\mathbf{r}')$ constant on $y$, we can integrate and arrive on Equation K6.

$$\int_{-\infty}^{\infty} \frac{e^{-ik|\mathbf{r}-\mathbf{r}'|}}{|\mathbf{r}-\mathbf{r}'|} dy' = -i\pi H_0^{(2)}(k|\rho - \rho'|) \quad (K6)$$

Where $H_0^{(2)}$ is the 0th order Hankel function of the second kind. And $\rho = (x, 0, z)$.

By computing the differentiation, we arrive on **Equation K7**.
<figure><img src="image_22.png" /><figcaption>**FIGURE J12** | Graph showing the difference between an ideal parabolic phase profile and the phase profile obtained using the z-component of the magnetic field. It can be seen the discontinuity near the origin caused by the discontinuity in the deflector geometrical parameters. The calculated separation is 390 nm, smaller than the pitch of the deflectors on this region. The phase mas was obtained from a 150 µm lens with 0.3 of NA. We collect the fields at a distance of 4λ of the metasurface and calculated the phase using the argument of the z-component of the magnetic field.</figcaption></figure>

$$U(x) = \frac{ikd}{2} \int_{-\infty}^{\infty} \frac{H_{1}^{(2)}(k|\rho - \rho'|)}{|\rho - \rho'|} U_{0}(\rho') dx' \quad (K7)$$

Where *d* is the propagation distance. We can use the asymptotic formula (for large |*z*|) for the Hankel function and substitute *k* and arrive on Equation K8.

$$U(x) = \frac{d}{i\sqrt{\lambda}} e^{\frac{i3\pi}{4}} \int_{-\infty}^{\infty} U_0(x') \frac{e^{-ik|\rho-\rho'|}}{|\rho-\rho'|^{3/2}} dx' \quad (K8)$$

With numerical integration, we can compute the field at any distance *d* from where we collected *U*₀.

The comparison between the utilization of the formula found and COMSOL Multiphysics FEM simulation (full simulation including the focal plane) is displayed in Figure K13. In our implementation, we use PyTorch to run the algorithm in the Graphics Processing Unit (GPU) by using a Simpson's rule numerical integration, rendering simulation times orders of magnitude faster compared to a full-wave FEM simulation. Furthermore, to compute the efficiencies of the metalens it is required to compute the field only at the focal plane and not in the whole simulation domain. By using our formula, we can compute the fields at the desired location without the need to perform a full FEM simulation of the entire domain, proving to be useful for large-scale simulations.

## L | Fabrication

We developed, fabricated, and characterized a 700 µm lens with focal length of, approximately, 860 µm, with a NA of 0.37. To design a practically implementable metalens, it is employed -2 order mode in the entire lens. Additionally, to demonstrate the necessity of the tandem network, the fabricated metalens was obtained by training the condenser instead of the tandem network, obtaining a model with multiple discontinuities on the parameters, as shown in Figure L14.

The metalens was fabricated using a combined process of EBL, sputtering of TiO₂, and liftoff in acetone. To obtain a reflective device, initially, a 100 nm Ag film was deposited on top of a silicon wafer. Then, to act as a buffer layer, 130 nm of SiO₂ was deposited on the sample by plasma enhanced chemical vapor deposition (PECVD). After coating the substrate with the electron-sensitive

<figure><img src="image_23.png" /><figcaption>**FIGURE K13** | Comparison between full-wave simulation of COMSOL Multiphysics FEM for solving Helmholtz equation and Rayleigh-Sommerfeld diffraction formula in 2D given by Equation K8. Showing the accuracy of the approximation.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**FIGURE L14** | Parameter distribution of the fabricated 700 µm lens. The parameters were obtained by training the condenser and not using the tandem neural network, obtaining a discontinuous distribution of parameters.</figcaption></figure>

resist (polymethyl methacrylate (PMMA) AR-P 679), the metalens pattern was directly transferred through EBL, with 10 kV and 200 µC/cm² and developed. Then, we formed a 140 nm-TiO₂ film on top of the patterned sample by sputtering. The sample was immersed in acetone to solubilize the resist and achieve
<table><caption>**TABLE N1** | Simulation and training times. RCWA simulations executed with python library S⁴ on a single thread of a 13th Gen Intel i7-13700H at 4.800GHz. Neural network training performed on a NVIDIA RTX A1000 6GB Laptop GPU. The total system memory is 32 GB.</caption><tr><td>Data generation (10<sup>5</sup> simulations)</td><td>31 min 2 seconds</td></tr><tr><td>LPA data generation (100 simulations)</td><td>2 seconds</td></tr><tr><td>Surrogate model training</td><td>30 min 13 seconds</td></tr><tr><td>MVANN Training</td><td>15 min 55 seconds</td></tr><tr><td>Condenser pre-training (100 epochs)</td><td>200 ms</td></tr><tr><td>Condenser training (4 batches of 200 deflectors)</td><td>40 Epochs per second</td></tr></table>

complete liftoff. Finally, the proposed metalens was ready for characterization.

## M | Characterization

To characterize the fabricated sample, we assembled the setup shown in Figure M15 (a). In short, the 633 nm polarized beam is incident on the metalens at 45° and the normal reflection is acquired with a charge-coupled device (CCD) camera, following a confocal configuration. By translating the objective stage alone while the remainder of the setup remains stationary, it is adjusted the z-distance depicted in Figure M15 (a), hence altering the acquisition distance at which the field intensities are recorded in order to analyze the focal distance of the fabricated device. Figure M15 (b) and (c) shows the characterization results of the fabricated lens, where (b) shows the picture taken at z = 0 µm, representing an optical image of the metalens and (c) illustrates the propagation along the z-axis, obtained by processing images taken by a CCD camera at different z-positions that was obtained by averaging the intensities along the y-axis, and manually aligning the x-position of each array to form the image. Figure M15 (d) and (e) shows the scanning electron microscopy image of the fabricated lens.

We can see, at 860 µm, a concentration of fields indicating a tendency to focus at this distance. While it agrees with the design, we also observe that it corresponds to a small portion of the metalens. This effect can be explained by the long distance (15 cm) between the tube lens and the objective, filtering peripheral waves, and decreasing the maximum wave propagation angle that is registrable by the optical system [37]. Nevertheless, it is still observed focusing at the designed focal distance.

An additional remark is about the discontinuities present in the model, as shown by Figure L14. It can be observed that the spots that divert light away from the focus are regions with parameters discontinuities, indicating that these regions decreases the efficiency of the lens.

## N | Simulation and training time

The graph displaying the training times and simulation times are displayed in Table N1.

## O | Experimental Section

Silicon substrates were organically cleaned in ultrasonic baths with acetone and then, isopropyl alcohol for 5 minutes. Both chemicals were purchased from Sigma-Aldrich. Cleaning of the substrates was followed by an O₂ plasma etch at 60 sccm, for 2 minutes to ensure complete removal of organic impurities. The Ag film formation was carried out via electron beam evaporation (apparatus from Aja Int.). This deposition resulted in 100 nm of

Ag on the silicon substrate, using a rate of 2.5 Å/s to 2.7 Å/s at a base pressure of 9 × 10⁻⁷ Torr and 115 mA. On top of the Ag layer, we used PECVD to form a SiO₂ film with thickness of 130 nm. The deposition rate for the SiO₂ was 3.5 Å/s to 4 Å/s with a current of 25 mA. Once the Ag and the SiO₂ films were properly deposited on silicon, the substrates were ready for pattern transfer using EBL. We utilized PMMA 679-04 (Allresist) as the electron sensitive resist for patterning with EBL. The resist was spin-coated at 4000 rpm, for 40 s and then, baked for 5 minutes at 180 °C. We transferred the metalens optimized design through EBL (Raith eLine Plus) with 10 kV and 200 µC/cm². The design exhibits large area of 0.7 mm × 0.7 mm, hence, it was split into writing fields of 100 µm × 100 µm with step sizes of 20 nm. After developing the patterned sample for one minute (AR-P 600 56, from Allresist), we performed a hard baking for 5 minutes at 180 °C. A film of 140 nm of TiO₂ was deposited via DC magnetron sputtering (Ulvac MCH-9000 system), using a gas mixture of 60% O₂ and Ar at deposition pressure of 10⁻³ Torr. In sequence, we performed liftoff of the sample by immersion in acetone.

For the characterization of the fabricated metalens, we assembled the aforementioned confocal experimental setup using the following components for the incidence of the beam: a 632.8 nm HeNe laser (Thorlabs HNL050L), a λ/2 waveplate (Thorlabs), a collimator (Newport LC-075) and a pinhole of 0.5 mm. For the infinity-corrected optical system, we used a long working distance objective (Mitutoyo Plan Apochromat MY5X-802B), with magnification of five times, a tube lens (Edmund Optics), with focal distance length of 200 mm, and a CCD camera (Point Grey Chameleon CMLN-13S2M-CS).
<figure><img src="image_25.png" /></figure>

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

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

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

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

**FIGURE M15** | (a) Schematic of the experimental setup for the characterization of the metalens, utilizing a 633 nm laser passing through a system for controlling its phase, polarization and coherence. Subsequently, the light is reflected by the metalens and is collected by an infinity-corrected optical system comprised of an objective and a tube lens. (b) The intensity profile obtained for z = 0 µm, resembling an optical image of the metalens. (c) Presents the normalized intensity for each z, demonstrating its capability to focus light at the designed distance of 860 µm. A rapidly varying color map was utilized to percept minute variations in intensity, facilitating the identification of the position and rotation of the metalens. (d) and (e) are scanning electron microscopy image of (d) whole lens, and (e) individual deflectors, indicating the shape of a supercell built with two resonators while optimizing the -2 order mode to deflect to the focal spot.