## Abstract

Topology-optimized metasurfaces are thin film optical devices that have received much interest because they support ultra-high diffraction efficiencies. An important design consideration is ensuring that devices are insensitive to imperfections arising from realistic fabrication processing. We show that topology-optimized metasurfaces can be made robust by incorporating the performance of geometrically eroded and dilated devices directly into the iterative optimization algorithm. We additionally apply topology optimization to refine conventional designs and make them more robust. Unexpectedly, we find that devices robust to systematic erosion and dilation variations are also robust to random periodic perturbations. An analysis of the optical modes of robust devices indicates that robustness is enforced via highly complex and non-intuitive interactions between these modes and cannot be enforced using simple design rules. These concepts can more generally address other fabrication imperfections, such as thickness and refractive index variation, and can extend to other diffractive and nanophotonic platforms.

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

## 1. Introduction

Metasurfaces are optical devices that can control light with capabilities exceeding those of bulk refractive and scalar diffractive optics [1–4]. They are powerful tools for wavefront control [5–7], imaging [8–10], sensing [11,12], and many other applications in photonics [13,14]. An emergent class of metasurface designs is based on topology optimization, which uses numerical adjoint-based methods to produce devices with unusual geometries and non-intuitive dynamics. This technique has been used to design photonic crystals [15], integrated photonics devices [16,17], and metamaterials [18–24] with superior performance metrics than those based on conventional designs.

A principal challenge in translating metasurface designs from theory to experiment is the fabrication of nanoscale features. Patterning [25–27], thin-film etching, lift-off, and additive manufacturing [28] all introduce inevitable geometric deviations between experimental devices and ideal, theoretical designs. Some of these deviations, such as sidewall roughness, are random. Others, such as sloped sidewalls, over- and under-dosing during patterning, and over- and under-etching, are more systematic. Figure 1(a) shows scanning electron microscope (SEM) images of a representative topology-optimized silicon device layout patterned using two different electron beam doses, which results in geometric variations. The device shown in Fig. 1(b) is more optimally dosed and features small holes approximately 50 nm in diameter, while the device in Fig. 1(c) is clearly underdosed and the small holes disappear. It is not immediately obvious whether these holes are critical to the operation of this device, but such variations are potentially of great concern given the utility of deeply subwavelength features in metasurfaces. It is therefore imperative that devices are designed to be robust and that they maintain high performance even with fabrication process variations.

In this study, we examine methods for the robust design of topology-optimized metasurfaces and the underlying physical principles of robustness. As a model system, we analyze periodic silicon deflectors, termed metagratings, for study. We show that robustness to systematic geometric deviations in the form of geometric erosion and dilation can be incorporated into the iterative optimization process, and we capture these distortions experimentally. This study builds on prior work in our group [18], where we applied such robustness criteria to the design of large angle metagrating deflectors, and it validates our use of erosion and dilation as metrics for robustness. We also show that these concepts can be used to introduce robustness into conventional metasurface designs and can apply to other forms of fabrication imperfections, such as thickness or refractive index variation.

## 2. Modeling fabrication errors

To capture geometric variations in a quantitative manner, we will analyze devices that are geometrically eroded and dilated relative to the ideal pattern. While these geometric perturbations do not capture all possible fabrication variations, they do mimic the effects of certain systematic errors, such as over- and under-dosing during patterning and over- and under-etching. The use of erosion and dilation as robustness metrics was previously demonstrated in the design of topology-optimized photonic crystals [29].

To specify the eroded and dilated geometries of an ideal pattern comprised of two materials with different refractive indices, we first define the refractive index distribution of the ideal pattern, $\rho (r)$. $r$represents position, and the device and background refractive index values are normalized to 1 and 0, respectively. A blurred version of the pattern, $\tilde{\rho}(r)$, is produced by convolving $\rho (r)$with the Gaussian distribution:

The spread of the Gaussian filter, $\sigma $, corresponds to the blur radius and is determined by a length scale in the fabrication process that corresponds with patterning resolution. In our subsequent analysis, we set$\sigma $to 20 nm, which is the nominal beam size in our electron beam lithography tool. The eroded and dilated patterns, $\overline{\rho}(r)$, are produced by thresholding the blurred pattern:

To quantify the amount of erosion and dilation with a single parameter, we define the edge deviation,$\Delta $, which is related to the blurring and thresholding parameters:

