Exciton-assisted electron tunnelling in van der Waals heterostructures

The control of elastic and inelastic electron tunnelling relies on materials with well-defined interfaces. Two-dimensional van der Waals materials are an excellent platform for such studies. Signatures of acoustic phonons and defect states have been observed in current-to-voltage measurements. These features can be explained by direct electron–phonon or electron–defect interactions. Here we use a tunnelling process that involves excitons in transition metal dichalcogenides (TMDs). We study tunnel junctions consisting of graphene and gold electrodes separated by hexagonal boron nitride with an adjacent TMD monolayer and observe prominent resonant features in current-to-voltage measurements appearing at bias voltages that correspond to TMD exciton energies. By placing the TMD outside of the tunnelling pathway, we demonstrate that this tunnelling process does not require any charge injection into the TMD. The appearance of such optical modes in electrical transport introduces additional functionality towards van der Waals material–based optoelectronic devices.


S1. Device details and I-V characteristics
Geometrical details of the four devices discussed in the main text are provided in Table S1.The device size is determined by the overlap region of graphene, the hBN tunnel barrier, and the Au electrode.The TMD coverage is defined as the percentage of the device area under a TMD monolayer.The thickness of the hBN tunnel barrier is ∼3-4 nm, as estimated from the combination of the optical contrast and AFM measurements.However, the exact number of hBN atomic layers is difficult to determine with our AFM setup.The measured I-V characteristics of the devices are plotted in Fig. S1.The exciton-assisted tunneling feature is less visible in the WSe 2 device than in the other three devices due to a smaller TMD coverage.Nevertheless, the exciton signature is clearly observed in the dI/dV spectrum of this device, as shown in the main text.The current of the MoSe 2 device at a given bias voltage is more than two orders of magnitude larger than that of the other three devices due to a thinner hBN tunnel barrier.An exciton-assisted light emission can be measured for this device and is shown in the main text, while no light emission is detected for the other three devices.Assuming a similar electron-to-photon efficiency for all of our TMD devices, the number of emitted photons is at least two orders of magnitude smaller in the other three devices than in the MoSe 2 device at a given bias voltage, and thus, the electroluminescence signal lies below the detection limit of our setup.

S2. Data for additional devices
To further corroborate our findings, we fabricated two more devices, from which the results shown in the main text are reproduced.The I-V and dI/dV curves of a second WS 2 device and a second WSe 2 device are displayed in Figs.S2 and S3, respectively.Resonance features can be clearly observed in the dI/dV spectra, at ∼2.05 V for WS 2 and ∼1.64 V for WSe 2 , which are essentially the same results as we find in the devices discussed in the main text.The reproducibility of the results further supports our exciton-assisted tunneling interpretation.

S3. Discrimination between actual features and measurement instabilities
In this section, we are presenting dI/dV curves for the four TMD devices (WS 2 , MoS 2 , WSe 2 , and MoSe 2 ) measured at 10 K after (Fig. S6a-d) and before (Fig. S6e-h) filtering noise by means of a Savitzky-Golay filter.We conclude that the features that appear at high voltages in the dI/dV curve of the WSe 2 device in Fig. S6c can be attributed to measurement instabilities and increased noise.This becomes apparent in the unfiltered version of the measurements in Fig. S6g, where increased noise and fluctuations are observed.Finally, the features appearing close to the exciton resonances for the rest of the devices (WS 2 , MoS 2 , and MoSe 2 ) are still observed in the unfiltered data.To further support this interpretation, we show in Fig. S7a the MoS 2 dI/dV curve at both 50 K and 10 K, with features appearing consistently in both of them, suggesting that they are physical.This is in contrast to the dI/dV curve of the MoS 2 device in Fig. S7b, where the high-energy features vary for different measurements at 50 K and 10 K, suggesting that they are due to instabilities.FIG.S6.Measured dI/dV curves at 10 K for the various TMDs studied in this work.Curves in panels (a-d) are treated with a Savitzky-Golay filter, whereas the unfiltered data are presented in (e-h).Features observed at V b >2.2 V appear to be artificial due to the excess noise at those voltages, as shown in (g).In (a), some features appear consistently in both measurements close to the exciton resonance, whereas in (b), features at higher energies fluctuate, suggesting that they could originate in measurement instabilities.

S4. Graphene/hBN/Au tunnel junction: comparison between theory and experiment
In the tunnel junctions presented in the main text, we explore a range of voltages up to 3 V.To further our understanding on whether we measure contributions from Fowler-Nordheim (F-N) tunneling or thermionic effects, we use a simple tunneling model with an infinitely high tunneling barrier to compare to our measurements.The phonon-induced tunneling current between graphene (Gr) and Au is calculated as [1] where ℏΩ is the phonon energy, T is the transmission matrix element (assumed to be constant in accordance with the DFT calculations presented in the supplementary material of Ref. [1]), ρ Au is the DOS of Au (also assumed to be constant for V b > 0 according to the DFT calculations in Fig. S9), and ρ Gr is the DFT-calculated DOS of graphene taken from Ref. [2] and shifted in energy according to the level of graphene doping (∼ 0.3 eV [1]).By taking the derivative of the calculated current as a function of voltage, we obtain the dI/dV curve shown in Fig. S8 (red dashed curve).We scale the calculation to fit the low-voltage region of the measured dI/dV curve.This is the region in which contributions of F-N or thermionic effects should be limited.Our comparison reveals that the calculated dI/dV curve closely follows our measurement up to 2 V, with significant deviations only above 2.5 V.These high-voltage contributions, which may be the result of F-N or thermionic effects, seem to become relevant at voltages outside the range in which we explore exciton energies (1.5 − 2.25 V).

S5. Excitonic light emission from a WSe2 device with a thinner hBN tunnel barrier
Here, we investigate another WSe 2 device that has a thinner hBN tunnel barrier (∼ 2.3 nm).The measured I-V profile is plotted in Fig. S4, showing a much larger tunnel current than in the devices with thicker tunnel barriers discussed above.In the dI/dV spectrum, the phonon-related feature shows up near zero bias voltage [3], but no clear exciton signatures can be observed.We attribute this result to the fact that the background arising from the larger total current overshadows any signature from the much smaller exciton-assisted tunneling current [4].However, exciton-assisted light emission is observed from this device, as shown in Fig. S5.The peak at ∼1.63 eV in the electroluminescence spectra matches well the WSe 2 A exciton ground state [5], thus pointing to a coupling of the tunneling electrons to excitons in the WSe 2 monolayer.The bump feature around 1.4 eV in the spectra can be attributed to the broad light emission resulting from the coupling of tunneling electrons to optical modes such as surface-plasmon polaritons in the device [1].

