## Abstract

We use inverse design to discover metalens structures that exhibit broadband, achromatic focusing across low, moderate, and high numerical apertures. We show that standard unit-cell approaches cannot achieve high-efficiency high-NA focusing, even at a single frequency, due to the incompleteness of the unit-cell basis, and we provide computational upper bounds on their maximum efficiencies. At low NA, our devices exhibit the highest theoretical efficiencies to date. At high NA—of 0.9 with translation-invariant films and of 0.99 with “freeform” structures—our designs are the first to exhibit achromatic high-NA focusing.

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

## 1. Introduction

Metasurfaces [1,2] are patterned optical thin films that offer the possibility of manipulating light, for applications from holography to lenses, with precision equal to or greater than their bulky conventional-optics counterparts. For metasurface lenses, i.e., metalenses, a now-standard approach of “stitching” together wavelength-scale resonators into a larger device has demonstrated the possibility of focusing [2–9], but has suffered from narrow-bandwidth operation, low-numerical-aperture restrictions, and/or low focusing efficiencies. In this paper, we prove that unit-cell designs *cannot* have high efficiency for high numerical apertures; conversely, we theoretically demonstrate that “inverse design,” a large-scale computational design technique optimizing all geometrical degrees of freedom [10–20] can discover high-numerical-aperture (“fast”) lenses that operate over visible bandwidths with best-in-class focusing efficiencies. Inverse design enables rapid computation of gradients with respect to arbitrarily many geometrical degrees of freedom; our implementation using a minimax formulation of the design criteria identifies fabrication-ready designs that can simultaneously achieve the large bending angles of high-NA lenses with the broad-bandwidth control that is necessary. Our results demonstrate the capabilities of adjoint-based approaches for superior design, and the emerging possibilities for combining multiple high-efficiency, hard-to-achieve functionalities in a single metasurface.

Metasurface designs and implementations have demonstrated the core functionality required for many applications: holography [22–24], retroreflection [25], antennas [26], flat lenses [5–9,27] and tunable optical components [28]. Yet that functionality is typically highly restricted in bandwidth, angular acceptance, and/or numerical aperture. Underlying these restrictions is the intuitive “unit cell” approach that pioneered initial metasurfaces designs [1,2,29], whereby large-area films (tens to thousands of wavelengths in diameter) are constructed from libraries of wavelength-scale unit cells whose outgoing-wave phases, under “locally periodic” boundary conditions, are optimized for a single or few frequencies. This approach has been a crucial first step for designing large-area structures with complex patterning. Yet in the case of metalenses, which require precise control of potentially rapidly varying wavefronts, no metalens designs to date have achieved broad bandwidth and high numerical aperture. And as we show here, the assumption of local periodicity is incompatible with the requirements of high-NA focusing, wherein the amplitude and phase of the outgoing field must undergo rapid variations. Using a basis-projection approach we show that an assumption of local periodicity necessarily entails large focusing-efficiency losses, even at a single frequency.

To circumvent the limitations of the resonator-based unit-cell approach, we instead use inverse design to discover full-wave structures that simultaneously achieve broad bandwidth, high numerical aperture, and high efficiency. Inverse design centers around computing gradients with respect to large numbers of structural degrees of freedom by “adjoint”-based methods. Adjoint-based sensitivities trace their roots to control theory [30–33], and have since been used for rapid, efficient optimization in circuit theory [34], aerodynamics [35], mechanics and elasticity [36], quantum dynamics [37,38], and deep learning [39–41], where it is known as “backpropagation.” More recently it has emerged as a promising design tool for nanophotonics [18,19,42] for applications including waveguide demultiplexers [17,43,44], beam deflectors [45,46], photonic bandgaps [47], solar cells [48], and many others. Preliminary studies have applied inverse design to metasurfaces [13,20,46,49], including large-area metasurfaces [15,50], albeit thus far limited to isolated frequencies or metrics other than lens focusing (beam deflection, polarizers, etc.).

In this work, we use inverse design to demonstrate broadband achromatic metalenses operating with relatively high efficiencies across the visible for both small and large numerical apertures as shown in Fig. 1. To fully explore the design space and the tradeoffs associated with numerical aperture, bandwidth, and efficiency, we primarily design two-dimensional dielectric profiles, while demonstrating that the methods scale to fully three-dimensional films. We make no periodic / unit-cell approximations, and indeed the designs that we discover often show rapidly varying spatial profiles; to ensure feasible computation times, we design devices in the 10–60$\lambda$ size range. In these exponentially large design spaces, one can never certify global optimality (unless there are known global bounds by some other means [51–57]), and we make no claims that our high-performance structures are globally optimal. It is possible that deep learning or related techniques may overcome the non-uniqueness from an inverse-problem perspective [58]. A key advantage of the adjoint-based approach is that computational complexity does not increase with the number of parameters, nor do we empirically find a number-of-parameters cost in the quality of the local optima. (This may be related to the “blessing of dimensionality” in the structure of high-dimensional data [59].) We do find that a combination of minimax optimizations and random parameter initializations (both discussed further below) enables avoidance of low-quality optima to a significant extent. With recent developments in fast electromagnetic solvers [60,61], one can anticipate applying this approach to devices orders of magnitude larger, with small and controllable errors, in the near future. The full-device optimizations enable us to overcome the efficiency losses and single-functionality limitations associated with assuming small, periodic unit cells, which are especially prominent at high NA. Our metalens designs demonstrate the capability for inverse design to lead the discovery of non-intuitive yet superior metalens devices.

## 2. Unit-cell-metalens efficiency limits

