# Nonlocality in photonic materials and metamaterials: roadmap

Francesco Monticone, N. Asger Mortensen, Antonio Fernández-Domínguez, Yu Luo, Xuezhi Zheng, Christos Tserkezis, Jacob Khurgin, Tigran Shahbazyan, André Chaves, Nuno Peres, et al.

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

► To cite this version:

Francesco Monticone, N. Asger Mortensen, Antonio Fernández-Domínguez, Yu Luo, Xuezhi Zheng, et al.. Nonlocality in photonic materials and metamaterials: roadmap. Optical Materials Express, 2025, 15 (7), pp.1544. 10.1364/OME.559374. hal-05170351

HAL Id: hal-05170351

https://hal.science/hal-05170351v1

Submitted on 19 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Nonlocality in photonic materials and metamaterials: roadmap

FRANCESCO MONTICONE,<fnref n="1" />,*, N. ASGER MORTENSEN,<fnref n="2" />,<fnref n="3" />,‡

ANTONIO I. FERNÁNDEZ-DOMÍNGUEZ,<fnref n="4" />,<fnref n="5" /> YU LUO,<fnref n="6" />,<fnref n="7" />

XUEZHI ZHENG,<fnref n="8" /> CHRISTOS TSERKEZIS,<fnref n="2" /> JACOB B. KHURGIN,⁹

TIGRAN V. SHAHBAZYAN,¹⁰ ANDRÉ J. CHAVES,<fnref n="2" />,¹¹

NUNO M. R. PERES,<fnref n="2" />,¹²,¹³ GINO WEGNER,¹⁴

KURT BUSCH,¹⁴,¹⁵ HUATIAN HU,¹⁶ FABIO DELLA SALA,¹⁶,¹⁷

PU ZHANG,¹⁸ CRISTIAN CIRACÌ,¹⁶ JAVIER AIZPURUA,¹⁹,²⁰,²¹

ANTTON BABAZE,²² ANDREI G. BORISOV,¹⁹,²³

XUE-WEN CHEN,¹⁸ THOMAS CHRISTENSEN,²⁴ WEI YAN,²⁵,²⁶

YI YANG,²⁷ ULRICH HOHENESTER,²⁸ LORENZ HUBER,²⁸

MARTIJN WUBS,²⁴,²⁹ SIMONE DE LIBERATO,³⁰,³¹,³²

P. A. D. GONÇALVES,<fnref n="2" />,³³ F. JAVIER GARCÍA DE ABAJO,³³,³⁴

ORTWIN HESS,³⁵ ILLYA TARASENKO,³⁵ JOEL D. COX,<fnref n="2" />,<fnref n="3" />

LINE JELVER,<fnref n="2" /> EDUARDO J. C. DIAS,<fnref n="2" />

MIGUEL SÁNCHEZ SÁNCHEZ,³⁶ DIONISIOS MARGETIS,³⁷

GUILLERMO GÓMEZ-SANTOS,³⁸ IGOR M. VASILEVSKIY,³⁶

TOBIAS STAUBER,³⁶ SERGEI TRETYAKOV,³⁹

CONSTANTIN SIMOVSKI,³⁹ SAMANEH PAKNIYAT,⁴⁰

J. SEBASTIÁN GÓMEZ-DÍAZ,⁴⁰ IGOR V. BONDAREV,⁴¹

SVEND-AGE BIEHS,⁴² ALEXANDRA BOLTASSEVA,⁴³,⁴⁴

VLADIMIR M. SHALAEV,⁴³ ALEXEY V. KRASAVIN,⁴⁵

ANATOLY V. ZAYATS,⁴⁵ ANDREA ALÙ,⁴⁶,⁴⁷

JUNG-HWAN SONG,⁴⁸ MARK L. BRONGERSMA,⁴⁹

URIEL LEVY,⁵⁰ OLIVIA Y. LONG,⁵¹ CHENG GUO,⁵²

SHANHUI FAN,⁵¹,⁵² SERGEY I. BOZHEVOLNYI,⁵³

ADAM OVERVIG,⁵⁴ FILIPA R. PRUDÊNCIO,⁵⁵,⁵⁶

MÁRIO G. SILVEIRINHA,⁵⁵ S. ALI HASSANI GANGARAJ,⁵⁷

CHRISTOS ARGYROPOULOS,⁵⁸ PALOMA A. HUIDOBRO,<fnref n="4" />,<fnref n="5" />

EMANUELE GALIFFI,⁴⁶ FAN YANG,⁵⁹ JOHN B. PENDRY,⁶⁰

AND DAVID A. B. MILLER⁵²

<footnote n="1">¹School of Electrical and Computer Engineering, Cornell University, Ithaca, New York 14853, USA</footnote>

<footnote n="2">²POLIMA—Center for Polariton-driven Light-Matter Interactions, University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark</footnote>

<footnote n="3">³D-IAS—Danish Institute for Advanced Study, University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark</footnote>

<footnote n="4">⁴Departamento de Física Teórica de la Materia Condensada, Universidad Autónoma de Madrid, E-28049 Madrid, Spain</footnote>

<footnote n="5">⁵Condensed Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid, E-28049 Madrid, Spain</footnote>

<footnote n="6">⁶School of Electrical and Electronic Engineering, Nanyang Technological University, Singapore 639798, Singapore</footnote>

<footnote n="7">⁷Key Laboratory of Radar Imaging and Microwave Photonics, Ministry of Education, College of Electronic and Information Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 211106, China</footnote>

<footnote n="8">⁸WaveCoRE Division, Department of Electrical Engineering, KU Leuven, Leuven B-3001, Belgium</footnote>
⁹Department of Electrical and Computer Engineering, Johns Hopkins University, Baltimore, Maryland 21218, USA

¹⁰Department of Physics, Jackson State University, Jackson, Mississippi 39217, USA

¹¹Department of Physics, Aeronautics Institute of Technology, 12228-900 So Jos dos Campos, SP, Brazil

¹²Centro de Física (CF-UM-UP) and Departamento de Física, Universidade do Minho, P-4710-057 Braga, Portugal

¹³International Iberian Nanotechnology Laboratory (INL), Av Mestre José Veiga, 4715-330 Braga, Portugal

¹⁴Max-Born-Institut, 12489 Berlin, Germany

¹⁵Humboldt-Universität zu Berlin, Institut für Physik, AG Theoretische Optik and Photonik, 12489 Berlin, Germany

¹⁶Center for Biomolecular Nanotechnologies, Istituto Italiano di Tecnologia, Via Barsanti 14, 73010 Arnesano, Italy

¹⁷Institute for Microelectronics and Microsystems (CNR-IMM), Via Monteroni, Campus Unisalento, 73100 Lecce, Italy

¹⁸School of Physics, Huazhong University of Science and Technology, Luoyu Road 1037, Wuhan 430074, China

¹⁹Donostia International Physics Center DIPC, 20018 Donostia-San Sebastián, Spain

²⁰Ikerbasque, Basque Foundation for Science, 48009 Bilbao, Spain

²¹Department of Electricity and Electronics, University of the Basque Country, 48940 Leioa, Spain

²²Department of Applied Physics, University of the Basque Country, 20018 Donostia, Spain

²³Université Paris-Saclay, CNRS, Institut des Sciences Moléculaires d'Orsay, Orsay, France

²⁴Department of Electrical and Photonics Engineering, Technical University of Denmark, Kongens Lyngby, Denmark

²⁵Key Laboratory of 3D Micro/Nano Fabrication and Characterization of Zhejiang Province, School of Engineering, Westlake University, Hangzhou, Zhejiang Province, China

²⁶Institute of Advanced Technology, Westlake Institute for Advanced Study, Hangzhou, Zhejiang Province, China

²⁷Department of Physics and HK Institute of Quantum Science and Technology, The University of Hong Kong, Pokfulam, Hong Kong SAR, China

²⁸Institute of Physics, University of Graz, Universitätsplatz 5, 8010 Graz, Austria

²⁹NanoPhoton – Center for Nanophotonics, Technical University of Denmark, 2800 Kongens Lyngby, Denmark

³⁰School of Physics and Astronomy, University of Southampton, Southampton SO17 1BJ, United Kingdom

³¹Sensorium Technological Laboratories, Nashville, Tennessee 37210, USA

³²Istituto di Fotonica e Nanotecnologie, Consiglio Nazionale delle Ricerche (CNR), 20133 Milano, Italy

³³ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain

³⁴ICREA-Institució Catalana de Recerca i Estudis Avançats, Passeig Lluís Companys 23, 08010 Barcelona, Spain

³⁵School of Physics and CRANN Institute, Trinity College Dublin, Dublin 2, Ireland

³⁶Instituto de Ciencia de Materiales de Madrid, CSIC, E-28049 Madrid, Spain

³⁷Department of Mathematics, and Institute for Physical Science and Technology, University of Maryland, College Park, Maryland 20742, USA

³⁸Departamento de Física de la Materia Condensada, Instituto Nicolás Cabrera and Condensed Matter Physics Center (IFIMAC), Universidad Autónoma de Madrid, E-28049 Madrid, Spain

³⁹Aalto University, PO Box 15500, FI-00076 Aalto, Finland

⁴⁰Electrical and Computer Engineering Department, University of California, Davis, Davis, California 95616, USA

⁴¹Department of Mathematics & Physics, North Carolina Central University, Durham, North Carolina 27707, USA

⁴²Institut für Physik, Carl von Ossietzky Universität, 26111, Oldenburg, Germany

⁴³Elmore Family School of Electrical and Computer Engineering, Purdue Quantum Science and Engineering Institute, and Birck Nanotechnology Center, West Lafayette, Indiana 47907, USA

⁴⁴School of Materials Engineering, Purdue University, West Lafayette, Indiana 47907, USA
<sup>45</sup>Department of Physics and London Centre for Nanotechnology, King's College London, London, WS2R 2LS, United Kingdom

<sup>46</sup>Photonics Initiative, Advanced Science Research Center, City University of New York, New York, New York 10031, USA

<sup>47</sup>Physics Program, Graduate Center of the City University of New York, New York, New York 10016, USA

<sup>48</sup>Department of Electrical and Computer Engineering, National University of Singapore, Singapore 117583, Singapore

<sup>49</sup>Geballe Laboratory for Advanced Materials, Stanford University, Stanford, California 94305, USA

<sup>50</sup>Department of Applied Physics, The Faculty of Science, The Hebrew University of Jerusalem, Jerusalem, Israel

<sup>51</sup>Department of Applied Physics, Stanford University, Stanford, California 94305, USA

<sup>52</sup>Ginzton Laboratory and Department of Electrical Engineering, Stanford University, Stanford, California 94305, USA

<sup>53</sup>Center for Nano Optics, University of Southern Denmark, Campusvej 55, DK-5230 Odense M, Denmark

<sup>54</sup>Department of Physics, Stevens Institute of Technology, Hoboken, New Jersey 07030, USA

<sup>55</sup>University of Lisbon—Instituto Superior Técnico and Instituto de Telecomunicações, Avenida Rovisco Pais 1, 1049-001 Lisbon, Portugal

<sup>56</sup>Instituto Universitário de Lisboa (ISCTE-IUL), Avenida das Forças Armadas 376, 1600-077 Lisbon, Portugal

<sup>57</sup>Optical Physics Division, Corning Research and Development, Corning, New York 14831, USA

<sup>58</sup>Department of Electrical Engineering, The Pennsylvania State University, University Park, Pennsylvania 16802, USA

<sup>59</sup>College of Physics and Key Laboratory of High Energy Density Physics and Technology of the Ministry of Education, Sichuan University, Chengdu, Sichuan 610065, China

<sup>60</sup>The Blackett Laboratory, Department of Physics, Imperial College London, London SW7 2AZ, United Kingdom

‡asger@mailaps.org

*francesco.monticone@cornell.edu

**Abstract:** Photonic technologies continue to drive the quest for new optical materials with unprecedented responses. A major frontier in this field is the exploration of nonlocal (spatially dispersive) materials, going beyond the local, wavevector-independent assumption traditionally adopted in optical material modeling. The growing interest in plasmonic, polaritonic, and quantum materials has revealed naturally occurring nonlocalities, emphasizing the need for more accurate models to predict and design their optical responses. This has major implications also for topological, nonreciprocal, and time-varying systems based on these material platforms. Beyond natural materials, artificially structured materials—metamaterials and metasurfaces—can provide even stronger and engineered nonlocal effects, emerging from long-range interactions or multipolar effects. This is a rapidly expanding area in the field of photonic metamaterials, with open frontiers yet to be explored. In metasurfaces, in particular, nonlocality engineering has emerged as a powerful tool for designing strongly wavevector-dependent responses, enabling enhanced wavefront control, spatial compression, multifunctional devices, and wave-based computing. Furthermore, nonlocality and related concepts play a critical role in defining the ultimate limits of what is possible in optics, photonics, and wave physics. This Roadmap aims to survey the most exciting developments in nonlocal photonic materials and metamaterials, highlight new opportunities and open challenges, and chart new pathways that will drive this emerging field forward—toward new scientific discoveries and technological advancements.
Contents

<table><tbody><tr><td>1. Introduction (Monticone & Mortensen)</td><td>1548</td></tr><tr><td><b>Part I—Fundamental frameworks</b></td><td>1550</td></tr><tr><td>2. Linear-response formalism for nonlocal electrodynamics (Mortensen & Monticone)</td><td>1550</td></tr><tr><td>3. Nonlocal hydrodynamics (Fernández-Domínguez et al.)</td><td>1554</td></tr><tr><td>4. Imaginary part of nonlocal permittivity: how it localizes absorption and imposes limits on field confinement and enhancement in plasmonics (Khurgin)</td><td>1562</td></tr><tr><td>5. Landau damping of surface plasmons in metal nanoparticles: the RPA approach (Shahbazyan)</td><td>1566</td></tr><tr><td>6. Madelung hydrodynamics, nonlocal plasmonics, and nonlinear optics (Chaves et al.)</td><td>1571</td></tr><tr><td><b>Part II—Nonlocal effects in free-electron metals and plasmonic nanostructures</b></td><td>1576</td></tr><tr><td>7. Continuum framework for plasmonic systems (Wegner & Busch)</td><td>1576</td></tr><tr><td>8. Quantum hydrodynamic model (Hu et al.)</td><td>1580</td></tr><tr><td>9. Jellium model within time-dependent density-functional theory to address quantum effects in nanoplasmonics (Aizpurua et al.)</td><td>1585</td></tr><tr><td>10. Atomic scale nonlocal plasmonics (Zhang & Chen)</td><td>1589</td></tr><tr><td>11. Nonlocal limitations to light-matter interactions (Tserkezis)</td><td>1593</td></tr><tr><td>12. Quantum surface-response functions for nanoplasmonics: Feibelman parameters (Christensen et al.)</td><td>1596</td></tr><tr><td>13. Computational electrodynamics with Feibelman surface-response formalism (Hohenester et al.)</td><td>1602</td></tr><tr><td><b>Part III—Nonlocal effects beyond free-electron bulk metals</b></td><td>1605</td></tr><tr><td>14. Nonlocal response in doped semiconductors (Wubs)</td><td>1605</td></tr><tr><td>15. Nonlocality in polar dielectrics (De Liberato)</td><td>1609</td></tr><tr><td>16. Nonlocal effects in graphene plasmonics (Gonçalves & García de Abajo)</td><td>1613</td></tr><tr><td>17. Nonlocal quantum gain, loss compensation, and plasmon amplification in graphene hyperbolic metamaterials (Hess & Tarasenko)</td><td>1618</td></tr><tr><td>18. Surface electrodynamics of crystalline noble metals (Cox & Mortensen)</td><td>1621</td></tr><tr><td>19. Concerted nonlinear and nonlocal electrodynamics in 2D materials (Jelver et al.)</td><td>1624</td></tr><tr><td>20. Nonlocal chirality in a twisted multilayer (Sánchez Sánchez et al.)</td><td>1628</td></tr><tr><td><b>Part IV—Engineered nonlocal responses in metamaterials and metasurfaces</b></td><td>1631</td></tr><tr><td>21. Strong spatial dispersion in metamaterials and metasurfaces (Tretyakov & Simovski)</td><td>1631</td></tr><tr><td>22. Nonlocal hyperbolic metamaterials and metasurfaces (Pakniyat & Gómez-Díaz)</td><td>1633</td></tr><tr><td>23. Nonlocal effects in transdimensional plasmonics (Bondarev et al.)</td><td>1638</td></tr><tr><td>24. Light management using structural nonlocality in nanorod metamaterials (Krasavin & Zayats)</td><td>1643</td></tr><tr><td>25. Nonlocal wave phenomena in spatio-temporal metamaterials (Alù)</td><td>1648</td></tr><tr><td>26. Spectral and angular control of light with nonlocal metasurfaces (Song & Brongersma)</td><td>1651</td></tr><tr><td>27. Nonlocal metalenses (Levy)</td><td>1655</td></tr><tr><td>28. Nonlocal metasurfaces for space compression (Long et al.)</td><td>1659</td></tr><tr><td>29. Electro-optic spatiotemporal nonlocal metasurfaces (Bozhevolnyi)</td><td>1663</td></tr><tr><td>30. Cascaded diffractive nonlocal metasurfaces for highly multifunctional meta-optics (Overvig)</td><td>1666</td></tr></tbody></table>
Part V—Nonlocality, topology and nonreciprocity, extreme geometries and fundamental limits

<table><tr><td>31. Role of nonlocality in topological metamaterials (Prudêncio & Silveirinha)</td><td>1671</td></tr></table>

<table><tr><td>32. Nonlocality in nonreciprocal plasmonics (Hassani Gangaraj & Argyropoulos)</td><td>1675</td></tr></table>

<table><tr><td>33. Nonlocal effects in singular plasmonic geometries unveiled by transformation optics (Huidobro et al.)</td><td>1679</td></tr></table>

<table><tr><td>34. The concept of overlapping nonlocality and the ultimate thickness limits of optics (Miller)</td><td>1683</td></tr></table>

Closing remarks and back matter

<table><tr><td>35. Conclusion and outlook (Monticone & Mortensen)</td><td>1687</td></tr><tr><td>Acknowledgments</td><td>1688</td></tr><tr><td>Disclosures</td><td>1688</td></tr><tr><td>References</td><td>1688</td></tr></table>

# 1. Introduction

FRANCESCO MONTICONE & N. ASGER MORTENSEN

In many ways, the remarkable progress in the field of optics over the past decades and centuries has been based on the ability to shape materials and progressively relax assumptions about their optical response. Consider, for example, a simple piece of glass. Its macroscopic electromagnetic properties—at least at the level relevant to applications such as lens design—are relatively simple and can be described by linear, time-invariant, isotropic, local, nonmagnetic, and passive constitutive parameters: specifically, a scalar frequency-dispersive electric permittivity and the free-space magnetic permeability. Significantly greater complexity can be unlocked and accessed by relaxing these assumptions.

For instance, breaking the assumption of linearity gave rise to the field of nonlinear optics, which gained prominence with the invention of the laser and has had profound implications for both fundamental science and technology. Relaxing the assumption of an isotropic response leads, among other things, to polarization-dependent functionalities, while more complex linear constitutive relations—incorporating effects such as negative constitutive parameters, magneto-electric coupling, chirality, nonreciprocity, and biaxisotropy—have been central to the field of metamaterials, enabling unprecedented control over electromagnetic waves. Furthermore, relaxing the assumption of time-invariance is currently one of the most active areas of research in applied electromagnetics, photonics, and wave physics, allowing for the breaking of frequency/energy conservation and time-reversal symmetry. These examples illustrate how efforts to challenge any of these fundamental assumptions have sparked and motivated entirely new research fields.

Yet, among these assumptions, spatial locality has arguably been the least explored. This is not to say that nonlocal effects have been overlooked until recently. Early studies date back to the mid-20th century with even earlier seminal contributions by various authors (for further details, see Section 2). By the early 1960s, the basic physical mechanisms underlying optical nonlocality in natural materials were well understood, as in the notable work by Hopfield and Thomas [1] (the former gained recognition also for his work on physical computing and artificial neural networks, for which he won the 2024 Nobel Prize in Physics). Part of the Introduction of this paper is worth quoting here for its clarity [1]: In [the local] approximation, the dielectric polarization P within a small volume of radius $r_0$ ($r_0 \ll \text{any wavelength involved}$) depends only on the value of the electric field inside this volume (at the present time and in the past) and is not explicitly dependent on the electric field or other parameters outside
the volume under consideration. The term “spatial dispersion” has been used to apply to dielectric behavior for which the local description is not valid. In general, spatial dispersion refers to the wave-vector dependence of the dielectric constant. Implicitly contained in the supposition of local dielectric behavior is the neglect of the transport of energy by any mechanism other than electromagnetic waves. When energy transport by other mechanisms must be considered anomalous (nonlocal) dielectric behavior results, often accompanied by new physical phenomena.

Although the fundamental ideas of nonlocality were already present in these early works, it was only with the advent of the field of plasmonics and nano-optics that their significance became clearly recognized. In particular, as research in metamaterials and plasmonics started uncovering extreme forms of light-matter interactions, such as ultra-strong field confinement and enhancement in plasmonic nanostructures, it became evident that the traditional local approximation was insufficient. More sophisticated models were required to regularize and accurately describe how these systems respond to tightly confined electromagnetic fields, which contain large wavevector components in their spatial spectrum. This need is even more evident today in the study of nonreciprocal, topological, and time-varying systems, especially those made of continuous media with negative constitutive parameters, whose overall response often depends on the large-wavevector behavior of the system [2]. As another example, nonlocal response is also important for quantum-fluctuation phenomena [3] where spatial correlations in nonequilibrium quantum systems—often overlooked in local thermal equilibrium approximations—significantly influence interactions like quantum friction and near-field radiative forces [4]. In many cases, local models can lead to incorrect or even nonphysical predictions, whereas nonlocal models correct the asymptotic response for large wavevectors, ensuring a gradual vanishing of the polarization response.

In this context, nonlocality is both a *curse and a blessing*. On one hand, it resolves the unphysical issues of local models and predicts new physical effects that would otherwise be missed by the local approximation. On the other hand, it also imposes fundamental limits on how strong and localized light-matter interactions can be.

