Abrupt current switching in graphene bilayer tunnel transistors enabled by van Hove singularities

In a continuous search for the energy-efficient electronic switches, a great attention is focused on tunnel field-effect transistors (TFETs) demonstrating an abrupt dependence of the source-drain current on the gate voltage. Among all TFETs, those based on one-dimensional (1D) semiconductors exhibit the steepest current switching due to the singular density of states near the band edges, though the current in 1D structures is pretty low. In this paper, we propose a TFET based on 2D graphene bilayer which demonstrates a record steep subthreshold slope enabled by van Hove singularities in the density of states near the edges of conduction and valence bands. Our simulations show the accessibility of 3.5 × 104 ON/OFF current ratio with 150 mV gate voltage swing, and a maximum subthreshold slope of (20 μV/dec)−1 just above the threshold. The high ON-state current of 0.8 mA/μm is enabled by a narrow (~0.3 eV) extrinsic band gap, while the smallness of the leakage current is due to an all-electrical doping of the source and drain contacts which suppresses the band tailing and trap-assisted tunneling.


Self-consistent evaluation of band gap and carrier density in gated bilayer
The intrinsic graphene bilayer (GBL) is a gapless semiconductor with symmetric electron-hole dispersion being parabolic near the band edges. However, the application of transverse electric field F ⊥ leads to the potential energy difference ∆ = eF ⊥ d between the layers comprising GBL (d = 3.35 A is the interlayer separation in GBL) and opening of the band gap. As the transverse electric field depends on the carrier density, which, in turn, depends on the band structure, a self-consistent procedure for the determination of those quantities in gated GBL is required.
In this section, we describe the procedure used to evaluate the carrier density, bandgap, and electric potential in GBL at given potentials of top and bottom gates (V t and V b ). The tight-binding Hamiltonian of graphene bilayer yields the following energy bands depicted in Fig. 1: 1 with the corresponding four-component eigen wave functions (non-normalized): p miñ (1) (2) (3) (4) is the electron momentum at the band edges The components of wave function determine the amplitudes of finding an electron on the atoms of type "A" and "B" of the unit cell on the top and bottom layers in the i-th band. In Eqs. (1) -(4), p is the electron quasimomentum, γ 1 ≈ 0.4 eV is the interlayer hopping integral, and v ≈ 10 6 m/s is the Fermi velocity in graphene. Given the eigenfunctions, the probabilities of finding an electron of the i-th energy band on the top and the bottom graphene layer can be calculated as follows: here ||...|| stands for the norm of the vector. The net charge densities at top and bottom graphene layers are obtained by summing the probabilities (5) timed by the distribution functions over the bands and momenta: Here, g = 4 is the spin-valley degeneracy factor, f (E, µ) = [e E−µ kT + 1] −1 is the Fermi-Dirac distribution function, µ is the chemical potential reckoned from the midgap, k is the Boltzmann constant, and T is the absolute temperature.
We further consider graphene bilayer placed between the top and bottom gates separated by gate dielectrics of thicknesses d t and d b and dielectric constants ε t and ε b , respectively (see Fig. 2).
With the aid of Gauss law, the potential energy difference between layers, ∆ = e(ϕ b − ϕ t ), is readily expressed via the gate potentials V t and V b : where ε 0 is the vacuum permittivity. Assuming that far from gates graphene is undoped (µ = 0), and the source contact is Near the grounded source contact there is no bandgap and GBL has zero potential, while under the "doping" gate GBL exhibits a bandgap∆ and has potential (ϕ t + ϕ b )/2, which leads to a corresponding midgap shift grounded, we obtain an expression for the chemical potential in the source region: Equations (2)-(10) form a system of equations used for a self-consistent evaluation of ∆ and µ in the source region. The system, which can be written in the symbolic form was solved numerically using an iterative method similar to the nonlinear successive over relaxation: 2 where i is the iteration number, and λ is the relaxation parameter adjusted to achieve convergence. The same approach was applied for the evaluation of ∆ and µ in the drain region (with an obvious shift of electric potentials by the drain voltage). Considering the GBL band structure under the control (middle) gate, we assumed that the conduction band electrons are in equilibrium with the drain contact, while the holes are in equilibrium with the source contact. This approach is justified as far as the tunneling probability is small and the electrons injected from the valence band do not distort the potential distribution. To this end, two different chemical potentials for electrons and holes were used in Eqs. (7), (8) when calculating the charge density under the control gate.