A. Preliminary considerations regarding energy-momentum conservation
We first note that the electronic density of states (DOS) of gold near the Fermi level is dominated by the 6sp-band, which is well-described by the free-electron gas (FEG) model up to energies ∼ 3 eV above the Fermi level (see Fig. S9).The contribution of the 5d band shows up as a sharp increase in the DOS below ∼ −2 eV, which explains the sharp jump in the tunneling current observed in the experiment at negative biases exceeding that value.For the remainder of this discussion, we focus however on the smooth region dominated by the sp band, which we describe in the FEG model.
In a simplified picture, the surface-momentum-projected sp band of gold defines a circle of radius given by the Fermi wave vector k Au F , which we compare to the conduction band of graphene in Fig. S10a.This plot shows that graphene electrons near its Fermi level are situated far from the Γ point, where gold conduction electrons are expected to extend further beyond the metal surface.A band projection on energy and parallel momentum along the ΓK g direction of n-type doped graphene (Fig. S10b) also illustrates how the electrons of this material need to bridge a large momentum mismatch to tunnel to gold states near the Γ point.This mismatch is smaller for transitions to gold states with high parallel wave vectors, but then, a huge jump in out-of-plane energy is encountered.In this respect, the landscape of the conduction band bottom (CBB) in the hexagonal boron nitride (hBN) layer [6] separating graphene from gold in our samples could play a role, whose analysis would require a realistic description of the electronic states in the system beyond the scope of the present theoretical model.However, a quantitative argument can be constructed upon , estimated from a one-dimensional potential barrier determined by the hBN gap at the Γ point [6] Eg ≈ 8.9 eV and the effective mass along the out-of-plane ΓA direction m * hBN ≈ 0.63 me.In (a) and (b), we consider two possible channels for graphene-to-gold tunneling: (1) one-step direct inelastic transitions to near-threshold gold states around the surface projected Γ point (coinciding with the Γ point of graphene) assisted by the creation of indirect excitons in a neighboring TMD monolayer (dark blue arrows); and (2) two-step quasi-elastic tunneling from graphene to gold states near the surface-projected Γ point assisted by large-momentum phonon excitation (purple arrows), followed by quasi-vertical inelastic decay to states in the vicinity of the Fermi level of gold mediated by direct excitons (green arrows).
examination of the exponential attenuation factor displayed by the wave functions of gold conduction electrons across a distance of 3 nm outside the metal in the interface with hBN.The corresponding plot (Fig. S10c) clearly indicates that tunneling should be exponentially attenuated for energies and wave vectors deviating from the conditions of both quasi-elastic tunneling (E − E Au F ≲ eV b = 3 eV in the figure) and electron momentum in the surface-projected gold Γ point.
Based on these arguments, we expect tunneling in the TMD/graphene/hBN/Au sandwich for positive bias of the graphene (as considered in the experiment) to involve a momentum transfer to graphene electrons such that they can transition to gold states near the Γ point (since graphene does not have any states near such point).The required wave vector transfer, k ∥ ∼ ΓK g = 4π/(3 √ 3a CC ) ≈ 17.02 nm −1 , dictated by the graphene lattice and the carbon-carbon bond distance a CC ≈ 0.1421 nm, is too high to be supplied by the creation of low-momentum optical modes such as direct excitons.For example, the nonlocal conductivities of the four studied TMDs show a strong attenuation for k ∥ ≳ 0.5 nm −1 (see Fig. S18 in Sec.S9).An additional source of momentum is thus required, for which we postulate two possible channels that can assist graphene-to-gold tunneling: • Direct one-step transitions assisted by indirect TMD excitons.Indirect excitons in the TMD monolayer can be formed by involving an electron and a hole of very different wave vectors.Upon examination of the possible indirect excitons connecting flat points in the band structure of TMD monolayers (i.e., those capable of producing Van Hove singularities), we find KQ TMD excitons (consisting of one hole at the K point and one electron at the Q point of the TMD) as good candidates to contribute to the observed tunneling current because their separation in wave vector is close to the ΓK g distance in graphene (see Table S2).In addition, our calculated energies for KQ TMD excitons are in relatively good agreement with the observed peaks in our conductance measurements, as shown in Table S2, where we also include the corresponding KQ TMD gap energies (again not too far from the observed peaks), as well as information on the A and B excitons.For reference, we show the first Brillouin zone (1BZ) of MoSe 2 (pink) in the following scheme, superimposed on the 1BZ of graphene and the surface-projected gold Fermi surface taken from Fig. S10a: In Table S2, we show values of the KQ TMD and K ′ Q TMD distances, indicated by double arrows in this scheme.
Although the relative azimuthal orientation between TMD and graphene lattices is not well defined in our samples, the set of KQ TMD distances spanned by the combinations of inequivalent K g and Q points is dense, and some of them lie close to the required graphene ΓK g momentum transfer, particularly considering that KQ TMD excitons span a wave vector distribution of finite extension in the ∼ 1 nm −1 range (see Fig. S19 in Sec.S9).
• Two-step processes assisted by phonon excitation and direct TMD excitons.Phonon excitation can bridge the graphene ΓK g momentum mismatch, as both hBN and graphene have relatively close lattice parameters, and in addition, the electronic and phononic Brillouin zones are identical within each of them.The involved phonon energies are ∼ 70 meV and ∼ 120 − 150 meV [7,8], much smaller than the exciton energies considered in this work.Apart from indirect excitons (see previous point), we cannot find a single source of momentum and energy that can simultaneously break the energy-momentum mismatch between graphene states and gold conduction electrons near the Fermi level at the Γ point.However, a two-step process can be invoked, involving a first quasi-elastic transition from graphene to gold near the Γ point, assisted by the excitation of a high-momentum phonon in graphene (or hBN, although proximity favors graphene phonons); and a subsequent inelastic decay from higher-to lower-energy states in gold (both of them near the Γ point), taking place by investing the transition energy in the creation of an exciton in the TMD monolayer placed on top of graphene.Inverting the order of these two steps would involve inelastic decay from graphene to an intermediate unoccupied gold state of large parallel momentum, and thus, this possibility should produce a negligible contribution because of the marginal spill out of such intermediate state near the Fermi level of gold (see Fig. S10c).For negative bias, with p-doped graphene, electrons could still tunnel from gold to graphene following a two-step process, but the intermediate and final states (both in graphene) have a small local DOS compared with gold, and thus, two-step inelastic tunneling should occur with lower probability in this direction.
We conclude these preliminary considerations by estimating the dependence on the graphene-gold bias voltage V b of the shift in the graphene Fermi energy E g F relative to the Dirac point energy E D (i.e., relative to the graphene Fermi energy at zero bias).Throughout this work, we refer all energies to the CBB of gold, so we write the noted shift as In addition, we assume hBN to fill the regions both below and above the TMD+graphene heterostructure, and ignore the effect of the TMD monolayer on graphene doping.A variational minimization of the total energy in the system shows that the graphene doping density n (electrons per unit area) is determined by the capacitor formed with gold as n = (ϵ hBN /4πe d)(V b − ∆E g F /e), where ϵ hBN ≈ 3.76 is the out-of-plane DC permittivity of hBN, d is the hBN spacing layer thickness separating graphene from gold, and the effective bias voltage needs to be   S2.Direct and indirect excitons in TMD monolayers.We present a compilation of energies for the measured KQ TMD gaps (from photoluminescence [9,16,18] and I-V [14] studies), the conductance peaks in our experiment, luminescence measurements from different sources, and our calculated A and B excitons.We show data for the four different TMD monolayer materials under consideration.All energies are in eV.We also provide wave vector distances between singular points of the 1BZ (four rightmost columns).The Q point is defined by an intermediate conduction band minimum along the ΓK TMD line.
We consider two of the four distances between K and Q points in the 1BZ of each TMD, whose values lie close to the ΓKg distance in graphene (17.02 nm −1 ).
reduced by ∆E g F /e to account for the change in energy needed to add/subtract every new electron to/from graphene as its Fermi energy is varied.Using the additional relation [21] πn for positive bias, where v g F ≈ 10 6 m/s is the graphene Fermi velocity, we finally obtain the Fermi energy shift as From this expression, we have, for example, ∆E g F ≈ 0.49 eV for d = 3 nm and eV b = 3 eV (≫ E 0 ≈ 0.09 eV).Incidentally, this value is slightly shifted with respect to the large-bias limit FIG. S11.Electron potential landscape.We model the out-of-plane electron dynamics through a one-dimensional potential energy U (z), here sketched for a graphene-gold positive bias potential energy eV b and defined by the conduction band bottom (CBB) of gold at z < 0, as well as the potential energies U0 and U1 of the CBB of hBN at z = 0 and z = d (the graphene-gold spacing), respectively.The graphene Fermi level E g F is separated from the gold Fermi level E Au F by the bias potential energy 45 eV as half of the hBN gap at the Γ point [6].