In parallel to the exploration of nonlocality in natural materials, the field of metamaterials has enabled the realization of strong *artificial* nonlocality, emerging—at the “effective medium” level—from long-range interactions across multiple unit cells or from the response of multipolar meta-atoms. In this case, nonlocality is clearly a “blessing”, offering a powerful new tool to engineer strongly wavevector-dependent functionalities. “Nonlocality engineering” adds wavevector as a new degree of freedom in the metamaterial design toolbox, enabling a new degree of control over the angular response of incident light, strong angular/frequency selectivity, space compression, multifunctional devices, and much more. Notably, nonlocal metamaterials or metasurfaces with a tailored spatial impulse response can even be used to perform mathematical operations on an input wavefront, effectively using artificial nonlocality for wave-based physical computing (an incidental but intriguing parallel to the career of one of the pioneers of this field, John J. Hopfield.)

The goal of this Roadmap is to review the current state of this emerging frontier in optics and photonics, highlighting both the challenges and opportunities of nonlocality—its *curse and blessing*. By surveying recent developments across various aspects of nonlocality, we also aim to define new directions that will drive this field forward, opening pathways for fundamental scientific discoveries and technological advancements.

Toward these goals, we have assembled a large collection of contributions from many of the leaders of this field. The Roadmap consists of thirty-four sections, structured into five parts:

* Part I starts with a general introduction of the concept of nonlocal response in the interaction of electromagnetic waves with matter (Section 2) and an introduction to semi-classical
hydrodynamic models (Section 3), followed by three sections that explore some of the
foundational aspects of this field in greater detail (Sections 4, 5, 6).

* Part II introduces the continuum framework for nonlocal plasmonics (Section 7), and continues with discussions of quantum hydrodynamics (Section 8). This is followed by contributions on density-functional theory approaches in the jellium model (Section 9) and in atomistic representations (Section 10). Implications of nonlocality for light-matter interactions are also discussed (Section 11). Finally, aspects of the surface-response formalism are discussed (Sections 12 and 13).

* Part III highlights important directions that go beyond the nonlocal response of free-electron bulk metals, ranging from the response of doped semiconductors (Section 14) and polar dielectrics (Section 15) to the nonlocal behavior of 2D plasmons in graphene and noble-metal surface-states (Sections 16, 17 and 18), nonlinear and nonlocal effects in 2D materials (Section 19) and nonlocal chirality in moiré multilayer systems (Section 20).

* Part IV explores the principles, design strategies, and emerging applications of nonlocal metamaterials and metasurfaces, highlighting their potential to transform modern optics and photonics. Part IV starts with a general discussion of strong spatial dispersion in metamaterials and metasurfaces (Section 21) and then explores the role of nonlocality and the new opportunities enabled by nonlocality engineering in a wide range of meta-structures for several applications, including hyperbolic metamaterials and metasurfaces (Section 22), “transdimensional” plasmonic thin films (Section 23), nanorod metamaterials for light management (Section 24), spatio-temporal metamaterials (Section 25), metasurfaces with high spectral/angular control (Section 26), nonlocal metalenses (Section 27), metasurfaces for space compression (Section 28), electro-optic spatiotemporal nonlocal metasurfaces (Section 29), and multi-functional meta-optics (Section 30).

* Part V delves into the implications of combining nonlocality with other advanced concepts to investigate new physics, extreme scenarios, and fundamental limits in photonics. Topics include the role of nonlocality in topological metamaterials (Section 31) and nonreciprocal plasmonics (Section 32). This section also features contributions on combining transformation optics and nonlocal models to study singular geometries (Section 33) and it explores the idea of “overlapping nonlocality”—a concept related to, but distinct from, standard nonlocality—and its implications for the ultimate thickness limits of optical systems (Section 34).

The Roadmap concludes with some brief remarks by the lead authors (Section 35).

## Part I—Fundamental frameworks

### 2. Linear-response formalism for nonlocal electrodynamics

N. ASGER MORTENSEN & FRANCESCO MONTICONE

#### *Overview*

The early fundamental contributions to the concepts of nonlocal response and spatial dispersion in condensed matter physics and electromagnetism arose primarily in the mid-20th century. These ideas expanded on classical theories of material response to electromagnetic fields by considering the fact that the response at a given point in a material might depend on the electromagnetic fields at other, nonlocal points [5]. However, this nonlocal dependence is limited to a vicinity of the initial point, with the range governed by the underlying microscopic dynamics of the system.
The numerous independent developments in this field make it challenging to present a chronological or cohesive review that fully honors the significant contributions across quantum physics, condensed-matter physics, plasma physics, and electrodynamics. The study of nonlocal phenomena is inherently tied to the hydrodynamic framework, including the Navier-Stokes equations of classical continuum mechanics. Notably, Madelung's seminal work established a vital connection between the Hamiltonian formalism of quantum mechanics and the hydrodynamic formalism of classical physics [6]. Equally important are Bloch's contributions to the hydrodynamic descriptions of electrons in solids [7], Landau's work on the damping of longitudinal waves in plasmas [8], Lindhard's quantum mechanical calculations of the nonlocal dielectric response function of electrons in solids [9], and the development of the Kubo formalism for analyzing the linear response of quantum systems [10]. Each of these contributions fundamentally addresses the nonlocal response of the underlying electron systems.

**The common implicit assumption.** As a general theme, these developments have largely relied on the common and often intuitive assumptions of a linear and local response of matter to light fields [11]. Mathematically, Maxwell's equations are augmented by linearized and spatially local constitutive relations between the relevant fields, such as $D(r) = \varepsilon(r)E(r)$ for the displacement field $D(r)$ arising from an applied electric field $E(r)$ in the spatial point $r$. The use of Ohm's law $J(r) = \sigma(r)E(r)$ is another such example where the induced current density $J(r)$ is arising from an applied electric field $E(r)$ in the spatial point $r$. However, there is no fundamental reason to assume that the material's response will be confined to the exact point of perturbation.

## Current status

**Beyond the local-response approximation.** Building on the concepts of Taylor expansions, the Volterra series offers a systematic approach to phenomenologically describe spatio-temporal nonlocal effects. Assuming a system response $\mathcal{L}$ associated with a system input $\mathcal{F}$, i.e. $\mathcal{L}[\mathcal{F}]$, the linear contribution in the Volterra series is formally expressed as [12,13]

$$ \mathcal{L}(t, \mathbf{r}) = \int_{-\infty}^{t} dt' \int d\mathbf{r}' \chi(t - t'; \mathbf{r}, \mathbf{r}') \mathcal{F}(t', \mathbf{r}'), \quad (1) $$

where $\chi(t - t'; \mathbf{r}, \mathbf{r}')$ is the spatio-temporal nonlocal response function, with the spatial and temporal integrals extending over all space and time, naturally subject to the restrictions by causality [12]. Figure 1 illustrates the spatio-temporal domain, where the response function is finite (shaded area), while being strictly zero in the noncausal part. The “time cone” is defined by the bounds on the speed $v$ of propagation in the medium, being naturally bounded by the speed of light $c$.
giving the common LRA result

$$ \mathcal{L}(t, \mathbf{r}) \approx \int_{-\infty}^{t} dt' \chi(t-t', \mathbf{r}) \mathcal{F}(t', \mathbf{r}). \qquad (2b) $$

The remaining temporal integral is associated with the frequency-dispersive response, subject to KK conditions imposed by causality.

<figure><img src="image_3.png" /><figcaption>**Fig. 1.** Schematic illustration (inspired by Ref. [12]) of the constraints of causality on the response function $\chi(t-t', \mathbf{r}, \mathbf{r}')$ and the spatio-temporal integrals in Eq. (1). The causal regime (gray-shaded) represents the region where the response function can take non-zero values, while the non-causal regime (orange-shaded) indicates where the response function is strictly zero. The “time cone” is determined by the characteristic speed $v$ of propagation in the medium, ensuring no response occurs from points $\mathbf{r}'$ located farther from the origin $(t, \mathbf{r})$ than $v|t' - t|$. The finite range of the spatially nonlocal response (associated with intrinsic dynamics of the matter) is denoted by $\xi_{\text{NL}}$.</figcaption></figure>

**Translational invariance.** For translationally invariant media, it is useful to Fourier transform Eq. (1) with respect to the temporal and spatial coordinates, resulting in a form that is particularly familiar to the communities of condensed-matter physics and solid-state spectroscopy [22]

$$ \mathcal{F}(\omega, \mathbf{q}) = \chi(\omega, \mathbf{q}) \mathcal{F}(\omega, \mathbf{q}), \qquad (3) $$

with the response function $\chi(\omega, \mathbf{q})$ now explicitly exhibiting both its temporal and spatial dispersion through its explicit dependence on both the temporal frequency $\omega$ and the spatial wave vector $\mathbf{q}$.

Exploiting light fields with the dispersion relation $\omega(q) = cq$, there has been a natural preference to account for light-matter interactions within the LRA, which in Fourier space amounts to the approximation $\chi(\omega, \mathbf{q}) \approx \chi(\omega, \mathbf{q} \to 0)$. This limit is well motivated by the fact that commonly $q\xi_{\text{NL}} = \omega\xi_{\text{NL}}/c \ll 1$ for the most encountered optical materials. In many cases, the momentum of optical photons is too low, and the nonlocal length scale is too short for nonlocal effects to have any significant impact. Alternatively, this can be understood as the speed of light $c$ being much greater than any intrinsic velocities $v$ characteristic for excitations in matter. Consequently, in most practical scenarios, common optical materials can be approximated as homogeneous and local in their response to light fields. This assumption is largely justified by the small size of the constituent atoms and their dense arrangement
within solids, where all microscopic length scales are much smaller than the wavelength of light.

**Beyond the low-q regime.** What would happen if translational invariance were broken or if we could artificially increase the separation between atoms? In the first scenario, where the system has small but finite dimensions $a$, the relaxation of Noether's conservation laws [23] would pragmatically allow probing at finite wavevectors $q \sim 1/a$. In the second scenario, introducing an artificial lattice constant $a$, the associated discrete translational invariance would enable *Umklapp* scattering processes [24] among photonic Bloch states [25], which also facilitate probing at finite wavevectors $q \sim 1/a$ (or higher multiples). Manipulating and engineering a finite $q$—response is no longer merely aspirational—the advent of photonic nanostructures and metamaterials has turned this idea into reality. In recent decades, advancements have enabled the fabrication of increasingly smaller nanostructures and intricate photonic crystals and metamaterials with engineered lattice constants, along with improved capabilities for exploring light fields at sub-wavelength scales, including various photon and electron-based near-field techniques [26–30]. As a result, we are now moving beyond the traditional paradigm of probing only the $q \to 0$ limit when investigating naturally occurring nanostructures and artificially engineered metamaterials. This progress has enabled a deeper fundamental understanding of nonlocal light-matter interactions, including potential quantum effects within the underlying matter systems. Simultaneously, it has opened entirely new possibilities for engineering the complex optical responses of nonlocal metamaterials.

**Homogeneous media versus metamaterials.** The nonlocal responses of both homogeneous materials and metamaterials involve spatial dispersion, where field variations influence the optical response. In homogeneous materials, this nonlocality arises from intrinsic microscopic and mesoscopic properties, leading to wavevector-dependent effects. In metamaterials, nonlocality can emerge from the artificial periodicity of the composite with strong interactions between neighboring resonant meta-units, likewise giving rise to wavevector-dependent effects. These two interpretations—material-induced and resonance-driven nonlocality—are often discussed separately, yet they share a common underlying principle: both modify wave propagation by introducing additional modes and spatial dispersion at different scales. As noted by Podolskiy during the 2015 Faraday Discussions on nonlocality [31]:

> "The drastic difference between nonlocality in homogeneous media studied previously and nonlocality in metamaterials lies in the origin of additional waves. In homogeneous materials, the nonlocal response is attributed to the spatial dispersion of the material itself. In metamaterials, however, nonlocality appears at the 'effective medium' level; the response of every component of the metamaterial may remain local, while the granularity of the composite leads to the spatial dispersion of the effective permittivity."

Recognizing this dual nature of nonlocality is essential for establishing a unified conceptual framework for nonlocal photonics, as explored throughout this Roadmap.

**Short-range approximation.** Returning to Eq. (1) we finally offer a short-range version, which corresponds to retaining terms up to order $q^2$ in Eq. (3). Initially motivated by Ref. [14], this is enabled by the introduction of real-space moments of the short-range response function [15,16]. For simplicity, we consider an isotropic scalar response, and we get the following spatial expression

$$ \mathcal{L}(\omega, \mathbf{r}) \approx [\chi(\omega, \mathbf{r}) + \xi_{\text{NL}}^2 \nabla^2] \mathcal{F}(\omega, \mathbf{r}), \quad (4) $$

where the first term on the right-hand side emerges as a result of the LRA in Eq. (2a) (corresponding to $q \to 0$), while in the second nonlocal term with the Laplacian operator, $\xi_{\text{NL}}$ is the range of the nonlocal response, being formally defined as the second spatial moment of the response function [15,16]. This scalar expression can be generalized to vectorial fields and non-scalar responses,
with the nonlocal hydrodynamic model serving as a key example, where a nonlocal Ohm’s law can be formulated in the form [32]

$$ \mathcal{L}(\omega, \mathbf{r}) \approx [\chi(\omega, \mathbf{r}) + \xi_{\text{NL}}^2 \nabla[\nabla \cdot \mathbf{r}]] \mathcal{F}(\omega, \mathbf{r}). \quad (5) $$

The change from a Laplacian to a gradient-of-divergence nonlocal correction term introduces the possibility for additional longitudinal waves [33,34], which is a main difference between nonlocality in homogeneous (conductive) media versus metamaterials (with local constituents).

## Concluding remarks

Having attempted to offer a general introduction to the concept of nonlocal response in the interaction of electromagnetic waves with matter, the remainder of this first part of the Roadmap explores some of these foundational aspects in greater detail, and for further insight, we direct readers to key textbooks [5,35], as well as reviews on the nonlocal electrodynamics of metals [16,34,36–38] and plasma [39–41].

**Funding.** Danmarks Grundforskningsfond (DNRF165); Air Force Office of Scientific Research (FA9550-22-1-0204); Office of Naval Research (N00014-22-1-2486).

# 3. Nonlocal hydrodynamics

ANTONIO I. FERNÁNDEZ-DOMÍNGUEZ ID, YU LUO ID, XUEZHI ZHENG & CHRISTOS TSERKEZIS ID

## Overview

**The hydrodynamic Drude model.** The hydrodynamic Drude model (HDM) constitutes the first and most straightforward realization of a nonlocal permittivity of metals, with a lifetime that is approaching a century [7]. It is the direct descendant of the standard Drude theory developed 30 years earlier [42], which describes free electrons in the metal as an ideal gas, whose atoms move freely until they collide with the heavier nuclei (or, with much smaller probability, with each other). The corresponding angular-frequency-dependent relative permittivity $\epsilon$ is then [24]

$$ \epsilon_m(\omega) = \epsilon_\infty(\omega) - \frac{\omega_p^2}{\omega(\omega + i\gamma)}, \quad (6) $$

where $\omega_p = \sqrt{ne^2/(m_e\epsilon_0)}$ is the plasma frequency (here $\epsilon_0$ is the vacuum permittivity, $e$ is the elementary charge, $m_e$ the mass of electrons and $n$ their density), while $\gamma$ is a phenomenological damping rate accounting for the aforementioned collisions. Any contribution to the permittivity from bound electrons or other mechanisms (e.g. interband transitions) is included via $\epsilon_\infty$, which can be obtained by subtracting the free-electron contribution from experimental data, or by modeling any additional known mechanism as a Lorentzian; for a purely free-electron metal (as in alkali), $\epsilon_\infty = 1$.

HDM emerged as a generalization of the above permittivity to reflect, in an effective manner, the quantum character of electrons. As quantum particles, they obey Pauli’s exclusion principle, which inevitably turns the electron gas into a compressible one [7], obeying an equation of motion of the form [34,36,37]

$$ \frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla) \mathbf{v} = -\frac{e}{m_e} (\mathbf{E} + \mathbf{v} \times \mathbf{B}) - \gamma \mathbf{v} - \frac{\beta^2}{n} \nabla n. \quad (7) $$

