Advanced Numerical Modeling Methods for the Characterization and Optimization of Metasurfaces

The last 10 years have witnessed an impressive amount of works aiming at the development of thin metamaterials for controlling the wavefront of light, and thus realize planar photonics also referred as flat optics or metaoptics. The concept of metasurface is at the heart of almost all the discoveries in this domain. Metasurfaces are arrays of subwavelength-spaced and optically thin optical elements, which enable new physics and phenomena that are distinctly different from those observed in three-dimensional bulk metamaterials. We present here our recent activities and achievements in relation with the design of metasurfaces, which are concerned with two topics: on one hand, we study numerical characterization approaches that are well suited to the multiscale nature of metasurfaces; on the other hand, we develop inverse design strategies for discovering non-classical metasurface configurations for a target optical functionality. These two topics are addressed in the context of a multidisciplinary collaborative project, which involve computational scientists and physicists. In particular, we apply the proposed numerical methodologies to the design of phase gradient metasurfaces and light front shaping metalenses. In some cases, the numerically designed metasurfaces have been frabricated and experimentally characterized to confirm their predicted performances.


Introduction
During the last decade, metasurfaces have been extensively studied due to their ability to precisely control the phase, amplitude, and wavefront of light. These light-matter interactions are mediated by ensembles of subwavelength metaatoms, made of plasmonic or high dielectric refractive index materials, which have thicknesses within the range of the operating wavelength. The design of metasurfaces is generally achieved using one of the following two approaches. In the first approach, which we will refer to as the direct design approach, a rigorous full wave electromagnetic solver is used to study different classes of meta-atoms. Large sets of key meta-atom parameters are parametrically swept to create meta-atom libraries (also referred as lookup tables), and ensembles of meta-atoms from these libraries are assembled together to create metasurfaces with given de-sired optical responses. This method works well for certain applications. However, it does not incorporate nearfield electromagnetic coupling effects between neighboring meta-atoms and does not generalize to large area, freeform devices. The second approach is kown as inverse design. With this approach, the desired optical response is defined as an objective cost function (for example, the deflection or focusing efficiency) and the inverse problem solves for the shape and dimensions of meta-atoms, and for the metasurface organization, in a manner that maximizes the cost function value. In principle, the inverse design approach can account for near-field interactions by optimizing relatively large metasurface regions at a time. However, it requires rigorous and computationally efficient optimization techniques that can deal with large parameter spaces and multiobjective cost functions. Here we report on our recent efforts to address these two goals by relying on and extending two state-of-the-art numerical bricks: a high order discontinuous finite element method for the simulation of nanoscale light-matter interactions, and a statistical learning-based global optimization technique.

Numerical optimization
The optimization of metasurface performance is a crucial concern as they impact a wide range of applications. Yet, current design techniques are mostly based on engineering knowledge, and may potentially be improved by resorting to inverse design. Recently, several numerical optimization methodologies including both local (gradient-based) and global search methods have been considered to tune the design of metasurfaces according to the desired applications. There are two main classes of optimization methods that have currently been used in inverse design of metasurfaces: gradient-based algorithms and gradient-free approaches.
Gradient-based methods depend on the initial guess of the solution and are efficient in finding local optima, even for a large number of design parameters. Gradient-based method require the knowledge of the derivatives of the cost function with respect to design parameters, which can be evaluated analytically sometimes, or approximated numerically in most cases [3]. These methods have a lengthy history in metasurface design and several descent methods, ranging from steepest-descent to quasi-Newton methods for both constrained and unconstrained problems, have been applied to metasurfaces. The most common methods from the recent the literature include the objective-first and topology optimization algorithms, which have proven to be efficient and rigorous design methodologies [2,13,16,17].
Gradient-free approaches are capable of capturing global optima, albeit for a limited number of parameters, thereby overcoming the local minimum trapping issue of gradientbased algorithms. In other words, with gradient-free global optimization techniques, the final optimized results are not influenced by the initialization of the optimizer. In addition, some of these algorithms can deal with discrete optimization parameters and non-differentiable objective functions, which are conditions that are generally not handled by gradient-based algorithms. Without gradients to help guide the optimization, convergence with global optimization algorithms is often considerably slower than convergence with gradient-based algorithms. To date, several global optimization techniques have been proposed in the context of metasurface design. To properly deal with the extensive (often discrete) parameter space and the existence of several local optima, the majority of inverse design methods of interest for metasurface design are stochastic and include genetic algorithms and evolutionary strategies [10,4].
In addition to the methods discussed above, emergent approaches, including methods based on artificial neural networks and Bayesian optimization methods, have the potential to uncover surprising new metasurface designs. We refer to our review article [8] for a more detailed discusison of the different optimization techniques with their advantages and disadvantages in the field of flat optics.
In our works, we use one of the most advanced optimization techniques based on a statistical learning-based approach, which is known as Efficient Global Optimization (EGO) [11]. The EGO algorithm is a global optimization algorithm that substitutes the complex and costly iterative electromagnetic evaluation process with a simpler and cheaper metamodel. EGO is related to the class of Bayesian optimization. Contrary to the traditional common global optimization strategies like genetic algorithms, EGO is not based on adaptive sampling but on a surrogate model, which is constructed on the basis of available objective function evaluations. This surrogate model utilizes a statistical learning criterion related to the optimization target (usually called merit function) in order to identify which design (set of parameters) should be tested in the next iteration that would provide better results close to the predefined goal.
In general, the EGO is based on two phases. The first one is the Design Of Experiment (DOE), in which an initial database is generated. In essence, a uniform sampling strategy (e.g. Latin Hypercube Sampling) is deployed in order to generate different designs in which the cost function is evaluated using an electromagnetic solver. In the second phase, using the data obtained from the DOE, a Gaussian Process (GP) model, is constructed to fit these data. This GP model allows us to predict the values of the cost function in the parameter space without the need to perform additional electromagnetic simulations. Once this GP model is determined, one can estimate at any point of the design space, the objective function (mean of the GP model) and an uncertainty value (variance of the GP model).The mean and the variance are used together to determine a statistical merit function. In our case, we rely on the expected improvement, which is a function whose maximum defines the next design parameters set to be evaluated. That is to say, in the search parameter space where this function is maximized, we extract the corresponding parameter values, and the corresponding design will be simulated using our electromagnetic solver. Then the database is updated accounting for this new observation (construction of a new GP model based on the updated database). We repeat this process until a predefined convergence criterion is reached, or when the expected improvement is sufficiently small.
We have developed an inverse design methodology that combines two ingredients: a high order DGTD (Discontinuous Galerkin Time-Domain) solver of the time-domain Maxwell equations coupled to a Generalized Dispersion Model [12] for the numerical characterization of a given metasurface design and the EGO algorithm for the optimization of the geometrical characteristics of constituting meta-atoms. Based on this methodology, we can distinguish three main achievements so far: -In [7], we have optimized GaN (Gallium Nitride) semiconductor phase gradient metasurfaces operating at visible wavelengths, and we have clarified the importance of near-field coupling for a simple supercell made of four nanopillars. Our numerical results reveal that rectangular and cylindrical nanopillar arrays can achieve more than respectively 88% and 85% of diffraction efficiency for TM polarization and both TM and TE polarization respectively, using only 150 fullwave simulations. To the best of our knowledge, this is the highest blazed diffraction efficiency reported so far at visible wavelength using such metasurface architectures.
-In [6], we presented for the first time to the metasurface community a multiobjective inverse design approach leveraging the EGO statistical learning-based global optimization method. We combined this approach with our in-house developed 3D high order DGTD fullwave solver in order to optimize 3D achromatic metalens with numerical aperture NA above 0.5.
Our metalens attempts to focus the RGB colours at the same focal plane with the maximum feasible focusing efficiency. In other words, we seek to get the most suitable compromise between the chromatic dispersion and the focusing efficiency for a given fixed focal distance. It is worth mentioning that the optimized designs have been fabricated and experimentally characterized by our physics collaborators at CRHEA.
-Finally, in [5], we addressed the optimization of metasurface designs taking into account uncertainties due to fabrication errors. Despite the advances of nanofabrication facilities developed during the past few years, systematic fabrication errors are usually nonnegligible, and are the source of efficiency degradation compared to modeling and simulation results. In order to discover designs that are robust to these fabrication uncertainties, we again rely on the statistical learningbased methodology that was considered in [6]. Hover, this time, we formulate a multiobjective optimization problem that aims at minimizing two statistical robustness measures defined on the parameterized metasurface geometry. We apply this technique to optimize and design realistic phase gradient metasurfaces in the visible regime. Our preliminary results reveal that we obtain robust designs, which are less sensitivity to the geometrical fabrication errors as compared to the one optimized without considering uncertainties.

