Introduction

Graphene, a two-dimensional material with unique mechanical, thermal and electronic transport properties1 is a promising candidate for nanodevices as it is deeply scaled in one dimension and the lithography offers scaling in the other two dimensions. Building devices based on graphene is, however, partially impeded by the lack of compatible insulating substrate. Hexagonal boron nitride (hBN) has an atomically smooth two dimensional (2D) layered structure with a lattice constant very similar to that of graphene (1.8% mismatch), sufficiently large electrical band gap (~4.7 eV) and excellent thermal and chemical stability2, allowing it to be stacked with graphene to build device structures with desired functionalities. Also, hBN reduces the surface roughness of graphene without degrading its giant mobility3,4. The nano-scale devices based on graphene employing atomically thin hBN with novel electrical and optical properties have recently been reported5,6,7,8,9,10,11,12,13,14,15. An appearance of negative differential resistance (NDR) in such devices further interests the researchers as it could potentially impact the number of applications such as high-speed IC circuits, signal generators, data storage and so on16.

NDR in double barrier resonant tunneling diodes (DB-RTD), appears when the quasi-bound levels can no longer enhance the tunneling resonantly17. Recent theoretical investigations report the appearance of NDR features in pure graphene based devices, involving nanoribbon superlattice18, doped junctions19,20,21,22,23, tunnel-FET24,25 and MOSFET structures26. 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-graphene27,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 memory31,32. The multilayer graphene-hBN-graphene heterostructure based electronic devices particularly attract the attention of engineers due to their relatively simpler fabrication33,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).

Figure 1
figure 1

(a) A schematic view of the heterostructure device. Nx and Ny represent the width and stacking length of the device respectively. Nz is the number of hBN layers sandwiched between the two AGNR ribbons. All the dimensions are in unit of atomic layers. (b) The average DOS versus Energy of hBN for device with Nx = 200, Ny = 32 and Nz = 3. This shows a 4.72 eV bandgap of atomically thin hBN material and a 1.38 eV valence band-offset between graphene and hBN stacking structure.

Full size image

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 parameters11,36: , , and , , , . Only the low energy pz 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 (Nz), the system width (Nx) and the length of the multilayer stacking region (Ny), where units of Nx and Ny are number of atoms, on the device performance is investigated. The nanostructure thickness is (Nz+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 = −eVb, where Vb 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 = −eVb 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 eVg, where Vg is the gate voltage. We choose Nx equal to 3n+2, where n is an arbitrary positive integer37,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) formalism39. 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 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,

where and and the current as a function of bias is given by,

where f(E) 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 Nx = 200, Ny = 32 and with three hBN layers Nz = 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 work40.

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 simulated41.

Mechanisms

Origin of Multiple NDR peaks (Mechanism 1)

Figure 2(a) presents the computed current-voltage characteristics of the heterostructure with a transverse width of Nx = 62 (7.6 nm), stacking length Ny = 32 (6.8 nm) and three hBN layers (1.4 nm) serving as the tunneling barrier. First we consider the highlighted curve with Vg = 0, in which case two NDR peaks emerge at Vb = 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.

Figure 2
figure 2

(a) Current versus drain voltage for a device with Nx = 62, Ny = 32 and Nz = 3. Vg varies from −45 V to 0. The black solid arrows in the four plots mark the current resonant peaks due to mechanism 2 and the empty arrows marks the NDR peak due to mechanism 1. The inset explains the resonant tunneling induced from mechanism 2. The difference between Fermi energy and Dirac point in bottom graphene is induced by gate potential. When Vb = −0.01 Vg, the electronic spectra of top and bottom electrodes are tuned into alignment, allowing the resonant tunneling. (b) Current versus drain voltage for large device with Nx = 200, Ny = 32 and Nz = 1. Here Vg varies from −45 V to +45 V. Inset shows an asymmetric PVR relationship with the applied vertical gate potential.

Full size image

To understand the multiple NDR peaks, we calculate the transmission and the average density of states (DOSg and DOShBN) at Vb = 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 DOSg curves (Fig. 3), the average DOS of bottom and top graphene layers are plotted separately. In particular, the huge peaks in DOSg are marked as peak S, which captures the edge states due to the zigzag-shaped cut-ends of graphene ribbons. The blue DOShBN 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 DOSg 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 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.

Figure 3
figure 3

Transmission and DOS plot at various biases for I-V curve (Vg = 0) in Fig. 2(a). (a)–(c) specifies the bias potential Vb = 0.3 V, 0.46 V and 0.66 V respectively. In the transmission plots, black dashed curves are transmission coefficient when Vb = 0. In DOSg 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.