Here, $\mathbf{v}(\omega, \mathbf{r})$ is the hydrodynamic velocity of the free-electron gas (in general, a function of position $\mathbf{r}$ and angular frequency, while $\mathbf{E}$ and $\mathbf{B}$ are the electric and magnetic field of the
incident electromagnetic (EM) field, respectively, which is the driving force for the motion of the electron gas. The hydrodynamic parameter $\beta$ is directly proportional to the Fermi velocity $v_F$ of the metal; in particular, for a three-dimensional (3D) electron gas, it takes the value $\beta^2 = v_F^2/3$ in the low-frequency limit ($\omega \ll \gamma$), while $\beta^2 = 3v_F^2/5$ in the high-frequency (optical) limit ($\omega \gg \gamma$) [44]. For more details about the derivation of these values, we direct the reader to [16,34,45]. The pressure term $\beta^2\nabla n/n$ in Eq. (7) can be interpreted as an additional force that tends to homogenize the electron density.

Using $J = en_0v$ to define the current density (where $n_0$ is the equilibrium electron density, while the driving EM field induces a perturbation $n_1$), the equation of motion Eq. (7), together with the continuity equation, result in the generalized Ohm's law for HDM,

$$ \frac{\beta^2}{\omega(\omega + i\gamma)} \nabla [\nabla \cdot J(\omega, \mathbf{r})] + J(\omega, \mathbf{r}) = \sigma(\omega)E(\omega, \mathbf{r}), \quad (8a) $$

where $\sigma$ is the Drude conductivity. This equation is to be solved coupled to the EM wave equation

$$ \nabla \times \nabla \times E(\omega, \mathbf{r}) = \frac{\omega^2}{c^2} \varepsilon_{\infty} E(\omega, \mathbf{r}) + i\omega\mu_0 J(\omega, \mathbf{r}), \quad (8b) $$

where $\mu_0$ is the vacuum permeability, and $c = 1/\sqrt{\varepsilon_0\mu_0}$ the speed of light in vacuum.

The aforementioned inclusion of quantum pressure has several immediate implications. First of all, it introduces a degree of spatial extent of induced charges, rather than being strictly localized at the metal-environment interface. But this means that the centroid of induced charge is pushed inwards (so-called electron spill-in); thinking in terms of classical harmonic oscillators, this reduced path increases the strength of the restoring force experienced by the electron cloud, leading thus to an increase in resonant frequency or, in terms of wavelength, a blueshift of the localized surface plasmons (LSPs) as compared to the predictions of LRA. This is shown in Fig. 2(a), comparing the black dotted line (LRA) with the light-blue solid one (HDM), for a sodium (Na) sphere of radius $R = 5$ nm, described by the Drude model in Eq. (6)

<figure><img src="image_4.png" /><figcaption>**Fig. 2.** Effect of nonlocality to the far-field spectra. A Normalized (to the geometrical cross section) extinction cross section for (a) Na sphere of radius $R = 5$ nm, as calculated within LRA (dotted black line), HDM (solid light-blue line), GNOR (dash-dotted dark-blue line), and SRF (dashed red line), in the energy window of the dipolar LSP mode. (b) Normalized extinction cross section (in logarithmic scale) for an infinitely-long Ag cylinder, within LRA (dashed black lines), HDM (solid green lines), and HDM with a curl-free approach (dotted red lines), for two different cylinder radii. The inset shows the blueshift of the main LSP as a function of radius. Adapted with permission from Ref. [43] (Copyright © 2011 American Physical Society).</figcaption></figure>
with $\hbar\omega_p = 5.89$ eV, $\hbar\gamma = 0.1$ eV, and $\varepsilon_\infty = 1$. As the nanoparticle (NP) size becomes smaller, the deviation between HDM and LRA predictions increases, see, e.g. Fig. 1(a) in Ref. [46]. On the one hand, the LRA spectra converge to the quasistatic approximation ($\omega_p/\sqrt{3}$ for a metal with $\varepsilon_\infty = 1$), since the NP behaves increasingly more like a classic point dipole; on the other hand, the nonlocal length and the shift of the centroid of the induced charge become more significant compared to the overall NP size, leading to steadily stronger resonance blueshifts.

The physical origin of the electron spill-in described above can be linked to the optical excitation of longitudinal waves (plasmons) in metal NPs modeled using the HDM. Longitudinal waves are not allowed in metals described by Eq. (6), as the permittivity vanishes only at the plasma frequency, well above the optical range of the EM spectrum. On the contrary, within the HDM, the convective motion of the electron gas induces a spatial dispersion on the longitudinal component of the metal permittivity, while leaving the transverse part as given by Eq. (6). Thus, the dispersion of the longitudinal waves within the metal acquires the form [47]

$$ \varepsilon_{mL}(\omega, k_{mL}) \equiv \varepsilon_{\infty}(\omega) - \frac{\omega_p^2}{\omega(\omega + i\gamma) - \beta^2 k_{mL}^2} = 0, \quad (9) $$

where $k_{mL}$ is the longitudinal wavenumber in the metal. Equation (9) shows that for $\omega<\omega_p$, the longitudinal waves are evanescent, decaying exponentially in the radial direction into the NP, giving rise to the aforementioned spill-in. In the limit $\omega \ll \omega_p$, the thickness of the spilled-in surface charges can be estimated as [48]

$$ \delta_{mL} \sim \frac{1}{k_{mL}} \approx \frac{\omega_p}{\sqrt{\varepsilon_\infty \beta}}, \quad (10) $$

which yields $\delta_{mL} \sim 0.1$ nm for realistic parameters (see below), in good agreement with the Thomas—Fermi wavelength in noble metals [49]. This spatial encoding of the HDM has been exploited for the development of effective, simplified models of nonlocal effects in the optical response of metal nanostructures [50], which take into account the spill-in of induced charges at their boundaries. Back to Eq. (9), at $\omega>\omega_p$, the longitudinal waves become propagating in all directions, leading to longitudinal bulk plasmons resonances, as shown by the green lines in the upper panel of Fig. 2(b). In this case, for a gold (Au) cylinder with $\omega_p = 8.81$ eV, $\hbar\gamma = 0.075$ eV, $\varepsilon_\infty = 1$ and $v_F = 1.39 \times 10^6$ m/s [43]. These bulk plasmon modes have been reported already since the 1970s [51,52], but their resulting resonances are substantially weaker (notice the logarithmic scale in the figure) and largely affected by material losses.

**Boundary conditions.** Finally, the addition of an extra differential equation, and the existence of longitudinal waves, directly leads to the need of an additional boundary condition, for the optical response of a nonlocal sphere to be calculated analytically or numerically. Boundary conditions have been discussed in detail in the past, with different choices often leading to different physics, as can be seen for example in the spectra of Fig. 2(b), and the striking differences between the so-called curl-free model [53] (red dotted line) and the implementation following the so-called Sauter boundary condition [54] (green solid line). This, more natural choice, is based on the clear, physical understanding that the induced-current component normal to the metal-environment interface must be equal to zero, i.e. there is no electron spill-out. For $\varepsilon_\infty = 1$ this boundary condition yields naturally the continuity of the normal component of the electric field as well, while for $\varepsilon_\infty \neq 1$, differences emerge between imposing current or electric field continuity, which can be directly related to the (to some degree) artificial distinction between bound and free electrons inherent to the Drude model and HDM [55].

The so-called hard-wall boundary condition introduced above, i.e. the premise that the work function of the metal is so high that the electron-density profile terminates abruptly at the assumed interface [56] (see Ref. [57] for a detailed discussion) has been the cornerstone of criticism
towards the traditional HDM, as it fails, by construction, to describe the response of metals where spill-out is considerable, such as Na [58]. Spill-out shifts the centroid of induced charge outwards, rather than inwards, thus leading to resonance redshifts, as shown in Fig. 2(a) with a dashed red line—here, the spectrum is calculated with the surface-response formalism (SRF) based on Feibelman parameters (see e.g. Reference [59] as well as Secs. 12 and 13 of this Roadmap). Ways to relax the hard-wall boundary condition within the framework of HDM have also been proposed, and are currently in use [60,61]. Finally, we should also mention that HDM, as described here, also fails to account for the well-known size-dependent broadening of the resonances, already measured experimentally back in the 1970s [62,63]. This is well reproduced by SRF, while in the context of HDM it is dealt with by the generalized nonlocal optical response (GNOR) model [64] [see dark-blue dashed-dotted line in Fig. 2(a)], which extends HDM by allowing $\beta$ to become complex, with the imaginary part describing electron diffusion (and thus increased collisions and damping), on top of their convective motion.

*Current status, challenges, and opportunities*

**HDM solutions and other theoretical models.** From the above discussion, it becomes clear that HDM and its generalizations introduce an additional dimension to the analytical or numerical treatment of electrodynamic problems. The existence of longitudinal waves, and the additional boundary condition that is needed, have to be implemented in any solver. For canonical architectures such as infinite planes, spheres, and infinitely long cylinders, analytic solutions that describe the reflection, scattering, and absorption of light when spatial nonlocality is taken into account can be derived, and this was done already in the 1970s [51,65,66]. But, in order to be of practical use to the nanophotonics community, this treatment needed to be generalized to account for more complex particle shapes and configurations, exploited in plasmonics designs. This triggered a rapid growth of nonlocal computational electrodynamics in the late 2000s and early 2010s. McMahon and collaborators introduced the simplified curl-free incarnation of HDM into finite-difference time-domain (FDTD) calculations, essentially by Fourier-transforming the corresponding version of Eq. (8a) and solving it coupled to the wave equation in the time domain [53,67]. Figure 3 displays the squared modulus of the displacement field at the brightest localized plasmon resonance sustained by cylindrical, square, and triangular Au nanowires with size 50 nm (see the inset sketching the incident fields). Top and bottom panels present nonlocal and local FDTD calculations, respectively, illustrating the optical excitation of spurious longitudinal plasmons below the plasma frequency in the former [53]. Toscano et al. proceeded directly in the frequency domain, using the finite-element method (FEM) as their tool for solving coupled differential equations [43,68] using a commercial solver, while Hiremath et al. presented the full details of a rigorous formulation based on Nédélec finite elements [69]. Figure 3(b) shows FEM calculations of the energy loss rate experienced by swift electrons interacting with Na (top) and Au nanospheres [34]. As discussed in the context of experiments, electron-energy-loss spectroscopy (EELS) allows the characterization of plasmonic resonances with excellent spatial and spectral precision, and thus, it is also an ideal testbed for nonlocal approaches, here comparing GNOR versus local predictions.

Combining the strengths of the two methods, HDM was also introduced in discontinuous-Galerkin time-domain codes [73,74] (see also Sec. 7 of this Roadmap). This method, like FDTD and FEM, is volume-based, which comes with an increased computational cost, since the entire simulation domain needs to be discretized. On the contrary, boundary-element methods (BEM) only need to describe the interface between two different media, and are thus more efficient in terms of computational time and resources. Figure 3(c) illustrates the implementations of HDM into BEM, presenting the extinction cross-section of ellipsoidal NPs with different eccentricity [70], revealing the strong dispersion of the bright and dark plasmonic modes that they support. Other BEM descriptions of quantum nonlocal effects in metallic nanostructures can be found
in Refs. [29,75,76]. Although most of the theoretical research in this context has focused on finite systems comprising nanometric-sized metallic elements (such as antennas or cavities), HDM approaches have also shed light into the optical performance of extended devices. Thus, Fig. 3(d) shows the comparison between local (top) and nonlocal (bottom) predictions for the near-field of a nanotip. These maps reveal that spatial dispersion effectively blurs the structure boundaries, reducing the sensitivity of EM fields to surface roughness with sizes of the order

<figure><img src="image_5.png" /><figcaption>**Fig. 3.** Nonlocal modeling of plasmonic devices: (a) Spurious longitudinal plasmons in NPs of different sizes are absent (present) in local (curl-free HDM) FDTD calculations in the top (bottom) maps. Adapted with permission from Ref. [53] (Copyright © 2010 American Physical Society). (b) EELS spectra obtained through the GNOR and local FEM models of metallic NPs. Adapted with permission from Ref. [34] (Copyright © 2015 Institute of Physics). (c) BEM-HDM calculations of the extinction spectra of ellipsoid NPs of increasing eccentricity, representing data similar to Ref. [70]. (d) FEM calculations that show the nonlocal robustness to surface roughness of the plasmonic modes propagating in metal tips. Reprinted (adapted) with permission from Ref. [71] (Copyright © 2012 American Chemical Society). (e) Transformation optics HDM description of the absorption spectrum of Au spheres separated by gaps between 10 and 0.01 nm. The left panel sketches the mapping procedure in the analytical calculations. Adapted with permission from Ref. [72] (Copyright © 2014 National Academy of Sciences).</figcaption></figure>
of $\delta_{mL}$ in Eq. (10). Thus, as a result of quantum nonlocality, the optical modes at the tip suffer lower scattering losses in their propagation towards its apex.

In the literature, the versatile and powerful numerical tools introduced above are complemented by analytical efforts to introduce nonlocal refinements to local EM models. One of these exploited the theoretical framework of transformation optics to describe the plasmonic fields sustained by nanostructures in the quasi-static regime [77]. This approach is illustrated in Fig. 3(e), which presents the calculation of the absorption cross-section of 30 nm radius Au spherical dimers separated by nano and subnanometric gaps [72]. The left sketch shows that, through the inversion of two parallel HDM metal surfaces, sphere dimers with a $\beta$ parameter that varies in space are obtained. The asymmetric orange shell renders $\Delta d' \approx \delta'_{mL} \propto \beta'(\mathbf{r}')$ [48], as defined in Eq. (10) (the primes indicate the dimer coordinate system). The correct treatment of these mapping effects, leveraging the geometric interpretation of HDM in Ref. [50], yields the absorption spectra in the right contour plot, which reproduces correctly the main features of the nonlocal phenomenology reported experimentally (see later discussion on experiments). Another analytically-inspired approach sheds light into the emergence of electron tunneling effects in these plasmonic gaps. This consists in the modeling of the medium between the metallic elements through a permittivity of the form of Eq. (6) with a very large $\gamma$, parametrized with electron transmission (ab initio) calculations [78].

**Experiments.** The investigation of nonlocal effects in metal nanostructures has attracted much experimental attention in recent years. Figure 4 presents a general perspective on the different optical characterization approaches and sample geometries explored, as well as the nonlocal phenomenology reported in each experimental configuration. EELS has been instrumental for our current understanding of nanophotonic systems [29,79]. It has also been a crucial tool enabling the assessment of surface quantum and nonlocal effects in metal NPs with morphologies presenting sub-nanometric features. Figure 4(a) shows EELS measurements of the resonant frequency of the dipolar plasmonic mode sustained by silver (Ag) spheres with decreasing diameter [80]. The EELS data demonstrate a modal blueshift of ~ 0.5 eV from 20 nm to 1 nm diameter, and a much lower effect on the bulk plasmon resonance at higher frequencies. The authors of this study linked these observations to quantum confinement effects in the NP conduction-electron wavefunctions. Their findings were also completely compatible with a hydrodynamic model of the Ag permittivity as described above, as demonstrated by a subsequent report [81]. Similar EELS experiments shed light into the disappearance of higher order plasmon modes in isolated Au NPs, and the occurrence of a different nonlocal phenomenology, modal redshift with decreasing junction cross-section, in connected Au bowties [82].

Apart from monolithic structures, nanometric gaps (at the encounter between two metallic elements) have been probably the platform where quantum nonlocality has been investigated most intensively. Although EELS has also been employed in this geometry [86,87], dark-field spectroscopy has enabled a deep far-field exploration of gap plasmonic modes. Note that the systematic far-field probing of nonlocal effects in nanometric NPs is not possible due to their diminutive optical cross-section [88]. Figure 4(b) plots the position of the lowest plasmonic resonance sustained by 30 nm radius Au nanoparticle-on-mirror (NPoM) samples in which the gap length was set through molecular spacers generated by two different techniques: layer-by-layer deposition (from 25 to 3 nm) and self-assembled monolayers (from 2 to 0.5 nm). These experiments revealed that, contrary to the continuous redshifting with decreasing gap predicted by local EM calculations, a saturation of the plasmon resonance was observed, in quantitative agreement with a nonlocal hydrodynamic model [83] with realistic parameters ($\beta \sim 10^6$ m/s). These experiments indirectly established the limits of plasmon-assisted field enhancement within plasmonic gaps, and proved that NPoM structures function as far-field rulers with Ångström-level precision [89].

The investigation of smaller and smaller plasmonic gaps revealed a regime where the quantum nature of conduction electrons in metals manifests in a completely different fashion. When the
gap size enters the Ångström range and becomes comparable to the electronic spill-out length, the conduction electrons can tunnel between NPs, giving rise to charge transfer phenomena that alter the plasmonic modes at the gap. Figure 4(c) presents experimental evidence of the onset in this tunneling regime between two Au atomic-force-microscopy tips [84]. The color plot renders dark field scattering versus incident wavelength and force applied to approach the nanotips. The lateral panel plots the electrical conductance through the gap, whose abrupt step reveals the point of contact between the tips. Through the comparison against full quantum and quantum-corrected [78] classical simulations, the authors of the study set the onset of the tunneling regime at $d_{QR} = 0.31$ nm (8 nN applied force), where the dispersion of the gap plasmonic modes changes abruptly. At larger gap sizes (lower forces), the modes redshift with decreasing distance in the same way as reported in Ref. [83] due to the emergence of strong spatial dispersion in the Au permittivity. At gaps smaller than $d_{QR}$, the modes blueshift with applied force, which was identified as the fingerprint of the occurrence of charge transfer between the tips. Apart from EELS and dark-field scattering, the impact of quantum tunneling in plasmonic gaps has been

<figure><img src="image_6.png" /><figcaption>**Fig. 4.** Experimental demonstration of quantum nonlocal effects in metal nanostructures: (a) Blueshift of the EELS-determined plasmon frequency with decreasing nanoparticle diameter, with the inset depicting bulk resonance energies. Adapted with permission from Ref. [80] (Copyright © 2012 Nature Springer). (b) Deviation from local predictions of the dark-field scattering resonances of a 30 nm radius NPoM with gap sizes in the nanometric range. Adapted with permission from Ref. [83] (Copyright © 2012 American Association for the Advancement of Science). (c) Scattering spectra as a function of the distance (measured through the force exerted) between two atomic force microscopy nanotips, revealing the onset of electron tunneling effects. Adapted with permission from Ref. [84] (Copyright © 2012 Nature Springer). (d) Imaginary versus real part of the effective index of the gap plasmon modes propagating between Au surfaces, revealing the impact of nonlocal damping. Adapted with permission from Ref. [85] (Copyright © 2022 Nature Springer).</figcaption></figure>
assessed using nonlinear optical processes. These range from surface enhanced Raman scattering (SERS) from thiophenol monolayers attached to the surface of Au nanodisk dimers [90] to third-harmonic generation from alkanethiol monolayers at the gap of Au NPoM cavities [91].

More recently, advances in nanofabrication techniques have allowed the scrutiny of quantum nonlocal effects in extended systems, which involve interfacing plasmonic surfaces within areas larger than the operating wavelength, orders of magnitude larger than the gap between NPs. Thus, the characterization of the dispersion characteristics of the extremely confined modes propagating within the nanometric gap, $t_d$, between two Au flakes could be carried out using scanning near-field optical microscopy [85]. Figure 4(d) plots the imaginary part of the effective mode index of these gap surface plasmons, $n_{GSP}$, versus its real part for $t_d$ between 2 and 20 nm and at an excitation wavelength of 1550 nm. The local-response approximation (LRA) predictions underestimate the ratio $Im\{n_{GSP}\}/Re\{n_{GSP}\}$. On the contrary, the GNOR (generalized nonlocal optical response) model, which incorporates a diffusive-like, imaginary part in the $\beta$ factor, reproduced the experimental observations accurately. These findings revealed that nonlocal damping, encoded in $Im\{n_{GSP}\}$, has a significant impact in propagating plasmon modes confined at the nanoscale.

Van der Waals heterostructures provide another platform for achieving extreme optical confinement. Notably, graphene supports surface plasmons at THz frequencies, with naturally transverse decay lengths on the order of a few nanometers. Furthermore, the dispersion characteristics of these plasmon modes can be tuned by adjusting the graphene carrier density. Coupling THz graphene plasmons to metal structures and other two-dimensional materials has enabled the characterization and electrical tailoring of quantum nonlocal effects in these systems [92,93]. Recently, the interface between inorganic semiconductors, such as indium antimonide (InSb) or gallium arsenide (GaAs), and metamaterial resonators has also been explored to investigate the emergence of nonlocal effects at THz frequencies [94,95]. While the influence of dielectric spatial dispersion in these systems differs significantly from that observed in nanostructured metals at optical frequencies, these experiments have demonstrated the universal nature of nonlocal EM phenomena.

## Future developments to address challenges

Given that in this Roadmap, the theoretical treatment of quantum nonlocal effects in nanophotonics systems will be discussed in great detail, we provide here only a brief overview of the broader prospects of the field, focusing on general insights beyond the specific topics explored in detail in the following sections. On the theoretical front, a comprehensive framework combining electrodynamics with ab-initio condensed matter physics is still lacking. Such framework will be essential to fully capture the quantum phenomenology found in the experiments briefly introduced above. These were conducted on samples of relevance from an optical perspective, in terms of size, material complexity and geometrical design. Another promising direction is the effective integration of nonlocal electromagnetics with quantum optical models, which would enable the unified description of the quantum nature of electrons and photons on the same footing. However, both of these avenues face significant challenges in terms of numerical complexity, demanding substantial efforts in the years ahead. On the experimental side, cleaner experiments are still needed to efficiently isolate quantum electronic effects in optical measurements, overcoming current constraints on sample design and detection sensitivity to provide unambiguous measurements of their strength in devices optimizing photonic functionalities ranging from light collection and concentration to light-matter interaction and quantum nano-optics.

## Concluding remarks

The hydrodynamic Drude model, presented here in its linearized form, remains a foundational framework for the investigation and interpretation of nonlocal phenomena. Secs. 8 and 7 delve
deeper into the hydrodynamic framework, addressing also physical and computational aspects of
self-consistency and nonlinearities in greater detail.

**Funding.** Danmarks Grundforskningsfond (DNRF165); Ministerio de Ciencia, Innovación y Universidades (PID2021-126964OB-I00); Otto Mønsteds Fond (24-12-1979).

4. Imaginary part of nonlocal permittivity: how it localizes absorption and imposes limits on field confinement and enhancement in plasmonics

JACOB B. KHURGIN ID

Overview

