The authors present an overview of the main theoretical approaches used to describe tunnel processes in graphene nanoelectronics. Two currently central theoretical methods of calculating tunnel current, the Bardeen tunneling Hamiltonian approach and the method of nonequilibrium Green's functions, are reviewed in a pedagogical fashion. Several examples are used to illustrate the specific features of the methods. An application of both methods to the analysis of current flow in graphene tunnel field-effect transistors is discussed.

## I. INTRODUCTION

Present-day progress in the field of microelectronics follows the way of improving the technique of creating electronic components with the aim of their miniaturization, which increases performance and reduces energy consumption of electronic circuits on their basis. Moreover, a transition from the usual silicon electronics to a radically new component basis is expected on the basis of modern nanomaterials. This is evidenced by the development and production of prototypes of a new generation of diodes, transistors, microchips, sensors, nanoelectromechanical devices, etc. It should be noted that of particular interest is a possibility to apply new materials with dimension less than three. Among these are 2D crystalline materials such as graphene, chalcogenides of transition metals, phosphorene, silicene, and others; 1D materials are presented by carbon nanotubes, metal nanowires, and linear polymers. Examples of zero-dimensional structures are fullerenes, metallic, and semiconductor clusters.

The unique properties of graphene make it possible to consider it as one of the basic materials for new electronics.^{1} Indeed, exceptionally high mobility of charge carriers assumes a possibility to increase quick operation of electronic devices. Monoatomic thickness in combination with high stability and flexibility provides good processibility of synthesis and laying of graphene on different substrates. Considerable progress in the production technique of the defectless, perfect crystal graphene allows one to expect creation of reproducible prototypes of effective electronic devices.

One of the priority directions is the development of nanodevices whose physical principle of operation is based on the use of the tunnel effect. We are mostly interested in technological achievements and theoretical ideas of fabrication of a graphene-based current switching tunnel device. The most well-known results in this direction are related to the elaboration of a vertical tunnel field-effect transistor (TFET) on the basis of heterostructures.^{2,3} This transistor has a configuration of a parallel-plate capacitor with an additional electrode playing the role of a gate [see Fig. 1(a)].

It should be emphasized that the concept of tunneling between parallel two-dimensional electron systems (2DEG) has been developing over the last 25 years. In the first experiments, the tunneling conductance between two 2DEGs in semiconductor heterostructures AlGaAs/GaAs was measured at different temperatures, carrier densities, and applied magnetic fields.^{4–8} The 2D-2D tunneling is a resonance process that is very sensitive to the change of the band structure under the action of external electric and magnetic fields as well as to the change of the relaxation time of charge carriers. This makes it possible both to obtain a noticeable switching effect and study specific features of the electronic spectrum in the material playing the role of coatings. For example, in Ref. 9, measuring quantum capacity of the graphene heterostructure revealed asymmetry between electrons and holes in graphene. This allowed one to determine the value of the overlap integral between distant neighbors and show that the experimental value significantly exceeds theoretical estimates. At the present time, active work is underway on the synthesis and investigation of heterostructures on the basis of graphene with the use of various materials as dielectric layers. One of the new interesting results is the realization of the light emitter on the graphene heterojunction base.^{10}

Theoretical studies of the tunnel heterostructure on the basis of graphene were carried out in Ref. 11, where a resonance nature of tunnel current was demonstrated and nonlinearity of current–voltage characteristics was revealed. Later, the concept of symmetric field tunnel transistor (SymFET) was proposed.^{12} In that device, the tunnel current flows from a graphene electrode of the *n*-type to an electrode of the *p*-type. It reaches a maximum value when the alignment of Dirac's points occurs at a certain voltage between the electrodes. The characteristics of the device were shown to slightly depend on the temperature, and the resonance peak of the tunnel current is well controlled by the doping level and the gate voltage.

An alternative configuration of thin-film electrodes is a planar arrangement where they are separated by a narrow gap (nanogap). One of the first experimental attempts to create a device on the basis of the nanogap is presented in Ref. 13. It is worth noting that the creation of a rather small gap for a direct electrode–electrode tunneling is a highly complicated task. In Ref. 13, a gap of 100 nm was created, which made it possible to observe only the effect of field emission that is weakly amendable to the change by the gate voltage. In Ref. 14, a tunnel contact between the sheets of the multilayer graphene in the planar configuration was produced. The edge of the adjacent sheets turned out to be uneven, and the gaps had closest separations on the order of 1–2 nm. In this case, tunneling proceeded between those protrusions on the contacts which appeared to be closer to each other.

An original theoretical concept of a switching device on the basis of a plain configuration of graphene electrodes is formulated in Ref. 15. A distinctive feature of this device is that the channel of a graphene FET is formed by hexagonal boron carbon-nitrogen (hBCN) domains. The possibility of influencing electronic properties of the channel with hBCN allows the development of very efficient graphene-based transistors.

An interesting version of the transistor in the planar configuration was proposed in Ref. 16. It should be noted that a common difficulty of using graphene without the vacuum or heterostructural barrier is caused by that the main charge carriers are massless quasirelativistic fermions with two quasispin states. This determines the presence of the Klein paradox in graphene^{1} which manifests itself in the fact that tunnel barriers inside graphene are transparent for electrons that are incident on them at a right angle. This fact makes the control of the current in graphene by the external electrostatic field difficult. The original idea of the authors of Ref. 16 consists in separating the source/drain electrodes by the saw-shaped geometry gate (in place of the conventional bar-shaped geometry).