The edge deviation is physically intuitive and represents the distance that a semi-infinite flat edge will shift. Devices with negative and positive edge deviations correspond to those that have been eroded and dilated, respectively. To visualize these geometric perturbations, Fig. 2 shows eroded and dilated versions of the device depicted in Fig. 1 with −10 nm and + 10 nm edge deviations, respectively.## 3. Robustness in topology-optimized metasurfaces

In our metasurface topology optimization method, we begin with a random distribution of refractive indices,$\rho (r)$, which comprise a continuum of values bounded by the device and background indices. Our goal is to maximize the Fig. of Merit (FoM), which for metagratings is the absolute efficiency of light deflected to the desired diffraction channel. Absolute efficiency is defined as the intensity of light transmitted to the desired diffraction order divided by the incident light intensity. For each iteration, we perform forward and adjoint electromagnetic simulations to compute the FoM gradient at each pixel of the device:

This gradient specifies perturbative changes in the refractive index at$\rho (r)$that improve the FoM. Constraints are added to the optimization to ensure that the final pattern converges to a binary design comprising only two refractive index values. Details on the optimization can be found in [18].Devices designed in this manner are not intrinsically robust and are highly sensitive to small geometric variations. With a representative topology-optimized metagrating, Figure 3(a) shows the absolute efficiency of the device as a function of edge deviation. The device consists of 325 nm tall, polycrystalline silicon patterns on a fused silica substrate, operates at a wavelength of 1050 nm, and deflects normally-incident TM-polarized light to 75°.

Historically, such high-angle metagratings have proven difficult to design with high efficiencies but are of great interest in areas such as spectroscopy and large-numeric-aperture lenses for imaging [18].

The device exhibits a peak theoretical efficiency of over 97% when ideally patterned (i.e., with zero edge deviation). However, the device efficiency greatly diminishes as the magnitude of edge deviation increases. We experimentally capture these trends in device efficiency as a function of edge deviation by patterning devices with different levels of electron beam lithography dosing. The devices are fabricated on a single sample so all other fabrications parameters such as etch characteristics or material thickness are the same between the devices. The experimental data most closely fit a theoretical edge deviation range from +0.8 to +8.2 nm, which was determined using a least-squares fit between the experimental and simulated efficiency trends. The trends in experimental efficiencies generally match well with dilated theoretical trends and indicate that this particular fabrication batch was slightly overdosed and/or overetched. Representative SEM images of devices patterned with the highest and lowest dosages are presented in Fig. 3(b) and agree well with theoretical models, shown as insets. Key features in the model with small edge deviation, such as gaps and holes, are faithfully reproduced in the experimental sample. Similarly, the model with large edge deviation accurately predicts the loss of small shapes, as found in the experimental device.

Such extreme sensitivity to small geometric deviations highlights the difficulty of reliably fabricating topology-optimized devices: a theoretical precision of ±1 nm is required for this particular metagrating to support efficiencies greater than 80%. Furthermore, we observe a decrease in experimental efficiencies compared to theoretical values for devices fabricated with edge deviations near zero. This decrease is likely due to the sensitivity of the device to other non-systematic geometric perturbations from fabrication. As such, we were unable to produce samples with over 80% experimental efficiency.

## 4. Robustness control

To design devices that are less sensitive to edge deviation, we modify our topology optimization technique to account for systemic geometric errors that may occur during fabrication. In versions of topology optimization that do not account for robustness, gradients to the FoM with respect to refractive index are computed for only the ideal device pattern. In our modified method incorporating robustness, we compute FoM gradients for ideal, eroded, and dilated patterns for each optimization iteration. With this approach, eroded and dilated patterns must be calculated for devices comprising a continuum of refractive index values. As such, we modify our thresholding function used to generate eroded and dilated versions of the device with the following:

A schematic illustrating our incorporation of robustness during a single topology optimization iteration is shown in Fig. 4(a). We begin with the base pattern and first apply blurring and thresholding filters to generate an eroded pattern,${\overline{\rho}}_{e}$, and dilated pattern,${\overline{\rho}}_{d}$, with a desired value of edge deviation. Next, we compute separate gradients,${\overline{G}}_{q}$, for each of the three patterns,${\overline{\rho}}_{q}$. Here, the three patterns are indexed by$q$. At this stage, the three gradients act to improve the FoM for three different patterns. To produce a single gradient of the ideal pattern from these three gradients, we project the gradients from the eroded and dilated patterns to the ideal pattern using the chain rule. The final gradient used to modify the ideal pattern is computed as:

Devices produced with this procedure exhibit much greater robustness to edge deviation variations than their non-robust counterparts. As a proof-of-concept, we designed a metagrating that is robust to edge deviations of ±5 nm with the same design parameters as the device in Fig. 3. The calculated efficiencies of this fully optimized device as a function of edge deviation are plotted in Fig. 4(b), together with those of the non-robust device from Fig. 3(a). The topology-optimized device designed with robustness maintains an efficiency above 90% for edge deviations ranging from −4 to +3 nm, and an efficiency above 80% for edge deviations ranging from −6 to +5 nm. Its peak efficiency at zero edge deviation is 95%, which is lower than the peak efficiency of the non-robust device. This difference in peak efficiency can be attributed to the fact that the incorporation of robustness serves as a design constraint and yields an inherent tradeoff between peak efficiency and robustness. Devices designed to be robust for larger ranges of edge deviation than those used here will be even more robust to geometric perturbations at the expense of further reductions in peak efficiency.

Figure 5(a) shows experimental efficiencies for the robust device in Fig. 4. Similar to the analysis presented in Fig. 3(a), a set of devices is fabricated with varying dosages. The experimental efficiency values match well with the trends of the simulated devices. SEM images of the devices patterned with highest and lowest dosages, corresponding to edge deviations of +1.7 to +8.9 nm, are shown in Fig. 5(b) and exhibit good qualitative agreement with the target, theoretical geometries. The robust device has measured efficiencies as high as 88%, which is 10% higher than the highest efficiency non-robust device. Importantly, devices possessing edge deviations near zero all exhibit high efficiencies with no observable drop-off in efficiency, indicating that these devices are tolerant to systematic errors incurred during fabrication.

## 5. Robust optimization of conventional metasurface designs

While robustness must be specifically engineered into topology-optimized metasurfaces, conventional phased-array metasurfaces that deflect light to modest angles already benefit from some measure of intrinsic robustness. The inset in Fig. 6(a) shows a metasurface designed following conventional design strategies [30]. The device consists of 580 nm tall polycrystalline silicon posts on a glass substrate and deflects 1050 nm, TM-polarized light to 40°. The curve in Fig. 6(a) shows that this device is very robust as a function of edge deviation. This result is not altogether unsurprising because these phase-array designs typically rely on an adiabatic sampling of the desired phase profile, where adjacent posts are similar in size.

However, such intrinsic robustness does not necessarily hold as the deflection angle increases and the size difference between adjacent posts becomes exaggerated. In this regime, we are no longer guaranteed intrinsic robustness to edge deviations. The red inset in Fig. 6(b) shows a metasurface designed to deflect light to 75° with the same operating wavelength and device thickness. The red curve in Fig. 6(b) shows the corresponding robustness characteristics of the device. While the device shows a nominal 81% peak efficiency, edge deviations of only ±3 nm cause the efficiency to quickly plummet below 40%. This device is not robust to erosion and dilation. In order to make this device more robust, we can apply topology optimization with robustness controls.

We begin by applying a blurring filter to the conventionally designed posts, which prevents the optimization from getting caught in the local optimum. Then, we apply one hundred steps of topology optimization with robustness enforcement to the blurred structure. Such a procedure requires around one-third to one-fourth of the computational resources as a full topology optimization. The resulting topology optimized structure is shown as the blue inset in Fig. 6(b). The blue curve shows the corresponding robustness characteristics of the refined structure. Not only is the resulting structure notably more robust to edge deviations, the peak efficiency of 85% is also higher than that of the original structure.

While our previous studies have shown that using conventional structures as starting points typically results in worse performance than fully randomized starting points [21], such a procedure is still of interest due to the high computational cost of topology optimization. With this framework, previously designed metasurfaces can be made more robust by applying an abridged topology optimization step.

## 6. Robustness under random local deformations

In this section, we examine how devices designed with robustness criteria respond to random, local perturbations. Such deformations mimic more stochastic fabrication errors, such as sidewall roughness stemming from imperfect liftoff, resist tail residue, or material resputtering during etching. For our analysis, we will examine metagratings with random perturbations added to the grating unit cell, meaning that these perturbations are periodic in nature. While these perturbations are not random across an entire macroscopic device, as would be expected in an experimental system, they serve as an approximation to realistic sidewall roughness and are computationally tractable to analyze. This type of analysis is similar to those previously performed on mechanical systems [31].

