Dispersive surface-response formalism to address nonlocality in extreme plasmonic field confinement

The surface-response formalism (SRF), where quantum surface-response corrections are incorporated into the classical electromagnetic theory via the Feibelman parameters, serves to address quantum effects in the optical response of metallic nanostructures. So far, the Feibelman parameters have been typically obtained from many-body calculations performed in the long-wavelength approximation, which neglects the nonlocality of the optical response in the direction parallel to the metal–dielectric interface, thus preventing to address the optical response of systems with extreme field confinement. To improve this approach, we introduce a dispersive SRF based on a general Feibelman parameter d⊥(ω, k‖), which is a function of both the excitation frequency, ω, and the wavenumber parallel to the planarmetal surface, k‖. An explicit comparison with time-dependent density functional theory (TDDFT) results shows that the dispersive SRF correctly describes the *Corresponding authors: Antton Babaze, Materials Physics Center CSIC-UPV/EHU, Paseo Manuel de Lardizabal 5, 20018, Donostia-San Sebastián, Spain; Donostia International Physics Center DIPC, Paseo Manuel de Lardizabal 4, 20018, Donostia-San Sebastián, Spain; andDepartment of Electricity and Electronics, FCT-ZTF, UPV-EHU, 48080 Bilbao, Spain, E-mail: anttonbabaze@dipc.org. https://orcid.org/0000-0002-9775062X; Javier Aizpurua, Materials Physics Center CSIC-UPV/EHU, Paseo Manuel de Lardizabal 5, 20018, Donostia-San Sebastián, Spain; Donostia International Physics Center DIPC, Paseo Manuel de Lardizabal 4, 20018, Donostia-San Sebastián, Spain, E-mail: aizpurua@ehu.eus. https://orcid. org/0000-0002-1444-7589; and Andrei G. Borisov, Institut des Sciences Moléculaires d’Orsay, UMR 8214 CNRS-Université Paris-Saclay, Bât. 520, 91405 Orsay Cedex, France, E-mail: andrei.borissov@universite-paris-saclay.fr. https://orcid.org/ 0000-0003-0819-5028 Tomáš Neuman, Institut des Sciences Moléculaires d’Orsay, UMR 8214 CNRS-Université Paris-Saclay, Bât. 520, 91405 Orsay Cedex, France, E-mail: neuman@fzu.cz. https://orcid.org/0000-0003-3089-5805 Ruben Esteban, Materials Physics Center CSIC-UPV/EHU, Paseo Manuel de Lardizabal 5, 20018, Donostia-San Sebastián, Spain; and Donostia International Physics Center DIPC, Paseo Manuel de Lardizabal 4, 20018, Donostia-San Sebastián, Spain, E-mail: ruben.esteban@ehu.eus. https:// orcid.org/0000-0002-9175-2878 plasmonic response of planar and nonplanar systems featuring extreme field confinement. This work thus significantly extends the applicability range of the SRF, contributing to the development of computationally efficient semiclassical descriptions of light–matter interaction that capture quantum effects.