The prototypical approach to metalens design is the unit-cell method [3,9,27], wherein the fields at the exit plane of the metalens are designed by breaking them into wavelength-scale “unit cells” that vary slowly and can be treated locally as periodic elements. This reduces the large-scale design problem into a large number of much smaller design problems, while simplifying the physics of the problem. In the prototypical case with a unit-cell period smaller than the wavelength, there is only a single outgoing diffracted-wave order for each cell, and the design problem simplifies to one of controlling the phase and amplitude of each cell’s outgoing wave. Such a formulation extends to broadband metalens operation, where all of the complexity of large-area, many-frequency operation can be simplified to control of phase and its first two derivatives—group velocity and group velocity dispersion—at a single frequency [8]. Yet there are complex tradeoffs associated with the unit-cell approach. For any fixed fabrication tolerance, a larger period allows for more complex unit cells, but as the period increases beyond the wavelength, higher-order modes are introduced, which if not controlled can lead to significant efficiency losses. Conversely, smaller periods reduce the number of geometrical degrees of freedom, particularly if the local-periodicity assumption is to hold.

We use a modal-decomposition analysis to derive bounds on focusing efficiencies for metalenses designed by a unit-cell approach, at any given frequency. We show that even without incorporating deviations from the uncontrolled assumption of periodic boundary conditions, a high numerical aperture *necessarily* incurs significant efficiency losses in the unit-cell approach. The key issue is that the ideal exit-plane fields are non-periodic over each unit cell, and cannot be represented by *periodic* diffraction orders without mismatch. This mismatch, integrated over the entire metasurface diameter, leads to imperfect focusing and efficiency losses. Mathematically, this mismatch is a physical manifestation of the *incompleteness* of the diffraction-order basis in a non-periodic system, an effect that is exacerbated in unit cells with only a small number of diffraction orders.

Consider a lens focusing light to a focal point. At the exit plane of the lens, the optimal fields $\boldsymbol {E}^{\textrm{ideal}}$ are the time-reversed conjugates of those radiated from a dipole at the focal point, with the dipole polarization selected for maximum efficiency. (It can be shown that this is the optimal excitation for field concentration from unit-normalized excitation fields; if sidelobe or spot-size constraints are required, there are known techniques for identifying the optimal exit-plane fields [62–64].) In a unit-cell approach, designing a metalens by stitching together library elements simulated with periodic boundary conditions (over period $\Lambda$) is equivalent to designing the field $\boldsymbol {E}^\textrm {UC,ideal}$ as a linear combination of basis functions $\boldsymbol {u}_{m,\Lambda }(\boldsymbol {x})$ that are the diffraction orders of each “periodic” unit cell:

Every unit-cell metasurface design [5,8,9,22,27,29] is implicitly using Eq. (1) to describe the outgoing fields. The true ideal field is the time-reversed field from a dipole at the focal point, which we denote $\boldsymbol {E}^\textrm {ideal}$. The best possible unit-cell design is the one that minimizes the difference between $\boldsymbol {E}^\textrm {UC,ideal}$ and $\boldsymbol {E}^\textrm {ideal}$; for example, the field that minimizes $\|\boldsymbol {E}^\textrm {ideal} - \boldsymbol {E}^\textrm {UC,ideal}\|_2^{2}$. (Often only the phase is optimized and the efficiency decreases further; since we are interested in upper bounds, we assume the ideal scenario.) By the orthogonality of the basis functions in Eq. (2), there is a set of unit-cell coefficients $c_{m,\Lambda }$ that minimize this least-squares quantity, given by (Appendix)

We can semi-analytically determine the closest approximation to Eq. (1) which *is* a valid solution of Maxwell’s equations. At the exit plane of the metalens, we can write the fields not as a superposition of unit-cell basis functions, but instead as a linear combination of plane waves:

Equation (5) can be evaluated for any unit-cell period $\Lambda$ and numerical aperture, with the resulting focusing-plane fields computed by Eq. (4). The intensity of the fields at the focal point, relative to those of the ideal field $\boldsymbol {E}^\textrm {ideal}$, thereby represents an upper bound on the maximum designable efficiency by a unit-cell approach. We plot the upper bound in Fig. 2(a), showing the steep dropoffs in maximum efficiency as numerical aperture increases. Figure 2(b) shows the ideal and unit-cell-ideal fields at the exit plane of the metasurface, with the latter unable to capture the necessarily rapid variations in phase and amplitude. The prediction of low maximal efficiencies at high NA might seem incompatible with the results of Ref. [27], where single-frequency, high-efficiency (approaching 90%), high-NA (0.8) devices are theoretically simulated and experimentally demonstrated. However, the efficiencies reported in Ref. [27] are polarization-conversion efficiencies across the full transmitted wave, *not* efficiencies of power concentration at the focal point, which are certainly significantly smaller, offering no contradiction with our theoretical predictions.

We emphasize that the efficiency upper bound accounts only for the incompleteness of the unit-cell basis (including effects such as insufficient spatial sampling of the phase [66,67]); it assumes no further losses due to multiple-scattering effects between non-identical neighboring cells (that violate the local-periodicity assumption), and thus almost certainly *over*estimates the maximum efficiency possible in a unit-cell approach. Even so, these modal-decomposition bounds predict significant efficiency losses for unit-cell-based high-numerical-aperture metalenses. An alternative is to design the entire device at once, a task of significant complexity where inverse design may be ideal.

## 3. Metalens inverse design framework

Inverse design requires specification of a figure of merit as well as the geometrical degrees of freedom. To determine the focusing properties of a given metalens geometry, we compute the fields at an exit plane of the metasurface and project them to the far field. The plane-wave decomposition at the exit plane is given by