B. Vertical potential landscape, graphene electrons, and electron Green function
We assume the tunneling electron to evolve in the potential energy landscape sketched in Fig. S11, as described by the Hamiltonian where a uniform isotropic electron mass is assumed for simplicity.We note that the strong variation of the electron mass and energy gap within the in-plane 1BZ of hBN [6] can produce substantial corrections to this model, which deserve a future investigation including a more detailed description of the electron states involved in the tunneling process.In addition, we approximate the wave function of each initial graphene electron i as a separable state relative to the CBB of gold and wave vector Q i relative to one of the K g points of momentum K g (see Figs. S10 and S11).Here, s i = ±1 refers to the upper (+) and lower (−) graphene Dirac cones, K g and Q i and both 2D vectors, the in-plane position coordinates are R = (x, y), and we introduce the graphene area A for normalization.We also note that E g F is referred to the CBB of gold.Incidentally, the two inequivalent Dirac cones within the 1BZ should contribute equally to the tunneling current, and therefore, we absorb valley and spin degeneracies in an overall factor of 4.
The evolution of the tunneling electron in the potential landscape of Fig. S11 can be conveniently described in terms of the three-dimensional (3D) electron Green function Due to the in-plane translational symmetry of the system, the Green function depends on R and R ′ only through the difference R − R ′ , so it can be written as ), which we correct in Eq. (S3) by subtracting the in-plane energy ℏ 2 Q 2 /2m e from ℏε.
In the calculations presented below, we only need the Green function evaluated at G 0 (0, z, ε) and G 0 (z, d, ε), as well as the derivative of the former with respect to the first argument.It is then useful to point out that the spatial arguments of these functions are interchangeable in virtue of the reciprocity property G 0 (z, z ′ , ε) = G 0 (z ′ , z, ε).In addition, we note that the required functions are solutions of (H 0 − ℏε)G 0 (z, z ′ , ε) = −δ(z − z ′ ) for z ′ = 0 and z ′ = d.This allows us to express them in terms of analytical functions within each of the three distinct regions of the potential U (z) in Fig. S11: where we consider outgoing waves with k = 2m e ε/ℏ and κ = 2m e (U 1 /ℏ − ε)/ℏ, which are suitable solutions decaying at |z| → ∞ for a source placed at The functions in the intermediate region can be written as in terms of the Airy functions Ai and Bi [22] with an argument θ = (2m e eV b d 2 /ℏ 2 ) 1/3 z/d + (U 0 − ℏε)/eV b .Finally, the coefficients in Eq. (S4) are obtained from the conditions In particular, Eqs.(S5a) and (S5b) guarantee continuity at z = 0 and z = d, respectively, whereas Eqs.(S5c) and (S5d) relate to the jump in the derivative (2m e /ℏ 2 ) produced by the inhomogeneous term −δ(z − z ′ ).More precisely, we need to set n 0 = 1 and n d = 0 for z ′ = 0; and n 0 = 0 and n d = −1 for z ′ = d.The required Green function values are then given by Eq. (S4) with coefficients A, B, C, and D determined by solving the 4 × 4 system of equations (S5).Likewise, we evaluate − by writing a solution similar to Eq. (S4), with coefficients determined by the equations This solution if obtained by first writing G 0 (z ′ , z, ε) (with z ′ < 0 regarded as a parameter) in terms of z-dependent outgoing waves within the z < z ′ and z > d regions, as well as Airy functions for 0 ≤ z ≤ d and e ±ikz waves for z ′ ≤ z < 0. We then take the derivative of the resulting coefficients with respect to z ′ and calculate the z ′ → 0 − limit (i.e., approaching z ′ = 0 from the left), which leads to Eqs. (S6).

