High-efficiency, large-area, topology-optimized metasurfaces

ABSTRACT Metasurfaces are ultrathin optical elements that are highly promising for constructing lightweight and compact optical systems. For their practical implementation, it is imperative

to maximize the metasurface efficiency. Topology optimization provides a pathway for pushing the limits of metasurface efficiency; however, topology optimization methods have been limited to

the design of microscale devices due to the extensive computational resources that are required. We introduce a new strategy for optimizing large-area metasurfaces in a computationally

efficient manner. By stitching together individually optimized sections of the metasurface, we can reduce the computational complexity of the optimization from high-polynomial to linear. As

a proof of concept, we design and experimentally demonstrate large-area, high-numerical-aperture silicon metasurface lenses with focusing efficiencies exceeding 90%. These concepts can be

generalized to the design of multifunctional, broadband diffractive optical devices and will enable the implementation of large-area, high-performance metasurfaces in practical optical

systems. SIMILAR CONTENT BEING VIEWED BY OTHERS ADVANCES IN OPTICAL METALENSES Article 22 December 2022 OFF-AXIS METASURFACES FOR FOLDED FLAT OPTICS Article Open access 12 September 2023

INVERSE DESIGN OF METASURFACES WITH NON-LOCAL INTERACTIONS Article Open access 05 August 2020 INTRODUCTION Metasurfaces are optical devices that utilize subwavelength-scale structuring to

shape and manipulate electromagnetic waves1. They are powerful complements to bulk refractive and scalar diffractive optics and have a wide range of potential applications in imaging2,

lithography3, sensing4, and computing platforms5. To date, a broad range of metasurface design concepts have emerged. The most widely used methods, which we will term “conventional methods,”

sample the desired phase profile using discrete phase-shifter elements to form a nanoscale phased array (Fig. 1a). These methods utilize a library of simple, physically intuitive building

blocks, including anisotropic waveguides6, Mie resonators7, plasmonic resonators8, and dielectric transmit arrays9, and can quickly produce macroscale device designs. However, these

approaches lack the necessary degrees of freedom for realizing high-efficiency in devices that are designed for large-angle deflections, multiple functions, and broadband responses10,

thereby preventing metasurfaces from being practically applied in many contexts. Topology optimization is an alternative design method that can overcome the efficiency limitations of

conventional approaches11. Inverse design methods, including objective-first and adjoint-based topology optimization12, have produced photonic crystals13,14, optical demultiplexers15, and

spectral splitters16 with outstanding performance resulting from the devices’ unusual geometries and nonintuitive optical dynamics. More recently, topology optimization has led to

high-performance metasurfaces that have a broad range of capabilities, such as high-efficiency light deflection with nearly arbitrary input and output angles17, spectral sorting of plane

waves to distinct diffraction orders18, and light focusing with field-curvature correction19. While devices that are designed using topology optimization support enhanced efficiencies, they

require extensive computational resources that dramatically increase with the size of the device. Hence, current topology-optimized devices are either microscale in dimension or limited to

periodic structures with microscale unit cells. We report a conceptually new approach to metasurface design that extends the high performance of topology-optimized devices to macroscopic

areas in a computationally efficient manner. The concept is schematically illustrated in Fig. 1b and is comprised of three parts. First, we discretize the desired phase profile into a series

of wavelength-scale, linear sections. Next, we use topology optimization to design metasurface elements that optimally scatter light as required by each section of the phase profile.

Finally, we stitch all the elements together to form a full metasurface. Our approach produces devices that operate more efficiently than conventional designs by accounting for and

optimizing the near-field optical coupling between neighboring nanostructures. Our design method is also substantially more computationally efficient than existing topology optimization

methods. RESULTS ASSESSMENT OF COMPUTATIONAL EFFICIENCY To benchmark the improvements in computational efficiency that are afforded by our approach, we perform adjoint-based topology

optimization on metagratings that are made of silicon ridges. These periodic metasurfaces are designed to deflect light of a specified wavelength λ to the +1 diffraction order. We perform

electromagnetic simulations for topology optimization using rigorous coupled wave analysis (RCWA), which is also known as the Fourier modal method20,21,22,23, on a personal computer and plot

the time that is required for optimizing a full device as a function of its width _L_ in Fig. 1c. The simulation time scales approximately as \({\cal{O}}\left( {L^{2.4}} \right)\), which

corresponds to the general scaling trend for electromagnetic solvers that utilize standard matrix multiplication and inversion algorithms24. These trends demonstrate that prohibitively large

