Fullwave Maxwell inverse design of axisymmetric, tunable, and multi-scale multi-wavelength metalenses

We demonstrate new axisymmetric inverse-design techniques that can solve problems radically different from traditional lenses, including \emph{reconfigurable} lenses (that shift a multi-frequency focal spot in response to refractive-index changes) and \emph{ultra-broadband} multi-wavelength lenses ($\lambda = 1\,\mu$m and $10\,\mu$m). We also present experimental validation for an axisymmetric inverse-designed monochrome lens in the near-infrared fabricated via two-photon polymerization. Axisymmetry allows fullwave Maxwell solvers to be scaled up to structures hundreds or even thousands of wavelengths in diameter before requiring domain-decomposition approximations, while multilayer topology optimization with $\sim 10^5$ degrees of freedom can tackle challenging design problems even when restricted to axisymmetric structures.


I. INTRODUCTION
In this paper, we demonstrate that axisymmetric metalenses can be designed with fullwave Maxwell simulations (as opposed to the scalar-diffraction [1] or domain-decomposition approximations [2] used in prior metasurface designs), for > 100 wavelengths (λ) in diameters, combined with multilayer variable-height topology optimization (TO) with > ∼ 10 4 degrees of freedom (∼ 10 per λ per layer) as shown in Fig. 1. The capability and flexibility of our design approach is demonstrated by solving two challenging new design problems with 10-layer metalenses. First (in Sec. III), we design a multi-scale metalens that simultaneously focuses λ = 1 µm and λ = 10 µm at the same diffraction-limited focal point (numerical aperature NA = 0.85, Strehl ratios 0.84 and 0.60, efficiencies 82% and 95%). Second (in Sec. IV), we design an active metalens that shifts its achromatic multi-wavelength (three λs over a 6% bandwidth) focal spot from NA = 0.7 to 0.8 as the index of the material (GSS4T1 [3]) is changed from n = 3.2 to n = 4.6 (thermally or optically), in contrast to previous work that showed only monochromatic reconfigurability [4]. As a proof of concept, we also show (in Sec. V) an experimental realization of single-layer axisymmetric TO-designed metalens for λ = 1550 nm, fabricated by two-photon polymerization 3D-printing (Nanoscribe Professional GT), demonstrating that such variable-height surfaces are manufacturable. As discussed in Sec. VI, our approach could easily be scaled to much larger diameters and number of layers, and the vast number of design degrees of freedom coupled with the lack of approximations makes it a uniquely attractive method for the most difficult metasurface inverse designs.
Flat-optics metalenses have received widespread attention due to their potential for achieving multiple functionalities within an ultra-compact form factor [5][6][7][8][9]. Prior work on metalens design has largely focused on exploiting local resonant conditions [5] under locally periodic approximation (LPA), using rather small unit cells ( < ∼ λ) [2,10]. Recently, it has been shown that the unit-cell approach with LPA can lead to fundamental limitations on metalens performance [11,12], whereas some of these limitations may be mitigated by using overlapping boundaries, perfectly matched layers or larger domains [13,14] some may not. Meanwhile, axisymmetric multi-level diffractive lenses (MDL) have been proposed as an alternative to metalenses for achieving enhanced functionalities [15,16]; however, MDL designs utilize scalar diffraction theory subject to locally uniform approximation, neglecting multiple scatterings or resonant phenomena, and, thus, have limited design complexity and physical behavior [17]. In contrast to previous works, our approach considers axisymmetric multilayered freeform metalenses which can be modelled by rigorous fullwave Maxwell equations without any uncontrolled approximations.
The prospect of fabricating single-layer metasurfaces with traditional single-step lithography processes is very promising for large-scale high-throughput integration [6,7]. However, single-layer metasurfaces have been limited in their functionality to narrow angular and spectral bandwidths of operation, with some progress towards chromatic [9] and geometrical aberration correction [18]. Achieving truly multifunctional metasurfaces requires more advanced designs, such as closely-packed multilayer structures [10,11,13,19]. Recently, there has been a surge in interest in fabricating multilayer metasurfaces [20,21], bolstered by advances in inverse-designed nanophotonics. However, these designs can only be fabricated with more advanced fabrication techniques, such as multi-step lithography, or multi-photon lithography. Multi-photon lithography/polymerization is a technique that enables the fabrication of sub-micron (down to ∼ 150 nm) arbitrary three-dimensional structures. Two-photon polymerization enabled the demonstration of three-dimensional chiral/helical structures [22,23] and was more recently applied to the fabrication of supercritical lenses [24,25] and to the demonstration of full three dimensional control of optical fields using a metasurface [26]. Nonetheless, the possibility of 3D printing inverse-designed metasurfaces with two-photon lithography processes remains largely unexplored. In this work, we realize a proof-of-concept experiment with an inverse-designed, freeform, single-layer metalens working at λ = 1550 nm fabricated via two-photon polymerization.
The use of TO as a tool for inverse design in nanophotonics has increased steadily over the last two decades [27,28] with a recent application to metalens design [19]. Our proposed multilayer metalens design framework utilizes density-based topology optimization [29] as the inverse design tool. Rather than allowing fully free-form 3D designs, not amenable to nano-scale fabrication, we propose using TO to control the radial height-variation of the N -layers constituting the lens. In addition, we propose using a filtering technique [30,31] combined with a threshold operation to limit the gradient of the height variations, making it possible to ensure that they comply with fabrication constraints.