C. Inelastic electron tunneling assisted by excitations in the heterostructure
The tunneling process can be assisted by the creation of an excitation in the heterostructure, which absorbs the energy and momentum change undergone by the transferred electron.In this section, we derive a general expression for the inelastic tunneling current that encapsulates such excitations n of energies ℏω n in the dielectric response function.Although we consider electron tunneling from graphene to gold in this work, the formalism presented in this section can be applied to any combination of two conductive materials separated by an insulating layer and subject to a relative bias.Starting from an electron in an initial graphene state i of wave function φ i (r) and energy ℏε i , the final electron wave function component of energy ℏ(ε i − ω n ) associated with an additional excitation n is given by within first-order perturbation in the interaction Hamiltonian H 1 of matrix elements sandwiched by the ground and excited states of the heterostructure, |0⟩ and |n⟩, respectively.Here, ρ n0 (r) = ⟨n| ρ(r) |0⟩ is the matrix element of the charge density operator ρ(r) = j q j δ(r − r j ), incorporating all electrons and nuclei of charges q j at positions r j in the system.The Hamiltonian H 1 thus describes the Coulomb interaction between such charges and the tunneling electron in a completely general fashion.The details of the structure are additionally captured by the unperturbed Green function G 0 (r, r ′ , ε), which remains general in this section, but we then specify it for the heterostructure under consideration by using the methods described in Sec.S6 B.
From the wave functions given by Eq. (S7), we obtain the tunneling current density (charge moving downwards and entering the gold surface per unit time and unit area) as where we sum over all initial occupied states i and material excitations n that lead to final unoccupied gold states of energy ℏ(ε i − ω n ), as imposed by using the Fermi-Dirac distribution f T at temperature T .The gold surface is placed at z = 0, so we evaluate the current right below that plane at z = 0 − in Eq. (S9).From Eqs. (S7)-(S9), we find Following a method similar to the derivation of the electron energy-loss probability in electron microscopy (see Ref. [23]), this expression can be written in terms of a response function by exploiting the identities [24] Im χ(r, r ′ , ω) = − π where χ(r, r ′ , ω) is the nonlocal susceptibility and W ind (r, r ′ , ω) is the induced part of the screened interaction.The latter is defined as the potential created at a position r by a time-dependent unit charge of magnitude e −iωt placed at r ′ (see Secs.S7 and S8 for explicit calculations of this quantity in the planar heterostructures under consideration).Using Eqs.(S11) to manipulate Eq. (S10), we obtain where we have substituted the induced part of the screened interaction by the full interaction because

D. Indirect exciton-assisted one-step tunneling
We now apply Eq. (S12) to discuss one-step transitions from graphene states i assisted by the creation of indirect excitons in the TMD layer.In what follows, the sum over i is understood to be restricted to occupied graphene states and we set T = 0, such that 1 in expressed in terms of the step function Θ and the gold Fermi energy E Au F .Accordingly, we rewrite Eq. (S12) as Momentum mismatch is then considered to be bridged by large wave vectors k ∥ in the screened interaction, contributed by indirect excitons in the TMD layer.
Assuming isotropy and translational invariance in the in-plane response, one can write the screened interaction as in terms of momentum components W (k ∥ , z, z ′ , ω).However, for the large values of k ∥ under consideration, the in-plane atomic lattice can play a substantial role, imprinting a lattice periodicity on the screened interaction, which can be written as This expression, which generalizes Eq. (S14) by considering W G,G ′ (k ∥ , z, z ′ , ω) components labeled by TMD reciprocal lattice vectors G and G ′ and restricting the k ∥ integral to the 1BZ, is general to describe the linear optical response of a TMD layer including atomic periodicity.In our analysis, we neglect lattice contributions to the Green function as a less relevant effect than in the excitonic optical response of the TMD material.Upon insertion of Eqs.(S3) and (S15) into Eq.(S13), noticing the plane-wave dependence of the wave functions on in-plane coordinates, we can readily carry out the integrals over R ′ and R ′′ .Further averaging the current over in-plane positions R, we obtain where K g is the wave vector at the graphene K g point and we find contributions from different Brillouin zones labeled by G.We note that the average over R has eliminated nondiagonal W GG ′ components (G ̸ = G ′ ) from this expression.Incidentally, the rightmost argument of the function G 0 (z, z ′ , ε) in Eq. (S16) refers to the out-of-plane electron energy ℏε relative to the gold CBB, and consequently, the sums over i and G need to be restricted to only yield positive values of such argument.We proceed by using the fact that the initial graphene states i are spatially localized near the z = d plane [i.e., the out-of-plane spatial extension of φ i (z) (∼ ±1 Å [25]) is small compared with the tunneling distance d ∼ 2 − 3 nm], and therefore, Eq. (S16) can be approximated by setting z ′ = z ′′ = d.In addition, the initial wave vectors Q i relative to the graphene K g point (see Sec. S6 B) are small compared with K g , and their energies ℏε i close to E g F , so we approximate The sum over i is then trivially yielding an overall constant.Putting these elements together, we transform Eq. (S16) into Finally, we argue that G 0 (0, d, ε) is a strongly decreasing function of ε, and therefore, the lowest values of the argument ε should contribute maximally to the tunneling current.Such values are encountered at ℏω = eV b in Eq. (S17).Using an analogous argument, the region near the point k ∥ = G − K g for each reciprocal lattice vector G should produce a dominant contribution to the k ∥ integral.In addition, because the relative orientation of the TMD and graphene lattices is undefined, we average over the azimuthal angle of K g , which we denote as φ Kg .By doing so, the tunneling current is found to depend on bias voltage V b roughly as The 1BZ integration region of Eq. (S18) is represented in k ∥ space in Fig. S12 for MoSe 2 , while similar results are found for the other TMDs.The six nearest-neighbor TMD Brillouin zones are indicated by hexagons in this figure, each of them centered around its corresponding G vector.We also show the circles defined by the respective conditions |G − k ∥ | = K g .It is clear that only the six smallest non-vanishing TMD reciprocal lattice vectors G can satisfy such conditions for k ∥ within the 1BZ.Furthermore, as we have averaged over the orientation of ΓK g in graphene, and in virtue of symmetry, all of those G's should make equal contributions.Consequently, we only need to consider a single , where a TMD is the in-plane lattice constant of the TMD.This allows us to rewrite Eq. (S18) as where , and we recall that K g ≈ 17.02 nm −1 is the ΓK g distance in graphene.The integration contour in Eq. (S19) is indicated by a thick circular segment inside the 1BZ of Fig. S12.The optical conductivity of the TMD, illustrated for MoSe 2 at ℏω = 2.1 eV in Fig. S12 through a color plot of Re{σ TMD,00 (k ∥ , ω)} (i.e., with G = G ′ = 0, see Secs.S8 and S9), shows that the response is dominated by localized momentum regions corresponding to the KQ TMD excitons (see Table S2), inside of which the wave vector component k y does not vary significantly, so we approximate the current as for the sake of the discussion of one-step tunneling processes in the main text.[Nevertheless, the results presented in Fig. S13 and Fig. 4c in the main text are obtained by evaluating Eq. (S19).]We remind that k ∥ is understood to be determined by the condition |G 0 x − k ∥ | = K g in Eq. (S20).We evaluate Eq. (S19) from the diagonal elements of the screened interaction W GG (k ∥ , d, d, ω), which are in turn calculated as explained in Sec.S8 below.A result for the I-V curve is presented in Fig. S13 for MoSe 2 with a hBN tunnel barrier distance d = 3 nm.Spectral features are clearly visible in the 1.9-2.2eV range, corresponding to indirect excitons that are also observed in our first-principles calculations of the TMD conductivity (see Fig. S19 in Sec.S9).The spectral range in which these features show up is in good correspondence with experiments, although the detailed spectral shape presents differences that are possibly connected to the approximations adopted within the present theoretical model.

