## Abstract

Nonlocal optical response is one of the emerging effects on the nanoscale for particles made of metals or doped semiconductors. Here we classify and compare both scalar and tensorial nonlocal response models. In the latter case the nonlocality can stem from either the longitudinal response, the transverse response, or both. In phenomenological scalar models the nonlocal response is described as a smearing out of the commonly assumed infinitely localized response, as characterized by a distribution with a finite width. Here we calculate explicitly whether and how tensorial models, such as the hydrodynamic Drude model and generalized nonlocal optical response theory, follow this phenomenological description. We find considerable differences, for example that nonlocal response functions, in contrast to simple distributions, assume negative and complex values. Moreover, nonlocal response regularizes some but not all diverging optical near fields. We identify the scalar model that comes closest to the hydrodynamic model. Interestingly, for the hydrodynamic Drude model we find that actually only one third (1/3) of the free-electron response is smeared out nonlocally. In that sense, nonlocal response is stronger for transverse and scalar nonlocal response models, where the smeared-out fractions are 2/3 and 3/3, respectively. The latter two models seem to predict novel plasmonic resonances also below the plasma frequency, in contrast to the hydrodynamic model that predicts standing pressure waves only above the plasma frequency.

© 2015 Optical Society of America

## 1. Introduction

In plasmonics research, conducting materials are described on many levels, from classically all the way to atomistically [1, 2]. Here we focus on effects of nonlocal response (also known as ‘spatial dispersion’) that for metals emerge on the few-nanometer scale, which is the intermediate regime where deviations from classical electrodynamics occur while atomistic descriptions are not needed yet. Nonlocal response means that the polarization field at a certain point depends not just on the electric field in that point, but rather on the electric field in its neighborhood. In Fourier space it is characterized by a wavevector-dependent dielectric response.

In microscopic descriptions of light-matter interactions, nonlocal optical response may ultimately occur due to the quantum mechanical position uncertainty of the particles or material excitations that the light interacts with. The optical properties of conductors is typically dominated by the free-carrier response, which is a collective effect (except for semiconductors at very low doping [3]) that allows semi-classical descriptions of nonlocal response, both for metals [4–6], for doped semiconductors [3, 7–9], and recently also for graphene [10, 11]. Indeed various nonlocal-response models for conductors have been proposed, predicting different phenomena. This article is not an attempt to review all of these models but rather to classify some of them, after introducing a few developments in this active subfield of plasmonics.

Probably the best known semiclassical model to describe nonlocal collective response is the hydrodynamic Drude model (HDM) that goes back to Bloch [4,12] who included the dynamics of the electron gas in the Thomas-Fermi approximation. Although the HDM has been around for a while, it is relevant for current research both because new experiments become possible and because new effects are predicted with the HDM. Three nonclassical phenomena that are well known to be explained by the HDM are (i) size-dependent blueshifts of the dipolar resonance of noble metal nanoparticles as their sizes get smaller than typically 10 nanometers, see for example Ref. [13], (ii) standing bulk-plasmon wave resonances above the plasma frequency, first observed as anomalous absorption in Ag thin films [14], and (iii) somewhat reduced plasmonic field enhancements [15, 16].

A natural extension of the HDM to include nonlocality due to surface diffusion is the Generalized Nonlocal Optical Response (GNOR) model [6, 17]. Recently it was recognized that the HDM predicts that surface-plasmon resonances of higher multipolar order are blueshifted more [18, 19], and in the GNOR model also broadened more [6]. The latter nonlocal broadening was used very recently to explain experimental electron-loss spectra of embedded silver spheres, in particular the disappearance of higher-order plasmon resonances for spheres of radii smaller than 4 nanometers [20].

Very actively studied nowadays are the optical properties of plasmonic dimer systems with (sub-)nanometer-size gaps. Again an important role is played by semi-classical models. Even though they cannot explain all phenomena, to get insight it is important to find out which phenomena can and which cannot be explained by these relatively simple models. In Ref. [16] the HDM was used to show that large dimers with small gaps respond non-classically even though classical electrodynamics suffices to describe the response of the individual monomers. An intermediate nonlocal regime of gap sizes was anticipated to be well described by semiclassical models, with gaps too small for classical electrodynamics to be accurate, yet too large to allow short-circuiting due to electronic spill-out at both interfaces [16]. In some experiments this intermediate nonlocal regime was not identified [21], while in others it was [22]. The observed broadening of optical spectra for smaller-gap dimers [21] is typically interpreted as due to quantum tunneling [23], so it is rather interesting that the semiclassical GNOR theory systematically neglects tunneling but still predicts an analogous spectral broadening of smaller-gap dimers as due to diffusive nonlocal response [6, 17, 24]. What is more, without further ingredients the GNOR theory also explains the observed spectral broadening for smaller monomers, which of course cannot be explained as due to quantum tunneling. Here we have the interesting situation that a semiclassical nonlocal-response theory suggests to reinterpret microscopic quantum calculations, because besides quantum tunneling there must be at least one other broadening mechanism at work in small-gap dimers. This is confirmed in Ref. [25] where it is shown that besides damping due to electron-hole pair creation in forward-scattering processes across the gap (i.e. tunneling) also backward scattering processes at metal-air interfaces are important, the latter processes both for monomers and dimers.

Another very recent generalization of the standard hydrodynamic Drude model is the self-consistent hydrodynamic model (with electron gas dynamics beyond the Thomas-Fermi approximation) that can describe electronic spill-out semiclassically [26]. It predicts nonlocal blueshifts for noble metals while for simple metals redshifts are predicted for smaller nanoparticles [26], both in agreement with experiments. By contrast, the usual (hard-wall) hydrodynamic model systematically excludes spill-out at the interfaces and consequently always predicts nonlocal blueshifts. Furthermore, the self-consistent hydrodynamic model is the first semiclassical theory to describe so-called multipole surface plasmons, also known as Bennett resonances, and which have been experimentally observed, as reviewed in [4].

As a final research trend, it deserves mentioning that effects of nonlocal response are being investigated for metamaterials [27–29], metasurfaces [11] and other periodic structures [30,31].

In this article we will classify several types of nonlocal-response models, both scalar and dyadic, as summarized in Table 1. We will derive their nonlocal response functions in real space, and compare the results with phenomenological scalar models for nonlocal response [32, 33]. Aim is to obtain a sharper intuition about the models considered and the phenomena that they describe, especially the analogies and differences between scalar and dyadic models.

## 2. Preliminaries: nonlocal response in real space and in wavevector space