*maximin*over the geometrical degrees of freedom and the frequency range of interest: We find that use of the minimax approach leads to nearly constant response over broad bandwidths. In all of the designs presented below, we use TiO$_2$ as the metalens material, incorporating its dispersion across visible wavelengths [68]. For the geometrical degrees of freedom, we allow the density of TiO$_2$ to vary between 0 and 1 at every point in the structure (a “topology optimization” approach), and then add penalty functions [12] to Eq. (8) to enforce a binary-material constraint (ultimately converging to densities of 0 and 1 at every point). In some cases we also enforce the constraint of constant permittivity in one ($z$) direction, akin to a slab that can be fabricated by conventional lithography techniques. The initial density parameters are determined based on random distribution within the range between 0.45 and 0.55.

The critical step in inverse design is the efficient computation of the gradient with respect to the arbitrarily many geometrical degrees of freedom, to discover efficient updates from one geometry to the next. Instead of *independently* simulating every geometrical perturbation, by reciprocity (or its generalization for nonreciprocal media [69]), one can use a single coherent set of dipole sources, with spatially dependent amplitudes determined by the form of the figure of merit, to compute an “adjoint” field that in one simulation provides information about all possible perturbations [70]. For a metric with the form of Eq. (7), the current sources for the adjoint simulation are given by $\boldsymbol {J}_\textrm {adj} = -i\omega \boldsymbol {P}_\textrm {adj} = -i\omega \partial \mathcal {F} / \partial \boldsymbol {E}$ (Ref. [14]), which on the exit plane of the metalens are given by:

The plane-wave-decomposition approach to metalens inverse design is depicted in Fig. (a). We simulate every geometry with Meep, a free-software implementation [71] of the finite-difference time-domain (FDTD) method [72]. We consider broad-bandwidth response for 450–700nm wavelengths, representing a large 43% relative bandwidth. We optimize the fields at 20 wavelengths in this range and verify that the optimal devices have a smooth spectral response at substantially higher resolution (120 wavelengths). (We found that optimizing with 10 frequency points or fewer may lead to highly oscillatory frequency response with wide deviations within the visible-frequency bandwidth at intermediate, non-optimized frequencies, as shown in the Appendix.) The maximin metric of Eq. (8) is not differentiable everywhere, caused by by frequency crossings whereby infinitesimal geometric perturbations cause the lowest-efficiency frequency to change. This can be handled by a standard optimization transformation to epigraph form [73], in which a dummy scalar variable is maximized, and the value of $\mathcal {F}$ in Eq. (7) at each frequency becomes an inequality constraint, but for metalenses that are many wavelengths in size, with many frequencies to be computed, this can require substantial data storage. Instead, for every geometry we select only the minimum-efficiency frequency’s gradient and use that as the gradient for the more general FOM, which is the correct gradient everywhere that the FOM is differentiable. In theory such an approach could lead to oscillations and slow convergence, but in practice we see relatively good convergence to high-quality bandwidth-averaged response (as shown in the Appendix. Figure 15). The dummy-variable approach may lead to even better performance, though at the expense of much greater memory requirements. With our approach, each iteration takes approximately 45 seconds on 25 cores in our computational cluster (Intel Xeon E5-2660 v4 3.2 GHz processors). The least-squares figure of merit rapidly improves then converges in about 200 iterations, after which the geometric penalization transforms the grayscale design to a binary one, which is a slow process (on the order of 2000 iterations) but which results in little-to-no efficiency losses. By further improving the penalization process, the total number of iterations could be reduced to a few hundred.

## 4. Achromatic metalens design

In this section, we design metalenses across a range of low to high numerical apertures, achieving relatively high efficiencies for each optimal device. We assume two-dimensional devices with thicknesses of 250 nm and widths of 12.5 $\mu$m, with the latter chosen to be large enough ($\approx 20\lambda$) to require large phase variations but small enough for relatively fast optimizations. In the Appendixwe demonstrate 2D designs with similar efficiencies for widths up to 60$\lambda$, as well as fully 3D metalens designs. For simplicity we do not incorporate a particular substrate, though again in the Appendix we demonstrate similar performance for an optimized device on a substrate. We use TiO$_2$ for the material, and incorporate its material dispersion [68] by fitting a Lorentz–Drude model to the susceptibility. A transverse magnetic wave is used for the incident field. For the geometric “pixels,” we enforce a 25 nm minimum size to avoid non-robust, highly sensitive designs that are difficult to fabricate (validation with a finer grid spacing shown in Fig. 16). We also verify robustness after optimization by simulating material imperfections in the optimized designs, which retain high efficiency even for moderate geometrical imperfections.

We have designed optimal metalenses at each of the 10 numerical apertures $0.1, \ldots , 0.9, 0.99$, with their efficiencies plotted in Fig. (d). We include further characteristic data for each in the Appendix; in the remainder of the section, we highlight and discuss three of the designs, at NA = 0.1, 0.9, and 0.99. All geometrical degrees of freedom for each design are included in the Appendix. The designs must achieve exquisite control of phase and amplitude over a broad bandwidth in a multiple-scattering physical system; accordingly, it is difficult to physically interpret the real-space profiles of the optimal structures.

#### 4.1 Low-NA, high-efficiency metalens

All broadband metalens designs to date [2–5,8,9,21] operate in the low-NA regime. Thus, we start with low-NA designs for a more direct comparison, and we demonstrate higher efficiencies than the current state-of-the-art.