II. MODEL AND DESIGN PROBLEM
The design problem is modelled in an axisymmetric domain, Ω, sketched in Fig. 1A(right). The model domain contains a designable region, Ω D (gray), where the metalens is placed on a solid material surface (dark gray). The sketch also indicates the plane, Γ FF (magenta line), used for computing the far-field transformation in Eq. 1, the focal plane(s), Γ FP (blue and red lines), and the focal spot(s), r FP (black dots), of the lens. Finally, the model domain is truncated using a perfectly matching layer in Ω PML [32,33]. The lens-design itself consists of N layers of material ( Fig. 1A(left)), each with a variable height controlled by the design field,ξ L (r). Each designable layer is separated from the next by a layer of air (light gray) and a layer of solid material (dark gray) of fixed thicknesses.
The physics is modelled in Ω using the Maxwell Equations [34] assuming time-harmonic behavior. Doing so, we capture the full wave-behaviour of the electromagnetic field without simplifying assumptions, thus enabling the exploitation of the full wave-behavior to design metalenses exerting precise control of the electromagnetic field. To make the model problem numerically tractable for large design domains, we assume an axisymmetric geometry. This enables the reduction of the full 3D Maxwell Equations to their 2D axisymmetric counterparts [35], therefore significantly reducing the computational effort required to solve the model problem at the cost of a geometric restriction on the design.
The model problem is solved using the scattered-field formulation, E = E b + E s , where the background field, E b , is taken to be a planewave propagating along the z-axis (decomposed into two counter-rotating circularly polarized waves). The model problem is discretized using the Finite Element Method (FEM) [35] and solved using COMSOL Multiphysics v5.5 [36].
A far-field transformation [37] may be used to compute the electric field at any point in space given knowledge of the field in the plane Γ FF above the lens, see Fig. 1. The use of this transformation removes the need for simulating the spatial domain between the lens and focal point, hereby significantly reducing computational cost. The far field transformation may be written as, Here E FF (r) denotes the electric far field at the point r, G E (r, r ) and G H (r, r ) denotes the electric and magnetic field Green's functions, respectively. Finally K(r ) and J(r ) denote the equivalent magnetic and electric surface currents computed from the electric and magnetic near field obtained by solving the model problem.
The figure of merit (FOM) used in the design process is the electric field-intensity at the focal point, r FP . The design problem is formulated as the following continuous constrained optimization problem, Here ξ L (r) denotes a radially-varying design field, which controls the thickness of the L'th layer of the metalens. The electric field at the focal point, E FF (r FP , ξ), is computed using the solution to the physical model problem and Eq. 1.
We propose using a standard PDE-filter [31] to limit the layer-thickness gradient, by applying it to ξ L (r) through the choice of filter radius, r f . After filtering we propose using ξ L (r) to control the layer height through the smoothed threshold operation [38] as,ξ Here z L denotes the spatial position inside each designable layer. The value z L = A L corresponds to the bottom of the designable region in the L'th layer, and the value z L = B L corresponds to the top of the designable region in the L'th layer. The threshold sharpness is controlled by β. In the limit of β → ∞ the fieldξ L takes the value 0 when z L > ξ L and the value 1 when z L < ξ L . A continuation approach may be used to gradually increase β during the inverse design process to enforce a 0/1 final design.
The fieldξ L is used to interpolate the relative permittivity, ε r (r), in space between the background material and the material constituting the metalens using a linear scheme, ε r (r) = ε r,bg +ξ L (ε r,lens − ε r,bg ).
Here ε r,bg (resp. ε r,lens ) denotes the relative permittivity of the background (resp. lens). The design problem, Eqs. (2-3), is solved using the Method of Moving Asymptotes [39], for which the sensitivites of the FOM are computed using adjoint sensitivity analysis [40]. Details regarding the modelling and optimization process as well as the parameter choices for each example are found in Appendix A and Appendix B, respectively.
The final designs are all evaluated numerically using a high resolution model by exciting the lens using a linearly polarized planewave decomposed into two counter-rotating circularly polarized waves introduced in the model using a first order scattering boundary condition. An example of a reconfigurable metalens operating at λ = 10 µm for two different refractive indices is shown in Fig. 1B.

