A universal macroscopic theory of surface plasma waves and their losses

Recently, we have revealed an intrinsic instability of metals due to surface plasma waves (SPWs) and raised the prospect of using it to create lossless SPWs. The counter-intuitive nature of this finding prompts one to ask, why had not this instability been disclosed before, given the long history of this subject? If this instability does exist, how far is it from reality? The present work is devoted to answering these questions. To this end, we derive a unified macroscopic theory of SPWs that applies to any type of electron dynamics, be they local or non-local, classical or quantum-mechanical. In light of this theory, we analyze the behaviors of SPWs according to several electron dynamics models, including the widely used local dielectric model (DM), the hydrodynamic model (HDM) and the specular reflection model (SRM), in addition to the less common semi-classical model (SCM). We find that, in order to unveil the instability, one must (i) self-consistently treat surface effects without any of the usually imposed auxiliary conditions and (ii) include translation symmetry breaking effects in electron dynamics. As far as we are concerned, none existing work had fulfilled both (i) and (ii). To assess the possibility of realizing the instability, we analyze two very important factors: the dielectric interfacing the metal and inter-band transitions, which both were ignored in our recent work. Whereas inter-band absorption -- together with Landau damping -- is shown adverse to the instability, a dielectric brings it closer to occurrence. One may even attain it in common plasmonic materials such as silver under not so tough conditions.


