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.

 

Fig. 1. Metalens inverse design. (a) The figure of merit for optical focusing is computed at an exit plane of the metalens (dashed line). (b) Schematic comparison of the unit-cell design approach, whereby a single (ideal) phase is fixed, using an incomplete basis periodic-boundary diffraction orders, versus large-area inverse design, whereby the full device scattering properties are incorporated into the design. (c) Compilation of broad-bandwidth inverse-designed metalenses using “freeform” (red) and constant-$z$ (“Const-z”) geometries, showing relatively high efficiencies for numerical apertures ranging from 0.1 to 0.99. The black circles indicate recently published results [8,9,21] with similar bandwidths, which operate only in the low-NA regime.

Download Full Size | PDF

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.

 

Fig. 2. (a) Upper bounds to the single-frequency focusing efficiency of a metalens designed by a unit-cell approach. (Shown for a 2D metalens with diameter $L = 100\lambda$.) For high numerical apertures, there is a significant efficiency loss incurred by the unit-cell approach, especially for unit-cell periods $\Lambda$ close to but smaller than the wavelength $\lambda$, where there is a single diffraction order ($N_d = 1$) and the best theoretical efficiency may only be 20%. Inset: variation of the maximum efficiency as a function of period, with abrupt discontinuities as new orders appear. (b) Optimal fields, real part (upper) and phase (lower), of the ideal metalens (blue) and the ideal unit-cell design (orange), showing that non-periodic nature of the focusing fields leads to amplitude and phase inaccuracies in the unit-cell approach.

Download Full Size | PDF

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 overestimates 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

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.

 

Fig. 3. NA = 0.1, high-efficiency achromatic metalens for visible wavelengths (450 – 700 nm) using a lithography-compatible “Const-z” geometry. (a) Metalens dimensions. (Not to scale.) (b) Focusing efficiency and focal lengths of the optimized metalens. The average efficiency is 65% over visible spectrum and minimum efficiency of 43%. The efficiency was calculated by power within the first minimum point divided by the incidence power. Calculated focal lengths remain within the depth of focus area at all wavelengths. (c) Full-width half-maximum (FWHM) of the optimized metalens and the corresponding Airy disk. (d) Optimized TiO$_2$-based metalens design. (e) Normalized intensity ($|E|^{2}$) profile. The curves are normalized intensity profiles at the focal plane.

Download Full Size | PDF

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).

 

Fig. 4. High-NA (NA = 0.9) achromatic metalens for visible wavelengths. (a) Metalens dimensions. (b) Focusing efficiency and focal lengths of the optimized metalens. The average efficiency is 23% over visible spectrum and minimum efficiency of 13%. (c) FWHM of the optimized metalens and the corresponding Airy disk. (d) Optimized metalens design. (e) Normalized intensity ($|E|^{2}$) profile everywhere (images) and at the focal plane (curves).

Download Full Size | PDF

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.

 

Fig. 5. Near-unity-NA (NA = 0.99) achromatic metalens achieved by inverse design of a freeform geometry. (a,b) Metalens dimensionsi, focusing efficiency, and focal lengths. (c) FWHM of the optimized metalens and the corresponding Airy disk. (d) Optimal design. (e) Normalized field-intensity profiles.

Download Full Size | PDF

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.

 

Fig. 6. Simulated and calculated reflection and transmission of 100-nm-thick TiO$_2$ dielectric slab. They match very well over the visible spectrum.

Download Full Size | PDF

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.

 

Fig. 7. Achromatic metalens for the visible wavelength (450 – 700 nm) with substrate (n = 1.5) (a) Schematic illustration of metalens dimension. (b) Focusing efficiency of the optimized metalens. The average efficiency is 32 % and minimum efficiency of 19% over the visible spectrum. The efficiency was calculated by the power within the first minimum point divided by the incidence power. (c) Full Width Half Maximum (FWHM) of the optimized metalens and the corresponding airy disk. (d) Optimized TiO$_2$ based metalens structure. The width and thickness of the metalens are 12.5 $\mu$m and 250 nm, respectively. 25 nm grid spacing was used. The maximum aspect ratio is 10. (e) Normalized intensity ($|$E$|^{2}$) profile.

Download Full Size | PDF

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:

(10)$$\min_{c_i} F(c_i) = \min_{c_i} \left\| \boldsymbol{E}_\textrm{tar} - \boldsymbol{E}(\boldsymbol{x}) \right\|_2^{2} = \min_{c_i} \left\| \boldsymbol{E}_\textrm{tar} - \sum_i c_i \boldsymbol{E}_i(\boldsymbol{x}) \right\|_2^{2}.$$

In the main text we stated that the optimal coefficients are given by

(11)$$c_i = \int {\boldsymbol{E}}^{*}_i \cdot \boldsymbol{E}_\textrm{tar}.$$

This statement is a well-known result for linear least squares problems. Here, we provide a short proof in this context.

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}^{*}$:

(12)$$\begin{aligned}\frac{\partial F}{\partial {c_i}^{*}} &={-} \int {\boldsymbol{E}}^{*}_i \cdot \left( \boldsymbol{E}_\textrm{tar} - \sum_j c_j \boldsymbol{E}_j \right)\\ &={-}\int {\boldsymbol{E}}^{*}_i \cdot \boldsymbol{E}_\textrm{tar} + c_i, \end{aligned}$$

where in the latter expression we have used the orthonormality of the basis functions, $\int {\boldsymbol {E}}^{*}_j \cdot \boldsymbol {E}_i = \delta _{ij}$. Thus the condition for zero first derivative requires Eq. (11) to hold; as the only extremum it must also be the minimum, thus proving that it is the solution to the minimization problem.

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.

 

Fig. 8. Optimized metalenses for various NA (0.1 – 0.9) using translation invariant geometry. The thickness and width are 250 nm and 12.5 $\mu$m, respectively. The minimum resolution is 25 nm with the maximum aspect ratio of 10.

Download Full Size | PDF

 

Fig. 9. Achieved efficiency over 450 – 700 nm wavelength range for various NA (0.1 – 0.9) using translation invariant geometry. The red curves are focusing efficiencies (as defined in the main text).

Download Full Size | PDF

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.

 

Fig. 10. Optimized metalenses for various NA (0.1 – 0.9) using freeform geometry. The thickness and width are 500 nm and 12.5 $\mu$m, respectively. The minimum resolution is 25 nm.

Download Full Size | PDF

 

Fig. 11. Achieved efficiency over 450 – 700 nm wavelength range for various NA (0.1 – 0.9) using freeform geometry. The red curves are focusing efficiencies (as defined in the main text).

Download Full Size | PDF

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.

 

Fig. 12. Average focusing efficiencies for different device sizes.

Download Full Size | PDF

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.

 

Fig. 13. Frequency sensitivity study for NA = 0.99 metalens. (a) The optimization was initially performed over 10 discrete wavelengths (blue circled data), then the performance of metalens was plotted over 30 wavelengths (red curve) to check the broadband operation where the efficiency dropped from 31 % (10 wavelengths calculation) to 18 % (30 wavelengths calculation). The minimum efficiency in the 10 wavelengths optimization is 20 % while the minimum efficiency in 30 wavelengths validation is 1 %. (b) Focusing efficiency of NA = 0.1, 0.9, 0.99 metalenses demonstrated in this work. The black dots indicate 20 target wavelengths within the visible spectrum while the solid lines represent evaluation at 120 wavelengths.

Download Full Size | PDF

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

 

Fig. 14. Optimized 3D metalens. Thickness of device is 0.25 um and lens diameter is 7.5 $\mu$m and NA is 0.25 (a) top view of the 3D metalens. (b) E-field intensity at the focal plane.

Download Full Size | PDF

6.9. Figure of merit convergence over inverse design iteration

 

Fig. 15. Figure of merit evolution by iteration over the course of the design process. The intensity is normalized against that of the Airy disc. The steepest descent method might cause FOM values oscillation after they are converged to the local optima.

Download Full Size | PDF

6.10. Mesh-convergence study

 

Fig. 16. Validation of mesh-size accuracy for NA = 0.6 achromatic metalens in freeform. Focusing efficiencies are plotted over visible spectrum for the three mesh sizes shown here, the average focusing efficiences are 63.4%, 63.0%, and 63.3% for $\lambda /20$, $\lambda /40$, and $\lambda /80$ mesh sizes, respectively.

Download Full Size | PDF

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).