computational resources are required for directly optimizing devices that are many times larger than the wavelength. If we instead divide the metagrating into sections of width _d_, we have

_L_/_d_ sections to optimize. The total computation time now scales as \({\cal{O}}\left( {d^{2.4} \cdot L/d} \right) = {\cal{O}}\left( {d^{1.4} \cdot L} \right)\), which is a linear

function of the total size. This expression also shows that the computation time decreases with decreasing section size. In practice, we found that using 3λ-wide sections minimizes the total

computation time, as the benefits from using smaller sections are outweighed by the computational overhead in the electromagnetic solver. We plot the computation time for a metasurface that

has been optimized using 3λ-wide sections in Fig. 1c. The observed reduction in the computation time enables us to optimize devices with dimensions that are many times larger than the

wavelength using realistic computational resources. For example, via our method, we can produce millimeter-scale topology-optimized metasurfaces that operate at near-infrared wavelengths in

less than 1 day using a personal computer. In contrast, optimizing the full device at once would take nearly 1 year and require intractably large amounts of memory. We realize additional

savings in computation time via our method when using multiple computing cores because individual segments can be optimized on different cores. Our design concept allows the segments to be

treated independently, thereby enabling the optimizations to be parallelized without issues concerning race conditions and synchronization25. With _N_ computing nodes, metasurfaces that are

subdivided into _N_ sections can be optimized in the same amount of time as it takes to optimize a single segment. These reductions in computational complexity also apply to fully

three-dimensional topology-optimized metasurfaces, which exhibit even more severe scaling trends: the time it takes to optimize a metasurface of size _L_ × _L_ all at once scales

approximately as \({\cal{O}}\left( {L^{4.8}} \right)\). This trend indicates the necessity of linear sectioning for these more intricate design problems, which would reduce the runtime to a

more reasonable \({\cal{O}}\left( {L^2} \right)\). SECTIONING OF A CURVILINEAR PHASE PROFILE The process of approximating a curvilinear phase profile with a series of linear sections

introduces wavefront error, which is denoted as ϵ (Fig. 1b, inset). However, this error has a negligible impact on the overall metasurface performance if the sections are sufficiently small.

To analyze the effect of linearizing a general curvilinear phase profile _ϕ_(_x_), we locally describe each section of _ϕ_(_x_) at location _x_0 using a 2nd-order Taylor series expansion as

follows: $$\phi \left( x \right) \approx \phi \left( {x_0} \right) + \phi \prime \left( {x_0} \right)\left( {x - x_0} \right) + \frac{1}{2}\phi \prime \prime \left( {x_0} \right)\left( {x -

x_0} \right)^2$$ (1) This section can be approximated as a line of slope \(\phi \prime \left( {x_0} \right)\) with a phase offset of _ϕ_(_x_0) + Δ_ϕ_, which incurs an error of

\({\it{\epsilon }}\left( x \right) = \phi \prime \prime \left( {x_0} \right)\left( {x - x_0} \right)^2/2 - {\mathrm{\Delta }}\phi\). Given a section of length \(d\) at position _x_0, the

root-mean-square (RMS) wavefront error is minimized when \({\mathrm{\Delta }}\phi = \phi \prime \prime \left( {x_0} \right)d^2/24\) and is expressed as follows: $${\it{\epsilon

}}_{{\mathrm{rms}}}\left( {{\mathrm{\Delta }}x} \right) = \frac{1}{{12\sqrt 5 }}\phi \prime \prime \left( {x_0} \right)d^2$$ (2) We use this result to analyze the impact of linearization on

a focusing cylindrical lens, which enables us to benchmark the device performance using well-established metrics in lens design. To quantify the performance of lenses that are constructed

with linear sections, we use the Strehl ratio26, which is a metric that compares the diffraction efficiency of our lenses to that of an ideal lens. For a lens that focuses normally incident

light, the ideal phase profile is \(\phi \left( x \right) = \left( {2\pi /\lambda } \right)( {f - \sqrt {f^2 + x^2} } )\)27, where _f_ is the focal length and _λ_ is the wavelength. A lens

with a Strehl ratio of 0.98, which corresponds to an RMS wavefront error of λ/50, can be realized if we use linear sections that are no larger than the following: $$d\, < \,0.73\sqrt

{f\lambda }$$ (3) This equation provides a practical and quantitative guide for linearizing phase profiles in a manner that minimizes the phase error. As an example, consider a cylindrical

lens with a focal length of _f_ = 36_λ_ and a numerical aperture (NA) of 0.7 (Fig. 2a). Using Eq. 3, we expect that linearizing the phase profile with segments that are smaller than 4.4λ