#### 2.1. From real space to wavevector space: Fourier relations

In general, nonlocal response means that the response of a medium at a certain point **r** depends on the exciting field in a neighborhood of that point, rather than on the exciting field at **r** only. Here for simplicity we only consider dielectric (nonmagnetic) response of a medium to light, and assume the response to be linear. Then the dielectric function (or rather: the 3 × 3 dielectric tensor ** ε**) depends on two positions rather than one, so that the displacement field at a certain point depends on the electric field in a neighborhood of that point:

**(**

*ε***r**

_{1},

**r**

_{2},

*ω*) that depends on two positions. One can also study its Fourier representation

**(**

*ε***k**

_{1},

**k**

_{2},

*ω*) but for arbitrarily shaped media this does not necessarily simplify the analysis. Although nonlocal response in finite geometries will be discussed in this article, we only calculate response functions for translationally invariant systems, where the nonlocal response can only depend on the spatial difference between the two position arguments,

**(**

*ε***r**

_{1},

**r**

_{2},

*ω*) =

**(**

*ε***r**

_{1}−

**r**

_{2},

*ω*) =

**(**

*ε***r**,

*ω*), where we defined the position difference vector

**r**≡

**r**

_{1}−

**r**

_{2}. Due to this translational invariance, the response assumes a simple form in

*k*-space,

**(**

*ε***k**

_{1},

**k**

_{2},

*ω*) = (2

*π*)

^{3}

*δ*

^{3}(

**k**

_{1}−

**k**

_{2})

**(**

*ε***k**

_{1},

*ω*), where the tensor function

**(**

*ε***k**

_{1},

*ω*) is the Fourier transform of

**(**

*ε***r**,

*ω*). For translationally invariant systems, it is usually simpler to first find the dielectric function in

*k*-space rather than in real space. By subsequent inverse Fourier transformation (see below in Sec. 2.3) one then finds the sought real-space form of the nonlocal response that features in the constitutive relation (1). The nonlocal response function of infinite systems so obtained is also useful for the analysis of finite systems, as an approximation for example [32], or as input for a phenomenological specular reflection model for conduction electrons at the metal surface [4].

#### 2.2. Transverse and longitudinal response

The dielectric function is written as a tensor, and indeed in a few important models that we consider below, the response function is tensorial (or dyadic, meaning 3 × 3 tensorial) rather than scalar. For isotropic bulk media there are only two fundamentally different directions in wavevector space, namely the longitudinal direction parallel to a given wavevector, and the transverse directions orthogonal to it. More specifically, given a wavevector **k** of length *k* = |**k**| in Fourier space, **k̂** ≡ **k**/*k* is the corresponding longitudinal unit vector, where ‘longitudinal’ means that it points in the direction of **k**. In addition, we can find two mutually orthogonal *transverse* unit vectors **e**_{1,2}(**k**) that by definition are also orthogonal to **k**. These three unit vectors together thus form an orthonormal basis for *k*-space. An arbitrary vector **v** has a projection of length **v** ·**k̂** in the longitudinal direction **k̂**. The longitudinal projection of the vector is thus given by (**v** ·**k̂**)**k̂**, which is the same as the inner product of **v** with the dyadic product **k̂** ⊗ **k̂**. Sometimes the symbol ⊗ that signifies dyadic product is omitted, but here we keep it. The 3 × 3 unit tensor **I** in *k*-space can then be expanded as

**I**in the {

**k̂**,

**e**

_{1}(

**k**),

**e**

_{2}(

**k**)} basis. For convenience we also introduce the longitudinal and transverse projectors

*k*-space will be different in the longitudinal and in the transverse directions. Even in isotropic media the dielectric response function need not be a scalar. The most general response of an isotropic bulk medium is rather of the form

*ε*

_{L}(

*k*,

*ω*) and

*ε*

_{T}(

*k*,

*ω*) describe the response in the longitudinal and transverse directions, respectively. Because of the assumed isotropy of the medium, the scalar functions

*ε*

_{L,T}(

*k*,

*ω*) depend only on the magnitude but not on the direction of the wavevector.

#### 2.3. Real-space nonlocal response: from k-space back to real space

For a given dielectric tensor in *k*-space, we can find the corresponding dielectric function in real space by inverse Fourier transformation:

*ε*_{L}=

*ε*_{L}(

**k**,

*ω*) back to real space, making use of its property that the scalar function

*ε*

_{L}(

*k*,

*ω*) depends only on the magnitude of the wavevector:

*π*-angular integral over wave vector directions

**k̂**

_{1}, and then the integral over the wave vector magnitude

*k*

_{1}. (Doing the magnitude integral first may also work in some cases.) The angular integral can be evaluated as follows (see the Appendix of Ref. [34]):

*k*

_{1}-integral and with

*ε*

_{L}(−

*k*

_{1},

*ω*) =

*ε*

_{L}(

*k*

_{1},

*ω*) obtain

*ε*

_{T}(

**r**,

*ω*) is the Fourier transform of the scalar function

*ε*

_{T}(

*k*,

*ω*).

We are interested in fully explicit forms of the real-space response functions *ε*_{L,T}(**r**, *ω*). The procedure to arrive at that in subsequent sections will be to assume specific model functions for *ε*_{L,T}(*k*_{1}, *ω*), to perform the *k*_{1}-integrals in Eqs. (9a) and (9b), and then to evaluate the double spatial derivatives. Here we make the general point that the result of doing the *k*_{1}-integral will be a radial function, i.e. depending on the length *r* rather than on the full position difference vector **r**. That makes doing the spatial derivatives easier, since they follow immediately from the dyadic identity

*δ*(

**k**

_{1}) in between the square brackets. The sum rule (11) will be a convenient consistency test of our results. On the right-hand side we see the

*k*-space response in the limit

*k*→ 0. This limit typically coincides with the local-response limit, and then the sum rule describes that the integral over the nonlocal response function equals that of the corresponding local-response model. In this sense, nonlocal-response models ‘smear out’ the response of the corresponding local model. There are infinitely many different functions that have the same

*k*→ 0 limit, and equally many ways to smear out the delta-function response of local models. Below we focus on specific response functions that in most cases are physically motivated.

#### 2.4. Relation between dyadic nonlocal response function and dyadic Green function

The dyadic Green function can conveniently be determined using the longitudinal and transverse projectors of Eq. (3). The retarded (dyadic) Green function for an infinite medium with nonlocal response is defined by the equation

**G**(

**r**,

*ω*) → 0 for

*r*→ ∞. The translation invariance means that it is best solved in

*k*-space. Using the vector identity

**k̂**×

**k̂**×

**A**= (

**k̂**·

**A**)

**k̂**− (

**k̂**·

**k̂**)

**A**= −

**T**·

_{k}**A**, the defining equation becomes

**, the Green function**

*ε***G**(

**k**,

*ω*) for isotropic translationally invariant media can be written as the sum of a transverse and a longitudinal part, where

**G**

_{L}(

**k**,

*ω*) equals

*G*

_{L}(

*k*,

*ω*)

**L**and

_{k}**G**

_{T}(

**k**,

*ω*) equals

*G*

_{L,T}(

*k*,

*ω*)

**T**, and where the

_{k}*G*

_{L,T}(

*k*,

*ω*) are scalar functions depending only on the magnitude of the wave vector. Inserting this expansion into Eq. (13) and using the properties of the projectors

**L**and

_{k}**T**of Eq. (4), we obtain two scalar identities for the transverse and the longitudinal parts of the resulting equation separately, whereby the scalar Green functions can be expressed in terms of the dielectric response functions as

_{k}*ε*

_{L}(

*k*,

*ω*) = 0 and the transverse modes to

*k*

^{2}= (

*ω*/

