## Abstract

We present a computational framework for efficient optimization-based “inverse design” of large-area “metasurfaces” (subwavelength-patterned surfaces) for applications such as multi-wavelength/multi-angle optimizations, and demultiplexers. To optimize surfaces that can be thousands of wavelengths in diameter, with thousands (or millions) of parameters, the key is a fast approximate solver for the scattered field. We employ a “locally periodic” approximation in which the scattering problem is approximated by a composition of periodic scattering problems from each unit cell of the surface, and validate it against brute-force Maxwell solutions. This is an extension of ideas in previous metasurface designs, but with greatly increased flexibility, e.g. to automatically balance tradeoffs between multiple frequencies or to optimize a photonic device given only partial information about the desired field. Our approach even extends beyond the metasurface regime to non-subwavelength structures where additional diffracted orders must be included (but the period is not large enough to apply scalar diffraction theory).

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

## 1. Introduction and motivation

In this paper, we present and validate a fast method for optimization-based “inverse design” of large (hundreds of wavelengths *λ*) aperiodic metasurfaces for wavefront shaping [1–5], incorporating both scattered amplitude and phase for multiple incident *λ* and angles. Previous methods either optimized the full Maxwell equations [6–10] (which is infeasible for large surfaces), were restricted to weakly coupled scatterers [11], or started with a desired scattered phase and tried to design a corresponding metasurface unit cell [12–21] (but if attainable unit cells fail to exactly match the desired *λ*-dependent phase there was no systematic way to choose the best compromise). In contrast, our approach starts with a family of manufacturable unit cells and directly optimizes an aperiodic composition for the desired field pattern by a fast approximate model, automatically finding the best compromise for the given constraints. Whereas phase-design methods typically assume that the desired scattered field is known everywhere [12–21], our approach allows one to specify the field objective only in regions of interest. As outlined in Fig. 1, given exact scattering calculations for small metasurface unit cells (Sec. 2.1 and Fig. 2), we build an approximate convolutional model of an arbitrary metasurface (Sec. 2.2) that can then be optimized (Sec. 3) rapidly (seconds to find the optimum for a 200-*λ* 2d aperiodic surface with hundreds of parameters) using two different objective functions. We validate the optimized design with a brute-force Maxwell solver (Sec. 3.1) and we find excellent quantitative agreement (Fig. 4) even for rapidly varying aperiodic surfaces that challenge the assumptions of our model. We present example designs (Sec. 4) for a multi-wavelength optimization (Fig. 5), a wavelength demultiplexer (Fig. 6), and an multi-angle optimization (Fig. 7). Our approach is not limited to true “metasurfaces” whose features are small enough to mimic effective-impedance [22–30] surfaces: we show that it even works well for large-period microstructures that scatter multiple diffracted beams. Indeed, our method is easily extensible to incorporate multiple diffraction coefficients, multi-layer/multi-parameter unit cells, multiple polarizations, and other complications (Sec. 5 and Fig. 8).

## 2. Locally periodic approximation

The key to “metasurface” design is to be able to quickly calculate the transmitted/reflected field for a large-area structure, possibly thousands of wavelengths in diameter—too large to solve the full Maxwell equations without some simplifying assumption. Similar to [12–20], the central approximation of our approach is to assume that the metasurface is *locally periodic*: the scattering in any small region is almost the same as the scattering from a periodic surface. The use of periodic calculations to compute the specular reflection phase only, typically discarding amplitudes and additional diffracted orders, has sometimes been called a “local phase approximation” [31,32]. (Contrast this with the regime of scalar diffraction theory [33,34], valid for period Λ ≫ wavelength *λ*, in which the surface is treated as locally *uniform*, separately computing the transmission coefficient at each *point* on the surface.) In a separate paper [35], we develop the rigorous foundations and convergence rates of a related approximation, along with higher-order corrections, but here (similar to previous authors [12–19]) we will simply perform brute-force Maxwell simulations at the end (Sec. 3) to validate our designs. (In fact, we will find in Fig. 4 that the locally periodic approximation gives excellent agreement with full simulations even for surfaces where the unit cell is rapidly varying in some regions.) Unlike previous authors who calculated only the scattered phase and not the total scattered field [12–19], we employ both amplitude and phase information to formulate a complete approximate solver (scattered field for any given incident field) that can be used to optimize the metasurface for arbitrary “objective” functions of the field. Not only does this approximate solver enable very general optimization, it also allows us to evaluate the optimized metasurface for different (non-optimized) incident fields (e.g. the wavelength sensitivity computed in Sec. 4 for the multi-wavelength lens). As discussed in Sec. 5, our approach is also easily extensible to a "non-metasurface" regime in which there are multiple diffracted beams from a large-period surface, as well as to computing near fields (via the scattering coefficients of the evanescent waves). In this section, we explain in detail how the locally periodic approximation allows us to compute the total scattered/reflected field (at any point in space) for any incident wave.

To make it easier to understand our approach, it is helpful to consider a specific example of a two-dimensional “metasurface” unit cell, based on [17]: TiO_{2} pillars on top of a silicon dioxide substrate, as shown in Fig. 2. The height of the pillar is fixed to 600 nm, the period *a* is fixed to 235 nm, and the pillar width varies: *p* ∈ [50, *a* −50] nm (imposing a minimum feature size of 50 nm for practical fabrication). One could easily add more parameters and/or constraints, as discussed in Sec. 6. Given this unit cell, an aperiodic metasurface is formed by taking a group of such unit cells with independent parameters and juxtaposing them next to each other.