Space-time duality in photonics has been recognized since the formulation of Maxwell’s equations, but only in recent years has it been rigorously studied both theoretically and experimentally [96]. The duality between beam diffraction in space and the dispersion of short pulses in fibers has been demonstrated and applied to phenomena such as temporal lensing. When combined with nonlinearity, this duality leads to analogies between spatial and temporal solitons. More recently, space-time duality has been further explored through the introduction of temporal reflection and photonic time crystals, which serve as full analogs to their spatial counterparts [97]. Simultaneously, nonlocality—whether in natural materials like metals or in artificial metamaterials—has become a subject of significant interest [98]. Nonlocality is the spatial counterpart to delayed (non-instantaneous) temporal responses, and in the Fourier domain, it manifests as wavevector-dependent dielectric permittivity $\epsilon(\omega, k)$—or susceptibility $\chi(\omega, k)$—in addition to its ubiquitous frequency dependence. The permittivity of a material is a complex quantity, with the real part describing retardation and the imaginary part associated with decay or absorption. It is well established that the dispersion of the real part of permittivity is linked to the dispersion of the imaginary part via the Kramers–Kronig (KK) relations. Strong, and even anomalous, dispersion occurs near absorptive resonances, as first noted by Kundt in 19th Century [99]. This suggests that the presence of dispersion implies the existence of absorption somewhere in the frequency domain. But what is the spatial analog of this relation? The KK relations arise from the causality of susceptibility $\chi(\tau<0) = 0$ which is manifested in the Fourier transform being analytic in the upper complex half-plane. In the spatial domain, however, the nonlocal susceptibility behaves differently, $\chi_{\omega}(-\rho) = \chi_{\omega}(\rho)$, and a direct application of the KK transform is not possible. Since the Fourier transform of a symmetric function is real, the real and imaginary parts of complex nonlocal susceptibility $\chi_{\omega}(\rho) = \chi'_{\omega}(\rho) + i\chi''_{\omega}(\rho)$ have their own separate spatial Fourier transforms $X'_{\omega}(k) = F[\chi'_{\omega}(\rho)]$ and $X''_{\omega}(k) = F[\chi''_{\omega}(\rho)]$. Nevertheless, because the susceptibility $X_{\omega}(k)$ still contains frequency dependence and for a given $k$ KK relations hold true in temporal domain, it becomes possible under certain conditions to establish KK-like relationships between the real and imaginary parts in the spatial domain. It follows then that manifestation of strong spatial dispersion in a material that is nominally transparent at small wavevectors is usually accompanied by strong absorption somewhere at higher $k$'s. While the nonlocality of the real part of permittivity has been extensively studied, and several approximations have been developed to describe it effectively, its impact remains relatively limited. Typically, it results in shifts of resonant frequencies in various plasmonic and nanophotonic structures. In contrast, the effects of nonlocality on the imaginary part of permittivity are far more significant. These changes not only contribute to enhanced losses but also impose a limit on the field confinement enhancement in the aforementioned plasmonic and nanophotonic structures, as will be briefly discussed below.

Let us begin by examining the simplest scenario of nonlocality in a material with a single resonance frequency $\omega_0(k) = \omega_0(0) + ak^2/2$, where $a$ can be positive, as in the case of a Wannier-Mott exciton [1], or negative, such as in the case of an optical phonon [94]. Since $2ak$ is
a velocity these examples highlight the fundamental physical origin of nonlocality: the motion of the medium interacting with light. This interaction causes polarization induced at one point $r - \rho$ to propagate to another point $r$ giving rise to nonlocal susceptibility $\chi_{\omega}(\rho)$.

Consider the scenario of a positive a as shown in Fig. 5. For the fixed wavevector $k_1$, the susceptibility exhibits a Lorentzian $\chi_{\omega}(k_1) \sim [\omega_0(k_1) - \omega - i\gamma]^{-1}$ profile with a broadening $\gamma$. The real and imaginary components of susceptibility, visualized on the left side of Fig. 5, are clearly interconnected by the KK relation. At the same time, for fixed frequency $\omega_1 > \omega_0(0)$, the spatial dispersion of susceptibility takes the form of

$$X_{\omega_1}(k) \sim [\omega_0(0) + ak^2 - \omega_1 - i\gamma]^{-1} \sim [k^2 - k_0^2(\omega_1) - i\gamma/a]^{-1}, \quad (11)$$

where the “resonance” wavevector is $k_0^2 = [\omega_1 - \omega_0(0)]/a$. The spatial dispersion of $X'(k)$ and $X''(k)$ is illustrated at the bottom of Fig. 5 and the KK-like relation between them is apparent. Note, however, that for $\omega_1 < \omega_0(0)$, the $X_{\omega_1}(k)$ remains real for all wavevectors and a KK-like relation does not exist. We may conclude that unlike the temporal domain, the presence of absorption in the spatial domain always implies nonlocal dispersion, but the converse is not necessarily true.

<figure><img src="image_7.png" /><figcaption>Fig. 5. Spatial dispersion of resonant frequency in nonlocal dielectric. On the left – temporal dispersion of susceptibility. At the bottom – spatial dispersion of susceptibility.</figcaption></figure>

## Current status

Consider now the case of a metal such as gold (Au) or silver (Ag), where in the absence of scattering, absorption is forbidden at small wavevectors due to momentum conservation. However, once the wavevector approaches $k_0(\omega) = \omega/v_F$, where $v_F$ is a Fermi velocity, absorption becomes permitted as indicated by the tilted solid arrow in Fig. 6(a). This phenomenon is called Landau damping (LD), being also further discussed in Sec. 5. In one dimension the dispersion of resonant frequency is $\omega_0(k) = v_F k$, and a KK-like relation similar to the one depicted in Fig. 5 can be established. In three dimensions, the resonance broadens as absorption is now allowed for all $k \ge k_0(\omega)$, and the integration of 11 is expected to show that KK-like relation will demonstrate its preservation. Indeed, according to Lindhard the spatial dispersion of dielectric permittivity is
given by [9]

$$ \varepsilon(\omega, \mathbf{k}) = \varepsilon_b + \frac{3\omega_p^2}{v_F^2 k^2} \left[ 1 - \frac{\omega}{2v_F k} \ln \frac{\omega + v_F k}{\omega - v_F k} \right], \quad (12) $$

where $\varepsilon_b$ is the dielectric constant due to bound electrons and $\omega_p$ is the plasma frequency. At $k=0$ the permittivity is $\varepsilon(\omega, 0) = \varepsilon_b - \omega_p^2/\omega^2$ and we can write for $\Delta\varepsilon(\omega, q) = \varepsilon(\omega, q) - \varepsilon(\omega, 0)$ that

$$ \Delta\varepsilon'(\omega, q) = \frac{\omega_p^2}{\omega^2} \left[ 1 + \frac{3}{q^2} - \frac{3}{2q^3} \ln \left| \frac{1+q}{1-q} \right| \right], \quad \Delta\varepsilon''(\omega, q \ge 1) = \frac{\omega_p^2}{\omega^2} \frac{3\pi}{2q^3}, \quad (13) $$

where $q = k/k_0 = v_F k/\omega$ as shown in Fig. 6(b,c). It is evident that a KK-like relation exists between the two parts, and the onset of LD at $q=1$ is a spatial analog to the onset of band edge absorption in temporal domain. One can interpret the spatial dispersion of the real part of the permittivity as the change in permittivity due to “virtual” LD, i.e., virtual indirect transitions in $k$-space where energy and momentum conservation cannot be simultaneously satisfied. This is illustrated by a dashed line in Fig. 6(a) where $q<1$ is not sufficient to achieve LD but still contributing to the change of real part of permittivity.

<figure><img src="image_8.png" /><figcaption>Fig. 6. (a) Real (solid arrow) and virtual (dashed arrow) spatially indirect transitions in metal. (b) Real and (c) imaginary parts of the change of permittivity induced by LD and their overlaps with modal power spectrum.</figcaption></figure>

## Challenges and opportunities

To assess the influence of LD on the properties of plasmonic modes, one can employ a straightforward technique. This involves calculating the power density spectrum of the longitudinal electric field component and subsequently determining the effective change in dielectric constant for the entire mode

$$ \Delta\varepsilon_{\text{eff}} = \frac{\int_{0}^{\infty} \Delta\varepsilon(q) |F_{\parallel}(q)|^2 d^3\mathbf{q}}{\int_{0}^{\infty} |F_{\parallel}(q)|^2 d^3\mathbf{q}}. \quad (14) $$

Figure 6(b,c) illustrate the overlap between the mode and the LD-induced permittivity change. As the mode’s power spectrum extends to larger wavevectors, this overlap increases. For visible light ($\omega \sim 3 \times 10^{15} \text{s}^{-1}$) in Au or Ag ($v_F \sim 1.4 \times 10^6 \text{ m/s}$) the value of LD offset wavevector $k_0 \sim 0.34 \text{ nm}^{-1}$ is significantly larger than the mode bandwidth of $|F_{\parallel}(q)|^2$. Therefore, most of
the overlap occurs at large wave vectors in the tail of the $|F_{\parallel}(q)|^2$ only due to presence of sharp boundary indicating that LD and nonlocality is largely surface phenomena. The impact of LD on the real part of the permittivity is relatively limited due to sign change at $q = 1$ primarily resulting in modifications to the resonant frequency, as extensively documented in the literature. In contrast, the impact of alterations in the imaginary component of the permittivity is substantially more pronounced. The increase in absorptive loss not only reduces field enhancement but also influences field confinement as the field distribution tends to avoid high-loss regions.

In fact, one can evaluate the effective damping rate due to LD as

$$ \gamma_{\text{LD}} = \frac{3\pi\omega \int_{q>1}^{\infty} q^{-3} |F_{\parallel}(q)|^2 d^3\mathbf{q}}{2 \int_{0}^{\infty} |F_{\parallel}(q)|^2 d^3\mathbf{q}} \quad (15) $$

and then simply add it to the bulk damping rate $\gamma_b$ in Drude permittivity

$$ \varepsilon_{\text{eff}}(\omega) = \varepsilon_b - \frac{\omega_p^2}{\omega(\omega + i\gamma_b + i\gamma_{\text{LD}})}. \quad (16) $$

As shown in Ref. [100] one can evaluate effective LD rate as $\gamma_s = \frac{3v_F}{8d_{\text{eff}}}$, where the effective surface to volume ratio is [100–102]:

$$ d_{\text{eff}}^{-1} = \frac{\int \varepsilon_0 E_{\perp}^2(\mathbf{r}) dS}{\int_{\text{metal}} \varepsilon_0 E^2(\mathbf{r}) dV}, \quad (17) $$

where $E_{\perp}^2$ is component of the electric field normal to the surface. This confirms that LD is a surface phenomenon. Its rate, $\gamma_{\text{LD}}$ can be interpreted as the inverse of the time it takes an average electron, moving along the electric field direction, to reach the surface. For plasmon polaritons propagating along the metal-dielectric interface, the result $\gamma_{\text{LD}} = \frac{3}{4}\alpha v_F$ where $\alpha$ is the exponential decay constant inside the metal. For spherical nanoparticles of diameter $a$, the result is $\gamma_{\text{LD}} = \frac{3}{4}v_F/a$ closely matching Kreibig's intuitive interpretation of surface damping [103]. Given that in noble metals bulk damping is on the scale of $\gamma_b \sim 10^{14} \text{s}^{-1}$, it follows that LD becomes the dominant damping process when the characteristic mode size is less than 10 nm.

The impact of LD and changes in imaginary part of permittivity is now clear. Concentrating electric fields in SPP modes increases damping, which alters the dielectric permittivity according to 16, affecting the mode shape and necessitating reevaluation of the damping constant. This iterative process leads to a self-consistent solution, and the achievable field concentration and enhancement are significantly lower than predicted by models neglecting nonlocal effects or considering only first-order real changes in the real part of permittivity (hydrodynamic models [37]). Notably, the mode shape is influenced both inside and outside the metal, with the latter being crucial for applications like nonlinear optics. This is evident in various structures, such as propagating SPPs [104], nanospheres [105], and nano-dimers [106], where field enhancement is often desired outside the metal itself. It's worth noting that even including the diffusion term [64] in the hydrodynamic model of nonlocality only redistributes the field within the metal, leaving the field outside unaffected.

Beyond its impact on field concentration and enhancement, LD absorption also plays a crucial role in hot carrier effects in plasmonics. As absorption occurs near the metal surface and the generated hot carriers are preferentially directed normally to the surface, it makes them more likely to be ejected from the metal. This process underpins plasmon-assisted photodetection and photocatalysis [100,107].

### Future research and concluding remarks

While this summary has highlighted the impact of loss-induced nonlocality on the imaginary part of the permittivity, there remain unexplored aspects of this phenomenon. One intriguing question
is how nonlocality affects absorption in the quantum regime, specifically whether the presence of a single photon at a point $r_1$ can influence the absorption of a second photon at another point $r_2$. Additionally, the “reverse LD”, or luminescence of metallic structures, warrants further investigation, particularly regarding the spatial and angular distribution of emitted photons and polaritons. It is hoped that future research will uncover more phenomena with potential practical applications.

## 5. Landau damping of surface plasmons in metal nanoparticles: the RPA approach

### TIGRAN V. SHAHBAZYAN

#### Overview

There has been a renewed interest in the role of nonlocal phenomena in optical response of metal-dielectric structures [14,64,108]. Landau damping (LD) of localized surface plasmons (LSP) is one of the earliest manifestations of nonlocal effects observed as broadening of the LSP resonance in optical spectra of small metal nanoparticles (NP) [103]. The optically excited LSP decays into single-particle excitations while momentum matching is provided by the electron scattering off the confining potential. For small NPs, this momentum relaxation mechanism can be incorporated, along with the bulk phonon and impurity scattering, into Drude’s dielectric function of the metal $\varepsilon(\omega) = \varepsilon_\infty - \omega_p^2/\omega(\omega + i\gamma)$, where $\omega_p$ is the plasma frequency, and $\gamma$ is the scattering rate. The latter is presented as the sum $\gamma = \gamma_0 + \gamma_s$ of bulk scattering rate $\gamma_0$ and of surface-induced rate

$$ \gamma_s = A \frac{v_F}{L}, \qquad (18) $$

where $v_F$ is the electron Fermi velocity, $L$ is NP’s characteristic size, and $A$ is a phenomenological constant in the range 0.3–1.5 accounting for surface-related effects [103].

#### Current status

The scattering rate $\gamma_s$ was initially associated with electron’s classical scattering (CS) time $\tau_{cs} = L/v_F$ across the NP [109,110], but later it was recognized as intrinsically nonlocal effect. An alternating electric field excites an electron-hole (e-h) pair with energy $\hbar\omega$ across the Fermi level $E_F$. For typical LSP energies $\hbar\omega \ll E_F$, such a process requires momentum transfer $q = \hbar\omega/v_F$, facilitated by surface scattering, which defines nonlocal length scale $\xi = \hbar/q = v_F/\omega$ [64]. For common plasmonic metals such as gold (Au) and silver (Ag), this scale is below 1 nm (e.g., for Au, $\xi \approx 0.5$ nm at wavelength $\lambda = 700$ nm), implying that e-h pair excitation takes place in the close proximity $\xi$ to the metal surface (see Fig. 7). For a NP with characteristic size $L$, the probability of such process occurring during the optical cycle is $\sim \omega\xi/L$, leading to Eq. (18). Despite subnanometer scale of $\xi$, the broadening of optical spectra associated with $\gamma_s$ has been observed for NPs of various shapes with sizes up to $\sim 10$ nm. In the presence of surface scattering, the LSP resonance quality factor is $Q = \omega/\gamma \approx Q_0/[1 + (\xi/L)Q_0]$, where $Q_0 = \omega/\gamma_0$, implying that nonlocal effects in plasmonics persist at much larger scale $L \sim \xi Q_0 \gg \xi$.

For spherical NPs in the size range transitioning from metal clusters to several nm, calculations within jellium model using time-dependent local-density approximation (TDLDA) highlighted the important role of electron confining potential and electron density spillover beyond the NP classical boundary [111,112]. At the same time, for larger NPs in the size range above several nanometers, TDLDA calculations revealed that the precise shape of confining potential is largely unimportant and the overall magnitude of $\gamma_s$, defined by the coefficient $A$, is determined by the electron spillover and dielectric environment effects [113,114]. The reasonable accuracy of Eq. (18) even for relatively large NPs indicates that, for $L \gg \xi$, the LD rate $\gamma_s$ can be obtained as
<figure><img src="image_9.png" /><figcaption>Fig. 7. Schematics for surface-assisted excitation of an e-h pair with energy $\hbar\omega$. (a) An external optical field incident on a metal nanostructure of characteristic size $L$, (b) excites LSP that decays into an e-h pair, (c) with momentum matching provided by carrier surface scattering in a region of size $v_F/\omega$.</figcaption></figure>

a correction to the Drude dielectric function within a standard analytical approach such as the random-phase approximation (RPA) and the Lindhard dielectric function [9].

RPA calculations of $\gamma_s$ have been carried out in the course of several decades [100–102, 115–119]. While in earlier studies, the scattering rate for spherical NPs with radius $a$ was derived as $\gamma_{sp} = 3v_F/4a$, the subsequent works focused on obtaining $\gamma_s$ for NPs of more complex shapes. The major challenge was to obtain the LD rate for NPs of arbitrary shape, even irregular one, for which no numerical calculations are feasible. The condition $L \gg \xi$ implies that the length scales for electronic and LSP excitations are well separated. For hard-wall confining potential, the electronic contribution to $\gamma_s$ can be effectively “integrated out” and the LD rate is obtained as a nonlocal correction to the Drude decay rate as $\gamma = \gamma_0 + \gamma_s$, where $\gamma_s$ has the form 18 with the characteristic size $L$ depending on the field distribution in the metal [100–102]:

$$L = \frac{\int dV |\mathbf{E}|^2}{\int dS |E_n|^2} \qquad (19)$$

Here, $\mathbf{E}$ is the electric field inside a NP of volume $V$, $E_n$ is the field component normal to the NP surface $S$, and $A = 3/4$. Importantly, $\gamma_s$ is sensitive to polarization of the LSP field that drives the electrons towards the interface. Such form of $L$ (and $\gamma_s$) is valid for NPs in the size range $\xi \ll L \ll \lambda$, where $\lambda$ is the optical wavelength, which includes most plasmonic systems used in the applications.

RPA approach to Landau damping in small metal nanoparticles. Here we briefly outline the RPA approach to LD of LSP in metal NPs following Refs. [101,102]. An external field excites an LSP in a NP which subsequently decays into an e-h pair by promoting, with its alternating field $Ee^{-i\omega t}$, a conduction band electron across the Fermi level, while the momentum matching is provided by carriers’ surface scattering (see Fig. 7). The full dissipated power $Q$ in a metal NP is given by

$$Q = \frac{\omega}{2} \mathrm{Im} \int dV \mathbf{E}^* \cdot \mathbf{P}, \qquad (20)$$

where $\mathbf{P}(\mathbf{r})$ is the electric polarization vector (the star stands for complex conjugation). The bulk contribution to $Q$ is obtained by relating, in the local limit, the polarization vector to the electric field as $\mathbf{P}(\mathbf{r}) = \mathbf{E}(\mathbf{r})[\varepsilon(\omega) - 1]/4\pi$, yielding the standard expression

$$Q = \frac{\omega \varepsilon''(\omega)}{8\pi} \int dV |\mathbf{E}|^2, \qquad (21)$$

where $\varepsilon''(\omega)$ is the imaginary part of metal dielectric function due to bulk relaxation processes. For small NPs, there is also a surface contribution $Q_s$ to the dissipated power arising from momentum relaxation due to surface scattering. The general expression for $Q_s$ is obtained by
relating $P(r)$ to nonlocal electron polarization operator $\Pi(\omega; r, r')$ via the induced charge density $\rho(r)$ as

$$ \nabla \cdot P(r) = -\rho(r) = -e \int dV' \Pi(\omega; r, r') \Phi(r'). \quad (22) $$

Here, the potential $\Phi(r)$ is defined as $eE(r) = -\nabla\Phi(r)$, where $e$ is the electron charge. Using Eq. (22), after integrating Eq. (20) by parts, the dissipated power takes the form

$$ Q_s = \frac{\omega}{2} \operatorname{Im} \int dV dV' \Phi^*(r) \Pi(\omega; r, r') \Phi(r'), \quad (23) $$

where $\Pi(\omega; r, r')$ includes only the electronic contribution. Within RPA, $\Pi(\omega; r, r')$ is replaced by the polarization operator for noninteracting electrons, yielding

$$ Q_s = \pi\omega \sum_{\alpha\beta} |M_{\alpha\beta}|^2 [f(\epsilon_{\alpha}) - f(\epsilon_{\beta})] \delta(\epsilon_{\alpha} - \epsilon_{\beta} + \hbar\omega), \quad (24) $$

where $M_{\alpha\beta} = \int dV \psi_{\alpha}^* \Phi \psi_{\beta}$ is the transition matrix element of potential $\Phi(r)$ calculated from the wave functions $\psi_{\alpha}(r)$ and $\psi_{\beta}(r)$ of electron states with energies $\epsilon_{\alpha}$ and $\epsilon_{\beta}$ separated by $\hbar\omega$, $f(\epsilon)$ is the Fermi distribution function, and spin degeneracy is included. For NPs of arbitrary shape, a direct numerical evaluation of $M_{\alpha\beta}$ is not possible due to complexity of the electron wave functions. However, for NPs with characteristic size $L \gg \xi$, this issue can be bypassed by extracting the surface contribution to the matrix element as [101]

$$ M_{\alpha\beta}^{s} = \frac{-e\hbar^4}{2m^2\epsilon_{\alpha\beta}^2} \int dS [\nabla_n\psi_{\alpha}(s)]^* E_n(s) \nabla_n\psi_{\beta}(s), \quad (25) $$

where $\nabla_n\psi_{\alpha}(s)$ is wave function's derivative normal to the surface, $E_n(s)$ is the corresponding normal field component, $\epsilon_{\alpha\beta} = \epsilon_{\alpha} - \epsilon_{\beta}$ is the $e-h$ pair excitation energy, and $m$ is the electron mass. The integration in Eq. (25) takes place over the NP surface $S$, while the volume contribution to the matrix element is negligibly small due to near-vanishing overlap of the electron wave-functions in the presence of slowly varying potential $\Phi$. Using this expression for the matrix element, the surface contribution to the dissipated power (24) can be recast as

$$ Q_s = \int \int dS dS' E_n(s) E_{n'}^*(s') F(s, s'), \quad (26) $$

where $F(s, s')$ is the $e-h$ correlation function, which is expressed via normal derivatives of the electron and hole Green functions in a hard-wall cavity (see Ref. [101]).

Evaluation of $Q_s$ hinges on the observation that excitation of an $e-h$ pair with energy $\hbar\omega$ takes place in a region of size $\sim \xi \ll L$. Then it can be shown that $F(s, s')$ peaks in the region $|s - s'| \sim \xi$ and rapidly oscillates outside of it. Since the electric field is relatively smooth on the scale $L$, the $e-h$ correlation function can be approximated by $F(s, s') = F_0 \delta(s - s')$, where the coefficient $F_0$ is evaluated using single-scattering approximation for the electron Green function as $F_0 = (3v_F\omega_p^2/32\pi\omega^2)$ [101]. The final expression for the surface-induced dissipated power has the form

$$ Q_s = \frac{3v_F \omega_p^2}{32\pi \omega^2} \int dS |E_n|^2. \quad (27) $$

Comparing $Q_s$ with the bulk expression (21), one observes that both contributions can be combined together by adding an imaginary nonlocal correction $\delta\epsilon_s$, where

$$ \delta\epsilon_s = i \frac{\omega_p^2 \gamma_s}{\omega^3}, \quad \gamma_s = \frac{3v_F}{4} \frac{\int dS |E_n|^2}{\int dV |E|^2}, \quad (28) $$

to the Drude bulk dielectric function. The above form for $\delta\epsilon_s$ implies that the scattering rate in the Drude dielectric function $\epsilon(\omega)$ should be modified as $\gamma = \gamma_0 + \gamma_s$. A similar expression for $\gamma_s$ was obtained in Ref. [100] using a different approach.
Although $\gamma_s$ is independent of the electric field's overall amplitude, it is highly sensitive to field's polarization that defines the electrons motion relative to the metal-dielectric interface. In the case of asymmetric NPs, such polarization dependence implies that the LD rate can vary significantly for different LSP modes excited in the same system. This point is illustrated in Fig. 8 for longitudinal (L) and transverse (T) modes in nanorods and nanodisks, modeled by the prolate (P) and oblate (O) spheroids, respectively. For such systems, explicit analytic expressions for the LD rate has been obtained in the form [101] $\gamma_s = \gamma_{sp}f_{L,T}$, where

$$f_L = \frac{3}{2 \tan^2(\alpha)} \left[ \frac{2\alpha}{\sin(2\alpha)} - 1 \right], \quad f_T = \frac{3}{4 \sin^2(\alpha)} \left[ 1 - \frac{2\alpha}{\tan(2\alpha)} \right]. \qquad (29)$$

Here, for a prolate spheroid (nanorod), $\alpha = \arccos(b/a)$, where a and b are semi-major and semi-minor axis respectively, while for an oblate spheroid (nanodisk), the LD rates have the same form 29 but with $\alpha = \operatorname{arccosh}(b/a)$. For comparison, the CS model rate $\gamma_{cs} = v_F S / 4V$, which is independent of mode polarization, is also plotted in Fig. 8. For visual convenience, all rates are normalized by the LD rate $\gamma_{sp} = 3v_F / 4a$ for a spherical NP of radius a. At the sphere point a = b, the normalized rates continuously transition into each other (e.g., PL to OL and PT to OT), but away from it, the rates for different modes exhibit dramatic difference in magnitude depending on the mode polarization.

## Discussion and further developments

The surface scattering rate Eq. (28) quantifies LD in “simple” NPs characterized by a single metal surface. Further developments involved more complex plasmonic systems such as thin films [120], core-shell hybrid structures [121–123] and NP dimers [106,124]. In thin films or metal nanoshells with dielectric core, the surface scattering that accompanies e-h pair excitation can take place from both the inner and outer metal boundaries. The interference between these processes can lead to coherent oscillations (quantum beats) of $\gamma_s$ with changing metal thickness d. Such oscillations were reported in TDLDA studies of thin Ag films [120] and later in RPA calculations of $\gamma_s$ for spherical metal nanoshells [121,122]. Specifically, for metal nanoshells with thickness d, the LD rate acquires an interference-induced contribution $\gamma_s^{int} \approx B \sin(d/\xi)$ where the coefficient B depends on the nanoshell's thickness and overall size [122].

<figure><img src="image_10.png" /><figcaption>**Fig. 8.** Normalized surface scattering rates for prolate and oblate spheroids along with the CS rate are plotted against aspect ratio $b/a$. Insets: Schematics of LSP modes.</figcaption></figure>
In plasmonic dimers, LD has been shown to play a critical role in limiting the field enhancement in the gap between closely spaced NPs [106,124]. Numerical calculations of the field intensity performed within RPA using the Lindhard dielectric function [106] and within the generalized nonlocal optical response (GNOR) model [124] revealed a significant increase in the LD rate that reduced the field enhancement factor in the gap by nearly two orders of magnitude. Qualitatively, the LD rate increase for narrow gaps can be inferred from the characteristic size *L*, given by Eq. (19), which depends on the field distribution inside the metal. Note that, in the absence of electron wave-function overlap between two NPs, excitation of an *e*-h pair can take place in either NP, so that the LD rate for a dimer is the sum of individual NP LD rates. In dimers, the field intensity is highest in the gap between the NPs and is much weaker in the metal, where it is concentrated in a relatively small volume around the gap. More precisely, the integrated LSP density of states (DOS) $\rho(\omega)$ can be split into contributions from the regions inside the metal and outside of it as $\rho(\omega) = \rho_{\text{in}}(\omega) + \rho_{\text{out}}(\omega)$, with their relative magnitude $\rho_{\text{out}}(\omega)/\rho_{\text{in}}(\omega) = |\varepsilon'(\omega)| \gg 1$, where $\varepsilon'(\omega)$ is the real part of metal dielectric function [125]. With decreasing NP separation, as the LSP resonance redshifts [106,124], this ratio increases with $|\varepsilon'(\omega)|$ implying that the electric field is further pushed out from the metal into the gap. As a result, the volume-integrated field intensity in the numerator of Eq. (19) decreases relative to the surface-integrated intensity in the denominator, leading to a reduced *L* and, hence, enhanced LD rate for narrow gaps.

The overall magnitude of the LD rate 18 is defined by the coefficient *A* that depends on the electron confining potential and charge density profile, which, in turn, determine the electric fields near the interface. Recent TDLDA calculations for relatively large (up to 10 nm) NPs [113,114] indicate that the main impact on the LSP resonance width comes from the electron density spillover and dielectric environment, rather than the electronic states in the cavity. Since the electron spillover is not sensitive to the overall NP shape, the TDLDA value *A* ≈ 0.32 for a nanosphere [113,114] is likely shape independent. Note, however, that substantially larger values in the range 0.3–1.5 were reported in the experiment depending on the surrounding dielectric medium [103]. On the other hand, the presence of d-band electrons with a nearly step-like density profile in noble-metals gives rise to a thin surface layer, in which the conduction electrons with extended density tail are no longer screened by the d-band electrons, which leads to the field enhancement near the interface [126]. In Ag NPs, this effect has been shown to cause a blueshift of the LSP resonance which, in fact, overcompensates the resonance redshift due to the electron density spillover [127]. One can expect that a similar competition between these two nonlocal mechanisms will take place for the LD rate as well and bring the coefficient *A* closer to the experiment. A related effect has been recently studied for Ag NPs coated with a thin dielectric shell [123]. In this case, the electron spillover into a medium with a much weaker field screening results in a noticeable LD rate enhancement.

## Concluding remarks

In this contribution, we tried to present a brief outlook on the Landau damping (LD) of surface plasmons in metal nanoparticles (NPs). This phenomenon has attracted a constant interest for over 50 years due to its prominent role in the optical spectra of small NPs. There is a very extensive literature on LD in NPs, and several theoretical and numerical approaches have been developed during this time span. We have focused on recent advances within the RPA approach which resolved a long-standing problem of describing, by means of a relatively simple analytical model, the LD of surface plasmons in metal NPs of arbitrary shape. On the applications side, LD is an efficient source of hot electrons widely used in photochemistry and light harvesting; there are excellent reviews on this topic the readers are referred to.

**Funding.** National Science Foundation (DMR-2000170, DMR-2301350, NSFPREM-2423854).
6. Madelung hydrodynamics, nonlocal plasmonics, and nonlinear optics

ANDRÉ J. CHAVES ID, N. ASGER MORTENSEN ID & NUNO M. R. PERES ID

Overview

Quantum–hydrodynamic models provide a robust framework for exploring the electrodynamic behavior of both three-dimensional (3D) and two-dimensional (2D) electron gases, particularly in the fields of nonlocal plasmonics and nonlinear optics. The Madelung equations, when coupled with Poisson’s equation, allow for the analysis of magnetoplasmon spectra, magneto-optical conductivity, and nonlocal Fermi pressure corrections, and are generally relevant for nonlinear and nonlocal light-matter interactions. Applications include plasmonics [55,61,128–130], nonlinear optics [131–134], transistor design [135], the analysis of two-dimensional materials [136–139], and semiconductor optics [140].

The model originates from early developments in quantum mechanics, incorporating both the statistical interpretation of the wave function and the hydrodynamic formulation. In 1926, Born introduced the statistical interpretation of the wave function [141–143]. Shortly thereafter, Madelung reformulated the single-particle time-dependent Schrödinger equation in hydrodynamic form [6,144], linking quantum mechanics to some kind of fluid dynamics.

The resulting Madelung equations [6] resemble the continuity and Euler equations in fluid dynamics but initially lacked mechanisms to account for many-body effects, such as statistical pressure [40]. Over time, the framework was expanded to include these effects, such as Fermi pressure and exchange-correlation terms [145]. The relationship between the hydrodynamic model and the Schrödinger equation is bidirectional; starting from the hydrodynamic formulation, an effective Schrödinger equation can also be derived [146–148]. This approach facilitates the study of complex quantum systems using established numerical methods.

When studying electronic fluid dynamics under electric and magnetic fields [149–153], the self-consistent determination of these fields becomes necessary, often requiring that the hydrodynamic model to be solved alongside Maxwell’s equations. In particularly within plasmonics [154], the non-retarded approximation [155] simplifies these calculations by replacing Maxwell’s equations with Poisson’s equation, applicable under specific conditions.

The original Madelung equations, however, do not account for many-body quantum phenomena [156]. This limitation can be addressed by either reformulating the approach using a many-body wave function [156,157] or by deriving the equations from the moments of the Wigner distribution function [146]. The latter approach naturally incorporates elements such as statistical pressure, the Bohm term, and additional contributions from the energy potential [146], thus extending the formalism to include many-body effects.

**Madelung hydrodynamics.** The time-dependent Schrödinger equation is given by

$$i \frac{\partial \Psi(\mathbf{r}, t)}{\partial t} = H \Psi(\mathbf{r}, t), \qquad (30)$$

where $H$ is the system’s Hamiltonian describing a particle of mass $m$ subjected to a potential $U$, and $\Psi(\mathbf{r}, t)$ is the complex-valued wave function. By expressing the wave function as $\Psi(\mathbf{r}, t) = \sqrt{n(\mathbf{r}, t)} \exp[iS(\mathbf{r}, t)]$, where $n(\mathbf{r}, t)$ and $S(\mathbf{r}, t)$ are real-valued functions, and substituting into Eq. (30), two equations emerge:

$$\frac{\partial n}{\partial t} + \nabla \cdot (n\mathbf{v}) = 0, \qquad (31a)$$

$$\frac{\partial \mathbf{v}}{\partial t} + \frac{1}{2} \nabla \mathbf{v}^2 = -\frac{1}{m} \nabla U + \frac{\hbar^2}{2m^2} \nabla \left( \frac{1}{\sqrt{n}} \nabla^2 \sqrt{n} \right). \qquad (31b)$$
Here, $v = \hbar \nabla S/m$ is the velocity field. Equation (31a) resembles the continuity equation, while (31b) has the form of the Euler equation for fluid motion. The term $-m^{-1} \nabla U$ accounts for external forces, and the final term, known as the quantum (or Bohm) potential, arises purely from quantum mechanics, vanishing as $\hbar \to 0$. This quantum potential introduces spatial dispersion, leading to nonlocal effects in systems like plasmons. When $U(r, t)$ is the self-consistent electrostatic potential $V_{sc}(r, t)$, the system is supplemented by Poisson's equation:

$$ \nabla^2 V_{sc}(r, t) = -\frac{q}{\epsilon_0} n(r, t), \qquad (32) $$

where $q$ is the particle charge and $\epsilon_0$ is the vacuum permittivity. In cases where retardation is significant, Maxwell's equations replace Poisson's equation.

**Many-body effects.** The derived Euler equation excludes many-body effects because it originates from the single-particle Schrödinger equation. Many-body effects can be incorporated by introducing a phenomenological term or by deriving equations from a many-body wave function.

In many-body systems, the Euler equation acquires an additional term, the Fermi force, given by

$$ \mathcal{F}_F = -\frac{\nabla p_F[n(r, t)]}{mn}, \qquad (33) $$

where $p_F$ is the Fermi pressure, which depends on the system's dimensionality and particle dispersion. For specific cases, it is given by

$$ p_F = \begin{cases} \frac{\hbar^2 n}{5m} (3\pi^2 n)^{2/3}, & (3D electron gas) \\ \pi \frac{\hbar^2 n^2}{2m}, & (2D electron gas) \\ \frac{\hbar v_F}{3\pi} (\pi n)^{3/2}, & (graphene) \end{cases} \qquad (34) $$

where $v_F$ is the Fermi velocity for graphene. Both the Bohm potential and the Fermi pressure contribute to the k-dependence of the dielectric function, differing in the powers of k, i.e., $k^2$ for Fermi pressure and $k^4$ for the Bohm potential.

**Wigner function approach.** The Madelung equations can also be derived from the Wigner distribution function, which preserves information about both position and momentum while maintaining the interference properties characteristic of quantum mechanics. The Wigner distribution approach assumes a distribution function of the following form (in 2D)

$$ f(\mathbf{r}, \mathbf{p}, t) = \frac{1}{(2\pi)^2 \hbar^2} \sum_{\alpha} P_{\alpha} \int d\mathbf{r}' \psi_{\alpha}^*(\mathbf{r} + \mathbf{r}'/2, t) \psi_{\alpha}(\mathbf{r} - \mathbf{r}'/2, t) e^{i\mathbf{p} \cdot \mathbf{r}' / \hbar}, \qquad (35) $$

where $\psi_\alpha(\mathbf{r}, t)$ is a single-particle wave function obeying the time-dependent Schrödinger Eq. (30). The system is assumed to be in a mixed state with a density matrix $\rho$ of the form $\rho(\mathbf{r}, \mathbf{r}', t) = \sum_\alpha P_\alpha \psi_\alpha^*(\mathbf{r}', t) \psi_\alpha(\mathbf{r}, t)$, where the probabilities $P_\alpha$ obey the sum rule $\sum_\alpha P_\alpha = 1$, where $N_s$ is the number of states in the mixture. The time derivative of $f(\mathbf{r}, \mathbf{p}, t)$ can be obtained using the time-dependent Schrödinger Eq. (30). Introducing the following two moments of the Wigner distribution

$$ n(r, t) = \int dp f(r, p, t), \quad n(r, t)v(r, t) = \frac{1}{m} \int dp p f(r, p, t), \qquad (36) $$

the continuity equation arises from the first moment, while Euler's equation arises from the second moment, with the latter incorporating both statistical and quantum (Bohm term) pressures.
Current status

**Magnetoplasmon spectrum and optical conductivity.** The hydrodynamic formulation of quantum mechanics has been effectively applied to problems such as second-harmonic generation in 2D and 3D electron gases and the calculation of plasmon dispersion, both in the presence and absence of a magnetic field (magnetoplasmons). Including a magnetic field in the Madelung's equation amounts to changing the right-hand side of Euler's equation by including the term

$$-\frac{q}{m} \nabla V - \frac{q}{m} \frac{\partial A}{\partial t} + \frac{q}{m} \mathbf{v} \times (\nabla \times \mathbf{A}), \qquad (37)$$

where V and A are the scalar and vector potential, respectively. This does not include the Fermi pressure, which can be incorporated in an ad hoc manner. Additionally, it is assumed that a weak magnetic field allows for the use of the previously provided expressions for the Fermi pressure. When an external electromagnetic field is introduced, the magnetic field component must be accounted for to ensure a consistent calculation of the system's nonlinear optical properties. With all terms included, Euler's equation reads

$$\frac{\partial \mathbf{v}}{\partial t} + (\mathbf{v} \cdot \nabla)\mathbf{v} = -\frac{\pi \hbar}{m^2} \nabla n - \frac{q}{m} \nabla V_{\text{sc}} + \frac{q}{m} \mathbf{v} \times \mathbf{B}, \qquad (38)$$

where $V_{sc}$ is the self-consistent potential. By linearizing the problem, the calculation of the magnetoplasmon spectrum arises from the simultaneous solution of the continuity, Euler's, and Poisson's equations, giving [158]

$$\omega = \Omega_k = \sqrt{\omega_c^2 + ak + \beta^2 k^2}. \qquad (39)$$

Here, $a = q^2 n_0 / (2m\epsilon_0)$, $\omega_c = |qB|/m$ is the cyclotron frequency, and $\beta^2 = v_F^2 / 2$, where $\beta$ corresponds to the speed of the first sound in a 2D electron gas [159]. For the 3D electron gas we have $\beta^2 = v_F^2 / 3$.

Applying an external electromagnetic field rather than a self-consistent electrostatic potential induces a charge current $J_1 = qn_0\mathbf{v}_1$, which leads to the optical conductivity tensor

$$\overleftrightarrow{\sigma} = -\frac{q^2 n_0}{D(k, \omega)} \begin{bmatrix} im\omega & -qB + ik_x k_y \beta^2 \\ qB + ik_x k_y \beta^2 & im\omega \end{bmatrix} \qquad (40)$$

where $D(k, \omega) = q^2 B^2 + m^2 \beta^2 k^2 - m^2 \omega^2$. When nonlocal corrections are neglected [160], this result agrees with calculations based on Boltzmann's kinetic equation (see also Refs. [161,162]). From the off-diagonal term of the conductivity sensor, we see that, due to nonlocal effects, there exists a transverse response even in the absence of a magnetic field.

**Nonlinear effects.** We can extend the formalism for the calculation of nonlinear optical properties of the electron gas [163–165]. Up to second order, the nonlinear current is defined as $J_{\text{NL}} = qn_0\mathbf{v}_2 + qn_1\mathbf{v}_1$, where $\mathbf{v}_2$ is the new term resulting from the expansion of the velocity field up to second order. Proceeding as in linear response, and assuming a zero static magnetic field, we find that

$$J_{\text{NL},x} = \frac{3}{2} k_x \frac{e^3 v_F^2}{4\pi \hbar^2 \omega^3} E_{\omega,x}^2 e^{-i2\omega t} - \frac{3}{2} k_x \frac{e^3 v_F^2}{4\pi \hbar^2 \omega^3} E_{\omega,x}^* E_{2\omega,x} e^{-i\omega t} + \text{c.c.} \qquad (41)$$

This result agrees with the literature [166,167] and where we have used that $q = -e < 0$, with e being the elementary charge. We emphasize that $E_\omega(r)$ corresponds to the total field and not to the external field. Note that we have assumed that the total electric field has a finite in-plane
momentum component $k_x$. For normal incidence, $k_x = 0$ and the second order response vanishes, as it should according to inversion-symmetry arguments.

**Plasmon-enhanced second-harmonic generation.** The calculation of second-harmonic generation assisted by plasmons [168] can be performed by including a potential of the form $V_{\text{ext}} + V_{\text{ind}}$ in Euler's equation, where $V_{\text{ext}}$ corresponds to the external potential and $V_{\text{ind}}$ represents the induced potential from plasmon excitation, which is accounted for via Poisson's equation.

Proceeding in the same manner as for the determination of the magnetoplasmon spectrum, we find, including nonlocal corrections, that [169]

$$ \phi_{2\omega}^{\text{ind}} = -\frac{e^3 k^3 n_0}{2\epsilon_0 m^2} \frac{3\omega^2 - \gamma^2 + \beta^2 k^2}{(\beta^2 k^2 - \omega(\omega + i\gamma))(4\beta^2 k^2 - 2\omega(2\omega + i\gamma))} \times \frac{1}{(\omega + i\gamma)^2} \frac{1}{\hat{\epsilon}(2\mathbf{k}, 2\omega)} \frac{1}{[\hat{\epsilon}(\mathbf{k}, \omega)]^2} (\phi_{\omega}^{\text{ext}})^2, \quad (42) $$

where $\phi_{\omega}^{\text{ext}}$ is the external field. The denominators of $\phi_{2\omega}^{\text{ind}}$ can be easily expressed in terms of the dielectric function with damping and reads as follows:

$$ \hat{\epsilon}(\mathbf{k}, \omega) = 1 - \frac{\omega_k^2}{\omega(\omega + i\gamma) - \beta^2 k^2}. \quad (43) $$

The longitudinal excitations of the gas are determined from $\hat{\epsilon}(\mathbf{k}, \omega) = 0$, implying the spectrum $\omega(k) = \sqrt{e^2 n_0 k / (2\epsilon_0 m) + \beta^2 k^2}$. Therefore, the SHG is enhanced by plasmon excitation.

As opposed to the case studied previously for zero magnetic field, where nonlocal corrections were ignored, the finite value of $\phi_{2\omega}^{\text{ind}}$ requires a finite in-plane wave vector.

In Fig. 9 we depict the function $\phi_{2\omega}^{\text{ind}}$ as function of the plasmons wave number and energy $\hbar\omega$. The two harmonics are clearly visible, whose dispersion follows the plasmon curve of the 2D electron gas.

<figure><img src="image_11.png" /><figcaption>**Fig. 9.** Second-harmonic generation assisted by plasmons. The two bright lines correspond to the fundamental and second harmonic. The parameters are $n_0 = 10^{12} \text{cm}^{-2}$ and $\hbar\gamma = 0.25 \text{ meV}$, and using for $m$ the free electron mass.</figcaption></figure>

## Challenges and opportunities

The primary advantage of the hydrodynamic model lies in its simplicity compared to other methods that incorporate quantum mechanical effects. This simplicity proves advantageous in
two distinct contexts: i) deriving straightforward analytical results for relevant problems and
ii) numerically solving complex systems where the hydrodynamic model effectively captures
quantum effects. Below, we briefly discuss how these features can motivate the application of the
hydrodynamic model to recent and intriguing problems.

