Simulation of Figures of Merit for Barristor Based on Graphene/Insulator Junction.

We investigated the tunneling of graphene/insulator/metal heterojunctions by revising the Tsu-Esaki model of Fowler-Nordheim tunneling and direct tunneling current. Notably, the revised equations for both tunneling currents are proportional to V3, which originates from the linear dispersion of graphene. We developed a simulation tool by adopting revised tunneling equations using MATLAB. Thereafter, we optimized the device performance of the field-emission barristor by engineering the barrier height and thickness to improve the delay time, cut-off frequency, and power-delay product.

1. Introduction

Graphene barristor [GB] has been introduced to break the limitation of the low ION/IOFF of the graphene field-effect transistor (GFET) [1]. Since the introduction of the barristor, not only Si, but either organic [2,3,4,5] or inorganic [6,7,8,9,10,11,12] materials, have been exploited to make the graphene-semiconductor junction and applied on photosensors, gas sensors, etc., with its superior ION/IOFF.

Extreme temperature changes or degradation due to the environment may affect the performance of GFET or GB [13]. However, graphene-insulator-junction barristors exhibit a more stable performance compared to semiconductor-based ones [14].

In this study, we revised the models for Fowler–Nordheim tunneling (FNT) and direct tunneling (DT) in graphene/insulator/metal (GIM) junctions using the Tsu–Esaki tunneling model to reflect the graphene’s linear band structure [15]. Compared to the traditional FNT equation—proportional to V2—the revised FNT equation—proportional to V3—fits better with the experimental data. Then, we simulated how delay time (τ), power-delay product (PDP), and cut-off frequency (fT) of the field-emission barristor (FEB) could be improved by varying the tunneling-barrier height (), and the thickness of the hexagonal boron nitride (hBN) (tTunnel) with the revised tunneling models. We also considered thermionic emission for the low barrier height and DT for the thin insulator channel. The figures of merit of the FEB were extracted from the experimentally measured I–V characteristics in [14]. Notably, we could improve the device performance by decreasing tTunnel. This is because , followed by the channel current (ID), decreased with tTunnel. However, the improvement in the device performance by increasing tTunnel has a limitation because the increase in tTunnel deteriorates the PDP. Therefore, to not only improve τ and fT but also PDP, both the tunneling barrier height and tTunnel should be decreased.

2. Materials and Methods

We used the polydimethylsiloxane (PDMS) stamping method to create the bottom structure (gate electrode/hBN). Few-layer hBN was prepared on PDMS (PF Film-X4-6.5 mil bought in Gel-Film®) by mechanical exfoliation, and it was attached to the slide glass upside down. The hBN was subsequently aligned on top of the target substrate using a 3-axes manipulator. The PDMS was heated to 60 °C while it was transferred onto the target. Finally, PDMS was mechanically peeled off after the temperature was decreased to 25 °C. The top structure (hBN/graphene) was then transferred to the metal/hBN structure using the PMMA transfer method. The drain and source electrodes were deposited using an e-beam evaporator. The processes were conducted at the Core Facility Center for Quantum Characterization/Analysis of Two-Dimensional Materials and Heterostructures.

The simulation data were calculated using MATLAB R2021b. The source code implements the equations derived in the manuscript. We utilized MATLAB to calculate the current by varying the barrier height and the electric field. Using the double for loops, we obtained the current ranging from the height up to 3.0 eV and the field up to 0.4 V/nm with the discrete tunneling thickness. The source code is provided in the Supplementary Materials.

3. Results

3.1. Tunneling Current Model in GIM Junction

We derived a revised model for tunneling in graphene/insulator/metal (GIM) junctions. We started with the Tsu–Esaki model and applied the band structure of graphene. The current density tunneling from graphene to the metal can be calculated as follows [16]: where is the elementary charge; is the velocity in the x-direction, which is perpendicular to the graphene plane (to the insulator); is the transmission coefficient; and are the wave vectors of the carrier parallel and perpendicular to the graphene plane (y-z plane), respectively; and is the density of states of the graphene with momentum . In addition, and indicate the Fermi–Dirac distribution of the filled states in the graphene and empty states in the metal, respectively. The density of states can be obtained as , where and is the number of states per unit cell in the 2D momentum space. The areas of the primitive cell (S1) in real space and reciprocal lattice space (S2) of graphene are and , respectively (Figure S1a) , considering the spin and valley degeneracy [17,18]. Thus, , as shown in Figure S2b. Since graphene has a linear dispersion relation around the K-point, so that or (inset of Figure S2b), and the density of states in energy space , where is the Fermi velocity and is the reduced Plank’s constant.