E. Two-step tunneling assisted by phonon and exciton creation
We now analyze two-step processes involving phonon excitation followed by the creation of nearly direct excitons as an additional channel contributing to electron tunneling from graphene to gold and also producing an excitonic signature in the dI/dV curves.S13.Calculated I-V curve for indirect exciton-assisted one-step tunneling.We plot the current obtained by using Eq.(S20) for MoSe2 (see Fig. 4 in the main text for other TMDs) with a hBN tunnel barrier distance d = 3 nm.The inset shows the corresponding dI/dV curve.Incidentally, the current drops beyond ∼ 2.5 eV because of the reduced number of bands used to calculate the response associated with indirect excitons (see Sec. S9).

Quasi-elastic phonon-assisted transitions
The phonon energies at the K point are large compared with k B T at room temperature T , and therefore, we can safely neglect phonon absorption processes.Starting from a graphene state i (see Sec. S6 B), within first-order perturbation theory, the excitation of a phonon p of wave vector K g + Q p (i.e., with Q i relative to the graphene K g point) and frequency ω p contributes with an electron wave function φip is expressed in terms of the electron-phonon coupling Hamiltonian H ph p (r ′ ) and the 3D electron Green function see Sec.S6 B for details of this function and the definitions of φ i (z), Q i , and A].Separating the in-and out-of-plane dependence of the phonon Hamiltonian as and noticing again that φ i (z) is tightly localized around z = d, we can approximate this expression as where the constant C p depends on the specific phonon mode p under consideration.In addition, as tunneling is expected to involve large phonon momentum transfers that place the electron near the Γ point, only a reduced region in the phonon bands contributes to the process, and therefore, we assume C 1 ≡ |C p | to be roughly independent of p.
The tunneling current per unit area traversing the hBN/Au interface is finally obtained as , which is evaluated at a position z = 0 − right inside the metal (see Fig. S11 and Sec.S6 C), with the step function restricting the sum over i and p such that only final unoccupied electron gold states above the Fermi level E Au F are included.

Phonon+exciton-assisted tunneling
In the two-step channel, after a phonon has been created, a direct exciton absorbs most of the energy lost by the tunneling electron.This second step is described by Eq. (S12) in Sec.S6 C, with φ i substituted by the phononscattered wave function φip given in Eq. (S22).Following a similar procedure as in Sec.S6 D, we insert Eqs. ( S3) and (S14) into Eq.(S12) and carry out the integrals over R, R ′ , and R ′′ to find where we have assumed translational invariance in the screened interaction [Eq.(S14)] and defined We note that the intermediate phonon-assisted state φip (z) is generated from the graphene region (at z = d), and therefore, we expect it to be exposed to further inelastic interactions from that region, including TMD excitons.Rather than propagating it elastically (once the phonon has been created) to tunnel into the metal (i.e., like in Sec.S6 E 1), we are interested in the immediate interaction with the excitons.Therefore, proceeding as we did in Sec.S6 D, we approximate Eq. (S23) by setting z ′ = z ′′ = d and encapsulating the dependence on the intermediate states in an overall constant factor, so we write We further simplify the evaluation of Eq. (S24) by arguing that the response function dies out quickly as k ∥ increases to values beyond 1 nm −1 (see Sec. S7 and Fig. S18), which is small compared with both K g ≈ 17.02 nm −1 and k Au F ≈ 12 nm −1 , so we can dismiss the k ∥ dependence in front of Q p in the last argument of the Green function G 0 .Likewise, the wave vectors of graphene electrons relative to the Dirac point are also small and can be neglected in G 0 when compared with Q p .We thus approximate In contrast, the phonon wave vector Q p is large, so it has to be retained.We replace the sum over phonons p by C 2 ∞ 0 QdQ, where Q ≡ Q p and the constant C 2 accounts for further details of the electron-phonon coupling.Analogously, we replace the sum over graphene electrons i by C 3 s=±1 Qs 0 Q ′ dQ ′ , where Q ′ ≡ Q i and the constant C 3 incorporates further graphene band details.We also sum over s = +1 and s = −1, referring to the upper (with Q ′ running up to the graphene Fermi wave vector Q 1 = k g F ) and lower (with Q −1 = ∞) Dirac cones, respectively.In addition, we assume ω p to be independent of p.These considerations allow us to reduce Eq.(S24) to d, d, ω) where ℏε 0 ≡ E D = E g F − ∆E g F is the Dirac point energy relative to the CBB of gold (see Figs. S10 and S11), and a cutoff wave vector is introduced to guarantee that the final out-of-plane electron energy is above the CBB of gold.
We are interested in examining the tunneling current for bias potential energies close to the TMD exciton energies, which are large compared with the graphene electron energies relative to the Dirac point.Therefore, we further approximate Eq. (S25) by setting and treating s Q ′ dQ ′ as an additional overall constant.This leads to where we need to redefine Finally, we argue that G 0 (0, d, ε) is a strongly decreasing function of ε, and therefore, the lowest values of the argument ε should contribute maximally to the tunneling current.Such values are encountered at ℏω = eV b according to Eq. (S26), from which the phonon+exciton-assisted tunneling current is found to depend on bias voltage V b roughly as This result relates the phonon+exciton-assisted tunneling current to the screened interaction in the heterostructure.The latter is studied in detail in Sec.S7.In Fig. S14, we plot the I-V curves obtained from Eq. (S27) for d = 3 nm, with the screened interaction obtained by describing gold in the specular-reflection model, hBN through a local anisotropic permittivity, graphene in the random-phase approximation (RPA), and the TMD monolayer through a nonlocal surface conductivity (see Sec. S9).The tunneling current shows dominant features associated with C excitons, while A and B excitons emerge as weak shoulders at lower biases.
Upon examination of the k ∥ dependence of Im − W (k ∥ , d, d, eV b /ℏ − ω p ) (see Fig. S18), we conclude that the integral in Eq. (S27) yields a result similar to Im − W (0, d, d, eV b /ℏ − ω p ) (i.e., at k ∥ = 0), so we write where we insert subindices indicating that this part involves G = 0 (i.e., no umklapp processes), for the sake of the discussion of two-step tunneling in the main text.[Nevertheless, the results presented in Fig. S14 and Fig. 4d in the main text are obtained by evaluating Eq. (S27).]

