## Abstract

We demonstrate optimization of optical metasurfaces over 10^{5}–10^{6} degrees of freedom in two and three dimensions, 100–1000+ wavelengths (*λ*) in diameter, with 100+ parameters per *λ*^{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 *λ*^{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.

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

## 1. Introduction

We present a method for large-area (≥ 100 wavelengths *λ*) topology optimization (TO) [1–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–13], TO allows us to consider thousands of parameters per unit cell, with much larger (∼ 5*λ*) 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 ≲ 10*λ* amenable to brute-force simulation [14,15], we exploit a locally periodic approximation [6–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 ∼ 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–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–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–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–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–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 ≲ 10*λ*. 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 **r _{0}** [16]:

*∊̄*(

**r**)} is related to the position-dependent dielectric profile via

*∊*(

**r**) = (

*∊*

_{st}−

*∊*

_{bg})

*∊̄*(

**r**) +

*∊*

_{bg}, where

*∊*

_{st (bg)}denotes the relative permittivity of the structural (background) dielectric material. Whereas

*∊̄*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.,

*∊̄*

_{optimal}∈ {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

**G**with the near-field

_{0}**E**at some reference plane

**r**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,

_{s}**E**is the steady-state solution to the frequency-domain Maxwell’s equation:

**J**at a frequency

*ω*where we set the magnetic permeability

*μ*= 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

**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 $\frac{\partial f}{\partial \overline{\u220a}}$ at every pixel point in the device region, which we obtain by an adjoint method [3,30]:

**Ẽ**is evaluated by solving an additional Maxwell’s equation:

#### 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 (≲ 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 {*ω _{i}*,

**J**

*}. In such cases, the optimization problem can be easily extended by the maximin formulation:*

_{i}#### 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 {*θ* (°) = 3*i* |*i* = 0, 1, ..., 10 the same focus. The lens consists of five layers of TiO_{2} (refractive index *n* ≈ 2.4), four of which are buried in silica (*n* ≈ 1.5). The diameter of the lens is *D* = 1200*λ*, where *λ* is the operational wavelength, and the focal length is *F* = 1000*λ*, 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*λ* long, for 3 × 10^{5} DoFs in total, while the entire computation is parallelized over 1200 CPUs. While we perform the optimization under the locally periodic approximation, we validate the optimized result with rigorous full-wave FDTD simulations [42] with a high spatial resolution of *λ*/50 (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 ∼ *λ* while the average transmission efficiency is found to be *T* ≈ 40% and an average focusing efficiency of ∼ 55% (we define the focusing efficiency as the integrated field intensity within a radius of 2 FWHMs around the focal spot divided by the total intensity over the image plane [43]). The capability to focus at discrete angles (over a finite angular bandwidth) is consistent with reciprocity and brightness theorems [44].

#### 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 {*λ*_{1}, *λ*_{2}, *λ*_{3}} with *λ*_{2} = *λ*_{1}/1.2 and *λ*_{3} = *λ*_{1}/1.38. In particular, one may choose *λ*_{1} = 650 nm, *λ*_{2} = 540 nm and *λ*_{3} = 470 nm, corresponding to red, green and blue (RGB) wavelengths. The lens consists of two layers of TiO_{2} (*n* ≈ 2.4), one of which is buried in silica (*n* ≈ 1.5). The diameter of the lens is *D* = 200*λ*_{1} and the focal length is *F* = 100*λ*_{1}, corresponding to NA = 0.71. During the optimization, the entire lens is broken up into 40 cells, each of which is 5*λ* long, for 1.6 × 10^{4} DoFs in total, while the entire computation is parallelized over just 20 CPUs. redDuring optimization, we employ RCWA with 100 Fourier orders to model each unit cell; at the end, a rigorous FDTD simulation of the entire optimized lens is performed with a high resolution of *λ*_{1}/80. 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 ≈ 0.71*λ*_{1}, 0.72*λ*_{2} and 0.76*λ*_{3} respectively and with an average transmission efficiency of ∼ 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–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 are found to be ∼ 51%, 50% and 45% while LPA predicts ∼ 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], by considering overlapping domains instead of sharply truncated boundaries at the unit cell level, 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:

*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}*
×

*L*=

_{y}*λ*×

*λ*. The cells are repeated along

*y*but are allowed to differ from each other along

*x*. Figure 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 × 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 *λ* × *λ* cells, corresponding to NA=0.37, with 6.4 × 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–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–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].

## Funding

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

## Acknowledgements

We thank Owen Miller for helpful discussions.

## References

**1. **J. S. Jensen and O. Sigmund, “Topology optimization
for nano-photonics,” Laser Photonics
Rev. **5**,
308–321
(2011). [CrossRef]

**2. **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]

