the unique morphologies of these structures and enhance their efficiency by precomputing light scattering by the individual scatterers that compose the metasurface.

In the widely used T-matrix method, the field outside the scatterer is expanded using spherical vectorial wave harmonics and a T-matrix is utilized to describe how each scatterer transforms incident into scattered multipolar fields. The method has been extensively explored over decades⁹¹¹ and is now available through multiple open-source implementations.⁹¹¹,⁹¹² Alternative methods account for multiple scattering using vectorial Green functions and treat scatterers as ensembles of equivalent electric and magnetic surface currents,⁹¹³,⁹¹⁴ or employ numerical non-local dipoles⁹¹⁵ obtained by numerically solving an inverse scattering problem. A detailed comparative analysis of these methods strengths and limitations is beyond this roadmap's scope and is deferred to a comprehensive future review.⁹¹⁶

Two approaches exist for the extrapolation step. In the field stitching approach, Maxwell equations are solved for numerous small square subdomains, each corresponding to an independent realization. These subdomains are then stitched together to artificially expand the metasurface area, with field discontinuities smoothed at subdomain boundaries to mitigate artifacts.⁹¹⁷ In the supercell approach, identical subdomains are “stitched” with pseudo-periodic boundary conditions. As the artificial period increases, both approaches generally yield stable numerical outcomes for specular and diffuse light.⁹¹⁸,⁹¹⁹ Although finite-size computations unavoidably introduce boundary effects,⁹¹⁹ the observed relative errors are often encouragingly below 10%.

However, the two approaches necessitate numerous computations to assess statistical convergence, and they offer limited physical insights. The challenge of effectively simulating infinite metasurfaces remains, highlighting the need for approximate models.

### 3. Approximate models

There are a few available models for predicting the BSDF of disordered metasurfaces. Similar to research conducted on metamaterials, advanced studies have primarily focused on homogenization, i.e., on the specular contribution to the BSDF.⁸⁸³,⁸⁸⁴,⁹²²,⁹²³ To our knowledge, there has been no research on diffuse light, despite its critical significance in our context. An exception to this is the independent scattering approximation (ISA), which assumes that each constituent particle scatters light independently, thus neglecting multiple scattering effects. Although the ISA may appear simplistic to readers, it holds substantial value. This approximation plays a significant role in investigations on the physics of waves in complex media. It also largely facilitates a transparent interpretation that is beneficial during the initial design phases and is amenable to more accurate models relying on mean fields.⁹¹⁶ For a monolayer comprising N identical particles randomly positioned at coordinates $\mathbf{r}_p = [x_p, y_p, z_p = 0]$, $p = 1, ..., N$, the ISA simply leads to

$$ \langle I_s \rangle (\hat{\mathbf{k}}_s, \hat{\mathbf{e}}_s, \hat{\mathbf{k}}_i, \hat{\mathbf{e}}_i, \omega) = I_0 N \frac{d\sigma_s}{d\Omega} (\hat{\mathbf{k}}_s, \hat{\mathbf{e}}_s, \hat{\mathbf{k}}_i, \hat{\mathbf{e}}_i, \omega) S(\mathbf{q}), \quad (7) $$

for the total radiant intensity $\langle I_s \rangle$ (measured in W sr⁻¹), encompassing both coherent and incoherent components, scattered by the particles. In Eq. (7), $d\sigma_s/d\Omega$ (measured in m² sr⁻¹) is the differential scattering cross section of individual particles—a quantity directly related to the so-called form factor in condensed matter physics. $I_0$ is the incident

intensity (measured in W m⁻²). $S(\mathbf{q})$ is the structure factor, which incorporates the effect of far-field interference between pairs of particles on the scattered intensity, and is defined as

$$ S(\mathbf{q}) = 1 + \frac{1}{N} \left\langle \sum_{m \neq n}^{N} \exp[i\mathbf{q} \cdot (\mathbf{r}_m - \mathbf{r}_n)] \right\rangle, \quad (8) $$

where $\mathbf{q} = \omega/cn_b(\hat{\mathbf{k}}_i - \hat{\mathbf{k}}_s)$ is the wavevector difference between the incident and scattered waves and $n_b$ is the refractive index of the background medium.

The structure factor can be readily evaluated through statistical averaging. On the other hand, accurate determination of the form factor involves iteratively solving Maxwell equations for multiple wavelengths, incidences, and polarizations, followed by applying a near-to-far-field transformation to obtain the radiation pattern. Alternatively, a quicker approach with deeper insights can be achieved by expanding the differential scattering cross section using the resonance modes of the nanoparticles.⁹²³

Exploiting Eq. (7), one can derive straightforward expressions for the BSDF of infinite metasurfaces with $N \to \infty$. However, due to the inherent limitations of the ISA, the model's predictions exhibit notable inaccuracies for high particle densities, grazing incident angles, or large angles. These inaccuracies prompt the need for further refinements and extensions. For readers seeking a more comprehensive understanding, let us note that advanced models for the diffuse contribution in the BSDF have been presented in Ref. 918. For the specular contribution, an overview of various models can be found in a recent book chapter⁹²⁰ and the references therein. Noteworthy is the remarkable efficacy demonstrated by the quasi-crystalline approximation.⁹²¹,⁹²⁵

### 4. Fabrication

Several techniques exist for fabricating disordered metasurfaces. Top-down patterning methods, which make use of electron-beam or optical lithography, facilitate the precise positioning of the meta-atoms with nanometer-level precision and a fine-tuned control over meta-atom size and shape. Furthermore, with the advent of nanoimprint lithography⁹²⁴ and soft lithography,⁹²⁵ they have opened up new avenues for cost-effective fabrication.

Conversely, bottom-up approaches prioritize scalability and cost-efficiency over precision. Although they may not match top-down techniques in shaping and arranging nano-objects, their utility remains unequaled.

Among metasurfaces fabricated by scalable and cost-effective techniques, semicontinuous colloidal films composed of metal nano-islands stand out. These films naturally arise from deposition techniques that interrupt metal film growth before percolation threshold is reached.⁹²⁶ Originally investigated for enhanced Raman scattering,⁹²⁷ these ultra-affordable films can be mass-produced by the glass industry on extensive scales, spanning tens of square meters. Because island dimensions and interspace distances are much smaller than visible wavelengths, nanoisland metasurfaces act as coherent layers, absorbing rather than scattering incident light. Alternatives, such as utilizing block copolymer templates for self-assembly,⁹²⁸,⁹²⁹ also facilitate monolayer creation with such tiny particles.

Producing disordered configurations with larger particles involves assembling pre-synthesized particles onto surfaces. Various techniques fit within this category:
<figure><img src="image_39.png" /><figcaption>**FIG. 31.** A few selected examples of disordered metasurfaces fabricated with colloidal chemistry. (a) Availability of almost monodisperse nanospheres in various dielectric materials, e.g., PMMA, high-index Si nanospheres in the present figure.<sup>930</sup> (b) Amorphous colloidal monolayers of polystyrene nanospheres with sizes of 170 and 270 nm with short-range order.<sup>931</sup> (c) Monolayer of Au nanocubes with correlated disorder on a 2.5 cm-diameter half sphere.<sup>932</sup> All images are reprinted (adapted) with permission from the respective Journal and copyright remains with the original publisher.</figcaption></figure>

* Colloidal lithography: This method hinges on nearly uniform in size colloids in solutions of materials like polystyrene, PMMA, or Si nanospheres<sup>930</sup> [Figs. 31(a) and 31(b)]. After drop-casting and drying, a monolayer forms,<sup>931</sup> subsequently serving as a mask to imprint patterns into underlying substrates. Since the sizing of polymer beads and strategies for self-assembly—utilizing repulsive and attractive forces for instance—are well-mastered [Fig. 31(b)], the transferred patterns, often arrays of nanodiscs or nanoholes in various materials, are fabricated on large scales with high fidelity and few aggregates.

* Self-assembly: With this method, colloidal nanoparticles, initially stabilized in solution, arrange on surfaces via methods like blade coating, dropcasting, and dip coating.<sup>932</sup> Diverse particle shapes, sizes and materials are commercially available, and controlled short-range arrangements can be achieved through nanoparticle and substrate surface functionalization [Fig. 31(c)].

For deeper insights, recent comprehensive reviews provide substantial information regarding the subject. They additionally cover deposition methods enabling the creation of densely packed monolayers. These monolayers exhibit high sensitivity to electromagnetic near-field interactions. They offer complex and possibly rich optical responses due to enhanced particle couplings and spatial dispersion. Future research directions should, therefore, include a more thorough investigation of these dense configurations, extending modeling approaches far beyond the limits of ISA and mean-field theories.<sup>99,933–936</sup>

### C. Future directions and outlook

Disordered optical metasurfaces have unveiled a plethora of applications, some of which can be traced back to historical contributions, such as porous film-based anti-reflective coatings by Fraunhofer in 1887,<sup>900</sup> and observations of ultra-strong Raman signals from semicontinuous nanoisland metal films in the 1970s.<sup>902</sup> Over the past two decades, the field of disordered optical metasurfaces has experienced remarkable activity, encompassing a diverse array of applications. These include wide-angle broadband absorption,<sup>932,937,938</sup> anti-counterfeiting measures,<sup>939</sup> anti-reflection

<figure><img src="image_40.png" /><figcaption>**FIG. 32.** (a) Harnessing the BRDF toward full manipulation of diffuse and specular light. In the sketchy example, the surface color varies abruptly, differently from the continuous color variation of classical iridescence, when the light source and/or the observer position move. (b) Hyperuniform metasurface for light trapping in solar cell application. After Ref. 968. (c) Rendered images showing the visual appearances of a car coated with disordered arrays of identical Ag nanospheres with a density of 4 µm⁻¹. Top: no structural correlation. Bottom: short-range correlated disorder. Adapted from Ref. 918. (d) Transparent display made of a bilayer metasurface which offer completely different properties in reflection and transmission. The rightmost images are camera pictures of a flower bouquet taken in reflection (top) or in transmission (bottom). After Ref. 945. All images are reprinted (adapted) with permission from the respective Journal and copyright remains with the original publisher.</figcaption></figure>
coating,⁹⁰⁰,⁹¹⁷,⁹⁴⁰ wavefront manipulation,⁷⁰¹,⁹⁴¹,⁹⁴² transparent displays,⁹⁴³,⁹⁴⁴ chiral films,⁹⁴⁵ fade-resistant structural coloring,⁹³⁸,⁹⁴⁶⁻⁹⁵⁰ visual effects,⁹¹⁸,⁹²³ plasmonic electronic paper,⁹⁵¹ solar radiation management in architectural and consumer products glazing,⁹⁵²⁻⁹⁵⁵ radiative cooling,⁹⁵⁶ photocatalysis,⁹⁵⁷ light trapping in solar cells,⁹³¹ and light extraction in LEDs.⁹⁵⁸⁻⁹⁶¹

As nanofabrication capabilities continue to advance, these applications are progressively moving closer to industrial-scale production. In the subsequent sections, we delve into emerging applications that are poised to inspire numerous studies and emphasize two prospective avenues for enhancing manufacturing flexibility and reducing manufacturing costs.

## 1. Prospective applications

Anticipating the forthcoming wave of research, we expect a concerted effort to tackle the enduring question of how to master the nanopatterning of surfaces for shaping the far-field radiation across both spatial and spectral frequency domains. The ultimate ambition is to attain comprehensive command over both specular and diffuse light. It encompasses manipulating independently these two components, sculpting the directional patterns of diffuse light, monitoring alterations in light color for varying viewing and illumination angles, and, in essence, molding the BSDF at will [Fig. 32(a)].

Two key physical parameters that play pivotal roles can be used for design [see Eq. (7)]. Through the form factor, an assortment of shapes and high-index materials may be used to manipulate the resonances of the meta-atoms at the nanoscale. These resonances may be further enriched through mode hybridization at the wavelength scale to achieve wavelength-selective control over the directionality of scattered light.¹⁷¹,⁹⁶²,⁹⁶³ Through the structure factor, we may engineer intricate mesoscale interferences. This engineering encompasses a vast landscape of arrangements,⁵²⁴,⁹⁶⁴ ranging from full periodicity to entirely uncorrelated disorders. This gamut includes quasicrystals, weakly defective crystals, hyperuniform configurations, and more. For example, enforcing greater inter-particle distances engenders a short-range amorphous-type correlation that mitigates the presence of diffuse light around the specular direction,⁹¹⁸,⁹⁶⁵ regardless of the angle of incidence.

This effect is even more pronounced within stealthy hyperuniform metasurfaces, which entirely suppress scattering within an angular range centered on the specular direction.⁵²⁴ Ongoing endeavors in the domain of hyperuniform disordered scattering structures, particularly for solar cells, are only just emerging.⁹⁶⁵⁻⁹⁶⁸ A particularly remarkable recent achievement is the 2.5-fold enhancement in absorption of a 1 µm-thick silicon slab, thanks to a hyperuniform metasurface [Fig. 32(b)] meticulously designed to trap incident light into the slab guided modes.⁹⁶⁸

Disordered optical metasurfaces also hold significant promise for applications in applied and fine arts. Metasurfaces offer the potential to manipulate the visual characteristics—such as color, glossiness, haze, transparency, and iridescence—of macroscopic objects of any conceivable shape.⁹¹⁸ A recent breakthrough in this domain introduced a multiscale numerical tool capable of generating lifelike images of diverse macroscopic objects covered with disordered metasurfaces [Fig. 32(c)]. This tool not only forecasts subtle effects resulting from the interplay between diffuse and specular light but also unveils entirely new visual effects.

<figure><img src="image_41.png" /><figcaption>FIG. 33. (a) Metasurfaces fabricated on curved substrates: air holes etched in pyramid-textured Si surfaces. Adapted from Ref. 975. (b) Production of vivid colors despite a strong polydispersity in the size and shape of Al nanoislands used to resonantly absorb light. The resonance is due to nanogap plasmons supported by the capped nanoislands obtained by dewetting on Al₂O₃/Al substrate (inset). The photography shows a black canvas painted with oil-based plasmonic paints obtained by dispersing metasurface flakes in commercial paint oils. Adapted from Ref. 950. All images are reprinted (adapted) with permission from the respective Journal and copyright remains with the original publisher.</figcaption></figure>

Future exploratory directions may capitalize on the monolithic integration of multiple metasurface layers. A compelling feature of these surfaces is their distinct behavior concerning diffuse light in transmission vs reflection. This contrast can be readily comprehended by examining the structure factor, S(q) in Eq. (8), in which q signifies the difference in wavevector between the scattered direction and incident directions. It becomes evident that $(k_i - k_s) \cdot (r_m - r_n)$ in the exponential term remains independent of the normal coordinate z around the specular transmitted direction (where $k_s \equiv k_t \approx k_i$). In contrast, this term does vary with the normal coordinate around the specular reflected direction (where $k_s \equiv k_r \approx -k_i$). This intrinsic property has recently been ingeniously harnessed to create a transparent display,⁹⁴⁴ that selectively scatters reflected light while preserving transmitted light integrity [Fig. 32(d)].

## 2. Metasurfaces on curved substrates

An overwhelming number of nanopatterning methods, including optical and electronic lithography, nanoimprint lithography,⁹²⁴ soft-lithography,⁹²⁵ DNA-templating,⁹⁶⁹ and functional block copolymers,⁹⁷⁰ are available nowadays. They often function on 2D and atomically flat surfaces.

In general, adapting methods tailored for flat surfaces to curved ones necessitates significant changes to mitigate curvature-induced variations in feature attributes like size, orientation, and spacing. While some fabrication methods can accommodate relatively simple curvatures, the adaptation challenge grows as curvature intensifies. This puts severe limitations onto the use of metasurfaces in a variety of every-day-use or advanced technologies, where materials with non-planar, rough, soft, mechanically agile, or chemically sensitive surfaces are used.⁹⁷¹⁻⁹⁷³

Recently, several transfer protocols have been proposed to overcome the limitations.⁹³²,⁹⁷³⁻⁹⁷⁵ Figure 33(a) shows the example of a pyramid-textured Si surfaces which has been patterned with air holes to lower reflection at the Si-air interfaces and enhance the total absorption. The hole pattern is obtained from a template manufactured with colloidal lithography (other bottom-up or top-down approaches may be used as well). In many instances, the key is the use of a 10 nm-thin carbon film deposited on a sacrificial layer for the template. Such a
thin film combines extreme strength with a large degree of flexibility that facilitates conformation of the transferred pattern to nearly any macroscopic or microscopic 3D object.

### 3. Striving for ultra-low costs with augmented resilience to fabrication imperfections

Optical metasurfaces, which harness resonant nanostructures at subwavelength scales, impose strong constraints on the nanofabrication precision. These constraints are growing in significance, as metasurfaces progress toward practical applications that require large sample areas. In this evolving landscape, the emergence of metasurfaces designed to exhibit resilience against imperfections becomes increasingly crucial. Recent studies have converged to provide compelling evidence that disordered arrays of nanoparticles, when placed on a transparent dielectric spacer above a metallic substrate, exhibit an extraordinary resilience to variations in material, size, and nanoparticle shape. One fundamental mechanism that imparts resilience is through the process of mode hybridization.⁹²³ In scenarios where plasmonic nanoparticles are situated on low-index substrates without undergoing hybridization, the electromagnetic field becomes profoundly confined around the individual particle. Consequently, even slight variations in material composition or shape induce significant shifts in resonance frequencies. For context, even a minor alteration of 10 nm in the dimensions of a metal nanoparticle leads to a resonance wavelength shift of 10–20 nm. In contrast, within hybridized modes, the field localization diminishes as it extends into the transparent layer. Consequently, hybridization inherently mitigates sensitivity to size and material polydispersity. This general attribute of hybrid modes has facilitated the discovery of a remarkable diffuse iridescence phenomenon, in which a mere two distinct and vivid colors are discernible across all angles of illumination and viewing, even when confronted with significant variations in nanoparticle size and shape across the metasurface.⁹²³

An alternative resilience mechanism emerges when nanoparticles are substantially smaller than the wavelength. This scenario is exemplified by semicontinuous nanoisland films on a dielectric spacer atop a metal substrate, where empirical observations reveal that the pivotal factor shaping the specular reflectance spectrum is the spacer thickness,⁹⁴⁷ rather than the typical expectation of the nanoparticle material choice. This material-independent response has been experimentally verified for various materials like Ag, Au, Pt, Pd, and numerically for Al, Cr, Cu, Ni, among others, under normal incidence.⁹⁴⁷ Thus commonly accessible materials such as Cu or Al could substantially mitigate costs related to nanostructure-based coloration.

This striking property can be attributed to the amalgamation of localized plasmonic resonances at the nanoisland scale, hot spots within the interstitial gaps, and delocalized modes, irrespective of the metal used for the nanoislands. This synergy results in broad absorption, further enhanced by the diversity in island size and shape. Consequently, the film uniformly absorbs the incident light, except for the Fabry–Pérot resonance of the thin transparent spacer, for which the field is nearly null within the semicontinuous film.

Yet another illustration of resilience against manufacturing defects is observed in a recent work on ultralight plasmonic structural-color paints⁹⁵⁰ with a notably low-cost morphology shown in the inset of Fig. 33(b). Impressively, vibrant colors are attainable despite the substantial variance in nanoisland sizes. While plasmon hybridization

leading to a confined gap-plasmon resonance within the alumina film stands as a logical explanation, it is important to acknowledge the influence of other mechanisms as well. Among these, the generation of subtractive color through the absorption of a “homogenized” film emerges as a noteworthy consideration, deviating from the diffusion-based coloration observed with larger nanoparticles.

### ACKNOWLEDGMENTS

P.L. and K.V. acknowledge financial support from the French National Agency for Research (ANR) under the project “NANO-APPEARANCE” (No. ANR-19-CE09-0014). A.D. acknowledges the Swedish Research Council for Sustainable Development (Formas) (Project No. 2021-01390). P.L. acknowledges financial support from the European Research Council Advanced grant (Project UNSEEN No. 101097856). C.R. and A.S. acknowledge financial support from the Deutsche Forschungsgemeinschaft (DFG) (Nos. 413644979, RO 3640/11-1, and WE 4051/26-1).

## XVII. FROM PRINCIPLES TO PERFORMANCE: EPSILON-NEAR-ZERO METASURFACES FOR NONLINEAR OPTICAL DEVICES

Jeremy Upham, M. Zahirul Alam, Israel De Leon, Robert W. Boyd, and Sebastian A. Schulz*

*sas35@st-andrews.ac.uk

Epsilon-near-zero (ENZ) materials, particularly thin films of transparent conductive oxides, exhibit a variety of intriguing optical phenomena, two of which should be of particular interest to the metasurface community: (1) ENZ materials of sub-wavelength thicknesses can exhibit large and very fast nonlinear optical responses, including significant nonlinear phase shifts; (2) the nanoantenna of a metasurface can strongly couple to long-range surface plasmon polariton modes of these low index thin films, making these nonlinear optical responses even stronger. In this article, we examine how these enormous nonlinearities are currently understood and highlight examples of how this functionality is already being exploited for applications in nonlinear optics, such as frequency conversion, self-phase modulation, time-refraction, and beyond. Looking forward, the current limitations on ENZ-metasurface performance appear to be predominantly engineering challenges rather than physical limitations. We argue that through creative metasurface design and intuitive materials science optimization, ENZ-metasurfaces are poised to transform from compelling proofs of principle to high-performance, all-optical devices.

### A. Introduction: The physics of ENZ

Exciting subsets of current metasurface research include the development of dynamic, temporally varying devices and the enhancement of nonlinear optical responses to expand metasurface functionality. Fortunately, this has proven to be achievable simply by incorporating thin film layers of a class of materials with vanishing permittivity known as epsilon-near-zero (ENZ) materials. A material is said to have ENZ properties if its dispersion curve features a spectral region where $-1 < \text{Re}(\varepsilon) < 1$, typically with the ENZ-crossing being the wavelength where the real part of the permittivity crosses from
<figure><img src="image_42.png" /><figcaption>FIG. 34. Dispersion of the real and imaginary components of the permittivity, as well as the real part of the refractive index, for a typical ENZ material. $\lambda_{ENZ}$ indicates the wavelength where the real part of the permittivity crosses zero. The shaded region indicates the ENZ spectral range where $-1 < \text{Re}(\varepsilon) < 1$. Because the real part of the refractive index also depends on the imaginary part of the permittivity, the refractive index minimum is near, but not necessarily within the ENZ spectral range.</figcaption></figure>

positive to negative. This effect leads to a wide variety of unconventional optical behaviors in both the linear and nonlinear regimes.⁹⁷⁷,⁹⁷⁸ For example, waveguides formed from ENZ materials can feature coupling through arbitrary shaped channels, phase-matching free nonlinear processes or manipulation of radiation from nearby or embedded emitters.⁹⁷⁷ In other platforms more relevant to metasurface technologies, ENZ films have shown enhancement of nonlinear effects such as harmonic generation⁹⁷⁹ and the nonlinear refractive index.¹⁸⁴,¹⁸⁵

