Plasmon generation through electron tunneling in twisted double-layer graphene and metal-insulator-graphene systems
Sandra de Vega F. Javier García de Abajo




ICFO-Institut de Ciencies Fotoniques, The Barcelona Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain

javier.garciadeabajo@nanophotonics.es
ICFO-Institut de Ciencies Fotoniques,  The Barcelona Institute of Science and Technology, 08860 Castelldefels (Barcelona), Spain
ICREA-Institució Catalana de Recerca i Estudis Avançats, Passeig Lluís Companys, 23, 08010 Barcelona, Spain




The generation of highly-confined plasmons through far-field optical illumination appears to be impractical for technological applications due to their large energy-momentum mismatch with external light. Electrical generation of plasmons offers a possible solution to this problem, although its performance depends on a careful choice of material and geometrical parameters. Here we theoretically investigate graphene-based structures and show in particular the very different performance between (i) two layers of graphene separated by a dielectric and (ii) metal|insulator|graphene sandwiches as generators of propagating plasmons assisted by inelastic electron tunneling. For double-layer graphene, we study the dependence on the relative tilt angle between the two sheets and show that the plasmon generation efficiency for 4^∘ twist angle drops to ∼20% of its maximum for perfect stacking. For metal|insulator|graphene sandwiches, the inelastic tunneling efficiency drops by several orders of magnitude relative to double-layer graphene, regardless of doping level, metal|graphene separation, choice of metal, and direction of tunneling (metal to or from graphene), a result that we attribute to the small fraction of the surface-projected metal Brillouin zone covered by the graphene Dirac cone. Our results reveal a reasonable tolerance to graphene lattice misalignment and a poor performance of structures involving metals, thus supporting the use of double-layer graphene as an optimum choice for electrical plasmon generation in tunneling devices.



Plasmon generation through electron tunneling in twisted double-layer graphene and metal-insulator-graphene systems
    F. Javier García de Abajo
    December 30, 2023
===================================================================================================================











§ INTRODUCTION


Surface plasmons –the collective excitations of conduction electrons at the surface of conducting materials– have been the focus of intense research over the past decades because of a their ability to interact strongly with light, producing huge concentration of electromagnetic energy down to deep-subwavelength regions, accompanied by strong enhancement of the optical field intensity.<cit.> Plasmons are relatively tolerant to defects in the fabrication of the structures supporting them, which are commonly relying on colloid chemistry<cit.> and lithography<cit.> methods, with application to optical sensing,<cit.> photocatalysis,<cit.> nonlinear optics,<cit.> and optical signal processing,<cit.> among other feats. However, the strong spatial confinement produced by optical plasmons has a negative side effect: the coupling cross section from far-field radiation to plasmons is small, typically comparable to or even smaller than the lateral size of the supporting plasmonic nanostructures. Plasmon excitation with light is thus inefficient, which represents a serious obstacle in the development of practical applications. Various methods have been devised to solve this problem, among them the use of funneling structures that focus light down to the spatial extension of the plasmon through the use of gratings<cit.> and tips.<cit.> Quantum emitters have also been employed as a source of plasmons,<cit.> and although they can function as a source of single plasmons, their efficiency is also low.

An alternative to the optical generation of plasmons is provided by electron beams, which have been used in pioneering studies of surface plasmons, including their discovery.<cit.> However, the production and control of electron beams is difficult to combine with integrated devices. Electrons can also inelastically undergo tunneling between neighboring conductors, losing energy that results in the emission of light.<cit.> Recently, plasmon emission produced by inelastic tunneling has been explored in different systems,<cit.> and this mechanism has also been theoretically explored in sandwiches formed between graphene layers,<cit.> where it is found to be particularly efficient.<cit.>

In this article, we extend a previous study<cit.> and present a comprehensive study of the plasmon-emission efficiency associated with inelastic tunneling between two graphene layers and compare the results with the emission in metal|insulator|graphene (MIG) configurations. For double-layer graphene (DLG), we assess the dependence of the emission efficiency on the relative graphene twisting angle  and find a reduction of the efficiency when the K points in the two layers are misaligned by more than a few degrees. For MIG configurations, we find plasmon generation rates way below that of DLG.




§ SPECTRALLY RESOLVED INELASTIC TUNNELING CURRENT



The large plasmon confinement and small thickness compared with the light wavelength in the structures under consideration allow us to work within the quasistatic limit. We thus calculate the probability for an electron to inelastically tunnel between initial and final states ψ_i and ψ_f, while transferring an amount of energy ħω to the materials (e.g., a plasmon), by using the frequency-resolved screened interaction W(,',ω), which describes the potential created at  by a point charge placed at ' and oscillating with frequency ω. More precisely, the probability can be written as<cit.>

	
    Γ(ω)= 42e^2ħ∑_i,f∫^3 ∫^3' ψ_i^†()·ψ_f()  ψ_f^†(') ·ψ_i(') Im{ - W(,',ω) }δ(ε_f - ε_i + ω)   f_i(ħε_i)   [1 - f_f(ħε_f)] ,


