A Theory of Dynamical Responses for Metal Films: Surface Roughness Effects

: A generic expression is derived for the dynamical response function of metal ﬁlms, with conductivity tensors as the only input. The semi-classical model is then used to provide an analytical expression for the conductivity tensor, thus establishing a kinetic theory for the response function. A major advantage of the theory is its ability to handle surface roughness effects through the use of the so-called specularity parameter. We applied the theory to study the properties of surface plasma waves. It is found that surface roughness does not affect the dispersion, but rather the decay rate of these waves. Furthermore, it signiﬁcantly affects the spectral weight carried by the SPW resonances, which diminishes toward zero as the specularity parameter approaches unity.

Electrodynamic responses can be quantified with the charge density-density response function, written as χ(x, x , t), which measures the amount of charges induced at point x by a transient electrostatic potential introduced a while t ago, localized at x .For infinite metals with full translational symmetry, χ depends on x − x rather than x and x individually, and its properties have been well known since the seminal work of Bohm and Pines in the 1950s [20][21][22][23].Their work provided a thorough understanding of plasma waves in metals, namely the charge density oscillations sustained by the long-range part of the Coulomb force between the charges.
For finite systems such as metal films, translational symmetry is partial, and the calculation of χ is much less straightforward [5,[24][25][26].Most earlier works were based on local responses, assuming that χ is significant only for x x , which is reasonable if the probing potential does not vary fast spatially.An important discovery from such work was that the charge density oscillations could also exist on metal surfaces, not just in the bulk [27][28][29][30].These oscillations, called surface plasma waves (SPWs) to be distinguished from the oscillations in the bulk (i.e., volume plasma waves, VPWs), are at the centre of on-going research in the name of plasmonics [31][32][33].
In order to go beyond local responses, the first and the simplest approach was to treat electrons in a metal like a fluid within the hydrodynamic model [34][35][36][37][38][39], which remains popular and has gained momentum in recent years [40][41][42][43].Another natural approach was to go microscopic and fully quantum mechanical [5,44], either employing the Greenwood-Kubo formalism [45][46][47] or the density functional theory (DFT) [44,[48][49][50][51][52][53][54][55].The Greenwood-Kubo community invariably works with an infinite barrier model [44], assuming that electrons are confined by an infinite potential barrier with wave function nodes at the surfaces, while the DFT way is in general computational, relying on supercomputers, and seems to be the favoured approach nowadays [44].DFT allows one to study microscopic effects for ideal surfaces such as the electron spill-out [48], but is not best suited to deal with ubiquitous features such as surface roughness [56].Yet, there is a semi-classical approach based on Boltzmann's kinetic equation, which was pursued by several authors in the 1970s but has since received undue attention [57][58][59].As a major advantage, this kinetic approach can deal with surface roughness analytically via a specularity parameter [60], and seems to be the only approach that has been used to account for the anomalous skin effect [61][62][63][64][65]. The kinetic approach is appropriate when the microscopic details of a surface are negligible in the problems of interest.
The main purpose of the present work is to derive generic analytical expressions for the density response function χ for metal films, and tailor it with the kinetic method.To date, as we are concerned, existing studies using this method have been exclusively on infinite or semi-infinite metals and none on films.As an innovation, our derivation treats the surfaces of a film as physical though negligibly thin regions rather than geometric separations [66][67][68][69].In this way we are able to avoid the imposition of additional boundary conditions permeating existing work [57][58][59].As an application, we apply the analytical expressions to study the properties of VPWs and SPWs in metal films.Although these waves have been well studied, little has been known of the impact from surface roughness, while revisiting well known results, we also see some unexpected features about these waves.
The rest of the paper is organized as follows.In the next section, we first derive a generic formalism for evaluating the response function with the conductivity tensor as the only input.We then discuss the macroscopic limit of this formalism.The result is illustrated with the local responses and the calculation of a screening function that describes how charges are induced in a film by exterior charges.Finally, we calculate the conductivity tensor using the kinetic equation and obtain the kinetic theory for the response function.It is shown that the response contains two types of resonances, one corresponding to the VPWs and the other to the SPWs.In Section 3, we analyse the properties of SPWs in light of the kinetic theory.We found that surface roughness, as characterized by the specularity parameter, has a strong effect on the spectral weight carried by the SPWs.Specifically, the spectral weight is reduced toward zero as the specularity parameter approaches unity.We also found that surface roughness has negligible effect on the SPW dispersion but considerably affects the decay rate.This can, in principle, be experimentally studied, for example in electron scattering.The paper is summarized in Section 4.

Analytical Expressions for the Dynamical Response Function
The system under consideration is a metal film placed in a vacuum, of thickness L with two surfaces located at z = 0 and z = L, respectively.The surfaces are macroscopically flat, but could be rough microscopically (i.e., on the scale of a Fermi wavelength) so that translational symmetry along the surfaces could be preserved macroscopically.The metal is treated in the spirit of the jellium model, by which the electrons are moving in a homogeneous and static background of positive charges that keeps the metal neutral as a whole.
In the presence of an external probing potential φ ext (x, t), electric charges of density ρ(x, t) = en(x, t) are induced in the metal.Here, e is the charge of an electron and n(x, t) is the deviation of the electron number density from the constant mean density n 0 .In the regime of a linear response, writing x = (r, z) with r = (x, y) being planar coordinates, or, equivalently, by a Fourier transform with respect to the planar coordinates and time, where k and ω represent the planar wave vector and frequency of the Fourier components, respectively.To avoid notation cluttering, we shall hereafter drop the scripts k and ω wherever no confusion is caused.In the above equation, both z and z lie inside the metal.
We may then introduce a cosine series for ρ(z) and φ ext (z).For example, where q n = nπ L are the wavenumbers and L n = L 2−δ n,0 , with δ n,0 being the Kronecker symbol.It may be noted that q n has a natural cutoff: ρ(z) cannot vary faster than allowed by the largest possible wavenumber in the first Brillouin zone of the underlying crystal lattice, which is comparable to ω p /v F .Thus, the cutoff in n is of the order of Lω p /v F .In other words, ρ n is negligible for n beyond this value.Now, Equation (1) can be rewritten as where φ ext,m is the cosine component of φ ext (z) and is the key quantity to be calculated.
In the rest of this section, we shall first describe a generic density response theory with the electrical conductivity tensor as the only input.Then, we apply the theory to the specific case for which the conductivity is obtained with Boltzmann's kinetic equation.We shall also derive the screening function describing screening effects.

Generic Formalism
We begin the construction of the theory with a discussion of the electric current flowing in the metal, as electric currents are the ultimate reason for induced charges.In the first place, we may distinguish between two types of electric currents [70]: one of them is driven by short-range thermal electronic collisions, while the other is driven by the long-range mean electric field generated by the probing potential and the induced charges.The electric current density is then written as j f ield + j col , where j f ield is due to the mean field and j col due to electronic collisions.As the collisions tend to drive the system toward thermodynamic equilibrium, the main effect of j col is to smear the non-equilibrium electronic distribution and cause any induced charges to decay.Within the relaxation-time approximation, we may assume that ∂ x • j col (x, t) = ρ(x,t) τ , where 1 τ gives the collision rate.This result is consistent with and can be directly derived from the kinetic theory [71].
Meanwhile, in the regime of linear responses, we may write j f ield = j + j ext , where j is driven by the electric field E of the induced charges and j ext by that of the probing potential.
The currents and electric fields are related by the electrical conductivity tensor Σ µν , where µ, ν = x, y, z label the components.Namely, we have, component-wise, Or, after a Fourier transform with respect to the planar coordinates and time, The scripts of wave vector and frequency have been dropped here as before.It should be noted that the conductivity Σ µν excludes the effects due to the long-range electronelectron Coulomb interaction supposedly included in the mean field E. The electric fields are obtained from the law of electrostatics, where ∇ = (ik, ∂ z ) and, with k = |k|, is the electrostatic potential of the induced charges.Equations ( 5)-( 7) are closed by the equation of continuity, Within the aforementioned relaxation-time approximation and after Fourier transform, the equation becomes where ω = ω + i/τ.Now, substituting the expressions of the currents and eliminating the electric fields in favour of the charge density, we find where is an integral kernel.Here ∇ = (ik, ∂ z ).Equation ( 9), which is linear, is formally solved by an inversion, yielding where Λ(z, z ) inverts the integral in the left-hand-side of Equation ( 9), satisfying From Equation ( 11) it follows that which gives the density response function in terms of the electrical conductivity as the only input.

Macroscopic Limit
The expression for χ(n, m), Equation ( 12) is generic and, in principle, exact within the mean field framework.However, to make progress, one must compute the conductivity tensor Σ µ,ν , which is not straightforward for a finite system.The difficulty, to be further explained below, lies with the surfaces.Here, we argue that this difficulty can be bypassed, largely at the cost of introducing some surface characteristics, in the macroscopic limit, whereby the surfaces can be considered infinitesimally thin.Needless to say, many theories, such as the currently popular hydrodynamic theory and the kinetic theory to be studied in the next section, are macroscopic.The macroscopic approach can be justified if the fields do not vary appreciably on the scale of the thickness of surfaces (usually a few atomic layers thick), as often implicitly presumed in the subject of electromagnetism [8].
In macroscopic theories, a surface is usually reduced as a geometric separation and its physical attributes have been habitually implemented with some boundary conditions [8].Such a macroscopic view tends to obscure the aforementioned difficulty, which, however, becomes conspicuous atomistically, e.g., in microscopic theories such as those based on DFT or the Greenwood-Kubo formula.Atomistically, a film normally divides itself into three regions: two inhomogeneous surface regions (SRs), one for each surface, and a homogeneous bulk region (BR).The SRs differ from the BR due to surface reconstruction, and the atomistic environment seen by electrons in the SRs, and therefore differs from that in the BR.In microscopic theories, such division is not needed, as the SRs and the BR are treated on equal footing from the outset and the exact Σ µν (z, z ) is computed, which applies to the entire system including both the SRs and the BR.Nevertheless, to perform the computation, one must provide a detailed atomistic profile of the whole system.In the BR, the profile is supposed to be the same as for an infinite system.In the SRs, however, the profile is expected to vary from sample to sample and is virtually unknown a priori.First-principles microscopic calculations often prescribe an ideal atomistic profile for the SRs and ignore realistic complications.Such microscopically computed Σ µν (z, z ), though applicable to the whole system, does not necessarily capture all the interesting aspects of realistic surfaces, in particular those due to surface roughness.
In macroscopic theories, on the other hand, the SRs are regarded as infinitesimally thin geometric separations and completely disregarded in the dynamical equations governing the conductivity tensor, and hence only the BR is directly (and often unwittingly) considered and the conductivity tensor obtained in macroscopic theories, denoted by σ µν (z, z ), applies only to the BR.In other words, σ µν (z, z ) and Σ µν (z, z ) are equal if both z and z lie in the BR, but they differ if z or z or both lie in the SRs.It follows that the actual electric current densities, j µ (z) and j ext,µ (z), determined by Σ µν (z, z ), are generally not the same as which are determined by σ µν (z, z ).Nonetheless, in the macroscopic limit where the electric fields do not vary appreciably over the SRs, j µ (z) J µ (z) and j ext,µ (z) J ext,µ (z) for z in the BR.Indeed, the BR the discrepancy between j µ (z) and J µ (z), and similarly between j ext,µ (z) and J ext,µ (z), is given by ∑ ν SR dz Σ µν (z, z ) − σ µν (z, z ) E ν (z ), which represents a higher-order correction ∼ d s /Λ that vanishes in the macroscopic limit.Here, d s is the SR thickness and Λ a macroscopic length.In the SRs, however, the discrepancy between j µ (z) and J µ (z) amounts to ∑ ν SR+BR Σ µν (z, z ) − σ µν (z, z ) E ν (z ), which cannot be neglected even in the macroscopic limit.In particular, it is worth pointing out that, at the microscopic division between the metal and the vacuum, Σ µν (z, z ) and j z (z) must vanish, as no electrons can flow into the vacuum (tunnelling ignored), but σ µν (z, z ) and J µ (z) do not need to, as these are valid only in the BR.To summarize, we may introduce some phenomenological functions w µ (z) so as to write j µ (z) = w µ (z)J µ (z) and j ext,µ (z) = w µ (z)J ext,µ (z).By definition, w µ (z) vanish at the microscopic divisions and cross over to unity as z goes over a SR into the BR.In the macroscopic limit, w µ (z) invariably reduces to the simple Heaviside step functions θ(z), i.e., w µ (z) → θ(z) − θ(z − L), all microscopic variations within the SRs thereby erased as they should in macroscopic theories.As such, in macroscopic theories, This simple expression dictates how j(z) and j ext (z), which are the actual current densities throughout the entire system, are related to J(z) and J ext (z), which give the current densities correctly only in the BR.The values of the latter in the SRs are extrapolations.In the literature, these two sorts of densities are often confused.For example, most authors impose that J z (z) vanishes at the surfaces, exemplifying the so-called additional boundary conditions first proposed by Pekar [72].Such an imposition is not justified in light of the above analysis; see also Refs.[66,69,73] and references therein.
As aforesaid, σ µν (z, z ) is governed by the same dynamical equations as for infinite systems.However, this does not mean that it is the same as the conductivity for an infinite system.Actually, σ µν (z, z ) contains an additional contribution standing for surface scattering effects, which in general carries some parameters characterizing such effects that can be regarded as fingerprints of a surface.In quantum mechanical theories (e.g., DFT based), the dynamical equation is the Schrödinger equation, and the parameters correspond to the scattering matrix elements with the surface as the scatterer.In the kinetic theory covered in the next section, the dynamical equation is Boltzmann's transport equation and the parameter is the specularity parameter (also called the Fuchs parameter), specifying the fraction of incident electrons that are specularly reflected by a rough surface; see details below.
With the relation (14), the density response function can be expressed in terms of σ µν (z, z ).After some straightforward manipulations, we obtain or compactly, treating χ as a matrix, where I is the identity matrix and H, G and A are square matrices with elements given by where 2 is a shorthand standing for ∑ µ,ν which is proportional to the potential generated by the m-th component of ρ(z).Note that Φ m (z) and cos(q m z) behave in the same way under reflection about the the mid-plane of the film, that is, As for A, we can split A = A (1) + A (2) with A Note that both G and A (2) originate from the derivatives taken in h(z, z ) and Equation (12) of the step functions in Equation ( 14).Had we mistaken σ µν (z, z ) for Σ µν (z, z ), these quantities would not have existed.They represent surface capacitive effects (viz., effects due to charges occupying the SRs) [67].
The wave equation is satisfied by the charge density, and Equation ( 9) can now be rewritten as a matrix equation, where ρ and φ ext are column vectors with elements ρ n and φ ext,n , respectively.The normal modes of the charge density waves (e.g., VPWs and SPWs) in a film are then determined from the secular equation It is instructive to note the differences between an infinite system and a film.There are two basic differences.Firstly, for an infinite system, G (and also A (2) ) vanishes as there is no surface capacitive effects.Secondly, as stressed before, the conductivity tensor for an infinite system, which we may denote by σ inf,µν (z − z ), differs from σ µν (z, z ) for a film: σ inf,µν (z − z ) does not bear any surface scattering effects, while σ µν (z, z ) does.As a result, H (and A (1) ) for a film differs from that for an infinite system, which we denote by H inf .The corresponding secular equation for an infinite system reads ω2 I − H inf = 0.In the so-called dielectric approximation [74], the difference between H and H inf is ignored.As will be seen later, this difference is crucial for obtaining the correct frequencies of SPWs in the kinetic theory.
Hereafter, we study films symmetric under reflection about its mid-plane z = L/2.The symmetry requires that σ µν (L − z, L − z ) = σ µν (z, z ) for (µ, ν) = (x, x), (y, y), (z, z), (x, y) and (y, x) whereas σ µν (L − z, L − z ) = −σ µν (z, z ) for (µ, ν) = (x, z), (z, x), (y, z) and (z, y).With this, one can show that H n,m vanishes identically unless n − m is even (i.e., n and m of the same parity).The same can be shown for G and A (both A (1) and A (2) ).Physically, this means that the symmetry of ρ(z) must be the same as that of φ ext (z): a symmetric (anti-symmetric) φ ext can only induce a symmetric (anti-symmetric) ρ(z).Therefore, we can decompose the quantities into two disconnected sectors, a symmetric one (indicated by '+') and an anti-symmetric one ('−'), obtaining and, similarly, G = G + ⊕ G − and A = A + ⊕ A − .Here, ρ + is a column vector with elements ρ + l = ρ 2l and ρ − is a column vector with elements ρ − l = ρ 2l+1 , where l = 0, 1, ..., H + is a square matrix with elements H + l,l = H 2l,2l while H − is a square matrix with elements H − l,l = H 2l+1,2l +1 , and similarly, For the sake of completeness, one may also note that where ρ ± (z) = ∑ ∞ l=0 ρ ± l cos(q ± l z) with q + l = q 2l and q − l = q 2l+1 .Further transformations can be performed by noting that, using the expression (17), and, with In matrix form, G ± = Z ± G ± , where Z ± is a column vector with elements Z ± l and G ± is a row vector with elements G ± l .In terms of G ± , we find where the function s is given by Note that the second term in Equation ( 24) does not exist for infinite systems and is unique to bounded systems (capacitive effects).The density response function of an infinite system, denoted by χ inf , has the same form as the first term in Equation ( 24), but with H ± and A ± replaced by the counterparts H inf and A inf , respectively.Namely, A inf .For an infinite system, there is no difference between the symmetric and the anti-symmetric sectors.
Equation (24) shows that χ ± possesses two types of resonances.One type of resonance occurs where In infinite systems, these resonances give rise to the VPWs.In films, depending on the dynamics, surface scattering effects might convert some of these resonances into surface localized ones.Within the kinetic theory, though, all of them represent VPWs, analogous to what happens in infinite systems.
The other type of resonance occurs where where ω p is the characteristic frequency of the metal.Upon using this expression, we obtain δ l,l .After some algebra, one gets the density response function χ ± (l, l ) for this model, which we write as Here, is the SPW frequency by this model, which is often quoted in the literature and textbooks.Clearly, χ ± LDM (l, l ) has a resonance at ω = ω p , which corresponds to the VPWs.It also has resonances at ω = ω + s0 and ω = ω − s0 , representing the symmetric and anti-symmetric SPWs, respectively.
Screening function.As another application, let us study the screening response of a metal film to a distribution of charges placed outside the film.We denote by ρ ext (z) the (ω, k)-th Fourier component of the density of these charges, which we assume all reside in the half space z < 0 without loss of generality, i.e., ρ ext (z > 0) ≡ 0. The potential generated by this component is given by ) n e −kL .Now, the corresponding Fourier component of the density of the charges induced in the film reads The screening function P n determines the response of the film to exterior charges.It plays a fundamental role in electron energy loss spectroscopy and electron scattering with metals.

Kinetic Theory
Here, we assume semi-classical dynamics for the electrons to obtain σ µν (z, z ) and then the corresponding χ ± through the formalism established in the preceding section.
The semi-classical distribution function for the electrons may be written as , where v is the velocity of electrons, f 0 (ε) is the equilibrium distribution (i.e., the Fermi-Dirac distribution) with ε = mv 2 /2 being the energy of the electrons of effective mass m, and g(v, z) is the non-equilibrium part of the distribution.As before, we have dropped the Fourier frequency and wave vector labels.No position dependence of f 0 has been assumed, as we are concerned with the electronic distribution in the BR, wherein g(v, z) satisfies the following Boltzmann equation, where is the electric field present in the BR.Following Fuchs [60,75,76], we assume that a fraction of electrons impinging on the surfaces of the film are specularly reflected.With for the surface at z = 0 and for the surface at z = L. Here, p 0 and p L are the Fuchs parameters, which give the fractions of specularly reflected electrons at the corresponding surfaces and measure microscopic surface roughness seen by electron waves resolving structural variations over a Fermi wavelength.Solving Equation ( 28) subjected to such considerations, we find Here, D = 1 1−p 0 p L e 2i ωL/|vz | expresses the effect due to electrons bouncing back and forth between the surfaces.Now, the current density can be calculated as usual and the conductivity tensor is obtained as where, as indicated by the greater symbol '>', the integral is restricted to v z ≥ 0, and the function Γ µν (z, z ) can be split in three parts, Here, we have assumed that p 0 = p L = p [and hence D = 1/(1 − p 2 exp(2i ωL/|v z |))], and the three parts Γ 0 µν (z − z ), Γ µν (z − z ) and Γ (2) µν (z, z ) are given as follows, which describes responses due to 'out-going' electrons (i.e., moving away from the perturbation caused by the electric field at z ) excluding specular surface scattering effects, and together with describe specular scattering effects.Γ 0 µν (z − z ) preserves translation symmetry, which alone would revisit the conductivity for an infinite system by Equation (30).Γ µν (z − z ) also preserves translation symmetry, which superimposes an out-going electron wave propagating in one direction on top of a specularly reflected wave travelling in the other direction.Γ (2) µν (z, z ) depends on z and z separately and breaks the symmetry.Now, the quantities, H ± , G ± and A ± in the density response function χ(n, m) given by Equation (15), are obtained by substituting the conductivity Equation (30) in Equations ( 16)- (19).Since Γ µν (z, z ) [and σ µν (z, z )] consist of three parts, we accordingly write each of these quantities in three parts.After some straightforward but tedious calculations, we find where stands for m 2πh 3 > d 3 v(−e 2 f 0 ), and Here, we have introduced Note that only the first term in Equation ( 36) contributes a diagonal part to H ± , while the rest represents scattering effects on the VPWs.In a previous work [67], these terms were ignored.
In addition, we find that where and In the above, As a sanity check, we have verified that G ± l vanishes identically for p = 1.In the limit L → ∞, one may show that G ± l ∝ 1 − p and revisits the expression given in Ref. [67].Finally, we have A ± l,l = A (1)± l,l + A (2)± l,l , where with and A (2)± l,l = Z ± l A ± l , where with In matrix form, we may write A (2)± = Z ± A ± , where A ± is a row vector with elements A ± l .Such a structure is already clear in Equation (19).As a sanity check, we have verified that A ± vanishes identically for p = 1.
Before proceeding to look at plasma waves, we remark on a few properties of the matrices H ± , G ± and A ± .In the first place, we may split where H ± (l, l ) represents scattering effects on the VPWs and with Q being the magnitude of Q.The diagonal part was derived from the first term in Equation (36).One may easily show that Ω 2 (Q, ω) ≈ ω 2 p for small Q, but it tends to zero as Qv F / ω increases.A peak (dip) develops around Qv F / ω ∼ 1 in the real (imaginary) part of Ω(Q, ω).This is numerically verified in Figure 1, where we also see that H ± amounts to a small perturbation with elements much less than ω 2 p .Despite its apparent smallness, H ± has a strong impact on the behaviours of SPWs and should not be ignored; see below.We have numerically computed G ± l as well.A typical result is displayed in Figure 2, where it is seen that G ± l features a peak around l ∼ L ω/v F and tends to zero for l beyond.These features are reminiscent of Ω(K ± l , ω).

Application to Plasma Waves
In this section, we examine the properties of VPWs and SPWs using the above derived kinetic theory.We first look at VPWs.These are represented by the resonances determined by ± b (k, ω) = ω2 − H ± ≈ 0. Their frequencies and decay rates are obtained from the eigenvalues of H ± .Since H ± is mostly diagonal with diagonal elements Ω 2 , the resonances are then located where Ω 2 (K ± l , ω) ≈ ω2 .In Figure 3, we map out log Det   As for SPWs, their dispersion and decay rates are determined by ± s (k, ω) = 0. To visualize them, we numerically compute the quantity and map it out in the ω-k plane, which features a peak at each SPW resonance for given k.The peak position gives the SPW frequency ω ± (k), while its half width at full maximum gives the decay rate ν ± (k).The maps are displayed in Figure 4 at various values of p, and example peaks are shown in Figure 5.The SPW resonances disperse along the bright lines in the map.It is seen that the symmetric SPWs and the anti-symmetric ones are nearly degenerate, i.e., ω + (k) ≈ ω − (k) for large kL 1, whereby the surfaces of the film are decoupled.As k decreases toward zero while still larger than 1/L, both ω + (k) and ω − (k) tend to the universal value ω p / √ 2 regardless of the value of p. Indeed, ω ± (k) exhibits virtually no dependence on p.For kL > 1, ω ± (k) increase linearly with increasing k.As well known, this stands in contrast with the LDM, by which the frequency would approach the universal constant.(L = 100v F /ω p ) alongside a thick one (L = 500v F /ω p ).One may note that, for a not-sosmall k, the FWFM increases linearly with k and the resulting ν ± can be well fitted by ν ± (k) ≈ ak + b.Here, both the coefficient a > 0 and the residual b vary with p.In Figure 6, one may see that a decreases with p, yet b increases with p.However, b seems to be bounded by 1/τ from below, i.e., ν ± (k) ≥ τ −1 .Such a linear increase is a characteristic of Landau damping [77,78].As k decreases toward zero, ν ± (k) starts to deviate from the linear dependence and rises up to a high value, which may be due to the strong coupling between the surfaces.For thicker films, such coupling occurs for smaller k and, hence, the linear dependence persists to a lower value of k, as can be seen by comparing the upper two rows of panels to the lower two in Figure 6.In previous works [67,71,73,79,80], one of us studied SPWs for a semi-infinite metal with a similar method.It was shown there that surface roughness could alter the SPW frequency, in disagreement with the present finding.This discrepancy occurs because in those works the term H ± in Equation ( 48) was neglected.The present work shows that H ± is essential to obtain the correct SPW frequency, even though its elements are small.Moreover, H ± represents surface scattering effects and contributes to the decay of SPWs.Without this term, which is indeed negligible for VPWs, the decay rate ν ± could go below 1/τ for a small k, as shown in those works.

Summary
In summary, we have provided a generic formalism for evaluating the charge densitydensity response function for a metal film with the conductivity tensor as the only input.We then applied the semi-classical model to calculate the conductivity, thereby establishing a kinetic theory for the response function.The main advantage of the kinetic theory is that it enables a simple way to study the effects of surface roughness through the Fuchs parameter p.
As an application, we employed the theory to analyse the properties of SPWs and VPWs.We found that the SPW dispersion does not depend much on p, whereas the decay rate of SPWs increases mildly as p increases.The most significant effect of surface roughness is seen in the spectral weight of the SPW resonances, which is drastically reduced as p

Figure 2 .
Figure 2. Numerically computed G + l (a) and G − l (b).Parameters are the same as for Figure 1.
2 (K ± l ,ω) in the ω-k plane at various values of the Fuchs parameter p.Since this is a sum of the contributions from all resonances with the same k, it is not possible to discern individual resonances unless L is extremely small so that the resonances are well separated.The VPW resonances are largely concentrated around ω = 0 and ω = ω p , since Ω(K ± l , ω) sits around ω p for small l but goes to zero for large l.