To clarify, while the ENZ condition typically has the real part of the permittivity cross from positive to negative, the imaginary component need not change sign, nor does the permeability. Therefore, the refractive index need not become negative as in the case of negative index metamaterials⁹⁸⁰,⁹⁸¹ or even reach zero as in near zero index materials.⁹⁸²,⁹⁸³ The ENZ condition often occurs naturally in Drude materials simply by the impact of the plasma frequency on the dispersion curve, delineating where a material crosses from metallic to dielectric behavior. One can also obtain ENZ conditions in a Lorentz-type material for a sufficiently large value of the oscillator strength. While a near-zero real part of the permittivity can also be achieved for values of the refractive index larger than one in such materials at absorption resonances, these are typically accompanied by high imaginary permittivity values. In these cases, such excessively high losses prohibit the observation of the desired ENZ effects. Therefore, we generally consider ENZ materials to be Drude-type materials where both the absolute real part of the refractive index and the permittivity are less than unity over some range of frequencies. The schematic of Fig. 34 shows the dispersion of permittivity and refractive index for a typical ENZ material. The spectral range where the real part of the permittivity has a magnitude less than one is the ENZ range. The real part of the refractive index can be below 1 in this same region, but this also depends on the imaginary component of the permittivity.

## 1. Materials with ENZ properties

Depending on what part of the EM spectrum is being considered, there are multiple ways to achieve this ENZ condition. At radio frequencies well-established metamaterials such as arrays of split ring resonators can easily achieve this state,⁹⁸⁴,⁹⁸⁵ and similarly, ENZ behavior has been shown in microwaves⁹⁸⁶ and through phononic resonances in the THz spectral region.⁹⁸⁷,⁹⁸⁸ Within the deep UV, many metals have ENZ crossings (e.g., silver near $\lambda = 320$ nm⁹⁸⁹) For the visible and near- to mid-infrared spectral regions, two solutions exist. The first is to engineer structures, either photonic crystal-like band structures⁹⁸² or multi-layer materials, with an effective index⁹⁹⁰⁻⁹⁹³ such that their effective permittivity is near zero. However, in the context of metasurfaces, the prevailing solution has been to use transparent conductive oxides (TCOs).⁹⁹⁴,⁹⁹⁵ These non-stoichiometric materials are typically transparent dielectrics such as indium oxide, aluminum oxide, or cadmium oxide that are degenerately doped so that the Fermi level moves to the conduction band. These TCOs hit a sweet spot for ENZ performance in the near to mid-infrared range: they feature sufficient free carriers to achieve a crossing of the real part of the permittivity from positive to negative, but still fewer than noble metals, and hence the crossing occurs at longer wavelengths. Furthermore, the concentration and mobility of the carriers (the key factors determining the plasma frequency) can be controlled during TCO deposition, e.g., by tuning the doping concentration.⁹⁹⁶⁻⁹⁹⁹ Together with the naturally slightly different plasma frequency depending on the dopant (e.g., indium, gallium, aluminum), this means that TCOs can place thisepsilon near zero condition anywhere within a wide range of the near- to mid-IR spectrum, covering important regions for optical communications and sensing.¹⁰⁰⁰

Certain TCOs, such as indium tin oxide (ITO) and aluminum zinc oxide (AZO), have become reasonably mature fabrication materials because they can be deposited as transparent yet conductive thin films that enable touchscreen and solar cell technology.¹⁰⁰¹,¹⁰⁰² Driven by these applications to be transparent to visible light (dielectric in the visible spectrum) but also as conductive as possible (lots of free carriers), they are readily available materials for thin films with a plasma frequency just past the visible spectrum in the near-infrared. As a result, these thin films typically exhibit ENZ behavior between 1 and 2 $\mu$m, with the additional benefit of a reasonably small absorptive loss.⁹⁹⁶ Because the real part of the permittivity crosses zero at a wavelength where the imaginary component is also small, the refractive index is less than 1 as well.

## 2. ENZ induced behavior in TCOs

A consequence of the refractive index falling below 1 is that the wavelength of light within that material is stretched compared to that of free space. Equivalently, as the index approaches zero, spatial phase variation within the ENZ medium also goes to zero. Because light in a zero index medium can only radiate away normal to the bounding interfaces, this leads to interesting linear optical effects,⁹⁷⁷ such as light propagating through arbitrary channels,¹⁰⁰³ or wavefront shaping at interfaces.¹⁰⁰⁴

However, from the point of the original motivation, i.e., the desire for large, dynamic modulation of metasurface behavior, the nonlinear enhancement in ENZ materials is of more interest. Many conventionally held equations for nonlinear optical phenomena, such as the Kerr
<figure><img src="image_43.png" /><figcaption>**FIG. 35.** First experimental demonstrations of large, optically driven changes of refractive index of TCOs in their ENZ spectral regions. (a) Pump-probe measurement of $\Delta n = 0.45$ for AZO at 1390 nm. Reproduced with permission from Caspani et al., Phys. Rev. Lett. **116**, 233901 (2016). Copyright 2023 American Physical Society.¹⁸⁵ (b) Z-scan measurements of ITO showed a $\Delta n = 0.72$ at 1240 nm.¹⁸⁴ All images are reprinted (adapted) with permission from the respective Journal and copyright remains with the original publisher.</figcaption></figure>

effect, are written to be inversely proportional to the real part of the linear refractive index, e.g., $n_2 = \frac{3\chi^{(3)}}{4\varepsilon_0 c n_0 R(n_0)}$ (where $n_2$ and $n_0$ are the nonlinear and linear refractive indices, $\chi^{(3)}$ is the third order susceptibility and $\varepsilon_0$ is the vacuum permittivity), suggesting that such phenomena would be enhanced or even diverge within an ENZ material.¹⁰⁰⁵ This led to the prediction¹⁰⁰⁶ and demonstration¹⁰⁰⁷ of ENZ materials boosting the nonlinear generation of harmonic frequencies.

In 2016, two studies on nonlinear optics in TCOs investigated the effect of all-optical nonlinear interactions on light propagation and reflection within ITO¹⁸⁴ and AZO films.¹⁸⁵ Both studies demonstrate very strong enhancement of the nonlinear refractive index, see Fig. 35, with values of $\Delta n = 0.72$ and $\Delta n = 0.45$, respectively, corresponding to up to 500% of the linear refractive index,¹⁸⁴,¹⁸⁵ reaching the regime where the nonlinear response exceeded the linear component of the refractive index. We thus note that the conventionally held equations of the Kerr and other nonlinear effects that are based on perturbative descriptions are no longer appropriate. ENZ materials feature such a strong nonlinear response that the equations of nonlinear optics must

be reconsidered, specifically the assumptions of $n \le 1$ and $n \ll \Delta n$ are no longer valid, making a compelling case for the use of ENZ materials in nonlinear optics.¹⁰⁰⁸

While other, non-perturbative nonlinear effects have been demonstrated before (e.g., the semiconductor-to-metallic light- or heat-induced phase transition in vanadium dioxide¹⁰⁰⁵) such strong effects are typically slow, as such thermo-optic or thermo-chroic responses occur on the ms timescale or longer. In contrast to this, the nonlinear effects in the TCO materials occur on an ultrafast timescale.¹⁸⁴,¹⁸⁵,¹⁰¹⁰ Often the onset seems to match the temporal width of an ultrafast pump pulse and relaxation times are on the scale of a few-hundreds femtosecond.¹⁸⁴,¹⁸⁵,¹⁰¹⁰ This holds potential for THz switching speed operation of ENZ-based devices. The underlying nonlinear response of the TCOs is dominated by the free carriers, resulting in a hot electron-dominated effect. When optically pumped, the available carriers in the conduction band experience intraband transitions that significantly change their distribution and consequently the plasma frequency of the material. Because the transitions are principally intraband, the electrons cool quickly, explaining the rapid onset and relaxation of the index change. These transitions are so fast that light witnessing changes of the refractive index will undergo an equivalent change of frequency through adiabatic frequency conversion, otherwise known as time refraction.⁸⁵²,¹⁰¹¹,¹⁰¹² discussed further in Sec. XV.

Interestingly, the magnitude of the observed index change stems from the fact that the response is free carrier dominated, even though these are present in low concentrations. Because of the significantly lower carrier concentrations in the conduction band (compared to conventional metals), the excitation of even a modest amount of carriers is enough to significantly alter the overall carrier distribution, the Fermi level, and consequently the effective carrier mass and plasma frequency. The Drude model can thus explain the change of the TCO's dispersion curve almost entirely as a change of the electron effective mass due to the changed effective electron temperature¹⁰¹³ (see Fig. 36). The observed response is all the more dramatic in the ENZ region, where the initially low linear index can be effectively tripled by the magnitude of the nonlinear index change.

While these rapid and robust index changes show significant potential for rapidly modulated optical components, the transitions are not efficient. As a low-index material and a conductor, the index contrast or impedance mismatch with optical modes in free space

<figure><img src="image_44.png" /><figcaption>**FIG. 36.** Modeling changes to the real part of the permittivity of ITO as a function of the electron temperature. The green curve (1240 nm) has a real permittivity of 0 when $T_e$ is at room temperature, but increases with electron temperature. When heated up to 4000 K, the black curve (1300 nm) becomes ENZ crossing wavelength. Reproduced with permission from Baxter et al., Adv. Photonics Res. 4(3), 2200280 (2023). Copyright 2023 Wiley.¹⁰¹³ All images are reprinted (adapted) with permission from the respective Journal and copyright remains with the original publisher.</figcaption></figure>
<figure><img src="image_45.png" /><figcaption>**FIG. 37.** Experimental demonstrations of plasmonic metasurfaces enhancing the nonlinear response of TCOs. (a) Plasmonic dimers on ITO show a large, fast (sub-picosecond) modulation of optical density when exciting antenna resonances, suggesting a strong, fast Kerr-type nonlinearity.<sup>1018</sup> (b) A sufficiently thin, bare layer of ENZ material will support a Berreman (LR-SPP) mode with high field enhancement, but cannot easily couple to light in free space.<sup>189</sup> (c) Nanoantennas with resonances in the ENZ region of a thin TCO layer exhibit mode splitting from strong coupling between the metasurface and the LR-SPP mode.<sup>1015</sup> (d) When a metasurface enables optical pumping of this LR-SPP mode, it not only shows an enhanced nonlinear response but can also extend this nonlinear enhancement well beyond the narrow spectral range where ε is near zero.<sup>184</sup> All images are reprinted (adapted) with permission from the respective Journal and copyright remains with the original publisher.</figcaption></figure>

makes coupling light into ENZ materials difficult. For this to be a truly useful optical component, we must overcome this challenge.

### 3. Incorporating metasurfaces with ENZ materials

Fortunately, metallic nanoantennas can effectively couple light between propagating and localized modes. If a plasmonic nanoantenna is placed on the surface of an ENZ material, the incident light can efficiently couple to the local surface plasmon resonance and then have the fields of this resonance penetrate the TCO and excite free carriers<sup>1014</sup> [Fig. 37(a)]. Excitation of free carriers is further enhanced if the ENZ material can support a Berreman mode, which is a Long Range Surface Plasmon Polariton (LR-SPP) mode traveling in the low index medium of a sufficiently thin ENZ material with field enhancement inversely proportional to its thickness<sup>189</sup> [Fig. 37(b)]. In fact, plasmonic metasurfaces of appropriately designed nanoantennas can even exhibit strong coupling with the LR-SPP mode of a thin film of TCO,<sup>1015,1016</sup> which greatly improves energy transfer into the TCO through near-field interaction [Fig. 37(c)]. Critically, this strong coupling between modes in the ENZ material and the nanoantenna enhances the nonlinear response of the metasurfaces well beyond what the bare ENZ material could achieve: Light couples through the antenna into the ENZ film, changing the latter's refractive index. This change modifies the dispersion of the TCO's LR-SPP modes, which in turn modifies the

coupled resonance spectrum. As a result, light of high intensity experiences a dramatically different refractive index compared with that experienced by low-intensity light<sup>154</sup> [Fig. 37(d)]. Therefore, pairing an array of plasmonic nanoantenna with an ENZ thin film will not only improve energy transfer to the free carriers within the material but can further increase the nonlinear phase shift of the metasurface by its shifted resonance. Significantly, the broadening of the split resonance also extends the spectral range over which a large nonlinear phase shift can be achieved to be well beyond the ENZ region of the bare TCO material.

This synergistic interaction enables a faster, larger, and potentially more flexible nonlinear optical response from a sub-wavelength ENZ metasurface than can be achieved otherwise, significantly outperforming thicker layers of TCOs or even highly nonlinear bulk materials.<sup>1017</sup> Having established that a metasurface not only gains a large, dynamic optical response from incorporating an ENZ material but in fact can magnify the size of this response through strong coupling to the LR-SPP modes, we next consider how this functionality can be put to use by metasurface design.

### 4. Metasurfaces exploiting the nonlinear ENZ response

Having motivated the potential of ENZ-based metamaterials for nonlinear optics, we next consider some recent, exciting examples of devices where this strategy has been put into practice.
<figure><img src="image_46.png" /><figcaption>FIG. 38. Compelling results in ENZ metasurfaces: (a) dynamically changing the direction and optical frequency of light diffracting from an ENZ-metasurface is different for diffraction orders.879 (b) SHG is enhanced by a factor of 10⁴ and is polarization selective.1019 (c) Enhanced THz generation from an ENZ-metasurface that also dictates polarization and wavefront.1020 (d) Four wave mixing enhanced by dynamic metasurface-ENZ strong coupling.851 (e) Enhanced adiabatic frequency conversion.875 (f) Strong coupling and dynamic modulation of a dielectric metasurface on a thin TCO layer, for both the electric or magnetic dipole Mie resonances.1021 All images are reprinted (adapted) with permission from the respective Journal and copyright remains with the original publisher.</figcaption></figure>

An example of time-refraction in an ENZ metasurface was realized in a system consisting of a plasmonic blazed diffraction grating on an ITO thin film with a pump-probe excitation. In this instance, the pump pulses cause the diffracted probe to undergo a change of central wavelength while interacting with the metasurface and therefore have the diffraction angle dictated by the probe's new wavelength. Interestingly, the light diffracting into the -1 diffraction order (the order suppressed by the blazing) appears to interact with the metasurface longer than light diffracting to the +1 order. Consequently, light in the -1 diffraction order experiences about twice as much wavelength shift and thus twice as much change to the diffraction angle as it's counterpart.879 [Fig. 38(a)].

Both the previously demonstrated large refractive index change, as well as the time refraction are demonstrations of the Kerr-type effect, i.e., an intensity-dependent refractive index change. However, ENZ-metasurfaces can enhance a wide range of other nonlinear effects. For example, Deng et al. used a hexagonal array of Y-shaped nanoantennas to enhance the second harmonic conversion from an ITO thin film.1019 In this work the antenna array once again overcomes the issue of coupling light into the ENZ film, resulting in an experimental observation of a 10 000-fold enhancement of the SHG efficiency. Furthermore, this work demonstrates the multi-functionality of ENZ metasurfaces, as the symmetry properties of the antennas are used to control the polarization of the second harmonic light. For example, the metasurfaces can be configured such that excitation with a right-hand circularly polarized pump results in left circularly polarized SHG

emission. SHG enhancement has also been demonstrated in other ENZ-metasurface geometries, such as through gap-plasmon enhancement, based on square-shaped antennas coupled to ENZ films.1022

ENZ metasurface enhancement has also been demonstrated for other frequency conversion processes. For example, split ring resonators (SRRs) on ITO metasurfaces have been shown to result in orders of magnitude enhancement in the THz generation efficiency, compared to the same SRRs on glass,1020,1023,1024 which allows the generation of single-cycle THz pulses. This enhancement can be directly linked to the hot-electron dynamics in the ENZ film, also altering the spectral properties of the emitted THz pulse, with spectral broadening being observed.1023

A wide range of other nonlinear processes are enhanced in such ENZ-metasurfaces, with the geometry of both the ENZ film (e.g., thin films featuring the LP-SPP mode vs thicker films not featuring this guided mode) as well as the antenna shape and dimensions together with the excitation scheme controlling the achieved nonlinear effect. For example, an array of rectangular antennas on an ENZ film was examined in a pump-probe configuration, however, with the pump at double the probe frequency.851 This results in a number of nonlinear optical processes, including four-wave-mixing and negative refraction.

Furthermore, ENZ-metasurfaces retain the design freedom of linear metasurfaces. For example, through rectangular antennas we can design the strong coupling—and hence nonlinear enhancement—to be present for only one polarization,154,851 or through the use of antennas with a 90°
rotation symmetry, e.g., + shaped antennas, we can obtain polarization independent nonlinear enhancement.⁸⁷⁵,¹⁰²⁵ In a polarization-dependent setting, we can use the dynamic nature of the device to obtain controllable polarization rotation.¹⁰²⁶ Here for a given incident polarization featuring components both parallel and perpendicular to the antenna—e.g., of axis linear polarization or circular/elliptical polarization—the two polarizations will experience a different delay due to their varying interaction with the metasurfaces. However, this interaction will change for the parallel polarization as the antenna is excited nonlinearly (and its resonance shifts), while the cross-polarization retains the same interaction as in the linear case. We thus have the ability to control the delay of one polarization component relative to the other.

While ENZ-metasurfaces based on plasmonic nanoantennas are used in the majority of nonlinear optical investigation efforts at this time, it is important to also consider how other nano-resonances could be integrated with ENZs as well. Dielectric metasurfaces and their corresponding Mie resonances do not provide the same field enhancement in the TCO by the simple fact that their mode volumes lie principally within the dielectric material of the resonator rather than at the edges. Nevertheless, experiments have shown strong coupling between an LR-SPP mode and either electric or magnetic dipole resonances,¹⁰²¹ as well as confirming that these resonances can be optically tuned in a pump-probe set-up with sub-picosecond recovery time. Just as with plasmonic metasurfaces, strong coupling between the Mie resonances and the ENZ mode of a TCO thin-film appears to enable large, rapidly tunable modulations of the optical modes and should thus facilitate a broad range of dynamic, nonlinear optical responses.

It is evident from these examples that ENZ-metasurfaces enable strong and fast nonlinear optical responses while still maintaining the tremendous design freedom available to metasurfaces. There is a significant opportunity for devices that rapidly manipulate the phase, direction, polarization, and even mode of propagating light, all within flat, sub-wavelength structures.

## 5. Future of applications

Going forward, there are still several outstanding research questions that need to be addressed for the nonlinear response of ENZ materials in general and optimizing the potential of ENZ-metasurfaces in particular.

a. *Speed of modulation.* For the experiments presented here, the nonlinear optical response of the ENZ-metasurfaces is clearly faster than the short (sub picosecond) pump pulses can resolve. Knowing it can be faster than 100 fs supports the working understanding that rapid heating and intraband transitions of the free carriers are driving the index changes in the TCO. However, a more direct investigation into the onset time of this response, how it saturates, the thermalization time as well as how each of these is affected by coupling to plasmonic nanoantennas could provide insight into devices that modulate at rates beyond GHz. The speed and magnitude of the nonlinear phase shifts are likely large enough to enable non-adiabatic transitions, which could have significant implications for applications in which mode degeneracy or spectrum engineering are important.

b. *Absorptive losses.* While the imaginary component of TCOs refractive indices are typically smaller than for properly conductive materials,⁹⁹⁶ there is still a significant amount of absorptive loss seen

through sub-wavelength layers and this is exacerbated by coupling with plasmonic resonances. Material science research into TCOs intends to strike a balance between strong nonlinear response and absorption¹⁰²⁷ and be able to produce materials tuned for specific applications. However, there must always be some absorption remaining for TCOs, making them ineffective for most low-light sensing or quantum optics applications where loss management is critical.

c. *Saturation and damage threshold.* The majority of the above examples of ENZ-metasurfaces are based on gold nanoantennas. It is well known¹⁰²⁸ that many metallic nanostructures and gold particles in particular, are very sensitive to absorption-induced heating and melting, with the melting point significantly lower than for the bulk metal. However, more thermally robust materials such as TiN have higher absorption losses. This therefore presents an interesting research area, where materials and antenna design can be optimized to allow higher optical powers and operating temperatures. This low damage threshold, together with the increased design freedom of the ENZ metasurfaces means that nonlinear saturation effects have not yet been studied extensively in these systems.

d. *Alternative modulation mechanisms.* In this paper we have focused on the ultra-fast, all-optical response of ENZ metasurfaces. However, the underlying tuning mechanism is much more general. We recall that the mechanism is based on a shift of the ENZ dispersion, induced by a change in the free carrier properties. Such changes can be induced in an alternative fashion, for example through electro-optic modulation¹⁰²⁹,¹⁰³⁰ or thermo-optic tuning,¹⁰³¹ creating the potential for dynamically reconfigurable metasurfaces. Here we note that the sign of the nonlinear response can differ compared to the hot-electron nonlinearities discussed earlier, as some of these tuning mechanisms, can include interband transitions, or extraction or addition of carriers, hence allowing either a blue or redshift of the plasma frequency.

e. *Pushing back technical limits.* Given the unprecedented large changes of refractive index observed in ENZ-metamaterials, a natural question is how much further such devices be engineered before reaching a new limit? Some examples still needing exploration include, can an ENZ-metamaterial be engineered to produce a rapidly modulating nonlinear phase shift of $\Delta\phi = \pi/2$? In the electro-optic modulation scheme, this full phase shift has been achieved, however, at the cost of a highly reduced efficiency¹⁰³⁰ and limited to the tuning speed of the control electronics. Reaching the $\Delta\phi = \pi/2$ limit all-optically would make it possible to create a sub-wavelength, THz-speed switch.

In adiabatic frequency conversion work, it's been shown that large changes of the carrier frequency can be achieved in TCOs¹⁰³¹ and that plasmonic metasurfaces can manipulate this response more efficiently.⁸⁷⁵,⁸⁷⁹ However, in time refraction the fractional wavelength shift is proportional to the refractive index change witnessed by the photons, which so far has been far less than the full index shift of the TCO ($\Delta n \sim 1.0$). If ENZ-metasurfaces could be designed to increase this interaction time, it should in principle be possible to shift the carrier frequency of light by an octave.

## B. Perspective on future work

We can imagine that ENZ metasurfaces could be designed for a much wider set of nonlinear interactions. For example, some of the
most widely used interactions in commercial devices featuring nonlin-
ear optics are the sum and difference frequency generation in optical
parametric amplifiers and oscillators. Similarly, the field of attosecond
science relies on efficient high-harmonic generation. As ENZ metasurfaces
have shown enhancement of second harmonic light it would be a natural
next step to investigate and extend this to higher nonlinear orders.
Furthermore, the broad response of ENZ metasurfaces could also be of
interest to other processes such as supercontinuum generation and cas-
caded nonlinearities. However, these applications typically require nonlin-
ear interactions over an extended propagation length. How this could be
replicated in a metasurface geometry requires further research.

As the combination of metasurface resonances with TCO thin films shows such potential, a natural next question is what functionality can arise from patterning the metasurface directly into the TCO thin film itself? Initial investigations suggest sufficiently good coupling and field enhancement to produce significant nonlinear optical effects and non-reciprocal radiation patterns,¹⁰³² so further experimentation could be fruitful.

The majority of the nonlinear optical work on ENZ metasurfaces
has been performed in the near-infrared (1100–2000 nm) spectral
region. Naturally, these same effects hold promise for other pump
wavelength regions, such as visible light, but also mid-infrared through
to THz spectral regions. Implementing these effects across different
wavelength regimes will likely require different ENZ materials or at
least wide tuning of the ENZ wavelength.