S7. Screened interaction under the in-plane isotropic and translational-invariant approximations
Tunneling electrons couple to excitations in the structure through the screened interaction W (r, r ′ , ω).Because the distances involved are small compared with the optical wavelength at frequency ω, we can safely neglect retardation effects and treat the interaction in the electrostatic limit.The structure under consideration is sketched in Fig. S15a and consists of a layer of 2D materials separated by a distance d from a gold surface.The region in between the two of them is filled with hBN, and the medium above the 2D materials is also hBN.
The layer of 2D materials here considered consists of a TMD monolayer on top of graphene.We describe it in the zero-thickness approximation through a nonlocal wave-vector-and frequency-dependent surface conductivity σ(k ∥ , ω).This approximation simplifies the algebra considerably and yields reasonable results in comparison to finite-thickness layers [26].We thus set σ(k ∥ , ω) = σ TMD (k ∥ , ω) + σ g (k ∥ , ω) as the sum of the surface conductivities of the TMD and graphene monolayers.A nonlocal version of the former is calculated from first principles (see Sec. S9), whereas the nonlocal conductivity of graphene is modeled in the RPA [27].The structure is covered by an upper hBN layer.Gold is described through a surface-corrected nonlocal dielectric function ϵs(k ∥ , ω), whereas hBN is modeled through its anisotropic frequency-dependent permittivity [28] of in-plane and out-of-plane components ϵx(ω) and ϵz(ω), respectively.The gold and upper hBN media are assumed to have an infinite thickness, while the 2D material is treated in the zero-thickness limit.(b,c) The response of the system is described in terms of Fresnel transmission and reflection coefficients for each of the two interfaces.
The interfaces in our structure are described by first defining Fresnel's reflection coefficients, as shown in Fig. S15b,c.More precisely, only p-polarization coefficients survive in the electrostatic limit.For convenience, we define in-plane coordinates R = (x, y) and work in the space of in-plane wave vectors k ∥ = (k x , k y ), with electrostatic quantities such as the electric potential written as ϕ(r, ω) = (2π We start by considering the electrostatic potential created by a point charge of time-dependent magnitude e −iωt placed at the origin in bulk hBN.The momentum components of such potential are given by [2π/k ∥ ε(ω)] e −q|z| , where ε(ω) = ϵ x (ω)ϵ z (ω) and q = k ∥ ϵ x (ω)/ϵ z (ω), with the square root sign chosen to yield positive real parts.Now, with the 2D material layer placed at z = d, scattering of an incident potential component e −q(z−d) coming from the z < d region gives rise to a transmitted potential t 2D e −q(z−d) at z > d and an incident+reflected potential e −q(z−d) − r 2D e q(z−d) at z < d. (Incidentally, note that the reflection coefficient is preceded by a − sign in this expression for the potential [26].)The corresponding Fresnel coefficients are determined by imposing the continuity of the potential at z = d, as well as the condition that the jump in the normal electric displacement is 4π times the surface charge produced through the conductivity σ(k ∥ , ω).Writing such charge in terms of the current through the continuity equation, and this in turn as the conductivity times the in-plane electric field, we are left with a set of two linear equations from which we obtain ) for the Fresnel coefficients of the 2D layer, which generalize previous results [29] to deal with an anisotropic host medium.We note that the z dependence of the potential components inside hBN is simply given by e ∓q(z−d) for upward (− sign) and downward (+ sign) components (see Fig. S15b).
Given the small thickness of the hBN tunnel layer, nonlocal effects in the response of gold can play an important role, and therefore, we use the nonlocal bulk permittivity of this material ϵ Au (k, ω).As a reasonably accurate prescription, we express ϵ Au (k, ω) as the sum of a conduction-electron component, treated in the RPA, and a local contribution associated with interband transitions, as described in detail elsewhere [30].In addition, we adopt the specular-reflection model [23,30,31] to incorporate surface effects through a surface response function where Incidentally, this quantity becomes ϵ s (k ∥ , z, ω) = e −k ∥ |z| /ϵ Au (ω) in the local limit (i.e., when dismissing the k dependence of the gold permittivity).We also define ϵ s (k ∥ , ω) ≡ ϵ s (k ∥ , 0, ω), in terms of which we find the Fresnel coefficients for the hBN/Au interface by imposing the continuity of both the potential and the normal electric displacement.Incidentally, the latter involves ϵ Au • ϵ s (in operator notation), which reduces to e −k ∥ |z| in (k ∥ , z) space.With the hBN/Au interface placed at z = 0, the z dependence of the potential indicated by the arrows in Fig. S15c is given by e ∓qz for z > 0. However, the dependence inside gold (z < 0) is more complex [23]: it is partially encapsulated in a factor ϵ s (k ∥ , z, ω)/ϵ s (k ∥ , ω) multiplying the transmission coefficients; in addition, an involved dependence is encountered for internal reflection from the gold side (see more details in Ref. [23]).
Finally, the screened interaction can also be decomposed in parallel wave vectors as shown in Eq. (S14), where ) is in turn separated into direct and reflected components.The former corresponds to the direct Coulomb interaction when z and z ′ lie in the same medium (see Fig. S15a): where the part inside gold incorporates nonlocal bulk corrections.The remaining surface contribution is given by for z ≥ z ′ , while the reciprocity property W dir (k ∥ , z, z ′ , ω) = W dir (k ∥ , z ′ , z, ω) can be used to compute it for z < z ′ .
In Eq. (S31), we use the Fabry-Perot factor ∆ = 1 − r 2D r BN,Au e −2qd −1 to account for multiple round trips within the 2D/hBN/Au cavity.We have explicitly verified that the sum of Eqs.(S30) and (S31) is continuous in both z and z ′ across the z = 0 and z = d interfaces, while the normal displacement is continuous at z = 0 and undergoes a jump by 2 ε k ∥ /(1 − 1/r 2D ) = iωε 2 /πσ(k ∥ , ω) at z = d due to the presence of the surface conductivity associated with the 2D materials.