Another concept for a graphene TFET was proposed in our papers.^{17–19} The source and drain electrodes were suggested to be made of the sheets of monolayer graphene on the insulator substrate in one plane oriented toward each other by crystallographically matching edges of the zigzag type and separated by a tunnel-transparent, for charge carriers, vacuum barrier [see Fig. 1(b)]. The switching principle of this device essentially differs from that suggested in Ref. 15. In our case, the main idea is based on the use of specific localized electronic states, edge states, which are resistant to external perturbations due to their being “topologically protected.” They provide increased density of states (DOS) near the Fermi energy. The gate voltage changes the position of the regions with increased DOS with respect to an admissible energy interval (the energy window) within which an electronic transport is possible, which ensures the switching effect.

For practical implementation of this kind of transistor, one needs to create a rather narrow gap between the electrodes whose edges preserve a correct crystallographic configuration of the zigzag type. The set of available technologies makes it possible to produce gaps of width of an order of tens of nanometers, and attempts are made at their further reduction.^{20} This is very important not only for the production of transistors but also other urgent devices, for example, such as DNA sequenator of a new generation.^{21,22} Of the advanced technological methods of creating a nanogap with the required parameters, one should distinguish a nanolithography with the help of the scanning tunnel microscopy, electromigration, local anode oxidation, nanofabrication with the use of transmission electron microscopy, and catalytic nanocutting.

The work is organized as follows. Section II presents briefly the necessary aspects of Bardeen's theory and Green's functions approach. Three-plate capacitor model is considered in Sec. III. The fact is that all known promising graphene tunneling devices have a gate, so that electrostatics of such structures is of importance. In Sec. IV, we present a solution of a number of interesting problems where either one or both approaches are used. In particular, examples in Secs. IV A and IV B illustrate the effect of changing the tunneling area from 0D to 2D. In Sec. IV C, the tunneling problem is examined with the example of armchair-edged graphene tunnel contacts where both approaches are applied. Sections IV D and IV E are devoted to the description of two examples of graphene TFETs based on graphene-insulator-graphene (GIG) structures: the vertical and the planar ones. The final section discusses the recent literature where both theoretical methods were effectively used to describe tunnel current in GIG structures.

## II. CALCULATION METHODS

Tunneling is one of the brightest manifestations of quantum effects. Josephson, Esaki, and Giaever were awarded the Nobel Prize in 1973 for the discoveries regarding tunneling phenomena in solids. The design of the scanning tunneling microscope (STM) made a revolution in the physics of meso- and nanoscale phenomena for which Binnig and Rohrer were awarded the Nobel Prize in 1986. We will discuss the key points of two basic approaches used to describe tunnel processes in various structures: the Bardeen transfer (tunneling) Hamiltonian approach (TH) and the method of nonequilibrium Green's functions (NEGF).

### A. Transfer Hamiltonian approach

The method of tunnel Hamiltonian^{23–25} was developed simultaneously with the experimental methods of tunnel spectroscopy.^{26} At present, this approach is actively used for theoretical description of tunneling in various 2D heterostructures including on the basis of graphene. The conventional TH theory uses the following expression for the tunnel current:^{27}

with **k** and **q** being the wave vectors of charge carriers in the left and right contacts, respectively, $ f ( \epsilon ) = [ 1 + exp ( \epsilon / k B T ) ] \u2212 1 $ the Fermi-Dirac distribution function, $ M ( k , q ) $ the matrix element describing the interaction between contacts, and *V _{b}* the bias voltage. Spectral densities can be written down as the Lorentz functions

where $ S = ( L , R ) , \u2009 E S ( k ) $ is the energy spectrum, $ \gamma S = \u210f / \tau S $, and *τ ^{S}* is the relaxation time of electrons in the contacts. Whenever inelastic processes in the tunnel contacts can be neglected, spectral functions take the form of delta-functions: $ lim \gamma \u2192 0 A ( k , \epsilon ) = 2 \pi \delta ( E ( k ) \u2212 \epsilon ) $.

Bardeen proposed the following expression for the matrix element:^{23,25}

where *ψ ^{L}* and

*ψ*are the electron wave functions in the contacts and integration is over an arbitrary surface separating contacts. Notice that the functions

^{R}*ψ*and

^{L}*ψ*can be used in various approximations and below we give examples of the use of WKB and local combination of atomic orbitals (LCAO) approximations and draw an analogy between them.

^{R}#### 1. WKB wavefunctions: Harrison's approach

Harrison was among the first to use this approach in solving the problem of tunneling between two electrodes separated by a thin dielectric film in the case of different effective mass of electrons on the right and left.^{24} In taking translation symmetry into account in the directions perpendicular to tunneling the problem reduces to the one-dimensional one. The main approximation used in Ref. 24 is the assumption of a smooth change of the parameters defining the band structure when passing from the left to the right contact. This makes it possible to use the WKB method to calculate the electron wave functions in the contacts. As a result, the square of the matrix element calculated by Eq. (3) takes the following form:

where $ \nu S = v g S / 2 l S $ determines the so-called attempt frequency,^{28}, $ v g S $ are group velocities of electrons in the contacts having an extension *l _{S}*. Equation (4) contains the exponent typical for the WKB approximation;

*x*and

_{L}*x*fix the width of the tunnel barrier. This result is physically understood: electron transmittance is proportional to the frequency of electron collision with the boundaries and the probability of tunneling through a barrier dividing them. This approach was further developed in Refs. 29 and 30. Notice that the method given in Ref. 24 allows one to illustrate qualitatively the role of the effective mass in tunnel phenomena.

_{R}#### 2. LCAO approximation

