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.

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 V 2 -the revised FNT equation-proportional to V 3 -fits better with the experimental data. Then, we simulated how delay time (τ), power-delay product (PDP), and cut-off frequency (f T ) of the field-emission barristor (FEB) could be improved by varying the tunneling-barrier height (∅ B ), and the thickness of the hexagonal boron nitride (hBN) (t Tunnel ) 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 t Tunnel . This is because ∅ B , followed by the channel current (I D ), decreased with t Tunnel . However, the improvement in the device performance by increasing t Tunnel has a limitation because the increase in t Tunnel deteriorates the PDP. Therefore, to not only improve τ and f T but also PDP, both the tunneling barrier height and t Tunnel should be decreased.

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 Nanomaterials 2022, 12, 3029 2 of 10 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.

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 q is the elementary charge; v x is the velocity in the x-direction, which is perpendicular to the graphene plane (to the insulator); T (k x ) is the transmission coefficient; k g and k x are the wave vectors of the carrier parallel and perpendicular to the graphene plane (y-z plane), respectively; and D g k g is the density of states of the graphene with momentum k g . In addition, f g (E) and 1 − f m (E) indicate the Fermi-Dirac distribution of the filled states in the graphene and empty states in the metal, respectively. The density of states D g k g can be obtained as D g k g dk g = D * g k y , k z dk y dk z k=k g , where k = k x 2 + k y 2 and D * g k y , k z is the number of states per unit cell in the 2D momentum space. The areas of the primitive cell (S 1 ) in real space and reciprocal lattice space (S 2 ) of graphene are S 1 = √ 3 2 a 2 and S 2 = 8π 2 √ 3a 2 , respectively (Figure S1a) D * g k y , k z = 2 * 2 * 1 S 1 * 1 S 2 = 1 π 2 , considering the spin and valley degeneracy [17,18]. Thus, D g k g dk g = 1 π 2 2πk g dk g = 2 π k g dk g , as shown in Figure S2b. Since graphene has a linear dispersion relation around the K-point, so that E g = v F k g or v F dk g = 1 dE g (inset of Figure S2b), and the density of states in energy space where v F is the Fermi velocity and is the reduced Plank's constant. Using the parabolic dispersion relation so that (1) can be re-expressed with k x changed to E x as follows: As the total energy E = E x + E g , Equation (2a) can be re-expressed as: Therefore, considering that E g > 0, 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 E x , E F,g 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 E F,m is the Fermi level of the metal. Therefore, J g→m 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 E F,g is the Fermi level of graphene.
Therefore, considering that > 0, 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 ∫ ( − ) , = 1 2 ( , − ) 2 , 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. 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

Fowler-Nordheim Tunneling Model in GIM Junction
In Equation (4) While the total current density is the difference between the current tunneling from graphene to metal (J g→m ) and that from metal to graphene (J m→g ), the total current density J of the GIM junction becomes J g→m because J m→g does not contribute to J. This is because

Fowler-Nordheim Tunneling Model in GIM Junction
In Equation (4), the transmission coefficient T(E x ) 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) , as shown in Figure 1b. Assuming that the work function of graphene and metal are the same, V = (E F,g − E F,m )/q, 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 x 1 , the tunneling coefficient can be described as follows: Then, Equation (4) can be rewritten for the FNT current when E x is between E F,m , and E F,g as follows: The exponential term in Equation (7) can be rewritten by using the following Taylor series as Then, by substituting E x − E F,g with E, the FNT current J is given by the following equation: The integral part can be simplified as Assuming that E F,g E F,m , 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 V 3 , whereas the FNT in the Tsu-Esaki model is proportional to V 2 .

Direct Tunneling Model in GIM Junction
T(E x ) 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 (q∅ B ) 3 2 of the Taylor series, the current density for DT can be rewritten as follows: where E x is located between E F,m and E F,g , qV q∅ B , and E x − E F,g 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 V 2 , whereas the current in the revised model is proportional to V 3 .

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 I D is the drain current; V D is the drain voltage, γ = ln 8πh∅ B d 2 m * 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 I D -V D 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 I D -V D 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].
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 ( ∅ ) 3 2 ⁄ of the Taylor series, the current density for DT can be rewritten as follows: where Ex 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 V 2 , whereas the current in the revised model is proportional to V 3 .

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, =  (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].  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 I D -V D curves obtained using the original (blue) and revised (red) FNT equations. The black curve indicates the measured I D -V D . 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.

