Inelastic resonant tunnelling through adjacent localised electronic states in van der Waals heterostructures

Introduction

The advent of 2D crystals and van der Waals heterostructures^1,2 has boosted the development of the field of nanoelectronics offering unprecedented opportunities for creation of functional devices based on novel physical principles or enhanced performance characteristics. Graphene and hBN are the pivotal components required in a variety of device architecture designs^3,4. Vertical tunnelling transistor devices present a special architecture⁵ with an operational principle based on tunnel effect. The resonant mechanisms in vertical tunnelling transistor devices can be broadly categorized into two types – twist-controlled elastic^6,7,8 and a variety of assisted inelastic events. The latter involves an energy loss during tunnelling process due to the interactions with quasiparticles^9,10 or, for instance, with moiré potential¹¹. The palette of assisted inelastic mechanisms can be further expanded to include excitons, magnons, photons and others when replacing the constituting graphene or hBN layers with 2D transition metal dichalcogenides or other representatives of 2D family^12,13,14,15. Notably, replacing the emitter and collector layers with the latter allows observation of twist-controlled and spin-conserved elastic tunnelling events¹⁶. Resonant tunnelling can also manifest without third quasiparticle’s involvement. Those resonances, however, are of an elastic nature and arise due to the sequential tunnelling through atomic defect based individual localised electronic states within the tunnel barrier^17,18,19.

The elastic sequential resonant tunnelling through a singular localised electronic state in hBN emerges when the chemical potential of collector or emitter layer aligns with it energetically. This results in the opening of a conductive channel with step (peak) resonant feature in current-voltage (conductance-voltage) characteristics of the device^17,18 and can be of great use for infrared photodetection applications²⁰. Notably, introducing several localised electronic states usually results in the emergence of parallel conductive channels¹⁹ exhibiting similar features. Nevertheless, late 20^th century studies of tunnelling through two localised electronic states highlighted the role of more complex dynamics at play showing that those may interact within the tunnel barriers²¹ and lead to the rise of correlated current behaviours²². Furthermore, more recent studies demonstrate that the sequential resonant tunnelling through localised electronic state in hBN may also manifest as with an assisted mechanism as an inelastic event when the corresponding atomic defects within the barrier are coupled to the lattice phonons²³.

It is worth noting, that atomic defects in hBN can serve as effective visible light single-photon emitters^24,25 offering numerous advantages, such as room-temperature operation and potential compatibility with photonic devices. Recent studies demonstrate that hBN can host a variety^26,27 of optically active atomic defects including carbon or oxygen impurities or boron and nitrogen vacancies. The former can be controllably introduced into the hBN crystal. Those defects exhibit strong photo- and electro-luminescent narrow emission lines^28,29, which are essential characteristics for designing next generation quantum light sources^30,31 based on single-photon emission mechanism. Despite this promising development, challenges remain in understanding the underlying mechanisms governing the emission properties of a variety of atomic defects in hBN.

In this work, we use dual gated vertical tunnelling transistors to study resonant electron tunnelling through two adjacent localised electronic states in few atomic layer hBN barriers sandwiched in-between moiré monolayer and AB Bernal bilayer graphene collector and emitter layers. We demonstrate that the resonant tunnelling through such two adjacent localised electronic states can be of an inelastic nature as a generic phenomenon in the case of wide tunnel barriers and short localization lengths. We also display that those inelastic tunnelling events dominate all over the other channels of elastic events through individual localised electronic states giving rise to explicitly strong resonant features. It allows their utilization in accurate resonant tunnelling spectroscopy of delicate features in electronic density of states of emitter and collector layers, such as second neutrality point valence-band bandgap of moiré monolayer and electric field induced bandgap of Bernal bilayer graphene. Our work enriches the understanding of interaction mechanisms among localised electronic states in hBN facilitating future investigations into their potential applications.

Inelastic resonant tunnelling through adjacent localised electronic states in hBN