The method of LCAO, which is closely related to the tight binding (TB) method, can be considered in a sense as an alternative to the foregoing approach. Within the Harrison approach, the tunnel contacts are simulated by potential wells and their real structure is considered only by introducing the effective mass which is smoothly changing from one well to another. In the framework of the LCAO approximation, the electron wave function is taken to be a superposition of orbitals of all atoms in the structure

Here, the information on the structure is contained in a set of coefficients *c _{n}*, atomic positions $ r n $, and in the form of orbital wave functions $\varphi $. With Eq. (5) taken into account, the matrix element Eq. (3) becomes

where only orbital wave functions remain under the integral. Let us determine the interaction between the orbitals, *t _{nm}*, through

To calculate the coefficients *t _{nm}*, it is necessary to define the explicit form of the orbital wave functions of interacting atoms. Various versions of the wave functions for

*s*and

*p*orbitals and, correspondingly, for

*t*were proposed in Refs. 31–34. It should be noted that in all the cases $ t n m \u223c \u2009 exp ( \u2212 \lambda | r m \u2212 r n | ) $, i.e.,

_{nm}*t*exponentially decreased with increasing distance between the atoms but the pre-exponential factors had characteristic differences.

_{nm}It is worth to note that the TH approach is well suited to a qualitative analysis and revealing of the basic laws in studying the tunnel current. The algorithm of its application can be described as follows: (1) Hamiltonians of tunnel contacts are formulated, and the problem of eigenvectors and eigenvalues is solved, (2) eigenvectors are used to calculate the square of the matrix element describing the interaction of contacts with each other, (3) the matrix element together with the spectrum of charge carriers in the contacts is used for calculation of the current according to Eq. (1). Among the drawbacks of this kind of approach is the difficulty of describing tunnel contacts which influence each other so strongly that their DOS spectra are changed. From a formal point of view, such “strongly interacting” contacts can be described within the given method only if the proper problem for the general Hamiltonian of two contacts is solved. This offers serious problems in practice. However, we show below that neglecting the interference of contacts one can obtain, in some cases, a result analytically.

### B. Nonequilibrium Green's function method

At present, the NEGF method is in great demand in calculating the tunnel current. The fact is that this approach is best suited to numerical calculations because the solution of the transport problem requires only a correct definition of the total Hamiltonian of a system. All further actions are simply formalized, and finally, the calculation of tunnel characteristics entirely depends on the computational power.

The total Hamiltonian of the tunnel problem can be represented as a sum of the Hamiltonians of the left and right contacts as well as an operator responsible for tunneling

In the second quantization representation for identical atoms, one has

where $ c i \u2020 $ and *c _{j}* are the electron creation and annihilation operators on the lattice sites

*i*and

*j*, $ \varphi L = e V b / 2 , \u2009 \varphi R = \u2212 e V b / 2 $,

*V*is the gate voltage (when the gate is present),

_{g}*t*and $ t i j n $ are the hopping integrals for atoms belonging to one and neighboring contacts, respectively.

_{ij}The tunnel current is defined as the rate of change of the particle number (N) in any contact, for example, in *R*

The Liouville equation with the Hamiltonian Eq. (8) gives

As a result, the current is expressed by the following known relation

Notice that each term in this sum defines the current between atoms *i* and *j* of different contacts. Therefore, Eq. (13) is convenient to use when describing spatial characteristics of the tunnel current. Moreover, it can be rewritten in terms of the Fourier transform of the Keldysh Green's function^{35}

where $ G i j K ( t ) = \u2212 i \u3008 [ c i ( t ) , c j \u2020 ( 0 ) ] \u3009 $, and in the last equality, the matrix notation is used. This equation was successfully used in the description of spatial current patterns in graphene nanoribbons.^{36} In the case of tunnel process, the Keldysh function can be expressed in terms of matrix *T*, Fermi functions *f _{s}*, the retarded $ G S + ( \epsilon ) $, and advanced $ G S \u2212 ( \epsilon ) $ Green's functions for the isolated contacts at thermal equilibrium through the solution of the corresponding Dyson's equation.

^{35,37}Finally, we get

Here, $ f S = f ( \epsilon + \varphi S ) $, and $1$ is the unit matrix. Due to the exponential decrease of the coefficients $ t i j n $ with distance, electron tunneling will proceed only from the nearest atoms in the contacts. Therefore, the matrices *G*, defining the total structures, can be replaced by the matrices *g*, corresponding to the local regions of the tunnel contacts of an adequate dimension. This type of Green's functions is called the surface ones (SGF). As a result, the expression for the tunnel current assumes the following form:

where $ A S = i ( g S + \u2212 g S \u2212 ) $ is the spectral density of electrons, and *t* is the interaction matrix for atoms in the local regions of the contacts. It should be noted that the second terms in parentheses of Eq. (15) take into account the effect of one contact on the electron states of the other. If the effect is weak, these terms can be neglected, and then (15) becomes identical to Eq. (1) within LCAO approximation.

A method for analyzing the quantum transport in the electrode-quantum dot (QD)-electrode system [see Fig. 2(a)] is well known. The current is written as^{38}

where $ G D \xb1 $ is the Fourier transform of corresponding Green's functions for a QD, $ \Gamma S = i ( \Sigma S \u2212 \Sigma S \u2020 ) $ is the tunneling rate from the left or right contact to the QD, $ \Sigma S = t S g S + t S \u2020 $ is the self-energy, *t _{S}* is the matrix describing the interaction of contacts with the QD, and $ g S + $ are the retarded SGF. Let the QD Hamiltonian be

*H*, then