Perhaps most intriguing, the field of photonic metasurfaces is
rapidly becoming more sophisticated, as is readily demonstrated by
this roadmap. As new materials, resonance structures and metasurface
principles are being introduced, there are opportunities to investigate
each of their interactions with ENZ thin films. For instance, incorpo-
rating bound states in the continuum, disordered metasurfaces or
active emitters with ENZ thin films could demonstrate the enhance-
ment or quenching of various light-matter interactions.

C. Conclusion

In this perspective, we have summarized the physics governing
nonlinear optical interactions in ENZ metasurfaces, consisting of
antenna arrays coupled to thin films of ENZ materials. The strong cou-
pling results in huge nonlinear enhancements, exceeding the perturba-
tive regime, on sub-picosecond timescales. We then highlighted some
recent key contributions to this rapidly growing field, well aware that
an even broader body of work exists. Considering the potential future
directions for this field and identifying opportunities across physics,
engineering, and materials science, there is clearly a strong potential
for novel nonlinear optical applications. In effect, almost any applica-
tion based on the mentioned nonlinear processes could benefit from
the enhancement in ENZ metasurfaces. Finally, as there remains an
even wider set of nonlinear optical interactions to which ENZ metasur-
faces have not yet been applied, we leave them as an exercise to the
reader to make their own mark on this exciting and growing field.

XVIII. DEEP LEARNING FOR PHOTONIC
METASURFACES

Willie J. Padilla\* and Jordan M. Malof

Deep learning (DL) is a specific type of artificial intelligence and has recently garnered tremendous attention within science and engineering, including for the design and study of photonic metasurfaces. However, there is a tendency to overvalue the current capabilities of DL, often attributing powerful reasoning skills to it. Here we aim to clarify the current strength of DL methods for photonic metasurface research and development. We highlight two major successful use-cases. First, DL has been highly effective at automatically constructing models that accurately predict the properties of complex metasurfaces based upon their experimental parameters (e.g., geometric structure, or illumination). Additionally, DL models has been used to build accurate models allowing users to predict the experimental parameters needed to yield some desired behavior, significantly accelerating metasurface design. Despite these successes it is important to recognize that DL primarily serves as a collection of advanced regression methods, which represent a relatively rudimentary form of intelligence compared to, for example, human scientists. The fundamentally limited reasoning capabilities of DL manifest as key practical limitations that undermine its current value and impact, such as DL's need for large quantities of training data, and their limited generalization to new scenarios that were not represented in training data. We discuss several promising acute and long-term strategies for addressing these limitations, which may ultimately lead to DL models that fulfill more of their current promise.

A. Introduction

Photonic metasurfaces have demonstrated great utility as a mate-
rials platform for the demonstration of numerous exotic electromag-
netic (EM) phenomena in both applied¹⁰³³ and fundamental¹⁰³⁴⁻¹⁰³⁶
settings. A common objective of photonic metasurfaces is to tailor
the EM field to achieve a desired configuration which may consti-
tute various forms, including free-space scattered states, surface
wave conditions,¹⁰³⁷ or localized field concentrations.⁶⁸² Moreover,
metasurfaces enable this control over diverse dimensions of
space,¹⁰³⁸ time,¹⁰³⁹ and frequency.¹⁰⁴⁰ Although design constraints
including material choice, physical restrictions, and operational
environmental limitations are unavoidable, there is otherwise great
design flexibility in the type of photonic metasurface utilized. This
rich history of EM design, evolving through distinct theoretical,
computational, and now artificial intelligence (AI)-driven eras,
informs our understanding of the interplay between design free-
dom and practical limitations in photonic metasurfaces, opening
new avenues for exploration and optimization.

In stark contrast to the vast EM design tools available today, at
the end of the 19th century scientists relied on manual calculations to
perform EM analysis¹⁰⁴¹ using the newly established Maxwell equa-
tions of 1864 (Ref. 1042) and the refined set in 1873.¹⁰⁴³ The initial
focus was on finding closed form solutions to Maxwell's equations
and by 1910 there were a dozen or so including those of
Sommerfeld,¹⁰⁴⁴,¹⁰⁴⁵ Lord Rayleigh,¹⁰⁴⁶ Mie,¹⁰⁴⁷ Hondros,¹⁰⁴⁸
Nicholson,¹⁰⁴⁹ Debye,¹⁰⁵⁰ and Hondros and Debye.¹⁰⁵¹ During the
first half of the 20th century more complex EM problems were solved
by generalization of these explicit solutions. EM problems involving
Sommerfeld integrals were more challenging to solve and benefited
from the development of asymptotic methods including stationary
phase and the method of steepest descent. Perturbation and

<footnote>\*willie.padilla@duke.edu</footnote>
<figure><img src="image_47.png" /><figcaption>**FIG. 39.** Illustration depicting a conventional scientific workflow wherein a human scientist acts as a reasoning engine, denoted (b); in this capacity, a scientist synthesizes empirical data, denoted (a), and existing scientific knowledge, in the form of models/theory, denoted (c), to develop novel scientific models, or perform design of physical systems, denoted (d). In many modern branches of science, empirical data are often collected through both physical experiment and electromagnetic (EM) simulations. The promise of deep learning is to automate the process of reasoning that is currently fulfilled by humans; however the reasoning capabilities of deep learning are currently far more limited. Note that, in general, the relationships between these components can be more complex than indicated here.</figcaption></figure>

approximation methods were also developed to solve more complex problems, and these are generally termed high frequency methods, to which geometrical and physical optics belongs.¹⁰⁵²

A key component of scientific inquiry is empirical observation, or *data*, which is crucial to guide theoretical investigations, and ultimately to validate them. Prior to the advent of computational EM solvers (CEMS), the primary means to collect data were the construction of physical experiments, which was relatively slow and costly, severely limiting progress. The computational era began in the 1940s with the development of the computer, but it was not until the 1960s, when various numerical methods emerged to solve Maxwell's equations, that the field really gained momentum. The CEMS developed were those of the finite-difference time-domain (FDTD) method,¹⁰⁵³ the method of moments (MOM),¹⁰⁵⁴ and the finite element method (FEM).¹⁰⁵⁵ In contrast to physical experiment, these approaches essentially allowed scientists to conduct experiments by writing and executing software, which is generally much faster and cheaper.

The evolution of CEMS has continued in earnest, and in the past few decades, we have witnessed the emergence of many accurate and efficient options. Today, we have commercial computational EM tools

available including full-wave approaches based on either integral equations (both time and frequency domain) or on differential equations (both time and frequency domain), or alternatively high frequency methods based on either EM fields (geometrical optics) or currents (physical optics). With modern commercial software, programs can be created much faster, and a single program can be used to collect data over a wide variety of physical settings.

The last decade has witnessed the emergence of deep learning (DL), which has received significant attention in the literature. In the modern academic milieu, it is reasonable to ask whether this attention is justified. To address this question, it is helpful to understand the fundamental role of DL in the EM scientific workflow: for example, how does DL relate to CEMS, human scientists, and the theoretical models they develop? DL is a specific type of artificial intelligence (AI) technique, and the long-term objective of AI is to enable machines to achieve human-like intelligence.¹⁰⁵⁶ In this regard DL occupies the same role in the scientific process as human scientists, which is to reason about empirical observation and existing knowledge (e.g., theory and models) to develop new knowledge, or to perform EM material or device design. This broad scientific workflow is illustrated in Fig. 39. The prospect of automating the role of humans—or even enhancing it—using machines is enticing, and understandably garners tremendous attention.

While DL represents an important step toward automated reasoning, at present, DL methods are quite limited in their reasoning capabilities. In their current form, the term “Deep Learning” most often refers to a collection of machine learning (ML) methods—namely, deep neural networks (DNNs)—which are broadly designed construct models that *approximate* the relationships present in observed data, e.g., CEMS.¹⁰⁵⁷,¹⁰⁵⁸ This capability is illustrated and described in Fig. 40. DL does, therefore, automate an important reasoning task that has historically been accomplished by human scientists—building models of the physical world from data. However, DNNs do not efficiently integrate knowledge of the physical world, or perform deductive reasoning. Nonetheless, DL does present a significant advance over prior machine learning approaches in nanophotonic metasurfaces and has enabled the analysis of far more complex materials than previously possible,¹⁰⁵⁹⁻¹⁰⁶² justifying some of the attention they have received.

**Figure 40** depicts the current state of EM design approaches, including (a) theory, (b) computational EMs (simulation), and (c) deep learning. A theoretical approach strives to represent observations mathematically, aiming for closed-form solutions or, less favorably, accurate approximations. In EMs, Maxwell's equations are derived from the experimental observations of Faraday, Ampère, Ørsted, Gauss, and others, and verify the soundness of this method. We depict this theoretical approach in Fig. 40(a) where the open green symbols are points of observations which have given rise to the exact theory, $y = f(x)$, represented by the solid red curve. We next consider an EM design problem where the theory is not known, i.e., no exact solution exists. This may arise due to the complexity of the problem or the incorporation of some newly discovered material. Although no exact solution to our problem exists, CEMS, which are based on Maxwell's equations, are able to model the problem accurately and we are able to determine $f(x)$ at some discrete values, shown as the black squares in Fig. 40(b). However, since CEMS are limited by computer memory and calculation time, the problem is usually restricted to some range, denoted by $x_m$ in Fig. 40(b). The deep learning approach to EM design permits the

<figure><img src="image_48.png" /><figcaption>**FIG. 40.** Illustration depicting the relationship between DL, theoretical findings, and experimental results (e.g., CEMS). Given some EM system we often have experimental parameters $x$, and observed experimental results $y$. (a) The true relationship underlying these data (red curve), denoted $y = f(x)$, can sometimes be uncovered via observation (open green symbols) and theoretical deduction. (b) Physical experiments and CEMS (open black squares) can be used to sample this function, providing pairs of data ($x$, $y$) for specific settings of $x$. (c) Deep learning methods can utilize empirical data (open black squares) to automatically construct models, denoted $\tilde{f}$ and shown as the dashed blue curve, that accurately approximate the true underlying relationships (solid red curve) over the domain of available data.</figcaption></figure>
<figure><img src="image_49.png" /><figcaption>FIG. 41. Illustration of forward modeling (bottom) and inverse modeling (top) illustrated using a metasurface comprised of repeating 2 × 2 supercells of cylindrical resonators. Each resonator in the supercell has an adjustable geometric structure, denoted x in the main text. In forward modeling we build some mathematical model or algorithm, denoted $\hat{f}$, for computing the behavior of the system, denoted y, given some setting of x. One widely used method for forward modeling is computational simulation, however, in some cases we can dramatically accelerate forward modeling using a neural network. Using this acceleration, one can perform design by rapidly evaluating $\hat{y} = \hat{f}(x)$ for different candidate settings of x, and then comparing $\hat{y}$ to a desired system behavior, denoted $y^*$. In inverse modeling we build a model, termed $\hat{g}$, that can predict a setting of x that would yield a desired set of properties, denoted $y^*$. Deep Inverse Models are a class of DNNs that have been shown capable of directly making accurate inverse predictions $\hat{x}$ for complex metasurfaces.</figcaption></figure>

determination of a surrogate model from data (often CEMS data), denoted $\hat{f}(x)$ which approximates the true function $f(x)$. In Fig. 40(c) the open black symbols represent the same x points used in (b) for simulation, but now they are used in the deep learning workflow to determine the surrogate function $\hat{f}$, represented by the dashed blue curve. The surrogate function $\hat{f}(x)$ is significantly faster to evaluate than performing CEMS. Past results have demonstrated that one evaluation of $\hat{f}(x)$ compared to one numerical simulation is nearly $10^6$ times faster.

Next we discuss in greater detail where DL has been successful in EM materials, and the roadmap for realizing the greater potential for DL to automate human reasoning.

## 1. State of the art

**Where is deep learning currently most successful?** Although DL has been applied for a variety of purposes in EM materials, its primary success has been as an instrument to dramatically accelerate forward and inverse modeling. The goal of forward modeling is to uncover some, yet unknown, mathematical relationship of the form $y = f(x)$, that relates the parameters of an experiment, denoted $x$, to the corresponding measurements obtained from that experiment, denoted $y$. In the context of metasurfaces, $x$ might quantify parameters such as a metasurface geometry (e.g., radii and height of cylinders), its bulk construction material, or the illumination of the material; similarly, the observations encoded in $y$ might include electric/magnetic field measurements, or absorption,¹⁰⁶³ reflection,¹⁰⁶⁴ transmission¹⁰⁶⁵ spectra, measured at a specific frequencies. On the other hand, the goal of inverse modeling is to construct some model of the form $x = \hat{g}(y)$, wherein we can provide some desired experimental observation, and then have the model efficiently produce the experimental parameters that would give rise to those observations. These two basic objectives are illustrated in Fig. 41. We next explain how forward and inverse modeling have been utilized EM metasurfaces, to important effect.

**Forward modeling.** DL accelerates forward modeling by automating the process of uncovering $f$ from observational data. In practice DL is usually implemented with some type of DNN which is essentially a parameterized function, denoted $\hat{f}_0$. The parameters of the DNN are

denoted $\theta$ and they influence the relationship between the input to the DNN and its output, so that $\hat{f}_0$ can be made to accurately approximate a wide variety of complex functions depending upon the settings of $\theta$. To approximate a forward function, DNNs assume the availability of $N$ empirical observations of the form $D = (x_i, y_i)_{i=1}^N$, which is typically procured by CEMS in EM applications, and where $N$ may vary substantially across applications. The DNN is then “trained” to approximate the relationship between $x$ and $y$ by iteratively adjusting the parameters $\theta$ in a way that tends to shrink the error between predictions made by the DNN, denoted $\hat{y} = \hat{f}(x)$ and the true setting of $y$ associated with the input setting provided to the DNN, $x$. In other words, $\hat{f}(x_i) \approx y_i$ for most of the $x_i \in D$. Given a large enough DNN (i.e., enough adjustable parameters), it has been shown that a DNN can approximate any underlying function to an arbitrary degree of accuracy.¹⁰⁶⁶,¹⁰⁶⁷ In practice however, the accuracy of the model obtained may vary depending upon many factors (e.g., the available training data, $N$, or the complexity of the true underlying function that we wish to approximate).

Once the model $\hat{f}$ is constructed, it can be used to make accurate predictions for settings of $x$ that were not present in $D$. This process is essentially regression analysis, for which there are many techniques. However, in contrast to previous approaches, DL has been demonstrated to often produce much more accurate approximations for challenging high-dimensional systems that exhibit complex (e.g., highly non-linear) relationships between $x$ and $y$. Furthermore, and crucially, DL models have been shown in some cases to make predictions for novel settings of $x$ at a speed that is several orders of magnitude faster than CEMS, e.g., by a factor of nearly $10^6$. Therefore, using $\hat{f}(x)$, it is possible to rapidly explore the behavior of the system under study, which can aid scientists with theoretical investigations.

In addition to scientific investigation, forward modeling can also be used to accelerate design. In contrast to modeling, the process of design involves the manipulation of experimental parameters to achieve a certain system behavior. In other words, we wish to find some setting of $x$ that will result in a particular desired setting of $y$. A common design strategy involves using a forward model to evaluate different candidate settings of $x$ until one is found that produces the desired behavior. For advanced materials CEMS is often the only available means of evaluating the properties of arbitrary settings of $x$. CEMS is often computationally expensive however, creating a major bottleneck to design. In such cases DL can be utilized to build a fast

<table><caption>TABLE I. A list of three major components for construction of a deep neural network: Data, Architecture, and Loss. The NN architectures shown are some of the most commonly employed in photonic metasurfaces and include: residual neural network (RNN), convolutional neural network (CNN), multi-layer perceptron (MLP), graph neural network (GNN). Here the Transformer is a broad class of neural networks with similar applicability as MLPs, and Generative Models are NNs that can be trained to create novel instances of data that resemble those from some collection of data.</caption><thead><tr><th>Data</th><th>Architecture</th><th>Loss function</th></tr></thead><tbody><tr><td>Experiments</td><td>RNNs</td><td>Error metric, e.g. MSE</td></tr><tr><td rowspan="4">Simulations</td><td>CNNs</td><td>Differential equation</td></tr><tr><td>FCNs</td><td>Algebraic equation</td></tr><tr><td>GNNs</td><td>Complexity penalties</td></tr><tr><td>Transformers</td><td></td></tr></tbody></table>
<table><caption>TABLE II. Examples of specialized DNN approaches to perform inverse modeling, often for the purpose of design.</caption><thead><tr><th>Model</th><th>Type</th><th>Applications to AEM problems</th></tr></thead><tbody><tr><td>Tandem (TD)</td><td>Deterministic</td><td>Liu et al.,<sup>522</sup> Ma et al.,<sup>1068</sup> Gao et al.,<sup>1069</sup> and Hou et al.<sup>1070</sup><br>So et al.,<sup>1071</sup> Long et al.,<sup>1072</sup> He et al.,<sup>1073</sup> and Xu et al.<sup>1074</sup><br>Ashalley et al.,<sup>1075</sup> Mall et al.,<sup>1076</sup> and Pilozzi et al.<sup>1077</sup><br>Phan et al.,<sup>1078</sup> Singh et al.,<sup>1079</sup> and Malkiel et al.<sup>1080</sup></td></tr><tr><td>Mixture Manifold Networks (MMN)</td><td>Deterministic</td><td>Spell et al.<sup>1081</sup></td></tr><tr><td>Genetic Algorithm (GA)</td><td>Optimization</td><td>Zhang et al.,<sup>312</sup> Johnson and Rahmat-Samii,<sup>1082</sup> and Deng et al.<sup>1061</sup><br>Forestiere et al.<sup>1083</sup> and Li et al.<sup>1084</sup></td></tr><tr><td>Neural Adjoint (NA)</td><td>Optimization</td><td>Deng et al.<sup>1085</sup> and Peurifoy et al.<sup>1086</sup><br>Asano et al.<sup>678</sup> and Miyatake et al.<sup>1087</sup></td></tr><tr><td>Variational Auto-encoder (VAE)</td><td>Probabilistic</td><td>Ma et al.,<sup>1088,1089</sup> Qiu et al.,<sup>1090</sup> and Kudyshev et al.<sup>1091,1092</sup><br>Shi et al.,<sup>1093</sup> Liu et al.,<sup>1094</sup> and Kiarashinejad et al.<sup>1095</sup></td></tr><tr><td>Mixture Density Network (MDN)</td><td>Probabilistic</td><td>Unni et al.<sup>1096,1097</sup></td></tr><tr><td>Invertible Neural Networks (INN)</td><td>Probabilistic</td><td>Ren et al.<sup>1098</sup></td></tr><tr><td>Generative Adversarial Networks (GAN)</td><td>Probabilistic</td><td>Wen et al.<sup>1099</sup></td></tr></tbody></table>

approximation to $f$, permitting a more efficient and effective design. Consequently, DL represents a powerful tool to accelerate this type of forward design procedure.<sup>1100</sup>

State-of-the-art DL most often takes the form of DNNs, which have a large number of free design choices, sometimes called hyperparameters, that must be chosen to suit the particular system under study. Some of the most important choices are the data utilized for training, the model architecture (i.e., what functional form does $\hat{f}$ or $\hat{g}$ take), and the loss function (i.e., what criteria is used to identify the “best” model during training). Fortunately, there are a number of highly effective choices for these hyperparameters depending upon properties of the problem, of which some are listed in Table I. The number of neural network hyperparameters is vast and we do not attempt to include an exhaustive list here. For the reader interested in learning more we refer to reviews on these topics.<sup>1101-1104</sup>

**Inverse modeling.** Despite the acceleration provided by DL, the forward design procedure can quickly become intractable for advanced photonic metasurfaces with many free design parameters (e.g., $x$ is high-dimensional), or the number of candidate settings for each design parameter is large. Another powerful capability of DNNs is so-called inverse modeling, whereby a DNN utilizes data to build a model of the form $x = \hat{g}(y)$ with the property that $f(\hat{g}(y)) = y$; in other words, the DNN approximates the inverse of $f$. In principle then, a scientist can specify a precise set of desired EM properties, denoted $y^*$, and the DNN will produce the design that will yield those desired properties—a transformative capability depicted in the top path of Fig. 41. This approach essentially uses DL for regression, similar to modeling and forward design, except with the roles of $x$ and $y$ reversed. However, this reversal of roles introduces some significant technical difficulties for regression, such as the so-called non-uniqueness problem, where several different settings of $x$ can produce the same desired setting of $y^*$. Non-uniqueness causes conventional DL methods to fail, however, substantial recent work has focused upon specialized DL models—termed deep inverse models (DIMs)—that can overcome these challenges. Recent studies have demonstrated that DIMs can be trained

successfully to rapidly identify designs that achieve customized settings of $y^*$, and do so for highly complex EM material systems.<sup>1098,1105</sup>

There are various DIMs which have been applied to photonic metasurfaces and in Table II we list some of these models, their type, and references to associated studies. DIMs fall into three basic categories—namely *Deterministic*, *Optimization*, and *Probabilistic*. Broadly, each class of model represents a trade-off between computational cost (e.g., speed) and design accuracy. Optimization approaches initialize a large population of randomly chosen candidate solutions and then sequentially improve all of them in parallel—a computationally intensive process, but one that often achieves the highest design accuracies. Like Optimization methods, Probabilistic DIMs typically initialize a large set of candidate solutions, but rather than incrementally optimizing them, they try to make the quality of these initial estimates higher by learning a probability distribution of where good solutions are likely to exist, given some target $y^*$. These approaches are relatively fast because they avoid sequential processing, but they typically yield less-accurate solutions because their solutions are still randomly drawn. Deterministic methods attempt to learn a direct mapping from $y^*$ to the appropriate design that will yield $y^*$. Consequently these approaches require just one evaluation of a DNN (e.g., for Tandem), or just a few evaluations in parallel (e.g., MMN), making them computationally fast. There are also hybrid methods that combine different types of DIMs, or combine DIMs with traditional optimization approaches (see review in Deng et al.<sup>1061</sup>) Further information about DIMs for photonic metasurfaces can be found in recent reviews of the topic.<sup>177,1098,1106</sup>

## 2. Future directions and outlook

Although DL has become a powerful tool to accelerate EM modeling and design, the long-term potential of DL is far greater. The future promise of DL, and AI more broadly, is to enhance and automate the role of human reasoning within the scientific process, as shown in Fig. 39. Currently however the reasoning capabilities of the most successful DL approaches—namely DNNs—are relatively limited
compared to humans. Although there are several limitations of DNNs compared to humans, we highlight two specific limitations that are especially influential: (i) world modeling, and (ii) deductive reasoning.

*World modeling* refers here to the ability of an intelligent agent (e.g., a DNN) to represent and utilize knowledge of objects in the world, and their known relationships (e.g., physical laws), when reasoning. For example, if we wish to build a model relating the position of a vehicle on a road, *y*, as a function of time, *x*, it is helpful to know that the vehicle is powered by a combustion engine, thereby limiting its maximum speed and rates of acceleration. This knowledge immediately and dramatically reduces the number of potential explanatory models, yet building a DNN that can represent and utilize such information when training (i.e., building a model) still requires substantial time of a human expert, and often still doesn't fully leverage this information. Consequently DNNs often lack world models altogether, or they are simplistic and limited in scope.