First, the connection between the hydrodynamic and Schrödinger equations, as demonstrated
in the previous section, enables their parallel use alongside Maxwell's or Poisson's equations.
This approach has been previously explored to study electron spill-out [148], but could be further
applied to various plasmonic problems using numerical methods from quantum mechanics.

The Wigner transformation has been applied in the context of nanoscale electronic transport [170], with the possibility of incorporating electron-phonon [171] and electron-impurity interactions [172]. Similar approaches could be employed to incorporate the effects of phonons and impurities in the hydrodynamic model, aiding in the understanding of their role in plasmon propagation, spectral broadening, and plasmon-phonon hybridization. Furthermore, the hydrodynamic equation can account for viscosity [173,174].

As nanoparticle size decreases, light-matter interactions deviate from their bulk counterparts
due to quantum confinement [175]. The hydrodynamic model can incorporate quantum effects,
such as the Fermi pressure term, while remaining compatible with most numerical approaches in
nano-optics, such as the finite-element method [55]. Thus, the hydrodynamical equation can be
used to investigate the plasmonic properties of nanostructures, including novel 2D van der Waals
materials.

Although this study focused on second-harmonic generation in graphene, the procedure can be extended to higher-order nonlinear responses, such as third-harmonic generation [91], four-wave mixing [176], plasmon self-modulation [177], and solitons [178]. In the case of self-modulation phenomena driven by nonlinearity, we expect a renormalization of the plasmon dispersion and an amplification of nonlocal effects due to nonlinearity. Using the hydrodynamic formalism, an analytic expression for the renormalized spectra can, in principle, be derived. Furthermore, rather than analyzing the Fourier decomposition of the fields for each order, one could first solve the hydrodynamic equations in the time domain and subsequently decompose the fields in the frequency domain.

This work considered the Schrödinger Hamiltonian corresponding to a Fermi liquid described
by an effective mass Hamiltonian. However, in the realm of quantum materials, ranging from
graphene to topological insulators [179], one can expect phenomena such as flat bands, strongly
correlated systems, and nontrivial Bloch curvature. Other material excitations, like excitons,
phonons, and Cooper pairs, strongly couple with light [180]. A relevant question, then, is whether
a hydrodynamic model akin to the one presented here can be developed for these excitations.

Future developments to address challenges

Generalizations of the quantum hydrodynamic model to address light-matter interactions in anisotropic systems, such as phosphorene, and strongly correlated systems, such as magnetic systems, are urgently needed. Finally, a formulation based on the Wigner distribution is a possible development to address the mentioned challenges.

Developing models that go beyond the hydrodynamic model while retaining similar simplicity
is a challenging yet rewarding endeavor. One such relevant problem is the study of plasmons in
magic-angle twisted bilayer graphene, which exhibits flat bands [181,182].

Concluding remarks

The quantum-hydrodynamic framework, derived from the Madelung equations, has evolved
into a versatile tool for analyzing nonlocal and nonlinear effects in photonic materials. Its
capability to incorporate many-body interactions, magnetic fields, and self-consistent electrostatic
potentials enables rich insights into electron dynamics. By bridging quantum mechanics and hydrodynamics, this approach paves the way for advancing applications of advanced materials.

**Funding.** Fundação para a Ciência e a Tecnologia (UIDB/04650/2020, COMPETE 2020, PORTUGAL 2020, FEDER, PTDC/FIS-MAC/2045/2021); Danmarks Grundforskningsfond (DNRF165); Conselho Nacional de Desenvolvimento Científico e Tecnológico (423423/2021-5, 408144/2022-0, 315408/2021-9); Fundação de Amparo à Pesquisa do Estado de São Paulo (2022/08086-0); Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (PrInt scholarship).

## Part II—Nonlocal effects in free-electron metals and plasmonic nanostructures

### 7. Continuum framework for plasmonic systems

GINO WEGNER & KURT BUSCH

#### Current status

The continuum framework for plasmonic systems deals with the underlying many-body dynamics while, at the same time, embracing only a small set of degrees of freedom, thereby providing access to a large class of analytically and numerically tractable observables. This approach is particularly suited to handle the collective conduction electron dynamics, i.e. the protagonists of plasmonic excitations. In fact, besides bare plasmons, these plasmon excitations also include hybrid light-plasmon excitations commonly referred to as plasmon polaritons. While plasmons occur both, in the bulk and near surfaces, the finite penetration depth of electromagnetic fields in metals confines plasmon polaritons to surfaces (surface plasmon polaritons, SPPs) or small nano-particles (localized surface plasmon polaritons, LSPs). Over the past few decades, these plasmonic excitations have been studied extensively both, theoretically and experimentally, and find numerous applications in nano-photonic devices and plasmon-driven chemistry.

The continuum framework rests on local averaging schemes (in the multi-scale modeling also known as coarse-graining), which yield balance equations for the bulk of the form

