Negative Differential Resistance in Boron Nitride Graphene Heterostructures: Physical Mechanisms and Size Scaling Analysis

Hexagonal boron nitride (hBN) is drawing increasing attention as an insulator and substrate material to develop next generation graphene-based electronic devices. In this paper, we investigate the quantum transport in heterostructures consisting of a few atomic layers thick hBN film sandwiched between graphene nanoribbon electrodes. We show a gate-controllable vertical transistor exhibiting strong negative differential resistance (NDR) effect with multiple resonant peaks, which stay pronounced for various device dimensions. We find two distinct mechanisms that are responsible for NDR, depending on the gate and applied biases, in the same device. The origin of first mechanism is a Fabry-Pérot like interference and that of the second mechanism is an in-plane wave vector matching when the Dirac points of the electrodes align. The hBN layers can induce an asymmetry in the current-voltage characteristics which can be further modulated by an applied bias. We find that the electron-phonon scattering suppresses the first mechanism whereas the second mechanism remains relatively unaffected. We also show that the NDR features are tunable by varying device dimensions. The NDR feature with multiple resonant peaks, combined with ultrafast tunneling speed provides prospect for the graphene-hBN-graphene heterostructure in the high-performance electronics.

NDR in double barrier resonant tunneling diodes (DB-RTD), appears when the quasi-bound levels can no longer enhance the tunneling resonantly 17 . Recent theoretical investigations report the appearance of NDR features in pure graphene based devices, involving nanoribbon superlattice 18 , doped junctions [19][20][21][22][23] , tunnel-FET 24,25 , and MOSFET structures 26 . These structures typically employ graphene with fine-tuned bandgap, such that graphene behaves more like a semiconductor. NDR effect is also being reported in a single as well as multilayer heterostructure of graphene-hBN-graphene [27][28][29] . Reference 30 models NDR peak in a near metallic bi-layer graphene device. Apart from this, such devices could also find applications in multi-valued memory 31,32 . The multilayer graphene-hBN-graphene heterostructure based electronic devices particularly attract the attention of engineers due to their relatively simpler fabrication [33][34][35] .
The NDR features in multilayer based devices are to be investigated by exploring current-voltage characteristics as a function of (i) number of hBN layers, (ii) lateral dimensions in determining both the voltage location of NDR peaks and the peak-to-valley ratio, which are essential in the device design, (iii) the role of the asymmetric band offset between hBN and graphene, and (iv) defects and scattering. This, however, necessitates further research to rationalize the underlying physics of the NDR effect and gain insight on how to control its critical properties mentioned above.
In this work, we focus on a prototypical multilayer device shown in Fig. 1(a), which consists of layers of graphene and hBN that are vertically stacked. The graphene layers serve as conducting electrodes with a unique band structure while the hBN layers are tunnel barriers. We model the electron transport in these devices by atomistic non-equilibrium Green's function (NEGF) method. A decoherence mechanism, the electron-phonon scattering is introduced and its impact on NDR effects is presented. Additionally, we demonstrate how the magnitude of current, locations of resonant peaks, and peak-to-valley ratio (PVR) values can be tuned by the device parameters. The modeled devices range from a small system with 6,000 atoms to experimentally feasible sizes up to 70,000 atoms (lateral dimensions 24.6 nm × 27 nm).
Next, section 2 defines our method by discussing the underlying Hamiltonians and the methodology for the computation of quantum transport. Section 3 discusses the results including the two underlying mechanisms with the role of electron-phonon scattering included. Section 4 presents the size scaling analysis followed by section 5 concluding the work.