_{D}^{38}

where $ \Sigma L , R $ describe the influence of the contacts on the QD.

The present method can be adapted for calculation of the current in the tunnel contact. As is shown in Fig. 2(b), the tunnel contact can be treated as an electrode-QD-electrode system if a small boundary part of one of the contacts is chosen as a QD. The resultant expression for the current will be identical to Eq. (15).

Thus, the calculation of the current is reduced to that of the SGF of the tunnel contacts. In deriving an equation for determining SGF, it is convenient to apply the following procedure. First, we divide the contact into the same recurring blocks uniformly interacting with each other [Fig. 2(c)]. The size and shape of the blocks depend on a concrete problem. The first block should include all atoms interacting with the QD. Let us denote this block by *H*_{00}, and the matrix of interaction of the blocks with each other by *H*_{01}. Now we consider the extreme block of the semi-infinite contact [Fig. 2(c)] as a QD standing separately. The SGF for it is $ g 00 + = [ ( \epsilon + i 0 + ) 1 \u2212 H 00 \u2212 \Sigma 01 ] \u2212 1 $. For $ \Sigma 01 $, the role of the interaction matrix will be played by *H*_{01}, and of the SGF by the function identical to *g*_{00}, as the characteristics of the edge of the semi-infinite structure will not be changed if it is lengthened or shortened to a finite number of blocks. As a result, we obtain the matrix equation for *g*_{00}

## III. THREE-PLATE CAPACITOR MODEL

All known concepts of graphene-based transistors include a gate electrode. A change of potential on the gate provides a switching effect. The applied voltage modifies the chemical potential of graphene *μ* which, in turn, affects the tunnel current. To take this into account let us consider a simple electrostatic model of a three parallel-plate capacitor device with the metal gate and a couple of graphene electrodes (see Fig. 3 where the corresponding band diagram is also shown).

Let us define that *μ* is positive/negative when the level is shifted to the conduction/valence band, respectively. The energy spectrum of graphene electrodes will be accordingly modified as $ E ( k ) \u2192 E ( k ) \u2212 \mu $. Our task is to calculate the values of *μ* with known *V _{g}* and

*V*.

_{b}We assume that the *x*-axis is directed normal to the plates. Applying Gauss's law around each plate, one obtains

where *F _{g}* and

*F*are the electric fields between plates 1–2 and 2–3, respectively. Keeping the charge neutrality condition $ \sigma 1 + \sigma 2 + \sigma 3 = 0 $ gives

_{b}On the other hand, the electric field in such GIG structure can be found by the potential difference across the plates and the distance between them: $ F b = \u2212 ( \varphi 3 \u2212 \varphi 2 ) / d 23 , \u2009 F g = \u2212 ( \varphi 2 \u2212 \varphi 1 ) / d 21 $. In turn, the electrode potential is the sum of the voltage and the electrochemical potential *μ*, so that

The electrochemical potential is determined by the charge density. In the case where the first electrode operates as a metal gate, *μ*_{1} may be taken to be zero. Thus, finally one obtains a system of equations

with *σ*_{2} and *σ*_{3} being unknown quantities. Notice that an alternative system with variables *μ*_{2} and *μ*_{3} can also be formulated. Using the approximate expression $ F g = \u2212 V g / d 21 $ reduces the problem to a single equation,^{3,41} while in Refs. 42 and 43, the system (23) was solved without any approximation.

The induced charge is connected with the electrochemical potential and a doping level as $ \sigma = e ( n ( \mu ) \u2212 p ( \mu ) \u2212 N D + N A ) $, where *N _{D}* and

*N*determine the

_{A}*n*and

*p*-type doping concentrations,

*n*and

*p*are 2D carrier densities for electrons and holes. For a single-layer graphene, one has

^{11}

where $ \epsilon D $ is the Dirac point. The density of states in graphene is defined as $ DOS ( \epsilon ) = \xb1 8 \pi ( \epsilon \u2212 \epsilon D ) / ( h v F ) 2 $. At zero temperature and in the absence of doping, one has $ \mu = sign ( \u2212 \sigma ) \u210f v F \pi | \sigma / e | $, where $ v F $ is the Fermi velocity.

Summing, the problem of finding the electrochemical potentials is reduced to a system of two transcendental equations that can be solved numerically. What is important, variation of *V _{g}* and

*V*allows one to independently handle both the electrochemical potentials and positions of charge-neutrality points.

_{b}## IV. APPLICATION EXAMPLES

For structures with reduced dimension such as thin films, graphene, nanowires, nanotubes, etc., the number of summations in Eq. (1) decreases due to dimensional quantization. In this case, the calculations are directly dependent on what dimension is inherent in the contacts and what configurations they form. Some possible variants of configurations are considered in Ref. 44.

### A. 1D tunneling

As an illustration of the TH method, let us consider a simple model which is useful in the analysis of tunneling from very small areas. In particular, this may correspond to a situation where tunneling is restricted to a small area of one electrode over to the same small area of the opposite electrode.^{45} Suppose that tunneling occurs only between an atom at one contact and the nearest atom at another one. In Eq. (6), the atoms are characterized by the coefficients *c _{R}* and

*c*, respectively. In this case, the sum in Eq. (6) reduces to a single term and one obtains

_{L}As is seen, the squared modulus of the matrix element is the product of probabilities for an electron to be on the selected atoms and the tunnel exponent. The latter is obvious and corresponds to the result Eq. (4) obtained within the WKB-approximation. Notice that in Eq. (4) the role of *c _{R}* and

*c*is played by the attempt frequencies.

