non-abrupt and smooth variation of the electron density near material surfaces (e.g., Friedel oscillations and electron spill-out), and (3) absorption processes involving non-negligible change of quasiparticle momentum (e.g., surface-enabled Landau damping). While such single-effect explorations have yielded valuable physical and conceptual insights, they are challenged by the reality that no single beyond-CEM effect strongly dominates the others in general and, moreover, that they may act in opposition. In extreme cases, including only one effect may decrease predictive power, as in the case of hydrodynamic accounts of nonlocality for simple metal surface plasmons. Observations of beyond-CEM deviations should therefore generally be expected as arising from the interplay of multiple, physically distinct effects.

Feibelman identified this challenge in the early eighties and introduced his eponymous $d_⊥$ and $d_‖$ parameters as a possible remedy [192,309], further developed with subsequent contributions from Apell, Liebsch, and others [310–312]. In particular, the CEM implies that induced charge density is accumulated only at interfaces, and consequently amounts to a zeroth-order (monopole) approximation of a multipole expansion of a continuously varying induced charge density. $d_⊥$ parameterizes the omitted first-order (dipole) term of this expansion and $d_‖$ parameterizes the discrepancy between the classical and quantum-mechanical monopole terms (i.e., total charge); together, they capture the leading-order corrections to CEM [230,294]. Considering a planar interface (residing at, say, $x = 0$), $d_⊥$ and $d_‖$ are the frequency-dependent centroids of the induced charge density $\rho(x, \omega)e^{iky}$ [Fig. 16(a)] and the normal derivative of the tangential current density $J_y(x, \omega)e^{iky}$ due to an incident plane wave propagating in the xy-plane at frequency $\omega$ [192,310,312]:

$$d_{\perp}(\omega) = \frac{\int x \rho(x, \omega) dx}{\int \rho(x, \omega) dx}, \qquad d_{\parallel}(\omega) = \frac{\int x \partial_x J_y(x, \omega) dx}{\int \partial_x J_y(x, \omega) dx}. \qquad (53)$$

As centroids of induced quantities, $d_⊥$ and $d_‖$ have units of length with complex amplitudes in the Ångström-range for typical plasmonic materials. They represent surface-response functions of the interface (i.e., intrinsic functions of the material composition across the interface) that describe the tendency of interfaces to polarize, analogously to how the bulk permittivity $\epsilon$ describes the tendency of a bulk material to polarize. The associated surface polarization is [258]:

$$\pi = \epsilon_0 d_{\perp} \llbracket E_{\perp} \rrbracket \hat{n} - \epsilon_0 d_{\parallel} \llbracket D_{\parallel} / \epsilon_0 \rrbracket, \qquad (54)$$

with $\llbracket f \rrbracket = f(r + 0^+ \hat{n}) - f(r + 0^- \hat{n})$ denoting the discontinuity of a field $f$ across interface with normal $\hat{n}$ and with $f_⊥ = \hat{n} \cdot f$ and $f_‖ = (I - \hat{n}\hat{n}^T)f$ denoting the field’s normal and tangential parts. Unlike the bulk polarization, sourced by the E-field amplitude, the surface polarization is sourced by the interface-discontinuities $\llbracket E_⊥ \rrbracket$ and $\llbracket D_⊥ \rrbracket$. Compellingly, the corrections to CEM are thus driven by CEM’s failure to achieve field continuity.

While Feibelman’s original formulation focused only on planar interfaces, the surface polarization formulation—and the recognition that any surface-curvature dependence of $d_⊥$ is negligible for leading-order beyond-CEM considerations [294,311]—generalizes the scope to arbitrary photonic structures [230,258]. Moreover, as discovered three and a half decades after the Feibelman parameters’ introduction, the surface polarization can in fact be abstracted away through a modification of the CEM boundary conditions [258]:

$$\llbracket D_{\perp} \rrbracket = d_{\parallel} \nabla_{\parallel} \cdot \llbracket D_{\parallel} \rrbracket, \qquad \llbracket B_{\perp} \rrbracket = 0, \qquad (55a)$$

$$\llbracket E_{\parallel} \rrbracket = -d_{\perp} \nabla_{\parallel} \llbracket E_{\perp} \rrbracket, \qquad \llbracket H_{\parallel} \rrbracket = i \omega d_{\parallel} \llbracket D_{\parallel} \rrbracket \times \hat{n}. \qquad (55b)$$

While Eqs. (55) are amenable to direct incorporation in full-wave computational frameworks, the general smallness of beyond-CEM corrections naturally suggests perturbative approaches.
Such approaches have been pursued in both nonretarded [294] and retarded [258] settings, suitable for deep-subwavelength resonators and multi-scale structures respectively. In the latter setting, a quasi-normal mode perturbation theory [314] gives the leading-order shift $\omega_1$ of the complex resonance frequency $\omega = \omega_0 + \omega_1 + ...$ from the CEM frequency $\omega_0$ [258]:

$$ \omega_1 = \omega_0(\kappa_{\perp}d_{\perp} + \kappa_{\parallel}d_{\parallel}), \qquad (56) $$

with mode-, shape- and scale-dependent perturbation amplitudes $\kappa_{\perp,\parallel}$ (inverse length units), obtainable from CEM quasi-normal modes, and with $d_{\perp,\parallel}$ evaluated at $\omega_0$. The perturbative expression cleanly separates the mode-, shape-, and scale-dependence ($\kappa_{\perp,\parallel}$) from the material-dependence ($d_{\perp,\parallel}$) of beyond-CEM shifts. In the nonretarded limit and for a pure plasma (plasma

<figure><img src="image_18.png" /><figcaption>**Fig. 16.** Measuring Feibelman parameters. (a) The $d_{\perp}$ and $d_{\parallel}$ parameters are the frequency-dependent centroids of the induced charge $\rho(x)$ (shown) and the normal derivative of the interface-tangential current $\partial_x J_y(x)$ (not shown). Both $\rho(x)$ and $J_y(x)$ are induced quantities due to an incident plane wave, appearing as perturbations to the static equilibrium density $n(x)$. (b-d) Comparing observations of large deviations between CEM and experiments in gapped film-coupled nanodisks (b,c) with perturbation theory, Eq. (56), can enable experimental inference of the Feibelman parameters (d), here of $d_{\perp}$ for a gold (Au) aluminum oxide (AlO$_x$) interface [258]. (e,f) The angle- ($\theta$) and frequency-dependent amplitude contrast $\Psi$ and phase difference $\Delta$ of the complex reflection ratio $r_p/r_s = e^{i\Delta\tan\Psi}$ between p- and s-polarized light are conventionally used to obtain a material's permittivity by fitting to the classical Fresnel reflection coefficients $r_{p,s}$. By using the Feibelman-corrected reflection coefficients instead, sufficiently precise measurements of $\Psi$ and $\Delta$ can enable inference of the interface's Feibelman parameters (f), here of $d_{\perp}$ at an Au(111)-air interface [313]. (g,h) An external gate voltage of 10 V, applied over an Au nanoresonator, leads to a frequency-dependent shift $\Delta S$ of the nanoresonator's scattering signal, presented relative to its ungated signal $S_0$ (b). Measurements over different nanoresonator lengths $l$ are labeled, shifted vertically, and differentiated by color. The shift of the scattering signal may be interpreted as due to a gate-tunable (interface-averaged) shift $\langle\Delta d_{\perp}\rangle$ of the $d_{\perp}$-parameter at the Au-air interface (h), obtained by fitting to experimental data [236]. Inset in (g) schematically illustrates the nanoresonator, its dimensions, and gating configuration. Panels (a,c,d) adapted with permission from Ref. [258] (Copyright © 2019 Springer Nature).; (e,f) adapted from Ref. [313]; (g,h) adapted from Ref. [236] (American Association for the Advancement of Science, CC BY 4.0).</figcaption></figure>
frequency $\omega_p$) adjacent to vacuum, analytical expressions can e.g., be obtained for the resonance frequency $\omega_k^{SP}$ of a surface plasmon wave vector $k$ at a planar interface ($\kappa_{\perp} = -\kappa_{\parallel} = -\frac{1}{2}k$) and for the lth order localized surface plasmon frequency $\omega_l^{LSP}$ of a nanosphere of radius $R$ ($\kappa_{\perp} = -\kappa_{\parallel} = -\frac{l+1}{2}R^{-1}$) [59,294,312]:

$$ \omega_k^{SP} = \frac{1}{\sqrt{2}} \omega_p (1 - \frac{1}{2} d_{eff} k), \quad \omega_l^{LSP} = \frac{1}{\sqrt{2+l^{-1}}} \omega_p (1 - \frac{l+1}{2} d_{eff} R^{-1}), \qquad (57) $$

with $d_{eff} = d_{\perp} - d_{\parallel}$ evaluated at $\omega_0$. The incorporation of Feibelman parameters lifts the characteristic scale-invariance of the nonretarded CEM by introducing the length-scale of $d_{\perp,\parallel}$ to compete with the geometric or modal length-scale $L$, resulting in corrections of order $O(d_{eff}L^{-1})$. These simple examples also reveal the physical interpretations of the complex parts of $d_{eff}$: (i) Re $d_{eff}$ contributes to a shift of the frequency center, (ii) the sign of Re $d_{eff}$ determines the shift's direction (blue-shifting for charge "spill-in", Re $d_{eff} < 0$, and red-shifting for "spill-out", Re $d_{eff} > 0$), and (iii) Im $d_{eff}$ contributes size- or momentum-dependent broadening.

## Challenges and opportunities

While the introduction Feibelman parameters now lies more than four decades in the past, the framework's application to modern nanoplasmonic-response questions remains in an early stage and several open, unsolved challenges remain. Below, we highlight opportunities and challenges in the experimental and theoretical determination of the Feibelman d-parameters.