For 12.5 $\mu$m width and 250 nm thickness, a 0.1-NA device has a 60 $\mu$m focal length, depicted in Fig. 3(a). As shown in Fig. 3(b), the designed device achieves an average focusing efficiency of 65% and a maximum efficiency of 78% (minimum efficiency of 43%) over visible wavelengths (450 – 700 nm). We use the same definition of focusing efficiency as Ref. [8]: the ratio of the power concentrated within the region out to the first electric-field minimum divided by the total incident power (*including* all reflection loss). The focal length remains nearly constant over wavelength. Any focal-length deviations remain within the depth of focus as measured by the full-width half-max of the field intensity. As shown in Fig. 3(c), the focused electric fields are nearly diffraction-limited, and are even narrower than the diffraction limit at shorter wavelengths. (Imperfect focusing efficiencies enable sub-diffraction-limited fields [64,74].) Figure 3(d) shows the optimized structure, comprising an alternating series of TiO$_2$ blocks and air holes. Considering the high-level geometric patterning, the density of TiO$_2$ increases from the lens edge to its center, as also seen in unit-cell-based designs [8,9]. On the other hand, very distinguishable from unit-cell designs, the small-scale geometrical variations are clearly non-periodic. Figure 3(e) shows the electric-field intensity profiles over many wavelengths, exhibiting clear high focusing efficiency; despite significant variations in the near-zone fields as the wavelength changes, each wavelength depicted (and all intermediate wavelengths not shown) produce focusing at the same focal plane. The images in Fig. 3 show the entire simulation region with no stray fields unaccounted for. Further increases in focusing efficiency would be possible with materials that have less dispersion than TiO$_2$. We enforced relatively large minimum features sizes (25 nm), for fabrication compatibility; smaller feature sizes may yield modest further efficiency improvements.

#### 4.2 High-NA metalens, NA = 0.9

In this subsection we design and demonstrate high-NA broadband achromatic metalenses. The optimal device, shown in Fig. 4(a), has the same dimensions as the NA=0.1 lens while the target focal length is now 3 $\mu$m. The optimal-device focal length, shown in Fig. 4(b), is nearly unchanged over visible wavelengths, and the focusing efficiency ranges from 13% to 32% with an average value of 23%. As expected, the focusing efficiency is smaller than for the low-NA device, due to the increased difficulty of focusing light to a closer point over the same frequency bandwidth. The FWHM of the optimized metalens is slightly larger than that of the diffraction limited lens. As shown in Fig. 4(d), the optimized device structure is again clearly non-periodic, with rapid changes in topology that accommodate the $\approx 12$ intervals of 0 to $2\pi$ phase change required at these dimensions. Normalized intensity ($|E|^{2}$) profiles are shown in Fig. 4(e) for 8 selected wavelengths. The 450-nm-wavelength intensity profile (blue) shows a primary focal spot at the desired length of 3$\mu$m as well as a second focal spot at 1.9$\mu$m, which arises due to the relatively low focusing efficiency at that wavelength. To further improve efficiency of the high NA achromatic metalens, as we discussed in the previous section, one could use a much finer resolution of the structure up to fabrication limit. This is illustrated by the fact that increasing the degrees of freedom with a “freeform” topology varying in the $z$ direction improves the average focusing efficiency from 23% to 41% for 0.9 NA metalenses, as shown in Fig. (d).

#### 4.3 Very-high-NA freeform metalens, NA=0.99

In this section, we probe the extreme limit of high-NA design, relaxing the constant-$z$ topology constraint and designing a freeform-topology device that achieves NA = 0.99. The optimized device has 12.5 $\mu$m width, 500 nm thickness, and a focal length now of 0.9 $\mu$m, as depicted in Fig. 5(a). The focal length shown in Fig. 4(b) is again nearly constant over visible frequencies and the average focusing efficiency is 27% and minimum efficiency of 23%. This result shows much higher efficiency compared to the same NA constant-$z$ structure metalens due to the increased degrees of freedom. The FWHM of the optimized metalens is slightly larger than that of the diffraction limited lens. As shown in Fig. 5(d), the optimized structure exhibits diagonal patterns that presumably help redirect the light towards the nearby focal spot. Normalized intensity ($|$E$|^{2}$) profile is shown in Fig. 5(e) for 8 selected wavelengths. It shows clear focal spots at the target focal plane (0.9 $\mu$m) over the visible spectrum.

## 5. Conclusions

In this work, we have demonstrated the capability for inverse design to discover high-efficiency achromatic metalenses across the visible spectrum. We focused on 2D devices, to enable rapid systematic design of devices for ten different numerical apertures and two topologies (const-$z$ and freeform), as well as the many optimizations with varying initial conditions and hyperparameter choices for each case. In the Appendix we demonstrate two generalizations: to larger diameters (up to 60$\lambda$) with very similar efficiencies, and to similarly sized 3D metalens devices. Our results corroborate those from the literature demonstrating that there is not a significant dropoff from 2D device designs [2,3,75] to their 3D counterparts [8,9,21,28].

Through inverse design, we have demonstrated the highest-efficiency low-NA achromatic metalenses to date, as well as the first theoretical demonstration of broadband high-NA structures. Breaking the local-periodicity assumptions of the unit-cell design approach results in significant enhancements in device efficiency and a new regime for performance. Scaling this approach to macroscopic length scales should be possible with emerging fast-solver techniques. A natural concern with increasing the size scale will be the scalability of inverse design, yet the efficiency of giga-scale adjoint-variable gradient-based optimization has been proven in two other domains: deep learning, where networks with billions of parameters are trained efficiently, and aerodynamic design [76], where optimal distributions of more than a billion voxels can be achieved. Macroscopic metalens inverse design would offer significant advantages for applications including wearable optical devices, microscopy and integrated optics.