Our goal is to compute the scattered (transmitted or reflected) field for such an aperiodic surface, for any given incident wave (e.g. a planewave or gaussian beam), given *only* the exact Maxwell solutions for scattering of planewaves by *periodic* surfaces of the different pillar widths. In this paper, we consider only incident propagating (not evanescent) waves, but in another paper [35] we show that a similar approach can be extended to evanescent fields as well. The key “locally periodic” assumption is that the pillar width (the unit cell) changes sufficiently slowly from one pillar to the next. (This assumption is rigorously quantified in [35], is validated numerically in Sec. 4, and it turns out that we even obtain good accuracy when there are sudden changes in pillar width at a few locations.) As mentioned above, this assumption is similar in spirit to other metasurface work [12–19], where it was found to work well for a wide variety of metasurface designs; the main contribution of this paper is to couple the locally periodic approximation to general optimization tools and near-to-far-field transformations. Of course, this approximation can break down for devices requiring extremely rapid surface variations such as diffraction to nearly glancing angles [35, 36], although it can be generalized by including a next-order correction [35], but this is not a problem for the moderate-NA lens-like applications considered in this paper.

#### 2.1. Periodic sub-problems

In Fig. 2(left) is shown the fundamental assumption of our approach. For each unit cell of the aperiodic structure, we approximate the field in a plane/line just above the unit cell by the solution for the equivalent periodic structure. Three examples are highlighted corresponding to three different parameters of the unit cell. When the period of the unit cell is subwavelength, the zeroth diffractive order is the only propagating wave [37]. Therefore, if we are interested only in the far field, we can make an additional approximation: we replace the scattered field by its zeroth Fourier component which is simply the average of the field on the plane just above the pillar. Given this approximate field just above the surface, in Sec. 2.2 we construct an approximate field everywhere above the surface. In Sec. 5, we go beyond the zeroth-order (specular) approximation by including additional diffractive orders.

We will consider periodic structures with hundreds of different pillar widths (with a fixed period), but we would like to avoid having to do hundreds of unit-cell calculations. We can take advantage of the fact that the scattered fields are smooth functions of the pillar width [38,39] by solving the scattering function for a few widths and then interpolating to any other widths.

Given a smooth function *f* (*p*) of some parameter *p*_{1} ≤ *p* ≤*p*_{2}, Chebyshev methods evaluate *f* (*p*) at a few special points *p* and construct a polynomial approximation $\tilde{f}(p)$ that can be used to rapidly evaluate the *f* (*p*) with exponentially good accuracy. This can be extended to multiple parameters using products of Chebyshev polynomials [40] or by more sophisticated methods such as sparse grids [41]. In this way, we only need to solve the unit-cell Maxwell problem a few times to obtain our polynomial approximant $\tilde{f}(p)$, which is then evaluated, along with its derivative, many times during optimization. In particular Fig. 2(right) shows the real part, the imaginary part, and the phase of the zeroth Fourier coefficient of the transmitted field versus the pillar width. We evaluate this coefficient for 21 different widths (at Chebyshev points [40]), and can then interpolate to high accuracy using Chebyshev polynomials [40]. Multiple parameters per unit cell could be interpolated using a product of Chebyshev polynomials [40] or, for more than 3–4 parameters, sparse-grid interpolation [41]. Here, we use the finite-difference frequency-domain (FDFD) method [42,43] with perfectly matched layer (PML) absorbing boundaries [44,45], but any other computational method for periodic Maxwell problems would work as well.

#### 2.2. Green’s functions and the equivalence principle

Once the fields are known in a plane above the metasurface, we can obtain the fields *everywhere* above the metasurface using the *principle of equivalence* [46,47], also known as a near-to-far-field transformation [48]): the fields in the *y* = *y*_{0} plane can be treated as equivalent current sources that generate the fields everywhere else. These equivalent electric (**J**) and magnetic (**K**) current densities are defined by [47]:

*y =y*

_{0}.

A further simplification is possible if we only care to compute the fields *above* the surface. The currents (1) produce the desired scattered fields above the plane and zero fields below the plane [46,47], and this means that the same fields above are produced if we add or subtract the mirror-image currents (which produce fields below and zero above). Subtracting the mirror-image sources, however, cancels the **J** term and leaves only the **K** current (a pseudovector under mirror flips [49]). This allows us to use only $\widehat{\mathbf{n}}\times \mathbf{E}$ sources arising from the electric field computed by the locally periodic approximation in section 2.1. As explained in section 2.1, we can further approximate the **E** field by its average in each unit-cell calculation for subwavelength periods, since this gives the far-field diffracted order.

Given these equivalent currents, or their approximation by far-field locally periodic calculations, the electric (or magnetic) fields at any point **x** above the surface can be computed by integrating along with the Maxwell Green’s function (the field at **x** from a source at **x′**) [46]. For our two-dimensional model problem (*xy* plane) with the *E _{z}* polarization, where we only have a current

*K*(

_{x}**x′**) = −

*E*(

_{z}**x′**)

*δ*(

*y*−

*y*

_{0}) and we let

*G*(

**x**,

**x′)**denote the relevant of tensor component the Green’s [

*E*(

_{z}**x**) from

*K*(

_{x}**x′**)], this integral takes the form

*G*is a Hankel function [50] $G(\mathbf{x},{\mathbf{x}}^{\prime})=-\frac{ik}{4}{H}_{1}^{(1)}(kr)\widehat{\mathbf{n}}\cdot \frac{\mathbf{r}}{r}$, where $k=\frac{2\pi}{\lambda}$,

**r**=

**x − x′**, and

*r*= |

**r**| For a finite metasurface with an infinite silica substrate, we use a standard “windowing” method to truncate this integral accurately to a finite region [51–53].