*Deductive reasoning* refers here to the ability of an intelligent agent to use existing knowledge to draw valid conclusions about the world via logic. In the context of physics, deductive reasoning might include, for example, the use of mathematics to deduce closed-form solutions to Maxwell's equations under specific conditions. Another example may be that field measurements from a EM material must obey Maxwell's equations, and therefore any exact model of the system [e.g., of the form $y=f(x)$] must also obey Maxwell's equations. Although DNNs sometimes incorporate the results of human deductive reasoning (e.g., so-called physics-informed neural networks, or Lorentzian Neural Networks) they do not explicitly perform deductive reasoning. Like world modeling, deductive reasoning can be used to dramatically reduce the number of valid models of observed data, and also ensure certain levels of generality.

**Near-term research directions.** The most acute consequence of DNN's limited reasoning capabilities is that it limits their *generalization* capabilities; that is, their ability to make inferences (e.g., predictions) for scenarios (e.g., inputs) that were not present in their training data. Consequently, DNNs typically require substantial quantities of training data to build accurate models, which is problematic because these data are typically generated by expensive and/or time-consuming CEMS, and it must be generated independently for each unique problem (i.e., physical system) under consideration. Furthermore, the quantity of data required to achieve a particular level of modeling accuracy grows rapidly with respect to the complexity of the problem under consideration (e.g., the dimensionality of *x* and *y*, and the non-linearity of *f*).1107 At the same time, CEMS is also often more computationally intensive for such problems, limiting the quantity of data that can be obtained. Therefore, data inefficiency represents a major limitation of DL, undermining its value even in applications where it has been highly successful. Importantly it also makes DL completely impracticable for many important open problems in EM, such as those involving multi-physics or multi-scale problems.1108 A variety of research is under way to improve the data efficiency of DNNs. These approaches are relatively superficial, in the sense that they do not significantly improve fundamental reasoning capabilities of DNNs, but they are also often more tractable, and have demonstrated some effectiveness:

*   **Physics-guided learning:** These approaches try to utilize known physical laws or principles to guide the interpretation of empirical data, ultimately reducing the amount of data needed. For

example, so-called physics-informed neural networks (PINNs) encourage any DL model to be consistent with Maxwell's equations.1109

*   **Transfer and incremental learning:** These approaches attempt to build upon knowledge gained from one problem—perhaps a problem with substantial quantities of data—and leverage it to aid the construction of a model for a new problem. Transfer learning is widely used in the broader DL community, and it has also recently been shown to dramatically reduce the quantity of data needed to solve novel metasurface problems as well.1110,1111

*   **Active learning:** These approaches attempt to identify data points that are most “informative,” in the sense that the model can improve its accuracy faster, and thereby use less total training data.1112

*   **Modeling ill-posed problems:** This research area is unique to inverse modeling, where non-uniqueness is a major obstacle to building inverse models. This is an active and important area of research for ultimately building effective deep inverse models.1098

**Long-term research directions.** Long-term DL research should aim to endow DL methods with more sophisticated reasoning capabilities, such as general World Modeling and Deduction. Achieving these goals represents a tremendous challenge; however, some progress is being made within the AI research community. Here we list three promising areas of current research that may endow DL with these two more sophisticated reasoning capabilities.

*   **Large language models (LLMs):**1113,1114 LLMs have recently demonstrated the capability of representing human language, and using it to perform some basic reasoning. If successful, LLMs represent one way to provide DNNs with knowledge about physical systems, and potentially perform deductive reasoning.

*   **Neuro-symbolic models:**1115–1117 These approaches attempt to combine the existing capabilities of DNNs (e.g., data-driven models) with powerful existing symbolic reasoning engines, thereby providing DNNs with the capability to apply deductive reasoning.

*   **Causal models**1118,1119 These approaches attempt to use deductive reasoning to help DNNs to distinguish between spurious correlations and causal relationships present in data.

Although we highlight these areas of DL research, it is noteworthy that there are many important active areas of DL research (e.g., interpretability)

**Outlook.** Since its emergence roughly a decade ago, DL has had a significant impact on many fields, including EM metasurfaces. Among their successes, DNNs have been especially useful for utilizing empirical data (in particular from CEMS) to build computationally efficient approximate models that can dramatically accelerate the modeling and design of complex photonic metasurfaces. The capabilities of DNNs are quite robust for these tasks, and they have quickly become pervasive within the EM research community. Consequently, we expect DL, serving in this capacity, to become an important tool for modeling and design of EM materials, alongside other tools such as theory and CEMS. Despite their impact however, DL represents a relatively rudimentary form of reasoning about empirical data. The long-term potential of DL, insofar as it is an AI technique, is tremendous. The long-term objective of AI is endow machines with human-level intelligence. Currently, DL fundamentally lacks important reasoning capabilities of
humans, for example, explicit world modeling and deductive reasoning. If DL or other AI techniques can be endowed with these more sophisticated forms of reasoning, it may result in much greater advances than has yet been provided by DL. Furthermore, DL possesses some advantages over humans; specifically, DL possesses substantially greater computation speed and memory capacity. Consequently, upon possession of human-like reasoning capabilities, it may also dramatically accelerate the rate at which such reasoning is done, considerably accelerating science. Although it is unclear how quickly DL and AI will progress, important inroads have been made toward tackling some of the aforementioned reasoning limitations it possesses (e.g., with LLMs, and neuro-symbolic models). Given the tremendous potential impact of advances in AI, and the commensurate investment in its development, we anticipate that it will be an increasingly important area of research within EM metasurfaces.

## ACKNOWLEDGMENTS

W.J.P. and J.M.M. acknowledge support from the Department of Energy under U.S. Department of Energy (DOE) (No. DESC0014372).

# XIX. COMPLEX-VALUED SINGULARITIES FOR METASURFACE DESIGNS

Aloke Jana, Zijin Yang, Rémi Colom, Qinghua Song, and Patrice Genevet

patrice.genevet@mines.edu

Optical singularities are points at which a physical quantity, such as the phase, polarization, polarization ellipticity, etc., is not defined, but around which this same physical quantity can take almost any value. This latter property is very beneficial for designing photonic devices with desired optical response. Moreover, singularities are ubiquitous in nanophotonics as they generically occur when varying at least two optical or geometrical parameters of the system. Here, we present approaches to design nanophotonic building blocks with desired optical responses via singularity position engineering. In particular, we show that the control of the positions of zeros and poles of optical response functions (e.g., reflection, transmission coefficient or scattering, and Jones matrices), the so-called singularities, can be used to design optical metasurfaces with desired amplitude and phase behaviors. Additionally, breaking the in-plane symmetries of the metasurface is shown to achieve a phase singularity for only one circular polarization state, while its orthogonal counterpart remains unmodulated, thereby opening new avenues for the information multiplexing with nanophotonic designs. Sharing similarities and inspiration from established works in signal processing and control theory, it is expected that the control of the position of complex-valued singularities in the “pole-zero complex-plane,” for example relying on symmetries, would provide new design methodologies in photonics.

## A. Introduction

Optical singularities refer to points or regions, either in the physical or parameter space, where certain optical or electromagnetic properties become undefined within the vicinity of steep gradients, as illustrated in Fig. 42. In ray optics, singularities manifest as points or

<figure><img src="image_50.png" /><figcaption>**FIG. 42.** Optical singularities and their topological features. (a) Depiction of a catacaustic formed by a circle with the source positioned on the circle. (b) Vortex beam associated with the phase singularity of the complex scalar wave, (c) polarization singularity of the vector beam, (d) optical skyrrelevant lattice, (e) exceptional points: depiction of the topological features of non-Hermitian systems linked with the degeneracy of the eigenvalue spectrum. Inset illustrates the real and imaginary parts of the eigenvalue spectrum in the parametric space. (f) Scattering singularities featuring a zero-pole pair of the complex reflection/transmission coefficients. Inset displays the phase profile of reflection/transmission coefficients, along with a simple illustration of the reflection/transmission zero condition of the scattering processes. Figure (a) reproduced with permission from “Catacaustics, Resultants and Kissing Conics — Greg Egan,” https://www.gregegan.net/SCIENCE/Catacaustics/Catacaustics.html (2012). Copyright Author(s);¹¹⁵³ (c) reproduced with permission from Copyright Hindawi Publishing Corporation;¹¹²⁹ (d) reproduced with permission from Tsesse et al., Science **361**, 993 (2018). Copyright 2018 AAAS;¹¹³⁰ (e) reproduced with permission from Miri and Alu, Science **363**, eaar7709 (2019). Copyright 2019 AAAS;³⁹ (e) reproduced with permission from Ozdemir et al., Nat. Mater. **18**, 783 (2019). Copyright 2019 Nature Publishing Group;¹³⁶ and (f) reproduced with permission from Colom et al., Laser and Photonics Rev. **17**(6), 2200976 (2023). Copyright 2023 Wiley-VCH.¹¹⁵⁴ Copyright remains with the respective publishers.</figcaption></figure>
lines of extreme light intensity (caustics) known as optical catastrophes. Meanwhile, phase singularities can be found in complex scalar wave fields as a form of vortices, wavefront dislocation, or nodal lines.<sup>1120-1124</sup> This ensemble of phenomena belongs to the realm of Singular optics. One very familiar example of a singular optics phenomenon is the existence of Laguerre Gaussian beams having zero intensity lines along the optical axis and a helical phase front associated with an intrinsic orbital angular momentum.<sup>1125-1127</sup> Incorporating the vector nature of EM waves gives rise to additional polarization singularities.<sup>1128</sup> These polarization singularities can manifest as singularities in one or more parameters of the polarization ellipse, such as the azimuthal angle and ellipticity angle, resulting in features like C points and L lines.<sup>1129</sup> Skyrmions are the stable topological defects of the three-dimensional vector field that exhibit promising applications in the domain of magnetic storage and spintronics. Similar possibilities, i.e., optical information, transfer, and storage arise in the realm of the photonic system by utilizing the topological properties of three-dimensional EM waves.<sup>1130</sup> Over the decades, singularities in optical systems have been studied extensively in different scenarios starting from the seminal work "Dislocation in wave trains" by Berry and Ney to "Observation of the optical mobius strips."<sup>1123,1131-1135</sup> In addition, recent advancements of spectral singularities associated with non-Hermitian systems have led down to fertile ground, revolutionizing various fundamental and technological aspects.<sup>39,1136-1139</sup> It's well known that a closed and conservative system governed by Hermitian Hamiltonian exhibits real eigenvalues with orthogonal eigenvectors due to its unitary evolution nature. Instead, open and non-conservative systems, described by non-Hermitian Hamiltonians, exhibit complex eigenvalues and non-orthogonal eigenvectors. Non-Hermitian systems possess many features with no counterpart in Hermitian systems. For example, exceptional points, a new type of degeneracies, may occur in non-Hermitian systems. They correspond to the degeneracy of two or more eigenvalues along with the coalescence of their associated eigenvectors. These EPs are linked with sharp transitions in the eigenvalue spectrum, capable of dramatically altering the response of the system. Exceptional points are very important in parity-time symmetric systems as they are associated with spontaneous symmetry breaking. They thus separate regions of the parameter space where the Hamiltonian spectrum is real from regions where the eigen-frequencies occur in complex conjugated pairs.<sup>112</sup> Interestingly, parity-time (PT) symmetric systems can display entirely real eigenvalue spectra despite being governed by non-Hermitian Hamiltonian. While EPs are initially discovered in the context of quantum mechanical Hamiltonians, realizing EPs experimentally in quantum systems poses significant challenges. However, optical and photonic systems have emerged as promising platforms for exploring non-Hermitian physics. This is primarily due to the controllable introduction of optical gain and loss—non-conservative elements in these systems. Given their design versatility, metasurfaces (MSs) represent a convenient approach to study and validate non-Hermitian concepts in photonics. MSs are two-dimensional artificial structures consisting of subwavelength unit cells which exhibit peculiar electromagnetic responses, not achievable in natural materials. Most importantly, the response of MSs building blocks, i.e., nanostructures, are characterized by the scattering matrices, which are inherently non-Hermitian. The non-Hermitian nature of metasurfaces greatly influences their scattering properties. This requires to look at the properties of coefficients of the scattering,

reflection or transmission coefficients in the complex-frequency plane. Phase singularities of these coefficients such as zeros and poles play a fundamental role. In addition to the position of single zeros, two or more of these quantities can merge in the complex plane, potentially exhibiting EP behavior. This characteristic also holds true for the conventional Jones matrix, which describes changes in the polarization state between input and output channels. Therefore, the non-homicity of the scattering/Jones matrices provides a ubiquitous tool to control and manipulate the properties of EM waves. By tuning the geometry and changing the composition of the individual unit cells, anomalous light propagation or scattering effect, including unidirectional propagation, polarization and diffraction control, maximum helical dichroism, ultrasensitive sensing and biosensors, imaging, and surface waves, can be studied.<sup>5,880,1140-1151</sup> Interestingly, even the most basic zeros of the system response functions, such as the complex reflection/transmission function, also exhibit phase singularities associated with non-trivial topological charge of $C = \frac{1}{2\pi} \oint_{\gamma} \nabla(\phi)d\omega$, where $\phi$ is the response phase function and the integration is carried out along a path $\gamma$ enclosing the singular zero response point in two-dimensional parameter space. It turns out that the latter reflection and transmission zeros play a leading, and often neglected, role in the design of photonic systems, including MSs.<sup>1152</sup>

This roadmap article focuses on these recent developments, drawing some of the perspectives and implications offered by the complex-valued zeros and EPs singularities of transfer functions for the design of innovative photonic systems. This manuscript is organized into three sections. In Sec. XIX C, we begin by introducing the basic notion of zeros in electromagnetic quantities and elucidate the impact of the position of singularities in the complex frequency plane, such as zeros and poles, emphasizing the influence of symmetries on altering system responses. We discuss two examples of MSs designs relying on the complex-valued singularity approach in reflection and transmission. In Sec. XIX D, we address a new wavefront engineering mechanism—topology-protected 2π phase accumulation for a specific polarization channel by encircling EPs in the parameter space. Furthermore, the construction of EP pairs is harnessed to achieve uni-directional full-color vectorial holography employing a wavelength-multiplexing scheme. Finally, in Sec. XIX E, we delve into the potential future prospects of complex-valued singularities particularly higher-order ones. We address a novel phase excursion mechanism associated with the zero mode conversion in parametric space and the designing of ingenious toolboxes that would exploit the control theory methodology for nanophotonics design.

## B. Designing metasurfaces according to the position of complex-plane zero singularities

MSs may be considered as linear systems which transform inputs into outputs. These transformations are mediated by transfer functions or matrices which generally depend on the frequencies and directions of the input signals. Depending on the situation and the type of MS to be designed, the transfer function can be the reflection or transmission coefficient or even the *S*− or the Jones *J*− matrix. Designing a MS often amounts to looking for the right set of geometrical and optical parameters allowing one to obtain the desired amplitude and phase of the transfer function for a given frequency or an interval of frequencies. One can therefore take inspiration from other fields concerned with the design of linear systems such as electronics or control
<figure><img src="image_51.png" /><figcaption>FIG. 43. (a). Behaviors of the amplitude and phase of a function possessing a zero and a pole in the complex frequency plane (top) and on the real axis (bottom). (b). Symmetries of transmissionless (top) and reflectionless states (bottom). (c). Silicon nanocone MS possessing a reflection zero and pole complex conjugated and whose amplitude remains large and constant while the phase varies from 0 to 2π. (d). (Top) trajectories of a couple of transmission zeros of a Sb₂S₃ nanocube MS as the height of the nanocubes h is varied. At two degeneracies of the two zeros, there is a transition between a situation where zeros are real and one where they are complex conjugated. (bottom) At h = 160 nm, one zero is in the upper part of the complex plane and almost complex-conjugated of a pole leading to a high amplitude and a variation of the phase between 0 and 2π for the transmission coefficient. Figure (b) reproduced with permission from Colom et al., Laser and Photonics Rev. 17(6), 2200976 (2023). Copyright 2023 Wiley-VCH;¹¹⁵⁴ (c) reproduced with permission from Mikheeva et al., Optica 10, 1287 (2023). Copyright 2023 OSA;¹⁴⁰ (d) reproduced with permission from Colom et al., Laser and Photonics Rev. 17(6), 2200976 (2023). Copyright 2023 Wiley-VCH.¹¹⁵⁴ Copyright remains with the respective publishers.</figcaption></figure>

engineering to describe and design MSs.¹¹⁵⁵ In these fields, it is common to look at the analytic continuation of response functions $F(\omega)$ of the systems in the complex-frequency plane. As discussed above, the necessity of analyzing the response at complex frequencies is rooted in the non-Hermitian nature of scattering processes used for the design of metasurfaces, particularly the resonant ones that suffer from scattering losses, and also sometime, of additional absorption losses.

Of particular importance are the zeros and poles of the response functions which generally occur in the complex-frequency plane. As displayed in Fig. 43(a), simple zeros and poles are a given type of topological defects at which the phase is singular and around which the phase behaves like a vortex with a $2\pi$ winding counterclockwise around zeros and clockwise around poles (more generally it winds $2n\pi$ times for nth order phase singularities). Their topological charge C, as defined earlier, is, respectively, +1 and -1 when choosing a contour $\gamma$ encircling the singularity in the complex frequency plane. In this study, the phase winding direction adheres to the right-hand rule, where $2\pi$

increase in the phase is defined as counterclockwise. The well-defined sense of rotation of the phase comes from the fact that we consider analytic functions.¹¹²⁶,¹¹⁵⁶ These phase singularities bring fundamental deeper topological aspects to the study and design of MSs. In this section, we will discuss approaches to shape the optical response based on the control of the relative positions of these phase singularities in the complex-frequency plane. We will also discuss the relationship between zeros of transmission, reflection, and S-matrix coefficients and the symmetries of MSs. Finally, we will provide two examples of MS designs that rely on this approach, one in reflection and another in transmission.

The analytic continuation in the complex-frequency plane of response functions $F(\omega)$ can be expanded according to the Weierstrass expansion¹¹⁵⁷,¹¹⁵⁸ depending on both zeros $\omega_{z,n}$ and poles $\omega_{p,n}$ of the response functions,

$$F(\omega) = A \exp(iB\omega) \prod_n \frac{\omega - \omega_{z,n}}{\omega - \omega_{p,n}}. \qquad (9)$$
It is well known that $\omega_{p,n}$ generally occur in the lower part of the complex plane [for an exp($-i\omega t$) time convention] in passive photonic systems because of causality. However, a special occurs for the bound state in the continuum (BIC), poles can occur at real frequencies. To preserve energy conservation, this can only occur when they merge with a zero.¹⁴¹

Zeros may on the contrary occur at any position in the whole complex plane. Depending on the positions of zeros, the variations of the response functions display different behaviors as a function of the input frequency. To illustrate this, let us first consider the case of a single zero-pole pair in the expansion (9) with $A=1$ and $B=0$, i.e. $F(\omega) = \frac{\omega - \omega_{RZ}}{\omega - \omega_p}$. It has been shown in Ref. 1156, that when the zero is in the lower part of the complex plane, there is generally a dip in the amplitude of the response function and the phase varies only on a small phase interval (usually smaller than $\pi$). When the zero is on the real axis, the response function cancels out for a certain frequency and as a result, there is a $\pi$ phase jump. Eventually, as illustrated in Fig. 43(a), when there is a zero in the upper part of the complex plane, the phase spans the full $[0, 2\pi]$ interval as the frequency is varied passing across the line separating the zero-pole pair. This property is better evidenced when the phase is unwrapped. Moreover, if the zero becomes complex-conjugate of the pole, one obtains a similar response as a Gires-Tournois etalon,⁴⁰,¹¹⁵⁹ i.e., a constant and unitary amplitude alongside with a $2\pi$ phase variation. This shows that one can use the position of zeros and poles of response functions as a design methodology to address MS with the desired response in both amplitude and phase.

To this end, it is necessary to have a good understanding of the rules constraining the positions of zeros in the complex-frequency plane. These constraints are related to the symmetries of the MS. It is for example well known that time-reversal symmetry imposes that zeros of the S-matrix coefficients are complex-conjugated of poles for lossless passive systems. Let us consider the symmetries governing the positions of transmission and reflection zeros. These symmetries are illustrated in Fig. 43(b) for normal incidence for transmission zeros (top) and reflection zeros (bottom). Delving into the scenario of a passive, lossless MS with a transmission zero at the frequency $\omega_{TZ}$, the time-reversal transformation leaves the MS invariant while causing a reversal of momenta in the incident and reflected waves. Consequently, the time-reversed counterpart of a transmissionless state is also a transmissionless state. As time-reversal implies complex conjugation of frequency, one can conclude that a lossless MS possessing a transmission zero $\omega_{TZ}$ also possesses a transmission zero at the frequency $\omega_{TZ}^*$. As a consequence, illuminating this device at normal incidence should display transmission zeros at real frequencies or pairs of complex-conjugated transmission zeros. This is reminiscent of eigenvalues of PT-symmetric Hamiltonian in the unbroken and broken phases.¹¹⁶⁰