*c*)

^{2}

*ε*

_{T}(

*k*,

*ω*).

The corresponding real-space representations of the Green functions (14) are found by inverse Fourier transformation, where we can use the isotropy of the medium, analogous to Eqs. (9a) and (9b) for the longitudinal and transverse response functions. The result is therefore

*G*

_{T}(

**r**,

*ω*) is the three-dimensional inverse Fourier transform of the scalar function

*G*

_{T}(

*k*,

*ω*). The total Green function in real space

**G**(

**r**,

*ω*) is given by the sum of

**G**

_{L}(

**r**,

*ω*) and

**G**

_{T}(

**r**,

*ω*). This shows how the total Green function can be computed once the scalar functions

*ε*

_{L,T}(

*k*,

*ω*) are given.

## 3. *k*-independent response [*ε*_{L}(*ω*) and *ε*_{T}(*ω*)]

#### 3.1. Local response [ε_{L}(ω) = ε_{T}(ω)]

When the scalar functions *ε*_{L,T}(*k*, *ω*) actually have no *k*-dependence, it may be expected that the corresponding response in real space will be local. In general this is incorrect, as shown below, but let us first consider the simplest case for which it is true: if the longitudinal and transverse responses are equal and independent of *k*, *ε*_{L}(*k*, *ω*) = *ε*_{T}(*k*, *ω*) = *ε*(*ω*), then ** ε**(

**k**

_{1},

*ω*) =

*ε*(

*ω*)

**I**. In this simple situation we do not need to consider longitudinal and transverse response separately since from Eqs. (2) and (6) it follows immediately that

**I**denotes the 3 × 3 unit matrix in three-dimensional

*k*-space, while the second

**I**is the 3 × 3 unit matrix in real space. Because of the delta function

*δ*

^{3}(

**r**), we see that the response in this simple case is indeed local.

Local response models are employed in most of plasmonics research, and different metals are described by different permittivities *ε*(*ω*). No need to elaborate much on this here. For later reference we mention the widely used Drude-like models, with scalar local response functions of the form

*χ*

_{D}(

*ω*) features the plasma frequency

*ω*

_{p}and Drude damping

*γ*

_{D}, while

*ε*

_{core}(

*ω*) (not further specified here) models the response from the bound ions and electrons, and includes effects of interband transitions.

#### 3.2. Nonlocal response model with local trace [ε_{L}(ω) and ε_{T}(ω)]

A slightly more complex situation occurs if the scalar transverse and longitudinal response functions are different, but still are independent of the wavevector, that is if we assume

*k*-independent matrix elements in the basis of Eq. (2). The real-space response for this model will nevertheless be nonlocal, because the basis of Eq. (2) itself does depend on the wavevector, on its direction to be precise. We now calculate the spatial dispersion explicitly, starting with

*ε*_{L}(

**r**,

*ω*):

*δ*^{‖}(

**r**) ≡ −∇ ⊗ ∇(1/(4

*πr*)), which is a dyadic, see the derivation in the Appendix and in Refs. [35,36]. Hence the longitudinal response does not solely give rise to a delta-function, but also to a “dipolar” nonlocal term scaling as 1/

*r*

^{3}. Similarly, for the transverse part we find

*δ*^{⊥}(

**r**) [35,36]. Again, besides the delta-function term that describes the local part of the response, there is a dipolar nonlocal term scaling with distance as 1/

*r*

^{3}. From the above two expressions it is clear that the adjectives “longitudinal” and “transverse” refer to being parallel or perpendicular to the wavevector

**k**in Fourier space, rather than parallel or perpendicular to the position vector

**r**in real space. The total response of this model is the sum of its longitudinal and transverse parts, and is given by

*ε*

_{T}(

*ω*) =

*ε*

_{L}(

*ω*), the nonlocal 1/

*r*

^{3}terms of the transverse and longitudinal responses cancel exactly. Then Eq. (21) reduces as it should to the scalar local response found earlier in Eq. (16). On the other hand, when

*ε*

_{T}(

*ω*) and

*ε*

_{L}(

*ω*) differ, then the dyadic tensor

**(**

*ε***r**,

*ω*) is a non-scalar dyadic tensor. What does this mean in a bulk system? It means that the response in the direction defined by

**r**differs from the response in directions perpendicular to

**r**. In particular, along

**r**and perpendicular to it we find the diagonal tensor components

**r**and a two-dimensional subspace perpendicular to it. All off-diagonal tensor components are identically zero. We can partly characterize this dyadic response by a scalar, namely its average over all directions, in other words by taking one third of the trace of Eq. (21). The result is [2

*ε*

_{T}(

*ω*) +

*ε*

_{L}(

*ω*)]

*δ*

^{3}(

**r**)/3, which is a local quantity. In other words, the nonlocality of the dyadic response of this model does not show up in its trace. Related to this vanishing angle average, the nonlocal term also disappears when taking the volume average ∫ d

**r**of the dielectric response in Eq. (21), because the 4

*π*solid-angle integral over (

**I**− 3

**r̂**⊗

**r̂**) vanishes. In general one could hope that dyadic nonlocal response functions allow an approximately correct scalar description by taking the angle-averaged response, but the present model shows a qualitative change where taking the average turns a nonlocal model into a local one.