This equivalent-currents formulation is exact if the true aperiodic *E _{z}* field is used for the

*K*source term, and in section 3 we find that it has excellent accuracy with the locally periodic approximation for typical metasurface designs. (A related approximation is made in scalar diffraction theory, where the locally

_{x}*uniform*approximate scattered fields can be thought of as sources producing fields everywhere else [33,34], further approximated in the far field by e.g. the Fraunhofer diffraction theory [33].)

## 3. Metasurface inverse design methods

The previous section gives us a fast way to solve the forward problem for the scattered field above a given metasurface. In this section, we see how we use it to solve the inverse problem, i.e. find the parameters of a metasurface to produce a desired scattered field. We solve it as an optimization problem: we minimize or maximize an *objective* function of the unit cell parameters subject to some constraints. Given an efficient way to compute any objective function and its gradient, there are a wide-variety of well-known optimization methods that can be applied; we use the “CCSA-MMA” algorithm [54] via a free-software implementation [55]. To avoid getting trapped in poor local optima, we use a technique called successive refinement [56–58]: We successively double the number of degrees of freedom, using the optimized coarser structures as starting points for optimization of the finer structures. (The result was not very sensitive to the starting parameter guess; we simply started each parameter in the middle of its allowed range.) What objective function should we optimize? In Sec. 3.1 we consider an objective function similar to previous work [15,21], which matches the field just above the metasurface with the desired field. In Sec. 3.2 we optimize more general functions of the scattered field, e.g. the intensity at a single focal point, which is more flexible when only partial information is known about the desired field. In both cases, in this section we will design a simple lens structure that we will validate using brute force simulations. In Sec. 4, we will consider more difficult design problems. In Sec. 3.3 we generalize our approach to multiple frequencies and angles of incidence via a maxmin formulation.

#### 3.1. Optimizing the wavefront

When the exact desired field is known everywhere above the metasurface, as in lens design [12–21] and other wavefront shaping problems, by the equivalence principle [46] it is sufficient to produce this field on a plane just above the surface. Since the approximate scattered fields *s*(*p*) just above the surface (the *E _{z}* produced for a given metasurface parameter

*p*) are given by the locally periodic approximation in Sec. 2.1, we can directly minimize the difference between this

*s*and the desired field

*a*(

*x*)

*e*

^{iϕ}^{(}

^{x}^{)}:

*s*

_{0}and

*ϕ*

_{0}are an unknown overall amplitude/phase and

*p*(

*x*) describes the metasurface parameters along the surface. This approach eliminates the need for any Green’s function integral (Sec. 2.2) to obtain the field elsewhere.

For a lens application, typically *a*(*x*) = 1 and all of the information is in the desired phase *ϕ*(*x*) [15]. A closely related approach was used for metalens design in several previous works [12–19, 21]. There, since both *a*(*x*) and the locally periodic far-field |*E _{z}|* were approximately constant, the amplitude was ignored and they simply attempted to match the desired phase. If this phase can be matched exactly in a given unit cell by tuning its parameter

*p*(

*x*) (e.g. pillar width), then no explicit optimization formulation is needed [12], but an optimization-based approach is more flexible at balancing tradeoffs in cases where the desired

*ae*cannot be exactly obtained, especially in multi-frequency problems (Sec. 4.1). A phase-based optimization approach was directly employed in [21] for topology optimization of a small area (no locally periodic approximation).

^{iϕ}For example, in Fig. 3 we minimize equation (3) for a single-frequency *λ* = 532 nm lens design problem: we focus an incident planewave at a 5-degree angle on a focal point 14.7 *µ*m from the surface, using the target phase *ϕ*(*x*) from [15]. We optimize over piecewise-constant parameters, effectively one parameter per unit cell, with a standard optimization algorithm [54] utilizing analytically computed gradients of the objective function with respect to the parameters. Starting the optimization from a constant-*p* initial guess was sufficient to obtain a local minimum with excellent performance shown in Fig. 3. (This 40 unit-cell optimization required < 100 ms on a laptop.) At the top is the |*E _{z}*|

^{2}intensity plot computed with by our approximate solver (Sec. 2.2). Below this is shown the 96% match between the desired and obtained fields (from the locally periodic approximation) just above the metasurface. At the bottom is shown the optimized metasurface geometry, which is mostly slowly varying but has sudden jumps in the pillar widths when the desired phase passes through 2

*π*.

In Fig. 4, the locally periodic approximate solver (left) is compared to a brute-force surface-integral equation (SIE) Maxwell solver [51–53] for this optimized solution, showing good quantitative agreement. More precisely, at right we compare the computed intensities |*E _{z}* (

*x, y*)|

^{2}for several separations

*y*from the surface. On the focal line, the mean squared difference between the solutions divided by the mean squared intensity is only 0.3%, validating our locally periodic approximation. The errors increase as one approaches the surface because of effects that decay with distance—scattered waves (intensity ~ 1/

*y*) from sudden jumps in the pillar width (which violate the locally periodic approximation) along with evanescent fields that we neglected in our far-field approximation—combined with the fact that small errors are more apparent in low-intensity regions far from the focal point.

#### 3.2. Optimizing arbitrary functions of the field

Alternatively, since Sec. 2.2 allows us to compute the approximate field anywhere above a metasurface, we can optimize *any* function of this field. This is especially useful if the desired field is only partially known: perhaps one cares about the field in some regions but not others, or is interested in amplitude but not phase. In particular, here we approach the lens-design problem by directly maximizing the intensity |*E _{z}* (

**x**)|

^{2}at a single focal point |

