Nastasi and Romano Journal of Mathematics in Industry (2022) 12:4

https://doi.org/10.1186/s13362-022-00120-3

RESEARCH Open Access

Drift-diffusion models for the simulation of a

graphene field effect transistor
Giovanni Nastasi1* and Vittorio Romano1
*
Correspondence:
giovanni.nastasi@unict.it Abstract
1
Department of Mathematics and
Computer Science, University of A field effect transistor having the active area made of monolayer graphene is
Catania, Viale A. Doria, 6, 95125 simulated by a drift-diffusion model coupled with the Poisson equation. The adopted
Catania, Italy 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.
Keywords: Graphene FET; Drift-diffusion model; Graphene

1 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 simu-
lation tools. In the case of GFETs several applications are already present in the literature
[4–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-

© The Author(s) 2022. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use,
sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original
author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other
third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a
copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Nastasi and Romano Journal of Mathematics in Industry (2022) 12:4 Page 2 of 11

tant issue is to address the importance of the degenerate statistics since often the standard
drift-diffusion model based on the Maxwell–Boltzman approximation of the Fermi distri-
bution is used. For such a reason, in this paper we compare the results obtained with the
drift-diffusion model in the degenerate case and in the non-degenerate case for the GFET
proposed in [1].
The plan of the paper is as follows. In Sect. 2 the mathematical model is described; in
Sect. 3 the numerical method is explained; in Sect. 4 numerical results and comparisons
are shown.

2 Mathematical model
The physical setting is represented by the device proposed in [1], here depicted in Fig. 1. It
is made of a sheet of graphene sandwiched between two layers of oxide. Source and drain
metallic contacts are placed along all the device thickness while gate contacts are located
in the central area of top and bottom oxide.
In order to simulate the current flowing in the channel of the device, we adopt the 1D
bipolar stationary drift-diffusion model

∂Jn ∂Jp
= eR, = –eR, (1)
∂x ∂x

e being the (positive) elementary charge. Jn (x) and Jp (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 [x1 , x4 ] 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.

Figure 1 Schematic representation of the simulated GFET. In the direction orthogonal to the xy-plane the
device is assumed infinitely long
Nastasi and Romano Journal of Mathematics in Industry (2022) 12:4 Page 3 of 11

2.1 Degenerate case

The electron and hole current densities are expressed in terms of the quasi-Fermi energies
(see for example [7])

(p)
∂εF(n) ∂εF
Jn = μn n , Jp = μp p , (2)
∂x ∂x

(p)
where n and p, εF(n) and εF , μn and μp are the densities, quasi-Fermi energies and mobil-
ities of electrons and holes, respectively.
The electron and hole densities, n(x) and p(x) respectively, are evaluated as

gs gv  
n(x) = fFD k, εF(n) dk, x ∈ [x1 , x4 ],
(2π)2 R2
 (3)
gs gv  (p) 
p(x) = fFD k, –εF dk, x ∈ [x1 , x4 ],
(2π)2 R2

with gs = 2 and gv = 2 the spin and valley degeneracy. fFD indicates the Fermi–Dirac distri-
bution
  
(ε(k) – εD ) – (εF – εD ) –1
fFD (k, εF ) = 1 + exp ,
kB T

where εD = –eφ(x, ygr ) is the Dirac energy and φ(x, y) is the electrical potential, here eval-
uated on y = ygr , ygr being the average y-coordinate of the graphene sheet (see Fig. 1). εF
denotes the Fermi level (in pristine graphene εF = εD ), ε(k) – εD = vF |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, vF the Fermi velocity,
kB 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 R2 .
We impose Dirichlet boundary conditions on the system (1), (2)

(p)
εF(n) (x1 ) – εD (x1 ) = εF (x1 ) – εD (x1 ) = αεF ,
(p)
εF(n) (x4 ) – εD (x4 ) = εF (x4 ) – εD (x4 ) = αεF .

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].

2.2 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])

∂n ∂φ ∂p ∂φ
Jn = Dn – enμn , Jp = –Dp – epμp , (4)
∂x ∂x ∂x ∂x
Nastasi and Romano Journal of Mathematics in Industry (2022) 12:4 Page 4 of 11

where Dn = μn kB T and Dp = μp kB T are the classical Einstein relations (strictly valid only
at equilibrium) for the electron and hole diffusion coefficients, and φ represents the elec-
trostatic potential along the graphene sheet.
In this case, we impose the following Dirichlet boundary conditions on the system (1),
(4)