To summarize, in a model in which neither the *k*-longitudinal nor the *k*-transverse dielectric functions depend on the wavevector, we find that in real space there is nevertheless a concomitant spatial dispersion that scales as 1/*r*^{3}, unless the two dielectric functions *ε*_{L}(*ω*) and *ε*_{T}(*ω*) are identical. This result also puts the commonly employed local response models in another perspective: in local-response models, neither the longitudinal nor the transverse responses are local; rather, the overall local response is the result of an exact cancellation between the longitudinal and the transverse nonlocal response functions.

#### 3.3. Examples of nonlocal response models with local trace?

We have not come across situations in which different *ε*_{L}(*ω*) and *ε*_{T}(*ω*) are employed, so in that sense the nonlocal response model with local trace can be considered a toy model that at least gave some insight. It is a rather special model in the sense that the nonlocality does not involve a new length scale. Below we will see instead examples where either *ε*_{L} or *ε*_{T} has an explicit *k*-dependence, and where in the small-*k* limit the two become equal. It is actually an interesting open question whether nonlocal response models with local trace may describe a physical state of affairs in some limit of other dyadic nonlocal response models. Toy models or not, below we will make use of the results obtained for this model in order to understand other nonlocal-response models that do have a clear physical motivation.

## 4. Longitudinal nonlocal response [*ε*_{L}(*k*, *ω*) and *ε*_{T}(*ω*)]

Although we have just seen that even without explicit *k*-dependence in the longitudinal and transverse scalar functions the real-space response is typically still nonlocal, we now look at models where there *is* a *k*-dependence in the scalar longitudinal response response function. Let’s call such models longitudinal nonlocal response models.

#### 4.1. Hydrodynamic Drude model

The hydrodynamic Drude model (HDM) in the Thomas-Fermi approximation in general describes a response of metals to light that is nonlinear [6, 12], but here we will only consider the linear-response limit. The linearized dynamics of the HDM can be described by two coupled equations of motion, in which the electric field is driven by the free-carrier current density and *vice versa*,

*β*is proportional to the Fermi velocity of the conductor. Our results will also apply to the mentioned generalized nonlocal optical response (GNOR) theory, a recent generalization of the HDM where due to diffusion the parameter

*β*becomes complex-valued [17]. Mathematically, by taking the limit |

*β*| → 0 the equation (23b) reduces to Ohm’s law

**J**=

*σ*

**E**. By inserting this into the first equation, we would obtain the Maxwell wave equation in a medium with local response as described by

*ε*

_{D}(

*ω*) of Eq. (17). From the non-classical term proportional to

*β*

^{2}∇ ⊗ ∇·

**J**in Eq. (23b) one can identify the nonlocal response in the HDM as due to non-classical pressure waves in the electron density; the observation of corresponding standing pressure waves in finite structures as in [14] is an important example of the predictive power of the HDM [19].

The above coupled equations indeed also describe finite systems when combined with the right boundary conditions. For infinite systems on the other hand, we can find the longitudinal and transverse response functions in *k*-space, where the above coupled equations turn into

**k**,

*ω*) were dropped for readability. Now both the electric field and the current density have longitudinal and transverse parts. By Eq. (4) we can decompose

**E**(

**k**,

*ω*) uniquely into the sum of

**E**

_{L}=

**L**·

_{k}**E**and

**E**

_{T}=

**T**·

_{k}**E**. This is essentially the Helmholtz decomposition carried out in k-space. Projection properties such as

**L**·

_{k}**E**

_{L}=

**E**

_{L}and

**T**·

_{k}**E**

_{L}=

**0**etcetera follow immediately. An analogous decomposition exists for

**J**(

**k**,

*ω*). By multiplying both equations in (24) from the left with

**T**, one obtains two coupled equations involving only the transverse fields

_{k}**E**

_{T}and

**J**

_{T}, the solution of which gives the dispersion relation

*k*

^{2}=

*ε*

_{T}(

*k*,

*ω*) (

*ω*/

*c*)

^{2}. Multiplying Eq. (24) from the left instead with

**L**gives the longitudinal dispersion relation

_{k}*ε*

_{L}(

*k*,

*ω*) = 0. This way one obtains these response functions in the HDM as

*β*is proportional to the Fermi velocity of the material. Since

*ε*

_{core}(

*ω*) shows up both in the longitudinal and in the transverse response, the total response of the bound electrons is local in this model, whereas the response due to the free carriers is nonlocal, as we shall see. From Eq. (25) one can appreciate that the response in the HDM is both a generalization of the Drude-like local response of Eq. (17) and a special case of the general nonlocal response model (5). We note in passing that the self-consistent hydrodynamic model [26] that can describe electronic spill-out also belongs to the class of longitudinal nonlocal response models.

#### 4.2. Transverse response in real space in hydrodynamic model

From Eq. (25) it can be seen that the transverse dielectric function has no dependence on *k*, exactly as in the result Eq. (20) for the “nonlocal model with local trace” considered in Sec. 3.2. In fact, the integral is identical, and so is the result, namely that in the hydrodynamic Drude model the spatial dependence of the transverse response function is given by

**(**

*ε***r**,

*ω*) that fall off as 1/

*r*

^{3}, originating from the transverse hydrodynamic response. Similar rather singular terms appeared in the total response of “nonlocal model with local trace”, and in the local-response limit they cancelled exactly against similar terms originating from the longitudinal response. It is thus interesting to see whether such cancellation also exists in the HDM, so let us now also consider the real-space version of the longitudinal response in the HDM.

#### 4.3. Longitudinal response in real space in hydrodynamic model

In order to obtain *ε*_{L}(**r**, *ω*) for the HDM, we insert into the integral Eq. (9a) the specific hydrodynamic form of *ε*_{L}(*k*, *ω*) as given in Eq. (25b), which has a bound-carrier part and a free-carrier part. For the first part we need the integral

*q*

^{2}≡

*ω*(

*ω*+

*iγ*)/

*β*

^{2}. It is interesting to see the familiar local Drude susceptibility emerging here, but modified by the spatially dependent factor (1 −

*e*). We can combine these two results and write the spatially dependent longitudinal response in the HDM as

^{iqr}*ε*

_{D}(

*ω*) and

*χ*

_{D}(

*ω*) as defined in Eq. (17).

#### 4.4. Total response in real space in hydrodynamic model

By adding up Eqs. (26) and (29), the same two terms cancel that cancel in local-response models, and we obtain the total real-space response of the HDM in the compact form

*r*contains a delta-function against which some of the local free-carrier response of the first term cancels. Shortly we can and will work out the spatial derivatives of the second term using Eq. (10), but that identity is valid everywhere except in