Distribution of electric potential, barrier transparency, and current
The calculation of the tunnel current requires the knowledge of the barrier transparency which, in turn, depends on the distribution of electric potential in the tunneling region. Various studies of the tunnel diodes confirm that it is the maximum field in the junction region that governs the barrier transparency. [3][4][5] Thus, for an estimate of the tunnel current, the full knowledge of potential distribution is not necessary.
To obtain the junction field, we have solved the Laplace equation in the region between two gate electrodes: the "doping gate" above the source (held at potential U S ) and the top control gate above the channel (held at potential V G ) treating them as two keen metal plates. 6 The lateral distance between the gates is denoted by d g . The effect of the bottom gate on the junction field is weak provided d b d t . The difference between dielectric constants of top and bottom dielectric layers, ε t = ε b , can be approximately taken into account by replacing the top gate oxide thickness with effective oxide thickness, d e f f = d t ε b /ε t . The distribution of electric potential in the geometry considered is obtained with conformal mapping technique, and the result is Figure 3. Distribution of the electric potential produced by two keen planes held at potentials U S and V G . Potential is encoded by colour, increasing from blue to red, while arrows represent field lines (see Fig. 3), where the origin of the coordinate grid is located between the two gates and lies in their plane.
The maximum electric field in the graphene plane (y = −d e f f ) is achieved at x = 0 and equals: Having obtained the junction field F, we calculate the WKB transparency of the potential barrier D (p ⊥ , E) separating the conduction and valence bands. Assuming that the potential energy depends on x -coordinate linearly, we find where p and p ⊥ are the components of electron quasimomentum that are parallel and perpendicular to the direction of junction field, respectively, E is the electron energy measured from the midgap in the source region, the integrals are taken between the classical turning points x 1 and x 2 . Generally, the interlayer asymmetry parameters in the source (∆ S ) and in the channel (∆ C ) are different. For simplicity, we perform a linear interpolation between ∆ S and ∆ C in the junction region, where ϕ S and ϕ C are the electric potentials in the source and channel regions (obtained in the previous section), and l is the length of the transition region between them. The tunneling current density J tun can be written as where   3) and (4)] which contribute to the current independently the transverse electron momentum p ⊥ and its total energy E which are conserved during the tunneling: , µ S and µ D are the chemical potentials in source and drain regions, respectively (we assume the electrons in the conduction band of the channel region are in equilibrium with the drain contact), p max (E) is the maximum transverse momentum allowed for given total electron energy E in the valence band of source and in the conduction band of the channel. More precisely, p max (E) = min {p c (E), p v (E)}, where p c (E) and p v (E) are the inverse functions of the electron dispersion law E(p) in the conduction band of the channel and the valence band of source, respectively. 7 Finally, we note the appearance of the factor of two before the whole expression for the tunnel current which is unique to the band structure of graphene bilayer. It appears due to the presence of two points with zero group velocity in the dispersion laws for electrons and holes, which implies that a carrier incident on a barrier attempts to pass through it twice (see the red arrow in Fig. 4).
Thermionic currents J e and J h carried by electrons and holes, respectively, are calculated in a similar way, except for the absence of the factor of two and barrier transparency equal to unity. The integration limits E in the expressions for the thermionic current are restricted to the energies above the barriers created by the doping gates: where E S c and E D v are the edges of conduction band in source region and valence band in drain region, respectively. Electrons from the "hat" (i.e. with p < p min ) do not contribute to the thermionic current because they can neither remain in the "hat" in all three regions (source, gate and drain) nor get out of it (it is possible only at band edges).

