Plasmon-assisted resonant tunneling in graphene-based heterostructures

TL;DR

This paper presents a theory of plasmon-assisted resonant tunneling in graphene heterostructures, revealing sharp peaks in tunneling current and plasmon emission spectra due to dynamic Coulomb screening and surface plasmon resonances.

Contribution

It introduces the concept of plasmaronic resonance in graphene heterostructures and analyzes its dependence on interlayer twist and plasmon lifetime.

Findings

Resonant peaks in tunneling current at specific voltages

Sharp plasmon emission spectra corresponding to surface plasmon energies

Resonance persists despite interlayer twist and varies with carrier density

Abstract

We develop a theory of electron tunneling accompanied by carrier-carrier scattering in graphene - insulator - graphene heterostructures. Due to the dynamic screening of Coulomb interaction, the scattering-aided tunneling is resonantly enhanced if the transferred energy and momentum correspond to those of surface plasmons. We reveal the possible experimental manifestations of such plasmon-assisted tunneling in current-voltage curves and plasmon emission spectra of graphene-based tunnel junctions. We find that inelastic current and plasmon emission rates have sharp peaks at voltages providing equal energies, momenta and group velocities of plasmons and interlayer single-particle excitations. The strength of this resonance, which we call plasmaronic resonance, is limited by interlayer twist and plasmon lifetime. The onset of plasmon-assisted tunneling can be also marked by a cusp in the junction $I(V)$-curve at low temperatures, and the threshold voltage for such tunneling weakly depends on carrier density and persists in the presence of interlayer twist.

Equations

$$I_{t\to b}=e{g}^2\sum _{\begin{array}{c}\mathbf{p}{\mathbf{p}}_1\mathbf{q}\\ ls{s}_1{s}^{\prime }{s}_1^{\prime }\end{array}}{W}_{fi}{f}_{\mathbf{p}}^{+s}[1-f_{{\mathbf{p}}^{\prime }}^{-s^{\prime }}]f_{{\mathbf{p}}_1}^{l{s}_1}[1-f_{{\mathbf{p}}_1^{\prime }}^{l{s}_1^{\prime }}],$$

$$\begin{gathered} I_{t\to b}=\frac{2e}{\pi }\underset{-\infty }{\overset{+\infty }{\int }}d\omega \sum _{\mathbf{q},l}{\Pi }_{+-}^{\prime \prime }(q,\omega ){\Pi }_{ll}^{\prime \prime }\left(q,\omega \right)\times \\ {\left|V_{+l,-l}\right|}^2{N}_{\omega -eV}\left[N_{\omega }+1\right], \end{gathered}$$

$${\Pi }_{l{l}^{\prime }}=\frac{g}{A}\sum _{\begin{array}{c}\mathbf{p}s{s}^{\prime }\\ {\mathbf{p}}^{\prime }=\mathbf{p}+\mathbf{q}\end{array}}\frac{1+s{s}^{\prime }\cos {\theta }_{\mathbf{p}{\mathbf{p}}^{\prime }}}{2}\frac{f_{\mathbf{p}}^{ls}-f_{{\mathbf{p}}^{\prime }}^{l^{\prime }{s}^{\prime }}}{\omega +i\delta -{ϵ}_{\mathbf{p}}^{ls}+{ϵ}_{{\mathbf{p}}^{\prime }}^{l^{\prime }{s}^{\prime }}}.$$

$$\begin{gathered} V_{+l,-l}=\frac{V_0(\mathbf{q})}{\epsilon (\mathbf{q},\omega )}\frac{z_{\pm }}{d}(1-e^{-qd})\times \\ \left[1-V_0(\mathbf{q}){\Pi }_{-l-l}(\mathbf{q},\omega )(1+e^{-qd})\right], \end{gathered}$$

$$I_{t\to b}^{\mathrm{pl},\mathrm{em}}=I_{t\to b}^{\mathrm{pl},\mathrm{em}}+I_{t\to b}^{\mathrm{pl},\mathrm{abs}},$$

$$I_{t\to b}^{\mathrm{pl},\mathrm{em}}=2\pi e\sum _{\mathbf{q}p}{\left|\frac{e{\phi }_{\mathbf{q}\pm }^p}{2}\right|}^2{\Pi }_{tb}^{\prime \prime }(\mathbf{q},{\omega }_{\mathbf{q}}^p)[N_{{\omega }_{\mathbf{q}}^p}+1]N_{{\omega }_{\mathbf{q}}^p-eV}.$$

$$\begin{gathered} \frac{(e{\phi }_{\mathbf{q}}^p{)}_{\pm }^2}{2}=V_0(\mathbf{q}){\left|\frac{z_{\pm }}{d}\right|}^2{\left(1-e^{-qd}\right)}^2\times \\ \times \frac{{\left[1-V_0(\mathbf{q}){\Pi }_{++}^{\prime }(\mathbf{q},{\omega }_q^p)(1+e^{-qd})\right]}^2}{\frac{\partial {\epsilon }^{\prime }}{\partial {\omega }_q^p}\left|1-V_0(\mathbf{q}){\Pi }_{++}^{\prime }(\mathbf{q},{\omega }_q^p)(1-e^{-2qd})\right|}. \end{gathered}$$

$${\Delta }^{\ast }+qv={\omega }_q^p,$$

$$\partial {\omega }_q^p/\partial q=v.$$

$$I_{t\to b}^{\mathrm{pl},\mathrm{em}}(\Delta )\approx I_0{\ln \left|\frac{q^2{\partial }^2{\omega }_q/\partial q^2}{2(\Delta -{\Delta }^{\ast })}\right||}_{q=q^{\ast }},$$

$$I_0={\left|\frac{e{\phi }_{\mathbf{q}\pm }^p}{2}\right|}^2{\frac{eq\left(N_{\omega }+1\right)N_{\omega -eV}{\tilde{\Pi }}_{\pm }(q)}{2\pi {\hslash }^2\sqrt{qv{\partial }^2{\omega }_q/\partial q^2}}|}_{\begin{array}{c}q=q^{\ast }\\ \omega ={\omega }_{q^{\ast }}\end{array}}$$