Full size image

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 (Vb < 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 Vb = 0.3 V, which induces a drop in current. At Vb = 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 Vb = 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 DOShBN plot at Vb = 0.3 V (Fig. 3(a) DOShBN 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), 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 DOSg 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 (Vg). Figure 2 shows the current-voltage characteristics for a family of Vg ranging from −45 V to + 45 V for two devices with different sizes. Take the current-voltage curve at Vg = −45V as an example; At Vb = 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 Vb = 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 conservation35. 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 (Vb = −0.01Vg) 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, Vg = −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 Vg 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 Vg. In addition, the PVR values of these NDR peaks are also Vg -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 Vg, but can become significantly stronger at large gate voltage, as shown in Fig. 2(a) when Vg = −45 V. When the device structure is enlarged to Nx = 200 and Ny = 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 Vg = 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 experiments35.

Effect of scattering

It is to be noted that the mechanism 2 induced single current peak only occurs at specific bias points, determined by Vg. Besides, the amplitude of this single peak is Vg -sensitive, which helps us to distinguish it from the multiple peaks induced by mechanism 1. However, in the experiment35 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 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 elsewhere42. The typical values of electron-phonon coupling constants and phonon energy are taken from Ref. 43: elastic deformation potential Del = 0.01 eV2, inelastic deformation potential Dinel = 0.07 eV2 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 Del and Dinel 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.

Figure 4
figure 4

Current-voltage curves at Vg = −45 V for devices with Nz = 1, Nx = 62 and Ny = 64 with consideration of electron-phonon scattering. The deformational potentials used in the calculation are Del = 0.01 eV2, Dinel = 0.07 eV2. The simulated electron mean free path in the multilayer system is about 1.48 μm.

Full size image

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 (Nz)

Representing the thickness of tunneling barrier, Nz 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 (Nz) is varied from 1 layer (0.6 nm) to 5 layers (2 nm), while both Nx and Ny are fixed. The magnitudes of current are scaled by a multiplicative factor to present results on the same plot for different values of Nz. The transmission versus energy (Fig. 5(b)) shows that while the magnitude of transmission depends strongly on Nz, the locations of peaks depend weakly on Nz. Note that the dependence of current magnitude on Nz 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.

Figure 5
figure 5

(a) Current-voltage curves for devices with different Nz, with fixed Nx = 62 and Ny = 32. Here the current value for cases when Nz = 3 and Nz = 5 are scaled by 1E2 and 1E4 respectively. The inset plots the low bias conductance of the three current-voltage curves. (b) Transmission relationship for devices with different Nz and fixed Nx = 62 and Ny = 32 at Vb = 0.3 V, corresponding to the first current peaks shown in (a). Again, the transmission coefficient value for cases when Nz = 3 and Nz = 5 are scaled by 1E2 and 1E4 respectively.

Full size image

Width of AGNR (Nx)

For the mechanism 1 (at Vg = 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 Nx = 200 are about three times smaller than that for the structure with Nx = 62. Therefore, a larger number of subbands contribute to current under lower biases, resulting in initial increase in current with Nx, as seen in Fig. 6(a). When the gate voltage is −45 V, the NDR peaks induced by the mechanism 2 are observed near Vb = 0.45 V for different Nx 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 Peak current and PVR values as a function of Nx for both mechanisms. (I-V curves from Fig. 6).
Full size table
Figure 6
figure 6

(a) Current versus drain voltage at Vg = 0 for devices with various Nx, with fixed Ny = 32 and Nz = 1. (b) Current versus drain voltage at Vg = − 45 V for devices with various Nx, with Ny = 32 and Nz = 1. Black arrows mark the NDR peaks due to mechanism 2.

Full size image

Table 1 exhibits a PVR value up to 13 for Nx = 62, which can be further increased to over 60 when Nx 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 Nx and minimization of decoherence.

Length of the heterostructure (Ny)

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 Ny 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 Vb =  0.38 V, 0.8 V for Ny = 32, which shifts to lower Vb i.e. at 0.2 V, 0.4 V respectively, when Ny = 64. This is because the resonant transmission appears at various energies, which vary inversely proportionally with the length of the interference cavity Ny. 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 Ny. We note that experiments with changing overlap dimensions have been performed in carbon nanotubes before44,45 and future experiments in BN-graphene heterostructures should be useful in studying these features.

Figure 7
figure 7

(a) Current versus drain voltage at Vg = 0 for devices with various Ny , with fixed Nx = 200 and Nz = 1. (b) Current versus drain voltage at Vg = −45 V for devices with various Ny, with Nx = 200 and Nz = 1. Black arrows mark the NDR peaks due to mechanism 2.

Full size image

With an ultrathin tunneling barrier (Nz = 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 Nz (Nz  ≥ 3, results not shown), the peak current increases with Ny rapidly without saturation. This is consistent with our previous discussion since a thicker hBN tunneling barrier greatly suppresses 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 (Nx), 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.

Additional Information

How to cite this article: Zhao, Y. et al. Negative Differential Resistance in Boron Nitride Graphene Heterostructures: Physical Mechanisms and Size Scaling Analysis. Sci. Rep. 5, 10712; doi: 10.1038/srep10712 (2015).