Derivation of approximate analytic expressions for tunnel and thermionic currents
Firstly, we need to obtain an analytic expression for barrier transparency. The integrand in Eq. (15) is too complicated to carry out the integration exactly. However, using the fact that the integrand has an arch-like shape (it turns to zero at x 1 and x 2 and reaches maximum somewhere between), we can approximate it with a half ellipse of the same width and height (see Fig. 5 A).
In the simplest case of constant ∆ = ∆ tun the maximum of the integrand is reached in the middle of the barrier and equals max Im p (x) = Im p 2 (0, ∆ tun ) − p 2 ⊥ = Im The main contribution to the tunnel current is from electrons with p ⊥ < p min , because for them barrier is the lowest and the shortest. Such electrons tunnel from the edge of the valence band to the edge of the conduction band, which gives the following barrier length: Equations (15), (21) and (22) yield the following expression for barrier transparency in the "semielliptical" approximation: where To avoid special functions, we extend the integration region up to infinite p ⊥ and obtain Now we move on to the calculation of the thermionic current. Near the band edges the dispersion law is approximately parabolic: where plus sign corresponds to the conduction band and minus sign corresponds to the valence band. Taking into account that in Eq. (19) the Fermi-Dirac distribution almost coincides with the Boltzmann distibution, we obtain where index "S" corresponds to source, and "D" corresponds to drain.

Comparison with TFETs based on 2d semiconductors with parabolic bands
In this section, we compare the dependencies of interband tunnel current on the band overlap in GBL and an equivalent 2D semiconductor without the 'Mexican hat' band structure. The equivalence of the materials implies equal values of band gap and barrier transparency in the electric field, while the band structure outside the gap is different. For the sake of analytical traceability, we assume the barrier transparency D 0 to be weakly dependent on electron energy and gate voltage. The low-energy band structure of GBL is given by Eqs. (2). The electron dispersion of a parabolic-band material is in the conduction band, in the valence band. Here m c and m v are the effective masses for in the respective bands. The interband tunneling occurs mainly between conduction and light hole bands as the barrier transparency is small for heavy carriers.
Similarly to the derivation of Eq. (2) of the main text, one can write down the current density in a parabolic-band 2d semiconductor

8/11
where p max (E) is the maximum transverse momentum of electron with energy E in the conduction and valence bands. The value of p max (E) is the minimum of electron momenta at given energy in the conduction and valence bands: Graphically, the current J t,par is proportional to the area under the inverse dispersion laws in the conduction and valence bands, this region is filled with blue in Fig. 7. The integral in (33) is evaluated analytically with the result In graphene bilayer, the tunneling current density is given by where an extra factor of g v = 2 accounts for the valley degeneracy, and the factor of two in front of the barrier transparency is due to the two points in electron dispersion E c (p) with zero group velocity where an incident electron attempts to tunnel from the valence band. Again, p max (E) = min {p c (E) , p v (E)}, where p c,v (E) are the dependencies of electron momentum on energy in the conduction and valence bands. Inverting the electron dispersion, one finds The energy integral (36) for GBL cannot be evaluated analytically. However, the area under the dependence p max (E) is naturally larger for graphene bilayer than for a typical narrow-gap semiconductor, not to mention an extra factor of four due to valley degeneracy and twofold tunneling. At small energies close to E c , p c (E) ≈ p min > 2m c (E − E c ), while at large energies p c (E) ≈ E/v 0 , which again lies above 2m c (E − E c ). In Fig. 5 of the main text, we numerically compare the current densities for GBL and 2d semiconductor with parabolic bands using Eqs. (36) and (35). At the operating voltage of 150 mV, the current in GBL is 15 times larger than the current in a parabolic band semiconductor, of which the factor of 2 is due to the extra valley degeneracy, the factor of 2 is due to the two-fold tunneling, and the factor of 3.5 is due to the 'Mexican hat' dispersion.