*r*= 0, and thus is blind for any delta-function contribution. But if for the moment it is only the delta-function contribution that we are interested in, then we can zoom in on ∇ ⊗ ∇[

*e*/(4

^{iqr}*πr*)] so close to the origin that |

*q*|

*r*≪ 1 and find that the sought delta-function contribution indeed exists and equals that of the longitudinal delta function of Eq. (19). Thus we find that the second term in Eq. (30) was hiding a delta-function term −

*χ*

_{D}(

*ω*)

*δ*(

**r**)

**I**/3, which partly cancels the local response described by the first term of Eq. (30). It is interesting that the HDM thereby apparently “smears out” only one third of the local free-electron response over a finite spatial region, while leaving two thirds of the free-carrier local response and (of course) the entire local core response unaffected. We arrive at one of our main results, the explicit spatial response of the HDM as

**I**/3, so that the integral over terms proportional to (

**I**− 3

**r̂**⊗

**r̂**) vanishes. The last term in Eq. (31), the one proportional to

**r̂**⊗

**r̂**, does contribute to the volume integral, and it does so in a neat way since it ensures that indeed ∫ d

**r**(

*ε***r**,

*ω*) equals

*ε*

_{D}(

*ω*)

**I**, the same outcome as for the local-response Drude-like model. This agrees with the sum rule (11) and it also implies that the nonlocal term in Eq. (30) had a vanishing volume integral. But most importantly, it entails that indeed in the HDM one third (no more and no less) of the free-carrier response is smeared out nonlocally. This interesting result is a direct consequence of the dyadic nature of the HDM. In

*k*-space the response is only

*k*-dependent in the longitudinal direction, not in the two other directions, which makes the HDM only one third non-local. This is a rather general argument to explain the factor 1/3 and we did not use any specifics of the hydrodynamic model other than that it belongs to the longitudinal nonlocal class. We therefore anticipate that the argument holds for the entire class, which for example includes models with a fourth-order polynomial of

*k*in the denominator of

*ε*

_{L}(

*k*,

*ω*) instead of the parabolic

*k*-dependence for the HDM in Eq. (25b).

Intuitively one might have expected that nonlocal response gives rise to a smooth real-space response function, but for the HDM the expression Eq. (31) shows several diverging near-field terms. In the near field (i.e. for *qr* ≪ 1), the nonlocal terms are proportional to the longitudinal delta function, thus scaling as *r*^{−3}, both in directions parallel and perpendicular to **r**. By contrast, in the far field (*qr* ≫ 1) the dyadic response along **r** scales with distance as *e ^{iqr}/r*, while perpendicular to

**r**the response falls off faster, namely as

*e*

^{iqr}/r^{2}. So in that sense, the

*k*-dependence of the

**k**-longitudinal response

*ε*

_{L}(

*k*,

*ω*) gives rise to a nonlocal response in real space that is also predominantly

**r**-longitudinal, stronger along

**r**than perpendicular to it.

The response function ** ε**(

**r**,

*ω*) in Eq. (31) is an integration kernel in Maxwell’s equations, in particular in the constitutive relation Eq. (1). Because of the diverging terms in the kernel, evaluating the integral may be numerically much more costly than in case of a Gaussian kernel, say. Instead of solving nonlocal-response problems as an integro-differential equation, it is therefore often simpler to solve coupled differential equations such as in Eq. (23). But one could try to approximate the kernel to make the integro-differential approach more manageable. A first approximation could be to neglect the dipole term out of the kernel (31)), since its volume integral vanishes. A further approximation would be to map the dyadic HDM response onto the scalar response function that captures it best. We propose to use the angle-averaged response, obtained by taking one third of the trace of Eq. (31) and multiplying by

**I**. The result is

*r*-dependence is more similar to the HDM response than a Gaussian function would be. The numerical accuracy of this tailor-made scalar function in capturing hydrodynamic nonlocal response in phenomenological approaches such as in Ref. [32] would be interesting for further study. For more scalar models see Sec. 6.1 and for a warning about using them we refer to the Discussion in Sec. 7.

## 5. Transverse nonlocal response [*ε*_{L}(*ω*) and *ε*_{T}(*k*, *ω*)]

While in the previous section we considered models with *k*-dependent longitudinal and *k*-independent transverse scalar response functions, here we consider just the opposite: models defined by a *ε*_{L}(*ω*) in combination with a function *ε*_{T}(*k*, *ω*) that has an explicit dependence on the wavevector. We discuss some physical predictions of such models in Sec. 7 below, but first let us specify one such a model and study its response in real space. A main question is whether the real-space response function is qualitatively different from the one of the HDM.

We study the model from the recent work of Ref. [8], which in our notation becomes

*χ*

_{free}(

*ω*) reduces to the Drude response

*χ*

_{D}(

*ω*) that we considered above. The corresponding real-space longitudinal response is simple and follows the longitudinal part of the “toy model”, with a spatial dependence of the longitudinal delta function, see Eq. (19). The real-space transverse response requires a bit more work, and after adding up the longitudinal and transverse contributions one finds the total response function

*β*

_{T}by

*β*this

*Q*coincides with

*q*as introduced earlier. As follows from the first line of this expression, in the transverse nonlocal model, two-thirds of the free-electron response is smeared out due to the nonlocality, whereas one third remains a delta-function response. So again we find that the nonlocal spreading is only partial, but this time twice as much as in the hydrodynamic model. This is directly related to the dyadic nature of this transverse nonlocal-response model, and since in

*k*-space the response is nonlocal in two out of three independent directions, namely the two transverse directions perpendicular to

**k**, it starts to become intuitive that only two-thirds of the free-electron response becomes nonlocal. It can be shown that the integral over all of space of Eq. (34) indeed gives the local-response result

**I**[

*ε*

_{core}(

*ω*) +

*χ*

_{free}(

*ω*)], in accordance with the sum rule of Eq. (11).

## 6. Both longitudinal and transverse nonlocal response [*ε*_{L}(*k*, *ω*) and *ε*_{T}(*k*, *ω*)]

#### 6.1. Scalar models [ε_{L}(k, ω) = ε_{T}(k, ω)]

Scalar models for nonlocal response have also been proposed in the literature [7, 37], and most phenomenological models are also scalar [32, 33]. In the light of the previous models, the specialty about scalar models is that their longitudinal and scalar response functions are identical,

*k*-dependence, such scalar models reduce to simple local-response models. Ref. [7] assumes

*ε*(

*k*,

*ω*) to have the functional form of

*ε*