Using the parabolic dispersion relation so that , , Equation (1) can be re-expressed with changed to as follows:

As the total energy , Equation (2a) can be re-expressed as:

Therefore, considering that ,

Notably, of the integrand originates from the linear dispersion relation of graphene, which does not exist in the original transport equation.

Equation (3) is a model of the tunneling current density in a graphene/insulator/metal junction. It comprises two integral parts: the first integral indicates the tunneling probability between graphene and the metal. Because it includes the tunneling coefficient T, which is related to the barrier height and width, it has different forms depending on the tunneling mechanisms such as FNT and DT [19,20,21,22]. The difference will be described in the next subsection. The second integral is related to the carriers supplied at the interface by an applied voltage. Thus, it has an identical form for both the FNT and DT. At the temperature T = 0, the integrand is finite only in the interval of because the carrier tunnels only from the filled state of graphene and the empty state of the metal, as shown in Figure 1b,c. Thus, the second integral can be simplified to , where is the Fermi level of the metal. Therefore, can be written as T = 0, as follows: where is the interval limit between the filled state of graphene and empty state of the metal, and is the Fermi level of graphene.

Figure 1

(a) A schematic of the field-emission barristor and schematic illustration of the graphene/insulator/metal (GIM) junction in Fowler–Nordheim tunneling; (b) an illustration of the FN tunneling. Electron tunnels through the insulator with the tunneling thickness of ; (c) an illustration of DT; electrons tunnel through the insulator with a width of d.

While the total current density is the difference between the current tunneling from graphene to metal () and that from metal to graphene (), the total current density of the GIM junction becomes because does not contribute to . This is because becomes zero because under .

3.2. Fowler–Nordheim Tunneling Model in GIM Junction

In Equation (4), the transmission coefficient for FNT in the Wentzel–Kramers–Brillouin approximation can be expressed as [16,19], where m∗ is the electron effective mass, the conduction band minimum (CBM) , and , as shown in Figure 1b. Assuming that the work function of graphene and metal are the same, , where V is the potential bias between the graphene and metal. Then, the transmission coefficient can be rewritten as follows:

Since the first term in the curly bracket disappears because of the value of x1, the tunneling coefficient can be described as follows:

Then, Equation (4) can be rewritten for the FNT current when E is between , and as follows:

The exponential term in Equation (7) can be rewritten by using the following Taylor series as .

Then, by substituting with E, the FNT current is given by the following equation:

The integral part can be simplified as where and .

Assuming that , the k is very small, as is Ak. Therefore, only the constant in Equation (8b) survives. Equation (8a) can then be rewritten as follows:

Equation (9) represents a new model for the FNT current in a graphene/insulator/metal junction [15]. Notably, in the revised model, the FNT is proportional to V3, whereas the FNT in the Tsu–Esaki model is proportional to V2.

3.3. Direct Tunneling Model in GIM Junction

T(Ex) for DT can be expressed from Equation (4) as follows [23]:

Note that the integration interval for DT is [0, d], as shown in Figure 1c. After integrating and taking the terms up to the order of of the Taylor series, the current density for DT can be rewritten as follows: where E is located between and , , and is replaced by E. Therefore, the DT current can be rewritten as follows:

Equation (12) is a revised DT model for a graphene/insulator/metal junction. Compared to the DT equation based on the Tsu–Esaki model in the metal/insulator/metal junctions [15], the tunneling current in the Tsu–Esaki model is proportional to V2, whereas the current in the revised model is proportional to V3.

3.4. FNT Barrier Height

