Drift-diffusion models for the simulation of a graphene field effect transistor

A field effect transistor having the active area made of monolayer graphene is simulated by a drift-diffusion model coupled with the Poisson equation. The adopted geometry, already proposed in (Nastasi and Romano in IEEE Trans. Electron. Devices 68:4729–4734, 2021, https://doi.org/10.1109/TED.2021.3096492), gives a good current-ON/current-OFF ratio as it is evident in the simulations. In this paper, we compare the numerical simulations of the standard (non-degenerate) drift-diffusion model, that includes the Einstein diffusion coefficient, with the degenerate case.


Introduction
Since its discovery, graphene has attracted scientists for its possible use in electronics due principally to the high electron mobility; moreover, the thickness of just one atom [2] could represent the ultimate miniaturization. In the last years, several graphene field effect transistors (GFETs) have been designed [3] with the aim to substitute the standard silicon metal oxide field effect transistors (Si MOSFETs). The main issue is due to the absence of gap in monolayer graphene, causing the ambipolar conduction in the device and, consequently, a very narrow current-off region. In addition, the physical features of the metal-graphene interface at the contacts are not completely understood because they are rather different with respect to the case of traditional semiconductors in which the Ohmic or Schottky conditions are usually appropriate. The considerations above make to conceive GFETs a very challenging target in nano-electronics.
Appropriate and accurate simulations of the current flowing in the channel of GFETs are crucial to help the designers in improving and optimizing such devices. The drift-diffusion models represent the simplest way of modeling and are employed in all the standard simulation tools. In the case of GFETs several applications are already present in the literature [4][5][6][7] but only compact models are employed for the electrostatic potential. Instead, in [8] a full 2D Poisson equation is included and coupled with the drift-diffusion equations. A new GFET geometry that seems to overcome the zero gap issue and presents a wide off region, robust enough for engineering purposes, has been introduced in [1]. An impor-

Mathematical model
The physical setting is represented by the device proposed in [1], here depicted in Fig. 1 In order to simulate the current flowing in the channel of the device, we adopt the 1D bipolar stationary drift-diffusion model e being the (positive) elementary charge. J n (x) and J p (x) are the electron and hole current respectively. R denotes the generation-recombination term. The previous equations are coupled with the 2D Poisson equation for the electrostatic potential in the whole section.
The system (1) is set in the interval [x 1 , x 4 ] and is augmented with Dirichlet boundary conditions for electron and hole densities as will be explained below.
In this paper, we adopt two different models for the expression of currents, one in the degenerate case and another one for the non-degenerate case.

Degenerate case
The electron and hole current densities are expressed in terms of the quasi-Fermi energies (see for example [7]) where n and p, ε (n) F and ε (p) F , μ n and μ p are the densities, quasi-Fermi energies and mobilities of electrons and holes, respectively.
The electron and hole densities, n(x) and p(x) respectively, are evaluated as with g s = 2 and g v = 2 the spin and valley degeneracy. f FD indicates the Fermi-Dirac distribution where ε D = -eφ(x, y gr ) is the Dirac energy and φ(x, y) is the electrical potential, here evaluated on y = y gr , y gr being the average y-coordinate of the graphene sheet (see Fig. 1). ε F denotes the Fermi level (in pristine graphene ε F = ε D ), ε(k)ε D = v F |k| represents the graphene dispersion relation (strictly valid around the Dirac points), which is the same for electrons and holes (see [2,9,10]), is the reduced Planck constant, v F the Fermi velocity, k B the Boltzmann constant, T the lattice temperature, kept at 300 K (room temperature). The crystal momentum of electrons and holes is assumed to vary over R 2 .
We impose Dirichlet boundary conditions on the system (1), (2) The quantity ε F is the difference of the work functions between metal and graphene. Its value depends on the material the contact is made of. However, only part of the electrons, which the metal can furnish, enters graphene because of quantum reasons (for further explanations see [1]). This gives rise to a kind of contact resistance which is modeled by multiplying ε F by the parameter 0 < α ≤ 1. In this paper, we use the value ε F = 0.25 eV and α = 1. Simulations with different values of α can be found in [1].

Non-degenerate case
In the presence of low density the Fermi-Dirac distribution can be approximated by the Maxwell-Boltzmann one and the current relations can be written as (see for example [11]) where D n = μ n k B T and D p = μ p k B T are the classical Einstein relations (strictly valid only at equilibrium) for the electron and hole diffusion coefficients, and φ represents the electrostatic potential along the graphene sheet.
In this case, we impose the following Dirichlet boundary conditions on the system (1), (4)

Generation-recombination and mobilities
In both cases, for the generation-recombination term R we adopt the Shockley-Read-Hall model, as suggested in [4] by analogy with standard semiconductors where n e and p e are the equilibrium electron and hole densities, that is obtained with zero source-drain voltage, and τ is the generation-recombination relaxation time. For the simulations, in the present paper we have set τ = 10 ps. The qualitative behavior is similar for other values of τ [1]. Indeed, in [4] an additional term has been included for describing tunnelling effects. It depends on several experimental parameters and its effect should be not very relevant [1]. Here, it is neglected. The mobility models μ n = μ n (n, E) and μ p = μ n (p, E) depend in general on the carrier densities n and p and on the electric field E, defined along the graphene sheet as They represent a crucial point for an accurate determination of currents. A mobility model obtained from experimental results has been introduced in [12]. In this paper, we adopt a model deduced by a fitting procedure on extensive simulation of the homogeneous Boltzmann equation solved by a deterministic numerical method [13,14]. It has been presented in [15,16] and in the general form, including also the presence of a substrate as in the situation adopted here, in [8].

Poisson's equation
The system (1) is coupled with the 2D Poisson equation for the electrostatic potential, solved in the section of the device where The dielectric function is given by ox , otherwise, with gr and ox dielectric constants in graphene and oxide respectively. n imp is the areal density of the impurity charges at the graphene/oxide interface. Since n and p are areal densities, the charge in the graphene layer is considered as distributed in the volume enclosed by the parallelepiped of base the area of the graphene sheet and height t gr . Equation (6) is augmented with the following boundary conditions ∇ ν φ = 0 at the remaining part of the boundary.
Here V b is the bias voltage, V G u is the upper gate-source potential, V G d is the down gatesource potential. V G u and V G d include the (subtracted) flat-band voltages. We have denoted by ∇ ν the external normal derivative.

Numerical method
In order to solve the system (1)-(2) we adopt a finite difference scheme. A uniform mesh of N x grid points and size x is introduced in [x 1 , x 3 ] and for y = y gr . We denote by u i the numerical approximation of a generic function u(x) at the node i. First we discretize the equations (1) at the midpoint Then we discretize the equations for currents (2) where μ i+ 1 2 n , μ i+ 1 2 p , n i+ 1 2 and p i+ 1 2 are obtained by interpolating the same quantities at i and i + 1 nodes in the midpoint. In this way, an implicit scheme is devised.
We remark that the quantities μ  (3). For such a reason, the numerical solution is obtained by the following iterative procedure.
We denote by u (r) the value of the function u at the step r.
where n i is the intrinsic graphene concentration, and get ϕ (1) solving the Poisson equation (6)-(7). 2. For each r > 1: (a) find (ε (n) F ) (r) and (ε (p) F ) (r) solving the drift-diffusion equations (1)-(2) with the electrostatic potential ϕ (r-1) ; (b) calculate n (r) and p (r) through (3); (c) solve the Poisson equation (6)-(7) with densities n (r) and p (r) getting ϕ (r) . 3. Iterate step 2 until convergence. In order to improve the convergence, Pulay's method has been employed as in [6]. The difference between two iterations is measured as the relative variation of the current where J (r) n and J (r) p are the electron and hole currents at the rth iteration. As stopping criterion e n and e p less than 10 -4 has been required. Here · denotes the supremum norm.
The system (1), (4) is solved by the Scharfetter-Gummel scheme in the stationary case with the aid of Pulay's method, even in this case. The discretization is similar to the time dependent one and the interested reader is referred to [8].
To solve the Poisson equation numerically a uniform two-dimensional mesh is adopted. It is discretized by finite differences at the internal points with the standard five points stencil [8].

Simulation results
In this section, we show the results of numerical simulations of the device of Fig. 1. We have set the length 100 nm, the width of the lower and upper oxide (SiO 2 ) 10 nm. The graphene thickness t gr , including also the double graphene/oxide interface, is set to 1 nm. The gate contacts are long 50 nm. The lateral (source and drain) contacts are long 21 nm. The two gate potentials are both set equal to the same value V G in the simulated cases.
We considered a mesh of 41 grid points along the x-direction and 23 grid points along the y-direction. In the graphene layer a single row of 39 nodes has been employed.
In Fig. 2 we show, for several drain-source voltages V b and gate voltages V G , the characteristic curves of the simulated GFET obtained with the model (1), (2). In Fig. 3 the same curves are reported adopting the model (1), (4). There are evident qualitative and quantitative differences, in particular at the lower bias voltages considered in the simulations. For example in the case V b = 0.1 V, when V G = 1.5 we have a difference of about 74%. The discrepancy reduces at higher voltage but it is clear that the model (1), (4) overestimates the current anyway.
A comparison of electron densities with both models is shown in Fig. 4. The simulations are obtained with V b = 0.2 V and V G = -1.5 and 1.5 V. In Fig. 5 a similar plot is presented for hole densities. The differences are less evident with respect to the currents. This seems to indicate that the crucial point is a reduction of the mean velocity when the degeneracy effects are taken into account.
For the sake of completeness also the electrostatic potential is depicted (Fig. 6) in the same cases. For the electrostatic potential no significant differences are obtained. On the left the adopted model is (1), (2); on the right it is (1), (4) In order to better understand the origin of the discrepancy between the degenerate and non-degenerate cases, we look closer to the models. By observing that dk with ϕ azimuthal angle of k, the electron and hole densities can be written in terms of the Fermi integrals of order k with (x) Euler gamma function, as follows Here η (n) (x) = ε (n) F (x)ε D (x) and η (p) (x) = ε (p) F (x)ε D (x) are the chemical energy for electrons and holes respectively. Therefore J n = μ n n∇ε (n) F = μ n n∇ η (n) + ε D = μ n n∇η (n)eμ n n∇φ.
The second term is the standard drift one. By observing that F 1 (z) is a strict monotone function and recalling the well known properties of the Fermi integral dF k (z) dz = F k-1 (z), k = 1, 2, . . . and F 0 (z) = ln(1 + exp z), the first term is given by μ n n∇η (n) = μ n n d η (n) dn ∇n = μ n n dn dη (n) ∇n from which the electron diffusion coefficient in the degenerate case Similarly the diffusion coefficient of holes in the degenerate case is related to that in the non-degenerate case by .
In Fig. 7 the ratios r n = D d n D n and r p = D d p D p are plotted along the graphene in the simulated cased considered above. For electrons r n is considerably different from one near the contacts while r p is in a significant way greater than one in the channel. Anyway there is a clear evidence that the classical Einstein relations for the diffusion coefficients are not adequate in the simulation of the considered device. This casts serious doubts on the use of the non-degenerate drift-diffusion models in GFET.

Conclusions
The degenerate and non-degenerate drift-diffusion models have been compared in the case of charge transport in graphene by simulating the GFET proposed in [1]. Major differences are observed regarding the characteristic curves making the drift-diffusion models based on the Maxwell-Boltzmann statistics not adequate in the considered device. A more in-depth analysis shows that the root of the failure is in the diffusion coefficient which is very different in the degenerate and non-degenerate case. The Einstein relation is no longer valid in GFET.