Method
A prototypical heterostructure consists of two semi-infinitely long monolayer armchair-edged graphene nanoribbon (AGNR) electrodes sandwiching an ultra-thin hBN film, with a vertically applied external gate electric field as shown in Fig 1(a). AGNR is employed because it can be engineered as an intrinsic conductor. This forms a vertical tunneling heterostructure with hBN acting as a potential barrier. The hBN film is sandwiched between a bottom and top AGNR, forming a central overlapping heterostructure/ multilayer region stacked in AB order (Bernal stacking). The lattice constant mismatch between hBN and graphene is negligibly small, 1.8%, therefore, we build the device structure with the uniform lattice constant of graphene (2.46Å) only. The system Hamiltonian is constructed using the nearest neighbor tight binding approximation, with the parameters 11 Only the low energy p z orbitals are considered here; so that the Hamiltonian has the same dimension as the total number of atoms simulated. The effect of number of tunneling hBN layers (N z ), the system width (N x ) and the length of the multilayer stacking region (N y ), where units of N x and N y are number of atoms, on the device performance is investigated. The nanostructure thickness is (N z + 2) in units of atomic layers, which includes the two monolayer graphene sheets at the ends.
The bias across the heterostructure is applied by rigidly shifting the electrostatic potential of the bottom graphene electrode by the amount equal to the applied bias as the metallic graphene layers have much higher conductivity than hBN. The electrostatic potential energy of the bottom layer is U = − eV b , where V b is the applied bias, and the electrostatic potential of the top graphene layer remains zero. The electrostatic potential at each of the sandwiched hBN layers are determined by linearly increasing/decreasing the potential from top to the bottom graphene layer, because the c-axis (out of plane) conductivity of hBN is orders of magnitude smaller than the in-plane conductivity of graphene. The chemical potential of contacts are controlled by bias voltage, namely μ B = − eV b and μ T = 0 for bottom and top graphene leads respectively. The gate voltage is modeled by shifting the electrostatic potential at the bottom electrode by Δ U = − 0.01 eV g , where V g is the gate voltage. We choose N x equal to 3n+ 2, where n is an arbitrary positive integer 37,38 , such that the AGNR have zero bandgap. All calculations were performed at 300 K.
We simulate the transport properties of the device by using the state of the art Green's function (NEGF) formalism 39 . After the generation of real-space Hamiltonian matrix, the retarded Green's function is computed by where E is the electrons energy, I denotes the identity matrix, H is the system Hamiltonian matrix, and L R r Σ / is the self-energy capturing the semi-infinite left (right) graphene sheets shown in Fig. 1. The local density of states at a given atom site i can be calculated subsequently by The transmission coefficient is then calculated using the expression, , and the current as a function of bias is given by, is the Fermi function and μ L/R represent the Fermi level of the left (top) and right (bottom) graphene monolayer electrode. In Fig. 1(b), we plot the DOS for a structure with a width of N x = 200, N y = 32, and with three hBN layers N z = 3, at zero bias. The bandgap of hBN is found to be around 4.72 eV, and the valence band offset between hBN and graphene is around 1.38 eV, which is consistent with the prior work 40 . Traditionally, the recursive Green's function (RGF) approach is applied to solve for the retarded Green's function. Here, an algorithm named HSC-extension is developed, which is more efficient and is based on the nested dissection method, so that the requisite large scale systems can be simulated 41 . Figure 2(a) presents the computed current-voltage characteristics of the heterostructure with a transverse width of N x = 62 (7.6 nm), stacking length N y = 32 (6.8 nm) and three hBN layers (1.4 nm) serving as the tunneling barrier. First we consider the highlighted curve with V g = 0, in which case two NDR peaks emerge at V b = 0.3 V and 0.66 V, respectively. We attribute the formation of these multiple NDR peaks to the Fabry-Pérot like interference in the multilayer region (mechanism 1), as rationalized below.

