Topology optimization of freeform large-area metasurfaces

TL;DR

This paper presents a scalable topology optimization method for large-area, multi-frequency metasurfaces, enabling the discovery of complex designs with unprecedented degrees of freedom and multi-layer configurations.

Contribution

It extends topology optimization to large metasurface areas using a locally periodic approximation, allowing for the design of both sub- and supra-wavelength patterns across multiple frequencies and angles.

Findings

Successfully optimized metasurfaces with over 10^6 degrees of freedom.

Discovered novel multi-layer, multi-frequency metasurface designs.

Achieved designs that control both near and far fields.

Abstract

We demonstrate optimization of optical metasurfaces over $10^5$--$10^6$ degrees of freedom in two and three dimensions, 100--1000+ wavelengths ($\lambda$) in diameter, with 100+ parameters per $\lambda^2$. In particular, we show how topology optimization, with one degree of freedom per high-resolution "pixel," can be extended to large areas with the help of a locally periodic approximation that was previously only used for a few parameters per $\lambda^2$. In this way, we can computationally discover completely unexpected metasurface designs for challenging multi-frequency, multi-angle problems, including designs for fully coupled multi-layer structures with arbitrary per-layer patterns. Unlike typical metasurface designs based on subwavelength unit cells, our approach can discover both sub- and supra-wavelength patterns and can obtain both the near and far fields.

Figures

Equations

$$\displaystyle\max_{\big{\{}\bar{\epsilon}(\mathbf{r})\big{\}}}f\left(\mathbf{E};\bar{\epsilon}\right)$$

$$0$$

$$\nabla \times {\mu }^{-1}\nabla \times \mathbf{E}-{\omega }^2ϵ\mathbf{E}=-i\omega \mathbf{J}$$

$${\nabla }_{\overline{ϵ}}f(\mathbf{r})$$

$$g$$

$$\nabla \times {\mu }^{-1}\nabla \times \tilde{\mathbf{E}}-{\omega }^2ϵ\tilde{\mathbf{E}}={\mathbf{G}}_0\left({\mathbf{r}}_0,\mathbf{r}\right)\delta (\mathbf{r}-{\mathbf{r}}_{\mathbf{s}}).$$

$$\displaystyle\max_{\bar{\epsilon}}\min_{i}\Big{\{}f_{\omega_{i},\mathbf{J}_{i}}\Big{\}},$$

$$\underset{\overline{ϵ},t}{\max }$$

$$t-f_{{\omega }_i,{\mathbf{J}}_i}$$

$$\displaystyle r(\epsilon)=\Big{|}\nabla\times\mu^{-1}\nabla\times\mathbf{E_{\text{LPA}}}-\omega^{2}\epsilon\mathbf{E_{\text{LPA}}}+i\omega\mathbf{J}\Big{|}^{2}$$

Full text

1 Introduction

We present a method for large-area ($\geq 100$ wavelengths $\lambda$) topology optimization (TO) [1, 2, 3] of optical “metasurfaces” [4, 5] with millions of degrees of freedom (DoF) determined automatically. Whereas previous metasurface design methods used only a few parameters per subwavelength unit cell [5, 6, 7, 8, 9, 10, 11, 12, 13], TO allows us to consider thousands of parameters per unit cell, with much larger ($\sim 5\lambda$) unit cells in which both sub-wavelength and supra-wavelength features are discovered by making every “pixel” a DoF. Whereas previous TO methods in optics were restricted to computational domains $\lesssim 10\lambda$ amenable to brute-force simulation [14, 15], we exploit a locally periodic approximation [6, 7, 8, 9, 10, 11, 12, 16] to combine many smaller simulations into a single optimization problem for a huge surface (Sec. 2). We validate example designs for high numerical-aperture (NA) multi-frequency/multi-angle lenses in 2d against full-wave simulations of the entire domain (Sec. 3.1 and 3.2). We also present example 3d designs for monochromatic high-NA lenses with $\sim 10^{6}$ DoF (Sec. 3.3), enabled by an efficient massively parallel implementation combining thousands of RCWA (rigorous coupled-wave analysis) unit-cell solutions [17, 18] and fast adjoint-method [1, 3] gradient computations. Because our method employs the full scattered field for each surface unit, as opposed to a single phase or amplitude in previous works [5, 6, 7, 8, 9, 10, 11, 12], we could also apply our method to near-field or guided-wave optimization [19] (Sec. 4).