Simulation for Barrier Height Engineering to Improve Delay Time and Cut-Off Frequency
To evaluate the performance of the FEB, we extracted τ and f T from the experimental I-V curve, where τ is a time delay required to charge the gate electrode with I ON , and f T 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 (I ON ) because the delay time and cut-off frequency depend on I ON . In the FEB, the tunneling barrier height should be decreased to improve I ON . 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 (E Field ) using the revised FNT equation. Figure 3a shows the FNT current simulated by varying ∅ B and E Field . To increase the J ON , a higher E Field and a lower ∅ B are required. However, Figure 3b describes that ∆J D /∆∅ B ratios (slope of lines) decreased with E Field . This indicates that the device requires more charge on graphene to modulate its work function. Increasing E Field to improve J ON requires more switching energy for the device. Therefore, the maximum E Field 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 k B is the Boltzmann constant and T is the temperature [27].
Because the thermionic emission current exponentially depends on ∅ B and the temperature, the total channel current under E Field lower than 0.1 V/nm is affected by the thermionic emission current, as shown in Figure 3b. Therefore, an E Field above 0.1 V/nm should be applied to avoid temperature dependence of the I D -V D characteristics. As shown in Figure 3b, these conditions for improving the delay time and cut-off frequency were satisfied when the E Field was 0.2 V/nm, J ON was 10 −4 A/µm 2 with ∅ B 0.5 eV, and J Off was 10 −10 A/µm 2 with ∅ B 1.055 eV. The required charge (Q) to decrease ∅ B from 1.055 eV to 0.505 eV was calculated by using the equation: where ∆w G is the work function shift of graphene [28]. The delay time and cut-off frequency were calculated using the following equations [29]: 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/μm 2 with ∅ 0.5 eV, and JOff was 10 −10 A/μm 2 with ∅ 1.055 eV. The required charge ( ) to decrease ∅ from 1.055 eV to 0.505 eV was calculated by using the equation: ∆ = ℎ 2 √ , 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 f T of 3.98 GHz that were better than 1.1 ns and 1.3 GHz in the graphene/Si barristor. Therefore, when the applied E Field was 0.2 V/nm and ∅ B was changed from 0.5 eV to 1.055 eV or vice versa, FEB could achieve the optimized delay time and cut-off frequency.

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/µm 2 from the I D -V G graph, which was 47 times greater than that of the graphene/Si barristor (7.5 fJ/µm 2 ). This high energy consumption for device switching is because the semiconductor-less device requires a high V D to control the FNT current. Therefore, to reduce the PDP, t Tunnel , where the FNT takes place, should be reduced. Figure 3d exhibits the DT current by varying t Tunnel . When t Tunnel 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 t Tunnel was 6 nm. Therefore, when t Tunnel was 6 nm, and E Field was 0.2 V/nm, we obtained a PDP of 21 fJ/µm 2 from the following equation: Although still greater than that of the graphene/Si barristor (7.4 fJ/µm 2 ), the PDP decreased to 6% of the measurement by engineering t Tunnel .

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 V 3 . We then extracted the tunneling barrier height in the graphene/insulator junction from the slope of the line in the I D -V D curve replotted with the axes of ln(I/V 3 ) and 1/V. The barrier height obtained using the revised model was not significantly different from that of the original model. However, the I D -V D curve estimated by the revised FNT fit better with the experimental data than the I D -V D curve simulated by the original FNT. Then, we simulated the τ, f T , and PDP of the FEB by varying ∅ B and t Tunnel by using the revised FNT equation. These significantly improved by decreasing ∅ B and increasing E Field . 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 f T of 3.98 GHz when E Field was 0.2 V/nm, and ∅ B was changed from 0.5 to 1.055 eV or vice versa. We improved the PDP by decreasing the t Tunnel . As DT exponentially increased as the thickness decreased, we obtained a lower boundary t Tunnel of 6 nm, and then the PDP decreased to 17 times lower than the experimental data.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/nano12173029/s1, Figure S1: An illustration of the reciprocal lattice space of graphene; Figure S2: The optical microscope image of the device; Supplementary Software Files: "S1 MATLAB code.zip" simulation MATLAB Code.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The fitting lines projected in Figure 2a,b were obtained by linear fitting of each plotted data. The interval of the x-axis for the linear regression ranged from 0.034 V −1 to 0.039 V −1 , corresponding to the V D of 25.6 V and 29.0 V, respectively. This limited interval was to distinguish the FNT current from that of the DT. Then, the currents in Figure 2c were extrapolated by using each current equation (revised or traditional) with the key parameters obtained by the linear regression. Figure 2c projects the data ranges from 20 to 30 V.