$$ \partial_t \mathcal{F}^{(\alpha)}(\mathbf{r}, t) = -\nabla \cdot \mathcal{F}^{(\alpha+1)}(\mathbf{r}, t) + \mathcal{S}^{(\alpha)}(\mathbf{r}, t), \quad (44) $$

thereby connecting consecutive-order moments $\mathcal{F}^{(\alpha)}$ of increasing tensorial rank while introducing different source/sink terms $\mathcal{S}^{(\alpha)}$. Here, a position $\mathbf{r}$ specifies a volume element $dV$ which is small compared to the whole continuum, but large compared to the typical distance of constituents [183]. Contemporary bulk models in plasmonics describe the dynamics of the charge density $\rho$, which, in the absence of external electromagnetic fields, is assumed to be spatially homogeneous. In turn, the charge density is coupled to the current density $\mathbf{J}$ by the continuity equation [11]

$$ \partial_t \rho(\mathbf{r}, t) = -\nabla \cdot \mathbf{J}(\mathbf{r}, t). \quad (45) $$

Similarly, the current balance reads

$$ \partial_t \mathbf{J}(\mathbf{r}, t) = -\nabla \cdot \mathbf{\Pi}(\mathbf{r}, t) + \mathcal{S}^{(1)}(\mathbf{r}, t), \quad (46) $$

and is tied to momentum balance by the ubiquitous relation $\mathbf{J} = (\rho/m)\mathbf{p}$ where $\mathbf{p}$ denotes the momentum.

Different continuum models differ in the form of $\mathbf{\Pi}$ which, in linear response, describes the stress tensor [183] and is typically split into pressure- and shear-like contributions, i.e. $\mathbf{\Pi} = \mathbf{\Pi}_p + \mathbf{\Pi}_{sh}$. The stress tensor describes locally the macroscopic response due to short-ranged (classical and quantum) interactions induced by deformations of the continuum [184] [see Fig. 10(a,b)]. In practice, it is often expressed by equations of state [185] via $\rho$ and/or $\mathbf{J}$-introducing nonlocality which originates in the conduction electron's Fermi-Dirac statistics and scales with the Fermi velocity $v_F$. In the Drude model [42], the charge density is assumed to be incompressible and
constant $\Pi \equiv 0$ and the current sink/source is given by a phenomenological damping term via a relaxation rate $\gamma$ and an external field, which adds to the internal Coulomb field. The ionic background is treated as a rigid continuum for overall charge-neutrality and a restoring force [24]. As a general feature, the Drude model allows for the existence of plasmons, SPPs, LSPs and non-radiative dissipation. Within the Euler-Drude model [45] (i.e., the linearized version of the HDM discussed in Sec. 3 of this Roadmap) attributed to Bloch [7], the Drude model is augmented by a hydrostatic pressure term, such that $\nabla \cdot \Pi_p = \beta_{LF}^2 \nabla \rho$, where $\beta_{LF} = v_F/\sqrt{3}$ is the characteristic velocity of compressional waves. When combined with hard-wall boundary conditions, this model is known to qualitatively reproduce the size-dependent nonlocal blueshift of LSPs in noble metals. In fact, the hard-wall boundary condition prohibits the spill-out of the charge density which is implied by the large values of these systems’ work function [16] (see also Sec. 3 of this Roadmap). Further, the viscoelastic model [173] introduces shear. It interpolates between the viscous shear at small frequencies of the incident field compared to the rate $\gamma$ and non-dissipative elastic shear [184] at large frequencies. The dissipative nature of the former limit [183,186] [corresponding to a Navier-Stokes(-Drude) model] is tied to a damping of the internal stresses on the time scale $\sim 1/\gamma$. The viscoelastic shear stress is reminiscent of highly viscous fluids and evolves according to [184,187]

$$[\partial_t + \gamma] \nabla \cdot \Pi_{\equiv \text{sh}}(\mathbf{r}, t) = -\frac{4}{3} \beta_{\text{el}}^2 \nabla [\nabla \cdot \mathbf{J}(\mathbf{r}, t)] + \beta_{\text{el}}^2 \nabla \times \nabla \times \mathbf{J}(\mathbf{r}, t), \quad (47)$$

where $\beta_{el} = v_F/\sqrt{5}$ denotes the characteristic velocity of transverse plasma waves due to elastic shear [188]. The double-curl-operator introduces transverse nonlocality, while the other term, together with the pressure embodies longitudinal nonlocality.

<figure><img src="image_12.png" /><figcaption>Fig. 10. Illustration of a continuum framework in plasmonics. (a) Connection of the micro- and macroscopic scales of the conduction electron (left) and the ionic background (right) continuum, where $\ell_{dV}$ and $L$ denote the extent of a continuum parcell and a linear dimension of the respective continuum, respectively. The microscopic scales $n_0^{-1/3}$ and $d$ correspond to the mean electronic distance and lattice constant, respectively. (b) Continuum deformations that are considered. (c) Selected bulk material models, their (hierarchical) interrelation in Fourier-space (with wavevector $k$ and frequency $\omega$), correspondingly added bulk dynamics and associated impact on LSPs. $\mathcal{L}$ denotes the linear projection operator, while $O$ refers to an arbitrary observable used to study the spectral properties of LSPs. Blue (red) lines denote nonlocal (local) realizations of the observable. Nonlocal models are combined with hard-wall conditions throughout.</figcaption></figure>

As a last model, we introduce a contribution due to Halevi [44], which describes the longitudinal projection ($\mathcal{L}$) of the viscoelastic model. As such, the nonlocal correction for plasmonic frequencies, where $\omega \gg \gamma$, is given by $\mathcal{L} \nabla \cdot \left[ \Pi_p + \Pi_{\equiv \text{sh}} \right](\mathbf{r}, \omega) \approx [\beta_{TF}^2 + (4/3)\beta_{el}^2] \nabla \rho(\mathbf{r}, \omega)$. Therein, the often introduced (squared) velocity $\beta_{HF}^2 = \beta_{TF}^2 + (4/3)\beta_{el}^2$ appears due to joint action of pressure and shear [45] (see also Sec. 3 of this Roadmap). It increases the expected LSP
blueshifts. In fact, the Halevi model gives rise to a unified (qualitative) description of line
shifts and broadenings of LSPs in noble metals [45]. In Fig. 10(c), we illustrate the hierarchical
interrelation in Fourier space of the above-described bulk continuum models along with the
correspondingly added dynamics and concomitant effects on LSPs. We emphasize, that only
two bulk parameters, $\omega_p$ and $\gamma$, have to be determined (e.g. via fits to experimental refractive
indices that exclude surface effects [16,45,189]). Specifically, the Fermi velocity represents a
third parameter which occurs in those models that feature a stress tensor and it can be derived
from the metal's equilibrium electron density $n_0$, using $\omega_p = \omega_p(n_0)$.

Besides a bulk material description, a complete material model for plasmonic scatterers
requires an appropriate treatment of the surface. In the simplest case, a description of the metal
by a local bulk response (e.g. Drude) is carried all the way up to a discontinuous transition to
a dielectric environment with spatially-constant bulk response (which includes vacuum) [16]
leading to a sharp boundary. A more general framework is based on a treatment of Maxwell's
equations and the introduction of a sub-wavelength transition layer (also called selvedge [190]).
This layer's response deviates from that of the corresponding bulk medium [191]. In the spirit
of continuum models, the charges, currents and fields in this microscopic surface layer are
integrated out (in surface-normal direction) to yield surface excess quantities. These govern the
(dis)continuities of the bulk fields at material interfaces [191], thereby leading to corrections of
the treatment of sharp boundaries above. The (dis)continuities represent a formal set of boundary
conditions including, as special cases, the classical Maxwell [11] as well as the mesoscopic
boundary conditions discussed in Secs. 3 and 13 of this Roadmap. In principle, these excesses
can be derived from the charges and currents induced in the surface layer of a scatterer with
nonlocal continuum bulk models, assuming the latter to be equipped with additional boundary
conditions (ABCs) [19]. In particular, the surface charges and currents using nonlocal continuum
models can then be used to derive surface-response functions (e.g., Feibelman parameters [192]),
for use in the new boundary conditions (such as, e.g., in the mesoscopic treatment as displayed in
Ref. [16]). In such scattering problems, the lower bound of resolved bulk scales is lifted, such,
that typically a local description of the metal's bulk is used. This scheme is analytically more
tractable and simplifies time-domain simulations [16], while allowing for a successively more
elaborate surface-treatment via continuum-derived excesses (see e.g. the Euler-Drude-derived
Feibelman parameter listed in Ref. [16]). In turn, the latter would enrich the desired database
of Feibelman parameters, discussed in Sec. 13 of this Roadmap. We would like to stress, that
the reduction to a local bulk model removes confined bulk plasmons, e.g. the standing waves in
spheres employing the Euler-Drude model [193].

## Challenges and opportunities

Naturally, with increasing complexity of a continuum model and/or scatterer geometry, the
underlying set of differential equations (completed with initial and boundary conditions) may
not anymore be amenable for analytical treatment so that numerical Maxwell solvers have to be
employed. In this context, discontinuous-Galerkin finite-element methods represent a natural
choice as they are tailored for the spatial discretization of balance equations in the presence
of material interfaces. In a nodal discontinuous-Galerkin time-domain (DGTD) approach
to Maxwell's equations [194,195], the electromagnetic and matter fields are discretized on
unstructured finite-element meshes for adaptive discretization of complex geometries, equipped
with element-local interpolating nodal basis functions of adjustable order for efficient sampling
of the relevant bulk length scales. In order to obtain a global solution, the fields in neighboring
elements are connected via numerical fluxes which facilitates a flexible incorporation of boundary
conditions. The spatial discretization leads to a set of coupled ordinary differential equations in
time which, owing to the element-local discretization, can be solved with high-order explicit time-
stepping schemes (e.g. low-storage Runge-Kutta schemes) that allow to take full advantage of the
efficient spatial discretization. Besides pulsed plane-wave excitations (e.g. for computing spectra and field distributions via on-the-fly Fourier transforms [45,196,197]), the DGTD approach has been extended to point-dipole sources (which allow, among others, to compute Green's tensors, e.g. for Casimir-Polder calculations [198]) and electron beams (for determining energy-electron loss and cathodoluminescence spectra [199–201]). In addition, the numerical scheme allows the treatment of cavities and open systems as well [195].

On the analytical level, a number of Mie solutions with successively refined continuum bulk models exist [16,45,188,202,203]. In fact, with increasing order of the spatial derivatives in Eq. (46) more field modes appear (satisfying separate Helmholtz equations [188]) each with its own wavenumbers $k$ via the implicit longitudinal and transverse dispersion relations [204], $\epsilon_L(k_L, \omega) = 0$ and $k_T^2 = \epsilon_T(k_T, \omega)\omega^2/c_0^2$, that result from the corresponding models' longitudinal and transverse nonlocal dielectric functions $\epsilon_L$ and $\epsilon_T$. With an increasing number of modes, the problem becomes analytically intractable. Independent of this and since the particular choice of ABCs significantly impacts the interplay of these modes, it is highly desirable to develop a general framework to derive ABCs for nonlocal continuum models (beyond typical hard-wall [16] and no-slip conditions [205]).

## Future developments to address challenges

In order to discuss future developments of the continuum framework for plasmonic systems it is, again, instructive to distinguish analytical and numerical pathways.

In the analytical realm, the focus lies on a systematic refinement of the bulk models, where connections to kinetic equations for the statistical dynamics in phase-space have been suggested [45,173,185,187] and relations to orbital-free approximations of dynamical equations appearing in TDDFT have been pointed out [206]. Eventually the refined bulk models should also consider the influence of the Fermi-Dirac statistics on electron-electron interaction. As an early precursor, we would like to mention the incorporation of the Dirac exchange term in Eq. (46) [207] as a function of the density gradient. When keeping same-spin electrons apart, this term counteracts the hydrostatic pressure, thereby reducing the LSP lineshifts in noble metals. However, since this represents an extension of the dynamical equation suggested by Bloch, the need for modifications towards plasmonic frequencies still has to be discussed. Moving beyond noble metals, refined bulk models should address the modifications to the Drude model that are required to accurately describe transition metals [208]. While all the above treatments of the bulk refer to the conduction electrons, a complete model of the optical response of metals also necessitates the (on the linear level additive [209,210]) incorporation of the response of bound electrons due to the significant absorption stemming from interband transitions (IBTs) as well as due to the (often d-band) screening of the conduction electrons away from the IBT range [211] using few parameters with clear physical meaning [189].

Regarding the surface response, derivations of Feibelman parameters based on TDDFT-informed continuum bulk models would be very interesting, since eventually Feibelman parameters are often determined from or checked against TDDFT simulations (see Secs. 12 and 13 of this Roadmap). Further, this would allow the separation of the approximations that are tied to bulk continuum models from those that are tied to use of Feibelman parameters. Further, we suggest to examine the applicability of the excess formalism also to balance equations beyond the level of charge conservation. This will allow the connection to the bulk nonlocal response and should lead to a generalized set of ABCs. Specifically, assuming a planar surface in the xy-plane, the continuity equation has led to the general result [191]

$$J_n(x, y, z \to 0^-) - J_n(x, y, z \to 0^+) = \partial_t \rho_s(x, y) + \nabla_{\|} \cdot J_{s, \|}(x, y) \quad (48)$$

where $\nabla_{\|} = (\partial_x, \partial_y, 0)$ is the gradient with surface-parallel components only and the subscript s denotes an excess.
In the numerical realm, the focus lies on the implementation of the analytical forms of the above-discussed bulk and surface response. On the bulk level, for DGTD a particular difficulty emerges with higher-order spatial derivatives related to the stress tensor such as when introducing the von Weizsäcker potential [212]. As a time-domain formulation, the DGTD allows to treat also nonlinear response. So far, the convective nonlinearity [183] of the full Euler-Drude model has been treated and this could serve as a starting point for similar developments regarding more refined models. The geometrical flexibility inherent in finite-element approaches such as boundary-element methods (see Sec. 13 of this Roadmap) or DGTD should further be leveraged by considering different boundary conditions in order to reflect on the subtleties of the transition layer also for irregular shapes, including e.g. surface roughness. Here we stress, that some surface-response functions such as the Feibelman parameters obtained from the Euler-Drude model with hard-wall boundary conditions [16], are usually derived in frequency domain. Consequently, their inclusion in boundary-element methods is quite advanced (see Sec. 13 of this Roadmap) but, so far, is largely absent for DGTD.

As regards the analytical and numerical studies as well, modern plasmonic studies are also focused on the interaction with external resonances, e.g. in core-shell particles, where excitons of dye-shells (like J-aggregate layers) couple to the plasmons in the metallic core. While aiming for more realistic descriptions, it has been realized that the choice of model for each bulk material determines the nature and also amount of resonances of the core-shell system (see, e.g., Ref. [213]). We also suggest to explore the role of metallic bulk nonlocality and surface roughness on the matching of plasmonic and excitonic resonance as well as the losses, determining whether strong coupling prevails. Further, the applicability of continuum bulk models in other materials, such as doped semiconductors [214] or meta-materials, has to be studied, similar to the suggestion on such attempts for surface response functions in Sec. 12 of this Roadmap.

## Concluding remarks

Overall, the continuum framework for plasmonic systems provides an efficient and extensible description of the collective dynamics of conduction electrons. Depending on the resolved scales and, thus, accuracy, either nonlocal bulk models are matched with additional boundary conditions or local models are equipped with mesoscopic boundary conditions, which can be derived from the charges and currents originating in the former case. The established plasmonic continuum models have demonstrated the necessity to distinguish between longitudinal and transverse bulk nonlocality, with analytical solutions being available for selected combinations of material models and geometries. To increase the class of tractable material models and geometries the DGTD has proven itself as a flexible and resource-efficient Maxwell solver.

**Funding.** Deutsche Forschungsgemeinschaft (398816777–Project A06, 510943930–Project A06).

## 8. Quantum hydrodynamic model

HUATIAN HU ID, FABIO DELLA SALA ID, PU ZHANG ID &
CRISTIAN CIRACI ID

### Overview

Optical spatial nonlocality, a spatially dispersive effect where local polarization depends on nearby fields, is an omnipresent phenomenon in natural materials. In fact, beyond the temporal nonlocality that leads to the well-known frequency dispersion, spatial dispersion depending on the momentum is often neglected in conventional photonics involving bulk materials. This effect is instead very relevant in the context of nanophotonics, especially plasmonics, where light can be concentrated into a deep-subwavelength region. In this case, the near field of the
light can resolve the medium’s microscopic spatial graininess, which violates the conception of homogeneous constitutive parameters in the bulk materials with an averaged field. Once the scale of light-matter interaction is confined down to nanometer-scale, quantum effects play a dominant role and the dielectric constant becomes nonlocal.

Current status

The Kohn–Sham (KS) density-functional theory (DFT) approach and its time-dependent extension (TD-DFT) is a powerful method for accurately describing quantum effects [215] (see also Sec. 9 of this Roadmap). However, its computational cost is prohibitively high for mesoscopic plasmonic nanostructures containing billions of electrons. Such systems are predominantly governed by the collective behavior of electrons rather than individual single-particle properties, allowing for a trade-off between accuracy and the computational demand of calculating each electron’s orbital. Orbital-free density functional theory (OF-DFT) methods [206,216,217], which approximate the kinetic energy functional by assuming an effective orbital, provide an efficient alternative for addressing these large-scale systems. They are less accurate but much faster than KS methods. Building on this principle, the quantum hydrodynamic theory (QHT) [60,61,129,206] provides a self-consistent framework that incorporates microscopic quantum effects—such as nonlocality, electron spill out, and quantum tunneling—while seamlessly integrating with full-vector Maxwell’s equations. This approach effectively predicts both near- and far-field properties of mesoscopic systems. It borrows the concept of electron density $n(\mathbf{r}, t)$ from DFT, together with another macroscopic variable—the electron fluid velocity $v(\mathbf{r}, t)$, enabling the dynamics of the electron fluid to be expressed as:

$$m_e \left( \frac{\partial}{\partial t} + \mathbf{v} \cdot \nabla + \gamma \right) \mathbf{v} = -e (\mathbf{E} + \mathbf{v} \times \mathbf{B}) - \nabla \frac{\delta G[n]}{\delta n}. \quad (49)$$

where all quantum effects are collected into the functional $G[n]$. The functional $G[n]$ can be expressed as

$$G[n] = T[n] + E_{XC}[n] \quad (50)$$

i.e. as the sum of non-interacting kinetic energy, $T[n]$, and exchange-correlation potential energy $E_{XC}[n]$, based on the level of approximation that is considered. Notably, if the energy functionals $G[n]$ are neglected, in the linear regime Eq. (49) will reduce to the Drude model, i.e., local-response approximation (LRA), as shown in Fig. 11.

The Fermi theory of a free-electron gas introduces electron-electron interactions leading to a collective repulsion of electrons due to the Pauli exclusion principle. As a result, when electrons are forced toward the material boundary by an external electric field, they spread (see the induced density $n_1$ of middle panel of Fig. 11). This approximation, denoted as $G[n] = T_{TF}[n]$, is referred to as Thomas–Fermi hydrodynamic theory (TF-HT). To properly describe the boundary behavior, additional conditions—so-called hard-wall boundary condition (HW-BC)—must be imposed, ensuring that the normal component of the polarization vanishes. Physically, the HW-BC imply that electrons cannot cross the boundary. TF-HT has been shown to successfully account for the blueshift of plasmonic resonances and the decrease of the field enhancement in confined regions [55,83].

The TF-HT framework is not always accurate because it often considers the background equilibrium density ($n_0$) to be constant and strictly confined within the material boundaries, neglecting the quantum effects such as electron spill out and tunneling. To address these limitations QHT [60,61,129] incorporates a more advanced KE approximation by employing a fraction ($\lambda$) of the von Weizsäcker (vW) kinetic energy functional, namely

$$T[n] = T_{TF}[n] + \lambda T_{vW}[n] \quad (51)$$

together with the exchange-correlation with local-density approximations (LDA). We note that XC functionals beyond LDA can also be employed. However, given the approximations
<figure><img src="image_13.png" /><figcaption>**Fig. 11.** A summary of three different response models of a free-electron gas. In the local-response model, the induced charge density is a Dirac delta function on the metal surface. The background equilibrium density $n_0$ is uniform and confined in the ion boundary; in the Thomas-Fermi hydrodynamic theory (TF-HT), $n_0$ is uniform as local-response model while the induced electrons accumulate near the metal surface without being able to escape (hard-wall condition); in the quantum hydrodynamic theory (QHT) the equilibrium density is exponentially decaying across the metal surface and the induced charges are smeared both inside and outside the metal-air interface. Different functionals used are summarized as below.</figcaption></figure>

inherent in the leading term—the kinetic energy (KE) functional—it is not particularly beneficial to use more sophisticated XC functionals in this context. The KE functional in Eq. (51) is also widely used for ground-state calculations of finite systems within OF-DFT. The OF-DFT ground-state calculation in the jellium model is equivalent to solve the zero-order equation on the QHT, which can be easily coupled to the Poisson’s equation to account for the electrostatic potentials [218]. Additionally, the near-field potential from the ionic lattice and the electron affinity of the interfacing dielectric become relevant when considering electron spill-out and its tunability at the interface between, for example, a noble metal and a dielectric [219]. Thus the spatially-varying equilibrium electron density $n_0(r)$ can be calculated [Fig. 11(right), for metallic structures of any shape. Structures such as dimers and other nanogap configurations exhibiting electronic quantum tunneling effects can be accurately modeled using this framework, eliminating the need for additional boundary conditions. More precisely, in the context of QHT, the concept of well-defined boundary conditions, such as HW-BC of TF-HT, is no longer applicable since electrons can cross boundaries and decay in space, resulting in continuous fields throughout the system. However, realistic simulations require the physical domain to have finite dimensions, and the choice and truncation of the simulation region might significantly influence the accuracy of the results. Greater accuracy and stability in QHT can be achieved by improving the KE functional, incorporating for example Laplacian-level kinetic energy contributions [130]. Moreover, numerical schemes for modeling the nonlinear optical response of QHT have been proposed [220], along with efficient simulation approaches that can handle full three-dimensional systems [221].

Opportunities and challenges

QHT is a versatile multiscale numerical tool well-suited for a broad range of applications in systems that require quantum-mechanical modeling of electrons alongside the ability to account for macroscopic optical interactions.

**Applicability beyond metals.** In recent decades, QHT has been successfully applied to describe optical properties of alkali, e.g., sodium (Na), and noble metals, e.g., gold (Au) and silver (Ag) [219]. The study of alkali metals holds fundamental significance as it validates QHT by benchmarking its predictions against DFT calculations. On the other hand, research on noble metals has its practical relevance as these materials are the most used in conventional plasmonic applications. However, extending the application of QHT to other material systems has
been limited. More recently, significant efforts have focused on exploring alternative materials that support plasmons, including degenerate electron systems such as transparent conduction oxides (TCO), e.g., indium tin oxide (ITO) in the near-infrared, graphene, and heavily doped semiconductors, e.g., indium gallium arsenide (InGaAs) and indium phosphide (InP), in the mid-infrared. These materials exhibit exceptional linear and nonlinear optical properties, as well as electrical tunability. TF-HT has been employed to characterize the linear and nonlinear optical responses of ITO [222] and doped semiconductors [223]. Intersubband transitions in heavily doped semiconductor systems, or metallic quantum wells (QWs), present particularly intriguing phenomena since their properties can be tuned externally, and they exhibit exceptionally strong optical nonlinearities [224]. When doping levels are sufficiently high, such that the Fermi level can cover the conduction band and a few QW states, the intersubband transitions become plasmonic since the potential has a much lower contribution to the Hamiltonian. Organized into nanometer-scale layers, these heavily doped semiconductor or metallic QWs can support remarkably strong optical nonlinearities [225,226]. Conventionally, QW properties are solved using Schrödinger–Poisson equations [227] and other phenomenological models [228], but their nonlinear susceptibilities remain complex to predict. Preliminary work based on TF-HT has been undertaken to address the strong nonlinearity of such QWs [229]. However, more precise predictions of QHT applied to these and other emerging material systems remains a promising yet largely unexplored avenue for future research.

**Nonlinear nano-optics empowered by QHT.** Among the various applications of QHT, those focused on nonlinear phenomena are particularly compelling. Contrarily to effective-response approaches, such as the use of Feibelman *d*-parameters [192] (see also Secs. 12 and 13 of this Roadmap), or quantum-corrected methods [78,230], which are often restricted to linear response regimes, QHT inherently provides a comprehensive, fully dynamical framework. Moreover, while nonlocal and electron spill-out effects are traditionally regarded as “unwanted” side effects that increase broadening and degrade field enhancement, from a nonlinear perspective, these effects can instead enhance the free-electron response [231]. Low-electron-density Drude materials such as TCOs, and heavily doped semiconductors can indeed support high optical nonlinearity, as their effective nonlinear susceptibilities scale inversely with some power (depending on the order) of the electron density [223] as well as their effective electron mass [232]. Although TF-HT has been successfully applied to describe heavily doped semiconductors [233], ITO [222], and silicon (Si), these approaches rely on the assumption that spill-out effects can be neglected. By contrast, QHT—accounting for spatially varying equilibrium electron densities—becomes particularly relevant for configurations where geometry dimensions are comparable with the typical scale of electron spill-out. Examples include QWs and atomic scale protrusions supporting extremely localized optical modes (see Sec. 10 in this Roadmap). Additionally, the local spatially varying electron density can indeed have a significant impact either directly on the nonlinear optical properties [231,234] or on the linear response in a nonlinear manner through externally controlled $n_0$ [235,236]. A self-consistent QHT proves very suitable for both scenarios.

**Longitudinal bulk plasmons.** A longitudinal bulk plasmon (LBP) is an electronic density wave within the bulk material, enabled by nonlocal effects in the medium's constitutive relations. LBPs have been extensively studied, both theoretically and experimentally, for decades in free-electron gas materials with TF-HT [51,204,214,237,238]. Yet, their application scenarios remain largely unexplored, primarily due to the challenges associated with efficiently exciting them in conventionally sized structures. Recent advances in material science and nanofabrication have however have begun to unlock new possibilities. Unlike localized surface plasmons, which occur at the material boundaries at lower energy, LBPs resonate above the plasma frequency in the bulk region of the material. This distinction leads to unique opportunities: Firstly, the resonance of LBPs lies above the plasma frequency, which correspond to the epsilon-near-zero (ENZ) region commonly associated with TCOs. This allows LBPs in TCOs to possibly operate
in the near-infrared or even visible spectrum, hereby providing a valuable complement to silicon photonics. Secondly, since the LBP optical mode profile reaches into the bulk regions, they offer larger active interaction volumes compared to surface-bound modes. This larger volume enables higher nonlinear efficiencies. Recent studies using TF-HT have revealed that LBPs in plasmon-QW hybrid structures can exhibit exceptionally high Kerr nonlinearities and low-power-threshold optical bistability [229]. Looking forward, the application of QHT is expected to unveil new physics in this domain. By incorporating a more accurate description of spatially varying electron densities and dynamical responses, QHT can provide deeper insights into LBP behavior and open pathways for innovative applications in nonlinear and quantum photonics.

**QHT as a constitutive description of nonlocal media.** QHT couples with Maxwell's equations to form a closed description for plasmonic responses, it effectively serves as a constitutive relation of spatially dispersive metals and semiconductors. Building upon the full-wave solution obtained by solving the coupled Maxwell-QHT equations, a more formal and comprehensive electromagnetic theory for nonlocal media emerges. This theory, grounded in QHT, is particularly useful for analyzing the optical responses and bridging gaps with other models. Notably, QHT, as a continuum theory, has recently been integrated into the generalized Lorentz model (GLM) [239]. Within the theoretical framework of GLM electromagnetic theorems such as Poynting theorem and Lorentz reciprocity are generalized to nonlocal media. This extension further enables the development of a QHT-based quasinormal mode (QNM) theory (see Sec. 10 of this Roadmap). The QNM theory not only provides modal analysis, as evidenced with the classical QNM theory, it promises a route to field quantization through the QNM expansion of Green's function [240]. The elevation of QHT to a quantized theory is yet still an unaccomplished feat.

*Future developments to address challenges*

While QHT can incorporate functionals to capture both linear and nonlinear optical properties of degenerate electron systems, there remain considerable challenges to be solved.

**Full self-consistency.** Due to the mathematical complexity of QHT equations, these are often solved using a perturbation approach where the equilibrium (ground-state) electron density $n_0$ is calculated independently by solving the zeroth-order equation [218]. Then, one assumes that the sum of all higher-order perturbations remains much smaller than $n_0$. However, this assumption becomes less straightforward near boundaries, where $n_0$ exhibits exponential decay, where a different asymptotic decay velocity might render this condition false for any arbitrary small perturbation. In the context of quantum systems, one must ensure that $n = n_0 + \delta n > 0$, while this is straightforward when $\delta n > 0$ it might be problematic when $\delta n < 0$. From a foundational perspective, this highlights the need for a time-dependent, fully self-consistent model to address QHT and eliminate the perturbative shortcomings [241, 242]. A non-perturbative, time-dependent treatment would account for saturation effects and allow for the full exploration of plasmon dynamics, including modulation effects driven by changes in $n_0$.

**Accuracy of QHT.** The accuracy of QHT is strongly influenced by the choice of KE functionals, far more so than by the exchange-correlation energy. A detailed discussion on the subject can be found in Ref. [206]. Two main pathways for enhancement can be identified:

The first approach involves better capturing the nonlocality of the KE functional. In the TF-HT framework the KE functional is entirely local, although the use of Eq. (51) introduces a semilocal description. Laplacian-level (LL) KE functionals, such as those described in Ref. [130], further increase nonlocality, but they also introduces significant complexity in practical implementations. This is particularly problematic for nonlinear applications, where expansions in LL-QHT become quite cumbersome. In solid-state physics, fully nonlocal KE functionals based on the Lindhard function [9], which describes the exact response of the
homogeneous electron gas, have been extensively studied [216,217]. Unlike LL-KE functionals, these nonlocal functionals are not derived from high-order derivatives of the electron density but instead use convolutions of the electron density with the Lindhard function. These nonlocal KE functionals are, however, typically defined in the reciprocal space and thus limited to periodic systems and/or finite systems inside a periodic cell. For applications in nano-optics, where far-field scattering calculations are of interest, KE functionals must be defined in the real-space. Recently, fully nonlocal functionals in real-space have emerged [243]. These functionals involve the calculation of the screened Coulomb potential and its functional derivatives, which can be implemented within the QHT framework for nano-optics. Thus, such functionals can represent an interesting direction to follow to increase the QHT accuracy for diverse systems.

The second pathway is to address the frequency dependence—or non-adiabaticity—of the KE functional, as the exact $G[n]$ functional also depends on time. In the frequency-domain, using Eq. (51), the KE contributions to the QHT are static and real-valued. However, the exact KE contributions are frequency-dependent and complex-valued [206], with the imaginary part accounting for the broadening [218]. This issue have been discussed in Refs. [244–246]. For just two electrons, the (static) vW KE functional is exact and the QHT linear response exactly coincides with the reference TD-DFT results [206]. For many electron systems, a frequency dependent KE is required [246]. However, to date only nonlocal approximations derived from the Lindhard function exists [244–247], which cannot be directly applied for real-space nano-optics QHT implementation. Finding a spatially nonlocal and frequency-dependent KE functional in real-space is a challenging interesting direction for future developments.

For both the above development directions machine-learning can be a valuable tool [248–251].

### Concluding remarks

QHT is a rapidly evolving and versatile tool for modeling quantum effects in nanophotonics. By incorporating quantum corrections into macroscopic optical models, QHT enables the accurate prediction of both linear and nonlinear optical responses, bridging the gap between quantum mechanics and classical electrodynamics. Despite its successes, significant challenges remain, particularly in expanding its applicability to nonlinear and more complex systems.

Future advancements in QHT hinge on two critical fronts: the development of more accurate and versatile KE functionals and the adoption of fully self-consistent, time-dependent numerical solvers. Addressing the nonlocality of KE functionals will be essential for extending QHT's utility to a broader range of materials and geometries.

QHT's ability to describe intricate quantum phenomena such as electron spill-out, quantum tunneling, and longitudinal bulk plasmons positions it as a cornerstone for the future of quantum-enabled nanophotonics. By addressing its current limitations and leveraging advancements in computational techniques, QHT can pave the way for transformative innovations in the study of light-matter interactions at the nanoscale.

**Funding.** HORIZON EUROPE European Innovation Council (NEHO, project No. 101046329); ICSC—Centro Nazionale di Ricerca in High Performance Computing, Big Data and Quantum Computing, funded by European Union—NextGenerationEU—PNRR.

## 9. Jellium model within time-dependent density-functional theory to address quantum effects in nanoplasmonics

JAVIER AIZPURUA ID, ANTTON BABAZE ID & ANDREI G. BORISOV ID

### Current status

The collective response of a free-electron gas confined by a metal surface or by the boundaries of a nanoparticle has been one of the main exponents of light control and manipulation at the
nanoscale during the last decades. The capacity of a finite structure to sustain localized surface plasmons provides one of the most effective mechanisms to localize light below the diffraction limit in extremely small effective mode volumes, associated with a very large electromagnetic near field. Classical approaches to calculate the interaction of plasmonic nanoparticles with external electromagnetic field are typically based on solving Maxwell's equations within the linear-response theory, where materials are described by local frequency-dependent dielectric functions. These approaches have been very successful to address optical properties in field-enhanced spectroscopy and microscopy. However, nowadays the building blocks in nanoplasmonics reach the realm of the atomic scale, such as in ultranarrow plasmonic gaps, in metallic clusters, in emitters closely attached to metallic surfaces, or in picocavity configurations. Strong quantum effects thus emerge, derived from the inhomogeneous electron density distribution at the boundaries as well as from nonlocal dynamical screening [58,252,253].

A powerful method to account for quantum effects in the optical response of metal surfaces or metal nano-objects relies on the time-dependent density-functional theory (TDDFT) [254]. Different flavors of TDDFT have been developed in condensed-matter physics and in quantum chemistry, depending on the properties to be addressed. The most widely used approach is based on the Kohn–Sham (KS) scheme. The electron density $n(r, t)$ is represented using effective single-particle orbitals $\Psi(r, t)$, which evolve in time under the action of the KS potential $v_{KS}(n; r, t)$. This potential comprises the Hartree, exchange-correlation, and external potential, as well as (pseudo) potentials representing the interaction of electrons with nuclei.

The solution of the KS equations within TDDFT in space and time allows for obtaining key quantities that characterize the linear and nonlinear optical response of a nanosystem, and is a state-of-the-art full *ab initio* approach at the forefront of modern theoretical capabilities to address such response. This approach fully accounts for quantum effects such as i) nonlocal dynamical screening, ii) finite-size effects and electron spill-out, and iii) charge transfer in close-contact nanoparticle surfaces [58,252,253]. In a metal–dielectric or metal–insulator–metal configuration, the TDDFT properly describes multiphoton and optical-field induced photoemission, as well as optically assisted electron tunneling from one metal to another. This allows for addressing challenging goals in nanoscale transport related to the engineering of petahertz (PHz) electronic devices [255].

