Lecture Notes in Mathematics

woundcallousΗμιαγωγοί

1 Νοε 2013 (πριν από 5 χρόνια και 9 μέρες)

483 εμφανίσεις

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 15-22,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
e-mail:anile@dmi.unict.it
Walter Allegretto
Department of Mathematical and
Statistical Sciences
Alberta University
Edmonton AB T6G 2G1
Canada
e-mail:retl@retl.math.ualberta.ca
Christian Ringhofer
Department of Mathematics
Arizona State University
Tempe,Arizona 85287-1804,USA
e-mail:ringhofer@asu.edu
Cataloging-in-Publication 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
ISSN0075-8434
ISBN3-540-40802-9 Springer-Verlag 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 fromSpringer-Verlag.Violations are
liable for prosecution under the German Copyright Law.
Springer-Verlag Berlin Heidelberg NewYork a member of BertelsmannSpringer
Science + Business Media GmbH
springer.de
c
￿
Springer-Verlag 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:Camera-ready T
E
X output by the authors
SPIN:10952481 41/3142/du - 543210 - Printed on acid-free 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-
diffusion 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
Wigner-Boltzmann-Poisson 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 difficult case).
Another problem of fundamental interest is that of the hydrodynamical limit
which one expects to be quite different from the Navier-Stokes-Fourier one,
since the collision operator is substantially different from the one in rarefied
gas case.
From the application viewpoint it is of great practical importance to study
efficient 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 post-docs
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 15-22 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 23-24 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 drift-diffusion 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:drift-diffusion 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 rarefied 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 (co-director
of the course) on the asymptotic limit for strong fieds,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 Wigner-Poisson system,and of
A.Marrocco on the mixed finite element discretization of the energy transport
model.
A.M.Anile
CIME’s activity is supported by:
Ministero dell’Universit`a Ricerca Scientifica 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 H-Theorem 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
Drift-Diffusion Equations and Applications
W.Allegretto....................................................57
1 The Classical Semiconductor Drift-Diffusion System...............57
1.1 Derivation................................................57
1.2 Existence.................................................58
1.3 Uniqueness and Asymptotics................................63
2 Other Drift-Diffusion 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 Effects.......................................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 Multi-Body Equations and Effective Single Electron Models........98
1.1 Effective 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 Diffusion Approximations to Kinetic Equations....................111
3.1 Diffusion Limits:The Hilbert Expansion.....................113
3.2 The Drift Diffusion 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 definition 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-field phenomena such as hot electron propagation,impact
ionization and heat generation.In fact the standard drift-diffusion models
cannot cope with high-field 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 flux 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-
diffusion model,is not restricted to quasi-stationary situations.
These models are,loosely speaking,called hydrodynamical models and
they are usually derived by suitable truncation procedures,from the infi-
nite hierarchy of the moment equations of the Boltzmann transport equation.
However,most of these suffer from serious theoretical drawbacks due to the
A.M.Anile,W.Allegretto,C.Ringhofer:LNM 1821,A.M.Anile (Ed.),pp.1–56,2003.
c￿Springer-Verlag 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 L-valleys have been considered.Therefore electrons in
the conduction band have been treated as a mixture of two fluids,one repre-
senting electrons in the Γ-valley and the other electrons in the four equivalent
L-valleys.
Both in the Si and in the GaAs case,the models comprise the balance
equations of electron density,energy density,velocity and energy flux.The
only difference 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 field 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 quasi-particles).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 simplified 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 first Brilloiun zone,
near the X points,which,for this reason,are termed the X-valleys.In GaAs
there is an absolute minimum at the center of the Brillouin zone,the Γ-point,
and local minima at the L-points along the Λ cristallographic orientations.
As mentioned above,in the simplified 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 effective electron mass in the A-valley and  the reduced Planck
constant.
A more appropriate analytical approximation,which takes into account
the non-parabolicity 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
effects 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 field,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 Boltzmann-Poisson 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 electron-electron scat-
terings and the interaction with stationary imperfections of the crystal as
vacancies,external and internal crystal boundaries.In many situations the
electron-electron collision term can be neglected since the electron density is
not too high.However in the case of high doping,electron-electron collisions
must be taken into account because they might produce sizable effects.Re-
taining the electron-electron 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 electron-electron 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 first 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 Fermi-Dirac
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 H-theorem and the null space of the collision operator
In [9,10,11] an H-theorem has been derived for the physical electron-
phonon operator in the homogeneous case without electric field.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 field,neglecting the electron-electron
interaction and assuming the electron gas sufficiently dilute to neglect the
degeneracy effects.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 so-called 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 Bose-Einstein 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 definition 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 Boltzmann-Poisson system
(5)-(7).The first 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 finding 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 electron-phonon
operator was tackled and solved in general in [11] where it was proved that
the equilibrium solutions are not only the Fermi-Dirac distributions but form
an infinite 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 differs 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 Fermi-Dirac 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 fluxes;
ii) the closure for the production terms.
As in gasdynamics [15],multiplying eq.(5) by a sufficiently regular func-
tion ψ(k) and integrating over B,the first 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 field 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


=

ψ(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 sufficiently 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 different 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 sufficiently 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
defined,does not arise here because the moments are obtained by integrating
over the first 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
finds the closure relations for the fluxes and the production terms appearing
in the balance equations.On account of the analytical difficulties 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 differential 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 first 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 flux
10 A.M.Anile,G.Mascali and V.Romano
and
g =

η

(f)C[f]dk entropy production.
The electric field contributes neither to the entropy production nor to the
entropy flux.
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
=


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


+

ψ
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 field 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 defined 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 defined 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 well-posed for the system (24),at least in the simple case
where the electric field is considered as an external field.
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);
• electron-phonon 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

effective 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
acoustic-phonon 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 Bose-Einstein 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

2
ρv
2
s
,K
α
=
Z
f α
(D
t
K)
2
α

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 final 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
first 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 fitting parameters.Their values depend on the approxima-
tion used for the energy bands,on the specific 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 flux,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 flux,
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 flux 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 flux of energy flux,
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 flux,
These moment equations do not constitute a set of closed relations because
of the fluxes 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 fluxes and
in [3] for the production terms up to second order in δ.However it has been
found in [26] that the first order model is sufficiently 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 sufficiently 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 finds 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 field which is chosen as x-direction
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 x-component 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 field are
reported in figs.1 (velocity),2 (energy) and 3 (energy flux).
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 field
as in figure 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 flux (eV cm/sec) versus time (ps) for the same values of the electric
field as in figure 1
The typical phenomena of overshoot and saturation velocity are both qual-
itatively and quantitatively well described (see [33] fig.3.22 for a comparison
with results obtained by MC simulations).
Similar results were reported in [3],but there a different 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,figs.4,5,6.The differences,especially in the energy,with respect to
the Kane case,confirm that the parabolic band is an oversimplification 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 field (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 flux (eV cm/sec) versus electric field (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 profiles 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-
flux.
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 profile 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 five),after a short transient with
wide oscillations.
As first 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 ( fig.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
(fig.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 flux S is not a field 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 justified [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 five).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 fields are
involved,as shown in the previous section for the n
+
−n −n
+
diode.
The density is plotted in figure 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 (figure 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 figure 14.The higher values of
the x-component are at the edges of the gate contact.This happens also for
the y-component,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 flux (figure 15) is qualitatively similar to that of
the velocity.
Very large tangential and normal components of the electric field (figure
16) are present again at the edges of the gate.For completeness the electric
potential is also presented in fig.17.
The results are qualitatively similar to those presented in [40,41] for all
the variables except the y-component 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 first 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 non-polar optical phonon scattering,
• the polar optical phonon scattering,
• the impurity scattering.
All these scattering mechanisms are intravalley,except the non-polar 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 non-polar 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 non-polar optical phonon Bose-Einstein dis-
tribution.
Z
AB
is the degeneracy of the final valley (Table 4),that is the valley
the electron reaches after the scattering,with respect to the initial one.In
GaAs,non-polar optical phonons contribute to the electron intervalley transfer
between two equivalent L-valleys and between Γ-valley and L-valleys.
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

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 electron-electron 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

Γ
effective electron mass in the Γ-valley
0.067 ×m
e
m

L
effective electron mass in the L-valley
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 L-valley
0.242 eV
−1

r
relative dielectric constant
12.90


relative dielectric constant at optical frequency range
10.92
Ξ
d
acoustic-phonon deformation potential
7 eV
D
t
K
non-polar optical phonon deformation potential
10
9
eV/cm
￿ω
np
non-polar optical phonon energy
0.03 eV
￿ω
p
polar optical phonon energy
0.03536 eV
E
Γ(0)
Γ-valley bottom energy
0 eV
E
L(0)
L-valley bottom energy
0.32 eV
Z
ΓL
degeneracy from Γ to L valleys
4
Z

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 L-valley.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-fluxes of both Γ and L-valley 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 flux,
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 flux,
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 flux of energy flux,
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 flux production.
All these quantities refer to electrons in the A-valley,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 above-written 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 fluxes 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 find 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 fluxes
and production terms up to the first 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 stiff) source terms,see
appendix.
The complete method is based on a second-order splitting technique which
separately solves the system with the source put equal to zero (convection
step) and the one with the flux 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 different 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 differential 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 satisfied with n equal to the value of the doping
concentration and Econstant.The motion is along the direction of the electric
field and,if we take this as the x-direction,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
Γ
=


Γ
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
=


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 x-components of V
A
and S
A
and G
A
is the xx-
component of G
ij
A
,A = Γ,L.In the evolution equations the x-component E
of the electric field 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 fields 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 different values of
the electric field 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 field
(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 L-valleys 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 field 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
L-valleys.The population inversion is observed at electric fields of about
1.8V/µm.
Fig.20 shows that the low field 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 profile 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 final time of t = 75ps,at which the solutions are
judged to reach a steady state.The results relative to the electric field,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 difference in the velocity,which can be seen in Fig.21-b is essentially
due to the difference between the effective masses of Silicon electrons and Γ-
electrons in GaAs.Indeed near the first junction the mean energy is low and
electrons are mainly in the Γ-valley where they have a small effective mass
and consequently high velocities.Instead close to the second junction,where
the electric field reaches its maximumvalue,the L-valley is almost exclusively
populated.Since here the effective mass is much greater than that in the Γ-
valley and comparable with the effective 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 field 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 profile as that
in [46] and [47].This is also coupled to an RLC tank circuit which stimulates
Gunn oscillatory effects.The one-dimensional diode has length L
d
= 2µm
and its doping profile 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 profile 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 finite difference 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 cross-sectional 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:L-electron density,dotted-dashed 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
+