gs gv
n(x1 ) = n(x4 ) = fFD (k, αεF ) dk,
(2π)2 R2

gs gv
p(x1 ) = p(x4 ) = fFD (k, –αεF ) dk.
(2π)2 R2

2.3 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

np – ne pe
R= √ , (5)
τ (n + p + 2 ne pe )

where ne and pe 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

∂φ
E=– .
∂x

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 Boltz-
mann 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 sit-
uation adopted here, in [8].

2.4 Poisson’s equation

The system (1) is coupled with the 2D Poisson equation for the electrostatic potential,
solved in the section of the device

∇ · ( ∇φ) = h(x, y), (6)

where
⎧
⎨e(n(x) – p(x) – n )/t , if (x, y) ∈ [x , x ] × [y , y ],
imp gr 1 4 2 3
h(x, y) = (7)
⎩0, otherwise.
Nastasi and Romano Journal of Mathematics in Industry (2022) 12:4 Page 5 of 11

The dielectric function is given by

⎧
⎨ , if y ∈ [y , y ],
gr 2 3
=
⎩ ox , otherwise,

with gr and ox dielectric constants in graphene and oxide respectively. nimp 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 en-
closed by the parallelepiped of base the area of the graphene sheet and height tgr .
Equation (6) is augmented with the following boundary conditions

φ=0 at y ∈ [y1 , y4 ], x = x1 ,

φ = Vb at y ∈ [y1 , y4 ], x = x4 ,

φ = VGd at y = y1 , x ∈ [x2 , x3 ],

φ = VGu at y = y4 , x ∈ [x2 , x3 ],

∇ν φ = 0 at the remaining part of

the boundary.

Here Vb is the bias voltage, VGu is the upper gate-source potential, VGd is the down gate-
source potential. VGu and VGd include the (subtracted) flat-band voltages. We have de-
noted by ∇ν the external normal derivative.

3 Numerical method
In order to solve the system (1)–(2) we adopt a finite difference scheme. A uniform mesh
of Nx grid points and size x is introduced in [x1 , x3 ] and for y = ygr . We denote by ui the
numerical approximation of a generic function u(x) at the node i. First we discretize the
equations (1) at the midpoint xi+1/2 of each interval [xi , xi+1 ], that is

i+ 12 i– 12 i+ 12 i– 12
Jn – Jn i Jp – Jp
= eR , = –eRi . (8)
x x

Then we discretize the equations for currents (2)

i+ 12 i+ 1 1 (εF(n) )j+1 – (εF(n) )j

Jn ≈ μn 2 ni+ 2 ,
x
(9)
(p) (p)
i+ 1 i+ 1 1 (εF )j+1 – (εF )j
Jp 2 ≈ μp 2 pi+ 2 ,
x