where the leading factor of 4 accounts for spin and valley degeneracies in graphene. This factor is the same for DLG and MIG structures, as they involve at least one graphene layer. Electrons undergo transitions from initial occupied states to final empty states of energies ħε_i and ħε_f, respectively. The occupation of these levels follows the Fermi-Dirac distributions f_i|f(ħε_i|f)=1/{1 + exp[(ħε_i|f-E_ F,i|f)/k_ BT]}, where E_ F,i and E_ F,f are the corresponding Fermi energies in the emitting and receiving materials, referenced to a common origin of energies. We work here in the T=0 limit, for which these distributions become step functions f_i(ħε_i)=Θ(E_ F,i-ħε_i) and 1-f_f(ħε_f)=Θ(ħε_f-E_ F,f). For initial or final states in graphene, we recast the sum over electron states as ∑_i|f→(2π)^-2A∫^2_i|f, while for metals it becomes ∑_i|f→(2π)^-3V∫^3 _i|f, where A and V are the corresponding normalization area and volume, respectively. In the evaluation of Eq. (<ref>), we separate the electron wave functions into parallel and perpendicular components as ψ()=φ^∥()φ^⊥(z), where =(x,y) are in-plane coordinates. Finally, we present some results below for the loss probability normalized per unit of film surface area J(ω)=Γ(ω)/A (i.e., the spectrally resolved inelastic tunneling current).






§ GRAPHENE|HBN|GRAPHENE STRUCTURES: DEPENDENCE ON DOPING LEVEL



We first discuss structures formed by two graphene layers separated by a hexagonal boron nitride (hBN) film, assuming perfect alignment between the two graphene reciprocal crystal lattices. We introduced this system in a previous publication,<cit.> but here we study the dependence on the doping level of the two graphene layers. The tunneling structure is sketched in Fig. <ref>(a): a sandwich of two perfectly stacked graphene layers, with different doping levels, conductivities σ_1 (bottom) and σ_2 (top), and connected with a bias voltage , are separated by a distance d≃ 1 nm (∼3 atomic layers of hBN), which is maintained throughout this work. The first Brillouin zone of the graphene reciprocal lattice is schematically plotted in Fig. <ref>(b). Around the K point [=(2π/3a) (1,1/√(3)), where a=1.42 Å  is the nearest-neighbor C-C distance], the electronic band structure is conical in  (the 2D wave vector relative to K) and we can express the electron wave functions as spinors with two components, each of them associated with one of the two carbon atoms in the unit cell.<cit.> This representation leads to a closed-form expression for the graphene surface conductivity in the random-phase approximation<cit.> (RPA), which we use throughout this work with a phenomenological inelastic electron lifetime of 66 fs (10 meV energy width, or equivalently, a moderate mean-free path of 66 nm for a Fermi velocity v_ F≈10^6m/s).

For the evaluation of Eq. (<ref>), we incorporate the following elements:


  * Screened interaction. We use an analytical expression for the screened interaction describing the graphene|hBN|graphene sandwich, given in Appendix <ref> as a generalization of the results presented in Ref. paper295 to an asymmetric environment. This expression incorporates the anisotropy of hBN through its parametrized dielectric function<cit.> for in- and out-of-plane directions, which accounts for optical phonons in this material.

  * Parallel electron wave functions. We use a Dirac-fermion description of the in-plane graphene electron wave functions as

    φ^∥_i|f() = 1√(2A) ^_i|f·^ϕ_i|f /2^-ϕ_i|f /2^·,

where ϕ_i|f is the azimuthal angle of _i|f and  is the wave vector at the K point relative to the Γ point [see Fig. <ref>(b)]. We disregard inter-valley scattering, so each of the two inequivalent K points in the first Brillouin zone produces an identical contribution.

  * Perpendicular electron wave functions. We describe the evolution of the electron along z through a 1D wave function trapped in the graphene layers by potential wells, which are vertically off-set due to . Full details of this wave function are given elsewhere for graphene,<cit.> while explicit expressions for metals are offered in Appendix <ref> by treating the surface in the one-electron step-potential approximation.


The hybridization of plasmons in the two neighboring layers produces a characteristic hybridization scheme, leading to two plasmon branches: optical and acoustic. Acoustic plasmons have lower energy and their out-of-plane electric-field profile looks antisymmetric. In contrast, optical plasmons possess higher energy and symmetric field profiles, so they are generally easier to excite and manipulate, and therefore, more suitable for photonic applications.

Tunneling requires ≠ 0 and ≠. We take >, and hence, favor tunneling from layer 1 (bottom) to layer 2 (top). We explore the effect of varying the doping-level difference between the two graphene layers by calculating the spectrally resolved inelastic tunneling current shown in Fig. <ref>(c), where we fix =0.3 eV and e=1.2 eV. Energy splitting due to phonons in the hBN film can be observed at low energy transfers (below 0.3 eV), losing strength as the Fermi energy difference Δ=- decreases.

Acoustic plasmons become apparent for Δ>0.2eV. These resonances have an average plasmon lifetime of ∼50-70 fs [Fig. <ref>(f), blue-dashed curve] and undergo a redshift with increasing Δ [Fig. <ref>(f), blue-solid curve]. In contrast, optical plasmons become more pronounced when > 2, they experience a milder redshift with increasing Δ [Fig. <ref>(f) green-solid curve], and their lifetime is rapidly increasing from ∼10 fs to ∼80 fs [Fig. <ref>(f) green-dashed curve].

We quantify the fraction of the total inelastic current invested into exciting plasmons by comparing the area under the whole spectrum in Fig. <ref>(g) for different values of  [Fig. <ref>(g), black curve] to the area under either the optical (green curve) or acoustic (blue curve) plasmon regions. The average generation efficiency for the selected values of  and  lies in the ∼20-60% range for optical plasmons and < 1% for acoustic modes.






§ GRAPHENE|HBN|GRAPHENE STRUCTURES: DEPENDENCE ON TWIST ANGLE