The tunnelling current resulting from inelastic resonant events between two adjacent localised electronic states can, under certain conditions, surpass the elastic resonant current associated with each state. As illustrated in Fig. 1a, a schematic band diagram depicts the distribution of localised electronic states (A) and (B) within the hBN tunnel barrier. Notably, the resonant feature in the tunnelling conductance as a function of bias voltage—corresponding to inelastic transitions between states (A) and (B) in series (denoted as (A) (to) (B) in Fig. 1b)—is approximately an order of magnitude greater than elastic tunnelling features observed through the individual electronic states ({A}) or ({B}). In this context, an electron in state (A) has the option to either directly traverse to the drain or take a shorter route to state (B). The latter pathway may incur the additional cost of emitting a phonon⁹ or a plasmon¹⁰ or emerge due to Auger-like process.

**Fig. 1: Inelastic resonant tunnelling through two adjacent localised states in hBN barrier, conceptualization.**

The inelastic tunnelling event’s rate contains extra small coupling constant (g) responsible for electron-phonon and electron-electron interactions. It exhibits a weak (non-exponential) dependence on the distance between the localised electronic states. When separating the elastic and inelastic components of the total (from source to drain) resonant tunnelling current in the balance equation (see Supplementary Note 1 for further details), it becomes evident that the inelastic component predominates if ({{|g|}}^{2} > {e}^{-2{{kappa }L}_{{SA}}}+{e}^{-2{{kappa }L}_{{BD}}}), where ({kappa }) represents the localization length of the electronic state, ({L}_{{SA}}) (({L}_{{BD}})) denotes the distance from the source ((B)) to localised state (A) (drain). Notably, inelastic resonant events consistently dominate when the tunnelling probability is low, such as in scenarios involving wide barriers or short localization lengths. The evaluated plots presenting the inelastic and elastic components of the tunnelling current as a function of the total attenuation between the source and drain layers, ({e}^{-2{kappa }d}), are illustrated in Fig. 1c. Our calculations incorporate the precise positions of localised electronic states within the hBN barrier (see Methods for modeling electrostatic parameters), which were determined through fitting resonant tunnelling features and will be elaborated upon in next section. Here, the inelastic current significantly surpasses the elastic current as the tunnelling weakens (({ sim 10}^{-7}) for ({{|g|}}^{2}=0.1)). This intriguing behaviour, where higher-order inelastic tunnelling events prevail over first-order elastic events, can be attributed to the selection of the easiest tunnelling path for the electron. If the coupling is weak, i.e. the barrier is wide, and is also high, the electron would always “choose” an inelastic path.

Electric field fine-tuning of resonant features arisen from inelastic tunnelling events through adjacent states

To investigate the resonant features arising from our inelastic tunnelling events, we conducted low-temperature differential tunnelling conductance measurements for our dual gated Device 1 with moiré monolayer and Bernal bilayer graphene layers used as bottom and top layers correspondingly (refer to Methods for further details and Supplementary Note 2 for the device schematics). The results are presented as a function of backgate and bias voltages for several fixed topgate voltages of ({V}_{{rm{tg}}}=2,3,4,{rm{V}}) in Figs. 2a, d, and g, respectively. The maps reveal that certain resonant features, manifested as peaks in the conductance, stand out prominently against a backdrop of various weaker features. This prominence can be attributed to the predominance of inelastic tunnelling events occurring through adjacent localised electronic states (A) and (B) within hBN barrier. Additionally, a range of other resonant features emerge from the tunnelling through individual localised electronic states, which are of a different origin and give rise to elastic resonant tunnelling. These include tunnelling features associated with localised electronic states (A) or (B), as well as contributions from numerous other states located in different layers of hBN barrier (see Supplementary Note 3 for more details). Notably, recent studies²³ display that some of weaker features observed can be also attributed to the defect-phonon coupling assisted tunnelling through hBN barriers. Here, we limit our discussion to the generic phenomenon describing conceptually the inelastic resonant tunnelling through two adjacent localised electronic states.