**x**|, which can be rapidly computed by a single integral (2) of the locally periodic surface fields. As in the previous sections, we used standard optimization techniques [54] with an analytically computed gradient (essentially via an adjoint method [59]), and the optimized structure for 40 unit cells was found in < 1 s on a laptop (whereas our brute-force solver was about 10

^{5}times slower). A comparison of the two methods when the period is not sub-wavelength appears in Sec. 5.

#### 3.3. Max–min multi-objective optimization

Many design problems involve a combination of multiple objectives: maximizing performance at different wavelengths, angles, and/or focal spots, for example. One common way to do this is a *max–min* formulation: we optimize the worst objective:

(For example, in Sec. 4.1, the “objective” function for an RGB lens is the intensity at the focal spot, and max–min optimization means that we try to maximize the *lowest* intensity across the three design wavelengths.) Although the expression [⋯] being maximized is no longer differentiable, which would make the most efficient high-dimensional optimization methods inapplicable, it can be transformed into an equivalent differentiable problem [60]

*t*∈ ℝ is a new “dummy” optimization parameter. Assuming that the original objective function is differentiable, we can now use a standard nonlinear constrained-optimization algorithm [54]. In particular, the CCSA-MMA algorithm [54] only requires us to supply the functions

*t*and

*t*− objective(parameters,

*λ*) and their gradients (with respect to

*t*and the parameters) in order to solve the local-optimization problem. Efficient gradient formulas for our cost functions from Sec. 3.1 and Sec. 3.2 are given in the Appendix.

We will show examples of such optimization problems in Sec. 4, where we will use max–min to optimize for multiple frequencies (Fig. 5 and Fig. 6) or angles (Fig. 7).

## 4. Applications: RGB lens, demultiplexer, and angle-insensitive lens

In this section, we show some larger and more interesting design problems that can be solved by our methods from the previous section. We still use the same TiO_{2} pillar unit cells as in Sec. 2, but now we consider metasurfaces consisting of 1000 unit cells, combining multiple frequencies and/or angles, and we could solve the resulting optimization problems in a few minutes on a laptop. In particular, we consider three applications: a lens which has the *same* focal spot for RGB (red, green, blue) wavelengths, a demultiplexer that focuses RGB wavelengths at three *different* focal spots, and a lens that focuses four incident *angles* at the same wavelength to the same focal spot. We will also show that our methods are suitable for sensitivity analyses with respect to wavelengths or angles, by evaluating designs at non-optimized inputs using our fast (locally periodic) solver.

#### 4.1. Max–min RGB (red, green, blue) focusing

Here, we use the max–min method of Sec. 3.3 to focus normally incident plane waves of three different wavelengths—480 nm (blue), 530 nm (green), and 650 nm (red)—on a *single* focal spot, by maximizing the minimum (worst) intensity at that spot for all three wavelengths. The diameter of the lens is 235 microns (1000 unit cells), and the focal length is 350.6 microns, which corresponds to a numerical aperture of 0.3.

At the bottom of Fig. 5 is shown the intensity on the focal line for all three wavelengths, demonstrating nearly diffraction-limited focusing (RGB half-maximum widths of 975, 997, and 850 nm, respectively). In Fig. 5(middle), we evaluate our optimized design along the focal axis (a fixed *x* = 0) versus distance *y* from the surface and versus wavelength across the visible spectrum, in order to show the wavelength sensitivity of our RGB design. This plot reveals that the optimized design is actually producing three different focal spots (local intensity maxima) on the focal axis for every wavelength, and at each of the RGB wavelengths a different spot is brought to the 350.6 *µ*m target. At this target focal spot, the intensity |*E _{z}*|

^{2}is plotted versus wavelength in Fig. 5(top), showing the narrowband nature of the RGB focus. The ability of our approximate solver to rapidly evaluate the performance of the design with many different (non-optimized) inputs (< 100 ms each) is a powerful tool for characterizing and understanding the metasurface.

#### 4.2. Demultiplexer

Here, we design a *demultiplexer* that focuses normally incident plane waves of three different wavelengths (RGB again) at three *different* points, which are sixty microns laterally (*x*) apart from each other on the same focal plane (again 350.6 *µ*m from the surface, a numerical aperture of 0.3). As above, we use the max–min formulation from Sec. 3.3 to maximize the worst case intensity at the focal spots.

In Fig. 6(top), we show the field intensities in the vicinity of the three focal spots for the RGB wavelengths, and in Fig. 6(bottom) we plot the corresponding intensities along the focal line *y* = 350.6 *µ*m. The focal spots for the two side focal points are tilted outward from the focal axis, which makes sense because they required off-axis focusing relative to the center of the metasurface. As in Sec. 4.1, we attain nearly diffraction-limited RGB foci half widths of 825, 785, and 795 nm, respectively.

#### 4.3. Max–min multi-angle focus

Our last application is a metasurface focusing incident plane waves coming at four different angles of incidence (normal 0°, 3°, 6°, and 9°) at the same focal point for the wavelength 532 nm, inspired by earlier topology-optimization work [21]. As in the previous sections, we target a focal length of 350.6 *µ*m (numerical aperture 0.3), and use the max–min formulation of Sec. 3.3 to maximize the worst-case focal-point intensity.