Mechanisms. Origin of Multiple NDR peaks (Mechanism 1).
To understand the multiple NDR peaks, we calculate the transmission and the average density of states (DOS g and DOS hBN ) at V b = 0.3 V, 0.46 V and 0.66 V, corresponding to the first peak, the first valley and the second peak in the current-voltage characteristics. In the DOS g curves (Fig. 3), the average DOS of bottom and top graphene layers are plotted separately. In particular, the huge peaks in DOS g are marked as peak S, which captures the edge states due to the zigzag-shaped cut-ends of graphene ribbons. The blue DOS hBN curves show the average DOS of the three hBN barrier layers. For the sake of comparison, the transmission coefficient at equilibrium is also plotted with the black dashed curves. The chemical potentials of the bottom and top AGNR are marked as vertical black lines (μ B and μ T ). The transmission and DOS g show a strong Fabry-Pérot like resonant feature in the low energy window. The semi-infinite top and bottom AGNRs couple with hBN at the central heterostructure (multilayer) region. The potential discontinuity caused by the interaction between the hBN cut-ends and the graphene layers    Fig. 2(a). (a)-(c) specifies the bias potential V b = 0.3 V, 0.46 V and 0.66 V respectively. In the transmission plots, black dashed curves are transmission coefficient when V b = 0. In DOS g plots, blue and red curves represent DOS of bottom and top graphene sheets respectively. Vertical dash-dot lines give the chemical potentials at both graphene ends μ B and μ T , which determines the bias window. S mark the DOS peaks resulting from the zigzag shaped edges of graphene cut ends. P mark the transmission peaks that mainly contribute to the current. T represent the tunneling peaks due to the energy alignment of subbands in top and bottom graphene contacts; they do not contribute significantly to current. Units of DOS are number of states per atom per eV. create a resonant cavity in the overlapping region at both the top and bottom graphene layer. When electrons transport across the boundaries between graphene monolayer and hBN multilayer regions, partial reflections occur at the interfaces. As a result, the transmission is oscillatory with peaks and valleys corresponding to constructive and destructive interferences.
The current is determined from the area enveloped by the transmission curve in the energy window bounded by the Fermi levels of two electrodes, μ B and μ T (black dash-dot lines in Fig. 3). The transmission peak (P) at E = − 0.3 eV in Fig. 3(a) mainly contributes to the tunneling current. At low bias regime (V b < 0.18 V), we find that this transmission peak P is enhanced, resulting in the increase of current with applied bias. The resonant tunneling occurs when the constructive quantum interference assists the tunneling of electrons from the top to bottom electrodes at specific energies. When the bias is further increased (until 0.46 V), the transmission peak P is reduced due to destructive interference despite the fact that the energy window for carrying current enlarges. This transmission reduction begins to dominate after V b = 0.3 V, which induces a drop in current. At V b = 0.46 V, there is a large suppression of transmission within the bias window, which creates a large tunneling gap, leading to the current valley as reflected in the highlighted curve of Fig. 2(a). Note that the density of states is large in both graphene electrodes even when the transmission is small as see in Fig. 3(b). Then, the transmission starts to increase again at around V b = 0.66 V due to the constructive interference.
An interesting feature of Fig. 3(a) is that only one transmission peak is observed at around E = − 0.3 eV (μ T ), while a symmetric peak at E = 0 eV (μ B ) is clearly absent. This is due to the fact that the presence of hBN layers break the symmetry. We could understand this from DOS hBN plot at V b = 0.3 V (Fig. 3(a) DOS hBN curve), which shows a peak near E = − 0.3 eV but a valley at E = 0 eV. This means that electrons at E = 0 eV see a stronger barrier when tunneling between AGNR layers, and suggests the breaking of the π− π * symmetry. This argument is tested by considering a symmetric tunnel barrier, where such an asymmetry in transmission does not exist. In Fig. 3(b),(c), sharp transmission peaks (marked as peak T) are observed. This significant tunneling enhancement results from the energy level alignment between the subbands of top and bottom graphene electrodes, as reflected in the corresponding DOS g features.
Gate induced NDR peak (Mechanism 2). We next discuss the second mechanism that leads to single intense NDR peak by investigating the operational behavior of the heterostructure in the presence of an external gate voltage (V g ). Figure 2 shows the current-voltage characteristics for a family of V g ranging from − 45 V to + 45 V for two devices with different sizes. Take the current-voltage curve at V g = − 45V as an example; At V b = 0, the negative gate voltage shifts the energy of Dirac point in the bottom AGNR electrode to U = 0.45 eV at equilibrium, while preserving the chemical potentials from two contacts at μ B = μ T = 0. At V b = 0.45 V (see Fig. 2(a) inset), the Dirac points of bottom and top AGNR electrodes are aligned. As a result, electrons can tunnel from the valence band of the top graphene layer to the conduction band of the bottom graphene layer owing to the in-plane wave vector conservation 35 . This particular mechanism (mechanism 2) induces the resonant transmission and results in the large current peaks marked by solid arrows in Fig. 2. It is noticeable that current peak positions are shifted from the theoretical prediction (V b = − 0.01V g ) based on mechanism 2 only. This occurs when the strength of mechanism 2 is comparable to that of mechanism 1, when the voltage at which the peak current occurs is influenced by the Fabry-Pérot like interference. The superposition of mechanisms 1 and 2 leads to the current peak displacement, which is larger at low gate voltage (for example, V g = − 15 V).
Gate response. Gate voltage has different impacts on NDR induced by two distinct mechanisms discussed above. In Fig. 2(a), the peak current (solid arrows) increases with V g as the peak current is proportional to the number of carriers between μ B and μ T when Dirac points of the top and bottom graphene align (see inset of Fig. 2(a)). In contrast to this, we find that in Fig. 2(a), the peak current of the NDR peaks induced by mechanism 1 (empty arrows) are relatively insensitive to V g . In addition, the PVR values of these NDR peaks are also V g -insensitive because the vertical gate potential tunes the resonant energies for constructive interference, but do not affect the number of tunneling carriers. Consequently, the amplitudes of current peaks for mechanism 2 is weaker than that for mechanism 1 at low V g , but can become significantly stronger at large gate voltage, as shown in Fig. 2(a) when V g = − 45 V. When the device structure is enlarged to N x = 200 and N y = 32, the current-voltage curves for various gate voltages ( Fig. 2(b)) show that the multiple NDR peaks stay clearly defined and their locations are strongly gate-controlled. We point out that the current-voltage curves are asymmetric for positive and negative biases even at V g = 0 as the hBN layers breaks the π− π * symmetry in the multilayer system. We also note that after the NDR peak, our calculations clearly show a trend of increase in current with increase in drain voltage, in a manner qualitatively similar to the experiments 35 .
Effect of scattering. It is to be noted that the mechanism 2 induced single current peak only occurs at specific bias points, determined by V g . Besides, the amplitude of this single peak is V g -sensitive, which helps us to distinguish it from the multiple peaks induced by mechanism 1. However, in the experiment 35 only strong resonant peaks due to mechanism 2 are observed. The disappearance of peaks originated from mechanism 1 merits discussion.
In experiments on large area devices such as Ref. 35, unavoidable decoherence mechanisms such as electron-phonon, electron-electron and electron-environment interactions are non-negligible. To model Scientific RepoRts | 5:10712 | DOi: 10.1038/srep10712 decoherence, here we next introduce electron-phonon scattering in top and bottom graphene electrodes and investigate its influence on NDR peaks. The detailed modeling methodology can be found elsewhere 42 . The typical values of electron-phonon coupling constants and phonon energy are taken from Ref. 43: elastic deformation potential D el = 0.01 eV 2 , inelastic deformation potential D inel = 0.07 eV 2 and phonon energy ħω = 180 meV. The deformation potential values phenomenologically represent the strength of electron-phonon coupling, in terms of the electron mean free path which can be experimentally measured. Physically, larger D el and D inel represent stronger electron-phonon interaction and thus shorter electron mean free path. The simulated electron mean free path corresponding to the above parameters is about 1.48 μm, which represents moderate scattering. The current-voltage curves are plotted in Fig. 4.
Apparently, from Fig. 4, electron-phonon scattering suppresses both NDR mechanisms and therefore, PVR values of these NDR peaks are reduced. However, the suppression of mechanism 2 is not as substantial as that of mechanism 1 as the quantum interference is more vulnerable to decoherence introduced by electron-phonon scattering. This might explain the absence of NDR peaks due to mechanism 1 in the experiments of Ref. 35. However, in an experiment with sufficiently smaller devices, both mechanisms leading to the multiple NDR peaks can occur.
Size scaling analysis. System dimensions are the key ingredients in engineering the device performance. In this particular multilayer heterostructure, for instance, the device width determines the number of subbands in ANGR electrodes and the heterostructure length determines the length of interference region. Based on the two distinct mechanisms responsible for the multiple NDR peaks, it is intuitive that the device dimensions have significant and non-trivial influence on the NDR features rather than simply tuning the current magnitude by following quantum mechanical rules or Ohm's law. In order to comprehend such influences, the scaling analysis of the device dimensions, namely the lateral (x,y) dimensions which defines the overlap area between hBN and graphene layers and the z-direction (number of hBN layers), is performed in this section. However, we do not consider the electron-phonon scattering effects during this analysis as they do not significantly alter the outcomes of the analysis.