**Fig. 2: The dominance of inelastic resonant tunnelling events and their electric field tuneability.**

To leverage our findings, we model the electrostatic parameters of our devices (see Methods for details) to elucidate the characteristics associated with inelastic tunnelling events. In line with³¹, our approach employs gated configuration featuring graphene electrodes being gated by metallic gates with an additional Au topgate (see Supplementary Note 2). To characterize the energetic spectrum of Bernal bilayer graphene, we adopt a four-band model with a self-consistent bandgap, as introduced in ref.³², alongside an effective screening approach for evaluating the electric field between the graphene layers. Next, we apply modified Landauer-Buttiker model to derive an expression for inelastic resonant tunnelling conductance (current) through our localised electronic states (A) and (B) (see Methods for further details). Our best fit to the experimental data for Device 1 indicates that the hBN tunnel barrier consists of six monolayers ((d)) with localised states of (A) and (B) positioned at ({z}_{{rm{A}}}=d/)6, and ({z}_{{rm{B}}}=) 2 (d/)3. Their zero electric field energies correspond to ({E}_{0{rm{A}}}=) 83 (,{rm{meV}}),({E}_{0{rm{B}}}=) 73 (,{rm{meV}}). As a result, when introducing the resonant conditions of ({mu }_{{rm{Gr}}} > {E}_{{rm{A}}} > {E}_{{rm{B}}} > {mu }_{{rm{BGr}}}) and ({mu }_{{rm{Gr}}} < {E}_{{rm{A}}} < {E}_{{rm{B}}} < {mu }_{{rm{BGr}}}), where ({mu }_{{rm{Gr}}}) (({mu }_{{rm{BGr}}})) denote the chemical potentials of moiré monolayer (Bernal bilayer) graphene, and ({E}_{{rm{A}}}) (({E}_{{rm{B}}})) represent the in-field energies of localised electronic state (A) ((B)), we generate differential tunnelling conductance maps. These maps, plotted as a function of backgate and bias voltages, are presented in Figs. 2b, e, h for fixed topgate voltages of ({V}_{{rm{tg}}}=) 2, 3, 4 (,{rm{V}}), respectively. Here, dashed white lines indicate deeps in the differential tunnelling conductance illustrating the alignment of chemical potential of moiré monolayer graphene with its neutrality point.

When we take a close look at the differential tunnelling conductance maps, we can see that inelastic resonant tunnelling events between two localised electronic states (A) and (B) are influenced significantly by the electric field, specifically, the topgate voltage applied to the device. Figure 2c presents the band diagram of an arbitrary point of the inelastic resonant tunnelling conductive channel indicated by white dot in Fig. 2b. Here, the changes in electrostatic field (due to the application of bias and backgate voltage) along the conductive channel shift the energy levels of localised electronic states ({A}) and (B) fulfilling the resonant conditions all over the channel. Figure 2 f, i present the band diagrams at more specific points also indicated by white dots in Fig. 2 e, h. At these points, the energy levels of states (A) and (B) align causing the tunnelling process to shift from inelastic to elastic. Interestingly, when topgate voltage of 4 V is applied, the four conductive channels converge at a single point. Notably, topgate voltage can be finely adjusted to tune (A) and (B) states in such a way to fail the resonant conditions along the inelastic sequential tunnelling conductive channel (see Supplementary Note 4 for additional details). Those regions are highlighted as a transparent shaded area in Fig. 2b, e. The fact that resonant inelastic tunnelling can take place already at very small bias voltages (see Fig. 2g, h) implies the absence of energy threshold for this process. Obviously, it does not necessarily exclude the resonant inelastic tunnelling associated with phonon emission from physical candidates as it also may take place for higher bias voltage cases, but predominantly displays that the origin of these inelastic transitions can be associated with Auger-like processes or, for instance, an emission of 2D plasmons possessing gapless energy spectrum.