Fig. 7(right) shows the field intensities in the vicinity of the target focal spot for the four angles, exhibiting an unsurprising “tilt” proportional to the angle of incidence. As in the previous sections, the spots are nearly diffraction limited (half widths of 787, 787, 807, and 724 nm). Fig. 7(left) shows the corresponding intensities on the focal plane *y* = 350.6 *µ*m versus *x*. This plot shows that, in addition to a peak at the target point *x* = 0, the metasurface produces three auxiliary side peaks. (Preliminary work indicates that, similar to [21], these auxiliary peaks can be mostly eliminated by redesigning the unit cell via additional parameters; we will address this in a future manuscript.) That is, much like in Fig. 5, the metasurface is creating four focal spots, such that at each angle of incidence a different focal spot is brought to the *x* = 0 target point. The complex surface design and resulting transmitted field here would be very difficult to reproduce without large-scale optimization.

## 5. Beyond subwavelength periods

The term “metasurface” should strictly apply only to deeply subwavelength structures that can be accurately described by an effective surface impedance/admittance or similar [22–30], and most previous work operated in a subwavelength regime [12–19]. Conversely, when the period is larger than the wavelength, additional diffracted waves appear in the far field [37] that cannot be described by a uniform effective medium or by a single Fourier coefficient. Nevertheless, if the unit cells are mostly slowly varying it should still be valid to describe the surface by a locally periodic approximation (analogous to the adiabatic theorem for propagation through nearly periodic media [61]) to approximate the field just above the surface and hence the field everywhere as in Sec. 2. When we solve the local periodic problems in non-subwavelength structures we can no longer retain only the 0th order Fourier coefficient, but instead we must retain either the full *E _{z}* field on the surface or, for far-field calculations, the Fourier coefficients corresponding to all of the non-evanescent diffracted orders [37].

In Fig. 8 we show a single-wavelength (*λ* = 532 nm) lens design for a period of 800 nm > *λ*, so that even a periodic surface produces two additional diffracted orders 1 in addition to the 0th-order “specular” transmission. (Other than the period, the structure is±the same TiO_{2} pillar geometry considered in the previous sections, we use normal incidence, and design for a focal length of 48.6 *µ*m with 40 unit cells similar to Sec. 2.) We considered both the wavefront and the intensity optimization approaches, and validated against a brute-force Maxwell solution as in Sec. 3. Since the additional diffracted orders propagate at oblique angles, they have little influence on the focal intensity if the lens is designed to focus the 0th-order (specular) transmitted wave sufficiently far from the surface, so we carry out the inverse design using only the 0th-order term in the approximate model.

The results in Fig. 8 show that the intensity method still produces an excellent (near diffraction-limited) focal spot with high intensity that agrees well with the brute-force validation, whereas the wavefront optimization produces a much weaker focus that agrees poorly with the validation. In both cases, the brute-force calculation and the approximate solver (which includes also the diffractive orders ±1) clearly show the additional diffracted orders scattering to oblique angles that have low amplitude at the focal spot. One major problem with the wavefront approach in this geometry is that varying the pillar width in this case changes the amplitude from 1 to 0.2, very different from the constant amplitude ≈ 1 in the subwavelength case. The best phase match corresponds to a weak efficiency, whereas the intensity method can compensate by utilizing both amplitude and phase variations. The resulting lens designs shown in Fig. 8(left) correspondingly have an average amplitude twice bigger for the intensity approach than for the wavefront approach (0.8 vs 0.4). Another challenge of non-subwavelength structures, which would become more acute for larger-aperture lenses, is that large-period gratings with only a small number of parameters per unit cell cannot easily implement the rapid variations in phase that are called for by large lenses. There are too few parameters to fit the complex intra-cell phase variation.

## 6. Concluding remarks

We believe that our locally periodic inverse-design approach represents a powerful extension to the ideas in previous work, allowing one to balance competing tradeoffs in wavefront design, optimize arbitrary functions of the scattered field (e.g. intensity in selected regions), evaluate parameter sensitivity, design for robustness to uncertainties, and to go beyond the regime of subwavelength structures and far-field designs. A similar max–min formulation can be used to implement a standard *robust* optimization method to account for manufacturing uncertainty [56, 62]. Our approximate solver remains orders of magnitude faster than optimization methods based on full Maxwell solvers, allowing it to scale to aperiodic structures hundreds or thousands of wavelengths in diameter while retaining acceptable accuracy for typical designs. We find that complex behaviors can be designed even from very simple unit cells without plasmonic resonances, and without operating in deeply subwavelength regimes.

This paper presented a proof of concept and validation of the approach, and opens up many future possibilities. We are currently working on extension to the design of 3d surfaces and vector fields, and believe such problems to be tractable with a few hours of computation (rather than the few minutes required here for 2d inverse problems). We can easily extend our inverse design from a single parameter per unit cell to multiple parameters per cell. With a few (≲ 10) parameters, one can use a similar library-based approach via multidimensional interpolation, for which the main limitation is the number *N* of unit-cell calculations that need to be solved beforehand in order to build the interpolation library. The simplest method is a tensor product of Chebyshev polynomials [40], which is practical for at most 2–3 parameters because *N* grows exponentially with the number of parameters. Polynomial scaling of *N* can be achieved by sparse-grid methods [41] *or neural networks* [63, 64]. To handle hundreds or thousands of parameters per unit cell for topology optimization [6–10,21,65], the library approach must be abandoned in favor of directly solving Maxwell’s equations in every metasurface unit cell for each optimization iteration (still via the locally periodic approximation). In this case, the cost is essentially independent of the number of parameters and scales linearly with the number of unit cells, which can be solved in parallel; we have successfully optimized metasurfaces with > 1000 parameters per unit cell in this way and are currently preparing a manuscript on those results. Multiple parameters per unit cell could describe more complicated surface patterns (e.g. the V-shaped antennas of [4]), but also includes the possibility of *multi-layer* patterns (e.g. stacked gratings). Additional degrees of freedom could prove crucial for obtaining truly wide-bandwidth devices, coupling multiple polarizations, minimizing unwanted reflections, and so on.

