Inverse design of large-area metasurfaces

We present a computational framework for efficient optimization-based"inverse design"of large-area"metasurfaces"(subwavelength-patterned surfaces) for applications such as multi-wavelength and 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).


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 [34,59,32,58,26], incorporating both scattered amplitude and phase for multiple incident λ and angles. Previous methods either optimized the full Maxwell equations [35,44,45,46,60] (which is infeasible for large surfaces), were restricted to weakly coupled scatterers [37], or started with a desired scattered phase and tried to design a corresponding metasurface unit cell [2,4,30,3,29,31,5,51,21,36] (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 [2,4,30,3,29,31,5,36,51,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 [1,18,19,23,24,33,43,53,54] 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). Figure 1: Schematic of our design method: exact Maxwell scattering solutions for a set of periodic unit cells (top) are composed into an approximate solution for an arbitrary aperiodic composition (right), and this approximation is then used for large-scale optimization to determine the metasurface parameters to maximize a given objective (e.g. the focal intensity, bottom). The performance of the final design can then be fed back into adjusting the design of the unit cell (left).

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 [2,4,30,3,29,31,5,51,21], 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" [55,14]. (Contrast this with the regime of scalar diffraction theory [6,40], 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 [42], we develop the rigorous foundations and convergence rates of a related approximation, along with higher-order corrections, but here (similar to previous authors [2,4,30,3,29,31,5,51]) 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 [2,4,30,3,29,31,5,51], 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 [31]: 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 [42] 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 [42], 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.)

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 Transmitted field 40

80
Half width (nm)  Figure 2: Left: an arbitrary aperiodic metasurface (top) is approximated by solving a set of periodic scattering problems (bottom), one for each unit cell, to obtain the scattered field just above the surface (horizontal line segments). Right: 0th diffracted-order amplitude (top) and phase (bottom) of periodic subproblems as a function of the pillar width. This is precomputed for several widths (markers) and interpolated as needed.
diffractive order is the only propagating wave [27]. 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 [16,17] 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 approximationf (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 [7] or by more sophisticated methods such as sparse grids [20]. In this way, we only need to solve the unit-cell Maxwell problem a few times to obtain our polynomial approximantf (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 [7]), and can then interpolate to high accuracy using Chebyshev polynomials [7]. Multiple parameters per unit cell could be interpolated using a product of Chebyshev polynomials [7] or, for more than 3-4 parameters, sparsegrid interpolation [20]. Here, we use the finite-difference frequency-domain (FDFD) method [12,47] with perfectly matched layer (PML) absorbing boundaries [48,38], but any other computational method for periodic Maxwell problems would work as well.

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 [22,41]: the fields in the 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 [41]: wheren is the surface unit-normal vector and the delta function implies that these are surface currents.
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 [22,41], 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 [25]). This allows us to use onlyn × 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 ) [22]. For our twodimensional model problem (xy plane) with the E z polarization, where we only have a current 1 (kr)n · r r , where k = 2π λ , 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 [10,9,11].
This equivalent-currents formulation is exact if the true aperiodic E z field is used for the K x 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 uniform approximate scattered fields can be thought of as sources producing fields everywhere else [6,40], further approximated in the far field by e.g. the Fraunhofer diffraction theory [6].)

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. We use a standard gradient-based optimization method [52]. To avoid getting trapped in poor local optima, we use a technique called successive refinement [39,13,15]: we successively double the number of degrees of freedom, using the optimized coarser structures as starting points for optimization of the finer structures. What objective function should we optimize? In Sec. 3.1 we consider an objective function similar to previous work [3,36], 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.

Optimizing the wavefront
When the exact desired field is known everywhere above the metasurface, as in lens design [2,4,30,3,29,31,5,21,51,36] and other wavefront shaping problems, by the equivalence principle [22] 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) : where 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) [3]. A closely related approach was used for metalens design in several previous works [2,4,30,3,29,31,5,51,36]. 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 [2], but an optimization-based approach is more flexible at balancing tradeoffs in cases where the desired ae iφ cannot be exactly obtained, especially in multi-frequency problems (Sec. 4.1). A phase-based optimization approach was directly employed in [36] for topology optimization of a small area (no locally periodic approximation). 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 [3]. We optimize over piecewise-constant parameters, effectively one parameter per unit cell, with a standard optimization algorithm [52] 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 surfaceintegral equation (SIE) Maxwell solver [10,9,11] 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.

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 [52] with an analytically computed gradient (essentially via an adjoint method [50]), 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.

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: max parameters min λ∈wavelengths objective(parameters, λ) .
(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 [8] max t,parameters t subject to t ≤ objective(parameters, λ) for λ ∈ wavelengths.
where t ∈ R is a new "dummy" optimization parameter. Assuming that the objective function is differentiable, we can now use a standard nonlinear constrained-optimization algorithm [52].
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).

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.

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.

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.

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 [36]. 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. λ=470 nm λ=530 nm λ=650 nm Figure 6: Bottom: focal lines for the three target wavelengths (blue, green and red) focus on points sixty microns apart. Top: the field produced by our design focuses on the desired foci, the high-intensity regions for blue (left) and red (right) are tilted because their foci are off-axis. 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 [36], 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.

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 [1,18,19,23,24,33,43,53,54], and most previous work operated in a subwavelength regime [2,4,30,3,29,31,5,51]. Conversely, when the period is larger than the wavelength, additional diffracted waves appear in the far field [27] 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 [28]) 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 [27].
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  Figure 7: Left: the focal lines for 0-degrees, 3-degrees, 6-degrees, and 9-degrees angles of incidence show four foci with the maximum intensity at the target focal spot (at x = 0), the other three peaks correspond to the other three target angles. Right: corresponding produced field around the foci, the focal spot becomes more tilted as the angle of incidence increases.
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 diffractionlimited) 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. Figure 8: Bottom: the geometry (left) from intensity optimization shows big variations in the width of the pillar, and produce good focusing when simulated with a brute force simulation (right), or our locally periodic solver (middle) which includes the diffractive orders ±1. Top: the geometry (left) from wavefront optimization shows poor focusing both using our locally periodic solver (middle) or a brute force calculation (right). All the intensity plots have the same color scale.

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 [39,49]. 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, either continuing with a library-based approach using multidimensional interpolation (perhaps based on sparse-grid methods [20] for more than 2-3 parameters), or switching to a more flexible approach in which the periodic calculations are repeated for each unit cell during optimization (which opens the possibility of large-scale topology optimization [57,36,35] of the unit cells). Multiple parameters per unit cell could describe more complicated surface patterns (e.g. the V-shaped antennas of [58]), 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 [6,40]) 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 [42], 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 [28], 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 nearlyperiodic surfaces, a version of which is solved in [42].

Funding
This work was supported in part by the U. S. Army Research Office through the Institute for Soldier Nanotechnologies at MIT under contract number W911NF-13-D-0001.