In contrast, within the weak coupling regime and at high barrier conditions, the tunnelling current resulting from inelastic resonant events between adjacent localised electronic states does not consistently prevail over elastic events. For Device 2, which comprises nine monolayers of hBN that form the tunnel barrier, that features several localised states distributed throughout the barrier, we observe only elastic resonant events under sequential tunnelling mechanisms through two localised states (refer to Device 2 data in Supplementary Note 5). This event manifests as a pronounced negative differential conductance as a function of bias voltage and is comparable in magnitude to the resonant events associated with individual states. Moreover, the tuning of the electric field can indeed shift the relative positions of these elastic resonant features in relation to the backgate and topgate, however, this adjustment does not result in the emergence of a forbidden tunnelling regions observed in Device 1.

Inelastic resonant tunnelling for the spectroscopy of delicate electronic density of state features

Studies of low-temperature differential tunnelling conductance at elevated backgate voltages reveal that our inelastic tunnelling events can serve as a powerful tool for probing subtle features in the electronic density of states of emitter or collector layers (see Fig. 3). This technique further boosts the defect-assisted tunnelling spectroscopy³³ enabling the visualization³⁴ of electronic bandgap at the second neutrality point of the bottom moiré monolayer graphene, which is aligned with underlying encapsulating hBN layer. Additionally, it facilitates the measurement of the electric field-tuneable bandgap of top Bernal bilayer graphene with an exceptional precision. The experimental differential tunnelling conductance contour plot, as a function of bias and backgate voltages at fixed topgate voltage of ({V}_{{rm{tg}}}=4,{rm{V}}), is presented in Fig. 3a. The evaluated differential tunnelling conductance map for our inelastic events is shown in Fig. 3b. In both cases, resonant features associated with inelastic events through two adjacent localised electronic states ({A}) and (B) manifest as a series of kinks and forbidden regions within a limited windows of bias and backgate voltages applied to the device. These regions are indicated by transparent squares in the experimental (see Fig. 3a) and theoretical (see Fig. 3b) maps and are further detailed in Figs. 3c, e (experimental) and Figs. 3d, f (theoretical) as zoomed-in maps. Figures 3c, d present the case of a reconstruction of energetic band structure at the valence band second neutrality point of moiré monolayer graphene. Figures 3e, f, on the other hand, illustrate the case of electric field induced bandgap opening of Bernal bilayer graphene. Notably, our experimental results are in an excellent agreement with theoretical evaluations.

**Fig. 3: Resonant inelastic tunnelling spectroscopy of delicate features in electronic density of states of emitter and collector layers.**

The bias and backgate voltage positions of these resonant features enable the determination of the twist angle in moiré monolayer graphene³⁵ measured as the distance between the two parallel dashed lines in Fig. 3b. Here, similar to the approach used for the extraction of second neutrality point energetic values from the gate dependence of field-effect moiré monolayer transport devices³⁶, we used Eq. (1) system of electrostatic equations to obtain the corresponding lines when the chemical potential of moiré monolayer nests at its first neutrality point and valence band second neutrality point (see white dashed lines in Fig. 2). Our evaluations showed that the latter corresponds to the energies of 220 meV. Next, it was directly converted to the twist angle of approximately 1° using the approach presented in ref.^37,38. This twist angle induces a bandgap at the valence band of the energetic spectrum^37,38,39, which starts to reconstruct at 220 meV. Here, the resonant tunnelling from electronic state (A) into the collector moiré graphene (see the inset in Fig. 3a) is effectively prohibited within a narrow range of bias and backgate voltages, particularly when its chemical potential nests at second neutrality point bandgap (see Fig. 3c), which we observe to be of ({varDelta }_{2{rm{NPmbox{-}Gr}}}=) 7 meV. The same twist angle of moiré graphene layer also influences the conduction band. However, in this case, it instead manifests as a subtle, smeared feature at the background of inelastic and elastic resonant features (see the loci distortions indicated by black arrows in Supplementary Fig. 4), which we attribute to the absence of conduction band second neutrality point bandgap. Instead, it represents a deep in the dependence of differential tunnelling conductance from bias or backgate voltages.