Returning to our topology optimized devices, we introduce local perturbations as random elastic deformations, using a formalism described in [32]. First, the coordinates of each pixel in the unit cell are displaced to a random nearby position, with the coordinate transformation represented by a displacement vector field. The displacement field is then convolved with a Gaussian filter to ensure that the resulting deformations are smooth. This process is illustrated in Fig. 7(a) using a semi-infinite edge as our ideal pattern. Box I shows a random displacement field with an average displacement of 6 nm and a Gaussian filter with spread of 10 nm. For the case of a semi-infinite edge, only the displacement field near the edge contributes to deformations, as shown in Box II. This coordinate transformation results in the elastic deformation shown in Box III. The dashed red lines indicate 6 nm displacement values from the ideal edge. To further visualize these elastic deformations, Fig. 7(b) displays circular shapes that have been transformed using differing displacement values. At high displacement values, new topological features, such as holes and islands, begin to appear.

To analyze the impact of local deformations on device performance, we apply random deformations to our non-robust device from Fig. 3 and our robust device from Fig. 5. We consider absolute displacement values ranging from 2 to 10 nm, and one hundred different random deformations are applied to each device for each displacement value. Histograms of the resulting device efficiencies are plotted in Fig. 7(c). The histograms of the non-robust devices indicate that this metagrating design is highly sensitive to random perturbations. For displacements of only 2 nm, the distribution already shows a broadening of efficiency values, albeit with values mostly above 80%. For larger displacements, the distribution broadens substantially, and for 8 nm displacements, the majority of devices have efficiencies below 80%.

Unexpectedly, we find that the device designed to be robust to erosion and dilation perturbations is also robust to random perturbations. For displacements of 2 and 4 nm, the efficiency distributions are narrowly peaked, and nearly 90% of the devices with 4 nm displacement have efficiencies over 90%, which is within 5% of the value of the ideal device. The efficiency histogram begins to visibly broaden for a displacement of 6 nm. This displacement value is near the edge deviation value used for robustness design. For these deformations, the median device efficiency is 91%, and 90% of the devices still have efficiencies over 80%. While the efficiency histograms continue to broaden for increasing displacement values, they still comprise a narrower distribution with higher values compared to those of the non-robust device. These histograms suggest that devices designed to be robust for a given range of edge deviations are also robust to random deformations, provided that the random deformation length scales are smaller than the edge deviation parameter.

## 7. Mechanisms for robustness

To provide insight into the mechanisms for robustness, we examine the optical modes of the robust metagrating design from Fig. 5 under different geometric perturbations. We utilize coupled Bloch mode analysis to solve for the optical modes of the device and their relative contributions to device efficiency. Details pertaining to this calculation can be found in [33,34]. To summarize, we treat the metagrating as a vertical Fabry-Perot cavity that supports *n* propagating Bloch modes. The efficiency contribution of the *i*^{th} mode (ordered in decreasing n_{eff}) is computed as the difference in diffractive efficiency of the device utilizing the first *i* and *i-*1 modes.

One may naively expect that robust devices have optical modes that are relatively insensitive to geometric perturbations. Our mode analysis shows that this is not the case. Figure 8(a) shows how the efficiency contribution from the different Bloch modes varies as a function of edge deviation, for our robust metagrating. It is evident that in the regime of robust operation (i.e., edge deviations from −4 to +3 nm), the efficiency contribution from each mode varies for differing edge deviations. In addition, n_{eff} of each mode varies as a function of edge deviation, in a manner that is highly mode dependent as shown in Fig. 8(b). This variation can be significant, as we observe that the mode with the lowest n_{eff} only exists for edge deviations larger than −3 nm and is otherwise cut off. These variations in mode characteristics as a function of edge deviation can be attributed to the strongly varying structure of the modes, shown in Fig. 8(c), which range from highly confined in the silicon to higher order and air-like. This analysis indicates that robust devices support multiple modes that possess special and optimized combinations of refractive indices and scattering properties as a function of edge deviation. While these modes evolve and interact in different ways for different edge deviations, they nonetheless enable high efficiency diffraction into the desired diffraction order.