**Experimental measurement and tabulation.** A key strength of the Feibelman parameters is that they provide a physics-agnostic framework for incorporating beyond-CEM effects. No specific physical effects are included or excluded: the physics in the underlying microscopic response quantities $\rho$ and $J$—equivalently, $\varepsilon(\mathbf{r}, \mathbf{r}')$—on which $d_{\perp,\parallel}$ depends is automatically incorporated. Specific models of $d_{\perp,\parallel}$ can therefore usually encompass (to leading order) existing “single-effect” approaches, including e.g., hydrodynamic accounts of nonlocality [192]. This freedom and versatility, however, brings a concomitant challenge: the values and frequency-dependence of $d_{\perp,\parallel}(\omega)$ are not prescribed by the framework itself. For real materials, they must be obtained by experimental measurements or other sophisticated theoretical modeling.

An analogous problem arises for the bulk permittivity $\varepsilon(\omega)$, which has been met extensive experimental tabulation efforts [315]. By comparing observations of large beyond-CEM deviations in $AlO_x$–Au thin-film-gapped Au-disks [Fig. 16(b,c)], Yang et al. [258] demonstrated that equivalent efforts are feasible for the Feibelman parameters, by inferring $d_{\perp}$ of an Au–$AlO_x$ interface by comparison with perturbation theory [Fig. 16(d)]. Near-field measurements of highly confined thin-film gap surface plasmons at an Au–$Al_2O_3$ interface [85] demonstrate the feasibility more broadly, as do measurements of beyond-CEM deviations in graphene–metal heterostructures [92,316].

While feasible, these approaches often require multiple parametrically varied geometries and are intrinsically limited to the frequency range of the resonant plasmonic feature. Measurements of the Feibelman parameters across broad frequency ranges is significantly more challenging, especially above the screened plasma frequency, contending with increased smallness of beyond-CEM shifts away from resonance and worsened signal-to-noise ratios. Recent efforts have made strides towards broadband characterizations [Fig. 16(e,f)] [313], by incorporating the Feibelman parameters into ellipsometric analysis via Feibelman-corrected Fresnel reflection coefficients and, e.g., demonstrating that while Au spills in at typical surface plasmon frequencies (Re $d_{\perp} < 0$), it spills out above the plasma frequency (Re $d_{\perp} < 0$).

**Theory and first-principles modeling.** State-of-the-art computational modeling of the Feibelman parameters arguably remains closely aligned with techniques first developed and applied in the 1980s [126,192]. These techniques leverage time-dependent density functional theory (TDDFT), accounting for the dynamics of the highest-lying conduction electrons and
treating the remaining electrons through an effective jellium, i.e., uniform potential (see Sec. 9 of this Roadmap). This captures the essential physics of simple metals well, since they exhibit only negligible screening from lower-lying electrons, and correctly predicts spill-out ($\text{Re } d_{\perp}>0$) and beyond-CEM red-shifts ($\omega_1<0$), consistently with measurements [317]. Noble metals, however, are strongly screened by bound electrons, and exhibit beyond-CEM blue-shifting ($\omega_1>0$) associated with spill-in ($\text{Re } d_{\perp}<0$), contrary to jellium predictions. Even semi-classically screened generalizations [126,294] do not substantially remedy this mismatch. These facts call for beyond-jellium calculations of the Feibelman parameters, explicitly incorporating both bound electrons and ionic potentials, e.g., by leveraging modern TDDFT software.

First steps have been made through the incorporation of non-uniform atomic potentials which has e.g., shed light on the influence of crystallographic facet terminations [318]. Many important questions remain unaddressed, however, including:

1. Influence of dielectric screening for metal–dielectric interfaces [319]. From semiclassical considerations, $d_{\perp}$ is expected to scale roughly inversely dielectric permittivity for noble-metal-dielectric heterostructures [258]. No similar expectation approximate rules arise for simple metals, however, since their free electrons spill into the dielectric cladding, and a general, predictive, and accurate modeling scheme has not been developed.

2. Impact of residual chemical elements or few-layer intercalates at interfaces, including the wetting layers that form under ambient conditions. Their potential influence is crucial for the interpretation and assignment of experimentally obtained Feibelman parameters.

3. While $d_{\|}$ vanishes in charge-neutral jellium calculations [310] and has been assumed small in the experimental treatments, the presence of bound-electron screening [126], surface roughness [320], and electronic interface states [318,321] can lift this requirement. An understanding of the relative amplitudes of $d_{\|}$ and $d_{\perp}$ in general settings remains absent.

4. Tunability of $d_{\perp,\|}$ by external means, e.g., applied voltages or strain. Recent experiments by Zurak et al. [236] used differential scattering response under voltage gating of a plasmonic nanoresonator [Fig. 16(g)] to estimate the associated tunability of the Feibelman parameters, finding sensitivities on the order of $10^{-1}$ pm/V for $d_{\perp}$ at Au–air interfaces [Fig. 16(h)]. A full theoretical understanding of these experimental results as well as whether larger tunability is attainable by other external controls are exciting open questions.

5. Jellium TDDFT’s failure to model $\text{Re } d_{\perp,\|}$ (frequency shift) also limits its ability to model $\text{Im } d_{\perp,\|}$ (frequency broadening) due to Kramers–Kronig relations. The impact of bound-electron screening on beyond-CEM broadening is thus still poorly understood. One open question is whether some plasmonic materials might be naturally “quantum low-loss” (i.e., small $\text{Im } d_{\perp,\|}$) or “quantum high-loss” (i.e., large $\text{Im } d_{\perp,\|}$), which may constrain outlooks for low-loss ultra-confined plasmons and hot-electron generation, respectively.

6. The electron–phonon coupling and phonon-assisted transitions play a significant role in the bulk losses of materials [322]. Surprisingly little, however, is known about the electron–phonon coupling’s impact on surface damping. E.g., the contribution to $\text{Im } d_{\perp}$ from phonon-assisted transitions is not currently understood, nor is its relative magnitude compared to intraband electron–hole pair transitions (surface-enabled Landau damping).

Answering these questions will likely necessitate new computational and algorithmic developments for TD-DFT approaches of thin-film and interfacial optical response.

In both theoretical and experimental pursuits of the Feibelman parameters’ values, two facts are helpful to validate obtained results and ensure physical consistency. First, since the Feibelman
parameters are response functions, they are causal and consequently obey Kramers–Kronig relations, which equip them with several useful sum rules [323]. Second, passivity—i.e., that linear passive materials do no net work on the electromagnetic fields—also constrains the Feibelman parameters, just as it constrains the local permittivity by $Im(\omega\varepsilon) \ge 0$. In particular, passivity requires that $Im(\omega d_{\perp}[[\varepsilon^{-1}]]) \ge 0$ and $Im(\omega d_{\parallel}[[\varepsilon]]) \le 0$ for all frequencies ($Im \omega > 0$) [313].

## Future developments to address challenges

While the noted opportunities largely fall within Feibelman's original scope, i.e., anchored in plasmonics, future developments may well be driven most efficiently by pursuing an expanded scope and generalizations of the underlying ideas of surface-response functions.

One particularly apt and simple example of such possible generalizations of scope is doped semiconductors, which are excellent candidates for realizing large beyond-CEM corrections due to long Fermi wavelengths and pronounced nonlocal response [324]. As are generalizations to synthetically nonlocal meta-materials. In particular, we anticipate that both phenomenological and semiclassical modeling of the Feibelman parameters may be sufficient to match experiments, since these materials afford a clearer separation between the involved beyond-CEM length-scales than metals. In turn, this may inform and revitalize phenomenological modeling approaches for Feibelman parameters of metals.

More broadly, several other material and quasiparticle platforms support ultra-confined optical modes, which may correspondingly host substantial beyond-CEM corrections. Polar materials such as hexagonal boron nitride, for instance, supports highly-confined surface phonon polaritons [325]. Even ordinary high-index dielectrics, such as silicon, can support ultra-confined, nanometer-scale optical modes in the gaps of bowtie nano-antennae [326]. Lacking mobile charge carriers, the physical mechanisms for beyond-CEM effects must involve phonons in either material platform, rather than electrons alone: accordingly, the momentum-dependence of the electron–phonon coupling appears likely to play a key role. Importantly, regardless of the underlying physical mechanism, the Feibelman-parameter framework will be applicable due to its physics-agnostic formulation. Recent observations of super-thermal evaporation at water surfaces [327,328], suggest that even broader material classes may exhibit beyond-CEM effects associated with nonzero Feibelman parameters, with physical mechanisms that are wholly unrelated to the Fermi liquid perspective of plasmons. Even the classical effects of surface roughness can be abstracted into the framework [320], a fact that is important to hold in view for experimental inference of Feibelman parameters.

The existence of nonlocality, even the mere existence of an interface, is also known to enable nonlinear effects otherwise forbidden in a local-response bulk. This includes second-harmonic generation in centrosymmetric materials. Generalizations of Feibelman parameters to the nonlinear domain, e.g., for capturing hydrodynamic accounts of nonlinearity [132], would consequently be of substantial interest for describing surface-enabled nonlinear phenomena. The fact that the second-harmonic nonlinear susceptibility of an interface contains terms proportional to the interface-discontinuity of the electric field [329] is highly encouraging in this regard.

Taken together, these outlooks suggest the possibility of more far-reaching generalizations to physical problems beyond electromagnetism. At their core, the Feibelman parameters amount to a leading-order remedy of the field discontinuity imposed by the simple but unphysical assumptions of piecewise constant, and local response. Such assumptions are widespread in physics, spanning e.g., heat transport, acoustics, and semiconductor physics. These facts suggest that many notions of surface- and interface-response functions might find generalization and unification in the ideas underlying the Feibelman parameters.
## Concluding remarks

The Feibelman parameters provide a natural and elegant framework for incorporating the salient beyond-classical effects in electromagnetic modeling. While the physics of Feibelman parameters for jellium plasma have been well-understood for decades, many questions remain open for the real materials of interest to contemporary polaritonics, e.g., metals, doped semiconductors, polar materials, and high-index dielectrics. In addition to these questions, there are broad opportunities for generalizations and extensions of the underlying ideas to physics beyond linear electromagnetism. We hope this overview provides an interesting and illustrative selection of such questions and opportunities, and that it may serve to inspire future contributions to the field.

**Funding.** Villum Fonden (42106); Ms. Belinda Hung, the Asian Young Scientist Fellowship; Croucher Foundation; New Cornerstone Science Foundation through the Xplorer Prize.

# 13. Computational electrodynamics with Feibelman surface-response formalism

ULRICH HOHENESTER ID, LORENZ HUBER & JAVIER AIZPURUA

## Current status

Metallic nanoparticles can focus light down to extreme sub-wavelength volumes through excitation of surface plasmons, these are coherent surface charge oscillations of electrons at the metal surface [330]. The question of how to account for quantum surface effects around sharp particle features, small gap regions, or for incoming fields with high spatial variations such as those produced by quantum emitters in the vicinity of a particle, has puzzled researchers for several decades. A particularly appealing solution is due to Feibelman [192] (see also Sec. 12 of this Roadmap), who devised an approach where all quantum modifications at the boundary are lumped into two dispersive so-called Feibelman parameters, which provide a measure of the distance over which the electric fields deviate from those for a sharp interface. Feibelman parameters were first brought to the field of plasmonics in [58] for the explanation of confinement-induced plasmon shifts, and further developed to include surface-induced nonlocal effects in the optical response of arbitrarily-shaped plasmonic nanoparticles [230]. More recently [59,258], it has been suggested that Feibelman parameters can be directly incorporated into classical Maxwell solvers by introducing modified “mesoscopic” boundary conditions. By comparing results from modified Mie theory including Feibelman parameters with those of *ab initio* simulations based on time-dependent density-functional theory (TDDFT) [253], it has been shown that for electromagnetic field variations on the nanometer scale one should rather consider Feibelman parameters that are nonlocal in the directions tangential to the boundary, in accordance to the original work of Feibelman who considered for a planar slab geometry a parameter dependence on the wavevector components parallel to the interface. Despite the strong recent interest in Feibelman parameters and mesoscopic boundary conditions, there are surprisingly few attempts to implement them into computational Maxwell solvers. In the seminal paper of [258] the authors proposed two methods for dealing with such boundary conditions in the COMSOL Multiphysics software based on the finite-element method (FEM): either by introducing an additional scalar auxiliary potential together with a surface current, or by computing the quasinormal modes of the system under study and accounting for the modifications due to the Feibelman parameters in lowest order of perturbation theory. Both approaches lead to simulations that are slower in comparison to those using normal boundary conditions. As an alternative approach, in Refs. [307,331] we have developed the methodology for directly introducing mesoscopic boundary conditions into a Maxwell solver based on the boundary element method (BEM), using either local or nonlocal Feibelman parameters, and have implemented them in our software toolbox nanobem
[332]. Agreement with analytic results based on a modified Mie theory has been demonstrated for single and coupled spheres, together with the expected computational convergence when increasing the number of boundary elements. Further developments are needed to bring our software to a form that can be used flexibly by other users, and to set up a database for Feibelman parameters for various material combinations.

## Challenges and opportunities

For computational Maxwell solvers incorporating local or nonlocal Feibelman parameters, the main challenges concern the methodology and the numerical implementation. To understand the challenges, we have to briefly dwell on some of the technicalities. In BEM, the Calderon identities, see Eq. (9.55) of Ref. [330], relate the tangential electromagnetic fields at the particle boundary to the incoming fields, for instance those produced by plane waves or quantum emitters placed in the vicinity of the particle. Together with the boundary conditions, one can use the identities to compute the unknown tangential fields and, in turn, the fields everywhere else. When considering Feibelman parameters, one can start from the same Calderon identities, and only the field matching at the boundary has to be modified. Here one has to be careful about the tangential derivative $\nabla_\| (d_\perp \circ [[E_\perp]])$ in the mesoscopic boundary conditions [258], where the electric field component normal to the boundary $E_\perp$ additionally has to be related in BEM to the tangential magnetic field $H_\|$ through Ampère's law. This leaves us with two tangential derivatives in the mesoscopic boundary condition. Fortunately, they can be handled in the same manner as the hypersingular operator in the Galerkin implementation of BEM, see Eq. (11.40) of Ref. [330], namely by swapping the derivatives through integration by parts to the Raviart-Thomas shape elements $f_\nu(r_\|)$ that are used for the representation of the tangential electromagnetic fields at the boundary. Furthermore, the derivative $\nabla \cdot \hat{n} \times f_\nu(r_\|)$ originating from this partial integration, with $\hat{n}$ being the outer surface normal, has to be related to $\nabla_\| \cdot f_\nu(r_\|)$ using a suitable translation procedure [307]. The important point for our discussion here is that through this translation not only neighbor but also distant shape elements become directly coupled when matching the tangential fields at the boundary, even for local Feibelman parameters, as can be seen in Fig. 17(b,b*) for the matching matrix before and after translation. Quite generally, dense matrices play a central role in BEM [330] and the inversion of an additional matrix constitutes no major bottleneck in the computational approach. The same conclusion holds for nonlocal Feibelman parameters, where from the outset distant shape elements are coupled, see Fig. 17(c), and the translation procedure leads to matrices that are equally dense than those for local Feibelman parameters. With the exception of the filling of these matrices, the consideration of nonlocal Feibelman parameters thus represents no computational drawback. The implementation of local and nonlocal Feibelman parameters to solve Maxwell's equations within the BEM requires the availability of a reliable set of these parameters associated with the boundaries separating two different media, which can be obtained from a quantum calculation of the optical response of a flat interface separating the same pair of media. To that end, when it comes to describe an interface between a dielectric and a free-electron gas of s-like electrons, such as in sodium (Na) or aluminum (Al), a TDDFT where the metallic medium is considered to be a jellium turns to be very adequate to describe the surface response [293]. Interfaces involving noble metals such as gold (Au) or silver (Ag), with a strong effect of d-like electrons on the optical response require more sophisticated approaches.

## Future developments to address challenges

In the future, the computation of Feibelman parameters and the development of Maxwell solvers incorporating mesoscopic boundary conditions should go hand in hand. Only in combination they provide a powerful machinery for considering quantum surface effects in classical Maxwell solvers, and can serve the community as a flexible simulation toolkit. The current BEM implementations of the mesoscopic boundary conditions [307,331] are somewhat experimental, and are neither
<figure><img src="image_19.png" /><figcaption>**Fig. 17.** Schematics of BEM simulations including Feibelman parameters. (a) We start by computing local or nonlocal Feibelman parameters within an *ab initio* approach, which can be submitted to the mesoscopic boundary conditions and implemented in a BEM solver. The latter provides us with *T*-matrices that could be used for the simulation of particle clusters or particle arrays. (b,c) The lower panels show the matrices used in our BEM approach for matching the electromagnetic fields at the boundary, with a representation of the tangential electromagnetic fields in terms of Raviart-Thomas shape elements (see inset) as global degrees of freedom (DOF). (b) For local Feibelman parameters the matrix is spare, however, due to the technical details described in the text and in Refs. [307,331], the *b*\* matrix actually used in the BEM approach is already relatively dense. (c,c\*) For nonlocal Feibelman parameters the matching matrices are equally dense. As the inversion of the matching matrices is computationally expensive, the performance of our BEM approach is not too different for local and nonlocal Feibelman parameters.</figcaption></figure>

optimized for flexibility nor speed. In contrast to state-of-the-art BEM solvers, such as nanobem [332], the additional computation of the matrices for field matching must properly account for the novel features of a spatial cutoff for the nonlocal Feibelman parameters $d(r_\|)$ and a sufficiently fine integration grid to resolve the spatial variations of $d(r_\|)$. A first step in this direction has been undertaken in Ref. [331], where we have suggested a separation of the Feibelman parameters into a local contribution and a remainder, accounting for the most salient short- and long-range features. Feibelman parameters are expected to be of importance for sharp particle features, gap structures, or quantum emitters nearby particles, thus calling for a boundary discretization that is fine in the vicinity of the critical points and coarser further away. This requires refined integration routines that can efficiently deal with boundary elements of unequal size. Since the consideration of the boundary matching matrices is computationally equally expensive for local and nonlocal Feibelman parameters, one must assure that the filling of these matrices does not become the slowest part of the BEM simulations. Nanobem is a solver particularly suited for small-scale problems. However, with the recent developments of the toolbox it is now possible to compute *T*-matrices, which can incorporate the Feibelman parameters, and to directly load them into *T*-matrix solvers such as SMUTHI or treams. With this it is possible to simulate clusters, arrays, or metasurfaces. Mesoscopic boundary conditions can be also implemented in other Maxwell solvers. The approach
developed in Refs. [307,331] can be probably adapted without too many modifications for finite element method (FEM) solvers using a Galerkin scheme. As for other solvers, such as those based on the finite-difference time domain (FDTD) or discrete-dipole approximation (DDA) method, it has to be clarified how the spatial derivatives in the tangential directions can be handled. Regarding the challenges in the calculation of Feibelman parameters, extensive efforts have been invested in developing jellium TDDFT simulations to obtain local and nonlocal parameters of interfaces involving a purely free electron gas, however, the implementation of such a nonlocal scheme of parameters for noble metal crystallographic interfaces is also a need in nowadays nano-plasmonics. Gold and silver interfaces, for instance, are commonly used in practical optical spectroscopy implementations, where strong band-structure effects are in place induced by the d-band of more localized-electrons, and a strong presence of electronic image states also modifies the spectral dependence of the parameters at lower energies. All of these effects enforce the implementation of atomistic *ab initio* calculations of the surface response in periodic supercells, capable to address these issues.

### Concluding remarks

Mesoscopic boundary conditions incorporating Feibelman parameters have been recently introduced as a versatile tool to account for quantum surface effects at metal-dielectric interfaces. However, computational Maxwell solvers incorporating Feibelman parameters are still scarce. One exception is the boundary element method approach, for which Feibelman parameters that are either local or nonlocal in the directions tangential to the interface have been implemented successfully. In the future, we plan to render our software suitable for an easy-to-use and flexible open-source simulation toolbox, and to supplement it with Feibelman parameters computed with state-of-the-art *ab initio* Schrödinger solvers for various material systems.

**Funding.** Austrian Science Fund (10.55776/P37150).

## Part III—Nonlocal effects beyond free-electron bulk metals

### 14. Nonlocal response in doped semiconductors

#### MARTIJN WUBS ID

##### Current status

Doped semiconductors are a vast and active research field, see Ref. [333] for a recent review. It is textbook knowledge that the Drude model [42] can describe free-electron systems, not only metals but also doped semiconductors. In doped semiconductors the response can be due to combinations of electrons, and light and heavy holes. Free electrons typically have much smaller reduced masses than holes, so they dominate the response if present. Electron reduced masses $m^*$ of small-bandgap semiconductors are smaller than in metals. Free-carrier densities in doped semiconductors are much smaller than in metals. The net effect is that semiconductor plasma frequencies $\omega_p = \sqrt{ne^2/(\varepsilon_0 m^*)}$ occur in the infrared to THz regimes, compared to the ultraviolet frequencies for metals.

For nanostructured metals, it is known that the Drude model ceases being accurate because of nonlocal response and because of “spill-out” across geometric boundaries [see Sec. 3 of this Roadmap]. Typical sizes (particle radii, thin-film thicknesses) should be well below ten nanometers to see nonlocal effects in metals. Signatures of nonlocal response in metallic nanoparticles are the spectral blueshift as particles get smaller, and longitudinal resonances above the plasma frequency arising in confined geometries, and hydrodynamic surface nonlinearities. It is therefore interesting whether similar nonlocal effects occur in semiconductors.

In Ref. [336], few-nm sized doped semiconductor nanocrystals are studied in a range of doping concentrations, showing features ranging from size-quantization to classical plasmonic behavior. In Ref. [324] it was argued that the hydrodynamic Drude model (HDM) should apply and nonlocal behavior should be observable in larger doped semiconductor particles, of tens of
nanometers in size (considerably larger than for metals), described by the longitudinal dielectric function $\epsilon(\omega, k) = \epsilon_{\infty} - \omega_p^2/[\omega^2 + i\gamma\omega - \beta^2k^2]$. The surface plasmon resonance of a nanosphere in air then has a nonlocal blueshift of $\Delta\omega = \beta/(\sqrt{2}R)$. In small-bandgap semiconductors [such as indium antimonide (InSb) with $E_g = 0.17$ eV] electrons are thermally excited to the conduction band, and at low temperatures the nonlocal parameter becomes $\beta = 3k_B T/m_e^*$ [324,334]. For n-doped semiconductors on the other hand, the $\beta$ is given by $\sqrt{3}/5(\hbar/m^*)(3\pi^2n)^{1/3}$, where n is the electron density [324,334]. Thus the same hydrodynamic Drude response function can describe microscopically quite different situations. In both cases, exciton effects can be neglected because of large screening. And in both cases, the relative nonlocal blueshift $\Delta\omega/\omega_p$ for semiconductors—doped gallium arsenide (GaAs) or thermally excited InSb—is much larger than for metals [324], e.g., silver (Ag), see Fig. 18(b). For the doped semiconductor, $\Delta\omega/\omega_p \propto n^{-1/6}(m^*)^{-1/2}/R$, which shows that, interestingly, both the lower carrier densities and smaller reduced masses in semiconductors contribute to larger relative nonlocal blueshifts.

<figure><img src="image_20.png" /><figcaption>**Fig. 18.** Examples of theory and experiments on nonlocal response in semiconductors. (a) Comparison of hydrodynamic and Boltzmann-Mermin theories with electron energy loss spectroscopy experiments on InSb. Reproduced with permission from Ref. [334] (Copyright © 2006 American Physical Society). (b) Prediction of large relative nonlocal blueshifts of nanospheres made of doped GaAs or thermally excited InSb, as compared to Ag. Reproduced with permission from Ref. [324] (Copyright © 2017 Europhysics Letters Association). (c) Observation of nonlocal blueshift and additional thickness-dependent resonances above the bulk plasma frequency in an In:CdO thin film (here: thickness 20 nm). Reproduced with permission from Ref. [214] (Copyright © 2018 Nature Springer). (d) Plasmonic dispersion of compact and tunable InSb-based THz surface plasmon cavities. Data points, obtained by thermal tuning, deviate from the local-response curve. Reproduced with permission from Ref. [95] (Copyright © 2023 Nature Springer). (e) Nonlocal effects and increased electron gas temperature tend to pull the transmitted SHG efficiency resonance(s) in opposite directions. Reproduced with permission from Ref. [335] (Copyright © 2020 AIP Publishing). (f) Surface charge depletion and corresponding free-electron third-harmonic efficiency $\eta$ for various static external fields (in V/µm). Reproduced with permission from Ref. [234] (Copyright © 2022 American Physical Society).</figcaption></figure>

In linear optics, convincing hydrodynamic nonlocal behavior in semiconductors has been observed in planar geometries: infrared optical nonlocalities due to conduction electrons have been observed in indium (In)-doped cadmium oxide (CdO) thin films [214], both a nonlocal
blueshift of the main plasmon resonance and two additional longitudinal resonances at higher energies, for film thicknesses of 20 nm, see Fig. 18(c). In a recent preprint, heavily n-doped indium arsenide antimonide (InAsSb) thin films are also reported to exhibit nonlocal resonances, even for film thicknesses up to 200 nm [337]. Also, compact plasmonic-based THz cavities, with thermal activation of carriers across the small bandgap of InSb, have been shown to operate until the emergence of nonlocality and Landau damping [95], see Fig. 18(d).

Besides this progress in nonlocality in linear optics, there are interesting developments in the nonlinear response of doped semiconductors as well. Reference [335] described both the second- and the third-harmonic response from an indium tin oxide (ITO) nanolayer, modeled as a combination of bulk and hydrodynamic nonlinear response, see Fig. 18(e). Reference [234] proposed to use surface modulation of the equilibrium charge density of heavily doped semiconductors via an external static potential $n_0$, see Fig. 18(f). It combines the insights that the hydrodynamic nonlinear response is a surface rather than a bulk response, and that the third-order nonlinear polarization is proportional to $1/n_0^2$, so that a reduction of $n_0$ at the surface in response to a static potential can make the nonlinear response orders of magnitude stronger. A recent preprint [233] attributes the strong third-harmonic response of n-doped indium gallium arsenide (InGaAs) thin layers to free-electron nonlocal response, with measurements agreeing much better with hydrodynamic than with local-response theory.

## Challenges and opportunities

**Quantum-confined versus hydrodynamic nonlocality.** While for planar geometries, accurate agreement of hydrodynamic theory and experiment have been found [214,337], for nanoparticles and for other geometries their linear optical response has been harder to identify as being hydrodynamic. There can be challenges to make spherical particles, and a more fundamental challenge is to distinguish hydrodynamic from size-quantization effects that also produce blueshifts as sizes are reduced [336]. And thus non-classical phenomena, as observed in Ref. [338] for example, may be due to either. It is difficult to draw a hard boundary, but generally in larger systems, single-particle size quantization effects become less likely and before we reach larger sizes where classical plasmonics applies, the hydrodynamic Drude model is expected to apply. While "spill-out" on the atomic scale is less important for the larger doped semiconductor structures than for metals that exhibit nonlocal response, charge depletion and inhomogeneous doping levels are novel possible features—and opportunities—for the semiconductors.

**Embedded eigenstates.** Recently, it was argued that taking nonlocal response into account does not destroy embedded eigenstates [339]. On the contrary, nonlocal nanospheres feature infinitely many such states, and the nonlocality relaxes the size restrictions of local theories. These do not couple to light, but could be observed instead with electron beams. As nonlocal effects show up in larger systems in semiconductors, as discussed above, doped semiconductors may be the preferred system to observe nonlocal embedded eigenstates, probed by electrons.

**Two-fluid hydrodynamic response.** The optical response of semiconductors may be produced by more than one type of charge carrier, for example by both heavy and light holes, each with their corresponding plasma frequency. In Ref. [340], a hybridization theory for two-component carrier plasmas was proposed, with an antibonding mode where the two components move in phase depends sensitively on the doping densities, and the low-energy bonding mode with opposite charge alignment. The antibonding mode agreed well with experiment. A hydrodynamic two-fluid theory was also proposed for doped semiconductors [140,341]. Observing the predicted low-energy hybridized hydrodynamic mode is still an open challenge. More generally, in semiconductors it can be a challenge that plasmon and phonon resonances have similar energies.

**Complex nonlocal parameter.** In the standard hydrodynamic Drude model, the $\beta$ parameter is real-valued with $\beta^2 = 3v_F^2/5$ and describes quantum pressure convection, while in the GNOR extension it is complex-valued [64], with $Im\beta$ interpreted as representing induced
charge diffusion and accounting for size-dependent broadening of arbitrary plasmonic structures, including semiconductors. In the theory for semiconductors of Ref. [214], an imaginary part of β is due to the viscosity of the electron gas. Likewise in the theory for semiconductors in Ref. [337], a second viscosity term appears exactly as the imaginary part of β² in the standard hydrodynamic model. A complex β is also obtained in Halevi's formula β²(ω) = (3ω/5 + iγ/3)v_F²/(ω + iγ) [44], where β becomes a running coupling constant interpolating between the ω ≪ γ and ω ≫ γ limits, see Refs. [44,45,334]. Each of these microscopic mechanisms can be the cause of size-dependent damping, i.e. the commonly observed broadening of spectra for smaller plasmonic structures. Yet size-dependent damping could also be caused by surface roughness, or by a combination of these causes. This makes identifying the actual causes of damping in a given experiment challenging. For example, in the intriguing quest for signatures of viscous electron plasmas, as reviewed in Ref. [342], many causes will need to be suppressed or ruled out before size-dependent plasmonic damping can be uniquely attributed to viscosity of the electron gas, in semiconductor plasmonics and elsewhere.

**Beyond hydrodynamic theory.** The hydrodynamic Drude model, or the simple extension provided by the GNOR model [64], are robust in the sense that their nonlocal dielectric function accurately describes many optical experiments where nonlocality shows up [234,337]. A more general dispersion relation for plasmons would be

$$ \omega^2 + i\omega \left( \gamma_0 + \sum_{n=1} a_{2n} k^{2n} \right) = \omega_0^2 + \sum_{n=1} b_{2n} k^{2n}, \qquad (58) $$

with spatial dispersion both in the frequency and in the damping constants and with constants a₂ₙ and b₂ₙ. This general form reduces to ω² = ωₚ²/ε∞ for the local Drude model and to ω² + iγω = ωₚ²/ε∞ + β²k² for the HDM and GNOR models, in each case found by putting the corresponding epsilon to zero. So in the hydrodynamic theory, both series in powers of k² truncate beyond k². One reason why higher-order terms would be masked in experiments is the damping, where the Drude damping plus quadratic damping ∝ k² may be strong enough to broaden and flatten plasmonic spectra before higher-order dispersion comes in sight for larger k. By contrast, in an electron energy loss spectroscopy experiment on InSb [334], i.e. with electrons rather than light, the measurements did not agree well with the hydrodynamic Drude model, but remarkably well with a Boltzmann-Mermin model. So challenges are to explore this boundary between theories, to find out whether other electron energy loss spectroscopy (EELS) experiments can be found where the HDM is accurate, and vice versa to find examples in linear and nonlinear optics where the HDM ceases to be accurate.

**Applications.** It has been shown that plasmonic cavities of a few microns in size of low-carrier density InSb can operate in the nonlocal regime of THz plasmonics [95]. These tunable cavities may become preferable to metallic ones, with useful modes of operation even in the nonlocal regime. The recent progress in nonlinear free-electron response of doped semiconductors [233,234,335], especially their tunable and efficient operation at low carrier densities, has put nonlocal hydrodynamic theory in the spotlight and may lead to novel schemes and compact semiconductor devices for second- and third-harmonic generation. Here further progress depends on material science and design based on demanding simulations of systems that combine both nonlocal and nonlinear response.

## Future developments to address challenges

**Combine linear and nonlinear optics.** Evidence is mounting that the nonlinear second- and third-harmonic response of doped semiconductors can largely be attributed to free-carrier hydrodynamic response. This makes doped semiconductors interesting materials for tunable on-chip nonlinear devices that may lead to novel applications. On the fundamental side, as
discussed above, the size-dependent linear optics effects of doped semiconductor nanoparticles have been explained either as due to quantum confinement or as hydrodynamic. In the future, combining nonlinear with linear measurements may help to identify the hydrodynamic nature (or not) of the linear response, by studying whether the nonlinear response shows a hydrodynamic nature.

**Comparing optical with conductivity measurements.** The relation between Drude conductivity and corresponding dielectric function is well known, and that dielectric functions can be studied both by optical and electron transport experiments. Some researchers in optics have attributed the imaginary part of the nonlocal parameter $\beta$ for semiconductors to viscosity of the electron gas [214,337], while solid-state physicists claim that viscous electron transport has been elusive and was first observed in graphene, see Ref. [342]. In the future, it would be interesting to do both electron transport and optical measurements on the same tunable plasmonic semiconductor systems, since the viscosity of the electron gas should affect both.

**Comparing optical with EELS measurements.** Above we discussed that the hydrodynamic Drude model could be used to explain nonlocal response in measurements on thin semiconductor films, while in some EELS measurements the HDM was inadequate, clearly not as accurate as a Boltzmann-Mermin model [334]. It would be worthwhile to explore the limits of the HDM, by doing combinations of optical and EELS measurements on the same planar structures.

## Concluding remarks

Here we have given an overview of some predictions and observations of nonlocal effects in plasmonic semiconductors, both in linear and nonlinear optics, and listed some challenges, opportunities and suggestions for future research.

**Funding.** Danmarks Grundforskningsfond (DNRF147).

# 15. Nonlocality in polar dielectrics

SIMONE DE LIBERATO ID

## Current status

Phonon polaritons, hybrid modes formed by the strong coupling of electromagnetic fields with optical phonons in polar crystals, have emerged as a powerful platform for sub-diffractional mid-infrared (mid-IR) photonics. Unlike plasmons, which rely on free carrier oscillations and can suffer from high losses, phonon polaritons benefit from the intrinsic low-loss nature of lattice vibrations. As such, they have found application in a range of areas including sensing, nonlinear optics and thermal emitters [343]. The standard approach to modeling polar crystals relies on a local dielectric approximation, which implicitly assumes a negligible momentum dependence of the phonon response. This local treatment works well at larger feature sizes, as in bulk transverse optical (TO) and longitudinal optical (LO) phonons have well-defined frequencies at optical wavelengths, separated by the Reststrahlen band. However, at the nanoscale, extreme light confinement can probe the phonon dispersion well inside the Brillouin zone, phonons acquire a non-negligible wavevector dependence, and their propagation within a confined geometry cannot be neglected. Under these conditions nonlocal effects, where the material response at a given point depends on fields in a finite surrounding region, become significant. Studying the electrodynamics of the system requires then taking into account the phononic, mechanical degrees of freedom on-par with the electromagnetic ones. This boils down to the Poynting vector, describing energy fluxes, acquiring a second term, proportional to the stress tensor of the lattice and describing energy propagating in the form of elastic energy in the solid [344]. Imposing the continuity of such a generalized Poynting vector allows one to determine boundary conditions on the mechanical fields that, together with the usual Maxwell boundary conditions, make it possible to
fully solve the coupled light-matter problem. One important feature of these boundary conditions is that they mix longitudinal and transverse degrees of freedom, leading to a modified optical response and the emergence of novel excitations known as longitudinal-transverse polaritons (LTPs), which mix both transverse and longitudinal character. While nonlocal effects are known to be present also in plasmonic systems, and some consequences of nonlocal phenomenology are similar in plasmonic and phononic systems, the origin and consequences of nonlocality in those two leading nanophotonic platforms is fundamentally different, as schematically shown in Fig. 19 [303]. Plasmons disperse toward the blue and longitudinal plasmons, which only exist above the plasma frequency, cannot thus hybridize with electromagnetically localized plasmon-polariton modes, which instead exist only below the plasma frequency where the dielectric function of the metal is negative. Their interaction is necessarily an evanescent phenomenon, leading to nonlocal phenomenology in plasmons being limited to an Ångström-thin skin depth. Optical phonons disperse instead toward the red, making the coupling between longitudinal and transverse degrees of freedom a resonant, propagative phenomenon, leading to the emergence of LTPs, and enhancing by orders of magnitude the length scales on which nonlocal phenomena can be observed and exploited. The existence of LTPs has been theoretically predicted and experimentally observed in silicon carbide (SiC) nanopillar arrays [345] and in crystal hybrids, polar superlattices with nanometric-sized layers in which standard dielectric modeling fails to reproduce even the qualitative features of experimental reflectance [303]. This success highlights that nonlocality in phonon polariton systems is not a theoretical curiosity but a practical concern that can affect device design and operation. Other predicted impacts of nonlocality on phonon polariton systems are an increased linewidth of Fröhlich resonances in dielectric nanoparticles [303] and the reduction of field confinement in nanogap resonators [346], as in both cases energy can leak out via emission of LO phonons.

## Challenges and opportunities

As the length scales of photonic devices approach the phonon propagation length, phonon nonlocality introduces several fundamental and practical challenges. The first and perhaps most direct challenge lies in accurate theoretical modeling. Traditional modeling approaches rely on local dielectric functions fitted to bulk optical constants, which fail to capture the momentum-dependent response of optical phonons. As a result, standard simulation tools and design methodologies yield inaccurate predictions of device resonance frequencies, confinement factors, and field enhancements. This discrepancy complicates the design of ultra-thin and deeply subwavelength structures where the effects of nonlocality are most pronounced. Multiple approaches have been explored, taking explicitly the phononic field into account either in finite-element method (FEM) simulations [346] or in scattering-matrix codes [347], relying on effective medium theories [348], or developing analytical mappings to well-known quantum optics models [349]. Still, these methods have been applied to relatively simple materials and interfaces, with the most complex case a single asymmetric Reststrahlen band in Ref. [347]. For more complex and potentially highly anisotropic materials, in which multiple phonon modes are present in overlapping Reststrahlen bands, the resulting set of coupled boundary conditions could be substantially more challenging to solve. Even for simple materials current approaches rely on approximations which limit the theory's overall domain of applicability. For a start, all the theoretical development cited rely on a quadratic approximation of the phonon dispersion. This means that the frequency of each optical phonon mode is written as a function of the wavevector as $\omega(k) = \sqrt{\omega^2(0) - \beta^2k^2}$, where $\omega(0)$ and $\beta$ are the zone-center frequency and velocity parameter of the mode. This can be problematic when considering interfaces between materials with Reststrahlen bands with broad overlaps, in which the frequencies provided by the quadratic approximation can diverge significantly from their physical values. While a completely numerical approach that takes as input a parametrization of the phonon dispersion across all
<figure><img src="image_21.png" /><figcaption>**Fig. 19.** Comparison of nonlocality in plasmon-polariton (top) and phonon-polariton (bottom) systems. On the left column we show a sketch of the surface polariton dispersions (yellow solid line) with highlighted the regions in which the Drude (top) and Lorentz (bottom) dielectric functions change sign. The red solid lines are the light line and the horizontal lines mark the plasma frequency (P) and the LO and TO phonon frequencies. The dashed blue lines sketch the dispersions of the longitudinal modes, written in the quadratic approximation in the central column for each case. On the right column we show the Fröhlich resonance, of frequency $\omega_F$, for a gold (top) and SiC (bottom) nanosphere, with the extinction cross section calculated using a local (red dashed line) and nonlocal (blue solid line) theory. The shaded region corresponds to the frequency region spanned by the dispersive mode for real values of the wavevector $k$. Only for phonon-polaritons, due to the red-dispersion of the LO phonon, such a region overlaps with the region of negative dielectric function in which the Fröhlich resonance can exist, leading to resonant energy exchange and to a substantially larger nonlocal length scale $l_{nl}$. Data adapted from Ref. [303].</figcaption></figure>

the Brillouin zone could in-principle be used to solve the problem in a more general way, this remains an open problem. Moreover, for extremely small nanolayers (1–2 nm), the microscopic crystal structure starts to play a role, and the thickness of the each layer, used as a parameter in the dielectric continuum modeling, is not precisely defined. This adds a material-dependent adjustable parameter to the theory which could be fixed (once each material interface) by ab initio simulations [303]. A second challenge is related to experimental characterization. Probing nonlocal phonon effects requires techniques capable of exploring few-nanometer scales, where these phenomena manifest. While various mature approaches to near-field microscopy exist, the dependency of the nonlocal effects on the exact dimensions of the objects requires samples with extremely high uniformity, or the capability to probe a single nano-object, which remains to-date challenging. Despite these hurdles, the opportunities that nonlocal phonon effects unlock are significant. Nonlocality offers a pathway to create novel types of hybrid modes like LTPs which combine the radiative nature of transverse fields and the strong interactions with electrical currents of longitudinal phonons, potentially bridging the gap between electrical excitations and far-field mid-IR emission. Such a capability could enable direct electrical pumping of mid-IR
phonon polariton modes, offering a route to electrically driven emitters without the need for engineered electronic transitions such as those in quantum cascade structures. In such a scheme, the Fröhlich interaction, responsible for LO-phonon emission in polar dielectrics, can become dressed by the strong interaction crating the LTPs. The scattering of a conduction electron would then resonantly create an LTP which could subsequently decay emitting mid-infrared light in the far-field. While initial theoretical analysis shows this mechanism to be potentially intense enough for practical applications [350], the intrinsically multi-scale nature of the problem makes it difficult to obtain precise results and an experimental verification of the emission mechanism is still missing. Only recently some more advanced non-equilibrium results obtained with a non-equilibrium Green function approach have appeared [351] and still only for simple structures and low values of the applied bias. Additionally, nonlocal responses could be used to enhance control over optical anisotropy and hyperbolicity in polar materials. By carefully designing nanostructured heterostructures, superlattices, and patterned nanophotonic elements, it may be possible to engineer the dispersion of phonons and thus the properties of the resulting polaritons. This could lead to enhanced sensing capabilities, ultra-confined modes with engineered dispersion, and new opportunities for nonlinear optical processes.

## Future developments to address challenges

To effectively harness phonon nonlocality, the next years will likely see a concerted effort along two main axes: theoretical modeling and experimental characterization. From a theoretical standpoint, developing comprehensive and user-friendly computational frameworks that include phonon dispersion and satisfy both electromagnetic and mechanical boundary conditions is a priority. Existing macroscopic models, some already demonstrated in simplified geometries, must be extended and made accessible to a broad community. Incorporating nonlocality into commercial electromagnetic simulation packages will be an important milestone, facilitating widespread exploration. Such models should be flexible enough to handle complex anisotropic materials, arbitrary nanostructures, and large-scale integrated photonic devices. On a fundamental level, more intuitive analytical approaches to interpret nonlocal behavior and guide rational design strategies will also prove valuable. Advanced multi-scale modeling of LTP-based electroluminescent devices will have to advance, allowing to better guide the design of prototype devices to gain a first experimental proof of this novel electroluminescence emission channel. Experimentally, there is a need for the characterization of more LTP materials, which will play the double role of both verifying and pushing forward theoretical developments. Spectroscopic characterization of monodisperse polar nanosphere of decreasing dimension could provide evidence for some of the yet unobserved predictions of the nonlocal theory as well as clarifying the behavior of materials below critical nanometer dimensions. Experiments are necessary to try and observe LTP electroluminescence. Only after a first unequivocal observation of this emission channel an effort to optimize the emission and collection efficiency could take place, which in turn will be necessary to ascertain whether LTP electroluminescence is at most an interesting curiosity or if it has the potential to empower a novel generation of mid-infrared optoelectronic devices. Finally, synergy with other fields—such as ultrafast optics, quantum photonics, and topological photonics—may lead to novel concepts and applications. For instance, the interplay of nonlocal phonon effects with coherent control schemes or coupling to quantum emitters (like defects in diamonds, molecular vibrations or quantum wells) could yield unprecedented control over mid-IR light-matter interactions. Similarly, integrating nonlocal phononic elements into topological platforms might provide robust, loss-resistant channels for infrared light. In short, the future developments needed to address the challenges posed by phonon nonlocality involve building a comprehensive toolkit encompassing theory, experiment, and device engineering. These advances will enable the community to fully realize the potential of nonlocal phonon polaritons in creating next-generation mid-IR photonic technologies.
## Concluding remarks

Phonon nonlocality represents both a theoretical and technological frontier in mid-infrared photonics. By moving beyond the local dielectric approximation, we acknowledge the true complexity of polar materials and unlock new phenomena such as longitudinal-transverse polaritons. Although this introduces modeling challenges and demands more refined experimental techniques, it could open up the possibility of creating more compact, efficient, and versatile mid-IR devices. The ability to couple electrical currents directly to radiative modes, leverage nonlocal effects for enhanced field confinement, and engineer the phonon dispersion through nanostructuring and heterostructures suggests a rich landscape of future innovations. As theoretical tools mature and integrate with experimental platforms, and as fabrication techniques continue to improve, phonon-based nonlocal photonics will likely play an increasingly important role. Ultimately, this may enable a new generation of mid-IR emitters, detectors, sensors, and nonlinear optical devices that harness the full power of phonon physics and broaden the scope of photonic materials and metamaterials research.

**Funding.** Leverhulme Trust (RPG-2022-037).

# 16. Nonlocal effects in graphene plasmonics

P. A. D. GONÇALVES & F. JAVIER GARCÍA DE ABAJO

## Current status

Plasmons in graphene and graphene-based nanostructures possess extraordinary properties, featuring deeply subwavelength field confinement ($\lambda_p \ll \lambda_0$, where $\lambda_p$ and $\lambda_0$ denote the graphene-plasmon and photon wavelengths, respectively), relatively low losses, and active tunability through electrostatic gating [352,353] and exposure to magnetic fields [354]. Graphene plasmonics evolved into a field of its own, fueled by the demonstration of gate-tunable plasmons that combine high confinement with low losses [355,356], an exceptional performance as ultrasensitive biochemical sensors [357], strong nonlinearities promoted by graphene plasmons [358], and potential for the design of quantum optoelectronic circuitry [359], just to mention a few directions. Additionally, and arguably just as important, the accelerated growth of graphene plasmonics played a pivotal role in paving the way for the broader field of two-dimensional (2D) polaritonics in atomically thin materials and related van der Waals (vdW) heterostructures [360].

In addition to their remarkable plasmonic properties, graphene plasmons (GPs) can exhibit strong signatures of nonlocal effects (i.e., arising from the material's *q*-dependent response, where *q* is the in-plane plasmon wave vector), in contrast to plasmons in three-dimensional (3D) or even ultrathin noble metals, where nonlocality is inherently weak and challenging to probe—often requiring nanostructures with characteristic length scales in the few-nanometer range [127, 252, 256], nanometric graphene-emitter separations [59], or the ability to measure large momentum transfers [361]. The unparalleled susceptibility of GPs to nonlocal effects can be largely attributed to two key properties:

1) *Two-dimensional character and unique band structure.* Graphene is a 2D semimetal in which charge carriers behave as massless Dirac fermions, following a linear energy–momentum dispersion. This enables on-demand control of the carrier density *n* (allowing for a complete on/off switch of graphene’s plasmonic response) and, consequently, control over the electronic Fermi wave vector via $k_F = \sqrt{\pi n}$. For moderate doping, the associated Fermi wavelength $\lambda_F = 2\pi/k_F$ can take values of tens of nanometers (Table 1), which also defines the length scale (e.g., for the plasmon wavelength) below which nonlocal effects are important. Moreover, in finite-sized nanostructures, $\lambda_F$ can become comparable to the system size *D*, making both nonlocal and atomistic effects appreciable.
2) *Large optical confinement.* Plasmons in graphene-based systems can achieve extremely large wave vectors $q = 2\pi/\lambda_p$, which can reach a substantial fraction of the Fermi wave vector $k_F$, ultimately approaching the $q/k_F \approx 1$ regime in graphene–dielectric–metal heterostructures, thereby enhancing the role of nonlocality.