Utilizing a similar approach, one may enable the determination of the electric field-tuneable electronic bandgap of Bernal bilayer graphene. Here, the emergence of the bandgap results in a distinct window of bias and gate voltages (at a fixed topgate voltage) characterized by forbidden states for tunnelling events from state (B) to the collector layer. This occurs when the chemical potential of bilayer graphene resides within its bandgap (see the inset in Fig. 3b). For topgate voltage of ({V}_{{rm{tg}}}=4,{rm{V}}), as illustrated in Fig. 3, we observe strong correlation between experimental and theoretical results, yielding to the bandgap of ({Delta }_{{rm{BGr}}}=42,{rm{meV}}) (see Figs. 3e, f). Our technique facilitates high precision extraction of the bandgap of Bernal bilayer graphene across various topgate voltages (see Supplementary Fig. 7). Furthermore, under perpendicular magnetic fields, the studies of differential tunnelling conductance displays that our resonant inelastic tunnelling events can also serve as a valuable tool for magneto-spectroscopy, like quantum-dot assisted spectroscopy⁴⁰, revealing subtle features, such as spin and valley associated gaps in quantized spectra of emitter/collector layers at moderate magnetic fields (see Supplementary Fig. 8).

In conclusion, our work highlights the dominance of inelastic resonant tunnelling events over elastic tunnelling events for two localised electronic states closely positioned in wide hBN barriers and exhibiting short localization lengths. Our findings reveal that inelastic tunnelling current resonances, which may be driven by interactions like electron-phonon, can exceed the elastic current resonances by an order of magnitude under the conditions of small tunnelling probabilities. We also demonstrate that the observed strong resonances associated with inelastic tunnelling through two adjacent localised electronic states can be of great use for the spectroscopy of delicate electronic density of state features in collector and emitter layers. This includes accurate spectroscopy of second neutrality point valence-band bandgap of moiré monolayer and electric field induced bandgap of Bernal bilayer graphene. Our observations not only deepen the understanding of underlying mechanisms governing the electron tunnelling in quantum systems, but also suggest potential pathways for designing unique device architectures for future electronic applications.

Methods

Device preparation

Tunnelling devices were fabricated using the dry transfer technique⁴¹ applied to the micro-mechanically cleaved layers of monolayer graphene, Bernal bilayer graphene, and hBN of various thicknesses (tunnelling layer, two encapsulating top and back gate electrode layers) from bulk graphite and hBN crystals. Those were assembled on top of each other into van der Waals heterostructures using PC (PolyBisphenol carbonate) stamps prepared on commercial PDMS (poly-dimethylsiloxane) films, and eventually, were deposited on top of 285 nm thick silicon dioxide/strongly (p)–doped silicon wafers. Twist angle of about 1° was introduced in-between the bottom encapsulating layer of hBN and bottom monolayer graphene during the fabrication procedure for Device 1. Next, Cr/Au edge contacts⁴² were made on the bottom monolayer and top Bernal bilayer graphene using electron-beam lithography followed by hBN etching, metal deposition and a lift-off process. hBN was etched in a reactive ion etching system using CHF₃ chemistry. Contacts were made in such a way to have a four-probe measurement geometry. The top encapsulating hBN layer was additionally covered by a Cr/Au pad at the cross-sectional area of monolayer and Bernal bilayer graphene of (Sapprox) 20 ({rm{mu }}{{rm{m}}}^{2}) ((Sapprox) 25 ({rm{mu }}{{rm{m}}}^{2})) for Device 1 (Device 2), which served as topgate electrode.

