HTML
-
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-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)$.
-
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 x0 using a 2nd-order Taylor series expansion as follows:
$$ \phi \left( x \right) \approx \phi \left( {x_0} \right) + \phi ' \left( {x_0} \right)\left( {x - x_0} \right) + \frac{1}{2}\phi '' \left( {x_0} \right)\left( {x - x_0} \right)^2 $$ (1) This section can be approximated as a line of slope $\phi ' \left({x_0} \right)$ with a phase offset of ϕ(x0) + Δϕ, which incurs an error of ${\it{\epsilon }}\left(x \right) = \phi '' \left({x_0} \right)\left({x - x_0} \right)^2/2 - {\mathrm{\Delta }}\phi$. Given a section of length $d$ at position x0, the root-mean-square (RMS) wavefront error is minimized when ${\mathrm{\Delta }}\phi = \phi '' \left({x_0} \right)d^2/24$ and is expressed as follows:
$$ {\it{\epsilon }}_{{{rms}}}\left( {{\mathrm{\Delta }}x} \right) = \frac{1}{{12\sqrt 5 }}\phi '' \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.
Fig. 2 Impact of phase profile linearization on the cylindrical lens performance.
a Cylindrical lenses of focal length 36λ and NA of 0.7 are constructed using linear sections of length d for various values of d. b Line scans of the field intensity at the focal planes of the lenses. Lenses that are linearized with sections that are smaller than 4λ have peak intensities that are within 1% of that of the ideal lensOur 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) -
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 Etgt 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 Efwd and phase ϕfwd. Near-to-far-field transformations from the forward simulations are used to evaluate Efwd 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 A1 and A2 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.
-
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-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 ${{\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:
$$ {{\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:
$$ {{\Phi }}_s\left( {x, y, z} \right) = \mathop {\sum }\limits_p B_pe^{ - ik\gamma _pz}\mathop {\sum }\limits_{nm} {{\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.
-
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.
Fig. 3 Metasurface element optimization.
a At the beginning of the optimization process, the initial dielectric distribution is a random dielectric continuum. The dashed lines indicate the desired scattering angle of 20° and phase of π/2. After 10 iterations, the scattering profile is already highly directional. After 100 iterations, the optimization is complete and the metasurface section is a binary structure of silicon and air that supports the desired scattering metrics. b An intensity plot that shows the scattered fields of the optimized metasurface elementAn 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.
-
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.
Fig. 4 Simulation results of optimized metalenses.
a Relative and absolute efficiencies for metalenses that were designed with various numerical apertures. b A full-field simulation of a metalens with an NA of 0.9The 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.
-
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.
Fig. 5 Experimental characterization of fabricated metalenses.
a Scanning electron micrographs of the metalenses with tilted and top-down views. b Relative and absolute efficiencies of the fabricated metalenses, along with their simulated values. c-e Intensity line scans at the focal planes of the metalenses with NAs of 0.2, 0.5, and 0.8, respectively, along with comparisons with the simulated lenses. f-h Efficiencies of the fabricated metalenses as a function of the wavelength, along with their simulated values. These metalenses have the same design but were fabricated on a different sample than those that were measured in (a-e)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.