Therefore, the combination of large Fermi wavelengths and highly confined optical modes renders graphene plasmons a unique platform for exploring the material’s nonlocal response. Beyond its significance for accurately modeling plasmonic resonances in graphene structures [353,362,363], nonlocal effects can also be harnessed to uncover many-body interactions encoded in the nonlocal optical conductivity of this carbon material [92].

<table><caption>**Table 1.** Key physical quantities underpinning the electronic properties of graphene (for three different carrier densities) and gold. (*)Taking the electron density from Table 1.1 in Ref. [24] and assuming an effective electron mass $m = 9.109 \times 10^{-31}$ kg.</caption><thead><tr><th></th><th>n</th><th>E<sub>F</sub> (eV)</th><th>k<sub>F</sub> (nm<sup>-1</sup>)</th><th>λ<sub>F</sub> (nm)</th></tr></thead><tbody><tr><td rowspan="3">Graphene</td><td>10<sup>11</sup> cm<sup>-2</sup></td><td>0.04</td><td>0.06</td><td>112</td></tr><tr><td>10<sup>12</sup> cm<sup>-2</sup></td><td>0.12</td><td>0.18</td><td>35.4</td></tr><tr><td>2 × 10<sup>13</sup> cm<sup>-2</sup></td><td>0.52</td><td>0.79</td><td>7.9</td></tr><tr><td>Gold(*)</td><td>5.9 × 10<sup>22</sup> cm<sup>-3</sup></td><td>5.53</td><td>12</td><td>0.52</td></tr></tbody></table>

Figure 20(a) shows the dispersion of plasmons in a free-standing graphene monolayer. In the long-wavelength limit (i.e., within the local-response approximation), graphene plasmons exhibit a $\omega_{GP} \propto \sqrt{q}$ dispersion, just like any homogeneous 2D electron gas (2DEG) (although with a different dependence on the carrier density, that is, $n^{1/4}$ instead of the $n^{1/2}$ scaling in conventional 2DEGs) [352,364]. As the plasmon wave vector approaches a significant fraction of $k_F$, nonlocal effects induce a blueshift in the dispersion. At very large wave vectors, the plasmon mode eventually enters the electron–hole continuum, enabling Landau damping [observed as a pronounced broadening of the plasmon band in Fig. 20(a)].

In extended graphene, the large optical confinement provided by graphene plasmons can be pushed to the ultimate limit in graphene–dielectric–metal (GDM) heterostructures [92,93,365]. In this configuration, screening from the metal at small graphene–metal separations significantly alters the plasmon spectrum, reshaping it into a nearly linear dispersion [Fig. 20(b)] [92,316,366]. As a result, these so-called acoustic graphene plasmons (AGPs) attain even larger wave vectors than those of conventional (i.e., unscreened) graphene plasmons at the same frequency. This, together with the substantial slow down of the plasmon group velocity, makes AGPs particularly susceptible to nonlocal effects [Fig. 20(b)]. Notably, Lundeberg et al. [92] measured the AGP dispersion using scattering-type scanning near-field optical microscopy (s-SNOM), revealing not only significant nonlocal effects but also electron–electron interactions, including Fermi-velocity renormalization and compressibility corrections. Furthermore, AGPs are capable of squeezing electromagnetic radiation down to one-atom-thick regions [93,365].

The significance of nonlocality in graphene plasmonics is not limited to spectroscopic measurements of nonlocal spectral shifts and broadening. The interaction between emitters and graphene plasmons, mediated by near-fields with large wave vectors, presents a promising avenue for probing nonlocal phenomena through their influence on emitter dynamics and emission spectra. In particular, graphene plasmons have been predicted to enhance otherwise “forbidden” light–matter interactions—e.g., multipolar transitions, two-plasmon spontaneous emission, and singlet-to-triplet transitions—by several orders of magnitude [367]. For instance, in the presence of large confinement factors ($\lambda_0/\lambda_p$ or $q/k_0$) achievable with AGPs, an atomic transition involving a change in angular momentum of $\Delta l = 5$ (a 5th-order electric multipole transition) can occur at a rate exceeding that of a dipole-allowed transition ($\Delta l = 1$) for the same atom in free space [367]. In this regime, accounting for nonlocal effects is not just important, but essential for
<figure><img src="image_22.png" /><figcaption>**Fig. 20.** Plasmons in extended graphene systems. (a) Plasmon dispersion in free-standing graphene obtained using different conductivity models, namely, the Drude approximation as well as the local and nonlocal versions of the random-phase approximation (RPA). (b) Local- and nonlocal-RPA dispersion of acoustic graphene plasmons in hBN-encapsulated graphene placed on a metal substrate (which behaves as a perfect conductor in the THz range). Landau damping is active in the hatched region ($\omega < qv_F$). The thickness of the top hBN slab is fixed at 10 nm, whereas the graphene–metal separation $t$ (controlled by the thickness of the bottom hBN film) is varied as indicated in the figure. hBN is modeled as a uniaxial crystal with $\epsilon_{xx}^{hBN} = \epsilon_{yy}^{hBN} = 6.7$ and $\epsilon_{zz}^{hBN} = 3.56$ (z along the normal to the graphene surface). For graphene (in both (a) and (b)), the carrier density and relaxation time are $n = 10^{12} \text{ cm}^{-2}$ (corresponding to $E_F = 0.12 \text{ eV}$) and $\tau = 500$ fs, respectively.</figcaption></figure>