## 6. Appendix

## 6.1. Material-dispersion simulation verification

In this work, Finite Difference Time Domain (FDTD) method is used to solve a full Maxwell equation [71,72]. FDTD can cover a wide range of frequencies with a single time domain simulation. However, it requires time-domain modeling of material dispersion. We model the dispersion of TiO$_2$ material with Drude model and validate the material modeling with a theoretical calculation of Fabry-Perot oscillations in a 1D dielectric slab as shown in Fig. 6.

## 6.2. Achromatic metalenses with a glass substrate

Metalenses are generally fabricated on transparent substrates. Here, we demonstrate achromatic metalenses with SiO$_2$ substrate (n = 1.5). The incidence power is calculated within the glass substrate. As shown in Fig. 7(a), the optimized device has 250 nm thickness, 12.5 $\mu$m width, and 8.2 $\mu$m focal length, which corresponds to NA = 0.6. As shown in Fig. 7(b), a calculated average focusing efficiency is 32 % and minimum efficiency is 19% over visible wavelengths (450 – 700 nm) while the highest efficiency is 39 %. The focusing efficiency is defined as the ratio of the E-field intensity within the first minimum point and the incidence power [8]. As shown in Fig. 7(c), the calculated full-width half-maximum (FWHM) over visible wavelength is very close to the diffraction limited FWHM except for the longer wavelength. Figure 7(d) shows the optimized structure which is feasible to existing lithography techniques [29,77]. Figure 7(e) shows E-field intensity profiles at the focal plane. Except for the 700nm wavelength, they are very close to diffraction limited focal spots.

## 6.3. Least-squares optimal coefficients

In the main text, we solve for optimal unit-cell coefficients to achieve a desired field distribution. In the setup of the problem, there is a given field $\boldsymbol {E}(\boldsymbol {x})$ that can be decomposed into an orthonormal basis of functions $\boldsymbol {E}_i(\boldsymbol {x})$ and coefficients $c_i$, and a target field $\boldsymbol {E}_\textrm {tar}(\boldsymbol {x})$ that one would like $\boldsymbol {E}(\boldsymbol {x})$ to most closely approximately, i.e., there is a functional $F$ of the coefficients $c_i$ to be minimized:

For any given $i$, we can compute the value of $c_i$ for which $\partial F / \partial c_i$ equals 0. In fact, we must carefully consider both the real and imaginary parts of $c_i$; instead, we can use the complex CR calculus [78], whereby $c_i$ and its conjugate ${c_i}^{*}$ are formally considered independent variables, and the derivatives of $F$ with respect to each must be zero. Since $F$ is real-valued the two derivatives contain redundant information, and only one must be set to zero; for simplicity in the ultimately derivation, we choose the derivative with respect to ${c_i}^{*}$:

## 6.4. Summary of optimized metalenses for various numerical apertures using translation invariant geometry

We demonstrated low to high NA metalenses in the main text. Here, we attach detailed optimized structures and efficiencies at the visible wavelength. As shown in Fig. 8, the optimized structures, in the macroscopic view, have a gradually varying TiO$_2$ filling ratio over device radial direction. The speed of this variation increases over increasing NA. Figure 9 shows achieved efficiencies for various NA (0.1 – 0.9) using translation invariant geometry.

## 6.5. Summary of optimized metalenses for various numerical apertures using freeform geometry

The numerical aperture (NA) versus efficiencies plot is shown in the main text. The optimized geometries and efficiencies data are shown in Figs. 10 and 11. It might be hard to figure out any physical insight from the optimized geometries. Again, in the macroscopic view, they seem to have a gradually varying TiO$_2$ filling ratio over device radial direction.

## 6.6. Analysis of larger device size

Large area design is a particular problem in the metalens community [67]. To address this issue, we optimize metalenses with different device widths from 10 $\lambda$ to 60 $\lambda$. NA is fixed to 0.3 and then the device width is adjusted to see the focusing efficiency variation. The optimized structure is translational invariant geometry with the 250-nm-thick TiO$_2$ material. As shown in Fig. 12, the average efficiencies are close to 40 % for different device sizes. There is no sign of significant efficiency drop for larger devices. Note that it was tested up to the 60 $\lambda$ long device where our computational capacity can handle.

## 6.7. Frequency sensitivity study

Early metalens works have demonstrated discrete operating wavelengths [2,3,29]. Then, in order to achieve a pure achromatic operation, all of phase profile, group delay, and group delay dispersion are modeled by the unit cell approach [8]. Also, densely selected wavelengths [4,9,21] could also realize a continuous achromatic metalenses. The former should work well for a certain bandwidth centered at the selected wavelength. However, many studies using the latter method do not check the performance over continuous bandwidth. Here, we optimize an achromatic metalens with 10 discrete wavelengths ranging from 450 nm to 700 nm for NA = 0.99, which corresponds to $\Delta \lambda =$ 35 nm. As shown in Fig. 13, the selected ten wavelengths show 31 % average efficiency and minimum efficiency of 20% while the average efficiency calculated with finer $\Delta \lambda$ is only about 18 % and minimum efficiency of 1%. This implies that many existing achromatic metalenses which were evaluated at well separated wavelengths may only work at those wavelengths.

## 6.8. Metalens design in 3D