_{T}(

*k*,

*ω*) in Eq. (33b). Here we focus instead on the interesting scalar model of Ref. [37] where

*ε*(

*k*,

*ω*) is assumed to be of the same hydrodynamic form as

*ε*

_{L}(

*k*,

*ω*) in Eq. (25b). The crucial difference with the real HDM is that in the scalar model also an explicit

*k*-dependence of the transverse response is assumed. The model can be derived from two coupled equations analogous to the two equations for the HDM, the only difference being that in the second equation, that is in Eq. (23b), the ∇ ⊗ ∇ is replaced by a Laplacian ∇

^{2}

**I**:

*ε*

_{L}(

*k*,

*ω*) and

*ε*

_{T}(

*k*,

*ω*) are equal so that

**(**

*ε***k**,

*ω*) is scalar and given by the right-hand side of Eq. (25b). This leads to the scalar real-space nonlocal response function

*r*

^{−3}for small

*r*as we found for the HDM, but a 1/

*r*divergence still remains, and on top of that

*ε*(

**r**,

*ω*) is again complex-valued. Unlike in the HDM, here

*all*free-carrier response is smeared out over a finite volume, because nonlocal response was assumed for all three independent directions in

*k*-space. The volume integral of Eq. (37) indeed equals

*ε*

_{D}(

*ω*)

**I**, in accordance with the general sum rule (11). The scalar model (37) that we just obtained can be contrasted with the traced-out HDM model of Eq. (32). Only the latter has the correct amount of smearing out of the local response, whereas Eq. (37) overestimates the amount of nonlocality by a factor of three, at least when assuming that the HDM is accurate.

#### 6.2. Non-scalar models with L- and T-nonlocal response [ε_{L}(k, ω) ≠ ε_{T}(k, ω)]

In the general case that both the transverse and the longitudinal scalar response functions *ε*_{L,T}(*k*, *ω*) are explicitly *k*-dependent but different, the total real-space response functions can be calculated by adding Eqs. (9a) and (9b). In this case all of the local free-space response will be smeared out, one third as determined by the typical length scales in *ε*_{L}(*k*, *ω*), and two thirds on the length scales of *ε*_{T}(*k*, *ω*). When assuming a hydrodynamic form of the response functions, this would correspond to different *β* parameters for the longitudinal and for the transverse scalar responses. Then again the total integral of ** ε**(

**r**,

*ω*) will give the local Drude

*ε*

_{D}(

*ω*)

**I**, because the sum rule (11) holds for longitudinal and transverse parts separately.

The microscopic RPA calculations for a free-electron gas by Lindhard [38] also lead to response functions *ε*_{L,T}(*k*, *ω*) that are both nonlocal (and different). The main difference with the hydrodynamic Drude model is obviously that the latter has a local transverse response. But also the wavevector dependence of the longitudinal response is different, with more pronounced differences for larger *k*. In the long-wavelength approximation on the other hand, Lindhard’s *ε*_{L,T}(*k*, *ω*) both tend to the same classical Drude form, so we do not obtain a nonlocal model with local trace (see Sec. 3.2) by taking the *k* → 0 limit. For applications of the Lindhard response functions also for finite systems, see Ref. [4].

## 7. Discussion

We have compared scalar, longitudinal, and transverse models for nonlocal response and already for translationally invariant systems found qualitatively different predictions, for example whether all local free-carrier response gets nonlocally distributed or not. If we suppose that longitudinal nonlocal models describe metal nanospheres accurately and transverse nonlocal models apply to low-doped semiconductors, can one then still call scalar nonlocal response models phenomenological? It all depends how well they describe physical phenomena.

For finite structures the HDM predicts novel resonances (standing pressure waves of the electron density) due to nonlocal response only above the plasma frequency, and indeed these have been observed experimentally [14]. It was found for metals that scalar nonlocal models with hard-wall boundary conditions exhibit novel nonlocal standing-wave resonances also *below* the plasma frequency, for example in Ref. [37] (and follow-up papers) where a scalar hydrodynamic-like response function was assumed that in real space becomes Eq. (37). Likewise in Ref. [32] with assumed Gaussian scalar nonlocal response functions, near the dipole resonance new standing-wave patterns were predicted on a scale determined by the strength of the nonlocality. Since experimental observations of such resonances in metallic nanoparticles to date have not been reported, the longitudinal hydrodynamic model seems to be the better simple model of nonlocal response in metals. Further discussion and graphical comparisons of the resonances in scalar models as compared to the HDM are given in Ref. [39]. It should be mentioned that non-classical plasmonic resonances of another kind *have* been observed below the plasma frequency for some metals, namely the so-called multipole surface plasmons (reviewed in [4]), but they are related to electronic spill-out at interfaces which classical electrodynamics and the standard HDM both neglect; the first semiclassical model to describe these non-classical resonances is the self-consistent hydrodynamic model [26].

In Eq. (32) we constructed the scalar nonlocal response function that comes closest to the standard HDM. It is an interesting question for future studies how well it captures the predictions of the HDM. However, one concern should be that ‘scalarizing’ the HDM response may again lead to unphysical predictions of novel nonlocal standing-wave resonances in nanostructures. If so, then one should perhaps give up attempts to model nonlocal response phenomena in metals by scalar functions. This would still not rule out transverse nonlocal response in metals, but only the assumption that the transverse response function equals the longitudinal one.

The for metals not-so-physical novel nonlocal resonances below the plasma frequency calculated in Refs. [32, 37] are related to the fact that also the transverse response was assumed to be nonlocal [39]. Let us assume that Ref. [8] is right in modeling the nonlocal response of ultra-low doped semiconductor particles with a transverse nonlocal model. (But here the issue is not settled either, Ref. [8] should be compared with the scalar nonlocal model for doped semiconductors in Ref. [7] and with the longitudinal nonlocal model in Ref. [9]). Then the very interesting question arises whether this transverse model gives rise to new resonances below the plasma frequency, just like in the scalar models of Refs. [32, 37]. If so, then it is to be expected that new nonlocal resonances will show up below the plasma frequency when miniaturizing doped semiconductor nanoparticles. That would constitute a qualitative difference in the nonlocal response of doped semiconductors as compared to metals. Indeed, Ref. [8] features spectral asymmetries and oscillatory behavior on the high-energy side that could very well be due to such novel nonlocal resonances below the plasma frequency.

## 8. Conclusions

