Plasmon Generation through Electron Tunneling in Graphene
PPlasmon Generation through Electron Tunneling in Graphene
Sandra de Vega and F. Javier Garc´ıa de Abajo
1, 2 ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute ofScience and Technology, 08860 Castelldefels (Barcelona), Spain ICREA-Instituci´o Catalana de Recerca i Estudis Avan¸cats,Passeig Llu´ıs Companys 23, 08010 Barcelona, Spain ∗ (Dated: November 10, 2018) The short wavelength of graphene plasmons relative to the light wavelength makesthem attractive for applications in optoelectronics and sensing. However, this propertylimits their coupling to external light and our ability to create and detect them. Moreefficient ways of generating plasmons are therefore desirable. Here we demonstratethrough realistic theoretical simulations that graphene plasmons can be efficiently ex-cited via electron tunneling in a sandwich structure formed by two graphene monolay-ers separated by a few atomic layers of hBN. We predict plasmon generation rates of ∼ − / s over an area of the squared plasmon wavelength for realistic values of thespacing and bias voltage, while the yield (plasmons per tunneled electron) has unityorder. Our results support electrical excitation of graphene plasmons in tunnelingdevices as a viable mechanism for the development of optics-free ultrathin plasmonicdevices.I. INTRODUCTION Plasmons offer the means to concentrate optical fieldsdown to nanometer-sized regions and enhance the inten-sity of externally incident light by several orders of mag-nitude [1, 2]. These properties have been extensively in-vestigated to develop applications in areas as diverse asoptical sensing [3–6], medical diagnosis and treatment [7–9], nonlinear optics [10, 11], catalysis [12–15], and pho-tovoltaics [16, 17]. In this context, plasmons in high-quality graphene offer the advantages of being electri-cally [18–23], magnetically [24, 25], and optically [26, 27]tunable, while exhibiting short wavelengths and long life-times [28]. Graphene plasmons are highly customizablethrough lateral patterning, stacking of 2D-crystal atomiclayers, and coupling to the dielectric environment [29],enabling applications such as light modulation [22, 30],control of thermal emission [31–34], spectral photome-try [35, 36], and optical sensing [37–40]. While thesephenomena and applications have been mainly exploredat mid-infrared and lower frequencies, prospects for theirextension to the visible and near-infrared spectral regionslook promising [26]. Actually, there is no fundamentalreason that limits the existence of graphene plasmons atsuch high frequencies, although in practice they requirelarge doping ( e.g. , an attainable E F ∼ <
10 nm. Although thelatter is challenging using top-down lithography, excitingpossibilities are opened by bottom-up approaches to thesynthesis of nanographenes with only a few nanometersin lateral size [42, 43], while recent experiments on self-assembled graphene nanodisks already reveal plasmonsclose to the near-infrared [44]. In fact, plasmon-like res- ∗ Electronic address: [email protected] onances have been measured in small aromatic hydro-carbons [45, 46], which can be considered as molecularversions of graphene.The strong confinement of graphene plasmons is clearlyindicated by the ratio of plasmon-to-light wavelengthsfor a self-standing carbon monolayer [26], λ p /λ =2 α ( E F / ¯ hω ) (cid:28)
1, where α ≈ /
137 is the fine struc-ture constant and ω is the light frequency. This rela-tion implies that graphene plasmons can be describedin the quasistatic limit. Unfortunately, it also meansthat their in/out-coupling to propagating light is weak.Let us emphasize that far-field coupling in homogeneousgraphene is forbidden because plasmons are evanescentsurface waves with energy and momentum outside thelight cone; however, we refer here to the limited abilityof graphene nanostructures to produce in/out plasmoncoupling, which is neatly illustrated by their optical ex-tinction and elastic-scattering cross-sections. Indeed, asimple extension of previous analytical theory [26] basedon the Drude conductivity indicates that the radiativedecay rate of a low-order plasmon in a graphene struc-ture ( e.g. , a disk) is κ r ∼ cAE F /λ E p , which for a plas-mon energy E p comparable with the Fermi energy E F ,a graphene area A ∼ λ ∼ − λ , and a characteristicphoton wavelength of 2 µ m, leads to κ r ∼ . / ns. This isorders of magnitude smaller than the total plasmon decayrate κ . Now, even for an unrealistically optimistic value κ = 1 / ps, the on-resonance plasmon-driven extinction( σ ext = (3 λ / π ) κ r /κ ∼ − λ ) and elastic-scattering( σ scat = ( κ r /κ ) σ ext ∼ − λ ) cross-sections are small,thus rendering plasmon coupling to propagating light as amere retardation correction. In order to circumvent thislimitation, many experimental studies in this field haverelied on near-field nanoscopy [47], which is based on theuse of sharp tips to enhance this coupling [19, 20, 27].However, more compact devices for the generation anddetection of graphene plasmons are needed in order to a r X i v : . [ c ond - m a t . m e s - h a ll ] A ug enable the design and development of applications in in-tegrated architectures.A promising approach consists in exploiting inelas-tic electron transitions to excite and detect plasmons.Understandably, focused beams in electron microscopeshave been the probe of choice to excite and map plasmonswith nanometer spatial resolution via the recorded elec-tron energy-loss and cathodoluminescence signals [48].Electron tunneling has been considered for a long time asa mechanism for the excitation of plasmons [49–57], whileelectron injection by tunneling from conducting tips intometallic structures has been also demonstrated to pro-duce efficient plasmon excitation [58–62]. More dedicateddesigns have incorporated emitting devices in which thegenerated light directly couples to the plasmon near-field[63, 64]. A recent theoretical study has shown that plas-mons can boost tunneling across an insulator separatingtwo graphene layers [65], with potential use as a plasmon-gain device [66], while a similar structure has been exper-imentally used for the electrical generation of THz radia-tion [67]. Additionally, evidence of radiative gap-plasmondecay has been experimentally obtained associated withhot electron tunneling under external illumination [68].In a related context, efforts to realize electrical detec-tion of plasmons have relied on near-field coupling tostructures consisting of organic diodes [69], supercon-ductors [70], nanowire semiconductors [71, 72], Schot-tky barriers [73–75], and 2D semiconductors [76]. Ther-moelectric detection has been recently demonstrated ingraphene [36], while a nanoscale detector has been pro-posed based on changes produced in the conductivity ofgraphene nanojunctions [35]. These technologies shouldbenefit from advances in the design, nanofabrication, andtheoretical modeling of transistors made from 2D mate-rial sandwiches [77–79].Here, we theoretically investigate a simple structurein which two closely spaced graphene sheets serve bothas gates to produce electron tunneling and as plasmon-supporting elements. We predict a generation efficiencythat approaches one plasmon per tunneled electron. Con-sidering attainable doping conditions and bias voltagesapplied to the graphene layers, we find a generation rateof ∼ − plasmons per second over an area of asquared plasmon wavelength for a 1 nm inter-graphenespacing filled with hexagonal boron nitride (hBN). Thistype of structure provides a key element for future optics-free integrated devices and could also be operated in re-versed mode to detect plasmons by decay into tunneledelectrons. II. PLASMON GENERATION THROUGHELECTRON TUNNELING
The structure under consideration (Fig. 1a) consists oftwo graphene layers separated by a hBN film and placedat different Fermi energies E F1 and E F2 . Upon applica-tion of a bias voltage V b between the two carbon layers, Diracpoint graphenegraphene plasmon BN (a)(b) BN electrons FIG. 1:
Electron Tunneling in Double-LayerGraphene . (a) We consider a graphene-hBN-graphenesandwich structure consisting of two crystallographically-aligned monolayer graphene sheets separated by a hBNfilm of thickness d . Electron tunneling takes place uponapplication of a bias voltage V b . The graphene layers are de-scribed through their wave-vector- and frequency-dependentconductivities σ ( k (cid:107) , ω ) and σ ( k (cid:107) , ω ). The electric-fieldintensity associated with the symmetric optical plasmonsustained by this structure is plotted schematically. (b) Weshow the electronic bands of the graphene layers and theirrelative energy positions due to the bias. The Fermi energies E F1 and E F2 and the work function Φ are indicated. Theelectron wave functions have a dependence along the filmnormal direction modeled through two quantum wells (blueprofile), while Dirac fermions describe their dependenceon parallel directions (red cones and parallel electron wavevector Q ). Electron tunneling (wiggly arrow) is limited byparallel momentum conservation. electrons can tunnel from one to the other, assisted bythe excitation of a propagating plasmon (inelastic tunnel-ing). We schematically depict the profile of a symmetricplasmon superimposed on the scheme ( i.e. , an excitationto which the tunneling electrons couple with high effi-ciency, see below).Elastic tunneling ( i.e. , direct tunneling of electronsfrom one layer to the other) is forbidden due to the ver-tical displacement of their respective Dirac cones associ-ated with the bias, which introduces a mismatch of par-allel wave vectors for all electron energies (see Fig. 1b).In contrast, the creation of a plasmon (wiggly arrow)can break this mismatch and enable inelastic tunnelingof the electron. Nonetheless, this process involves only aspecific energy in the graphene electron bands for eachinelastic frequency ω , which is controlled through theapplied bias voltage and the Fermi energies of the twocarbon layers. The theoretical investigation of plasmonexcitation during electron tunneling has been pioneeredin the context of junctions between semiconductors us-ing the self-energy [52]. A recent theoretical study[65]discusses in great detail plasmon generation associatedwith electron tunneling using a second-quantization for-malism. We follow instead a more direct approach basedupon methods borrowed from the analysis of electronenergy-loss spectroscopy [48].The elements involved in the theoretical description ofinelastic electron tunneling are (i) the conduction elec-tron wave functions; (ii) the screened Coulomb poten-tial W that mediates the interaction between those elec-trons; and (iii) the evaluation of the transition proba-bility through a well-established linear-response formal-ism [48]. We present a a self-contained derivation inthe Appendix, under the assumption of crystallographi-cally aligned graphene sheets ( i.e. , their Dirac cones areon top of each other in 2D momentum space). In par-ticular, electrons are assumed to have their wave func-tions factorized as the product of components parallel(Dirac fermions [80]) and perpendicular to the surface.The latter accounts for the transversal confinement tothe graphene layers, which we describe through a simplepotential-well model (blue profile in Fig. 1b), with pa-rameters fitted to match the work function and electronspill out of p z orbitals in the material (see Appendix).Actually, spill out is a crucial element because it deter-mines the overlap of electronic wave functions betweenthe two graphene layers. For simplicity, we neglect therepulsive electron potential in the hBN region, whichshould produce a correction to this overlap. The inelas-tic tunneling current density is then expressed as an in-tegral over parallel electron wave vectors (see Eq. (A2)),which involves the Fermi-Dirac occupation distributionsof the graphene bands in the two layers (we assume atemperature of 300 K). The screened Coulomb interac-tion enters this integral and incorporates the optical re-sponses of the graphene layers ( i.e. , the correspondingconductivities) and the hBN film (see Appendix). Weuse the random-phase approximation [81–83] (RPA) todescribe the frequency- and wave-vector-dependence ofthe graphene conductivities, thus incorporating nonlocaleffects that are important because both interlayer dielec-tric coupling and electron tunneling involve large valuesof the transferred parallel momentum for the short sepa-rations under consideration. The conductivities dependon the Fermi energies of the graphene layers, as well ason the electron lifetime, which we set to a conservativevalue τ = 66 fs. We remark that our results for the in- tegral of the tunneling current over the plasmon widths( ∼ τ − ) is rather independent of the choice of τ ( i.e. , thisproperty is inherited from the negligible dependence ofthe frequency integral of the screened interaction Im { W } over the plasmon width). III. RESULTS AND DISCUSSIONA. Plasmons in the Double-Layer GrapheneStructure
The optical response of the sandwich structure inFig. 1a is dominated by graphene plasmons and opti-cal phonons of hBN. The frequency- and wave-vector-dependence of these excitations is illustrated in the dis-persion diagram of Fig. 2a, which shows the imaginarypart of the Fresnel coefficient for p polarization and in-cidence from layer 1 (see Eq. (D1)) with a representa-tive choice of interlayer distance d = 1 nm and grapheneFermi energies E F1 = 1 eV and E F2 = 0 . E F1 (cid:54) = E F2 . In thesmaller gap region (1 eV maximum energy), the opticaland acoustic plasmons are both well defined, with theirlifetimes roughly equal to the assumed intrinsic value of τ (see above). The acoustic plasmon then fades away whenits energy increases and it enters the region of interbandtransitions of the low-doping layer (Landau damping).A similar behavior is observed for the optical plasmons,although the described interband absorption effect is sig-nificantly weaker. Examination of the field intensity pro-files associated with both plasmons explains this differ-ent behavior (Fig. 2a, upper inset): the acoustic plasmonhas larger weight on the low-doping layer (2, at z = d ),while the optical plasmon receives substantial contribu-tions from both layers; therefore, the effect of interbandLandau damping produced by layer 1 (at z = 0) actsmore strongly on the acoustic plasmon. Additionally, op- (nm –2 ) -10 -7 (a) (b) (c) o p t i c a l p l a s m o n BN-phonons
AB C a c o u s t i c p l a s m o n FIG. 2:
Energy- and momentum-resolved electron tunneling . (a) Optical dispersion diagram of a graphene-hBN-graphene structure (Fig. 1a), plotted through the imaginary part of the Fresnel reflection coefficient for p polarization andincidence on layer 1 as a function of optical energy ω and parallel wave-vector k (cid:107) . The upper inset shows the electric intensityprofiles associated with acoustic and optical plasmons at the energies plotted in (c). The energy and wave-vector axes arenormalized to the Fermi energy and wave vector of layer 1, E F1 and k F1 = E F1 / ¯ hv F , respectively. (b) Contribution of differentpoints in the dispersion diagram to inelastic electron tunneling J ( k (cid:107) , ω ) (Eq. (A2)). (c) Frequency-dependent momentum-integrated tunneling probability per unit area (horizontal scale) as a function of inelastically transferred energy ¯ hω (verticalscale). We consider a graphene spacing d = 1 nm, Fermi energies E F1 = 1 eV and E F2 = 0 . eV b = 1 . tical plasmons have larger Drude weights than acousticplasmons ( i.e. , they are comparatively more weighted onlayer 1, which has higher doping electron density), andthis contributes to reduce the relative effect produced bycoupling to the interband transitions of layer 2. B. Inelastic Tunneling Probability
The tunneling current density can be decomposed intothe contributions of different parallel-wave-vector andfrequency transfers J ( k (cid:107) , ω ) as J = (cid:90) ∞ dω (cid:90) ∞ dk (cid:107) J ( k (cid:107) , ω ) . We find that J ( k (cid:107) , ω ) (Eq. (A2)) exhibits features thatfollow the hybrid modes of the system ( cf. Figure 2a and2b). The relative strengths of the different excitationsobviously depend on the applied bias voltage. Addition-ally, the plot of J ( k (cid:107) , ω ) is dominated by the conditionimposed by energy and momentum conservation (see Eq.(A2) in the Appendix), which renders a zero contributionoutside the region eV b / ¯ h + v F k (cid:107) < ω < eV b / ¯ h − v F k (cid:107) .This is more clearly illustrated in Figure 3, where weshow | J ( k (cid:107) , ω ) | for different bias voltages, plotted now inlinear scale. The momentum-integrated spectral decom-position of the tunneling current J ( ω ) = (cid:82) ∞ dk (cid:107) J ( k (cid:107) , ω )(Fig. 2c) is thus peaked at frequencies ω ≈ eV b / ¯ h − v F k (cid:107) for values of k (cid:107) determined by the dispersion relations ofthe hBN phonon and the graphene plasmons.The above analysis indicates that the peak frequencyof the excited plasmons can be controlled through theapplied voltage. This is corroborated by Fig. 4a, in whichthe peak intensities are in excellent agreement with thecombination of threshold and mode dispersion relations,as we argue above. Importantly, the bias also affects themagnitude of the tunneling current (Fig. 4a).We are interested in producing a substantial tunnelingcurrent, but ultimately, we need that a sizable fraction ofthe tunneled electrons really generate plasmons. With asuitable choice of the bias voltage, the tunneling currentis dominated by the production of either acoustic or op-tical plasmons (see Fig. 4d). As expected, for small volt-ages, below the thresholds for generation of both acous-tic and optical plasmons, phonons dominate the inelasticcurrent, while for large voltages the plasmon contribu-tion becomes small compared with other channels of in-elastic tunneling associated with interband electron tran-sitions. Therefore, there is an optimum voltage range inwhich the generation yield (plasmons per electron) hasunity order; tunneling in that range is actually domi-nated by the excitation of optical plasmons (see regionnear eV b ∼ E F1 = 1 eV in Fig. 4d).Additionally, we are interested in producing long-livedplasmons. Under the conditions of Fig. 4, for a real-istic bias voltage < / nm, the lifetimes are boosted (a) (b) (c) FIG. 3:
Dispersion diagram of tunneling current components . Contribution of different points in the dispersion diagramto the inelastic electron tunneling J ( k (cid:107) , ω ) (Eq. (A2)), plotted in linear scale for different bias voltages V b (see labels) underthe same conditions as in Fig. 2 ( d = 1 nm, E F1 = 1 eV, and E F2 = 0 . when the plasmons enter the gap region of the low-dopinggraphene layer 2 (Fig. 4c, solid curves), where they reachvalues of the order of the assumed intrinsic lifetime τ ,which depends on material quality and is ultimately lim-ited by acoustic phonons [28].Similar qualitative results are obtained when consid-ering other values of the Fermi energies or when thegraphene layers are separated by a vacuum gap, as shownin the additional plots presented in the Appendix. IV. CONCLUDING REMARKS
In order to place the calculated tunneling current den-sities in context, it is pertinent to consider the currentproduced over an area of a squared plasmon wavelength.As shown in Fig. 4c (right scale), this quantity reachesvalues of ∼ − / s for a realistic graphene spacingof 1 nm, corresponding to roughly three hBN atomic lay-ers. In this respect, the maximum applied voltage beforebreakdown takes place is an important factor to consider,as well as the threshold voltage for plasmon generation,which depends on the Fermi energies of the two graphenelayers. It is important to stress that the generation effi-ciency reaches one plasmon per tunneled electron, whichshould enable the generation of single and few plasmonswhen combined with quantum electron transport setups[85].Electrical detection of graphene plasmons can also beaccomplished within the sandwich structures under con-sideration. A plasmon propagating through the structureor laterally confined in a finite-size island can decay bytransferring its energy to a tunneling electron under theappropriate bias conditions. Indeed, using high-quality graphene and considering plasmons of frequency and mo-mentum in the gap region of both graphene layers, inelas-tic electron tunneling should be the dominant channel ofplasmon decay, leaving an electron in the low-potentiallayer and a hole in the high-potential one. In this way,charge separation is naturally accomplished, thus facil-itating the electrical detection of the plasmon. Thereare two favorable factors that support this possibility ingraphene: the noted momentum mismatch that preventsdirect elastic tunneling (obviously, graphene quality is animportant factor to prevent defect-assisted elastic tunnel-ing); and the high spatial concentration of the plasmons,which result in large coupling to inelastic transitions com-pared with coupling to photons in tunneling-based in-elastic light emission measurements [55]. As an alterna-tive approach, electrical detection could be performed bya separate graphene-based structure via thermoelectricmeasurements [36] or through thermo-optically activatednanoscale junctions [35].These concepts are equally applicable to other van derWaals crystals supporting 2D polaritons ( e.g. , plasmonsin black phosphorous [86, 87], or even optical phononsin hBN). Then, our formalism can be easily adapted tomaterials other than graphene by correcting the electronenergies and matrix elements of Eq. (A2) ( e.g. , usingtight-binding electron wave functions [88]), as well as theconductivities of Eq (C1).The sandwich structure under consideration is conve-nient for the design of integrated plasmonic circuitry.We envision plasmon generation and detection in patcheswith a lateral size of the order of the plasmon wavelength( ∼ −
100 nm for the plasmon energies under considera-tion). The proposed source generates plasmons peaked atbias-dependent frequencies, but further frequency selec-
Bias (eV) e ! / ( s ! n m " ) -10 -8 -6 Bias (eV) e ! / ( s ! n m " ) P l a s m on li f e t i m e ( f s ) La y e r SP f r a c t i on (eV) (a)(b) (eV) (nm) (d) scaling (eV) (c) total A BC
C BA E l e c t r on c u rr en t ( s –1 n m –2 ) optical plasmonacoustic plasmon BN phonons E l e c t r on c u rr en t ( s –1 –2 ) –10 –8 acoustic ( n m –2 ) CB E l e c t r on c u rr en t ( s –1 n m –2 ) –6 optical FIG. 4:
Probability of plasmon-generation by electron tunneling . (a) Spectrally resolved inelastic tunneling probabilityfor different bias voltages. (b)
Bias voltage dependence of the total inelastic tunneling current density for various graphenespacings (see labels), multiplied by the indicated factors for the sake of clarity. (c)
Plasmon lifetimes (solid curves, left scale,assuming an intrinsic RPA input lifetime of 66 fs) and fractions of the plasmon weight in layer 2 (dashed curves, right scale,showing the fraction of squared plasmon-induced charge density in that layer). (d)
Contribution of the BN main optical phonon(A), the acoustic graphene plasmon (B), and the optical graphene plasmon (C) to the total tunneling current density (blackcurve), corresponding to the spectral features labeled in (a), plotted per square nanometer (solid curves, left axis) and persquare plasmon wavelength (dashed curves, right axis). We consider Fermi energies of 1 eV and 0.5 eV for the graphene layers1 and 2 in all cases, while the graphene spacing and bias voltage are indicated by axis labels and text insets. tion could be performed upon transmission of the plas-mons through engineered waveguides. These elementsshould grant us access into optics-free devices, in whichplasmons can perform different operations, for exampleby interacting with quantum dots and localized molec-ular excitations, truly relying on their very nature ascollective electron excitations, without the mediation ofphotons whatsoever. Plasmon-based sensors could thenconsist of a plasmon generation stage, a transmission rib-bon exposed to the analyte, and a plasmon detectionstage, with some degree of spectral resolution possiblyachieved by playing with the graphene Fermi levels andthe applied bias voltages. The versatility and excellentperformance predicted for the double-graphene sandwich structure opens exciting prospects for the design of mod-ular integrated devices with multiple functionalities.
Appendix A: Calculation of Inelastic ElectronTunneling Probabilities
We consider electron transitions between initial i andfinal f states assisted by the creation of inelastic exci-tations in the system. The initial and final states aretaken to be in the graphene layers 1 and 2, repectively.The transition probability Γ can be decomposed intothe contribution of different inelastic frequencies ω asΓ = (cid:82) ∞ dω Γ( ω ), where the spectrally resolve probabilityadmits the expression [48]Γ( ω ) = 2 e ¯ h (cid:88) i,f (cid:90) d r (cid:90) d r (cid:48) ψ ∗ i ( r ) ψ f ( r ) ψ ∗ f ( r (cid:48) ) ψ i ( r (cid:48) )Im {− W ( r , r (cid:48) , ω ) } (A1) × δ ( ε f − ε i + ω ) f (¯ hε i ) [1 − f (¯ hε f )] . Here, ψ i and ψ f are initial and final electron wave func-tions of energies ¯ hε i and ¯ hε f , respectively, f j ( E ) = { E − E F j ) /k B T ] } − is the Fermi-Dirac electrondistribution of the graphene layer j (we set T = 300 K),and W ( r , r (cid:48) , ω ) is the screened interaction, defined asthe electric potential produced at r by a unit chargeoscillating with frequency ω at the position r (cid:48) . Thisresult is rigorous up to first-order perturbation in thescreened interaction W . We assume that the electronwave functions of the system depicted in Fig. 1a canbe factorized as the product of decoupled in-plane andout-of-plane components. The latter, which we approx-imate by a real quantum-well wave function ϕ ⊥ ( z ) (seebelow), is shared by all conduction electrons. The re- maining in-plane wave functions are Dirac fermions [80]characterized by their 2D wave vectors Q . Addition-ally, the translational symmetry of the system allowsus to write the screened interaction as W ( r , r (cid:48) , ω ) =(2 π ) − (cid:82) d k (cid:107) exp[i k (cid:107) · ( R − R (cid:48) )] W ( k (cid:107) , z, z (cid:48) , ω ) in terms ofparallel wave vector components k (cid:107) . Explicit expresionsfor W ( k (cid:107) , z, z (cid:48) , ω ) are offered below. Likewise, the transi-tion rate can be separated as Γ = (cid:82) ∞ dω (cid:82) ∞ dk (cid:107) Γ( k (cid:107) , ω ).We then divide by the sandwich area A to obtain thetunneling current J = Γ /A . As we show in the detailedderivation provided in Sec. F, Eq. (F1) allows us to writethe wave-vector- and frequency-resolved tunneling cur-rent as J ( k (cid:107) , ω ) = e k (cid:107) π ¯ h (cid:90) d Q i (cid:18) Q i · ( Q i − k (cid:107) ) Q i | Q i − k (cid:107) | (cid:19) f (¯ hv F Q i ) [1 − f (¯ hv F | Q i − k (cid:107) | )] (A2) × (cid:90) dz (cid:90) dz (cid:48) ϕ ⊥ ( z ) ϕ ⊥ ( z − d ) ϕ ⊥ ( z (cid:48) ) ϕ ⊥ ( z (cid:48) − d ) Im (cid:8) − W ( k (cid:107) , z, z (cid:48) , ω ) (cid:9) × δ ( v F ( | Q i − k (cid:107) | − Q i ) + ω − eV b / ¯ h ) , where v F ≈ m / s is the graphene Fermi velocity. Thisexpression, which is independent of the direction of k (cid:107) ,includes an integral over initial-state parallel wave vec-tors Q i , as well as a factor of 4 accounting for spinand valley degeneracies. Additionally, the δ function im-poses energy conservation and limits the spectral rangefor which J ( k (cid:107) , ω ) is nonzero to eV b / ¯ h − v F k (cid:107) < ω We approximate the dependence of the graphene elec-tron wave function ϕ ⊥ ( z ) on the out-of-plane direction z by one of the states of a quantum well. More pre-cisely, we consider the first excited well state, in order tomimic the antisymmetry of the p z orbitals in grapheneconduction band. For simplicity, we assume a squarewell potential of depth V and width a . These parame-ters are fitted in such a way that (i) the binding energy ofthe well state coincides with the graphene work functionΦ = 4 . (ii) the centroid of the electron- density spill-out is the same in the well state and in the p z orbital (i.e., (cid:82) dz | z | | ϕ ⊥ ( z ) | = (cid:82) d r | z | | ϕ p z ( r ) | ). Us-ing the atomic carbon 2 p orbital for ϕ p z [90], we find V = 45 eV and a = 0 . 12 nm. This value of a is closeto (but smaller than) the interlayer spacing in graphite(0.335 nm), which is reasonable due to the spill out of ϕ ⊥ outside the well. A detailed derivation of these results isgiven in Sec. E. Incidentally, the graphene-electron spill-out toward the gap must depend on the details of theelectron band structure, which should also vary with thefilling material. In particular, hBN presents a band gapof ∼ Appendix C: Screened Interaction We consider a sandwich formed by two graphene layersof conductivities σ and σ , separated by an anisotropicfilm of dielectric permittivity (cid:15) z and (cid:15) x for directions par-allel ( z ) and perpendicular ( x , y ) to the film normal, re-spectively, and surrounded by vacuum above and belowthe structure. In the electrostatic limit here assumed, thesurface-normal wave vector inside the dielectric can be written as i q with q = k (cid:107) (cid:112) (cid:15) x /(cid:15) z , where we take the signof the square root such that Re { q } > 0. It is also conve-nient to define the effective permittivity (cid:15) = √ (cid:15) x (cid:15) z , takento have positive imaginary part. The screened interactionadmits the closed-form expression (see Appendix) W ( k (cid:107) , z, z (cid:48) , ω ) = W dir ( k (cid:107) , z, z (cid:48) , ω ) + W ref ( k (cid:107) , z, z (cid:48) , ω ) , (C1)where W dir ( k (cid:107) , z, z (cid:48) , ω ) = 2 πk (cid:107) × e − k (cid:107) | z − z (cid:48) | , z, z (cid:48) < z, z (cid:48) > d(cid:15) − e − q | z − z (cid:48) | , < z, z (cid:48) < d , otherwiseis the interaction between charges placed inside the homogeneous vacuum or dielectric parts of the structure, while W ref ( k (cid:107) , z, z (cid:48) , ω ) = (2 π/k (cid:107) )1 − A (cid:48) A (cid:48) e − qd × e k (cid:107) (2 d − z − z (cid:48) ) (cid:2) A + A (cid:48) ( A + B (cid:48) ) e − qd (cid:3) , d < z d < z (cid:48) B e k (cid:107) ( d − z ) (cid:104) e − q ( d − z (cid:48) ) + A (cid:48) e − q ( d + z (cid:48) ) (cid:105) , d < z, < z (cid:48) < d(cid:15)B B e k (cid:107) ( d − z + z (cid:48) ) e − qd , d < z, z (cid:48) < B e k (cid:107) ( d − z (cid:48) ) (cid:2) e − q ( d − z ) + A (cid:48) e − q ( d + z ) (cid:3) , < z < d, d < z (cid:48) (cid:15) − (cid:8) A (cid:48) e − q ( z + z (cid:48) ) + A (cid:48) e − q (2 d − z − z (cid:48) ) + A (cid:48) A (cid:48) (cid:104) e − q (2 d + z − z (cid:48) ) + e − q (2 d − z + z (cid:48) ) (cid:105) (cid:9) , < z < d, < z (cid:48) < dB e k (cid:107) z (cid:48) (cid:2) e − qz + A (cid:48) e − q (2 d − z ) (cid:3) , < z < d, z (cid:48) < (cid:15)B B e k (cid:107) ( d + z − z (cid:48) ) e − qd , z < , d < z (cid:48) B e k (cid:107) z (cid:104) e − qz (cid:48) + A (cid:48) e − q (2 d − z (cid:48) ) (cid:105) , z < , < z (cid:48) < d e k (cid:107) ( z + z (cid:48) ) (cid:2) A + A (cid:48) ( A + B (cid:48) ) e − qd (cid:3) , z < , z (cid:48) < B j = 2 / (1 + (cid:15) + β j ), B (cid:48) j = (cid:15)B j , A j = B j − A (cid:48) j = B (cid:48) j − 1, where β j = 4 π i k (cid:107) σ j /ω . Reassur-ingly, Eq. (C1) explicitly satisfies the reciprocity con-dition W ( k (cid:107) , z, z (cid:48) , ω ) = W ( k (cid:107) , z (cid:48) , z, ω ). In our calcula-tions, we consider two graphene layers separated by ei-ther vacuum ( (cid:15) x = (cid:15) z = 1) or hBN. For hBN, we ap-proximate optical phonons as Lorentzians and write [92] (cid:15) (cid:96) ( ω ) ≈ (cid:15) ∞ ,(cid:96) + (cid:80) i =1 , s (cid:96)i / [ ω (cid:96)i − ω ( ω + i γ (cid:96)i )], with pa-rameters (cid:15) ∞ ,z = 4 . s z = 70 . ω z = 97 . γ z = 0 . 99 meV, s z = 126 meV, ω z = 187 meV, and γ z = 9 . 92 meV for the out-of-plane component ( (cid:96) = z );and (cid:15) ∞ ,x = 4 . s x = 232 meV, ω x = 170 meV, and γ x = 3 . s x = 43 . ω x = 95 . γ x = 3 . (cid:96) = x ). Thegraphene conductivities σ j ( k (cid:107) , ω ) are calculated in theRPA [81–83], including nonlocal effects through their k (cid:107) dependence. We assume a graphene plasmon lifetime τ = 66 fs ( i.e. , ¯ hτ − = 10 meV) in all cases. Appendix D: Fresnel’s Reflection Coefficient Within the quasistatic approximation, the Fresnel re-flection coefficient of the graphene-hBN-graphene sand-wich vanishes for s polarization, whereas it can be readilyextracted from Eq. (C1) for p polarization. Indeed, con-sidering a distant point source in the z < ∝ (2 π/k (cid:107) )e − k (cid:107) z and ∝ − r p (2 π/k (cid:107) )e k (cid:107) z , whichallow us to identify r p = − A + A (cid:48) ( A + B (cid:48) )e − qd − A (cid:48) A (cid:48) e − qd (D1)by comparison to Eq. (C1) for light incident on layer 1. Appendix E: Effective quantum-well description ofthe out-of-plane graphene electron wave function We describe the out-of-plane electron wave functionof conduction electrons in graphene as a quantum wellstate. For simplicity, we assume a square well potentialof depth V and width a . We focus on the first excitedstate, which is antisymmetric along the normal direction z , just like the p z orbital in graphene, as shown in thefollowing scheme:We are therefore interested in the red wave function,which we denote ϕ ⊥ ( z ). In this model, the parameters V and a are used to fit (i) the state binding energy tothe graphene work function [89] Φ = 4 . (ii) thecentroid of the electron density outside the z = 0 planeto the value obtained for the carbon 2 p orbital.We find the quantum well state ϕ ⊥ ( z ) by solving theSchr¨odinger equation ( − ¯ h / m ) d ϕ ⊥ ( z ) /dz + [ V ( z ) − E ] ϕ ⊥ ( z ) = 0, where V ( z ) is the potential shown inthe scheme above. The wave function of the first ex-cited state can be written A sin( k in z ) inside the well,with k in = (cid:112) m ( E + V ) / ¯ h , while it decays evanes-cently in the outer region as sign( z ) B e − k out | z | , with k out = √− mE/ ¯ h . The electron energy E < z = ± a/ − k in /k out = tan( k in a/ E of asymmetricstates. Combining these conditions with normalization( (cid:82) dz | ϕ ⊥ ( z ) | = 1), we find | A | = ( a/ /k out ) − and | B | = | A | e k out a sin ( k in a/ z = 0 plane is calculated as (cid:82) dz | z | | ϕ ⊥ ( z ) | and com- pared with the centroid of the p z orbital of graphene (cid:82) d r | z | | ϕ p z ( r ) | . We approximate the latter by us-ing a tabulated 2 p atomic carbon wave function [90], ϕ p ( r ) = z (cid:80) j β j e − α j r , where the parameters α j and β j are expressed in the following table in atomic units: j α j β j p electron proba-bility density integrated over parallel ( x, y ) directions( (cid:82) dx (cid:82) dy | ϕ p ( r ) | , solid curve), compared with thefitted well state ( | ϕ ⊥ ( z ) | , dashed curve): P r obab ili t y den s i t y well2p orbital The agreement between the two probability densities isexcellent using fitted values V = 45 eV and a = 0 . 12 nm. Appendix F: Derivation of Eq. (A2) In this section, we start from Eq. (F1) for the inelasticelectron transition probability,Γ( ω ) = 2 e ¯ h (cid:88) i,f (cid:90) d r (cid:90) d r (cid:48) ψ † i ( r ) · ψ f ( r ) ψ † f ( r (cid:48) ) · ψ i ( r (cid:48) ) Im {− W ( r , r (cid:48) , ω ) } (F1) × δ ( ε f − ε i + ω ) f (¯ hε i ) [1 − f (¯ hε f )](see definitions of different elements in the Appendix),and specify it for the sandwich structure depicted inthe following scheme, consisting of two graphene layers separated by a film of thickness d and permittivity (cid:15) :0We factorize the electron wave functions as the productof in-plane ϕ (cid:107) and out-of-plane ϕ ⊥ components. Thelatter is described in Sec. E and is common to all con-duction electrons, so we can write ψ i ( r ) = ϕ (cid:107) i ( R ) ϕ ⊥ ( z )for initial states in layer 1, centered at z = 0, and ψ f ( r ) = ϕ (cid:107) f ( R ) ϕ ⊥ ( z − d ) for final states in layer 2, cen-tered at z = d . Here, we use the notation r = ( R , z ),with R = ( x, y ).The in-plane wave functions are Dirac fermions charac-terized by their parallel wave vector Q = ( Q x , Q y ), spin,and valley (K or K’ points). The transition probabilitymust be independent of spin and valley, so we performthe calculation for only one combination of these degreesof freedom and multiply the result by a factor of 4. Wefurther consider negative doping and a bias such thatthe initial and final wave functions lie within the upper(conduction) band of their respective layers. The Diracfermions admit the expression [80] ϕ (cid:107) ( R ) = 1 √ A e i Q · R (cid:18) e i φ Q / e − i φ Q / (cid:19) , (F2) where A is the normalization area, φ Q = tan − ( Q y /Q x )is the azimuthal angle of Q , and the upper and lowercomponents refer to the value of the wave function ineach of the two graphene carbon sublattices. The cor-responding electron energy relative to the Dirac point is¯ hε Q = ¯ hv F Q .We have adapted Eq. (F1) from a previous derivation[48] in order to incorporate the sum over both carbonsublattices, which is accounted for through the indicatedspinor products. Before inserting Eq. (F2) into Eq. (F1),we recast the sum over i as (cid:80) i → ( A/ π ) (cid:82) d Q i , andsimilarly for the sum over f . Additionally, as a con-sequence of translational invariance, the integrand in-side (cid:82) d r should be independent of R , so we can re-place (cid:82) d R → A . Now, we express the screened in-teraction as W ( r , r (cid:48) , ω ) = (2 π ) − (cid:82) d k (cid:107) exp[i k (cid:107) · ( R − R (cid:48) )] W ( k (cid:107) , z, z (cid:48) , ω ) (see Sec. G), which allows us to carryout the integral over R (cid:48) analytically to yield a δ functionfor conservation of parallel momentum, δ ( Q f − Q i + k (cid:107) ),and this in turn can be used to perform the integral over Q f . Putting these elements together, Eq. (F1) becomesΓ( ω ) = e A π ¯ h (cid:90) d k (cid:107) (cid:90) d Q i (cid:12)(cid:12)(cid:12)(cid:12) (e − i φ Q f / e i φ Q f / ) · (cid:18) e i φ Q i / e − i φ Q i / (cid:19)(cid:12)(cid:12)(cid:12)(cid:12) f (¯ hv F Q i ) [1 − f (¯ hv F | Q i − k (cid:107) | )] × (cid:90) dz (cid:90) dz (cid:48) ϕ ⊥ ( z ) ϕ ⊥ ( z − d ) ϕ ⊥ ( z (cid:48) ) ϕ ⊥ ( z (cid:48) − d ) Im (cid:8) − W ( k (cid:107) , z, z (cid:48) , ω ) (cid:9) × δ ( v F ( | Q i − k (cid:107) | − Q i ) + ω − eV b / ¯ h ) . (F3)Notice that the Fermi-Dirac distributions are referredto the Dirac point of their respective graphene lay-ers. However, the electron energy in layer 2 isshifted by the bias energy − eV b relative to layer 1(last term inside the δ function). We also note thatthe spinor product yields (cid:12)(cid:12) φ Q i − φ Q f )] (cid:12)(cid:12) =2 (cid:2) Q i · ( Q i − k (cid:107) ) Q − i | Q i − k (cid:107) | − (cid:3) , where we have ex-pressed the angle between Q i and Q f = Q i − k (cid:107) interms of the inner product of these two vectors. Fi-nally, inserting this expression into Eq. (F3), and notic-ing that the result is independent of the direction of k (cid:107) once the Q i integral has been carried out, we can makethe substitution (cid:82) dφ k (cid:107) → π for the azimuthal integraland divide the result by the graphene area A to read-ily obtain Eq. (2) of the main text. Equation 2 is theexpression that we use in our numerical simulations of the tunneling current, in which the azimuthal φ Q i in-tegral is carried out analytically by using the relation δ [ F ( φ Q i )] = (cid:80) j δ ( φ Q i − q j ) / | F (cid:48) ( q j ) | for the δ function(notice that the poles of F ( φ Q i ) are of first order). Appendix G: Derivation of Eq. (C1) The screened interaction W ( r , r (cid:48) , ω ) is defined as thescalar potential produced at r by a charge placed at r (cid:48) and oscillating with frequency ω . Translational invari-ance allows us to write W ( r , r (cid:48) , ω ) = (cid:90) d k (cid:107) (2 π ) e i k (cid:107) · ( R − R (cid:48) ) W ( k (cid:107) , z, z (cid:48) , ω ) , k (cid:107) space ( i.e. , we assume anoverall e i k (cid:107) · ( R − R (cid:48) ) dependence). A point charge placedat z (cid:48) produces a direct scalar potential (2 π/k (cid:107) )e − k (cid:107) | z − z (cid:48) | in vacuum ( i.e. , this is the direct Coulomb interactionterm in Eq. (3)). Additionally, inside the bulk of ananisotropic dielectric (permittivity (cid:15) z along z , and (cid:15) x along x and y ), the Poisson equation ∇ · (cid:15) ∇ φ = 0 hassolutions φ = e ± i qz , where q = k (cid:107) (cid:112) (cid:15) x /(cid:15) z and we take thesquare root to yield Im { q } > 0; this allows us to write thepoint-charge potential as [2 π/ ( (cid:15) z q )]e − q | z − z (cid:48) | inside thatmedium. Now, the induced potential has the form A e k (cid:107) z below the sandwich ( z < D e k (cid:107) ( d − z ) above it ( z > d ),and B e − qz + C e − q ( d − z ) inside the dielectric (0 < z < d ).Here, the coefficients A , B , C , and D are used to satisfythe boundary conditions, namely: (1) the continuity ofthe potential at each graphene layer j = 1 , 2; and (2)the jump of normal displacement is equal to 4 π timesthe induced charge. From the continuity equation, theinduced charge can be expressed as the divergence of thecurrent, and this in turn as the product of the conductiv-ity σ j times the in-plane electric field. The jump of nor-mal displacement at layer j is then given by − π i k (cid:107) σ j /ω times the potential. Solving the resulting system of fourequations for each position of the external charge z (cid:48) , weobtain, after some tedious but straightforward algebra,the expression for the screened interaction W ( k (cid:107) , z, z (cid:48) , ω )presented in Eq. (3) of the main text. Alternatively, a more direct Fabry-Perot-like derivation can be made interms of the transmission and reflection coefficients ofthe dielectric/graphene/vacuum interface defined in themain text. Appendix H: Additional numerical simulations We present additional simulations of the dispersion di-agrams and tunneling currents for various combinationsof the graphene Fermi energies in Figs. 5-9. We also showin Figs. 10 and 11 calculations similar to those of the Figs.1-4 for graphene layers separated by vacuum instead ofhBN. We assume a conservative graphene plasmon life-time of 66 fs in all cases. Acknowledgments This work has been supported in part by theSpanish MINECO (MAT2014-59096-P and SEV2015-0522), the European Commission (Graphene Flag-ship CNECT-ICT-604391 and FP7-ICT-2013-613024-GRASP), Ag`encia de Gesti´o d’Ajuts Universitaris ide Recerca (AGAUR) (2014-SGR-1400), and Fundaci´oPrivada Cellex. SdV acknowledges financial supportthrough the FPU program from the Spanish MECD. [1] K. R. Li, M. I. Stockman, and D. J. Bergman, Phys. Rev.Lett. , 227402 (2003).[2] R. A. ´Alvarez-Puebla, L. M. Liz-Marz´an, and F. J. Garc´ıade Abajo, J. Phys. Chem. Lett. , 2428 (2010).[3] B. Liedberg, C. Nylander, and I. Lunstr¨om, Sens. Actu-ators , 299 (1983).[4] J. N. Anker, W. P. Hall, O. Lyandres, N. C. Shah,J. Zhao, and R. P. Van Duyne, Nat. Mater. , 442 (2008).[5] S. Zeng, D. Baillargeat, H.-P. Hod, and K.-T. Yong,Chem. Soc. Rev. , 3426 (2014).[6] J. Reguera, D. Jimenez de Aberasturi, J. Langer, andL. M. Liz-Marz´an, Chem. Soc. Rev. in press , DOI:10.1039/c7cs00158d (2017).[7] D. P. O’Neal, L. R. Hirsch, N. J. Halas, J. D. Payne, andJ. L. West, Cancer Lett. , 171 (2004).[8] P. K. Jain, I. H. El-Sayed, and M. A. El-Sayed, NanoToday , 18 (2007).[9] Y. L. Luo, Y. S. Shiao, and Y. F. Huang, ACS Nano ,7796 (2011).[10] S. Palomba and L. Novotny, Phys. Rev. Lett. , 056802(2008).[11] S. Palomba, M. Danckwerts, and L. Novotny, J. Opt. A:Pure Appl. Opt. , 114030 (2009).[12] J. Shen, Y. Zhu, X. Yang, and C. Li, Chem. Commun. , 3686 (2012).[13] C. Wang, O. Ranasingha, S. Natesakhawat, P. R. Ohod-nicki, M. Andio, J. P. Lewis, and C. Matranga, Nanoscale , 6968 (2013).[14] C. Clavero, Nat. Photon. , 95 (2014). [15] J. Y. Park, L. R. Baker, and G. A. Somorjai, Chem. Rev. , 2781 (2015).[16] K. R. Catchpole and A. Polman, Opt. Express , 21793(2008).[17] H. A. Atwater and A. Polman, Nat. Mater. , 205 (2010).[18] Z. Fei, G. O. Andreev, W. Bao, L. M. Zhang,A. S. McLeod, C. Wang, M. K. Stewart, Z. Zhao,G. Dominguez, M. Thiemens, et al., Nano Lett. , 4701(2011).[19] Z. Fei, A. S. Rodin, G. O. Andreev, W. Bao,A. S. McLeod, M. Wagner, L. M. Zhang, Z. Zhao,M. Thiemens, G. Dominguez, et al., Nature , 82(2012).[20] J. Chen, M. Badioli, P. Alonso-Gonz´alez, S. Thongrat-tanasiri, F. Huth, J. Osmond, M. Spasenovi´c, A. Cen-teno, A. Pesquera, P. Godignon, et al., Nature , 77(2012).[21] Z. Fang, S. Thongrattanasiri, A. Schlather, Z. Liu, L. Ma,Y. Wang, P. M. Ajayan, P. Nordlander, N. J. Halas, andF. J. Garc´ıa de Abajo, ACS Nano , 2388 (2013).[22] V. W. Brar, M. S. Jang, M. Sherrott, J. J. Lopez, andH. A. Atwater, Nano Lett. , 2541 (2013).[23] V. W. Brar, M. Seok Jang, M. C. Sherrott, S. Kim, J. J.Lopez, L. B. Kim, M. Choi, and H. A. Atwater, NanoLett. , 3876 (2014).[24] H. Yan, Z. Li, X. Li, W. Zhu, P. Avouris, and F. Xia,Nano Lett. , 3766 (2012).[25] N. Kumada, P. Roulleau, B. Roche, M. Hashisaka, H. Hi-bino, I. Petkovi´c, and D. C. Glattli, Phys. Rev. Lett. , (a) (b)(c) (d)(e) (f) FIG. 5: Additional dispersion diagrams. Same as Fig. 2afor various combinations of the gap distance and the grapheneFermi energies (see labels).266601 (2014).[26] F. J. Garc´ıa de Abajo, ACS Photon. , 135 (2014).[27] G. X. Ni, L. Wang, M. D. Goldflam, M. Wagner, Z. Fei,A. S. McLeod, M. K. Liu, F. Keilmann, B. ¨Ozyilmaz,A. H. C. Neto, et al., Nat. Photon. , 244 (2016).[28] A. Woessner, M. B. Lundeberg, Y. Gao, A. Prin-cipi, P. Alonso-Gonz´alez, M. Carrega, K. Watanabe,T. Taniguchi, G. Vignale, M. Polini, et al., Nat. Mater. , 421 (2015).[29] D. N. Basov, M. M. Fogler, and F. J. Garc´ıa de Abajo,Science , aag1992 (2016).[30] H. Yan, T. Low, W. Zhu, Y. Wu, M. Freitag, X. Li,F. Guinea, P. Avouris, and F. Xia, Nat. Photon. , 394(2013).[31] M. Freitag, H.-Y. Chiu, M. Steiner, V. Perebeinos, andP. Avouris, Nat. Nanotech. , 497 (2010).[32] O. Ilic, M. Jablan, J. D. Joannopoulos, I. Celanovic,H. Buljan, and M. Soljaˇci´c, Phys. Rev. B , 155422(2012).[33] A. Manjavacas, S. Thongrattanasiri, J.-J. Greffet, andF. J. Garc´ıa de Abajo, Appl. Phys. Lett. , 211102(2014).[34] V. W. Brar, M. C. Sherrott, M. S. Jang, S. Kim, L. Kim,M. Choi, L. A. Sweatlock, and H. A. Atwater, Nat. Com-mun. , 7032 (2015).[35] R. Yu and F. J. Garc´ıa de Abajo, ACS Nano , 8045 (a) (b)(c) (d) FIG. 6: Additional calculations of the energy- andmomentum-resolved tunneling current. Same as Fig.2b for fixed gap distance d = 2 nm and various combinationsof the graphene Fermi energies (see labels).(2016).[36] M. B. Lundeberg, Y. Gao, A. Woessner, C. Tan,P. Alonso-Gonz´alez, K. Watanabe, T. Taniguchi,J. Hone, R. Hillenbrand, and F. H. L. Koppens, Nat.Mater. , 204 (2016).[37] Y. Li, H. Yan, D. B. Farmer, X. Meng, W. Zhu, R. M.Osgood, T. F. Heinz, and P. Avouris, Nano Lett. ,1573 (2014).[38] A. Marini, I. Silveiro, and F. J. Garc´ıa de Abajo, ACSPhoton. , 876 (2015).[39] D. Rodrigo, O. Limaj, D. Janner, D. Etezadi, F. J. Garc´ıade Abajo, V. Pruneri, and H. Altug, Science , 165(2015).[40] D. B. Farmer, P. Avouris, Y. Li, T. F. Heinz, and S.-J.Han, ACS Photon. , 553 (2016).[41] C. F. Chen, C. H. Park, B. W. Boudouris, J. Horng,B. Geng, C. Girit, A. Zettl, M. F. Crommie, R. A. Segal-man, S. G. Louie, et al., Nature , 617 (2011).[42] R. Ye, C. Xiang, J. Lin, Z. Peng, K. Huang, Z. Yan, N. P.Cook, E. L. Samuel, C.-C. Hwang, G. Ruan, et al., Nat.Commun. , 2943 (2013).[43] K. M¨ullen, ACS Nano , 6531 (2014).[44] Z. Wang, T. Li, K. Almdal, N. A. Mortensen, S. Xiao,and S. Ndoni, Opt. Lett. , 5345 (2016).[45] A. Manjavacas, F. Marchesin, S. Thongrattanasiri,P. Koval, P. Nordlander, D. S´anchez-Portal, and F. J.Garc´ıa de Abajo, ACS Nano , 3635 (2013).[46] A. Lauchner, A. Schlather, A. Manjavacas, Y. Cui, M. J.McClain, G. J. Stec, F. J. Garc´ıa de Abajo, P. Nordlan-der, and N. J. Halas, Nano Lett. , 6208 (2015).[47] Z. Fei, G.-X. Ni, B.-Y. Jiang, M. M. Fogler, andD. N. Basov, ACS Photon. p. DOI: 10.1021/acsphoton-ics.7b00477 (2017). (a) (b) (c) FIG. 7: Additional dispersion current diagrams. Same as Fig. 3 with | J ( k (cid:107) , ω ) | plotted in logarithmic scale. Energy (eV) -11 -10 -9 -8 -7 -6 ! ( " ) ( / n m ! ) Bias (eV) Ef1 = 0.75 eVEf2 = 0.5 eVd = 1 nm Energy (eV) -11 -10 -9 -8 -7 -6 ! ( " ) ( / n m ! ) Bias (eV) Ef1 = 0.75 eVEf2 = 0.25 eVd = 1 nm Energy (eV) -11 -10 -9 -8 -7 -6 -5 ! ( " ) ( / n m ! ) Bias (eV) Ef1 = 1.0 eVEf2 = 0.25 eVd = 1 nm Energy (eV) -11 -10 -9 -8 -7 -6 -5 ! ( " ) ( / n m ! ) Bias (eV) Ef1 = 1.0 eVEf2 = 0.5 eVd = 1 nm (a) (b)(c) (d) –7 –6 –10 –10 FIG. 8: Additional calculations of the spectrally re-solved tunneling current. Same as Fig. 4a for fixed gapdistance d = 1 nm and various combinations of the grapheneFermi energies (see labels).[48] F. J. Garc´ıa de Abajo, Rev. Mod. Phys. , 209 (2010).[49] D. C. Tsui, Phys. Rev. Lett. , 293 (1969).[50] D. C. Tsui and A. S. Barker Jr, Phys. Rev. , 590(1969).[51] C. B. Duke, Phys. Rev. , 588 (1969).[52] C. B. Duke, M. J. Rice, and F. Steinrisser, Phys. Rev. , 733 (1969).[53] E. N. Economou and K. L. Ngai, Phys. Rev. B , 4105(1971).[54] P. Johansson, R. Monreal, and P. Apell, Phys. Rev. Lett. , 9210 (1990).[55] R. Berndt, J. K. Gimzewski, and P. Johansson, Phys.Rev. Lett. , 3796 (1991).[56] P. Rai, N. Hartmann, J. Berthelot, J. Arocas, G. Colas Bias (eV) e ! / ( s ! n m " ) Ef1(eV) Ef2(eV) d = 1 nm E l e c t r on c u rr en t ( s –1 n m –2 ) FIG. 9: Additional calculations of the tunneling cur-rent. Same as Fig. 4b for fixed gap distance d = 1 nm andvarious combinations of the graphene Fermi energies (see la-bels).des Francs, A. Hartschuh, and A. Bouhelier, Phys. Rev.Lett. , 026804 (2013).[57] K. J. A. Ooi, H. S. Chu, C. Y. Hsieh, D. T. H. Tan, andL. K. Ang, Phys. Rev. Applied , 054001 (2015).[58] N. L. Schneider, G. Schull, and R. Berndt, Phys. Rev.Lett. , 026601 (2010).[59] P. Bharadwaj, A. Bouhelier, and L. Novotny, Phys. Rev.Lett. , 226802 (2011).[60] C. Gro β e, A. Kabakchiev, T. Lutz, R. Froidevaux,F. Schramm, M. Ruben, M. Etzkorn, U. Schlickum,K. Kuhnke, , et al., Nano Lett. , 5693 (2014).[61] J. Kern, R. Kullock, J. Prangsma, M. Emmerling,M. Kamp, and B. Hecht, Nat. Photon. , 582 (2015).[62] A. V. Uskov, J. B. Khurgin, M. Buret, A. Bouhelier, (nm –2 ) -11 -8 (a) (b) (c) o p t i c a l p l a s m o n BC a c o u s t i c p l a s m o n FIG. 10: Energy- and momentum-resolved electron tunneling with a vacuum gap. Same as Fig. 2 with the hBNfilm replaced by vacuum. (cid:15514) ! (eV) -14 -12 -10 -8 Bias (eV) e ! / ( s ! n m " ) (a) (b) ( n m –2 ) (eV) B C E l e c t r on c u rr en t ( s –1 n m –2 ) (nm) scaling FIG. 11: Probability of plasmon-generation by electron tunneling with a vacuum gap. Same as Fig. 4a,b with thehBN film replaced by vacuum.I. V. Smetanin, and I. E. Protsenko, ACS Photon. ,1501 (2017).[63] D. M. Koller, A. Hohenau, H. Ditlbacher, N. Galler,F. Reil, F. R. Aussenegg, A. Leitner, E. J. W. List, andJ. R. Krenn, Nat. Photon. , 684 (2008).[64] R. J. Walters, R. V. A. van Loon, I. Brunets, J. Schmitz,and A. Polman, Nat. Mater. , 21 (2010).[65] V. Enaldiev, A. Bylinkin, and D. Svintsov, 0 ,arXiv:1706.05216v1 (2017).[66] D. Svintsov, Z. Devizorova, T. Otsuji, and V. Ryzhii,Phys. Rev. B , 115301 (2016).[67] D. Yadav, S. B. Tombet, T. Watanabe, S. Arnold,V. Ryzhii, and T. Otsuji, 2D Mater. , 045009 (2016).[68] X. Wang, K. Braun, D. Zhang, H. Peisert, H. Adler,T. Chass´e, and A. J. Meixner, ACS Nano , 8176 (2015).[69] H. D. F. R. Aussenegg, J. R. Krenn, L. G. Jakopic, andG. Leising, Appl. Phys. Lett. , 161101 (2006). [70] R. W. Heeres, S. N. Dorenbos, B. Koene, G. S. Solomon,L. P. Kouwenhoven, and V. Zwiller, Nano Lett. , 661(2009).[71] P. Neutens, P. Van Dorpe, I. De Vlaminck, L. Lagae, andG. Borghs, Nat. Photon. , 283 (2009).[72] A. L. Falk, F. H. Koppens, L. Y. Chun, K. Kang,N. de Leon Snapp, A. V. Akimov, M.-H. Jo, M. D. Lukin,and H. Park, Nature Physics , 475 (2009).[73] T. Dufaux, J. Dorfm¨uller, R. Vogelgesang, M. Burghard,and K. Kern, Appl. Phys. Lett. , 161110 (2010).[74] I. Goykhman, B. Desiatov, J. Khurgin, J. Shappir, andU. Levy, Nano Lett. , 2219 (2011).[75] M. W. Knight, H. Sobhani, P. Nordlander, and N. J.Halas, Science , 702 (2011).[76] K. M. Goodfellow, C. Chakraborty, R. Beams,L. Novotny, and A. N. Vamivakas, Nano Lett. , 5477(2015). [77] L. Britnell, R. V. Gorbachev, A. K. Geim, L. A. Pono-marenko, A. Mishchenko, M. T. Greenaway, T. M.Fromhold, K. S. Novoselov, and L. Eaves, Nat. Commun. , 1794 (2013).[78] H. Jeong, H. M. Oh, S. Bang, H. J. Jeong, S.-J. An, G. H.Han, H. Kim, S. J. Yun, K. K. Kim, J. C. Park, et al.,Nano Lett. , 1858 (2016).[79] V. L. Katkov and V. A. Osipov, J. Vac. Sci. Technol. B , 050801 (2017).[80] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S.Novoselov, and A. K. Geim, Rev. Mod. Phys. , 109(2009).[81] P. R. Wallace, Phys. Rev. , 622 (1947).[82] B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New J.Phys. , 318 (2006).[83] E. H. Hwang and S. Das Sarma, Phys. Rev. B , 205418(2007).[84] R. B. Pettit, J. Silcox, and R. Vincent, Phys. Rev. B ,3116 (1975).[85] T. Ihn, J. G¨uttinger, F. Molitor, S. Schnez, E. Schurten- berger, A. Jacobsen, S. Hellm¨uller, T. Frey, S. Dr¨oscher,C. Stampfer, et al., Mater. Today , 44 (2010).[86] T. Low, R. Rold´an, H. Wang, F. Xia, P. Avouris,L. Mart´ın Moreno, and F. Guinea, Phys. Rev. Lett. ,106802 (2014).[87] M. A. Huber, F. Mooshammer, M. Plankl, L. Viti,F. Sandner, L. Z. Kastner, T. Frank, J. Fabian, M. S.Vitiello, T. L. Cocker, et al., Nat. Nanotech. , 207(2017).[88] G.-B. Liu, W.-Y. Shan, Y. Yao, W. Yao, and D. Xiao,Phys. Rev. B , 085433 (2013).[89] Y.-J. Yu, Y. Zhao, S. Ryu, L. E. Brus, K. S. Kim, andP. Kim, Nano Lett. , 3430 (2009).[90] E. Clementi and C. Roetti, At. Data Nucl. Data Tables , 177 (1974).[91] K. Watanabe, T. Taniguchi, and H. Kanda, Nat. Mater. , 404 (2004).[92] R. Geick, C. H. Perry, and G. Rupprecht, Phys. Rev.146