S8. Screened interaction including the atomic lattice periodicity
A more rigorous description of the screened interaction incorporating the in-plane atomic corrugation of the TMD layer may become necessary when large wave vectors are involved, such as in the current experiment (i.e., for momentum transfers ∼ ΓK g in graphene).This is the case of the one-step tunneling channel mediated by the creation of indirect excitons (see Sec. S6 D), as described by Eq. (S17), in which the screened interaction enters through the components W GG ′ (k ∥ , d, d, ω) [see Eq. (S15)], with the source and probing positions both lying on the z = d plane (i.e., where the TMD-graphene heterostructure is located, assuming again a zero-thickness limit in the description of its optical response).
We obtain W GG ′ (k ∥ , d, d, ω) directly from a Fabry-Perot (FP) description of the pathways followed by the electric potential generated by a time-dependent unit charge of magnitude e −iωt placed at r ′ with z ′ = d + taken right above the 2D material layer.Using the notation and methods of Sec.S7, the direct potential in an infinite hBN host medium can be written in a form analogous to Eq. (S15): where and we define q G = |k ∥ + G| ϵ x (ω)/ϵ z (ω) similar to q in Sec.
where ∆ is a tensor of inverse components . The reflection and transmission tensors of the 2D layer entering Eq. (S33) can be in turn obtained from the conductivity tensor σ GG ′ (k ∥ , ω) by applying the same boundary conditions as invoked in the derivation of Eqs.(S28), but taking into account the expansion of different quantities in reciprocal-lattice-vector components.A detailed analysis leads to the result where I GG ′ = δ GG ′ is the identity matrix and we introduce a tensor of components Obviously, the G = G ′ = 0 elements r 2D,00 and t 2D,00 reproduce Eqs.(S28).Finally, we express the conductivity tensor σ GG ′ (k ∥ , ω) in terms of the non-interacting susceptibility tensor χ 0 GG ′ (k ∥ , ω) by writing the charge ρ ind induced in the 2D layer in response to a total potential ϕ expressed in (r, ω) space as ρ ind = (iω) −1 ∇ • j ind = (i/ω)∇ • σ∇ * ϕ, where j ind is the induced current, given in turn by the conductivity times the electric field −∇ϕ.Here, the dot stands for the scalar product of two 2D vectors and the star indicates a functional linear operation.By definition, the non-interacting susceptibility relates the total potential and the induced charge as ρ ind = χ 0 * ϕ, so we have the relation χ 0 = (i/ω)∇ • σ.Projecting on k ∥ (within the 1BZ) and G components, we obtain Incidentally, using Eq.(S36), we can write the matrix in Eq. (S35) as As explained in Sec.S7, the conductivity of the TMD-graphene 2D heterostructure is taken as the sum σ(k ∥ , ω) = σ TMD (k ∥ , ω)+σ g (k ∥ , ω).We apply Eq. (S36) to obtain the conductivity of the TMD monolayer from the susceptibility calculated from first principles in Sec.S9.In addition, because the excitons involved in the tunneling process are supported by the TMD layer, and the 1BZ of graphene is substantially larger than that of the TMD, we treat the conductivity of graphene in the isotropic and translational-invariant approximations [i.e., σ g (k ∥ , ω) is calculated in the RPA [27], as described in Sec.S7], with tensor components given by σ g,GG ′ (k ∥ , ω) = δ GG ′ σ g (|k ∥ + G|, ω) for k ∥ within the 1BZ of the TMD and the G vectors running over the reciprocal lattice of this material.
An example of the obtained response function is shown in Fig. S16 for a heterostructure containing MoSe 2 .We remark that k ∥ refers to the optical wave vector, which is associated with the difference between electron and hole wave vectors of the contributing excitons.The response is dominated by nearly direct excitons at k ∥ ≈ 0 as well as weaker features corresponding to ΓQ TMD excitons away from that region.

S9. First-principles nonlocal optical response of TMD monolayers
The exciton-dominated nonlocal optical conductivities of monolayer TMDs constitute a central ingredient in our theory.We perform first-principles calculations in the framework of many-body perturbation theory to appropriately describe the excited-state properties of the considered materials.The ground-state wave functions and oneelectron eigenenergies are obtained by performing density-functional-theory (DFT) calculations with the Perdew-Burke-Ernzerhof (PBE) [32] exchange-correlation functional as implemented in the Quantum Espresso code [33].In addition, the spin-orbit interaction is included using optimized norm-conserving pseudopotentials [34,35] with a kinetic energy cutoff of 80 Ry.We employ the Coulomb-cutoff technique [36] at the edges of the unit cells in the out-of-plane direction with a vacuum layer of 37 atomic units.
Using the DFT output, we perform Bethe-Salpeter-equation (BSE) calculations within the YAMBO code [37,38].The G 0 W 0 and plasmon-pole approximations [39,40] are employed to correct the Kohn-Sham eigenenergies.In Fig. S17, we present the so-obtained quasiparticle electronic band dispersions for the TMDs under consideration.The direct electronic band gap at the K point is found to be 2.32, 1.95, 2.61, and 2.41 eV for MoSe 2 , WSe 2 , MoS 2 , and WS 2 , respectively.We solve the BSE [41,42] to calculate the eigenenergies and eigenvectors of the excitons supported by these materials, from which we obtain the dielectric function.More precisely, we write the BSE as where a combination of the index j and the in-plane wave vector k ∥ are used to label each exciton state, ℏε j k ∥ is the corresponding exciton energy, ℏε vq and ℏε cq−k ∥ denote valence-and conduction-band quasiparticle energies obtained from G 0 W 0 simulations, |vq; cq − k ∥ ⟩ are pure electron-hole excitonic states, A j vcqk ∥ are the exciton eigenvector coefficients, K eh is the electron-hole interaction kernel, and q are 3D one-electron wave vectors.The generalized oscillator strengths of an excitonic transition can be defined as Ω j k ∥ G = vcq A j vcqk ∥ ⟨vq|e −i(k ∥ +G)•R |cq − k ∥ ⟩, from which we calculate the non-interacting susceptibility as Here, we set f s = 1 (2) for spin-polarized (-unpolarized) calculations, N q is the number of q points in the 1BZ, V denotes the 3D unit-cell volume, and we multiply by the out-of-plane lattice constant of the simulation unit cell t to obtain a 2D surface response function.
For the two-step tunneling process, we describe direct excitonic transitions using the dipole matrix elements ⟨vq|e −ik ∥ •R |cq − k ∥ ⟩, calculated from the YAMBO code for small values of k ∥ , in which the in-plane response is isotropic.To obtain an optical spectrum comparable to experimental results, we use a dense and uniform grid of electron wave vectors q of size 72×72×1 in the G 0 W 0 and BSE calculations.We ignore umklapp processes (i.e., we only use the G = G ′ = 0 element of the χ 0 GG ′ matrix of components labeled by the 2D reciprocal lattice vectors G and G ′ [37]).Transitions between the four highest-valence and the four lowest-conduction bands are incorporated in the BSE calculations.We employ the Scalable Library for Eigenvalue Problem Computations (SLEPC) [43] to obtain detailed information about each BSE eigenmode and validate our results with the full Lanczos-Haydock solution [44].Since the matrix that we need to diagonalize is very large, we only extract the BSE eigenmodes for transitions below 3 eV.The 2D conductivities σ TMD (k ∥ , ω) of the studied TMD monolayers are obtained from the 2D non-interacting susceptibility by using the relation σ TMD (k ∥ , ω) = i ωχ 0 00 (k ∥ , ω)/k 2 ∥ .The calculated conductivities are presented in Fig. S18 as a function of photon energy and parallel wave vector.
In the one-step tunneling process, we consider indirect excitonic transitions by calculating the matrix elements ⟨cq − k ∥ |e i(k ∥ +G)•r |vq⟩ for k ∥ running over the full 1BZ.We also include the first six neighboring Brouillin zones corresponding to the six smallest non-vanishing reciprocal lattice vectors G.To reduce the increased computational burden for such a large set of k ∥ and G vectors, we decrease the q-point grid size to 30×30×1 and do not include spin-orbit coupling in this part of our calculations.In addition, since we focus on momentum-transfer-induced changes in the excitonic spectra, we solve the BSE only for the degenerate highest-valence and lowest-conduction bands.We present examples of the so-obtained conductivity of MoSe 2 in Fig. S19, showing features associated with indirect excitons of high momentum within the 1BZ in optical wave vector space k ∥ .
FIG. S2.Measured I-V characteristics (a) and corresponding dI/dV curve (b) of a second WS2 device at room temperature.