We also perform a similar analysis for the robust metagrating undergoing random, local deformations. Figure 9(a) shows the efficiency contribution from the different Bloch modes for one hundred different locally-deformed devices. These devices are taken from the previous section and have an average displacement of 6 nm. This plot indicates that while the total device efficiency is high for most devices, there exist variations in mode contribution that are a function of the detailed geometric perturbations. The scope of random, geometric variation can be visualized in Fig. 9(b), which shows unit cell images of the top twelve highest efficiency devices from the distribution. This analysis suggests that these devices are robust to local geometric deformations for similar reasons they are robust to systematic geometric deformations: while the modes evolve and interact in different ways for different local perturbations, they work together in a manner that enables high efficiency diffraction.

## 8. Additional forms of robust optimization

In all the above analyses, we show that topology optimization allows us to incorporate robustness to specific geometric variations in the form of erosions and dilations of the pattern. More generally, topology optimization allows us to design gradients that enable robustness against other forms of fabrication imperfections. As an immediate extension, robustness against device layer thickness variations can be optimized using the following gradient:

_{${\overline{G}}_{t}$}are the gradients for devices with smaller, ideal, and larger thicknesses and${w}_{t}$represent the corresponding weights.

Additionally, this technique can be applied to robustness beyond purely geometric fabrication variations. As an example, we explore robustness to variations in the refractive index of the dielectric material, which can arise from variations in material density, porosity, and crystallinity. For silicon, these dielectric properties strongly vary as the material ranges from amorphous to polycrystalline and monocrystalline [7]. The corresponding gradient is:

In this case, the various${\overline{G}}_{i}$ give the gradients for devices with refractive indices${n}_{i}$that are lower, ideal, and higher than the target index${n}_{0}$.As a point of reference, we analyze the performance of our prior topology-optimized devices undergoing ±10% changes in refractive index. The target refractive index for our devices is ${n}_{0}$ = 3.57 for polycrystalline silicon at a wavelength of 1050 nm, yielding an index range from 3.2 to 3.9. Given the well-characterized nature of silicon deposition techniques, this range of index variation is not typically a concern during device processing. In other common dielectric material choices such as TiO_{2}, however, there is less control in materials processing, and variations in refractive index as large as 20% are observed [35].

Figure 10(a) plots the absolute efficiencies as a function of refractive index for the two topology-optimized metagrating deflectors from Figs. 3 and 5. Both devices have limited resilience to changes in refractive index, with comparable sensitivities to refractive index variations.

Now, we use topology optimization to design silicon metagratings that are robust to refractive index variations. Here, we compute gradients for three devices with the same physical pattern, but with changes of −5, 0 and +5% in refractive index, corresponding to n = 3.39, 3.57 and 3.75, respectively. The device was optimized for 350 iterations from a random starting point. The performance of the final device as a function of refractive index is given in Fig. 10(b). The device exhibits an absolute efficiency exceeding 80% across an index range of n = 3.3 to 3.8, doubling the range of devices designed without refractive index robustness.

Edge, thickness, and refractive index deviations are just three examples of fabrication imperfections. The true power of topology optimization is that we can utilize it to design robustness against a broad range of parameters, so long as a proper gradient can be defined. Furthermore, topology optimization is not limited to only one robustness objective at a time, as robustness to multiple objectives can be generally achieved by summing up the individual gradients of each objective. As an example, a three-parameter robust gradient would take the following form:

*n*robustness parameters, this optimization would require simulations on 3

*different patterns. Sufficient parallelization and/or judicious sampling of the parameter space could enable large-scale robust optimizations.*

^{n}## 9. Conclusion

We present a general method for designing robust metasurfaces using topology optimization. Robustness is incorporated directly into the iterative optimization process by simulating geometrically eroded and dilated devices during each iteration. We find that robust metasurfaces are relatively insensitive to edge deviations compared to devices designed without robustness criteria. We show that conventional metasurface designs can be made robust using additional topology optimization as refinement. We also find that robust metasurfaces are relatively insensitive to small scale, random geometric perturbations. An analysis of the optical modes of these devices indicate that robustness arises from the non-trivial interaction of these modes, and that such behavior relies on numerical optimization and cannot be achieved using simple design rules. Finally, we show that topology optimization can be applied to design robustness against refractive index variations as well as other parameters of interest.