Estimate of the band tails in graphene bilayer
The abruptness of tunnel switching is limited by the tunneling from the localized states in the band tails. 8 This mechanism can be alternatively interpreted as tunneling from the band tails, produced by the random dopant-induced potential fluctuations. In this section, we evaluate the density of states in graphene bilayer due to these fluctuations. In the quasi-classical limit, the 'true' density of states ρ(E) is evaluated by integrating the bare density of states shifted by eV , ρ 0 (E − eV ), timed by the probability of voltage fluctuation of magnitude eV In the above equation, V 2 is the root-mean square voltage fluctuation, where V 0 (q) is the Fourier transform of the potential produced by a single impurity and n i is the impurity density. To obtain a non-divergent result for fluctuations in two dimension, a finite distance from the charged impurities to the 2d plane has to be considered, 9 in this case the Fourier transform of potential is Hole leakage

Dielectric
Dielectric Metal Graphene d t Figure 8. A model band diagram of graphene bilayer -dielectric -metal system used for the evaluation of gate leakage current. Graphene bilayer is modelled as a delta-well with potential U(z) = −αδ (z + d t ), where α is chosen to provide a correct work function from graphene to the insulating material.
Above, κ is the background dielectric constant and ε(q) is the dielectric function of graphene bilayer itself. For analytical traceability, the latter is taken in the Thomas-Fermi screening approximation: 10 where the screening wave vector is calculated as The integral (40) is evaluated analytically with the result where K (x) = −e 2x (2x + 1)Ei(−2x) − 1. Equations (39) and (44) were used to calculate the density of states in the Fig. 5 of the main text numerically. The evaluation of DoS is, actually, an iterative procedure as the screening wave vector depends on the DoS itself. However, it converges quickly and the first approximation already is close to the final result.

Modeling of the gate leakage
To estimate the leakage current from graphene bilayer to the metal gate, one has to multiply the charge density in graphene by the area of the gate and by the characteristic 'tunneling frequency' ν t . The latter can be derived from a lifetime of electron state bound to graphene bilayer in the presence of continuum of states representing the metal gate (Fig. 8).
The highest voltages are applied to the 'doping' gates, and the largest leakage is expected therein. An example of the band diagram corresponding to source doping gate is shown in Fig. 8. Graphene bilayer is modelled as a delta-well U(z) = −αδ (z + d t ). The parameter α = 2h U b /2m * is chosen to provide a correct value of the work function from graphene bilayer to the surrounding insulator material, m * is the electron effective mass in the insulator.
A simple wave function matching procedure for the potential profile of Fig. 8 leads to the dispersion equation for the quasi-bound states in graphene bilayer where κ B = √ 2mU b /h is the decay constant of the wave function, k F is the electron wave vector in metal (approximately equal to the Fermi wave vector), and κ = √ −2mE/h is related to the sought-for state energy E < 0. The solutions of Eq. (45) are

10/11
complex, E ≈ −U b + ihν t /2, where For non-rectangular barriers, the exponent e −2κ B d should be replaced with the quasi-classical transparency of the barrier separating graphene bilayer and the gate. The full leakage current density (per unit width of the channel) is evaluated as J g = eL g ν t (n e + n h ), where n h is the density of holes induced at the source doping gate (4.5 × 10 12 cm −2 for U S = −0.6 V) and n e is the electron density at the drain doping gate (8 × 10 12 cm −2 for U D = 0.25 V). To provide minimum tunneling frequency, we choose zirconium oxide as a gate dielectric due to its relatively large effective (tunneling) mass m * = 0.3m 0 and small electron affinity χ ZrO 2 = 1.6 eV leading to a quite large value of U b = χ Gr − χ ZrO 2 = 2.9 eV.