**3. **S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, “Inverse design in nanophotonics,” Nat. Photonics **12**, 659 (2018). [CrossRef]

**4. **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).

**5. **N. Yu and F. Capasso, “Flat optics with designer metasurfaces,” Nat. Mater. **13**, 139 (2014). [CrossRef] [PubMed]

**6. **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]

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

**8. **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]

**9. **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]

**10. **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]

**11. **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–632 (2017). [CrossRef]

**12. **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]

**13. **M. Khorasaninejad, W. T. Chen, R. C. Devlin, J. Oh, A. Y. Zhu, and F. Capasso, “Metalenses at visible wavelengths: Diffraction-limited focusing and subwavelength resolution imaging,” Science **352**, 1190–1194 (2016). [CrossRef] [PubMed]

**14. **D. Sell, J. Yang, S. Doshay, R. Yang, and J. A. Fan, “Large-angle,
multifunctional metagratings based on freeform multimode
geometries,” Nano Lett. **17**, 3752–3757
(2017). [CrossRef] [PubMed]

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

**16. **R. Pestourie, C. Pérez-Arancibia, Z. Lin, W. Shin, F. Capasso, and S. G. Johnson, “Inverse design of
large-area metasurfaces,” Opt.
Express **26**,
33732–33747
(2018). [CrossRef]

**17. **M. Moharam and T. Gaylord, “Rigorous coupled-wave
analysis of planar-grating diffraction,”
J. Opt. Soc. Am. **71**,
811–818
(1981). [CrossRef]

**18. **V. Liu and S. Fan, “S^{4} : A free electromagnetic solver for layered periodic structures,” Comput. Phys. Commun. **183**, 2233–2244 (2012). [CrossRef]

**19. **C. Pérez-Arancibia, R. Pestourie, and S. G. Johnson, “Sideways adiabaticity: beyond ray optics for slowly varying metasurfaces,” Opt. Express **26**, 30202–30230 (2018). [CrossRef] [PubMed]

**20. **F. Capasso, “The future and promise of flat optics: a personal perspective,” Nanophotonics **7**, 953–957 (2018). [CrossRef]

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

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

**23. **G. Kim, J. A. Domínguez-Caballero, and R. Menon, “Design and analysis of
multi-wavelength diffractive optics,”
Opt. Express **20**,
2814–2823
(2012). [CrossRef] [PubMed]

**24. **W. T. Chen, A. Y. Zhu, V. Sanjeev, M. Khorasaninejad, Z. Shi, E. Lee, and F. Capasso, “A broadband achromatic metalens for focusing and imaging in the visible,” Nat. Nanotechnol. **13**, 220 (2018). [CrossRef] [PubMed]

**25. **S. Wang, P. C. Wu, V.-C. Su, Y.-C. Lai, M.-K. Chen, H. Y. Kuo, B. H. Chen, Y. H. Chen, T.-T. Huang, J.-H. Wang, R.-M. Lin, C.-H. Kuan, T. Li, Z. Wang, S. Zhu, and D. P. Tsai, “A broadband achromatic metalens in the visible,” Nat. Nanotechnol. **13**, 227 (2018). [CrossRef] [PubMed]

**26. **S. Shrestha, A. C. Overvig, M. Lu, A. Stein, and N. Yu, “Broadband achromatic
dielectric metalenses,” Light. Sci.
Appl. **7**, 85
(2018). [CrossRef]

**27. **S. M. Kamali, E. Arbabi, A. Arbabi, Y. Horie, M. Faraji-Dana, and A. Faraon, “Angle-multiplexed metasurfaces: Encoding independent wavefronts in a single metasurface under different illumination angles,” Phys. Rev. X **7**, 041056 (2017).

**28. **M. Mitchell, *An Introduction to Genetic Algorithms* (MIT, 1998).

**29. **J. Kennedy, “Particle swarm optimization,” in *Encyclopedia of Machine Learning*, (Springer, 2011), pp. 760–766.

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

**31. **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. Photonics **9**, 374–377 (2015). [CrossRef]

**32. **M. Zhou, B. S. Lazarov, F. Wang, and O. Sigmund, “Minimum length scale in topology optimization by geometric constraints,” Comput. Methods Appl. Mech. Eng. **293**, 266–282 (2015). [CrossRef]

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

**34. **A. Taflove, A. Oskooi, and S. G. Johnson, “Electromagnetic wave
source conditions,” in *Advances in
FDTD Computational Electrodynamics: Photonics and
Nanotechnology*, (Artech
House, 2013), pp.
65–100.

**35. **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]

**36. **W. D. Gropp, W. Gropp, E. Lusk, A. Skjellum, and A. D. F. E. E. Lusk, *Using MPI: Portable Parallel Programming with the Message-Passing Interface*, vol. 1 (MIT, 1999).