will have a negligible impact on the performance. To verify, we simulate lenses that are linearized with various segment lengths and calculate the field intensities at the focal plane. These

field intensity profiles are plotted in Fig. 2b, according to which the lenses that are linearized with section lengths of 2λ and 4λ are nearly indistinguishable from the ideal lens. Our

approach to sectioning readily extends to three-dimensional phase profiles, which can be approximated as series of planar tiles. This approximation is discussed in detail in the

Supplementary Section. For a hyperboloid that is discretized into tiles with dimensions _d_ × _d_, an RMS wavefront error of λ/50 can be realized if: $$d\, < \,0.61\sqrt {f\lambda }$$ (4)

TOPOLOGY OPTIMIZATION FOR FINITE-SIZED, ISOLATED DEVICES To design metasurface elements that have the desired linear phase profiles, we utilize adjoint-based topology optimization28.

Adjoint-based optimization is an iterative algorithm that modifies the device’s dielectric constant distribution, namely, _ε_(_x_), to maximize a figure of merit (FoM). Our objective is to

optimize a device that scatters normally incident electromagnetic waves in a desired direction with electric field amplitude _E_tgt and phase _ϕ_tgt. To compute the FoM, we run a forward

simulation in which waves that are incident onto the metasurface element scatter in the desired direction with field amplitude _E_fwd and phase _ϕ_fwd. Near-to-far-field transformations from

the forward simulations are used to evaluate _E_fwd and _ϕ_fwd29. The FoM describes the difference between the current and desired responses and has the following form: $$FoM = - A_1\left[

{\left| {\boldsymbol{E}}_{\text{tgt}} \right|^2 - \left| {\boldsymbol{E}}_{\text{fwd}} \right|^2} \right]^2 - A_2\left[ {\arg \left( {e^{i\left( {\phi _{\text{tgt}} - \phi _{\text{fwd}}}

\right)}} \right)} \right]^2$$ (5) The terms _A_1 and _A_2 are weights that balance how strongly the FoM is biased toward optimizing the amplitude and phase, respectively. To determine how

_ε_(_x_) should be modified to improve the FoM each iteration, we perform a pair of forward and adjoint simulations and record the electric fields in the device for each excitation

condition. These fields are used to calculate _δFoM_, which is the gradient of the FoM with respect to the dielectric constant at each position _x_: $$\begin{array}{l}\delta FoM =

{2}A_{1}\left( {\left| {\boldsymbol{E}}_{\text{tgt}} \right|^{2} - \left| {\boldsymbol{E}}_{\text{fwd}} \right|^2} \right){\cal{R}}e\left\{ {{\boldsymbol{E}}_{\text{fwd}} \cdot \delta

{\boldsymbol{E}}^\ast } \right\}\\ - {2}A_{2}\left( \phi _{\text{tgt}} - \phi _{\text{fwd}} \right)\frac{1}{{\left| {\boldsymbol{E}}_{\text{fwd}} \right|^2}}{\cal{I}}m\left\{

{\boldsymbol{E}}_{\text{fwd}} \cdot \delta {\boldsymbol{E}^ \ast } \right\}\end{array}$$ (6) where \(\delta {\boldsymbol{E}}\) is a function of the adjoint field and represents the

variations of the field in the target direction in response to variations of the refractive index within the device. A more detailed discussion of the adjoint optimization method that is

applied here is provided in the Supplementary Section. APERIODIC FOURIER MODAL METHOD To apply these concepts to the design of isolated, finite-sized device elements, we have developed an

aperiodic Fourier modal method (AFMM), which is a hybrid method that combines a solver for periodic systems with perfectly matched layers (PMLs). The key challenge of implementing PMLs

involves describing both the periodic incident plane wave (the input field) and the aperiodic scattered field (the output field) of the isolated device within the same formalism. To address

this challenge, we introduce a hybrid method that combines a Fourier basis, Maxwell’s equations in complex coordinates, and the Stratton-Chu integral formalism30,31,32,33,34. A metasurface

is typically composed of a single layer of patterned material. The patterned material can be expressed as a distribution of the relative permeability, namely, _ϵ_(_x, y_), and the

permeability, namely, _μ_(_x, y_), on the _xy_-plane. Along the thickness of the device in the _z_-direction, the device cross-section is constant. In this case, it can be shown from

Maxwell’s equations that the transverse electric fields satisfy the following eigenvalue equation: $$- \gamma ^2\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{E}}_x} \\ {{\boldsymbol{E}}_y}