_{L}Neglecting inelastic processes, and then substituting Eq. (27) into Eq. (1) allows one to express the tunnel current as

where

What is important, the relation (28) assumes that in-plane momentum is not conserved in the tunneling process. An indication of such momentum nonconservation was obtained in a number of experiments with GIG structures.^{3,41} One possible explanation was presented in Refs. 45 and 46 where it was shown that this behavior corresponds to a lateral coherence length for the tunneling of about 1 nm in the TH approach.

Notice that a similar to Eq. (28) expression is known in the theory of STM.^{31,32} In particular, it shows that measurements by means of the STM can give information on the distribution of the local DOS in the sample. The application of Bardeen's theory to different STM problems is discussed in detail in Ref. 47.

### B. Specificity of 2D-2D resonance tunneling

In what follows we will discuss the resonant nature of the tunneling in the 2D-2D tunnel contact having the form of a parallel-plate capacitor. Suppose that the current is directed normal to the plates (in the *x* direction). In this case, two of six sums in Eq. (1) disappear due to the effect of dimensional quantization in the *x* direction. After simple transformations, one gets

where $ k | | = ( k y , k z ) , \u2009 q | | = ( q y , q z ) $. This expression was used, for example, in Ref. 48 in analyzing the tunnel current between twisted graphene layers. For plates of infinite size, the wavefunction can be factored into the product of an *x*-dependent function and a 2D plane wave $ exp ( i k y y + i k z z ) $. Substitution of such a wave function into Eq. (6) gives two Kronecker's symbols $ \delta k y , q y $ and $ \delta k z , q z $ which remove two sums in Eq. (30). As a result, at zero temperature, we have

where the contacts are assumed to have the same band structure. When the energy spectrum is defined by a single-valued function of momentum, the resulting current will entirely be due to the overlap of the tails of the spectral functions whose maxima are spaced at *eV _{b}*. As a result, if the effects of inelastic scattering of electrons inside the contacts are neglected, tunneling turns out to be forbidden. However, in the case of a nonsingle valued spectrum, the current may flow even in the absence of inelastic scattering of electrons in the contacts. In particular, this situation takes place in graphene.

^{11}

### C. Armchair-edged graphene tunnel contacts

#### 1. TH method

Let us consider tunnel contacts on the basis of graphene layers with armchair edges (see Fig. 4). In the tight-binding approximation, eigenvectors were obtained in the form^{49}

where *n* and *m* number the vertical and horizontal rows, respectively (see Fig. 4), $ k x i = 2 \pi i / [ 3 a 0 ( N x + 1 ) ] , \u2009 k y j = 2 \pi j / ( 3 a 0 N y ) $, *a*_{0} is the distance between neighboring atoms, $ E ( k ) = \xb1 t 0 | f ( k ) | , \u2009 f ( k ) = \u2212 t 0 ( 2 e i k y a 0 / 2 \u2009 cos ( 3 a 0 / 2 ) k x + e \u2212 i k y a 0 ) \u2248 3 t 0 a 0 / 2 [ ( k x \u2212 K x ) + i ( k y \u2212 K y ) ] $. An approximate equality is valid near the point $ K = ( 4 \pi / ( 3 3 a 0 ) , 0 ) $. Notice that in the case of macroscopic contacts, both *N _{x}* and

*N*tend to infinity.

_{y}Taking into account Eq. (7) and the fact that transport proceeds near the *K*-point, one obtains for the first rows of contacts (*n* = 1)

The Kronecker delta allows one to remove summation over *q _{y}* in Eq. (1). Summation over

*k*and

_{x}*q*can be replaced by integration: $ \u2211 k x = ( \Delta k x ) \u2212 1 \u222b d k x $, where $ \Delta k x = k x i \u2212 k x i \u2212 1 $. The introduction of new variables

_{x}*E*provides additional factors $ ( \u2202 E S ( k x ) / \u2202 k x ) \u2212 1 $. The one-dimensional DOS reads

^{S}In the case when inelastic scattering of electrons in the contacts can be neglected, *δ*-functions make it possible to remove double integration in Eq. (1) with the substitution $ E L \u2192 \epsilon , \u2009 E R \u2192 \epsilon + e V b $. The result should be doubled since tunneling occurs between a pair of points *K* and $ K \u2032 $ in each contact. Finally, the tunnel current is written as

To determine the DOS let us take into account the explicit form of the energy dispersion near the *K*-point in graphene: $ E ( k ) = \xb1 ( 3 / 2 ) t 0 a 0 | k | $. One obtains

where $ v F = 3 t 0 a 0 / ( 2 \u210f ) $. Within the zero temperature approximation, integration over *k _{y}* gives

where $ z = \epsilon / ( \epsilon + e V b ) $, *K*(*z*) is the elliptic integral of the first kind, *V _{b}* is assumed to be low, and $ e V b > 0 $. At $ e V b < 0 $ the current will have the opposite sign. It should be noted that the resulting

*quadratic*I–V dependence is unique and is directly related to the semimetallic nature of graphene.

#### 2. NEGF method