Metasurfaces offer ultra-compact flat-optics alternatives for traditional bulky systems used in lensing, beam-shaping, holography and beyond [20]. The power and versatility of a metasurface resides in aperiodically patterned nanophotonic features covering hundreds or thousands of wavelengths in diameter. The term “meta" often refers to extremely subwavelength features that may be modeled by effective surface impedances [21], but in this work we consider both sub- and supra-wavelength features and do not employ any effective-medium approximation. Understandably, a meta-device poses a more challenging design problem than a traditional optical system since it requires an understanding and control of electromagnetic interactions within a vast number of nanostructures. Typically, such interactions are handled by a locally periodic approximation (LPA) in which the metasurface is divided into computationally tractable independent unit cells with periodic boundary conditions [5, 6, 7, 8, 9, 10, 11, 12]. One extreme limit of LPA is scalar diffraction theory [22], in which the surface is treated as locally uniform; this is often used to design diffractive surfaces [23] but is unlikely to accurately model subwavelength patterns. Meanwhile, the geometric library which provides the “meta-elements" is made up of primitive shapes such as subwavelength cylindrical or rectangular pillars which can be rapidly modeled (under LPA) by an electromagnetic solver such as RCWA [5, 6, 7, 8, 9, 10, 11, 12]. However, it may become increasingly difficult and ultimately infeasible to leverage only simple primitive geometries for more complex problems such as broadband achromatic focusing [24, 25, 26] or controlled angular dispersion [27] which impose stringent demands on the local phase shifts provided by the meta-elements.

Topology optimization (TO) is a large-scale computational technique that can handle an extensive design space, considering the dielectric permittivity at every spatial point as a DoF [1, 2, 3]. In contrast to the heuristic search routines regularly employed by the photonic community such as genetic algorithms [28] or particle swarm methods [29], TO employs gradient-based optimization techniques to explore hundreds to billions of continuous DoFs. Such a capability is made possible by a rapid computation of gradients (with respect to all the DoFs) via adjoint methods [30], which, in the context of electromagnetic inverse design, involve just one additional solution of Maxwell’s equations [3]. In fact, such techniques have been gaining traction lately in the field of photonic integrated circuits, especially for the design of compact modal multiplexers and converters [31]. Only recently, TO-based inverse design methods have been extended to metasurfaces, particularly in the context of freeform meta-gratings [14] and metalenses with angular phase control [15]. While such applications reveal TO as a very promising tool for realizing increasingly sophisticated meta-devices, they have been limited to small computational domains $\lesssim 10\lambda$. In this paper, we present a combination of TO and LPA to efficiently design large-area metasurfaces with enhanced functionalities.

2 Formulation

Because we employ a general optimization framework, we have great flexibility in what function of the electromagnetic fields to optimize (our “objective” function $f$) [16]. For many light-focusing problems, it is convenient to simply maximize the electric-field intensity at a focal point $\mathbf{r_{0}}$ [16]:

[TABLE]

Here, the set of DoFs \big{\{}\bar{\epsilon}(\mathbf{r})\big{\}} is related to the position-dependent dielectric profile via $\epsilon(\mathbf{r})=\left(\epsilon_{\text{st}}-\epsilon_{\text{bg}}\right)\bar{\epsilon}\left(\mathbf{r}\right)+\epsilon_{\text{bg}}$, where $\epsilon_{\text{st (bg)}}$ denotes the relative permittivity of the structural (background) dielectric material. Whereas $\bar{\epsilon}$ is a continuous DoF allowed to have intermediate values between 0 and 1, we employ Heaviside projection filters [1] to ensure that the final optimal design is binary, i.e., $\bar{\epsilon}_{\text{optimal}}\in\{0,1\}$. In this work, we consider only the binary filter to theoretically illustrate our design method. More generally, a variety of additional geometric filters and constraints can be straightforwardly incorporated into the formulation, including regularization filters and curvature constraints [1, 32, 33] to impose design robustness and minimum feature sizes conducive to fabrication. The objective function $f$ represents the far-field intensity obtained by convolving the free-space Green’s function $\mathbf{G_{0}}$ with the near-field $\mathbf{E}$ at some reference plane $\mathbf{r_{s}}$ over the entire metasurface. (To be precise, the free-space Green’s function is convolved with equivalent surface currents which are given by the tangential field components [34].) Here, $\mathbf{E}$ is the steady-state solution to the frequency-domain Maxwell’s equation:

[TABLE]

in response to the incident current $\mathbf{J}$ at a frequency $\omega$ where we set the magnetic permeability $\mu=1$ for typical optical materials. To efficiently model a large device area, we use the locally periodic approximation [16], in which the metasurface is broken up into multiple smaller cells, each of which is independently simulated with periodic boundary conditions (Fig. 1). The total near-field $\mathbf{E}$ is, therefore, approximated by the set of disjoint intra-cell electric fields which must be engineered such that they all constructively add up at the focal point. The optimization typically involves thousands of DoFs which must be efficiently handled by a numerical algorithm. Within the scope of TO, this requires efficient calculations of the derivatives ${\partial f\over\partial\bar{\epsilon}}$ at every pixel point in the device region, which we obtain by an adjoint method [3, 30]:

[TABLE]

where the adjoint field $\mathbf{\tilde{E}}$ is evaluated by solving an additional Maxwell’s equation:

[TABLE]

We emphasize that our framework is entirely different from the common approach where an ideal phase profile to be realized is approximately fitted with zeroth-order phase shifts allowed by primitive scattering elements, each residing within a sub-wavelength unit cell [5, 6, 7, 8, 9, 10, 11, 12]. In contrast, our approach greatly broadens the structural design space by considering non-intuitive shapes and forms as well as supra-wavelength unit cells within which the amplitude and phase of the electric field may vary considerably and must be fully taken into account during optimization.

2.1 Numerical Implementation

To ensure a speedy computation, we numerically implement the optimization problem (1) in C. Each of the cells is either solved by a C implementation of the finite-difference frequency-domain method (FDFD) [35] or the rigorous coupled-wave analysis (RCWA) [18]. To model a large surface area rapidly, all the cells are independently simulated in an “embarrassingly parallel” fashion using the message-passing interface (MPI) library [36, 37]. The simulated results are then consolidated to compute the global objective and gradients. The structural update during optimization is provided by a standard gradient-based nonlinear-optimization algorithm [38, 39] with a free-software implementation [40]. We note that the discretization error in the gradients arising from the discrepancy between the formal adjoint method and the numerical Maxwell’s solver (such as RCWA) is practically negligible ($\lesssim 0.01\%$).

3 Applications

The major advantage of a large-scale computational approach consists in its flexibility for handling multi-objective problems, such as those prescribing the behavior of a single metasurface for a set of frequencies and incident current conditions $\{\omega_{i},\mathbf{J}_{i}\}$. In such cases, the optimization problem can be easily extended by the maximin formulation:

[TABLE]

which can be recast into an equivalent differentiable form [41]:

[TABLE]

3.1 Metalens Concentrator

To demonstrate the efficacy of our design technique, we optimize a multi-layered 2D metalens (Fig. 2) that can concentrate incoming light at 11 incident angles \{\theta~{}(^{\circ})=3i~{}\big{|}~{}i=0,1,...,10\} to the same focus. The lens consists of five layers of TiO2 (refractive index $n\approx 2.4$), four of which are buried in silica ($n\approx 1.5$). The diameter of the lens is $D=1200\lambda$, where $\lambda$ is the operational wavelength, and the focal length is $F=1000\lambda$, corresponding to a rather high NA of 0.51. For simplicity, we consider the electric field to be polarized along the $y$ axis. During the optimization, the entire lens is broken up into 240 cells, each of which is $5\lambda$ long, for $3\times 10^{5}$ DoFs in total, while the entire computation is parallelized over 1200 CPUs. While we perform the optimization with the locally periodic approximation, we validate the optimized result with rigorous full-wave FDTD simulations [42] (Fig. 2) in which we compute the actual electric fields over the entire metasurface without any uncontrolled approximation. The lens is shown to exhibit diffraction-limited focusing at all the optimized angles with the full-width-at-maximum (FWHM) of $\sim\lambda$ while the average transmission efficiency is found to be $T\approx 40\%$ . The average efficiency of $40\%$ seems to violate Lorentz reciprocity and, in particular, the concentration bound $\leq{1\over 11}$ predicted by [43]. However, there is no such violation because, although the output field profiles have the same prominent peak at the optical axis, they are never entirely alike for any two input angles due to the weak but distinct diffracting patterns away from the axis.