\end{array}} \right] = {\cal{L}}_{EH}{\cal{L}}_{HE}\left[ {\begin{array}{*{20}{c}} {{\boldsymbol{E}}_x} \\ {{\boldsymbol{E}}_y} \end{array}} \right]$$ (7) where \({\cal{L}}_{EH}\) and

\({\cal{L}}_{HE}\) are differential operators that are based on _ε_ and _μ_ and the electric field solution can be expressed as \({\mathrm{\Phi }}\left( {x,y} \right)e^{ - ik\gamma z}\). The

Fourier modal method can easily be used to expand the input fields, which are assumed to be periodic, into a Fourier basis: $${\mathrm{\Phi }}_{in}\left( {x,y,z} \right) = \mathop {\sum

}\limits_p A_pe^{ - ik\gamma _pz}\mathop {\sum }\limits_{nm} \delta _{nm,p}e^{ - ik\alpha _nx}e^{ - ik\beta _my}$$ (8) Next, the PML boundary conditions are described similarly. The PMLs can

be introduced via a change of coordinates \(\left( {x,y,z} \right) \to (\tilde x,\tilde y,\tilde z)\)30,31: $$\left\{ {\begin{array}{*{20}{c}} {\tilde x\left( x \right) = \left( {\chi _x -

i\eta _x} \right)\left( {x - x_0} \right) + x_0} \\ {\tilde y\left( y \right) = \left( {\chi _y - i\eta _y} \right)\left( {y - y_0} \right) + y_0} \end{array}} \right\}$$ (9) Here, the

parameter _χ_ controls the scaling of the PML layers, while the parameter _η_ controls the PML attenuation. This transformation is also useful for expressing the output scattered fields,

which are computed as eigenfunctions of Maxwell’s equations in complex coordinates: $${\mathrm{\Phi }}_s\left( {x,y,z} \right) = \mathop {\sum }\limits_p B_pe^{ - ik\gamma _pz}\mathop {\sum

}\limits_{nm} {\mathrm{\Phi }}_{nm,p}e^{ - ik\alpha _nx}e^{ - ik\beta _my}$$ (10) Finally, the Stratton-Chu integral equation computes the radiated field in all of space34. In this study, we

design devices that are comprised of nanoridges, which can be described in this analysis by a one-dimensional Fourier basis. SECTION OPTIMIZATION DEMONSTRATION To demonstrate adjoint-based

optimization in the design of metasurface elements, we design a 2.5λ-wide element that scatters incident TM-polarized light at a 20° angle with a phase response of π/2. The dielectric

distribution and the amplitude and phase scattering profiles of the element at successive iterations in the optimization process are shown in Fig. 3a. The optimization begins with a random

dielectric continuum with values between air and silicon. After a few iterations, the continuum begins to strongly scatter light at the desired angle. The final metasurface section is a

binary structure of silicon in air that possesses a peak scattering amplitude and phase response matching with the targeted values. The full width at half maximum of the scattering peak is

consistent with that expected from light that is diffracting from a 2.5λ-wide aperture; hence, the element is performing directional scattering near its physical limits. An intensity plot of

the scattered fields from the fully optimized section (Fig. 3b) shows strong near-field coupling between neighboring nanostructures; hence, optimal near-field coupling is responsible for

mediating strong scattering in the desired direction. As a method of gradient descent, adjoint-based topology optimization is a local optimizer and is sensitive to the initial dielectric

distribution10. To obtain high-performance elements for a desired scattering angle and phase target, we perform ten optimizations with various initial dielectric distributions and select the

best result. OPTIMIZED METALENS SIMULATION RESULTS Combining everything into a proof-of-concept demonstration, we stitch together optimized metasurface elements to construct cylindrical

metalenses. We design these metalenses to focus TM-polarized light at a wavelength of 640 nm. To enable device operation at visible wavelengths, we use 250-nm-thick crystalline silicon,

which has relatively low absorption compared to polycrystalline and amorphous silicon, but much higher index contrast than materials such as titanium dioxide35,36. First, we design and

simulate 64-μm-wide metalenses with NAs that range from 0.2 to 0.9. We divide the metalenses into sections that are 2 μm wide, which is below the phase error limit of Eq. 3 and near the

optimal size for efficient computation (Fig. 1c). Further reductions in section size lead to degradation of the device performance. The reason can be traced to our design of each metasurface

element, which is optimized in isolation with PML boundary conditions. When the elements are stitched together to produce a device, the optical fields that are guided by a single element