rigorously describing plasmon–emitter interactions mediated by strongly confined graphene plasmons [367,368].

Additionally, the unconventional Dirac-cone electronic dispersion of graphene gives rise to a large intrinsic nonlinear response. Combined with the strong field confinement and enhancement associated with graphene plasmons, this makes graphene an attractive platform for nonlinear optics [358]. This potential can be leveraged to develop all-optical nanophotonic devices and implement quantum gates [359,369]. Incidentally, although graphene is a centrosymmetric material and, thus, not expected to exhibit second-order nonlinearities, the tightly confined graphene plasmons can give rise to spatially nonlocal nonlinear optical interactions. These interactions make the second-order response finite, with the corresponding coupling strength scaling as $\sim (q/k_F)^{7/4}$ [166].

In graphene nanostructures—such as ribbons and disks—nonlocal effects can be particularly pronounced due to the additional localization provided by the finite lateral size (i.e., the structure acts as a resonant cavity). This results in plasmon resonances with wavelength comparable to the characteristic system size $\lambda_p \sim D$, and introduces an *effective* wave vector $q \sim D^{-1}$. As mentioned earlier, and as shown in Table 1, the Fermi wavelength can satisfy $\lambda_F \sim \{\lambda_p, D\}$ for commonly used carrier densities. At low carrier concentrations, where $\lambda_F \gtrsim 100$ nm, nonlocal effects are observable even in relatively large graphene structures. In small nanostructures with significant edge-to-bulk ratios, nonlocality and quantum finite-size effects become closely intertwined, as localized plasmon resonances can probe the length scales associated with electronic edge states. Specifically, the type of edge termination—zigzag (ZZ) or armchair (AC)—leads to distinct plasmonic responses in nanostructured graphene (Fig. 21) [362]. Figures 21(b) and 21(c) show the plasmon energies and linewidths obtained through RPA calculations using tight-binding eigenstates for different ribbon widths. The classically predicted plasmon energies agree well with those obtained from quantum-mechanical calculations down to $D \approx 10$ nm. For smaller
ribbon widths, the atomistic theory predicts an increasing blueshift of the plasmon resonances as the width decreases. Note, however, that while the resonance position is almost identical for both ZZ- and AC-terminated ribbons (except for extremely narrow widths), the corresponding linewidths are significantly different: the plasmon linewidths in AC ribbons closely follow the classical result, whereas those in ZZ ribbons exhibit increased broadening for $D \lesssim 17$ nm. This broadening arises from edge states forming a weakly dispersive band near the Dirac point, enabling plasmon decay into electron-hole pair excitations for plasmon energies above $E_F$ [362].

Unlike nanoribbons, graphene nanodisks cannot be constructed with purely ZZ or AC terminations, inevitably resulting in mixed edge types (see inset of Fig. 21). Consequently, their plasmonic response becomes polarization-dependent, as the termination varies along the perimeter. This is clearly shown in Fig. 21(c), together with the observed blueshift of the plasmon resonances when compared with the classical description. Additionally, the figure reveals several weaker resonances, attributed to higher-order multipoles as well as hybridizations with interband transitions (involving electronic edge states from ZZ-terminated regions).

For systems comprising more than one nanostructure, the Coulomb interaction between neighboring elements leads to a richer plasmonic response [370,371].

In passing, we emphasize that, while in Fig. 21 the classical electromagnetic description accurately predicts plasmon resonances in graphene ribbons and disks with $D \gtrsim 10–17$ nm, this threshold is not universal, as it strongly depends on carrier density. For the cases considered in Fig. 21, with $E_F = 0.4$ eV ($n \approx 1.2 \times 10^{13}$ cm$^{-2}$), the corresponding Fermi wavelength is

<figure><img src="image_23.png" /><figcaption>**Fig. 21.** Plasmons in graphene nanostructures. (a) Illustration of the two standard types of edge termination in graphene ribbons. (b) Dipole plasmon resonances in graphene nanoribbons featuring either armchair or zigzag edge terminations (for normally incident light polarized across the ribbon). The circles denote atomistic calculations (see Ref. [362] for details), while the solid line indicates the classical result using the bulk conductivity of homogeneous, bulk graphene. (c) Linewidths associated with the plasmon resonances in (b). (d) Plasmon resonances in graphene nanodisks (see inset). The circles correspond to quantum-mechanical atomistic calculations for incident plane waves with different polarizations (see labels) under normal illumination. The circles are scaled according to the strength of the resonances. The solid line indicates the dipole resonance energy obtained using a classical electromagnetic framework. In all cases, the graphene nanostructures are considered to be free-standing, the Fermi energy is $E_F = 0.4$ eV, and the relaxation time is $\tau = 411$ fs. Reprinted (adapted) with permission from Ref. [362] (Copyright © 2012 American Chemical Society).</figcaption></figure>
$\lambda_F \approx 10$ nm, which aligns well with this threshold. Hence, at lower carrier densities, nonlocal and quantum finite-size effects can manifest in larger graphene nanostructures.

## Challenges, opportunities, and future directions

Graphene plasmonics is now a well-consolidated field from both theoretical and experimental standpoints. However, nonlocal and quantum size effects in graphene nanoplasmonics remain less explored, despite their significance for exploring the fundamental limits of plasmon-enhanced light-matter interactions in graphene and for designing truly nanometer-scale graphene devices.

On the experimental side, many theoretical proposals have yet to be realized, including the unambiguous observation of quantum size effects in small graphene nanostructures (down to “molecular-sized” graphene [372]) and the exploitation of ultraconfined graphene plasmons to infer quantum nonlocal effects in metals [316] and in strongly correlated materials [373]. While the latter could be already pursued using current state-of-the-art s-SNOM techniques [92,356,374], the fabrication of high-quality, pristine nanographenes on dielectric substrates suitable for plasmonics remains a challenge, with stable systems still limited to crystallographically well-defined metallic surfaces.

On the theoretical front, several promising directions remain open. In extended graphene, further investigation is needed in the hydrodynamic-transport regime—the opposite of the commonly studied quasi-collisionless scenario—where charge carriers behave as a viscous fluid, giving rise to intriguing transport phenomena such as negative nonlocal resistance and current whirlpools [375,376]. The implications of these effects for graphene plasmons remain largely unexplored. Moreover, conventional plasmon modes cannot propagate at velocities below $v_F$ [i.e., the plasmon dispersion cannot enter the electron-hole continuum, as shown in Fig. 20(b)]. However, in the hydrodynamic regime, plasmon-like energy waves (“demons”) can exist and propagate with velocities as low as $v_F/\sqrt{2}$ [377,378], further enhancing nonlocal effects. In finite-sized graphene, many open questions remain in relation to the influence of edge reconstructions—which can stabilize the edges both energetically and mechanically [379]—as well as disorder and edge-passivating adatoms, all of which constitute uncharted territory. Additionally, graphene nanoribbons can exhibit spin-polarized edges and edge magnetism [380], suggesting the exploration of spin-polarized plasmons [381].

In parallel, the growing interest in twisted moiré quantum materials is expected to extend into plasmonics in moiré graphene. This field is already gaining attention, with the first experimental observations reported [382] alongside theoretical advancements [383,384]. Here too, nonlocal effects could not only shed light into the mechanisms governing the complex phases arising in moiré materials but also open new coupling pathways where quasiparticles can exchange momentum with the twist-induced superlattice reciprocal vectors.

## Concluding remarks

Nonlocal effects in graphene plasmonics arise from the unique interplay between highly confined plasmons and large Fermi wavelengths in doped graphene, making this material an exceptional platform for exploring nonlocal phenomena. While often considered detrimental due to their association with increased broadening, nonlocal effects can also be leveraged to access otherwise inaccessible phenomena—such as nonlocality-induced nonlinearities and new coupling pathways—or to provide additional insights into electronic correlations intrinsically linked to nonlocal properties. Given the advancements and future prospects outlined in this section, the outlook for nonlocal graphene plasmonics is promising.

**Funding.** European Research Council (101141220-QUEFES); Horizon 2020 Framework Programme (101017720 FET-Proactive EBEAM, 964591-SMART-electron); Ministerio de Ciencia, Innovación y Universidades (PID2021-126964OB-I00, CEX2019-000910-S); Centres de Recerca de Catalunya.
# 17. Nonlocal quantum gain, loss compensation, and plasmon amplification in graphene hyperbolic metamaterials

ORTWIN HESS ID & ILLYA TARASENKO ID

## Current status

Insight into nonlocal quantum gain represents a pivotal advancement in the development of graphene-based hyperbolic metamaterials (GHMMs), addressing longstanding challenges in loss compensation and plasmon amplification. GHMMs are a class of hyperbolic metamaterials (HMMs) that leverage the exceptional properties of graphene, a two-dimensional material known for its high carrier mobility, tunable optical conductivity, and strong confinement of electromagnetic fields [385–389]. Traditional HMMs use metal layers to achieve hyperbolic dispersion [Fig. 22(a)], but these materials suffer from significant Ohmic losses and limited tunability. In contrast, GHMMs replace metallic layers with single or multiple graphene sheets [Fig. 22(b)], enabling dynamic control over optical properties while minimizing losses [389,390].

<figure><img src="image_24.png" /><figcaption>**Fig. 22.** Graphene-based hyperbolic metamaterials (GHMMs). (a) Traditional hyperbolic metamaterials (HMMs) use conductive sheets of finite thickness, where an incident transverse magnetic (TM) wave couples to carrier oscillations in both the x- and z-directions. (b) Replacing these metal layers with single sheets of graphene results in graphene-based hyperbolic metamaterials, where carrier oscillations are confined strictly to the plane of the graphene layer. (c) Density plot, showing the imaginary part of the Bloch wave vector for an infinite Graphene-dielectric periodic structure, shown in (b). Here graphene is optically pumped to create inversion. Positive (red) values correspond to amplification, negative (blue) values correspond to attenuation. Black curve bounds the hyperbolic regime of operation, where $\epsilon_{x-eff} < 0$. Green line is the single sheet plasmon dispersion. Reproduced with permission from Ref. [385] (Copyright © 2019 American Physical Society).</figcaption></figure>

Hyperbolic dispersion is a key feature of these materials, allowing propagation of high-wave-vector modes and enabling subwavelength confinement of light. These properties are essential for applications in imaging, sensing, and nanolithography [364,391]. However, conventional models for graphene conductivity, such as Drude and Kubo formulas, assume locality and linearity. These approximations break down for excitations with large wave vectors, particularly near the Dirac cone, leading to inaccuracies in describing plasmonic phenomena. The introduction of nonlocal quantum conductivity models, derived using the random-phase approximation (RPA), has addressed these limitations. These models account for both interband and intraband transitions, providing an accurate description of the optical response of graphene [392–394]. This enables the calculation of regions of plasmonic amplification of waves, propagating in active GHMMs, as illustrated in [Fig. 22(c)].

The chapter highlights the transformative potential of nonlocal quantum gain in GHMMs. To theoretically grasp loss compensation and amplification in GHMM requires a nonlocal quantum
conductivity model [385,386,395,396]

$$ \sigma_s(k_x, \omega) = i\omega e^2 \frac{\Pi(k_x, \omega)}{k_x^2}. \qquad (59) $$

Here, $\Pi(k_x, \omega)$ is the propagator of the electron-hole pairs. In RPA first order and expressed in Lindhard format [9], the polarization function of the arbitrary nonequilibrium graphene sheet can, in turn, be expressed in the form [385,395,396]

$$ \Pi(n) = \Pi|_{T=0,\mu=0} + \int_{0}^{\infty} d\epsilon \left[ \frac{\partial \Pi(k_x, \omega)|_{T=0,\mu=\epsilon} n(\epsilon)}{\partial \epsilon} \right] + \int_{0}^{\infty} d\epsilon \left[ \frac{\partial \Pi(k_x, \omega)|_{T=0,\mu=\epsilon} \bar{n}(\epsilon)}{\partial \epsilon} \right], \quad (60) $$

where $n(\epsilon)$ and $\bar{n}(\epsilon) = 1 - n(-\epsilon)$ are the distribution functions for electrons and holes, respectively.

Optical pumping of graphene introduces gain by creating an inverted electron-hole plasma, facilitating stimulated recombination that amplifies surface plasmons [360,395–398]. This gain mechanism not only compensates for inherent losses but also enables a shift from loss-dominated to gain-dominated regimes of operation. Such dynamic tunability opens up new possibilities for GHMM applications, including sub-diffraction imaging, efficient light-matter interaction platforms, and quantum photonics [399,400].

Recent studies have demonstrated the potential of GHMMs to achieve low-loss imaging and high-resolution performance in nanophotonic devices. For instance, GHMMs have been integrated into subwavelength optical filters, tunable modulators, and plasmonic waveguides. Their ability to dynamically switch between hyperbolic and elliptic regimes further enhances their versatility [401,402]. Despite these advances, challenges such as gain stability, fabrication precision, and scalability continue to hinder the practical implementation of GHMMs. The field is now focused on overcoming these limitations to fully exploit the capabilities of nonlocal quantum gain.

## Challenges and opportunities

While nonlocal quantum gain has unlocked new capabilities in GHMMs, significant challenges remain that must be addressed to realize their full potential. One of the primary challenges is achieving stable gain mechanisms. Optical pumping of graphene introduces nonlocal quantum gain by creating an inverted electron-hole plasma, but this process can lead to instability, noise, and mode competition [399]. These effects arise when the gain exceeds a critical threshold, resulting in uncontrolled amplification. Stabilizing gain requires precise control over doping levels, chemical potentials, and excitation regimes. Without careful optimization, gain mechanisms may inadvertently degrade device performance.

Thermal and collision losses also present a significant challenge. Graphene’s Fermi level smearing at elevated temperatures and carrier scattering at high doping levels introduce energy dissipation, limiting the effectiveness of loss compensation. These effects are particularly detrimental in applications requiring high field confinement or long propagation lengths. Mitigating these losses will require innovative material design and system integration [403–405].

Fabrication challenges further complicate the development of GHMMs. The precise layering of graphene and dielectric materials at nanoscale dimensions is technically demanding and prone to inconsistencies. Current fabrication techniques, such as chemical vapor deposition and atomic-layer deposition, show promise but require further refinement to ensure reproducibility and scalability. Achieving uniformity across large-scale devices is essential for translating laboratory advances into practical applications.

Despite these challenges, the opportunities presented by nonlocal quantum gain are immense. Amplified plasmonic modes in GHMMs enable enhanced light-matter interactions, which are critical for applications in quantum optics, sensing, and nanophotonics [406]. For example,
GHMMs can serve as platforms for developing efficient quantum light sources, entangled photon generation, and nanoscale waveguides. The ability to dynamically switch between hyperbolic and elliptic regimes also makes GHMMs ideal for reconfigurable optical devices, including modulators, filters, and frequency-selective surfaces.

Another exciting opportunity lies in extending the operational bandwidth of GHMMs. Nonlocal quantum gain enables the propagation of high-wave-vector modes, which can be harnessed for broadband plasmonic devices. This is particularly relevant in applications such as super-resolution imaging, where GHMMs can overcome the diffraction limit and achieve unparalleled spatial resolution. Moreover, hybrid systems combining GHMMs with complementary materials, such as transition metal dichalcogenides (TMDs) or hexagonal boron nitride (hBN), can further enhance device performance and introduce new functionalities [360,400].

## Future developments to address challenges

To address the challenges facing GHMMs, significant advancements are needed in material engineering, device design, and fabrication techniques. One of the most promising avenues is the optimization of gain mechanisms. This involves tailoring doping levels, chemical potentials, and excitation regimes to maximize electron-hole recombination while minimizing noise and instability [407]. Advances in laser excitation techniques and chemical functionalization could enable precise control over carrier dynamics, enhancing the stability and efficiency of plasmon amplification.

Thermal and collision losses can be mitigated by exploring new materials and hybrid structures. For example, integrating graphene with high-thermal-conductivity materials, such as hBN, can help dissipate heat more effectively. Similarly, introducing nanostructured dielectric layers or low-scattering materials into GHMMs can reduce energy dissipation and enhance performance under realistic conditions [389,398].

Innovations in fabrication techniques will also be crucial. Emerging technologies such as 3D nanoscale printing, machine learning-assisted design, and layer-by-layer assembly hold promise for improving precision, reproducibility, and scalability. These approaches could enable the mass production of GHMMs with consistent optical and electronic properties, facilitating their integration into existing photonic and electronic platforms [400].

Hybrid systems represent another promising direction for future development. Combining GHMMs with complementary materials, such as perovskites or TMDs, can address losses while introducing new functionalities, such as enhanced light absorption or tailored optical responses. These hybrid systems could leverage the unique properties of each material to create devices with superior efficiency and versatility [360,399].

In the longer term, GHMMs are expected to play a transformative role in quantum technologies. By leveraging nonlocal quantum gain, researchers can develop advanced platforms for single-photon generation, entangled photon sources, and low-loss quantum waveguides. These capabilities will be critical for advancing quantum communication, sensing, and computing.

Ultimately, the goal is to develop energy-efficient, ambient-compatible plasmonic devices that can operate reliably across diverse applications. Overcoming the current limitations will require a multidisciplinary approach, combining advances in materials science, device engineering, and theoretical modeling. By addressing these challenges, GHMMs could become foundational components of next-generation photonic and optoelectronic technologies [30,385].

## Concluding remarks

The advancements in our understanding leading to non-local quantum gain have established graphene-based hyperbolic metamaterials (GHMMs) as a promising solution to overcome the intrinsic limitations of traditional hyperbolic metamaterials. By leveraging the unique properties
of graphene and incorporating precise models of non-local quantum conductivity, GHMMs have demonstrated their potential to address critical challenges such as loss compensation and plasmon amplification. These innovations have enabled significant breakthroughs in nanophotonics, quantum technologies, and advanced imaging systems [360,385,388].

The key takeaway from the discussed work is the role of nonlocal quantum gain in mitigating losses and amplifying plasmonic modes, thereby enhancing the overall performance of GHMMs. This is achieved through the optical pumping of graphene, which facilitates stimulated electron-hole recombination and unlocks new operational regimes. As a result, GHMMs offer tunable and dynamically reconfigurable properties, making them suitable for a wide range of applications, including quantum emitters, nanoscale waveguides, and low-loss photonic devices [391,400].

Despite these advances, the field faces several challenges, such as achieving stable gain, reducing thermal and collision losses, and overcoming fabrication limitations. Addressing these issues will require multidisciplinary efforts in materials science, device engineering, and theoretical modeling. The development of hybrid systems, integration with complementary materials, and the use of advanced fabrication techniques are promising pathways to mitigate these challenges [389,399,402].

Looking ahead, GHMMs are poised to play a transformative role in next-generation technologies. Their ability to dynamically manipulate light at the nanoscale opens new frontiers in quantum communication, high-resolution imaging, and reconfigurable photonic systems. Moreover, the integration of GHMMs into existing photonic and electronic platforms will enable the development of energy-efficient, scalable, and versatile devices [398,401].

In conclusion, the realization of nonlocal quantum gain in GHMMs has paved the way for substantial advancements in the field of nanophotonics. While challenges remain, the opportunities presented by these materials far outweigh the limitations. With continued research and development, GHMMs are set to redefine the boundaries of photonic and optoelectronic technologies, offering solutions to some of the most pressing demands of modern science and engineering.

**Funding.** Research Ireland via the METAQUANT Research Professorship Programme (18/RP/6236).

# 18. Surface electrodynamics of crystalline noble metals

JOEL D. COX & N. ASGER MORTENSEN

## Current status