In practical devices, the alignment of the two graphene layers can be a challenge, so we examine the effect of a finite twist angle between their respective lattices. Rotations in real space result in rotations around the Γ point in momentum space. In Fig. <ref>(a) we depict two reciprocal lattices corresponding to the bottom (1, black) and top (2, red) layers for a finite rotation angle θ. The circles surrounding their respective K points represent the projections of their Fermi surfaces on 2D momentum space for Fermi energies =1 eV and =0.5 eV. Inelastic tunneling requires the two circles to be either overlapping or closer to each other than the plasmon momentum, which is small compared with the Fermi momenta (k_ F=/ħ v_ F, i.e., the radii of the plotted circles), leading to a cutoff angle ∼10^∘.
 
We introduce the twist angle θ in our formalism through the electron wave function of the final states, so that we maintain the initial states as in Eq. (<ref>) with =_i, but write the final wave function as

    φ^∥_f() = 1√(2A) ^_f·^ (ϕ_f - θ) /2^- (ϕ_f - θ) /2^_f·.

Notice that _i and _f are the K points corresponding to the centers of the circles in Fig. <ref>(a). The spinor product in Eq. (<ref>) involving the parallel wave functions thus becomes

    φ^∥_i^†()·φ^∥_f()  φ^∥_f^†(') ·φ^∥_i(') 
       = 12A^2 ^ (_f-_i+_f-_i)·( - ') [1 + cos(ϕ_i - ϕ_f + θ)],

with _f - _i = (2π/3a) (cosθ - sinθ/√(3) - 1, sinθ+cosθ/√(3)-1/√(3)). Introducing this expression together with the parallel-wave-vector decomposition of W [Eq. (<ref>)] into Eq. (<ref>) and performing the  and ' integrals, we find the condition +_f-_i+_f-_i=0, which guarantees momentum conservation. We find it convenient to separate the contribution of the perpendicular wave functions as

    I_1(,ω) = - ∫ z∫ z'    φ^⊥_i^†(z)φ^⊥_f(z)  φ^⊥_f^†(z')φ^⊥_i(z') 
       ×Im{  W(,z,z',ω)},

where φ^⊥_i|f are solutions corresponding to the two quantum wells that we use to model the graphene layers.<cit.> Putting these elements together, and considering the initial electron to tunnel from the conduction band of layer 1, we can write from Eq. (<ref>) the spectrally-resolved tunneling current density

    J(ω)   =Γ(ω)/A= e^24π^4∫^2_i ∫^2  I_1(,ω) 
       ×δ(ħ Q_f -ħ Q_i +-- e+ ħω) 
       × [1+cos(ϕ_i-ϕ'_f-θ)] 
       ×Θ( -ħ Q_i) Θ(ħ Q_f -),

where =_f-_i+, ϕ'_f=tan^-1[(Q_yi-_y)/(Q_xi-_x)], and Q_f=|_i-|. We note that the Dirac δ-function ensures energy conservation in Eq. (<ref>), limiting the spectral range to -+e/ħ<ω<+e/ħ. For θ=0 (no twist), this expression reduces to that of Ref. paper295, where incidentally, a spurious factor of 2π was accidentally introduced in the numerical results.

We represent in Fig. <ref>(b) the tunneling current obtained from Eq. (<ref>) for different rotation angles when fixing the bias voltage to =1.2 eV and the Fermi energies to =1 eV and =0.5 eV. When plotting the maximum of J(ω) associated with the optical plasmon (ħω_p^opt≃0.61 eV), we find a sharp decay with increasing twist angle [Fig. <ref>(c)], although we still maintain ∼20% of the maximum value for θ=4^∘, further indicating a reasonable tolerance of this type of device against unintended misalignments below the 1^∘ level.

So far, we have considered electron doping in both of the graphene layers. We can straightforwardly repeat the above analysis to find expressions that apply to situations in which one or both of the graphene layers is doped with holes. Furthermore, for sufficiently high , electrons can tunnel from the valence band of layer 1 even when both layers have electron doping, which leads to an additional term in the integrand of Eq. (<ref>); although we give the resulting expression below these lines, we do not consider such high voltages in this study. We find

    J(ω)= e^24π^4   ∫^2_i ∫^2 I_1(,ω) 
       × [1+cos(ϕ_i-ϕ'_f-θ)]  Δ(_i,,ω),

where


    Δ(_i,,ω) =
    {[      δ(ħ Q_f -ħ Q_i +-- e+ ħω)
                                               Θ( -ħ Q_i)  Θ(ħ Q_f -)                                                                                     (e-e doping);                                                          +δ(ħ Q_f +ħ Q_i +-- e+ ħω)  Θ(ħ Q_f -),                                                                                                 ;                                                                                                                                                                                                  ; δ(-ħ Q_f -ħ Q_i + + - e+ ħω)
                                               Θ( - ħ Q_i)  Θ( - ħ Q_f)                                                                                     (e-h doping);                         +δ(ħ Q_f  -ħ Q_i + + - e+ħω)
                                    Θ( - ħ Q_i)                                                                                                 ;                        +δ(-ħ Q_f +ħ Q_i + + - e+ ħω)
                                    Θ( - ħ Q_f)                                                                                                 ;                                                                    +δ(ħ Q_f  +ħ Q_i + + - e+ħω),                                                                                                 ;                                                                                                                                                                                                  ;        δ(ħ Q_f   +ħ Q_i - - - e+ ħω)
                                      Θ(ħ Q_i - )  Θ(ħ Q_f - ),                                                                                     (h-e doping);                                                                                                                                                                                                  ;                                          δ(-ħ Q_f   + Q_i - + - e+ ħω)  Θ(ħ Q_i - )  Θ( - ħ Q_f)                                                                                     (h-h doping);                                                     +δ(ħ Q_f   +ħ Q_i - + - e+ ħω)  Θ(ħ Q_i - ),                                                                                                  ].


the labels e (electron) and h (hole) indicate the type of doping in the first and second layers, respectively, and several δ-functions appear in each expression depending on whether the electron originates in the valence or conduction bands (for electron doping in layer 1) and whether the electron tunnels to the valence or conduction bands (for hole doping in layer 2). We define  and  as positive quantities, although it is understood that they represent a lowering of the Fermi energy relative to the Dirac point when doping with holes.

In Figs. <ref>(d,h-j) we show results for two perfectly stacked graphene layers with different types of doping. Arrows in the insets indicate the different tunneling channels. Additionally, we compare the results to a reference spectrally-resolved tunnelling probability associated with conduction-to-conduction transitions in two electron-doped graphene layers, which we conclude in fact to be the most effective mechanism to generate plasmons under the conditions of the figure.




§ PLASMON GENERATION IN METAL|INSULATOR|GRAPHENE TUNNELING STRUCTURES


We now compare the performance of DLG and MIG structures. In the latter, one of the graphene layers is substituted by a semi-inifinite metal medium. When the tunneling electrons go from the metal to the graphene sheet, not all metals can comply with momentum conservation. In Fig. <ref>(a) we represent the surface projection of the Fermi spheres (assuming for simplicity an independent-electron description of the metal Fermi sea) of different plasmonic metals (Au, Cu, and Al) compared with the first Brillouin zone of graphene, where we further show the graphene Fermi surface for a doping level =1 eV. When comparing the Fermi wave vectors of these metals (=12.1 nm^-1, =12.0 nm^-1, =13.6 nm^-1, =17.5 nm^-1) with the wave vector at the graphene K point (K=17.0 nm^-1), considering the small radius of the graphene Fermi circle (=/ħ=1.52 nm^-1), we find that only the Al Fermi sea overlaps the graphene Dirac cone, so this is the only one among the good plasmonic metals in which electrons can tunnel to graphene (without the mediation of phonons or defects). Then, it is reasonable to consider the structure represented in Fig. <ref>(b), consisting of a graphene layer on top of an aluminum surface coated with a 1 nm layer of oxide (Al|Al_2O_3|graphene). A gate voltage is then introduced to make metal electrons tunnel into the doped graphene sheet and excite plasmons in the MIG structure.

In contrast, when electrons tunnel from graphene to the metal, the situation is reversed, so it is favorable to have the occupied Dirac-cone region outside the projected metal Fermi sea. This situation is encountered with different choices of metal, and in particular with gold, separated from graphene by 3 atomic layers of hBN, as depicted in Fig. <ref>(c) (Au|hBN|graphene).

For simplicity, we calculate the screened interaction (see Appendix <ref>) assimilating the metal to a perfect electric conductor (|ϵ|→∞ is a good approximation for Au and Al within plasmonic energy range under consideration). We further approximate the response of aluminum oxide by using a constant isotropic permittivity ϵ_Al_2O_3=3 (the measured permittivity<cit.> only changes by ∼4% within the energy range under consideration). Finally, we use the RPA conductivity for graphene<cit.> (see Sec. <ref>) and the anisotropic permittivity described in 	Appendix <ref> for hBN.<cit.>

We calculate the spectrally-resolved tunneling probability using Eq. (<ref>) and the formalism described in Sec. <ref>, with the metal wave functions (either as initial or final states, depending on the configuration) described as

    ψ_i|f()= φ^⊥_i|f(z,k_i|f z) 1√(A) ^_i|f^∥·.
	
The perpendicular component of these wave functions depend on the incident wave vector along z (i.e., k_i|f z) according to the explicit expressions derived in Appendix <ref> under the assumption of a step potential to represent the metal|insulator interface.

As graphene has two sublattices, a small phase difference ^_i^∥· a has to be introduced to account in the coupling to the metal wave functions. Additionally, depending on the applied bias, the graphene doping level, and the type of doping (i.e., electrons or holes), only one or the two graphene Dirac cones can be engaged in the tunneling process. Taking these elements into account, the expression for the inelastic tunneling current calculated from Eq. (<ref>) needs to be specified for each of the following MIG configurations:



  * Al→Al_2O_3→electron-doped-graphene [see Fig. <ref>(a)]

    J(ω) = ħ e^24π^5   ∫^2_f ∫^3_i   I_2(,k_iz,ω)  Θ(  - ħ^2 k_i^2/2)[ 1 + cos(ϕ_f + _i^∥·) ]  
       ×δ(ħ Q_f - e -  - ħ^2 k_i^2/2 +    + ħω)  Θ(ħ Q_f - ),



  * Al→Al_2O_3→hole-doped-graphene [see Fig. <ref>(e)]

    J(ω) = ħ e^24π^5   ∫^2_f ∫^3_i   I_2(,k_iz,ω)  Θ(  - ħ^2 k_i^2/2)[ 1 + cos(ϕ_f + _i^∥·) ]  
       ×[ δ(-ħ Q_f - e +  - ħ^2 k_i^2/2 +    + ħω)  Θ (-ħ Q_f) 
       + δ(ħ Q_f - e +  - ħ^2 k_i^2/2 +    + ħω)],



  * electron-doped-graphene→hBN→Au [see Fig. <ref>(f)]

    J(ω) = ħ e^24π^5   ∫^2_i ∫^3_f   I_3(,k_fz,ω)  Θ( ħ^2 k_f^2/2 - )[ 1 + cos(ϕ_i + _f^∥·) ]  
       ×[ δ(-ħ Q_i - e +  + ħ^2 k_f^2/2 -    + ħω)  Θ (-ħ Q_i) 
       + δ(ħ Q_i - e +  + ħ^2 k_f^2/2 -    + ħω)],



  * hole-doped-graphene→hBN→Au [see Fig. <ref>(g)]

    J(ω) = ħ e^24π^5   ∫^2_i ∫^3_f   I_3(,k_fz,ω)  Θ( ħ^2 k_f^2/2 -)[ 1 + cos(ϕ_i + _f^∥·) ]  
       ×δ(ħ Q_i - e -   + ħ^2 k_f^2/2 -    + ħω)  Θ(ħ Q_i -  ),


where the arrows indicate the direction of electron tunneling, and we have defined

    I_2(,k_iz,ω) = -∫ z∫ z' φ^⊥_i^†(z,k_iz)φ^⊥_f(z)  φ^⊥_f^†(z')φ^⊥_i(z',k_iz) Im{W(,z,z',ω)}

for metal→graphene tunneling and

    I_3(,k_fz,ω) = -∫ z∫ z' φ^⊥_i^†(z)φ^⊥_f(z,k_fz)  φ^⊥_f^†(z',k_fz)φ^⊥_i(z') Im{W(,z,z',ω)}

for graphene→metal. We note that energy conservation can only be fulfilled for energies ħω<e in all of the these structures.


The contour plots in Fig. <ref>(d,e) portray Eqs. (<ref>) and (<ref>) for fixed values of k_iz and k_fz in the Al|Al_2O_3|graphene and Au|hBN|graphene structures, respectively, with a graphene Fermi energy fixed to =1 eV. They clearly reveal a plasmon mode arising from the poles of the screened potential W. Additionally, Fig. <ref>(e) shows low-energy features associated with hBN phonons in the Au|hBN|graphene structure. 



In order to quantify the amount of inelastic current associated with plasmon generation, we separate the screened interaction W=W^dir+W^ref into the sum of external (W^dir) and reflected (W^ref) components (see Appendix <ref>). Plasmons arise from the poles of W^ref, originating in a denominator of the form η()=1+4πσ/ω+ϵ̃_2+^2qd(ϵ̃_2-1-4π σ/ω), where ϵ̃_2=√(ϵ_2xϵ_2z), q=√(ϵ_2x/ϵ_2z), ϵ_2x and ϵ_2z are the permittivities of the insulator (Al_2O_3 or hBN) along in-plane (x) and out-of-plane (z) directions, and σ is the graphene conductivity. We now isolate this pole and write W^ref=η^-1W̃^ref. In the vicinity of the plasmon pole =k_ p, we can Taylor-expand η ()=η'(k_p)(-k_ p) to first order and approximate the plasmon contribution to Eqs. (<ref>) and (<ref>) as


    I_2^pl(,k_iz,ω)    = -∫ z∫ z' φ^⊥_i^†(z,k_iz)φ^⊥_f(z)  φ^⊥_f^†(z')φ^⊥_i(z',k_iz)Im{W̃(,z,z',ω)η'(k_p)( - k_p)}, 
    
    I_3^pl(,k_fz,ω)    = -∫ z∫ z' φ^⊥_i^†(z)φ^⊥_f(z,k_fz)  φ^⊥_f^†(z',k_fz)φ^⊥_i(z')Im{W̃(,z,z',ω)η'(k_p)( - k_p)}.


By using either Eqs. (<ref>)-(<ref>) or Eqs. (<ref>)-(<ref>) into Eqs. (<ref>)-(<ref>), we obtain the total spectrally-resolved inelastic tunneling current J(ω) or the contribution arising from plasmon generation J^pl(ω), respectively.



We show the calculated spectrally-resolved currents J(ω) and J^pl(ω) in Fig. <ref>(b,c) for Al→Al_2O_3→electron-doped-graphene inelastic tunneling with a bias voltage in the 0.1 eV <e<1.2 eV range and =1 eV. A significant spectral broadening in the total current [Fig. <ref>(b)] originates in the availability of multiple inelastic channels associated for example with electron-hole pair generation in the metal. Indeed, plasmons make a relatively moderate contribution  [see Fig. <ref>(c)], as we conclude by comparing the total to the plasmon-based ω-integrated inelastic currents [Fig. <ref>(d)]. The plasmon generation efficiency (number of plasmons per tunneling electron, obtained from the ratio J^pl/J) varies from <10^-6 for <0.3eV to 0.1 for voltages close to 1.2 eV, where the plasmon generation rate is ∼10^-2 nm^-2s^-1. Figure <ref>(e-g) represents the total and plasmon-based inelastic currents for the configurations depicted below the graphs, from which we conslude that Al→Al_2O_3→hole-doped-graphene exhibits better performance, with plasmon generation rates approaching ∼1nm^-2s^-1, while for Au|hBN|graphene structures the currents are 16 orders of magnitude smaller. Overall, comparing the performance of these MIG systems to DLG, we find the latter to be much more efficient, with inelastic currents reaching 10^8 nm^-2s^-1 (see Fig. <ref> and Ref. paper295).




§ CONCLUSIONS


Our simulations reveal the suitability of DLG heterostructures as plasmon sources even under moderately unfavorable conditions produced by a finite twist angle between the two graphene sheets, with >70% efficiency for a twist angle as large as 4^∘. Additionally, we find the relation =2 between the Fermi energies of the two layers to be an optimum choice to maximize the plasmon emission rate of both acoustic and optical plasmons, although the efficiency is still in the >10% range for order-one variations in the Fermi energies. Furthermore, our study of MIG structures leads to efficiencies that are orders of magnitude lower than DLG. 

From an intuitive viewpoint, these conclusions can be understood using the following argument. We are trying to project electrons from one electrode into the other, which requires conservation of parallel momentum, differing by just the plasmon momentum, which is small compared with the size of the Brillouin zone. For DLG, the matching is most efficient with perfect alignment, involving large overlaps of their respective Dirac cones, while the tolerance against twisting mentioned above can be roughly quantified by the angle required to produce total mismatch between the Dirac cones, of the order of a few degrees for doping levels of 0.5-1 eV. Unfortunately, when one of the electrodes is a metal, although the conduction electron density at the surface is 2-3 orders of magnitude larger than that of charge carriers in graphene, they are distributed over a larger momentum-space region, rendering the overlap with the Dirac cone smaller; this effect, together with a weaker plasmon strength in MIG compared with DLG, results in much poorer plasmon generation rates.

In conclusion, DLG offers an optimum choice for the electrical generation of plasmons based upon inelastic electron tunneling, which is moderately robust against twisting of the layers (i.e., misalignment of their respective K points in the first Brillouin zone) and takes place for broad ranges of doping (i.e., it is also tolerant with respect to variations in doping).






§ ELECTRON WAVE FUNCTIONS IN METAL|INSULATOR INTERFACES



We approximate the electron wave functions in the metal along the z-direction as those of a quantum step-potential, assuming the configuration shown in the following scheme (metal for z<0 and insulator for z>0):

 
    < g r a p h i c s >
 

The Fermi energies of the metals here considered are =11.7eV and =5.53eV. The parameter E_C accounts for the band gap of the insulators, for which we take E_C^Al_2O_3=3.5eV<cit.> and E_C^hBN=2.6eV (i.e., half of their gap energies, under the assumption that their Fermi levels are in the center of the gap, aligned with that of the metal). Moreover, T_1 and R_1 are the transmission and reflection coefficients at the interface, respectively. By solving the time independent Schrödinger equation (-ħ^2/2) d^2φ^⊥(z)/dz^2 + V(z) φ^⊥(z) = E φ^⊥(z) with V(z)=-V_0 θ(-z) and V_0=E_C+, an electron incident on the interface from the metal side has the wave function

    φ^⊥_i(z,k_iz)={[ 1√(2π)(^ k_iz z + R_1^- k_iz z),                              z<0;             1√(2π) T_1 ^- k_2 z,                              z>0;                                  ].

where k_iz is the electron wave vector along the interface normal in the metal side, k_2 = √(2 m V_0/ħ^2 - k_iz^2), R_1 = (k_iz -  k_2)/(k_iz +  k_2) and T_1=2k_iz / (k_iz + k_2).




§ SCREENED INTERACTION POTENTIAL



Because of the 2D translational invariance of the system, we can write the screened interaction as

    W(,',ω) = ∫^2(2π)^2 ^·(-')  W(,z,z',ω)
 
in terms of its parallel wave-vector components. We consider a multilayer structure with the permittivities ϵ_1|2|3 and surface conductivities σ_1|2 as defined by the following scheme:

 
    < g r a p h i c s >
 

This structure can describe the two configurations considered in this work: MIG (with on of the σ_j's set to zero) and DLG. An expression for W was presented in a previous work<cit.> for a symmetric environment (ϵ_1=ϵ_3=1), which we generalize here for the asymmetric configuration sketched above. Using the methods explained in that work, we separate

    W(,z,z',ω)=W^ dir(,z,z',ω)+W^ ref(,z,z',ω),

into direct and surface-reflection contributions, where


    W^ dir(,z,z',ω)=2π/×{[   (1/ϵ_1)^-|z-z'|,             z,z'<0;   (1/ϵ_3)^-|z-z'|,             z,z'>0; (1/ϵ̃_2)^-q|z-z'|,           0<z,z'<d;                 0,          otherwise ].

and

    W^ ref(,z,z',ω)=(2π/)/1-A'_1A'_2^-2qd
       ×{[ (1/ϵ_3)^(2d-z-z')[A_2+A'_1(A_2+B'_2) ^-2qd],                                         d<z,                                         d<z'; (1/ϵ_3)B_2 ^(d-z)[^-q(d-z')+A'_1 ^-q(d+z')],                                         d<z,                                       0<z'<d;            (ϵ̃_2/ϵ_1) B_1B_2 ^ (d-z+z')^-qd,                                         d<z,                                         z'<0;  (1/ϵ_3)B_2 ^(d-z')[^-q(d-z)+A'_1 ^-q(d+z)],                                       0<z<d,                                         d<z';   (1/ϵ̃_2){
    A'_1 ^-q(z+z')+A'_2 ^-q(2d-z-z')                                                                                          ;      +A'_1A'_2[^-q(2d+z-z')+^-q(2d-z+z')]
    },                                       0<z<d,                                       0<z'<d;        (1/ϵ_1)B_1 ^ z'[^-qz+A'_2 ^-q(2d-z)],                                       0<z<d,                                         z'<0;             (ϵ̃_2/ϵ_3) B_1B_2 ^(d+z-z')^-qd,                                         z<0,                                         d<z';      (1/ϵ_1)B_1 ^ z[^-qz'
    +A'_2 ^-q(2d-z')],                                         z<0,                                       0<z'<d;    (1/ϵ_1)^(z+z')[A_1+A'_2(A_1+B'_1) ^-2qd],                                         z<0,                                         z'<0 ].


This expression is applicable when the material in medium 2 is anisotropic (e.g., hBN) by defining ϵ̃_2=√(ϵ_2xϵ_2z) as the geometrical average of the in-plane (x) and out-of-plane (z) permittivities, with the square root chosen to yield a positive imaginary part. Also, q=√(ϵ_2x/ϵ_2z) (with Re{q}>0) is the effective out-of-plane wave vector in that medium, and we have defined the coefficients

    B_1   =2ϵ_1/(ϵ_1+ϵ̃_2+β_1), 
    
    B'_1   =(ϵ̃_2/ϵ_1) B_1, 
    
    B_2   =2ϵ_3/(ϵ_3+ϵ̃_2+β_2), 
    
    B'_2   =(ϵ̃_2/ϵ_3) B_2, 
    
    A_j   =B_j-1, 
    
    A'_j   =B'_j-1, 
    β_j   =4πσ_j/ω,

for j=1,2. In this work, we approximate the permittivity of alumina as a constant value ϵ̃_2=ϵ_Al_2O_3=3, we use the perfect-conductor limit for metals (|ϵ_1|→∞), we describe the graphene conductivity in the RPA,<cit.> and we use the expression <cit.>

    ϵ_2ℓ(ω)≈ϵ_∞,ℓ+∑_i=1,2s_ℓ i^2/[ω_ℓ i^2-ω(ω+γ_ℓ i)],

for the conductivity of hBN, where
ϵ_∞,z=4.10,
s_z 1=70.8meV,
ω_z 1=97.1meV,
γ_z 1=0.99meV,
s_z 2=126meV,
ω_z 2=187meV,
and
γ_z 2=9.92meV
for the z component (ℓ=z); and
ϵ_∞,x=4.95,
s_x 1=232meV,
ω_x 1=170meV,
and
γ_x 1=3.6meV,
s_x 2=43.5meV,
ω_x 2=95.1meV,
and
γ_x 2=3.4meV
for the x component (ℓ=x). The latter incorporates hBN phonons as Lorentzians terms.

This work has been supported in part by the Spanish MINECO (MAT2017-88492-R and SEV2015- 0522), the ERC (Advanced Grant 789104-eNANO), the European Commission (Graphene Flagship 696656), the Catalan CERCA Program, and the Fundació Privada Cellex. SdV acknowledges financial support through the FPU program from the Spanish MECD.




31
natexlab#1#1bibnamefont#1#1bibfnamefont#1#1citenamefont#1#1url<#>1urlprefixURL 


[Schuller et al.(2010)Schuller, Barnard,
  Cai, Jun, White, and Brongersma]SBC10
authorJ. A. Schuller,
  authorE. S. Barnard,
  authorW. Cai,
  authorY. C. Jun,
  authorJ. S. White, and
  authorM. L. Brongersma,
  journalNat. Mater. volume9,
  pages193 (year2010).

[Halas et al.(2011)Halas, Lal, Chang,
  Link, and Nordlander]HLC11
authorN. J. Halas,
  authorS. Lal,
  authorW. Chang,
  authorS. Link, and
  authorP. Nordlander,
  journalChem. Rev. volume111,
  pages3913 (year2011).

[Liz-Marzán(2006)]L06
authorL. M. Liz-Marzán,
  journalLangmuir volume22,
  pages32 (year2006).

[Nagpal et al.(2009)Nagpal, Lindquist,
  Oh, and Norris]NLO09
authorP. Nagpal,
  authorN. C. Lindquist,
  authorS.-H. Oh, and
  authorD. J. Norris,
  journalScience volume325,
  pages594 (year2009).

[Zeng et al.(2014)Zeng, Baillargeat, Hod,
  and Yong]ZBH14
authorS. Zeng,
  authorD. Baillargeat,
  authorH.-P. Hod, and
  authorK.-T. Yong,
  journalChem. Soc. Rev. volume43,
  pages3426 (year2014).

[Clavero(2014)]C14
authorC. Clavero,
  journalNat. Photon. volume8,
  pages95 (year2014).

[Danckwerts and Novotny(2007)]DN07
authorM. Danckwerts and
  authorL. Novotny,
  journalPhys. Rev. Lett. volume98,
  pages026104 (year2007).

[Zia et al.(2006)Zia, Schuller, Chandran,
  and Brongersma]ZSC06
authorR. Zia,
  authorJ. A. Schuller,
  authorA. Chandran, and
  authorM. L. Brongersma,
  journalMater. Today volume9,
  pages20 (year2006).

[Ropers et al.(2007)Ropers, Neacsu,
  Elsaesser, Albrecht, Raschke, and Lienau]RNE07
authorC. Ropers,
  authorC. C. Neacsu,
  authorT. Elsaesser,
  authorM. Albrecht,
  authorM. B. Raschke,
  and authorC. Lienau,
  journalNano Lett. volume7,
  pages2784 (year2007).

[Alcaraz Iranzo et al.(2018)Alcaraz
  Iranzo, Nanot, Dias, Epstein, Peng, Efetov, Lundeberg, Parret, Osmond, Hong
  et al.]AND18
authorD. Alcaraz Iranzo,
  authorS. Nanot,
  authorE. J. C. Dias,
  authorI. Epstein,
  authorC. Peng,
  authorD. K. Efetov,
  authorM. B. Lundeberg,
  authorR. Parret,
  authorJ. Osmond,
  authorJ.-Y. Hong,
  et al., journalScience
  volume360, pages291 (year2018).

[Fei et al.(2012)Fei, Rodin, Andreev,
  Bao, McLeod, Wagner, Zhang, Zhao, Thiemens, Dominguez et al.]FRA12
authorZ. Fei,
  authorA. S. Rodin,
  authorG. O. Andreev,
  authorW. Bao,
  authorA. S. McLeod,
  authorM. Wagner,
  authorL. M. Zhang,
  authorZ. Zhao,
  authorM. Thiemens,
  authorG. Dominguez,
  et al., journalNature
  volume487, pages82 (year2012).

[Chen et al.(2012)Chen, Badioli,
  Alonso-González, Thongrattanasiri, Huth, Osmond, Spasenović, Centeno,
  Pesquera, Godignon et al.]paper196
authorJ. Chen,
  authorM. Badioli,
  authorP. Alonso-González,
  authorS. Thongrattanasiri,
  authorF. Huth,
  authorJ. Osmond,
  authorM. Spasenović,
  authorA. Centeno,
  authorA. Pesquera,
  authorP. Godignon,
  et al., journalNature
  volume487, pages77 (year2012).

[Bharadwaj et al.(2011)Bharadwaj,
  Bouhelier, and Novotny]BBN11
authorP. Bharadwaj,
  authorA. Bouhelier,
  and authorL. Novotny,
  journalPhys. Rev. Lett. volume106,
  pages226802 (year2011).

[Nikitin et al.(2011)Nikitin, Guinea,
  García-Vidal, and Martín-Moreno]NGG11_2
authorA. Y. Nikitin,
  authorF. Guinea,
  authorF. J. García-Vidal,
  and
  authorL. Martín-Moreno,
  journalPhys. Rev. B volume84,
  pages195446 (year2011).

[Ritchie(1957)]R1957
authorR. H. Ritchie,
  journalPhys. Rev. volume106,
  pages874 (year1957).

[Powell and Swan(1959)]PS1959
authorC. J. Powell and
  authorJ. B. Swan,
  journalPhys. Rev. volume115,
  pages869 (year1959).

[García de Abajo(2010)]paper149
authorF. J. García de Abajo,
  journalRev. Mod. Phys. volume82,
  pages209 (year2010).

[Persson and Baratoff(1992)]PB92
authorB. N. J. Persson
  and authorA. Baratoff,
  journalPhys. Rev. Lett. volume68,
  pages3224 (year1992).

[Parzefall et al.(2015)Parzefall,
  Bharadwaj, Jain, Taniguchi, Watanabe, and Novotny]PPJ15
authorM. Parzefall,
  authorP. Bharadwaj,
  authorA. Jain,
  authorT. Taniguchi,
  authorK. Watanabe, and
  authorL. Novotny,
  journalNat. Nanotech. volume0,
  pagesDOI: 10.1038/NNANO.2015.203 (year2015).

[Walters et al.(2010)Walters, van Loon,
  Brunets, Schmitz, and Polman]WVB10
authorR. J. Walters,
  authorR. V. A. van Loon,
  authorI. Brunets,
  authorJ. Schmitz, and
  authorA. Polman,
  journalNat. Mater. volume9,
  pages21 (year2010).

[Zhang et al.(2013)Zhang, Allen, and
  DeCorby]ZAD13
authorM. C. Zhang,
  authorT. W. Allen, and
  authorR. G. DeCorby,
  journalAppl. Phys. Lett. volume103,
  pages071109 (year2013).

[Svintson et al.(2016)Svintson,
  Devizorova, Otsuji, and Ryzhii]SDO16
authorD. Svintson,
  authorZ. Devizorova,
  authorT. Otsuji, and
  authorV. Ryzhii,
  journalPhys. Rev. B volume94,
  pages115301 (year2016).

[Svintsov et al.(2016)Svintsov,
  Devizorova, Otsuji, and Ryzhii]SDR16
authorD. Svintsov,
  authorZ. Devizorova,
  authorT. Otsuji, and
  authorV. Ryzhii,
  journalPhys. Rev. B volume94,
  pages115301 (year2016).

[Woessner et al.(2017)Woessner, Misra,
  Cao, Torre, Mishchenko, Lundeberg, Watanabe, Taniguchi, Polini, Novoselov
  et al.]WMC17
authorA. Woessner,
  authorA. Misra,
  authorY. Cao,
  authorI. Torre,
  authorA. Mishchenko,
  authorM. B. Lundeberg,
  authorK. Watanabe,
  authorT. Taniguchi,
  authorM. Polini,
  authorK. S. Novoselov,
  et al., journalACS Photon.
  volume4, pages1461 (year2017).

[de Vega and García de Abajo(2017)]paper295
authorS. de Vega and
  authorF. J. García de Abajo,
  journalACS Photon. volume4,
  pages2367 (year2017).

[Castro Neto et al.(2009)Castro Neto,
  Guinea, Peres, Novoselov, and Geim]CGP09
authorA. H. Castro Neto,
  authorF. Guinea,
  authorN. M. R. Peres,
  authorK. S. Novoselov,
  and authorA. K. Geim,
  journalRev. Mod. Phys. volume81,
  pages109 (year2009).

[Wunsch et al.(2006)Wunsch, Stauber,
  Sols, and Guinea]WSS06
authorB. Wunsch,
  authorT. Stauber,
  authorF. Sols, and
  authorF. Guinea,
  journalNew J. Phys. volume8,
  pages318 (year2006).

[Hwang and Das Sarma(2007)]HD07
authorE. H. Hwang and
  authorS. Das Sarma,
  journalPhys. Rev. B volume75,
  pages205418 (year2007).

[Geick et al.(1966)Geick, Perry, and
  Rupprecht]GPR1966
authorR. Geick,
  authorC. H. Perry, and
  authorG. Rupprecht,
  journalPhys. Rev. volume146,
  pages543 (year1966).

[Palik(1985)]P1985
authorE. D. Palik,
  titleHandbook of Optical Constants of Solids
  (publisherAcademic Press, addressSan Diego,
  year1985).

[Filatova and Konashuk(2015)]FK15
authorE. O. Filatova and
  authorA. S. Konashuk,
  journalJ. Phys. Chem. C volume119,
  pages20755 (year2015).