We considered different classes of scalar and dyadic nonlocal response models, and in particular their real-space behavior. We were inspired by the intuitive scalar models of Refs. [32, 33] that stressed the smearing out of the local response onto a finite region of space. That is what we also found here, both for scalar and for dyadic models. The simple sum rule Eq. (11) describes exactly this: the spatial integral of the nonlocal response functions equals that of the corresponding local-response model. So whatever delta-function response gets ‘missing’ due to nonlocality is balanced by the sum of the response elsewhere. This equality holds for longitudinal and transverse response functions separately.

What the sum rule does not describe is *how much* of the local response is smeared out in this way. The phenomenological models Refs. [32, 33] assumed implicitly that all free-carrier response is redistributed. That is indeed what we find for scalar models, but for longitudinal nonlocal models including the hydrodynamic Drude model we found to our initial surprise that only one third of the free-electron response is redistributed nonlocally, while for the class of transverse nonlocal models we find two thirds. In this sense scalar models are more nonlocal than transverse ones, which in turn are more nonlocal than longitudinal ones. We explained these properties of real-space response functions back in *k*-space, where longitudinal nonlocal response modified only one out of three directions (i.e. only along **k**), and transverse nonlocal response changed only the two other directions (only normal to **k**). These general arguments apply to the whole classes of models.

Another clear difference with phenomenological models is that the nonlocal response functions ** ε**(

**r**

_{1}−

**r**

_{2}) that we found for some important dyadic models are interesting and rather non-Gaussian: no simple non-negative distributions but instead complex-valued functions that diverge as

*r*

^{−3}for

*r*→ 0. Even for the scalar version of the hydrodynamic model a 1/

*r*divergence remains. Although these divergencies are integrable, as guaranteed by the sum rule Eq. (11), it does explain that solving these nonlocal response problems as integro-differential equations can be slow, and slower for dyadic models such as the HDM than for scalar ones.

We finally suggested that novel resonances below the plasma frequency may arise in transverse nonlocal models. If so, and if the transverse nonlocal model of Ref. [8] indeed applies to doped semiconductors, then effects of novel collective nonlocal resonances may be observed in semiconductor nanoparticles that have no counterpart in metal nanoparticles.

## A. Derivation of Eq. (10) and connection with spherical harmonics