Among the semiclassical approaches, the surfaceresponse formalism (SRF) [23,[56][57][58][59] has prompted great interest in the nanophotonics community over the last few years. The SRF is based on the theoretical framework proposed by Peter Feibelman in the eighties [60] to describe the interaction of electromagnetic radiation with planar metal surfaces. This formalism introduces quantum surface-response corrections into the classical Maxwell's theory [61] through the so-called Feibelman parameters. As a consequence, quantum effects such as surface-enabled Landau damping [62][63][64], spill-out of the induced electron density [34,51,54], and nonlocal dynamical screening [27,28,57,65] are accounted for. The Feibelman parameters, commonly denoted as d ⊥ and d ‖ , are complex-valued functions characterizing the dynamically induced charges (d ⊥ ) and the parallel-to-the-surface component of the induced currents (d ‖ ) at the metal-dielectric interface [32].
In the literature, the practical implementation of the SRF within the Feibelman theory typically involves two approximations. First, d ‖ is assumed to be zero, which is a reasonable approach given that d ‖ strictly vanishes for charge-neutral planar interfaces [32,56,66]. Second, the nonlocality of the optical response in the direction parallel to the metal-dielectric interface is neglected. Within this long-wavelength approximation, it is considered that the characteristic wavenumber k ‖ parallel to the metal surface is negligible, so that the Feibelman parameter d ⊥ ( ) is exclusively a function of the excitation frequency, . Considering the long-wavelength limit with k ‖ = 0 reduces the computational effort in obtaining d ⊥ ( ) from quantum calculations, and simplifies the implementation of the SRF in existing numerical tools that solve Maxwell's equations in nanophotonics, as employed in a number of recent studies [59,[67][68][69][70][71][72][73][74]. In what follows, we refer to d ⊥ ( ) as the nondispersive Feibelman parameter.
Using the nondispersive Feibelman parameter is reasonable when the nonlocality of the optical response in the direction perpendicular to the metal surface dominates, i.e., when the characteristic length scale of the optical field variation along the surface is relatively large, 2 ∕k ‖ ≫ F , where F is the Fermi wavelength of the metal. This is the case of e.g. typical individual nanoparticles subjected to plane-wave illumination. However, for localized probes in proximity of metal nanoparticle surfaces, where high-order plasmonic modes can be excited, the long-wavelength approximation to d ⊥ ( ) is compromised [75]. Indeed, in such situations, multipolar plasmon modes characterized by localized surface charges that rapidly vary along the nanoparticle surface are involved. This is analogous to exciting surface plasmons with large transverse wavenumber k ‖ at planar metal-dielectric interfaces, and therefore requires going beyond the long-wavelength limit of d ⊥ .
In this work, we demonstrate that considering the nonlocality of the optical response in the direction parallel to the metal surface significantly improves the performance of the SRF. While earlier theory addresses the role of nonlocality by introducing an ad hoc dipole layer [58], we base our approach on the implementation of a dispersive Feibelman parameter d ⊥ ( , k ‖ ) that is a function of and k ‖ , providing a direct connection of the theory to the microscopic quantum aspects of screening. The numerical values of this parameter are obtained from linear-response frequencydomain TDDFT calculations using a planar free-electron metal slab. We then show that the dispersive SRF based on the same set of d ⊥ ( , k ‖ ) can be used to describe the plasmon modes and optical response not only for planar surfaces, but also for nonplanar geometries. In particular, the method proposed in this work is relevant for situations where plasmon resonances with large transverse momenta (small plasmon wavelength) are excited. We test the validity of the dispersive SRF by using quantum many-body TDDFT simulations as a reference for a planar metal surface ( Figure 1a) as well as for canonical plasmonic nanostructures such as cylindrical metallic nanowires, small spherical metal nanoparticles, and plasmonic gaps as exemplified by a spherical dimer (Figure 1b). In all the systems studied, the dispersive SRF reproduces the TDDFT results, demonstrating its validity to capture the effects linked to nonlocality of dynamical screening along the metal surface. The theoretical framework proposed here provides significant conceptual advances toward the implementation of efficient semiclassical approaches that adequately account for quantum effects in the optical response of plasmonic systems.

Methodology
To obtain the dispersive Feibelman parameter, d ⊥ ( , k ‖ ), we employ linear-response frequency-domain TDDFT calculations for a planar metal slab. The adiabatic local-density approximation (ALDA) with the exchange-correlation potential given by Gunnarsson and Lunqvist [76] is used. The metal slab is infinite along the (x, y)-plane and has a finite thickness L in the z-direction. The metal surface is located in the (x, y, z = 0) plane, with vacuum at z > 0 (see Figure 1a).
The electronic structure of the slab is described within the jellium model of free-electron metals [77][78][79][80][81] using a Wigner-Seitz radius r s = 4 a 0 that corresponds to sodium (a 0 = 0.0529 nm is the Bohr radius). The bulk plasma frequency is p = 5.89 eV, and the nonretarded n(z, , k ‖ ) relative to the geometrical metal surface, and the imaginary part is related to surface-enabled Landau damping. The dispersive Feibelman parameter d ⊥ ( , k ‖ ) is a function of the excitation frequency, , and the wavenumber parallel to the planar metal surface, k ‖ = 2 ∕ ‖ (with ‖ the wavelength of the nonretarded surface plasmon). The nondispersive SRF assumes that k ‖ = 0 ( ‖ → ∞), so that d ⊥ is a function of exclusively. n 0 (z) represents the equilibrium electron density. (b) Sketch of nonplanar structures considered in this work. Left: cylindrical metallic nanowire with radius R c , infinite along the z-axis, with the azimuth angle. The system possesses rotational and translational symmetry with respect to the z-axis. The equivalence between k ‖ and |m|∕R c is used to implement the dispersive SRF in the cylindrical nanowire, k m ‖ = |m|∕R c (with m the magnetic quantum number of a particular plasmonic mode). Right: a dimer of two identical spherical metallic nanoparticles of radius a. An effective wavenumber k ‖ = √ ( + 1)∕a is assigned to implement the dispersive SRF in the spherical nanoparticles (with the angular momentum quantum number of the plasmonic mode). All the structures are considered to be made of sodium, characterized using a Wigner-Seitz radius r s = 4 a 0 (with a 0 = 0.0529 nm the Bohr radius).
surface plasmon frequency for k ‖ = 0 is SP = p ∕ √ 2 = 4.16 eV. The choice of this material allows for a direct comparison between our results and those obtained in recent works using the nondispersive SRF [67,70,75]. Further details can be found in ref. [32] and in Section S1 of Supplementary Material. The methodology employed in this work is well suited for simple free-electron metals. For noble metals, interband transitions involving bound d-band electrons contribute to the screening of the electromagnetic field and thus affect the optical response. The present approach could be extended to these situations by mimicking the effect of d-electrons with the inclusion of a polarizable background [65,71,82].
Since the nonlocality of the metal response in the direction parallel to the surface is important in geometries confining fields at typical length scales comparable to that of the Fermi wavelength of electrons [58,75] (i.e., within the nm range), we use the nonretarded approximation within the linear-response theory. The optical response of the metal slab to a time-dependent external potential V ext (z, r ‖ , t) = V ext (z, , k ‖ )e ik ‖ r ‖ −i t oscillating at frequency is determined from the induced electron density n(z, r ‖ , t) = n(z, , k ‖ )e ik ‖ r ‖ −i t , where r ‖ is the component of the position vector parallel to the metal surface, and k ‖ is the parallel wavevector Within the many-body formalism, n(z, , k ‖ ) is given by where ( is the many-body response function (see further details in Section S1 of Supplementary Material). Here we use the Fourier component of the external potential given by which exponentially decays within the metal (z < 0). Note that, up to a multiplying factor, Eq. (2) represents the asymptotic behavior of the potential created by any charge distribution located in vacuum far from the surface, since the plane-wave decomposition of such potential can be obtained from the Green's function G(r − r ′ ) corresponding to the potential of a point charge located at r ′ : In this respect, within the linear-response regime the induced electron density n(z, , k ‖ ) given by Eq. (1) is independent of the specific form of the external potential that excites the system.
is therefore a characteristic of the surface response inherent to the specific metal and the surrounding material [32]. In contrast to the nondispersive parameter [59,67,[70][71][72][73], the dispersive d ⊥ ( , k ‖ ) introduced in Eq. (4) depends on the wavenumber k ‖ parallel to the surface, and thus explicitly accounts for the nonlocality of the many-body response and screening in this direction. Other definitions different from Eq. (4) that express the Feibelman parameter in terms of the electromagnetic fields and nonlocal dielectric functions have also been used. A discussion on the equivalence between Eq. (4) and other definitions can be found in refs. [32,60,83]. As follows from Eq. (2), the z-extension of the external potential within the metal slab increases with decreasing k ‖ . In our calculations (where k ‖ ≠ 0), we use a slab with thickness L = 870a 0 (≈ 46 nm), large enough to avoid the interaction between the charges induced at the opposite surfaces of the slab. For k ‖ = 0 the slab geometry is not applicable, and we rely for this case on the nondispersive parameter d ⊥ ( ) calculated by Christensen et al. [67] for the semi-infinite metal characterized by the same r s as that used here. ). The imaginary part of the Feibelman parameter is related to the surface loss function and thus to the energy absorption by electronic excitations at the surface [32,84]. These electronic excitations mainly consist of the Bennett plasmon [42,85], characterized by the presence of induced electron density of opposite sign across a single metal-vacuum interface (see Figure 2 in ref. [51]), and electron-hole pair excitations involved in the Landau damping mechanism of plasmons [62,86,87]. A detailed discussion of the frequency dependence of the Feibelman parameter within the long-wavelength limit k ‖ = 0 can be found in ref. [32]. Therefore, we only discuss briefly the main features of

Dispersive Feibelman parameter
For small k ‖ , the TDDFT results in Figure 2 are consistent with earlier findings for k ‖ = 0 [67], with d ⊥ ( , k ‖ ) featuring a complex Lorentzian-like resonance at ∼ 4.6 eV.
This resonance is associated with the excitation of the Bennett plasmon at ∼ 0.8 p [42,85,88]. With increasing k ‖ , the Bennett plasmon resonance broadens, blueshifts, and gradually disappears.
} is then dominated for large k ‖ by electron-hole pair excitations that extend over the whole frequency range addressed here. At low frequencies, ≲ 1 eV, the optical response of the system is close to that of an ideal metal, and becomes nearly independent of k ‖ and . In this situation, lim 1.2 a 0 determines the position of the electrostatic image plane of the free-electron metal with respect to its jellium edge.
Importantly, in the frequency range ∼ 2 − 4.5 eV relevant for plasmon excitations in different structures (see below), the following trends can be observed with increas- } decreases and even changes the sign from positive to negative. This indicates that the centroid of the induced electron density shifts inwards the metal for large k ‖ .
In general, the results in Figure 2 demonstrate a strong dependence of the dispersive Feibelman parameter d ⊥ ( , k ‖ ) on k ‖ due to nonlocal dynamical screening. Thus, one can expect that the dispersive and nondispersive SRF yield very different results in situations that involve a rapid variation of plasmon-induced charges along the metalvacuum interface. In the following sections, we analyze the limitations of the nondispersive SRF and illustrate the good performance of the dispersive SRF using the example of canonical plasmonic nanostructures described semiclassically within the SRF and quantum mechanically within TDDFT. Our results show that the dispersive SRF correctly captures quantum surface effects within a computationally efficient semiclassical framework.