The ability to design non-subwavelength surface patterns (but still far from the ≫ *λ* regime of scalar diffraction theory [33, 34]) could prove useful for a variety of applications, starting with designs for short wavelengths (e.g. near UV) where subwavelength fabrication is difficult. The additional diffracted orders of large-period structures may also become useful for near-field focusing and related design problems or for focusing a single incident beam at multiple spots.

Another interesting direction to explore would be further development of the theory of nearly periodic structures and locally periodic approximations. In a companion work [35], we develop a rigorous theory of slowly varying (nearly *uniform*) structures, and show that the analogous “locally uniform” approximation appears as the 0th-order term in a convergent series of integral corrections. A corresponding rigorous theory of higher-order corrections to the locally periodic approximation, analogous to coupled-mode expansions for propagation *through* nearly periodic media [61], along with efficient numerical methods to obtain corrections, is an important goal for the theory of metasurfaces. A closely related problem is coupling radiation to and from guided modes by nearly-periodic surfaces, a version of which is solved in [35]. In another paper [35], we have recently shown that similar local approximations can indeed be used to compute both near fields and coupling to guided waves.

## Appendix

To use standard high-dimensional optimization algorithms, one needs to provide an efficient computation of both the objective (cost) function and its gradient. There is a well-known technique called an *adjoint method* [59] that can be used to efficiently compute the gradient for any number of parameters with a cost comparable to evaluating the objective function at most *twice*, which is commonly used in topology optimization [6–10,21,65]. In the case of the two objectives presented in Sec. 3.1 and Sec. 3.2, the gradient is especially simple to evaluate as described in this Appendix.

In Sec. 3.1, $f|(p,{s}_{0},{\varphi}_{0})={\displaystyle \int {\left|s(p(x))\u2013{s}_{0}a(x){e}^{i\varphi (x)+i{\varphi}_{0}}\right|}^{2}dx}$, so the gradient is:

*∂f*/

*∂p*denotes the functional derivative [66] with respect to the parameter function

*p x*and ∗ denotes complex conjugation. Notice that the computation of the gradient requires only the evaluation of a few simple integrals, comparable to the cost of evaluating

*f*. Similarly, in Sec. 3.2, $g(p,\mathbf{x})={\left|{E}_{z}(\mathbf{x})\right|}^{2}={\left|{\displaystyle {\int}_{\mathrm{y}={\mathrm{y}}_{0}}G(\mathbf{x},({x}^{\prime},0))s(p({x}^{\prime}))d{x}^{\prime}}\right|}^{2}$, and so its gradient is:

## Funding

U. S. Army Research Office through the Institute for Soldier Nanotechnologies (W911NF-13-D-0001).

## References

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

**2. **N. Yu, P. Genevet, M. A. Kats, F. Aieta, J.-P. Tetienne, F. Capasso, and Z. Gaburro, “Light propagation with phase discontinuities: Generalized laws of reflection and refraction,” Science **334**, 333–337 (2011). [CrossRef] [PubMed]

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

**4. **N. Yu, P. Genevet, F. Aieta, M. A. Kats, R. Blanchard, G. Aoust, J.-P. Tetienne, Z. Gaburro, and F. Capasso, “Flat optics: Controlling wavefronts with optical antenna metasurfaces,” IEEE J. Sel. Top. Quantum Electron. **19**, 4700423 (2013). [CrossRef]

**5. **S. Jahani and Z. Jacob, “All-dielectric metamaterials,” Nat. Nanotech. **11**, 23–36 (2016). [CrossRef]

**6. **C. M. Lalau-Keraly, S. Bhargava, O. D. Miller, and E. Yablonovitch, “Adjoint shape optimization applied to electromagnetic design,” Opt. express **21**, 21693–21701 (2013). [CrossRef] [PubMed]

**7. **A. Y. Piggott, J. Lu, K. G. Lagoudakis, J. Petykiewicz, T. M. Babinec, and J. Vučković, “Inverse design and demonstration of a compact and broadband on-chip wavelength demultiplexer,” Nat. Photon. **9**, 374–377 (2015). [CrossRef]

**8. **A. Y. Piggott, J. Petykiewicz, L. Su, and J. Vučković, “Fabrication-constrained nanophotonic inverse design,” Sci. Rep. **7**, 1786 (2017). [CrossRef] [PubMed]

**9. **D. Sell, J. Yang, S. Doshay, and J. A. Fan, “Periodic dielectric metasurfaces with high-efficiency, multiwavelength functionalities,” Adv. Opt. Mater. **5**, 1700645 (2017). [CrossRef]

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

**11. **K. H. Matlack, M. Serra-Garcia, A. Palermo, S. D. Huber, and C. Daraio, “Designing perturbative metamaterials from discrete models,” Nat. Mater. **17**, 323 (2018). [CrossRef] [PubMed]

**12. **F. Aieta, P. Genevet, M. A. Kats, N. Yu, R. Blanchard, Z. Gaburro, and F. Capasso, “Aberration-free ultrathin flat lenses and axicons at telecom wavelengths based on plasmonic metasurfaces,” Nano Lett. **12**, 4932–4936 (2012). [CrossRef] [PubMed]

**13. **A. Arbabi, E. Arbabi, Y. Horie, S. M. Kamali, and A. Faraon, “Planar metasurface retroreflector,” Nat. Photon. **11**, 415–420 (2017). [CrossRef]