By the chain rule, one obtains the identity for radial functions ∇*f*(*r*) = *f′*(*r*)∇*r*, and by repeating this chain rule once more one obtains ∇ ⊗ [∇*f*(*r*)] = ∇ ⊗ [*f′*(*r*)∇*r*] = *f″*(*r*)∇*r* ⊗ ∇*r* + *f′*(*r*)∇ ⊗ ∇*r*. By writing out in Cartesian coordinates it is easy to show that ∇*r* = **r**/*r* = *r̂*, and after only a little more work that ∇ ⊗ ∇*r* = (**I** − **r̂** ⊗ **r̂**)/*r*. From the previous identities it then follows that

*f*(

*r*) = −1/(4

*πr*) gives

*r*

^{−3}follow from Eq. (38) that was derived for

**r**≠

**0**, thereby missing that a delta-function term must be added, as shown in a few lines in Ref. [35]. To make the connection with spherical harmonics, notice that Eq. (38) can be rewritten as

*f*(

*r*) = sin(

*k*

_{1}

*r*)/(

*k*

_{1}

*r*) this leads to the identity presented in Ref. [34],

*j*that are defined by

_{n}## Acknowledgments

Stimulating discussions in the Structural Electromagnetic Materials group at DTU Fotonik are gratefully acknowledged, in particular with Asger Mortensen, Thomas Christensen, Søren Raza, Nicolas Stenger, Giuseppe Toscano, and Wei Yan. Funding is acknowledged from the Danish Council for Independent Research, Grant No. 1323-00087. The Center for Nanostructured Graphene is sponsored by the Danish National Research Foundation, Project DNRF58.

## References and links

**1. **M. Barbry, P. Koval, F. Marchesin, R. Esteban, A. G. Borisov, J. Aizpurua, and D. Sánchez-Portal, “Atomistic near-field nanoplasmonics: reaching atomic-scale resolution in nanooptics,” Nano Lett. **15**(5), 3410–3419 (2015). [CrossRef] [PubMed]

**2. **X. Chen, J. E. Moore, M. Zekarias, and L. Jensen, “Atomistic electrodynamics simulations of bare and ligand-coated nanoparticles in the quantum size regime,” Nature Commun. **6**, 8921 (2015). [CrossRef]

**3. **H. Zhang, V. Kulkarni, E. Prodan, P. Nordlander, and A. O. Govorov, “Theory of quantum plasmon resonances in doped semiconductor nanocrystals,” J. Phys. Chem. C **118**, 16035–16042 (2014). [CrossRef]

**4. **E. Forstmann and R. R. Gerhardts, *Metal Optics Near the Plasma Frequency* (Springer, 1986). [CrossRef]

**5. **O. Keller, “Screened electromagnetic propagators in nonlocal metal optics,” Phys. Rev. B **34**(6), 3883–3899 (1986). [CrossRef]

**6. **S. Raza, S. I. Bozhevolnyi, M. Wubs, and N. A. Mortensen, “Nonlocal optical response in metallic nanostructures,” J. Phys.: Condens. Matter **27**, 183204 (2015).

**7. **R. Ruppin, Optical properties of a spatially dispersive cylinder,” J. Opt. Soc. Am. B **6**(8), 1559–1563 (1989). [CrossRef]

**8. **R. Carmina Monreal, T. J. Antosiewicz, and S. P. Apell, “Diffuse surface scattering in the plasmonic resonances of ultralow electron density nanospheres,” J. Phys. Chem. Lett. **6**, 1847–1853 (2015). [CrossRef] [PubMed]

**9. **A. Moradi, “Infrared absorption spectra of a spatially dispersive polar semiconductor nanowire,” Solid State Commun. **212**, 10–13 (2015). [CrossRef]

**10. **T. Christensen, W. Wang, A.-P. Jauho, M. Wubs, and N. A. Mortensen, “Classical and quantum plasmonics in graphene nanodisks: role of edge states,” Phys. Rev. B **90**(24), 241414(R) (2014). [CrossRef]

**11. **D. Correas-Serrano, J. S. Gomez-Diaz, M. Tymchenko, and A. Alù, “Nonlocal response of hyperbolic metasurfaces,” Opt. Express **23**(23), 29434 (2015). [CrossRef]

**12. **F. Bloch, “Bremsvermögen von Atomen mit mehreren Elektronen,” Z. Physik A **81**, 363–376 (1933). [CrossRef]

**13. **S. Raza, N. Stenger, S. Kadkhodazadeh, S. V. Fischer, N. Kostesha, A.-P. Jauho, A. Burrows, M. Wubs, and N. A. Mortensen, “Blueshift of the surface plasmon resonance in silver nanoparticles studied with EELS,” Nanophotonics **2**(3), 161–166 (2013). [CrossRef]

**14. **I. Lindau and P. O. Nilsson, “Experimental verification of optically excited longitudinal plasmons,” Phys. Scr. **3**, 87–92 (1971). [CrossRef]

**15. **A. I. Fernández-Domínguez, A. Wiener, F. J. García-Vidal, S. A. Maier, and J. B. Pendry, “Transformation-optics description of nonlocal effects in plasmonic nanostructures,” Phys. Rev. Lett. **108**(10), 106802 (2012). [CrossRef] [PubMed]

**16. **G. Toscano, S. Raza, A.-P. Jauho, N. A. Mortensen, and M. Wubs, “Modified field enhancement and extinction in plasmonic nanowire dimers due to nonlocal response,” Opt. Express **20**(4), 4176–4188 (2012). [CrossRef] [PubMed]

**17. **N. A. Mortensen, S. Raza, M. Wubs, T. Søndergaard, and S. I. Bozhevolnyi, “A generalized non-local optical response theory for plasmonic nanostructures,” Nature Commun. **5**, 3809 (2014). [CrossRef]

**18. **W. Yan, N. A. Mortensen, and M. Wubs, “Green’s function surface-integral method for nonlocal response of plasmonic nanowires in arbitrary dielectric environments,” Phys. Rev. B **88**(15), 155414 (2013). [CrossRef]

**19. **T. Christensen, W. Yan, S. Raza, A.-P. Jauho, N. A. Mortensen, and M. Wubs, “Nonlocal response of metallic nanospheres probed by light, electrons, and atoms,” ACS Nano **8**(2), 1745–1758 (2014). [CrossRef] [PubMed]

**20. **S. Raza, S. Kadkhodazadeh, T. Christensen, M. Di Vece, M. Wubs, N. A. Mortensen, and N. Stenger, “Multipole plasmons and their disappearance in few-nanometer silver nanoparticles,” Nature Commun. **6**, 8788 (2015). [CrossRef]

**21. **K. J. Savage, M. M. Hawkeye, R. Esteban, A. G. Borisov, J. Aizpurua, and J. J. Baumberg, “Revealing the quantum regime in tunnelling plasmonics,” Nature **491**, 574–577 (2012). [CrossRef] [PubMed]

**22. **H. Jung, H. Cha, D. Lee, and S. Yoon, “Bridging the nanogap with light: continuous tuning of plasmon coupling between gold nanoparticles,” ACS Nano, Article ASAP (2015). [CrossRef]

**23. **R. Esteban, A. G. Borisov, P. Nordlander, and J. Aizpurua, “Bridging quantum and classical plasmonics with a quantum-corrected model,” Nature Commun. **3**, 825 (2012). [CrossRef]

**24. **S. Raza, M. Wubs, S. I. Bozhevolnyi, and N. A. Mortensen, “Nonlocal study of ultimate plasmon hybridization”, Opt. Lett. **40**(5), 839–842 (2015). [CrossRef] [PubMed]

**25. **W. Yan, M. Wubs, and N. A. Mortensen, “Projected-dipole model for quantum plasmonics,” Phys. Rev. Lett. **115**(13), 137403 (2015). [CrossRef] [PubMed]

**26. **G. Toscano, J. Straubel, A. Kwiatkowski, C. Rockstuhl, F. Evers, H. Xu, N. A. Mortensen, and M. Wubs, “Resonance shifts and spill-out effects in self-consistent hydrodynamic nanoplasmonics,” Nature Commun. **6**, 7132 (2015). [CrossRef]

**27. **W. Yan, M. Wubs, and N. A. Mortensen, “Hyperbolic metamaterials: nonlocal response regularizes broadband supersingularity,” Phys. Rev. B **86**(20), 205429 (2012). [CrossRef]

**28. **W. Yan, N. A. Mortensen, and M. Wubs, “Hyperbolic metamaterial lens with hydrodynamic nonlocal response,” Opt. Express **21**(12), 15026 (2013). [CrossRef] [PubMed]

**29. **A. Moradi, Maxwell-Garnett effective medium theory: quantum nonlocal effects,” Phys. of Plasmas **22**(4), 042105 (2015). [CrossRef]

**30. **A. Yanai, N. A. Mortensen, and U. Levy, “Absorption and eigenmode calculation for one-dimensional periodic metallic structures using the hydrodynamic approximation,” Phys. Rev. B **88**(20), 205120 (2013). [CrossRef]

**31. **F. Intravaia and K. Busch, “Fluorescence in nonlocal dissipative periodic structures,” Phys. Rev. A **91**(5), 053836 (2015). [CrossRef]

**32. **P. Ginzburg and A. Zayats, “Localized surface plasmon resonances in spatially dispersive nano-objects: phenomenological treatise,” ACS Nano **7**(5), 43344342 (2013). [CrossRef] [PubMed]

**33. **N. A. Mortensen, “Nonlocal formalism for nanoplasmonics: phenomenological and semi-classical considerations,” Phot. Nanostr. **11**(4), 303 (2013). [CrossRef]

**34. **P. de Vries, D. V. van Coevorden, and A. Lagendijk, “Point scatterers for classical waves,” Rev. Mod. Phys. **70**(2), 447 (1998). [CrossRef]

**35. **D. P. Craig and T. Thirunamachandran, *Molecular Quantum Electrodynamics. An Introduction to Radiation Molecule Interactions* (Academic, 1984); see Sec. 3.1: Transverse and longitudinal δ-dyadics.

**36. **C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, *Photons and Atoms. Introduction to Quantum Electrodynamics* (Wiley, 1989); see pp. 36–42.

**37. **J. M. McMahon, S. K. Gray, and G. C. Schatz, “Nonlocal optical response of metal nanostructures with arbitrary shape,” Phys. Rev. Lett. **103**(9), 097403 (2009). [CrossRef] [PubMed]

**38. **J. Lindhard, “On the properties of a gas of charged particles,” Kgl. Dansk Videnskab. Selskab, Mat. Fys. Medd. **28**(8), 1–57 (1954).

**39. **S. Raza, G. Toscano, A.-P. Jauho, M. Wubs, and N. A. Mortensen, “Unusual resonances in nanoplasmonic structures due to nonlocal response,” Phys. Rev. B **84**(12), 121412(R) (2011). [CrossRef]