The emergence of crystalline metals in plasmonics [408], which are often considered superior to their polycrystalline counterparts [409], has stimulated substantial efforts in the synthesis [410–412] and optical characterization [413,414] of these pristine samples to enable demanding experiments [85,415–417]. However, crystalline metals are also turning notable in relation to their surfaces. In particular, surface science has long established [418] that (111) noble-metal surfaces [Fig. 23(a)] host in-plane conductive surface states—often referred to as Tamm-Shockley (TS) surface states [419,420]—governed by a free-electron-like parabolic energy-momentum dispersion relation [Fig. 23(b)] and hosting a two-dimensional electron gas (2DEG) localized to the very surface. The TS surface state is formally characterized by a surface conductivity $\sigma_{\parallel}(\omega, \mathbf{q}_{\parallel})$, exhibiting both frequency and spatial dispersion. The impact of the TS 2DEG can be seamlessly incorporated into electrodynamics by modifying the boundary conditions to account for the surface conductivity. This approach establishes a connection to the Feibelman surface-response function (see also Sec. 12 in this Roadmap), as $d_{\parallel} \propto \sigma_{\parallel}$ [192].

The surfaces of crystalline noble-metal particles are typically comprised of multiple facets corresponding to the intrinsic crystal planes of the solid [421]. The facets are in turn the
<figure><img src="image_25.png" /><figcaption>**Fig. 23.** (a) SEM image of gold crystalline metal flake with indication of the (111) surface facet. Reproduced with permission from Ref. [413] (Copyright © 2018 Optica Publishing Group). (b) Schematic energy-momentum electron dispersion relation for the TS surface states and the corresponding valence band in Au(111). (c) Plot of the loss function at an Au(111) surface and d acoustic-like surface-plasmon energy-momentum dispersion relation associated with the TS surface states on (111) facets of noble metals. Reproduced with permission from Ref. [318] (Copyright © 2021 Optica Publishing Group).</figcaption></figure>

cause for intriguing morphology-dependent and polarization-dependent plasmonic resonances [318,422–424]. In contrast, large planar flakes are predominantly characterized by (111) surfaces [Fig. 23(a)], although their edges may naturally exhibit other facets [413]. The conductive properties of the (111) facet in crystalline flakes are now renewing interest in TS surface states, particularly in the realm of their electrodynamic response.

## Challenges and opportunities

The electrodynamic response of the homogeneous 2DEG (possibly embedded in general dielectric media), and in particular its supported plasmon excitations exhibiting a square-root dependence of the plasmon energy on the wave vector [425], are theoretically well explored. A straightforward approach is to derive $\sigma_{\|}(\omega, \mathbf{q}_{\|} \rightarrow 0)$ using a Drude model within the phenomenological relaxation-time approximation, employing a 2D electron density consistent with the common Fermi energy of the bulk and TS surface states. A key challenge remains in gaining a deeper understanding of the relaxation rate in the 2DEG and how it may differ from its bulk counterpart.

Turning to the plasmon excitations of the TS surface states, strong screening by the adjacent 3DEG of the bulk influences the plasmon dispersion [426], resulting in an acoustic-like relation $\omega(\mathbf{q}_{\|}) \approx v_{\phi}\mathbf{q}_{\|}$ [Fig. 23(c)]. This is analogous to the 2D mirror plasmons observed in graphene-on-metal structures [92], but represents the extreme case where the separation between the 2DEG and the metallic “mirror” approaches zero. This makes the phase velocity $v_{\phi}$ even lower, leading to longer wave vectors at the same frequency [Fig. 23(c)], posing a significant challenge for experimental observation in near-field setups. Angle-resolved electron energy loss spectroscopy (EELS) has identified 2D plasmons on similar surfaces [427,428], but with even lower phase velocities $v_{\phi}$. This makes them even more tightly confined to the surface than graphene plasmons are to a graphene monolayer [92,429].

In the context of the plasmonic response of the bulk 3D electron gas (3DEG), the surface-response functions are modified by the presence of the TS 2DEG. In addition to the Feibelman perpendicular parameter $d_{\perp}(\omega)$, there is now also a non-vanishing parallel component $d_{\|}(\omega) \propto \sigma_{\|}(\omega, \mathbf{q}_{\|} \rightarrow 0)$ [192]. Thus, contrary to any other facet where $d_{\|} = 0$ [192,312], the (111) surface of noble metals is characterized by $d_{\|} \neq 0$ [318]. A compelling question arises: are there any plasmonic signatures that could differentiate (111) surfaces [430,431] from other surfaces of the same metal? For faceted noble-metal nanoparticles, the $d_{\|}$ thus varies among its facets, which is predicted to reflect polarization-dependent responses for even the dipole resonances [318]. The definitive experimental investigation of this phenomenon remains challenging, requiring precise
control and characterization of particle morphology and orientation on the surface, along with accurate measurement of its polarization response, for instance, using dark-field microscopy.

Probing surface-response functions remains an experimental challenge [258], as their extraction often relies on comparing computational spectra with and without surface effects, which in turn critically depends on accurate experimental knowledge of morphology on the nanoscale. An alternative approach involves electrically gating a single nanostructure, where strong bulk screening confines added electrons to the surface, effectively altering the surface response while leaving the bulk response unchanged. This method recently enabled the first modulation of surface response in a single plasmonic nanoresonator [236]. While simple models qualitatively capture these results [236,432], a more precise theoretical and computational description remains elusive. In addition to rigorous determination of $d_{\perp}(\omega)$ and $d_{\parallel}(\omega)$ Feibelman parameters for charge-neutral (111) surfaces and other facets with the aid of *ab initio* methods that account for d-band screening, additional key challenges include understanding the mutual equilibrium distribution of added electrons between TS surface states on (111) facets and nearby bulk states on all facets of a nanostructure. While atomistic models and the random-phase approximation (RPA) have been used to study surface-response functions of (111) surfaces [318], further investigations with advanced methods of charged systems (contrary to the common charge-neutral scenario), such as time-dependent density-functional theory, are needed. Consequently, the experiments by Zurak *et al.* [236] outpace the current theoretical foundation.

## Future developments to address challenges

To address the outlined challenges, future developments could focus on several key areas. First, there is a need to refine theoretical models to gain a more precise understanding of the relaxation rates for TS states and their differences from bulk counterparts. This potentially involves the study of the coupling to phonons in surface-terminated systems. By enhancing these theoretical frameworks, researchers can better capture the complexities of TS surface-state dynamics and interactions with bulk states.

Second, advancing experimental techniques for probing surface-response functions is essential. This includes innovating experimental techniques that improve the accuracy of measurements. The highly confined acoustic-like plasmons of the TS 2DEG can be probed using scanning near-field optical microscopy (SNOM) and electron spectroscopy techniques, though both methods present technical challenges, particularly due to low signal-to-noise ratios. In this context, recent ellipsometry-based characterization of the surface response offers an interesting alternative [313]. The excitation of acoustic plasmons based on nonlinear wave mixing constitutes another appealing possibility that exploits nonlocal and nonlinear light-matter interactions (see Sec. 19 in this Roadmap). Here, the energy- and momentum-matching of impinging light to highly confined polaritons could be achieved by interfering counter-propagating beams in a second-order difference-frequency generation scheme, such as that used to probe graphene plasmons from free space by Constant *et al.* [433]. Additional optical momentum can be provided by coupling light from dielectric waveguides to highly confined acoustic plasmon modes. The polarization dependence of the second-order nonlinear response at (111) surfaces [414] and plasmonic field enhancement [434] can provide additional insight to characterize the TS 2DEG. For nanoparticles, achieving greater control over nanoparticle morphology, orientation, and polarization response characterization will be critical for obtaining reliable data on the effects of surface response and in particular the electrodynamic effects of TS surface states.

Another crucial area is the integration of computational and experimental insights. Future research should focus on designing computational models that can explain experimental results with high fidelity, particularly in the context of gated plasmonic nanostructures. These models and *ab initio* descriptions must address the equilibrium distribution of added electrons between
TS surface states and nearby bulk states, potentially allowing theoretical predictions to better align with observed phenomena.

Additionally, further development of electrically gated single-nanostructure experiments will enhance the understanding of TS surface responses. Exploring methods to dynamically modulate and tune plasmonic properties in these nanostructures will yield new insights into the behavior of plasmonic systems under varying conditions. In this context, the combination of electrical gating with the AC lock-in technique employed by Zurak *et al.* [236] is a promising method to extract weak spectral changes otherwise buried in noise.

Moreover, investigating the plasmonic signatures unique to (111) noble-metal surfaces compared to other facets is vital. Identifying these distinct polarization-dependent resonances can provide valuable information about the fundamental interactions at play in these materials. Incidentally, for delafossite metals, recent reports of nonlocal electrodynamic effects in the hexagonally-symmetric conducting palladium (Pd) planes of palladium cobalt oxide (PdCoO₂) reveal the importance of facets in the Fermi surface, suggesting analogous dependencies on Fermi-surface facets in reciprocal space [435].

Furthermore, there is a significant opportunity to explore transdimensional plasmonics [436] (see also Sec. 23 in this Roadmap), particularly the transition from three-dimensional (3D) to two-dimensional (2D) plasmonic behavior in few-atom thin films subjected to out-of-plane quantum confinement. In these circumstances, where distinguishing between surface and bulk properties becomes less straightforward, it may be interesting to revisit Sipe's selvage considerations [190].

Lastly, the similarities between graphene plasmonics and TS surface states could be advantageous, as graphene plasmons may serve as a tool to probe and reveal quantum phenomena in the 2DEG linked to TS surface states. For instance, Gonçalves *et al.* [316] proposed extending graphene-on-metal experiments [92] to investigate the interplay between nonlocal 2D plasmons (from both graphene and TS types) and their screening by the adjacent 3D electron gas of the bulk metal, which also responds nonlocally [437,438]. While impressive experiments have utilized polycrystalline metal films [92,93], the change to crystalline flakes with (111) surfaces may present a new paradigm when interfaced with graphene. The hybridization of graphene and TS-type plasmons could pave the way for new explorations of quantum effects in highly confined electron gases.

### Concluding remarks

Crystalline noble metals, particularly their (111) facets with Tamm–Shockley surface states, present opportunities to explore surface-specific phenomena, including nonlocal response and refined surface-response formalisms. Progress requires advancing theoretical models for relaxation rates, nonlocal effects, and their distinctions from bulk properties, as well as improving experimental methods to probe surface responses with high precision. Greater control over nanoparticle morphology and the integration of computational and experimental approaches, especially in gated nanostructures, are key. Investigating surface-state plasmonic properties, including their nonlocal behavior, and hybrid interactions with graphene plasmons could pave the way for deeper insights into surface electrodynamics.

**Funding.** Danmarks Frie Forskningsfond (0165-00051B); Danmarks Grundforskningsfond (DNRF165).

## 19. Concerted nonlinear and nonlocal electrodynamics in 2D materials

LINE JELVER ID, EDUARDO J. C. DIAS ID & JOEL D. COX ID

### Current status

Nonlocal electrodynamic phenomena in nanophotonics have been extensively studied in the regime of linear optics, where a one-to-one correspondence is maintained between the energy
and momentum carried by photons and the polarization they induce in matter. However, spatially-
extended nonlocal light-matter interactions become inherently more complex in the realm of
nonlinear optics, where the self-interaction of intense optical fields comprised of many photons
is mediated by their mutual strong coupling with matter. For instance, sharp variations in the
spatial distribution of near-fields can significantly impact the nonlinear response of individual
nanostructures [439], while engineered long-range interactions among subwavelength scatterers
in metasurfaces can give rise to a large collective nonlinear optical response [440–442].

The regimes of nonlocal and nonlinear electrodynamics overlap when high-intensity light
is confined on extremely small length scales, as summarized in Fig. 24(a). This situation is
uniquely embodied by polaritons—quasiparticles formed by the hybridization of light with
collective dipole-carrying matter excitations—in two-dimensional (2D) atomically thin or van der
Waals (vdW) materials. Here, reduced screening due to the vanishing thickness of 2D materials
pushes the near-field confinement of their supported polaritons to extreme levels [360,366,401],
leading to enhanced spatial dispersion and an intrinsically nonlocal response. In tandem, the
concentrated electromagnetic near-fields associated with 2D polaritons can effectively drive
nonlinear light-matter interactions that exhibit an involved dependence on both frequency and
spatial dispersion (see Fig. 24).

<figure><img src="image_26.png" /><figcaption>**Fig. 24.** (a) Map of nonlocal and nonlinear light-matter interaction regimes according to interaction strength and field confinement. (b) Examples of local nonlinear optical processes, where electronic transitions are independent of photon momentum (left), contrasted with momentum-dependent nonlinear light-matter interactions (right). (c) Illustration of plasmon-driven second harmonic generation in graphene nanoribbon heterostructures enabled by passive (nanostructuring, left) and active (charge carrier doping, right) mechanisms that shape near-fields to break inversion symmetry. (d) Proposed scheme to excite highly confined polaritons, e.g., image polaritons in a 2D vdW material/metal heterostructure (inset), from free space through nonlinear wave mixing.</figcaption></figure>

Although significant effort has been dedicated to studying the nonlocal behavior of polaritons
in 2D materials (see Sec. 16 in this Roadmap), such as through scanning near-field optical
microscopy (SNOM) imaging [92,356,443,444], most investigations have focused on the linear
response regime. On the other hand, studies of nonlinear polariton-driven light-matter interactions
in 2D materials have thus far been mainly restricted to graphene plasmonics [166,433,445–449].
Owing to its linear electronic band structure and atomic thickness, graphene plasmon resonances
that strongly confine electromagnetic near-fields can be induced and actively tuned by electrostatic
gating. The Dirac cone electronic dispersion relation of graphene also ensures that its free charge
carriers undergo highly anharmonic motion, leading to a large intrinsic nonlinear optical response.
These properties combined have motivated intensive research efforts in nonlinear graphene plasmonics to leverage the highly confined and actively tunable 2D plasmons for applications in nonlinear nanophotonics [358]. However, graphene plasmon resonances are generally limited by achievable charge doping levels to terahertz (THz) and infrared (IR) frequencies, and, albeit longer-lived than their noble metal counterparts, suffer from substantial losses [356].

Beyond isolated 2D materials, nonlocal light-matter interactions become particularly important in vdW heterostructures. In particular, interfacing a vdW material with a metal film behaving as a perfect conductor gives rise to *image polaritons*, formed by the self-interaction of 2D polaritons with their mirror image in the metal [450]. The extreme confinement of the associated near-fields is characterized by an acoustic-like plasmon dispersion relation that extends well beyond the light cone, making these so-called acoustic plasmon modes extremely sensitive to nonlocal effects. In the case of metal/graphene hybrids, the deviation of their experimentally observed behavior from predictions based on a local description allows these polaritons to probe the nonlocal response of the combined system [437]. Although the strong near-fields associated with image polaritons hold great promise for nonlinear optics, studies of these excitations have mainly focused on the linear response regime.

### Challenges and opportunities

In centrosymmetric crystals, the leading contribution to the second-order nonlinear optical response is proportional to the gradient of the electromagnetic field that breaks inversion symmetry, and thus is inherently nonlocal [166,451]. While this symmetry dependence impedes second-order nonlinear optical effects (e.g., second-harmonic generation, sum- or difference-frequency generation, and optical rectification) in bulk crystals, it presents an opportunity to control the nonlinear response by engineering spatially-extended light-matter interactions. This concept has been theoretically explored in graphene plasmonics, where the resonant excitation of plasmons in inversion-symmetric ensembles of graphene nanoribbons is predicted to drive a large second-order nonlinear response when the inversion symmetry is broken by tuning the ribbon doping charge densities independently [452], as illustrated in Fig. 24(c).

Meanwhile, the extreme near-field confinement associated with image polaritons in graphene-metal heterostructures should give rise to an intense optical nonlinearity that is highly sensitive to nonlocal effects. Preliminary results have been achieved with heterostructures of graphene interfacing noble metal ribbons, where the broken translational symmetry of the metal ribbons facilitates in- and out-coupling to the far-field in experiments on harmonic generation [448,453]. However, the large volume of metal in these systems presents challenges when attempting to distinguish the metal and graphene contributions, and is problematic due to the large Ohmic losses in noble metals.

Nonlinear optical effects can also facilitate the excitation of propagating polaritons in extended systems. For instance, the difference frequency and wave vector produced by the second-order mixing of intense counter-propagating electromagnetic fields can be tailored to couple with the polariton dispersion extending beyond the light cone, as schematically represented in Fig. 24(d). This concept has been experimentally demonstrated for graphene plasmons, enabling the excitation of plasmons by free space illumination [433], as well as the optoelectronic tuning of multiple waveguided plasmon modes sustained by a graphene sheet embedded in an optical waveguide [449].

The combined nonlocal and nonlinear interactions of waveguided polaritons offers exciting prospects for developing integrated quantum photonic devices. In particular, propagating quantized plasmons in graphene nanoribbons are predicted to exhibit unity-order coherent nonlinear interactions, including two-plasmon absorption [359] and $\pi$-phase shifts induced by the collision of counter-propagating plasmons [369]. However, these schemes impose demanding
structural and electronic quality requirements on 2D sheets and nanostructures, presenting a significant challenge to their experimental realization. Although recent advancements have enabled atomically precise control over the width and edge configurations of graphene nanoribbons [454–456], as well as in the fabrication of nanoscale transition metal dichalcogenide (TMD) ribbons and tubes [457,458], scalable and robust production methods for polariton-based devices remains a critical barrier. Addressing this challenge will require further development of scalable fabrication techniques to fully unlock the potential of 2D materials in optoelectronic applications.

### Future developments to address challenges

The unique optoelectronic properties of atomically thin materials open numerous avenues for applications in nanophotonics, while the strong field confinement of polaritons in 2D materials beyond graphene can help overcome current limitations in graphene plasmonics. For instance, the lack of crystal inversion symmetry in monolayer TMDs allows even-ordered nonlinear optical processes to occur in extended samples, which can be enhanced by excitons with valley-dependent degrees of freedom [459]. Although the low group velocities exhibited by exciton-polaritons limit their ability to propagate optical signals, nonlinear and nonlocal light-matter interactions in engineered heterostructures could alleviate this constraint while leveraging the high exciton coherence for emerging technologies [460]. The directional dependence of propagating polaritons in anisotropic 2D materials—such as plasmons in doped phosphorene or in emerging vdW materials like palladium phosphoselenide (PdPSe) and arsenic trisulfide ($As_2S_3$)—offers additional spatial control that can influence polariton-driven nonlinear optical phenomena. Notably, $As_2S_3$ features a record high birefringence [461], further expanding the design possibilities for advanced nanophotonic systems.

The greatest challenges for the realization of devices that harness concerted nonlinear and nonlocal effects lie in the need for highly precise fabrication methods that produce defect-free and well-aligned nanostructures while ensuring seamless integration with existing technologies. On the theory side, advanced condensed matter and electrodynamic models that incorporate spatial dispersion in the nonlinear response functions of extended 2D materials must be developed to accurately capture nonlocal effects [462]. In mesoscopic 2D systems, atomistic simulations offer a powerful framework to simultaneously incorporate nonlocal and quantum-mechanical effects in the polariton-driven nonlinear optical response [463,464].

The creation of heterostructures that integrate graphene with other 2D materials offers another promising direction for next-generation photonic and optoelectronic technologies. Such heterostructures offer highly confined polaritons, exhibiting enhanced in-plane dispersion and strong associated near-fields, and can leverage the high intrinsic nonlinear response and optoelectronic tunability of graphene to explore nonlocal and nonlinear optical phenomena. Future research could focus on graphene-hBN heterostructures, where hybrid plasmon-phonon modes enable ultrafast nonlinear switching and low-loss waveguiding in the mid-IR [465]. The exploration of graphene-phosphorene heterostructures will assess their potential for reconfigurable and anisotropic photonics in broadband signal processing [466]. Furthermore, hybrid structures comprised of hyperbolic 2D materials open new paths to engineer polariton dispersion, unlocking multifunctional photonic platforms [467,468]. Finally, advances in vertical integration and interlayer coupling could lead to more efficient nonlinear optical devices with dynamic control of both spectral and spatial responses.

### Concluding remarks

While low-dimensional materials are ubiquitous in nanophotonics, studies of spatially-extended nonlocal optical phenomena have largely been limited to the linear response regime. The intense near-field concentration associated with polaritons in 2D materials presents unique opportunities
to explore the interplay of nonlocal and nonlinear light-matter interactions. To this end, novel theory frameworks that incorporate optical nonlocality in the response functions of 2D materials are needed to describe the complex energy- and momentum-dependent nonlinear interactions of their supported polaritons. On the experimental front, advanced nanofabrication techniques are required to create high-quality 2D material heterostructures capable of probing and leveraging spatial and frequency dispersion in integrated nonlinear photonic devices.

**Funding.** Danmarks Frie Forskningsfond (0165-00051B); Danmarks Grundforskningsfond (DNRF165).

