Table Of Content

A universal self-amplificationchannel forsurfaceplasma waves Hai-Yao Deng1,2, Katsunori Wakabayashi1,3, and Chi-Hang Lam4 ∗ † ‡ 1Department of Nanotechnology forSustainable Energy, School of Science and Technology, Kwansei Gakuin University, Gakuen 2-1, Sanda 669-1337, Japan 2DepartmentofPhysicsandAstronomy,UniversityofExeter,StokerRoadEX44QLExeter,UnitedKingdom 3National InstituteforMaterialsScience (NIMS), Namiki 1-1, Tsukuba 305-0044, Japan and 4Department ofAppliedPhysics, TheHongKongPolytechnicUniversity, HungHum, Hong Kong Wepresent atheoryofsurfaceplasmawaves(SPWs)inmetalswitharbitraryelectroniccollisionrateτ 1. 7 − 1 Weshowthatthereexistsauniversalintrinsicamplificationchannelforthesewaves,subsequenttotheinterplay 0 betweenballisticmotionsandthemetalsurface. Weevaluate thecorresponding amplificationrateγ0, which 2 turnsouttobeasizablefractionoftheSPWfrequencyωs.Wealsofindthatthevalueofωsdependsonsurface scatteringproperties,incontrastwiththeconventionaltheory. n a J I. INTRODUCTION interplaybetweenballisticelectronicmotionsandboundaries 4 resultsinanelectricalcurrentthatallowsanetamountofen- ] Collective electronicoscillations on the surface of metals, ergytobedrawnfromtheelectronsbySPWs,whichwouldin- ll dubbed surface plasma waves (SPWs) [1–3], have emerged evitablyself-amplifyinthecollision-lesslimit,therebydesta- a asapivotalplayerinnanoscopicmanipulationoflight[4–8]. bilizing the system. The hereby predicted self-amplification h Thefunctionalityofmanyprototypicalnanophotonicdevices of collision-less SPWs is analogous to the Landau damping - s critically relies on the distance SPWs can travel before they in collisionless bulk plasma waves [28–31]. While the latter e aredampedoutduetoenergylossesviaseveralchannels[7– is caused by slowly moving electrons that strip energy from m 12]. SPWs canlose energydue to Joule heat, inter-bandab- the wave, the former is attributable to ballistically moving t. sorption,radiationemissionandindividualelectronicmotions electronsimpartingenergyto the wave when boundariesare a (Landaudamping). Most of the losses can be efficiently but present. m nottotallyreducedunderappropriatecircumstances. Ampli- Inthenextsection,wedescribethesystemunderconsider- - fiers have been contrived to compensate for the losses [13– ation,stateourmainresultsandconceivetheirpossibleexper- d 21],whichareallextrinsicandrequireexternalagentssuchas imentalsignature.InSec.III,theSPWtheoryispresented.In n o adipolargainmediumtosupplytheenergy. Sec.IV,wedetailthemethodsusedinthetheory.Wediscuss c ThestandardtheoryofSPWswasformulatedshortlyafter the results and conclude the paper in Sec. V. Some calcula- [ theirdiscoveryin1957[1]andhasbeencomprehensivelydis- tionsandargumentsnotcoveredinSec.IIIandIVaregiven 1 coursed in many textbooks [3, 5, 6, 22, 23]. In this theory, inappendicesAandB. v theelectricalpropertiesofmetalsareeffectedbyafrequency- 0 dependent dielectric function ǫ. To analytically treat ǫ, the 6 simple Drude model or the slightly more involved hydrody- II. RESULTS 0 namicmodel[22,24,25,27]isinvariablyinvoked.Foreither 1 model to be valid, electronic collisions must be sufficiently 0 frequentsothattheelectronicmeanfreepath,l =v τ,where System - We consider a prototypical system, namely, a . 0 F 1 v is the Fermi velocity and τ the relaxation time, is much semi-infinitemetal(SIM)occupyingthehalfspacez 0and F ≥ 0 shorter than the SPW wavelength or the typical length scale interfacing with the vacuum at a geometrical surface z = 0. 7 of the system. The general case with arbitrary τ, especially ASIMisnotsheerlyofacademicinterest: itcanberegarded 1 the collision-less limit, where τ , defies these models asonehalfofathickmetalfilm. Inthespiritoftheso-called : v andhasyettobeentertained.Othe→rm∞odelsbasedonabinitio jelliummodel,[32–34]themetalistreatedasafreeelectron Xi quantummechanicalcomputations[22,26]arehelpfulinun- gasembeddedinastaticbackgroundofuniformlydistributed derstandingthecomplexityofrealmaterialsbutfallsshortin positivecharges. Onthewholethesystemisneutral. Theki- ar providinganintuitiveandsystematicpictureofSPWsunder- netic energy of the electrons is written ε(v) = 1mv2, where 2 pinnedbyelectronsexperiencinglessfrequentcollisions. m and v denote the mass and velocity of the electrons, re- InthepresentworkweemployBoltzmann’stransportequa- spectively. Inter-bandtransitionsareaccordinglyignoredbut tiontoderiveatheoryofSPWsinmetalswitharbitraryτ. Our their effects will be discussed in Sec. V. The surface is of a analysisrevealsauniversalintrinsicamplificationchannelfor hard-walltypeandpreventselectronsfromleakingoutofthe SPWs. Theexistenceofthischanneldoesnotdependonτbut system. Throughoutwe reserver = (x,y) for planarcoordi- iswarrantedbyageneralprinciple. Weshowthattheunique nates while x = (r,z) denotes the complete position vector. Weneglectretardationeffectsintotal. Key results - We find that the SPW frequency ω and its s amplificationrateγcanbewritteninthefollowingform, [email protected] ∗ †[email protected] 1 ‡[email protected] ωs =ωs0β(p), γ=γ0− τ, γ0 >0. (1) 2 Here ω = ω /√2 is the SPW frequency in the hydro- own electric field E(x,t). Here t denotes the time. Thanks s0 p dynamic/Drude model, with ω being the characteristic fre- to the linearity and symmetries of the system, we can write p quencyofthemetal,p [0,1]isaparameterthataccountsfor ρ(x,t)=Re ρ(z)ei(kx ωt) andE(x,t)=Re E(z)ei(kx ωt) with- ∈ − − surfacescatteringeffects,β(p) 1weaklydependson pand out loss of ghenerality. Hiere Re/Im takeshthe real/imaiginary ∼ γ0 includesLandaudampingandisexactlyindependentofτ. partofaquantity,k 0denotesthe wavenumberwhile ωis Both ωs and γ0 are determined by the secular equation (18) the SPW eigen-frequ≥ency to be determined by the equation giveninSec.IV.Analyticalexpressionscanbefoundforβ(p) of motion. In general, ω can be complex. Instead of ρ(z), it and γ0 under certain approximations. Exact ωs and γ0 have provesmoreconvenienttoworkdirectlywith beencomputednumericallyandaredisplayedinFig.1,where weEoqbusaetrivoent(h1a)tfγu0rn/ωishs0es≈aαn(i1nt+rinps),icwaimthpαlifi≈ca0t.i1o.nchannelγ0 ρq =Z0∞dzρ(z) cos(qz). for SPWs. It also shows that, in contrast to the hydrody- In order for the jellium model to be valid, we must impose namic/Drude model, ω relies on surface properties via the s that ρ = 0, where q is a cut-off. q 1 roughly gives the parameter p. In particular, as seen in Fig. 1 (a), ωs is sig- q>qc c −c microscopiclatticeconstantofthemetal.Wemaychooseq nificantly – more than 10% – larger than ω , the value ex- c pectedfromthehydrodynamic/Drudemodel.s0Inotherwords, kF = mvF/~ ∼ n01/3 ∼ ks = ωs0/vF, where kF is the Ferm∼i β(p)isunityinthesemodelswhereasitcouldreachupto1.2 wavelengthandn0 themeanelectrondensity. Theresultsdo in our theory. Such a great contrast would be ideal for ex- notdependontheexactvalueofqcaslongasitissufficiently perimentallyverifyingthetheory. Unfortunately,incommon large.Bydefinition, materialssuchasnoblemetals,duetopronouncedinter-band 2 transitions,thereisnosimplerelationbetweenω andω . ρ(z)= ∞dqρ cos(qz). s p πZ q Experimentalsignature - We conceiveone possible reper- 0 cussion of the self-amplification channel, noting that the net Intermsofρ ,wecaneasilyfind,bysolvingLaplace’sequa- q amplification rate γ = γ0 1/τ and its temperature depen- tion,that − dence can be directlymeasured in variousways, e.g. by ex- 4kρ aminingtheSPWpropagationdistanceorthewidthoftheen- E (z)= i ∞dq q 2cos(qz) e kz , (2) ergyloss peaksin electronenergyloss spectroscopy(EELS) x − Z0 k2+q2 h − − i orthequalityofthereflectancedipsintheKretschmann-Otto 4kρ E (z)= ∞dq q 2(q/k)sin(qz) e kz . (3) configuration. In this paper, our calculation of γ0 is done at z Z k2+q2 − − 0 h i zero temperature. However, arguably γ could bear a dif- 0 AsnapshotofE(x,t)foratypicalSPWisdisplayedinFig.3 ferent – probably weak – temperature dependence than 1/τ. In sufficiently pure samples, in which the residual scattering (a)andthemagnitudeofE(z)isplottedin (b). Of course,in makingtheseplotsρ hasbeendeterminedbytheequationof is small enough, there might exist a critical temperature T , q ∗ motiontobesetupshortly. above which γ < 0 while below it γ > 0, i.e. T marks a ∗ Current densities - Under E(x,t), an electrical current, of transition of the system from the Fermi sea to another state, seeSec.Vforfurtherdiscussions.InFig.2,wesketchthesit- densityJ(x,t) = Re J(z)ei(kx−ωt) , willflow. Itcan becalcu- uationdescribedhere,assuming1/τisdominatedbyphonon latedusingBoltzmanhn’stransporitequation,seeSec.IV.Asa andimpurityscattering. crucialobservation,wefindthatJ(x,t)canalwaysbewritten To experimentally investigate the self-amplification chan- in two disparate contributions, J(z) = JD(z) + JB(z). What nelusingthemostexperimentedmaterialslikegoldandsilver, setsthemapartistheirdistinctrelationstoE(z),asillustrated wemustclarifytheeffectsofinter-bandtransitions,whichnot in Fig. 4. It turnsout that JD(z) follows E(z) almost locally, onlybringaboutadditionalenergylossesbutalsostronglyaf- asintheconventionalhydrodynamic/Drudemodel–whichis fectthevalueofω .WediscusstheseeffectsfurtherinSec.V. applicableonlywhenelectronsexecutediffusivemotions,i.e. s τisverysmall. Explicitly,J (z)canbewritten D i ω2 III. THEORY J (z)= pE(z)+J(z), (4) D ω¯ 4π ′ Inthissection, wegiveasystematicexpositionofthethe- where ω¯ = ω+i/τ, ω = 4πn e2/m is the characteristic p 0 orythatsupportsequation(1). Sometechnicalaspectsareleft plasmafrequencyofthemetapland to Sec. IV. The overallobjectiveis to set up the equationof motion for the charge density, analyze its internal structure J(z)= ∞dq8ρqF(k,q;ω¯)cos(qz), (5) andsolveit. We firstdiscusstherelationbetweenthecharge ′ Z k2+q2 0 density andthe electric field, andthen prescribesthe current densityandshowthatitpossessesaspectacularstructure.We signifies a non-local contribution that generates dispersive proceedthencetotheequationofmotion,extractingitssolu- plasmawaves,with tionsandunveilingthepropertiesofSPWs. Electrostatics-Weaimforestablishingtheequationofmo- F(k,q;ω¯)= m 3 d3v( e2f )v ∞ kvx+qvz l. tion for the charge density ρ(x,t) under the influence of its (cid:18)2π~(cid:19) Z − 0′ Xl=2 ω¯ ! 3 (a) (b) 0.95 0.2 Numerical (exact) 0.18 0.9 Fitting by 0.09+0.1 p 0.16 0.85 ωp ωp / / 0.14 ωs γ0 0.8 0.12 Numerical (exact) 0.75 0.1 Equation (1) 0.7 0.08 0 0.2 0.4 0.6 0.8 1.0 0 0.2 0.4 0.6 0.8 1.0 p p FIG.1. NumericallycalculatedSPWfrequencyω anditsintrinsicamplificationrateγ asafunctionoftheFuchsparameter p,byequation s 0 (18).k/k =0.1andq /k =1.5.Theapproximateanalyticalsolutionsarealsoshown. s c s whereJ (z)iscontributedbyelectronsdirectlyemerging B,emg 1/τ at the surface while J (z) by reflected electrons. Here p B,ref is the Fuchsparameter, whichgaugesthe probabilitythat an electron impinging on the surface gets reflected back. Both contributionsarisefromelectronsmovingawayfromthesur- face,v 0. Infact, z γ ≥ 0 m 3 JB,emg/ref(z)=(cid:18)2π~(cid:19) Z d3vΘ(vz)eveiωv˜zzgB,emg/ref(v). (6) whereω˜ =ω¯ kv and x − 4ρ L (v ,v ,k,q;ω¯) residual scattering rate g (v)= ef ∞dq q emg/ref x z ,(7) B,emg/ref − 0′Z k2+q2 0 T* with Temperature k(iv v ) 2(q2v2 ω˜kv ) L (v ,v ,k,q;ω¯)= x− z z ± x . emg/ref x z kv iω˜ − ω˜2 q2v2 z± − z FIG.2. Possibleinstabilitycausedbytheself-amplificationchannel uponcoolingdownthesystem. AtthecriticaltemperatureT∗,γ0 = For p ≈ 1 andsmallkvF/ω¯, it iseasy to show that JB,x(z) ≈ 1/τ. Herethesketchassumesthat1/τisdominatedbyphononand 0, implying that JB(z) points normal to the surface. This is impurityscattering. illustratedinFig.4(b). Positiveness of γ - Equation (6) strongly constrains the 0 value of ω¯. Actually, it requiresthat Im(ω¯) 0; otherwise, Here f0(ε) denotes the Fermi-Dirac distribution and f0′ its theintegralinthisequationwoulddivergefor≥anyz,because dpearpievra.tivOen,lywhteicrmhsiswtiothboedtdakleinnatthezesreoriteesmcpoenrtartiubruetei.nOthnies eiωv˜zz divergesexponentially for vz → 0. In Sec. IV and Ap- pendixA,thispointisfurtherelaborated.Since can show that J (0) 0. The factor headingE(z) in Eq. (4) z′ ≡ is recognized as the Drude conductivity. As this conductiv- 1 γ=Im(ω)=γ , γ =Im(ω¯), (8) ityisnearlyimaginary,JD(x,t)pointsalmostperpendicularto 0− τ 0 E(x,t),asseeninFig.4(a). In contrast, J (z) – the ’surface-ballistic’ current– has no the fact that Im(ω¯) 0 suggests an intrinsic amplification simple relationBwith E(z) and reflects genuine surface ef- channelγ0 competing≥withthelosschannelτ−1. Inwhatfol- fects: JB(z) would totally disappear if the surface were ab- lows, we show thatω¯ is independentof τ−1 by meansof the sent. While J (z) is essentially a bulk property, J (z) orig- equationofmotion. D B inates from the interplay between ballistic motions and the Equation of motion - The equation of motion can be ob- surface.Itconsistsoftwoparts, tainedbyrelatingρ(z)andJ(z)viatheequationofcontinuity, J (z)=J (z)+pJ (z), iω¯ρ(z)+ J(z)= J (0)δ(z), (9) B B,emg B,ref z − ∇· − 4 (a) x (in units of k-1) (b) s 20 30 40 50 60 70 40 14 0.8 12 Metal 30 0.4 -1k)s 10 10 of 8 0 s 20 nit 6 u -0.4 (in 4 5 z 10 -0.8 2 EEzx 0 0 0 0 0.3 0.6 0.9 -6 0 6 Vacuum ρ(z) FIG.3. SketchofSPWssupportedonthesurfaceofasemi-infinitemetal. k/k =0.1and p=1. (a)Chargedensitymap(color)andelectric s fieldmap(arrows). (b)Plotsofρ(z)andE (z),whereµ = x,z. ρ(z)iscalculatedbyEq.(16)whileE (z)by(2)and(3). Bothρ(z)andE (z) µ µ µ havebeennormalized. (a) (b) 14 Diffusive current density 14 Surface - ballistic current density -1k )s 12 -1k )s 12 of 10 of 10 s 8 s 8 nit 6 nit 6 n u 4 n u 4 z (i 02 z (i 02 20 30 40 50 60 70 20 30 40 50 60 70 x (in units of ks-1 ) x (in units of ks-1 ) FIG.4.Snapshotsof(a)diffusivecurrentdensityJ (x,t)and(b)surface-ballisticcurrentdensityJ (x,t)inSPWs.Thesetwocurrentdensities D B arenotdiscriminatedbythevalueofτ,butbytheirdependenceonthesurface.J signifiesgenuinesurfaceeffectsandwouldtotallydisappear B withoutthesurface,whileJ isabulkproperty.k/k =0.1andp=1.Im(ω¯)hasbeenneglectedtoemphasizethedifferences. D s which is derived and discussed in the next section. Here from J (z) with elements (q,q;ω¯) given in Sec. IV. B ′ = (ik,∂ ,∂ ). Substituting J(z) into this equation and In Ap∇pe·ndix B, we show thatM(ω¯) M Z, where M y z 0 0 F∇ouriertransformingit,wearriveat kv /ω isaconstantandZistheMunitym≈atrixwillallelement∼s F p being one. It makes only a minor correction to the diagonal ∞dq (q,q;ω¯) ω¯2δ(q q) ρ =S(ω¯), (10) matrixwithelementsΩ(k,q;ω¯)δ(q q′). ′ ′ ′ q − Z H − − ′ Theequationofmotion(10) cannowbere-castin a com- 0 h i pactmatrixform, whereS(ω¯)=iω¯J (0)countsasasourcetermand z (ω¯) ω¯2I ρ=S(ω¯)E, (12) (q,q;ω¯)=Ω2(k,q;ω¯)δ(q q)+ (q,q;ω¯) H − H ′ − ′ M ′ whereIistheidenhtitymatrix,ρi isacolumnvectorcollecting aretheelementsofamatrixdenotedby (ω¯),with all ρ and E is a column vectorwith all elements beingone. H q Forlateruse,thesourcetermcanberewritten 4πω¯ k F(k,q;ω¯) Ω2(k,q;ω¯)=ω2 + · , k=(k,q). (11) p k·k S(ω¯)=Gt(ω¯)ρ=Z ∞dqGqρq, Gq = 4Gk2(k+,qq;2ω¯), (13) See that Ω(k,q;ω¯) is an even function of ω¯ and depends on 0 thelengthbutnotthedirectionofk. Thematrix (ω¯)arises where (ω¯) is another column vector, t takes the transpose, M G 5 and which determinesω¯ and hence the SPW eigen-frequencyω. Uponomitting (ω¯)from (ω¯),thisequationbecomes M H G(k,q;ω¯)=G (k)+G (k,q;ω¯), D B 4G(k,q;ω¯) 1 1= ∞dq . (18) with Z k2+q2 Ω2(k,q;ω¯) ω¯2 0 − ω2 Let us write the solution as ω¯ = ω + iγ and hence ω = p s 0 GD(k)=k 4π ωs +iγ with γ = γ0 1/τ. See thatthe solutions ωs +iγ0 − ± occur together, in accord with the fact that ρ(x,t) is a real- arisingfromJ (0)and valuedfield. InEq.(18),Ω(k,q;ω¯)isgenerallycomplexand D,z γ automaticallyincludesLandaudamping. 0 GB(k,q;ω¯)=GB,emg(k,q;ω¯)+pGB,ref(k,q;ω¯) Hydrodynamic/Drudelimits - The hydrodynamicmodelis revisited if we replace in Eq. (18) G(k,q;ω¯) and Ω(k,q;ω¯) fromJB,z(0). Here withGD(k) and ωb(k) given by Eq. (15), respectively. If we furtherdisregardthedispersioninω (k),theDrudemodelis b m 3 GB,emg/ref(k,q;ω¯)=iω¯ (cid:18)2π~(cid:19) tIhnebnortehcmovoedreelds,,iω¯niwshriecahl.caseweimmediatelyfindωs = ωs0. d3vΘ(v )( e2f )v L (v ,v ,k,q;ω¯). (14) Approximate solutions - We can solve Eq. (18) approxi- ×Z z − 0′ z emg/ref x z mately. Tothelowestorderinγ /ω ,whichisassumedtobe 0 s small,wemaydetermineω byapproximatingtherealpartof s Equations(10)-(14)areexactanddonotexplicitlyinvolveτ, (18)asfollows therebyconcludingthatω¯ isindependentofτ. Bulkmodesandlocalizedmodes-Withoutthe surfacewe 4Re G(k,q;ω ) 1 have (q,q;ω¯) 0andS(ω¯)=0. Equation(10)reducesto 1 ∞dq s . (19) Ω2(k,Mq;ω¯) ′ω¯2 =≡ 0 forany ρ , 0, which all are extended ≈Z0 (cid:2)k2+q2 (cid:3) Ω2(k,q;ωs)−ω2s q − modes and therefore describe cosine bulk plasma waves. If Substitutingtheso-obtainedω intheimaginarypartof(18), s wediscretizeqinstepdq,intotalthereare Nc = qc/dqsuch wefind modes. Solvingthisequation,wefindthedispersionrelation iωimnbtae(kgg)rinafnaordrytihpneatsrhetedmuinoetdeteogsr.aaGlpieonnlveeo,rlalvolelcyda,tiΩend(Fka,(tqkω¯;,qω¯=;)ω¯ck)ov,uxgl+divpqinovgsz,sreoissfsethatnoe ωγ0s ≈−12 ∞dqR04∞Rdeq[G4(kI,mq;ω[kG2s)+](kq,2q;ωs)]1Ω2(k,q;1ωs)−ω2sω2s , (20) thecelebratedLandaudamping. Neglectingthis,oneobtains R0 k2+q2 Ω2(k,q;ωs)−ω2s Ω2(k,q;ωs)−ω2s forsmall k which can be brought into a rather simple form if we take | | Ω(k,q;ω ) ω2 andω ω . Namely, 3k kv2 s ≈ p s ≈ s0 ω2(k) ω2 1+ · F ω2, (15) whichwaswellbknow≈nfpromthe5hydωr2podyn≈amicp/Drudemodel. ωγ0s ≈−12R0∞∞ddqq RIme[[kGG2((+kkq,,qq2;;ωωss))]] = 12IRme(cid:2)JJzz((00))(cid:3). (21) 0 k2+q2 In the presence of the surface, equation (10) admits not R (cid:2) (cid:3) only extended modes, for which S(ω¯) = 0, but also local- In this approximation, Ω(k,q;ω ) could have an imaginary s izedmodes,forwhichS(ω¯) , 0. Thenumberoftotalmodes partonlyifq > ω /v . Landaudampingwouldbeexcluded c s F cannot change and is still Nc. The extended modes again fromγ0otherwise. represent bulk waves and for them equation (12) reduces to Toproceed,weneedtoevaluateG (k,q;ω¯). To thelinear B (ω¯) ω¯2I ρ = 0. As shown in Appendix B, because of orderink,wefind H − thheconstrainitthatS(ω¯) = 0,thereareintotalatmostN 1 c− m 3 solutionstothisequation. Thenumberoftotalbulkmodesis G (k,q;ω¯) iω¯e2 thenreducedbyonefromthatwithoutthesurface. Theirdis- B ≈ (cid:18)2π~(cid:19) bpyerTωshiboe(nkmr)e.islasitniognbiuslnkemgloigdiebhlyasaffbeecetnedcobnyveMrte(ωd¯)inatnodasltoilclagliivzeedn ×Z d3vΘ(vz)f0′vz2qω¯22v2z−(1q2+v2zp) + kvz(i1ω¯−p). (22) modesatisfyingS(ω¯) , 0andrepresentingSPWs,forwhich Carryingouttheintegralgives equation(12)yields ω2 k(p 1) 3(1+p)v G (k,q;ω¯) p − i Fq2 . (23) ρ=S(ω¯) (ω¯) ω¯2I −1E. (16) B ≈ 4π " 2 − 4 ω¯ # H − h i Thus,Re G(k,q;ω ) k(ω2/4π)(1+ p)/2. Insertingthisin PluggingthisinEq.(13),weobtain s ≈ p Eq. (19),(cid:2)we get ω(cid:3) ω 2 p+1. We see that ω gen- 1= t(ω¯) (ω¯) ω¯2I −1E, (17) erally depends on susrf≈ace ss0caqtter−ing2effects, in contrasst with G H − h i 6 whatis expectedof the hydrodynamic/Drudemodel. Analo- mann’sequationcanbe uniquelydeterminedupto some pa- gously, using Eq. (21) and q v ω , we find γ αω , rameters,whichsummarizetheeffectsof–whilewithoutac- c F s0 0 s0 ∼ ≈ withα = 3/2π. Landaudampinghasbeenexcludedhere,as tually knowing – φ (x). With translational symmetry along s theapproximationonlytakestherealpartofΩ(k,q;ω ). the surface, only one such parameter, namely, the so-called s Numerical solutions - Equation (18) can also be exactly Fuchs parameter p, is needed for the simple specular scat- solved numerically. A comparisonwith the approximateso- tering picture. Physically, p measures the probability that lutionis displayedin Fig. 1. Theagreementin the matterof an electronimpinginguponthe surfaceis bouncedback. As ω is satisfactory, while that for γ is not. The discrepancy usual,wewrite f(x,v,t)= f (ε(v))+g(x,v,t),whereg(x,v,t) s 0 0 mightbe because the approximatesolution excludesLandau isthenon-equilibriumpartofthedistributionfunctionin the damping. It should be emphasized that, our numericalsolu- presence of E(x,t). The current density can then be calcu- tionsdonotdependonthechoiceofq ,providedthelatteris lated by J(x,t) = (m/2π~)2 d3v ev g(x,v,t), where e de- c bigenough,i.e. qc ks. notes the electron charge. IRt is worth pointing out that, as ≥ g(x,v,t) is a distribution only for the bulk, the charge density is not given by ρ˜(x,t) = (m/2π~)2 d3v e g(x,v,t), i.e. IV. METHODS ρ(x,t) , ρ˜(x,t). Actually,ItiseasytoseRethatρ˜(x,t)satisfies (∂ +1/τ)ρ˜(x,t)+∂ J(x,t) = 0ratherthanEq.(24). Obvi- t x · Thissectiondiscussesfurthersometechnicalaspectsofthe ously,whatismissingfromρ˜(x,t)isexactlythechargesjust theory. We set outwith a discussion of the equationof con- localizedonthesurface. tinuity in the presence of surfaces. Then we describe Boltz- Electronic distribution function - Let us write g(x,v,t) = mann’sequationandsolveittoobtaintheelectronicdistribu- Re g(v,z)ei(kx ωt) . Intheregimeoflinearresponses, Boltz- − tionfunctions.Thencewederivethedynamicequationforthe manhn’sequationreiads chargedensity. ∂g(v,z) v E(z) tinEuiqtuy,ati∂otn+o1f/cτonρti(nxu,itt)y+-W∂xesjt(axr,ttw)i=th0th,ewehqiucahtiroenlaotefscothne- ∂z +λ−1g(v,z)+ef0′(ε) ·vz =0, (25) · chargedensityρ(x,t)andthecurrentdensityj(x,t)inageneric (cid:0) (cid:1) where λ = iv /ω˜. In this equation, the velocityv is more of way. Herethedampingterm ρ(x,t)/τisincludedtoaccount z − a parameter than an argument and can be used to tag elec- for electronic collisions that would drive the system toward tronbeams. Itis straightforwardto solve the equationunder equilibrium. In the jellium model, ρ(x,t) appears when the appropriateboundaryconditions(AppendixA).Naturally,we electrondensityisdisturbedawayfromitsequilibriumvalue have n . Thesurfacepreventsanyelectronsfromleakingoutofthe 0 metal. Explicitly,wewritej(x,t) = Θ(z)J(x,t),whereΘ(z)is g(v,z)=g (v,z)+g (v,z), bulk surface theHeavisidestepfunction. Insodoing,wehavetreatedthe surfaceasahardwallandconsideredthefactthatJ(x,t)may where the bulk term would exist even in the absence of sur- notvanishevenintheimmediateneighborhoodofthesurface faceswhereasthesurfacetermwouldnot. Usingtheexpres- – which is obviouslythe case with the hydrodynamic/Drude sionsforE(z),weobtain model.Theequationofcontinuitybecomes 4ρ kv +qv g (v,z)= ef ∞dq q x z eiqz, (26) ∂t+ 1τ!ρ(x,t)+∂x·J(x,t)=−Jz(x0,t)δ(z), (24) bulk − 0′Z−∞ k2+q2ω¯ −(kvx+qvz) whichhasthesameformasonewouldexpectforabulksys- where x0 = (r,0) denotes a point on the surface and δ(z) is tem. Herewehavedefinedρq<0 :=ρ q. theDiracfunctionpeakedonthesurface.Physically,theright As for gsurface(v,z), we find it w−ith a subtle structure: it handsideofEq.(24)meansthat,chargesmustpileuponthe can be written as a sum of two contributions, one of which, surfaceiftheydonotcometoahaltbeforetheyreachit. gD,surface(v,z),hasasingleformforallelectronsirrespective Equation(24)reducestoequation(9)uponusingtheplane oftheirvelocitieswhiletheother,gB,surface(v,z),doesnot.Ex- waveformforρ(x,t)andJ(x,t). plicitly,wehave Boltzmann’sapproach-WeemployBoltzmann’sequation 4ρ k(v iv ) stococpailcculelvateel,Jw(xe,mt)aaysinatrroedspuocensaestuorfEac(xe,pt)o.tenOtinalthφes(xm)icinroto- gD,surface(v,z)=−ef0′Z0∞dq k2+qq2 kvzz−+iω˜x e−kz. (27) theequationtoaccountforsurfacescatteringeffects.Thecor- Wemaycombineg (v)andg (v,z)inasingleterm, bulk D,surface responding surface field E (x) = ∂ φ (x) is peaked on the s x s − surfaceand,complyingwiththehard-wallpicture,mayhave g (v,z)=g (v,z)+g (v,z), D bulk D,surface aninfinitesimalspreadd 0 . Unfortunately,asφ (x)can s + s hardlybeknownandvarie→sfromonesampletoanother,this inordertoseparatethemfrom microscopicmethodisimpracticalandfutile. g (v,z):=g (v,z). Alternatively, surface scattering effects can be dealt with B B,surface using boundaryconditions[32, 34–36]. This is possible be- Insodoing,wehavedecomposed cause E (x) acts only on the surface. In the bulk, the solu- s tions–theelectronicdistributionfunction f(x,v,t)–toBoltz- g(v,z)=g (v,z)+g (v,z) D B 7 ina diffusiveanda ballisticcomponent. Itisunderlinedthat electroniccollisionrate1/τ. Asakeyconsequenceofthethe- g (v,z) arises only when the surfaces are present. For bulk ory, we find thatthere exists a self-amplificationchannelfor B systems without surfaces, it does not exist even if the elec- SPWs,whichwouldcausethelattertospontaneouslyamplify tronicmotionsaretotallyballistic, i.e. τ . Assuch, we at a rate γ if not for electronic collisions. Surprisingly, the 0 → ∞ callit’surface-ballistic’. valueofγ turnsouttobeindependentofτ. Thepresenceof 0 g (v,z) exists only for electrons leaving the surface, i.e. thischannelisguaranteedbythecausalityprinciple. B vz 0. Those electronscoulddirectly emergefromthe sur- Whetherthesystemcouldactuallyamplifyornotdepends ≥ face or be those incident on but subsequently get reflected on the competition between γ and 1/τ. If γ > 1/τ, SPWs 0 0 bythesurface. WhatfundamentallysetsgB,surface apartfrom willamplifyandthesystemwillbecomeunstable.Inourthe- gD(v,z)restwithitssimplezdependence.Actually,wehave ory,thenon-equilibriumdeviationg(v,z)referstotheFermi- Dirac distribution f (ε); as such, the instability is one of the gB(v,z)=Θ(vz)eiωv˜zz gB,emg(v)+ pgB,ref(v) , Fermisea. Needless0 tosay,theinstabilitywillbeterminated h i once the system deviatesfar enoughfromthe Fermisea and wheretheFuchsparameter pgaugesthefractionofreflected settlesinastablestate. Clarifyingthenatureofthedestination electrons and g (v) has been given in Eq. (7). As al- B,emg/ref stateisasubjectofcrucialimportanceforfuturestudy. ready remarked, since gB(v,z) ∝ eiω˜z/vz ∝ e−γ0z/vz, we must Onecentralfeatureofourtheoryistheclassificationofcur- haveγ 0;otherwise,itwoulddivergeeitherwhenz 0 ≥ → ∞ rentdensitiesintoadiffusivecomponentJ (z)andasurface- or for small v . In Appendix A, we argue that this result is D z ballisticcomponentJ (z). Thisclassificationisnotbasedon aconsequenceofthecausalityprinciple,whichstatesthatthe B thevalueofτbutaccordingtowhetherthecomponentobeys numberofreflectedelectronsisdeterminedbythatofincident the (generalized) Ohm’s law or not. Apart from this, these electrons,ratherthanotherwise. componentsarealsodiscriminatedinotherways.Firstly,they Divergence of the current densities - The current density J(z)=(m/2π~)2 d3vevg(v,z)issplitintwoparts:J (z)and arecontrolledbydifferentlengthscales. Asitlargelyfollows D the local electric field E(z), the characteristic length associ- J (z), where J R(z) = (m/2π~)2 d3v ev g (v,z). Using B D/B D/B atedwithJ (z)isk 1. Ontheotherhand,thelengthforJ (z) D − B the expressions for gD/B(v,z), it isR straightforward to obtain isv /γ ,becauseofitssimplez-dependence. Secondly,they F 0 J (z) given by Eq. (4) and J (z) by Eq. (6). To obtain the D B are oriented disparately. J (z) is largely oriented normal to D equation of motion (10) from the equation of continuity (9), E(z)locallywhereasJ (z)normaltothesurfaces–especially B theFouriertransformof J (z)isneeded.Wefind ∇· D/B for pclosetounity. Thirdly,JD(z)isabulkpropertyandex- Z ∞dz cos(qz)∇·JD(z)= ω¯i Ω2(k,q;ω¯)ρq, (28) tirsutserseugrafardcleeessffeocftsthaendsuirtfwacoeu;ldOdnisthapepceoanrtwrairtyh,oJuBt(szu)rfraecfleesc.ts 0 If we replace the vacuum by a dielectric ǫ, γ may be re- 0 withΩ(k,q;ω¯)givenbyEq.(11). Additionally, duced by an order of 1/ǫ due to weakened E (z) and J (0). z z Roughlyspeaking,ρ(z)presentonthemetalsideinducesmir- i ∞dzcos(qz) J (z)= ∞dq (q,q)ρ , (29) ror charges amounting to ρ(z) = ρ( z)(ǫ 1)/(ǫ + 1) on Z0 ∇· B ω¯ Z0 ′M ′ q′ thedielectricside. Ifρ(z) is′highly−loc−alized−abouttheinter- where (q,q′)isamatrixgivenby face, as is with SPWs, ρ(z) and ρ′(z) will combine to give a M netchargeof2ρ(z)/(ǫ +1). Asa result, E and J (0)willbe z z 4ω¯2 m 3 reducedbyafactorof2/(ǫ+1). Thisleadstosmallerγ and (q,q)= 0 M ′ k2+q′2 (cid:18)2π~(cid:19) smallerωs0 = ωp/√ǫ+1. Studyingdielectriceffectsmaybe iω˜v importantinapplications. d3vΘ(v ) e2f z L(v ,v ,k,q;ω¯),(30) ×Z z − 0′ ω˜2 q2v2 x z ′ Anotherproblemthatneedstobeaddressedforexperimen- (cid:16) (cid:17) − z tal studies is concerned with the effects of inter-band tran- where sitions. In the most experimented materials, such as silver L(v ,v ,k,q;ω¯)= L (v ,v ,k,q;ω¯)+ pL (v ,v ,k,q;ω¯). and gold, these transitions are known to have dramatic ef- x z emg x z ref x z fects. They not only open a loss channel due to inter-band InAppendixB,weshowthat (q,q′)generallymakesami- absorption, but also significantly shift the SPW frequency. norcorrection,oftheorderofMkvF/ωs,toΩ2(k,q;ω¯)δ(q q′). Including them in our formalism consists of a simple gen- − Ultimately, the insignificance of this correction may be as- eralization: in addition to J (z) and J (z), the total cur- D B cribed to the factor eiω˜z/vz in gB(v,z), which is extremelyos- rentdensityJ(z)mustnowalso havea componentJint(z) ac- cillatoryoverzwithaperiodlessthan vF/ωs. counting for inter-band transitions. The equation of motion ∼ SubstitutingEqs.28and(29)intheequationofcontinuity, is obtained by substituting J(z) in Eq. (9). One may write weimmediatelyobtaintheequationofmotion(10). J (z)= dz σ (z,z;ω)E (z),whereµ,ν= x,y,zand int,µ ν ′ µν ′ ν ′ theinter-bPandRconductivityσµν caninprinciplebecalculated usingGreenwood-Kuboformula. Inpractice,calculatingσ µν V. DISCUSSIONSANDCONCLUSIONS could be a formidable task even for the imaginablysimplest surfaces. Nevertheless, one may argue that J (z) primarily int Thus, on the basis of Boltzmann’s equation, we have es- affectsthe propertiesofbulkwaves, namely,Ω(k,q;ω¯). The tablishedarigoroustheoryforSPWsinmetalswitharbitrary causalityprincipleshouldstillprotecttheamplificationchan- 8 nel. Asystematicanalysiswillbepresentedelsewhere. z for any v. For electrons movingaway from the sur- → ∞ We remarkthat, γ can also be calculatedby studyingthe face, i.e. v > 0, thisconditionis fulfilledfor anyC(v). For 0 z temporalevolutionoftheelectrostaticpotentialenergyofthe electrons moving toward the surface, i.e. v < 0, we may z system. In particular, equation (21) can be directly derived choose inthisway. Detailedcalculationsalongthislinewillbepub- e∂ f lishedinaseparatepaper. C(v)= v 0 ∞ dz′e−iω˜vzz′E(z′), vz <0, (A2) To conclude,we havepresenteda theoryforSPWs taking mvz ·Z0 intoaccounttheuniqueinterplaybetweenballisticelectronic whichleadsto motions and boundary effects, from which it emerges a uexnpiveecrtseadltshealtf-tahmepsltiufidcyatiwonillcbheaanrnefalrf-oreracthheinseg wpraavcetisc.alIatnids g(v,z)= em∂vvfz0 ·Zz∞dz′eiω˜(zv−zz′)E(z′), vz <0. (A3) fundamentalconsequences,tobeexploredinthefuture. TodetermineC(v)forv >0,theboundaryconditionatz=0 z Acknowledgement – HYD acknowledges the Interna- hasto be used, which dependson surface properties. In this tional Research Fellowship of the Japan Society for the paper, we assume that a fraction p (Fuchs parameter)of the Promotion of Science (JSPS). This work is supported by electrons impinging on the surface are bounced back in the JSPS KAKENHI Grant Nos. 15K13507and 15K21722and absenceofE (z),i.e. z MEXTKAKENHIGrantNo. 25107005. g (v ,v ,v >0),z=0 = pg (v ,v , v ),z=0 x y z x y z − (cid:16) (cid:17) (cid:16) (cid:17) AppendixA:Electronicdistributionfunctions evaluatedatEz(z)=0. Thenweget ThegeneralsolutiontoEq.(25)isgivenby C(v)= p e∂vf0 ∞dz′eiω˜(zv+zz′)E(z′), vz >0. (A4) − mvz ·Z0 e∂ f z g(v,z)=eiωv˜zz C(v) v 0 dz′e−iω˜vzz′E(z′) , (A1) Bydefinition, pvariesfromzerotounity.Thus,weobtain − mvz ·Z0 ! g(v,z)=Θ(v )g (v,z)+Θ( v )g (v,z), z > z < where C(v) is an arbitrary integration constant to be deter- − mined by boundary conditions. We require g(v,z) = 0 at where g>(v,z)= zeiω˜|z−vzz′|+p ∞eiω˜|z+vzz′| ev·E(z′)∂f0 dz′, g<(v,z)= ∞eiω˜|z−vzz′|ev·E(z′)∂f0 dz′. (A5) −"Z0 Z0 # |vz| ∂ε −Zz |vz| ∂ε Utilizingequations(2)and(3)forE(z)andcarryingouttheintegraloverz,wefind ′ ∂f ε(v) e 4kρ g (v,z)= 0 ∞dq q g (v,z;k,q), (A6) >/< − ∂(cid:0)ε (cid:1) vz Z0 k2+q2 >/< (cid:12) (cid:12) (cid:12) (cid:12) where (cid:12) (cid:12) g>(v,z;k,q)= kv+z−iω˜iv/vx e−kz+ 2((cid:18)ω˜q/vvzqk)2+ω˜qvv2xz(cid:19) cos(qz)+i 2(ω˜(cid:16)q/vvx)+2 ω˜qqk2(cid:17) sin(qz) z z z − − 2 ω˜ vx +qv q 2 ω˜ vx qv q and +kiv+xi−ω˜/vvzz − ((cid:18)ω˜/vvzz)2−qz2k(cid:19)+ pkiv−xi−ω˜/vvzz + ((cid:18)ω˜/vvzz)−2−qz2k(cid:19) eiωv˜zz, (A7) g<(v,z;k,q)= ki−vxiω−˜/vvzz e−kz+ 2(cid:18)(qω˜(cid:12)(cid:12)(cid:12)v/zv(cid:12)(cid:12)(cid:12)zqk)2+−ω˜q|2vvzx|(cid:19) cos(qz)−i 2(ω˜(cid:16)q/vvzx)+2−ω˜qqk2(cid:17) sin(qz). (A8) (cid:12) (cid:12) (cid:12) (cid:12) Now we can combine the terms with cos(q(cid:12)z)(cid:12)and those with sin(qz) into a single term called g (v,z), while the rest into bulk g (v,z).TheirexpressionshavebeengiveninSec.IV. surface Ifthesurfaceisnotuniform,weshouldhaveafunction p(r)insteadofaconstant p. Theboundaryconditionshouldthenbe writteng(x ,(v ,v ,v >0),t)= p(r)g(x ,(v ,v , v ),t). Insuchcase,thetranslationalsymmetryalongthesurfaceisgenerally 0 x y z 0 x y z − 9 lost and one cannotwork with a single k-componentany more, i.e. there is scattering effects and differentk-componentsare mixed.Wedonotconsiderthisinthispaper. Causalityprinciple-InapplyingtheboundaryconditionstoobtainC(v),wehaveimplicitlyassumedIm(ω˜) 0;otherwise, ≥ wewouldfindunphysicalsolutionsthatviolatetheprincipleofcausality,whichstatesthatthenumberofout-goingelectronsis determinedbythenumberofin-comingelectrons,nototherwise. Itiseasytoshowthat,hadweassumedIm(ω˜)<0,wewould havefoundtheopposite: thenumberofreflectedelectronswouldbefixedwhilethenumberofincidentelectronswouldgoto infinityas p 0,whichisunphysical. → AppendixB:Thematrix (ω¯) M Inthefirstplace,weshowthat (ω¯)/ω2 ikv /ω¯ +...,wheretheellipsisstandsforhigherordertermsinkv /ω¯. Forthis M p ∝ F F purpose,letusdiscretizeqinstepdq = 2π/d,i.e. q = l(2π/d),wherekd 1andl = 0,1,2,...takesintegervalues. Equation l ≫ (10)canthenbebroughtintothefollowingform Ω2(k,q;ω¯) ω¯2 ρ + (ω¯)ρ =S(ω¯), (B1) l l l,l l − M ′ ′ (cid:16) (cid:17) Xl ′ whereMl,l′(ω¯) = (2π/d)M(ql,ql′;ω¯)andρl = ρql. Writing d3vΘ(vz) = 02πdϕ 0π/2dθsinθ 0∞dv2(v/2)andintegratingover v,wefind R R R R (ω¯) 3 2π π/2 ω¯2cosθ ω¯kv sinθcosθcosϕ kω¯/v Ml,l′ =i dϕ dθsinθ − F F L(v sinθcosϕ,v cosθ,k,q ,ω¯).(B2) ω2 2πkd Z Z (ω¯ kv sinθcosϕ)2 q2v2 cos2θk2+q2 F F l′ p 0 0 − F − l F l′ Tothelowestorderinkv /ω¯,weonlyneedtoretainL(0) intheexpansionLsym = L(m)(kv /ω¯)m. Thus, F ∞m=0 F P 2q2v2 cos2θ L(v sinθcosϕ,v cosθ,k,q,ω¯) F 1+p . (B3) F F ≈ ω¯2 q2v2 cos2θ − F (cid:0) (cid:1) Substitutingthisbackin(B2)andapproximating ω¯2cosθ ω¯kv sinθcosθcosϕ ω¯2cosθ q2 − F , 1, (B4) (ω¯ kv sinθcosϕ)2 q2v2 cos2θ ≈ ω¯2 q2v2 cos2θ q2+k2 ≈ − F − F − F wearriveat (ω¯) 3(1+p) kv 2π π/2 cos3θ 1 Ml,l′ = i F dϕ dθsinθ . (B5) ω2p − πkd ω¯ !Z0 Z0 1−(qlvF/ω¯)2cos2θ1−(ql′vF/ω¯)2cos2θ Obviously,wehave (ω¯)/ω2 kv /ω¯,asstated. M p ∼ F WemayproceedfurtherIfwetake 3 2π π/2 cos3θ 1 3 2π π/2 dϕ dθsinθ dϕ dθsinθ cos3θ 1, (B6) πZ0 Z0 1−(qlvF/ω¯)2cos2θ1−(ql′vF/ω¯)2cos2θ ≈ πZ0 Z0 ∼ from which it follows that (ω¯) M = iω2(1/kd)(kv /ω¯), which is a constant. We then write (ω¯) M Z, where Z = 1 constitutesa unityMmatrixl,l.′ T≈o sim0 plif−ythpe followingF discussions, let us take Ω(k,q;ω¯) = ω . MEquatio≈n(B01) is now l,l p ′ written (ω2 ω¯2)I+M Z ρ=S(ω¯)E. (B7) p− 0 h i For bulk modes, S(ω¯) = 0. Note that Z has only one non-vanishingeigenvalueamounting to its dimension N = q d/2π c c (ω /v )(kd/2π). Thecorrespondingeigenvectorisρ E, whichobviouslydoesnotsatisfyS(ω¯) = 0andisnotabulkmode∼. s F WethereforeconcludethatthereareintotalatmostN∝ 1bulkmodes. Forthesemodes,theeigenvaluesofZareallzeroand c − therefore (ω¯)hasnoimpactonbulkmodes. WecanMalsoshowthat (ω¯)isnegligibleforthelocalizedmode,forwhichS(ω¯),0. Tothisend,wewrite M 1 1 ω2 ω¯2 I+M Z − =U 1 ω2 ω¯2 I+M Z˜ − U, (B8) (cid:20)(cid:16) p− (cid:17) 0 (cid:21) − (cid:20)(cid:16) p− (cid:17) 0 (cid:21) 10 whereU is a similarity transformationthatdiagnolizesZ. We have useda tilde to indicatethe transformedmatrices, e.g. we writeZ˜ = UZU 1. Assaid,Z˜ hasonlyonenon-vanishingelement. Letitbethel -thelement. ThenZ˜ = N δ δ . Seethat M 1/N and −˜ (ω¯) i(ω2/2π)δ δ . Introducing ˜t = tU 1 andE˜ = U0E,wecanrewritethel,ls′eculacrel,lq0ula′,tli0onforthe 0 ∼ Ml,l′ ∼ − p l,l0 l′,l0 G G − localizedmodeas 1 1= ˜t ω2 ω¯2 I+M Z˜ − E˜. (B9) G (cid:20)(cid:16) p− (cid:17) 0 (cid:21) Explicitly, 1 1 1 1 1 1= ˜t E˜ = t E + ˜t E˜ t E t E. Glω2 ω¯2+M Z˜ l Glω2 ω¯2 l Gl0ω2(1+i/2π) ω¯2 l0 −Gl0ω2 ω¯2 l0≈ Glω2 ω¯2 l Xl p− 0 l,l Xl p−  p − p−  Xl p− Theterminthesquarebracketmakesonlyacontributionoftheorderof 1/N andcanbeneglectedforlargeN . c c ∼ [1] R.H.Ritchie,Phys.Rev.106,874(1957). [20] S. Keńa-Cohen, P. N. Stavrinou, D. D. C. Bradley and S. A. [2] R.A.Ferrell,Phys.Rev.111,1214(1958). Maier,NanoLetters13,1323(2013). [3] H. Raether, Surface plasmons on smooth and rough surfaces [21] C. M. Aryal, B. Y.-K.Huand A.-P.Jauho, arXiv:1508.01271 andongratings(SpringerBerlinHeidelberg,1988). (2015). [4] B.RothenhaüslerandK.Wolfgang,Nature332,615(1988). [22] J.M.Pitarke,V.M.Silkin,E.V.ChulkovandP.M.Echenique, [5] A.V.Zayats,I.S.IgorandA.A.Maradudin, Phys.Rep.408, Rep.Prog.Phys.70,1(2007). 131(2005). [23] D. Sarid and W. Cgallener, Modern introduction to surface [6] S. A. Maier Plasmonics: fundamentals and applications plasmons: theory, mathematic modeling and applications (SpringerScience&BusinessMedia,2007). (CambridgeUniversityPress,Cambridge,UK,2010). [7] E.Ozbay,Science311:189(2006). [24] J.Harris,Phys.Rev.B4,1022(1971). [8] M.L.BrongersmaandP.G.Kik,Surfaceplasmonnanophoton- [25] A. L. Fetter, Ann. Physics 81, 367 (1973); Phys. Rev. B 33, ics(Springer,2007). 3717(1986). [9] W.L.Barnes,A.Dereux,andT.W.Ebbesen,Nature424,824 [26] P.J.Feibelman,Prog.Surf.Science.12,287(1982). (2003). [27] O.Schnitzer,V.Giannini,S.A.MaierandR.V.Craster,Proc. [10] W.L.Barnes,J.ofOpticsA8,S87(2006). R.Soc.A472,20160258(2016). [11] T.W.Ebbesen,C.GenetandS.I.Bozhevolnyi,PhysicsToday [28] L.Landau,J.Phys.USSRX,25(1946). 61,44(2008). [29] A. Y. Wong, R. W. Motley and N. D’angelo, Phys. Rev. 133, [12] D. M. Giuliana, Y. Sonnefraud, S. Keńa-Cohen, M. Tame, S. A436(1964). K.O¨zdmir,M.S.KimandS.A.Maier,NanoLetters12,2504 [30] R.BinghamandJ.T.MendoncaandJ.M.Dawson,Phys.Rev. (2012). Lett.78,247(1997). [13] D.J.BergmanandM.I.Stockman,Phys.Rev.Lett.90,027402 [31] S.S.NatuandR.M.Wilson,Phys.Rev.A88,063638(2013). (2003). [32] J. M. Ziman, Electronsand Phonons: the theory of transport [14] J.Seidel,S.FrafstronandL.Eng, Phys.Rev.Lett.94,177401 phenomenainsolids(OxfordUniversityPress,2001). (2005). [33] D. Pines, Elementary excitations in solids (W. A. Benjamin, [15] I.DeLeonandP.Berini,Phys.Rev.B78,161401(R)(2008). NewYork,1963) [16] I.DeLeonandP.Berini,Nat.Photonics4,382(2010). [34] A.A.Abrikosov,Fundamentalsofthetheoryofmetals(Elsiver [17] D. Y. Fedyanin and A. Y. Arsenin, Opt. Express 19, 12524 SciencePublishersB.V.,North-Holland,1988) (2011). [35] G.E.H.ReuterandE.H.Sondheimer, Proc.R.Soc.Lond.A [18] P.BeriniandI.DeLeon,Nat.Photonics6,16(2012). 195,338(1948). [19] D.Y.Fedyanin, A.V.ArseninandA.V.Zayats,NanoLetters [36] M.I.Kaganov,G.Y.LyubarskiyandA.G.Mitina,Phys.Rep. 12,2459(2012). 288,291(1997).