$$I_{t\to b}^{\mathrm{pl},\mathrm{em}}({\Delta }^{\ast })\approx I_0{\ln \left|\frac{q^2{\partial }^2{\omega }_q/\partial q^2}{{\epsilon }^{\prime \prime }/[\partial {\epsilon }^{\prime }/\partial \omega ]}\right||}_{q=q^{\ast }}.$$

$${\Pi }_{\pm }^{(T)}(\mathbf{q},\omega )=\frac{1}{3}\sum _{i=1}^3{\Pi }_{\pm }(\mathbf{q}+\Delta {\mathbf{q}}_i,\omega ).$$

$$I_{t\to b}^{\mathrm{pl},\mathrm{em}}({\Delta }^{\ast })\approx \frac{I_0}{2\pi }\sqrt{\frac{q_0^2{\frac{{\partial }^2\omega }{\partial q^2}|}_{q_0}}{v\Delta q}}\ln \left|\frac{q_c^2{\partial }^2{\omega }_q/\partial q_c^2}{8v\Delta q}\right|.$$

$$\frac{d{N}_{pl}}{dt}=\frac{1}{e}\left[I_{t\to b}^{\mathrm{pl},\mathrm{em}}+I_{b\to t}^{\mathrm{pl},\mathrm{em}}-I_{t\to b}^{\mathrm{pl},\mathrm{abs}}-I_{b\to t}^{\mathrm{pl},\mathrm{abs}}\right].$$

$$\frac{d{N}_{ph}}{dt}=\frac{{\gamma }_{rad}}{{\gamma }_{rad}+{\gamma }_{sc}}\frac{d{N}_{pl}}{dt}.$$

$${\gamma }_{rad}\approx \frac{8LdW{\omega }^4}{3{\pi }^3{c}^3}{\left[2+\frac{3}{K}+\frac{1}{K^2}\right]}^{-1},$$

$$Z_{rad}{\sigma }^{\prime \prime }({\omega }_{\mathbf{q}})W/L\sim 1,$$

$${\gamma }_{rad,\max }\sim {\omega }_{\mathbf{q}}^{\prime },$$

$$J(V,\Delta )=\sum _{\mathbf{p}\mathbf{q}s{s}^{\prime }}\theta (eV-{\omega }_{\mathbf{q}})\delta ({ϵ}_{\mathbf{p}+}^{s+}-{ϵ}_{\mathbf{p}+}^{s^{\prime }-}-{\omega }_{\mathbf{q}}),$$

$$e{V}_{th}={\omega }_q,$$

$${\omega }_q={\Delta }_{th}-(q+\Delta q)v.$$

$$e{V}_{th}=\frac{\Delta -v\Delta q}{1+v/s}.$$

$${\Delta }_{th}=\left(1+\frac{s}{v}\right)({\epsilon }_{F-}-{\epsilon }_{F+})-s\Delta q.$$

$$\phi \left(\mathbf{r}\right)=\frac{1}{2}\sum _{\mathbf{q},p}{\phi }_q^p(z)e^{i\mathbf{qr}-i\omega t}+c.c.$$

$$I_{t\to b}^{pl,em}=2\pi e\sum _{\mathbf{p}\mathbf{q}s{s}^{\prime }}{f}_{\mathbf{p}}^{+s}(1-f_{\mathbf{p}-\mathbf{q}}^{-s^{\prime }})|e{\phi }_{tb}^{em}{|}^2\delta ({\epsilon }_{\mathbf{p}}^{+s}-{\epsilon }_{\mathbf{p}-\mathbf{q}}^{-s^{\prime }}-{\omega }_q)$$

$$I_{t\to b}^{pl,abs}=2\pi e\sum _{\mathbf{p}\mathbf{q}s{s}^{\prime }}{f}_{\mathbf{p}}^{+s}(1-f_{\mathbf{p}+\mathbf{q}}^{-s^{\prime }})|e{\phi }_{tb}^{abs}{|}^2\delta ({\epsilon }_{\mathbf{p}}^{+s}-{\epsilon }_{\mathbf{p}+\mathbf{q}}^{-s^{\prime }}+{\omega }_q)$$

$$|e{\phi }_{tb}^{em}{|}^2=[N_{{\omega }_q^p}+1]|e{\phi }_{tb}{|}^2,$$

$$|e{\phi }_{tb}^{abs}{|}^2=N_{{\omega }_q^p}|e{\phi }_{tb}{|}^2.$$

$$e{\phi }_{tb}={\int }_{-\infty }^{+\infty }dz{\psi }_t^{\ast }(z)e{\phi }_{\mathbf{q}}(z){\psi }_b(z),$$

$${\phi }_{\mathbf{q}}(z)={\phi }_0{S}_q(z),$$