i+ 1 i+ 1 1 1
where μn 2 , μp 2 , ni+ 2 and pi+ 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.
i+ 1 i+ 1 1 1
We remark that the quantities μn 2 , μp 2 depend on ni+ 2 and pi+ 2 and the latters depend
(p)
also on εF(n) and εF through (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.
Nastasi and Romano Journal of Mathematics in Industry (2022) 12:4 Page 6 of 11

1. n(1) = p(1) = ni for each x ∈ [x1 , x4 ], where ni is the intrinsic graphene concentration,
and get ϕ (1) solving the Poisson equation (6)–(7).
2. For each r > 1:
(p)
(a) find (εF(n) )(r) and (ε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

Jn(r+1) – Jn(r)  Jp(r+1) – Jp(r) 

en = , ep = , (10)
Jn(r+1)  Jp(r+1) 

where Jn(r) and Jp(r) are the electron and hole currents at the rth iteration. As stopping cri-
terion en and ep 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].

4 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 (SiO2 ) 10 nm. The
graphene thickness tgr , 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 VG in the simulated cases.
In Eqs. (6)–(7) we adopt the values gr = 3.3 0 and ox = 3.6 0 , where 0 is the dielectric
constant in the vacuum, and nimp = 2.5 · 103 μm–2 .
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 Vb and gate voltages VG , the charac-
teristic 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 quan-
titative differences, in particular at the lower bias voltages considered in the simulations.
For example in the case Vb = 0.1 V, when VG = 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 Vb = 0.2 V and VG = –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.
Nastasi and Romano Journal of Mathematics in Industry (2022) 12:4 Page 7 of 11

Figure 2 Characteristic curves of the simulated GFET in linear scale (left) and logarithmic one (right). The
adopted model is (1), (2)

Figure 3 Characteristic curves of the simulated GFET in linear scale (left) and logarithmic one (right). The
adopted model is (1), (4)

Figure 4 Electron density of the simulated GFET in the case Vb = 0.2 V and VG = –1.5 and 1.5 V. 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 = (ε–εD(v
) d(ε–εD ) dϕ
2 ,
F)
with ϕ azimuthal angle of k, the electron and hole densities can be written in terms of the
Nastasi and Romano Journal of Mathematics in Industry (2022) 12:4 Page 8 of 11

Figure 5 Hole density of the simulated GFET in the case Vb = 0.2 V and VG = –1.5 and 1.5 V. On the left the
adopted model is (1), (2); on the right it is (1), (4)

Figure 6 Electrostatic potential of the simulated GFET in the case Vb = 0.2 V and VG = –1.5 V (left) and
VG = 1.5 V (right)

Fermi integrals of order k

 +∞
1 χk
Fk (z) = dχ,
(k + 1) 0 1 + eχ–z

with (x) Euler gamma function, as follows

 2  
2 εF(n) (x) – εD (x)
kB T
n(x) = F1 , (11)
π vF kB T
   (p) 
2 kB T 2 ε (x) – εD (x)
p(x) = F1 – F . (12)
π vF kB T

(p)
Here η(n) (x) = εF(n) (x) – εD (x) and η(p) (x) = εF (x) – εD (x) are the chemical energy for elec-
trons and holes respectively. Therefore

 
Jn = μn n∇εF(n) = μn n∇ η(n) + εD = μn n∇η(n) – eμn n∇φ.
Nastasi and Romano Journal of Mathematics in Industry (2022) 12:4 Page 9 of 11

The second term is the standard drift one. By observing that F1 (z) is a strict monotone
function and recalling the well known properties of the Fermi integral

dFk (z)
= Fk–1 (z), k = 1, 2, . . . and F0 (z) = ln(1 + exp z),
dz

the first term is given by

d η(n) μn n
μn n∇η(n) = μn n ∇n = dn ∇n
dn (n) dη

from which the electron diffusion coefficient in the degenerate case

(n) (n)
μn n F1 ( kηB T ) F1 ( kηB T )
Ddn = dn
= μn kB T (n)
= Dn (n)
. (13)
dη(n) ln[1 + exp( kηB T )] ln[1 + exp( kηB T )]

Similarly the diffusion coefficient of holes in the degenerate case is related to that in the
non-degenerate case by

(p)
F1 ( –η
kB T
)
Ddp = Dp (p)
.
ln[1 + exp( –η
kB T
)]

Dd Dd
In Fig. 7 the ratios rn = Dnn and rp = Dpp are plotted along the graphene in the simulated
cased considered above. For electrons rn is considerably different from one near the con-
tacts while rp 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 ad-
equate in the simulation of the considered device. This casts serious doubts on the use of
the non-degenerate drift-diffusion models in GFET.

Figure 7 Ratios of the degenerate and non-degenerate diffusion coefficients in the case Vb = 0.2 V and
VG = –1.5 V (left) and VG = 1.5 V (right)
Nastasi and Romano Journal of Mathematics in Industry (2022) 12:4 Page 10 of 11

5 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 dif-
ferences are observed regarding the characteristic curves making the drift-diffusion mod-
els 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.
A further step forward could be the inclusion of the thermal effects (see for example
[17–20]) or quantum correction [21–23].

Acknowledgements
Not applicable.

Funding
The authors acknowledge the support from INdAM (GNFM) and from Università degli Studi di Catania, Piano della Ricerca
2020/2022 Linea di intervento 2 “QICT”. GN acknowledges the financial support from Progetto Giovani GNFM 2020
“Trasporto di cariche e fononi in strutture a bassa dimensione”.

Abbreviations
FET, field effect transistor; GFET, graphene field effect transistor; MOSFET, metal oxide semiconductor field effect transistor.

Availability of data and materials

Not applicable.

Declarations
Competing interests
The authors declare that they have no competing interests.

Authors’ contributions
Both authors contributed equally to the writing of the manuscript. All authors read and approved the final manuscript.

Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Received: 9 November 2021 Accepted: 13 January 2022

References
1. Nastasi G, Romano V. An efficient GFET structure. IEEE Trans Electron Devices. 2021;68:4729–34.
https://doi.org/10.1109/TED.2021.3096492.
2. Castro Neto AH, Guinea F, Peres NMR, Novoselov KS, Geim AK. The electronic properties of graphene. Rev Mod Phys.
2009;81:109–62. https://doi.org/10.1103/RevModPhys.81.109.
3. Schwierz F. Graphene transistors. Nat Nanotechnol. 2010;5:487–96. https://doi.org/10.1038/nnano.2010.89.
4. Ancona MG. Electron transport in graphene from a diffusion-drift perspective. IEEE Trans Electron Devices.
2010;57:681–9. https://doi.org/10.1109/TED.2009.2038644.
5. Jiménez D, Moldovan O. Explicit drain-current model of graphene field-effect transistors targeting analog and
radio-frequency applications. IEEE Trans Electron Devices. 2011;58:4049–52.
https://doi.org/10.1109/TED.2011.2163517.
6. Feijoo PC, Jimńez D, Cartoixà X. Short channel effects in graphene-based field effect transistors targeting
radio-frequency applications. 2D Mater. 2016;3:025036. https://doi.org/10.1088/2053-1583/3/2/025036.
7. Champlain JG. A first principles theoretical examination of graphene-based field effect transistors. J Appl Phys.
2011;109:084515. https://doi.org/10.1063/1.3573517.
8. Nastasi G, Romano V. A full coupled drift-diffusion-Poisson simulation of a GFET. Commun Nonlinear Sci Numer
Simul. 2020;87:105300. https://doi.org/10.1016/j.cnsns.2020.105300.
9. Jacoboni C. Theory of electron transport in semiconductors. 1st ed. Berlin: Springer; 2010.
10. Kittel C. Introduction to solid state physics. 7th ed. Hoboken: Wiley; 2005.
11. Selberherr G. Analysis and simulation of semiconductor devices. Vienna: Springer; 1984.
12. Dorgan VE, Bae M-H, Pop E. Mobility and saturation velocity in graphene on SiO2 . Appl Phys Lett. 2010;97:082112.
https://doi.org/10.1063/1.3483130.
13. Coco M, Majorana A, Romano V. Cross validation of discontinuous Galerkin method and Monte Carlo simulations of
charge transport in graphene on substrate. Ric Mat. 2017;66:201–20. https://doi.org/10.1007/s11587-016-0298-4.
Nastasi and Romano Journal of Mathematics in Industry (2022) 12:4 Page 11 of 11

14. Majorana A, Nastasi G, Romano V. Simulation of bipolar charge transport in graphene by using a discontinuous
Galerkin method. Commun Comput Phys. 2019;26:114–34. https://doi.org/10.4208/cicp.OA-2018-0052.
15. Majorana A, Mascali G, Romano V. Charge transport and mobility in monolayer graphene. J Math Ind. 2016;7:4.
https://doi.org/10.1186/s13362-016-0027-3.
16. Nastasi G, Romano V. Improved mobility models for charge transport in graphene. Commun Appl Ind Math.
2019;10:41–52. https://doi.org/10.1515/caim-2019-0011.
17. Coco M, Romano V. Simulation of electron–phonon coupling and heating dynamics in suspended monolayer
graphene including all the phonon branches. J Heat Transf. 2018;140:092404. https://doi.org/10.1115/1.4040082.
18. Mascali G. A hydrodynamic model for silicon semiconductors including crystal heating. Eur J Appl Math.
2015;26:447–96. https://doi.org/10.1017/S0956792515000157.
19. Mascali G, Romano V. Charge transport in graphene including thermal effects. SIAM J Appl Math. 2017;77:593–613.
https://doi.org/10.1137/15M1052573.
20. Mascali G, Romano V. Exploitation of the maximum entropy principle in mathematical modeling of charge transport
in semiconductors. Entropy. 2017;19:36. https://doi.org/10.3390/e19010036.
21. Luca L, Romano V. Quantum corrected hydrodynamic models for charge transport in graphene. Ann Phys.
2019;406:30–53. https://doi.org/10.1016/j.aop.2019.03.018.
22. Barletti L, Cintolesi C. Derivation of isothermal quantum fluid equations with Fermi–Dirac and Bose–Einstein statistics.
J Stat Phys. 2012;148:353–86. https://doi.org/10.1007/s10955-012-0535-5.
23. Camiola VD, Luca L, Equilibrium RV. Wigner function for fermions and bosons in the case of a general energy
dispersion relation. Entropy. 2020;22:1023. https://doi.org/10.3390/e22091023.