III. MULTI-SCALE MULTI-WAVELENGTH MULTILAYER METALENS
As the first example of our framework we tailor a 10-layer silicon (n = 3.46) in air metalens to focus λ 1 = 1 µm light ( Fig. 2A) and λ 2 = 10 µm light (Fig. 2B) simultaneously at the same focal spot (NA= 0.85). The lens is 100 µm in diameter and has a thickness of 10 µm. The inversely-designed lens is presented in Fig. 2E with the insert showing an example of the layer-height variations. E) 3D rendering of the metalens design.
From Fig. 2A-2B it is clear that the lens exhibits the desired numerical aperture (green line). The focusing capability of the lens is found to reach the diffraction-limit for both wavelengths, when measured in terms of the Full Width at Half Maximum (FWHM) of the main lobe in the focal plane ( Fig. 2C-D). The Strehl ratio (SR) at the two targeted wavelengths, λ 1 = 1 µm and λ 2 = 10 µm, is computed to SR =≈ 0.60 and SR =≈ 0.84, respectively. The SR is computed based on the power flow through the focal plane (blue lines) and the corresponding Airy discs (dashed red lines), shown in Fig. 2C-2D. The absolute power transmission from the substrate of silicon through the lens is computed to T A,λ1 ≈ 82% and T A,λ1 ≈ 95%, relative to the incident power in the silicon substrate within the lens diameter. Appendix C includes an additional design example targeting NA= 0.65 rather than NA=0.85 while keeping all other parameters fixed, demonstrating the methods versatility. For that second example we also achieve diffraction-limited focusing and attain Strehl ratios of SR≈ 0.66 and SR≈ 0.99 for λ 1 = 1 µm and λ 2 = 10 µm, respectively.

IV. TUNABLE MULTI-WAVELENGTH MULTILAYER METALENS
As a second example of our framework, we design of a 10-layer tunable three-wavelength metalens (see Fig. 1C) capable of shifting the numerical aperture of the lens from NA= 0.7 (see Fig. 3[Left Column]) to NA= 0.8 (see Fig. 3[Right Column]) by changing the refractive index of the active material (GST41T1 [3]) from n = 3.2 to n = 4.6.

