Plasmon–emitter interactions at the nanoscale

Plasmon–emitter interactions are of central importance in modern nanoplasmonics and are generally maximal at short emitter–surface separations. However, when the separation falls below 10–20 nm, the classical theory deteriorates progressively due to its neglect of quantum effects such as nonlocality, electronic spill-out, and Landau damping. Here we show how this neglect can be remedied in a unified theoretical treatment of mesoscopic electrodynamics incorporating Feibelman \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d$$\end{document}d-parameters. Our approach incorporates nonclassical resonance shifts and surface-enabled Landau damping—a nonlocal damping effect—which have a dramatic impact on the amplitude and spectral distribution of plasmon–emitter interactions. We consider a broad array of plasmon–emitter interactions ranging from dipolar and multipolar spontaneous emission enhancement, to plasmon-assisted energy transfer and enhancement of two-photon transitions. The formalism gives a complete account of both plasmons and plasmon–emitter interactions at the nanoscale, constituting a simple yet rigorous platform to include nonclassical effects in plasmon-enabled nanophotonic phenomena.

. Illustration of a planar dielectric-metal interface, defined by the z = 0 plane. In the asymptotic regions, i.e., for |z| ≥ |z1,2|, the cladding dielectric is characterized by a dielectric constant, d, whereas the electromagnetic properties of the metal are accounted for by its (local) bulk dielectric function, m(ω). In the surface region, quantum mechanical features lead to a nonuniform equilibrium electron density, n0(z), and an induced charge density, ρ(z) (that arises due to the system's response to an external perturbation).
As indicated schematically in the figure, the Feibelman d ⊥ parameter can be regarded as the position of the centroid of the induced charge density with respect to the positive background edge.
where it is implicit that the interface runs along the x-direction and the z-direction lies normal to it, as depicted in Supplementary Figure 1. It is apparent from Eqs. (1) that d ⊥ is a length scale defined by the position of the centroid of the induced charge density (see Supplementary Figure 1) whereas d can be regarded as the (integrated) imbalance between the uniform ionic background and n 0 [3]. Hence, the latter (to leading-order) is identically zero for charge-neutral interfaces [4]. Once the microscopic variations of the quantum mechanical charge and current densities across the dielectric-metal boundary (here defined by the so-called jellium edge) are computed, e.g., using TDDFT, one can derive the corresponding nonclassical reflection and transmission coefficients. These quantities entail all the knowledge about the system's near-and far-field electromagnetic response; therefore, they provide us with, for instance, the structure's plasmon resonances, Purcell enhancement, and other nanophotonic phenomena. Wigner-Seitz radius, rs = {2, 4} (in units of the Bohr radius), and for silver (Ag). We note that a jellium with rs = {2, 4} is resemblant of aluminum (Al) and sodium (Na), respectively. For silver (Ag), we normalize the frequencies to the screened plasma frequency, ω * p 3.81 eV. The circles represent data obtained from TDDFT calculations by Christensen et al. [3] and the solid lines are the corresponding interpolations. The dielectric medium interfacing the metal is assumed to have d = 1.  3 . It should be noted that the d-parameters are microscopic surface-response functions that obey Kramers-Kronig relations and specific sum-rules [11], just like the traditional bulk response functions do [12]. The TDDFT data presented here is for air-metal (or vacuum-metal) interfaces (i.e., for d = 1); we note that the specific values of the d-parameters depend on the dielectric medium next to the metal [4; 13], since it changes the screened interaction.
Lastly, we emphasize that although in Supplementary Figure 2 we have d 0 for silver (last panel in Supplementary Figure 2) in this case silver is of course still charge-neutral. The finiteness of d here is merely a consequence of the methodology used to incorporate d-band screening in the TDDFT calculation [3] (see Refs. 3 and 4 for additional details).

SUPPLEMENTARY NOTE 2. MODIFIED BOUNDARY CONDITIONS WITH FEIBELMAN d-PARAMETER SURFACE TERMS
Here we summarize the recently introduced approach to incorporate the Feibelman d-parameters as boundary conditions applied to the conventional macroscopic Maxwell equations [14]. The boundary conditions facilitate a convenient and straightforward d-parameter generalization of conventional scattering coefficients-e.g., the reflection and transmission coefficients of a planar interface, or the Mie coefficients of a spherical interface-thereby providing the necessary ingredients to investigate the influence of quantum mechanical effects on the surface plasmon polariton (SPP) dispersion, electromagnetic local density of states (LDOS), etc.
In the absence of external charges and currents, the classical boundary conditions of macroscopic electrodynamics at an interface separating two media read [15]: wheren is a unit vector normal to the interface ∂Ω (pointing from medium 1 to medium 2). These govern the tangential (or parallel) components of the electromagnetic fields, imposing their continuity across the interface. However, as can be clearly seen in a nonretarded formulation [3; 4], the introduction of the Feibelman d-parameters renders the above boundary conditions inapplicable. Hence, a revision set of boundary conditions which reflect the presence of the d-parameters is required. The incorporation can be inferred by considering their impact as initially encoded through (normally-oriented) surface polarization and (tangentially-oriented) current terms, P(r) = π(r)δ(r − r ∂Ω ) and J(r) = K(r)δ(r − r ∂Ω ), respectively, which are nonzero only at the interface, i.e., at r = r ∂Ω . The surface polarization and current densities are each proportional to a d-parameter and are driven by the discontinuities of the field components [14] π After some manipulations, these surface polarizations and currents can be self-consistently absorbed into a revision of the conventional boundary conditions 2 [14] n × ( where ∇ is the surface nabla operator ∇ ≡ (1 −nn T )∇. The boundary conditions can be equivalently stated via the parallel and perpendicular components (with respect to the interface) of the fields as Clearly, the incorporation of Feibelman d-parameters introduces discontinuities in the parallel components of the electric and magnetic fields. The magnitude of these discontinuities is naturally proportional to the Feibelman d-parameters. 2 While the boundary conditions on the tangential parts of E and H is fully sufficient to uniquely couple solutions across the interface, the complementary set of boundary conditions-on the normal components of B and D-is occasionally more convenient. The boundary condition for the normal components of B, is unchanged from its classical counterpart, i.e.,n · (B 2 − B 1 ) = 0; for the normal components of D, it isn · (

SUPPLEMENTARY NOTE 3. NONCLASSICAL REFLECTION AND TRANSMISSION COEFFICIENTS FOR A PLANAR DIELECTRIC-METAL INTERFACE
In possession of the boundary conditions personified by Eqs. (5), we now have all the necessary ingredients to determine the nonclassical equivalents of Fresnel's reflection and transmission coefficients for a single dielectric-metal interface. In what follows, we adopt the same coordinate system as in Supplementary Figure 1.

Nonclassical reflection and transmission coefficients for TM waves (p-polarization)
We seek transverse magnetic (TM) solutions, in the dielectric and metal half-spaces, of the form: in the dielectric half-space (z > 0), and in the metal half-space (z < 0). Here k z, j = j k 2 0 − q 2 for j ∈ {d, m}. The explicit form of the amplitudes E x/z, j (z) readily follows from Maxwell's equations. Making use of the modified boundary conditions in Eqs. (5), one obtains the following linear system: After some algebraic manipulations, one finally arrives at the nonclassical Feibelman reflection and transmission coefficients for p-polarization (or TM polarization), which after changing to the Fresnel-like coefficients (i.e., defined in terms of the electric field), that is, r tm ≡ r (H) p and t tm ≡ d k z,m m k z,d t (H) p , read [9; 13]: where only terms up to linear order in qd ⊥, have been included (for consistency with the assumptions made when computing the d-parameters [3]), and where k z, j = j k 2 Scattering coefficients for TM polarization in the nonretarded limit. In the nonretarded limit (c → ∞) one has k z, j → iq and the reflection and transmission amplitudes reduce to their nonretarded forms [3] The derivation of the nonclassical reflection and transmission coefficient for TE solutions (s-polarization) follows the same lines as the one detailed above for the TM polarization. In particular, the former are given by

SUPPLEMENTARY NOTE 4. NONCLASSICAL SURFACE PLASMON POLARITONS: PLANAR INTERFACE
Equipped with the reflection coefficient for TM-polarization in Eq. (9a), the resulting surface plasmonpolariton (SPP) dispersion can be readily obtained from its poles. Hence, the implicit condition for the SPP dispersion relation including quantum surface corrections via Feibelman d-parameters is given by We stress that here m ≡ m (ω) is still the local bulk dielectric function of the metal, and therefore the nonclassical corrections (e.g., spill-out, nonlocality, electron-hole pair generation) are encoded and enter via the d-parameters alone. Naturally, it is apparent from Eq. (12) that the well-know classical result, Nonretarded surface plasmon dispersion for a planar interface. Taking the nonretarded limit of Eq. (12) leads to the following condition for the nonretarded surface plasmon dispersion: Clearly, contrary to the corresponding classical case, i.e., m = − d , the surface plasmon spectrum is not dispersionless anymore, even in the nonretarded regime. Indeed, one now may also write instead: Yet, writing an explicit solution in closed-form for the surface plasmon frequency as a function of wavevector, i.e., ω vs. q, is not possible due to the implicit frequency dependence of the d-parameters. Nevertheless, starting from Eq. (13), one may work towards a perturbative solution for the SPP (complex) eigenfrequencies.
In particular, within the nonretarded limit, and for a Drude-like dielectric function m (ω) = ∞ − ω 2 p /(ω 2 + iωγ), approximate expressions for the real and imaginary parts of the first-order spectral correction, ω = ω cl + ω (1) + O ω − ω cl 2 , can be obtained. Then, assuming that Re ω γ and qd ⊥, 1, we have where ω (0) sp ≡ Re ω cl = ω p / √ ∞ + d is the classical nonretarded surface plasmon frequency. Also, sp is a pole-like approximation (necessary here owing to the explicit frequency dependence of the d-parameters). Equations (15b) underscore that the direction (magnitude) of the nonclassical shift of the SPP resonance depends on the sign (magnitude) of the real part of d eff ≡ d ⊥ − d ; and, that the nonclassical SPP resonance broadening is proportional to the imaginary part of d eff .
Notably, for the instructive case of a homogeneous three-dimensional electron gas (3DEG) (i.e., in which ∞ = 1 and d = 0) with d = 1, Eqs. (15b) take a particularly compact form, reading Re ω ω (0) which reproduces the result obtained by Feibelman [9], and, more recently, by Christensen et al. [3]. At this point, the reader may appreciate that, despite being only approximate, Eqs. (15b) deliver a clear message: to lowest order, the direction of the frequency shift due to nonlocal quantum surface phenomena ultimately depends on the sign of Re d ⊥ − d . Specifically, for a fixed frequency, these nonclassical corrections lead to redshift (blueshift) of the surface plasmon resonance-relative to the classical case-if Re d ⊥ − d > 0 (< 0). Furthermore, for charge-neutral surfaces so that d = 0, and recalling the interpretation of d ⊥ as the centroid of the induced charge density, then one may associate the positiveness or negativeness of d ⊥ , respectively, to a "spill-out" or a "spill-in" of the induced electron density with respect to the fixed ionic background of the metal.
SPP dispersion for rs = 2 jellium and Ag. Figure 3 depicts the quantum mechanical SPP dispersion diagram for a planar interface between air and a jellium metal with r s = 2 (representative of Al), together with the SPP dispersion for a flat air-silver interface. The corresponding classical spectral properties are also SPPs in a air-Ag interface. We model silver with a bulk dielectric function m(ω) = ∞(ω) − ω 2 p /(ω 2 + iωγ) and use the d-parameters data outlined in Supplementary Figure 2 to account for the quantum mechanical effects discussed in the text. The background permittivity of silver is obtained from Johnson and Christy's experimental data [16] via ∞(ω) = exp(ω) + ω 2 p /(ω 2 + iωγ). We use ωp = 9.02 eV and γ = 22 meV (which fits well the experimental data). The screened plasma frequency is ω * p = 3.81 eV.
shown to facilitate the comparison and interpretation of the nonclassical results. As the figure indicates, we predict a nonclassical redshift (blueshift) of the SPP resonance, consistent with the fact that Re d eff > 0 (Re d eff < 0) in the relevant frequency ranges. Thus, the observed redshift of the SPP frequency for jellium metals is a consistent with the of the quantum mechanical spill-out of the induced electron density beyond the classically-forbidden jellium edge [4; 17-20]. Conversely, for noble metals like silver, Re d eff < 0 around ω cl due to d-band screening [3; 4]. Therefore, and contrary to the jellium case, here we predict a blueshift of the SPP resonances, in line with experimental observations [20][21][22][23][24]. This can be readily seen upon inspection of Supplementary Figure 3b which shows a blueshift of the SPP eigenfrequencies for large wavevectors.

SUPPLEMENTARY NOTE 5. GENERALIZED MIE THEORY WITH FEIBELMAN d-PARAMETERS: SPHERICAL METAL PARTICLES
In what follows we develop a generalized Mie theory that includes nonlocal and quantum mechanical effects incorporated via Feibelman d-parameters. As in the planar case, the possession of such quantum Mie coefficients is all that is required in order to readily determine the sphere's localized plasmon resonances, scattering cross-sections, Purcell enhancements, etc.

Theoretical framework
Wave equation in spherically symmetric systems. In a linear, isotropic and homogeneous medium the electric and magnetic fields must satisfy their respective vector wave equations, where k j ≡ √ j k 0 . The electric and magnetic fields in spherical coordinates can be constructed in terms of vector harmonics, M ν (r) and N ν (r), which in turn are defined in terms of a pilot wave c and a generating scalar function ψ ν (r) that obeys the (scalar) Helmholtz equation ∇ 2 ψ ν (r) + k 2 j ψ ν (r) = 0 [25]. Strictly speaking, the pilot vector c should be a constant vector. However, and rather fortunately, in spherical coordinates and in this case alone, it turns out that we can take c = rr; this enables us to construct a tangential solution M ν (r) (though N ν (r) is not purely normal) [26]. This is desirable because it allows us to associate each of the vector harmonics M ν (r) and N ν (r) with TE and TM waves, respectively. Therefore, the vector harmonics M ν (r) and N ν (r)-which are two independent solutions of the vector wave equation in spherical coordinates-may be written as [25; 26] M e/o lm (r) ≡ ∇ ×r rψ e/o lm (r), (18a) where the generation function is and where the radial part z l (k j r) is either a spherical Bessel or Hankel function of the first kind, j l (k j r) and h (1) l (k j r), respectively representing incoming and outgoing waves. Furthermore, the "quantum numbers" l and m are integers with values in the range l ∈ [1, ∞[ and m ∈ [0, l], and P m l denote the associated Legendre polynomials. The explicit form of the solenoidal (i.e., divergence-free) vector waves M e/o lm (r) and N e/o lm (r), obtained via Eqs. (18), reads 3 where ρ j ≡ k j r and µ ≡ cos θ (for the sake of brevity of notation), and the prime denotes differentiation with respect to ρ j .
The spherical vector harmonics M e/o lm (r) and N e/o lm (r) may be regarded as fundamental solutions of the vector wave equation in spherical coordinates. Hence, the electric and magnetic fields can be constructed from an expansion in terms of such vector waves.
Incident plane-wave in terms of vector spherical harmonics. Let us assume an incident plane-wave traveling along the positive z-direction and polarized alongx, namely 4 We now want to cast the incident field, not as in Eq. (20), but rather in terms of the vector spherical waves introduced earlier in the text, that is where the expansion coefficients A tm,σ lm and A te,σ lm follow from [25] A tm,σ Carrying out the explicit calculations and exploiting the orthogonality relations of the sine and cosines functions, one finds that A tm,o lm = 0 and A tm,e lm ∝ δ m1 , and, similarly, that A te,e lm = 0 and A te,o lm ∝ δ m1 . The calculation of these coefficients then leads to where the extra superscript "[d]" indicates that the radial function is z l (ρ d ) ≡ j l (ρ d ) because the impinging field has to remain finite at the origin, and we have also defined E l 0 = i l 2l+1 l(l+1) E 0 . From here, the corresponding incident magnetic field follows from Faraday's law, yielding where the properties of the vector spherical harmonics have been used [25; 26].
Equations (23) and (24) thereby reflect the expansion of the incident plane-wave in terms of vector spherical waves.
Internal and scattered fields. In possession of the incident field expansion, Eqs. (23) and (24), we now express in a similar fashion the internal and scattered fields using vector spherical waves. In this spirit, we write the fields inside the sphere as where now the superscript "[m]" highlights that the radial function is z l (ρ m ) ≡ j l (ρ m ) because it needs to remain finite at the origin.
Outside the sphere, both spherical Bessel functions j l and y l are permitted and therefore we shall use the spherical Hankel function of the first kind for the radial part of the generating function. Specifically, we choose z l (ρ d ) ≡ h (1) l (ρ d ) owing to its appropriate asymptotic behavior at large ρ d [25], i.e., corresponding to an outgoing spherical wave (which is what we expect for the asymptotic scattered field). Within this reasoning, we therefore express the scattered fields as where this time the superscript "[s]" indicates that the radial function is h (1) l (ρ d ) as we have just mentioned.
The amplitudes a tm l , a te l , c tm l , and c te l , are the so-called Mie coefficients, of the scattered and internal (i.e. transmitted) kind, in TM and TE flavors, respectively. These can be determined once the adequate boundary conditions are imposed (at the sphere's surface).

Nonclassical Mie coefficients
The Mie coefficients essentially entail all the relevant physics describing the electromagnetic response of a sphere of arbitrary size and material constitution. In the following, we derive closed-form expressions for the generalized (i.e., nonclassical) Mie coefficients within the formalism of Feibelman d-parameters. In order to proceed, we need to invoke the appropriate boundary conditions, reading (cf. SUPPLEMENT-ARY NOTE 2), in spherical coordinates, with Ω = {θ, φ} denoting the angular components (i.e., the components that are tangential to the sphere's surface). Equations (27), after some algebra, produce the following set of equations: with x j = k j R. Additionally, the Riccati-Bessel functions Ψ l (x) = x j l (x) and ξ l (x) = xh (1) l (x) have been employed, and the prime denotes the functions' derivatives with respect to their argument. Lastly, for shorthand notation, we have also introduced the dimensionless quantities (29a) Solving the system of equations posed by Eqs. (28), we obtain the mesoscopic generalization of the Mie coefficients incorporating the Feibelman d-parameters: for the coefficients associated with the scattered outgoing wave, and for the coefficients associated with the internal fields. As in the planar case, we have only kept terms up to linear order in the d-parameters. We have also added tm and te as superscripts in the Mie coefficients to highlight the character of the associated waves.
Notice that the quantum mechanical Mie coefficients a te l and c te l , which describe TE waves, are independent of the perpendicular Feibelman d-parameter, d ⊥ . This should not really constitute a surprise, since for a TE wave the electric field is purely tangential to the particle's surface and therefore has no perpendicular component capable of inducing a displacement of the charge density along the perpendicular to the sphere's surface. On the other hand, the Mie coefficients a tm l and c tm l , associated to TM waves, and thus the electric field in general possesses both radial and tangential components, and consequently both d ⊥ and d emerge in these coefficients.
To the best of our knowledge, such quantum mechanical Mie coefficients with the incorporation of nonlocal and quantum surface corrections via Feibelman d-parameters, given above by Eqs. (30) and (31), have not been derived before. Thus, these results could open new perspectives for the investigation of nonlocal and quantum effects in the context of light-matter interactions in nanophotonics with metal spheres. With this in mind, they constitute a springboard to study phenomena beyond the classical regime, namely the influence of quantum nonlocal effects on the cross-sections, localized surface plasmon (LSP) resonances, electromagnetic LDOS, etc, for spherical metal particles of arbitrary size.
Nonretarded multipolar polarizability and localized surface plasmons. In the nonretarded limit, the outgoing Mie coefficients reduce to a single quantity: the nonretarded multipolar polarizability, α l . Since in the nonretarded regime we have ωR/c → 0, then one may perform small-argument expansions of j l (z), h (1) l (z), and of the derivatives of the Riccati-Bessel functions, thereby yielding the following relation between the Mie coefficient a tm l and the multipolar polarizability: where Then, the nonretarded localized surface plasmon resonances immediately stem from the poles of α l , reading In the limit of vanishing d-parameters, Eqs. (33) and (34) appropriately reduce to their well-known classical equivalents, that is, α cl l = 4πR 2l+1 m − d m + d (l + 1)/l and m + l + 1 l d = 0, respectively.

Optical response of plasmonic spheres
A particularly alluring feature of the theoretical formalism set forth here, it that the mathematical structure of the usual expressions for the cross-sections, LDOS, etc, is unchanged; that is, the nonclassical optical response is obtained using the very same textbook expressions, simply by replacing the classical Mie coefficients by the quantum mechanical Mie coefficients presented in Eqs. (30). Hence, for instance, the extinction cross-section for a metal sphere illuminated by a monochromatic plane-wave is therefore given by For relatively small metal spheres, the optical response is dominated by a series of peaks originating from the excitation of LSP resonances. The latter occur at frequencies for which a tm l has pole, as illustrated in Supplementary Figure 4. A LSP resonance is said to be of electric dipole character for l = 1, of electric quadrupole character for l = 2, and so on. It is also apparent from the figure that the dipole resonance (l = 1) contributes the most for the cross-section, with the successively higher-order resonances becoming increasingly negligible (notice the logarithmic scale). The TE Mie coefficients a te l are essentially featureless and their contribution for the extinction cross-section is several orders of magnitude smaller than the one associated with TM Mie coefficients a tm l of the same order. The impact of quantum nonlocal effects, i.e., a redshift owing to spill-out (Re d ⊥ > 0) and broadening due to Landau damping, can be clearly observed in Supplementary Figure 4 as well. It is also interesting to note that the TE Mie coefficients are unchanged by such effects, which can be understood by the fact that here d = 0 [recall Eq. (30b)]. In the nonretarded regime, the TM Mie coefficient a tm l reduces to lim ωR/c→0 and σ ext is completely characterized by the nonretarded multipolar polarizability, α l .
In Supplementary Figure 5 we plot the normalized extinction cross-sections, σ ext /(πR 2 ), for r s = 4 jellium spheres of different radii. As expected, although the predictions from classical electrodynamics suffice for large spheres, they rapidly decline for the spheres of smaller radii. Moreover, the figure also underscores the importance of retardation effects for metal spheres that are small, but not too small. We end the discussion of the far-field optical response of plasmonic spheres by comparing the extinction cross-section of a r s = 2 jellium sphere with the one of a silver nanoparticle. The results are depicted in Supplementary  Figure 6, and show a nonclassical redshift for the former and a blueshift for the latter, consistent with the fact that Re d ⊥ > 0 for jellium metals and that Re d eff < 0 for silver [4].   7) of the main text [the deviation here is due to the pole-approximation performed owing to the explicit spectral dependence of the d-parameters-see discussion after the aforementioned Eqs. (7)]. Figure 8 shows the relative nonclassical shifts and corresponding resonance widths for the two plasmonic structures considered in this work, a planar dielectric-metal interface and metal spheres or arbitrary radii. We emphasize that in Supplementary Figure 8 we employ the full solution-that is, including both retardation and the frequency dependence of the d-parameters-in order to obtain quantitatively accurate predictions. The usefulness of the approximate first-order spectral corrections given by Eqs. (7)  main text hinges upon their ability to provide physical intuition and to elucidate in a qualitative fashion the scaling (with either q or R −1 ) of the plasmonic spectral features.

Planar dielectric-metal interface
The (electric) local density of states (LDOS) experienced by an emitter with dipole moment µ and located at a height h above a planar metal interface (defined by z = 0) is given by [27; 28] where we have normalized the LDOS to its value in a homogeneous medium with relative permittivity d (i.e., in the absence of the metal half-space). In the above, r tm ≡ r tm (u, ω) and r te ≡ r te (u, ω) are the reflection coefficients for a planar interface (cf. Supplementary Note SUPPLEMENTARY NOTE 3), and we have introduced the dimensionless variable u ≡ q/k d . Considering separately the cases of an ideal dipole emitter with perpendicular and parallel dipole moments (with the respect to the interface), we have and respectively. The corresponding orientation-averaged LDOS is then ρ E = 1 3 ρ E ⊥ + 2 3 ρ E , which can be particularly useful when dealing with an ensemble of randomly-oriented dipoles (e.g., in dyes).
It should now be clear that the work in Supplementary Note SUPPLEMENTARY NOTE 3 pays off also for the calculation of the LDOS taking quantum mechanical effects into account, since one can simply employ the corresponding nonclassical reflection coefficients, derived in Supplementary Note SUPPLE-MENTARY NOTE 3, in Eqs. (37) and (38) (with their classical equivalents naturally defining the classical LDOS). The results of such computations are summarized in Supplementary Figure 9.

Spherical metal sphere
The expression for the LDOS felt by an ideal electric dipole in the neighborhood of a spherical nanoparticle can be derived by following the same steps that lead to the LDOS in the planar case [29]. In particular, the LDOS for a radial (⊥) and for a tangent ( ) dipole near a metal sphere is given by 5 [29; 30] and where we have defined y = k d (R + h) = x d (1 + h/R) for the sake of clarity. Clearly, by comparing Eqs. (38) and (39) one realizes that the Mie coefficients in the latter play a similar role as the reflection coefficients in the former, and that the integration over in-plane momentum now appears as a sum over angular momenta instead, reflecting the (spherical) symmetry of the problem. Furthermore, and as before, the nonclassical LDOS possesses the same mathematical structure as its classical version. Hence, the quantum mechanical LDOS is immediately obtained upon substituting the classical Mie coefficients by their quantum mechanical counterparts. The outcome of Eqs. (39) is shown in Supplementary Figure 10.

SUPPLEMENTARY NOTE 8. ENHANCEMENT OF MULTIPOLAR DECAY RATES: INFLUENCE OF QUANTUM NONLOCAL EFFECTS
Atom-field Hamiltonian. In the context of quantum electrodynamics (QED) we consider atom-field interactions whose dynamics are governed by the following nonrelativistic Hamiltonian: where H a is the atom Hamiltonian (or, more generically, of an emitter, e.g., a quantum dot, etc), H EM is the Hamiltonian of the electromagnetic (EM) field, and H int stands for the Hamiltonian governing the atom-field interaction. As written above, the matter Hamiltonian H a is treated at the Hartree-Fock level (consisting of a kinetic part together with the Coulomb potential), that is, electron-electron interactions, spin-orbit coupling, etc, will be neglected in what follows. Naturally, in any case these may be readily incorporated here. Moreover, the Hamiltonian akin to the EM field is just the like the one of a typical quantum harmonic oscillator, where the operator f † j (r, ω) ( f j (r, ω)) creates (annihilates) a field excitation at position r with frequency ω and oriented along the direction j. For the part of the Hamiltonian describing the interaction between the atom and the field, H int , we take the so-called minimal-coupling interaction Hamiltonian 6 . Finally, A(r) is the vector potential associated with the EM field, from which both the electric and magnetic fields stem, in particular, E(r) = − ∂ ∂t A(r) and H(r) = µ −1 0 ∇ × A(r) [15; 27]. Notice that we have implicitly assumed that the scalar potential is identically zero (in our gauge). It should be noted that such vector and scalar potentials are simply useful mathematical entities/abstractions and, unlike E(r) and B(r), do not represent per se real physical quantities.
Below, we will be dealing with (linear) dispersive and dissipative media and therefore the common canonical quantization based on modal expansion is unsuitable. Ergo, we employ instead the formalism of macroscopic quantum electrodynamics [31; 32].
From here, and noting that A(r) = dω A(r, ω) + H.c., we now write the components of the vector potential operator featuring in the interaction Hamiltonian, Eq. (40d), as where summation over repeated indices is implicitly assumed (Einstein's summation notation).
Fermi's Golden Rule and first-order processes. In what follows we consider electronic first-order processes alone. This means, for instance, that we neglect the second term of the minimal-coupling Hamiltonian in Eq. (40d) [i.e., the term proportional to A 2 ≡ A · A]. Moreover, for definiteness, we assume that the emitter is a Hydrogen atom. Hence, the atom-field interaction Hamiltonian then becomes We want to determine the transition rate between an initial (excited) state |i = |e, 0 ≡ |e ⊗ |0 and a final (not necessarily the lowest) state | f = |g, 1 x jω 0 = f † j (x, ω 0 ) |g, 0 . The transition rate for such a first-order process calculated within first-order perturbation theory follows from Fermi's Golden Rule The matrix element entering in Eq. (51) reads 7 where the commutation relations of the bosonic field operators have been used. Carrying out a similar calculation for i| H int | f , and making use of the identity in Eq. (44c), one may write (taking the wavefunctions as real) Γ = 2π 2 e 2 π 0 m 2 c 2 dr dr ψ e (r )ψ e (r) Im G i j (r, r ; ω 0 )(p i ψ g (r))(p * j ψ g (r )).
Now, the electromagnetic Green's function admits an analytical expression of the form (for z > 0) Since deviations of the LDOS from classicality owing to quantum surface corrections incur at relatively short emitter-metal distances, i.e., at relatively large wavevectors, we here assume the nonretarded limit, i.e., q k 0 (except for the reflection coefficient, which is still assumed to be the retarded version in order to obtain the plasmon pole in accurately), such that k z → iq and ↔ M s → 0 as c → ∞: where withˆ q denoting the polarization vector defined viaˆ q = (q + iẑ)/ It is useful to rewrite Eq. (57) as which allows the evaluation of the matrix element with the emitter at the origin.
For simplicity, in what follows we restrict ourselves to electric multipolar transitions whose initial and final states posses m f = m i = 0. Then, we may choose q along the x-direction, and, writing the momentum operator as p = −i ∇, we obtain where we have introduced the dimensionless variable u ≡ q/k 0 . Since for an arbitrary nth electric transition, En, the matrix element in Eq. (59) is dominated by only one of the terms coming from the expansion of the exponential, 8 we get e| e ik 0 ux−k 0 u(z−z 0 )ˆ q · ∇ |g 2 ∝ u 2(n−1) .
Finally, this leads to the following expression of the transition rate of a given hydrogenic En transition: where a b denotes the Bohr radius and where the quantity Ξ originates from the matrix element, To extrapolate the argument-range of the d ⊥ parameter beyond ω p , we perform a multi-Lorentzian fitting of previous data [3] that cover d ⊥ across 0 < ω < ω p . By enforcing sum-rules and appropriate asymptotics [35; 36], we anticipate that the fit maintains a reasonable validity in the extrapolating domain ω > ω p .

Fitting scheme
We seek a fit of the jellium data of Ref. 3 to a sum of Lorentzians: with f n ∈ C, γ n ∈ R, and ω n ∈ R. In addition, we impose the following basic requirements on each of the oscillator parameters: Amplitudes f n . To strictly ensure the large-ω behavior Im d ⊥ (ω ω p ) ∼ ω −3 , we explicitly enforce that N n=1 Im f n = 0 by setting Im f N = − N−1 n=1 Im f n .
Resonance widths γ n . To ensure causality, γ n is restricted to positive real values, i.e., γ n > 0. This is enforced via a penalty term.
Resonance frequencies ω n . In the expectation that all resonant behavior occurs below ω p , we restrict ω 2 n ∈ [0, ω 2 p ]. This is enforced via penalty terms. To strictly ensure that Im d ⊥ (ω → 0) = 0, we explicitly enforce that N n=1 Im f n /ω 2 n = 0 by setting ω 2 In addition, we enforce the following sum-rules 9 [35] via penalty terms: Finally, we enforce that Im d ⊥ (ω > ω p ) < 0 to ensure that power is absorbed rather than gained at these frequencies [the opposite is the case for frequencies below ω p since the surface power absorption is proportional to ε(ω)d ⊥ (ω)].
The fitting itself is performed by minimizing the absolute square of residuals plus the penalty-terms noted above. The minimization itself is performed by the Julia package Optim.jl [37], using the Nealder-Mead algorithm.

Results
We've carried out this fitting procedure for r s = 4. We obtain a reasonable fit and small relative deviation from the sum-rules with five oscillator terms. The results are summarized in Supplementary Table 1. 10   9 There are additional sum-rules, that we nevertheless do not enforce, e.g.,

Energy transfer between two dipoles
The total energy transfer rate from a donor (D) to an acceptor (A), characterized by a dipole moment µ D and µ A , respectively, follows from [27; 38-40] where the kernel w ET (r D , r A , ω) is the energy transfer amplitude, given by Notice that all the information about the nanophotonic environment is embodied in the Green's dyadic. Therefore, if the Green's function is either known or can be calculated for the system under consideration, then the rate of energy transfer can be determined by evaluating Eq. (66). For the half-plane and sphere plasmonic structures, the nonclassical energy transfer rates are then straightforwardly obtained by simply substituting the classical scattering coefficients by their corresponding nonclassical ones (namely, Eqs. (9) and (11) for the planar interface, and Eqs. (30)) and (31) for the sphere).
At this point it is also instructive to note that, contrary to the computation of the LDOS where the (imaginary part of the) Green's function is evaluated at the same point (r = r = r 0 ), the rate of energy transfer depends on the (absolute square of the) Green's function taken simultaneously at the positions of the donor and of the acceptor. Furthermore, for the calculation of the total (broadband) energy transfer rate one has to perform an integration over frequency-albeit typically restricted to the overlap between the emission spectrum of the donor, f em D (ω), and the absorption spectrum of the acceptor, f abs A (ω). Figure 11. Illustration of two dipoles in the vicinity of a planar metal surface. The dipolesurface separation is given by z D(A) for the donor (acceptor), and its dipole moment is characterized by µ D(A) . Without loss of generality, we choose our coordinate system so that R = |xD − xA|.

R z
We now consider two electric point-dipoles (e.g., each a generic two-level system, or dye molecules, or quantum dots, etc) above a planar metal surface, as illustrated in Supplementary Figure 11. Therefore, according to Eqs. (66) the enhancement (with respect to the energy transfer rate in same medium but without the interface) of the energy transfer rate between the donor and the acceptor can be calculated as where ↔ G 0 (r D , r A ; ω) and ↔ G ref (r D , r A ; ω) are the Green's dyadics associated with the homogeneous medium and with the reflected part owing to the metal surface, correspondingly. Here,n A andn D are unit vectors denoting, respectively, the orientation of the acceptor and donor dipole moments.
Example I: two vertical dipoles. Forn D =n A =ẑ, the relevant element of the ↔ G ref (r D , r A ; ω) dyadic is given by with k z ≡ (k 2 0 − q 2 ) 1/2 , and where J n is the nth order Bessel function of the first kind. This is the example highlighted in the main text.
Example II: two horizontal dipoles. Forn D =n A =x, the relevant entry of the (69)

Second-order processes: two-photon emission
The theory of two-photon emission in arbitrary macroscopic media was derived in Ref. 33 and 41. Here, we summarize the main results, and defer to these references for additional details. From the standpoint of macroscopic quantum electrodynamics, two-photon emission can be represented as a transition from an initial state |i ≡ |e, 0 (an excited matter state e and zero photons) to a final state | f ≡ |g, {rωk, r ω k } (a lower-energy matter state g and two photons at positions r and r , frequencies ω and ω , and orientations k and k ). The energy difference between matter states is ω 0 ≡ ω e − ω g . The rate of two-photon emission is then governed by Fermi's Golden rule at second-order in quantum electrodynamics, which takes the form Γ TPE = 2π 2 1 2 dr dr dω dω where |i 1 are intermediate states containing both the atom and field degrees of freedom. The sum is understood to be a sum over discrete degrees of freedom and an integral over continuous ones. The factor of 1/2 comes from the fact that when we sum over all {rωk, r ω k } pairs, each pair of excitations appears twice. As we will calculate two-photon emission within the dipole approximation, V = −d · E is the dipole Hamiltonian. After a series of manipulations, it is found that the second-order differential decay rate dΓ TPE /dω-i.e., the decay rate per unit frequency, defined such that Γ TPE = ω 0 0 dω (dΓ TPE /dω)-is given by: Im G ri (r, r, ω) Im G s j (r, r, ω 0 − ω)T i j (ω)T * rs (ω), where T i j is an atomic structure factor given by where x gn j ≡ g|x j |n is the transition matrix-element between matter states |g and |n along the jth Cartesian direction. We note that this derivation makes no assumption about the Green's function ↔ G (with elements G i j ), and therefore may be used to describe nonclassical corrections to two-photon emission via d-parameters.
The expression for the two-photon spectrum is quite general. Now, we specialize it both for the photonic structure and the emitter. For the photon, we consider a planar system (with in-plane translational symmetry). For such a translationally invariant system (and taking a nonretarded, high-Purcell enhancement approximation; both well-realized in our considerations), one has Im G i j (r, r, ω) = ω 6πc F p (r, ω)D i j , where D ≡ diag(1/2, 1/2, 1) and F p (ω) ≡ ρ E ⊥ (r, ω)/ρ E 0 (ω) is the Purcell factor for one-photon emission for a perpendicularly-oriented dipole (frequency ω and position r) near this material (and ρ E 0 (ω) ≡ ω 2 /π 2 c 3 ). For the matter, we consider a hydrogenic atom undergoing a transition between s-states. Under these assumptions, the two-photon differential decay rate is given by dΓ TPE dω = 2 3πc 4 α 2 ω 3 (ω 0 − ω) 3 F p (ω)F p (ω 0 − ω)|T zz | 2 . (73) Comparing this to the free-space two-photon differential decay rate for an s → s transition (see Ref. 41), one finally finds that the enhancement of the two-photon emission spectrum is given by

Framework assumptions
The assumptions underlying the d-parameter framework are (see also Ref. 14): Dipolar nature: The d-parameters emerge from an interface-centered multipole expansion of the quantum mechanical charge and current density. The monopole term gives the classical framework; the dipole term produces d ⊥ and d ; the general nth order multipole term is of order ∼ (k eff x) n , with x a length scale (e.g., d ⊥ and d at n = 1) and k eff an effective modal wavevector (e.g., for a planar interface, the plasmon momentum k; for the dipole resonance of a sphere, the inverse radius R −1 ). Truncation at the dipole term thus produces a leading-order formalism, i.e., we require that {|kd ⊥, |, |R −1 d ⊥, |} 1.
Since |d ⊥, | is Ångström-scale, this is generally very well-satisfied.
Local curvature: The use of d-parameters-whose properties derive from planar interfaces-at curved interfaces implies an assumption of "local flatness". For a sphere, this is equivalent to requiring R |d ⊥, |, i.e., imposes no additional restrictions relative to the dipole-expansion itself.
Surface-centric: Since the d-parameters are surface-response quantities they cannot describe finitevolume effects such as the so-called quantum-size effect [42; 43] (fragmentation of the bulk electronic band structure into discrete levels) nor quantized volume plasmons (the splitting of the bulk plasmon dispersion ω p (k) into a set of discrete "levels"). 11 The assumptions underlying our treatment of the emitter are: Point-emitter: We use macroscopic QED to couple the states of the emitter and the electromagnetic modes of the plasmonic object (Supplementary Note SUPPLEMENTARY NOTE 8). While macroscopic