A Universal Self-Amplification Channel for Surface Plasma Waves
aa r X i v : . [ c ond - m a t . m e s - h a ll ] J a n A universal self-amplification channel for surface plasma waves
Hai-Yao Deng , , ∗ Katsunori Wakabayashi , , † and Chi-Hang Lam ‡ Department of Nanotechnology for Sustainable Energy, School of Science and Technology,Kwansei Gakuin University, Gakuen 2-1, Sanda 669-1337, Japan Department of Physics and Astronomy, University of Exeter, Stoker Road EX4 4QL Exeter, United Kingdom National Institute for Materials Science (NIMS), Namiki 1-1, Tsukuba 305-0044, Japan and Department of Applied Physics, The Hong Kong Polytechnic University, Hung Hum, Hong Kong
We present a theory of surface plasma waves (SPWs) in metals with arbitrary electronic collision rate τ − .We show that there exists a universal intrinsic amplification channel for these waves, subsequent to the interplaybetween ballistic motions and the metal surface. We evaluate the corresponding amplification rate γ , whichturns out to be a sizable fraction of the SPW frequency ω s . We also find that the value of ω s depends on surfacescattering properties, in contrast with the conventional theory. I. INTRODUCTION
Collective electronic oscillations on the surface of metals,dubbed surface plasma waves (SPWs) [1–3], have emergedas a pivotal player in nano scopic manipulation of light [4–8].The functionality of many prototypical nanophotonic devicescritically relies on the distance SPWs can travel before theyare damped out due to energy losses via several channels [7–12]. SPWs can lose energy due to Joule heat, inter-band ab-sorption, radiation emission and individual electronic motions(Landau damping). Most of the losses can be e ffi ciently butnot totally reduced under appropriate circumstances. Ampli-fiers have been contrived to compensate for the losses [13–21], which are all extrinsic and require external agents such asa dipolar gain medium to supply the energy.The standard theory of SPWs was formulated shortly aftertheir discovery in 1957 [1] and has been comprehensively dis-coursed in many textbooks [3, 5, 6, 22, 23]. In this theory,the electrical properties of metals are e ff ected by a frequency-dependent dielectric function ǫ . To analytically treat ǫ , thesimple Drude model or the slightly more involved hydrody-namic model [22, 24, 25, 27] is invariably invoked. For eithermodel to be valid, electronic collisions must be su ffi cientlyfrequent so that the electronic mean free path, l = v F τ , where v F is the Fermi velocity and τ the relaxation time, is muchshorter than the SPW wavelength or the typical length scaleof the system. The general case with arbitrary τ , especiallythe collision-less limit, where τ → ∞ , defies these modelsand has yet to be entertained. Other models based on ab initio quantum mechanical computations [22, 26] are helpful in un-derstanding the complexity of real materials but falls short inproviding an intuitive and systematic picture of SPWs under-pinned by electrons experiencing less frequent collisions.In the present work we employ Boltzmann’s transport equa-tion to derive a theory of SPWs in metals with arbitrary τ . Ouranalysis reveals a universal intrinsic amplification channel forSPWs. The existence of this channel does not depend on τ butis warranted by a general principle. We show that the unique ∗ [email protected] † [email protected] ‡ [email protected] interplay between ballistic electronic motions and boundariesresults in an electrical current that allows a net amount of en-ergy to be drawn from the electrons by SPWs, which would in-evitably self-amplify in the collision-less limit, thereby desta-bilizing the system. The hereby predicted self-amplificationof collision-less SPWs is analogous to the Landau dampingin collisionless bulk plasma waves [28–31]. While the latteris caused by slowly moving electrons that strip energy fromthe wave, the former is attributable to ballistically movingelectrons imparting energy to the wave when boundaries arepresent.In the next section, we describe the system under consider-ation, state our main results and conceive their possible exper-imental signature. In Sec. III, the SPW theory is presented. InSec. IV, we detail the methods used in the theory. We discussthe results and conclude the paper in Sec. V. Some calcula-tions and arguments not covered in Sec. III and IV are givenin appendices A and B. II. RESULTS
System - We consider a prototypical system, namely, asemi-infinite metal (SIM) occupying the half space z ≥ z = ε ( v ) = m v , where m and v denote the mass and velocity of the electrons, re-spectively. Inter-band transitions are accordingly ignored buttheir e ff ects will be discussed in Sec. V. The surface is of ahard-wall type and prevents electrons from leaking out of thesystem. Throughout we reserve r = ( x , y ) for planar coordi-nates while x = ( r , z ) denotes the complete position vector.We neglect retardation e ff ects in total. Key results - We find that the SPW frequency ω s and itsamplification rate γ can be written in the following form, ω s = ω s β ( p ) , γ = γ − τ , γ > . (1)Here ω s = ω p / √ / Drude model, with ω p being the characteristic fre-quency of the metal, p ∈ [0 ,
1] is a parameter that accounts forsurface scattering e ff ects, β ( p ) ∼ p and γ includes Landau damping and is exactly independent of τ .Both ω s and γ are determined by the secular equation (18)given in Sec. IV. Analytical expressions can be found for β ( p )and γ under certain approximations. Exact ω s and γ havebeen computed numerically and are displayed in Fig. 1, wherewe observe that γ /ω s ≈ α (1 + p ), with α ≈ . γ for SPWs. It also shows that, in contrast to the hydrody-namic / Drude model, ω s relies on surface properties via theparameter p . In particular, as seen in Fig. 1 (a), ω s is sig-nificantly – more than 10% – larger than ω s , the value ex-pected from the hydrodynamic / Drude model. In other words, β ( p ) is unity in these models whereas it could reach up to 1 . ω s and ω p . Experimental signature - We conceive one possible reper-cussion of the self-amplification channel, noting that the netamplification rate γ = γ − /τ and its temperature depen-dence can be directly measured in various ways, e.g. by ex-amining the SPW propagation distance or the width of the en-ergy loss peaks in electron energy loss spectroscopy (EELS)or the quality of the reflectance dips in the Kretschmann-Ottoconfiguration. In this paper, our calculation of γ is done atzero temperature. However, arguably γ could bear a dif-ferent – probably weak – temperature dependence than 1 /τ .In su ffi ciently pure samples, in which the residual scatteringis small enough, there might exist a critical temperature T ∗ ,above which γ < γ >
0, i.e. T ∗ marks atransition of the system from the Fermi sea to another state,see Sec. V for further discussions. In Fig. 2, we sketch the sit-uation described here, assuming 1 /τ is dominated by phononand impurity scattering.To experimentally investigate the self-amplification chan-nel using the most experimented materials like gold and silver,we must clarify the e ff ects of inter-band transitions, which notonly bring about additional energy losses but also strongly af-fect the value of ω s . We discuss these e ff ects further in Sec. V. III. THEORY
In this section, we give a systematic exposition of the the-ory that supports equation (1). Some technical aspects are leftto Sec. IV. The overall objective is to set up the equation ofmotion for the charge density, analyze its internal structureand solve it. We first discuss the relation between the chargedensity and the electric field, and then prescribes the currentdensity and show that it possesses a spectacular structure. Weproceed thence to the equation of motion, extracting its solu-tions and unveiling the properties of SPWs.
Electrostatics - We aim for establishing the equation of mo-tion for the charge density ρ ( x , t ) under the influence of its own electric field E ( x , t ). Here t denotes the time. Thanksto the linearity and symmetries of the system, we can write ρ ( x , t ) = Re h ρ ( z ) e i ( kx − ω t ) i and E ( x , t ) = Re h E ( z ) e i ( kx − ω t ) i with-out loss of generality. Here Re / Im takes the real / imaginarypart of a quantity, k ≥ ω isthe SPW eigen-frequency to be determined by the equationof motion. In general, ω can be complex. Instead of ρ ( z ), itproves more convenient to work directly with ρ q = Z ∞ dz ρ ( z ) cos( qz ) . In order for the jellium model to be valid, we must imposethat ρ q > q c =
0, where q c is a cut-o ff . q − c roughly gives themicroscopic lattice constant of the metal. We may choose q c ∼ k F = mv F / ~ ∼ n / ∼ k s = ω s / v F , where k F is the Fermiwavelength and n the mean electron density. The results donot depend on the exact value of q c as long as it is su ffi cientlylarge. By definition, ρ ( z ) = π Z ∞ dq ρ q cos( qz ) . In terms of ρ q , we can easily find, by solving Laplace’s equa-tion, that E x ( z ) = − i Z ∞ dq k ρ q k + q h qz ) − e − kz i , (2) E z ( z ) = Z ∞ dq k ρ q k + q h q / k ) sin( qz ) − e − kz i . (3)A snapshot of E ( x , t ) for a typical SPW is displayed in Fig. 3(a) and the magnitude of E ( z ) is plotted in (b). Of course, inmaking these plots ρ q has been determined by the equation ofmotion to be set up shortly. Current densities - Under E ( x , t ), an electrical current, ofdensity J ( x , t ) = Re h J ( z ) e i ( kx − ω t ) i , will flow. It can be calcu-lated using Boltzmann’s transport equation, see Sec. IV. As acrucial observation, we find that J ( x , t ) can always be writtenin two disparate contributions, J ( z ) = J D ( z ) + J B ( z ). Whatsets them apart is their distinct relations to E ( z ), as illustratedin Fig. 4. It turns out that J D ( z ) follows E ( z ) almost locally,as in the conventional hydrodynamic / Drude model – which isapplicable only when electrons execute di ff usive motions, i.e. τ is very small. Explicitly, J D ( z ) can be written J D ( z ) = i ¯ ω ω p π E ( z ) + J ′ ( z ) , (4)where ¯ ω = ω + i /τ , ω p = p π n e / m is the characteristicplasma frequency of the metal and J ′ ( z ) = Z ∞ dq ρ q F ( k , q ; ¯ ω ) k + q cos( qz ) , (5)signifies a non-local contribution that generates dispersiveplasma waves, with F ( k , q ; ¯ ω ) = (cid:18) m π ~ (cid:19) Z d v ( − e f ′ ) v ∞ X l = kv x + qv z ¯ ω ! l . p0.70.750.80.850.90.95 Numerical (exact)Equation (1)
Numerical (exact)Fitting by 0.09+0.1 p (a) (b) ω s / ω p γ / ω p0 FIG. 1. Numerically calculated SPW frequency ω s and its intrinsic amplification rate γ as a function of the Fuchs parameter p , by equation(18). k / k s = . q c / k s = .
5. The approximate analytical solutions are also shown. γ residual scattering rate TemperatureT*
FIG. 2. Possible instability caused by the self-amplification channelupon cooling down the system. At the critical temperature T ∗ , γ = /τ . Here the sketch assumes that 1 /τ is dominated by phonon andimpurity scattering. Here f ( ε ) denotes the Fermi-Dirac distribution and f ′ itsderivative, which is to be taken at zero temperature in thispaper. Only terms with odd l in the series contribute. Onecan show that J ′ z (0) ≡
0. The factor heading E ( z ) in Eq. (4)is recognized as the Drude conductivity. As this conductiv-ity is nearly imaginary, J D ( x , t ) points almost perpendicular to E ( x , t ), as seen in Fig. 4 (a).In contrast, J B ( z ) – the ’surface-ballistic’ current – has nosimple relation with E ( z ) and reflects genuine surface ef-fects: J B ( z ) would totally disappear if the surface were ab-sent. While J D ( z ) is essentially a bulk property, J B ( z ) orig-inates from the interplay between ballistic motions and thesurface. It consists of two parts, J B ( z ) = J B , emg ( z ) + p J B , re f ( z ) , where J B , emg (z) is contributed by electrons directly emergingat the surface while J B , re f ( z ) by reflected electrons. Here p is the Fuchs parameter, which gauges the probability that anelectron impinging on the surface gets reflected back. Bothcontributions arise from electrons moving away from the sur-face, v z ≥
0. In fact, J B , emg / re f ( z ) = (cid:18) m π ~ (cid:19) Z d v Θ ( v z ) e v e i ˜ ω zvz g B , emg / re f ( v ) . (6)where ˜ ω = ¯ ω − kv x and g B , emg / re f ( v ) = − e f ′ Z ∞ dq ρ q L emg / re f ( v x , v z , k , q ; ¯ ω ) k + q , (7)with L emg / re f ( v x , v z , k , q ; ¯ ω ) = k ( iv x − v z ) kv z ± i ˜ ω − q v z ± ˜ ω kv x )˜ ω − q v z . For p ≈ kv F / ¯ ω , it is easy to show that J B , x ( z ) ≈
0, implying that J B ( z ) points normal to the surface. This isillustrated in Fig. 4 (b). Positiveness of γ - Equation (6) strongly constrains thevalue of ¯ ω . Actually, it requires that Im( ¯ ω ) ≥
0; otherwise,the integral in this equation would diverge for any z , because e i ˜ ω zvz diverges exponentially for v z →
0. In Sec. IV and Ap-pendix A, this point is further elaborated. Since γ = Im( ω ) = γ − τ , γ = Im( ¯ ω ) , (8)the fact that Im( ¯ ω ) ≥ γ competing with the loss channel τ − . In what fol-lows, we show that ¯ ω is independent of τ − by means of theequation of motion. Equation of motion - The equation of motion can be ob-tained by relating ρ ( z ) and J ( z ) via the equation of continuity, − i ¯ ωρ ( z ) + ∇ · J ( z ) = − J z (0) δ ( z ) , (9) z ( i n un i t s o f k s ) ρ( z ) x (in units of k s ) -1 - VacuumMetal -6 0 6
EzEx (a) (b) FIG. 3. Sketch of SPWs supported on the surface of a semi-infinite metal. k / k s = . p =
1. (a) Charge density map (color) and electricfield map (arrows). (b) Plots of ρ ( z ) and E µ ( z ), where µ = x , z . ρ ( z ) is calculated by Eq. (16) while E µ ( z ) by (2) and (3). Both ρ ( z ) and E µ ( z )have been normalized. z ( i n un i t s o f k s ) - (b) 20 30 40 50 60 70x (in units of k s ) -1 z ( i n un i t s o f k s ) - (a) 20 30 40 50 60 70x (in units of k s ) -1 Diffusive current density Surface - ballistic current density
FIG. 4. Snapshots of (a) di ff usive current density J D ( x , t ) and (b) surface-ballistic current density J B ( x , t ) in SPWs. These two current densitiesare not discriminated by the value of τ , but by their dependence on the surface. J B signifies genuine surface e ff ects and would totally disappearwithout the surface, while J D is a bulk property. k / k s = . p =
1. Im( ¯ ω ) has been neglected to emphasize the di ff erences. which is derived and discussed in the next section. Here ∇ = ( ik , ∂ y , ∂ z ). Substituting J ( z ) into this equation andFourier transforming it, we arrive at Z ∞ dq ′ h H ( q , q ′ ; ¯ ω ) − ¯ ω δ ( q − q ′ ) i ρ q ′ = S ( ¯ ω ) , (10)where S ( ¯ ω ) = i ¯ ω J z (0) counts as a source term and H ( q , q ′ ; ¯ ω ) = Ω ( k , q ; ¯ ω ) δ ( q − q ′ ) + M ( q , q ′ ; ¯ ω )are the elements of a matrix denoted by H ( ¯ ω ), with Ω ( k , q ; ¯ ω ) = ω p + π ¯ ω k · F ( k , q ; ¯ ω ) k · k , k = ( k , q ) . (11)See that Ω ( k , q ; ¯ ω ) is an even function of ¯ ω and depends onthe length but not the direction of k . The matrix M ( ¯ ω ) arises from ∇ · J B ( z ) with elements M ( q , q ′ ; ¯ ω ) given in Sec. IV.In Appendix B, we show that M ( ¯ ω ) ≈ M Z , where M ∼ kv F /ω p is a constant and Z is the unity matrix will all elementsbeing one. It makes only a minor correction to the diagonalmatrix with elements Ω ( k , q ; ¯ ω ) δ ( q − q ′ ).The equation of motion (10) can now be re-cast in a com-pact matrix form, h H ( ¯ ω ) − ¯ ω I i ρ = S ( ¯ ω ) E , (12)where I is the identity matrix, ρ is a column vector collectingall ρ q and E is a column vector with all elements being one.For later use, the source term can be rewritten S ( ¯ ω ) = G t ( ¯ ω ) ρ = Z ∞ dq G q ρ q , G q = G ( k , q ; ¯ ω ) k + q , (13)where G ( ¯ ω ) is another column vector, t takes the transpose,and G ( k , q ; ¯ ω ) = G D ( k ) + G B ( k , q ; ¯ ω ) , with G D ( k ) = k ω p π arising from J D , z (0) and G B ( k , q ; ¯ ω ) = G B , emg ( k , q ; ¯ ω ) + pG B , re f ( k , q ; ¯ ω )from J B , z (0). Here G B , emg / re f ( k , q ; ¯ ω ) = i ¯ ω (cid:18) m π ~ (cid:19) × Z d v Θ ( v z )( − e f ′ ) v z L emg / re f ( v x , v z , k , q ; ¯ ω ) . (14)Equations (10) - (14) are exact and do not explicitly involve τ ,thereby concluding that ¯ ω is independent of τ . Bulk modes and localized modes - Without the surface wehave M ( q , q ′ ; ¯ ω ) ≡ S ( ¯ ω ) =
0. Equation (10) reduces to Ω ( k , q ; ¯ ω ) − ¯ ω = ρ q ,
0, which all are extendedmodes and therefore describe cosine bulk plasma waves. Ifwe discretize q in step dq , in total there are N c = q c / dq suchmodes. Solving this equation, we find the dispersion relation ω b ( k ) for these modes. Generally, Ω ( k , q ; ¯ ω ) could possess animaginary part due to a pole, located at ¯ ω = kv x + qv z , of theintegrand in the integral involved in F ( k , q ; ¯ ω ), giving rise tothe celebrated Landau damping. Neglecting this, one obtainsfor small | k | ω b ( k ) ≈ ω p + k · k v F ω p ≈ ω p , (15)which was well known from the hydrodynamic / Drude model.In the presence of the surface, equation (10) admits notonly extended modes, for which S ( ¯ ω ) =
0, but also local-ized modes, for which S ( ¯ ω ) ,
0. The number of total modescannot change and is still N c . The extended modes againrepresent bulk waves and for them equation (12) reduces to h H ( ¯ ω ) − ¯ ω I i ρ =
0. As shown in Appendix B, because ofthe constraint that S ( ¯ ω ) =
0, there are in total at most N c − ff ected by M ( ¯ ω ) and still givenby ω b ( k ).The missing bulk mode has been converted into a localizedmode satisfying S ( ¯ ω ) , ρ = S ( ¯ ω ) h H ( ¯ ω ) − ¯ ω I i − E . (16)Plugging this in Eq. (13), we obtain1 = G t ( ¯ ω ) h H ( ¯ ω ) − ¯ ω I i − E , (17) which determines ¯ ω and hence the SPW eigen-frequency ω .Upon omitting M ( ¯ ω ) from H ( ¯ ω ), this equation becomes1 = Z ∞ dq G ( k , q ; ¯ ω ) k + q Ω ( k , q ; ¯ ω ) − ¯ ω . (18)Let us write the solution as ¯ ω = ω s + i γ and hence ω = ω s + i γ with γ = γ − /τ . See that the solutions ± ω s + i γ occur together, in accord with the fact that ρ ( x , t ) is a real-valued field. In Eq. (18), Ω ( k , q ; ¯ ω ) is generally complex and γ automatically includes Landau damping. Hydrodynamic / Drude limits - The hydrodynamic model isrevisited if we replace in Eq. (18) G ( k , q ; ¯ ω ) and Ω ( k , q ; ¯ ω )with G D ( k ) and ω b ( k ) given by Eq. (15), respectively. If wefurther disregard the dispersion in ω b ( k ), the Drude model isthen recovered, in which case we immediately find ω s = ω s .In both models, ¯ ω is real. Approximate solutions - We can solve Eq. (18) approxi-mately. To the lowest order in γ /ω s , which is assumed to besmall, we may determine ω s by approximating the real part of(18) as follows1 ≈ Z ∞ dq (cid:2) G ( k , q ; ω s ) (cid:3) k + q Ω ( k , q ; ω s ) − ω s . (19)Substituting the so-obtained ω s in the imaginary part of (18),we find γ ω s ≈ − R ∞ dq Im[ G ( k , q ; ω s ) ] k + q Ω ( k , q ; ω s ) − ω s R ∞ dq Re[ G ( k , q ; ω s ) ] k + q Ω ( k , q ; ω s ) − ω s ω s Ω ( k , q ; ω s ) − ω s , (20)which can be brought into a rather simple form if we take Ω ( k , q ; ω s ) ≈ ω p and ω s ≈ ω s . Namely, γ ω s ≈ − R ∞ dq Im[ G ( k , q ; ω s ) ] k + q R ∞ dq Re[ G ( k , q ; ω s ) ] k + q =
12 Re (cid:2) J z (0) (cid:3) Im (cid:2) J z (0) (cid:3) . (21)In this approximation, Ω ( k , q ; ω s ) could have an imaginarypart only if q c > ω s / v F . Landau damping would be excludedfrom γ otherwise.To proceed, we need to evaluate G B ( k , q ; ¯ ω ). To the linearorder in k , we find G B ( k , q ; ¯ ω ) ≈ i ¯ ω e (cid:18) m π ~ (cid:19) × Z d v Θ ( v z ) f ′ v z q v z (1 + p )¯ ω − q v z + kv z (1 − p ) i ¯ ω . (22)Carrying out the integral gives G B ( k , q ; ¯ ω ) ≈ ω p π " k ( p − − i + p )4 v F ¯ ω q . (23)Thus, Re (cid:2) G ( k , q ; ω s ) (cid:3) ≈ k ( ω p / π )(1 + p ) /
2. Inserting this inEq. (19), we get ω s ≈ ω s q − p + . We see that ω s gen-erally depends on surface scattering e ff ects, in contrast withwhat is expected of the hydrodynamic / Drude model. Analo-gously, using Eq. (21) and q c v F ∼ ω s , we find γ ≈ αω s ,with α = / π . Landau damping has been excluded here, asthe approximation only takes the real part of Ω ( k , q ; ω s ). Numerical solutions - Equation (18) can also be exactlysolved numerically. A comparison with the approximate so-lution is displayed in Fig. 1. The agreement in the matter of ω s is satisfactory, while that for γ is not. The discrepancymight be because the approximate solution excludes Landaudamping. It should be emphasized that, our numerical solu-tions do not depend on the choice of q c , provided the latter isbig enough, i.e. q c ≥ k s . IV. METHODS
This section discusses further some technical aspects of thetheory. We set out with a discussion of the equation of con-tinuity in the presence of surfaces. Then we describe Boltz-mann’s equation and solve it to obtain the electronic distribu-tion functions. Thence we derive the dynamic equation for thecharge density.
Equation of continuity - We start with the equation of con-tinuity, (cid:0) ∂ t + /τ (cid:1) ρ ( x , t ) + ∂ x · j ( x , t ) =
0, which relates thecharge density ρ ( x , t ) and the current density j ( x , t ) in a genericway. Here the damping term − ρ ( x , t ) /τ is included to accountfor electronic collisions that would drive the system towardequilibrium. In the jellium model, ρ ( x , t ) appears when theelectron density is disturbed away from its equilibrium value n . The surface prevents any electrons from leaking out of themetal. Explicitly, we write j ( x , t ) = Θ ( z ) J ( x , t ), where Θ ( z ) isthe Heaviside step function. In so doing, we have treated thesurface as a hard wall and considered the fact that J ( x , t ) maynot vanish even in the immediate neighborhood of the surface– which is obviously the case with the hydrodynamic / Drudemodel. The equation of continuity becomes ∂ t + τ ! ρ ( x , t ) + ∂ x · J ( x , t ) = − J z ( x , t ) δ ( z ) , (24)where x = ( r ,
0) denotes a point on the surface and δ ( z ) isthe Dirac function peaked on the surface. Physically, the righthand side of Eq. (24) means that, charges must pile up on thesurface if they do not come to a halt before they reach it.Equation (24) reduces to equation (9) upon using the planewave form for ρ ( x , t ) and J ( x , t ). Boltzmann’s approach - We employ Boltzmann’s equationto calculate J ( x , t ) as a response to E ( x , t ). On the micro-scopic level, we may introduce a surface potential φ s ( x ) intothe equation to account for surface scattering e ff ects. The cor-responding surface field E s ( x ) = − ∂ x φ s ( x ) is peaked on thesurface and, complying with the hard-wall picture, may havean infinitesimal spread d s → + . Unfortunately, as φ s ( x ) canhardly be known and varies from one sample to another, thismicroscopic method is impractical and futile.Alternatively, surface scattering e ff ects can be dealt withusing boundary conditions [32, 34–36]. This is possible be-cause E s ( x ) acts only on the surface. In the bulk, the solu-tions – the electronic distribution function f ( x , v , t ) – to Boltz- mann’s equation can be uniquely determined up to some pa-rameters, which summarize the e ff ects of – while without ac-tually knowing – φ s ( x ). With translational symmetry alongthe surface, only one such parameter, namely, the so-called Fuchs parameter p , is needed for the simple specular scat-tering picture. Physically, p measures the probability thatan electron impinging upon the surface is bounced back. Asusual, we write f ( x , v , t ) = f ( ε ( v )) + g ( x , v , t ), where g ( x , v , t )is the non-equilibrium part of the distribution function in thepresence of E ( x , t ). The current density can then be calcu-lated by J ( x , t ) = ( m / π ~ ) R d v e v g ( x , v , t ), where e de-notes the electron charge. It is worth pointing out that, as g ( x , v , t ) is a distribution only for the bulk, the charge den-sity is not given by ˜ ρ ( x , t ) = ( m / π ~ ) R d v e g ( x , v , t ), i.e. ρ ( x , t ) , ˜ ρ ( x , t ). Actually, It is easy to see that ˜ ρ ( x , t ) satisfies( ∂ t + /τ ) ˜ ρ ( x , t ) + ∂ x · J ( x , t ) = ρ ( x , t ) is exactly the charges justlocalized on the surface. Electronic distribution function - Let us write g ( x , v , t ) = Re h g ( v , z ) e i ( kx − ω t ) i . In the regime of linear responses, Boltz-mann’s equation reads ∂ g ( v , z ) ∂ z + λ − g ( v , z ) + e f ′ ( ε ) v · E ( z ) v z = , (25)where λ = iv z / ˜ ω . In this equation, the velocity v is more ofa parameter than an argument and can be used to tag elec-tron beams. It is straightforward to solve the equation underappropriate boundary conditions (Appendix A). Naturally, wehave g ( v , z ) = g bulk ( v , z ) + g sur face ( v , z ) , where the bulk term would exist even in the absence of sur-faces whereas the surface term would not. Using the expres-sions for E ( z ), we obtain g bulk ( v , z ) = − e f ′ Z ∞−∞ dq ρ q k + q kv x + qv z ¯ ω − ( kv x + qv z ) e iqz , (26)which has the same form as one would expect for a bulk sys-tem. Here we have defined ρ q < : = ρ − q .As for g sur face ( v , z ), we find it with a subtle structure: itcan be written as a sum of two contributions, one of which, g D , sur face ( v , z ), has a single form for all electrons irrespectiveof their velocities while the other, g B , sur face ( v , z ), does not. Ex-plicitly, we have g D , sur face ( v , z ) = − e f ′ Z ∞ dq ρ q k + q k ( v z − iv x ) kv z + i ˜ ω e − kz . (27)We may combine g bulk ( v ) and g D , sur face ( v , z ) in a single term, g D ( v , z ) = g bulk ( v , z ) + g D , sur face ( v , z ) , in order to separate them from g B ( v , z ) : = g B , sur face ( v , z ) . In so doing, we have decomposed g ( v , z ) = g D ( v , z ) + g B ( v , z )in a di ff usive and a ballistic component. It is underlined that g B ( v , z ) arises only when the surfaces are present. For bulksystems without surfaces, it does not exist even if the elec-tronic motions are totally ballistic, i.e. τ → ∞ . As such, wecall it ’surface-ballistic’. g B ( v , z ) exists only for electrons leaving the surface, i.e. v z ≥
0. Those electrons could directly emerge from the sur-face or be those incident on but subsequently get reflectedby the surface. What fundamentally sets g B , sur face apart from g D ( v , z ) rest with its simple z dependence. Actually, we have g B ( v , z ) = Θ ( v z ) e i ˜ ω zvz h g B , emg ( v ) + p g B , re f ( v ) i , where the Fuchs parameter p gauges the fraction of reflectedelectrons and g B , emg / re f ( v ) has been given in Eq. (7). As al-ready remarked, since g B ( v , z ) ∝ e i ˜ ω z / v z ∝ e − γ z / v z , we musthave γ ≥
0; otherwise, it would diverge either when z → ∞ or for small v z . In Appendix A, we argue that this result isa consequence of the causality principle, which states that thenumber of reflected electrons is determined by that of incidentelectrons, rather than otherwise. Divergence of the current densities - The current density J ( z ) = ( m / π ~ ) R d v e v g ( v , z ) is split in two parts: J D ( z ) and J B ( z ), where J D / B ( z ) = ( m / π ~ ) R d v e v g D / B ( v , z ). Usingthe expressions for g D / B ( v , z ), it is straightforward to obtain J D ( z ) given by Eq. (4) and J B ( z ) by Eq. (6). To obtain theequation of motion (10) from the equation of continuity (9),the Fourier transform of ∇ · J D / B ( z ) is needed. We find Z ∞ dz cos( qz ) ∇ · J D ( z ) = i ¯ ω Ω ( k , q ; ¯ ω ) ρ q , (28)with Ω ( k , q ; ¯ ω ) given by Eq. (11). Additionally, Z ∞ dz cos( qz ) ∇ · J B ( z ) = i ¯ ω Z ∞ dq ′ M ( q , q ′ ) ρ q ′ , (29)where M ( q , q ′ ) is a matrix given by M ( q , q ′ ) = ω k + q ′ (cid:18) m π ~ (cid:19) × Z d v Θ ( v z ) (cid:16) − e f ′ (cid:17) i ˜ ω v z ˜ ω − q v z L ( v x , v z , k , q ′ ; ¯ ω ) , (30)where L ( v x , v z , k , q ; ¯ ω ) = L emg ( v x , v z , k , q ; ¯ ω ) + pL re f ( v x , v z , k , q ; ¯ ω ) . In Appendix B, we show that M ( q , q ′ ) generally makes a mi-nor correction, of the order of kv F /ω s , to Ω ( k , q ; ¯ ω ) δ ( q − q ′ ).Ultimately, the insignificance of this correction may be as-cribed to the factor e i ˜ ω z / v z in g B ( v , z ), which is extremely os-cillatory over z with a period less than ∼ v F /ω s .Substituting Eqs. 28 and (29) in the equation of continuity,we immediately obtain the equation of motion (10). V. DISCUSSIONS AND CONCLUSIONS
Thus, on the basis of Boltzmann’s equation, we have es-tablished a rigorous theory for SPWs in metals with arbitrary electronic collision rate 1 /τ . As a key consequence of the the-ory, we find that there exists a self-amplification channel forSPWs, which would cause the latter to spontaneously amplifyat a rate γ if not for electronic collisions. Surprisingly, thevalue of γ turns out to be independent of τ . The presence ofthis channel is guaranteed by the causality principle.Whether the system could actually amplify or not dependson the competition between γ and 1 /τ . If γ > /τ , SPWswill amplify and the system will become unstable. In our the-ory, the non-equilibrium deviation g ( v , z ) refers to the Fermi-Dirac distribution f ( ε ); as such, the instability is one of theFermi sea. Needless to say, the instability will be terminatedonce the system deviates far enough from the Fermi sea andsettles in a stable state. Clarifying the nature of the destinationstate is a subject of crucial importance for future study.One central feature of our theory is the classification of cur-rent densities into a di ff usive component J D ( z ) and a surface-ballistic component J B ( z ). This classification is not based onthe value of τ but according to whether the component obeysthe (generalized) Ohm’s law or not. Apart from this, thesecomponents are also discriminated in other ways. Firstly, theyare controlled by di ff erent length scales. As it largely followsthe local electric field E ( z ), the characteristic length associ-ated with J D ( z ) is k − . On the other hand, the length for J B ( z )is v F /γ , because of its simple z -dependence. Secondly, theyare oriented disparately. J D ( z ) is largely oriented normal to E ( z ) locally whereas J B ( z ) normal to the surfaces – especiallyfor p close to unity. Thirdly, J D ( z ) is a bulk property and ex-ists regardless of the surface; On the contrary, J B ( z ) reflectstrue surface e ff ects and it would disappear without surfaces.If we replace the vacuum by a dielectric ǫ , γ may be re-duced by an order of 1 /ǫ due to weakened E z ( z ) and J z (0).Roughly speaking, ρ ( z ) present on the metal side induces mir-ror charges amounting to ρ ′ ( z ) = − ρ ( − z )( ǫ − / ( ǫ +
1) onthe dielectric side. If ρ ( z ) is highly localized about the inter-face, as is with SPWs, ρ ( z ) and ρ ′ ( z ) will combine to give anet charge of 2 ρ ( z ) / ( ǫ + E z and J z (0) will bereduced by a factor of 2 / ( ǫ + γ andsmaller ω s = ω p / √ ǫ +
1. Studying dielectric e ff ects may beimportant in applications.Another problem that needs to be addressed for experimen-tal studies is concerned with the e ff ects of inter-band tran-sitions. In the most experimented materials, such as silverand gold, these transitions are known to have dramatic ef-fects. They not only open a loss channel due to inter-bandabsorption, but also significantly shift the SPW frequency.Including them in our formalism consists of a simple gen-eralization: in addition to J D ( z ) and J B ( z ), the total cur-rent density J ( z ) must now also have a component J int ( z ) ac-counting for inter-band transitions. The equation of motionis obtained by substituting J ( z ) in Eq. (9). One may write J int ,µ ( z ) = P ν R dz ′ σ µν ( z , z ′ ; ω ) E ν ( z ′ ), where µ, ν = x , y , z andthe inter-band conductivity σ µν can in principle be calculatedusing Greenwood-Kubo formula. In practice, calculating σ µν could be a formidable task even for the imaginably simplestsurfaces. Nevertheless, one may argue that J int ( z ) primarilya ff ects the properties of bulk waves, namely, Ω ( k , q ; ¯ ω ). Thecausality principle should still protect the amplification chan-nel. A systematic analysis will be presented elsewhere.We remark that, γ can also be calculated by studying thetemporal evolution of the electrostatic potential energy of thesystem. In particular, equation (21) can be directly derivedin this way. Detailed calculations along this line will be pub-lished in a separate paper.To conclude, we have presented a theory for SPWs takinginto account the unique interplay between ballistic electronicmotions and boundary e ff ects, from which it emerges auniversal self-amplification channel for these waves. It isexpected that the study will bear far-reaching practical andfundamental consequences, to be explored in the future. Acknowledgement – HYD acknowledges the Interna-tional Research Fellowship of the Japan Society for thePromotion of Science (JSPS). This work is supported byJSPS KAKENHI Grant Nos. 15K13507 and 15K21722 andMEXT KAKENHI Grant No. 25107005.
Appendix A: Electronic distribution functions
The general solution to Eq. (25) is given by g ( v , z ) = e i ˜ ω zvz C ( v ) − e ∂ v f mv z · Z z dz ′ e − i ˜ ω z ′ vz E ( z ′ ) ! , (A1)where C ( v ) is an arbitrary integration constant to be deter-mined by boundary conditions. We require g ( v , z ) = z → ∞ for any v . For electrons moving away from the sur-face, i.e. v z >
0, this condition is fulfilled for any C ( v ). Forelectrons moving toward the surface, i.e. v z <
0, we maychoose C ( v ) = e ∂ v f mv z · Z ∞ dz ′ e − i ˜ ω z ′ vz E ( z ′ ) , v z < , (A2)which leads to g ( v , z ) = e ∂ v f mv z · Z ∞ z dz ′ e i ˜ ω ( z − z ′ ) vz E ( z ′ ) , v z < . (A3)To determine C ( v ) for v z >
0, the boundary condition at z = p ( Fuchs parameter) of theelectrons impinging on the surface are bounced back in theabsence of E z ( z ), i.e. g (cid:16) ( v x , v y , v z > , z = (cid:17) = p g (cid:16) ( v x , v y , − v z ) , z = (cid:17) evaluated at E z ( z ) =
0. Then we get C ( v ) = − p e ∂ v f mv z · Z ∞ dz ′ e i ˜ ω ( z + z ′ ) vz E ( z ′ ) , v z > . (A4)By definition, p varies from zero to unity. Thus, we obtain g ( v , z ) = Θ ( v z ) g > ( v , z ) + Θ ( − v z ) g < ( v , z ) , where g > ( v , z ) = − "Z z e i ˜ ω | z − z ′ vz | + p Z ∞ e i ˜ ω | z + z ′ vz | e v · E ( z ′ ) | v z | ∂ f ∂ε dz ′ , g < ( v , z ) = − Z ∞ z e i ˜ ω | z − z ′ vz | e v · E ( z ′ ) | v z | ∂ f ∂ε dz ′ . (A5)Utilizing equations (2) and (3) for E ( z ) and carrying out the integral over z ′ , we find g >/< ( v , z ) = − ∂ f (cid:0) ε ( v ) (cid:1) ∂ε e (cid:12)(cid:12)(cid:12) v z (cid:12)(cid:12)(cid:12) Z ∞ dq k ρ q k + q g >/< ( v , z ; k , q ) , (A6)where g > ( v , z ; k , q ) = v z − iv x k + i ˜ ω/ v z e − kz + (cid:18) qv z qk + ˜ ω v x v z (cid:19) ( ˜ ω/ v z ) − q cos( qz ) + i (cid:16) qv x + ˜ ω qk (cid:17) ( ˜ ω/ v z ) − q sin( qz ) + iv x − v z k + i ˜ ω/ v z − (cid:18) ˜ ω v x v z + qv z qk (cid:19) ( ˜ ω/ v z ) − q + p iv x − v z k − i ˜ ω/ v z + (cid:18) ˜ ω v x v z − qv z qk (cid:19) ( ˜ ω/ v z ) − q e i ˜ ω zvz , (A7)and g < ( v , z ; k , q ) = iv x − v z k − i ˜ ω/ (cid:12)(cid:12)(cid:12) v z (cid:12)(cid:12)(cid:12) e − kz + (cid:18) q (cid:12)(cid:12)(cid:12) v z (cid:12)(cid:12)(cid:12) qk + ˜ ω v x | v z | (cid:19) ( ˜ ω/ v z ) − q cos( qz ) − i (cid:16) qv x + ˜ ω qk (cid:17) ( ˜ ω/ v z ) − q sin( qz ) . (A8)Now we can combine the terms with cos( qz ) and those with sin( qz ) into a single term called g bulk ( v , z ), while the rest into g sur face ( v , z ). Their expressions have been given in Sec. IV.If the surface is not uniform, we should have a function p ( r ) instead of a constant p . The boundary condition should then bewritten g ( x , ( v x , v y , v z > , t ) = p ( r ) g ( x , ( v x , v y , − v z ) , t ). In such case, the translational symmetry along the surface is generallylost and one cannot work with a single k -component any more, i.e. there is scattering e ff ects and di ff erent k -components aremixed. We do not consider this in this paper. Causality principle - In applying the boundary conditions to obtain C ( v ), we have implicitly assumed Im( ˜ ω ) ≥
0; otherwise,we would find unphysical solutions that violate the principle of causality, which states that the number of out-going electrons isdetermined by the number of in-coming electrons, not otherwise. It is easy to show that, had we assumed Im( ˜ ω ) <
0, we wouldhave found the opposite: the number of reflected electrons would be fixed while the number of incident electrons would go toinfinity as p →
0, which is unphysical.
Appendix B: The matrix M ( ¯ ω ) In the first place, we show that M ( ¯ ω ) /ω p ∝ ikv F / ¯ ω + ... , where the ellipsis stands for higher order terms in kv F / ¯ ω . For thispurpose, let us discretize q in step dq = π/ d , i.e. q l = l (2 π/ d ), where kd ≫ l = , , , ... takes integer values. Equation(10) can then be brought into the following form (cid:16) Ω ( k , q l ; ¯ ω ) − ¯ ω (cid:17) ρ l + X l ′ M l , l ′ ( ¯ ω ) ρ l ′ = S ( ¯ ω ) , (B1)where M l , l ′ ( ¯ ω ) = (2 π/ d ) M ( q l , q l ′ ; ¯ ω ) and ρ l = ρ q l . Writing R d v Θ ( v z ) = R π d ϕ R π/ d θ sin θ R ∞ dv ( v /
2) and integrating over v , we find M l , l ′ ( ¯ ω ) ω p = i π kd Z π d ϕ Z π/ d θ sin θ ¯ ω cos θ − ¯ ω kv F sin θ cos θ cos ϕ ( ¯ ω − kv F sin θ cos ϕ ) − q l v F cos θ k ¯ ω/ v F k + q l ′ L ( v F sin θ cos ϕ, v F cos θ, k , q l ′ , ¯ ω ) . (B2)To the lowest order in kv F / ¯ ω , we only need to retain L (0) in the expansion L sym = P ∞ m = L ( m ) ( kv F / ¯ ω ) m . Thus, L ( v F sin θ cos ϕ, v F cos θ, k , q , ¯ ω ) ≈ q v F cos θ ¯ ω − q v F cos θ (cid:0) + p (cid:1) . (B3)Substituting this back in (B2) and approximating¯ ω cos θ − ¯ ω kv F sin θ cos θ cos ϕ ( ¯ ω − kv F sin θ cos ϕ ) − q v F cos θ ≈ ¯ ω cos θ ¯ ω − q v F cos θ , q q + k ≈ , (B4)we arrive at M l , l ′ ( ¯ ω ) ω p = − i + p ) π kd kv F ¯ ω ! Z π d ϕ Z π/ d θ sin θ cos θ − ( q l v F / ¯ ω ) cos θ − ( q l ′ v F / ¯ ω ) cos θ . (B5)Obviously, we have M ( ¯ ω ) /ω p ∼ kv F / ¯ ω , as stated.We may proceed further If we take3 π Z π d ϕ Z π/ d θ sin θ cos θ − ( q l v F / ¯ ω ) cos θ − ( q l ′ v F / ¯ ω ) cos θ ≈ π Z π d ϕ Z π/ d θ sin θ cos θ ∼ , (B6)from which it follows that M ( ¯ ω ) l , l ′ ≈ M = − i ω p (1 / kd )( kv F / ¯ ω ), which is a constant. We then write M ( ¯ ω ) ≈ M Z , where Z l , l ′ = Ω ( k , q ; ¯ ω ) = ω p . Equation (B1) is nowwritten h ( ω p − ¯ ω ) I + M Z i ρ = S ( ¯ ω ) E . (B7)For bulk modes, S ( ¯ ω ) =
0. Note that Z has only one non-vanishing eigenvalue amounting to its dimension N c = q c d / π ∼ ( ω s / v F )( kd / π ). The corresponding eigenvector is ρ ∝ E , which obviously does not satisfy S ( ¯ ω ) = N c − Z are all zero andtherefore M ( ¯ ω ) has no impact on bulk modes.We can also show that M ( ¯ ω ) is negligible for the localized mode, for which S ( ¯ ω ) ,
0. To this end, we write (cid:20)(cid:16) ω p − ¯ ω (cid:17) I + M Z (cid:21) − = U − (cid:20)(cid:16) ω p − ¯ ω (cid:17) I + M ˜ Z (cid:21) − U , (B8)0where U is a similarity transformation that diagnolizes Z . We have used a tilde to indicate the transformed matrices, e.g. wewrite ˜ Z = U Z U − . As said, ˜ Z has only one non-vanishing element. Let it be the l -th element. Then ˜ Z l , l ′ = N c δ l , l δ l ′ , l . See that M ∼ / N and ˜ M l , l ′ ( ¯ ω ) ∼ − i ( ω p / π ) δ l , l δ l ′ , l . Introducing ˜ G t = G t U − and ˜ E = U E , we can rewrite the secular equation for thelocalized mode as 1 = ˜ G t (cid:20)(cid:16) ω p − ¯ ω (cid:17) I + M ˜ Z (cid:21) − ˜ E . (B9)Explicitly,1 = X l ˜ G tl ω p − ¯ ω + M ˜ Z l , l ˜ E l = X l G tl ω p − ¯ ω E l + ˜ G tl ω p (1 + i / π ) − ¯ ω ˜ E l − G tl ω p − ¯ ω E l ≈ X l G tl ω p − ¯ ω E l . The term in the square bracket makes only a contribution of the order of ∼ / N c and can be neglected for large N c . [1] R. H. Ritchie, Phys. Rev. , 874 (1957).[2] R. A. Ferrell, Phys. Rev. , 1214 (1958).[3] H. Raether, Surface plasmons on smooth and rough surfacesand on gratings (Springer Berlin Heidelberg,1988).[4] B. Rothenh¨ausler and K. Wolfgang, Nature , 615 (1988).[5] A. V. Zayats, I. S. Igor and A. A. Maradudin, Phys. Rep. ,131 (2005).[6] S. A. Maier
Plasmonics: fundamentals and applications (Springer Science & Business Media, 2007).[7] E. Ozbay, Science : 189 (2006).[8] M. L. Brongersma and P. G. Kik,
Surface plasmon nanophoton-ics (Springer, 2007).[9] W. L. Barnes, A. Dereux, and T. W. Ebbesen, Nature , 824(2003).[10] W. L. Barnes, J. of Optics A , S87 (2006).[11] T. W. Ebbesen, C. Genet and S. I. Bozhevolnyi, Physics Today , 44 (2008).[12] D. M. Giuliana, Y. Sonnefraud, S. K´ena-Cohen, M. Tame, S.K. ¨Ozdmir, M. S. Kim and S. A. Maier, Nano Letters , 2504(2012).[13] D. J. Bergman and M. I. Stockman, Phys. Rev. Lett. , 027402(2003).[14] J. Seidel, S. Frafstron and L. Eng, Phys. Rev. Lett. , 177401(2005).[15] I. De Leon and P. Berini, Phys. Rev. B , 161401(R) (2008).[16] I. De Leon and P. Berini, Nat. Photonics , 382 (2010).[17] D. Y. Fedyanin and A. Y. Arsenin, Opt. Express , 12524(2011).[18] P. Berini and I. De Leon, Nat. Photonics ,16 (2012).[19] D. Y. Fedyanin, A. V. Arsenin and A. V. Zayats, Nano Letters , 2459 (2012). [20] S. K´ena-Cohen, P. N. Stavrinou, D. D. C. Bradley and S. A.Maier, Nano Letters , 1323 (2013).[21] C. M. Aryal, B. Y.-K. Hu and A.-P. Jauho, arXiv:1508.01271(2015).[22] J. M. Pitarke, V. M. Silkin, E. V. Chulkov and P. M. Echenique,Rep. Prog. Phys. , 1 (2007).[23] D. Sarid and W. Cgallener, Modern introduction to surfaceplasmons: theory, mathematic modeling and applications (Cambridge University Press, Cambridge, UK, 2010).[24] J. Harris, Phys. Rev. B , 1022 (1971).[25] A. L. Fetter, Ann. Physics , 367 (1973); Phys. Rev. B ,3717 (1986).[26] P. J. Feibelman, Prog. Surf. Science. , 287 (1982).[27] O. Schnitzer, V. Giannini, S. A. Maier and R. V. Craster, Proc.R. Soc. A 472, 20160258 (2016).[28] L. Landau, J. Phys. USSR X , 25 (1946).[29] A. Y. Wong, R. W. Motley and N. D’angelo, Phys. Rev. ,A436 (1964).[30] R. Bingham and J. T. Mendonca and J. M. Dawson, Phys. Rev.Lett. , 247 (1997).[31] S. S. Natu and R. M. Wilson, Phys. Rev. A , 063638 (2013).[32] J. M. Ziman, Electrons and Phonons: the theory of transportphenomena in solids (Oxford University Press, 2001).[33] D. Pines,
Elementary excitations in solids (W. A. Benjamin,New York, 1963)[34] A. A. Abrikosov,
Fundamentals of the theory of metals (ElsiverScience Publishers B. V., North-Holland, 1988)[35] G. E. H. Reuter and E. H. Sondheimer, Proc. R. Soc. Lond. A , 338 (1948).[36] M. I. Kaganov, G. Y. Lyubarskiy and A. G. Mitina, Phys. Rep.288