To calculate the barrier height using the FNT equation for the graphene/insulator junction, the FNT equation can be rewritten as follows: where ID is the drain current; VD is the drain voltage, and . In the new equation, γ is replaced with from the original FNT equation because they are extracted from different second integrations in Equation (3), which are related to the density of states in the metal and graphene, respectively. However, β does not change because it is extracted from the tunneling coefficient, which is the first integration in Equation (3). Consequently, the results of the barrier height calculation using the original and revised FNT equations were not significantly different.

Figure 2a shows the replotted ID–VD curve of the FEB using the revised FNT equation. The graph fits well to the straight line, and its slope was estimated to be −683 V. Figure 2b shows a replotted ID–VD curve obtained using the traditional FNT equation. The slope of the straight line was estimated to be −689 V. The tunneling barrier heights extracted from the slopes were 2.10 eV (by the revised FNT equation) and 2.11 eV (by the traditional FNT equation), respectively. Therefore, the results of the barrier height calculation using the original and revised FNT equation were similar because they used an identical β [19,20,21,22,23,24].

Figure 2

Fitting of FNT to calculate the tunneling barrier height with the revised and traditional equations. (a) Straight line (red) fitted to FNT current by the revised FNT equation. Its slope is estimated to be −683 V. A barrier height of 2.10 eV was extracted from the slope. (b) Straight line (blue) fitted to the FNT current by the traditional FNT equation. Its slope was estimated to be −689 V. A barrier height of 2.11 eV was extracted from the slope. (c) The experimental ID–VD curve (black) of FEB consisting of graphene and hBN, the simulated I–V curve by the revised FNT equation (red), and simulated ID–VD curve by the original FNT equation (blue). Further information of fitting method is explained in detail on Appendix A.

However, the discrepancy between the two models became apparent when we estimated the FNT current in the graphene/insulator junctions. In contrast to the barrier calculation, α and γ affected the tunneling current estimated using each FNT equation. Figure 2c shows the experimental data and simulated ID–VD curves obtained using the original (blue) and revised (red) FNT equations. The black curve indicates the measured ID–VD. The red line was estimated using the revised FNT equation: . The blue line was obtained by the original FNT equation: [25]. The barrier heights of 2.10 eV and 2.11 eV were applied to the revised and traditional FNT equations, respectively. Although the barrier heights were similar, the simulated currents were significantly different because of the and . As shown in Figure 2c, the red curve simulated by the revised FNT equation was better fitted to the experimental data. Therefore, we used the revised FNT equation to simulate the figures of merit for the FEB.

3.5. Simulation for Barrier Height Engineering to Improve Delay Time and Cut-Off Frequency

To evaluate the performance of the FEB, we extracted τ and fT from the experimental I–V curve, where τ is a time delay required to charge the gate electrode with ION, and fT is a maximum frequency up to which the current of the transistor could be amplified [14]. We obtained a delay time of 154 ns and a cut-off frequency of 13.8 MHz for the FEB. Compared to the performance of the graphene/Si barristor, which was simulated using NanoTCAD ViDES (Device simulator) [26], the FEB’s delay time was 140 times slower than that of the graphene/Si barristor (1.1 ns), and the cut-off frequency was 92 times lower than that of the graphene/Si barristor (1.3 GHz). The low performance of the FEB originated from the low ON current (ION) because the delay time and cut-off frequency depend on ION. In the FEB, the tunneling barrier height should be decreased to improve ION. Therefore, to obtain the minimum barrier height in our device, we estimated the drain current by varying the barrier height and electric field strength between the source and drain electrodes (EField) using the revised FNT equation.

Figure 3a shows the FNT current simulated by varying and EField. To increase the JON, a higher EField and a lower are required. However, Figure 3b describes that ratios (slope of lines) decreased with EField. This indicates that the device requires more charge on graphene to modulate its work function. Increasing EField to improve JON requires more switching energy for the device. Therefore, the maximum EField should be determined by considering both the on-state current and the energy consumption for switching. Likewise, decreasing the tunneling barrier height is limited by thermionic emission. The thermionic emission current was estimated using the following equation: where kB is the Boltzmann constant and T is the temperature [27].

Figure 3

