INTRODUCTION
============

Designing optical elements based on arrays of discrete dielectric scatterers is currently a vibrant research topic ([@R1]). These discrete scatterers can control the local amplitude and phase of incident optical fields ([@R2]). Research on this class of optical elements resulted in the ultimate miniaturization of traditional optical elements such as lenses ([@R3]--[@R6]), free-form optics ([@R7], [@R8]), and retroreflectors ([@R9]). In addition, novel multifunctional optical elements based on polarization or wavelength multiplexing ([@R10], [@R11]) and engineered point spread functions (PSFs) ([@R12]) have been demonstrated by arrays of dielectric scatterers. The periodicity of these arrays has been pushed to wavelength, and even subwavelength, scales, giving rise to the field of dielectric metasurfaces ([@R13]). To date, this body of research has been primarily conducted using intuition-driven forward design methods. These forward design methods implement a desired phase profile using precompiled libraries of discrete scatterers. Properties of these scatterers are computed under periodic boundary conditions, and a single scatterer is assumed to behave like a periodic array of the scatterers. This approximation, generally known as the local phase approximation, neglects the differences in interactions between adjacent scatterers in the designed optical element and may not be valid for a phase profile with steep gradient ([@R14]). Moreover, in the absence of a known phase profile, such an approach cannot be used to create an optical element.

Inverse design methods, in contrast to forward design, begin by specifying a figure of merit (FOM) as a function of adjustable scatterer parameters that attempts to accurately encapsulate the performance of the optical element. Rather than picking the scatterer configuration by trial and error, the algorithm calculates both the FOM and the gradient of the FOM for a given configuration of the dielectric scatterers. Then, the algorithm iteratively moves toward a configuration using the gradient as an update direction to optimize the FOM. In this manner, inverse design methods can explore unintuitive scatterer configurations that would otherwise be ignored by intuitive forward design methods. In particular, inverse design methods excel for arbitrary designer optical functions that do not correspond to a known, analytical phase profile. These gradient-based inverse design methods were initially explored with great success in the context of plasmonic nanostructures ([@R15]) and dielectric planar integrated photonic elements ([@R16]--[@R19]). More recently, they have been used to design discrete scatterer--based free-space two-dimensional optical elements. The resulting devices have shown increased beam deflection efficiency at large angle ([@R20]) and multilayered lenses ([@R21], [@R22]). These high-performance devices with unorthodox forms clearly show the use of gradient-based optimization algorithms for optical element design. All these methods solve Maxwell's equations in a meshed design space; the dielectric permittivity at each mesh point is allowed to vary during each iteration. While such an approach is broadly applicable to any electromagnetics problem, it is memory intensive and does not scale well to large systems with small feature sizes. This is why all the inverse design methods have been limited to either a two-dimensional element design ([@R21]--[@R23]) or a small three-dimensional unit cell ([@R20]) designed for a periodic array. Recently, the large-scale inverse design of a two-dimensional metasurface was reported ([@R23]), although the design relies on the local phase approximation. In addition, neural network--based inverse design methods have shown great promise and impressive speedups, but thus far have only been used for tailoring individual scatterers ([@R24], [@R25]) and not for large-scale design.

Here, we use a gradient-based inverse design method for arrays of spherical Mie scatterers to engineer a PSF in three dimensions ([@R5]). Specifically, we aim to create an optical element that produces discrete focal points on a helical path along the optical axis. This type of optical element with a rotating PSF can be used for depth sensing, as the PSF rotates deterministically about the optical axis for different values of defocus ([@R26]). Our method solves the electromagnetic scattering problem using the generalized multisphere Mie theory (GMMT), which has been very successful in simulating the properties of large arrays of spheres ([@R27], [@R28]). The GMMT provides an analytical description of the scattering problem using smooth and continuous functions, making it ideal for a gradient-based optimization method. Instead of using the dielectric permittivity distribution as the variable, the method optimizes over the radius distribution of the spheres in an array. While this restricts the overall design space compared with methods solving Maxwell's equations directly, the GMMT is much less memory intensive for simulating large systems.

RESULTS
=======