## 20. Nonlocal chirality in a twisted multilayer

MIGUEL SÁNCHEZ SÁNCHEZ, DIONISIOS MARGETIS, GUILLERMO GÓMEZ-SANTOS, IGOR M. VASILEVSKIY & TOBIAS STAUBER

### Current status

Etymologically, the word chirality comes from the Greek word for “hand” and is thus related to a three-dimensional object whose mirror image does not match the original by any translation or rotation in 3D space. In reverse, this means that purely two-dimensional systems cannot be genuinely chiral, since the object itself would serve as a mirror plane.

Chiral objects define two enantiomers which are identical in all scalar properties such as density or eigenfrequencies. Nevertheless, in their interaction with other chiral objects, opposite enantiomers can be detected. Circularly polarized (CP) light can thus distinguish left or right-handed samples by measuring the different absorption cross-sections of left and right CP light giving rise to circular dichroism (CD) [469]. Since linearly polarized light is the superposition of left and right CP light, there will also be a rotation of the polarization plane. This has been observed in neutral twisted bilayer graphene (TBG) for certain resonant frequencies related to transitions around the van Hove singularities [470].

Polarization rotation and CD can usually be observed by breaking either time-reversal or rotational/mirror symmetries [471]. In TBG, both symmetries are conserved, which makes the observation of CD remarkable since Maxwell’s equations are reciprocal in unbiased two-dimensional systems and there should be no dependence on the direction of the CP light. However, there is a slight breaking of the two-dimensionality of the system as the two layers are separated by a distance $a = 3.4Å$. In fact, the chiral response in TBG is due to the non-local correlation between the current-density of layer 1 in say x-direction and the current density of layer 2 in the (transverse) y-direction:

$$ \chi^{\text{chiral}}(\omega) = -\frac{i}{\hbar} \int_{0}^{\infty} dt \, e^{i\omega t} \langle [j_{x}^{1}(t), j_{y}^{2}(0)] \rangle, \quad (61) $$

where $j_{y}^{\ell}(t)$ is the v-directed current operator ($\nu = x, y$) in layer $\ell$ in the interaction picture, and $\langle \cdot \rangle$ is the equilibrium average. Equation (61) provides the average current in layer 1 due to the electric field in layer 2, and could be interpreted as the system response to the “gradient” of the electric field along the third dimension [5].

The chiral response only depends indirectly on the distance between the two layers, $a$, via the equilibrium average. However, typical chiral observables are directly proportional to $a$ and thus depend on $a\chi^{\text{chiral}}$ [472].

The nonlocal response between perpendicular current-directions of different layers is not the only requirement to observe CD, because for electron-hole symmetric systems the chiral response is identically zero at charge neutrality. Only due to small electron-hole symmetry breaking, electron-like and hole-like transitions are not compensated and the enhanced density of states
around van Hove singularities can give rise to an observable CD [473,474]. In fact, two main resonances are observed in Ref. [470], related to transitions around two different van Hove singularities.

### Challenges and opportunities

To calculate the CD, one needs to go beyond the dipole approximation $e^{ikz} \approx 1 + ikz$, again emphasizing that chirality is a three-dimensional phenomenon. This introduces the dimensionless scale $ka$, where $a$ is the typical extension of the chiral object and $k = 2\pi/\lambda$ the wave number of light. This scale is usually small making chiral effects hard to use in typical nanoscale devices.

Field-theoretically, the chirality of an electromagnetic field with electric field $E$ and magnetic field $B$ in a dielectric medium with relative dielectric permittivity $\epsilon$ and magnetic permeability $\mu$ is defined by

$$C = \frac{\epsilon\epsilon_0}{2} E \cdot (\nabla \times E) + \frac{1}{2\mu\mu_0} B \cdot (\nabla \times B), \qquad (62)$$

which can further be related to a chiral flux via a continuity equation [469]. For CP light, the chirality is proportional to its frequency $\omega$ and intensity $|E_0|^2$. For fixed wavelengths, the chirality can thus only be modified by changing the amplitude which can be achieved by confining light. In fact, decay into evanescent modes gives rise to fluorescence quenching which has been used to detect chiral molecules with the help of TBG [475].

Also, surface-plasmon polaritons could enhance the chiral coupling and may enhance the chirality up to 4000 times the one of corresponding propagating CP light [476]. Cavities composed of TBG should thus offer tremendous opportunities as plasmon-induced near-field chiralities may enhance photocatalytic processes for enantiomer-selectivity.

Quasi two-dimensional moiré structures can host plasmons which are inherently chiral [472]. Especially, many organic molecules are intrinsically chiral mainly due to the chemical structure of carbon with its four valence electrons. In fact, all amino acids are chiral as only glycine is achiral. Moreover, the RNA aptamer or Spiegelmer is chiral as well as hexahelicene and pentahelicene. All these may be synthesized in an enantiomer-selective way giving rise to new functionalities.

### Future developments to address challenges

The crucial requirement for the chiral response in moiré systems is the spatial separation between the layers and thus the nonlocal correlations between the sheet currents. However, increasing the layer distance would result in weakening the quantum mechanical coupling. Hence, there is a trade-off between optimal coupling and large layer separation, probably best met by transition-metal dichalcogenides such as tungsten diselenide (WSe$_2$).

Another route to enhance the inherent chirality is to stack several van der Waals materials on top of each other [477,478]. The chiral response of multilayers with width $d$ should thus be increased by the geometric factor $d/a = N$, where $N$ is the number of layers. However, the moiré unit cell grows exponentially with the number of layers even for commensurable twist angles, and for incommensurable structures it is difficult to obtain the electromagnetic response [479].

For the special case of alternating twisted multilayers, the analysis can be considerably simplified as it can be mapped to independent twisted bilayers and one single layer in case of an odd number of layers [480]. In Fig. 25, the multilayer moiré systems with $N = 2, N = 3$, and $N = 4$ layers are shown. The structures with an even layer number are inherently chiral, whereas those with an odd layer number are achiral as they display a mirror plane at the central layer. Nevertheless, even for $N = 3$ the nonlocal “chiral” correlations that involve the perpendicular current densities of adjacent layers give rise to vertical gradients of the magnetic moment as a response to an electric in-plane field [481]. This layer nonlocality in nominally “achiral” systems could be detected by layer-discriminated contacts, as shown in Fig. 25(b) and also discussed in Ref. [482].
<figure><img src="image_27.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>**Fig. 25.** Multilayer moiré systems with (a) *N* = 2, (b) *N* = 3, and (c) *N* = 4 layers. The structures with alternating twist angles and an even layer number are inherently chiral, whereas those with an odd layer number are achiral as they display a mirror plane at the central layer.</figcaption></figure>

Let us note that we can estimate the chiral response for moiré systems with *arbitrary* twist angles between the layers using the decoupling procedure proposed in Ref. [480]. For that, we rely on the fact that decreasing (increasing) twist angles can be mimicked by properly increasing (decreasing) the interlayer tunnel amplitude [483]. Neglecting effects of non-commensurability, we can thus map any multilayer moiré system onto a set of independent twisted bilayers and one single layer in case of an odd number of layers.

## Concluding remarks

The nonlocal coupling between sheet-currents flowing the perpendicular directions gives rise to a chiral response for neutral multilayer moiré systems without a mirror plane and slight electron-hole asymmetry. Also for multilayer moiré systems with a mirror plane, nonlocal “chiral” correlations can be detected by layer-discriminating contacts. These effects can be considerably enhanced by non-radiative (plasmonic) modes leading the path to chiral cavities composed of multilayer moiré systems with improved enantiomer-selective capabilities. The functionality can be further optimized by tuning the twist angles and doping levels.

**Funding.** Ministerio de Ciencia, Innovación y Universidades (PID2021-126964OB-I00, PID2020-113164GBI00, PID2023-146461NB-I00); Consejo Superior de Investigaciones Científicas (PTI-001); Spanish Ministry of Science,
Innovation and Universities through the “María de Maeztu” Programme for Units of Excellence in R&D (CEX2023-001316-M).

# Part IV—Engineered nonlocal responses in metamaterials and metasurfaces

## 21. Strong spatial dispersion in metamaterials and metasurfaces

### SERGEI TRETYAKOV & CONSTANTIN SIMOVSKI ID

#### *Overview*

Metamaterials (MMs) are commonly defined as arrangements of artificial structural elements, designed to achieve advantageous and unusual electromagnetic properties [484]. This definition implies that metamaterials behave as effectively homogeneous media. This property distinguishes MMs from photonic crystals and other arrays that cannot be homogenized and modeled by constitutive relations for macroscopic fields. Structural elements of MMs, usually called *meta-atoms* (MAs), as a rule, are small compared with the wavelength $λ$ and are arranged with the interparticle distance which is also much smaller than $λ$. However, some arrays of large or even infinitely extent MAs (thin wires in a dielectric ambient or thin metal layers alternating with thin dielectric layers), can be classified as MMs (*hyperbolic MMs*) because these arrays can be homogenized if their MAs are arranged with a subwavelength period [484].

Among many parameters of meta-atoms that can be varied to design and optimize metamaterials and metasurfaces, the key role is played by their sizes and shapes. Shape effects on electromagnetic response of objects become significant when the object size becomes comparable (not negligibly small) with the wavelength. Dependence of the material response on the wavelength is called *spatial dispersion*. SD means non-locality of the polarization response. Understanding SD and using its implications is in the core of MM science and technologies. The same refers to metasurfaces (MSs) which are the 2D analogues of MMs – thin composite layers [485]. Instead of bulk material parameters, MSs are characterized either by their surface impedance, or by effective electric, magnetic, and magneto-electric collective polarizabilities, or susceptibilities [485], which are in most cases spatially dispersive. SD, being a spatial non-locality, is manifested in both MMs and MSs as the dependence of the material parameters (bulk or surface ones) on the wavevector $k$. In the absence of spatial dispersion, metasurfaces cannot realize such important functionalities as, for example, full absorption or reflection phase control, so that the possible applications are basically limited to frequency-selective surfaces and some polarizers.

Weak SD is inherent to MMs composed by small complex-shaped MAs. It corresponds to weak $k$-dependence of the permittivity tensor $\varepsilon(\omega, k)$ which obviously depends also on the frequency $\omega$. The first-order terms in the Taylor expansion of $\varepsilon(\omega, k)$, linear in $k$, define such effects of weak SD as chirality and omega coupling [484,486]. Among the effects of the second-order SD the most important one is artificial magnetism [484,486]. Both first- and second-order SD effects are observed also in MSs [485]. For MMs artificial magnetism results in the relative permeability being different from unity in the absence of natural magnetics. For MSs this effect results in the nonzero magnetic surface susceptibility and (in the alternative model) in magnetic collective polarizability. Chirality in both MMs and MSs leads to polarization transformation and conversion [484–486]. Omega coupling results in asymmetry of reflection (for both MS and MM layers in homogeneous space) [486]. More information about weak SD in MMs can be found in Refs. [484,486,487]. Weak SD in MSs is reviewed in Refs. [488,489]. Strong SD corresponds to the case when the $k$-dependence of $\varepsilon(\omega, k)$ is resonant [490], or its Taylor expansion poorly converges, including too many terms [484]. Strong SD is inherent to hyperbolic MMs [484,490,491]. MSs with strong SD are specially engineered so that their nonlocal response allows the manipulation of electromagnetic field spatial distributions and desired functionalities [488,489].
## Current status of research

Currently, physics and applications of MMs and MSs with weak SD is a mature field. In particular, there are well-established means to create and engineer resonant chirality for control and optimization of optical activity dichroism. In addition to the use of chiral shapes, these means include also the use of arrays of negligibly thin composite layers that are non-chiral in the geometric sense (pseudochirality [486,488,489], also called extrinsic chirality). If an electrically thin metasurface does not contain natural magnetic materials (like ferrites, for example), the only mean to use it for full control of reflection and transmission is to invoke spatial dispersion. This allows us to complement electric polarization of thin arrays with artificial magnetic response of spatially dispersive structures and realize non-reflecting (Huygens') and more general arrays for control of reflection and transmission [485,488,489]. Omega coupling plays the key role in realizing perfect anomalous refraction in metasurfaces [488]. Also, all electrically thin absorbers covering highly reflecting bodies are spatially dispersive (omega) layers. Strong SD effects in MMs and MSs have been also studied recently. It was found that strong SD in MMs offers long-distance transmission of subwavelength-resolution images, super-Planckian thermal emission and absorption and super-Planckian radiative heat transfer (see in Refs. [490,491]). For MSs, strong SD grants strong coupling of propagating and evanescent waves and offers a unique opportunity to engineer diffraction patterns [489]. Furthermore, full control of reflection to realize anomalous reflection or focusing by MSs requires strong SD. Most recent studies begin exploring new electromagnetic phenomena in time-varying media with spatial dispersion. For example, phenomena at time interfaces in chiral media were studied in Ref. [492], and effects in time-varying media with omega coupling are discussed in Ref. [493].

## Challenges and opportunities

To design and optimize metamaterials and metasurfaces, appropriate modeling tools are needed. As already discussed above, if SD is weak, there are well developed effective-medium models based on the bianisotropic material relations. However, even the second-order effects are modeled by these relations only partially, in form of artificial magnetism. Other second-order terms as well as all higher-order terms enter the material equations with spatial derivatives of the macroscopic fields [494,495]. Solving boundary problems for MMs and MSs with strong SD is even more difficult. In this case, one needs to derive additional boundary conditions (ABCs) for each particular MM microstructure [487]. In Ref. [494], ABCs were derived within a higher-order effective-medium model. In Refs. [496,497], ABCs were derived for hyperbolic MMs. For media with strong SD difficulties arise not only in solving boundary problems. While for low-loss bianisotropic media with artificial magnetism there is a general formula for the energy density (see e.g. in Refs. [484,486]), for media and MSs with strong SD such general formulas are absent. In Ref. [498], the energy density in wire media (one of two known classes of hyperbolic MMs) was deduced from microscopic consideration (telegrapher equations), and the results were shown to be consistent with the effective-medium model (k-dependent permittivity tensor). The main goals of introducing effective medium models are to enable simple solutions of extremely involved electromagnetic problems of fields in complex composite media and, most importantly, get physical insights into the computed or measured results. Unfortunately, higher-order effective material models become more and more complex, so that both goals become difficult to reach. In view of these difficulties, it appears that direct numerical optimization of metamaterial and metasurface structures remains the only available tool for development of advanced metadevices. Significant efforts have been devoted recently to development of tunable and reconfigurable MSs. One of the main envisaged applications is in future wireless communication technologies, where tunable MS can be potentially used to optimize propagation channels in wireless communications. As discussed above, only spatially dispersive MSs can offer full control of anomalous reflection
and transmission, which means that we need to develop effective and practical means to tune
non-local responses of MAs and their arrays or clusters electrically or optically.

*Future development to address challenges*

At present, research efforts in this field are focused on the design and optimization of MSs for control of reflection and transmission, including focusing. In the microwave and terahertz ranges the main envisaged applications are in future wireless communications (reconfigurable intelligent surfaces). In optics, researchers envisage, for example, extremely thin MS lenses. As discussed above, for any advanced control of waves, MSs must exhibit SD. Moreover, realizations of the optimal performance often demand specially engineered strong SD. For MSs, this usually means careful control of excitation and surface distribution of evanescent (surface) waves. Their fields exhibit fast variations at the wavelength scale, which means that MA sizes should be smaller than the wavelength scale. To realize such designs, technologies to create small, simple, and cheap MSs need to be developed. To make reconfigurable and adaptable devices, the properties of these MAs should be electrically or optically tunable in such ways that will allow optimization of meta-atom couplings. To optimize parameters, it is necessary to further advance numerical optimization methods, because fully reconfigurable devices cannot be realized as periodic arrays, so that a global optimization of huge sets of interacting meta-atoms will be needed. One possible route can be creation of standardized units (for each frequency range) and pre-computation of basic electromagnetic parameters (such as impedance matrices and effective antenna heights), after which fast purely arithmetic optimizations become possible [499]. Perhaps, some lessons from recent fast developments of artificial intelligence can be used in this field. Here, concerted efforts of experts in electromagnetics and computer science will be needed. To advance understanding and use of time-modulated (4D) MMs and MSs, very significant theoretic research is needed. So far, only initial studies of time-varying media have been made, mostly completely neglecting frequency and spatial dispersion. Studies of time-varying spatially dispersive media considered only weak spatial dispersion, where it was possible to use the bianisotropic effective medium model [493]. Moreover, always-present frequency dispersion of spatially dispersive materials was so far modeled only in the low-frequency approximation. Furthermore, no experimental studies of spatially dispersive media whose properties quickly change in time have been made so far. Here, we would like to remind that a huge majority of practically relevant MSs have significant SD! Many challenges are expected on this way: First, time-varying spatially dispersive media need to be studied in time domain. While the frequency-domain methods are well developed, this is not the case of the time-domain ones. For example, to solve for fields at time interfaces, time boundary conditions are needed, and to establish them, one needs to consider the microstructures of spatially dispersive materials. This means that no universal boundary condition (like continuity of fields in the frequency-domain theories) can be found. In summary, we expect seeing interesting and very much needed results of researchers working in this area.

**Funding.** European Research Council (PULSE, project No. 101099313); Research Council of Finland (project No. 365679).

**22. Nonlocal hyperbolic metamaterials and metasurfaces**

SAMANEH PAKNIYAT & J. SEBASTIÁN GÓMEZ-DÍAZ

*Overview*

Hyperbolic media, including hyperbolic metamaterials (HMMs) and metasurfaces (HMSs), have driven groundbreaking advancements in photonics as well as in related fields like acoustics and thermodynamics. Hyperbolic media display remarkable electromagnetic properties arising from
ith extreme anisotropy, hyperbolic dispersion, and field confinement. Such exciting properties also entail that the media response depends on the properties of the waves with which they interact, leading to stronger nonlocal effects than those found in other media. Nonlocal responses grow stronger with the confinement of the supported modes and appear due to the granularity of the artificial media and the spatial dispersion of the composing materials—for instance, the finite velocity of the electrons flowing in metals. To date, nonlocality has been mostly considered a hindrance in hyperbolic media, and thus it has not been exploited to enable new functionalities or to enhance device performance. This Roadmap explores the interplay of nonlocality and hyperbolic dispersion, highlighting both challenges and untapped opportunities. The limitations of conventional local effective medium theory are reviewed, emphasizing the need for sophisticated numerical tools able to solve general electromagnetic structures composed of nonlocal metallic elements. These tools will guide the development of hyperbolic devices in which nonlocality will be tailored to become an extra degree of freedom for wave control. Future directions include merging drift-current nonreciprocal plasmonic, moiré physics, and reconfigurability to realize tunable, nonlocal, and hyperbolic platforms with exciting applications in multidisciplinary fields, ranging from sensing, signal processing and computing to acoustics and thermal management.

## Current status

Hyperbolic media [500–502] exhibits extreme anisotropy enabled by a material response that changes sign as a function of the polarization of the electric and/or magnetic field and offers unique hyperbolic dispersion together with an ideally infinite local density of states and wave confinement [Fig. 26(a)]. Such features have driven advancements in ultra-high-resolution imaging, negative refraction, beam shaping, optical sensing, and spontaneous or coherent emission enhancement, among other applications. The electromagnetic response of hyperbolic media is limited in practice by two mechanisms, namely the presence of loss and nonlocality, which close their otherwise open hyperbolic dispersion relation and impose a wavenumber cutoff to the supported modes. The lossy nature of frequency-dispersive materials is a direct consequence of causality, and thus unavoidable in passive devices [5]. Nonlocality appears due to the granularity of artificial structures and to the intrinsic spatial dispersion of the materials that compose them [503,504], for instance, the finite velocity of electrons moving within metallic elements. Nonlocal media is usually described as wavenumber-dependent media and thus the material response depends on the properties of the waves with which it interacts. In addition to capturing the fundamental response of hyperbolic media and imposing upper limits to the field confinement and local density of states, nonlocal effects must be accounted for in the design of realistic hyperbolic devices and exploited to enhance their performance.

Hyperbolic responses appear naturally in certain materials and crystals, such as silicon nitrate, hexagonal boron nitride, α-phase molybdenum trioxide (MoO₃) bilayers, or β-phase gallium trioxide (Ga₂O₃) crystals, among others [505,506] and, as illustrated in Fig. 26(b), can also be implemented by arranging subwavelength unit-cells taking advantage of metamaterial concepts [500,501]. Hyperbolic media are usually characterized using local effective medium theory (EMT) based on the Maxwell-Garnett approximation coupled with a linear response model of the composing materials [500,501]. EMT describes the hyperbolic structure as a homogeneous medium with effective parameters and assumes that the structural elements are much smaller than the operation wavelength.

Interestingly, the response of hyperbolic media goes well beyond the predictions made by EMT. To illustrate nonlocal mechanisms, let us consider that a hyperbolic surface composed of graphene strips with periodicity *L* supports electromagnetic modes characterized with a momentum *k* and a wavelength *λ*. First, EMT remains valid only when the artificial structure can be considered homogeneous, i.e., when the condition *L* ≪ *λ* is satisfied. As the supported modes become
more confined, their momentum *k* increases and their wavelength *λ* decreases to a point in which they are comparable to the strips size (*L* ≈ *λ*). There, as shown in the left panel of Fig. 26(c), granularity becomes the main nonlocal mechanism that closes the isofrequency contour of the supported modes and imposes a cutoff wavenumber at *k*<sub>c</sub> ≈ π/L. Second, the response of metallic materials is usually modeled with local models (i.e., Drude local response approximation) that do not account for the finite speed of electrons. As the wavenumber of the supported modes increases, electrons might not be able to travel fast enough to follow spatial field variations. This imposes a wavenumber cutoff at *k*<sub>c</sub> ≈ (c/v<sub>F</sub>)*k*<sub>0</sub>, where v<sub>F</sub> is the Fermi velocity of electrons in the metal, *c* is the speed of light, and *k*<sub>0</sub> is the free-space wavenumber. The right panel of Fig. 26(c) illustrates this scenario, confirming that the intrinsic nonlocal response of graphene becomes dominant and closes the isofrequency contour when granularity is negligible. In a more general

<figure><img src="image_30.png" /><figcaption>**Fig. 26.** Hyperbolic media and nonlocal effects. (a) Isofrequency contour (i.e., dispersion relation at a fixed ω) found in a local, electric, and hyperbolic media type I (left) and II (right). ε⊥(k⊥) and ε‖(k‖) are the out-of-plane and in-plane effective permittivity (supported mode's wave-vector) components. (b) Examples of hyperbolic metamaterials (top) and metasurfaces (bottom) constructed by properly arranging subwavelength elements. (c) Isofrequency contour of the modes supported by an array of graphene strips that behaves as a hyperbolic surface [504]. Results are calculated using local EMT and rigorous full-wave numerical simulations considering local and nonlocal models of graphene. Left panel shows how granularity becomes dominant when the periodicity of the unit-cells is not much smaller than the wavelength of the supported modes. Right panel illustrates how the finite velocity of electrons traveling in graphene dominates the response when granularity is negligible. (d) Spontaneous emission rate of a z-oriented infinitesimal dipole versus its distance above the graphene-based hyperbolic surface [504]. Reproduced panels (c) and (d) with permission from Ref. [504] (Copyright © 2015 Optica Publishing Group).</figcaption></figure>
case, the interplay between these two nonlocal mechanisms will tailor the medium’s response. In general, EMT is just an approximation and thus local predictions regarding infinite local density of states, wave confinement, and spontaneous emission rate are incorrect. Figure 26(d) illustrates how the presence of nonlocality limits the maximum spontaneous emission rate of an emitter located near the graphene-based hyperbolic surface. It should also be mentioned that EMT might not be accurate to describe a certain type of artificial hyperbolic media, for instance in terms of predicting the number of supported modes and their properties or accounting for polarization-dependent inter-element coupling among the unit-cells that compose the structure.

## Challenges and opportunities

Accounting for nonlocal effects in hyperbolic media is not an easy task. To accurately consider the influence of the elements that compose the media, the most straightforward approach is to solve Maxwell’s equations using full-wave simulations that account for both geometrical details and the optical response of the materials involved. This powerful approach faces important complications, including *i*) inability to incorporate nonlocal material behavior in most scenarios; *ii*) lack of physical insight and analytical models that guide device design; and *iii*) requirement of large computational resources. For these reasons, dedicated formulations have been developed to model canonical nonlocal hyperbolic structures, including multilayered [407,507] and wire [508] bulk configurations and planar metasurfaces [501,509]. Even though these efficient approaches lead to closed-form expressions in certain scenarios, they are valid only under well-defined operation conditions and cannot be generalized for arbitrary structures.