have evanescent tails and can couple to a neighboring element in a parasitic manner. Smaller section sizes require more elements to be stitched together to produce a desired metasurface,

thereby resulting in more boundaries and more parasitic coupling. Below a section size of 2λ, the device performance begins to degrade and below a section size of 1λ, the aperiodic boundary

conditions in our optimizer are no longer valid. There are a few approaches for addressing the issue of stitching error. One is to keep the section size relatively large compared to the

wavelength. Another is to perform boundary optimization on the stitched regions to eliminate the stitching error. A third approach, which we use here, is to separate silicon structures from

other sections by a gap of at least 0.2λ, thereby reducing the near-field coupling between sections. To ensure a reduction in the stitching error with this scheme, we simulate stitched

sections to check for spurious diffraction and redesign the sections in the event of excess error. The efficiencies of our simulated metalenses are summarized in Fig. 4a. The absolute

efficiency is defined as the amount of power that is contained in the principal lobe of the focus compared to that of an ideal lens with 100% transmission. The relative efficiency, or

focusing efficiency, compares the power in the principal lobe to that of an ideal lens that transmits the same amount of power as the device37. This efficiency corresponds to the efficiency

of the diffraction process, as it removes the effects of absorption from the material and reflection at the metalens interface. The efficiency plots show that the relative efficiencies are

consistently high, namely, above 93%, with minimal drop-off in performance as the NA increases. This trend is unlike that of conventional metalenses, where the efficiency decreases with

increasing NA because conventional designs cannot efficiently deflect light at large angles27,38. The absolute efficiencies of the metalenses all exceed 75%, with approximately 10% of the

light reflected from the metalens and 10% absorbed by the silicon. Reflection losses can be reduced via the use of more intricate three-dimensional silicon nanostructures, while absorption

losses can be minimized by designing silicon-based devices for longer wavelength operation28. A simulated field profile of a metalens with an NA of 0.9 is shown in Fig. 4b, which

demonstrates that the lens focuses strongly with minimal spurious diffraction. PERFORMANCE OF FABRICATED METALENSES We design, fabricate, and characterize 200-μm-wide metalenses with NAs of

0.2, 0.5, and 0.8. To prepare crystalline silicon thin films on glass, we use hydrogen silsesquioxane to bond silicon-on-insulator wafers onto Pyrex wafers under high temperature and

pressure36,39. After removing the silicon handle wafer and the buried oxide layer, we pattern and etch the devices via standard electron beam lithography and dry etching techniques. We

characterize the metalenses by collimating polarized, monochromatic light from a tunable white-light laser onto the devices and imaging the light at the focal plane with a ×100 objective (NA

 = 0.9) and a CCD sensor. Scanning electron microscope images of the center of a representative device are shown in Fig. 5a and show silicon nanostructures that exhibit smooth and vertical

sidewalls. The metalenses all have relative efficiencies that exceed 89% and absolute efficiencies that exceed 67% for all NAs, which are within 10% of the simulated values (Fig. 5b). All

the metalenses exhibit diffraction-limited performance, as shown by the theoretical and experimental intensity plots in Fig. 5c–e. The device with an NA of 0.8 can focus light at a

wavelength of 640 nm to a spot with a beam waist of 340 nm. The central lobes of the foci are all much stronger than the side lobes; hence, the focusing efficiency is high. The metalenses

maintain reasonably high efficiencies for wavelengths that range from 580 to 700 nm, as shown in Fig. 5f–h. As these lenses are not designed to be achromatic, the focal length shifts with

the wavelength. The simulated shifts are shown in the Supplementary section in Fig. S3. Future work will focus on generalizing our design approach to include achromatic functionality, which

can be addressed by modifying the optimizer’s figure of merit to include multiple wavelengths18. The figure of merit can be specified so that each metasurface section deflects all

wavelengths in the same direction and realizes the correct dispersion for ensuring constructive interference at the focus40,41. DISCUSSION In summary, we present a computationally efficient

method for designing large-area, topology-optimized metasurfaces. Our approach of approximating a desired phase profile with a series of linear segments renders the design problem

computationally tractable while minimally impacting the device performance. Wavelength-scale scattering elements that are designed via topology optimization can strongly scatter light at a

desired angle and phase, thereby serving as high-performance building blocks for larger metasurfaces. The metalenses that are demonstrated here have limited functionality and other heuristic

design methods are available that can produce high-performance, high-numerical-aperture metalenses42. What differentiates our method is that our design platform can be generalized

systematically to multifunctional, broadband, and even multilayer metasurfaces with large areas and high efficiencies. To practically extend topology optimization in these directions,