V. EXPERIMENTAL VALIDATION OF SINGLE-LAYER VARIABLE-HEIGHT METALENS
Finally, as a proof of concept, we demonstrate experimentally that the proposed method can be used to design variable-height metasurfaces for given fabrication specifications (details about the fabrication and experiment are given in Appendix D and Appendix E). Figure 4G shows a 3D rendering of the designed single-layer varying-height metalens. The metalens is fabricated via 3D two-photon polymerization in IP-Dip, a low-refractive-index polymer [41] that can be printed in voxel sizes with in-plane feature sizes ∼ 100 nm and fixed voxel aspect ratio of ∼ 1 to 3. This example is not aimed at designing the largest area lens   Fig. 3.
possible nor at achieving the highest possible numerical performance, but at designing a metalens that complies with fabrication constraints. In this respect, the design is restricted to a diameter of 200 µm with a 300 nm radial pixel size and a varying height with a maximum height of 900 nm, restricted to height-variations in 100 nm increments. The height of the individual radial pixel is allowed to vary independently of its neighbors (i.e. no filtering is applied to ξ L ). The lens is designed to focus λ = 1550 nm light at normal incidence with a numerical aperture of 0.4. The numerically computed electric-field intensity at 1550 nm for planewave illumination of the lens at normal incidence is shown in Fig. 4A, clearly showing that the targeted numerical aperture (green line) is achieved. Numerically the lens achieves near diffractionlimited focusing in terms of the FWHM of main lobe of the power flow in the z-direction through the focal plane. A FWHM of ≈ 1000 nm is computed numerically, corresponding to ≈ 3.2% above the diffraction limit (Using the theoretical limit λ 2NA ≈ 969 nm).
The absolute power transmission from the substrate of IP-Dip through the lens is computed at T A ≈ 93%, relative to the incident power in the IP-Dip substrate within the lens diameter. A Strehl ratio of SR ≈ 0.29 is computed by numerical integration of the power flow over the focal plane. This SR value reveals that a significant fraction of the power is not flowing through the focal point. From a design point of view, the Strehl ratio is easy to improve using our framework by increasing the design freedom, either by changing the metasurface material; by decreasing the radial pixel size; by increasing the number of height increments; by increasing the total height of the lens and/or by introducing multiple-layers in the lens. All of these were demonstrated in the two previous examples.
Experimentally the Strehl ratio is estimated to be ≈ 0.64 by integrating the power flow over an 8 µm×8.5 µm region centered at the focal spot. Computing the SR numerically using the same integration area we obtain SR ≈ 1.0 showing that a majority of the power transmitted through the lens is not focused at the focal spot but flows through the focal plane outside this area. The discrepancy between the experimentally measured and numerically computed SR suggests that the experiment overestimates the SR, due to the camera's limited field of view. This is supported by the relatively low measured absolute focusing efficiency of ≈ 5%. The measured focal spot (Fig. 4(D-F)) exhibits FWHMs of 2.28 ± 0.16 µm (resp. 2.22 ± 0.17 µm) along the horizontal (resp. vertical) direction, corresponding to 18 ± 8% (resp. 15 ± 9%) above the diffraction limit. These experimental results validate the feasability of freeform axisymmetric metasurfaces experimentally. While this proof-of-concept experiment was limited to a single-layer metasurface, the radially-varying height of the structure can, to the authors knowledge, only be implemented with fabrication techniques such as 2.5D lithography or multi-photon polymerization. This is a first step towards realizing the full potential of the freeform axisymmetric inverse design technique presented in this work.

VI. CONCLUSION
In this paper, we demonstrated that fullwave Maxwell Equation based inverse design of axisymmetric structures can tackle challenging new design problems involving radically different wavelengths or active materials. We believe that the proposed design framework opens the way to many new applications whose functionality goes far beyond traditional lenses, such as end-to-end design [42], hyperspectral imaging [43], depth sensing [44] and nonlinear imaging [45].
Computationally, there are many ways to scale our algorithm to much larger designs. The simplest improvement would be to incorporate near-to-farfield transformations [37] to omit the homogeneous region above the lens from the computation, which would allow us to increase the radial size by a factor of ∼ 10. Approximate domain-decomposition could be used to partition a larger lens into overlapping subdomains solved in parallel (but optimized together) [13]. To increase design freedom, the axisymmetry could be relaxed to various forms of N -fold or other rotational symmetries. One could also explore fully free-form topology optimization for 3D-printed structures with manufacturability constraints [46][47][48][49].
Experimentally, we have shown a proof-of-concept fabricated structure using a two-photon 3D-lithography process. The inverse-designed metasurface achieved focusing at the telecommunication wavelength of 1550 nm, close to the diffraction limit, with a numerical aperture of 0.4. In the future, we will develop multilayer fabrication of these structures, in order to realize the full potential of the design technique developed in this work. A key challenge is to realize mechanically stable multilayered structures from which unpolymerized resist can be extracted. Application-specific two-photon polymerization setups [24] can achieve more height levels and some control over the voxel aspect ratio. For devices operating at shorter wavelengths, thus requiring proportionately smaller feature sizes, the design process would shift to multilayer structures with piecewise-constant cross-section [10,50]. Conversely, at longer wavelengths such as for microwave wavefront shaping, multilayer structures could be straightforwardly fabricated, for instance, by stacking multiple stacks of 3D-printed resins or drilled materials [51].