Tunneling barrier thickness (N z ).
Representing the thickness of tunneling barrier, N z homogeneously modifies the current magnitude at different applied voltages, whereas has little effects on the peak positions and corresponding PVR. In Fig. 5(a), the hBN thickness (N z ) is varied from 1 layer (0.6 nm) to 5 layers (2 nm), while both N x and N y are fixed. The magnitudes of current are scaled by a multiplicative factor to present results on the same plot for different values of N z . The transmission versus energy (Fig. 5(b)) shows that while the magnitude of transmission depends strongly on N z , the locations of peaks depend weakly on N z . Note that the dependence of current magnitude on N z lose its validity in the case when all the incident modes can tunnel through a thin barrier. This is because a thinner barrier only increases the tunneling probability of electrons without affecting the number of incident modes. Width of AGNR (N x ). For the mechanism 1 (at V g = 0), the density of subbands for the monolayer AGNR electrodes and the heterostructure region depends on the graphene nanoribbon width, i.e. the energy intervals between subbands for the structure with N x = 200 are about three times smaller than that for the structure with N x = 62. Therefore, a larger number of subbands contribute to current under lower biases, resulting in initial increase in current with N x , as seen in Fig. 6(a). When the gate voltage is − 45 V, the NDR peaks induced by the mechanism 2 are observed near V b = 0.45 V for different N x in Fig. 6(b). The heights of these peaks increase with device widths because the number of subbands carrying current between μ B and μ T grows with the width of the AGNR electrode (inset of Fig. 2(a)).  We summarize the peak currents and PVR values for both mechanisms in Table 1. Although the peak current is larger for the wider device, a rapidly decreasing PVR value can be observed. This is because of the stronger band-to-band tunneling between two AGNR contacts with a larger width, arising from the smaller subband spacing. Table 1 exhibits a PVR value up to 13 for N x = 62, which can be further increased to over 60 when N x shrinks to 14, showing the potential for the heterostructure to be utilized in both digital logic and memory. However, in reality to achieve the large PVR values will require a downscaling of N x and minimization of decoherence.
Length of the heterostructure (N y ). The length of the central multilayer region determines the number of incident carriers, and also characterizes the size of the Fabry-Pérot like interference cavity. For mechanism 1, when the heterostructure length N y changes from 16 (3.4 nm) to 64 (13.6 nm), the number of transmission peaks increase as shown in Fig. 7(a). The NDR peaks appear at V b = 0.38 V, 0.8 V for N y = 32, which shifts to lower V b i.e. at 0.2 V, 0.4 V respectively, when N y = 64. This is because the resonant transmission appears at various energies, which vary inversely proportionally with the length of the interference cavity N y . Experiments where the overlap between two graphene nanoribbons are altered should be able to reveal the differences in oscillations of I-V characteristics as a function of N y . We note that experiments with changing overlap dimensions have been performed in carbon nanotubes before 44,45 and future experiments in BN-graphene heterostructures should be useful in studying these features.
With an ultrathin tunneling barrier (N z = 1), electrons have high tunneling probabilities and thus the current is mainly limited by the number of modes incident within the energy window. Graphene layer has low DOS near Dirac point, yielding the saturation of peak current at large. It is also observed that for larger N z (N z ≥ 3, results not shown), the peak current increases with N y rapidly without saturation. This is consistent with our previous discussion since a thicker hBN tunneling barrier greatly suppresses  Fig. 6). Black arrows mark the NDR peaks due to mechanism 2. the overall tunneling transport probability and therefore the peak current magnitude is a strong function of the number of carriers in graphene.