In the main text, we demonstrated 2D achromatic metalenses due to the limit of our computational resource. Here, we present a relatively small 3D achromatic Metalens to validate our methodology in 3D. The diameter of the lens is 7.5 $\mu$m and thickness is 250 nm. The focal length is 15.0 $\mu$m which corresponds to NA = 0.25. Figure 14(a) shows the optimized 3D achromatic metalens in top view. The black region indicates TiO$_2$ and the white region indicates air. The focal length remains nearly constant at three target frequencies (450, 650, 850 nm). The electric-field intensity profiles at the focal plane (z = 15.0 $\mu$m) is shown in Fig. 14(b). Note that this is a preliminary result of designing 3D achromatic metalens. To further improve this result, it may require (1) smaller feature sizes (2) a better angular resolution in order decomposition (3) a good guess on initial geometry parameters. The two formers can be resolved with a greater computing resource while the latter may need theoretical study [57].

## 6.9. Figure of merit convergence over inverse design iteration

## 6.10. Mesh-convergence study

## Funding

Air Force Office of Scientific Research (FA9550-17-1-0093).

## Disclosures

The authors declare no conflicts of interest.

## References

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

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

**3. **M. Khorasaninejad, F. Aieta, P. Kanhaiya, M. A. Kats, P. Genevet, D. Rousso, and F. Capasso, “Achromatic metasurface lens at telecommunication wavelengths,” Nano Lett. **15**(8), 5358–5362 (2015). [CrossRef]

**4. **O. Avayu, E. Almeida, Y. Prior, and T. Ellenbogen, “Composite functional metasurfaces for multispectral achromatic optics,” Nat. Commun. **8**(1), 14992 (2017). [CrossRef]

**5. **S. Shrestha, A. C. Overvig, M. Lu, A. Stein, and N. Yu, “Broadband achromatic dielectric metalenses,” Light: Sci. Appl. **7**(1), 85 (2018). [CrossRef]

**6. **B. H. Chen, P. C. Wu, V.-C. Su, Y.-C. Lai, C. H. Chu, I. C. Lee, J.-W. Chen, Y. H. Chen, Y.-C. Lan, C.-H. Kuan, and D. P. Tsai, “GaN metalens for pixel-level full-color routing at visible light,” Nano Lett. **17**(10), 6345–6352 (2017). [CrossRef]

**7. **R. Paniagua-Dominguez, Y. F. Yu, E. Khaidarov, S. Choi, V. Leong, R. M. Bakker, X. Liang, Y. H. Fu, V. Valuckas, L. A. Krivitsky, and A. I. Kuznetsov, “A metalens with a near-unity numerical aperture,” Nano Lett. **18**(3), 2124–2132 (2018). [CrossRef]

**8. **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]

**9. **S. Wang, P. C. Wu, V.-C. Su, Y.-C. Lai, M.-K. Chen, H. Y. Kuo, B. H. Chen, Y. H. Chen, T.-T. Huang, J.-H. Wang, R.-M. Lin, C.-H. Kuan, T. Li, Z. Wang, S. Zhu, and D. P. Tsai, “A broadband achromatic metalens in the visible,” Nat. Nanotechnol. **13**(3), 227–232 (2018). [CrossRef]

**10. **M. P. Bendsøe and O. Sigmund, “Material interpolation schemes in topology optimization,” Arch. applied mechanics **69**(9-10), 635–654 (1999). [CrossRef]

**11. **M. P. Bendsøe, “Topology optimization,” in * Encyclopedia of Optimization*, (Springer, 2001), pp. 2636–2638.

**12. **M. M. Neves, O. Sigmund, and M. P. Bendsøe, “Topology optimization of periodic microstructures with a penalization of highly localized buckling modes,” Int. J. for Numer. Methods Eng. **54**(6), 809–834 (2002). [CrossRef]

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

**14. **O. D. Miller, “*Photonic design: From fundamental solar cell physics to computational inverse design*,” Ph.D. thesis (University of California, Berkeley, 2012).

**15. **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]

**16. **P. Camayd-Munoz and A. Faraon, “Scaling laws for inverse-designed metadevices,” in * Conference on Lasers and Electro-Optics*, (Optical Society of America, 2018), pp. FF3C–7.

**17. **A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vučković, “Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer,” Nat. Photonics **9**(6), 374–377 (2015). [CrossRef]

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

**19. **L. Yang, A. V. Lavrinenko, J. M. Hvam, and O. Sigmund, “Design of one-dimensional optical pulse-shaping filters by time-domain topology optimization,” Appl. Phys. Lett. **95**(26), 261101 (2009). [CrossRef]

**20. **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]

**21. **N. Mohammad, M. Meem, B. Shen, P. Wang, and R. Menon, “Broadband imaging with one planar diffractive lens,” Sci. Rep. **8**(1), 2799 (2018). [CrossRef]

**22. **G. Zheng, H. Mühlenbernd, M. Kenney, G. Li, T. Zentgraf, and S. Zhang, “Metasurface holograms reaching 80% efficiency,” Nat. Nanotechnol. **10**(4), 308–312 (2015). [CrossRef]

**23. **X. Ni, A. V. Kildishev, and V. M. Shalaev, “Metasurface holograms for visible light,” Nat. Commun. **4**(1), 2807 (2013). [CrossRef]

**24. **W. Zhao, B. Liu, H. Jiang, J. Song, Y. Pei, and Y. Jiang, “Full-color hologram using spatial multiplexing of dielectric metasurface,” Opt. Lett. **41**(1), 147–150 (2016). [CrossRef]

**25. **A. Arbabi, E. Arbabi, Y. Horie, S. M. Kamali, and A. Faraon, “Planar metasurface retroreflector,” Nat. Photonics **11**(7), 415–420 (2017). [CrossRef]

**26. **R. Palmeri, M. T. Bevacqua, A. F. Morabito, and T. Isernia, “Design of Artificial-Material-Based Antennas Using Inverse Scattering Techniques,” IEEE Trans. Antennas Propag. **66**(12), 7076–7090 (2018). [CrossRef]