continued improvements in electromagnetic solver speed will be required, particularly for fully three-dimensional devices. For multilayer devices, new innovations will be needed to address

the stitching error, which can compound due to interactions between neighboring elements. We envision that these design methods will enable the realization of compact and monolithic

electromagnetic devices that exhibit high efficiencies and practical functionality across the electromagnetic spectrum. MATERIALS AND METHODS SAMPLE FABRICATION Crystalline silicon on glass

wafers was prepared via an approach that is similar to that done by Sell, et al.36. A Pyrex and silicon-on-insulator (SOI) wafer were bonded under high temperature and pressure using

hydrogen silsesquioxane. The silicon handle wafer was etched away using SF6 plasma and the buried oxide layer was etched away using HF solution. The wafer was diced into pieces for later

use. The metasurfaces were patterned onto a piece in the AR-P-6200 resist via electron beam lithography. Aluminum oxide was evaporated onto a piece and removed using liftoff to form a hard

mask. The silicon was etched using a Cl2 and HBr plasma. The aluminum oxide was removed using a solution of HCl at 50 °C. Additional details can be found in the supplementary information.

MEASUREMENT PROCEDURE A collimated beam of light from a tunable white light laser is filtered with a longpass filter and polarized with a linear polarizer. Then, the beam passes through a

weak cylindrical lens to focus it onto the metasurface aperture without significantly changing the angle of incidence. The metasurface focuses the light. Next, the focal plane is imaged onto

a CMOS detector using a 0.9 NA, 100x objective and a tube lens. To evaluate the lens efficiencies, a line scan of the focal plane image is obtained. The lens relative efficiency is

calculated by integrating over a small aperture around the main focal lobe and comparing it to what would be expected for an ideal lens of the same numerical aperture. The lens absolute

efficiency is calculated by multiplying the relative efficiency by the measured transmission. A schematic diagram of the optical setup and additional characterization details can be found in

the supplementary information. SIMULATIONS We used Reticolo, a rigorous coupled-wave analysis solver, as a basis to perform optimizations21. The final metasurface lens designs were

simulated using Lumerical FDTD. The efficiencies were calculated by determining the optical power that is contained in the main lobe of the focus and comparing it to that expected from an

ideal lens. REFERENCES * Yu, N. F. & Capasso, F. Flat optics with designer metasurfaces. _Nat. Mater._ 13, 139–150 (2014). Article  ADS  Google Scholar  * Pahlevaninezhad, H. et al.

Nano-optic endoscope for high-resolution optical coherence tomography in vivo. _Nat. Photonics_ 12, 540–547 (2018). Article  ADS  Google Scholar  * Chen, K. et al. Moiré nanosphere

lithography. _ACS Nano_ 9, 6031–6040 (2015). Article  Google Scholar  * Tittl, A. et al. Imaging-based molecular barcoding with pixelated dielectric metasurfaces. _Science_ 360, 1105–1109

(2018). Article  ADS  MathSciNet  Google Scholar  * Silva, A. et al. Performing mathematical operations with metamaterials. _Science_ 343, 160–163 (2014). Article  ADS  MathSciNet  Google

Scholar  * Mueller, J. P. B. et al. Metasurface polarization optics: independent phase control of arbitrary orthogonal states of polarization. _Phys. Rev. Lett._ 118, 113901 (2017). Article

  ADS  Google Scholar  * Iyer, P. P., Pendharkar, M. & Schuller, J. A. Electrically reconfigurable metasurfaces using heterojunction resonators. _Adv. Opt. Mater._ 4, 1582–1588 (2016).

Article  Google Scholar  * Ding, F. et al. Broadband high-efficiency half-wave plate: a supercell-based plasmonic metasurface approach. _ACS Nano_ 9, 4111–4119 (2015). Article  Google

Scholar  * Arbabi, A. et al. Dielectric metasurfaces for complete control of phase and polarization with subwavelength spatial resolution and high transmission. _Nat. Nanotechnol._ 10,

937–943 (2015). Article  ADS  Google Scholar  * Yang, J. & Fan, J. A. Topology-optimized metasurfaces: impact of initial geometric layout. _Opt. Lett._ 42, 3161–3164 (2017). Article  ADS

  Google Scholar  * Lalau-Keraly, C. M. et al. Adjoint shape optimization applied to electromagnetic design. _Opt. Express_ 21, 21693–21701 (2013). Article  ADS  Google Scholar  * Molesky,