Conclusions
We have systematically investigated the charge transport properties of a three-terminal graphene-hBN-graphene multilayer heterostructure device as a function of device dimensions, so as to further understand the underlying mechanisms for negative differential resistance. The prototypical graphene-hBN-graphene multilayer heterostructure has two distinct mechanisms that can introduce NDR behavior. The first mechanism involves a Fabry-Pérot like resonant feature due to interference in the multilayer heterostructure region, which can produce multiple current peaks. In the presence of an external gate, resonant tunneling can also occur when the electronic spectrum (Dirac points) of the top and bottom graphene electrodes align, which leads to a second mechanism for resonant tunneling. Both mechanisms respond to gate voltage distinctly. Gate voltage only controls the locations of NDR peaks from mechanism 1 while can tune both PVR and locations of NDR peaks due to mechanism 2. In the presence of electron-phonon scattering decoherence, mechanism 1 is more intensively suppressed compared to mechanism 2, which may be relevant to the disappearance of mechanism 1 in the experiments.
Size scaling analysis provides insight into the device physics that determines the number of NDR peaks, the variation of peak current and PVR value with change in device dimensions: (1) The hBN thickness exponentially controls the magnitude of current without significantly affecting the NDR features. (2) For devices with larger widths (N x ), the multiple current peaks preserve but with decreasing PVR values for both mechanisms. (3) For mechanism 1, the bias voltages at which multiple current peaks, occur depend on the length, and the number of peaks increase with length. In contrast to this, the location of the single peak from mechanism 2 is independent of the length. We believe that the negative differential resistance's sensitivity to the system dimensions will provide additional insights for future theoretical and experimental investigations.