APPENDIX B. STUDY PARAMETERS
The parameters used in setting up the models and associated optimization problems for the three examples follow here.

B. 1. Multi-scale multi-wavelength multilayer metalens
For the problem treated in Sec. III the following parameter values are used: The axisymmetric model domain Ω has a width of 57 µm in the r-direction and a height of 82 µm in the z-direction. Ω is surrounded on three of four sides by a perfectly matched layer with a depth of 1500 nm (Fig. 1). The metalens design domain Ω D is taken to have a radius of 50 µm and a height of 10 µm and is separated into ten layers of equal height. Each layer has a total height of 1 µm with the designable region having a height of 600 nm and the fixed air and silicon regions each having heights of 200 nm. It is placed on a slab of material of 2 µm thickness placed at the bottom edge of the model domain.
The radial design pixel size is restricted to a minimum of 200 nm and the height-variation is restricted to 25 nm increments.
The two wavelengths of the incident field are taken to be λ 1 = 1 µm and λ 2 = 10 µm. The lens is taken to be made of silicon in an air background. The refractive index of air are taken to be n air = 1.0. The refractive index of silicon is taken to be n si = 3.46 at both operating wavelengths. The speed of light is taken to be c = 3·10 8 m/s. The numerical aperture is taken to be NA= 0.65.
The initial guess for the design field is ξ L,initial (r) = 0.5 ∀ r ∈ Ω D for all 10 layers. A filter radius of r f = 400 nm is used to limit the gradient of the heigh variation in each layer to avoid rapid pixel-by-pixel oscillations in the design. The value of the thresholding sharpness parameters is β = 40.

B. 2. Tunable multi-wavelength multilayer metalens
For the problem treated in Sec. IV the following parameter values are used: The axisymmetric model domain Ω has a width of 342.5 µm in the r-direction and a height of 380 µm in the z-direction. Ω is surrounded on three of four sides by a perfectly matched layer with a depth of 15 µm (Fig. 1). The metalens design domain Ω D is taken to have a radius of 312.5 µm and a height of 25 µm and is separated into ten layers of equal height with a 2000 nm designable region and 250 nm fixed air region and 250 nm fixed solid region. It is placed on a slab of material of 5 µm thickness placed at the bottom edge of the model domain.
The radial design pixel size is restricted to a minimum of 600 nm and the height-variation is restricted to 100 nm increments.
The three wavelengths of the incident field are taken to be λ 1 = 9.7 µm, λ 1 = 10 µm and λ 2 = 10.3 µm. The lens is taken to be made of GST41T1 in an air background. The refractive index of air are taken to be n air = 1.0. The refractive index of the active material is taken to be n GST,1 = 3.2 in the first configuration and n GST,2 = 4.6 in the second at all operating wavelengths. The speed of light is taken to be c = 3 · 10 8 m/s. The numerical aperture of the lens is taken to be NA= 0.7 in the first configuration and NA= 0.8 in the second.
The initial guess for the design field is ξ L,initial (r) = 0.5 ∀ r ∈ Ω D for all 10 layers. A filter radius of r f = 3 µm is used to limit the gradient of the height-variation in each layer (see the insert in Fig. 1C). The value of the thresholding sharpness parameters is β = 40.