Another challenge is to develop nonlocal optical theories to accurately describe the material response. These wavenumber-dependent models include phenomena such as quantum effects, electronic transitions, and induced charge diffusion kinetics, among others, and have been applied to metals [64], and particularized for materials such as graphene [510], and black-phosphorus [511]. In this context, the hydrodynamic Drude model within the Thomas-Fermi approximation [503,512] is a common nonlocal approach to characterize thin metals by employing a combination of longitudinal and transverse waves that account for the electron dynamics.

An exciting opportunity in the field of nonlocal hyperbolic media is the development of theoretical frameworks able to simultaneously account for the medium granularity and the intrinsic nonlocal response of its constitutive materials. The main challenge lies in incorporating complex wavenumber-dependent material responses within electromagnetic approaches. Even though some studies have already been put forward in this area [504,513] tools able to predict, model, and tailor hyperbolic wave and plasmonics propagation in realistic structures would drastically accelerate the application of this technology in practice. To date, hyperbolic devices have been mostly implemented by relegating nonlocality to a mere hindrance. This strategy not only misses the rich potential for tailored wave-control offered by nonlocality in hyperbolic media but leads to complex implementations with sub-optimal responses.

## Future developments to address challenges

Reconfigurable hyperbolic media have already been demonstrated in the local regime using phase-change and correlated materials, revealing exciting potential for programmable nanophotonic applications [518]. Tunable nonlocality has been explored in the context of current-driven plasmonics, implemented using materials such as graphene [519] or ultrathin metals [520], in which nonlocality appears because the supported waves effectively experience a different media when propagating along or against the drifting electrons. Exploiting this phenomena, drift-biased hyperbolic metasurfaces have been theoretically proposed [513]. In these structures, the interplay between nonlocality and hyperbolic dispersion breaks the symmetry of the supported modes in the *k*-space [Fig. 27(a)], leading to strong isolation and immunity to backscattering. The experimental demonstration of such platform would provide nonlocal, tunable, and nonreciprocal
<figure><img src="image_31.png" /><figcaption>**Fig. 27.** Applications of nonlocal hyperbolic metasurfaces. (a) Graphene-based drift-biased metasurfaces [513]. Drifting electrons break the symmetry of the supported modes, leading to nonreciprocal hyperbolic plasmons and immunity to backscattering. (b) Stacked moiré van der Waals crystals (top) and electric field distribution of a supported mode [514]. Reprinted (adapted) with permission from Ref. [514] (Copyright © 2020 American Chemical Society). (c) Moiré hyperbolic metasurfaces coupled through evanescent fields [515]. Spin-orbit interactions allow a rotating dipole to excite a directional beam. Reprinted (adapted) with permission from Ref. [515] (Copyright © 2020 American Chemical Society). (d) Thermal emission with hyperbolic devices. Reproduced with permission from Ref. [516] (Copyright © 2024 Nature Springer). (e) Nonlocal acoustic hyperbolic metasurfaces. Reproduced with permission from Ref. [517] (Copyright © 2024 Wiley). (f) In-plane optical computing with tunable nonlocal hyperbolic metasurfaces.</figcaption></figure>

responses useful to control, collimate, and route waves with a performance well beyond to the one found in common hyperbolic structures.

Recent developments in the field of twistronics and moiré physics have also enriched the possibilities offered by hyperbolic media. For instance, it is possible to manipulate the twist angle between stacked van der Waals crystal [Fig. 27(b)] to control their photonic and electronic properties and engineer hyperbolic responses [514]. At the mesoscopic scale, such response can be obtained by rotating hyperbolic metasurfaces that are located close to each other and evanescently coupled [Fig. 27(c)] [514]. Even though these platforms have been considered mostly in the local regime, nonlocality controls and bounds their fundamental response. Merging drift-current media with moiré physics will open new and unexpected opportunities for tunable, nonlocal, nonreciprocal, and hyperbolic wave control in which tailored nonlocal responses can be spatially and temporally manipulated to enable applications in sensing, communication systems,
and computing, among others. Additionally, it would also open the door to the generation of high *k*-waves [521] and enhanced nonlinear responses [442]—enabled by relaxed phase-matching conditions provided by nonlocality. These developments are inherently multidisciplinary and can be translated to other fields. For instance, thermally engineered devices [Fig. 27(d)] and nonlocal acoustic moiré hyperbolic metasurfaces [Fig. 27(e)] have recently been demonstrated [516,517,522]. Arguably, the next step in these areas is the tunable manipulation of nonlocality in both space and time, which would significantly expand the possibilities to radiate, receive, and process acoustic and thermal waves.

In a related context, nonlocal metasurfaces have become transformative in the field of optical computing, leveraging their ability to manipulate optical wavefronts in the momentum domain with unparalleled efficiency [523]. Nonlocal metasurfaces have been mostly designed to interact with light propagating in free-space and to provide tailored transfer functions in the momentum domain (i.e., different responses versus the angle of incidence of the incoming waves) that perform complex mathematical operations such as spatial differentiation, edge detection, or Fourier filtering. A paradigm shift would be to exploit high *k*-waves (evanescent spectrum) for optical computing. In such platform, illustrated in Fig. 27(f), input signals will be provided by near-field emitters such as quantum dots or by capturing light coming from the far-field using scatters and/or gratings. Nonlocal hyperbolic media provides unprecedented possibilities to tailor their in-plane transfer function versus the incident wavenumber. This paradigm offers new and rich possibilities for optical computing and multi-functional devices enabled by nonlocal processing of high *k*-waves that are not accessible in metasurfaces interacting with light propagating in free-space. Many technological challenges remain to be solved before this technology becomes a reality, mostly in terms of demonstrating tunable and nonlocal hyperbolic field control. Moving beyond, the integration of nonlocal hyperbolic media with optical neural networks and neuromorphic computing could pave the way for advanced and ultra-efficient signal processing capabilities over miniaturized devices.

### Concluding remarks

While nonlocality has usually imposed constraints on features such as wave propagation and the local density of states offered by hyperbolic media, it also provides an opportunity to develop innovative strategies for tailored wave control. By advancing theoretical frameworks, improving nanofabrication techniques, and leveraging emerging concepts like drift-current media and moiré physics, the interplay between nonlocality and hyperbolic dispersion can be exploited to unlock the full potential of hyperbolic media and construct nonreciprocal and programmable transfer functions in the momentum domain. These advancements will address existing challenges but will also pave the way for novel applications in a variety of fields, with emphasis on optical computation using high *k*-waves, sensing, acoustics, and thermal engineering. Exploiting nonlocality as an extra degree of freedom in hyperbolic designs marks a critical step toward realizing next-generation photonic and multidisciplinary technologies with unprecedented capabilities.

**Funding.** W. M. Keck Foundation.

## 23. Nonlocal effects in transdimensional plasmonics

IGOR V. BONDAREV, SVEND-AGE BIEHS, ALEXANDRA BOLTASSEVA & VLADIMIR M. SHALAEV

### Overview

Plasmonic transdimensional (TD) materials are atomically thin metal, semimetal or doped semiconductor films of precisely controlled countable number of monolayers [436]. Due to
current progress in nanofabrication techniques [524–527], such materials can be reproducibly grown and offer high tailorability of their electronic and optical properties not only by altering their chemical and/or electronic composition (stoichiometry, doping) or strain but also by merely varying their thickness. While ultrathin films have been extensively studied, the optical properties of TD materials remain somewhat underexplored. So far, the focus has largely been on either purely two-dimensional (2D) structures including metal-dielectric interfaces and novel 2D materials [528], or on conventional bulk materials, where the dimensionality and composition are the primary factors influencing the optoelectronic response. Recently it was proposed that ultrathin TD nanostructures provide a new regime—transdimensional, in between 3D and 2D, turning into 2D as the film thickness tends to zero [525–527]. In this regime, the strong vertical quantum confinement makes the linear electromagnetic (EM) response of the film nonlocal (spatially dispersive), and the degree of nonlocality can be controlled by the film thickness [529–531]. This makes plasmonic TD materials indispensable for studies of the nonlocal light-matter interactions at the nanoscale [530–537], where they exhibit extraordinary tailorability including the capabilities of active tuning of their EM response and thus enabling new and unique light-matter coupling phenomena [367]. In addition to strong dependencies on structural parameters such as the type of substrate/superstrate and strain [538], plasmonic properties of ultrathin TD structures show high sensitivity to external optical and electrical stimuli [539–541], which could facilitate the realization of dynamically tunable ultrathin plasmonic devices. The remarkable opportunities for tuning their highly confined in-plane plasma modes open access to novel quantum, nonlocal and nonlinear optical effects [367,533–541] in both the near-infrared and the visible ranges [539,542]. Over the past years, TD plasmonics has grown into a promising research direction aiming to greatly enhance plasmonic functionalities through unparalleled tailorability of confinement-induced nonlocal EM response effects [527,543].

Current status

The nonlocal optical properties of TD plasmonic films can be understood in terms of the confinement-induced nonlocal EM theory built using the Keldysh–Rytova (KR) electron interaction potential [529,530]. In optically dense thin films of metals, semimetals and doped semiconductors, the electrostatic Coulomb field produced by vertically confined remote charge carriers outside of their confinement region starts playing a perceptible role with the confinement size reduction [544,545]. The Coulomb interaction of such confined charges is typically stronger than that in a homogeneous medium with the same dielectric permittivity. This is due to the increased field contribution from outside dielectric environment with lower dielectric permittivity. In the case of TD plasmonic films in a typical “sandwich” geometry (Fig. 28, top panel), where the film of thickness *d* in region 2 is surrounded by semi-infinite dielectrics of constant permittivities $\varepsilon_1$ (top panel region 1) and $\varepsilon_3$ (bottom region 3), this turns the electron-electron Coulomb repulsion potential in the film into a much stronger *d*-dependent KR interaction potential [544,545]. This in turn leads to the *d*-dependent electron plasma oscillation frequency (see Fig. 28, top panel), resulting in the nonlocal Drude-like in-plane EM response of the TD film [529]:

$$ \varepsilon_2(\omega, k) = \varepsilon_b \left[ 1 - \frac{\omega_p^2(k)}{\omega(\omega + i\Gamma_D)} \right], \quad \omega_p(k) = \frac{\omega_p^{3D}}{\sqrt{1 + 1/(\varepsilon kd)}}, \quad \tilde{\varepsilon} = \frac{\varepsilon_b}{\varepsilon_1 + \varepsilon_3} \quad (63) $$

Here, $\varepsilon_b(>\varepsilon_1, \varepsilon_3)$ is the nearly constant permittivity due to bound electrons, $\Gamma_D$ is the damping constant, and $\omega_p(k)$ is the spatially dispersive (nonlocal) plasma frequency as a function of the in-plane momentum *k* and film thickness *d*; $\omega_p^{3D}$ is the bulk plasma frequency for the material the TD film is made of.

The nonlocal theoretical model for Drude-like EM response (aka the nonlocal KR model [534]) is verified experimentally in a variety of settings [525,527,535,546] (Fig. 28). The theory
<figure><img src="image_32.png" /><figcaption>**Fig. 28.** Recent progress in theoretical and experimental studies of the nonlocal EM effects in TD plasmonic systems. Top: Experimental demonstration of the theoretically predicted thickness dependence of the Drude plasma frequency in TD TiN films [Reprinted (adapted) with permission from Ref. [525] – Copyright © 2022 American Chemical Society]. Middle: (a) Thickness-dependent epsilon-near-zero modes of TiN film with variable thickness [Reproduced with permission from Ref. [530] – Copyright © 2022 American Physical Society]; (b) near-surface distance-/thickness-dependent spontaneous emission rate relative to vacuum [531] for a quantum dipole emitter coupled to the TD film plasma modes in (a) [Reproduced with permission from Ref. [531] – Copyright © 2023 Wiley]; (c) Casimir attraction force calculated using the nonlocal KR model for two free-standing metallic slabs of 10 (black), 20 (blue) and 200 nm (green) in thickness [Ref. [536] – Copyright © 2023 Royal Society of Chemistry]. Bottom left: Far-field heat fluxes predicted by local Drude (dashed lines) and nonlocal KR (solid lines) model for various TD metals as compared to black body (BB) radiation [534]. Bottom (a) and (b): Near-field heat conductance measured for ultrathin Pt films (inset) as compared to local Drude and nonlocal KR model predictions [535].</figcaption></figure>

accounts for the vertical electron confinement due to the presence of the substrate and superstrate of dielectric permittivities, which are smaller than that of the TD film. Importantly, the thickness of the film becomes a parameter that controls its generalized linear magneto-optical response [547]. The nonlocal KR model covers not only ultrathin films with thickness smaller than the half-wavelength of the incoming light but also conventional films with the thickness exceeding the optical wavelengths [530]. The nonlocal EM response of TD plasmonic systems has been shown to enable a variety of new effects, such as thickness-controlled plasma frequency redshift [525], low-temperature plasma frequency drop-off [546], plasma mode degeneracy lifting [530], and other quantum-optical [531,533], nonlocal magneto-optical effects [547], thermal and vacuum field fluctuation effects responsible for the near-/far-field radiative heat transfer [534,535] and Casimir interaction phenomena [536,537] (Fig. 28).

## Challenges and opportunities

The confinement-induced thickness-dependent nonlocal EM response [Eq. (63)] of the TD plasmonic systems provides the *k*-infinitesimal (first) order nonlocality type in the ultrathin
regime. It is dominant in the region of small in-plane momenta $k \le k_c$, where $k_c$ is the in-plane plasmon cut-off wave vector, as compared to the much smaller standard $k^2$-infinitesimal (second) order type nonlocality. The latter is due to the degenerated electron gas pressure in bulk metals described by hydrodynamical Drude plasma models [43]. The remarkable opportunity to tune the first-order confinement-induced nonlocality through thickness-dependent plasma frequency $\omega_p(k)$ in Eq. (63) is advantageous for applications of TD plasmonics in the fields of surface and quantum optics and with regard to the development of extra sensitive optical sensor technologies. To illustrate this, we refer to the TD-plasmonic-film-enhanced Goos-Hänchen (GH) and Fedorov-Imbert (FI) shifts of a light beam incident on a metallic surface [548] [see Fig. 29(top) (left)]. Both effects are known to depend on the incident beam polarization [548]. The GH and FI shifts are representative of TM/TE ($p/s$) linearly polarized and circularly (or elliptically) polarized incident waves, respectively. Below we discuss the GH shift in more detail [542]. Originating from the spatial dispersion of the reflection or transmission coefficients due to the finite transverse size of the beam (and so nonlocal), the GH effect was observed in a variety of systems including plasmonic metamaterials [549] and graphene [550]. In the reflection configuration (Fig. 29, bottom left), the lateral $\Delta_{GH}$ and angular $\Theta_{GH}$ GH shifts of TM polarized incoming waves are given by [542]

$$ \Delta_{GH} = n_1 \cos \theta_i \frac{\partial \phi_p}{\partial k}, \quad \Theta_{GH} = -\frac{\theta_0^2}{2} k_0 n_1 \frac{\cos \theta_i}{|R_p|} \frac{\partial |R_p|}{\partial k} \quad (64) $$

Here, $k_0 = \omega/c$, $n_1$ is the refractive index for the medium that the finite cross-section (Gaussian) light beam comes from; $\theta_0 = 2/(k_0 w_0 n_1)$ with $w_0$ representing the beam waist; and the p-wave reflection coefficient is written as $R_p = |R_p| \exp(i\phi_p)$ in the complex exponential form. It can be seen that the spatial dispersion makes $\Delta_{GH}$ sensitive to reflectivity phase jumps and $\Theta_{GH}$ to zero reflection itself so that large effects are highly likely for zero-reflection modes in the system. Phase jumps and singularities make the phase ill-defined, and the reflection coefficient's absolute value must be zero as required by the causality.

The GH shifts in Eq. (64) are normally thought of as being due to the spatial dispersion of the incoming light beam itself while the material EM response nonlocality is usually ignored. For bulk materials or macroscopically thick films this is indeed the case. However, it is not the case for the ultrathin TD plasmonic films where the strong vertical confinement makes the in-plane EM response nonlocal as described by Eq. (63) of the KR model. We use TD titanium nitride (TiN) as an example here, which is known for its exceptional plasmonic properties and high crystallinity down to thickness as small as 1 nm [525,551]. Then, there is an additional contribution to be taken into account in Eq. (64) that is proportional to

$$ \frac{\partial \varepsilon_{\text{TiN}}(\omega, k)}{\partial k} = - \frac{\varepsilon_b \omega_p^2(k)}{k(1 + \varepsilon kd)\omega(\omega + i\Gamma_D)} = [\varepsilon_{\text{TiN}}(\omega, k) - \varepsilon_b] \frac{\varepsilon d}{(1 + \varepsilon kd)^2}, \quad (65) $$