Let us illustrate the Green's function technique by applying it to the above-discussed system because it makes it possible to compare the results. This method allows us to consider the configuration with two semi-infinite contacts of finite width shown in Fig. 5. As a first step, one has to calculate the edge Green's function $ g ( \epsilon ) $ of the contact (the SGF for this configuration). For this purpose, we choose the repeating block shown in Fig. 5 [cf. Fig. 2(c)] taking into account that its size significantly affects the calculation time. After that one has to formulate the TB Hamiltonian of the system. In our case, the chosen unit block contains four atoms, and therefore, the size of the Hamiltonian matrix is (4 × 4). The corresponding matrix elements are of the form $ H i j = \u3008 i | H | j \u3009 $, where $H$ is the electronic Hamiltonian and $ i , j = 1 , 2 , 3 , 4 $. For clarity, we will consider the simplest first-neighbor approximation. In this case, $ \u3008 i | H | i + 1 \u3009 = t 0 , \u2009 H i j = 0 $ for next neighbors, and $ H i i $ defines the Fermi energy which is taken to be zero. Finally, the TB Hamiltonian is the sum of matrices *H*_{00} (for the unit block), *H*_{01} describing the interaction between the blocks, and *T* characterizing the interaction between the contacts (see Fig. 5). Explicitly,

In accordance with Eq. (18), one can now calculate the SGF by using a simple iteration scheme

where the initial matrix *g*_{0} can be taken in an arbitrary form. The iteration procedure ends at step *N* if the difference $ g N \u2212 g N \u2212 1 $ becomes negligible.

A faster iterative method was proposed in Ref. 40. The algorithm is as follows:

Here, the initial values are taken to be $ \u03f5 0 = \u03f5 0 s = H 00 , \u2009 \alpha 0 = H 01 , \u2009 \beta 0 = H 01 \u2020 $, and the previous criterion of convergence can be used.

A slightly different approach can be used in the case of an infinitely long edge in the *y* direction. We choose a unit block in the *x* direction and divide it into periodic unit cells in the *y* direction (see Fig. 6). This makes it possible to use the Bloch's theorem, which leads to the appearance of a phase factor $ exp ( i k y y i ) $ for atomic orbitals $ \u3008 i | $ with *k _{y}* being the wavevector in the

*y*-direction and

*y*the coordinate of the

_{i}*i*th atom. The Hamiltonian matrices can be calculated as before, but now for the nearest neighbors $ \u3008 i | H | j \u3009 = t 0 \u2009 exp ( \u2212 i k y ( y i \u2212 y j ) ) $. One obtains

and *T* is the same as in Eq. (37).

Now one can start calculating the edge Green's function using the above-described iteration scheme. Let us take the parameter values as $ t 0 = 1 $ and $ t = 0.06 t 0 $, which approximately corresponds to the interatomic distance $ r 0 = 2 a 0 $. Notice that the values of *t*, characterizing the tunneling between atoms of different contacts, were calculated by

where *r*_{0} is the distance between the nearest atoms and *p _{z}* orbitals of carbon atoms were taken from the tables.

^{50}

The total DOS per atom on the edge of finite width is written as

where $ g n , n $ are the diagonal elements of the matrix *g* for the atoms on the edge. In the case of a repeating unit cell, one can use the modified variant of Eq. (42)

where *N _{c}* is the number of edge atoms in the unit block and $ a c = 3 a 0 $ is the translation vector (see Fig. 6). When performing numerical calculations, integration should be replaced by summation. A similar modification should be done in Eq. (15).

Notice that the tunnel current can also be calculated using Eq. (16). To this end, one has to consider a part of any edge as a dot lying between contacts, like that shown in Fig. 2(c). Further procedure is as follows: (1) we define the Hamiltonian of the dot *H _{D}* and the matrix

*t*describing the interaction of contacts with the dot, (2) calculate SGF for contacts as in the previous case, (3) calculate Σ

_{S}_{S}and Γ

_{S}and find Green's function for the dot $ G D + $ by Eq. (17), and (4) substitute the obtained Γ

_{S}and $ G D + $ in Eq. (16).

It is reasonable to expect that with increasing width of the tunneling contact shown in Fig. 5, both the DOS and the current will tend to their values for infinitely wide contacts. To illustrate this, we performed calculations for contacts with a width of 1, 5, 10, 50, and 200 times larger than the initial one (shown in Fig. 5) as well as for the infinitely wide contact (see Fig. 7 where the corresponding I–V curves are shown).

Notice that the behavior of I–V curves is directly connected to the existence of peaks in the DOS near the Fermi level. As is known, such peaks are present at zigzag-type edges, while there are no localized edge states in the armchair configuration. In our case, corner atoms belong to both zigzag and armchair edges. This provides an increased local DOS not only on them, but also on neighboring atoms of the armchair-type edge. For the very narrow structure shown in Fig. 5, the peaks are absent, and therefore, the current is small and monotonically increases with voltage. For wider structures, such peaks make the main contribution to the current (see, e.g., structures with a five- and tenfold increased width in Fig. 7). The further increase in the width leads to a decrease of the contribution from the corner atoms; however, additional peaks appear on the I–V curves due to van Hove singularities. It should be finally noted that the limiting curve (for the widest contacts) almost coincides with that for infinitely wide contacts and exactly follows the analytical result in Eq. (36).

### D. Vertical graphene TFET

The first experimental investigations of the vertical GIG tunneling structure were reported in Refs. 3, 41, and 46. The test samples exhibited room-temperature switching ratios of about 50 and 10,000 with atomically thin boron nitride or molybdenum disulfide acting as a transport barrier, respectively.^{3} Both types of devices did not show the effect of a negative differential resistance (NDR); however, it has been revealed in further experiments.^{2} It is important to note that the very appearance of NDR in these structures is directly related to the conservation of momentum.^{11,45}

The vertical GIG is a typical configuration where the 2D-2D tunneling occurs. Let us note several differences from the previously considered cases: (1) there is a gate electrode, so that varying the gate potential will change the chemical potential in the contacts, (2) graphene electrodes can be initially doped by electrons or holes, (3) we are dealing with the parallel-plate capacitor with a fairly large capacity, and therefore, the charging effects become significant, (4) the specificity of both the wave functions and the electronic spectrum of graphene as well as the fact that graphene layers are separated by insulators should be taken into account.

