Lecture Notes in Mathematics 1823
Editors:
J.M.Morel,Cachan
F.Takens,Groningen
B.Teissier,Paris
Subseries:
Fondazione C.I.M.E.,Firenze
Adviser:Pietro Zecca
3
Berlin
Heidelberg
NewYork
Hong Kong
London
Milan
Paris
Tokyo
A.M.Anile W.Allegretto C.Ringhofer
Mathematical
Problems in
Semiconductor Physics
Lectures given at the
C.I.M.E.Summer School
held in Cetraro,Italy,
July 1522,1998
With the collaboration of
G.Mascali and V.Romano
Editor:A.M.Anile
1 3
Authors and Editors
Angelo Marcello Anile
Dipartimento di Matematica
Universit`a di Catania
Viale A.Doria 6
95125 Catania,Italy
email:anile@dmi.unict.it
Walter Allegretto
Department of Mathematical and
Statistical Sciences
Alberta University
Edmonton AB T6G 2G1
Canada
email:retl@retl.math.ualberta.ca
Christian Ringhofer
Department of Mathematics
Arizona State University
Tempe,Arizona 852871804,USA
email:ringhofer@asu.edu
CataloginginPublication Data applied for
Bibliographic information published by Die Deutsche Bibliothek
Die Deutsche Bibliothek lists this publication in the Deutsche Nationalbibliografie;
detailed bibliographic data is available in the Internet at http://dnb.ddb.de
Mathematics Subject Classification (2000):82D37,80A17,65Z05
ISSN00758434
ISBN3540408029 SpringerVerlag Berlin Heidelberg NewYork
This work is subject to copyright.All rights are reserved,whether the whole or part of the material is
concerned,specif ically the rights of translation,reprinting,reuse of illustrations,recitation,broadcasting,
reproduction on microf ilmor in any other way,and storage in data banks.Duplication of this publication
or parts thereof is permitted only under the provisions of the German Copyright Lawof September 9,1965,
in its current version,and permission for use must always be obtained fromSpringerVerlag.Violations are
liable for prosecution under the German Copyright Law.
SpringerVerlag Berlin Heidelberg NewYork a member of BertelsmannSpringer
Science + Business Media GmbH
springer.de
c
SpringerVerlag Berlin Heidelberg 2003
Printed in Germany
The use of general descriptive names,registered names,trademarks,etc.in this publication does not imply,
even in the absence of a specif ic statement,that such names are exempt fromthe relevant protective laws
and regulations and therefore free for general use.
Typesetting:Cameraready T
E
X output by the authors
SPIN:10952481 41/3142/du  543210  Printed on acidfree paper
Preface
The increasing demand on ultra miniturized electronic devices for ever im
proving performances has led to the necessity of a deep and detailed under
standing of the mathematical theory of charge transport in semiconductors.
Because of their very short dimensions of charge transport,these devices must
be described in terms of the semiclassical Boltzmann equation coupled with
the Poisson equation (or some phenomenological consequences of these equa
tions) because the standard approach,which is based on the celebrated drift
diﬀusion equations,leads to very inaccurate results whenever the dimensions
of the devices approach the carrier mean free path.
In some cases,such as for very abrupt heterojunctions in which tunneling
occurs it is even necessary to resort to quantum transport models (e.g.the
WignerBoltzmannPoisson system or equivalent descriptions).
These sophisticated physical models require an appropriate mathematical
framework for a proper understanding of their mathematical structure as well
as for the correct choice of the numerical algorithms employed for computa
tional simulations.
The resulting mathematical problems have a broad spectrum of theoretical
and practical conceptually interesting aspects.
From the theoretical point of view,it is of paramount interest to investigate
wellposedness problems for the semiclassical Boltzmann equation (and also for
the quantum transport equation,although this is a much more diﬃcult case).
Another problem of fundamental interest is that of the hydrodynamical limit
which one expects to be quite diﬀerent from the NavierStokesFourier one,
since the collision operator is substantially diﬀerent from the one in rareﬁed
gas case.
From the application viewpoint it is of great practical importance to study
eﬃcient numerical algorithms for the numerical solution of the semiclassical
Boltzmann transport equation (e.g spherical harmonics expansions,Monte
Carlo method,method of moments,etc.) because such investigations could
have a great impact on the performance of industrial simulation codes for
VI Preface
TCAD (Technology Computer Aided Design) in the microelectronics indus
try.
The CIME summer course entitled MATHEMATICAL PROBLEMS
IN SEMICONDUCTOR PHYSICS dealt with this and related ques
tions.It was addressed to researchers (either PhD students,young postdocs
or mature researchers from other areas of applied mathematics) with a strong
interest in a deep involvement in the mathematical aspects of the theory of
carrier transport in semiconductor devices.
The course took place in the period 1522 July 1998 on the premises of the
Grand Hotel San Michele di Cetraro (Cosenza),located at a beach of astound
ing beauty in the Magna Graecia part of southern Italy.The Hotel facilities
were more than adequate for an optimal functioning of the course.About 50
“students”,mainly from various parts of Europe,participated in the course.
At the end of the course,in the period 2324 July 1998,a related workshop of
the European Union TMR(Training and Mobility of Researchers) on “Asymp
totic Methods in Kinetic Theory” was held in the same place and several of
the participants stayed for both meetings.Furthermore the CIME course was
considered by the TMR as one of the regular training schools for the young
researchers belonging to the network.
The course developed as follows:
• W.Allegretto delivered 6 lectures on analytical and numerical problems
for the driftdiﬀusion equations and also on some recent results concerning
the electrothermal model.In particular he highlighted the relationship
with integrated sensor modeling and the relevant industrial applications,
inducing a considerable interest in the audience.
• F.Poupaud delivered 6 lectures on the rigorous derivation of the quan
tum transport equation in semiconductors,utilizing recent developments
on Wigner measures introduced by G´erard,in order to obtain the semi
classical limit.His lectures,in the French style of pure mathematics,were
very clear,comprehensive and of advanced formal rigour.The lectures were
particularly helpful to the young researchers with a strong background in
Analysis because they highlighted the analytical problems arising fromthe
rigorous treatment of the semiclassical limit.
• C.Ringhofer delivered 6 lectures which consisted of an overview of the
state of the art on the models and methods developed in order to study
the semiclassical Boltzmann equation for simulating semiconductor de
vices.He started his lectures by recalling the fundamentals of semicon
ductor physics then introduced the methods of asymptotic analysis in or
der to obtain a hierarchy of models,including:driftdiﬀusion equations,
energy transport equations,hydrodynamical models (both classical and
quantum),spherical harmonics and other kinds of expansions.His lec
tures provided comprehensive review of the modeling aspects of carrier
transport in semiconductors.
Preface VII
d) D.Levermore delivered 6 lectures on the mathematical foundations and
applications of the moment methods.He presented in detail and depth the
concepts of exponential closures and of the principle of maximum entropy.
In his lectures he gave several physical examples of great interest arising
from rareﬁed gas dynamics,and pointed out how the method could also
be applied to the semiclassical Boltzmann equation.He highlighted the re
lationships between the method of moments and the mathematical theory
of hyperbolic systems of conservation laws.
During the course several seminars on specialized topics were given by lead
ing specialists.Of particular interest were these of P.Markowich (codirector
of the course) on the asymptotic limit for strong ﬁeds,of P.Pietra on the
numerical solution of the quantum hydrodynamical model,of A.Jungel on
the entropy formulation of the energy transport model,of O.Muscato on the
Monte Carlo validation of hydrodynamical models,of C.Schmeiser on ex
tended moment methods,of A.Arnold on the WignerPoisson system,and of
A.Marrocco on the mixed ﬁnite element discretization of the energy transport
model.
A.M.Anile
CIME’s activity is supported by:
Ministero dell’Universit`a Ricerca Scientiﬁca e Tecnologica;
Consiglio Nazionale delle Ricerche;E.U.under the Training and Mobility of
Researchers Programme.
Contents
Recent Developments in Hydrodynamical Modeling of
Semiconductors
A.M.Anile,G.Mascali and V.Romano............................1
1 Introduction..................................................1
2 General Transport Properties in Semiconductors..................2
3 HTheorem and the Null Space of the Collision Operator...........5
4 Macroscopic Models...........................................7
4.1 Moment Equations........................................7
4.2 The Maximum Entropy Principle............................8
5 Application of MEP to Silicon..................................11
5.1 Collision Term in Silicon...................................11
5.2 Balance Equations and Closure Relations.....................13
5.3 Simulations in Bulk Silicon.................................15
5.4 Simulation of a n
+
−n −n
+
Silicon Diode...................21
5.5 Simulation of a Silicon MESFET............................26
6 Application of MEP to GaAs...................................34
6.1 Collision Term in GaAs....................................34
6.2 Balance Equations and Closure Relations.....................36
6.3 Simulations in Bulk GaAs..................................38
6.4 Simulation a GaAs n
+
−n −n
+
Diode.......................43
6.5 Gunn Oscillations.........................................45
References......................................................54
DriftDiﬀusion Equations and Applications
W.Allegretto....................................................57
1 The Classical Semiconductor DriftDiﬀusion System...............57
1.1 Derivation................................................57
1.2 Existence.................................................58
1.3 Uniqueness and Asymptotics................................63
2 Other DriftDiﬀusion Equations................................66
2.1 Small Devices.............................................66
X Contents
2.2 C
α,α/2
Solutions and the Amorphous Silicon System...........68
2.3 Avalanche Generation......................................70
3 Degenerate Systems...........................................70
3.1 Degenerate Problems:Limit Case of the Hydrodynamic Models.70
3.2 Temperature Eﬀects.......................................73
3.3 Degenerate Problems:Thermistor Equations and
Micromachined Structures..................................74
4 Related Problems.............................................80
5 Approximations,Numerical Results and Applications..............82
References......................................................89
Kinetic and Gas – Dynamic Models for Semiconductor
Transport
Christian Ringhofer..............................................97
1 MultiBody Equations and Eﬀective Single Electron Models........98
1.1 Eﬀective Single Particle Models – The BBGKY Hierarchy......101
1.2 The Relation Between Classical and Quantum Mechanical Models104
2 Collisions and the Boltzmann Equation..........................107
3 Diﬀusion Approximations to Kinetic Equations....................111
3.1 Diﬀusion Limits:The Hilbert Expansion.....................113
3.2 The Drift Diﬀusion Equations:..............................114
3.3 The Energy Equations:.....................................115
3.4 The Energy Transport – or SHE Model......................116
3.5 Parabolicity..............................................118
4 Moment Methods and Hydrodynamic Models....................120
References......................................................130
Recent Developments in Hydrodynamical
Modeling of Semiconductors
A.M.Anile,G.Mascali and V.Romano
Dipartimento di Matematica e Informatica,
Universit`a di Catania
viale A.Doria 6  95125 Catania,Italy
anile@dmi.unict.it,mascali@dmi.unict.it,romano@dmi.unict.it
Summary.We present a review of recent developments in hydrodynamical mod
eling of charge transport in semiconductors.We focus our attention on the models
for Si and GaAs based on the maximum entropy principle which,in the framework
of extended thermodynamics,leads to the deﬁnition of closed systems of moment
equations starting from the Boltzmann transport equation for semiconductors.
Both the theoretical and application issues are examined.
1 Introduction
Enhanced functional integration in modern electron devices requires an in
creasingly accurate modeling of energy transport in semiconductors in order
to describe highﬁeld phenomena such as hot electron propagation,impact
ionization and heat generation.In fact the standard driftdiﬀusion models
cannot cope with highﬁeld phenomena since they do not comprise energy as
a dynamical variable.
Furthermore,for many applications in optoelectronics it is necessary to
describe the transient interaction of electromagnetic radiation with carriers in
complex semiconductor materials.Since the characteristic times are of order of
the electron momentumor energy ﬂux relaxation times,some higher moments
of the carrier distribution function must be necessarily involved.These are the
main reasons why more general models have been sought which incorporate
energy as a dynamical variable and whose validity,at variance with the drift
diﬀusion model,is not restricted to quasistationary situations.
These models are,loosely speaking,called hydrodynamical models and
they are usually derived by suitable truncation procedures,from the inﬁ
nite hierarchy of the moment equations of the Boltzmann transport equation.
However,most of these suﬀer from serious theoretical drawbacks due to the
A.M.Anile,W.Allegretto,C.Ringhofer:LNM 1821,A.M.Anile (Ed.),pp.1–56,2003.
cSpringerVerlag Berlin Heidelberg 2003
2 A.M.Anile,G.Mascali and V.Romano
ad hoc treatment of the closure problem [1].Recently,in the case of silicon
semiconductors,a moment approach has been introduced [2,3] (see also [4] for
a complete review) in which the closure procedure is based on the maximum
entropy principle,while the conduction bands in the proximity of the local
minima are described by the Kane dispersion relation.Later on,[5,6],the
same approach has been employed for GaAs.In this case both the Γvalley
and the four equivalent Lvalleys have been considered.Therefore electrons in
the conduction band have been treated as a mixture of two ﬂuids,one repre
senting electrons in the Γvalley and the other electrons in the four equivalent
Lvalleys.
Both in the Si and in the GaAs case,the models comprise the balance
equations of electron density,energy density,velocity and energy ﬂux.The
only diﬀerence is that for GaAs both electron populations are taken into
account.These equations are coupled to the Poisson equation for the electric
potential.Apart from the Poisson equation,the system is hyperbolic in the
physically relevant region of the ﬁeld variables.
In this paper we present a general overview of the theory underlying hydro
dynamical models.In particular we investigate in depth the closure problem
and present various applications both to bulk materials and to electron de
vices.
The considerations and the results reported in the paper are exclusively
concerned with silicon and gallium arsenide.
2 General transport properties in semiconductors
Semiconductors are characterised by a sizable energy gap between the va
lence and the conduction bands.Upon thermal excitation,electrons from the
valence band can jump to the conduction band leaving behind holes (in the
language of quasiparticles).Therefore the transport of charge is achieved both
through negatively charged (electrons) and positively charged (holes) carri
ers.The conductivity is enhanced by doping the semiconductor with donor or
acceptor materials,which respectively increase the number of electrons in the
conduction band or that of holes in the valence band.Therefore it is clear why
the energy band structure plays a very important role in the determination of
the electrical properties of the material.The energy band structure of crys
tals can be obtained at the cost of intensive numerical calculations (and also
semiphenomenologically) by the quantum theory of solids [7].However,for
most applications,a simpliﬁed description,based on simple analytical mod
els,is adopted to describe charge transport.In this paper we will be essentially
concerned with unipolar devices in which the current is due to electrons (semi
conductors doped with donor materials).Electrons which mainly contribute
to the charge transport are those with energy in the neighborhoods of the
lowest conduction band minima,each neighborhood being called a valley.In
Recent Developments in Hydrodynamical Modeling of Semiconductors 3
silicon,there are six equivalent ellipsoidal valleys along the main crystallo
graphic directions ∆ at about 85 % from the center of the ﬁrst Brilloiun zone,
near the X points,which,for this reason,are termed the Xvalleys.In GaAs
there is an absolute minimum at the center of the Brillouin zone,the Γpoint,
and local minima at the Lpoints along the Λ cristallographic orientations.
As mentioned above,in the simpliﬁed description employed,the energy
in each valley is represented by analytical approximations.Among these,the
most common are the parabolic and the Kane dispersion relation.
In the isotropic parabolic band approximation,the energy E
A
of the A
valley,measured from the bottom of the valley
E
A
,has an expression similar
to that of a classical free particle
E
A
(k
A
) =
2
k
A

2
2m
∗
A
.(1)
In this approximation k
A
,the electron wave vector,is assumed to vary in all
R
3
,m
∗
A
is the eﬀective electron mass in the Avalley and the reduced Planck
constant.
A more appropriate analytical approximation,which takes into account
the nonparabolicity at high energy,is given by the Kane dispersion relation
E
A
(k
A
) [1 +α
A
E
A
(k
A
)] =
2
k
2
2m
∗
A
,k ∈ R,(2)
where α
A
is the non parabolicity parameter.
The electron velocity v(k)
1
in a generic band or valley depends on the
energy E by the relation
v(k) =
1
∇
k
E.
Explicitly we get for parabolic band
v
i
=
k
i
m
∗
,(3)
while in the approximation of the Kane dispersion relation
v
i
=
k
i
m
∗
[1 +2αE(k)]
.(4)
In the semiclassical kinetic approach the charge transport in semiconduc
tors is described by the Boltzmann equation.For electrons in the conduction
band it reads
∂f
∂t
+v
i
(k)
∂f
∂x
i
−
eE
i
∂f
∂k
i
= C[f],(5)
1
the valley index has been omitted for simplicity
4 A.M.Anile,G.Mascali and V.Romano
where f(x,k,t) is the electron distribution function and C[f] represents the
eﬀects due to scattering with phonons,impurities and with other electrons.
In a multivalley description one has to consider a transport equation for each
valley.
The electric ﬁeld,E,is calculated by solving the Poisson equation for the
electric potential φ
E
i
= −
∂φ
∂x
i
,(6)
∇(∇φ) = −e(N
+
−N
−
−n),(7)
N
+
and N
−
denote the donor and acceptor density respectively (which depend
only on the position), the dielectric constant and n the electron number
density
n =
fdk.
The equations (5)(7) constitute the BoltzmannPoisson system that is the
basic semiclassical model of electron transport in semiconductors.
The main scattering mechanisms in a semiconductor are the electron
phonon interaction,the interaction with impurities,the electronelectron scat
terings and the interaction with stationary imperfections of the crystal as
vacancies,external and internal crystal boundaries.In many situations the
electronelectron collision term can be neglected since the electron density is
not too high.However in the case of high doping,electronelectron collisions
must be taken into account because they might produce sizable eﬀects.Re
taining the electronelectron collision term greatly increases the complexity
of the collision operator on the RHS of the semiclassical Boltzmann equation.
In fact the collision operator for the electronelectron scattering is a highly
nonlinear one,being quartic in the distribution function.
After a collision the electron can remain in the same valley (intravalley
scattering) or be drawn in another valley (intervalley scattering).
The general form of the collision operator C[f] for each type of scattering
mechanism is
C[f] =
P(k
,k)f(k
)
1 −4π
3
f(k)
−P(k,k
)f(k)
1 −4π
3
f(k
)
dk.(8)
The ﬁrst term in (8) represents the gain and the second one the loss.The
terms 1 − 4π
3
f(k) account for the Pauli exclusion principle.P(k,k
) is the
transition probability from the state k to the state k
.
Under the assumption that the electron gas is dilute,the collision operator
can be linearized with respect to f and becomes
C[f] =
[P(k
,k)f(k
) −P(k,k
)f(k)] dk.(9)
As we shall see at equilibrium the electron distribution must obey the Fermi
Dirac statistics
Recent Developments in Hydrodynamical Modeling of Semiconductors 5
f
eq
=
exp
−
E −µ
k
B
T
L
+1
−1
,
k
B
being the Boltzmannn constant,µ the chemical potential and T
L
the lattice
temperature which will be taken as constant.
In the dilute case,one can consider the maxwellian limit of the FermiDirac
distribution
f
eq
≈ exp
−
E −µ
k
B
T
L
.
In both cases from the principle of detailed balance [8],it follows that
P(k
,k) = P(k,k
) exp
−
E −E
k
B
T
L
,(10)
where E = E(k) and E
= E(k
).
3 Htheorem and the null space of the collision operator
In [9,10,11] an Htheorem has been derived for the physical electron
phonon operator in the homogeneous case without electric ﬁeld.The same
problem has also been discussed in [12] in the parabolic case.
Here we review the question in the case of an arbitrary form of the energy
band and in the presence of an electric ﬁeld,neglecting the electronelectron
interaction and assuming the electron gas suﬃciently dilute to neglect the
degeneracy eﬀects.By following [13] a physical interpretation of the results is
suggested.
The transition probability from the state k to the state k
has the general
form [14]
P(k,k
) = G(k,k
) [(N
B
+1)δ(E
−E +ω
q
) +N
B
δ(E
−E −ω
q
)] (11)
where δ(x) is the Dirac distribution and G(k,k
) is the socalled overlap factor
which depends on the band structure and the particular type of interaction
[14] and enjoys the properties
G(k,k
) = G(k
,k) and G(k,k
) ≥ 0.
N
B
is the phonon distribution which obeys the BoseEinstein statistics
N
B
=
1
exp(ω
q
/k
B
T
L
) −1
,(12)
where ω
q
is the phonon energy.
Given an arbitary function ψ(k) for which the following integrals exist,
the chain of identities [9,10,11]
6 A.M.Anile,G.Mascali and V.Romano
C[f]ψ(k)dk =
[P(k
,k)f(k
) −P(k,k
)f(k)] ψ(k)dkdk
=
B
2
P(k,k
)f(k) (ψ(k
) −ψ(k)) dkdk
=
B
2
G(k,k
) [(N
B
+1)δ(E
−E +ω
q
) +N
B
δ(E
−E −ω
q
)] ×
f(k) (ψ(k
) −ψ(k)) dkdk
=
B
2
G(k,k
)δ(E
−E−ω
q
) [(N
B
+1)f(k
)−N
B
f(k)] (ψ(k)−ψ(k
)) dkdk
holds.By following [11] if we set without loss of generality
f(k) = h(k) exp
−
E
k
B
T
L
,
and in analogy with the case of a simple gas we take
ψ(k) = k
B
log h(k),
by using the deﬁnition of δ(x),one has
k
B
C[f] log h(k)dk = k
B
G(k,k
)δ(E
−E −ω
q
)N
B
exp
−
E
k
B
T
L
(h(k
) −h(k)) (log h(k) −log h(k
)) dkdk
≤ 0.(13)
Therefore along the characteristics of eq.(5)
−
log h(k)
df
dt
dk = −
C[f] log h(k)dk ≥ 0.
This implies that
Ψ = k
B
log h(k) df
dk = k
B
f log f −f +
E
k
B
T
L
f
dk.(14)
can be considered as a Liapunov function for the BoltzmannPoisson system
(5)(7).The ﬁrst two terms are equal to the opposite of the entropy arising in
the classical limit of a Fermi gas,while the last term is due to the presence of
the phonons.Ψ represents the nonequilibrium counterpart of the equilibrium
Helmholtz free energy,divided by the lattice temperature.It is well known in
thermostatics that for a body kept at constant temperature and mechanically
insulated,the equilibrium states are minima for Ψ.
A strictly related problem is the one of determining the null space of the
collision operator.It consists in ﬁnding the solutions of the equation C(f) =
0.The resulting distribution functions represent the equilibrium solutions.
Physically one expects that,asymptotically in time,the solution to a given
initial value problem will tend to such a solution.
Recent Developments in Hydrodynamical Modeling of Semiconductors 7
The problemof determining the null space for the physical electronphonon
operator was tackled and solved in general in [11] where it was proved that
the equilibrium solutions are not only the FermiDirac distributions but form
an inﬁnite sequence of functions of the kind
f(k) =
1
1 +h(k) expE(k)/k
B
T
L
(15)
where h(E) = h(E + ω
q
) is a periodic function of period ω
q
/n,n ∈ N.
This property implies a numerable set of collisional invariants and hence of
conservation laws.The physical meaning is that the density of electrons whose
energy E diﬀers froma given value u by a multiple of ω
q
is constant.However
if there are several types of phonons,as in the real physical cases,and their
frequencies are not commensurable,the kernel of the collision operator reduces
to the FermiDirac distribution.
4 Macroscopic models
4.1 Moment equations
Macroscopic models are obtained by taking the moments of the Boltzmann
transport equation.In principle,all the hierarchy of the moment equations
should be retained,but for practical purposes it is necessary to truncate it at
a suitable order N.Such a truncation introduces two main problems due to
the fact that the number of unknowns exceeds that of the equations:these are
i) the closure for higher order ﬂuxes;
ii) the closure for the production terms.
As in gasdynamics [15],multiplying eq.(5) by a suﬃciently regular func
tion ψ(k) and integrating over B,the ﬁrst Brillouin zone,one obtains the
generic moment equation
∂M
ψ
∂t
+
ψ(k)v
i
(k)
∂f
∂x
i
dk −
e
E
j
ψ(k)
∂
∂k
j
fdk =
ψ(k)C[f]dk,(16)
with
M
ψ
=
ψ(k)fdk,
the moment relative to the weight function ψ.
Since
ψ(k)
∂f
∂k
j
dk =
∂B
ψ(k)fndσ −
f
∂ψ(k)
∂k
j
dk,
with n outward unit normal ﬁeld on the boundary ∂B of the domain B and
dσ surface element of ∂B,eq.(16) becomes
8 A.M.Anile,G.Mascali and V.Romano
∂M
ψ
∂t
+
∂
∂x
i
fψ(k)v
i
(k)dk +
e
E
j
f
∂ψ(k)
∂k
j
dk −
∂B
ψ(k)fn
j
dσ
=
ψ(k)C(f)dk.(17)
The term
∂B
ψ(k)fndσ
vanishes both when B is expanded to R
3
,as in the parabolic and Kane ap
proximations,( because in order to guarantee the integrability condition f
must tend to zero suﬃciently fast as k
→ ∞ ) and when B is compact and
ψ(k) is periodic and continuous on ∂B.This latter condition is a consequence
of the periodicity of f on B and the symmetry of B with respect to the origin.
Various models employ diﬀerent expressions of ψ(k) and number of mo
ments.
4.2 The maximum entropy principle
The maximum entropy principle (hereafter MEP) leads to a systematic
way of obtaining constitutive relations on the basis of information theory (see
[16,17,18,19] for a review).
According to MEP if a given number of moments M
A
,A = 1,...,N,are
known,the distribution function which can be used to evaluate the unknown
moments of f,corresponds to the extremal,f
ME
,of the entropy functional
under the constraints that it yields exactly the known moments M
A
ψ
A
f
ME
dk = M
A
.(18)
Since the electrons interact with the phonons describing the thermal vibrations
of the ions placed at the points of the crystal lattice,in principle we should
deal with a two component system (electrons and phonons).However,if one
considers the phonon gas as a thermal bath at constant temperature T
L
,only
the electron component of the entropy must be maximized.Moreover,by
considering the electron gas as suﬃciently dilute,one can take the expression
of the entropy obtained as limiting case of that arising in the Fermi statistics
s = −k
B
(f log f −f) dk.(19)
If we introduce the lagrangian multipliers Λ
A
,the problem of maximizing
s under the constraints (18) is equivalent to maximizing
˜s = Λ
A
M
A
−
ψ
A
fdk
−s,
Recent Developments in Hydrodynamical Modeling of Semiconductors 9
the Legendre transform of s,without constraints,
δ˜s = 0.
This gives
log f +
Λ
A
ψ
A
k
B
δf = 0.
Since the latter relation must hold for arbitrary δf,it follows
f
ME
= exp
−
1
k
B
Λ
A
ψ
A
.(20)
We stress that at variance with the monatomic gas,the integrability prob
lem due to the fact that the sign of the argument in the exponential is not
deﬁned,does not arise here because the moments are obtained by integrating
over the ﬁrst Brillouin zone,which is a compact set of R
3
.
In order to get the dependence of the Λ
A
’s on the M
A
’s,one has to invert
the constraints (18).Then by taking the moments of f
ME
and C[f
ME
],one
ﬁnds the closure relations for the ﬂuxes and the production terms appearing
in the balance equations.On account of the analytical diﬃculties this,in
general,can be achieved only with a numerical procedure.However,apart
from the computational problems,the balance equations are now a closed set
of partial diﬀerential equations and with standard considerations in extended
thermodynamics [16],it is easy to showthat they forma quasilinear hyperbolic
system.
Let us set
η(f) = −k
B
(f log f −f).
The entropy balance equation is obtained multiplying the equation (5) by
η
(f) = ∂
f
η(f) and afterwards integrating with respect to k,one has
∂
∂t
η(f)dk +
∂
∂x
i
η(f)v
i
dk −
e
E
i
η
(f)
∂
∂k
i
fdk =
η
(f)C[f]dk.
By taking into account the periodicity condition of f on the ﬁrst Brillouin
zone,the integral
η
(f)
∂
∂k
i
fdk =
∂η(f)
∂k
i
dk =
∂B
η(f)n
i
dk
vanishes and the entropy balance equations assumes the usual form
∂s
∂t
+
∂ϕ
i
∂x
i
= g,
with
ϕ
i
=
η(f)v
i
dk entropy ﬂux
10 A.M.Anile,G.Mascali and V.Romano
and
g =
η
(f)C[f]dk entropy production.
The electric ﬁeld contributes neither to the entropy production nor to the
entropy ﬂux.
Let us now rewrite the balance equations in the form
∂M
A
∂t
+
∂F
i
A
∂x
i
= G
A
(M
B
,E),(21)
where
F
i
A
=
fψ
A
(k)v
i
(k)dk,
G
A
(M
B
,E) =−
e
E
j
f
∂ψ
A
(k)
∂k
j
dk−
∂B
ψ
A
(k)fn
j
dσ
+
ψ
A
(k)C(f)dk.
It is easy to prove,by multiplying (21) by Λ
A
and taking the sum over
A,that the entropy balance equation is a consequence of the solution of (21)
with the MEP closure relations.More in particular the ﬁeld equations and
the entropy balance equations are related by the condition that
∂s
∂t
+
∂ϕ
i
∂x
i
−g −Λ
A
∂M
A
∂t
+
∂F
i
A
∂x
i
−G
A
(M
B
,E)
= 0,(22)
with the lagrangian multipliers Λ
A
being the same as those arising by em
ploying the maximum entropy principle.This implies [16] that
M
A
=
∂s
∂Λ
A
,F
i
A
=
∂ϕ
i
∂Λ
A
,(23)
with
s
= Λ
A
M
A
−s and ϕ
i
= Λ
A
F
i
A
−ϕ.
If
∂
2
s
∂Λ
A
∂Λ
B
is deﬁned in sign,one can globally invert [20] and express the moments M
A
as function of the lagrangian multipliers Λ
B
.As shown in [21] the previous
condition is equivalent to require that
f
ME
=
∂
∂χ
f
ME
(χ)
is deﬁned in sign,with χ = Λ
A
ψ
A
/k
B
.One trivally gets f
ME
< 0 and there
fore the balance equations (21) can be rewritten in terms of the lagrangian
multipliers as
Recent Developments in Hydrodynamical Modeling of Semiconductors 11
∂
2
s
∂Λ
A
∂Λ
B
∂Λ
B
∂t
+
∂
2
ϕ
i
∂Λ
A
Λ
B
∂Λ
B
∂x
i
= G
A
∂s
∂Λ
C
,E
.(24)
In this formit is immediate to recognize that the balance equations constitute
a symmetric quasilinear hyperbolic system [22].The main consequence of this
property is that according to a theorem due to Fisher and Marsden [23] the
Cauchy problem is wellposed for the system (24),at least in the simple case
where the electric ﬁeld is considered as an external ﬁeld.
5 Application of MEP to silicon
5.1 Collision term in Silicon
In Silicon the electrons which give contribution to the charge transport are
those in the six equivalent valleys around the six minima of the conduction
band.One assumes that the electron energy in each valley is approximated
by the Kane dispersion relation.Concerning the collision term,the electron
phonon scatterings which occur can be summarized as follows:
• scattering with intravalley acoustic phonons (approximately elastic);
• electronphonon intervalley inelastic scatterings,for which there are six
contributions:the three g
1
,g
2
,g
3
and the three f
1
,f
2
,f
3
optical and acous
tical intervalley scatterings [24]
m
e
electron rest mass
9.109510
−10
g
m
∗
eﬀective electron mass
0.32 m
e
T
L
lattice temperature
300
o
K
ρ
density
2330 g/cm
3
v
s
longitudinal sound speed
9.18 10
5
cm/sec
Ξ
d
acousticphonon deformation potential
9 eV
α
non parabolicity factor
0.5 eV
−1
r
relative dielectric constant
11.7
0
vacuum dieletric constant
1.24 × 10
−22
C/( eV cm)
Table 1.Values of the physical parameters used for silicon
12 A.M.Anile,G.Mascali and V.Romano
α
Z
f
ω(meV )
(D
t
K) (10
8
eV/cm)
g
1
1
12
0.5
g
2
1
18.5
0.8
g
3
1
61.2
11
f
1
4
19
0.3
f
2
4
47.4
2
f
3
4
59
2
Table 2.Coupling constants and phonon energies for the intervalley scatterings in
Si.
In the elastic case
P
(ac)
(k,k
) = K
ac
δ(E −E
),(25)
while for inelastic scatterings
P
α
(k,k
) = K
α
(N
(α)
B
+1)δ(E
−E +ω
α
) +N
(α)
B
δ(E
−E −ω
α
)
,(26)
where α = g
1
,g
2
,g
3
,f
1
,f
2
,f
3
,N
(α)
B
is the phonon equilibrium distribution
according to the BoseEinstein statistics (12) and ω
α
the phonon energy.
Besides these interactions,we will also consider the scattering with impu
rities,which is an elastic mechanism of interaction.Its transition rate reads
P
(imp)
(k
A
,k
A
) =
K
imp
[k
A
−k
A

2
+β
2
]
2
δ(E
A
−E
A
),(27)
where β is the inverse Debye length.
The parameters that appear in the scatterings rates can be expressed in
terms of physical quantities characteristic of the considered material
K
ac
=
k
B
T
L
Ξ
2
d
4π
2
ρv
2
s
,K
α
=
Z
f α
(D
t
K)
2
α
8π
2
ρω
α
,
K
imp
=
N
I
Z
2
q
4
4π
2
,β =
q
2
N
I
k
B
T
L
1/2
Recent Developments in Hydrodynamical Modeling of Semiconductors 13
where Ξ
d
is the deformation potential of acoustic phonons,ρ the mass
density of the semiconductor,v
s
the sound velocity of the longitudinal acoustic
mode,(D
t
K)
α
the deformation potential realtive to the interaction with the
α intervalley phonon and Z
f α
the number of ﬁnal equivalent valleys for the
considered intervalley scattering.N
I
and Z q are respectively the impurity
concentration and charge.
The deformation potentials are not known from calculations by means of
ﬁrst principles because the perturbation theory employed to evaluate the tran
sition probabilities is not able to calculate them from the quantum theory of
scattering.In all the simulators,even the Monte Carlo ones,these quantities
are considered as ﬁtting parameters.Their values depend on the approxima
tion used for the energy bands,on the speciﬁc characteristics of the material
and on the energy range of interest in the applications.
The values of all these quantities as well as the silicon bulk constants are
given in [25].For the sake of completeness we summarize them in tables 1 and
2.
5.2 Balance equations and closure relations
As already stated,the macroscopic balance equations are deduced by taking
the moments of the Boltzmann transport equation for electrons in semicon
ductors [8].We will consider the balance equations for density,momentum,
energy and energy ﬂux,which correspond to the weight functions 1,k,E,Ev
∂n
∂t
+
∂(nV
i
)
∂x
i
= 0,(28)
∂(nP
i
)
∂t
+
∂(nU
ij
)
∂x
j
+neE
i
= nC
P
i
,(29)
∂(nW)
∂t
+
∂(nS
j
)
∂x
j
+neV
k
E
k
= nC
W
,(30)
∂(nS
i
)
∂t
+
∂(nF
ij
)
∂x
j
+neE
j
G
ij
= nC
S
i
.(31)
This system is coupled to Poisson’s equation
E = −∇φ,∇(∇φ) = −e(N
+
−N
−
−n) (32)
Since we consider the unipolar case,the hole concentration will be not
included.The macrocopic quantities involved in the balance equations are
related to the one particle distribution function of electrons f(x,k,t) as follows
14 A.M.Anile,G.Mascali and V.Romano
n =
R
3
fdk is the electron density,
V
i
=
1
n
R
3
fv
i
dk is the average electron velocity,
W =
1
n
R
3
E(k)fdk is the average electron energy,
S
i
=
1
n
R
3
fv
i
E(k)dk is the energy ﬂux,
P
i
=
1
n
R
3
fk
i
dk = m
∗
V
i
+2αS
i
is the average crystal momentum,
U
ij
=
1
n
R
3
fv
i
k
j
dk is the ﬂux of crystal momentum,
G
ij
=
1
n
R
3
1
f
∂
∂k
j
(Ev
i
)dk,
F
ij
=
1
n
R
3
fv
i
v
j
E(k)dk is the ﬂux of energy ﬂux,
C
P
i
=
1
n
R
3
C[f]k
i
dkis the production of crystal momentum,
C
W
=
1
n
R
3
C[f]E(k)dk is the energy production,
C
S
i
=
1
n
R
3
C[f]v
i
E(k)dk is the production of the energy ﬂux,
These moment equations do not constitute a set of closed relations because
of the ﬂuxes and production terms.Therefore constitutive assumptions must
be prescribed.
If we assume as fundamental variables n,V
i
,W and S
i
,which have a
direct physical interpretation,the closure problem consists of expressing P
i
,
U
ij
,F
ij
,G
ij
and the moments of the collision term C
P
i
,C
W
and C
S
i
as
functions of n,V
i
,W and S
i
.
If we use MEP to get the closure relations,we face the problemof inverting
the constraints (18) with ψ
A
= 1,v,E,Ev.
This problem has been overcome in [2,3] with the ansatz of small
anisotropy for f
ME
,since Monte Carlo simulations for electron transport in
Si show that the anisotropy of f is small even far from equilibrium.
Formally a small anisotropy parameter δ has been introduced and explicit
constitutive equations have been obtained in [2] for the higher order ﬂuxes and
in [3] for the production terms up to second order in δ.However it has been
found in [26] that the ﬁrst order model is suﬃciently accurate for numerical
applications and avoids some irregularities due to nonlinearities as occur in
the parabolic band case [27].
Since the closure relations are an approximation of the exact MEP ones,
the hyperbolicity is not guaranteed,but must be checked.As proved in [26],
Recent Developments in Hydrodynamical Modeling of Semiconductors 15
the system (28)(31) is hyperbolic in the region n > 0,W ≥ W
0
,with W
0
=
3
2
k
B
T
L
equilibrium energy.
For the numerical integration we use the scheme developed in [28,29,30],
see appendix,which is based on the Nessyhau and Tadmor scheme [31,32].
5.3 Simulations in bulk silicon
The physical situation is represented by a silicon semiconductor with a
uniform doping concentration,which we assume suﬃciently low so that the
scatterings with impurities can be neglected.On account of the symmetry with
respect to translations,the solution does not depend on the spatial variables.
The continuity equation gives n = constant and from the Poisson equation
one ﬁnds that E is also constant.Therefore the remaining balance equations
reduce to the following set of ODEs for the motion along the direction of the
electric ﬁeld which is chosen as xdirection
d
dt
V = −
eE
m
∗
+
2αeEG
m
∗
+
c
11
m
∗
−2αc
21
V +
c
12
m
∗
−2αc
22
S,(33)
d
dt
W = −eV E +C
W
,(34)
d
dt
S = −eEG
(0)
+c
21
V +c
22
S,(35)
where V and S are the xcomponent of V and S,G = G
11
and the c
ij
,
i,j = 1,2 are production terms whose expressions can be found in [3].
As initial conditions for (33)(35) we take
V (0) = 0,(36)
W(0) =
3
2
k
B
T
L
,(37)
S(0) = 0.(38)
The stationary regime is reached in a few picoseconds.
The solutions of (33)(35) for several values of the applied electric ﬁeld are
reported in ﬁgs.1 (velocity),2 (energy) and 3 (energy ﬂux).
16 A.M.Anile,G.Mascali and V.Romano
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
time (ps)
velocity (107 cm/sec)
Fig.1.velocity (cm/sec) versus time (ps) for E = 10 kV/cm,30 kV/cm,50 kV/cm,
70 kV/cm,100 kV/cm,120 kV/cm,150 kV/cm
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
time (ps)
energy (eV)
Fig.2.energy (eV) versus time (ps) for the same values of the of the electric ﬁeld
as in ﬁgure 1
Recent Developments in Hydrodynamical Modeling of Semiconductors 17
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
0
0.2
0.4
0.6
0.8
1
1.2
1.4
time (ps)
energy flux (107 eV cm/sec)
Fig.3.energy ﬂux (eV cm/sec) versus time (ps) for the same values of the electric
ﬁeld as in ﬁgure 1
The typical phenomena of overshoot and saturation velocity are both qual
itatively and quantitatively well described (see [33] ﬁg.3.22 for a comparison
with results obtained by MC simulations).
Similar results were reported in [3],but there a diﬀerent modeling of the
collision terms has been considered and,moreover,instead of taking into ac
count all the intervalley and intravalley scatterings,mean values of the cou
pling constant Ξ and D
t
K have been introduced.The inclusion of all the
scattering (intervalley and intravalley) mechanisms notably improves the re
sults.
For the sake of completeness,the parabolic band case has been also inte
grated,ﬁgs.4,5,6.The diﬀerences,especially in the energy,with respect to
the Kane case,conﬁrm that the parabolic band is an oversimpliﬁcation of the
real band structure.
18 A.M.Anile,G.Mascali and V.Romano
Fig.4.velocity (cm/sec) versus time (ps) in the parabolic band case (dashed line)
and for the Kane dispersion relation
Recent Developments in Hydrodynamical Modeling of Semiconductors 19
0
0.5
1
1.5
2
0.035
0.04
0.045
0.05
0.055
0.06
0.065
time (ps)
energy (eV)
Electric field 10 kV/cm
0
0.5
1
1.5
2
0.02
0.04
0.06
0.08
0.1
0.12
0.14
0.16
0.18
time (ps)
energy (eV)
Electric field 30 kV/cm
0
0.5
1
1.5
2
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
time (ps)
energy (eV)
Electric field 50 kV/cm
0
0.5
1
1.5
2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
time (ps)
energy (eV)
Electric field 70 kV/cm
Fig.5.energy (eV) versus electric ﬁeld (kV/cm) in the parabolic band case (dashed
line) and for the Kane dispersion relation
20 A.M.Anile,G.Mascali and V.Romano
0
0.5
1
1.5
2
0
0.01
0.02
0.03
0.04
0.05
0.06
0.07
time (ps)
energy flux (107 eV cm/sec)
Electric field 10 kV/cm
0
0.5
1
1.5
2
0
0.05
0.1
0.15
0.2
0.25
time (ps)
energy flux (107 eV cm/sec)
Electric field 30 kV/cm
0
0.5
1
1.5
2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
time (ps)
energy flux (107 eV cm/sec)
Electric field 50 kV/cm
0
0.5
1
1.5
2
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
time (ps)
energy flux (107 eV cm/sec)
Electric field 70 kV/cm
Fig.6.energy ﬂux (eV cm/sec) versus electric ﬁeld (kV/cm) in the parabolic band
case (dashed line) and for the Kane dispersion relation
Recent Developments in Hydrodynamical Modeling of Semiconductors 21
5.4 Simulation of a n
+
−n −n
+
silicon diode
Here we present the simulation of a ballistic n
+
− n − n
+
silicon diode as
follos [26] (see also [34] for a comparison with MC data).The n
+
regions are
Fig.7.Schematic representation of a n
+
−n −n
+
diode
0.1µm long,while the various lengths of the channel are taken into account.
Moreover several doping proﬁles will be considered as reported in table 3.
Channel length
N
+
+
N
C
+
V
b
Test#
L
c
(µm)
(×10
17
cm
−3
)
(×10
17
cm
−3
)
(Volt)
1
0.4
5
0.02
2
2
0.3
10
0.1
1
3
0.2
10
0.1
1
Table 3.Length of the channel,doping concentration (respectively in the n
+
and
n regions) and applied voltage in the test cases for the diode
Initially the electron energy is that of the lattice in thermal equilibrium at
the temperature T
L
,the charges are averagely at rest and the density is equal
to the doping concentration
n(x,0) = N
+
(x),W(x,0) =
3
2
k
B
T
L
,V (x,0) = 0,S(x,0) = 0,
where V and S are the only relevant component of velocity and energy
ﬂux.
22 A.M.Anile,G.Mascali and V.Romano
Regarding the boundary conditions,in principle the number of indepen
dent conditions on each boundary should be equal to the number of char
acteristics entering the domain.However we impose,in analogy with similar
cases [26,27,35,36,37,38],a double number of boundary conditions.More
precisely,we give conditions for all the variables in each boundary,located at
x = 0 and x = L,the total length of the device,
n(0,t) = n(L,t) = N
+
+
(39)
∂
∂x
W(0,t) =
∂
∂x
W(L,t) = 0,(40)
∂
∂x
V (0,t) =
∂
∂x
V (L,t) = 0,(41)
∂
∂x
S(0,t) =
∂
∂x
S(L,t) = 0,(42)
φ(0) = 0 and φ(L) = V
b
,(43)
where V
b
is the applied bias voltage.In all the numerical solutions there
is no sign of spurious oscillations near the boundary,indicating that the con
ditions (39)(42) are in fact compatible with the solution of the problem.
The doping proﬁle is regularized according to the function
N
+
(x) = N
+
+
−d
0
tanh
x −x
1
s
−tanh
x −x
2
s
,
where s = 0.01µm,d
0
= N
+
+
(1 −N
C
+
/N
+
+
)/2,x
1
= 0.1µm,and x
2
= x
1
+
L
c
,with L
c
channel length.The total length of the device is L = L
c
+0.2µm.
A grid with 400 spatial nodes has been used.The stationary solution is
reached within a few picoseconds (about ﬁve),after a short transient with
wide oscillations.
As ﬁrst case we consider the test problem 1 (length of the channel 0.4
micron) with V
b
= 2 Volts.
Recent Developments in Hydrodynamical Modeling of Semiconductors 23
0
0.1
0.2
0.3
0.4
0.5
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
micron
velocity (10
7 cm/sec)
0
0.1
0.2
0.3
0.4
0.5
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
0.4
0.45
0.5
micron
energy (eV)
0
0.1
0.2
0.3
0.4
0.5
0
0.5
1
1.5
micron
energy flux (107 eV cm/sec)
0
0.1
0.2
0.3
0.4
0.5
−12
−10
−8
−6
−4
−2
0
2
4
micron
electric field (Volt/micron)
Fig.8.numerical results of the test case 1 after 5 picoseconds in the parabolic band
case (dashed line) and for the Kane dispersion relation (continuous line)
24 A.M.Anile,G.Mascali and V.Romano
0
0.1
0.2
0.3
0.4
0.5
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
micron
velocity (10
7 cm/sec)
0.1
0.2
0.3
0.4
0.5
0
0.05
0.1
0.15
0.2
0.25
micron
energy (eV)
0
0.1
0.2
0.3
0.4
0.5
0
0.1
0.2
0.3
0.4
0.5
0.6
micron
energy flux (107 eV cm/sec)
0
0.1
0.2
0.3
0.4
0.5
−8
−6
−4
−2
0
2
4
micron
electric field (Volt/micron)
Fig.9.numerical results of the test case 2 after 5 picoseconds in the parabolic band
case (dashed line) and for the Kane dispersion relation (continuous line)
Recent Developments in Hydrodynamical Modeling of Semiconductors 25
0
0.1
0.2
0.3
0.4
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
micron
velocity (107 cm/sec)
0
0.1
0.2
0.3
0.4
0
0.05
0.1
0.15
0.2
0.25
0.3
micron
energy (eV)
0.1
0.2
0.3
0.4
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
micron
energy flux (107 eV cm/sec)
0
0.1
0.2
0.3
0.4
−12
−10
−8
−6
−4
−2
0
2
4
micron
electric field (Volt/micron)
Fig.10.numerical results of the test case 3 after 5 picoseconds in the parabolic
band case (dashed line) and for the Kane dispersion relation (continuous line)
The simulation for the parabolic band approximation is also shown ( ﬁg.8
dashed line),but it is evident,like in the bulk case,that the results are rather
poor.
The other test cases have been numerically integrated with V
b
= 1 Volt
(ﬁg.s 9,10).
Recent Developments in Hydrodynamical Modeling of Semiconductors 27
Γ
D
= {(x,y):0 ≤ x ≤ 0.1,0.2 ≤ x ≤ 0.4,0.5 ≤ x ≤ 0.6,y = 0.2,}.
The other part of ∂Ω is labeled as Γ
N
.The boundary conditions are as
signed as follows:
n =
n
+
at source and drain
n
g
at gate
(44)
φ =
0 at the source
φ
g
at the gate
φ
b
at the drain
(45)
W = W
0
,V· t = 0,
n · ∇(V· n) = 0,S =
5
3
W
0
V,
on Γ
D
,(46)
n · ∇n = 0,n · ∇W = 0,n · ∇φ = 0,
n · ∇V
i
= 0,S =
5
3
WV
i = 1,2 on Γ
N
.(47)
Here ∇is the bidimensional spatial gradient operator while n and t are the
unit outward normal vector and the unit tangent vector to ∂Ω respectively.
n
+
is the doping concentration in the n
+
region and n
g
is the density at the
gate,which is considered to be a Schottky contact [42],
n
g
= 3.9 ×10
5
cm
−3
.
φ
b
is the bias voltage and φ
g
is the gate voltage.In all the simulations we
set φ
g
= −0.8V while Φ
b
= 1V.
In the standard hydrodynamical model considered in the literature (e.g.
[43,44]),the energy ﬂux S is not a ﬁeld variable and it is not necessary to
prescribe boundary conditions for it.The relations (46)
4
and (47)
5
are not
based on the microscopic boundary conditions for the distribution function,
but they may be justiﬁed [30] in a heuristic way with the same approach
followed in [45].The numerical scheme can be found in [30].
We start the simulation with the following initial conditions:
28 A.M.Anile,G.Mascali and V.Romano
n(x,y,0) = N
+
(x,y) −N
−
(x,y),W = W
0
=
3
2
k
B
T
L
,
V
i
= 0,S
i
= 0 i = 1,2
The main numerical problems arise fromthe discontinuous doping and the
boundary conditions at the Schottky barrier which give rise there to sharp
changes in the density of several orders of magnitude.The use of a shock
capturing scheme is almost mandatory for this problem.The stationary solu
tion is reached in a few picoseconds (less than ﬁve).The code takes about 9
minutes and 10 seconds in a PC with 1 Ghz PentiumIII microprocessor.After
the initial behaviour,the solution becomes smooth and no signs of spurious
oscillations are present.The numerical scheme seems suitably robust and it is
able to capture the main features of the solution.Only the Kane dispersion
relation will be considered here because the results obtained in the parabolic
band approximation are rather unsatisfactory when high electric ﬁelds are
involved,as shown in the previous section for the n
+
−n −n
+
diode.
The density is plotted in ﬁgure 12.As expected there is a depletion region
beneath the gate.Moreover one can see that the drain is less populated than
the source.
Concerning the energy (ﬁgure 13) there are steep variations near the gate
edges.The mean energy of the electrons reaches a maximum value of about
0.35 eV in the part of the gate closest to the source.
The results for the velocity are shown in ﬁgure 14.The higher values of
the xcomponent are at the edges of the gate contact.This happens also for
the ycomponent,but with a huge peak at the gate edge closest to the source.
The behaviour seems to indicate that there is a loss of regularity at the edge
of the gate.
The shape of the energy ﬂux (ﬁgure 15) is qualitatively similar to that of
the velocity.
Very large tangential and normal components of the electric ﬁeld (ﬁgure
16) are present again at the edges of the gate.For completeness the electric
potential is also presented in ﬁg.17.
The results are qualitatively similar to those presented in [40,41] for all
the variables except the ycomponent of the velocity in account of the huge
peak at the edge of the gate.
Recent Developments in Hydrodynamical Modeling of Semiconductors 29
30 A.M.Anile,G.Mascali and V.Romano
Recent Developments in Hydrodynamical Modeling of Semiconductors 31
32 A.M.Anile,G.Mascali and V.Romano
Recent Developments in Hydrodynamical Modeling of Semiconductors 33
34 A.M.Anile,G.Mascali and V.Romano
Fig.17.stationary solution (after 5 picoseconds) for the electric potential for Φ
b
=
1 Volt.
6 Application of MEP to GaAs
6.1 Collision term in GaAs
In GaAs as well,some interactions are intravalley transitions and others are
intervalley transitions,but at variance with silicon the Γvalley and the L
valleys are not equivalent.The collision term may be split according to
C
A
[f] = C
A
[f
A
,f
A
] +
B
C
A
[f
A
,f
B
].
The ﬁrst term represents the intravalley scatterings,the second one takes into
account the intervalley scatterings.As said,the main scattering mechanisms,
to which electrons are subjected in semiconductors,are the ones with phonons,
impurities,other electrons and stationary imperfections of the crystal such as
vacancies and external and internal crystal boundaries.
Since GaAs is a compound semiconductor,electrons interact with the
phonons not only because of the deformation of the crystal but also through
the polarization waves.
The coupling through the polarization waves is due to the permanent
electric dipole moment of the constituent ions in a compound material.This
Recent Developments in Hydrodynamical Modeling of Semiconductors 35
coupling can be mediated through both the optical and the acoustic branches,
but the latter contribution is marginal at room temperature in very pure
semiconductors.
In particular we consider the following scattering mechanisms:
• the acoustic phonon scattering,
• the nonpolar optical phonon scattering,
• the polar optical phonon scattering,
• the impurity scattering.
All these scattering mechanisms are intravalley,except the nonpolar op
tical phonon scattering.
In the case of acoustic phonon scattering in its elastic approximation,valid
when the thermal energy is much greater than that of the phonon involved in
the scattering,we have
P
(ac)
(k
A
,k
A
) = K
ac
δ(E
A
−E
A
),(48)
where δ is the Dirac delta function and K
ac
a physical parameter.
For the nonpolar optical phonon interaction,the transition rate is given
by the sum of an absorption and an emission term
P
(np)
(k
A
,k
B
) = Z
AB
K
np
N
(np)
δ(E
B
−E
A
−ω
np
)
+ (N
(np)
+1)δ(E
B
−E
A
+ω
np
))
,(49)
where K
np
is a physical parameter,ω
np
is the longitudinal optical phonon
energy,N
(np)
is the equilibrium nonpolar optical phonon BoseEinstein dis
tribution.
Z
AB
is the degeneracy of the ﬁnal valley (Table 4),that is the valley
the electron reaches after the scattering,with respect to the initial one.In
GaAs,nonpolar optical phonons contribute to the electron intervalley transfer
between two equivalent Lvalleys and between Γvalley and Lvalleys.
The polar optical phonon scattering plays a very important role in com
pound semiconductors.It is an intravalley inelastic process whose transition
rate is given by
P
(p)
(k
A
,k
A
) =
K
p
8π
2
k
A
−k
A

2
G(k
A
,k
A
)
N
(p)
δ(E
A
−E
A
−ω
p
)+
+(N
(p)
+1)δ(E
A
−E
A
+ω
p
)
,(50)
where K
p
is a physical parameter,ω
p
is the polar optical phonon energy,
N
(p)
the thermal equilibrium polar optical phonon number and the overlap
factor G is given by [14]
G(k
A
,k
A
) = (a
A
a
A
+c
A
c
A
e
A
· e
A
)
2
,(51)
36 A.M.Anile,G.Mascali and V.Romano
where e
A
and e
A
are unit vectors pointing in the directions of k
A
and k
A
and
a
A
=
1 +α
A
E
A
1 +2α
A
E
A
,a
A
=
1 +α
A
E
A
1 +2α
A
E
A
c
A
=
α
A
E
A
1 +2α
A
E
A
,c
A
=
α
A
E
A
1 +2α
A
E
A
.
At last,the transition rate for the scattering with impurities,is similar to that
for Si,see 27.The electronelectron scattering is not considered as well as the
scattering with imperfections
The values of the physical parameters appearing in the scattering rates
are reported in table 4.
m
∗
Γ
eﬀective electron mass in the Γvalley
0.067 ×m
e
m
∗
L
eﬀective electron mass in the Lvalley
0.35 ×m
e
T
L
lattice temperature
300 K
ρ
0
density
5.360 g/cm
3
v
s
longitudinal sound speed
5.24 × 10
5
cm/sec
α
Γ
non parabolicity factor in the Γvalley
0.611 eV
−1
α
L
non parabolicity factor in the Lvalley
0.242 eV
−1
r
relative dielectric constant
12.90
∞
relative dielectric constant at optical frequency range
10.92
Ξ
d
acousticphonon deformation potential
7 eV
D
t
K
nonpolar optical phonon deformation potential
10
9
eV/cm
ω
np
nonpolar optical phonon energy
0.03 eV
ω
p
polar optical phonon energy
0.03536 eV
E
Γ(0)
Γvalley bottom energy
0 eV
E
L(0)
Lvalley bottom energy
0.32 eV
Z
ΓL
degeneracy from Γ to L valleys
4
Z
LΓ
degeneracy from L to Γ valley
1
Z
LL
degeneracy from L to L valleys
3
Table 4.Values of the physical parameters used for GaAs
6.2 Balance equations and closure relations
One assumes,as for silicon,that each valley in the conduction band [7] is
described by the Kane dispersion approximation.
However in this case,at a kinetic level,the system is described by two
Boltzmann equations,one for the Γvalley and the other for one Lvalley.The
macroscopic balance equations are deduced,as usual,by taking the moments
of the Boltzmann transport equations.We consider the set of weight functions
necessary to get the macroscopic balance equations for densities,momenta,
Recent Developments in Hydrodynamical Modeling of Semiconductors 37
energies and energyﬂuxes of both Γ and Lvalley electrons,and the resulting
equations read
∂n
A
∂t
+
∂(n
A
V
i
A
)
∂x
i
= n
A
C
n
A
,(52)
∂(n
A
P
i
A
)
∂t
+
∂(n
A
U
ij
A
)
∂x
j
+n
A
eE
i
= n
A
C
P
i
A
,(53)
∂(n
A
W
A
)
∂t
+
∂(n
A
S
j
A
)
∂x
j
+n
A
eV
j
A
E
j
= n
A
C
W
A
,(54)
∂(n
A
S
i
A
)
∂t
+
∂(n
A
F
ij
A
)
∂x
j
+n
A
eE
j
G
ij
A
= n
A
C
S
i
A
,(55)
where
n
A
=
R
3
f
A
dk
A
is the electron density,
V
i
A
=
1
n
A
R
3
v
i
A
f
A
dk
A
the average electron velocity,
W
A
=
1
n
A
R
3
E
A
(k
A
)f
A
dk
A
the average electron energy,
S
i
A
=
1
n
A
R
3
v
i
A
E
A
(k
A
)f
A
dk
A
the average energy ﬂux,
P
i
A
=
1
n
A
R
3
k
i
A
f
A
dk
A
=m
∗
A
V
i
A
+2α
A
S
i
A
the average crystal momentum,
U
ij
A
=
1
n
A
R
3
v
i
A
k
j
A
f
A
dk
A
the average crystal momentum ﬂux,
G
ij
A
=
1
n
A
R
3
∂
∂k
j
A
1
v
i
A
E
A
(k)
f
A
dk
A
,
F
ij
A
=
1
n
A
R
3
v
i
A
v
j
A
E
A
(k
A
)f
A
dk
A
the average ﬂux of energy ﬂux,
C
n
A
=
1
n
A
R
3
C
A
[f]dk
A
the density production,
C
P
i
A
=
1
n
A
R
3
k
i
A
C
A
[f]dk
A
the crystal momentum production,
C
W
A
=
1
n
A
R
3
E
A
(k
A
)C
A
[f]dk
A
the energy production,
C
S
i
A
=
1
n
A
R
3
v
i
A
E
A
(k
A
)C
A
[f]dk
A
the energy ﬂux production.
All these quantities refer to electrons in the Avalley,A = Γ,L.
C
A
[f] are the scattering operators which appear in the electron transport
equations.
38 A.M.Anile,G.Mascali and V.Romano
The abovewritten system is obviously coupled to the Poisson equation,
which,in this case,becomes
E = −∇φ,∇(∇φ) = −e(N
+
−N
−
−n
Γ
−4n
L
) (56)
where Φ is the electric potential,N
+
and N
−
the donor and acceptor densities
respectively,and the dielectric constant.
These moment equations do not constitute a set of closed relations be
cause of the additional ﬂuxes and production terms.Therefore constitutive
assumptions must be prescribed.
If we assume as fundamental variables n
A
,V
i
A
,W
A
and S
i
A
,A = Γ,L,
which have a direct physical interpretation,the closure problem consists of
expressing U
ij
A
,F
ij
A
and G
ij
A
and the moments of the collision terms C
n
A
,C
P
i
A
,
C
W
A
and C
S
i
A
as functions of these variables.
We again resort to the MaximumEntropy Principle in order to obtain these
constitutive relations.In this case,one has to ﬁnd the distribution functions
which maximize the electron entropy
s [f
Γ
,f
L
] = −k
B
3
(f
Γ
logf
Γ
−f
Γ
) dk
Γ
+4
3
(f
L
logf
L
−f
L
) dk
L
.
After that,the inversion of the constraint relations is realized assuming,
as in the case of silicon,that the anisotropy of the distribution functions is
small.In this way explicit constitutive equations have been obtained for ﬂuxes
and production terms up to the ﬁrst order in δ.Here we do not report the
constitutive relations.The interested reader is referred to [5].
One can prove that the system(52)(55),closed with the maximumentropy
principle,is hyperbolic in the physically relevant region of the dependent
variables [6].
For the numerical integration we use the extension of the scheme developed
in [31,32] for homogeneous hyperbolic systems,which has been adapted in
[26,27,28,29,30,35] for balance laws with (possibly stiﬀ) source terms,see
appendix.
The complete method is based on a secondorder splitting technique which
separately solves the system with the source put equal to zero (convection
step) and the one with the ﬂux vector put equal to zero (relaxation step).For
more details see [6].
6.3 Simulations in bulk GaAs.
In this section we test the model in the case of a uniformly doped GaAs semi
conductor.Two diﬀerent impurity concentrations,N
+
=10
14
cm
−3
,10
17
cm
−3
are considered,and all the relevant scattering mechanisms are taken into ac
count.
Since the problem is homogeneous,in the evolution equations we can drop
the spatial dependence and the balance equations reduce to the following set
of ordinary diﬀerential equations
Recent Developments in Hydrodynamical Modeling of Semiconductors 39
d
dt
n
Γ
= n
L
C
(np)+
n ΓL
(W
L
) −n
Γ
C
(np)−
n ΓL
(W
Γ
),(57)
m
∗
Γ
d
dt
n
Γ
(V
i
Γ
+2α
Γ
S
i
Γ
) = −en
Γ
E
i
+c
Γ
11
n
Γ
V
i
Γ
+c
Γ
12
n
Γ
S
i
Γ
,(58)
d
dt
n
Γ
W
Γ
= −en
Γ
V
k
Γ
E
k
+n
L
C
(np)+
W ΓL
(W
L
) +
n
Γ
C
(p)
W
Γ
(W
Γ
) −C
(np)−
W ΓL
(W
Γ
)
,(59)
d
dt
n
Γ
S
i
Γ
= −e n
Γ
E
k
G
ik
Γ
+c
Γ
21
n
Γ
V
i
Γ
+c
Γ
22
n
Γ
S
i
Γ
,(60)
d
dt
n
L
= −n
L
C
(np)−
n LΓ
(W
L
) +n
Γ
C
(np)+
n LΓ
(W
Γ
),(61)
m
∗
L
d
dt
n
L
(V
i
L
+2α
L
S
i
L
) = −en
L
E
i
+c
L
11
n
L
V
i
L
+c
L
12
n
L
S
i
L
,(62)
d
dt
n
L
W
L
= −en
L
V
k
L
E
k
+n
Γ
C
(np)+
W LΓ
(W
Γ
) +n
L
C
(p)
W
L
(W
L
) +C
(np)
W LL
(W
L
)
−C
(np)−
W LΓ
(W
L
)
,(63)
d
dt
n
L
S
i
L
= −e n
L
E
k
G
ik
L
+c
L
21
n
L
V
i
L
+c
L
22
n
L
S
i
L
,(64)
where the meaning of the various productions terms which appear in the
equations can be found in [5].
From (57),(61) and the expressions of the production terms,one has that
n = n
Γ
+4n
L
= const,so that the total electron number is conserved,as it
must be.In cases when a constant bias voltage is applied to the semiconductor,
the Poisson equation is satisﬁed with n equal to the value of the doping
concentration and Econstant.The motion is along the direction of the electric
ﬁeld and,if we take this as the xdirection,the system (57)(64) reads
40 A.M.Anile,G.Mascali and V.Romano
d
dt
n
Γ
= n
L
C
(np)+
n ΓL
(W
L
) −n
Γ
C
(np)−
n ΓL
(W
Γ
),(65)
d
dt
n
Γ
V
Γ
=
2α
Γ
G
Γ
−
1
m
∗
Γ
e En
Γ
+
c
Γ
11
m
∗
Γ
−2 α
Γ
c
Γ
21
n
Γ
V
Γ
+
+
c
Γ
12
m
∗
Γ
−2 α
Γ
c
Γ
22
n
Γ
S
Γ
,(66)
d
dt
n
Γ
W
Γ
= −en
Γ
V
Γ
E +n
L
C
(np)+
W ΓL
(W
L
)
+n
Γ
C
(p)
W
Γ
(W
Γ
) −C
(np)−
W ΓL
(W
Γ
)
,(67)
d
dt
n
Γ
S
Γ
= −e n
Γ
EG
Γ
+c
Γ
21
n
Γ
V
Γ
+c
Γ
22
n
Γ
S
Γ
,(68)
d
dt
(n
Γ
+4n
L
) = 0,(69)
d
dt
n
L
V
L
=
2α
L
G
L
−
1
m
∗
L
e En
L
+
c
L
11
m
∗
L
−2 α
L
c
L
21
n
L
V
L
+
+
c
L
12
m
∗
L
−2 α
L
c
L
22
n
L
S
L
,(70)
d
dt
n
L
W
L
= −en
L
V
L
E +n
Γ
C
(np)+
W LΓ
(W
Γ
) +n
L
C
(p)
W
L
(W
L
) +C
(np)
W LL
(W
L
)
−C
(np)−
W LΓ
(W
L
)
,(71)
d
dt
n
L
S
L
= −e n
L
EG
L
+c
L
21
n
L
V
L
+c
L
22
n
L
S
L
,(72)
where V
A
and S
A
are the xcomponents of V
A
and S
A
and G
A
is the xx
component of G
ij
A
,A = Γ,L.In the evolution equations the xcomponent E
of the electric ﬁeld enters as a parameter.
As initial conditions for (65)(72),we take,in suitable units,
n
Γ
(0) +4n
L
(0) = 1,
n
Γ
(0)
n
L
(0)
=
R
3
f
(eq)
Γ
dk
Γ
R
3
f
(eq)
L
dk
L
=
m
∗
Γ
m
∗
L
3/2
d
Γ
0
(
1
k
B
T
L
)
d
L
0
(
1
k
B
T
L
)
,
V
A
(0) = 0,
W
A
(0) = g
A
1
k
B
T
L
,
S
A
(0) = 0,A = Γ,L.
The crystal temperature,T
L
,is assumed to be 300 K.The expressions of d
Γ
0
,
d
L
0
and g
A
are given in [5].
The solutions of (65)(72) for electric ﬁelds respectively equal to 0.2,0.5,1,
2 and 6V/µm,are reported in Fig.18.
Recent Developments in Hydrodynamical Modeling of Semiconductors 41
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
1
2
3
4
5
6
7
8
time(ps)
velocity (107 cm/s)
6 V/mum
2 V/mum
1 V/mum
0.5 V/mum
0.2 V/mum
N
+
=10
14
/cm
3
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
5
0
1
2
3
4
5
6
7
8
time(ps)
velocity (107 cm/s)
6 V/µm
2 V/µm
1 V/µm
0.5 V/µm
0.2 V/µm
N
+
=10
17
/cm
3
Fig.18.The time evolution of the electron average velocity for diﬀerent values of
the electric ﬁeld and for N
+
= 10
14
/cm
3
and 10
17
/cm
3
respectively.
The stationary regime is reached in a few picoseconds.
The typical phenomena of overshoot and saturation of the velocity are
both qualitatively and quantitatively well described.
We also report the curves representing the electron valley occupancy,
n
Γ
n
Γ
+4n
L
and
4n
L
n
Γ
+4n
L
,and average velocity as functions of the electric ﬁeld
(see Figs.19 and 20) for the above impurity concentrations.
42 A.M.Anile,G.Mascali and V.Romano
0
1
2
3
4
5
6
7
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
electric field (V/µm)
occupancy
L−valleys
Γ−valley
N
+
=10
14
/cm
3
0
1
2
3
4
5
6
7
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
electric field (V/µm)
occupancy
L−valleys
Γ−valley
N
+
=10
17
/cm
3
Fig.19.Electron occupancy in the Γ and Lvalleys for N
+
= 10
14
/cm
3
and
10
17
/cm
3
respectively
Recent Developments in Hydrodynamical Modeling of Semiconductors 43
0
1
2
3
4
5
6
7
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
electric field (V/µm)
velocity (10
7cm/s)
N
+
=10
14
/cm
3
N
+
=10
17
/cm
3
Fig.20.Electron velocity versus electric ﬁeld characteristics for N
+
= 10
14
/cm
3
and 10
17
/cm
3
respectively
Fig.19 clearly shows the electron transfer from the Γvalley to the
Lvalleys.The population inversion is observed at electric ﬁelds of about
1.8V/µm.
Fig.20 shows that the low ﬁeld mobility and the electron peak velocity
decrease with the increase of the impurity concentration.
6.4 Simulation a GaAs n
+
−n −n
+
diode
In this section we consider the case of a n
+
−n−n
+
diode which models the
channel of a MOSFET.The device is made of GaAs and its temperature is
300 K.The n
+
regions are 0.1 µm long,while the channel length is 0.4 µm.
The doping proﬁle is
N
+
=
2 ×10
18
cm
−3
in the n
+
regions,
10
16
cm
−3
in the n region,
An external voltage of 2 V is applied.To avoid initiating too much a com
plicated transient behaviour in the device,the voltage is slowly raised to this
value.The following initial and boundary conditions are considered
n
Γ
(x,0) =
r
(0)
4 +r
(0)
N
+
(x),n
L
(x,0) =
N
+
(x) −n
Γ
(x,0)
4
,
W
A
(x) = W
(0)
A
,V
A
(x,0) = S
A
(x,0) = 0,
44 A.M.Anile,G.Mascali and V.Romano
n
Γ
(a,t) =
r
(0)
4 +r
(0)
N
+
(a),n
L
(a) =
N
+
(a) −n
Γ
(a)
4
,
W
A
(a) = W
(0)
A
,
∂V
A
∂x

a
=
∂S
A
∂x

a
= 0,A = Γ,L,
a = 0,0.6 µm.
where r
(0)
is the equilibrium ratio between n
Γ
and n
L
,and W
(0)
Γ
and W
(0)
L
are the equilibrium energies at the lattice temperature.Boundary conditions
for Poisson’s equation are imposed by specifying the electric potential at the
device contacts:
φ(0) = 0V φ(0.6) = 2V.
Simulations are run to a ﬁnal time of t = 75ps,at which the solutions are
judged to reach a steady state.The results relative to the electric ﬁeld,the
total electron average velocity and energy are represented in Figs 21,22,where
a comparison with analogous results for a Si n
+
−n−n
+
diode is also shown.
The great diﬀerence in the velocity,which can be seen in Fig.21b is essentially
due to the diﬀerence between the eﬀective masses of Silicon electrons and Γ
electrons in GaAs.Indeed near the ﬁrst junction the mean energy is low and
electrons are mainly in the Γvalley where they have a small eﬀective mass
and consequently high velocities.Instead close to the second junction,where
the electric ﬁeld reaches its maximumvalue,the Lvalley is almost exclusively
populated.Since here the eﬀective mass is much greater than that in the Γ
valley and comparable with the eﬀective mass in silicon,the mean velocity is
lower.
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
−20
−15
−10
−5
0
5
10
micron
electric field (Volt/ µm)
____ GaAs
****** Si
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
−0.5
0
0.5
1
1.5
2
2.5
3
3.5
4
4.5
micron
mean velocity (10
7 cm/sec)
____ GaAs
****** Si
Fig.21.Acomparison between a Si diode and a GaAs one,electric ﬁeld and velocity
vs position.
Recent Developments in Hydrodynamical Modeling of Semiconductors 45
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
micron
mean energy (eV)
____ GaAs
***** Si
Fig.22.A comparison between a Si diode and a GaAs one,energy vs position.
6.5 Gunn oscillations
In this section we consider a GaAs diode with the same doping proﬁle as that
in [46] and [47].This is also coupled to an RLC tank circuit which stimulates
Gunn oscillatory eﬀects.The onedimensional diode has length L
d
= 2µm
and its doping proﬁle is
N
+
(x) =
10
17
for x < 0.125 µm,
10
16
for 0.125 µm < x < 0.15 µm,
0.5 ×10
16
for 0.15 µm < x < 0.1875 µm,(donors/cm
3
)
10
16
for 0.1875 µm < x < 1.875 µm,
10
17
for 1.875 µm < x.
(73)
The transitions in the doping proﬁle at the device junctions are discontinuous
(in contrast to [46]).The same initial and boundary conditions as before have
been used except for the potential.Now φ(L
d
) may be either equal to 2V or
is determined by coupling the device to a system of ODE which models the
circuit.These equations read
dV
d
dt
=
1
C
I −I
d
−
V
d
R
,
dI
dt
=
1
Λ
(V
b
−V
d
),(74)
where V
d
is the voltage through the device,V
b
= 2V the bias voltage of the
circuit,and I
d
,the particle current in the device,is calculated as
I
d
= −
q A
L
d
L
d
0
(n
Γ
v
Γ
+4n
L
v
L
) dx.(75)
46 A.M.Anile,G.Mascali and V.Romano
Simple ﬁnite diﬀerence versions of equations (74) and (75) allow the diode
voltage to be updated at each simulation time step.The values used for the
capacitance,C,resistance,R,and inductance,Λ,of the circuit are
C =
A/L
d
+0.82 ×10
−12
F,R = 25 ohm,Λ = 3.5 ×10
−12
henry,
where the crosssectional area,A,of the diode is equal to 1.0 ×10
−3
cm
2
.
The oscillator equations (74) are given the initial state
V
d
(t
0
) = 2V,I(t
0
) = 0.(76)
The circuit is engaged at time t
0
= 75ps when the Ga As diode is judged to
have reached the steady state illustrated in Figs 23,24,25.
Fig.23.Eletron density and velocity vs position.Fig.7a:Continuus line:Γelectron
density,dotted line:Lelectron density,dotteddashed line:total electron density
Recent Developments in Hydrodynamical Modeling of Semiconductors 47
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
0
0.05
0.1
0.15
0.2
0.25
0.3
0.35
micron
mean energy (eV)
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
0.15
0.155
0.16
0.165
0.17
0.175
0.18
0.185
micron
current (1024 cm−2 sec−1)
Fig.24.Eletron energy and current vs position.
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
1.8
2
0.15
0.155
0.16
0.165
0.17
0.175
0.18
0.185
micron
current (1024 cm−2 sec−1)
Fig.25.Potential vs position.
48 A.M.Anile,G.Mascali and V.Romano
One observes,Fig.26,that there are some initial oscillations that smooth
out and become negligible after about 200 ps.The qualitative behaviour,at
variance with other hydrodynamical models,as those presented in [47],is very
similar to the MC simulation reported in [47],even though this latter has been
obtained in the parabolic band approximation.
0
50
100
150
200
250
300
350
1.85
1.9
1.95
2
2.05
2.1
2.15
2.2
time (ps)
potential (Volt)
Gunn diode
Fig.26.The potential V
d
versus time for the Gunn diode.
Acknowledgments
This work has been partially supported by MURST,ex fondi 60%,and
by CNR grants n.98.01041.CT01,n.99.01714.01,n.00.00128.ST74 and n.
CNRG00DB7 (program Agenzia 2000).
Appendix:numerical method
The left hand side of the moment equations represents a quasilinear hyperbolic
operator,while the right hand side contains relaxation and drift terms.We
make use of a splitting scheme,based on the following decomposition.Without
loss of geenrality,let us consider a scalar one space dimensional system of the
form
∂u
∂t
+
Enter the password to open this PDF file:
File name:

File size:

Title:

Author:

Subject:

Keywords:

Creation Date:

Modification Date:

Creator:

PDF Producer:

PDF Version:

Page Count:

Preparing document for printing…
0%
Comments 0
Log in to post a comment