The simulated FNT current as a function of the tunneling barrier height and drain electric field. (a) A 3D plot of the FNT equation. (b) The channel current under a low drain electric field (0.1 V/nm) can be affected by the thermionic emission current. The ΔJD/Δ ratios decreased with increasing drain electric field. The on-state current Jon was 10−4 A/μm2 when was 0.5 eV, and the off-state current Joff was 10−10 A/μm2 when was 1.055 eV. (c) A 3D plot of the DT equation for different thicknesses of the insulator. The DT current exponentially increased with a decrease in the thickness. (d) The simulated FNT current and DT as a function of the tunneling barrier height for different thicknesses of the insulator. The tunneling current under EField = 0.2 V/nm was simulated by using the revised FNT and revised DT equations. The DT current increased with the decreasing thickness of the insulator. Jon and Joff indicate the on-state current and off-state current, respectively.

Because the thermionic emission current exponentially depends on and the temperature, the total channel current under EField lower than 0.1 V/nm is affected by the thermionic emission current, as shown in Figure 3b. Therefore, an EField above 0.1 V/nm should be applied to avoid temperature dependence of the ID–VD characteristics. As shown in Figure 3b, these conditions for improving the delay time and cut-off frequency were satisfied when the EField was 0.2 V/nm, JON was 10−4 A/μm2 with 0.5 eV, and JOff was 10−10 A/μm2 with 1.055 eV. The required charge () to decrease from 1.055 eV to 0.505 eV was calculated by using the equation: , where is the work function shift of graphene [28]. The delay time and cut-off frequency were calculated using the following equations [29]:

We obtained a delay time τ of 0.18 ns and cut-off frequency fT of 3.98 GHz that were better than 1.1 ns and 1.3 GHz in the graphene/Si barristor. Therefore, when the applied EField was 0.2 V/nm and was changed from 0.5 eV to 1.055 eV or vice versa, FEB could achieve the optimized delay time and cut-off frequency.

3.6. Simulation for Insulator Thickness Engineering to Improve Power-Delay Product

The power delay product (PDP) refers to the energy consumed during device switching. We estimated the PDP of the FEB to be 355 fJ/μm2 from the ID–VG graph, which was 47 times greater than that of the graphene/Si barristor (7.5 fJ/μm2). This high energy consumption for device switching is because the semiconductor-less device requires a high VD to control the FNT current. Therefore, to reduce the PDP, tTunnel, where the FNT takes place, should be reduced.

Figure 3d exhibits the DT current by varying tTunnel. When tTunnel was thinner than 2 nm, the FNT could not control the channel current because the DT current dominated the FNT current. In contrast, DT was suppressed below the off-state current in the FNT regime when tTunnel was 6 nm. Therefore, when tTunnel was 6 nm, and EField was 0.2 V/nm, we obtained a PDP of 21 fJ/μm2 from the following equation:

Although still greater than that of the graphene/Si barristor (7.4 fJ/μm2), the PDP decreased to 6% of the measurement by engineering tTunnel.

4. Conclusions

In conclusion, we introduce a new model for the FNT and DT currents of graphene/insulator/metal (GIM) heterojunctions. We obtained new models by revising the supply function of the Tsu–Esaki model. Notably, the tunneling current in the revised FNT equation was proportional to V3. We then extracted the tunneling barrier height in the graphene/insulator junction from the slope of the line in the ID–VD curve replotted with the axes of ln(I/V3) and 1/V. The barrier height obtained using the revised model was not significantly different from that of the original model. However, the ID–VD curve estimated by the revised FNT fit better with the experimental data than the ID–VD curve simulated by the original FNT. Then, we simulated the τ, fT, and PDP of the FEB by varying and tTunnel by using the revised FNT equation. These significantly improved by decreasing and increasing EField. By considering the thermionic emission at low barrier height and energy consumption at the high electric field, we obtained a τ of 0.18 ns and fT of 3.98 GHz when EField was 0.2 V/nm, and was changed from 0.5 to 1.055 eV or vice versa. We improved the PDP by decreasing the tTunnel. As DT exponentially increased as the thickness decreased, we obtained a lower boundary tTunnel of 6 nm, and then the PDP decreased to 17 times lower than the experimental data.