**14. **M. Khorasaninejad, W. T. Chen, A. Y. Zhu, J. Oh, R. C. Devlin, C. Roques-Carmes, I. Mishra, and F. Capasso, “Visible wavelength planar metalenses based on titanium dioxide,” IEEE J. Sel. Top Quantum Electron. **23**, 43–58 (2017). [CrossRef]

**15. **F. Aieta, M. A. Kats, P. Genevet, and F. Capasso, “Multiwavelength achromatic metasurfaces by dispersive phase compensation,” Science **347**, 1342–1345 (2015). [CrossRef] [PubMed]

**16. **M. Khorasaninejad, F. Aieta, P. Kanhaiya, M. A. Kats, P. Genevet, D. Rousso, and F. Capasso, “Achromatic metasurface lens at telecommunication wavelengths,” Nano Lett. **15**, 5358–5362 (2015). [CrossRef] [PubMed]

**17. **M. Khorasaninejad, Z. Shi, A. Y. Zhu, W. T. Chen, V. Sanjeev, A. Zaidi, and F. Capasso, “Achromatic metalens over 60 nm bandwidth in the visible and metalens with reverse chromatic dispersion,” Nano Lett. **17**, 1819–1824 (2017). [CrossRef] [PubMed]

**18. **E. Arbabi, A. Arbabi, S. M. Kamali, Y. Horie, and A. Faraon, “Controlling the sign of chromatic dispersion in diffractive optics with dielectric metasurfaces,” Optica **4**, 625 (2017). [CrossRef]

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

**20. **B. Groever, C. Roques-Carmes, S. J. Byrnes, and F. Capasso, “Substrate aberration and correction for meta-lens imaging: an analytical approach,” Appl. Opt. **57**, 2973–2980 (2018). [CrossRef] [PubMed]

**21. **Z. Lin, B. Groever, F. Capasso, A. W. Rodriguez, and M. Lončar, “Topology-optimized multilayered metaoptics,” Phys. Rev. Appl. **9**, 044030 (2018). [CrossRef]

**22. **K. Achouri, M. A. Salem, and C. Caloz, “General metasurface synthesis based on susceptibility tensors,” IEEE Trans. Antennas Propag. **63**, 2977–2991 (2015). [CrossRef]

**23. **A. Epstein and G. V. Eleftheriades, “Floquet-Bloch analysis of refracting Huygens metasurfaces,” Phys. Rev. B **90**, 235127 (2014). [CrossRef]

**24. **A. Epstein and G. V. Eleftheriades, “Passive lossless Huygens metasurfaces for conversion of arbitrary source field to directive radiation,” IEEE Trans. Antennas Propag. **62**, 5680–5695 (2014). [CrossRef]

**25. **C. L. Holloway, A. Dienstfrey, E. F. Kuester, J. F. O’Hara, A. K. Azad, and A. J. Taylor, “A discussion on the interpretation and characterization of metafilms/metasurfaces: The two-dimensional equivalent of metamaterials,” Metamaterials **3**, 100–112 (2009). [CrossRef]

**26. **C. L. Holloway, E. F. Kuester, J. A. Gordon, J. O’Hara, J. Booth, and D. R. Smith, “An overview of the theory and applications of metasurfaces: The two-dimensional equivalents of metamaterials,” IEEE Antennas Propag. Mag. **54**, 10–35 (2012). [CrossRef]

**27. **E. F. Kuester, M. A. Mohamed, M. Piket-May, and C. L. Holloway, “Averaged transition conditions for electromagnetic fields at a metafilm,” IEEE Trans. Antennas Propag. **51**, 2641–2651 (2003). [CrossRef]

**28. **C. Pfeiffer and A. Grbic, “Metamaterial Huygens'surfaces: Tailoring wave fronts with reflectionless sheets,” Phys. Rev. Lett. **110**, 197401 (2013). [CrossRef]

**29. **S. Tcvetkova, D.-H. Kwon, A. Díaz-Rubio, and S. Tretyakov, “Near-perfect conversion of a propagating plane wave into a surface wave using metasurfaces,” Phys. Rev. B **97**, 115447 (2018). [CrossRef]

**30. **S. Tretyakov, “Metasurfaces for general transformations of electromagnetic fields,” Phil. Trans. R. Soc. A **373**, 20140362 (2015). [CrossRef] [PubMed]

**31. **L. Verslegers, P. B. Catrysse, Z. Yu, W. Shin, Z. Ruan, and S. Fan, “Phase front design with metallic pillar arrays,” Opt. Lett. **35**, 844–846 (2010). [CrossRef] [PubMed]

**32. **J. Cheng, S. Inampudi, and H. Mosallaei, “Optimization-based dielectric metasurfaces for angle-selective multifunctional beam deflection,” Sci. Rep. **7**, 12228 (2017). [CrossRef] [PubMed]

**33. **M. Born and E. Wolf, *Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light* (Elsevier, 2013).

**34. **D. C. O’Shea, T. J. Suleski, A. D. Kathman, and D. W. Prather, *Diffractive Optics: Design, Fabrication, and Test*, vol. 62 (Spie, 2004).

**35. **C. Pérez-Arancibia, R. Pestourie, and S. G. Johnson, “Sideways adiabaticity: Beyond ray optics for slowly varying metasurfaces,” arXiv:1806.05330 (2018).

**36. **M. Kim and G. V. Eleftheriades, “Design and demonstration of impedance-matched dual-band chiral metasurfaces,” Sci. Rep. **8**, 3449 (2018). [CrossRef] [PubMed]

**37. **J. D. Joannopoulos, S. G. Johnson, J. N. Winn, and R. D. Meade, *Photonic Crystals: Molding the Flow of Light* (Princeton University, 2008), 2nd ed.