S. et al. Inverse design in nanophotonics. _Nat. Photonics_ 12, 659–670 (2018). Article  ADS  Google Scholar  * Frandsen, L. H. et al. Broadband photonic crystal waveguide 60° bend obtained

utilizing topology optimization. _Opt. Express_ 12, 5916–5921 (2004). Article  ADS  Google Scholar  * Borel, P. I. et al. Topology optimization and fabrication of photonic crystal

structures. _Opt. Express_ 12, 1996–2001 (2004). Article  ADS  Google Scholar  * Piggot, A. Y. et al. Inverse design and demonstration of a compact and broadband on-chip wavelength

demultiplexer. _Nat. Photonics_ 9, 374–377 (2015). Article  ADS  Google Scholar  * Xiao, T. P. et al. Diffractive spectral-splitting optical element designed by adjoint-based electromagnetic

optimization and fabricated by femtosecond 3D direct laser writing. _ACS Photonics_ 3, 886–894 (2016). Article  Google Scholar  * Sell, D. et al. Ultra-high-efficiency anomalous refraction

with dielectric metasurfaces. _ACS Photonics_ 5, 2402–2407 (2018). Article  Google Scholar  * Sell, D. et al. Periodic dielectric metasurfaces with high-efficiency, multiwavelength

functionalities. _Adv. Opt. Mater._ 5, 1700645 (2017). Article  Google Scholar  * Lin, Z. et al. Topology-optimized multilayered metaoptics. _Phys. Rev. Appl._ 9, 044030 (2018). Article  ADS

  Google Scholar  * Lalanne, P. & Morris, G. M. Highly improved convergence of the coupled-wave method for TM polarization. _J. Opt. Soc. Am. A_ 13, 779–784 (1996). Article  ADS  Google

Scholar  * Hugonin, J. P., & Lalanne, P. Reticolo Software for Grating Analysis. (Institut d’Optique, 2005). * Knop, K. Rigorous diffraction theory for transmission phase gratings with

deep rectangular grooves. _J. Opt. Soc. Am._ 68, 1206–1210 (1978). Article  ADS  Google Scholar  * Li, L. F. Use of Fourier series in the analysis of discontinuous periodic structures. _J.

Opt. Soc. Am. A_ 13, 1870–1876 (1996). Article  ADS  Google Scholar  * Coppersmith, D. & Winograd, S. Matrix multiplication via arithmetic progressions. _J. Symb. Comput._ 9, 251–280

(1990). Article  MathSciNet  Google Scholar  * Patterson, D. A. & Hennessy, J. L. Computer Organization and Design. 5th edn. (Morgan Kaufmann, Waltham: MK, 2014). * Strehl, K.

Aplanatische und fehlerhafte abbildung im fernrohr. _Z. für Instrum._ 15, 362–370 (1895). Google Scholar  * Arbabi, A. et al. Subwavelength-thick lenses with high numerical apertures and

large efficiency based on high-contrast transmitarrays. _Nat. Commun._ 6, 7069 (2015). Article  Google Scholar  * Sell, D. et al. Large-angle, multifunctional metagratings based on freeform

multimode geometries. _Nano Lett._ 17, 3752–3757 (2017). Article  ADS  Google Scholar  * Yang, J. J., Hugonin, J. P. & Lalanne, P. Near-to-far field transformations for radiative and

guided waves. _ACS Photonics_ 3, 395–402 (2016). Article  Google Scholar  * Chew, W. C. & Weedon, W. H. A 3D perfectly matched medium from modified Maxwell’s equations with stretched

coordinates. _Microw. Opt. Technol. Lett._ 7, 599–604 (1994). Article  ADS  Google Scholar  * Chew, W. C., Jin, J. M. & Michielssen, E. Complex coordinate stretching as a generalized

absorbing boundary condition. _Microw. Opt. Technol. Lett._ 15, 363–369 (1997). Article  Google Scholar  * Teixeira, F. L. & Chew, W. C. Differential forms, metrics, and the

reflectionless absorption of electromagnetic waves. _J. Electromagn. Waves Appl._ 13, 665–686 (1999). Article  MathSciNet  Google Scholar  * Edee, K., Granet, G. & Plumey, J. P. Complex

coordinate implementation in the curvilinear coordinate method: application to plane-wave diffraction by nonperiodic rough surfaces. _J. Opt. Soc. Am. A_ 24, 1097–1102 (2007). Article  ADS 

MathSciNet  Google Scholar  * Stratton, J. A. & Chu, L. J. Diffraction theory of electromagnetic waves. _Phys. Rev._ 56, 99–107 (1939). Article  ADS  Google Scholar  * Yang, J. J. &