**37. **P. Pacheco, *An Introduction to Parallel Programming* (Elsevier, 2011).

**38. **K. Svanberg, “The method of moving
asymptotes — a new method for structural
optimization,” Int. J. Numer. Meth.
Engng. **24**,
359–373
(1987). [CrossRef]

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

**40. **S. G. Johnson, “The NLOpt nonlinear-optimization package,” http://github.com/stevengj/nlopt (2014).

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

**42. **A. F. Oskooi, D. Roundy, M. Ibanescu, P. Bermel, J. D. Joannopoulos, and S. G. Johnson, “MEEP: A flexible free-software package for electromagnetic simulations by the FDTD method,” Comput. Phys. Commun. **181**, 687–702 (2010). [CrossRef]

**43. **M. Khorasaninejad, A. Y. Zhu, C. Roques-Carmes, W. T. Chen, J. Oh, I. Mishra, R. C. Devlin, and F. Capasso,
“Polarization-insensitive metalenses at visible
wavelengths,” Nano Lett. **16**, 7229–7234
(2016). [CrossRef] [PubMed]

**44. **H. Zhang, C. W. Hsu, and O. D. Miller, “Scattering
concentration bounds: Brightness theorems for wave
scattering,” arXiv 1810.02727
(2018).

**45. **C. V. Poulton, A. Yaacobi, D. B. Cole, M. J. Byrd, M. Raval, D. Vermeulen, and M. R. Watts, “Coherent solid-state LIDAR with silicon photonic optical phased arrays,” Opt. Lett. **42**, 4091–4094 (2017). [CrossRef] [PubMed]

**46. **A. Pick, B. Zhen, O. D. Miller, C. W. Hsu, F. Hernandez, A. W. Rodriguez, M. Soljačić, and S. G. Johnson, “General theory of
spontaneous emission near exceptional points,”
Opt. Express **25**,
12325–12348
(2017). [CrossRef] [PubMed]

**47. **Z. Lin, A. Pick, M. Lončar, and A. W. Rodriguez, “Enhanced spontaneous
emission at third-order dirac exceptional points in inverse-designed
photonic crystals,” Phys. Rev.
Lett. **117**, 107402
(2016). [CrossRef] [PubMed]

**48. **N. Rivera, I. Kaminer, B. Zhen, J. D. Joannopoulos, and M. Soljačić, “Shrinking light to allow forbidden transitions on the atomic scale,” Science **353**, 263–269 (2016). [CrossRef] [PubMed]

**49. **Z. Lin, X. Liang, M. Lončar, S. G. Johnson, and A. W. Rodriguez, “Cavity-enhanced second-harmonic generation via nonlinear-overlap optimization,” Optica **3**, 233–238 (2016). [CrossRef]

**50. **P. S. Venkataram, J. Hermann, A. Tkatchenko, and A. W. Rodriguez, “Unifying microscopic
and continuum treatments of van der Waals and Casimir
interactions,” Phys. Rev.
Lett. **118**, 266802
(2017). [CrossRef] [PubMed]

**51. **A. W. Rodriguez, O. Ilic, P. Bermel, I. Celanovic, J. D. Joannopoulos, M. Soljačić, and S. G. Johnson, “Frequency-selective
near-field radiative heat transfer between photonic crystal slabs: a
computational approach for arbitrary geometries and
materials,” Phys. Rev. Lett. **107**, 114302
(2011). [CrossRef] [PubMed]

**52. **M. Reid, O. Miller, A. Polimeridis, A. Rodriguez, E. Tomlinson, and S. Johnson, “Photon torpedoes and
Rytov pinwheels: Integral-equation modeling of non-equilibrium
fluctuation-induced forces and torques on
nanoparticles,” arXiv 1708.01985
(2017).

**53. **Y. E. Lee, O. D. Miller, M. H. Reid, S. G. Johnson, and N. X. Fang, “Computational inverse
design of non-intuitive illumination patterns to maximize optical
force or torque,” Opt.
Express **25**,
6757–6766
(2017). [CrossRef] [PubMed]

**54. **T. F. Chan and T. P. Mathew, “Domain decomposition algorithms,” Acta Numer. **3**, 61–143 (1994). [CrossRef]

**55. **N. Zhao, S. Verweij, W. Shin, and S. Fan, “Accelerating
convergence of an iterative solution of finite difference frequency
domain problems via schur complement domain
decomposition,” Opt. Express **26**, 16925–16939
(2018). [CrossRef] [PubMed]

**56. **S. Tao, J. Cheng, and H. Mosallaei, “An integral equation based domain decomposition method for solving large-size substrate-supported aperiodic plasmonic array platforms,” MRS Commun. **6**, 105–115 (2016). [CrossRef]