Adiabatic amplification of plasmons and demons in 2D systems Zhiyuan Sun,1 D. N. Basov,1,2 and M. M. Fogler1 1 Department of Physics, University of California San Diego, 9500 Gilman Drive, La Jolla, California 92093 2 Department of Physics, Columbia University, 538 West 120th Street, New York, New York 10027 (Dated: November2,2016) We theoretically investigate charged collective modes in a two-dimensional conductor with hot electrons where the instantaneous mode frequencies gradually increase or decrease with time. We show that the loss compensation or even amplification of the modes may occur. We apply our theory to two types of collective modes in graphene, the plasmons and the energy waves, which can be probed in optical pump- 6 probeexperiments. 1 0 2 Introduction. Plasmons in metals, semiconductors, and (a) (b) 𝑋(𝑡) g other solid-state systems have been a topic of intensive re- u search for over half a century [1]. Plasmonics has found a A 𝜅>𝛾 Time number of technological applications in chemical sensing, e 9 ltiagthiotnmbayniuplutrlaastihoonr,talansderinpfuolrsmesat[i2o–n4]pirsoocensesionfgt.hPehmoteotehxocdis- dutilpm 𝜅=𝛾 (c) 𝑋(𝑡) ] togenerateplasmons. Whenthepulsedexcitationisofhigh A l 𝜅<𝛾 l enough power, it can modify material properties of either a Time h the plasmonic medium or its electromagnetic environment, Time - whichistheprincipleunderlyingtheemergingfieldofactive s e plasmonics [2, 5–7]. For example, photoexcitation-induced FIG. 1. (a) A schematic showing the amplitude of a plasmon as m populationinversionmaypermitplasmonlosscompensation a function of time t for different relations between the mode fre- . oramplification[5,6]. Moreoften,plasmonlifetimeremains quency decay rate κ and the damping rate γ. (b) The canonical t a quite short, e.g., tens of femtoseconds (fs) in noble metals, coordinate X as a function of t in the κ>γ case. The amplitude m which is an obstacle to applications. In experiments using growsasthefrequencydrops. (c)X(t)forthecasewhereamplifica- - ultrafastopticalpulses,theplasmonfrequencychangeswith tion occurs while frequency increases, as in the tunneling process d sketchedinFig.3(a)below. time as the system relaxes back to equilibrium. However, n because of high damping, it has been customary to treat o c plasmonicresponseofthesystemasquasi-stationaryduring [ theplasmonlifetime. tionally,weshowthatthesameconceptappliestotheenergy 2 Recently, graphene has emerged as a new plasmonic waveingraphene[20,21],whichisacollectivemodesimilar v medium distinguished by record-high tunability and con- to acoustic plasmons (or “demons” [1]) in metals and semi- 2 finement [8, 9]. Combating damping remains a chal- conductors [22–24] and also to “cosmic sound” in the early 72 lenge; however, plasmon quality factors as high as Q ∼30 universe[25]. 2 have been demonstrated [4, 10] for graphene encapsulated Qualitativepicture. Tomodelanonequilibriumsystemun- 0 in hexagonal boron nitride. A new scientific frontier in derintensephotoexcitationweassumethatitselectrontem- 1. grapheneplasmonicsisnonlinear[11,12]andnonequilibrium perature T is much larger than the lattice temperature Tl. 0 dynamics probed in ultrafast optical experiments [13, 14]. Suchahot-electronstatetypicallyformsinmetalsandsemi- 6 Plasmon amplification through stimulated emission [15, 16] conductorafewtensoffsafteropticalpumping. Thisrapid 1 has been proposed theoretically and plasmon switching by thermalization(thatis,relaxationoftheelectrondistribution : v opticalpumpinghasbeendemonstratedexperimentally[4]. to the Fermi-Dirac form with the temperature T) is due to Xi These encouraging developments motivate us to study stronginteractionsofelectronswitheachotherandwithop- the regime where the plasmon lifetime is comparable or tical phonons. Subsequently, T gradually decreases toward r a longer than the characteristic relaxation time in a mate- T atamuchslower“cooling”ratemeasuredinpicoseconds l rial. Although this regime may or may not be realizable (ps), predominantly due to emission of acoustic phonons. in graphene, we consider this as a theoretical possibility. If the plasmon dispersion depends on T, plasmons propa- Previously, collective modes in media undergoing adiabatic gating in this transient state would have a slowly changing evolutionhavebeendiscussedintheoreticalastrophysics[17], frequency. We will show that such an adiabatic change of plasma physics [18], and general relativity [19]. In this Let- theplasmonfrequencycouldinduceadiabaticamplification ter, we apply similar ideas to solid-state materials, which oftheplasmonamplitude. are better suited for controlled experiments. Our key find- Adiabaticchangeofparametershasbeenpreviouslycon- ingisthatlosscompensationorevenamplificationcanbea sidered in the context of plasmon-polariton focusing in ta- natural outcome of the transient plasmon dynamics. Addi- pered waveguides [26]. As plasmon approaches the narrow 2 end of the waveguide, its group velocity decreases and its This model is motivated by a popular physical picture (see, electric field increases. In this situation the change of pa- for example, Refs. 35) where the current damping occurs rameters occurs in space. The mechanism we study relies because the “density of photoexcited carriers” decays with instead on having parameters changing in time. To explain the rate 2κ and because of additionally, these carriers ex- ourkeyidealetustreattheplasmonasaharmonicoscillator perience momentum relaxation with the rate Γ, see [36] for withtheequationofmotion furtherdiscussion. Let us focus for now on the case of undoped graphene (cid:161)∂2t+γ(t)∂t+ω2(t)(cid:162)X =0 (1) wheretheDrudeweightD(t)isproportionaltotheelectron foritscanonicalcoordinate X(t)(e.g.,chargedensity). Here temperature [30] D(t)=2ln2ħe22T(t) and where the under- γ(t)isthedampingrateandω(t)istheinstantaneousmode damped plasmons exist at frequencies τ−1(cid:191)ω(cid:191)T/ħ. The ee frequency. Suppose ω(t) changes monotonically with the lowerlimitissetbytheelectron-electronscatteringrateτ−1; ee decay rate κ≡−∂ lnω, which is slow enough, ω(cid:192)κ, then the upper limit is imposed by the Landau damping due to t the Wentzel-Kramers-Brillouin (WKB) approximation to the the interband transitions, see Fig. 2. In particular, the di- solutionofEq.(1)isvalid: mensionless Landau damping rate of the thermal plasmons is given by [34] Γ/ω=(π/16ln2)(ħω/T)2, which is small if X(t)=A(t)e−iS(t), (2) ħω(cid:191)T. Note also that the assumption of scalar D can t be justified if κ and γ are much smaller than the electron- 1 1(cid:90) A(t)=(cid:112)ω(t)exp−2 γ(t0)dt0, (3) etrliebcutrtoionnre(ilnaxthateioanbsreanteceτ−eoe1fsaoptrhoabtea)nisismoatrionptaicineeldec.trondis- 0 It is straightforward to show that the equation of motion t (cid:90) S(t)= ω(t )dt . (4) for the plasmon has the same form as Eq. (1) with X equal 0 0 to ρ , the Fourier harmonic of the charge density, and with 0 thedqissipationrateequalto If both γ and κ are constant, the time-dependent plasmon amplitudehastheform γ=2κ+Γ. (8) A(t)=e21(κ−γ)t. (5) Unfortunately, the condition κ>γ seems impossible to sat- isfy since Γ>0 and κ>0. In other words, the amplification Clearly,thefrequencydecayrateκcompeteswiththedamp- cannotoccurduetotheplasmondampingratebeinglarger ingrateγ. Iftheconditionκ>γismet,thentheoscillation than the frequency decay rate, see also the supplementary amplitudeincreaseswithtime,asshowninFigs.1(a)and(b). material. Although the adiabatic principle appears simple and Suppose,however,thattheDrudeweightisgrowing,κ<0. straightforward, its application to actual solid-state systems In this case the criterion for amplification κ<−Γ can be may require sorting out some important details. In the re- met if the growth rate is fast enough, see Fig. 2(c). Under mainderofthisLetterwedosoontheexamplesoftwotypes of collective modes: the plasmons and the energy waves in 𝜔 graphene. Plasmons in two-dimensional (2D) materials. 2D materi- 𝑇/ℏ als are very promising for active plasmonics because they Plasmon 𝜔 ∼ 𝑞 are not affected by a finite penetration length of optical beams and are much more tunable than bulk metals. It is well known [27–29] that such plasmons have a char- acteristic square-root dispersion with momentum (Fig. 2), 𝜏𝑒−𝑒1 𝑣𝐹 𝑞 Energy wave (cid:113) 2 ω = 2Dq,where(cid:178)isthepermittivityoftheenvironment q (cid:178) and D is the Drude weight (see below). Our goal is to show 𝜏−1 𝜔 ∼ 𝑞 that the time-dependence of D may give rise to adiabatic 𝑝ℎ amplificationofplasmons. 𝑙−1 𝑙−1 𝑞 If the system has the spatial translational symmetry, dif- 𝑝ℎ 𝑒𝑒 ferent momenta are decoupled. For a given q, in the linear- FIG. 2. Dispersion of the plasmon [30] and the energy wave [20, response regime, the plasmon dynamics is determined by 21, 31] in a weakly doped graphene with hot electrons (schemati- theelectricalconductivityoperatorwiththekernelσq(t,t0). cally). The plasmon (energy wave) exists at ω above (below) τ−ee1; Considerthefollowingmodelfortheconductivitykernel: otherwise, it is overdamped, as indicated by the fainting ends of σ (t,t )= 1D(t)e−Γ(t−t0)θ(t−t ), (6) tahnedcau∼rv4es[.21F,o3r2T–3(cid:192)4].µThanedplαas(cid:191)mo1n,iτs−eae1ls=oaoαve2rTdawmipthedαa=tωe2>∼/(cid:178)Tħ/vħF, q 0 π 0 whiletheenergywaveisdampedbyelectron-phononanddisorder D(t)=D(0)e−2(cid:82)0tκ(t(cid:48))dt(cid:48). (7) scatteringcharacterizedbytherateτ−ph1. 3 what conditions can this scenario be realized? One possi- tion is the energy wave in graphene. This mode is pre- bility is to leverage the dependence of the Drude weight on dicted [20, 21, 31] to exist in the hydrodynamic regime of the carrier density or effective mass, which is another com- frequencies that are lower than the electron-electron colli- mon attribute of ultrafast pump-probe experiments [2, 7]. sionrateτ−1,seeFig.2. Inthisregime,onlycollectivevari- ee We speculate that the plasmon amplification may be pos- ables immune to interparticle collisions, i.e, the zero modes sible by exploiting tunneling in a vertical semiconduc- of the collision integral are important: the local tempera- tor/insulator/graphene heterostructure, see Fig. 3(a). The ture T(r), chemical potential µ(r), and drift velocity u(r). semiconductor could be, e.g., a transition-metal dichalco- Theirdynamicsisdescribedbyasetofhydrodynamicequa- genide and the inslulator could be hexagonal boron nitride tions[21,33,36]. Theenergywaveisthepropagatinglongi- (hBN), as in recent experiments [37]. With a suitable bias tudinal mode resulting from this set of equations. Consider voltage applied, the initial state with a lower electrochemi- a weakly doped graphene, µ(cid:191)T. The dispersion relation calpotentialingraphenecanbemaintainedastheinsulator oftheenergywave,neglectingdissipation,is bandgap would prevent electron tunneling in any direction. (cid:115) However, once they are heated to energies close or above 1 4πe2 n2 ω =v q2+ q, (10) the insulator’s band edge, the electrons in the semiconduc- q F 2 3 (cid:178) n E tor layer would tunnel to graphene. (This is similar to a hot-electron doping effect [38] whereas in [37] the tunneling where n is the average electron density and nE ≡〈ε〉 is the average kinetic energy density (relative to the zero-doping, was in the opposite direction.) For tunneling to be rapid the insulator must be thin, which implies that the charges zero-temperature state). The latter behaves as [36] nE ∝T3 andcurrentinthetwolayerswouldalsobecoupledelectro- in the regime we consider, T (cid:192)|µ|. For q(cid:192) e2 n2 ≡q , the (cid:178) nE c magnetically. Therefore, the plasmons are the modes of the dispersion of Eq. (35) approaches ω= (cid:112)1 v q. This collec- F 2 combined system. If the effective carrier mass in the semi- tivemodeisneutralbecauseelectronsandholesoscillatein conductor is larger than that in graphene, then the initial phase. It is similar to acoustic plasmons observed in semi- Drude weight is low but as a result of tunneling, the com- conductors [23, 24]. Incidentally, the plasmons in Ref. [41] binedDrudeweightofthecarriersinthesystem(andhence, were referred to as acoustic because their dispersion was the electriccurrent) wouldincrease. Theupper limit forthe changedfromthesquare-roottoalinearoneduetoscreen- amplificationfactorcanbeestimatedbycompletelyneglect- ing by a nearby gate. This is unlike the original meaning ing the damping, Γ→0, in the expression for the charge of the term acoustic plasmon (or “demon”) introduced for a densityamplitude system where the screening is by electrons from a different bandofthesamematerial[1]. ρq(t)∝e−21(κ+Γ)t, κ<0. (9) For q(cid:191)q ,thesecondterminthesquarerootofEq.(35) c (cid:112) The amplification is proportional to the square root of the dominates, so ω ∝ q. In this case the energy wave is no q plasmon frequency, or the fourth root of the Drude weight. longer neutral: it involves both energy and charge density If the increase of the latter comes from the decrease of the oscillations. However,itisdifferentfromtheplasmon. First, effective mass by, say, a factor of two, then plasmon ampli- the energy wave is in the hydrodynamic regime ω(cid:191)τ−1 ee fication by as much as ∼20% may be possible. For more while the plasmon is in the high frequency regime ω(cid:192)τ−1. ee elaborate estimates, the carrier dynamics beyond the sim- ple Drude approximation would need to be included in the model (see, for example, Ref. 39 and the theory references (a) Pump (b) 𝜔 plasmon citedtherein). SC The tunneling time of hot electrons across ultrathin hBN hBN layers can as short as 7fs [40], which would correspond to κ perhaps as high as several tens of ps−1. In compari- G son, the damping rate in hBN-encapsulated graphene was 𝑞 (c) 𝜔 found to be Γ∼2 and 20ps−1 before and after the opti- demon cal pump, respectively [4]. Hence, fulfilling the condition κ>γ may be feasible. Since the semiconductor would par- tially absorb the pump pulse, graphene may remain rela- tively cool, which may help reduce the plasmon damping SC G 𝑞 duetoelectron-phononscattering[10]. Experimental investigation of the frequency, amplitude, FIG.3. (a)Sketchoftheheterostructuremadeofparallelultrathin and spatial interference patterns of the amplified plasmons layers of semiconductor (SC), insulator, and graphene (G) where as a function of time may be possible by far- and near-field plasmon amplification can occur when hot electrons tunnel from pump-probeopticaltechniques[2,4,9]. the semiconductor to graphene. (b) A qualitative change of the Energy wave (demon) in graphene. Our second example plasmondispersionduringthetunneling. (c)Aqualitativechange ofthedemondispersionwithincreasingtemperature. of a collective mode that may exhibit adiabatic amplifica- 4 (In practice, the range of admissible q is also limited from other. For example, for graphene plasmons we find T (t)∼ m below by the inverse mean-free path l due to electron- ω (t)e−Γt ∼[T(t)]1/2e−Γt. For graphene energy wave, the ph q phonon and disorder scattering, see Fig. 2.) Second, the sameargumentyieldsT (t)∼ω (t)e−Γt∼[T(t)]−3/2e−Γt. m q frequency of the plasmon increases with electron tempera- Insummary,weproposedtheconceptofadiabaticampli- tureT [similartowhatisshowninFig.3(b)]whilethatofthe ficationofcollectivemodesinnonequilibriumsystemsunder demonsdecreases[Fig.3(c)]. Wewillfocusonthesmallq re- photoexcitationandsuggestedtwopossibleroutesforitsex- gionoftheenergywavewhereitsfrequencycanbeefficiently perimentalrealizationin2Dmaterials. Althoughwefocused controlled by T. For µ∼40meV≈500K and T =3000K, on systems with hot electrons, the concept of adiabatic am- which is the regime probed in a recent experiment [4], the plificationisalsoapplicabletosystemswith“cold”electrons, wavelength corresponding to momentum qc is about 1µm. for example, superconducting films [42, 43], where plasmon The change of temperature causes the change of the energy dampingcanbeevensmallerthaningraphene. densityn ,whichinturnaffectsthefrequencyoftheenergy E ThisworkissupportedbytheU.S.DepartmentofEnergy wave through Eq. (35). We assume that this change is adia- underGrantDE-SC0012592. WethankM.K.LiuandG.Ni batic, in other words, that the decay rate κ=−1∂ lnn is 2 t E for discussions and also R. D. Averitt, B. N. Narozhny, and a small parameter. Keeping only the leading terms in the M.I.Stockmanforcommentsonthemanuscript. hydrodynamic equations, we get the WKB solutions for the chargedensityn andtheenergydensity[36] q n (t)∝[n (t)]−1/4e−iS(t), (11) q E n (t)∝[n (t)]3/4e−iS(t), (12) [1] D.Pines,Can.J.Phys.34,1379(1956). Eq E [2] K. F. MacDonald, Z. L. Sámson, M. I. Stockman, and N. I. with S(t) given by Eq. (4). Therefore, if we want to in- Zheludev,Nat.Photonics3,55(2008). crease the energy density oscillations, we need to increase [3] H. Baida, D. Mongin, D. Christofilos, G. Bachelier, A. Crut, the average energy density n ; in other words, we need to P. Maioli, N. Del Fatti, and F. Vallée, Phys. Rev. Lett. 107, E heat up graphene. This can be done using, for example, a 057402(2011). moderate intensity laser source that heats the sample faster [4] G. X. Ni, L. Wang, M. D. Goldflam, M. Wagner, Z. Fei, A. S. McLeod, M. K. Liu, F. Keilmann, B. Özyilmaz, A. H. Castro than the characteristic time τ of electron scattering by ph Neto,J.Hone,M.M.Fogler, andD.N.Basov,Nat.Photonics phonons and disorder. According to Eq. (12), the naive up- 10,244(2016). perboundfortheamplificationfactor(neglectinganydamp- [5] O. Hess, J. B. Pendry, S. A. Maier, R. F. Oulton, J. M. Hamm, ing) is (T/T )9/4 where T is the electron temperature after andK.L.Tsakmakidis,Nat.Mater.11,573(2012). l the photoexcitation. The amplification is only possible if [6] M. I. Stockman, in Active Plasmonics And Tuneable Plasmonic graphene is slightly doped, in which case the energy mode Metamaterials,Vol.32,editedbyA.V.ZayatsandS.A.Maier is not purely neutral. Hence, it can also be probed by opti- (Wiley,2013)pp.1–39. [7] H. R. Seren, J. Zhang, G. R. Keiser, S. J. Maddox, X. Zhao, cal pump-probe spectroscopy, at THz frequencies. Alterna- K. Fan, S. R. Bank, X. Zhang, and R. D. Averitt, Light Sci tively,itmaybepossibletoexploitcouplingofthismodeto Appl.5,e16078(2016). phononsandprobeitbyinelasticlightscattering,similarto [8] J. Chen, M. Badioli, P. Alonso-González, S. Thongrattanasiri, acousticplasmonsinsemiconductors[23,24]. F.Huth,J.Osmond,M.Spasenovic´,A.Centeno,A.Pesquera, Three-temperature state. A solid-state system exhibiting P. Godignon, A. Z. Elorza, N. Camara, F. J. García de Abajo, adiabatic amplification of collective modes would have an- R.Hillenbrand, andF.H.L.Koppens,Nature487,77(2012). otherinterestingnonequilibriumproperty. Itwouldhavenot [9] Z. Fei, A. S. Rodin, G. O. Andreev, W. Bao, A. S. McLeod, M. Wagner, L. M. Zhang, Z. Zhao, M. Thiemens, twobutthreedifferenteffectivetemperatures. Inadditionto G. Dominguez, M. M. Fogler, A. H. Castro Neto, C. N. Lau, the lattice temperature Tl and the electron temperature T, F.Keilmann, andD.N.Basov,Nature487,82(2012). itwouldalsohavethecollectivemodetemperatureTm. The [10] A.Woessner,M.B.Lundeberg,Y.Gao,A.Principi,P.Alonso- temperatureT characterizesthemodescreatedbyrandom González,M.Carrega,K.Watanabe,T.Taniguchi,G.Vignale, m thermal fluctuations rather than those induced by an ex- M.Polini,J.Hone,R.Hillenbrand, andF.H.L.Koppens,Nat. ternal probe pulse. In the absence of damping, Γ=0, the Mater.14,421(2015). time evolution of T can be deduced from the principle of [11] J.D.CoxandF.JavierGarcíadeAbajo,Nat.Commun.5,5725 m (2014). entropy conservation in an adiabatic process. The entropy [12] S.A.Mikhailov,Phys.Rev.B90,241301(2014). of an ensemble of identical harmonic oscillators depends [13] M. Wagner, Z. Fei, A. S. McLeod, A. S. Rodin, W. Bao, only on the ratio of temperature and their mode frequency. E. G. Iwinski, Z. Zhao, M. Goldflam, M. Liu, G. Dominguez, The damping introduces an additional factor e−Γt. There- M. Thiemens, M. M. Fogler, A. H. Castro Neto, C. N. Lau, fore, we expect Tm(t) ≈ ωq(t)e−Γt. [Here we assume that T S.Amarie,F.Keilmann, andD.N.Basov,NanoLett.14,894 Tm(0) ωq(0) m (2014). is still much larger than the final equilibrium temperature [14] D. Brida, A. Tomadin, C. Manzoni, Y. J. Kim, A. Lombardo, Tm(t =∞)=Tl.] When the hot-electron state is just cre- S.Milana,R.R.Nair,K.S.Novoselov,A.C.Ferrari,G.Cerullo, ated,T =T (0)andT =T(0)shouldpresumablybeofthe andM.Polini,Nat.Commun.4,1987(2013). m m same order. Thereafter, they would diverge from one an- [15] F.Rana,IEEETrans.Nanotechnol.7,91(2008). 5 [16] V.ApalkovandM.I.Stockman,LightSci.Appl.3,e191(2014). [47] B. N. Narozhny, I. V. Gornyi, M. Titov, M. Schütt, and A. D. [17] R.J.Morton,A.W.Hood, andR.Erdélyi,Astron.Astrophys. Mirlin,Phys.Rev.B91,035414(2015). 512,A23(2010). [18] P.F.Schmit,I.Y.Dodin, andN.J.Fisch,Phys.Rev.Lett.105, 175003(2010). [19] L.P.Grishchuk,Sov.Phys.JETP67,825(1975). [20] T. V. Phan, J. C. W. Song, and L. S. Levitov, “Ballistic Heat Transfer and Energy Waves in an Electron System,” arXiv:1306.4972(unpublished). [21] U. Briskot, M. Schütt, I. V. Gornyi, M. Titov, B. N. Narozhny, andA.D.Mirlin,Phys.Rev.B92,115426(2015). [22] J.Ruvalds,Adv.Phys.30,677(1981). [23] A.Pinczuk,J.Shah, andP.A.Wolff,Phys.Rev.Lett.47,1487 (1981). [24] P. Padmanabhan, S. M. Young, M. Henstridge, S. Bhowmick, P.K.Bhattacharya, andR.Merlin,Phys.Rev.Lett.113,027402 (2014). [25] R.A.SunyaevandY.B.Zeldovich,Astrophys.SpaceSci.7,3 (1970). [26] M.I.Stockman,Phys.Rev.Lett.93,137404(2004). [27] A. N. Grigorenko, M. Polini, and K. S. Novoselov, Nat. Pho- tonics6,749(2012). [28] D.Basov,M.Fogler,A.Lanzara,F.Wang, andY.Zhang,Rev. Mod.Phys.86,959(2014). [29] F.J.GarcíadeAbajo,ACSPhotonics1,135(2014). [30] O.Vafek,Phys.Rev.Lett.97,266406(2006). [31] D.Svintsov,V.Vyurkov,S.Yurchenko,T.Otsuji, andV.Ryzhii, J.Appl.Phys.111,083715(2012). [32] A.B.Kashuba,Phys.Rev.B78,085415(2008). [33] M. Müller, L. Fritz, and S. Sachdev, Phys. Rev. B 78, 115406 (2008). [34] M. Schütt, P. M. Ostrovsky, I. V. Gornyi, and A. D. Mirlin, Phys.Rev.B83,155441(2011). [35] J.OrensteinandJ.S.Dodge,Phys.Rev.B92,134507(2015). [36] See Supplemental Material at [URL to be inserted by pub- lisher]fortechnicaldetails. [37] M. Massicotte, P. Schmidt, F. Vialla, K. Watanabe, T.Taniguchi,K.J.Tielrooij, andF.H.L.Koppens,NatCom- mun7,(2016). [38] Z.Fang, Y.Wang, Z.Liu, A.Schlather, P.M.Ajayan, F.H.L. Koppens,P.Nordlander, andN.J.Halas,ACSNano6,10222 (2012). [39] R.Huber,F.Tauser,A.Brodschelm,M.Bichler,G.Abstreiter, andA.Leitenstorfer,Nature414,286(2001). [40] Q.Ma,T.I.Andersen,N.L.Nair,N.M.Gabor,M.Massicotte, C. H. Lui, A. F. Young, W. Fang, K. Watanabe, T. Taniguchi, J. Kong, N. Gedik, F. H. L. Koppens, and P. Jarillo-Herrero, NaturePhys.12,455(2016). [41] P.Alonso-González,A.Y.Nikitin,Y.Gao,A.Woessner,M.B. Lundeberg, A. Principi, N. Forcellini, W. Yan, S. Velez, A. J. Huber, K. Watanabe, T. Taniguchi, L. E. Hueso, M. Polini, J.Hone,F.H.L.Koppens, andR.Hillenbrand,“Ultra-confined acoustic THz graphene plasmons revealed by photocurrent nanoscopy,”arXiv:1601.05753(unpublished). [42] O.Buisson,P.Xavier, andJ.Richard,Phys.Rev.Lett.73,3153 (1994). [43] H. T. Stinson, J. S. Wu, B. Y. Jiang, Z. Fei, A. S. Rodin, B. C. Chapler, A. S. McLeod, A. Castro Neto, Y. S. Lee, M. M. Fogler, andD.N.Basov,Phys.Rev.B90,014502(2014). [44] SeeSupplementalMaterial[URLtobeinsertedbypublisher], whichincludesRefs.[45-47]. [45] M.MüllerandS.Sachdev,Phys.Rev.B78,115419(2008). [46] M. Müller, J. Schmalian, and L. Fritz, Phys. Rev. Lett. 103, 025301(2009). 1 Supplementarymaterialfor“Adiabaticamplificationof We define the instantaneous Drude weight D(t ) from the 0 plasmonsanddemonsin2Dsystems” condition that the current generated at time t +0 by the 0 electric field E(t ) equals 1D(t )E(t ). Combined with the 0 π 0 0 definition(Eq.(7)ofthemaintext)ofκ(t),thisimplies CONDUCTIVITYKERNEL D(t0,t0)=D(t0)=D(0)e−2(cid:82)0t0κ(t(cid:48))dt(cid:48). (8) The electric current induced in a system by a weak ex- ternal field can be derived from the Boltzmann transport Let us first examine a single-layer system, N =1, where equation. Consider a general 2D system consisting of N theaveragequasiparticlevelocityv(k,t)issimplyv(1)(k)and layers. Let f(i)(k,t) be the electron distribution function doesnotchangewith t. Supposethissystemiscoolingafter of layer i. The total distribution function f(k) is the sum photoexcitation by emitting acoustic phonons. A common f(k,t)=(cid:80)N f(i)(k,t). Function f(k) can be decomposed wisdom is that the phonon emission is much more effective i=1 into partial waves of different angular momenta. When inrelaxingthemomentumdistributionoftheelectronsthan the external electric field is uniform and the unperturbed incoolingthem. Indeed,atT (cid:192)µ,Pauliblockingofthefinal Hamiltonianisisotropic,onlythes-andp-wavescontribute: electronstatesisunimportantandthetypicalmomentumof f(k)=f (k)+f (k). Theangle-independentpart f contains anemittedphononisoftheorderoftheelectronmomentum s p s the information about the nonequlibrum process, e.g., the k∼T/ħv . Ontheotherhand,theenergyofsuchaphonon F change of electron temperature T(t) and chemical potential Tv /v is much smaller than the typical electron energy ph F µ(t). The p-wave part fp determines the total electric cur- T because the sound velocity vph is much smaller than the rent of the system. In the linear-response regime fp is a Fermi velocity vF. Accordingly, γ should be significantly small correction to fs. The linearized Boltzmann equation larger than the cooling rate 2κ. It is then sensible to split γ hastheform asfollows: ∂tf −eE∂kf =−γ0[fs−f0(t)]−γ(t)fp, (1) γ=Γ+2κ. (9) whereγ istherelaxationrateofthes-wavepart f ,function Wecantransformthetwo-timeDrudeweightto 0 s f (t) is the Fermi-Dirac distribution defined by T(t), µ(t), 0 and γ(t) is the relaxation rate of the p-wave part due to D(t,t )=D(t )e−s(t,t0) 0 0 disorder, electron-phonon, and electron-electron scattering. (10) Keeping the leading-order terms in the external field, we =D(t0)e−2(cid:82)tt0κ(t(cid:48))dt(cid:48)−(cid:82)tt0Γ(t(cid:48))dt(cid:48). obtainseparateequationsforthetwopartialwaves: D(t,t0)=D(t)e−(cid:82)tt0Γ(t(cid:48))dt(cid:48), (11) ∂tfs=−γ0[fs−f0(t)], (2) whichyieldsEq.(6)ofthemaintext. ∂ f =eE∂ f −γ(t)f . (3) Consider next a multilayer system where electrons can t p k s p tunnelbetweenadjacentlayers. Ifthetunnelingisthemajor If γ is the fastest rate in the problem, the approximate 0 mechanism affecting the electron distribution, and electron solutionsoftheseequationsare fs=f0(t)and momentumkisconservedintunneling,thenitismorenat- ural to set γ=Γ instead of Eq. (9). The total distribution t f (t)= (cid:90) eE(t )∂ f (t )e−s(t,t0)dt , (4) function f0 =(cid:80)i f0(i) does not depend on t0. However, if p 0 k 0 0 0 the effective carrier mass of the layers is different, the layer- −∞ averaged quasiparticle velocity v(k,t) [Eq. (7)] changes with where s(t,t )=(cid:82)t γ(t(cid:48))dt(cid:48) is the accumulated damping ex- t. TheresultforD(t,t0)canbewrittenintheformofEq.(11). ponent. The0 currte0ntisgivenby In theory, there is still another case. When the electron system is heated due to absorption of photons, its energy ∞ (or temperature) increases while its momentum remains un- (cid:90) j(t)=−e(cid:88)v(k,t)f (k,t)≡ σ(t,t )E(t )dt , changed. If the situation is possible where this photoexcita- p 0 0 0 tion does not activate other degrees of freedom (for exam- k −∞ ple, phonons) that cause rapid momentum relaxation, then where we may set γ=Γ, as in the case of tunneling. However, the 1 two-time Drude weight D(t,t0) increases with time due to σ(t,t0)=πD(t,t0)θ(t−t0), (5) f0(t0), as in the very first case considered. As a result, we π get D(t,t )=− (cid:88)e2v(k,t)∂ f (t )e−s(t,t0), (6) 0 k 0 0 2 k N D(t,t0)=D(t0)e−(cid:82)tt0Γ(t(cid:48))dt(cid:48), (12) v(k,t)=[∂ f (k,t)]−1(cid:88)v(i)(k)∂ f(i)(k,t). (7) k 0 i=1 k 0 incontrasttoEq.(11). 2 EQUATIONOFMOTIONFORPLASMONS where j ≡ 〈ε(k)v(k)〉 = (n +P)u is the energy current, E E P = 1n is the pressure, η and ζ are the shear and bulk 2 E If the system has the spatial translational symmetry, dif- viscosities, σ is the conductivity, and the angular brackets Q ferent momenta are decoupled. For a given q, in the mean the integral of a quantity over electron momenta k linear-responseregime,theplasmondynamicsisdetermined withtheweightequaltotheshiftedFermidistributionfunc- by the electrical conductivity operator σˆq with the kernel tion f =f0(µ,T,(cid:178)−ku). Inequilibrium[u≡0;µ(r,t),T(r,t)= σq(t,t0), which relates the Fourier components of the cur- const] the electron concentration n and energy density nE rentj andtheelectricfieldE : havethefollowinganalyticalform q q ∞ ∞ (cid:90) (cid:90) jq=σˆqEq≡ σq(t,t0)Eq(t0)dt0. (13) n= (cid:163)f0(µ,T,(cid:178))−f0(0,0,(cid:178))(cid:164)g((cid:178))d(cid:178) −∞ −∞ In turn, the charge density ρq, the electric potential Φq, = 2 T2 (cid:104)π2+1µ2 +2Li (cid:161)−e−µ/T(cid:162)(cid:105), (25) and the current j are related by the Coulomb law and the πħ2v2 6 2T2 2 q F continuityequation: ∞ (cid:90) Φ =v ρ ,E =−iqΦ ,∂ ρ=−iqj , (14) nE= (cid:163)f0(µ,T,(cid:178))−f0(0,0,(cid:178))(cid:164)(cid:178)g((cid:178))d(cid:178) q q q q q t q −∞ whereallthequantitiesarefunctionsoftimeand vq=2π/(cid:178)q (15) =π2ħT2v32(cid:104)π32 Tµ+13 Tµ33−4Li3(cid:161)−e−µ/T(cid:162)(cid:105). (26) F is the Fourier transform of the bare Coulomb potential In these expression, g((cid:178))=(2/π)(|(cid:178)|/ħ2v2) is the electron screened by the dielectric environment. (Note that screen- F density of states of graphene, Li (x) is the polylogarithm ing by the electrons themselves should not be included in (cid:178) z function, and f (0,0,(cid:178))=Θ(−(cid:178)), where Θ(x) is the unit step because Eq is the total electric field.) Equations (13) and (14) function. For w0eakly doped graphene (or for high T), T (cid:192) entail |µ|,onefinds ∂tρq+q2vqσˆqρq=0, (16) 4ln2 µT whichisthesameas n(cid:39) π ħ2v2 , (27) F ∞ (cid:90) 6ζ(3) T3 ∂tρq+q2vq σq(t,t0)ρq(t0)dt0=0. (17) nE(cid:39) π ħ2v2 , (28) −∞ F Withthefollowingmodelfortheconductivitykernel: whereζ(3)=1.202istheRiemannzeta-function. Notethatif n is fixed, which is usually the case in the experiment, then σ (t,t )= 1D(t)e−Γ(t−t0)θ(t−t ), (18) Eq.(25)implicitlydefinesµisafunctionofT. Thisfunction q 0 π 0 µ=µ(n,T) can be found by solving Eq. (25) numerically or D(t)=D(0)e−2(cid:82)0tκ(t(cid:48))dt(cid:48), (19) (totheleadingorder)Eq.(27)analytically, and taking the time derivative of Eq. (17), and combining it π n withEq.(18),weget µ(cid:39)4ln2 T ħ2vF2, (29) (cid:104) (cid:105) ∂2+γ(t)∂ +ω2(t) ρ =0, (20) t t q q with γ=2κ+Γ, (21) whichistheresultannouncedinthemaintext. EQUATIONOFMOTIONFORENERGYWAVES (DEMONS) Thelinearizedhydrodynamicequations[20,21,33,45–47] (inthenotationsofRef.21)are 1 (cid:181) 1 µ(cid:182) FIG. 1. Dispersion of the energy wave (demon) in graphene at ∂tn+∇(nu)=− σQ∇ E+ T∇ , (22) differentT calculatedfromEqs.(25),(26),and(35). Thereddashed e e T ∂ n +∇j =0, (23) line is the infinite temperature limit, ω= (cid:112)12vFq. The electron t E E concentrationisn=2.0×1012cm−2 andthedielectricconstantis ∂tjE+vF2∇P=−envF2E+η∇2u+ζ∇(∇u), (24) (cid:178)=1. 3 Having obtained µ, one can use Eq. (26) to compute the Therefore, the instantaneous frequency of the energy wave (cid:112) energydensityn =n (n,T)fromEqs.(26)or(28). (or“demon”)isω = c. Amoreaccurateexpression[20]is E E q If the electronic temperature T(t) is uniform but slowly obtainedifthenext-orderin q termsareretained: changing, the linearized equations for the Fourier harmon- (cid:115) ics of the concentration, energy density, and drift velocity 1 4πe2 n2 ω =v q2+ q, (35) become q F 2 3 (cid:178) n E which is Eq. (10) of the main text. Formula (35) predicts ∂tnq+iqnuq=0, (30) the crossover from (cid:112)q to linear in q behavior, which was 3 discussed therein. Representative plots of the energy wave ∂ n + iqn u =0, (31) t Eq 2 E q dispersion illustrating its dependence on T are shown in iqv2(cid:179)nEq +e2nv n (cid:180)+3(n ∂ +∂ n )u =0. (32) Fig.1. F 2 q q 2 E t t E q If n is slowly time-dependent, so are the coefficients b E and c in Eq. (33). This differential equation can be solved In these equations we neglected the dissipative terms be- within the WKB approximation, which yields Eqs. (10) and cause they are quadratic in the small parameter q. Simi- (11)ofthemaintext. larly,wekeptonlytheleadingtermsintheadiabaticchang- Compared with previous work, our results are in full ing rate κ. From these equations we can get the third-order agreement with those of Ref. 20. An equation similar to differentialequationforn alone: Eq. (35) is Eq. (43) of Ref. 21, however the coefficient for the q term linear in q differs from that in Eq. (35) by a factor of (cid:161)∂3+b∂2+c∂ (cid:162)n (t)=0, (33) 2π. Acollectivemodewiththeacousticdispersionwasalso t t t q discussedinRef.31. Althoughthestartingequationsofthat 2 n2q2 workaremathematicallyequivalenttoourEqs.(22)–(24),the b=2∂ lnn ,c= e2v2v . (34) (cid:112) t E 3 F q nE finalresultforthevelocityis0.6vF insteadof vF/ 2.