## Appendix A Wavelength and incidence angle sensitivity

While all of the above devices were designed only for the targeted wavelength and incidence angle, we compute their performance at under non-ideal conditions. Figure 11(a) shows the deflectors’ performance under different wavelengths. It is interesting to note that the robust device exhibits better bandwidth than the non-robust device despite wavelength no such constraint having been placed during optimization. Figure 11(b) shows the absolute efficiency of both deflectors as a function of incidence angle. Both robust and non-robust metagratings perform well under non-normal incidence. The +2 degree angle is a limit that arises from the cutoff of the diffracted beam at larger incidence angles.

## Funding

U.S. Air Force (FA9550-18-1-0070); the Office of Naval Research (N00014-16-1-2630); the Packard Fellowship Foundation; Stanford Graduate Fellowship; National Science Foundation (NSF) through the NSF Graduate Research Fellowship.

## Acknowledgements

The authors thank S. Hoyer for a discussion of elastic deformation concepts.

This work was performed in part at the Stanford Nanofabrication Facility (SNF) and the Stanford Nano Shared Facilities (SNSF), which are supported by the National Science Foundation as part of the National Nanotechnology Coordinated Infrastructure under award ECCS-1542152.

## Disclosures

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

## References

**1. **P. Genevet, F. Capasso, F. Aieta, M. Khorasaninejad, and R. Devlin, “Recent advances in planar optics: from plasmonic to dielectric metasurfaces,” Optica **4**(1), 139–152 (2017). [CrossRef]

**2. **A. V. Kildishev, A. Boltasseva, and V. M. Shalaev, “Planar photonics with metasurfaces,” Science **339**(6125), 1232009 (2013). [CrossRef] [PubMed]

**3. **H. T. Chen, A. J. Taylor, and N. Yu, “A review of metasurfaces: physics and applications,” Rep. Prog. Phys. **79**(7), 076401 (2016). [CrossRef] [PubMed]

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

**5. **L. Wang, S. Kruk, H. Tang, T. Li, I. Kravchenko, D. N. Neshev, and Y. S. Kivshar, “Grayscale transparent metasurface holograms,” Optica **3**(12), 1504–1505 (2016). [CrossRef] [PubMed]

**6. **P. Lalanne, S. Astilean, P. Chavel, E. Cambril, and H. Launois, “Design and fabrication of blazed binary diffractive elements with sampling periods smaller than the structural cutoff,” J. Opt. Soc. Am. A **16**(5), 1143–1156 (1999). [CrossRef]

**7. **D. Sell, J. Yang, S. Doshay, K. Zhang, and J. A. Fan, “Visible light metasurfaces based on single-crystal silicon,” ACS Photonics **3**(10), 1919–1925 (2016). [CrossRef]

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

**9. **S. Colburn, A. Zhan, and A. Majumdar, “Metasurface optics for full-color computational imaging,” Sci. Adv. **4**(2), r2114 (2018). [CrossRef] [PubMed]

**10. **E. Arbabi, S. M. Kamali, A. Arbabi, and A. Faraon, “Full-Stokes Imaging Polarimetry Using Dielectric Metasurfaces,” ACS Photonics **5**(8), 3132–3140 (2018). [CrossRef]

**11. **A. Pors, M. G. Nielsen, and S. I. Bozhevolnyi, “Plasmonic metagratings for simultaneous determination of Stokes parameters,” Optica **2**(8), 716–723 (2015). [CrossRef]

**12. **A. Tittl, A. Leitis, M. Liu, F. Yesilkoy, D. Y. Choi, D. N. Neshev, Y. S. Kivshar, and H. Altug, “Imaging-based molecular barcoding with pixelated dielectric metasurfaces,” Science **360**(6393), 1105–1109 (2018). [CrossRef] [PubMed]

**13. **A. Arbabi, Y. Horie, M. Bagheri, and A. Faraon, “Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission,” Nat. Nanotechnol. **10**(11), 937–943 (2015). [CrossRef] [PubMed]

**14. **C. Pfeiffer and A. Grbic, “Cascaded metasurfaces for complete phase and polarization control,” Appl. Phys. Lett. **102**(23), 231116 (2013). [CrossRef]

**15. **J. Jensen and O. Sigmund, “Topology optimization for nanophotonics,” Laser Photonics Rev. **5**(2), 308–321 (2011). [CrossRef]

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