The three-plate capacitor model considered in Sec. III allows us to take into account the first three factors. The remaining one is convenient to consider within the Bardeen's approach. This requires knowledge of the explicit form of the electron wave functions in graphene contacts. As was shown in Refs. 2 and 11, they can be written as a product of 2D-Bloch functions and 1D-function that exponentially decays in the region of the potential barrier. For large electrode areas with no misorientation between the electrodes, the matrix element was found to be

where $ | g ( k | | , q | | ) | $ is a factor of order unity describing the overlap of periodic part of the lateral wavefunctions and *l* is approximately equal to an interplanar separation in graphite. This fairly simple expression was obtained for the GIG tunnel junction in Ref. 11 where it was noted that the factor $ | g ( k | | , q | | ) | $ has only a relatively small influence on the final results for the tunnel current. Comparison with Eq. (4) shows that $ \u210f k \u22a5 / m $ plays the role of the velocity. The type of homogeneous insulator specifies the height of the tunnel barrier and the effective mass of tunneling electrons. This approximation is valid at high misalignment between the graphene layers and the hexagonal boron nitride (hBN) barrier structures. In this case, any tunneling process involving scattering by hBN reciprocal lattice vectors is unable to scatter graphene electrons between the vicinity of the first Brillouin zone corners on the two layers.^{43} The method of accounting an arbitrary orientation of hBN is developed in Ref. 51. Provided that momentum is conserved, the effect of the band structure on the tunneling process can be taken into account by means of Eq. (31). Neglecting inelastic scattering, one has

Here, as before, $ E S ( k | | ) $ is the energy spectrum of carriers in contacts. Therefore, tunneling is allowed only for electrons with such momentum at which the argument of *δ*-function vanishes. With Eq. (21) taken into account, one obtains the following condition:

where the sign ± corresponds to the electron/hole part of the energy spectrum of contacts, respectively. There are two types of solutions to Eq. (46): (1) the signs are chosen to be opposite (tunneling occurs between the conduction and valence bands), which gives $ p | | = | F b d / ( 2 v F ) | $, and (2) the signs are the same, so that *F _{b}* = 0. In the last case, Eq. (46) is fulfilled for any $ p | | $, that is all charge carriers located in the energy window are involved in tunneling. This leads to the appearance of a resonance peak in the tunneling current. The corresponding relation between

*V*and

_{b}*V*can be found from Eq. (23) taking into account that $ \sigma 3 = 0 $. In particular, at zero temperature and without doping one obtains

_{g}When *V _{g}* is fixed and

*V*grows from zero, the current tends to a resonance region, where Eq. (47) is fulfilled and

_{b}*F*vanishes. Further increase in voltage leads to disruption of the resonance condition. As a result, the current decreases, thus revealing a region with the NDR.

_{b}It should be stressed that the direct use of Eq. (45) at resonance leads to the divergence of the current due to the *δ*-function term. Two approaches have been proposed to solve this problem. The simplest way is to assume the presence of inelastic processes.^{43,52} In this case, the *δ*-function can be replaced by the Lorentzian containing an energy broadening parameter giving the energy bands a finite width. Another idea is to admit the limitation of lateral momentum arising through some limited area *L* over which the coherent tunneling occurs.^{2,11,45} In this case, one has to return to Eq. (44) where the Kronecker's symbol can be replaced by the Lorentzian function with a characteristic width $ q c \u2248 2 \pi / L $.

### E. Planar graphene TFET

Another interesting example of the GIG structure is based on the flat configuration where the tunnel barrier is governed by a graphene nanogap. Specificity of current flow through graphene tunnel contacts lying in the same plane [see Fig. 1(b)] was studied in Refs. 17 and 18. Simulations were performed using the NEGF formalism. Like in the previous case, this device shows NDR with respect to the drain voltage while the current on/off ratio can reach a value of the order of 10^{3}. What is important, the behavior of I–V curves is directly related to energy distribution of the DOS in the graphene contacts. The most promising for application are contacts with sharp peaks near the Fermi level in the DOS.

To clarify the importance of this finding, let us return to the above-discussed graphene nanogap with armchair edges. As is known, graphene edges have a significant effect on the low-energy electron spectrum.^{1} For contacts fabricated of graphene or bilayer graphene with the armchair edge, the DOS near the Fermi level is identical to that in the bulk. The analysis shows that for bilayer graphene contacts the drain current obeys Ohm's law and the conductivity increases with the gate potential *V _{g}*. The situation for tunnel contacts of monolayer graphene is similar to the previous one with one significant difference: the current depends quadratically on voltage at small

*V*in accordance with Eq. (36).