3.2 Multi-wavelength Focusing

Another example we consider is a partially achromatic 2D metalens (Fig. 3) with discrete chromatic aberration corrections at three wavelengths $\{\lambda_{1},\lambda_{2},\lambda_{3}\}$ with $\lambda_{2}=\lambda_{1}/1.2$ and $\lambda_{3}=\lambda_{1}/1.38$. In particular, one may choose $\lambda_{1}=650~{}\mathrm{nm}$, $\lambda_{2}=540~{}\mathrm{nm}$ and $\lambda_{3}=470~{}\mathrm{nm}$, corresponding to red, green and blue (RGB) wavelengths. The lens consists of two layers of TiO2 ($n\approx 2.4$), one of which is buried in silica ($n\approx 1.5$). The diameter of the lens is $D=200\lambda_{1}$ and the focal length is $F=100\lambda_{1}$, corresponding to $\mathrm{NA}=0.71$. During the optimization, the entire lens is broken up into 40 cells, each of which is $5\lambda$ long, for $1.6\times 10^{4}$ DoFs in total, while the entire computation is parallelized over just 20 CPUs. At the desired frequencies, the lens is shown to focus a normally incident $y$-polarized beam at or near the diffraction limit with the FWHMs of $\approx 0.71\lambda_{1},~{}0.72\lambda_{2}$ and $0.76\lambda_{3}$ respectively and with an average transmission efficiency of $\sim 71\%$. Recently, broadband achromatic metalenses have been gaining widespread attention for full color imaging; although the state-of-the-art designs continuously cover the full visible spectrum, they are mostly limited to smaller lens diameters and/or low NA [24, 25, 26]. Our results suggest that multi-layer designs enabled by TO could alleviate many of these limitations by allowing more complex geometries to be explored.

Generally, a high-NA multi-wavelength design is a challenge for LPA because the design rapidly varies across the surface, resulting in the scattering “noise” visible in a full-wave simulation of the entire structure but not captured by the conjoined LPA solution (Fig. 3), in contrast to the multi-angle case in Fig. 2. One consequence is that the focusing efficiencies at the three wavelengths (defined as the integrated field intensity around the focal spot divided by the total intensity over the image plane [44]) are found to be $\sim 51\%,~{}50\%$ and $45\%$ while LPA predicts $\sim 78\%,~{}55\%$ and $60\%$ respectively. Note that our objective function (1) merely maximizes the intensity at the focal spot but does not necessarily suppress stray diffractions elsewhere which become apparent even in LPA predictions (Fig. 3). If desired, the optimal design can be further refined by improving the LPA wavefront by minimizing a phase error [15] or wavefront error objective [16]. At the same time, deviations from LPA may be mitigated by considering larger metasurfaces where LPA is more appropriate, by including higher-order corrections to LPA [19], or by incorporating optimization constraints designed to enforce the validity of LPA. For example, the optimization process may be augmented with constraints restricting the structural variations between neighboring cells or specifying a bound on the physical residual:

[TABLE]

where the double curl operator is applied over the whole surface without LPA. Note that evaluating the residual $r$ requires just a single matrix-vector multiplication, albeit involving a very large and sparse matrix, and can be made relatively cheap and fast when parallelized over several cores.

3.3 Three-dimensional Examples