which is not only nonzero at finite $d$ but can also be both positive and negative, depending on the light frequency. It disappears when $d$ goes to infinity—as it should to make the EM response of thick films local in accord with the standard Drude model.

It was recently shown theoretically [542] that the confinement-induced nonlocality of the TD films leads to topologically protected singularities of the nonlocal reflection coefficient. Such singularities are shown to result in giant lateral and angular GH shifts in the millimeter and milliradian ranges, respectively, which exceed greatly those reported for light beams of finite transverse extent with no material-induced nonlocality [548–551]. The singularities appear in the TD film systems with broken in-plane reflection symmetry (substrate and superstrate of different dielectric permittivities), where due to the strong vertical confinement the eigenmode degeneracy is lifted creating the points of the topological darkness in the visible range that do not exist in the
<figure><img src="image_33.png" /><figcaption>**Fig. 29.** Top left: Schematic illustrations to demonstrate the Goos-Hänchen (top) and Fedorov-Imbert (bottom) effects for the reflection/refraction of the finite cross-section light beam as it is shifted/deflected in the plane of incidence and out of the plane of incidence, respectively (reproduced with permission from Ref. [548] – Copyright © 2013 IOP Publishing). Bottom left: Sketch of the lateral and angular GH shifts of the light beam incident on a TiN plasmonic slab of broken in-plane reflection symmetry. Right panel – top, (a) and (b): Reflection phase $\phi_p/\pi$ and absolute values of the lateral and angular GH shifts, respectively, calculated for a 40 nm thick TiN slab with $\varepsilon_1 = 1$ and $\varepsilon_3 = 3$ (MgO substrate) neglecting dissipation. Lines show the eigenmodes available in the system. Circles indicate the topological singularity points at the eigenmode intersections. See Ref. [542] for details.</figcaption></figure>

usual (local) Drude materials. The right panel in Fig. 29 shows the topological darkness points (green circles) that appear at the intersections of the non-degenerate EM eigenmodes of the TD film system. The full analysis and the classification of the singularity points can be found in Ref. [542]. The remarkable opportunity to bring the GH effect to the visible range comes from the fact that the nonlocal plasma frequency in Eq. (63) can be red-shifted not only by thickness but also by in-plane momentum reduction through the incident angle change. Lateral and angular GH shifts as large as ~0.4 mm and ~40 mrad, respectively, are predicted theoretically for typical 40 nm thick TD plasmonic films, using for example the conventional helium-neon (He-Ne) laser light (Fig. 29, right panel).

## Future developments to address challenges

While the experimental progress in TD materials for nanophotonics has earlier been impeded by challenges in producing atomically thin films of noble metals, TD films of emerging plasmonic materials such as transition metal nitrides (TiN, ZrN, HfN, etc. [551]) can be grown as epitaxial-quality films with thicknesses down to 1–2 nm (5–10 atomic layers) [525,527]. A variety of
transition metal nitrides, their ability to grow as high-quality, ultrathin epitaxial films and the sensitivity of their optical and electronic properties to the material/structural/geometrical parameters provide a rich playground for the realization of the confinement-induced nonlocal effects. This includes using transition metal nitrides in TD plasmonics to control strong electron correlations by precise variation of the film thickness to study fundamental solid-state physics phenomena such 2D electron Wigner crystallization [543] and metal-insulator transitions [527]. Another important avenue to tailor the nonlocal optical properties of plasmonic TD films includes strain engineering. Below a critical thickness, an epitaxial thin film is expected to retain the strain induced by the substrate. This can be achieved experimentally by growing strained ultrathin films on lattice-mismatched substrates. It has been theoretically demonstrated that by varying the in-plane lattice parameter of an ultrathin film its nonlocal optical response can be tuned [538]. In conjunction with the thickness dependence of the nonlocal EM response discussed above, the strain engineering offers an additional way to tailor the optical response of TD plasmonic materials.

### Concluding remarks

Plasmonic TD materials hold great potential for enabling a new quantum material platform that utilizes strong confinement and nonlocal effects and offers new opportunities for quantum optics, quantum computing and sensing applications. A crucial step for such development is the realization of high quality, ultrathin metallic films of precisely controlled thickness with reduced surface roughness that exhibit new material functionalities and unique light-matter interactions. The quantum effects that arise in metallic TD materials along with their strong tunability may pave the way to new optical phenomena and novel, atomically thin, dynamically tunable nanophotonic devices.

**Funding.** Army Research Office (W911NF2310206).

## 24. Light management using structural nonlocality in nanorod metamaterials

ALEXEY V. KRASAVIN & ANATOLY V. ZAYATS

<figure><img src="image_34.png" /></figure>

### Overview

Engineering the shape of elemental metamaterial components, meta-atoms, as well as their arrangement on a subwavelength scale results in artificial materials with pre-designed optical responses both in linear and nonlinear regimes. The optical properties of such discrete media can be understood, in principle, through conventional local effective medium models. Comparing these models with full-vectorial numerical simulations, however, shows intricate differences, which are relatively small in general, but become crucial for certain metamaterial parameters, for which a local effective medium theory (EMT) breaks down, even in the range of its applicability. This behavior is related to the so-called nonlocal, spatial dispersion effects which are not taken into account in conventional, local EMT models (numerical simulations of a discrete structure of a metamaterial, however, includes these effects automatically).

Nonlocality is a peculiar property of optical materials, when the optical response of the medium at a given point depends not only on the driving electromagnetic field at this location, but also on the values of the field in other points of the surrounding region. This property is not so common, and it is clear why. For natural condensed non-metallic materials, liquids or solids, the optical response of an atom (or a molecule) is influenced by the optical responses of its counterparts located only in its close vicinity, particularly at the distances of the order of a lattice constant (or a molecular size), which is much smaller than the wavelength of light. Consequently, there is no variation of the field across this region and it is a valid assumption that the optical response of
this atom is defined by the value of the field at its location. This can be applied to any atom of the material, and therefore the optical response of the material is local. The situation changes when large quasiparticles are present in the material, such as excitons in semiconductors, resulting in spatial dispersion which strongly influences the optical response at low temperatures and leads to the appearance of additional waves [20]. The nonlocality of metals important for the nanostructures with the size below few nanometers, associated with electron spill-out effects and wavevector-dependent dielectric permittivity, is discussed in Parts I and II of this Roadmap, including Secs. 3, 4, 5, 7, and 8. A different situation can occur in metamaterials, in which both the size of the meta-atoms and their spacing, although being smaller than the wavelength in order to satisfy effective medium requirements, are noticeable fractions of it. This results in *i*) the field variation across the meta-atom, which can lead to the excitation of higher-order multipoles, thus affecting its optical response, and *ii*) the field at the positions of other meta-atoms affecting the optical response of the considered meta-atom starts to depend on their positions [552]. Thus, the optical response of a metamaterial at a given location can be noticeably affected by the field at other points, naturally producing a physical mechanism of the metamaterial structural nonlocality. In the reciprocal-space representation, this means that the optical properties of the metamaterial depend on the wave vector, i.e. the material is spatially dispersive.

With subwavelength structuring and the consequent absence of the diffractive effects, a metamaterial can be considered as a uniform medium, which optical properties are defined by the metamaterial design and can be described using EMTs. In the first approximation the metamaterial optical response can be considered to be local, but for its precise description development of nonlocal EMTs is needed. This is especially important for optical effects depending on the local fields inside a metamaterial, such as photoluminescence or nonlinearity, as nonlocality can create highly nonuniform local fields due to the introduction of an additional wave. The development of nonlocal EMTs is not an easy task. Generally, it requires deriving a microscopic optical response of the metamaterial by solving the electromagnetic problem for the exact nanostructure designs. Within two common approaches, the effective optical parameters, depending on the wave vector, can be found by averaging the derived microscopic fields over the metamaterial unit cell, or setting trial (or Taylor-expanded in the *k*-space) expressions for the effective permittivity and/or permeability and matching the dispersion of the corresponding optical modes to that found in the microscopic description [553].

## Current status

There are numerous experimental demonstrations of optical phenomena, in which the nonlocal response of a metamaterial plays a crucial role. As a prominent example, a metamaterial formed by an array of aligned plasmonic nanorods represents a uniaxial medium possessing record-high optical anisotropy, hyperbolic dispersion related to different signs of dielectric permittivity along the optical axes, and an epsilon-near-zero (ENZ) response related to $Re(ε_∥) ≈ 0$ (the permittivity component along the optical axis parallel to the nanorods) [554] [Fig. 30(a,b)]. The presence of structural nonlocality, which can be derived using various EMT approaches and approximations [553,555,557] fundamentally modifies these characteristics. In the local approximation, optical extinction of the nanorod metamaterial layer features two peaks [located for the typical metamaterial parameters at around 500 and 650 nm in Fig. 30(c)]. The short-wavelength peak, existing for both TE and TM polarized illumination, is related to the excitation of the transverse plasmonic resonance of the nanorods, slightly shifted due to inter-rod coupling. Its long-wavelength counterpart, existing only for TM polarized illumination, corresponds to a metamaterial opacity region near the ENZ spectral point for $ε_∥$. The magnitude of this peak increases with the angle of incidence, but in the local approximation its spectral position is angle independent. However, for low ohmic losses of the nanorod plasmonic metal, the behavior of the ENZ peak reveals an intricate angle dependence [555] [Fig. 30(c)]. Instead of a single peak, two ENZ extinction peaks showing a spectral dynamics characteristic to an anti-crossing of optical
modes are observed. These peaks are related to two TM-polarized modes, which result from the hybridization of the ordinary TM mode of light transmitted through the metamaterial in the local approximation with an additional longitudinal mode supported solely due to the metamaterial structural nonlocality [555]. The latter, due to its dispersion characteristics and field distribution is optically accessible only in the ENZ region. The presence of two TM-polarized waves in the nonlocal case can be clearly seen in the plots of isofrequency surfaces [Fig. 30(d)] [556]. The nature of this nonlocality can be traced to the coupled cylindrical surface plasmon polaritons supported by the nanorods in the optical regime [557]. The interference of the hybridized TM modes and, therefore, the extinction of the metamaterial in the ENZ region have been shown to be highly dependent on the optical properties of the constituting materials. Active control of losses in plasmonic nanorods, e.g., through the excitation of hot electrons upon light absorption or temperature effects [558] results in effective switching the nonlocal response on/off, leading to ultrafast all-optical switching [559] and polarization control [560] in the nonlocal regime.

<figure><img src="image_35.png" /><figcaption>**Fig. 30.** (a) Schematics of a nanorod metamaterial with indicated characteristic geometrical parameters and materials. (b) Spectra of the components of the metamaterial permittivity tensor (*d* = 25 nm, *s* = 60 nm) calculated using a local EMT [554]. (c) Experimental extinction spectra of a 300-nm-thick nanorod metamaterial layer [Au/Al₂O₃, geometrical parameters the same as in (b)] under TM-polarized illumination for various angles of incidence. Reproduced with permission from Ref. [555] (Copyright © 2009 American Chemical Society). (d) Isofrequency surfaces of two hybridized TM-polarized waves plotted in the elliptical dispersion regime (λ = 550 nm, *d* = 50 nm, *s* = 100 nm). Reproduced with permission from Ref. [556] (Copyright © 2017 Nature Springer).</figcaption></figure>

The hybridization of the main and additional modes in a nonlocal regime in the ENZ spectral region creates very large effective frequency dispersion for optical pulses transmitted through a metamaterial layer, which has been used to manage their propagation and temporal characteristics [561]. Sub- and super-luminal, as well as backward propagation of femtosecond pulses have been demonstrated, with the switching of the group dispersion achieved in the same metamaterial by just changing the angle of illumination [Fig. 31(a)]. Combining this with nonlinear functionalities discussed above, one can potentially realize all-optical control of the pulse characteristics.
<figure><img src="image_36.png" /><figcaption>**Fig. 31.** (a) Transient extinction spectra of a plasmonic nanorod metamaterial (Au nanorods in alumina template, *d* = 20 nm, *s* = 70 nm, *h* = 400 nm) at a 20° incident angle for various pump fluences at *λ*<sub>p</sub> = 465 nm and *Δ*(*t*) = 130 fs. Reproduced with permission from Ref. [559] (Copyright © 2011 Nature Springer). (b) Group delay experienced by an optical pulse while propagating through a 525-nm-thick metamaterial layer (*d* = 25 nm, *s* = 60 nm). Reproduced with permission from Ref. [561] (Copyright © 2023 Wiley). (c) Enhancement of the spontaneous emission rate of quantum emitters inside a nanorod metamaterial (*d* = 50 nm, *s* = 100 nm): spectral dependences for experimentally measured, numerically simulated (with the exact microscopic geometry) for a 250-nm-thick metamaterial layer (yellow band shows the range of the values for various locations inside the metamaterial layer), and calculated using local and nonlocal EMTs for the bulk metamaterial are shown. Reproduced with permission from Ref. [556] (Copyright © 2017 Nature Springer). (d) Purcell factor of a dipole with averaged orientation located 10 nm above a multilayer metamaterial (layer thicknesses: *a*<sub>1</sub> = 12 nm Ag, *a*<sub>2</sub> = 83 nm SiO<sub>2</sub>) surface calculated using a nonlocal EMT (black solid line) and finite element numerical simulations (color lines, for various metamaterial unit cell sizes *a* = *a*<sub>1</sub> + *a*<sub>2</sub>, keeping the same ratio between Ag and SiO<sub>2</sub> thicknesses *a*<sub>1</sub>/(*a*<sub>2</sub>), together with experimental dependences measured using three types of quantum dots (scatter plots). Reprinted (adapted) with permission from Ref. [562] (Copyright © 2017 American Chemical Society).</figcaption></figure>

Crucially, the structural nonlocality has a strong effect on the emission of quantum emitters placed inside the metamaterial. The nonlocal description has been shown to provide a more adequate prediction for the emitter spontaneous decay rate, contrary to its almost singular enhancement given by the local EMT [556] [Fig. 31(b)]. Experiments, together with the microscopic numerical simulations have shown that considering the exact metamaterial structure (as opposed to implementation of EMT theories) and the exact position of the emitters with respect to the meta-atoms is important for correct understanding of the emission process.

Structural nonlocality has been also studied in multilayer metamaterials, produced by alternating metallic and dielectric layers with nanoscale thicknesses, which share major optical characteristics
with the nanorod metamaterial, such as extremely high anisotropy, an ENZ behavior and hyperbolic dispersion. Nonlocality in these metamaterials is underlined by strong variation of the fields across individual layers (i.e. across meta-atoms) [563]. It affects the metamaterial transmission/reflection properties, the spectral position of the ENZ region and the emission of quantum emitters placed in the vicinity or inside the metamaterial. Particularly, Fig. 31(d) shows that the Purcell factor predicted by a nonlocal EMT matches the numerically simulated counterpart in the limit of very small unit cell sizes, while the numerical predictions agree well with the experimentally measured Purcell factor values for realistic geometrical parameters [562].

## Challenges and opportunities

From the theoretical perspective, various approaches have been developed to derive nonlocal EMTs for metamaterials, starting from that assuming ideally conducting wires (relevant for microwave frequencies) and wires having large negative real part of the permittivity (appropriate for wavelength in mid and far-IR) to the models with no specific restriction made on the epsilon of the rods (apart from the implicit assumption that they are metallic [555,557]). Commonly, they return the same general behavior, featuring local effective $\epsilon_{\perp}$ and nonlocal $\epsilon_{\parallel}$ components, depending on the wave vector along the optical axis $k_{\parallel}$. At the same time, frequently models derived under similar assumptions result in different expressions for the nonlocal effective permittivity, which motivates for comparison and generalization of the nonlocal EMT models [553,564]. A related question here is determining the additional boundary conditions required for the additional wave, characteristic of the spatial dispersion. Substantial benefit would be the derivation of new or application of existing models to define the range of metamaterial geometrical parameters which maximize the influence of structural nonlocality and/or maximize the magnitude of nonlocality-assisted phenomena, which would provide a valuable guidance for experimental studies. Of course, all the above research directions are very challenging, as any broadly applicable nonlocal EMT should ideally satisfy several important conditions: i) provide a correct dispersion relation for the metamaterial modes, ii) give correct boundary conditions (note that in the case of nonlocality additional boundary conditions are needed to take into account the existence of the additional optical mode [557]), and iii) be valid beyond the case of plane wave excitation [565]. At the same time, this means that there is a large scope for the future theoretical research.

Experimentally, since the nonlocal response depends on losses of metallic components in plasmonic metamaterials, the realization of low loss material platforms (including 2D materials, e.g. graphene) together with a broader range of geometrical parameters will open opportunities for the realization and exploitation of stronger spatial dispersion contributions in a desired spectral range for controlling molecular emission, with potential applications in bio-imaging and quantum technologies, and enhancement of optical nonlinearities, important for laser technologies and optical information processing.

## Future developments to address challenges

Nonlocal metamaterials have been shown to provide multiple functionalities, bio and chemical sensing being among the most prominent. Usually, the sensitivity of resonant excitation of guided modes supported by a metamaterial layer to the changes of the dielectric environment is used for this purpose [554]. At the same time, nonlocal effects can also make the extinction of the metamaterial in the ENZ region highly sensitive to such changes, through the same physical mechanism as was illustrated above for the realization of optical switching. To demonstrate this, metamaterials with maximized structural nonlocality, achieved by engineering the design and reducing the optical losses in metal, may be needed.

Apart from incoherent Kerr-type nonlinearity which exploits optically-induced heating of the electron gas discussed above, nanorod metamaterials possess its coherent counterpart based on
a complex electron dynamic coherent with the excitation wave [134]. Particularly, it has been shown that engineering of the modal structure of a nanorod metamaterial layer can enhance second harmonic generation [566], while alternative approaches based on back-propagating waves have also been proposed [567]. In this respect, developing nonlinear EMTs including phenomena related to structural nonlocality, as proposed e.g. in [568], would provide a valuable guidance to the experimental research. Possessing extreme anisotropy and nonlocality which affect only certain components of the effective permittivity tensor, nanorod and multilayer metamaterials present a particular interest for engineering light-matter interaction with complex vectorial beams. In this case, the interaction of the nonlocal $ε_‖$ component in the nanorod metamaterials with co-directed longitudinal fields, present in radially polarized cylindrical beams of various orders, can be used for controlling optical spin-orbit coupling and related effects. Another research direction could be the study of chiroptical effects, which have been shown to be inherent to the artificial media with effective structural spatial dispersion [569]. Going beyond uniaxial metamaterials to truly 3D metamaterial architectures may provide additional flexibility for exploitation of spatial dispersion effects.

### Concluding remarks

Spatial dispersion in metamaterials arising from sub-wavelength structuring results in nonlocal effects with the electromagnetic response depending not only on the field at a given point but also on its variations across several unit cells. It is essential for accurately describing wave propagation in complex media, including ultrashort pulses, and in particular the phenomena depending on local fields inside the metamaterial. It has a much stronger impact on optical properties of metamaterials than in semiconductors where it is related to the excitonic excitation and observable at low temperatures. Recognizing additional flexibility provided by the spatial dispersion in the design of the metamaterial optical response is important in advancing nanophotonic technologies such as high-resolution imaging, bio- and gas-sensing, ultrashort pulse control, nonlinear optics, optical information processing and quantum technologies.

**Funding.** European Research Council (Advanced Grant 789340–ICOMM); Engineering and Physical Sciences Research Council (EPSRC project EP/Y015673/1).

## 25. Nonlocal wave phenomena in spatio-temporal metamaterials

ANDREA ALÙ ID

### Current status

The last few years have witnessed impressive developments in the context of metamaterials and metasurfaces, with remarkable progress in controlling electromagnetic wavefronts, and in the application of these phenomena within a wide range of technologically relevant contexts. One important breakthrough, detailed in several sections of this Roadmap, has been the discovery that nonlocal wave phenomena are not necessarily just a nuisance, a correction over well-understood phenomena in light-matter interactions, or a modeling headache, but rather they can be controlled, enhanced and engineered in metamaterial platforms, offering new degrees of freedom for enhanced wave control both in real and in reciprocal space. In this context, the emerging field of nonlocal metasurfaces [570] has been flourishing, with remarkable demonstrations in three relevant directions: i) spatial shaping of lattice resonances for enhanced metasurface control [571], ii) analog-based image processing and optical computing at the nanoscale [572,573], and iii) enhanced coherence control over light generation [522].

Beyond nonlocality, in a parallel research direction the dimension of time has recently offered new opportunities to enhance the degrees of freedom available for wave control in metamaterials. Beyond the three spatial dimensions, by engineering the arrow of time through tailored temporal