$$S_q(z)=\{\begin{array}{l}{S}_{+}{e}^{-q(z-d/2)},z>d/2\\ \frac{S_{+}\sinh [q(z+\frac{d}{2})]-S_{-}\sinh [q(z-\frac{d}{2})]}{\sinh qd},|z|<d/2,\\ S_{-}{e}^{q(z+d/2)},z<-d/2,\end{array}$$

$$\begin{gathered} W=\int d^{3}r\frac{\kappa \mathbf{E}{\mathbf{E}}^{\ast }}{16\pi }-\frac{A}{4}\sum _{l=\pm }{\frac{\partial {\sigma }_l^{\prime \prime }}{\partial \omega }|}_{{\omega }_q}{{\mathbf{E}}_{||}{\mathbf{E}}_{||}^{\ast }|}_{z=ld/2}= \\ =\frac{A{\phi }_0^{2}q{\omega }_q{S}_{-}}{8\pi }{\frac{\partial \epsilon }{\partial \omega }|}_{{\omega }_q}, \end{gathered}$$

Full text

I Introduction

The plasmons supported by two-dimensional (2d) electron systems, including graphene, can be confined at the distance by two orders of magnitude smaller than free-space light wavelength Woessner et al. (2015). High field confinement in graphene plasmon modes leads to an enhanced electron-plasmon interaction Koppens et al. (2011) and a number of interesting phenomena. The formation of new quasiparticles of bound holes and plasmons, called plasmarons Bostwick et al. (2010); perfect light absorption Thongrattanasiri et al. (2012); ultrafast recombination of photoexcited carriers mediated by plasmon emission Rana (2008); Watanabe et al. (2013) are among the brightest manifestations of electron-plasmon coupling. The strength of light-matter interaction is the highest for acoustic plasmon modes supported by 2d bilayers and gated 2d systems Ryzhii et al. (2007); Alonso-Gonzalez et al. (2017); thereat the electromagnetic energy is concentrated at the length scale of interlayer spacing $d$. Such bilayers are actively investigated both as building blocks of novel resonant tunneling diodes Gaskell et al. (2015), transistors Britnell et al. (2013), and as a polygon for the fundamental studies of tunneling in the presence of chirality Greenaway et al. (2015) and field-controlled interlayer twist Wallbank et al. (2016).

What new effects can stem from strong electron-plasmon interaction in tunnel-coupled graphene layers? One might naively expect the emergence of steps in the dependence of tunnel conductivity on interlayer bias Vdovin et al. (2016). These steps commonly occur at voltages $V=\hbar\omega_{i}/e$ Esaki et al. (1968), where $\hbar\omega_{i}$ are the characteristic energies of collective excitations. However, the plasmon dispersion is soft in two dimensions, i.e. its frequency tends to zero in the long-wavelength limit. Therefore, generally there is no preferable frequency for plasmon-assisted tunneling in 2d, and plasmonic fingerprints in tunneling can appear only at very specific conditions Zhang et al. (1994); Belenov et al. (1987), e.g., at the anticrossing of different plasmon modes Kempa et al. (2001). In this paper, we find that in aligned graphene double layers the plasmon-assisted tunneling can surprisingly lead to a resonant enhancement of tunnel current. The interlayer bias providing the resonance corresponds to the energy, momentum, and group velocity matching between plasmons and interlayer single-particle excitations. This effect is inherent to the graphene’s linear band structure. In some sense, it is similar to the formation of plasmarons a single graphene layer due to the consonance between plasmon and electron motion Bostwick et al. (2010) – for this reason, we call this effect plasmaronic resonance.

The plasmaronic resonance can manifest itself also in the spectra of the tunnel junctions’ electroluminescence Lambe and McCarthy (1976); Parzefall et al. (2017). We show that in aligned layers the integrated luminescence can demonstrate a strong spike at the resonant voltage, and be considerably larger than the recently observed luminescence from twisted layers Yadav et al. (2016). In the presence of twist, an only remainder of plasmons is the fine structure of low-temperature $I(V)$-curves due to switch-on of plasmon-assisted tunneling.

So far, a considerable work has been done in the studies of plasmon spectra and damping/gain in graphene double layers Hwang and Das Sarma (2009); Svintsov et al. (2013), including the tunneling effects Svintsov et al. (2016); Sensale-Rodriguez (2013). On the other hand, there exists a number of theoretical works on various tunneling channels in these systems: elastic Feenstra et al. (2012); Vasko (2013); Brey (2014), including the effects of finite carrier lifetime Guerrero-Becerra et al. (2016), phonon-Amorim et al. (2016) and photon-assisted Ryzhii et al. (2013). However, no attempt has been made to calculate the plasmon-assisted tunnel current. The present theory is based on the calculation of electron interlayer transition rates due to the dynamically screened Coulomb interaction with other carriers. By extracting the plasmon-pole contribution from the scattering-assisted tunnel current, we are able to find the current due to plasmons, and plasmon emission rates as well. As an added benefit, our formalism allows to calculate the full inelastic tunnel current, including the non-resonant scattering contributions. However, as we find, the largest fraction of Coulomb-scattering-assisted current is due to emission of plasmons. We also present a simplified approach for plasmon-assisted current based on the calculation of electron tunneling due to fluctuating electric fields, and evaluation of the magnitude of fluctuations with second quantization procedure. These two approaches coincide in the limit of weak electromagnetic dissipation.

II Theory of many-particle and plasmon-assisted tunneling

Electron states in coupled graphene layers with small interlayer twist can be labelled by in-plane momentum ${\bf p}$, the band index $s=\pm 1$ and the index $l=\pm 1$ governing the vertical localization of electron wave function Vasko (2013). The respective energies are $\epsilon^{ls}_{\bf p}=spv_{0}+l\sqrt{\Delta^{2}/4+\Omega^{2}}$, where $v_{0}$ is the Fermi velocity, $\Delta$ is the voltage-induced splitting of bands in neighboring layers (band offset, see Fig. 1), and $\Omega$ is the tunnel splitting. For strong bias and/or weak tunneling, $\Delta\gg\Omega$, the state $l=+1$ can be regarded as belonging to the top layer and $l=-1$ – to the bottom one. The Coulomb interaction couples the states in neighboring layers and induces inelastic tunneling current. The current from the forward biased (top) layer can be presented as

[TABLE]

where $g=4$ is the spin-valley degeneracy factor, $f^{ls}_{\bf p}$ is the occupancy of the quantum state with respective set of quantum numbers, and $W_{fi}$ is the golden rule transition probability between two-particle states $\ket{+{\bf p}s,l{\bf p}_{1}s_{1}}$ and $\ket{-{\bf p}^{\prime}s^{\prime},l^{\prime}{\bf p^{\prime}}_{1}s^{\prime}_{1}}$ due to Coulomb interaction. The full current including the reverse component is $I=I_{t\rightarrow b}[1-e^{-eV/kT}]$, where $V$ is the interlayer voltage. We note that $eV=\Delta+\varepsilon_{F+}-\varepsilon_{F-}$, where $\varepsilon_{F\pm}$ are the Fermi energies in respective layers. The quantities $\Delta$ and $\varepsilon_{F\pm}$ can be controlled independently with doping or extra gates.

A sequence of transformations common in the theory of Coulomb scattering phenomena Zheng and Das Sarma (1996) leads us to a particularly simple form of the tunneling current (we set $\hbar\equiv 1$)

[TABLE]

where $N_{\omega}=[e^{\omega/T}-1]^{-1}$ is the Bose distribution and $\Pi_{ll^{\prime}}=\Pi^{\prime}_{ll^{\prime}}+i\Pi^{\prime\prime}_{ll^{\prime}}$ is the intra- ($l=l^{\prime}$) or inter-layer ($l\neq l^{\prime}$) polarizability:

[TABLE]

The amplitude $V_{+l,-l}$ describes the tunneling of electron from top to bottom layer upon Coulomb interaction with a carrier in the $l$-th layer. The latter depends on frequency and wave vector, primarily, due to the dynamic screening of Coulomb interaction. In the dipole approximation, the transition amplitude becomes

[TABLE]

where $z_{\pm}=\int_{-\infty}^{+\infty}dz\psi_{+}(z)z\psi_{-}(z)$ is the coordinate matrix element between initial and final states, $V_{0}({\bf q})=2\pi e^{2}/\kappa|{\bf q}|$ is the bare Coulomb interaction in two dimensions, $\kappa$ is the background dielectric constant, and $\varepsilon({\bf q},\omega)$ is the dynamic screening function of the double layer Hwang and Das Sarma (2009). The factor $1-e^{-qd}$ in Eq. (4) is due to the long range of Coulomb interaction. Indeed, the potential created by an electron is almost identical on both layers. Hence, the tunneling is weakened for long wavelengths (or slowly varying potentials) $q\rightarrow 0$. This prefactor cancels the low-angle scattering singularity in bare Coulomb interaction, making the full expression for inelastic current convergent even in the absence of screening. Still, the screening has an important impact on interlayer tunneling, which is seen from the resonant enhancement of matrix element at $\varepsilon({\bf q},\omega)\rightarrow 0$. This is nothing but tunneling accompanied by emission of surface plasmons.

To extract the plasmonic contribution, one can expand the dielectric function in the vicinity of plasmon poles $\varepsilon({\bf q},\omega)\approx\left[\partial\varepsilon^{\prime}/\partial\omega\right](\omega-\omega^{p}_{\bf q})+i\varepsilon^{\prime\prime}$ (here $p=+1$ ($-1$) corresponds to optical (acoustic) modes, for details see Appendix B). Assuming the electromagnetic dissipation to be small, $|\varepsilon^{\prime\prime}/\varepsilon^{\prime}|\ll 1$, we arrive at the expression for plasmon-assisted component of the net tunneling current $I^{\rm pl}_{t\rightarrow b}$. This can be conveniently split into the emission ($I^{\rm pl,em}_{t\rightarrow b}$) and absorption ($I^{\rm pl,abs}_{t\rightarrow b}$) contributions

[TABLE]

The quantity $(e\varphi^{p}_{\bf q})_{\pm}$ can be viewed as a matrix element of electron interaction with zero-point field of plasmon:

[TABLE]

The same result could be obtained by calculating the electron transition rates due to the interaction with zero-point and thermal longitudinal fluctuations of electromagnetic field. The magnitude of these fluctuations can be found from quantum-classical correspondence Belenov et al. (1987), i.e. by equating the classical energy of electromagnetic field in the dispersive medium to $N_{\omega}\hbar\omega$. This procedure is described in detail in Appendix A.

III Manifestations of plasmon-assisted tunneling

With the general formalism of calculation developed, we start discussing the possible experimental manifestations of plasmon-assisted tunneling. We consider three such effects: (1) resonant enhancement of tunnel current due to group velocity coincidence between plasmons and inter-layer single particle excitations, which we call plasmaronic resonance (2) plasmonic electroluminescence of graphene tunnel junctions (3) fine structure of the low-temperature $I(V)$-curves due to onset of the plasmon emission.

III.1 Plasmaronic resonance in tunnel current

The imaginary part of interlayer polarizability $\Pi^{\prime\prime}_{\pm}({\bf q},\omega)$ in aligned layers has a square-root singularity at the threshold of interlayer excitations, $\Pi^{\prime\prime}_{\pm}({\bf q},\omega)\propto|(\Delta\pm qv)^{2}-\omega^{2}|^{-1/2}$. This singularity can be explained as resulting from prolonged interaction between electron and hole with collinear momenta and, hence, equal velocities. Similar singularities exist in the polarizability of a single graphene layer at the threshold of Landau dampingHwang and Das Sarma (2007) $\Pi^{\prime\prime}_{ll}({\bf q},\omega)\propto|q^{2}v^{2}-\omega^{2}|^{-1/2}$, and were recently assessed experimentally Lundeberg et al. (2017a). A pronounced effect of such ”collinear singularities” is that the plasmon phase velocity always lies above the Fermi velocity, and the Landau damping of plasmons is absent Ryzhii et al. (2007).

The group velocity of plasmons can be, however, equal to or below the Fermi velocity. This applies both to the modes supported by a single layer and graphene double layer as well, and is shown in Fig. 2 B. When the line of interlayer tunneling singularities $\omega=\Delta+qv$ approaches the tangent with plasmon dispersion, the plasmon-assisted tunneling current is resonantly enhanced. The resonant interlayer band offset $\Delta^{*}$ is determined from

[TABLE]

In the vicinity of resonance the plasmon-assisted contribution grows as

[TABLE]

where the large logarithm is evaluated at $q=q^{*}$ which is the momentum of plasmons in resonance with interlayer excitations. The characteristic current in Eq. (10) is

[TABLE]

where $\widetilde{\Pi}_{\pm}$ is non-singular part of $\Pi_{\pm}$ (i.e. $\Pi_{\pm}$ without the square-root singularity). The logarithmic growth of the current at the resonance is limited by plasmon damping. The latter was assumed to be infinitesimal in Eq. (5) but is automatically taken into account in many-particle formalism, Eq. (2), where both real and imaginary parts of dielectric function contribute to screening. It is possible to show that the damping-limited resonant value of current is, roughly

[TABLE]

It is possible to tune the structure parameters to achieve low damping of plasmons with $q=q^{*}$ by blocking the interband transitions, therefore making the resonant contribution very large. The dependence of full inelastic current on band offset $\Delta$ at fixed carrier densities in the layers is shown in Fig. 2. For the sake of generality, we normalize the current by $z^{2}_{\pm}/d^{2}$, therefore getting rid of material-dependent tunneling exponent $z_{\pm}\propto e^{-\varkappa d}$, where $\varkappa$ is the decay length of electron wave function. The two peaks in Fig. 2 correspond to the plasmaronic resonances on acoustic and optical modes. It is worth noting that the resonance with optical mode is possible only for different conductivities of the two layers. Otherwise, the average mode field between layers is zero and the tunneling matrix elements turn to zero as well.

In the perfectly aligned layers, the inelastic tunneling current is readily seen (and even surpasses) the elastic current, as shown in Fig. 5. The elastic current scales as $I_{\max}\gamma/[\gamma^{2}+\Delta^{2}]$, where $\gamma$ is the collisional broadening of the resonance. The mentioned effect is, however, highly sensitive to the twist between graphene layers Mishchenko et al. (2014). With the neglect of emerging weak tunneling between dissimilar sublattices, the general expression for inelastic current (2) still holds, but the interlayer polarizability is now angle-dependent. Denoting the wave vectors connecting the $K$-points in the neighboring layers as $\Delta{\bf q}_{i}$ ($i=1...3$), we can write the polarizability in the presence of twist $\Pi^{(T)}_{\pm}({\bf q},\omega)$ as

[TABLE]

When the twist wave vector is small compared to the plasmon wave vector at the resonance, $\Delta q\ll q^{*}$, the twist-limited contribution to the tunnel current can be estimated as

[TABLE]

III.2 Plasmon emission and junction electroluminescence

The emission of surface plasmons upon resonant tunneling can be detected not only implicitly by analyzing the features of inelastic current. Recent advances in near-field optical microscopy and electromagnetic sensing at the nanoscale Woessner et al. (2015); Lundeberg et al. (2017b) allow a direct measurement of plasmon emission rates. This emission rate, $dN_{pl}/dt$, is obtained by a simple rearrangement of terms in the expression for plasmon-assisted tunnel current

[TABLE]

Naturally, the bias dependence of integrated plasmon emission rate inherits all resonant features of plasmon-assisted current. This is shown in Fig. 3, where the characteristic peaks are the plasmaronic resonances discussed above. At large band offsets, the emission will be overwhelmed by absorption, and the double layer structure can operate as a resonant plasmonic photodetector Goodfellow et al. (2015); Ryzhii et al. (2014). By appropriate choice of doping, it is possible to achieve the plasmaronic resonance both in emission and absorption rates.

It is notable that the energy spectrum of plasmon emission, $dN_{pl,\omega}/dt$, possesses specific resonances at any bias. These resonances reflect the singular nature of interlayer polarizability $\Pi_{\pm}({\bf q,\omega})$ and occur at frequencies satisfying $\omega=\Delta\pm q^{p}(\omega)$, where $q^{p}(\omega)$ is the inverse of the dispersion law for the $p$-th mode Svintsov et al. (2016). The singularities in the frequency spectrum are integrable except for the case of merging singularities at the plasmaronic resonance.

The spectra of emitted surface plasmons can be studied via analysis of tunnel junction electroluminescence. Such electroluminescence is commonly a two-step process including the excitation of surface plasmon upon inelastic tunneling and its subsequent radiative decay into free-space modes. The direct emission of photon is less probable due to the large spatial extent of photonic modes and small photonic density of states as well. The radiative decay of plasmon with characteristic rate $\gamma_{rad}$ generally competes with damping due to the free carrier scattering and interband absorption having the rate $\gamma_{sc}$. When both decay rates are smaller than plasmon frequency, the plasmon-to-photon conversion rate can be presented as

[TABLE]

The main obstacle to the observation of plasmonic electroluminescence is the smallness of plasmon radiative damping. Eventually, this is due to the smallness of dipole moment of tunnel junction, and smallness of plasmon velocity $s$ compared to the speed of light. For the tunnel junction of width $W$ and length $L$ tuned to the fundamental mode of surface plasmon $qL=\pi$, the radiative decay rate can be estimated in the dipole approximation to be

[TABLE]

where $K=4\alpha_{c}k_{F}d$, $k_{F}$ is the Fermi wave vector and $\alpha_{c}=e^{2}/\hbar\kappa v$ is the coupling constant. For the acoustic mode of energy $\hbar\omega=200$ meV (corresponding to the resonance A in Fig. 2), one can estimate $\gamma_{rad}\approx 2\times 10^{8}$ s*-1*. This is quite a low value compared to the absorption rate due to electron scattering $\gamma_{sc}\approx 10^{12}$ s*-1* Woessner et al. (2015). The resulting efficiency of plasmon-to-photon conversion is about $\gamma_{rad}/\gamma_{sc}\approx 2\times 10^{-4}$.

This efficiency can be increased dramatically by coupling a resonant-tunneling structure with a nanoscale antenna Parzefall et al. (2017). An estimate for the radiative decay rate in this case can be obtained by solving the plasmon dispersion equation in a bounded double-layer structure loaded by an antenna impedance $Z_{rad}$. The maximum decay rate is achieved when the load impedance matches the impedance of a double-layer structure at the resonant frequency,

[TABLE]

the definite value of numerical factor depends on the shape and biasing of the double layer and is evaluated in Appendix B. Under optimal loading condition, the radiative decay rate becomes the same order as the eigenfrequency

[TABLE]

a numerical estimate for the structure with equal layer doping yields $\gamma_{rad,\max}=0.04\omega^{\prime}_{\bf q}$ for the lowest mode. As the scattering rate is well below the resonant frequency, the plasmon-to-photon conversion efficiency can be made close to unity.

III.3 Fine structure of the low-temperature $I(V)$-curves

At low temperatures, the tunneling with plasmon absorption if frozen out, while emission aided tunneling is possible only for $\omega_{q}<eV$ by the virtue of Pauli blocking. The combination of Pauli restriction and energy-momentum conservation results in suppression of inelastic current for certain range of voltages $V<V_{th}$ and band offset $\Delta<\Delta_{th}$, while above the threshold value the inelastic current switches on in a threshold-like manner. The threshold structure of the plasmon-aided current repeats that of an integral

[TABLE]

where ${\bf p}_{\pm}={\bf p}\pm({\bf q}+\Delta{\bf q})/2$ are the momenta of initial and final electrons, the Heaviside theta function is responsible for Pauli blocking, and the delta-function for the energy conservation. The intraband tunneling transitions are stronger, and therefore we consider the terms with $s=s^{\prime}$ only. Evaluation of integral (20) leads us to the following threshold condition

[TABLE]

The latter admits a simple geometrical interpretation shown in the inset of Fig. (4). The minimal frequency of plasmon in the domain of intraband tunneling $|\omega-\Delta|<(q+\Delta q)v$ (orange filled region) should lie below the line of Pauli blocking $\omega=eV$ (blue filled region). From this analysis we also see that finite interlayer twist $\Delta q$ reduces the threshold of plasmon emission upon tunneling.

For acoustic plasmon in graphene double layer with linear dispersion $\omega=sq$, the threshold condition (21) can be solved analytically to yield

[TABLE]

If the band offset $\Delta$ is fixed, the threshold voltage (23) weakly depends on carrier density because the plasmon velocity $s$ tends to the Fermi velocity at small interlayer distance $d$ Ryzhii et al. (2007); Principi et al. (2011); Svintsov et al. (2016), and the density dependence of $s$ is weak. This contrasts to the case of plasmon-assisted tunneling in bulk metal-insulator-metal junctions Lambe and McCarthy (1976) where the threshold voltage equals the plasmon energy, the latter scaling as square root of density.

If the carrier densities in graphene layers are fixed while band offset is swept, the threshold condition can be presented in an alternative form

[TABLE]

The step-like switch-on of tunnel current upon increase in band offset $\Delta$ is shown in Fig. (4) for fixed carrier densities and zero temperature. These cusps in the $I(V)$-curves become broadened very quickly with the temperature increase, as the broadening is governed both by smearing of Fermi distributions and interband damping of graphene plasmons.

IV Discussion and conclusions

We have theoretically identified several manifestations of plasmon-assisted tunneling in graphene-insulator-graphene junctions, the most striking of them being the plasmaronic resonance in tunnel current. The origin of this resonance is the enhanced interaction between plasmons and interlayer single-particle excitations due to the group velocity matching. The relation between discussed resonance and formation of plasmarons in a single graphene layer is elucidated as follows. Plasmaron is formed off the mass shell $\epsilon_{p}=pv_{0}$ at some energy separation $\delta\epsilon$ equal to the energy of plasmon quantum $\omega_{q}$. Contrary to 3D systems, the plasmons in two-dimensions have soft spectrum with energy tending to zero at long wavelength. A natural question arises: which plasmon wave vector $q^{*}$ provides the strongest interaction with electrons? The answer is that such a plasmon should have group velocity equal to the carrier velocityBostwick et al. (2010). In the case of interlayer tunneling, we are dealing essentially with on-shell electrons, however, the energy of interlayer single-particle excitations is tuned by interlayer bias $\Delta$. At some bias $\Delta^{*}$, the energies, momenta and group velocities of interlayer excitations and plasmons coincide. This bias corresponds to the resonantly large generation of surface plasmons by interlayer tunneling.

We note that the mentioned resonance is closely related to the presence of square-root singularities in the interlayer polarizability of double layer. Such singularities are inherent to the linear carrier dispersion in graphene Hwang and Das Sarma (2007) and are absent in systems with parabolic bands Stern (1967). Similar singularities in the polarizability of a single layer manifest themselves in a number of interesting phenomena, including ultrafast collinear scattering of photoexcited carriers König-Otto et al. (2016) and the absence of Landau damping for plasmons Ryzhii et al. (2007). These singularities can be smeared by electron-electron corrections to carrier dispersion and/or by vertex corrections Gangadharaiah et al. (2008). Though our original derivation of inelastic current was based on the scattering of noninteracting particles, the interaction effects can be conveniently included in the transformed equation (2) by replacing the bare polarizabilities $\Pi_{ll^{\prime}}$ with interacting ones. Here we just note that logarithmic interaction corrections to linear bands are most relevant for quasi-neutral graphene, while for doped samples the e-e interactions just enhance the band velocity under the Fermi surface Hwang et al. (2007). In most realistic situations, the broadening of predicted plasmaronic resonance would be governed by finite plasmon lifetime, Eq. (12), and interlayer twist, Eq. (14).

In the present calculation, we assumed the interlayer tunneling to be weak, so that the dielectric function of the double layer was not renormalized by tunneling. Such renormalization can be done Svintsov et al. (2016), and it would enhance the plasmon-assisted current. The reason for enhancement is the partial plasmon loss compensation by stimulated plasmon emission upon tunneling. At some critical strength of tunneling, corresponding to the complete undamping of plasmon modes, the current (2) would diverge. This effect is analogous to the enhancement of scattering rates in photoexcited semiconductors due to plasmon undamping Scott et al. (1992). The divergence would signalize on the onset of surface plasmon lasing, at this point one has to solve the coupled kinetic equations for electrons and plasmons for evaluation of tunnel current.

The present theory demonstrates the prospect of graphene heterostructures for resonant and voltage tunable light emission in the far infrared. Compared to the light sources based on tunneling injection and recombination of excitons Withers et al. (2015), the proposed structures can be scaled down to the nanometre dimensions due to high confinement of plasmons, and integrated in photonic and plasmonic waveguides. The proposed process of plasmon and photon generation is also different from interband recombination of electrons injected upon resonant tunneling considered in Khorasani (2014). The difference between these processes is the same as the difference of quantum-cascade lasing with vertical and diagonal radiative transitions. Importantly, the emission spectrum for diagonal transitions is resonant due to the singularities in the joint density of states between layers, while the emission spectrum for interband recombination of injected carriers repeats the smooth spectrum of emission in a single layer.

In conclusion, we have developed a theoretical formalism for the calculation of tunneling current accompanied by carrier-carrier scattering in graphene-insulator-graphene heterostructures. Our calculation shows that the main contribution of inelastic scattering-assisted current comes from emission of surface plasmons. The plasmon-assisted current can be resonantly enhanced if the energy, momentum and group velocity of interlayer excitations and plasmons coincide. This effect, which we call plasmaronic resonance, can also manifest itself in enhanced plasmon emission and electroluminescence of graphene-based junctions.

Acknowledgements.

This work was supported by the grants # 16-37-60110/16 and # 16-29-03402/16 of the Russian Foundation for Basic Research and by the Grant of the President of the Russian Federation. The authors are grateful to G. Alymov for valuable discussions.

Appendix A Second-quantization approach to plasmon-assisted tunneling

The plasmon-assisted tunneling current can be calculated as a current induced by random fluctuations of electric potential

[TABLE]

Using the Fermi golden rule, we present this current as

[TABLE]

for emission contribution, and

[TABLE]

for absorptive contribution. For classical field, the matrix elements are the same for emission and absorption. For quantized field, one has

[TABLE]

The tunneling matrix element is calculated as

[TABLE]

where $\psi^{*}_{l}(z)$ is the $z$-component of the wave function belonging to the $l$-th layer, and $\varphi^{p}_{q}(z)$ is the normalized electric potential in the plasmon mode with wave vector ${\bf q}$ and polarization $p$. The field magnitude can be parametrized by a single quantity $\varphi_{0}$, while the spatial dependence of electric potential is given by:

[TABLE]

with $S_{+}=e^{-qd}$, $S_{-}=1+i2\pi q\sigma^{\prime\prime}_{+}({\bf q},\omega_{q}^{p})\left(1-e^{-2qd}\right)/\kappa\omega_{q}$,

$\sigma_{+/-}({\bf q},\omega_{q}^{p})=i\omega_{q}e^{2}\Pi_{++/--}\left({\bf q},\omega_{q}\right)/q^{2}$ is the conductivity of top/bottom graphene layer. To determine the amplitude of plasmon field $\varphi_{0}$ we equate the classical field energy in a dispersive medium to the quantized energy $\hbar\omega^{p}_{\bf q}$. The classical energy of electromagnetic field is

[TABLE]

where ${\bf E}=\left({\bf E}_{||},E_{z}\right)=-\left(i{\bf q},\partial_{z}\right)\varphi_{q}(z)$ is the electric field in plasmon mode and $A$ is the sample area. This leads us to the final formula for $\varphi_{0}$:

[TABLE]

here we have introduced the permittivity of the double layer structure:

[TABLE]

Appendix B Plasmon poles in scattering-assisted tunneling

The current accompanied by emission of plasmons can be also derived by extracting the contribution to the integral (2) due to the poles of screening function $\varepsilon^{-1}({\bf q},\omega)$. Assuming the dissipation of electromagnetic energy to be small, one can determine the plasmon frequency $\omega_{q}^{p}$ from:

[TABLE]

If the frequency $\omega$ and momentum ${\bf q}$ satisfy the disperion relation (36), the transition amplitudes $V_{++,-+}$ with $V_{+-,--}$ are related as follows:

[TABLE]

With the help of Eqs. (36),(37) we can write down the tunneling current in the following way (keeping in mind $\omega\approx\omega_{q}$):

[TABLE]

In the limit $|\varepsilon^{\prime\prime}/\varepsilon^{\prime}|\ll 1$ one can make a substitution

[TABLE]

in the Eq. (B) after the last equality and obtain the resulting expression for the tunneling current:

[TABLE]

where

[TABLE]

is a dimensionless potential profile which equals the plasmon potential $S_{q}(z)$ (31).

Appendix C Tunneling matrix elements

Electron states in coupled graphene layers with small interlayer twist can be described using the Hamiltonian

[TABLE]

where the blocks ${{\hat{H}}}_{G\pm}$ describe isolated graphene layers, and $\hat{\mathcal{T}}$ describes tunnel hopping:

[TABLE]

Here $\Omega$ is the tunneling overlap integral and $\Delta{\bf q}_{j}$ are the wave vectors connecting the respective edges of hexagonal Brillouin zones in the layers. In the absence of interlayer twist ($\delta{\bf q}=0$), the tunneling matrix is diagonal, ${\hat{\mathcal{T}}}=\Omega{\hat{I}}$, where ${\hat{I}}$ is the identity matrix. In this case, the band and layer degrees of freedom are decoupled.

We proceed now to the evaluation of matrix element $\Omega$. The physical meaning of $\Omega$ is half the energy splitting between electron states in coupled graphene layers, as can be seen from diagonalization of Hamiltonian (42). On the other hand, this splitting can be estimated from a continuum model, where each graphene layer is represented by a delta-well Svintsov et al. (2016). The delta-well potential is

[TABLE]

where $U_{0}$ is the work function from graphene to the barrier material, and $m^{*}$ is the effective mass in the barrier. For boron nitride, $U_{0}\approx 1.5$ eV and $m^{*}\approx 0.5m_{0}$. The eigen functions in this potential are symmetric and anti-symmetric ones. The energy difference between these states is

[TABLE]

where $\varkappa=\sqrt{2m^{*}U_{0}}/\hbar$ is the decay constant of electron wave function.

The plasmon-assisted current is proportional to the matrix element of electric potential energy $e\varphi_{\pm}$. Its evaluation generally requires the knowledge of electron wave function inside the barrier layer. This evaluation can be, however, simplified in the dipole approximation. We write the potential distribution in the plasmon mode as $\varphi_{q}(z)={\bar{\varphi}}+(\varphi_{+}-\varphi_{-})z/d$, thus the potential matrix element becomes:

[TABLE]

It appears that the coordinate matrix element $z_{\pm}$ and the tunnel splitting $\Omega$ are bound by a simple relation. We consider two methods for calculation of current between states $\ket{+}$ and $\ket{-}$. From one hand, this can be expressed through velocity operator in the transverse direction:

[TABLE]

From the other hand, it can be found by evaluating the derivative of particle number in the state $\ket{+}$

[TABLE]

Comparing these two expressions, we find

[TABLE]

Appendix D Comparison of elastic and inelastic currents

The elastic current can be evaluated by considering the tunneling matrix elements in Hamiltonian (42) as small perturbations. This leads to the following formula Brey (2014)

[TABLE]

here $A_{ls}({\bf k},E)=-2{\rm Im}G^{R}_{ls}({\bf k},E)$ is the spectral function in the $l$-th layer and $s$-th band, and $G^{R}_{ls}({\bf k},E)$ is the retarded Green’s function in graphene in the band representation:

[TABLE]

$\Sigma^{ls}({\bf k},E)$ is the electron self energy. In the simplest approximation, the self energy can be treated as a constant $\Sigma^{ls}({\bf k},E)\approx\gamma$. In this case, one can approximate

[TABLE]

where $\delta_{\gamma}(x)=(\gamma/\pi)/(\gamma^{2}+x^{2})$ is the ”broadened” delta-function. If, in addition, the layers are aligned ($\Delta{\bf q}\rightarrow 0$), the integration is performed trivially yielding the electron ($n$) and hole ($p$) densities in the layers

[TABLE]

Appendix E Radiative decay of plasmon modes

The present section is aimed at the estimate of plasmon-to-photon conversion efficiency. The plasmon emitted upon tunneling can decay either radiatively (with the rate $\gamma_{rad}$) or it can be re-absorbed due to the Drude or intraband absorption in a single layer (the corresponding rate is $\gamma_{abs}$). Considering these competing channels of plasmon decay, we can estimate the plasmon-to-photon conversion probability as $\gamma_{rad}/(\gamma_{rad}+\gamma_{abs})$.

The radiative decay rate of plasmon $\gamma_{\rm rad}$ can be estimated as

[TABLE]

where $P_{\rm rad}=(4\omega^{4}/3c^{3})|{\bf d}_{\omega}|^{2}$ is the power of dipole radiation. We consider a double graphene layer sample of length $L$ and width $W$, so that $L$ corresponds to the fundamental plasmon mode $qL=\pi$. As example, we consider an acoustic plasmon mode with linear dispersion $\omega=sq$. Evaluating the average dipole moment, we find the radiated power

[TABLE]

The mode energy is

[TABLE]

where the first term is due to field and the second one is due to the particle motion. Further estimates will be done in the quasi-classical limit, $qv\ll\varepsilon_{F}$, $\omega\ll\varepsilon_{F}$. In this approximation, the in-plane conductivity is essentially intraband:

[TABLE]

The solution for acoustic plasmon dispersion law with conductivity (57) leads to $\omega=sq$ with the velocity

[TABLE]

where $K=4\alpha_{c}k_{F}d$, $\alpha_{c}=e^{2}/\kappa\hbar v$ is the coupling constant and $k_{F}$ is the Fermi wave vector. Evaluating the frequency derivative of conductivity (57) at the dispersion curve $\omega=sq$, we arrive at the following expression for mode energy in the long-wavelength limit

[TABLE]

Combining Eqs. (55) and (59), we find the radiative decay rate

[TABLE]

We note that the radiative decay rate of plasmons is generally small due to two reasons: (1) smallness of the dipole moment of double layer structure (the factor of $d$ in the numerator) (2) smallness of plasmon velocity compared to the velocity of light (the $c^{3}$-term in the denominator).

Appendix F Antenna coupling of plasmons

The plasmon-to-photon conversion efficiency can be markedly increased if the double layer device is loaded with an antenna. To model the plasmon decay in this situation, we consider the two graphene layers connected via the radiative resistance $Z_{rad}$. We shall solve the dispersion equation for plasmons in this structure and find their decay rate due to radiation. An electric potential is sought for as a superposition of forward and backward optical and acoustic waves

[TABLE]

where $q_{-}$ and $q_{+}$ are the wave vectors of acoustic and optical modes. The boundary conditions for the schematic in Fig. 6 are

[TABLE]

where $I$ is the current induced in external circuit. Solving Eq. (61) with boundary conditions (62), we obtain the following dispersion relation

[TABLE]

It is possible to estimate the solutions analytically in the limit $q_{+}/q_{-}\ll 1$, i.e. the wavelength of optical plasmon much exceeds that of acoustic one. This is generally fulfilled as the acoustic plasmons have linear dispersion while the optical have a square-root one (see also Fig. 2). In this limit, the general dispersion equation (64) is decoupled into two, neither depending on $q_{+}$:

[TABLE]

Only the solutions of the second equations are affected by the radiative decay. It is convenient to rewrite it introducing the dimensionless frequency $u=q_{-}L/2=\omega L/2s$, and the dimensionless radiative resistance

[TABLE]

where $\omega_{pl}=\pi s/L$. The dispersion equation becomes

[TABLE]

A general feature of its solutions is that the imaginary part of frequency has an extermum as a function of $\tilde{Z}$. This is illustrated in Fig. (7) which shows the decay rate of five lowest plasmon modes vs. radiative resistance $Z_{rad}$ calculated with numerical solution of Eq. (66).

There exists an optimal value of antenna resistance providing the maximum radiative decay rate. Decoupling the solutions into real and imaginary parts, $u=u^{\prime}+iu^{\prime\prime}$, we find that the maximum of decay rate is achieved if

[TABLE]

For the two lowest modes we have obtained $u^{\prime}=3.4$, $u^{\prime\prime}=-0.14$ and $u^{\prime}=6.4$, $u^{\prime\prime}=-0.08$, respectively. We note that the condition of maximum radiative decay ${\tilde{Z}}=u^{\prime}\sim 1$ represents the matching of antenna impedance and impedance of graphene layer at the resonant plasmon frequency. For the lowest mode, the maximum decay rate is $\gamma_{rad}\approx 0.04\omega$. This greatly exceeds the decay rate due to dipole radiation into free space.

Finally, we estimate the rate of plasmon absorption due to the Drude loss. It is worthwhile noting that $\gamma_{abs}$ is not just the inverse of electron momentum relaxation rate due to the non-negligible spatial dispersion of conductivity Principi et al. (2014). To account for the spatial dispersion and electron scattering simultaneously, one can solve the kinetic equation for electrons with particle-conserving collision integral Svintsov et al. (2016). This leads to the modification of collisionless conductivity (57) according to

[TABLE]

where $x=(\omega+i\nu)/qv$, and $\nu$ is the electron scattering rate. Considering electron-phonon and electron-impurity collisions as the dominant scattering sources, we evaluate $\nu$ as Vasko and Ryzhii (2007); Adam et al. (2009)

[TABLE]

Here $D\approx 30$ eV is the deformation potential in graphene, $\rho=7.6\times 10^{-7}$ kg/m2 is its mass density, and $c_{s}=2\times 10^{4}$ m/s is the sound velocity, $N_{i}$ is the impurity density, and $J(\alpha_{c})$ is the dimensionless integral

[TABLE]

The results of scattering rate and conversion efficiency calculation are shown in Figs. (8) and (9), respectively. An increase in scattering rate with reducing the Fermi energy in Fig. (8) is due to the impurity scattering contribution to plasmon damping which scales as $\varepsilon_{F}^{-1}$. At lower impurity density, the scattering will be dominated by phonons and an increase in Fermi energy would increase the scattering rates.