FIG. S5 .
FIG. S5.Measured spectra of the light emitted by the WSe2 device of Fig. S4 at different bias voltages.

2 FIG. S7 .
FIG. S7.Comparison between 50 K and 10 K measurements for the dI/dV curves of (a) MoS2 and (b) WSe2 devices.In (a), some features appear consistently in both measurements close to the exciton resonance, whereas in (b), features at higher energies fluctuate, suggesting that they could originate in measurement instabilities.
FIG. S8.Comparison between dI/dV measurements and a simple tunneling model with an infinitely high tunneling barrier.The model follows the measurement up to 2 V, with significant deviations only appearing above 2.5 V.The device structure under consideration is illustrated in the inset.
FIG. S9.Density of states (DOS) of gold.We compare density-functional-theory (DFT) calculations (including a Gaussian energy smearing of full width at half maximum σ = 0.3 eV) with the free-electron gas (FEG) model.In the latter, the DOS is ρE = (m * Au /π 2 ℏ 3 ) 2m * Au E and the Fermi energy is E Au F = (12π 2 ) 2/3 (ℏ 2 /2m * Au a 2 ) ≈ 5.52 eV, where a = 4.08 Å is the lattice constant and m * Au ≈ me is the effective mass.

FIG
FIG. S10.Electron tunneling pathways.(a) On-scale representation of the in-plane electron momentum distribution in graphene (light blue) and gold (orange), including their respective Fermi wave vectors.Graphene is taken to be negatively doped to a Fermi energy of 0.5 eV relative to the Dirac points.(b) Energy-parallel-momentum diagram corresponding to the configuration in (a) for a graphene-gold bias energy eV b = 3 eV.The downward dashed parabola shows the out-of-plane electron energy in gold at the Fermi level, E ⊥ F = E Au F − ℏ 2 k 2 ∥ /2m * Au , as a function of parallel wave vector k ∥ .(c) Exponential attenuation of gold conduction electrons across the interface with an hBN film (thickness d = 3 nm) as a function of energy and in-plane wave vector, as determined by the spill-out distance LE,k ∥ ≈ ℏ [2m * hBN (E Au F + Eg/2 − E + ℏ 2 k 2 ∥ /2m * Au )] −1/2, estimated from a one-dimensional potential barrier determined by the hBN gap at the Γ point[6] Eg ≈ 8.9 eV and the effective mass along the out-of-plane ΓA direction m * hBN ≈ 0.63 me.In (a) and (b), we consider two possible channels for graphene-to-gold tunneling: (1) one-step direct inelastic transitions to near-threshold gold states around the surface projected Γ point (coinciding with the Γ point of graphene) assisted by the creation of indirect excitons in a neighboring TMD monolayer (dark blue arrows); and (2) two-step quasi-elastic tunneling from graphene to gold states near the surface-projected Γ point assisted by large-momentum phonon excitation (purple arrows), followed by quasi-vertical inelastic decay to states in the vicinity of the Fermi level of gold mediated by direct excitons (green arrows).
FIG. S12.Leading exciton contributions to the tunneling current.We plot Re{σ TMD,00 (k ∥ , ω)} (i.e., the real part of the optical conductivity for G = G ′ = 0, see Secs.S8 and S9) within the 1BZ of MoSe2 at a photon energy ℏω = 2.1 eV, along with the six nearest-neighbor Brillouin zones (hexagons), each of them centered around a G vector.We also show the corresponding circles defined by the conditions |G − k ∥ | = Kg (i.e., the ΓKg distance in graphene).
FIG.S14.Calculated I-V curves for phonon+exciton-assisted two-step tunneling.We plot the current obtained by using Eq.(S27) for MoSe2 (see Fig.4in the main text for other TMDs) with a hBN tunnel barrier distance d = 3 nm.The inset shows the corresponding dI/dV curve.

1 FIG
FIG. S15.Analysis of the dielectric response of the studied multilayer structure.(a) We consider a 2D material described by a nonlocal surface conductivity σ(k ∥ , ω), separated from a gold surface by an hBN tunnel barrier of thickness d.The structure is covered by an upper hBN layer.Gold is described through a surface-corrected nonlocal dielectric function ϵs(k ∥ , ω), whereas hBN is modeled through its anisotropic frequency-dependent permittivity[28] of in-plane and out-of-plane components ϵx(ω) and ϵz(ω), respectively.The gold and upper hBN media are assumed to have an infinite thickness, while the 2D material is treated in the zero-thickness limit.(b,c) The response of the system is described in terms of Fresnel transmission and reflection coefficients for each of the two interfaces.
FIG. S16.Screened interaction in MoSe2.We plot Im{−WGG(k ∥ , d, d, ω)} for a heterostructure like that in Fig. S15 with d = 3 nm and the 2D layer consisting of monolayer MoSe2 and graphene.The upper plot shows the dependence on photon energy and momentum along the indicated excursion within the 1BZ of the TMD for G = 0.The lower plots show the k ∥ dependence within the 1BZ for ℏω = 2.2 eV and two values of G [0 and G0 x, with G0 = 4π/( √ 3a TMD ), where a TMD is the in-plane lattice constant of the TMD].The color scale is shared by all plots and saturated outside the displayed range.

TABLE S1 .
Sizes and TMD coverages of the devices.
S7, but with k ∥ substituted by |k ∥ + G|.We now substitute the reflection and transmission coefficients r 2D and t 2D of the 2D material layer by reflection and transmission tensors of components r 2D,GG ′ and t 2D,GG ′ , labelled by incident (G ′ ) and reflected/transmitted (G) reciprocal lattice vectors (see below).In addition, the reflection coefficient at the BN-Au interface is still described by Eq. (S29a), assuming isotropy and translational invariance (i.e., k ∥ and G are both conserved upon reflection at this interface), but introducing a G dependence by also replacing k ∥ by |k ∥ + G| in such an equation, which we indicate by writing the corresponding reflection coefficient as r BN,Au,G .Using these elements, a straightforward FP analysis allows us to write