I. INTRODUCTION
Electron density ripples propagating along metal surfaces, known as surface plasma waves (SPWs), have been intensively pursued as a promising enabler of nano photonics [1][2][3][4] . A fundamental issue hampering further progress is concerned with energy losses 5 . SPWs lose energy fast due to Joule heating, Landau damping and inter-band absorption alongside others such as radiation losses 5 . It has been suggested that these losses -especially those due to Landau damping -are intrinsic and cannot be significantly reduced without the addition of external gain medium [5][6][7][8][9] . They would ultimately handicap the functionalities of SPWs 5,6 .
Our recent work 10-12 challenged the above view and found that, thanks to an incipient instability, the losses may well be reduced to any level without taking energy from outside. A critical point was shown to exist at γ 0 τ = 1, where τ −1 is the thermal electronic collision rate and γ 0 is a positive-definite quantity. At the critical point, the energy released from the Fermi sea -at rate γ 0 -just compensates for the dissipation due to electronic collisions and SPWs become lossless 12 . This instability was revealed through a semi-classical model (SCM, based on Boltzmann's transport equation) of SPWs. Given the highly counter-intuitive nature of this possibility and the long history on this subject dated back to late 1950s 13 , the following questions must then be answered: Why has the criticality not been discovered in earlier work and what is missing in those work? While it is not widely used in the study of SPWs, the SCM had certainly been considered before 14,15 . Yet nobody had claimed lossless SPWs.
The main purpose of the present paper is to answer the above question. Our strategy is to construct a universal theory of SPWs that is applicable to any type of electron dynamics, be they local or non-local, classical or quantum-mechanical. In light of this theory, the properties of SPWs are then analyzed with several common electron dynamics models, including the local dielectric model (DM) 13 , the classical hydrodynamic model (HDM) [16][17][18][19] and the specular reflection model (SRM) [20][21][22] in addition to the SCM. The chief result obtained through this analysis is that, in order to unveil the instability, (i) a SPW theory must self-consistently deal with physical surfaces, and (ii) translation symmetry breaking effects must be included in the electron dynamics. Amongst these models, only the SCM is capable of (ii), but previous work based on the SCM had all failed with respect to (i). This explains why the existing work had failed to hit upon the possibility of lossless SPWs.
In existing work based on non-local models, auxiliary conditions have usually been imposed [23][24][25][26] , mostly assuming no normal current on surfaces, as in the hard-wall picture adopted in computational studies [27][28][29][30] . From a conceptual point of view, however, these conditions are obviously incompatible with local models and arise only due to an incomplete description of physical surfaces 31 . We will show that the conventional treatment of SPWs via auxiliary conditions does not recover the standard textbook SPWs in the local models. We prescribe a simple yet complete macroscopic description of physical surfaces, which remedies the conceptual deficiency and allows to derive a universal SPW theory accounting for surface effects self-consistently.
As a secondary purpose, we wish to address an experimentally interesting issue, that is, how far is the instability from reality and how can it be achieved? In general, τ −1 is comparable to the characteristic plasma frequency ω p , even in defectfree materials and at zero temperature 32 . This is because τ −1 is the collision rate at the SPW frequency, which is in the order of a few eVs in metals and thus effects as an effective high temperature of tens of thousands of Kelvins opening up a large phase volume for electron scattering. One then expects γ 0 τ < 1 usually. In order to devise a practically useful method of enhancing γ 0 , we employ our theory to analyze two important factors affecting γ 0 : a dielectric interfacing the metal supporting SPWs and inter-band transition effects, which were ignored in Refs. [10][11][12] . We find that, inter-band absorption (as well as Landau damping) strongly reduces γ 0 , whereas topping a dielectric can significantly increase γ 0  instability may well be in reach in a single crystal of silver. This paper is organized as follows. We develop a universal macroscopic SPW theory in the next section and apply it to the DM, the HDM and the SRM in Sec IV. Section III includes a discussion of a phenomenological approach to interband transition effects to make the paper self-contained. Section V is devoted to a thorough treatment of SPWs within the SCM. Dielectric and inter-band effects are analyzed. We assess the possibility of creating lossless SPWs by the SCM and conclude the paper in Sec. VII. In Appendices A and B, we discuss some historical misconceptions regarding the SPWs in the HDM and the SRM, respectively.

II. UNIVERSAL THEORY OF SPWS
In this section, we present a general theory of charge density waves in semi-infinite metals (SIMs). Extension to systems of other geometries is straightforward and will be considered elsewhere. The theory is formulated in terms of universal physical concepts and makes no resort to the particulars of electron dynamics. It is a macroscopic theory, thus valid as long as the thickness (d s , typically a few lattice constants) of the microscopic surface layer, which forms between the vacuum and the bulk metal, is much smaller than the SPW wavelength Λ. In this paper, by 'macroscopic' we always mean d s /Λ ≪ 1, with no regard to electron dynamics.
Retardation effects are neglected throughout this paper, which is reasonable provided the SPW phase velocity is much less than the speed of light in vacuum.

A. Macroscopic description of surfaces
The SIM under consideration is assumed to occupy the region z ≥ 0 and bounded by a flat interface/surface macroscopically located at z = 0. The other half space is either the vacuum or a dielectric with dielectric constant ǫ d , see Fig. 1. We shall use the vector x = (r, z) to denote a point in space, where r = (x, y) is the planar projection of x. Further, x 0 = (r, 0) denotes a point on the surface and t denotes time.
In a macroscopic description of SPWs, one usually considers the metal as a medium for an electromagnetic field and seeks surface localized (polariton) solutions of the governing Maxwell's equations 33 . The procedure is to write down the waves for (an infinite metal) on the metal side and those (for an infinite dielectric) on the dielectric side, and then invoke conditions to match those waves at the surface. With local dynamics, the usual Maxwell's boundary conditions would do the job. With non-local dynamics, however, they are insufficient. Historically, auxiliary conditions -mostly assuming zero normal current at the surface 14,15,17,25 -have been invented as a remedy since 1950s 23 , which however treat the symptoms not the cause. The cause is a conceptual deficiency in the knowledge of surfaces 31 . Considering that a real microscopic surface can hardly be specified even for the simplest material due to preparation procedures, one might deem it hopeless to have a complete description. The following elementary analysis suggests otherwise.
Let us imagine bringing two materials (A and B) in contact, and an interfacial layer of thickness d s -in the order of a few lattice constants -shall form in between. We may characterize this layer by a surface potential φ s , which should quickly decay to zero in the bulk regions outside the interfacial layer. The exact microscopic profile of the layer varies from one case to another and can hardly be known a priori. Despite this, we may still write down a generic form for the electric current density j(x, t) in the whole system including the interfacial layer. To this end, we observe that in the bulk regions where φ s vanishes, the form of j(x, t) can be completely specified with the respective dynamic equations for the infinite materials, apart from some parameters (such as the Fuchs parameter, see Sec. V) that encode the effects of surface scattering on the electron waves. Let us denote by J A/B (x, t) the forms of j(x, t) in the bulk region of A/B. Microscopically, j evolves from J A in the bulk region of A, through a rapid variation in the interfacial layer, to J B in the bulk region of B. Formally, we can write for the µ-th component of the current where the profile functions w α (z) approach unity in the bulk region of A and zero in that of B. The exact profile of w µ (z) depends on the microscopic details of the interfacial layer. On the scale of λ, however, the interfacial layer appears infinitely thin and w α (z) reduce to the Heaviside step function Θ(z), where Θ(z ≥ 0) = 1 and Θ(z < 0) = 0, in a macroscopic theory. One thus ends up with This holds valid for any w µ (z) and is thus a general and complete macroscopic description of a physical interface, as long as a perturbation on one side does not cause significant responses on the other. To recapitulate, equation (1) elegantly captures two importance physical consequences of an interface: the rapid variation of the current density through the step function Θ(z) and the surface scattering effects on electron dynamics through the parameters contained in the bulk forms. These scattering effects -including the symmetry breaking effects -have been ignored in most models except for the SCM. In general J A and J B are not equal on the interface, as is certainly the case for local dynamics models, and charges accumulate in the interfacial layer. Such capacitive effects would be mistakenly erased under auxiliary conditions, which often dictate continuity of current across the interface, e.g. the vanishing of normal current at the metal-vacuum interface.

B. Charge density waves
Now we derive the equations of motion for the charge density ρ(x, t) in the SIM. Our starting point is the equation of continuity. Specifying Eq. (1) to the SIM, we have j(x, t) = Θ(z)J(x, t), where J(x, t) is the current density in the metal. The equation of continuity then reads Here we have included a global relaxation term −ρ(x, t)/τ to account for the thermal relaxation of non-equilibrium charges due to microscopic electronic collisions driving the system toward thermodynamic equilibrium. In terms of J(x, t), this equation becomes where Θ ′ (z) = ∂ z Θ(z). The right-hand term of this equation signifies the capacitive effects, which is critical in the energy conversion process but had been overlooked until our recent work 12 . This term was also noticed by A. L. Fetter in studying the edge plasmon in two-dimension systems 19 . Without loss of generality, we may assume a quasi-plane wave for the fields and write ρ(x, t) = ρ(z)e i(k·r−ωt) and similarly for J(x, t) and the electric field E(x, t), where the wave vector k = (k ≥ 0, 0) is directed along the positive x-axis for the sake of definiteness. In general the frequency ω = ω s − iγ is complex. The relation between ω and k is determined by the wave equations to be established in what follows. Equation (2) thus becomes whereω = ω + i/τ and J z (x 0 , t) = J z (0)e i(k·r−ωt) . We shall writeω = ω s + iγ 0 , so that by definition. We shall see that γ 0 is negative in all models except for the SCM, in which it is positive-definite thanks to a fundamental physical reason.
For linear responses, J(z) can be related to E(z) as follows where σ µν (z, z ′ , ω) is the conductivity tensor. Considering that E(z) linearly depends on ρ(z), we can define a linear operator H with kernel H(z, z ′ ) so that It is easy to show that Note that S does not depend on z. Since ρ(z) is defined only for z ≥ 0, we can introduce a cosine Fourier transform In terms of ρ q , Eq. (7) can be rewritten as Here is the matrix element between the cosine waves. Finally, we close Eq. (8) by the fact that J z (0) and hence S are also linear functionals of ρ q . We can thus write where K = (k, 0, q), K 2 = k 2 + q 2 and G is a kernel given by As a key result of this paper, Eqs. (8) and (9) constitute a complete description of self-sustained charge density waves in SIMs. Their basic structures are valid regardless of the underlying electron dynamics encoded in H and G.

C. Universal secular equation for SPWs
Equation (8) is an inhomogeneous linear equation with a source term. Equations of this type generally admit of two classes of solutions, depending on whether S vanishes or not. Solutions with S = 0 represent nothing but the bulk plasma waves, while those with S 0 represent localized waves, i.e. SPWs in our system. The spectra of these two classes generally do not overlap. Had we imposed auxiliary boundary conditions as is usually done in the literature, SPWs would not exist at all.
To find the 'secular' equation for SPWs, we wish to make a simplification, which is not necessary but will make the resulting expressions more transparent. To this end, let us look at some general properties of σ µν (z, z ′ , ω). For an infinite system without boundaries, the translational symmetry is preserved along z-axis as well as along the surface plane, in which case σ µν only depends on z − z ′ . For SIMs, however, that symmetry is broken and σ µν must in general depend on z and z ′ separately. It shall prove useful to decompose σ µν into two parts, σ b,µν and σ s,µν , where σ b,µν is defined to be that of the infinite system and depends only on z − z ′ while σ s,µν signifies pure surface effects. Namely, Accordingly, H and G also each contain two parts arising from σ b,µν and σ s,µν , respectively. Let us then write H = H b + H s . Now that it reflects on the properties of plasma waves of an infinite system, where Ω(K, ω) is a frequency that only depends on K. By definition, one can easily show that the dielectric function of an infinite system is given by As for H s , it gives rise to scattering of plasma waves. Nevertheless, for bulk waves the scattering due to a surface should be insignificant and may be treated perturbatively, as we did in Refs. 10,11 . To the lowest order in this perturbation, we may simply put This expression contains complete information of bulk plasma waves. As expected, the zeros of ǫ(K, ω) give the dispersion ǫ b (K) of these waves. We analogously write G = G b + G s . Note that via Eq. (9) G b determines J z at the point z = 0 in a infinite system, for which, however this point is none more special than any other points. Thus, G b can only be a constant that does not depend on q and it is therefore purely due to the local part of σ b,µν (z − z ′ ). This part must be isotropic for the jellium model and it can be written as δ µν δ(z − z ′ )σ(ω), where δ µν is the Kroneckle symbol and σ(ω) is to be discussed further. It follows that G b = −4iωσ(ω)kβ, where we have taken into account the effect of the dielectric via the factor β = 2ǫ d /(1 + ǫ d ); see details below. We now obtain We see that G s describes translation symmetry breaking effects. As to be seen later, it is missing from the DM, the HDM and the SRM. The positive-definiteness of γ 0 arises solely through this term and therefore cannot be captured in these models. For SPWs, S 0. Inserting Eq. (11) in (8) and using (9), we find which determines the SPW frequency ω s and damping rate γ as a function of k. This is another key result of this paper. The SPW charge density is obtained as .
The resulting ρ(z) is peaked about the surface. Equations (13) and (14) constitute a universal description of SPWs. They are valid irrespective of the underlying electron dynamics, which enter only through ǫ(K, ω) and G(K, ω).

D. Dielectric effects
Let us place a dielectric with dielectric constant ǫ d -which can be complex -on the side z < 0. The electrostatic fields are affected by this dielectric, which can be calculated with the method of mirror charges. Let the mirror charge density in the dielectric be ρ d ( The electrostatic potential φ( In the metal, where z ≥ 0, the electric field and where β has been given in the above. See that E z (0) and hence G b are enhanced by the factor β, which is why this factor appears in Eq. (12). Physically this is because the surface sits between the charges in the metal and those in the dielectric and spatially separates them.

III. INTER-BAND TRANSITION EFFECTS
The dynamics models to be considered in the next two sections only describe the currents from electrons in the conduction band. Realistically, valence electrons can also contribute by virtual inter-band transitions. In this paper, for the sake of simplicity, we account for these inter-band transition effects by a phenomenological approach [35][36][37] . The observation is that, valence electrons are usually tightly held to their host atoms and the energy bands are largely non-dispersion, and thus their electrical responses are mostly local and not susceptible to the presence of boundaries. We may then describe this response by a local conductivity function σ p (ω), so that the electrical current density due to the valence electrons is given by Usually σ p (ω) may be modeled in the Lorentz form. It is related to the inter-band dielectric function by ǫ p (ω) = 4πiσ p (ω)/ω, which can be measured for example by means of ellipsometry or computed by density functional theory. ǫ p (ω) contains a real part ǫ pr (ω) and an imaginary part ǫ pi (ω). While ǫ pr (ω) acts to shield the conduction electrons, ǫ pi (ω) -which is always positive -leads to inter-band absorption. Being basically an atomic property, ǫ p is not sensitive to temperature.

IV. SPWS IN DM, HDM AND SRM
While the theory established in Sec. II is universal, the behaviors of SPWs do depend on electron dynamics through ǫ(K, ω) and G(K, ω). In this section, we apply the theory to examine SPWs within several common electron dynamics models: the DM, the HDM and the SRM. SPWs in the SCM will be thoroughly treated in the next section. We show that the usually quoted SPW solution in the HDM is false and clarify the origin of SPWs in the SRM. A summary of this analysis is tabulated in Table I A. Local dielectric model (DM) This is the standard model for SPWs. Unlike other models, it does not require and is incompatible with any auxiliary conditions. Here we reproduce by our theory the well-known properties of SPWs in this model.
According to the DM, the current density due to conduction electrons is given by where ω p = 4πn 0 e 2 /m is the characteristic plasma frequency of the metal with m and e being the effective mass and charge of an electron, respectively. The total current density J(z) then is yielding σ(ω) = σ DM (ω) + σ p (ω). From Eq. (6) one finds Ω 2 dispersionless, given as Bulk plasma waves. Equating Ω 2 0 withω 2 yields the frequency ω b,DM and damping rate γ b,DM for bulk plasma waves. For ω p τ → ∞ and assuming ǫ p (ω) independent of ω, they are given by Note that bulk waves bear no dielectric effects, i.e. no dependence on ǫ d . The damping rate γ b,DM arises due to thermal collisions and inter-band absorption.
SPWs. The electrical conductivity is purely local and thus G s = 0. As such, G(K, ω) becomes . Substituting this into Eq. (13), we immediately arrive at the often quoted frequency ω s,DM and damping rate γ s,DM for SPWs. Neglecting absorption in the dielectric, i.e. assuming real ǫ d , we find In this model, the SPW charge density is completely localized on the surface, gives the areal charge density. See that the dielectric tends to reduce the SPW damping rate.
Traditionally 38 , the above results have been obtained by treating the metal as a simple dielectric with dielectric constant ǫ(ω) = 1 − Ω 2 0 /ω 2 . Then exponentially decaying electromagnetic (EM) waves (or electrostatic potentials in the quasistatic limit) are written down on the metal and the dielectric sides, and Maxwell's boundary conditions are used to match the waves to obtain the above frequency and damping rate of SPWs. Our theory works directly with charge density rather than EM waves. The two approaches are equivalent.

B. Hydrodynamic model (HDM)
The DM assumes a purely local relation between the current density and the electric field. The HDM extends the DM by inclusion of leading-order non-local corrections. Recently, this model has attracted lots of attention in plasmonics and quantum forces 39,40 . It has also been synergized with density functional theory in the quantum hydrodynamic model [41][42][43][44] to study local plasmon resonances on metal particles.
In the HDM, Ω 2 includes leading-order dependence on K and is given by where v 0 is a parameter. The dielectric function 45 then reads ǫ HDM = 1 − Ω 2 HDM /ω 2 . The bulk waves are similar to those in the DM except for some dispersion with K.
where ω p denotes the characteristic plasma frequency of the metal. F = (m/2π ) 3 where v is the velocity of electrons of mass m. ǫ(K, ω) is the dielectric function of the infinite system, whose zeros give the dispersion of bulk plasma waves. SPWs are determined by this secular equation: The solution to this equation is written as ω = ω s − iγ. In the SCM, an extra contribution γ s arises due to symmetry breaking effects contained in G s . Very interestingly, γ 0 = γ s − γ Landau − γ interband cannot be negative.

Quantity
DM HDM SRM SCM Symmetry breaking effects As in the DM, no symmetry breaking effects are included in the HDM. G(K, ω) is then the same as that for the DM. Now Eq. (13) transforms into the following which determines the SPWs in the HDM by our theory. Note that the integrand contains no poles or resonances near the solutions, as the SPW spectra are always gapped from the bulk wave spectra gratifying ǫ HDM = 0. For v 0 = 0, the HDM reduces to the DM and so do the SPWs, as expected. With ǫ d = 0 and ǫ p = 0, the equation simplifies to where we have used the fact that ǫ HDM is even in q and that ω 2 p /ω 2 ≈ 2 for the solutions to this equation. As shown in Appendix A, it leads to a linear dispersion of ω s versus k. As in the DM, the SPW damps due to thermal collisions and inter-band absorption, at rate γ s,HDM ≈ γ s,DM apart from some dispersion effects.
In the literature, the condition that J z (0) = 0 is usually imposed in the HDM 15,17 . This would mean S = 0 and therefore would exclude any SPWs according to our theory. Nevertheless, SPWs have been claimed to exist under this condition.
In what follows we briefly show that this claim is false, more details to be found in Appendix A.
For illustration, we take ǫ p = 0 and ǫ d = 0. Impose S = 0 and the wave equation becomes Ω 2 HDM −ω 2 ρ q = 0, or equivalently in the real space The claimed SPW solution is then sought of the form ρ(z) = ρ 0 e −κz . Substituting this in the equation leads toω 2 = ω 2 p + v 2 0 (k 2 − κ 2 ). Imposing J z (0) = 0 gives another relation, , which expresses the balance between the electronic pressure and the electric force at the surface. Here Those two relations specify the solution. Combined, they lead toω 2 ≈ ω 2 0 + ω 0 βk. Nevertheless, this solution does not reduce in the limit v 0 = 0 to the SPWs found by Ritchie with the DM. Actually, κ diverges in this limit, yielding ρ(z)dz ∼ κ −1 = 0, i.e. the solution is empty of charges. This false solution is also what was observed in Refs. 14,15,47 . It is plausible that existing ab initio calculations have only observed this false solution as well 28 . A comprehensive account may merit a future study.
Although the false solution and the correct solution are conceptually disparate, their dispersion relations are quite similar, as shown in Appendix A.

C. Specular reflection model (SRM)
A natural step to go beyond the HDM is to use the full form of Ω rather than the approximation Ω HDM . Equation (22) then becomes This is exactly the equation established by Marusak and Ritchie in 1966 for the SRM 20 . Our derivation makes it clear that the SRM is just an extension of the HDM. From this point of view, one may also conclude that the usually claimed SPWs in the HDM are false, because they are not solutions of Eq. 24 in the HDM limit.
In contrast with the DM and the HDM, SPWs in the SRM can also decay via Landau damping, because of an imaginary part in Ω associated with electron-hole excitations. Thus, the SPW damping rate is γ s,SRM ≈ γ s,DM + γ Landau , see the next section for further discussion on this.
We wish to point out a logical inconsistency in the original contrivance of the SRM. There are two elements in this contrivance: (i) as nominally expected, electron waves impinging on the surface are assumed to be specularly reflected back, and (ii) a sheet of 'fictitious' charges exactly localized on the surface. Element (i) would mean J z (0) = 0 and hence, by our theory, no SPWs would materialize. Then how do those waves come about? The answer rests with element (ii). In Appendix B, we show that the fictitious charge sheet reinstates the capacitive effects lost under element (i). Actually, SPWs appear as a pole of this fictitious charge density.
As with the HDM and the DM, the SRM contains no symmetry breaking effects, i.e. G s is absent from these models. To account for these effects, further improvement is required, leading to the SCM.

V. SPWS IN THE SCM
In the SCM one calculates the electric currents in terms of a distribution function f (x, v, t) defined in the single-particle phase space. As usual, we write it as a sum of the equilibrium part f 0 (ε(v)) and the non-equilibrium part g(x, v, t). f 0 (ε) is taken to be the Fermi-Dirac function at zero temperature. ε(v) = mv 2 /2 is the dispersion of the conduction band. Within the relaxation time approximation and the regime of linear responses, g(x, v, t) = g(v, z)e i(k·r−ωt) satisfies the following Boltzmann's equation where m is the electron effective mass and v F the Fermi velocity.
With γ 0 ≥ 0, the general solution to Eq. (25) can be written as where C(v) = g(v, 0) is the non-equilibrium deviation on the surface to be determined by boundary conditions. We require g(v, z) = 0 distant from the surface, i.e. z → ∞. For electrons moving away from the surface, v z > 0, this condition is automatically fulfilled. For electrons moving toward the surface, v z < 0, it leads to It follows that To determine C(v) for v z > 0, the boundary condition at z = 0 has to be used, which, whoever, depends on surface properties. We adopt a simple picture first conceived by Fuchs 48 and then widely used in the study of for instance anomalous skin effect [49][50][51][52] . According to this picture a fraction p i.e. the Fuchs parameter varying between zero and unity, of the electrons impinging on the surface are specularly reflected back, i.e.
This condition is identical with the condition used in Ref. 12 but differs from that in Refs. 10,11 except for p = 0. It follows that Equations (26) -(30) fully specify the distribution function for the electrons. The electrical current density due to the conduction electrons is then calculated as Note that the charge density is not given bỹ The reason is simple: the as-obtained g(v, z) is for the bulk region and not valid on the surface, because Eq. (3.1) involves no surface potentials, see Sec. II A. Actually, J c (x, t) andρ(x, t) obey the equation rather than the equation of continuity [c.f. Eq. (2)], automatically embodying the condition that J z (0) = 0. This underlies the incorrect conclusion drawn by Harris and others 14 .
A. The positive-definiteness of γ 0 We substitute the expression of E(z) given by Eqs. (17) and (18) into (26) -(30) and perform the integration over z ′ . We find it instructive to split g(v, z) into two parts, one denoted by g b (v, z) and the other by g s (v, z). They are given by where we have introduced the following functions, Note that F +/− is an even/odd function of v z . In addition, (34) Moreover, we have One may show that g b (v, z) can also be obtained directly by Boltzmann's equation for an infinite system. Thus, this part contains exactly the responses of an infinite system. It is independent of surface properties, i.e. showing no dependence on the Fuchs parameter p, and the electrons incident on the surface (i.e. with v z < 0) and those departing it (i.e. with v z > 0) appear on equal footing in its expression. On the other hand, g s (v, z) signifies pure surface effects: it exists only for departing electrons, as indicated by the Heaviside function Θ(v z ) in its expression, and it depends on p thus reflecting on surface roughness. If we keep only g b (v, z), the SRM will be revisited, making it evident that the SRM does not correspond to the limit of p = 1 (specularly reflecting surface), in contrast with its nominal meaning.
Another important feature of g s (v, z) lies in its dependence on z, i.e. g s (v, z) ∝ e iωz/v z ∝ e −γ 0 z/v z , which implies that γ 0 ≥ 0 in accord with causality [see also preceding remarks above Eq. (26)]. Otherwise, it would diverge far away from the surface.

B. Ω and G in the SCM
The conduction current density is also written in two parts, z). For small kv F /ω, we may keep only the first term in the series of F 0 (k,ω, v). It is then straightforward to show that where J SRM (z) is responsible for the extension made through the SRM beyond the DM. It is given by 53 where we have defined a short-hand together with these functions See that J SRM,z (0) ≡ 0, which makes no contribution to G. Now the total current density becomes By definitions (6) and (11), we find where F is an odd function ofω and given by Partially performing the integral, we obtain The real part of this expression approximates ω 2 p 1 + 3 5 K 2 /k 2 p in the long wavelength limit, corresponding to the HDM limit Obviously, the second term in Eq. (39) generally contains an imaginary part even in the collisionless limit where τ −1 is vanishingly small, due to a pole atω = K · v in the integrand of the integral in F. This part gives rise to Landau damping. For γ 0 /ω p ≪ 1, we find forω = ω s + iγ 0 Here P takes the principal value. The sign of the second line depends on the sign of γ 0 . Only for γ 0 > 0, it is negative leading to damping, in line with causality. Equation (42) shows that, for bulk waves Landau damping exists only for sufficiently large K. For SPWs, however, Landau damping always exists, because q runs over all positive values in the secular equation (13). A major improvement of the SCM over the SRM comes through the quantity G(K, ω). In the SRM and its descendents, G = G 0 contains no symmetry breaking effects. In the SCM, one finds G = G 0 + G s , where G s stems from J s (0) and is given by where ′ > ′ indicates that the integral is restricted to departing electrons, i.e. v z ≥ 0. Note that G s depends linearly on p.

C. The frequency and losses of SPWs
With Ω and G, we now proceed to solve Eq. (13) to find the frequency ω s and damping rate γ of SPWs in the SCM.
Analytical analysis. As said before, the most significant improvement of the SCM over the SRM is through the term G s . We wish to do an approximate analysis to explicitly demonstrate how G s = G ′ s + iG ′′ s affects ω s and γ. For this purpose, let us write G = G ′ +iG ′′ and Ω = Ω ′ +iΩ ′′ as well as ǫ = ǫ ′ + iǫ ′′ , and assume that γ 0 /ω s , G ′′ /G ′ and Ω ′′ /Ω ′ are all small quantities. Then, one can show that the real part of Eq. (13) gives which determines the SPW frequency ω s . The imaginary part of Eq. (13) yields γ 0 as where the contribution stems directly from the imaginary part G ′′ -which equals G ′′ s -of G, and comes directly from the imaginary part Ω ′′ of Ω. For stable systems, η 0 must be non-negative.
Note that Ω ′′ signifies Landau damping and inter-band absorption, as is clear from Eq. (42). As such, we may further split η 0 = η Landau + η interband , so that the SPW damping rate becomes where γ Landau,interband,s = ω s η Landau,interband,s , the first term represents Joule heating, the second and the third stand for Landau damping and inter-band absorption, respectively, while the last one is completely new due to G ′′ = G ′′ s . For models where G s = 0, e.g. the DM, the HDM and the SRM, this new term is absent.
For the SCM, however, G s is finite. Retaining only the first term in the series of F 0 , we find The second line here approximately corresponds to iG ′′ s . In general G ′′ s ≤ 0 and hence η s ≥ 0. This implies that symmetry breaking effects tend to counteract the conventional damping and destabilize the metal. Our numerical solutions shall demonstrate that γ 0 is non-negative, in accord with the general argument given in Sec. V A.
In the long wavelength limit k ∼ 0, one has k/K 2 ≈ πδ(q). Under this approximation, Eq. (44) becomes Here K 0 = (k, 0). For models with G s = 0, one immediately recovers from this equation the relation that ǫ d + ǫ ′ = 0 as expected. For the SCM, however, one finds instead This leads to This result differs considerably from what is expected of other models. It shows that the value of the SPW frequency depends on surface conditions.
If we replace in Eq. (46) ǫ ′ (K, ω s ) with its non-dispersive part, as is reasonable for small k, then we find which implies that inter-band transitions have little effect on η s , whereas a dielectric can enhance it by as much as 200%. This is to be borne out in numerical analysis in what follows.
Numerical analysis. We solve Eq. (13) numerically to find ω s and γ 0 and how they vary with k and ǫ d . We present the results for diffuse surfaces with p = 0 only. In Eq. (13), the integral over q extends to infinity. In numerical calculations, we impose a cutoff q c . Roughly, q c ∼ a −1 , where a is a lattice constant. For metals, this means q c ∼ k F ∼ k p . In all the numerical results presented here, we have chosen q c = 1.5k p . For q c beyond this value, both ω s and γ 0 quickly converge, confirming that the results are independent of the choice of q c 54 . Our results should be taken with a grain of salt for very small k ≪ ω s /c ≈ k p v F /c ∼ 0.01k p , where c is the speed of light in vacuum, because of retardation effects neglected in our theory.
The results are displayed in Figs. 2 and 3. In Fig. 2, we showω = ω s + iγ 0 as a function of k for various ǫ d but without inter-band transition effects (ǫ p = 0). As seen in Fig. 2 (a), in agreement with the analytical expression of ω s , increasing ǫ d leads to smaller ω s . Note that ω s is considerably larger than what would be obtained by other models due to surface effects. Meanwhile, γ 0 increases with increasing ǫ d , as seen in Fig. 2 (b), in accord with Eq. (52). This increase comes from the factor β, i.e. the presence of a dielectric enhances the electric field at the surface; see Eqs. (17) and (18).
The effects of inter-band absorption and Landau damping are illustrated in Fig. 3. Here we plot γ 0 for ǫ d = 1 under several circumstances as described in the figure. We see that inter-band transitions can strongly diminish γ 0 in two ways, as can be deduced from Eqs. (51) and (52). Firstly, there is the screening effect (the curve with ǫ p = 5). This leads to smaller ω s and hence smaller γ 0 , while leaving η 0 unaffected. Secondly, inter-band absorption further reduces γ 0 (see the curve with ǫ p = 5 + 0.5i). As for Landau damping, it is sizable and generally increases with k; see Refs. 5,12 . As such, γ 0 decreases as k increases. In relation to this feature, we should mention a size effect 11 : in films of thickness d, γ 0 is strongly suppressed and quickly diminishes to zero when the wavelength becomes longer than d. Echoing this, one can show that γ 0 ∼ kv F vanishes for k ∼ 0 in the SIM 46 .

VI. POSSIBLE LOSSLESS SPWS BY THE SCM
According to the SCM, the SPW loss rate is γ = τ −1 −γ 0 . In view of energy conversion, the expression implies a competition between the loss due to thermalization and the gain due to energy transferred from the electrons to the waves 12 . Should the condition γ 0 τ = 1 be fulfilled, lossless SPWs may be produced. In this section, we discuss this possibility against two common plasmonic metals: silver and aluminum.
Note that τ −1 is the collision rate at the SPW frequency ω s . Even at zero temperature and in defect-free samples, there is sufficient phase space -available due to the effective temperature ω s /k B -for electronic scattering and thus τ −1 is comparable to ω p . Up to our knowledge, there is virtually no direct data on τ −1 for any materials. We then opt to estimate it by the following formula 32 , where γ residual is the residual rate given by Here the first contribution comes from phonon scattering, the second from electron-electron scattering and the third due to impurity scattering to be neglected hereafter. The second term generally underestimates the electron-electron rate in noble metals and Al by a few times 55 . The coefficient ν 0 ∼ k B T D / may be determined from the slope of the phononic part of the D.C. conductivity at temperatures higher than the Debye temperature T D . Equation (54) shows that γ residual may well amount to quite a few percentages of ω p for metals. Silver. In silver electronic transitions involving the 5s and 4d bands have a dramatic effect on the properties of SPWs 35,56,57 , leading to ω s ≈ 3.69 eV and ω b ≈ 3.92 eV at long wavelengths, both lying far below the characteristic frequency ω p = 9.48 eV. Experimentally 32,58 , it was found that ν 0 ∼ 0.1 eV ∼ 0.01 ω p , in consistency with the value of T D ≈ 220 K for Ag. The electron-electron scattering rate according to Eq. (54) would be less than one percent of ω p , while experimental measurements and more accurate expressions 55 place it about 0.02ω p . As such, we may reasonably take γ residual ∼ 0.03ω p as an estimate.
The damping rate γ can be read out from the line shape of the electron energy loss spectra (EELS). The temperature dependence of γ has been recorded on a high-quality Ag single crystal by means of EELS 59,60 . The data indicates that at low temperature γ amounts to less than one percent of ω p , a value a few times smaller than the as-projected γ residual . If Landau damping and inter-band absorption are also counted, the discrepancy can be more dramatic. In light of the present theory, this discrepancy gives an estimate of γ 0 . The values suggest that γ 0 has substantially compensated for the collision losses, i.e. γ 0 ∼ γ residual , as borne out in the following analysis.
To evaluate γ 0 , ǫ p (ω) needs to be supplied. Combining a semi-quantum model and ellipsometry as well as transmittance-reflectance measurements, Rakić et al. 36 employed K-K analysis and prescribed a parametrized dielectric function -ǫ (b) r (ω) in their notation -for the inter-band contribution. Here, we use their fitting as an input for ǫ p (ω) but with a number of caveats. Firstly, their dielectric function was deduced from measurements assuming the conventional electromagnetic responses without any surface effects considered in the present paper. These effects, however, should be considered when analyzing ellipsometry and reflectance spectra. A future study will be made to address this issue. Secondly, their function very poorly reproduces the electron energy loss spectra and the reflectance spectra, especially near the SPW frequency of interest. Thirdly, their function does not give an accurate partition into inter-band and intra-band contributions, e.g. ω p ≈ 8eV was used rather than the widely agreed 9.48eV 57 , which may overestimate the inter-band transition effects. Finally, their function is defined only for real frequency, and thus in general not suitable for ǫ p (ω), where ω is complex. To remedy this point, we substituteǫ (b) r [Re(ω)] for ǫ p (ω), which should be reasonable if γ/ω s ≪ 1.
The results are displayed in Fig. 4 (a), where the computed γ 0 is exhibited as a function of k at various values of ǫ d . While γ 0 for SPWs supported on a pristine Ag surface (i.e. ǫ d = 1) is negligibly small and way below γ residual , by using a dielectric it can be significantly enhanced at long wavelengths beyond γ residual . This trend is consistent with Eq. (52). The fact that γ 0 can be made higher than γ residual suggests the possibility of compensating for the plasmonic losses completely in Ag. The situation is shown in Fig. 5. By interfacing the metal with a lossless dielectric of ǫ d = 5 and cooling it down toward a critical temperature T * ∼ 120K, one can diminish the net losses as much as desired. It should be mentioned that, ǫ d is the constant at the frequency ω s as well.
Aluminum. Inter-band transitions in Al are widely considered less pronounced than in Ag. Nevertheless, their presence can still be felt, e.g. in the difference between the values of ω p ≈ 12.6eV and ω b ≈ 15.3eV. These numbers were obtained by density functional theory 61  tal fitting 37 . In addition, ω s ≈ 10.7eV 62,63 . γ residual may be deduced from the experimental measurements performed by Sinvani et al. 64 and others 65 . These authors measured the low temperature dependence of the d.c. resistivity ρ of Al. Their data shows that ρ ≈ ρ 0 + AT 2 + B(T/T D ) 5 , where ρ 0 stems from impurity and lattice dislocation scattering while A and B are constants characterizing electron-electron scattering and electron-phonon scattering, respectively. Analyzing the data, the authors found that A ≈ 0.21pΩ·cm/K 2 (a lower bound) and B ≈ 4.9 · 10 4 µΩ·cm for T D = 430K. From this we obtain ν 0 ≈ 0.18ω p and the residual electronelectron scattering rate -a few times larger than what would be obtained with Eq. (54) 32 -approximating 0.14ω p , yielding γ residual ≈ 0.21ω p . It is noted that this value is comparable to the width (∼ 1.5eV, nearly 0.12ω p ) of the electron energy loss peak near the frequency of the bulk plasma waves in Al 62 . The as-obtained ν 0 (and hence γ residual ) represents probably an overestimate 66 .
To compute γ 0 , we again resort to the fitting function constructed by Rakić et al. 37 and have it in place of ǫ p (ω) in our theory, in the same way as we did in the case of Ag. It goes without saying that the same caveats should be kept in mind. The results are shown in Fig. 4 (b). As expected, γ 0 is comparable to that in the absence of inter-band transitions [see Fig. 2 (b)], as these transitions are weak in Al. As is with Ag, γ 0 of SPWs in Al can also be fortified -but to a lesser extentby a dielectric. Nevertheless, the enhanced γ 0 still falls short of γ residual unless for very long wavelengths where retardation effects need to be properly accounted for.
The calculations reported in the above have assumed p = 0. For surfaces strongly reflecting electrons (i.e. p > 0), γ 0 could be much lower 12 . We also point out that additional losses such as due to SPWs converted into radiation are not considered. They can be absorbed in the definition of τ −1 .

VII. CONCLUSIONS
In order to answer the two questions posed at the beginning of this paper, i.e. (i) why had not previous work hit at the possibility of lossless SPWs and (ii) how far is the latter from reality, we have derived a universal macroscopic theory of SPWs that applies to any electron dynamics. In light of the theory, our answer to question (i) is simple: lossless waves are possible only within a self-consistent description of physical surfaces that takes care of translation symmetry breaking effects, a condition not met in existing work. As for question (ii), we can only suggest an optimistic prospect rather than an answer due to various uncertainties in inter-band transition effects: our estimate shows that lossless waves may well be within the reach in some materials.
Our results reveal two contradictory views regarding SPW losses, as compared in Table I. According to the conventional framework, as exemplified by the DM, the HDM and the SRM, and thus the SPW loss rate cannot be smaller than either of τ −1 and γ Landau . On the other hand, within the SCM, a totally different picture emerges, giving which suggests that the loss rate is always smaller than τ −1 .
Here γ s = η 0 ω s , see Eq. (51). To directly contrast these two views, one has to measure separately τ −1 and γ. While the latter can be measured in many ways, the former is difficult to be directly measured. In what follows, we mention some indirect observations defying γ (1) but supporting γ (2) .
Firstly, we note that γ (2) naturally resolves a long-standing puzzle, that of the apparent insignificance of Landau damping even at very short wavelengths 38 . For example, in single crystal Ag, the loss rate measured by EELS 59,60 is ∼ 1%ω p even for k ∼ 1nm −1 , whereas γ Landau ∼ kv F ∼ 10%ω p . This discrepancy is inexplicable by γ (1) , but easily comprehensible by γ (2) , i.e. Landau damping has been overcompensated by γ s .
Secondly, we note that the loss rate of bulk plasma waves differs from that of SPWs primarily because of the absence of γ 0 in the former, as least in materials where inter-band effects are not important. In those materials, bulk waves should be generally much more lossy than SPWs, an observation that seems in consistency with experience. For example 67 , the loss rates for the SPWs and the bulk waves in potassium are 0.1eV −1 and 0.24eV −1 , respectively, while those in cesium are 0.23eV −1 and 0.75eV −1 , respectively. In spite of these, the general situation is obviously unclear at this stage.
Finally, we mention an experiment performed on a van der Waals structure by Iranzo et al. 68 . These authors were able to confine propagating plasmon between a graphene layer and a metal array to the atomic limit without sacrificing its lifetime, which obviously beats the limit set by Landau damping. From an energy conversion point of view 12 , the plasmon in such a structure is not much different from the surface plasmon on a metal surface. Their result is compatible with γ (2) : in ultimate confinement γ 0 tends to zero due to increase of Landau damping (Fig. 3), leaving the loss rate saturating at τ −1 , as observed. We anticipate a similar trend for the losses of local plasmon resonance on metal particles.
In the SCM we have assumed that the ground state of the underlying metal be simply the Fermi sea. The fact that γ can be made negative means an instability of the Fermi sea. Upon entering such circumstances, the metals are expected to undergo a transition into a different stable state, of which the electrical responses cannot be captured by our current SCM calculations. We will clarify the nature of this transition in the future.
The results reported in this work should be of broad interest to the researchers working in plasmonics, surface science and condensed matter physics. We hope that experimentalists will find the results fascinating enough to put their hands on them.
Acknowledgement -The author enjoyed the hospitality during his stay with K. Wakabayashi's group at Kwansei Gakuin University, Japan, where part of the writing was undertaken. This work is not supported by any funding bodies.

Appendix A: SPWs in the HDM
Here we show that the usually claimed SPWs in the HDM are incompatible with the waves in the DM. We put ǫ p = 0 and ǫ d = 1 for simplicity. In the HDM, the electrons are treated as a fluid described by two field quantities: the velocity field v(x, t) and the electron density field n 0 + n(x, t), where n(x, t) denotes the deviation from the mean density n 0 . The charge density is then ρ(x, t) = n(x, t)e and the current density is J(x, t) = e(n 0 + n(x, t))v(x, t), which in the linear responses regime becomes J(x, t) = n 0 ev(x, t).
A small fluid element of volume δV feels a force consisting of two portions: the electric force n 0 eδVE(x, t) and the pressure due to density variation −mδVv 2 0 ∂ x n(x, t). Now the laws of mechanics states that in the linear regime one has Here shear viscosity effects have been ignored. Now assuming n(x, t) = n(x)e −iωt and similarly for v(x, t) and other field quantities, we obtain the current density as the divergence of which is then given by Combining this equation with the equation of continuity, one obtains the wave equation for the charge density.
The usually quoted SPWs. In the standard but erroneous prescription for SPWs in this model, one takes J z (x 0 ) ≡ 0, or equivalently v z (x 0 ) ≡ 0.
Here x 0 denotes a point on the surface. The continuity equation reads In conjunction with Eq. (A3), one finds which is Eq. (23) given in Sec. II C. One then seeks solutions of the form ρ(x) = ρ(z)e ikx and similarly for other quantities.
With this it follows from Eq. (A8) that Combining this relation with Eq. (A6), we arrive at This is the usually quoted dispersion relation claimed for the SPWs in the HDM. As briefly captured in Sec. II C, this claim is plainly false: in the limit v 0 → 0, κ ∼ v −1 0 diverges and hence ∞ 0 dz ρ(z) = ρ 0 κ −1 ∼ ρ 0 v 0 vanishes, thus in contradiction with the DM. This would also erroneously imply that E z (z) did not change sign across the surface charge layer. It was essentially this erroneous solution that had been identified by Harris 14 and Garcia et al 15 in their study based on Boltzmann's equation. It is plausible that this is also so with works employing a more microscopic approach such as the density functional theory, at least those using the so-called 'infinite barrier' model for mimicking the surface 30 . For example, Feibelman identified a solution of uniform potential and hence empty of charges but with frequency ω p / √ 2 in the long wavelength limit 47 , exactly in this kind.