Modeling of metasurfaces based on GSTCs
The inverse design methodology discussed in the previous section is able to produce optimal designs with a reasonable amount of fullwave solver calls, however each of these evaluations with our 3D high order DGTD solver remains a costly process when considering complex configurations. This computational overhead is partly mitigated by the ability of the solver to leverage parallel processing but, as a general rule, there is still a need to investigate modeling approaches that allow to substantially reduce the computational effort. To avoid a direct modeling of many small geometric features, homogenisation methods that approximate the complexity of an inhomogeneous material by its effective medium response, i.e., homogeneous artificial material, have been proposed. For the case of metasurfaces, consisting of a surfacic two-dimensional arrangement of nanostructures, equivalent transition conditions linking the values of the macroscopic field quantities on both sides of a thin homogenized layer have been derived, called Generalized Sheet Transition Conditions (GSTCs). GSTCs were originally derived in optics in the 90s in the seminal paper of Idemen [9]. These transition conditions conceal in a tensorial form the equivalent response occurring on reflecting and transmitting fields at complex interfaces. We investigate GSTC modeling and our achivements so far have been the following: -We applied the (surfacic) quasi-periodic homogeneization theory to derive equivalent transmission condition of deeply subwavelength microstructured metasurfaces. Our approach provides a direct link between the coefficients of the susceptibility tensors driving surfacic polarization fields and the structural design of a microstructured metasurface. Moreover, in view of a time-domain implementation of GSTCs, we revisited the classical zero-thickness GSTCs and have shown that it is more convenient to use an enlarged formulation of these conditions in which the original metastructure is replaced by GSTCs that exclude the layer from the physical or computational domain. While negative constant susceptibilities appear in the classical zero-thickness GSTCs, their values in the enlarged formulation are always positive which ensure the stability of the effective problem [14]. The implementation of this GSTC variant in our high order DGTD framework is underway. -We studied theoretically the problem of electromagnetic field transition conditions at conformal interfaces using so called Conformal-Generalized Sheet Transition Conditions (C-GSTCs) to achieve surface topographydependent transmitted and reflected fields. Our analysis, supported by 2D and 3D finite element simulations in the frequency-domain regime, provides a solid theoretical framework to design metasurfaces for cloaking, wearable optics and next generation of freeform imaging systems. An illustration for a cloacking problem in the 2D case is shown in Fig. 3. A paper describing our work and results is under review [15]. In the future, we plan to enable this GSTC-based modeling in our high order HDGFD (Hybridizable Discontinuous Galerkin Frequency-Domain) [1]. Figure 3. 2D FEM simulation (frequency-domain) of a cloaking problem: the cat structure is coated with surfacic susceptibilities so that the field scattered by this structure is the same as the one scattered by the mouse structure.