Classification of scalar and dyadic nonlocal optical response models

: 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 inﬁnitely localized response, as characterized by a distribution with a ﬁnite 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 ﬁnd 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 ﬁelds. We identify the scalar model that comes closest to the hydrodynamic model. Interestingly, for the hydrodynamic Drude model we ﬁnd 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.


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][5][6], for doped semiconductors [3,[7][8][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 in-termediate 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 selfconsistent 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].
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,

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: where the dielectric function has become a nonlocal response function ε ε ε(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].

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 directionk. 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 productk ⊗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 where on the right-hand side we wrote the matrix representation of I in the k , e 1 (k), e 2 (k) basis. For convenience we also introduce the longitudinal and transverse projectors L k ≡k ⊗k and They are indeed projectors and together they project on the entire Fourier space because of their properties In general the dyadic dielectric response in 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 where the scalar functions ε 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.

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: We can do this inverse transformation for the longitudinal and transverse response separately. Let us first transform the longitudinal response ε ε ε 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: So, because of the isotropic response, we can evaluate the three-dimensional integral in two steps: first do the 4π-angular integral over wave vector directionsk 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]): By working out the spatial derivatives, one can show that this angular integral can be written in terms of spherical Bessel functions, see Appendix A and Ref. [34], but here we prefer instead to take the derivatives also out of the k 1 -integral and with ε L (−k Following analogous steps, the transverse response function in real space becomes where ε 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 which is derived and related to spherical harmonics in Appendix A. There is a sum rule for the real-space response functions that we will derive, namely that their integrals over all of space should give dr ε ε ε L,T (r, ω) = d 3 k 1 1 (2π) 3 dr e ik 1 ·r ε ε ε L,T (k 1 , ω) = ε ε ε L,T (k = 0, ω), (11) where in the first equality we used Eq. (6) and in the second one we identified δ (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.

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 plus the boundary condition that G(r, ω) → 0 for r → ∞. The translation invariance means that it is best solved in k-space. Using the vector identityk Now analogous to Eq. (5) for ε ε ε, 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 k and G T (k, ω) equals G L,T (k, ω)T k , and where the 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 k and T k 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 The poles of these Green functions give the dispersion relations for the excitations in the medium. Thus the longitudinal modes correspond to ε L (k, ω) = 0 and the transverse modes to 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 where 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.
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, 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 where the first 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 Here the free-carrier susceptibility χ 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.

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 ε ε ε( which has 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, ω): in terms of the so-called longitudinal delta function δ δ δ (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 in terms of the dyadic transverse delta function δ δ δ ⊥ (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 Notice that in the special case of scalar response, i.e. if ε 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 componentsr so the three-dimensional space is subdivided into a one-dimensional subspace parallel to 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 dr of the dielectric response in Eq. (21), because the 4π solid-angle integral over (I − 3r ⊗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.

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.

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.

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, where β 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 2 T k · E + (ω/c) 2 ε core (ω)E = −iµ 0 ωJ and − β 2 k 2 ω(ω + iγ D ) L k · J + J = −iε 0 ω χ D E, (24) where the variables (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 k , one obtains two coupled equations involving only the transverse fields 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 k gives the longitudinal dispersion relation ε L (k, ω) = 0. This way one obtains these response functions in the HDM as In other words, the transverse response is the familiar Drude-like response well-known from local response models, while the longitudinal response is modified by a k-dependent term in the denominator of the free-electron response. The parameter β 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 selfconsistent hydrodynamic model [26] that can describe electronic spill-out also belongs to the class of longitudinal nonlocal response models.

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 featuring the well-known Drude-like dielectric function of Eq. (17). Thus in the HDM we find terms in the spatial dependence of ε ε ε(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.

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 while the integral over the hydrodynamic Drude part becomes where 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 iqr . We can combine these two results and write the spatially dependent longitudinal response in the HDM as with ε D (ω) and χ D (ω) as defined in Eq. (17).

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 Here the first term on the right-hand side is the usual scalar local Drude-like response of Eq. (17). All effects of spatial dispersion are therefore described by the other term in Eq. (30). Sometimes nonlocal response is said to "smear out" the local delta-function response over a finite volume. We can now determine whether that is indeed what happens for the HDM, the archetypical model with nonlocal response. "Smearing out" would only occur if the second term of Eq. (30) besides smooth functions of 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 iqr /(4π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 ε ε ε(r, ω) = ε core (ω) We have not yet shown that the "one third nonlocal smearing out" story indeed holds, because we have only shown that one third of the local free-carrier response was removed, but we have not checked whether exactly this part was smeared out elsewhere. This balance can be tested by taking the volume integral of Eq. (31). Now the solid-angle integral ∂ Ωr ⊗r/(4π) equals I/3, so that the integral over terms proportional to (I − 3r ⊗r) vanishes. The last term in Eq. (31), the one proportional tor ⊗r, does contribute to the volume integral, and it does so in a neat way since it ensures that indeed dr ε ε ε(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 ε ε ε av. (r, ω) = The first thing to appreciate is that tracing out the HDM response does give a nonlocal response, unlike for the nonlocal model in Sec. 3.2 that had a local trace. The advantages of this scalar model in capturing the HDM response would be that by construction its spatial average is the same as for the HDM, and secondly that here also only one third of the usual local free-electron response is smeared out, and in the third place the 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 kindependent 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 and which was proposed as a physical model of the nonlocal response in semiconductor spheres with very low doping. The parameter ∆ is a resonance frequency and in the limit ∆ → 0, the free-electron susceptibility χ 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 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.

Non-scalar models with
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].

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.

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 trans-