## Abstract

We demonstrate new axisymmetric inverse-design techniques that can solve problems radically different from traditional lenses, including reconfigurable lenses (that shift a multi-frequency focal spot in response to refractive-index changes) and widely separated multi-wavelength lenses (*λ* = 1 *µ*m and 10 *µ*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 ∼10^{5} degrees of freedom can tackle challenging design problems even when restricted to axisymmetric structures.

© 2020 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. 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 ($\lambda$) in diameters, combined with multilayer variable-height topology optimization (TO) with $\gtrsim 10^4$ degrees of freedom ($\sim 10$ per $\lambda$ 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. 3), we design a *multi-scale* metalens that *simultaneously* focuses $\lambda = 1\,\mu$m and $\lambda =10\,\mu$m at the same diffraction-limited focal point (numerical aperature NA $=0.85$, Strehl ratios of $0.60$ and $0.84$ and efficiencies of $82$% and $95$%, respectively). Second (in Sec. 4), we design an *active* metalens that shifts its achromatic multi-wavelength (three $\lambda$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. 5) an experimental realization of single-layer axisymmetric TO-designed metalens for $\lambda = 1550$ nm, fabricated by two-photon polymerization 3D-printing (Nanoscribe Professional GT), demonstrating that such variable-height surfaces are manufacturable. As discussed in Sec. 6, 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–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 ($\lesssim \lambda$) [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 $\sim$ 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 $\lambda = 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 $\mathcal {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.

## 2. Model and design problem

The design problem is modelled in an axisymmetric domain, $\Omega$, sketched in Fig. 1(A) (right). The model domain contains a designable region, $\Omega _D$ (gray), where the metalens is placed on a solid material surface (dark gray). The sketch also indicates the plane, $\Gamma _{\mathrm {FF}}$ (magenta line), used for computing the far-field transformation in Eq. (1), the focal plane(s), $\Gamma _{\mathrm {FP}}$ (blue and red lines), and the focal spot(s), $\textbf {r}_{\mathrm {FP}}$ (black dots), of the lens. Finally, the model domain is truncated using a perfectly matching layer in $\Omega _{\mathrm {PML}}$ [32,33]. The lens-design itself consists of $\mathcal {N}$ layers of material (Fig. 1(A) (left)), each with a variable height controlled by the design field, $\bar {\xi }_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 $\Omega$ using the Maxwell equations [34] assuming time-harmonic behavior. Doing so, we capture the full-wave behavior 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, $\textbf {E} = \textbf {E}_{\mathrm {b}} + \textbf {E}_{\mathrm {s}}$, where the background field, $\textbf {E}_{\mathrm {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 $\Gamma _{\mathrm {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,

The figure of merit (FOM) used in the design process is the electric field-intensity at the focal point, $\textbf {r}_{\mathrm {FP}}$. The design problem is formulated as the following continuous constrained optimization problem,

We propose using a standard PDE-filter [31] to limit the layer-thickness gradient, by applying it to $\xi _L(r)$ through the choice of filter radius, $r_f$. After filtering we propose using $\xi _L(r)$ to control the layer height through the smoothed threshold operation [38] as,

The field $\bar {\xi }_L$ is used to interpolate the relative permittivity, $\varepsilon _r(\textbf {r})$, in space between the background material and the material constituting the metalens using a linear scheme,

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 $\lambda = 10~\mu$m for two different refractive indices is shown in Fig. 1(B).

## 3. 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 $\lambda _1 = 1 \ \mu$m light (Fig. 2(A)) and $\lambda _2 = 10 \ \mu$m light (Fig. 2(B)) simultaneously at the same focal spot (NA$=0.85$). The lens is 100 $\mu$m in diameter and has a thickness of 10 $\mu$m. The inverse-designed lens is presented in Fig. 2(E) with the insert showing an example of the layer-height variations.

From Figs. 2(A)–2(B) 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 (Figs. 2(C) and 2(D)). The Strehl ratio (SR) at the two targeted wavelengths, $\lambda _1 = 1 \ \mu$m and $\lambda _2 = 10 \ \mu$m, is computed to be SR $\approx 0.60$ and SR $\approx 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 Figs. 2(C)–2(D). The absolute power transmission from the substrate of silicon through the lens is computed to be $\mathrm {T}_{\mathrm {A},\lambda _1} \approx 82\%$ and $\mathrm {T}_{\mathrm {A},\lambda _2} \approx 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$\approx 0.66$ and SR$\approx 0.99$ for $\lambda _1 = 1 \ \mu$m and $\lambda _2 = 10 \ \mu$m, respectively.

To illustrate the benefit of the proposed multi-layer metalens over a single-layer lens we consider a simple single-layer reference design (Fig. 2(F)). The single-layer design is optimized using our proposed approach, with all parameters used in the example held constant, except for the number of layers. For this single-layer reference design we obtain a Strehl Ratio of $\approx 0.37$ and $\approx 0.09$ and an absolute transmission efficiency of $\mathrm {T}_{\mathrm {A},\lambda _1} \approx 57\%$ and $\mathrm {T}_{\mathrm {A},\lambda _2} \approx 74\%$ at $\lambda _1$ and $\lambda _2$, respectively. Comparing the Strehl ratios and absolute transmission efficiencies to those obtained for the ten layer metalens design, the benefit of multi-layer metalens designs over single-layer designs for this example is clear.

## 4. Tunable multi-wavelength multilayer metalens

As a second example of our framework, we design of a 10-layer tunable three-wavelength metalens (see Fig. 1(C)) 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 + 0.01 \mathrm {i}$.

The lens is 625 $\mu$m in diameter and has a thickness of 25 $\mu$m. The lens is designed to operate in the mid-infrared region at wavelengths, $\lambda _1 = 9.7 \ \mu$m (Fig. 3 [Row 1]), $\lambda _1 = 10 \ \mu$m (Fig. 3 [Row 2]) and $\lambda _1 = 10.3 \ \mu$m (Fig. 3 [Row 3]). From Fig. 3 it is observed that the lens exhibits the desired numerical aperture at all three wavelengths for both values of the refractive index. The Strehl ratio, absolute power transmission and FWHM of the main lobe at the focal point for the three targeted wavelengths and two refractive indices are presented in Table 1. In brief, a Strehl ratio of approximately $0.5$ is achieved across all six cases with the spatial focusing being at most $11\%$ from the diffraction limit. Finally, a $T_A$ of $\approx 0.3$ for $n = 3.2$ and of $\approx 0.2$ for $n = 4.6 + 0.01 \mathrm {i}$ is achieved.

## 5. Experimental validation of a 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 4(G) 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 $\sim 100$ nm and fixed voxel aspect ratio of $\sim 1$ to 3. This example is not aimed at designing the largest area lens 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 $\mu$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 $\xi _{L}$).

The lens is designed to focus $\lambda = 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. 4(A), clearly showing that the targeted numerical aperture (green line) is achieved. Numerically the lens achieves near diffraction-limited focusing in terms of the FWHM of main lobe of the power flow in the z-direction through the focal plane. A FWHM of $\approx 1000$ nm is computed numerically, corresponding to $\approx 3.2\%$ above the diffraction limit (Using the theoretical limit $\frac {\lambda }{2\mathrm {NA}}\approx 969$ nm).

The absolute power transmission from the substrate of IP-Dip through the lens is computed at $\mathrm {T}_{\mathrm {A}} \approx 93\%$, relative to the incident power in the IP-Dip substrate within the lens diameter. A Strehl ratio of $\mathrm {SR} \approx 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 $\approx 0.64$ by integrating the power flow over an $8~\mu \mathrm {m} \times 8.5~\mu \mathrm {m}$ region centered at the focal spot. Computing the SR numerically using the same integration area we obtain $\mathrm {SR} \approx 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 $\approx 5 \%$. The measured focal spot (Figs. 4(D)–4(F)) exhibits FWHMs of $2.28 \pm 0.16 ~\mu m$ (resp. $2.22\pm 0.17~\mu m$) along the horizontal (resp. vertical) direction, corresponding to $18 \pm 8 \%$ (resp. $15 \pm 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.

Achieving true multi-layer closely-packed metasurfaces presents additional challenges, such as the accurate positioning and alignment of each layer. Yet another challenge – which is specific to two-photon polymerization – is to design structures that allow unpolymerized material to be extracted, a constraint that could be included in further refinements of our theory.

## 6. 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]. While we expect a small-angle paraxial regime to be valid for our lens designs, which may be used for imaging over a narrow field of view, we will consider, in a future work, thorough corrections of off-axis as well as chromatic aberrations in a single-piece axisymmetric metalens design.

An example of the significant performance benefits that can be attained by designing multi-layer metalenses, compared to single-layer metalenses, was given in Sec. 3. The relationship between the targeted number of layers and the performance of the lens has not been investigated in detail and such a study for different metalens applications is likely to provide valuable information and is thus interesting to pursue.

Computationally, there are several ways to scale our algorithm to much larger designs. The simplest would be to utilize near-to-farfield transformations [37] to omit simulation of the homogeneous region above the lens from the computation, which would allow us to increase the radial size by a factor of $\sim 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–49].

When employing the proposed approach for materials with a large non-zero extinction coefficient, $\kappa$, i.e. a complex refractive index $\tilde {n} = n + \mathrm {i} \kappa$, it is possible that one needs to consider a different material interpolation scheme, to achieve high quality results from the inverse design process [50].

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,51]. 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 [52].

## Appendix A. Optimization and numerical modelling

The physics is modelled in COMSOL Multiphysics [36] and the optimization problem is solved using the Globally Convergent Method of Moving Asymptotes (GCMMA) [39].

In the design process $\Omega _{\mathrm {D}}$ and the solid material regions in $\Omega$ are discretized using a structured quadrilateral mesh, while the surrounding air regions are discretized using an unstructured triangular mesh, both of which uses $\geq 10$ elements per $\lambda / n$. The finite element method with a linear Lagrangian basis is used to discretize the physics [35].

The following stopping criterion is used to terminate the iterative solution of the optimization problem:

Here $i$ denotes the current optimization iteration, $i_{\mathrm {min}} = 70$ denotes the minimum number of design iterations taken. $\Phi _{i}$ denotes the objective function value at the $i$’th iteration and $n \in \mathbb {N}^{+}$.

## 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. 3 the following parameter values are used:

The axisymmetric model domain $\boldsymbol {\Omega }$ has a width of 57 $\mu$m in the r-direction and a height of 82 $\mu$m in the z-direction. $\boldsymbol {\Omega }$ is surrounded on three of four sides by a perfectly matched layer with a depth of 1500 nm (Fig. 1). The metalens design domain $\boldsymbol {\Omega }_{\mathrm {D}}$ is taken to have a radius of 50 $\mu$m and a height of 10 $\mu$m and is separated into ten layers of equal height. Each layer has a total height of 1 $\mu$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 $\mu$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 $\lambda _1 = 1$ $\mu$m and $\lambda _2 = 10$ $\mu$m. The lens is taken to be made of silicon in an air background. The refractive index of air are taken to be $n_{\mathrm {air}} = 1.0$. The refractive index of silicon is taken to be $n_{\mathrm {si}} = 3.46$ at both operating wavelengths. The speed of light is taken to be $c = 3 \cdot 10^8$ m/s. The numerical aperture is taken to be NA$= 0.65$.

The initial guess for the design field is $\xi _{L,\mathrm {initial}}(\textbf {r}) = 0.5 \ \forall \ \textbf {r} \in \boldsymbol {\Omega }_{\mathrm {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 $\beta = 40$.

## B. 2. Tunable multi-wavelength multilayer metalens

For the problem treated in Sec. 4 the following parameter values are used:

The axisymmetric model domain $\boldsymbol {\Omega }$ has a width of 342.5 $\mu$m in the r-direction and a height of 380 $\mu$m in the z-direction. $\boldsymbol {\Omega }$ is surrounded on three of four sides by a perfectly matched layer with a depth of 15 $\mu$m (Fig. 1). The metalens design domain $\boldsymbol {\Omega }_{\mathrm {D}}$ is taken to have a radius of 312.5 $\mu$m and a height of 25 $\mu$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 $\mu$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 $\lambda _1 = 9.7$ $\mu$m, $\lambda _1 = 10$ $\mu$m and $\lambda _2 = 10.3$ $\mu$m. The lens is taken to be made of GST41T1 in an air background. The refractive index of air are taken to be $n_{\mathrm {air}} = 1.0$. The refractive index of the active material is taken to be $n_{\mathrm {GST,1}} = 3.2$ in the first configuration and $n_{\mathrm {GST,2}} = 4.6 + 0.01 \mathrm {i}$ in the second at all operating wavelengths. The speed of light is taken to be $c = 3 \cdot 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 $\xi _{L,\mathrm {initial}}(\textbf {r}) = 0.5 \ \forall \ \textbf {r} \in \boldsymbol {\Omega }_{\mathrm {D}}$ for all 10 layers. A filter radius of $r_f = 3~\mu$m is used to limit the gradient of the height-variation in each layer (see the insert in Fig. 1(C)). The value of the thresholding sharpness parameters is $\beta = 40$.

## B. 3. Single-layer variable-height metalens

For the problem treated in Sec. 5 the following parameter values are used:

The axisymmetric model domain $\boldsymbol {\Omega }$ has a width of 106 $\mu$m in the r-direction and a height of 301.8 $\mu$m in the z-direction. $\boldsymbol {\Omega }$ is surrounded on three of four sides by a perfectly matched layer with a depth of 3 $\mu$m (Fig. 1). The metalens design domain $\boldsymbol {\Omega }_{\mathrm {D}}$ is taken to have a radius of 100 $\mu$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 $\lambda = 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_{\mathrm {air}} = 1.0$. The refractive index of IP-Dip is taken to be $n_{\mathrm {si}} = 1.507$ at both operating wavelengths. The speed of light is taken to be $c = 3 \cdot 10^8$ m/s. The numerical aperture is taken to be NA$= 0.4$.

The initial guess for the design field is $\xi _{L,\mathrm {initial}}(\textbf {r}) = 0.5 \ \forall \ \textbf {r} \in \boldsymbol {\Omega }_{\mathrm {D}}$. No smoothing filter is applied. The value of the thresholding sharpness parameters is $\beta = 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 $\lambda _1 = 1 \ \mu$m light (Fig. 2(A)) and $\lambda _2 = 10 \ \mu$m light (Fig. 2(B)) simultaneously at the same focal spot (NA$=0.65$). The lens has identical dimensions and design resolution as the lens in Sec. 3. The final lens design is presented in Fig. 5(E) with the insert showing an example of the layer-height variations.Figs. 2(A) and2(B) 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, $\lambda _1 = 1 \ \mu$m and $\lambda _2 = 10 \ \mu$m, is computed to be SR $=\approx 0.66$ and SR $=\approx 0.99$, respectively, from the data in Figs. 2(C) and 2(D).

## 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 $~ 100 nm$. 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 = 200$mm 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~\mu$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 = K V_s^{1/g}$, 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 \sim 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

## Funding

Villum Fonden (8692); Danmarks Grundforskningsfond (DNRF147); Army Research Office (W911NF-18-2-0048); MIT-IBM Watson AI Laboratory (2415); Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (P2EZP2_188091).

## Acknowledgements

Y.S. acknowledges the Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (SNSF) through Project No. P2EZP2_188091.

## Disclosures

The authors declare that there are no conflicts of interest related to this article.

## References

**1. **M. Meem, S. Banerji, C. Pies, T. Oberbiermann, A. Majumder, B. Sensale-Rodriguez, and R. Menon, “Large-area, high-numerical-aperture multi-level diffractive lens via inverse design,” Optica **7**(3), 252–253 (2020). [CrossRef]

**2. **R. Pestourie, C. Pérez-Arancibia, Z. Lin, W. Shin, F. Capasso, and S. G. Johnson, “Inverse design of large-area metasurfaces,” Opt. Express **26**(26), 33732–33747 (2018). [CrossRef]

**3. **Y. Zhang, J. B. Chou, J. Li, H. Li, Q. Du, A. Yadav, S. Zhou, M. Y. Shalaginov, Z. Fang, H. Zhong, C. Roberts, P. Robinson, B. Bohlin, C. Ríos, H. Lin, M. Kang, T. Gu, J. Warner, V. Liberman, K. Richardson, and J. Hu, “Broadband transparent optical phase change materials for high-performance nonvolatile photonics,” Nat. Commun. **10**(1), 4279 (2019). [CrossRef]

**4. **M. Y. Shalaginov, S. An, Y. Zhang, F. Yang, P. Su, V. Liberman, J. B. Chou, C. M. Roberts, M. Kang, C. Rios, Q. Du, C. Fowler, A. Agarwal, K. Richardson, C. Rivero-Baleine, H. Zhang, J. Hu, and T. Gu, “Reconfigurable all-dielectric metalens with diffraction limited performance,” arXiv:1911.12970 (2019).

**5. **N. Yu and F. Capasso, “Flat optics with designer metasurfaces,” Nat. Mater. **13**(2), 139–150 (2014). [CrossRef]

**6. **M. Khorasaninejad, W. T. Chen, R. C. Devlin, J. Oh, A. Y. Zhu, and F. Capasso, “Metalenses at visible wavelengths: Diffraction-limited focusing and subwavelength resolution imaging,” Science **352**(6290), 1190–1194 (2016). [CrossRef]

**7. **M. Khorasaninejad and F. Capasso, “Metalenses: Versatile multifunctional photonic components,” Science **358**(6367), eaam8100 (2017). [CrossRef]

**8. **F. Aieta, M. A. Kats, P. Genevet, and F. Capasso, “Multiwavelength achromatic metasurfaces by dispersive phase compensation,” Science **347**(6228), 1342–1345 (2015). [CrossRef]

**9. **W. T. Chen, A. Y. Zhu, V. Sanjeev, M. Khorasaninejad, Z. Shi, E. Lee, and F. Capasso, “A broadband achromatic metalens for focusing and imaging in the visible,” Nat. Nanotechnol. **13**(3), 220–226 (2018). [CrossRef]

**10. **Z. Lin, V. Liu, R. Pestourie, and S. G. Johnson, “Topology optimization of freeform large-area metasurfaces,” Opt. Express **27**(11), 15765–15775 (2019). [CrossRef]

**11. **H. Chung and O. D. Miller, “High-NA achromatic metalenses by inverse design,” Opt. Express **28**(5), 6945–6965 (2020). [CrossRef]

**12. **F. Presutti and F. Monticone, “Focusing on bandwidth: achromatic metalens limits,” Optica **7**(6), 624–631 (2020). [CrossRef]

**13. **Z. Lin and S. G. Johnson, “Overlapping domains for topology optimization of large-area metasurfaces,” Opt. Express **27**(22), 32445–32453 (2019). [CrossRef]

**14. **T. Phan, D. Sell, E. W. Wang, S. Doshay, K. Edee, J. Yang, and J. A. Fan, “High-efficiency, large-area, topology-optimized metasurfaces,” Light: Sci. Appl. **8**(1), 48 (2019). [CrossRef]

**15. **S. Banerji, M. Meem, A. Majumder, F. G. Vasquez, B. Sensale-Rodriguez, and R. Menon, “Imaging with flat optics: metalenses or diffractive lenses?” Optica **6**(6), 805–810 (2019). [CrossRef]

**16. **M. Meem, S. Banerji, A. Majumder, J. C. Garcia, P. W. Hon, B. Sensale-Rodriguez, and R. Menon, “Imaging from the visible to the longwave infrared wavelengths via an inverse-designed flat lens,” arXiv:2001.03684 (2020).

**17. **J. Engelberg and U. Levy, “The advantages of metalenses over diffractive lenses,” Nat. Commun. **11**(1), 1991 (2020). [CrossRef]

**18. **B. Groever, W. T. Chen, and F. Capasso, “Meta-lens doublet in the visible region,” Nano Lett. **17**(8), 4902–4907 (2017). [CrossRef]

**19. **Z. Lin, B. Groever, F. Capasso, A. W. Rodriguez, and M. Lončar, “Topology-optimized multilayered metaoptics,” Phys. Rev. Appl. **9**(4), 044030 (2018). [CrossRef]

**20. **M. Mansouree, H. Kwon, E. Arbabi, A. McClung, A. Faraon, and A. Arbabi, “Multifunctional 2.5 d metastructures enabled by adjoint optimization,” Optica **7**(1), 77–84 (2020). [CrossRef]

**21. **Y. Zhou, I. I. Kravchenko, H. Wang, J. R. Nolen, G. Gu, and J. Valentine, “Multilayer noninteracting dielectric metasurfaces for multiwavelength metaoptics,” Nano Lett. **18**(12), 7529–7537 (2018). [CrossRef]

**22. **J. K. Gansel, M. Thiel, M. S. Rill, M. Decker, K. Bade, V. Saile, G. von Freymann, S. Linden, and M. Wegener, “Gold helix photonic metamaterial as broadband circular polarizer,” Science **325**(5947), 1513–1515 (2009). [CrossRef]

**23. **S. Wu, S. Xu, T. L. Zinenko, V. V. Yachin, S. L. Prosvirnin, and V. R. Tuz, “3d-printed chiral metasurface as a dichroic dual-band polarization converter,” Opt. Lett. **44**(4), 1056–1059 (2019). [CrossRef]

**24. **W. Fang, J. Lei, P. Zhang, F. Qin, M. Jiang, X. Zhu, D. Hu, Y. Cao, and X. Li, “Multilevel phase supercritical lens fabricated by synergistic optical lithography,” Nanophotonics **9**(6), 1469–1477 (2020). [CrossRef]

**25. **J. Coompson, M. Liang, C. Auginash, A. Gin, M. Yang, Z. Qu, I. B. Djordjevic, and H. Xin, “3d-printed phase controlled focusing metalens at 1550 nm wavelength,” in 2018 IEEE International Symposium on Antennas and Propagation & USNC/URSI National Radio Science Meeting, (IEEE, 20181685–1686

**26. **A. Zhan, R. Gibson, J. Whitehead, E. Smith, J. R. Hendrickson, and A. Majumdar, “Controlling three-dimensional optical fields via inverse mie scattering,” Sci. Adv. **5**(10), eaax4769 (2019). [CrossRef]

**27. **J. S. Jensen and O. Sigmund, “Topology optimization for nano-photonics,” Laser Photonics Rev. **5**(2), 308–321 (2011). [CrossRef]

**28. **S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vuckovic, and A. W. Rodriguez, “Inverse design in nanophotonics,” Nat. Photonics **12**(11), 659–670 (2018). [CrossRef]

**29. **M. P. Bendsøe and O. Sigmund, * Topology Optimization* (Springer, 2003).

**30. **B. Bourdin, “Filters in topology optimization,” Int. J. Numer. Meth. Eng. **50**(9), 2143–2158 (2001). [CrossRef]

**31. **B. S. Lazarov and O. Sigmund, “Filters in topology optimization based on helmholtz-type differential equations,” Int. J. Numer. Meth. Eng. **86**(6), 765–781 (2011). [CrossRef]

**32. **J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” J. Comput. Phys. **114**(2), 185–200 (1994). [CrossRef]

**33. **W. C. Chew, J. M. Jin, and E. Michielssen, “Complex coordinate stretching as a generalized absorping boundary condition,” Microw. Opt. Technol. Lett. **15**(6), 363–369 (1997). [CrossRef]

**34. **D. J. Griffiths, * Introduction to Electrodymanics - Fourth Edition* (Pearson Education Limited, 2014).

**35. **J.-M. Jin, * The Finite Element Method in Electromagnetics - Third Edition* (Wiley-IEEE, 2014).

**36. **“Comsol multiphysics® v. 5.5. www.comsol.com,”.

**37. **A. Taflove and S. C. Hagness, * Computational Electrodynamics: the Finite-Difference Time-Domain Method* (Artech house, 2005).

**38. **F. Wang, B. S. Lazarov, and O. Sigmund, “On projection methods, convergence and robust formulations in topology optimization,” Struct. Multidisc. Optim. **43**(6), 767–784 (2011). [CrossRef]

**39. **K. Svanberg, “A class of globally convergent optimization methods based on conservative convex separable approximations,” SIAM J. Optim. **12**(2), 555–573 (2002). [CrossRef]

**40. **D. A. Tortorelli and P. Michaleris, “Design sensitivity analysis: Overview and review,” Inverse Probl. Sci. Eng. **1**(1), 71–105 (1994). [CrossRef]

**41. **T. Gissibl, S. Wagner, J. Sykora, M. Schmid, and H. Giessen, “Refractive index measurements of photo-resists for three-dimensional direct laser writing,” Opt. Mater. Express **7**(7), 2293–2298 (2017). [CrossRef]

**42. **Z. Lin, C. Roques-Carmes, R. Pestourie, M. Soljačić, A. Majumdar, and S. G. Johnson, “End-to-end inverse design for inverse scattering via freeform metastructures,” arXiv:2006.09145 (2020).

**43. **K. Monakhova, K. Yanny, N. Aggarwal, and L. Waller, “Spectral DiffuserCam: Lensless snapshot hyperspectral imaging with a spectral filter array,” arXiv:2006.08565 (2020).

**44. **Q. Guo, Z. Shi, Y.-W. Huang, E. Alexander, C.-W. Qiu, F. Capasso, and T. Zickler, “Compact single-shot metalens depth sensors inspired by eyes of jumping spiders,” P. Natl. Acad. Sci. U.S.A. **116**(46), 22959–22965 (2019). [CrossRef]

**45. **C. Schlickriede, S. S. Kruk, L. Wang, B. Sain, Y. Kivshar, and T. Zentgraf, “Nonlinear imaging with all-dielectric metasurfaces,” Nano Lett. **20**(6), 4370–4376 (2020). [CrossRef]

**46. **R. E. Christiansen, B. S. Lazarov, J. S. Jensen, and O. Sigmund, “Creating geometrically robust designs for highly sensitive problems using topology optimization - acoustic cavity design,” Struct. Multidisc. Optim. **52**(4), 737–754 (2015). [CrossRef]

**47. **M. Zhou, B. S. Lazarov, F. Wang, and O. Sigmund, “Minimum length scale in topology optimization by geometric constraints,” Comput. Method. Appl. Mech. Eng. **293**, 266–282 (2015). [CrossRef]

**48. **F. Wang, R. E. Christiansen, Y. Yu, J. Mørk, and O. Sigmund, “Maximizing the quality factor to mode volume ratio for ultra-small photonic crystal cavities,” Appl. Phys. Lett. **113**(24), 241101 (2018). [CrossRef]

**49. **Q. Li, W. Chen, S. Liu, and L. Tong, “Structural topology optimization considering connectivity constraint,” Struct. Multidisc. Optim. **54**(4), 971–984 (2016). [CrossRef]

**50. **R. E. Christiansen, J. Vester-Petersen, S. P. Madsen, and O. Sigmund, “A non-linear material interpolation for design of metallic nano-particles using topology optimization,” Comput. Method. Appl. Mech. Eng. **343**, 23–39 (2019). [CrossRef]

**51. **A. Y. Piggott, J. Petykiewicz, L. Su, and J. Vučković, “Fabrication-constrained nanophotonic inverse design,” Sci. Rep. **7**(1), 1786 (2017). [CrossRef]

**52. **P. Camayd-Muñoz, C. Ballew, G. Roberts, and A. Faraon, “Multifunctional volumetric meta-optics for color and polarization image sensors,” Optica **7**(4), 280–283 (2020). [CrossRef]