The description and application of the theory to our specific problem are provided in section S1. We fabricated and experimentally tested the designed optical elements for the wavelengths of 1.55 and 3 μm. The designed elements cover areas of approximately 144 by 144 μm^2^ and 200 by 200 μm^2^, the largest area to date, designed using inverse electromagnetic design.

For devices at both wavelengths, the designed intensity profile is chosen to be a discrete helix specified by eight individual focal spots along the optical axis. Each focal spot is located in a distinct focal plane, separated along the optical axis by 28 μm (57 μm) for the 1.55 μm (3 μm) wavelength designs. The focal spots are arranged on circles of radius 12 μm (20 μm) for the 1.55 μm (3 μm) wavelength. A detailed description of the FOM and optimization method is provided in the Supplementary Materials. We compute the gradient using the adjoint method ([@R17], [@R19], [@R21], [@R29], [@R30]), and our formulation of the gradient can be found in the Supplementary Materials.

We chose to use the Nanoscribe GT two-photon lithography system to fabricate our devices to properly resolve the spherical scatterers. We fabricate the structures using the 63× objective paired with the IP-Dip resist (*n* \~ 1.47 at our design wavelengths), giving us the highest possible fabrication resolution (\~200 nm). For both devices, our initial condition is a set of identical spheres on a square grid in vacuum, as we found that the addition of a substrate affects the device performance only slightly ([@R30]). The specific optimization parameters are presented in the Supplementary Materials. For the 1.55 and 3 μm devices, we chose periodicities of 2.42 μm (superwavelength) and 2.9 μm (subwavelength), respectively. Schematic and scanning electron micrographs (SEMs) of the gold-coated 1.55 μm device are shown in [Fig. 1](#F1){ref-type="fig"}. As shown in the SEMs ([Fig. 1, C and D](#F1){ref-type="fig"}), the scatterers are not completely spherical and have rough edges. These fabrication imperfections and the substrate are not accounted for in our simulation and optimization process.

![SEMs of the 1.55 μm device.\
(**A**) Schematic of the sphere layout. (**B**) SEM of the 1.55 μm device coated in gold showing the entire device from a top-down view, respectively. Zoomed-in SEMs showing the fabrication imperfections from an angled (**C**) and top-down (**D**) view of the spheres.](aax4769-F1){#F1}

To test our devices, we used two separate microscopes for each of the two wavelengths. The 1.55 μm device's experimental and simulated performance is shown in [Fig. 2](#F2){ref-type="fig"}. In the experiment, we find that the optical element produces a clear high-intensity spot that matches simulation predictions. However, the experiment shows a lower contrast between the focal spot and the background compared with the simulation. In addition, in [Fig. 2](#F2){ref-type="fig"} (E, F, and N), we can identify smaller hotspots of substantial intensity that are not present in simulation. The actual focal spots are produced at the correct spatial location in plane, but there is a slight offset between the expected locations of the focal spot along the optical axis for the focal spots shown in [Fig. 2](#F2){ref-type="fig"} (M and N). We attribute these discrepancies between the simulation and experiment primarily to fabrication errors, which have been shown to be an important factor in performance. Results on a previous fabrication run with comparatively poor results are shown in the Supplementary Materials. Both the experimental and simulation results show a gradual increase in the focal spot size with increasing focal plane distance.

![Simulated and experimental device performance.\
(**A** to **D** and **I** to **L**) Simulated intensity profiles produced at specific distances from the device surface showing the focal spot rotating in the *x*-*y* plane. (**E** to **H** and **M** to **P**) Experimentally acquired intensity profiles at specific distances from the device surface. Experimental images share the same linear intensity scale, and simulated images share the same linear intensity scale. All solid white scale bars are 10 μm, and the window size is 80 μm by 80 μm.](aax4769-F2){#F2}

We compared the in-plane locations of the focal spots between the simulation and experiment. The simulated (red) spots shown in [Fig. 3A](#F3){ref-type="fig"} lie on the dashed black circle with a 12-μm radius. The experimental (blue) spots largely trace out the same shape and qualitatively behave according to the defined FOM. However, there are some deviations from the simulated spots, e.g., the first and last spots do not lie perfectly on top of each other. The position differences between the simulated and experimentally measured focal spots are shown in [Fig. 3B](#F3){ref-type="fig"}. These discrepancies can be attributed to alignment error during optical characterization, in addition to the fabrication imperfections, as noted earlier.

DISCUSSION
==========

![Focal spot location and error.\
(**A**) Comparison of the simulated focal spot positions and the experimental positions based on the location of maximum intensity. Simulated (experimental) data are plotted in red (blue). The dashed black line is a circle of radius 12 μm serving as an eye guide. (**B**) Relative positional error of each of the focal spots. Numbers correspond to the order in which spots appear, with 1 being the closest focal plane (100 μm) and 8 being the furthest (300 μm).](aax4769-F3){#F3}

Currently, the main drawbacks associated with this design method are its restriction to spherical scatterers, the low refractive index of the available Nanoscribe GT resists, and the nonconvex nature of the optimization. The GMMT can be extended into the T-matrix theory (TMM), which supports particles of arbitrary geometry such as cylinders and ellipsoids ([@R31]--[@R33]). In particular, extending the method to cylindrical scatterers would allow the method to be applied to conventional binary semiconductor manufacturing processes. Furthermore, a plane wave coupling method is needed to improve the accuracy and use of the TMM to be compatible with cylindrical scatterers and substrates as the TMM becomes unstable for the closely packed ensembles of nonspherical scatterers ubiquitous in discrete scatterer optics ([@R34], [@R35]). In addition to the plane wave coupling method, it is also possible to use spheroidal basis functions to allow for closer scatterer packing ([@R36]). Besides the T-matrix approach, the PSF of an optical element could be designed using a unit cell--based optimization procedure. While these methods currently rely on the local phase approximation, they can potentially enable finer control over the geometry of individual scatterers ([@R20], [@R23], [@R37]). In addition, materials with high refractive indices have been shown to be an important factor determining the optical element's quality ([@R38]). Unfortunately, the resists compatible with the Nanoscribe currently all have a relatively low refractive index when compared with materials commonly used for metasurface optics. Thus, expanding the applicability of this method by using the TMM in addition to the GMMT not only increases the library of shapes available for optimization but also indirectly allows a wider selection of high-index optical materials. Last, the nonconvex nature of this optimization problem only guarantees convergence of gradient methods to a local optimum. This is a well-documented limitation of these methods and can be ameliorated by FOM design or suitable initial condition choices ([@R39]). Despite this limitation, gradient-based inverse design methods have been shown to demonstrate high-performance optical devices in many different applications ([@R16]--[@R23]).

In conclusion, we demonstrate a design method using inverse Mie scattering that allows specifying the optical fields in three dimensions and is well suited for large arrays of discrete dielectric scatterers. This design method efficiently takes advantage of the large number of degrees of freedom available in these scatterer arrays and tunes individual scatterer properties to optimize their performance. We validate this method by fabricating and testing a novel optical element that produces a discrete helical optical focusing pattern. Such an element creating lens-like PSF that varies predictably with different values of defocus would find applications in imaging systems concerned with extracting features at specific depth planes without performing a deconvolution operation. Some example applications of the reported device include depth-sensing cameras for autonomous navigation ([@R40]) and bright-field imaging of larger biological structures ([@R26]). The demonstrated optical elements have no analog in traditional optics and are impractical to design via intuition. We note that this method can be readily extended to different optical elements in a straightforward manner by including different geometries and suitable FOM. This work is a step toward enabling a flexible model of "designer" optics where optical elements can be exquisitely tailored to a user-defined functionality.

MATERIALS AND METHODS
=====================

Fabrication
-----------

We fabricated all devices in the Washington Nanofabrication Facility at the University of Washington, Seattle. We used the Nanoscribe GT and IP-Dip with the 63× objective for the fabrication of all samples. We converted the devices from .stl files using Nanoscribe's proprietary DeScribe software and then fabricated them on high-resolution glass substrates provided by Nanoscribe. In general, exposure took around 20 to 30 min for each device measuring about 144 by 144 μm^2^. Following exposure, we developed the samples for 20 min using the MicroChem SU-8 developer. We then rinsed them in isopropyl alcohol and deionized water.

Measurement
-----------

We measured the devices with two different setups corresponding to two design wavelengths. We measured the 1.55 μm device using a Santec TSL-510 centered at 1.55 μm emitting 15-mW continuous wave. We excited the device with normally incident illumination which was collected by a microscope with a movable objective mounted on a stage with a micrometer with 5-μm step size. The microscope was composed of an infinity-corrected 40× Nikon Plan Fluor and a 20-cm focal length lens. The detector used was a Xenics Bobcat-6583 infrared camera with a 320 × 256--pixel array. The camera exposure time was set to be 10 μs.

For the 3 μm device, the midwave infrared idler of a Ti-sapphire pumped optical parametric oscillator (OPO), M Squared Firefly IR, was used to illuminate the metalens at normal incidence. A 0.56 numerical aperture GeSbSe antireflective-coated infrared asphere (effective focal length, 4 mm; working distance, 3.05 mm) was translated in steps of 0.01 μm along the optical axis behind the metasurface. The chiral focal spots of the metasurface at λ = 3 μm were imaged directly onto an InSb 640 × 512 pixel focal plane array (FPA), FLIR SC6700, with a 15-μm pitch cooled to 76 K. This resulted in a magnification of ≈56.2. The exposure time of the camera was set to 0.9 ms, integrating 135 pulses from the 150.6 kHz (\<10-ns pulses, 105-mW average power) from the OPO per frame. To maximize the 14-bit dynamic range of the FPA, a neutral density 1 filter was used prior to the metasurface.

Supplementary Material
======================

###### http://advances.sciencemag.org/cgi/content/full/5/10/eaax4769/DC1

###### Download PDF

###### Controlling three-dimensional optical fields via inverse Mie scattering

**Funding:** The research is supported by NSF grant number 1825308, Samsung GRO grant, the UW Reality Lab, Facebook, Google, and Huawei. Part of this work was conducted at the Washington Nanofabrication Facility/Molecular Analysis Facility, a National Nanotechnology Coordinated Infrastructure (NNCI) site at the University of Washington, which is supported, in part, by funds from the NSF (awards NNCI-1542101, 1337840, and 0335765), the NIH, the Molecular Engineering & Sciences Institute, the Clean Energy Institute, the Washington Research Foundation, the M. J. Murdock Charitable Trust, Altatech, ClassOne Technology, GCE Market, Google, and SPTS. The Nanoscribe GT used in this work is provided via the NSF-MRI-1624513 grant. J.R.H. acknowledges support from the Air Force Office of Scientific Research under contract number FA9550-15RYCOR159. **Author contributions:** A.Z., A.M., and J.R.H. conceived of the project. A.Z. implemented the method and performed all simulations and design. A.Z. and J.W. fabricated all of the structures. A.Z., R.G., and E.S. performed all of the experiments. A.Z. wrote the manuscript with comments from all the other authors. **Competing interests:** The authors declare that they have no competing interests. **Data and materials availability:** All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper and the MATLAB implementation of the method may be requested from the authors.

Supplementary material for this article is available at <http://advances.sciencemag.org/cgi/content/full/5/10/eaax4769/DC1>

Section S1. The GMMT and specification of FOM

Section S2. Derivation of the derivative

Section S3. Optimization parameters

Section S4. Mie coefficients

Section S5. Effect of fabrication defects

Section S6. The 3 μm experimental and simulation data

Section S7. Computing resources

Fig. S1. Optimization scheme.

Fig. S2. Mie coefficients.

Fig. S3. SEM of the initial device.

Fig. S4. Simulated and experimental device performance.

Fig. S5. Extracted focal spots and locations.

Fig. S6. In-plane focal spot comparison.

Fig. S7. The 3 μm simulated and experimental device performance.

Fig. S8. The 3 μm device SEM.

Table S1. Benefits and drawbacks of the design approach.

Table S2. Optimization parameters.

References ([@R41], [@R42])