Measurement technique

Tunnelling current-voltage curves were measured using a 2636B source meter from Keithley Instruments. The differential tunnelling conductance measurements were performed using an ({rm{AC}})–({rm{DC}}) mixing technique with the help of symmetric voltage divider, Keithley Instruments 2636B source meter (for ({rm{DC}}) voltage), and Stanford Research SR860 lock-in amplifier (for output ({rm{AC}}) sine voltage at fixed frequency of (f=) 7 (,{rm{Hz}})). The measurements were performed in a two–probe configuration considering that the contact and lateral resistances of electrode layers were several orders of magnitude smaller than the tunnelling resistance at high biases, backgated and topgate voltages. To obtain high-resolution data no additional device was introduced to the measurement circuits. The samples were held in a closed-loop cryostat at fixed temperature of (T=) 2 (,{rm{K}}).

Theoretical simulations

An electrostatic model of the device under the study with an hBN tunnel barrier of width of (d), the gap between the two graphene layers of ({d}^{{prime} }) in Bernal bilayer graphene and with distances to the top ({d}_{{rm{tg}}}) and back gates ({d}_{{rm{bg}}}) can be given by the following system of equations:

$$left{begin{array}{c}{{eV}}_{{rm{b}}}={mu }_{{rm{Gr}}}-{mu }_{{rm{BGr}}}-{{rm{e}}{dF}}_{{rm{b}}}\ {{eV}}_{{rm{bg}}}={mu }_{{rm{Gr}}}+{{ed}}_{{rm{bg}}}{F}_{{rm{g}}}\ {{eV}}_{{rm{tg}}}={mu }_{{rm{Gr}}}-{{ed}}_{{rm{tg}}}{F}_{{rm{t}}}-{{edF}}_{{rm{b}}}\ {varDelta }_{{rm{BGr}}}=4{rm{pi }}{e}^{2}{d}^{{prime} }({n}_{{rm{Gr}}}left({mu }_{{rm{Gr}}}right)-{n}_{{rm{BGr}}}^{{prime} }({mu }_{{rm{BGr}}},{varDelta }_{{rm{BGr}}}))end{array}right.,$$

(1)

where ({V}_{{rm{tg}}}), ({V}_{{rm{bg}}}) and ({V}_{{rm{b}}}) are the backgate, topgate and bias voltages, ({F}_{{rm{g}}}) is the electric field between the monolayer graphene and Si backgate electrode, ({F}_{{rm{t}}}) is the electric field between the Bernal bilayer graphene and Au topgate electrode, ({F}_{{rm{b}}}) is the electric field within the tunnel barrier region, ({mu }_{{rm{Gr}}}({mu }_{{rm{BGr}}})) is the chemical potential of monolayer (Bernal bilayer) graphene with respect to its neutrality point (or the middle of the bilayer’s bandgap, in case it is open), ({n}_{{rm{Gr}}}left({mu }_{{rm{Gr}}}right)={mu }_{{rm{Gr}}}^{2}/pi {hslash }^{2}{v}_{0}^{2}) is the carrier density in monolayer graphene, ({n}_{{rm{BGr}}}({mu }_{{rm{BGr}}},,{varDelta }_{{rm{BGr}}})=m/pi {Re}[{Sqrt}(2{mu }_{{rm{BGr}}}^{2}-{varDelta }_{{rm{BGr}}}^{2})]) is the carrier density in Bernal bilayer graphene, and ({n}_{{rm{BGr}}}^{{prime} }({mu }_{{rm{BGr}}},{varDelta }_{{rm{BGr}}})) (=) ({n}_{{rm{BGr}}}({mu }_{{rm{BGr}}},{varDelta }_{{rm{BGr}}})/2) (+) (varLambda {n}_{{rm{S}}}({mu }_{{rm{BGr}}},{varDelta }_{{rm{BGr}}})) is the carrier density in one particular layer of Bernal bilayer, where ({n}_{{rm{S}}}({mu }_{{rm{BGr}}},{varDelta }_{{rm{BGr}}})) is the screening term, and (varLambda approx 1) is the dimensionless parameter describing the effectiveness of the screening in the bilayer graphene. Note, that (varLambda =0) describes poor screening when the density on each layer is equal to ({n}_{{rm{BGr}}}({mu }_{{rm{BGr}}},{varDelta }_{{rm{BGr}}})/2).

Similar to¹⁷, an inelastic sequential tunnelling current ((J)) can be presented by phenomenological model based on a modification of the Landauer-Buttiker’s model using the following equations:

$$begin{array}{c}{J}_{{rm{A}}- > {rm{B}}}=-xi int {dE}{L}_{{rm{A}}}left(E-{E}_{{rm{A}}}right){varGamma }_{{rm{B}}}left(E-{E}_{{rm{B}}}right)times frac{{gamma }_{{rm{b}}}{gamma }_{{rm{AB}}}}{{gamma }_{{rm{b}}}+{gamma }_{{rm{AB}}}}{f}_{{rm{BGr}}}left(Eright)left(1-{f}_{{rm{Gr}}}left({E}_{{rm{A}}}right)right),\ {J}_{{rm{B}}- > {rm{A}}}=-xi int {dE}{L}_{{rm{B}}}left(E-{E}_{{rm{B}}}right){varGamma }_{{rm{A}}}left(E-{E}_{{rm{A}}}right)times frac{{gamma }_{{rm{t}}}{gamma }_{{rm{AB}}}}{{gamma }_{{rm{t}}}+{gamma }_{{rm{AB}}}}{f}_{{rm{Gr}}}left(Eright)left(1-{f}_{{rm{BGr}}}left({E}_{{rm{B}}}right)right),end{array}$$

(2)

where ({gamma }_{{rm{t}}}({gamma }_{{rm{b}}})) and ({gamma }_{{rm{AB}}}) are the electronic tunnelling rates from a localised state (A)((B)) into the top (bottom) electrode, and in-between (A) and (B) localised states, respectively, ({f}_{{rm{Gr}},{rm{BGr}}}=)1(/()1(+exp ((E-{mu }_{{rm{Gr}},{rm{BGr}}})/{k}_{{rm{B}}}T)) are the corresponding Fermi functions with a temperature (T) and Boltzmann constant ({k}_{{rm{B}}}), (varGamma) is a singly peaked function with full width half maximum (FWHM) of (gamma), and ({E}_{0}={E}_{0{rm{i}}}+{{rm{e}}F}_{{rm{b}}}{z}_{{rm{i}}}) (({E}_{0{rm{i}}}) corresponds to the unperturbed energies of (A) and ({B}) localised states, ({z}_{{rm{i}}}), corresponds to their positions in the tunnelling barrier). Here, the tunnelling is possible only in cases when ({mu }_{{rm{Gr}}} > {E}_{{rm{A}}} > {E}_{{rm{B}}} > {mu }_{{rm{BGr}}}) and ({mu }_{{rm{Gr}}} < {E}_{{rm{A}}} < {E}_{{rm{B}}} < {mu }_{{rm{BGr}}}). We used the experimental resonant features presented in Fig. 2 of the main text to get the theoretical curves that are obtained using the electrostatic model closest to them by varying the model parameters, i. e., the energies of localised states at zero electric field and their locations within the hBN barrier). For each localised state, the zero-field energy, and the location within the hBN barrier at which the best match with the experimentally observed maps is achieved were successively chosen (({E}_{0{rm{A}}}=) 83 (,{rm{meV}}), ({E}_{0{rm{B}}}=) 73 (,{rm{meV}}), ({z}_{{rm{A}}}=d/) 6, and(,{z}_{{rm{B}}}=) 2(d/)3()).