B. 3. Single-layer variable-height metalens
For the problem treated in Sec. V the following parameter values are used: The axisymmetric model domain Ω has a width of 106 µm in the r-direction and a height of 301.8 µm in the z-direction. Ω is surrounded on three of four sides by a perfectly matched layer with a depth of 3 µm (Fig. 1). The metalens design domain Ω D is taken to have a radius of 100 µm and a height of 900 nm and comprises a single layer constituting the designable region. The design domain is placed on a slab of material of 500 nm thickness placed at the bottom edge of the model domain.
The design is discretized into 300 nm radial increments and 100 nm height increments.
The wavelength of the incident field is taken to be λ = 1550 nm. The lens is taken to be made of IP-Dip in an air background. The refractive index of air are taken to be n air = 1.0. The refractive index of IP-Dip is taken to be n si = 1.507 at both operating wavelengths. The speed of light is taken to be c = 3 · 10 8 m/s. The numerical aperture is taken to be NA= 0.4.
The initial guess for the design field is ξ L,initial (r) = 0.5 ∀ r ∈ Ω D . No smoothing filter is applied. The value of the thresholding sharpness parameters is β = 40.

APPENDIX C. SECOND EXAMPLE OF A MULTI-SCALE MULTI-WAVELENGTH MULTILAYER METALENS DESIGN
We tailor a 10-layer silicon (n = 3.46) in air metalens to focus λ 1 = 1 µm light ( Fig. 2A) and λ 2 = 10 µm light (Fig. 2B) simultaneously at the same focal spot (NA= 0.65). The lens has identical dimensions and design resolution as the lens in Sec. III. The final lens design is presented in Fig. 5E with the insert showing an example of the layer-height variations. Figures 2A-2B show that the lens exhibits the desired numerical aperture at both wavelengths (green line). Further, the focusing capability of the lens is diffraction-limited for both wavelengths. The Strehl ratio (SR) at the two targeted wavelengths, λ 1 = 1 µm and λ 2 = 10 µm, is computed to SR =≈ 0.66 and SR =≈ 0.99, respectively, from the data in Fig. 2C-2D.

APPENDIX D. FABRICATION
The metalens was fabricated using a commercial two-photon polymerization system (Nanoscribe Photonic Professional GT) on a 700-micron-thick fused silica substrate, where the structures are written in circles with height increments of 100nm. For this purpose, piezo actuators move the sample in the out-of-plane direction after fabricating each layer. Geometrical parameters and dose (scanning speed and laser power) are optimized with a dose test on this specific machine. In the in-plane direction the laser beam is guided by galvanometric mirrors parallel to the substrate. After printing, the structures are put in a developer bath (PGMEA 5 min) and dried in IPA with a critical point dryer Auto Samdri 815 Series A.

APPENDIX E. EXPERIMENT
For the proof-of-concept experimental results presented in Fig. 4, we used the imaging setup shown in Fig. 6(a). A Ando AQ4321D Tunable Laser Source produces a fiber-coupled output at 1550 nm. The fiber output is collimated with a set of lenses. In the measuring configuration Fig. 6, the collimated beam is focused by the metasurface, and the focal spot is imaged by an objective -tube lens -IR imaging camera system. For this measurement, we used a 100X Mitutoyo Plan Apo NIR HR Infinity Corrected Objective, a ThorLabs f = 200mm tube lens, and a EC MicronViewer 7290A. The imaging setup was first calibrated using the configuration shown in Fig. 6(b), where the equivalent pixel size on the detector is evaluated by imaging a USAF1951 target. To evaluate the efficiency of the metasurface, we measured the equivalent power going through a 200 µm diameter pinhole with the configuration shown in Fig. 6(c).
To estimate the metasurface efficiency, we use the intensity-voltage relation of the NIR camera provided by the vendor. It has the form I = KV 1/g s , where I is the incident optical power on a pixel, V s the generated voltage at that pixel, and g the characteristic nonlinear slope of the intensity-voltage relation, which is given to be g ∼ 0.7. We first calibrate the proportionality constant K by measuring the signal produced by the camera of a known beam power. This allows us to translate the measured voltage on a pixel to an incident power (in W). We also measure the incident intensity on the metasurface area with the experimental configuration shown in Fig. 6(c). The efficiency is then calculated as where L is the estimated optical loss through the objective and tube lens, which is 0.55 (objective) × 0.88 (tube lens). We typically remove the background from the measured focal spot in order to estimate the metasurface efficiency. The objective-tube lens-camera assembly can be translated along the optic axis. The setup is first calibrated by imaging a USAF1951 target (b). To calibrate our power estimates, we measured the power going through a pinhole with the setup shown in (c).