Next, we turn to proof-of-concept 3D examples. First, we consider a collection of 3D cells, each of which has an area $L_{x}\times L_{y}=\lambda\times\lambda$. The cells are repeated along $y$ but are allowed to differ from each other along $x$. Fig. 4 shows the corresponding “cylindrical” lens (NA=0.71) comprising such cells which focuses a normally-incident $y$-polarized beam to a focal line running parallel to the $y$-axis. The lens, though periodic along $y$, shows complex geometric variations along both $x$ and $y$ within each cell. A rigorous FDTD simulation of the whole lens shows the focusing behavior with a transmission efficiency of 75%. Although the design of a monochromatic lens, no matter the geometric complexity, is only useful to illustrate the method, we note that our framework enables a relatively cheap and fairly quick computation using as few as 25 CPUs over a span of a few minutes while considering the entire structure with as many as 200 cells and $4\times 10^{4}$ DoFs in total.

Lastly, we consider the full 3D lens (Fig. 5) with cells varying along both $x$ and $y$. The lens has 6400 $\lambda\times\lambda$ cells, corresponding to NA=0.37, with $6.4\times 10^{5}$ DoFs in total. Again, our framework enables a cheap and speedy design of the entire lens, utilizing 40 CPUs and lasting only a few minutes. A major computational hurdle actually arises after the optimization: given the enormous size, the lens cannot be readily simulated via FDTD using a few CPUs. Instead, we provide an approximate far-field profile computed from the LPA combination of the unit-cell simulations, which has been shown to give excellent accuracy for monochromatic designs [13] and good accuracy even for the complex multi-frequency designs in Sec. 3.2.

These preliminary results suggest promising ways to scale up the 3D metasurface design. In particular, the computational workload can be reduced by orders of magnitude via rotational symmetry, specifically, by considering just a single row of cells and rotating it to fully cover the entire surface. Such an approach is, in fact, a more common practice in metalens design [5, 6, 7, 8, 9, 10, 11, 12, 13] rather than considering all the cells as we did here. Also, polarization insensitivity could be ensured by further imposing appropriate symmetries, such as $C_{4v}$ within each cell. Combined with such symmetry-related reductions, an access to several hundreds or low thousands of CPUs (readily available on institutional supercomputers or commercial cloud computing services) amply allow for the design of high-NA wide-area metalenses with greatly enhanced functionalities including broadband achromaticity and large-angle aberration corrections.

4 Conclusion and Outlook

We have provided an efficient computational framework for the inverse design of large-area freeform metasurfaces. We have also presented various examples including multi-angle and multi-wavelength metalens designs that demonstrate the versatility and power of our method. We expect that a large-scale optimization technique like ours may become indispensable for tackling challenging problems which inherently call for a large design space with multiple layers. In particular, our method may greatly benefit the design of a polarization-insensitive large-area high-NA single-piece metalens with chromatic and achromatic aberration corrections over the entire visible spectrum and over a large field of view. On the other hand, our framework need not be limited to metalens design. In fact, the particular objective we presented above might be more naturally suited for three-dimensional sculpting of the far-field intensity (or thin-film holography). Another intriguing application particularly suitable for the RCWA-based optimization framework is to inverse-design the far field between two meta-devices separated by thick (many-wavelengths) spacing layers. Ultimately, one would like to address the most challenging problem of electromagnetic design: to intimately manipulate the near-field physics of large aperiodic structures in which ultra-rich modal interactions involving both extended and localized resonances may offer a superb playground for exploring novel phenomena and devising next-generation technologies such as solid-state Lidar systems [45], strongly enhanced light-matter interactions and nonlinear processes [46, 47, 48, 49, 50], near-field radiative heat transfer [51], and opto-nanomechanical engineering [52, 53]. While LPA offers asymptotic advantages in the limit of “adiabatic” nearly periodic surfaces, our optimization framework also provides a promising pathway for handling a wider variety of large aperiodic systems by incorporating additional constraints (Sec. 3.2) or by augmenting LPA with domain decomposition concepts [54] including using the conjoined LPA solution as a preconditioner, using other boundary conditions (e.g. mirror boundaries), and using overlapping domains [55, 56].