## Challenges and opportunities

Successful applications of TDDFT within a full atomistic description of nanoparticles have paved the way towards understanding the ultimate limits of field localization in nanoplasmonics, as well as describing the nonlocal and band structure effects related to the size dependence of plasmon resonances in (noble) metal clusters [256,257]. However, the high computational cost of the method prevents studies on large-scale plasmonic systems with a wide parameter variation, often necessary to gain physical intuition and “nail down” the main nonlocal plasmonic effects. This challenge is met with the free-electron (jellium) model of TDDFT (JM-TDDFT). In this method, the ionic cores of the metal are not treated explicitly, but represented by a uniform charge density, $n_+ = (3/4\pi)r_s^{-3}$, where $r_s$ is the Wigner–Seitz radius. Despite its simplicity, JM-TDDFT has proven to be a reliable and adequate many-body description of the dynamics of the conduction band electrons in a metal under optical excitation. Results of this approach have been pivotal for the prediction of quantum phenomena in plasmonics, such as the effect of electron tunneling across a plasmonic gap and electronic coupling between a quantum emitter and a plasmonic nanoantenna. Moreover, the nonlocal optical response of metallic nanoparticles and nanogaps obtained within the JM-TDDFT has been a key reference for benchmarking other theories [58,252,253,255].

The results of the properties of plasmon modes of metal nanoantennas within the JM-TDDFT, as well as those regarding the Green’s function of a point-like dipolar emitter in proximity to a
plasmonic dimer highlight the good performance of the dispersive surface-response formalism which incorporates nonlocality also in the direction parallel to the metal-vacuum interface [253]. This formalism is based on the implementation of Feibelman parameters [192] (see also Secs. 12 and 13 of this Roadmap), $d_{\perp}(\omega, k_{\parallel})$ and $d_{\parallel}(\omega, k_{\parallel})$, into modified boundary conditions of Maxwell's equations [59,258] to capture extreme nonlocality at interfaces. The real part of the Feibelman parameter $d_{\perp}(\omega, k_{\parallel})$ accounts for the position of the centroid of the induced charge at a given frequency $\omega$ and at a given parallel momentum transfer $k_{\parallel}$, and the imaginary part determines the probability of creating electron-hole pairs in the surface region. Crucially, it has been demonstrated that using Feibelman parameters that explicitly depend on $\omega$ and $k_{\parallel}$ is significantly more accurate than using Feibelman parameters obtained within the long-wavelength approximation [253]. These Feibelman parameters $d_{\perp}(\omega, k_{\parallel})$ can be obtained using JM-TDDFT for a flat metal surface [Fig. 12(a)].

<figure><img src="image_14.png" /><figcaption>**Fig. 12.** Areas of development of nonlocal nanoplasmonics within the JM-TDDFT approach. (a) The calculation of Feibelman parameters requires the description of the metal to be extended to include the effect of localized electron bands and atomistic crystallographic facets. (b) The jellium description of the electron gas allows for performing first-principles calculations of its intrinsic nonlinear optical response together with the effect of the nanoantenna geometrical boundaries, as in symmetry-induced nanoplasmonic response to structured light. (c) Ultrafast single cycle optical pulses induce photocurrents in nanogaps, determined by the dynamics of sub-fs electron bursts, as described by jellium-based electron models.</figcaption></figure>

One of the main challenges of the surface-response formalism based on Feibelman parameters is the limited availability of quantum-informed parameters for the different pairs of metal-dielectric interfaces commonly used in experimental realizations of plasmonics. So far, the jellium model of metals has been implemented to generate quantum-informed Feibelman parameters, which allows for addressing nonlocality of free-electron plasmonic nanoparticles [aluminum (Al) and alkali metals such as sodium (Na)] in vacuum. Such a description is not adequate for silver (Ag) and gold (Au), routinely used in spectroscopy and microscopy experiments. In these metals, the localized d-band electrons strongly impact the optical response. Furthermore, the dielectric environment also affects the Feibelman parameters by changing the spill-out of the electron
density and the dynamical screening. To address these effects, a full-atomistic *ab initio* TDDFT approach beyond JM-TDDFT needs to be applied (see also Sec. 10 of this Roadmap), which makes the calculations extremely challenging numerically.

On the other hand, the lighter and more versatile JM-TDDFT opens the door to tackle other nonlocal and, more generally, quantum effects of light-matter interaction in nanoplasmonics, where the underlying physics is determined by the conduction electrons. Thus, the intrinsic *nonlinear electromagnetic response* of conduction electrons leading to frequency conversion in nanoplasmonic configurations can be directly obtained within the jellium scheme [Fig. 12(b)], enabling the identification of different sources of nonlinearity at the microscopic level. Furthermore, JM-TDDFT allows one to describe the field-induced electron emission from metallic surfaces and tips triggered by single-cycle optical pulses. In particular, when a metallic surface is facing another one in close proximity [Fig. 12(c)], the bursts of photoemitted electrons exhibit an intriguing and complex *femtosecond dynamics* that depends on the waveform of the incident transient pulse [255]. This opens up perspectives for designing nanoscale PHz optoelectronic devices.

## Future developments to address challenges

Both short- and long-term developments will be required in the implementation of many-body calculations of light-matter interaction in the context of nanoplasmonics. The calculation of Feibelman parameters for interfaces of non-free-electron materials will require the development of atomistic *ab initio* calculations of dispersive Feibelman parameters. This will enable the direct implementation of nonlocality in a set of materials such as gold (Au), silver (Ag), copper (Cu), platinum (Pt), palladium (Pd), etc., interfacing dielectrics, of high interest in field-enhanced catalysis and chemical reactivity. The role of the crystallographic facets and of surface-localized electronic states could thus be addressed in this way.

In nonlinear nano-optics, specific combinations of the intrinsic nonlinear response of the electron gas with specific geometrical symmetries shown by certain nanoantennas configurations could be exploited in nonlinear chiral nanoplasmonics, with the full nonlocal and nonlinear response directly accounted for by the jellium description of the electron gas. Finally, the fs photocurrents produced by intense single-cycle optical pulses can be studied in the context of electron transport in STM junctions, with the ultimate goal of providing theoretical support to understand and design devices that combine atomic-scale spatial and sub-fs temporal resolution.

## Concluding remarks

The jellium model (JM-TDDFT) has been a powerful tool for addressing conduction electron dynamics and nonlocal effects in metallic nanoparticles during the last years, offering a computationally efficient yet accurate framework for exploring fundamental quantum phenomena in nanoplasmonics. Moreover, it is evident that this methodology holds significant potential for future advancements, particularly in deepening our understanding of nonlinear nano-optics and ultrafast optoelectronics. At the same time, the development of more sophisticated atomistic descriptions beyond the jellium model will be crucial for advancing in the development of quantum-informed nonlocality. Together, these developments will not only expand our fundamental knowledge of light-matter interaction but also contribute to the design of applications in quantum technologies and ultrafast photonics.

**Funding.** Spanish Ministry of Science and Innovation (PID2022-139579NB-I00); Basque Government program for groups of the University of the Basque Country (IT1526-22); Deutsche Forschungsgemeinschaft (510943930–Project A06).
# 10. Atomic scale nonlocal plasmonics

PU ZHANG ID & XUE-WEN CHEN ID

## Overview

Nanoplasmonics is continuously pursuing ever tighter light confinement. With the conceptual innovations and advances in nanofabrication, the compressed plasmonic field inevitably encounters the eventual graininess of material, or reaches the atomic scale. Studying the impact of the atomic scale structure on the plasmonic response and related applications constitutes atomic scale plasmonics. Holding the promise of achieving the ultimate light confinement and enormous enhancement of light-matter interaction, this research direction opens up an avenue for the plasmonics realm. In this contribution we start with a brief account of atomic scale plasmonics from its commencement to the current status, necessarily reflecting the personal view of the authors. Then we identify some opportunities and challenges faced by this flourishing field. Emphasis is placed on the nonlocal effects and theoretical analysis of the extremely localized field.

## Current status