Fan, J. A. Analysis of material selection on dielectric metasurface performance. _Opt. Express_ 25, 23899–23909 (2017). Article  ADS  Google Scholar  * Sell, D. et al. Visible light

metasurfaces based on single-crystal silicon. _ACS Photonics_ 3, 1919–1925 (2016). Article  Google Scholar  * Palmer, C. & Loewen, E. Diffraction Grating Handbook. 6th edn. (New York:

Newport Corporation, Rochester, 2005). * Yang, J. J., Sell, D. & Fan, J. A. Freeform metagratings based on complex light scattering dynamics for extreme, high efficiency beam steering.

_Ann. der Phys._ 530, 1700302 (2018). Article  ADS  MathSciNet  Google Scholar  * Chung, J. W. et al. Seamless on-wafer integration of Si(100) MOSFETs and GaN HEMTs. _IEEE Electron Device

Lett._ 30, 1015–1017 (2009). Article  ADS  Google Scholar  * Shrestha, S. et al. Broadband achromatic dielectric metalenses. _Light.: Sci. Appl._ 7, 85 (2018). Article  ADS  Google Scholar 

* Chen, W. T. et al. A broadband achromatic metalens for focusing and imaging in the visible. _Nat. Nanotechnol._ 13, 220–226 (2018). Article  ADS  Google Scholar  * Paniagua-Dominguez, R.

et al. A metalens with a near-unity numerical aperture. _Nano Lett._ 18, 2124–2132 (2018). Article  ADS  Google Scholar  Download references ACKNOWLEDGEMENTS This work was supported by the

U.S. Air Force under Award Number FA9550-18-1-0070, the Office of Naval Research under Award Number N00014-16-1-2630, and the David and Lucile Packard Foundation. DS was supported by the

National Science Foundation (NSF) through an NSF Graduate Research Fellowship. SD was supported by the Department of Defense (DoD) through an NDSEG Research Fellowship. 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. AUTHOR INFORMATION AUTHORS AND AFFILIATIONS * Department of Electrical Engineering, Stanford University,

Stanford, CA, 94305, USA Thaibao Phan, Evan W. Wang, Jianji Yang & Jonathan A. Fan * Department of Applied Physics, Stanford University, Stanford, CA, 94305, USA David Sell & Sage

Doshay * Université Clermont Auvergne, Institut Pascal, BP 10448, F-63000, Clermont-Ferrand, France Kofi Edee * CNRS, UMR 6602, Institut Pascal, F-63177, Aubière, France Kofi Edee Authors *

Thaibao Phan View author publications You can also search for this author inPubMed Google Scholar * David Sell View author publications You can also search for this author inPubMed Google

Scholar * Evan W. Wang View author publications You can also search for this author inPubMed Google Scholar * Sage Doshay View author publications You can also search for this author

inPubMed Google Scholar * Kofi Edee View author publications You can also search for this author inPubMed Google Scholar * Jianji Yang View author publications You can also search for this

author inPubMed Google Scholar * Jonathan A. Fan View author publications You can also search for this author inPubMed Google Scholar CONTRIBUTIONS T.P., D.S., K.E. and J.Y. formulated the

design algorithms. T.P., D.S., E.W. and S.D. fabricated the devices. D.S. performed the measurements and analyzed the data. T.P. and J.F. wrote the manuscript with input from all authors.

CORRESPONDING AUTHOR Correspondence to Jonathan A. Fan. ETHICS DECLARATIONS CONFLICT OF INTEREST The authors declare that they have no conflict of interest. SUPPLEMENTARY INFORMATION

SUPPLEMENTARY INFORMATION. RIGHTS AND PERMISSIONS OPEN ACCESS This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation,

distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and

indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to

the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will

need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/. Reprints and permissions ABOUT THIS ARTICLE

CITE THIS ARTICLE Phan, T., Sell, D., Wang, E.W. _et al._ High-efficiency, large-area, topology-optimized metasurfaces. _Light Sci Appl_ 8, 48 (2019).

https://doi.org/10.1038/s41377-019-0159-5 Download citation * Received: 09 December 2018 * Revised: 30 April 2019 * Accepted: 05 May 2019 * Published: 29 May 2019 * DOI:

https://doi.org/10.1038/s41377-019-0159-5 SHARE THIS ARTICLE Anyone you share the following link with will be able to read this content: Get shareable link Sorry, a shareable link is not

currently available for this article. Copy to clipboard Provided by the Springer Nature SharedIt content-sharing initiative