**17. **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] [PubMed]

**18. **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] [PubMed]

**19. **D. Sell, J. Yang, S. Doshay, and J. A. Fan, “Periodic Dielectric Metasurfaces with High-Efficiency, Multiwavelength Functionalities,” Adv. Opt. Mater. **5**(23), 1700645 (2017). [CrossRef]

**20. **D. Sell, J. Yang, E. W. Wang, T. Phan, S. Doshay, and J. A. Fan, “Ultra-High-Efficiency Anomalous Refraction with Dielectric Metasurfaces,” ACS Photonics **5**(6), 2402–2407 (2018). [CrossRef]

**21. **J. Yang and J. A. Fan, “Analysis of material selection on dielectric metasurface performance,” Opt. Express **25**(20), 23899–23909 (2017). [CrossRef] [PubMed]

**22. **Z. Lin, B. Groever, F. Capasso, A. W. Rodriguez, and M. Loncar, “Topology-Optimized Multilayered Metaoptics,” Phys. Rev. Appl. **9**(4), 044030 (2018). [CrossRef]

**23. **A. Zhan, T. K. Fryett, S. Colburn, and A. Majumdar, “Inverse design of optical elements based on arrays of dielectric spheres,” Appl. Opt. **57**(6), 1437–1446 (2018). [CrossRef] [PubMed]

**24. **R. Pestourie, C. Perez-Arancibia, Z. Lin, W. Shin, F. Capasso, and S. Johnson, “Inverse design of large-area metasurfaces,” arXiv preprint arXiv:1808.04215 (2018). [CrossRef]

**25. **V.-C. Su, C. H. Chu, G. Sun, and D. P. Tsai, “Advances in optical metasurfaces: fabrication and applications [Invited],” Opt. Express **26**(10), 13148–13182 (2018). [CrossRef] [PubMed]

**26. **D. Lin, P. Fan, E. Hasman, and M. L. Brongersma, “Dielectric gradient metasurface optical elements,” Science **345**(6194), 298–302 (2014). [CrossRef] [PubMed]

**27. **S. Doshay, D. Sell, J. Yang, R. Yang, and J. A. Fan, “High-performance axicon lenses based on high-contrast, multilayer gratings,” APL Photonics **3**(1), 011302 (2018). [CrossRef]

**28. **R. C. Devlin, M. Khorasaninejad, W. T. Chen, J. Oh, and F. Capasso, “Broadband high-efficiency dielectric metasurfaces for the visible spectrum,” Proc. Natl. Acad. Sci. U.S.A. **113**(38), 10473–10478 (2016). [CrossRef] [PubMed]

**29. **F. Wang, J. Jensen, and O. Sigmund, “Robust topology optimization of photonic crystal waveguides with tailored dispersion properties,” J. Opt. Soc. Am. B **28**(3), 387–397 (2011). [CrossRef]

**30. **A. Arbabi, Y. Horie, A. J. Ball, M. Bagheri, and A. Faraon, “Subwavelength-thick lenses with high numerical apertures and large efficiency based on high-contrast transmitarrays,” Nat. Commun. **6**(6), 7069 (2015). [CrossRef] [PubMed]

**31. **M. Schevenels, B. S. Lazaroz, and O. Sigmund, “Robust topology optimization accounting for spatially varying manufacturing errors,” Comput. Methods Appl. Math. **200**, 3613–3627 (2011).

**32. **P. Simard, D. Steinkraud, and J. C. Platt, “Best Practices for Convolutional Neural Networks Applied to Visual Document Analysis,” in Proc. of the International Conference on Document Analysis and Recognition, 2003. [CrossRef]

**33. **J. Yang, D. Sell, and J. A. Fan, “Freeform metagratings based on complex light scattering dynamics for extreme, high efficiency beam steering,” Ann. Phys. **530**(1), 1700302 (2018). [CrossRef]

**34. **J. Yang and J. A. Fan, “Topology-optimized metasurfaces: impact of initial geometric layout,” Opt. Lett. **42**(16), 3161–3164 (2017). [CrossRef] [PubMed]

**35. **J. Aarik, A. Aidla, A.-A. Kiisler, T. Uustare, and V. Sammelselg, “Effect of crystal structure on optical properties of TiO_{2} films grown by atomic layer deposition,” Thin Solid Films **305**(1–2), 270–273 (1997). [CrossRef]