**27. **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]

**28. **S. M. Kamali, E. Arbabi, A. Arbabi, Y. Horie, and A. Faraon, “Highly tunable elastic dielectric metasurface lenses,” Laser Photonics Rev. **10**(6), 1002–1008 (2016). [CrossRef]

**29. **M. Khorasaninejad, A. Y. Zhu, C. Roques-Carmes, W. T. Chen, J. Oh, I. Mishra, R. C. Devlin, and F. Capasso, “Polarization-insensitive metalenses at visible wavelengths,” Nano Lett. **16**(11), 7229–7234 (2016). [CrossRef]

**30. **L. S. Pontryagin, V. G. Boltyanskii, R. V. Gamkrelidze, and E. F. Mishechenko, * The Mathematical Theory of Optimal Processes* (John Wiley & Sons, 1962).

**31. **J. Céa, A. Gioan, and J. Michel, “Quelques résultats sur l’identification de domaines,” Calcolo **10**(3-4), 207–232 (1973). [CrossRef]

**32. **O. Pironneau, * Optimal Shape Design for Elliptic Systems* (Springer-Verlag, 2012).

**33. **M. B. Giles and N. A. Pierce, “An Introduction to the Adjoint Approach to Design,” Flow, Turbul. Combust. **65**(3/4), 393–415 (2000). [CrossRef]

**34. **S. W. Director and R. A. Rohrer, “The Generalized Adjoint Network and Network Sensitivities,” IEEE Trans. Circuit Theory **16**(3), 318–323 (1969). [CrossRef]

**35. **A. Jameson, “Aerodynamic design via control theory,” J. Sci. Comput. **3**(3), 233–260 (1988). [CrossRef]

**36. **M. P. Bendsoe and O. Sigmund, * Topology optimization: theory, methods, and applications* (Springer Science & Business Media, 2013).

**37. **M. Demiralp and H. Rabitz, “Optimally controlled quantum molecular dynamics: A perturbation formulation and the existence of multiple solutions,” Phys. Rev. A **47**(2), 809–816 (1993). [CrossRef]

**38. **H. A. Rabitz, M. M. Hsieh, and C. M. Rosenthal, “Quantum Optimally Controlled Transition Landscapes,” Science **303**(5666), 1998–2001 (2004). [CrossRef]

**39. **P. J. Werbos, * The Roots of Backpropagation* (John Wiley & Sons, Inc., 1994).

**40. **D. E. Rumelhart, G. E. Hinton, and R. J. Williams, “Learning representations by back-propagating errors,” Nature **323**(6088), 533–536 (1986). [CrossRef]

**41. **Y. LeCun, B. Boser, J. S. Denker, D. Henderson, R. E. Howard, W. Hubbard, and L. D. Jackel, “Backpropagation applied to handwritten zip code recognition,” Neural Comput. **1**(4), 541–551 (1989). [CrossRef]

**42. **L. H. Frandsen, Y. Elesin, L. F. Frellsen, M. Mitrovic, Y. Ding, O. Sigmund, and K. Yvind, “Topology optimized mode conversion in a photonic crystal waveguide fabricated in silicon-on-insulator material,” Opt. Express **22**(7), 8525–8532 (2014). [CrossRef]

**43. **C. M. Lalau-Keraly, S. Bhargava, O. D. Miller, and E. Yablonovitch, “Adjoint shape optimization applied to electromagnetic design,” Opt. Express **21**(18), 21693–21701 (2013). [CrossRef]

**44. **L. Su, A. Y. Piggott, N. V. Sapra, J. Petykiewicz, and J. Vuckovic, “Inverse design and demonstration of a compact on-chip narrowband three-channel wavelength demultiplexer,” ACS Photonics **5**(2), 301–305 (2018). [CrossRef]

**45. **D. Sell, J. Yang, S. Doshay, R. Yang, and J. A. Fan, “Large-angle, multifunctional metagratings based on freeform multimode geometries,” Nano Lett. **17**(6), 3752–3757 (2017). [CrossRef]

**46. **F. Callewaert, V. Velev, P. Kumar, A. Sahakian, and K. Aydin, “Inverse-designed broadband all-dielectric electromagnetic metadevices,” Sci. Rep. **8**(1), 1358 (2018). [CrossRef]

**47. **H. Men, K. Y. Lee, R. M. Freund, J. Peraire, and S. G. Johnson, “Robust topology optimization of three-dimensional photonic-crystal band-gap structures,” Opt. Express **22**(19), 22632–22648 (2014). [CrossRef]

**48. **V. Ganapati, O. D. Miller, and E. Yablonovitch, “Light trapping textures designed by electromagnetic optimization for subwavelength thick solar cells,” IEEE J. Photovoltaics **4**(1), 175–182 (2014). [CrossRef]

**49. **B. Shen, P. Wang, R. Polson, and R. Menon, “Ultra-high-efficiency metamaterial polarizer,” Optica **1**(5), 356–360 (2014). [CrossRef]

**50. **Z. Lin, V. Liu, R. Pestourie, and S. G. Johnson, “Topology optimization of freeform large-area metasurfaces,” arXiv preprint arXiv:1902.03179 (2019).

**51. **C. Sohl, M. Gustafsson, and G. Kristensson, “Physical limitations on broadband scattering by heterogeneous obstacles,” J. Phys. A: Math. Theor. **40**(36), 11165–11182 (2007). [CrossRef]