_{b}The situation changes significantly in the case of tunnel contacts formed by two sheets of graphene with the zigzag edges. It has been found that the singular electronic states arise at the Fermi level, whose wave functions are mainly localized on the zigzag edge. This results in remarkably sharp peaks in the DOS near the Fermi level. Of prime importance is the fact that these states appear due to nontrivial topology of the electron band structure or more precisely, to the so-called Zack phase (a variety of Berry's phase) collected by the wave function in bypassing the Brillouin zone.^{53} Thus, the edge states are topologically protected in the sense that their existence is based on a topological property and remain robust to weak external perturbations. It was established that these states considerably affect the electron transfer in graphene ribbons with the zigzag edge.^{36} We managed to show that the unique properties of the edge states critically influence the tunnel processes too. As was mentioned in Ref. 18, three conditions must be met: (1) both graphene electrodes must be oriented toward each other with crystallographically matching zigzag edges, (2) in order to reach the regime of tunnel current, the gap between contacts must be narrow, and (3) the gate has influence on both electrodes.

In order to explain the working principle of the planar graphene TFET, let us consider schematic electronic DOS diagrams shown in Fig. 8. All states below the Fermi level are filled, and all above are empty. A dotted rectangle denotes a transport window, i.e., the energy interval in which a directed tunneling current is possible. At zero gate voltage, the picture is symmetric [Fig. 8(a)], and tunnel current is absent. For finite bias voltage, the bands are shifted and the transport window is open, which leads to the appearance of tunnel current. The maximum tunnel current takes place in the configuration shown in Fig. 8(b) when the filled states of the DOS peak in the left electrode are situated exactly opposite to the empty states of the DOS peak in the right electrode: the transistor is turned on. Further increase in bias voltage leads to a larger shift of the peaks [see Fig. 8(c)] so that the DOS is reduced and tunnel current comes down. This gives a region with NDR. Applied gate voltage will cause a shift of the Fermi level to the region of the low DOS in both contacts, and, as a result, tunnel current turns out to be strongly suppressed [see Figs. 8(d) and 8(e)] until one of the DOS peaks meets the energy window [Fig. 8(f)].

As a final remark, the suggested TFET has a rather simple configuration consisting of two facing each other in-plane electrodes with a nanogap between them. While tunneling has not been measured so far, current progress in STM lithography shows the feasibility for the fabrication of such device in the near future.^{20,54}

## V. CONCLUSION

It seems useful to give an overview of the papers where one of the above approaches was effectively used to describe tunnel current in the vertical GIG structure. As far as we know, the first application of the TH method to this problem has been reported in Ref. 11. The suggested theory was successfully applied in the calculation of I–V characteristics of the SymFET at room temperature.^{12} A similar description was used in the interpretation of the experimental results in Refs. 2, 3, 55, and 56 (see a detailed comparison of both approaches in Ref. 45). Calculation of the tunneling current in a G/hBN/G heterostructure as a function of the twisting between the crystals revealed a regime of coherent tunneling and NDR.^{51} Here, the TH method was adapted to analyze the role of the interfacial insulator layer (hBN) at small twisting angles. Twist-controlled resonant tunneling between monolayer and bilayer graphene was investigated in Ref. 43. A theory of the tunneling characteristics between misaligned graphene layers which takes into account the spectral properties of the tunneling electrons has been presented in Ref. 57. It was shown that, in the regime of resonant tunneling, vertical transport in G/hBN/G heterostructures carries precious information on electron–electron interactions and the quasiparticle spectral function of the two-dimensional electron system in graphene. Let us also mention the experimental papers^{42,52} on twist-controlled resonant tunneling in G/hBN/G heterostructures where different variants of the TH theory were applied. One can state that this method proved to be very productive and allowed to predict and explain many interesting phenomena.

There are only a few examples of applying the NEGF formalism to study the vertical GIG structure. An interesting article is^{58} where a self-consistent capacitance model to obtain the electrochemical potential profile across the device was used, and then the current and conductance across the device were computed within the NEGF method. The problem of NDR in G/hBN/G heterostructures with different number of hBN layers was investigated in Ref. 59. The effect of the rotational alignment between hBN slab and the graphene layers, as well as inelastic phonon assisted tunneling in the vertical current of a G/hBN/G device, was studied in Ref. 60.

A combination of both methods was used in the description of G/rotated hBN/G heterostructure.^{61} It was found that the NEGF formalism is effective in the case of large misorientation angles of the hBN with respect to the graphene layers. As the rotation angles become smaller, the commensurate unit cells become very large. As a result, NEGF calculations with the large tight-binding Hamiltonians become computationally challenging. In order to better understand the physics governing the interlayer transport at small rotation angles, an effective continuum model was suggested.

As a brief conclusion, we have discussed two approaches which seem to be most suitable in analyzing tunnel processes in different systems. As it was noted, the Bardeen tunnel Hamiltonian method provides a detailed analysis of the problem of weakly interacting contacts. The obvious advantage of the method is a possibility to obtain analytical solution in some cases, which gives physical insight into the processes and the role of different factors. The method of nonequilibrium Green's functions is more universal and well describes the cases of noticeable mutual influence of contacts. Another advantage of the method is its being adapted to numerical calculation in modeling an atomistic structure of electrodes or layers in GIG devices, yet it is a drawback since calculations of this type require large computing power with increasing sizes. One illustrative example is given in Sec. IV C. The use of TH method enabled us to obtain explicitly the quadratic I–V dependence and explain it by the linear energy dispersion in graphene contacts at low voltage. In turn, the application of the NEGF method gave accurate quantitative values of tunnel current in a wide range of voltages. At the same time, it should be clearly understood that both methods are tightly connected and only use a different calculation technique based on the same approximations.

We would like to emphasize once more that creation of devices based on the use of the tunnel effect is one of the most promising directions of modern nanoelectronics. It is anticipated that in the coming few years, the progress of silicon technology will face some fundamental difficulties related to critical miniaturization, namely, at characteristic sizes less than 10 nm the role of quantum tunnel effects becomes crucial. This makes it difficult to create electronic devices working on old principles. Owing to the synthesis of new materials with unique physical properties, it is possible not only to overcome this difficulty but also to use these materials for the production of nanoscale diodes, transistors, and other devices.