Spatial confinement of light is one of the most important strategies to boost light-matter interaction, which lies at the core of the majority, if not all, of light technologies. Conducting material-based nanostructures with abundant and freely moving electrons serve as the driving force for light confinement, as they are able to leverage both the non-resonant lightning rod effect from morphology induced electron accumulation and the resonant augment from the localized surface plasmon resonance. With the typical nanoplasmonic constructs such as sharp tips, nanogap structures and cascaded assemblies, the localized optical field reaches the nanoscale and even smaller. Going a step forward nanoplasmonics encounters the graininess of material. Atomic scale plasmonics was initiated with theoretical exploration. From the theoretical perspective, classical electromagnetic theory with the local response approximation, or the description of material by the permittivity $\epsilon(\omega)$, is normally applied in plasmonics, but starts to fail for plasmonic field squeezed into nanoscale spaces because of the emerging nonlocal effects. The alternatives include continuum nonlocal theories (see Sec. 3 of this Roadmap), which effectively provide a description based on a spatially dispersive permittivity $\epsilon(\omega, k)$, and *ab initio* approaches. Being a representative of the latter, time-dependent density-functional theory (TD-DFT) encompasses all the non-classical effects of electrons, and is originally developed to treat atomic structures. Using TD-DFT beyond the jellium approximation (see Sec. 9 of this Roadmap), the impact of the atomic structure on the plasmonic near field was first investigated and highlighted in a study on the sodium (Na) cluster dimers consisting of over 600 atoms in total [259]. As shown in Fig. 13(a), the near-field pattern clearly has the imprint of the detailed atomic structure. The atoms at the cluster surface shape the localized field into irregular distributions. Intense field even appears in the atomic voids and crevices. The hot spot in the plasmonic near field is of particular importance. Much attention was focused on the formation of hot spots in similar atomic cluster dimers. In particular, an asymmetric dimer in the tip-to-facet configuration illustrated in Fig. 13(b) was successfully predicted with TD-DFT to generate hot spots of down to $0.4 \text{ nm}^2$ lateral area [256]. The underlying mechanism is recognized as the non-resonant lightning rod effect. The first experimental observation relevant to atomic scale plasmonics was reported in low-temperature time-series surface-enhanced Raman spectroscopy (SERS) [260]. Fluctuational bright flashes were spotted in the real-time spectra besides persistent signals. The experiment was interpreted with the introduction of the “picocavity” in Fig. 13(c), where the plasmonic hot spot at a gold (Au) adatom activated by laser irradiation couples with an individual molecular bond. The interpretation was supported by the experimental observations including the extreme enhancement of the optomechanical coupling and alteration of SERS selection
rules by strong optical field gradient. The picocavity has a Raman localization volume below 1 nm³. Apart from the enormously enhanced Raman scattering, this highly localized hot spot boasts great potentiality in various types of ultrasensitive detection and sensing. For instance, a silver (Ag) tip decorated with an Ag cluster has been shown able to realize sub-nanometer resolution single-molecule photoluminescence imaging [261]. In this case the atomic structure, or the Ag cluster, is transferred to the tip apex during tip indentations, instead of activated by laser stimulation. Nevertheless, neither methods can deterministically produce stable atomic protrusions. The beyond-single-adatom stable atomic protrusion sketched in Fig. 13(d) is being actively sought for [262,263].

<figure><img src="image_15.png" /><figcaption>**Fig. 13.** The representative developments and trend of atomic scale plasmonics. (a) An early study on the plasmonic response of atomic clusters via an ab-initio approach reveals the impact of the atomic structure on the local field pattern and strength [259]. (b) The atomic cluster dimer in the tip-to-facet configuration is predicted with TD-DFT to give rise to atomic scale lightning rod effect [256]. (c) Experimental report on “picocavities” based on the observation of bright fluctuational flashes in the real-time SERS spectra [260]. (d) The beyond-single-adatom atomic structure is a preferred configuration for atomic scale plasmonics owing to the stability at room temperature.</figcaption></figure>

## Challenges and opportunities

The vibrantly progressing research area of atomic scale plasmonics certainly faces prominent challenges, which in the meanwhile signify great opportunities. In this section we name a few outstanding ones in hope that they could inspire more research interest in the area. The most fundamental question would be “What is the ultimate limit of light confinement?”. Before any meaningful discussion on the topic, what needs to be resolved in the first place is how to characterize quantitatively the degree of light confinement with a consistent standard. While mode volume is routinely used for the purpose, the way it is evaluated varies a lot from case to case. For instance, Ref. [260] obtained the 1 nm³ Raman localization volume referring to the Raman scattering enhancement by the picocavity. A more widely accepted approach is
the mode volume according to the mode theory of optical microcavity. However, the nonlocal effects including electron nonlocality, electron spill out and nonlocal damping, manifested in atomic scale plasmonics, call for a non-classical theory. Quantum hydrodynamic theory (QHT) elaborated on in Sec. 8 is a viable option. The quasinormal mode theory (QNM) has recently been generalized to QHT, so that a mode theory-based mode volume can be defined in atomic scale plasmonics [239,264]. Therein the nonlocal effects are encapsulated into the restoring force $\hat{\Theta}$ and damping $\hat{\Gamma}$ operators in the generalized Lorentz model (GLM) [239]

$$i\rho\epsilon_0\omega_p^2\mathbf{E} - i\rho\hat{\Theta}\mathbf{P} - i\rho\hat{\Gamma}\mathbf{J} = \omega\mathbf{J} \qquad (52)$$

which constitutes a general theoretical framework for any continuum linear response theory and facilitates the establishment of the corresponding QNM theory. According to the QHT-QNM theory, picocavity working under the non-resonant lightning rod mechanism has a several tens of nm³ mode volume. Plotted in logarithmic scale the local field still spreads across the whole nanogap as illustrated in Fig. 14(b). In contrast, the QHT-QNM theory helps discover an extremely localized mode (ELM) supported by the atomic protrusion [264]. At the resonant frequency its mode volume turns out about 1 nm³, with the modal field in Fig. 14(c) extremely localized at the protrusion. Moreover, the ELM can be made very radiative with a 30% efficiency when coupled with a host antenna. The resonance of the protrusion was also pointed out through the local response calculations [265,266]. Nevertheless, whether 1 nm³ mode volume represents the ultimate light confinement is still an open question. When an optical field is localized on the atomic scale, or commensurate with the extent of the electron wave function, it is naturally expected susceptible to the ambient conditions. Apart from serving for ultrasensitive detection, tuning the optical response with another physical factor in a controlled manner is a promising avenue. At an atomic scale protrusion, the minute number of electrons responsible for the optical response implies that the tuning could be ultrafast and energy efficient. In fact, by using a QHT tamed for noble metals and the corresponding QNM theory, electro-optic modulation of the ELM has been proposed [219, 235]. A 20 fs response time and below 100 aJ per-bit energy consumption were predicted. We note electrical control is far from the only way to modulate the atomic scale plasmon. It is tempting to explore, both theoretically and experimentally, such modulation via e.g. thermal and mechanical stimuli in multi-physics systems. Since the debut of atomic scale plasmonics, the interactions with quantum processes like molecular vibration [260], photoluminescence [261] and exciton [267] are of particular importance for fundamental interest

<figure><img src="image_16.png" /><figcaption>Fig. 14. Field confinement at an atomic scale protrusion. (a) The schematic of the protrusion in the nanoparticle-on-mirror (NPoM) configuration. The logarithmic-scale normalized magnitude of the modal electric field around the protrusion at the (b) non-resonant and (c) resonant frequencies calculated with quantum hydrodynamic theory based quasinormal mode theory.</figcaption></figure>
and potential applications. The tightly confined optical field is able to remarkably expedite the processes and as well to unlock originally forbidden ones by breaking the selection rules with the strong field gradient [260]. There's still much room left to tap the ELM for boosting light-matter interaction, but it is accompanied with intricate issues. Currently the typical atomic scale protrusions are generated dynamically through laser activation, and become destroyed equally fast at room temperature. The probabilistic nature severely limits the practical application of picocavity enhanced processes. In addition, the lack of controllable realization of an atomic scale protrusion prevents matching the protrusion's resonant frequency with the intrinsic frequencies of the quantum processes. Almost all the existing studies only leverage the non-resonant lightning rod effect of the protrusion.

## Future developments to address challenges

Facing the challenges we then discuss future developments hopeful for answering the questions. Light confinement needs to be characterized with the consideration of the working wavelength. In the visible and near infrared range, light confinement has been reported theoretically to reach 1 nm³, but the limit was not addressed. Since the electron response to optical excitation is at the bottom of the mechanisms of light confinement, the length scale of the ultimate confinement, if it exists, should be related to that of the electrons in the material, or the extent of the electron wave function. What role the latter plays in determining the confinement is to be explored. Full quantum or some sophisticated approximate atomistic approaches would be required for the task. On the other hand, there are other known strategies for light localization than those leading to the ELM. Especially applying them to the judicious design of the host of the atomic scale structure is expected to further compress the field. To exploit the potential of using atomic scale plasmon for modulation, multi-physics scenarios become highly relevant. The existing theories and numerical models for plasmonic response must coherently incorporate the controlling physical factor of interest to characterize its nonlinear coupling with the plasmonic field. In this regard, the QHT is a convenient theoretical platform, as it directly deals with the foundational quantities of electron and current densities. The influences of the controlling factors on the densities are usually easy to introduce in a first-principles fashion as in electro-optic modulation [235]. A comprehensive modulation study necessarily involves inspecting the dynamic tuning process, and thus a time-domain theory is indispensable. At least for QHT, however, the time-domain implementation is challenging (see Sec. 8 of this Roadmap), let alone one generalized to multi-physical systems. Aside from used for information processing, modulating the atomic scale plasmon enables tuning its resonance to match with quantum processes, so that the enhancement is maximized. Deterministic preparation of a stable atomic scale structure poses an obstacle to the practical application. While beyond-single-adatom stable atomic scale structures are alternatives and have been attempted experimentally in literature, they are still too tiny to be mass fabricated. An effective approach to circumvent the difficulty is to work with highly doped semiconductors instead of metals. Except the working frequency is shifted to the infrared range, almost all the intriguing features are realizable with structures correspondingly scaled to decananometers, which are readily feasible with the current nanofabrication techniques. The extremely localized field is retained considering the mode volume normalized with respect to the wavelength cubed. Dramatic plasmon modulation is also within reach owing to the relatively large screening length in semiconductors. Adapting the non-classical theories to semiconductors is however in need for the new expedition.

## Concluding remarks

Atomic scale plasmonics is a natural development of nanoplasmonics. In this contribution we briefly discussed the current status of this frontier, the challenges and opportunities, and envisioned the future development. Although the picocavity and extremely localized mode
represent actively studied topics of atomic scale plasmonics, more diverse plasmonic field localized in atomic scale structures are worthwhile to look into. Plasmonic field squeezed into atomic scale slots was introduced to explain transient “flare” emission observed in plasmonic cavities [268]. “Picophotonics” was also coined for the study on the microscopic optical waves in the atomic lattice [269].

**Funding.** National Natural Science Foundation of China (12274160).

# 11. Nonlocal limitations to light-matter interactions

## CHRISTOS TSERKEZIS

### Current status

Probably the biggest appeal of plasmonic nanostructures is related to the impact they can have on enhancing light-matter interactions. Individual plasmonic nanoparticles (NPs), and even more so plasmonic cavities, provide a dramatic increase in the local density of optical states (LDOS), allowing to control the properties of quantum emitters (QEs) placed near them. This potential manifests in both the so-called weak- and strong-coupling regimes. Under weak coupling conditions, the plasmonic nanostructure is treated as a dielectric environment that alters both the excitation and the emission rate of QEs via the Purcell effect, allowing e.g. to control molecular fluorescence, at the single-molecule level [270–272]. In the strong-coupling regime, the QE hybridizes with the optical modes of the cavity, and the two components enter a coherent and reversible exchange of energy, that results in the formation of half-light-half-matter entities, called polaritons [180,273,274]. In the case of cavities, the LRA predicts a diverging field enhancement as the gap between individual components shrinks [275]. This unphysical situation reflects the structural weaknesses of LRA. Nonlocality imposes more realistic limits to the performance of nanostructures, especially plasmonic ones; the spatial extent of induced charges, together with the increased Landau damping as a result of collisions with the nanostructure boundaries, regularize the plasmonic near-field enhancement, as discussed already in Ref. [83]. Such limitations have been verified ever since for individual NPs [76] and coupled structures [276], probed with a variety of sources including swift electrons and point dipoles [82,193,277].

**Weak coupling.** The hydrodynamic Drude model (HDM) was employed already in the 1980s to study the emission of a point dipole near a flat metallic surface [278], while NPs were more actively considered in the late 2000s [279]. Detailed studies of fluorescence enhancement followed a few years later [124,280–284], and Purcell factors have been calculated for NPs of various shapes [285,286], or for stratified media [287], employing both HDM and the more recent generalized nonlocal optical response (GNOR) theory [64]. These studies highlighted that fluorescence can be affected in two additive ways by nonlocality. At first, as sketched in Fig. 15(a), the QE has to make the transition from the ground, $S_0$, to the excited, $S_2$ state, a process that is governed by the $p_d \cdot E$ interaction, where $p_d$ is the dipole moment of the QE, while $E$ is the total field experienced by it, maximized at plasmonic resonances. The frequency shifts predicted by HDM become thus important for the efficient excitation of the QE, as they detune plasmonic resonances from the transition energy. In the second process, after being excited to $S_2$, the QE quickly relaxes to a metastable state $S_1$ at a close energy, from which it then decays spontaneously to the ground state after a characteristic time $\tau$. When placed near plasmonic NPs, the QE can couple to the radiative dipolar mode, but also to higher-order multipoles, which are non-radiative. Depending on the QE-NP separation, one of the two couplings prevails, as quantified by the ratio of radiative over total decay rates. An optimum QE-NP separation can thus be identified, and it is strongly modified once nonlocal damping is taken into account (see e.g. Fig. 15(b) for GNOR calculations; Landau damping dominates the process, much of the energy is absorbed by the NP, and the QE needs to be placed farther away from the NP surface for an efficient radiative coupling
[281]). This situation is slightly different when the QE is placed within plasmonic cavities, because the hybrid plasmon modes they sustain are more radiative. Different types of plasmonic cavities have been considered within the context of HDM and GNOR [124], with the main conclusion being that nonlocality tends to regularize the maximum enhancement that is possible for decreasing gaps, bringing theoretical predictions closer to experimental measurements [272].

**Strong coupling.** In the strong-coupling regime, single QEs, or QE collections supporting collective excitonic modes, interact with a cavity mode, leading to the formation of hybrid polaritonic states [288–290]; the energy difference between the hybrid modes is known as the Rabi splitting. The main question regarding nonlocality is: can it affect the coupling strength, e.g. by weakening the interacting EM fields, or by detuning the plasmonic mode? In a first attempt to answer these questions, Ref. [291] considered plasmon–exciton coupled systems in a spherical core–shell geometry. Taking into account that in any experiment the plasmonic resonances will be characterized by the measured frequencies and linewidths, with any nonlocal effects intrinsically included, anticrossing spectra were plotted renormalized to the plasmonic resonances predicted within each model, for different nonlocal models, showing an almost perfect superposition of anticrossings; both hybrid modes shift, as a result of nonlocality, by the same amount, leading to negligible differences in the width of the Rabi splitting. The damping introduced by the nonlocal models is small enough to ensure that no significant broadening occurs, and the main strong-coupling criterion of the Rabi splitting being larger than the average damping rate of the uncoupled modes [273] still holds. A similar response has been observed with electron-beam excitation as well [292]. With a different starting point, Ref. [282] placed QEs with adjustable oscillator strength inside different plasmonic cavities, and observed that the threshold for entering the strong-coupling regime would shift depending on the applied model (see Fig. 15(c), for three

<figure><img src="image_17.png" /><figcaption>**Fig. 15.** Nonlocality in QE-plasmon interactions. (a) Fluorescence near a plasmonic NP; a QE with dipole moment $p_d$ is placed at distance $r_d$ from the surface of the NP, and both are illuminated by an external plane wave. The QE goes from the ground state $S_0$ to an excited state $S_2$, quickly relaxes to a state $S_1$, and then decays to the ground state after time $\tau$. (b) Modification of fluorescence enhancement near gold (Au) or silver (Ag) nanospheres and nanoshells as a function of distance $r_d$, within LRA, HDM, and GNOR. Figure adapted with permission from Ref. [281] (Copyright © 2016 Royal Society of Chemistry). (c) Modification of strong coupling of a single QE in different Ag plasmonic cavities, as a function of the dipole oscillator strength $f$, within LRA, HDM, and GNOR. Reprinted (adapted) with permission from Ref. [282] (Copyright © 2017 American Chemical Society).</figcaption></figure>
different cavities, within LRA, HDM, and GNOR). The main difference between Refs. [291] and [282] is exactly in the oscillator strength of the excitonic material; while the former assumed a collective emitter layer that would always have a large enough dipole moment to enter the strong-coupling regime, the latter explored the conditions for a single QE to achieve this.

## Challenges and opportunities

The above discussion is more relevant in noble-metal NPs and cavities, with high work functions and negligible electron spill-out. The picture is expected to change if spill-out dominates, as in alkali metals. This has only recently started being considered; for instance, in Ref. [293] the self-interaction Green's function was calculated for sodium nanospheres and nanosphere dimers, both within the surface-response formalism (SRF) based on Feibelman parameters, and within time-dependent density functional theory (TDDFT), showing that the results deviate considerably from the predictions of LRA, already for QE-NP distances or interparticle gaps as large as 2.5 nm. Interestingly, the semiclassical SRF treatment was in excellent agreement with TDDFT down to ~ 0.7 nm. Nevertheless, deviations were observed for even smaller distances; these were attributed to the fact that the standard SRF is based on *local* Feibelman parameters, usually obtained for flat metal–dielectric interfaces. Eventually the curvature of the structure needs to be included [230,294]; this was done by the introduction of *nonlocal* Feibelman parameters [253], which is, however, still a rather challenging task, both in terms of obtaining those parameters, and implementing them in computational electrodynamics. NP dimers were studied in Ref. [295], focusing on the role of the electron tunneling that comes as a result of spill-out in two sodium nanospheres separated by sub-nm gaps. This study, performed with the quantum hydrodynamic (QHT) model [61], also known as self-consistent hydrodynamic Drude model (SC-HDM) [60]—a generalization of HDM that relaxes the hard-wall boundary condition—predicted considerable fluorescence quenching inside the cavity, with results that are qualitatively and quantitatively different from those of LRA and HDM; similarly to the case of noble-metal NPs, nonlocality seems to provide the upper limits to the system's performance. In the strong-coupling regime, Ref. [296] found that spill-out, and the plasmonic redshifts that accompany it, do not drastically change Rabi splittings. In a similar study for flat sodium-air interfaces, strong coupling was indeed retained, also accompanied by a drastic effect on the QE dynamics, which acquires a strong non-Markovian character [297]. On the other hand, Ref. [298] observed that the oscillator-strength threshold for achieving strong coupling changes within each model. In particular, using QHT to account for electron spill-out, a considerable increase in the oscillator threshold was predicted, also related to the fact that coupling to higher-order multipoles occurs already for larger QE-NP distances due to the larger effective NP size. Finally, it is worth mentioning a similar investigation related to the Mollow triplets in QE-NL coupled systems [299]; using GNOR as the nonlocal model, Ref. [300] observed that the side bands become stronger and with narrower linewidths, compared to the predictions of LRA.

Historically, nonlocality in plasmonics has been associated with “bad news”, i.e. limitations to the performance of an architecture or a device. While such studies are tremendously important from a fundamental point of view, to establish the physical mechanisms involved, they also tend to become slightly unattractive, by predicting less impressive results than those of LRA. There have only been few situations where nonlocality promised to improve a performance, as e.g. when it was shown to smooth the surface roughness of a metallic tip [71]. In the context of light-matter interactions, a recent study [59] showed that nonlocality and electron spill-out have the potential to enhance dipole-forbidden transitions [301] or plasmon-mediated energy transfer, eventually showing that nonlocality is not always a limiting factor, but a possibility of harnessing it for practical applications does exist.
## Future developments to address challenges

From the above discussion, it becomes clear that nonlocality can have a strong, and sometimes unpredictable impact on light-matter interactions. Nevertheless, current literature is largely restricted to plasmonics, and almost exclusively on a theoretical level. Furthermore, most studies have only focused on exploring the role of corrections in the response of the plasmonic component, while the QEs themselves are not considered. There is therefore still plenty of room for future explorations. First of all, it will be important to consider, both in theory and in experiment, the impact of nonlocality on other types of light-matter interacting architectures. For instance, Ref. [302] has made a first effort to address semiconductors, focusing on the radiative and nonradiative decay rates of multi-subband plasmons in quantum wells. Magnetoplasmonic systems and Landau polaritons were studied in Ref. [94], which managed to set quantitative limits to the miniaturization of polaritonic devices and to the enhancement of polariton gaps. Nonlocality in phonon-polariton resonators was studied in Ref. [303], albeit only on the level of far-field extinction and reflection properties.

At the same time, one is justified to ask, to what extent it is appropriate to continue refining the description of photonic nanostructures, while retaining the simple point-dipole approximation for the QE. It is almost inevitable that the next generation of nonlocal studies will have to introduce both a quantum-electrodynamic treatment of the problem [38,304], and also more elaborate models for the QEs, either entering the domain of computational quantum chemistry [305], or with semiclassical models as suggested e.g. in Ref. [298]. Strong coupling of a QE with a film described by a hydrodynamic or a Lindhard permittivity has been explored in Ref. [306], but retaining again some level of approximation. And this only makes sense, as quantization of the “cavity modes” for open, dispersive and lossy cavities, can be tricky [274]. Macroscopic cavity quantum electrodynamics (QED) appears to be the most promising route right now, but obviously no treatment can be seen as the unique, conclusive answer yet. This will require a combination of further theoretical work, possibly with the introduction of new numerical methods for mesoscopic electrodynamics [307,308], but more importantly, realization of fine experiments that can probe the small modal shifts and damping involved, and unambiguously attribute them to nonlocality.

## Concluding remarks

Nonlocality plays an important role in light-matter interactions, by introducing detuning and additional damping that can modify the coupling strength between plasmonic environments and QEs. It is, however, important to broaden the horizon of such studies by also departing from the point-dipole description of the QE, while more efficient computational methods might still be needed. On the other hand, nonlocality is much less explored when QE environments other than plasmonic are concerned, and this promises to be an active direction in the near future.

**Funding.** Danmarks Grundforskningsfond (DNRF165).

# 12. Quantum surface-response functions for nanoplasmonics: Feibelman parameters

THOMAS CHRISTENSEN ID, WEI YAN & YI YANG ID

## Current status

The study of quantum-response effects in nanoplasmonics has often centered on the incorporation of specific physical effects otherwise omitted in the piecewise-constant, local-response outset that underlies conventional (or “classical”) electromagnetic modeling (CEM). Focus has especially centered on the incorporation of one of the three key omissions in CEM, namely (1) bulk nonlocal response (i.e., the finite range of the dielectric tensor $\varepsilon(\mathbf{r}, \mathbf{r}')$ with respect to $|\mathbf{r} - \mathbf{r}'|$), (2) the