**52. **D.-H. Kwon and D. Pozar, “Optimal characteristics of an arbitrary receive antenna,” IEEE Trans. Antennas Propag. **57**(12), 3720–3727 (2009). [CrossRef]

**53. **O. D. Miller, S. G. Johnson, and A. W. Rodriguez, “Shape-independent limits to near-field radiative heat transfer,” Phys. Rev. Lett. **115**(20), 204302 (2015). [CrossRef]

**54. **O. D. Miller, A. G. Polimeridis, M. T. H. Reid, C. W. Hsu, B. G. Delacy, J. D. Joannopoulos, M. Soljačic, and S. G. Johnson, “Fundamental limits to optical response in absorptive systems,” Opt. Express **24**(4), 3329 (2016). [CrossRef]

**55. **O. D. Miller, O. Ilic, T. Christensen, M. T. H. Reid, H. A. Atwater, J. D. Joannopoulos, M. Soljačic, and S. G. Johnson, “Limits to the optical response of graphene and two-dimensional materials,” Nano Lett. **17**(9), 5408–5415 (2017). [CrossRef]

**56. **H. Shim, L. Fan, S. G. Johnson, and O. D. Miller, “Fundamental Limits to Near-Field Optical Response over Any Bandwidth,” Phys. Rev. X **9**(1), 011043 (2019). [CrossRef]

**57. **G. Angeris, J. Vuckovic, and S. P. Boyd, “Computational bounds for photonic design,” ACS Photonics **6**(5), 1232–1239 (2019). [CrossRef]

**58. **Y. Kiarashinejad, M. Zandehshahvar, S. Abdollahramezani, O. Hemmatyar, R. Pourabolghasem, and A. Adibi, “Knowledge discovery in nanophotonics using geometric deep learning,” Advanced Intelligent Systems1900132 (2019). [CrossRef]

**59. **A. N. Gorban and I. Y. Tyukin, “Blessing of dimensionality: Mathematical foundations of the statistical physics of data,” Philos. Trans. R. Soc., A **376**(2118), 20170237 (2018). [CrossRef]

**60. **Y. Liu and A. H. Barnett, “Efficient numerical solution of acoustic scattering from doubly-periodic arrays of axisymmetric objects,” J. Comput. Phys. **324**, 226–245 (2016). [CrossRef]

**61. **O. P. Bruno and M. Maas, “Shifted equivalent sources and FFT acceleration for periodic scattering problems, including Wood anomalies,” J. Comput. Phys. **378**, 548–572 (2019). [CrossRef]

**62. **T. Isernia and G. Panariello, “Optimal focusing of scalar fields subject to arbitrary upper bounds,” Electron. Lett. **34**(2), 162–164 (1998). [CrossRef]

**63. **H. Liu, Y. Yan, Q. Tan, and G. Jin, “Theories for the design of diffractive superresolution elements and limits of optical superresolution,” J. Opt. Soc. Am. A **19**(11), 2185–2193 (2002). [CrossRef]

**64. **H. Shim, H. Chung, and O. D. Miller, “Maximal free-space concentration of light,” arXiv preprint arXiv:1905.10500 (2019).

**65. **L. N. Trefethen and D. Bau III, * Numerical linear algebra* (SIAM, 1997).

**66. **G. Swanson, “Binary optics technology: the theory and design of multilevel diffractive optical elements,” MIT Linc. Lab. Tech. Rep. 854 (1989).

**67. **S. M. Kamali, E. Arbabi, A. Arbabi, and A. Faraon, “A review of dielectric optical metasurfaces for wavefront control,” Nanophotonics **7**(6), 1041–1068 (2018). [CrossRef]

**68. **E. D. Palik, * Handbook of optical constants of solids*, vol. 3 (Academic press, 1998).

**69. **J. A. Kong, * Theory of electromagnetic waves* (Wiley-Interscience, New York 1975. 348 p. 1975).

**70. **S. G. Johnson, M. Ibanescu, M. Skorobogatiy, O. Weisberg, J. Joannopoulos, and Y. Fink, “Perturbation theory for Maxwell’s equations with shifting material boundaries,” Phys. Rev. E **65**(6), 066611 (2002). [CrossRef]

**71. **A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. D. Joannopoulos, and S. G. Johnson, “Meep: A flexible free-software package for electromagnetic simulations by the fdtd method,” Comput. Phys. Commun. **181**(3), 687–702 (2010). [CrossRef]

**72. **A. Taflove, A. Oskooi, and S. G. Johnson, * Advances in FDTD computational electrodynamics: photonics and nanotechnology* (Artech house, 2013).

**73. **S. Boyd and L. Vandenberghe, * Convex optimization* (Cambridge university press, 2004).

**74. **P. J. Ferreira and A. Kempf, “Superoscillations: Faster than the Nyquist rate,” IEEE Transactions on Signal Process. **54**(10), 3732–3740 (2006). [CrossRef]

**75. **P. Wang, N. Mohammad, and R. Menon, “Chromatic-aberration-corrected diffractive lenses for ultra-broadband focusing,” Sci. Rep. **6**(1), 21545 (2016). [CrossRef]

**76. **N. Aage, E. Andreassen, B. S. Lazarov, and O. Sigmund, “Giga-voxel computational morphogenesis for structural design,” Nature **550**(7674), 84–86 (2017). [CrossRef]

**77. **S. Owa and H. Nagasaka, “Advantage and feasibility of immersion lithography,” J. Micro/Nanolithogr., MEMS, MOEMS **3**(1), 97–104 (2004). [CrossRef]

**78. **K. Kreutz-Delgado, “The complex gradient operator and the cr-calculus,” arXiv preprint arXiv:0906.4835 (2009).