Regarding reflection zeros, their positions depend not only on the time-reversal symmetry but also on the mirror symmetry with respect to the median plane of the MS which will be called Parity (P-symmetry in the following discussion. Let us consider a lossless MS possessing a reflection zero $\omega_{RZ}$. Time-reversal transformation amounts to reversing the momenta of the incident and transmitted waves, and as a result, this MS also possesses a reflection zero for the light impinging from the other side at the frequency $\omega_{RZ}^*$. If the MS also preserves P-symmetry, it possesses

two reflection zeros for light impinging from the same side at the frequencies $\omega_{RZ}$ and $\omega_{RZ}^*$. Just like transmission zeros, reflection zeros are consequently either real or occur in complex-conjugated pairs. On the other hand, when P-symmetry is broken, symmetry constrained are lifted and $\omega_{RZ}$ can be complex-valued. In fact, by breaking the MS P-symmetry, one can control the coupling coefficients of its modes to the top $\gamma_{Top}$ and the bottom $\gamma_{Bottom}$. From the temporal coupled mode theory, one finds that a zero and pole of the reflection coefficient of a lossless MS illuminated from the top with a resonance frequency whose real part is $\omega_0$ take the following expressions:⁴⁰

$$ \begin{aligned} \omega_{RZ} &= \omega_0 + i\gamma_{Top} - i\gamma_{Bottom}, \\ \omega_{P} &= \omega_0 - i\gamma_{Top} - i\gamma_{Bottom}. \end{aligned} \qquad (10) $$

It is thus clear that P-symmetry breaking allows control of the imaginary part of reflection zeros and their relative positions to the poles.

Let us illustrate the usefulness of this insight for designing MSs with a desired optical response. To create optical components such as meta-deflectors or metalenses, it is necessary to use MSs whose response function amplitude is large to ensure a good efficiency of the component. It is also important that the phase of the response functions spans the full $[0, 2\pi]$ to be able to generate arbitrary phase gradients. As explained before, this condition on the phase requires to design structures to possess a zero in the upper part of the complex plane while the condition on the amplitude requires that this zero is complex conjugated of the pole. Let us show how to fulfill these conditions for a MS working in reflection and another working in transmission.

For the reflection case, it appears from the expression of reflection zeros and poles (10) that it is necessary to cancel out the coupling of the mode to the bottom channel to obtain complex conjugated zeros and poles. This condition is quite intuitive as the energy transmitted can be considered as loss altering the reflection-based MS efficiency. To achieve this condition, we considered in Ref. 40 a MS consisting of an array of silicon 600 nm high nanopillar with a period of 800 nm embedded in a medium with a homogeneous refractive index equal to 1.5. The symmetry was broken by reducing the diameter at the top $D_{Top}$ of the pillar compared to the diameter $D_{Bottom}$ at the bottom which is equal to 500 nm leading to a truncated cone geometry as shown in the top left of Fig. 43(c). For a ratio $D_{Top}/D_{Bottom} = 0.84$, $\gamma_{Bottom}$ is almost zero for the resonance we consider and the zero and pole are consequently almost complex conjugated as can be seen in the right panels of Fig. 43(c). As a result, the amplitude of the reflection coefficient is almost unitary and the phase spans the full $2\pi$ interval as shown in the bottom left figure of Fig. 43(c).

The design of MS in transmission for normal incidence fulfilling the condition of zero-pole pairs which are complex conjugated is more challenging since the position of zeros does not depend on the spatial symmetries of the structure (as long as there is no polarization conversion). Yet, a solution exists for MSs working in the regime where zeros are complex-conjugated. This is in fact the underlying mechanism describing the physics of dielectric Huygens MSs.³⁶²,¹¹⁵⁴ To illustrate this, we consider a MS consisting of an array of nanocubes of $Sb_2S_3$. The top figures in Fig. 43(d) show the trajectories of the real and imaginary parts of a couple of transmission zeros as the height of the nanodisk is varied. It appears that there are two distinct heights ($\approx 140$ and
≈180 nm) at which transmission zeros are degenerated, this represents a specific type of exceptional points, previously dubbed scattering EPs.¹¹⁵⁴,¹¹⁶¹ In the height interval between these two degeneracies, the two zeros are complex conjugated while they are real outside this interval. As a result, one zero gets almost complex conjugated of a pole. The amplitude thus remains large while the phase varies from 0 to 2π as seen in the bottom plots in Fig. 43(d). Interestingly and unexpectedly, the nature of the 2π-phase modulation as a function of the real frequency is due to the crossing of the branch-cut connecting only one zero-pole pair, and it is not due to the sum of two resonances, as previously stated. The full phase modulation occurring here results from a spontaneous symmetry-breaking behavior, as previously discussed in the framework of topological photonics.¹¹⁶²,¹¹⁶³ As shown herein,

avoided crossing in non-Hermitian systems involves pairs of complex singularities. This behavior contrasts drastically with the more conventional avoided crossing of resonances in Hermitian systems.

### C. Asymmetric full 2π phase engineering based on exceptional points in metasurfaces

EPs may manifest in more complex-valued linear operators that represent input-output conversion relationships, with the scattering matrix and Jones matrix being the most typical among them. Taking the second-order Jones matrix as an example, particularly under the circular polarization basis, EPs correspond to conditions for which only one of the off diagonal terms cancel out while the other remains

<figure><img src="image_52.png" /><figcaption>**FIG. 44.** Exceptional points manifested in Jones matrix systems for asymmetric phase engineering. (a) Schematic diagram and results of achieving topologically protected 2π phase by encircling exceptional points. (b) Construction of a pair of exceptional points with orthogonal polarization state using mirror-symmetry strategy. (c) Illustration of asymmetric full-color vectorial holography using exceptional point pairs. (a) Reproduced with permission from Song *et al.*, Science **373**, 1133 (2021). Copyright 2021 AAAS;³⁷ (b) reproduced with permission from Yang *et al.*, Nat. Commun. **15**, 232 (2024). Copyright 2024 Nature Publishing Group;¹¹⁶⁵ and (c) reproduced with permission from Yang *et al.*, Nano Lett., **24**, 844–851 (2024). Copyright 2024 ACS.¹¹⁶⁸ Copyright remains with the respective publishers.</figcaption></figure>
not null, i.e., $T_{+-} = 0$ and $T_{-+}$ different from 0. This happens when the two eigen-polarizations of the Jones matrix coalesce and both correspond to the right or left handed circular polarization, where + and − represent left circular polarization (LCP) to right circular polarization (RCP), respectively. This signifies the suppression of the conversion channel from LCP to RCP, resulting in a zero or singularity point in the $T_{+-}$ channel. There is consequently a reflection or transmission zero for one conversion channel but not the other which will be useful for information multiplexing. The construction and observation of such EPs have been presented in previous papers.¹¹³¹,¹¹⁶⁴ Building upon these works, some of the authors of the present paper have implemented a new type of functional device by harnessing the robust topological phase, namely, the exceptional topological (ET) phase, by encircling the EPs. Using a meta-structure as shown in Fig. 44(a), by manipulating the dimensions of the structure in two orthogonal directions, we can attain desired x- and y-polarized resonances, facilitating the construction of an EP of the Jones matrix associated with the cancelation of $T_{+-}$ channel. Importantly, upon achieving a closed loop encircling the EPs in the parameter space, a phase accumulation of $\varphi = 2\pi$ will occur, which is attributed to the phase vortex existing around the singularity in the complex plane. Considering any point $\mathbf{R}$ on the complex plane, there is

$$P(\mathbf{R}) = Re(\mathbf{R}) + i Im(\mathbf{R}), \qquad (11)$$

where $P(\mathbf{R})$ represents the position function and $Re(\mathbf{R})$ and $Im(\mathbf{R})$ represent the real and imaginary parts of $\mathbf{R}$, respectively. Therefore, the phase of $P(\mathbf{R})$ can be described as

$$\varphi[P(\mathbf{R})] = \arctan \frac{Im(\mathbf{R})}{Re(\mathbf{R})}. \qquad (12)$$

When $\mathbf{R}$ moves continuously in parameter space from $P_1$ to $P_2$, the resulting accumulated phase can be given by the integral form:

$$\phi = \int_{P_1}^{P_2} \varphi[P(\mathbf{R})] d\varphi. \qquad (13)$$

According to Fig. 44(a) and Eq. (13), it is evident that a phase accumulation of $2\pi$ solely occurs when the loop encircles the zero point (or singularity) while failing to achieve full-phase modulation when the singular point moves away from the loop. The above analysis underscores the distinctive nature of channel $T_{+-}$, differing from its cross-polarization channel $T_{-+}$, as the latter's nearly constant amplitude distribution lacks the topological protection of the singular points. Leveraging this unique property, we successfully achieved the decoupling of two circularly polarized lights. Specifically, upon incidence of the RCP beam, a holographic image of the letter “C” appears at the designed 30° angle, while under LCP incidence, no corresponding image appears. Furthermore, by combining the Pancharatnam-Berry (PB) phase and ET phase, more sophisticated information multiplexing can be realized. This new phase modulation mechanism, together with resonant scattering, guided mode propagation phase, and orientation-related phase retardation, enriches the optical control methods, while its ease of implementation and compatibility with other phase-control mechanisms also open the avenue for the application of topological phenomena in the optical wavefront modulation. However, due to the inherent asymmetric behavior caused by the handedness of EP, current EP systems are confined to operation within a specific circular polarization state and lack the versatility to engineer arbitrary states. A PB phase,

equal but opposite in phase for both circular polarization states, is generally picked up by rotating polarization-converting optics. Operating specifically at an EP is shown to suppress one of the polarization-conversion channels, thus ensuring that PB-phase is encoded on only one of the polarization-conversion channels asymmetrically. To expand this specific asymmetric response of the EP photonic systems, researchers obtained the enantiomer $S^n$ of the L-shaped meta-structure (S) through a general mirror-symmetry strategy.¹¹⁶⁵ As shown in Fig. 44(b), upon $E_x$ incidence, the current distribution of the vertical arm in the meta-structure $S^n$ switches the flow direction relative to its counterpart in S, resulting in the degeneracy point of the eigenstates flipping from RCP to LCP, corresponding to the north and south poles of the Poincaré sphere, respectively. Therefore, in the parameter space defined by the two arms of the meta-structure ($L_1, L_3$), a pair of singularity points induced by EPs are obtained, manifesting as zero conversion coefficients for the two channels of $M_{LR}$ and $M_{RL}^m$ at ($L_1, L_3$) = (52 nm, 119 nm). Interestingly, these two singularities give rise to phase vortices characterized by opposite topological charges in the parameter space. By judicious arrangements of structures S and $S^n$ on the same MS, the resultant output polarization can be expressed as³⁷

$$\langle S|S \rangle = A_+ e^{i\varphi_+} |+\rangle + A_- e^{i\varphi_-} |-\rangle, \qquad (14)$$

where A represents amplitude and $\varphi$ represents the corresponding phase. According to Eq. (14), coupling a pair of EPs breaks the limitation of the circular output polarization state. By adjusting the amplitude ratio and phase difference of the two rows of meta-structures, a full-polarization state covering the Poincaré sphere can be constructed while maintaining the unidirectional asymmetric imaging properties of the EP system. In brief, this work extends the application of EP to any arbitrary polarization state, surpassing the inherent limitations of the asymmetric response brought by EP to the system, and enhancing the potential for topological wavefront shaping applications. Furthermore, based on the construction of the aforementioned EP pairs, they have achieved a unique unidirectional full-color vectorial hologram by introducing a wavelength-multiplexing strategy,¹¹⁶⁶ as shown in Fig. 44(c). Since the EPs can maintain the asymmetry of conversion in a wide wavelength range, the designed full-color vectorial hologram can still achieve optimized energy distribution. This asymmetric full-color holographic strategy based on EP pairs is expected to be widely used in more complex application systems such as information security and virtual reality.

## D. Concluding remarks and future perspective

In conclusion, we have presented methods to design MSs to achieve various functionalities by engineering the positions of zeros of R, T, S, and J via parameter tuning. Wavefront engineering relies on strategically placing complex-valued singularities in the parameter space while exploiting exceptional points (EPs) of the Jones matrix enabling asymmetric full $2\pi$ phase modulation and full-color vectorial holography. However, it is crucial to clarify a common misconception about the EPs of the Jones matrix as discussed in Sec. XIX D. Although they exhibit singular behavior in the parametric space, they are not the actual singular points of the linear system itself. Recent studies highlight that the degeneracy of 2 zeros of the S-matrix is associated with EPs, particularly those linked to coherent perfect absorption (CPA), with phase-sensitive single-mode absorption.¹¹⁶⁷ The attainment of resonant, absorbing, and CPA EPs in linear systems provides a versatile platform for controlling scattering responses. Additionally, the
potential for full 2π phase modulation with nearly 100% transmittance has been demonstrated using a glide-symmetric bilayer MS.¹¹⁶⁸ Recognizing the sensitivity of higher-order singularities to the perturbations, consideration of these in photonic systems can enhance sensitivity and robustness to their locations in the complex frequency plane. Notably, the geometric phase associated with spatial mode transformation would allow innovative implementation of the zeros in the mode conversion process. Akin to the exceptional ET phase in polarization channels, a novel phase-addressing mechanism can be realized by encircling the singular points of the conversion polarization coefficients in the parametric space. Drawing inspiration from the principles of control theory, wherein the placement of poles and zeros is crucial for shaping a system's behavior, including its stability, frequency response, and transient dynamic. This principle is similarly critical in nanophotonic device design, where the location of poles and zeros in the complex frequency plane profoundly influences properties such as transmission, reflection, and phase modulation, making it feasible to define the real-axis response and tailor nano-geometry based on the strategic placement of these topological singularities. Adopting control theory methodologies in photonics could thus pave the way to a new paradigm for designing nanophotonics devices.

## XX. MULTIPOLAR DESIGN METHODS FOR METASURFACES

Carsten Rockstuhl, Karim Achouri, Andrey B. Evlyukhin, Ulrich Lemmer, and Ivan Fernandez-Corbaton

carsten.rockstuhl@kit.edu

An efficient, fast, and reliable framework for the theoretical description and the numerical analysis of the optical properties of metasurfaces is vital to their exploration from a fundamental science perspective and for their integration into applications. Here, we review one such framework that resides on the broader framework of a semi-analytical scattering theory. In that description, the incident and scattered fields are expanded into a multipolar series, and their relation is expressed by a matrix-vector product where the matrix is called the transition or T-matrix. Besides summarizing the state-of-the-art of the method, we elaborate on current research questions and oversee future developments.

### A. Introduction

Computational means are nowadays indispensable for exploring the optical properties of metasurfaces.⁹⁴,¹¹⁶⁹ On the one hand, we can design them in silico according to predefined needs, achieving an optical response on-demand before any fabrication. In addition to being essential for the design, computational means are useful for obtaining insights into aspects notoriously challenging to measure. An example of the latter would be the vectorial details of the near-field distribution, an important quantity for nonlinear metasurfaces or those used for light management in solar cells. On the other hand, we can simulate the optical response of fabricated devices to explain measured quantities. Particularly for complicated metasurfaces, e.g., stacked ones, the measured response often does not correspond to what had been predicted before because of imperfections in the fabrication

process. Here, simulations might help disclose the deviations and suggest what to improve in a further fabrication round.

The numerical analysis of optical metasurfaces is particularly successful because multiple numerical tools that solve the governing equations, i.e., the Maxwell equations, have been developed and equipped with suitable constitutive relations. These tools are mainly derived from the general field of numerical mathematics and rely on discretizing the Maxwell equations by finite differences or finite elements in the time or frequency domain.¹¹⁷⁰,¹¹⁷¹ Also, for strictly periodic metasurfaces, various grating algorithms to describe diffractive optical elements and photonic crystals¹¹⁷² were developed. However, while indisputably helpful, the direct numerical solution of Maxwell equations often complicates obtaining insights into the governing physics of metasurfaces. Especially for metasurfaces, assembled from what is frequently called the meta-atoms and analyzed in terms of the multipolar response induced in such meta-atoms, a discussion that is closer to scattering theory is suggested to be often superior.

This Perspective article aims to give an overview of the recent developments of such a framework for the multipolar description of metasurfaces. The framework tremendously relies on the notion of scattering theory.¹¹⁷³ Initially, it puts the description of the optical response of the individual meta-atom into focus. It derives from that a set of insights into the optical functionality of entire metasurfaces. We emphasize that the optical response of the individual meta-atom is typically expressed by a T-matrix or generalized polarizabilities.¹¹⁷⁴ While both expressions can be made interchangeable, they are different. For example, the T-matrix relies on a problem formulation in spherical coordinates, while the polarizabilities use Cartesian coordinates. Here, we start from the multipolar description using the T-matrix and establish the link to polarizabilities at a later stage.

### B. State of the art

#### 1. Theoretical framework of a multipolar description

The multipolar description of a metasurface has to start with the multipolar description of its meta-atom. In the framework of scattering theory in the frequency domain, we initially seek elementary solutions to Maxwell equations that arise from their separability in spherical coordinates.¹¹⁷⁵ These vector spherical harmonics expand the incident and the scattered field from the meta-atom as

$$ \mathbf{E}^{\text{inc}}(\mathbf{r}, \omega) = \sum_{j=1}^{\infty} \sum_{m=-j}^{j} q_{jm}^{e}(\omega) \mathbf{N}_{jm}^{(1)}(\mathbf{r}, \omega) + q_{jm}^{m}(\omega) \mathbf{M}_{jm}^{(1)}(\mathbf{r}, \omega), \quad (15) $$

$$ \mathbf{E}^{\text{sca}}(\mathbf{r}, \omega) = \sum_{j=1}^{\infty} \sum_{m=-j}^{j} b_{jm}^{e}(\omega) \mathbf{N}_{jm}^{(3)}(\mathbf{r}, \omega) + b_{jm}^{m}(\omega) \mathbf{M}_{jm}^{(3)}(\mathbf{r}, \omega), $$

where $\mathbf{N}_{jm}^{(n)}(\mathbf{r}, \omega)$ and $\mathbf{M}_{jm}^{(n)}(\mathbf{r}, \omega)$ are the vector spherical harmonics. These modes have well-defined frequency $\omega$, angular momentum squared $j$, and angular momentum along a specific axis $m$. Physically, the fields with $j=1$ are dipolar fields, with $j=2$ quadrupolar fields, etc. We distinguish modes with electric parity $\mathbf{N}$ and magnetic parity $\mathbf{M}$. Finally, we distinguish between those vector spherical harmonics expanding the incident and the scattered field using the index $n=1$ and $n=3$, respectively. Instead of parity, a helicity base can also be chosen, especially if chiral light-matter interactions are at stake. The
relation between the concatenated amplitudes expanding the incident field **q** and scattered field **b** from a meta-atom in a linear response is written as

$$ \mathbf{b}(\omega) = \mathbf{T}(\omega)\mathbf{q}(\omega), \qquad (16) $$

where the T-matrix T(ω) had been introduced, which entirely characterizes the optical response of the meta-atom.¹¹⁷⁶ For a sphere, the T-matrix can be built analytically as a diagonal matrix containing the Mie-coefficients. For an arbitrary shaped object, the T-matrix can be obtained numerically by multiple means.¹¹⁷⁷ They all rely on exciting the meta-atom with different illuminations, computing and expanding the scattered fields into vector spherical harmonics, and reconstructing the T-matrix that links all the different incident and scattered field amplitudes determined numerically.

In the following, we discuss metasurfaces made from a periodic arrangement of identical meta-atoms.¹¹⁷⁸ Disordered metasurfaces or spatially inhomogeneous metasurfaces can also be considered. However, they are covered in other perspective articles in this APL Roadmap. While considering an arbitrary chosen meta-atom to be the central one, we can write the field illuminating the particle as a superposition of the external illumination and the scattered field from all the other particles. By dropping the frequency argument, this can be written as

$$ \mathbf{b} = \mathbf{T} \left( \mathbf{q} + \sum_{\mathbf{R}} \mathbf{C}^{(3)}(-\mathbf{R}) \mathbf{b}_{\mathbf{R}} \right), \qquad (17) $$

where C(3)(-R) are the lattice coupling matrices that convert the scattered field from meta-atom i into an incident field on the central meta-atom. The prime in the sum indicates that the central particle is

<figure><img src="image_53.png" /><figcaption>FIG. 45. Illustration of some phenomena and applications explored in the context of metasurfaces and which can be discussed extremely convenient using the multipolar description. Starting at midnight and going clockwise, the phenomena are lattice resonances (figure adapted from Ref. 1179). Then we have diffusers as discussed in this road-map also in the context of disordered metasurfaces (figure adapted from Ref. 1180). Next, we have holograms. Then, we have bound-states-in-the-continuum (BIC) that require a multipolar analysis in their design (figure adapted from Ref. 1181). As the following phenomena, we have mentioned here pseudochirality and then angular scattering. Finally, we also mention as a possible domain of exploration molecular metasurfaces. At last, we have as an example the wide field of Kerker effects (figure adapted from Ref. 1182). All figures reproduced (adapted) with permission from the respective Journals and copyright remains with the original publishers.</figcaption></figure>
excluded. Considering plane wave illumination characterized by a transverse wave vector that expresses the phase variation of the incident field across the metasurface and exploiting the Bloch theorem, we can write the scattered field amplitude from each meta-atom forming the metasurface as

$$ \mathbf{b}_R = e^{i\mathbf{k}_{\parallel} \cdot \mathbf{R}} \mathbf{b}. \qquad (18) $$

Substituting this expression into Eq. (17) and solving it for the scattered field amplitude vector, we obtain

$$ \mathbf{b} = \left( \mathbb{I} - \mathbf{T} \sum_{\mathbf{R}} \mathbf{C}^{(3)}(-\mathbf{R}) e^{i\mathbf{k}_{\parallel} \cdot \mathbf{R}} \right)^{-1} \mathbf{T} \mathbf{q}. \qquad (19) $$

This equation is intuitive and says that the external illumination induces multipole moments that expand the scattered field. The response is mediated by the T-matrix of the meta-atom renormalized by the interaction with all the other meta-atoms forming the lattice. Their amplitudes can be expressed by summing up the scattered field from all the particles forming the lattice and projecting it onto the plane wave constituting the discrete diffraction orders. These amplitudes of the diffraction orders are the primary quantity of the solution to the problem, and many further properties of interest can be extracted from them.

Examples are reflectance and transmittance spectra or their phases, polarization-dependent quantities, or near- and far-fields, to mention only a few selected ones. This is at that heart of many optical phenomena that can be explored. We mention some of these phenomena and applications thereof in Fig. 45 and elaborate on some of them further below.

Also, scattering matrices can be defined that express reflection and transmission upon illuminating the metasurface from both directions. These scattering matrices constitute the primary object of obtaining the response from stacked metasurfaces. Layer-doubling techniques are very efficient for getting the response from stacked metasurfaces.1183

## 2. Numerical implementation

The aforementioned theoretical framework is highly susceptible to numerical implementation. The algebraic formulation reduces the computational complexity to solving matrix-vector products. With the analytical availability of T-matrices for spheres, it is no surprise that the first numerical implementations of the above-sketched framework were done nearly 30 years ago.1184 At that time, the algorithms were mainly used to explore photonic crystals made from spheres and gratings from high-index dielectric or metallic spheres, which we call metasurfaces today.

<figure><img src="image_54.png" /><figcaption>FIG. 46. A computational example to highlight the strength of the multipolar description of metasurfaces. (a) Geometry of the considered example: The layer sequence corresponds to ITO (orange), amorphous silicon p-doped (green) and intrinsic (red), and crystalline silicon (purple). A metasurface made from titanium dioxide disks with a height of 100 nm and a radius of 150 nm lies on top of that layer stack. The lattice and the period would be subject to variations. (b) Example of a hexagonal lattice with a square root of the unit cell area equal to $\sqrt{A} = 520$ nm, we show the simulated reflectance at normal incidence and the reference reflectance from a sample without the metasurface. The figure also shows an AM1.5D spectrum. (c) Reflectance of the sample as a function of the wavelength and the lattice size for two different types of lattices, i.e., once for a hexagonal lattice and once for a rectangular lattice with an aspect ratio of 1.3. (d) When weighting the reflectance with the AM1.5D spectrum and integrating between 300 and 1100 nm, we obtain a measure for the current loss due to reflectance. That quantity is shown for different lattices and different lattice sizes. An optimum can be identified from that FIG. Thanks to the effective multipolar description, the computation is extremely fast. The T-matrix of the disk needs to be calculated only once for each frequency. And the reflectance analysis from the lattice are eventually matrix-vector-products of a small size, which can be done very efficiently. FIG adapted from Ref. 1178. All figures reproduced (adapted) with permission from the respective Journals and copyright remains with the original publishers.</figcaption></figure>
With the renewed interest in expressing the response of meta-atoms as building blocks of metamaterials on a multipolar basis, various groups developed numerical tools to study the optical response from metasurfaces. The key aspect of all implementations is the evaluation of the lattice couplings in Eq. (19). Here, Ewald summation techniques are exploited that split the badly converging sums into a part evaluated in real space and a part evaluated in reciprocal space.¹¹⁸⁵ Current development revolves around the computations of lattice sums along different dimensions and considering multiple meta-atoms per unit cell.¹¹⁸⁶ That latter would be needed when studying bound states in the continuum or topological metasurfaces. The proper consideration of particles closely space above substrates or the consideration of particles that penetrate their respective circumscribing spheres is also explored.¹¹⁸⁷ Finally, it remains to be mentioned that various codes have been made publicly available.¹¹⁸⁸⁻¹¹⁹⁰ A representative analysis of the optical response of a metasurface for an application as a light management structure is discussed in Fig. 46. All the simulations shown in that figure were computed within 4.6 h on an ordinary desktop PC. That includes the computation of the T-matrix of the considered disk. Performing the same simulations using a finite-element code would have taken 2.6 years. The analysis highlights the strength of the multipolar description of the metasurface.¹¹⁷⁸

However, it needs to be clearly said that the approach is less preferred in specific situations. For example, when only the specific response to a specific excitation of a fixed metasurface is needed, the additional overhead of computing the T-matrix first that allows us to express the response to any given illumination is not justified. This should be the largest weakness of the method. Moreover, other methods, including both full-wave Maxwell solvers and semi-analytical techniques, are equally applicable to model metasurfaces. Common full-wave solvers utilize finite element or finite difference methods in either time or frequency domain. Alternatively, semi-analytical approaches such as the rigorous coupled wave analysis (RCWA) can be effectively used. Initially developed for dielectric gratings, the RCWA has been adapted for photonic crystals, metamaterials, and metasurfaces with a strictly periodic arrangement of identical unit cells.

In the RCWA method, one first solves for Bloch-periodic eigenmodes of Maxwell's equations within a metasurface assumed to be invariant in the principal propagation direction. These eigenmodes are represented in a plane wave basis, with interface conditions subsequently applied to determine the amplitudes of modes excited within the metasurface and the adjacent spaces. This method is particularly suited for metasurfaces with discrete heights, treating the metasurface as a single layer. However, for metasurfaces with continuously varying surface profiles, the profile must be segmented into different layers. While RCWA rapidly predicts the optical response of a metasurface, it does not provide extensive physical insights. Here, the multipolar scattering framework is preferable, as it highlights the influence of individual unit cells on the optical response.

Moreover, although RCWA is efficient for metasurfaces with discrete heights, the analysis becomes more complex and time-consuming for continuously varying profiles or metasurfaces composed of free-form particles. Each layer requires solving a dedicated eigenvalue problem, and if the unit cell varies within a supercell, the expansion must include many plane waves, slowing the computation significantly. The complexity of these computations typically scales cubically with the number of plane waves, imposing practical limits on

the size of the unit cell that can be effectively analyzed. So, if a simulation of the response from a reasonably sized, uniformly heighted unit cell of a metasurface is required, a method like RCWA is a viable option. For all other configurations, the T-matrix-based scattering approach is more appropriate, provided that the additional computational effort to determine the T-matrix of the scatterer is worth it.

### 3. Analytical investigations

While the multipolar description constitutes a powerful tool to explore the functionality of metasurfaces numerically, it unfolds its forte when considered from an analytical perspective.¹¹⁹¹ Two aspects are important.

First, explicit analytical expressions for many quantities of interest can be derived under some restricted settings. When isotropic scatterers are considered, i.e., spheres, the T-matrices reduce to a diagonal one, and only a few coefficients are non-zero. Moreover, while theoretically, the T-matrix is infinite in size and numerically truncated to a sufficiently large size to ensure convergence behavior, it is analytically feasible to consider only a few lowest-order multipolar orders. Furthermore, high-symmetry lattices, e.g., square and subwavelength lattices, can be assumed to lower the complexity of the setting. Then, eventually, the system gets reduced in complexity, and explicit analytical equations can be expressed for most quantities of interest. For example, the reflection and transmission amplitude at normal or oblique incidence can be easily expressed. This dramatically simplifies the analysis and allows for parametric investigations of effects such as bound states in the continuum, lattice resonances, polarization filters, or diffraction gratings. We wish to make clear that the formulated assumptions are only made for convenience to reach manageable expressions. It ultimately is up to the convenience and efforts of the person who uses the method to decide up to which degree of complexity they would like to work analytically. One can easily derive analytical expressions for any observable quantity by lifting the assumption on isotropic scatterer and a square and subwavelength lattice, and one can also consider arbitrary higher multipolar orders, it is just that the expressions get longer and are not that handy anymore. In cases with too many unknowns, the problem can always be treated numerically. It is just that closed-form analytical expressions are no longer reasonable. We leave it up to the researcher to decide up to which order they wish to consider on analytical grounds.

Second, the entries of the T-matrix themselves do not attain arbitrary values. For example, for spheres, they can be parametrically expressed using a single angle $\theta_j$ bound between $-\pi/2$ and $\pi/2$ for lossless particles. For absorbing particles, a second angle $\theta'_j$ bound between 0 and $\pi/2$ parametrizes the dissipation.¹¹⁹² The possible value of a Mie coefficient reads then

$$b_j = \frac{1}{1 - i \tan \theta_j + \tan \theta'_j} \qquad (20)$$

This constitutes a significant simplification, as a desired response can be adjusted by identifying the Mie angles of all involved entries of the T-matrix. Combined with the analytical expressions, the necessary properties of the meta-atom in terms of Mie angles can be identified that provide a metasurface with a response on demand can be identified. Identifying suitable particles with the specified properties is a secondary problem that inverse-design techniques can solve.
C. Future directions and outlook

1. A description using polarizabilities

So far, the description has focused on an analysis using T-matrices that express the scattering response using basis functions in a spherical coordinate system. Sometimes, that is inconvenient, and a discussion in multipole moments linked to specific directions in a Cartesian coordinate system may be more intuitive. Then, one speaks of multipolar polarizabilities that are induced instead of the amplitudes of the vector spherical harmonics.¹¹⁹³ Instead of speaking of a T-matrix, we speak of the Cartesian multipole moments and polarizabilities. When care is taken to eliminate the trace and the redundant components of the Cartesian tensors, the change between Cartesian and spherical multipoles is merely a change of basis. It has no implication on the physics as long as no additional approximations are made.

Expressions exist to convert between these amplitudes of the spherical harmonics and the Cartesian multipole moments. For example, the two lowest-order electric coefficients can be converted as¹¹⁷⁴

$$
\begin{align}
\epsilon_d^{-1} \mathbf{p} &= \frac{\zeta_1}{ik_d^3} \mathbf{M}_{\text{Dipole}} \mathbf{b}_1^e, \tag{21} \\
k_d \epsilon_d^{-1} \mathbf{Q} &= \frac{\zeta_2}{ik_d^3} \mathbf{M}_{\text{Quadrupole}} \mathbf{b}_2^e,
\end{align}
$$

where the matrices that connect the electric multipole moments in a
spherical basis ($\mathbf{b}_1^e$ and $\mathbf{b}_w^e$) to a Cartesian basis ($\mathbf{p}$ and $\mathbf{Q}$) are

$$
\begin{align}
\mathbf{M}_{\text{Dipole}} &= \frac{1}{\sqrt{2}} \begin{bmatrix} 1 & 0 & -1 \\ -i & 0 & -i \\ 0 & \sqrt{2} & 0 \end{bmatrix}, \nonumber \\
\mathbf{M}_{\text{Quadrupole}} &= \frac{1}{\sqrt{2}} \begin{bmatrix} -i & 0 & 0 & 0 & i \\ 0 & -i & 0 & -i & 0 \\ 0 & 0 & \sqrt{2} & 0 & 0 \\ 0 & 1 & 0 & -1 & 0 \\ 1 & 0 & 0 & 0 & 1 \end{bmatrix}. \tag{22}
\end{align}
$$

Similar matrices also exist for higher-order terms.¹¹⁹¹ In these expres-
sions, $k_d$ is the wavenumber in the background medium with a permit-
tivity $\epsilon_d$ and $\zeta_j = \sqrt{(2j + 1)!\pi}$. Using that framework, T-matrices can
be converted into polarizabilities with the same simplicity. For exam-
ple, the electric dipolar part of the T-matrix $T_{ee}$ can be converted into
the electric dipolar polarizability matrix $\alpha_{ee}$ as

$$
\mathbf{T}_{ee} = - \frac{Z_d k_d^3 c}{6\pi i} \mathbf{M}_{\text{Dipole}} \boldsymbol{\alpha}_{ee} \mathbf{M}_{\text{Dipole}}^{-1}, \quad (23)
$$

with $Z_d$ the impedance of the background medium. Then, all the
details described above can be expressed in polarizabilities and
Cartesian multipole moments. Besides the numerical analysis, semi-
analytical expressions for observables that provide valuable insights
can also be obtained.

In the framework of the analytical treatment discussed above, for
example, for metasurfaces composed of identical isotropic meta-
atoms, the electric field coefficients of reflection *r* and transmission *t* at
the condition of the normal incidence of an external linear *x*-polarized
plane wave can be represented as¹¹⁹⁴

$$
\begin{align}
t &= 1 + \frac{ik_d}{E_0 2S_L \epsilon_0 \epsilon_d} \left( p_x + \frac{m_y}{c_d} - \frac{ik_d Q_{xz}}{6} - \frac{ik_d M_{yz}}{2c_d} + \dots \right), \tag{24} \\
r &= \frac{ik_d}{E_0 2S_L \epsilon_0 \epsilon_d} \left( p_x - \frac{m_y}{c_d} + \frac{ik_d Q_{xz}}{6} - \frac{ik_d M_{yz}}{2c_d} + \dots \right),
\end{align}
$$

where $E_0$ is the incident electric field at the center of meta-atoms, $S_L$ is the surface area of the metasurface elementary cell, and $c_d$ is the light speed in the surrounding medium. For the above equations, we assume that the metasurface plane coincides with the xy-plane of the Cartesian coordinate system, and the incident wave of the frequency $\omega$ is $(E_0, 0, 0)e^{i(k_d z - \omega t)}$. The considered Cartesian multipole moments of the meta-atoms ($\mathbf{p}$ and $\mathbf{m}$ are the electric and magnetic dipole moments, and $\mathbf{Q}$ and $\mathbf{M}$ are the electric and magnetic quadrupole moments, respectively) enter through the components that provide radiation in a direction perpendicular to the metasurface. Therefore, there is a useful opportunity to combine multipole moments with the same parity symmetry (with respect to inversion)¹¹⁹⁵ into one group and introduce the corresponding effective multipole $M^{(p)}$ and its effective polarizability $\alpha^{(p)}$, for the even (odd) multipole moments $p \to e$ ($p \to o$).¹¹⁹⁷ So that for the case (24) one can write

$$
t = 1 + \frac{ik_d}{2S_L} (\alpha^{(o)} + \alpha^{(e)}), \quad (25)
$$

$$
r = \frac{ik_d}{2S_L} (\alpha^{(o)} - \alpha^{(e)}), \qquad (26)
$$

where $\alpha^{(o)} = M^{(o)} / (E_0 \epsilon_0 \epsilon_d) = (p_x - ik_d M_{yz} / 2c_d + \dots) / (E_0 \epsilon_0 \epsilon_d)$ and $\alpha^{(e)} = M^{(e)} / (E_0 \epsilon_0 \epsilon_d) = (m_y / c_d - ik_d Q_{xz} / 6 + \dots) / (E_0 \epsilon_0 \epsilon_d)$. From the representation (25) and (26) one can see that the lattice anapole effect, when the light wave transmits through the metasurface without changes in amplitude and phase, is realized if $\alpha^{(o)} = 0$ and $\alpha^{(e)} = 0$ independently, that is due to destructive superposition between contributions of the multipole moments of the same inverse parity symmetry. In contrast, the lattice Kerker effect (only the suppression of reflection) is realized for $\alpha^{(o)} = \alpha^{(e)}$ due to the destructive superposition between the multipole moments with different inverse parity symmetry. The case of a bianisotropic response (when multipole coupling exists between the multipole moments of different inverse parity) can also be involved in this approach.¹¹⁹⁶,¹¹⁹⁷

Note that we have been talking about the effective polarizability of meta-atoms inside a metasurface so far. However, significant advantages can be obtained if the polarizability of single meta-atoms and the so-called lattice sums are considered.⁷⁰,¹¹⁹⁸,¹¹⁹⁹ Recently, a strategy has been proposed to tune a symmetry-protected bound state in the continuum in the desired spectral range by directly determining the period of a metasurface consisting of building blocks with known electric or magnetic dipole polarizabilities. This can be explicitly obtained from the point of intersection (coincidence) between the real part of the corresponding inverse dipole polarizability of an individual building block (meta-atom) and the real part of the corresponding lattice sum, taking into account inter-particle interaction in the metasurface. In particular, it is shown that to obtain a bound state in the continuum at the resonant wavelength $\lambda_R$ of an individual meta-atoms, the period of a metasurface with a square cell must be equal to $\approx 0.71\lambda_R$.¹²⁰⁰ In the case of including the resonant contribution of multipoles of various orders into the metasurface model, the approximation of effective polarizabilities and lattice sums makes it possible to study the lattice-induced
interaction between these multipoles in a semi-analytical form.¹¹⁹⁹ In some cases, such coupling can lead to the suppression of contributions from certain groups of multipoles and the implementation of anapole states or quasi-bound states in the continuum.

## 2. Multipolar interface conditions

In parallel to the elaboration of the T-matrix mentioned above method, an alternative metasurface modeling technique has been developed based on generalized sheet transition conditions (GSTC).¹²⁰¹⁻¹²⁰⁵ This modeling technique relies on the premise that a metasurface acts as a discontinuity that may be modeled as a homogeneous polarizable sheet whose response is described in terms of effective material parameters. For this purpose, the concept of polarizabilities is generally replaced by that of susceptibilities, as the latter is more commonly used to describe the response of a material from the perspective of a Cartesian lattice. The main difference between these two modeling techniques is that the T-matrix approach focuses on the scattering response of individual meta-atoms to build up to the response of the overall metasurface by including the coupling between the meta-atoms. In contrast, the GSTC method ignores the granular response of the meta-atoms and rather directly relates the fields interacting with the metasurface to its susceptibilities, which intrinsically includes the coupling between the meta-atoms.¹²⁰²⁻¹²⁰⁵

Fundamentally, the GSTC models how the fields and their spatial derivatives interact with a metasurface. This model is derived by expanding these fields, along with the metasurface polarization densities, in terms of a series of derivatives of the Dirac delta distribution. This rigorous foundation makes the GSTC ideal for modeling conventional dipolar metasurfaces and opens the door for including more complex interactions. As such, they now include second-order nonlinear polarizations,¹²⁰⁶ time-varying modulations,¹²⁰⁵ multipolar contributions, and nonlocal responses.¹²⁰⁷,¹²⁰⁸ Another strength of the GSTC is that they can be used for both synthesis, i.e., finding the metasurface susceptibilities in terms of specified fields, and for analysis, i.e., computing the fields scattered by a metasurface with known susceptibilities. For instance, it has been demonstrated using the GSTC that nonlocal responses are required to synthesize perfectly refractive metasurfaces.¹¹⁶⁹ Several numerical schemes have been developed for the analysis, such as finite-difference frequency- and time-domain techniques, a finite element method, and a spectral-domain integral-equation technique, that can compute the fields scattered by a metasurface under arbitrary illumination conditions.¹²⁰⁹

Finally, connecting the susceptibilities to the geometry of actual meta-atoms is typically achieved using equations that relate these susceptibilities to the metasurface full-wave simulated scattering parameters.¹²⁰⁵ While effective, this approach remains slow and cumbersome. For this reason, a methodology to connect the susceptibilities directly to the T-matrix of meta-atoms, which should significantly simplify the metasurface practical implementation procedure, is being actively developed. Note that it is currently already possible to map the dipolar polarizabilities of a metasurface with its corresponding dipolar susceptibilities.¹²⁰⁵

## 3. Further developments

Finally, we wish to reflect on selected advanced topics currently under development in the context of the multipolar description of

metasurfaces. Generally, the multipolar decomposition of the optical response of meta-atoms allows one to understand and engineer notable and beneficial effects, such as directional radiation. As understood by Kerker et al.,¹²¹⁰ the interference of the fields radiated by an electric dipole and a magnetic dipole can result in highly directional far-fields. The effect has been engineered and demonstrated in different frequency bands for samples under far-field illumination⁹⁶²,¹²¹¹ and for dipolar emitters near the nanostructures.¹²¹² Understanding these effects is a major challenge and many open questions remain, partially mentioned in this contribution to the roadmap but partially also in other contributions.

Moreover, the formalism introduced here cannot just be applied to artificially structured macroscopic materials but also to other (nano-) photonic systems. A major stream of efforts goes into developing metasurfaces with molecular building blocks.¹²¹³ From a gigantically large chemical space, the individual molecules can be chosen that are periodically arranged using different techniques, and materials can be provided that inherit their linear, bianisotropic, or nonlinear properties from the individual molecules. In a multi-scale effort, the framework described here can be used to express the properties of these materials from first principles and to study their interaction with light. To reach such a description, quantum-chemical tools, mostly time-dependent density functional theory, compute the T-matrix of an individual molecule or an aggregate. Afterward, periodically arranged molecules accessible, e.g., with surface-anchored metallic organic frameworks can be studied in detail. The detailed and accurate description of the molecules that surpasses a phenomenologically motivated description of a molecule as a polarizable two-level system provides unprecedented insights. Furthermore, hybrid systems can be explored where molecular and traditional scatterers are seamlessly blended.

Moreover, the metasurfaces can be homogenized in some cases, and effective properties can be assigned to a film so that the optical response from the metasurface and the effectively homogenous medium agree.¹²¹⁴ This unlocks novel possibilities when integrating metasurfaces into the design of more traditional optical elements. While the theory is well-established in dipole approximation, nonlocal metasurfaces are within reach and are currently being discussed. In these metasurfaces, higher-order multipolar components have a significant contribution. The complexity in their description advances, and more sophisticated properties are in reach. The nonlocal character, generally, leads to an angle dependent response, that can be used for image processing purposes. In addition, it unlocks the opportunity to multiplex the response depending on, e.g., the angle of incidence. Such a multiplexing can add enhanced functionality to the devices. Additionally, deriving generalized impedance, admittance, or even susceptibility tensors from the T-matrices or polarizabilities renormalized by the two-dimensional lattice interactions can be of practical benefit. It allows us to study and exploit the renormalization of the properties of the scatterer due to the lattice interaction and enables shaping the optical response actually not just by choosing a specific meta-atom, but rather by combining the meta-atom and the lattice. Exploiting the renormalization can lead to effects such as a giant magnetic dipolar moments, which may in turn find applications when perceiving impedance-matched optical devices with vanishing reflection losses.

Finally, significant efforts are currently in place to solve the inverse problem, i.e., instead of quantifying the optical response of a
given metasurface, we aim to identify the metasurface that provides an optical response on demand.¹¹⁰⁶ We have touched on that issue already in Sec. XXC3. Yet/However it can also be done with established methods from inverse design, such as topology optimization or Bayesian optimization, tools that effectively exploit gradient information, genetic algorithms, or particle swarm optimization. Moreover, techniques from the field of machine learning are explored. However, for individual design problems, the efforts in generating data used for training artificial neural networks (ANNs) usually exceed by far the efforts to solve the inverse problem directly. Yet/However independent of that, the community actively explores scenarios where using ANNs is appealing.

## D. Conclusions

This mini-review gave an overview of a highly flexible approach to theoretically and computationally explore the optical response of metasurfaces. It resides on the language of scattering theory, where initially, the interaction of light with a localized object, i.e., the meta-atom, is studied on semi-analytical grounds. That can be done with the T-matrix or the polarizability matrix, depending on whether a spherical or Cartesian coordinate system is considered. Thanks to the efficient evaluation of lattice sums, the interaction among the meta-atoms forming the metasurface is assessed, and the optical response in terms of reflection and transmission from metasurfaces can be evaluated. It needs to be emphasized that the computational approach is highly efficient, and the design and analysis of metasurfaces is tremendously accelerated. Yet/however it is also a framework that can be reduced to obtain handy analytical expressions to discuss numerical observables on analytical grounds. Multiple current developments were sketched.

What remains is to elaborate shortly on open issues. One aspect concerns densely packed metasurfaces made from asymmetric objects. The scattering formalism assumes that the circumscribing spheres of adjacent objects do not penetrate, as otherwise, the Rayleigh hypothesis breaks down. That restriction limits the types of metasurfaces that can be considered. The problem can be circumvented with concepts such as distributed T-matrices, where multiple origins are considered to expand the field into vector spherical harmonics, but the approach is yet premature. Moreover, combining the considered framework with methods from inverse design is ongoing, and we can expect much progress in this field. One aspect needs to be highlighted here, and it concerns the simulation of metasurfaces made from a spatially inhomogeneous arrangement of meta-atoms, which can be handled by considering rather large supercells. Spatially inhomogeneous metasurfaces are needed for many applications, e.g., meta-lenses, meta-holograms, and other optical meta-devices. The design of these devices usually follows a three-step procedure. First, the complex transmission function at each spatial location is identified by suitable means. Second, the optical response from a periodically arranged unit cell with varying geometrical parameters is simulated and stored in some lookup table. Last, the particles are fetched from the lookup and placed to the position of the unit cell where they shall provide the desired complex response. However, this approach does not consider the fact that the response of the meta-atom in an inhomogenous environment is not the same as that of the same meta-atom in a periodic arrangement. In consequence, the metasurface shows a degraded efficiency. Optimizing the spatial details of the unit cell by methods from the field of inverse design is potentially the key to reach devices of very high

efficiency, as these coupling effects need to be considered. The general framework described here can handle these effects, but efforts are needed to formulate everything in a differentiable form and to use it in an inverse design framework.

Furthermore, the combination of the optical properties of the metasurface with other physical phenomena would tremendously widen the range of use cases. One example would be photothermal behavior, where the dissipation of energy changes the temperature of the medium, which in turn changes the material properties, which changes the optical response. In general, many other effects could be considered to work toward switchable and tunable metasurfaces.

Finally, the extension of the formalism to treat nonlinear effects can safely be foreseen as a direction for future developments. Many of the other contributions in this roadmap will also serve as a guide for open issues that can be tackled with dedicated improvements to the framework sketched in this contribution.

## ACKNOWLEDGMENTS

C.R. and U.L. acknowledge support through the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy via the Excellence Cluster 3D Matter Made to Order (EXC-2082/1, Grant No. 390761711). A.B.E. acknowledges support through the Cluster of Excellence PhoenixD (EXC 2122, Project ID No. 390833453). I.F.-C. and C.R. acknowledge support through the CRC Waves: Analysis and Numerics (SFB 1173, Grant No. 258734477). K.A. acknowledges funding from the Swiss National Science Foundation (Project No. PZ00P2_193221). We all wish to acknowledge contributions to this work from our past and current group members, for which we are deeply grateful.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

S.A.S., R.F.O., and M.K. conceived and organized this roadmap, including selecting the section topics and lead authors. These three authors contributed equally and may claim first authorship. The following authors led individual sections (in order of the sections): Andrea Alu, Isabelle Staude, Rupert F. Oulton, Alasdair W. Clark, Maria Kafesaki, Anatoly V. Zayats, Stefan Maier, Xianzhong Chen, Andrea Di Falco, Maria V. Chekhova, Ren-Min Ma, Rachel Grange, Mohsen Rahmani, Andrei Faraon, R. Sapienza, Philippe Lalanne, Sebastian A. Schulz, Willie J. Padilla, Patrice Genevet and Carsten Rockstuhl. These authors supervised the writing of their relevant sections, selected their coauthors and determined the overall direction of their perspectives. All other authors contributed to the perspective articles under which they are listed. The author list is in order of appearance in the document.

Sebastian A. Schulz, Rupert. F. Oulton, Mitchell Kenney conceived and coordinated the roadmap together. All three act as joint corresponding and joint first authors for the full document.

**Sebastian A. Schulz:** Conceptualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Rupert F. Oulton:**
Conceptualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Mitchell Kenney:** Conceptualization (equal); Writing – original draft (equal); Writing – review & editing (equal). **Andrea Alù:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Isabelle Staude:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Ayesheh Bashiri:** Writing – original draft (equal); Writing – review & editing (equal). **Zlata Fedorova:** Writing – original draft (equal); Writing – review & editing (equal). **Radosław Kolkowski:** Writing – original draft (equal); Writing – review & editing (equal). **A. Femius Koenderink:** Writing – original draft (equal); Writing – review & editing (equal). **Xiaofei Xiao:** Writing – original draft (equal); Writing – review & editing (equal). **John Yang:** Writing – original draft (equal); Writing – review & editing (equal). **William J. Peveler:** Writing – original draft (equal); Writing – review & editing (equal). **Alasdair W. Clark:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **George Perrakis:** Writing – original draft (equal); Writing – review & editing (equal). **Anna C. Tasolamprou:** Writing – original draft (equal); Writing – review & editing (equal). **Maria Kafesaki:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Anastasiia Zaleska:** Writing – original draft (equal); Writing – review & editing (equal). **Wayne Dickson:** Writing – original draft (equal); Writing – review & editing (equal). **David Richards:** Writing – original draft (equal); Writing – review & editing (equal). **Anatoly Zayats:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Haoran Ren:** Writing – original draft (equal); Writing – review & editing (equal). **Yuri Kivshar:** Writing – original draft (equal); Writing – review & editing (equal). **Stefan Maier:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Xianzhong Chen:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Muhammad Afnan-Ansari:** Writing – original draft (equal); Writing – review & editing (equal). **Yuhui Gan:** Writing – original draft (equal); Writing – review & editing (equal). **Arseny Alexeev:** Writing – original draft (equal); Writing – review & editing (equal). **Thomas F. Krauss:** Writing – original draft (equal); Writing – review & editing (equal). **Andrea di Falco:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Sylvain D. Gennaro:** Writing – original draft (equal); Writing – review & editing (equal). **Tomás Santiago-Cruz:** Writing – original draft (equal); Writing – review & editing (equal). **Igal Brener:** Writing – original draft (equal); Writing – review & editing (equal). **Maria V. Chekhova:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Ren-Min Ma:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Viola V. Vogler-Neuling:** Writing – original draft (equal); Writing – review & editing (equal). **Helena C. Weigand:** Writing – original draft (equal); Writing – review & editing (equal). **Ülle-Linda Talts:** Writing – original draft (equal); Writing – review & editing (equal). **Irene Occhiodori:** Writing – original draft (equal); Writing – review & editing (equal). **Rachel Grange:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Mohsen Rahmani:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Lei Xu:** Writing – original

draft (equal); Writing – review & editing (equal). **S. M. Kamali:** Writing – original draft (equal); Writing – review & editing (equal). **E. Arbabi:** Writing – original draft (equal); Writing – review & editing (equal). **Andrei Faraon:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **A. C. Harwood:** Writing – original draft (equal); Writing – review & editing (equal). **S. Vezzoli:** Writing – original draft (equal); Writing – review & editing (equal). **Riccardo Sapienza:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Philippe Lalanne:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Alexandre Dmitriev:** Writing – original draft (equal); Writing – review & editing (equal). **Carsten Rockstuhl:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Alexander Sprafke:** Writing – original draft (equal); Writing – review & editing (equal). **Kevin Vynck:** Writing – original draft (equal); Writing – review & editing (equal). **Jeremy Upham:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **M. Zahirul Alam:** Writing – original draft (equal); Writing – review & editing (equal). **Israel De Leon:** Writing – original draft (equal); Writing – review & editing (equal). **Robert W. Boyd:** Writing – original draft (equal); Writing – review & editing (equal). **Willie J. Padilla:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Jordan M. Malof:** Writing – original draft (equal); Writing – review & editing (equal). **Aloke Jana:** Writing – original draft (equal); Writing – review & editing (equal). **Zijin Yang:** Writing – original draft (equal); Writing – review & editing (equal). **Rémi Colom:** Writing – original draft (equal); Writing – review & editing (equal). **Patrice Genevet:** Conceptualization (supporting); Writing – original draft (equal); Writing – review & editing (equal). **Karim Achouri:** Writing – original draft (equal); Writing – review & editing (equal). **Andrey B. Evlyukhin:** Writing – original draft (equal); Writing – review & editing (equal). **Ulrich Lemmer:** Writing – original draft (equal); Writing – review & editing (equal). **Ivan Fernandez-Corbaton:** Writing – original draft (equal); Writing – review & editing (equal).

DATA AVAILABILITY

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

REFERENCES

¹N. Meinzer, W. Barnes, and I. Hooper, Nat. Photonics **8**, 889 (2014).

²N. Yu and F. Capasso, Nat. Mater. **13**, 139 (2014).

³S. M. Choudhury, D. Wang, K. Chaudhuri, C. DeVault, A. V. Kildishev, A. Boltasseva, and V. M. Shalaev, Nanophotonics **7**, 959 (2018).

⁴A. I. Kuznetsov, M. L. Brongersma, J. Yao, M. K. Chen, U. Levy, D. P. Tsai, N. I. Zheludev, A. Faraon, A. Arbabi, N. Yu, D. Chanda, K. B. Crozier, A. V. Kildishev, H. Wang, J. K. W. Yang, J. G. Valentine, P. Genevet, J. A. Fan, O. D. Miller, A. Majumdar, J. E. Fröch, D. Brady, F. Heide, A. Veeraraghavan, N. Engheta, A. Alù, A. Polman, H. A. Atwater, P. Thureja, R. Paniagua-Dominguez, S. T. Ha, A. I. Barreda, J. A. Schuller, I. Staude, G. Grinblat, Y. Kivshar, S. Peana, S. F. Yelin, A. Senichev, V. M. Shalaev, S. Saha, A. Boltasseva, J. Rho, D. K. Oh, J. Kim, J. Park, R. Devlin, and R. A. Pala, ACS Photonics **11**, 816 (2024).

⁵N. Yu, P. Genevet, M. Kats, F. Aieta, J. Tetienne, F. Capasso, and Z. Gaburro, Science **334**, 333 (2011).
<footnote n="6">⁶X. Ni, N. K. Emani, A. V. Kildishev, A. Boltasseva, and V. M. Shalaev, *Science* **335**, 427 (2012).</footnote>

<footnote n="7">⁷B. A. Munk, *Frequency Selective Surfaces: Theory and Design* (Wiley, New York, 2000).</footnote>

<footnote n="8">⁸D. Pozar, S. Targonski, and H. Syrigos, *IEEE Trans. Antennas Propagat.* **45**, 287 (1997).</footnote>

<footnote n="9">⁹E. Hasman, V. Kleiner, G. Biener, and A. Niv, *Appl. Phys. Lett.* **82**, 328 (2003).</footnote>

<footnote n="10">¹⁰P. Lalanne, S. Astilean, P. Chavel, E. Cambril, and H. Launois, *Opt. Lett.* **23**, 1081 (1998).</footnote>

<footnote n="11">¹¹F. Monticone, N. M. Estakhri, and A. Alù, *Phys. Rev. Lett.* **110**, 203903 (2013).</footnote>

<footnote n="12">¹²C. Pfeiffer and A. Grbic, *Phys. Rev. Lett.* **110**, 197401 (2013).</footnote>

<footnote n="13">¹³H.-T. Chen, A. Taylor, and N. Yu, *Rep. Prog. Phys.* **79**, 076401 (2016).</footnote>

<footnote n="14">¹⁴N. Mohammadi Estakhri and A. Alù, *Phys. Rev. X* **6**, 041008 (2016).</footnote>

<footnote n="15">¹⁵D. N. Neshev and A. E. Miroshnichenko, *Nat. Photonics* **17**, 26 (2023).</footnote>

<footnote n="16">¹⁶A. Arbabi, Y. Horie, M. Bagheri, and A. Faraon, *Nat. Nanotechnol.* **10**, 937 (2015).</footnote>

<footnote n="17">¹⁷Metalenz, see https://metalenz.com/</footnote>

<footnote n="18">¹⁸Meta materials, Inc., see https://metamaterial.com/</footnote>

<footnote n="19">¹⁹Y. Liu, X. Liu, X. Mu, T. Hou, J. Xu, M. D. Renzo, and N. Al-Dhahir, *IEEE Commun. Surv. Tutorials* **23**, 1546 (2021).</footnote>

<footnote n="20">²⁰G. Zheng, H. Mühlenbernd, M. Kenney, G. Li, T. Zentgraf, and S. Zhang, *Nat. Nanotechnol.* **10**, 308 (2015).</footnote>

<footnote n="21">²¹E. Arbabi, A. Arbabi, S. M. Kamali, Y. Horie, and A. Faraon, *Optica* **3**, 628 (2016).</footnote>

<footnote n="22">²²Y. Wang, Q. Chen, W. Yang, Z. Ji, L. Jin, X. Ma, Q. Song, A. Boltasseva, J. Han, V. M. Shalaev, and S. Xiao, *Nat. Commun.* **15**, 5560 (2021).</footnote>

<footnote n="23">²³A. Epstein and G. V. Eleftheriades, *J. Opt. Soc. Am. B* **33**, A31 (2016).</footnote>

<footnote n="24">²⁴C. P. Jisha, S. Nolte, and A. Alberucci, *Laser Photonics Rev.* **15**, 2100003 (2021).</footnote>

<footnote n="25">²⁵F. Presutti and F. Monticone, *Optica* **7**, 624 (2020).</footnote>

<footnote n="26">²⁶D. A. B. Miller, *Science* **379**, 41 (2023).</footnote>

<footnote n="27">²⁷Y. Ra'di, D. L. Sounas, and A. Alù, *Phys. Rev. Lett.* **119**, 067404 (2017).</footnote>

<footnote n="28">²⁸Y. Raadi and A. Alu, *IEEE Photonics J.* **14**, 1–13 (2022).</footnote>

<footnote n="29">²⁹A. Overvig and A. Alù, *Laser Photonics Rev.* **16**, 2100633 (2022).</footnote>

<footnote n="30">³⁰A. Overvig and A. Alù, *Adv. Photonics* **3**, 026002 (2021).</footnote>

<footnote n="31">³¹S. C. Malek, A. C. Overvig, A. Alù, and N. Yu, *Light: Sci. Appl.* **11**, 246 (2022).</footnote>

<footnote n="32">³²H. Ren, G. Briere, X. Fang, P. Ni, R. Sawant, S. Héron, S. Chenot, S. Vézian, B. Damilano, V. Brändli, S. Maier, and P. Genevet, *Nat. Commun.* **10**, 2986 (2019).</footnote>

<footnote n="33">³³J.-H. Song, J. van de Groep, S. J. Kim, and M. L. Brongersma, *Nat. Nanotechnol.* **16**, 1224 (2021).</footnote>

<footnote n="34">³⁴E. Cortés, F. Wendisch, L. Sortino, A. Mancini, S. Ezendam, S. Saris, L. De, S. Menezes, A. Tittl, H. Ren, and S. Maier, *Chem. Rev.* **122**, 15082 (2022).</footnote>

<footnote n="35">³⁵G.-Y. Lee, J.-Y. Hong, S. Hwang, S. Moon, H. Kang, S. Jeon, H. Kim, J.-H. Jeong, and B. Lee, *Nat. Commun.* **9**, 4562 (2018).</footnote>

<footnote n="36">³⁶K. Wang, M. Chekhova, and Y. Kivshar, *Phys. Today* **75**(8), 38 (2022).</footnote>

<footnote n="37">³⁷Q. Song, M. Odeh, J. Zúñiga-Pérez, B. Kanté, and P. Genevet, *Science* **373**, 1133 (2021).</footnote>

<footnote n="38">³⁸A. Krasnok, D. Baranov, H. Li, M.-A. Miri, F. Monticone, and A. Alù, *Adv. Opt. Photonics* **11**, 892 (2019).</footnote>

<footnote n="39">³⁹M.-A. Miri and A. Alu, *Science* **363**, eaar7709 (2019).</footnote>

<footnote n="40">⁴⁰E. Mikheeva, R. Colom, K. Achouri, A. Overvig, F. Binkowski, J.-Y. Duboz, S. Cueff, S. Fan, S. Burger, A. Alù et al., *Optica* **10**, 1287 (2023).</footnote>

<footnote n="41">⁴¹H. Kwon, D. Sounas, A. Cordaro, A. Polman, and A. Alù, *Phys. Rev. Lett.* **121**, 173004 (2018).</footnote>

<footnote n="42">⁴²Y. Zhou, H. Zheng, I. I. Kravchenko, and J. Valentine, *Nat. Photonics* **14**, 316 (2020).</footnote>

<footnote n="43">⁴³M. Cotrufo, A. Arora, S. Singh, and A. Alù, *Nat. Commun.* **14**, 7078 (2023).</footnote>

<footnote n="44">⁴⁴A. Silva, F. Monticone, G. Castaldi, V. Galdi, A. Alù, and N. Engheta, *Science* **343**, 160 (2014).</footnote>

<footnote n="45">⁴⁵H. Zheng, Q. Liu, I. I. Kravchenko, X. Zhang, Y. Huo, and J. G. Valentine, “Multichannel meta-imagers for accelerating machine vision,” *Nat. Nanotechnol.* (unpublished) (2024).</footnote>

<footnote n="46">⁴⁶A. Cordaro, B. Edwards, V. Nikkhah, A. Alù, N. Engheta, and A. Polman, *Nat. Nanotechnol.* **18**, 365 (2023).</footnote>

<footnote n="47">⁴⁷A. C. Overvig, S. A. Mann, and A. Alù, *Phys. Rev. X* **11**, 021050 (2021).</footnote>

<footnote n="48">⁴⁸P. Iyer, R. Decrescent, Y. Mohtashami, G. Lheureux, N. Butakov, A. Alhassan, C. Weisbuch, S. Nakamura, S. Denbaars, and J. Schuller, *Nat. Photonics* **14**, 543 (2020).</footnote>

<footnote n="49">⁴⁹N. Muhammad, Y. Chen, C.-W. Qiu, and G. P. Wang, *Nano Lett.* **21**, 967 (2021).</footnote>

<footnote n="50">⁵⁰A. M. Shaltout, V. M. Shalaev, and M. L. Brongersma, *Science* **364**, eaat3100 (2019).</footnote>

<footnote n="51">⁵¹Y. Li, M. R. Krissnamurthi, W. Luo, A. K. Swan, X. Ling, and R. Paiella, *Opt. Mater. Express* **12**, 4528 (2022).</footnote>

<footnote n="52">⁵²J. Lee, M. Tymchenko, C. Argyropoulos, P.-Y. Chen, F. Lu, F. Demmerle, G. Boehm, M.-C. Amann, A. Alu, and M. A. Belkin, *Nature* **511**, 65 (2014).</footnote>

<footnote n="53">⁵³G. Li, S. Zhang, and T. Zentgraf, *Nat. Rev. Mater.* **2**, 17010 (2017).</footnote>

<footnote n="54">⁵⁴R. Kolkowski, S. Kovaios, and A. F. Koenderink, *Phys. Rev. Res.* **3**, 023185 (2021).</footnote>

<footnote n="55">⁵⁵S.-Q. Li, X. Xu, R. Maruthiyodan Veetil, V. Valuckas, R. Paniagua-Dominguez, and A. I. Kuznetsov, *Science* **364**, 1087 (2019).</footnote>

<footnote n="56">⁵⁶D. L. Sounas and A. Alù, *Nat. Photonics* **11**, 774 (2017).</footnote>

<footnote n="57">⁵⁷A. Vaskin, R. Kolkowski, A. F. Koenderink, and I. Staude, *Nanophotonics* **8**, 1151 (2019).</footnote>

<footnote n="58">⁵⁸S. Yuan, X. Qiu, C. Cui, L. Zhu, Y. Wang, Y. Li, J. Song, Q. Huang, and J. Xia, *ACS Nano* **11**, 10704 (2017).</footnote>

<footnote n="59">⁵⁹S. Liu, A. Vaskin, S. Addamane, B. Leung, M. Tsai, Y. Yang, P. Vabishchevich, G. Keeler, G. Wang, X. He, Y. Kim, N. Hartmann, H. Htoon, S. Doorn, M. Zilk, T. Pertsch, G. Balakrishnan, M. Sinclair, I. Staude, and I. Brener, *Nano Lett.* **18**, 6906 (2018).</footnote>

<footnote n="60">⁶⁰Y. Mohtashami, R. A. DeCrescent, L. K. Heki, P. P. Iyer, N. A. Butakov, M. S. Wong, A. Alhassan, W. J. Mitchell, S. Nakamura, S. P. DenBaars, and J. A. Schuller, *Nat. Commun.* **12**, 3591 (2021).</footnote>

<footnote n="61">⁶¹Y.-J. Lu, R. Sokhoyan, W.-H. Cheng, G. Kafaie Shirmanesh, A. R. Davoyan, R. A. Pala, K. Thyagarajan, and H. A. Atwater, *Nat. Commun.* **8**, 1631 (2017).</footnote>

<footnote n="62">⁶²C.-Y. Chang, C.-L. Yu, C.-A. Lin, H.-T. Lin, A. B. Lee, Z.-Z. Chen, L.-S. Lu, W.-H. Chang, H.-C. Kuo, and M.-H. Shih, *ACS Appl. Nano Mater.* **3**, 6855 (2020).</footnote>

<footnote n="63">⁶³L. Liu, W. Liu, F. Wang, H. Cheng, D.-Y. Choi, J. Tian, Y. Cai, and S. Chen, *Nano Lett.* **22**, 6342 (2022).</footnote>

<footnote n="64">⁶⁴W.-J. Joo, J. Kyoung, M. Esfandyarpour, S.-H. Lee, H. Koo, S. Song, Y.-N. Kwon, S. H. Song, J. C. Bae, A. Jo, M.-J. Kwon, S. H. Han, S.-H. Kim, S. Hwang, and M. L. Brongersma, *Science* **370**, 459 (2020).</footnote>

<footnote n="65">⁶⁵S. Cao, Y. Jin, H. Dong, T. Guo, J. He, and S. He, *J. Phys.: Mater.* **4**, 035001 (2021).</footnote>

<footnote n="66">⁶⁶G. Lozano, D. J. Louwers, S. R. Rodríguez, S. Murai, O. T. Jansen, M. A. Verschuuren, and J. Gómez Rivas, *Light: Sci. Appl.* **2**, e66 (2013).</footnote>

<footnote n="67">⁶⁷I. Staude, T. Pertsch, and Y. S. Kivshar, *ACS Photonics* **6**, 802 (2019).</footnote>

<footnote n="68">⁶⁸R. Mupparapu, T. Bucher, and I. Staude, *Adv. Phys.: X* **5**, 1734083 (2020).</footnote>

<footnote n="69">⁶⁹Y. Liu, S. C. Lau, W.-H. Cheng, A. Johnson, Q. Li, E. Simmerman, O. Karni, J. Hu, F. Liu, M. L. Brongersma, T. F. Heinz, and J. A. Dionne, *Nano Lett.* **23**, 6124 (2023).</footnote>

<footnote n="70">⁷⁰V. G. Kravets, A. V. Kabashin, W. L. Barnes, and A. N. Grigorenko, *Chem. Rev.* **118**, 5912 (2018).</footnote>

<footnote n="71">⁷¹S. I. Azzam and A. V. Kildishev, *Adv. Opt. Mater.* **9**, 2001469 (2021).</footnote>

<footnote n="72">⁷²A. Bashiri, A. Vaskin, K. Tanaka, M. Steinert, T. Pertsch, and I. Staude, *ACS Nano* **18**, 506 (2023).</footnote>

<footnote n="73">⁷³K. Koshelev, S. Lepeshov, M. Liu, A. Bogdanov, and Y. Kivshar, *Phys. Rev. Lett.* **121**, 193903 (2018).</footnote>

<footnote n="74">⁷⁴A. H. Schokker and A. F. Koenderink, *Optica* **3**, 686 (2016).</footnote>

<footnote n="75">⁷⁵P. Törmä and W. L. Barnes, *Rep. Prog. Phys.* **78**, 013901 (2014).</footnote>

<footnote n="76">⁷⁶D. Wang, W. Wang, M. P. Knudson, G. C. Schatz, and T. W. Odom, *Chem. Rev.* **118**, 2865 (2017).</footnote>

<footnote n="77">⁷⁷T. Gu, H. J. Kim, C. Rivero-Baleine, and J. Hu, *Nat. Photonics* **17**, 48 (2023).</footnote>

<footnote n="78">⁷⁸A. F. Koenderink, *ACS Photonics* **4**, 710 (2017).</footnote>

<footnote n="79">⁷⁹A. Krasnok, S. Glybovski, M. Petrov, S. Makarov, R. Savelev, P. Belov, C. Simovski, and Y. S. Kivshar, *Appl. Phys. Lett.* **108**, 211105 (2016).</footnote>
80K. Tanaka, E. Plum, J. Y. Ou, T. Uchino, and N. I. Zheludev, Phys. Rev. Lett. 105, 227403 (2010).

81L. Zhu, S. Yuan, C. Zeng, and J. Xia, Adv. Opt. Mater. 8, 1901830 (2020).

82M. Ramezani, A. Halpin, A. I. Fernández-Dominguez, J. Feist, S. R.-K. Rodriguez, F. J. García-Vidal, and J. Gómez Rivas, Optica 4, 31 (2017).

83A. Vaskin, J. Bohn, K. E. Chong, T. Bucher, M. Zilk, D.-Y. Choi, D. N. Neshev, Y. S. Kivshar, T. Pertsch, and I. Staude, ACS Photonics 5, 1359 (2018).

84J. Lin, Q. Sun, W. Feng, S. Guo, Z. Liu, H. Liang, and J. Li, Adv. Photonics Res. 2, 2000145 (2021).

85Y.-Y. Xie, P.-N. Ni, Q.-H. Wang, Q. Kan, G. Briere, P.-P. Chen, Z.-Z. Zhao, A. Delga, H.-R. Ren, H.-D. Chen, C. Xu, and P. Genevet, Nat. Nanotechnol. 15, 125 (2020).

86C. Huang, C. Zhang, S. Xiao, Y. Wang, Y. Fan, Y. Liu, N. Zhang, G. Qu, H. Ji, J. Han, L. Ge, Y. S. Kivshar, and Q. Song, Science 367, 1018 (2020).

87X. Zhang, Y. Liu, J. Han, Y. S. Kivshar, and Q. Song, Science 377, 1215 (2022).

88P. P. Iyer, N. Karl, S. Addamane, S. D. Gennaro, M. B. Sinclair, and I. Brener, Nat. Photonics 17, 588 (2023).

89J. Bohn, T. Bucher, K. E. Chong, A. Komar, D.-Y. Choi, D. N. Neshev, Y. S. Kivshar, T. Pertsch, and I. Staude, Nano Lett. 18, 3461 (2018).

90Y. Chang, J. Wei, and C. Lee, Nanophotonics 9, 3049 (2020).

91Y. Chen, J. Feng, Y. Huang, W. Chen, R. Su, S. Ghosh, Y. Hou, Q. Xiong, and C.-W. Qiu, Nat. Mater. 22, 1065 (2023).

92T. Bucher, A. Vaskin, R. Mupparapu, F. J. Löchner, A. George, K. E. Chong, S. Fasold, C. Neumann, D.-Y. Choi, F. Eilenberger, F. Setzpfandt, Y. S. Kivshar, T. Pertsch, A. Turchanin, and I. Staude, ACS Photonics 6, 1002 (2019).

93S. Vadia, J. Scherzer, K. Watanabe, T. Taniguchi, and A. Högele, Nano Lett. 23, 614 (2023).

94B. Gallinet, J. Butet, and O. J. Martin, Laser Photonics Rev. 9, 577 (2015).

95G. Yoon and J. Rho, Comput. Phys. Commun. 264, 107846 (2021).

96Y. Chen, Y. Zhang, and A. Femius Koenderink, Opt. Express 25, 21358 (2017).

97F. Capolino, D. R. Jackson, D. R. Wilton, and L. B. Felsen, IEEE Trans. Antennas Propagat. 55, 1644 (2007).

98S. Zhang, E. R. Martins, A. G. Diyaf, J. I. Wilson, G. A. Turnbull, and I. D. Samuel, Synth. Met. 205, 127 (2015).

99V. Su, C. Chu, G. Sun, and D. Tsai, Opt. Express 26, 13148 (2018).

100A. Vaskin, S. Liu, S. Addamane, P. P. Vabishchevich, Y. Yang, G. Balarishnan, M. B. Sinclair, T. Pertsch, I. Brener, and I. Staude, Opt. Express 29, 5567 (2021).

101I. Staude, V. V. Khardikov, N. T. Tofang, S. Liu, M. Decker, D. N. Neshev, T. S. Luk, I. Brener, and Y. S. Kivshar, ACS Photonics 2, 172 (2015).

102Z. F. Sadrieva, I. S. Sinev, K. L. Koshelev, A. Samusev, I. V. Iorsh, O. Takayama, R. Mureanu, A. A. Bogdanov, and A. V. Lavrinenko, ACS Photonics 4, 723 (2017).

103N. Li, Z. Xu, Y. Dong, T. Hu, Q. Zhong, Y. H. Fu, S. Zhu, and N. Singh, Nanophotonics 9, 3071 (2020).

104G. Brière, P. Ni, S. Héron, S. Chenot, S. Vézian, V. Brändli, B. Damilano, J.-Y. Duboz, M. Iwanaga, and P. Genevet, Adv. Opt. Mater. 7, 1801271 (2019).

105A. I. Kuznetsov, A. B. Evlyukhin, M. R. Goncalves, C. Reinhardt, A. Koroleva, M. L. Arnedillo, R. Kiyan, O. Marti, and B. N. Chichkov, ACS Nano 5, 4843 (2011).

106M. Vellaichamy, M. Škarabot, and I. Mušević, Liq. Cryst. 50, 935 (2023).

107S. Karaveli, S. Wang, G. Xiao, and R. Zia, ACS Nano 7, 7165 (2013).

108R. A. DeCrescent, N. R. Venkatesan, C. J. Dahlman, R. M. Kennard, X. Zhang, W. Li, X. Du, M. L. Chabinyc, R. Zia, and J. A. Schuller, Sci. Adv. 6, eaay4900 (2020).

109R. Kolkowski and A. Shevchenko, Nanophotonics 12, 3443 (2023).

110M. Yoshida, S. Katsuno, T. Inoue, J. Gelleta, K. Izumi, M. De Zoysa, K. Ishizaki, and S. Noda, Nature 618, 727 (2023).

111J. W. You, Z. Lan, Q. Ma, Z. Gao, Y. Yang, F. Gao, M. Xiao, and T. J. Cui, Photonics Res. 11, B65 (2023).

112R. El-Ganainy, K. G. Makris, M. Khajavikhan, Z. H. Musslimani, S. Rotter, and D. N. Christodoulides, Nat. Phys. 14, 11–19 (2018).

113Y. Zeng, U. Chattopadhyay, B. Zhu, B. Qiang, J. Li, Y. Jin, L. Li, A. G. Davies, E. H. Linfield, B. Zhang, Y. Chong, and Q. J. Wang, Nature 578, 246 (2020).

114M. Lyubarov, Y. Lumer, A. Dikopoltsev, E. Lustig, Y. Sharabi, and M. Segev, Science 377, 425 (2022).

115R. W. Boyd, A. L. Gaeta, and E. Giese, Springer Handbook of Atomic, Molecular, and Optical Physics (Springer, 2008), pp. 1097–1110.

116T. Santiago-Cruz, V. Sultanov, H. Zhang, L. Krivitsky, and M. Chekhova, Opt. Lett. 46, 653 (2021).

117G. B. Lemos, V. Borish, G. D. Cole, S. Ramelow, R. Lapkiewicz, and A. Zeilinger, Nature 512, 409 (2014).

118G. Brida, M. Genovese, and I. R. Berchera, Nat. Photonics 4, 227 (2010).

119G. Triginer Garces, H. M. Chrzanowski, S. Daryanoosh, V. Thiel, A. L. Marchant, R. B. Patel, P. C. Humphreys, A. Datta, and I. A. Walmsley, Appl. Phys. Lett. 117, 024002 (2020).

120C. A. Casacio, L. S. Madsen, A. Terrasson, M. Waleed, K. Barnscheidt, B. Hage, M. A. Taylor, and W. P. Bowen, Nature 594, 201 (2021).

121E. Pearce, N. R. Gemmell, J. Flórez, J. Ding, R. F. Oulton, A. S. Clark, and C. C. Phillips, Opt. Continuum 2, 2386 (2023).

122S. R. Marder, D. N. Beratan, and L.-T. Cheng, Science 252, 103 (1991).

123S. R. Marder, L.-T. Cheng, B. G. Tiemann, A. C. Friedli, M. Blanchard-Desce, J. W. Perry, and J. Skindhøj, Science 263, 511 (1994).

124M. P. Nielsen, X. Shi, P. Dichtl, S. A. Maier, and R. F. Oulton, Science 358, 1179 (2017).

125A. Autere, H. Jussila, Y. Dai, Y. Wang, H. Lipsanen, and Z. Sun, Adv. Mater. 30, 1705963 (2018).

126H.-W. Guo, Z. Hu, Z.-B. Liu, and J.-G. Tian, Adv. Funct. Mater. 31, 2007810 (2021).

127Z. Sun, A. Martinez, and F. Wang, Nat. Photonics 10, 227 (2016).

128B. Guo, Q.-l. Xiao, S.-h. Wang, and H. Zhang, Laser Photonics Rev. 13, 1800327 (2019).

129K. O'Brien, H. Suchowski, J. Rho, A. Salandrino, B. Kante, X. Yin, and X. Zhang, Nat. Mater. 14, 379 (2015).

130T. Pertsch and G. Kivshar, MRS Bull. 45, 210 (2020).

131V. Giannini, A. I. Fernández-Domínguez, S. C. Heck, and S. A. Maier, Chem. Rev. 111, 3888 (2011).

132D. Smirnova and Y. S. Kivshar, Optica 3, 1241 (2016).

133P. Genevet, F. Capasso, F. Aieta, M. Khorasaninejad, and R. Devlin, Optica 4, 139 (2017).

134M. F. Limonov, M. V. Rybin, A. N. Poddubny, and Y. S. Kivshar, Nat. Photonics 11, 543 (2017).

135B. Sain, C. Meier, and T. Zentgraf, Adv. Photonics 1, 024002 (2019).

136G. Quaranta, G. Basset, O. J. Martin, and B. Gallinet, Laser Photonics Rev. 12, 1800017 (2018).

137J. Zhang, J. Ma, M. Parry, M. Cai, R. Camacho-Morales, L. Xu, D. Neshev, and A. Sukhorukov, Sci. Adv. 8, 4240 (2022).

138L. Michaeli, S. Keren-Zur, O. Avayu, H. Suchowski, and T. Ellenbogen, Phys. Rev. Lett. 118, 243904 (2017).

139R. Czaplicki, A. Kiviniemi, J. Laukkanen, J. Lehtolahti, M. Kuittinen, and M. Kauranen, Opt. Lett. 41, 2684 (2016).

140K. Koshelev, A. Bogdanov, and Y. Kivshar, Sci. Bull. 64, 836 (2019).

141C. W. Hsu, B. Zhen, A. D. Stone, J. D. Joannopoulos, and M. Soljačić, Nat. Rev. Mater. 1, 16048 (2016).

142Y. Zhao, Y. Yang, and H.-B. Sun, PhotonX 2, 3 (2021).

143B. Metzger, L. Gui, J. Fuchs, D. Floess, M. Hentschel, and H. Giessen, Nano Lett. 15, 3917 (2015).

144G. Li, S. Chen, N. Pholchai, B. Reineke, P. W. H. Wong, E. Y. B. Pun, K. W. Cheah, T. Zentgraf, and S. Zhang, Nat. Mater. 14, 607 (2015).

145M. R. Shcherbakov, D. N. Neshev, B. Hopkins, A. S. Shorokhov, I. Staude, E. V. Melik-Gaykazyan, M. Decker, A. A. Ezhov, A. E. Miroshnichenko, I. Brener, A. A. Fedyanin, and Y. S. Kivshar, Nano Lett. 14, 6488 (2014).

146G. Grinblat, Y. Li, M. P. Nielsen, R. F. Oulton, and S. A. Maier, Nano Lett. 16, 4635 (2016).

147M. R. Shcherbakov, A. S. Shorokhov, D. N. Neshev, B. Hopkins, I. Staude, E. V. Melik-Gaykazyan, A. A. Ezhov, A. E. Miroshnichenko, I. Brener, A. A. Fedyanin, and Y. S. Kivshar, ACS Photonics 2, 578 (2015).

148T. Shibanuma, G. Grinblat, P. Albella, and S. A. Maier, Nano Lett. 17, 2647 (2017).

149H. Liu, C. Guo, G. Vampa, J. L. Zhang, T. Sarmiento, M. Xiao, P. H. Bucksbaum, J. Vučković, S. Fan, and D. A. Reis, Nat. Phys. 14, 1006 (2018).
¹⁵⁰B. Hopkins, D. S. Filonov, A. E. Miroshnichenko, F. Monticone, A. Alù, and Y. S. Kivshar, ACS Photonics 2, 724 (2015).

¹⁵¹A. P. Anthur, H. Zhang, R. Paniagua-Dominguez, D. A. Kalashnikov, S. T. Ha, T. W. Maß, A. I. Kuznetsov, and L. Krivitsky, Nano Lett. 20, 8745 (2020).

¹⁵²K. Koshelev, S. Kruk, E. Melik-Gaykazyan, J.-H. Choi, A. Bogdanov, H.-G. Park, and Y. Kivshar, Science 367, 288 (2020).

¹⁵³Y. Yang, J. Lu, A. Manjavacas, T. S. Luk, H. Liu, K. Kelley, J.-P. Maria, E. L. Runnerstrom, M. B. Sinclair, S. Ghimire et al., Nat. Phys. 15, 1022 (2019).

¹⁵⁴M. Z. Alam, S. A. Schulz, J. Upham, I. D. Leon, and R. W. Boyd, Nat. Photonics 12, 79 (2018).

¹⁵⁵A. Xomalis, X. Zheng, R. Chikkaraddy, Z. Koczor-Benda, E. Miele, E. Rosta, G. A. Vandenbosch, A. Martínez, and J. J. Baumberg, Science 374, 1268 (2021).

¹⁵⁶L. Kühner, L. Sortino, R. Berté, J. Wang, H. Ren, S. A. Maier, Y. Kivshar, and A. Tittl, Nat. Commun. 13, 4992 (2022).

¹⁵⁷C. Gigli and G. Leo, Opto-Electron. Adv. 5, 210093 (2022).

¹⁵⁸S. Gennaro, M. Rahmani, V. Giannini, H. Aouani, T. Sidiropoulos, M. Navarro-Cía, S. Maier, and R. Oulton, Nano Lett. 16, 5278 (2016).

¹⁵⁹E. Rahimi and R. Gordon, Adv. Opt. Mater. 6, 1800274 (2018).

¹⁶⁰P.-Y. Chen, C. Argyropoulos, and A. Alù, Nanophotonics 1, 221 (2012).

¹⁶¹L. Wang, S. Kruk, K. Koshelev, I. Kravchenko, B. Luther-Davies, and Y. Kivshar, Nano Lett. 18, 3978 (2018).

¹⁶²M. Celebrano, X. Wu, M. Baselli, S. Großmann, P. Biagioni, A. Locatelli, C. De Angelis, G. Cerullo, R. Osellame, B. Hecht et al., Nat. Nanotechnol. 10, 412 (2015).

¹⁶³S. Kruk and Y. Kivshar, ACS Photonics 4, 2638 (2017).

¹⁶⁴M. W. Klein, C. Enkrich, M. Wegener, and S. Linden, Science 313, 502 (2006).

¹⁶⁵Y. Yang and S. I. Bozhevolnyi, Nanotechnology 30, 204001 (2019).

¹⁶⁶G. Grinblat, Y. Li, M. P. Nielsen, R. F. Oulton, and S. A. Maier, ACS Photonics 4, 2144 (2017).

¹⁶⁷Y. Yang, V. A. Zenin, and S. I. Bozhevolnyi, ACS Photonics 5, 1960 (2018).

¹⁶⁸G. F. Walsh and L. Dal Negro, Nano Lett. 13, 3111 (2013).

¹⁶⁹K. Thyagarajan, J. Butet, and O. J. Martin, Nano Lett. 13, 1847 (2013).

¹⁷⁰L. Wang, S. Kruk, L. Xu, M. Rahmani, D. Smirnova, A. Solntsev, I. Kravchenko, D. Neshev, and Y. Kivshar, Nanoscale 9, 2201 (2017).

¹⁷¹T. Liu, S. Xiao, B. Li, M. Gu, H. Luan, and X. Fang, Front. Nanotechnol. 4, 891892 (2022).

¹⁷²S. Kruk, M. Weismann, A. Y. Bykov, E. A. Mamonov, I. A. Kolmychek, T. Murzina, N. C. Panoiu, D. N. Neshev, and Y. S. Kivshar, ACS Photonics 2, 1007 (2015).

¹⁷³L. Wang, A. S. Shorokhov, P. N. Melentiev, S. Kruk, M. Decker, C. Helgert, F. Setzpfandt, A. A. Fedyanin, Y. S. Kivshar, and D. N. Neshev, ACS Photonics 3, 1494 (2016).

¹⁷⁴S.-D. Liu, E. S. P. Leong, G.-C. Li, Y. Hou, J. Deng, J. H. Teng, H. C. Ong, and D. Y. Lei, ACS Nano 10, 1442 (2016).

¹⁷⁵F. H. Stillinger and D. R. Herrick, Phys. Rev. A 11, 446 (1975).

¹⁷⁶D. Marinica, A. Borisov, and S. Shabanov, Phys. Rev. Lett. 100, 183902 (2008).

¹⁷⁷Z. Liu, J. Wang, B. Chen, Y. Wei, W. Liu, and J. Liu, Nano Lett. 21, 7405 (2021).

¹⁷⁸N. Bernhardt, K. Koshelev, S. J. White, K. W. C. Meng, J. E. Froch, S. Kim, T. Tran, D.-Y. Choi, Y. Kivshar, and A. S. Solntsev, Nano Lett. 20, 5309 (2020).

¹⁷⁹M. J. Huttunen, P. Rasekh, R. W. Boyd, and K. Dolgaleva, Phys. Rev. A 97, 053817 (2018).

¹⁸⁰S. Keren-Zur, M. Tal, S. Fleischer, D. M. Mittleman, and T. Ellenbogen, Nat. Commun. 10, 1778 (2019).

¹⁸¹E. Minerbi, S. Keren-Zur, and T. Ellenbogen, Nano Lett. 19, 6072 (2019).

¹⁸²D. Polyushkin, E. Hendry, E. Stone, and W. Barnes, Nano Lett. 11, 4718 (2011).

¹⁸³L. Luo, I. Chatzakis, J. Wang, F. B. Niesler, M. Wegener, T. Koschny, and C. M. Soukoulis, Nat. Commun. 5, 3055 (2014).

¹⁸⁴M. Z. Alam, I. De Leon, and R. W. Boyd, Science 352, 795 (2016).

¹⁸⁵L. Caspani, R. Kaipurath, M. Clerici, M. Ferrera, T. Roger, J. Kim, N. Kinsey, M. Pietrzyk, A. Di Falco, V. Shalaev, A. Boltasseva, and D. Faccio, Phys. Rev. Lett. 116, 233901 (2016).

¹⁸⁶N. Engheta, Science 340, 286 (2013).

¹⁸⁷T. S. Luk, D. De Ceglia, S. Liu, G. A. Keeler, R. P. Prasankumar, M. A. Vincenti, M. Scalora, M. B. Sinclair, and S. Campione, Appl. Phys. Lett. 106, 151103 (2015).

¹⁸⁸X. Niu, X. Hu, S. Chu, and Q. Gong, Adv. Opt. Mater. 6, 1701292 (2018).

¹⁸⁹S. Campione, I. Brener, and F. Marquier, Phys. Rev. B 91, 121408 (2015).

¹⁹⁰L. Li, Z. Liu, X. Ren, S. Wang, V.-C. Su, M.-K. Chen, C. Chu, H. Kuo, B. Liu, W. Zang, G. Guo, L. Zhang, Z. Wang, S. Zhu, and D. Tsai, Science 368, 1487 (2020).

¹⁹¹H. Aouani, M. Rahmani, M. Navarro-Cía, and S. A. Maier, Nat. Nanotechnol. 9, 290 (2014).

¹⁹²A. Melikyan, L. Alloatti, A. Muslija, D. Hillerkuss, P. C. Schindler, J. Li, R. Palmer, D. Korn, S. Muehlbrandt, D. Van Thourhout et al., Nat. Photonics 8, 229 (2014).

¹⁹³T. J. Duffin, M. P. Nielsen, F. Diaz, S. Palomba, S. A. Maier, and R. F. Oulton, Opt. Lett. 41, 155 (2016).

¹⁹⁴W. Chen, P. Roelli, H. Hu, S. Verlekar, S. P. Amirtharaj, A. I. Barreda, T. J. Kippenberg, M. Kovylina, E. Verhagen, A. Martínez et al., Science 374, 1264 (2021).

¹⁹⁵C. Gigli, T. Wu, G. Marino, A. Borne, G. Leo, and P. Lalanne, ACS Photonics 7, 1197 (2020).

¹⁹⁶K. Frizyuk, I. Volkovskaya, D. Smirnova, A. Poddubny, and M. Petrov, Phys. Rev. B 99, 075425 (2019).

¹⁹⁷R. Camacho-Morales, M. Rahmani, S. Kruk, L. Wang, L. Xu, D. A. Smirnova, A. S. Solntsev, A. Miroshnichenko, H. H. Tan, F. Karouta et al., Nano Lett. 16, 7191 (2016).

¹⁹⁸V. F. Gili, L. Carletti, F. Chouchane, G. Wang, C. Ricolleau, D. Rocco, A. Lemaître, I. Favero, L. Ghirardini, M. Finazzi, M. Celebrano, C. De Angelis, and G. Leo, Nanophotonics 7, 517 (2017).

¹⁹⁹A. V. Kildishev, A. Boltasseva, and V. M. Shalaev, Science 339, 1232009 (2013).

²⁰⁰S. Keren-Zur, O. Avayu, L. Michaeli, and T. Ellenbogen, ACS Photonics 3, 117 (2016).

²⁰¹S. Chen, G. Li, F. Zeuner, W. H. Wong, E. Y. B. Pun, T. Zentgraf, K. W. Cheah, and S. Zhang, Phys. Rev. Lett. 113, 033901 (2014).

²⁰²W. Ye, F. Zeuner, X. Li, B. Reineke, S. He, C.-W. Qiu, J. Liu, Y. Wang, S. Zhang, and T. Zentgraf, Nat. Commun. 7, 11930 (2016).

²⁰³B. Liu, B. Sain, B. Reineke, R. Zhao, C. Meier, L. Huang, Y. Jiang, and T. Zentgraf, Adv. Opt. Mater. 8, 1902050 (2020).

²⁰⁴M. Jin, Z. Hu, X. Liu, and G. Li, Laser Photonics Rev. 17, 2300174 (2023).

²⁰⁵N. Segal, S. Keren-Zur, N. Hendler, and T. Ellenbogen, Nat. Photonics 9, 180 (2015).

²⁰⁶O. Wolf, S. Campione, A. Benz, A. P. Ravikumar, S. Liu, T. S. Luk, E. A. Kadlec, E. A. Shaner, J. F. Klem, M. B. Sinclair et al., Nat. Commun. 6, 7667 (2015).

²⁰⁷E. Almeida, O. Bitton, and Y. Prior, Nat. Commun. 7, 12533 (2016).

²⁰⁸H. Ren, X. Fang, J. Jang, J. Bürger, J. Rho, and S. A. Maier, Nat. Nanotechnol. 15, 948 (2020).

²⁰⁹T. Cui, B. Bai, H.-B. Sun, T. Cui, B. Bai, and H. Sun, Adv. Funct. Mater. 29, 1806692 (2019).

²¹⁰W. Cai, A. P. Vasudev, and M. L. Brongersma, Science 333, 1720 (2011).

²¹¹W. Ding, L. Zhou, and S. Y. Chou, Nano Lett. 14, 2822 (2014).

²¹²L. Kang, Y. Cui, S. Lan, S. P. Rodrigues, M. L. Brongersma, and W. Cai, Nat. Commun. 5, 4680 (2014).

²¹³K.-T. Lee, M. Taghinejad, J. Yan, A. S. Kim, L. Raju, D. K. Brown, and W. Cai, ACS Photonics 6, 2663 (2019).

²¹⁴J. Yu, S. Park, I. Hwang, D. Kim, F. Demmerle, G. Boehm, M.-C. Amann, M. A. Belkin, and J. Lee, Nat. Photonics 16, 72 (2022).

²¹⁵S. Kim, J. Jin, Y.-J. Kim, I.-Y. Park, Y. Kim, and S.-W. Kim, Nature 453, 757 (2008).

²¹⁶I.-Y. Park, S. Kim, J. Choi, D.-H. Lee, Y.-J. Kim, M. F. Kling, M. I. Stockman, and S.-W. Kim, Nat. Photonics 5, 677 (2011).

²¹⁷M. Krüger, M. Schenk, and P. Hommelhoff, Nature 475, 78 (2011).

²¹⁸G. Herink, D. R. Solli, M. Gulde, and C. Ropers, Nature 483, 190 (2012).

²¹⁹S. Sederberg, F. Kong, F. Hufnagel, C. Zhang, E. Karimi, and P. B. Corkum, Nat. Photonics 14, 680 (2020).

²²⁰S. Han, H. Kim, Y. W. Kim, Y.-J. Kim, S. Kim, I.-Y. Park, and S.-W. Kim, Nat. Commun. 7, 13105 (2016).