**38. **M. Costabel and F. Le Louër, “Shape derivatives of boundary integral operators in electromagnetic scattering. part i: Shape differentiability of pseudo-homogeneous boundary integral operators,” IEOT **72**, 509–535 (2012).

**39. **M. Costabel and F. Le Louër, “Shape derivatives of boundary integral operators in electromagnetic scattering. part ii: Application to scattering by a homogeneous dielectric obstacle,” IEOT **73**, 17–48 (2012).

**40. **J. P. Boyd, *Chebyshev and Fourier Spectral Methods.* (Dover Publications, 2001).

**41. **T. Gerstner and M. Griebel, “Numerical integration using sparse grids,” Numer. Algorithms **18**, 209 (1998). [CrossRef]

**42. **N. J. Champagne II, J. G. Berryman, and H. M. Buettner, “FDFD: A 3D finite-difference frequency-domain code for electromagnetic induction tomography,” J. Comput. Phys. **170**, 830–848 (2001). [CrossRef]

**43. **W. Shin, “MaxwellFDFD Webpage,” (2015). https://github.com/wsshin/maxwellfdfd.

**44. **W. Shin and S. Fan, “Choice of the perfectly matched layer boundary condition for frequency-domain Maxwell's equations solvers,” J. Comput. Phys. **231**, 3406–3431 (2012). [CrossRef]

**45. **R. Mittra and U. Pekel, “A new look at the perfectly matched layer (PML) concept for the reflectionless absorption of electromagnetic waves,” IEEE Microw. Wirel. Compon. Lett. **5**, 84–86 (1995).

**46. **R. F. Harrington, *Time-Harmonic Electromagnetic Fields* (Wiley-IEEE, 2001), 2nd ed. [CrossRef]

**47. **A. Oskooi and S. G. Johnson, “Electromagnetic wave source conditions,” arXiv preprint arXiv:1301.5366 (2013).

**48. **A. Taflove and S. C. Hagness, *Computational electrodynamics: the finite-difference time-domain method* (Artech house, 2005).

**49. **J. D. Jackson, *Classical Electrodynamics* (Wiley, 1999), 3rd ed.

**50. **K. Watanabe, *Integral Transform Techniques for Green’s Function*, vol. 71 of Lecture Notes in Applied and Computational Mechanics (Springer International Publishing, 2014). DOI: . [CrossRef]

**51. **O. P. Bruno, M. Lyon, C. Pérez-Arancibia, and C. Turc, “Windowed Green function method for layered-media scattering,” SIAM J. Appl. Math. **76**, 1871–1898 (2016). [CrossRef]

**52. **O. P. Bruno, E. Garza, and C. Pérez-Arancibia, “Windowed Green function method for nonuniform open-waveguide problems,” IEEE Trans. Antennas Propag. **65**, 4684–4692 (2017). [CrossRef]

**53. **O. P. Bruno and C. Pérez-Arancibia, “Windowed Green function method for the helmholtz equation in the presence of multiply layered media,” Proc. R. Soc. A **473**, 20170161 (2017). [CrossRef] [PubMed]

**54. **K. Svanberg, “A class of globally convergent optimization methods based on conservative convex separable approximations,” SIAM J. Optim. **12**, 555–573 (2002). [CrossRef]

**55. **S. G. Johnson, “The nlopt nonlinear-optimization package,” http://ab-initio.mit.edu/nlopt.

**56. **A. Mutapcic, S. Boyd, A. Farjadpour, S. G. Johnson, and Y. Avniel, “Robust design of slow-light tapers in periodic waveguides,” Eng. Optim. **41**, 365–384 (2009). [CrossRef]

**57. **T. F. Chan, J. Cong, T. Kong, and J. R. Shinnerl, “Multilevel optimization for large-scale circuit placement,” in Proceedings of the 2000 IEEE/ACM International Conference on Computer-Aided Design, (IEEE, 2000), pp. 171–176.

**58. **K. Chun and J. B. Ra, “Fast block-matching algorithm by successive refinement of matching criterion,” in *Visual Communications and Image Processing’92*, vol. 1818 (International Society for Optics and Photonics, 1992), pp. 552–561.

**59. **G. Strang, *Computational Science and Engineering*, vol. 791 (Wellesley-Cambridge, 2007).

**60. **S. Boyd and L. Vandenberghe, *Convex Optimization* (Cambridge University, 2004). [CrossRef]

**61. **S. G. Johnson, P. Bienstman, M. A. Skorobogatiy, M. Ibanescu, E. Lidorikis, and J. D. Joannopoulos, “Adiabatic theorem and continuous coupled-mode theory for efficient taper transitions in photonic crystals,” Phys. Rev. E **66**, 066608 (2002). [CrossRef]

**62. **M. Sniedovich, “From statistical decision theory to robust optimization: A maximin perspective on robust decision-making,” in *Robustness Analysis in Decision Aiding, Optimization, and Analytics*, (Springer, 2016), pp. 59–87.

**63. **D. Liu, Y. Tan, E. Khoram, and Z. Yu, “Training deep neural networks for the inverse design of nanophotonic structures,” ACS Photonics **5**, 1365–1369 (2018). [CrossRef]

**64. **J. Peurifoy, Y. Shen, L. Jing, Y. Yang, F. Cano-Renteria, B. G. DeLacy, J. D. Joannopoulos, M. Tegmark, and M. Soljačić, “Nanophotonic particle simulation and inverse design using artificial neural networks,” Sci. Adv. **4**, eaar4206 (2018). [CrossRef] [PubMed]

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

**66. **R. Courant and J. Moser, *Calculus of variations and supplementary notes and exercises, 1945-1946* (New York University, 1957).