Nonlocal dynamical screening of a planar metal surface
As a first application of the SRF, we consider the plasmon propagating at the interface between a planar metal surface and vacuum. The surface plasmon frequency dispersion s (k ‖ ) with parallel wavenumber can be obtained within TDDFT from the surface loss function Im { g( , k ‖ )) } , which reveals the rate at which the electron-hole pair excitations at a given frequency are produced at the surface (including both single-particle and collective excitations) [32,84].
Here, the so-called surface response function g( , k ‖ ) is given by where the induced electron density n(z, , k ‖ ) is obtained within TDDFT considering the slab geometry introduced in the previous section. Figure 3a shows the frequency dependence of the surface loss function, Im  [32,84,85,[92][93][94][95]. In particular, with increasing k ‖ , the frequency of the surface plasmon resonance (defined from the maximum of Im { g( , k ‖ ) } ) first redshifts at low k ‖ and, after reaching a minimum value, it continuously blueshifts for higher k ‖ . The width of the plasmon resonance increases with increasing k ‖ because of the enhancement of surface-enabled Landau damping [32,84,86,91,93].
The frequency dispersion s (k ‖ ) featured by the plasmon resonances in Figure 3a is intimately related to the k ‖ -dependence of the position of the centroid of the induced electron density (given by Re{d ⊥ ( , k ‖ )}). Indeed, for low k ‖ the induced electron density is shifted outward from the geometrical metal surface (Re{d ⊥ ( , k ‖ )} > 0 in Figure 2a), i.e., toward the region of lower electron density, which produces the redshift of s (k ‖ ) with respect to SP = p ∕ √ 2 = 4.16 eV (see Eq. (6) below) [32]. In contrast, for high k ‖ , the plasmon-induced electron density is shifted inward the geometrical surface (Re{d ⊥ ( , k ‖ )} < 0), thus producing a blueshift of s (k ‖ ) relative to SP .
We continue the quantitative analysis of the plasmon dispersion in Figure 3b, where we show the resonant plasmon frequency s (k ‖ ) obtained from TDDFT calculations (blue solid line) and compare it with dispersive (red dashed line) and nondispersive (green dotted line) SRF results. To find the s (k ‖ ) dispersion using the SRF based on the Feibelman parameters, we describe the dielectric function of the metal with a Drude model, and self-consistently solve the following transcendental equation [32,60,65,84,[96][97][98], where SP = p ∕ √ 2 = 4.16 eV, and |k ‖ Re{d ⊥ ( s , k ‖ )}| ≪ 1 is assumed (see details in Section S6 of Supplementary Material). Throughout this work, we use Re{d ⊥ ( , k ‖ )} as obtained from the TDDFT results reported in Figure 2a for the dispersive SRF, while for the nondispersive model we use d ⊥ ( ) as calculated by Christensen et al. [67]. Figure 3b shows that the dispersive SRF and TDDFT results are in good agreement within the broad range of k ‖ values considered here. The dispersive SRF reproduces the redshift followed by blueshift of s (k ‖ ) with increasing k ‖ as predicted by TDDFT. For high k ‖ values, the quantitative agreement between the two sets of results worsens. This may be attributed to the short lifetime of plasmon modes at high k ‖ . The resonant structures in the surface loss function become broad and asymmetric so that the resonant frequencies are ill defined. Additionally, higherorder terms of the induced electron density near the surface beyond the leading dipolar contribution accounted for with d ⊥ ( , k ‖ ) may gain in importance in such situations. In sheer contrast, the nondispersive SRF fails to describe the entire s (k ‖ ) dependence, since it predicts a continuous redshift of s (k ‖ ) with increasing k ‖ . Therefore, the nondispersive SRF is only accurate for small values of k ‖ ≲ 0.06 a −1 0 , where the long-wavelength approximation is well justified. Figure 3b thus demonstrates that using the dispersive Feibelman parameter d ⊥ ( , k ‖ ) within the SRF allows us to correctly capture the main nonlocal effect associated with the dependence of the dynamical screening on k ‖ .

Localized multipolar plasmon resonances in a cylindrical nanowire
We study next the applicability of the dispersive SRF to address nonlocal effects in the optical response and plasmon resonances of nonplanar plasmonic nanostructures used in a variety of applications in nanophotonics. To this end, we first consider localized multipolar plasmons (LMP) sustained by cylindrical metallic nanowires of radii R c = 75 − 150 a 0 (≈ 4 − 8 nm), infinite along the z-axis (see Figure 1b). The nanowire is described within the same jellium model as that used for the metal slab. The reference quantum calculations are performed using real-time ALDA-TDDFT as implemented in prior works [27,28,[99][100][101][102] (see Section S2 of Supplementary Material for further details).
The excitation of LMPs evolving in the (x, y)-plane determines the optical response of the nanowire for an incident electromagnetic wave polarized such that its electric field is perpendicular to the nanowire z-axis. The system is translationally invariant along the z-axis, and under the illumination conditions considered here, the wavevector describing the plasmon propagation along the nanowire is To support the equivalence between k ‖ and |m|∕R c , we show in Figure 3b the plasmon frequency m of cylinders with different radii as a function of k m ‖ = |m|∕R c . The results obtained for R c = 150 a 0 (blue filled circles) and R c = 100 a 0 (blue hollow circles) fall on a universal curve very close Owing to the equivalence k m ‖ = |m|∕R c , the dispersive SRF can be straightforwardly applied to describe the nonlocal optical response of the cylindrical nanowire. The nonretarded multipolar polarizability m ( ) of a cylindrical nanowire per unit length along the z-axis can be obtained within the SRF from (see Section S4 in Supplementary Material) where ( ) is the dielectric function of the metal described here with a Drude model (bulk plasma frequency p = 5.89 eV and intrinsic damping parameter p = 0.1 eV [75]).
In Eq. (7), d ⊥ is the dispersive Feibelman parameter calculated in this work for a planar metal surface ( Figure 2). For the sake of comparison, we also calcu- Note that the LMP resonance frequencies within the SRF can be obtained from the poles of Eq. (7), which are exactly given by Eq. (6) using k ‖ = k m ‖ = |m|∕R c . Figure 4b and c show the result of Im{ m ( )} for a Na nanowire, as obtained from Eq. (7) using the nondispersive SRF and the dispersive SRF, respectively. Again, the nondispersive SRF (panel b, green filled curves) erroneously predicts a continuous redshift of the localized multipolar plasmon of order m in the nanowire with increasing m. On the other hand, the dispersive SRF (panel c, red filled curves) correctly reproduces the general behavior of plasmon resonances when compared to the TDDFT results (redshift with increasing m for low m, and blueshift for larger m). The results shown in Figure 4 thus confirm that the dispersive SRF can be used to model the nonlocal optical response of a cylindrical nanowire of small radius where the field localization and oscillation in space can be extreme.

Individual spherical nanoparticles and nanoparticle dimers
We finally show that the dispersive SRF employed in this work can also be used to account for quantum surface effects and nonlocality in the optical response of plasmonic nanostructures with finite extension along the three spatial dimensions. As a first canonical plasmonic nanostructure, we consider a spherical Na metal nanoparticle (MNP) of radius a = 65.83 a 0 (≈ 3.5 nm). On the one hand, this MNP size is sufficiently small so that we can perform TDDFT calculations as a reference and, on the other hand; it is large enough to ensure well-developed plasmon modes [104,105]. Due to the small size of the nanoparticle, quantum effects and, in particular, nonlocality is expected to clearly affect the optical response [23,30,70,75,106]. The spectrum of LMPs in spherical MNPS is determined by the nonretarded multipolar polarizability ( ), which is given within the SRF by [56,70] We compare in Figure 5a the imaginary part of ( ), Im{ ( )}, obtained for different values of from TDDFT (blue lines), nondispersive SRF (left, green filled curves) and dispersive SRF (right, red filled curves). Details on the TDDFT calculations are given in ref. [75]. The LMPs of order appear as resonances in Im{ } associated with an increased absorption.
The TDDFT results of Im{ ( )} in Figure 5a show that the LMP of order continuously blueshifts with increasing in the considered range = 1 − 10. The resonance broadens as increases due to the enhancement of surfaceenabled Landau damping. In this situation, the coupling of the plasmon to electron-hole pair excitations and quantum finite-size effects lead to the emergence of discrete spectral features [42,62] (notice that these features merge into a smooth resonant profile in the metallic nanowires studied above because of their larger radius). As already stated in ref. [75], the nondispersive SRF accurately reproduces the TDDFT results for low values of ∼ 1 − 4, but it cannot describe the correct behavior for larger angular momenta. Indeed, within the nondispersive SRF the plasmon resonances start to redshift with increasing for ≥ 5 (left-hand side panel in Figure 5a). On the other hand, when accounting for the dependence of the Feibelman parameter d ⊥ on the effective wavenumber k ‖ = √ ( + 1)∕a, the dispersive SRF correctly captures the frequency blueshift of LMP resonances in Im{ ( )} (righthand side panel in Figure 5a). Although slight quantitative differences are present for large multipolar order = 7 − 10 (partially linked to the non-Lorentzian shape of Im{ ( )} within TDDFT due to finite-size effects), qualitatively there is a good agreement between TDDFT and the dispersive SRF results over the entire range of -s considered here, and the improvement over the nondispersive SFR is evident.
Finally, we address another canonical plasmonic system: a dimer of spherical MNPs forming a plasmonic gap. Specifically, we study the case of a point-dipole emitter located at the center of the gap of size D formed by two identical spherical MNPs of radius a = 65.83 a 0 (≈ 3.5 nm).
The emitter is oriented along the axis of the MNPs dimer (z-axis). The gap separation distance, D, is in the nanometer scale, and thus a strong impact of the nonlocality on the optical response of the system is expected. We obtain the decay-rate enhancement (Purcell factor) and the change of the resonant frequency (Lamb shift) of the emitter due to its interaction with the MNPs dimer by calculating the imaginary and real parts of the electric field in the middle of the gap in response to the excitation of the pointdipole emitter [107]. The TDDFT results are obtained from the time evolution of the electron density in response to an impulsive potential created by a point dipole. The SRF results are obtained from the solution of Laplace's equation by considering the electromagnetic coupling between different multipoles of the two individual MNPs. For the dispersive model we assign again k ‖ = √ ( + 1)∕a. Further details on the TDDFT and SRF calculations are provided in ref. [75]. The Lamb shift is calculated considering a transition dipole moment = 0.1 e nm (with e the electron charge). Figure 5b shows the Lamb shift Δ 0 (left) and Purcell factor F P (right) obtained for a gap separation D = 2.33 nm, as calculated within the three models tested here (TDDFT, dispersive SRF, and nondispersive SRF). The three approximations show qualitatively good agreement. For this relatively large gap, the excitation of low-LMP resonances dominates the response of the MNPs dimer [75], validating the long-wavelength approximation behind the nondispersive SRF results for D = 2.33 nm. Note, nonetheless, that the results obtained within the dispersive SRF are more accurate when comparing to those obtained within TDDFT.
The significant improvement introduced by the dispersive SRF to describe the emitter-dimer electromagnetic interaction is more evident when considering a smaller gap, where field localization and thus higher multipolar activation occurs. Figure 5c shows the Lamb shift and Purcell factor obtained for a gap separation D = 1.06 nm. In this situation, because of the larger spatial confinement of the charges induced at the metallic surfaces across the left-hand side panel shows the comparison between TDDFT (blue lines) and nondispersive SRF (green filled curves) results, whereas the right-hand side panel shows the comparison between TDDFT and dispersive SRF (red filled curves) results. Each spectrum of Im{ ( )} is normalized to 1 at its maximum, and vertically offset for clarity. (b, c) Lamb shift Δ 0 (left-hand side panels) and Purcell factor F P (right-hand side panels) obtained for a point-dipole emitter located at the center of a dimer of spherical MNPs of radius a = 65.83 a 0 (≈ 3.5 nm). The dipole is oriented along the dimer axis, and its transition dipole moment is = 0.1 e nm (with e the electron charge). In (b), the gap separation is D = 2.33 nm. In (c), D = 1.06 nm.
gap, plasmon modes with high multipolar order become important. These high-modes show overlapping resonant frequencies and thus contribute to a single broad peak (referred to as pseudomode [108,109]) at ∼ 3.4 eV as revealed by the TDDFT calculations (blue lines). Since the nondispersive model does not describe accurately highmultipolar modes for the individual MNP (Figure 5a), it clearly fails to reproduce the frequency and the width of the plasmon pseudomode obtained within TDDFT. As a consequence, the nondispersive SRF strongly overestimates the overlap between the pseudomode and the bonding dipolar plasmon resonance at ∼ 2.75 − 3 eV (green dotted lines), as observed in the Purcell factor and the Lamb shift. On the other hand, the dispersive SRF (red dashed lines) provides accurate results even for this small gap separation, further illustrating its ability to account for nonlocality in situations where plasmon-induced charges characterized by a rapid variation in the direction parallel to the metal surface are excited.

Summary and conclusions
We have introduced a dispersive surface-response formalism (SRF) that incorporates quantum surface effects and nonlocality into the optical response of plasmonic nanostructures via the so-called Feibelman parameters. While the nondispersive SRF typically implemented in the literature is based on the Feibelman parameter d ⊥ ( ) obtained within the long-wavelength limit (k ‖ = 0), the dispersive SRF proposed here is based on a dispersive Feibelman parameter d ⊥ ( , k ‖ ) that explicitly depends on the wavenumber parallel to the metal surface, k ‖ . We obtain the values of d ⊥ ( , k ‖ ) from quantum many-body calculations for a planar metal-vacuum interface, and use a recent formulation of the electromagnetic boundary conditions [59] to account for the nonlocality of the dynamical screening in the direction parallel to the metal surface in various nanostructures.
Using TDDFT calculations as a reference, we have demonstrated that the dispersive SRF is significantly more accurate than the nondispersive SRF in describing plasmonic systems characterized by extremely confined induced fields. We have first shown that the dispersive SRF correctly describes the nonlocal dynamical screening of a planar metal surface and provides with the correct surface plasmon frequency dispersion with k ‖ . Further, we have demonstrated that the Feibelman parameter d ⊥ ( , k ‖ ) calculated in this work using a planar metal surface can also be used to address the nonlocal optical response of nonplanar nanostructures relevant in plasmonics. As examples, we have considered infinite cylindrical nanowires, spherical metallic nanoparticles, and nanoparticle dimers forming gaps described within the free-electron metal approximation. The symmetry of these nanostructures allowed us to introduce geometry-dependent effective parallel wavenumbers, k m ‖ and k ‖ , in order to efficiently implement the dispersive SRF.
In all these systems, the dispersive SRF reproduces the TDDFT results that naturally incorporate quantum effects such as nonlocal dynamical screening, surface-enabled Landau damping, and the finite spatial extension of the induced electron density at the metal surface. This work thus establishes a milestone for the development of a theoretical model that captures quantum surface effects and nonlocality in extreme situations of field localization, while keeping the numerical efficiency and easy implementation of classical (local) electromagnetic theories [59,110,111]. Among others, our findings can be important for describing the optical response of metallic nanostructures interacting with fast electrons, nanoantennas coupled to molecules in close proximity, metallic nanoparticles with geometrical features characterized by small radii of curvature (such as picocavities), or nanoparticle ensembles with narrow gaps, all of them relevant situations in practical configurations of nowadays Nanophotonics.
Author contributions: All authors have accepted responsibility for the entire content of this manuscript and approved its submission. Conflict of interest statement: Authors state no conflicts of interest. Data availability: The datasets generated and/or analysed during the current study are available from the corresponding author upon reasonable request.