NUMERICAL METHODS FOR THE SEMICONDUCTOR BOLTZMANN EQUATION BASED ON SPHERICAL HARMONICS EXPANSIONS AND ENTROPY DISCRETIZATIONS

bentgalaxyΗμιαγωγοί

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

161 εμφανίσεις

NUMERICAL METHODS FOR THE
SEMICONDUCTOR BOLTZMANN
EQUATION BASED ON SPHERICAL
HARMONICS EXPANSIONS AND
ENTROPY DISCRETIZATIONS
Christian Ringhofer
Department of Mathematics,
Arizona State University,Tempe,AZ 85287-1804
E-mail:ringhofer@asu.edu
ABSTRACT
We present a class of numerical methods for the semiconductor
Boltzmann Poisson problem in the case of spherical band
energies.The methods are based on spherical harmonics
expansions in the wave vector and difference discretizations
in space–time.The resulting class of approximate solutions
dissipate a certain type of entropy.
Key Words:Boltzmann equation;Galerkin methods;Finite
differences
AMS Classification:65N35;65N05
431
DOI:10.1081/TT-120015508 0041-1450 (Print);1532-2424 (Online)
Copyright & 2002 by Marcel Dekker,Inc.www.dekker.com
TRANSPORT THEORY AND STATISTICAL PHYSICS
Vol.31,Nos.4–6,pp.431–452,2002
+ [28.8.2002–7:05am] [431–452] [Page No.431] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
1.INTRODUCTION
This paper is concerned with the space–time discretization of a class of
conservation laws which arise fromthe approximation of the semiconductor
Boltzmann–Poisson problem via an spherical harmonics expansion in the
wave vector direction.We consider the Boltzmann–Poisson problem
ðaÞ @
t
f þLð f Þ þ
1

Qð f Þ ¼ 0
ðbÞ Lð f Þ:¼ r
x
 ½r
k
"ðkÞ f  qr
k
 ½r
x
VðxÞ f 
ðcÞ 
x
V þq½D
dop
ðxÞ  ¼ 0,ðx,tÞ:¼
Z
f ðx,k,tÞ dk ð1Þ
where the kinetic density function f ðx,k,tÞ depends on the spatial variable
x 2 R
d
with d¼1,2,or 3,the wave vector k 2 R
3
and the time t.The
acceleration of particles is given by the electric field qr
x
V which,using
a mean field approximation,is determined by the Poisson Eq.(1c).Here,
q denotes the charge of one particle,so q < 0 holds for electrons and q > 0
holds if the Boltzmann equation describes the transport of holes.The
function D
dop
denotes a charge background,produced by the implantation
of different species atoms into the crystal. is the dielectricity constant of
the material.The dimensionless constant  in Eq.(1a) denotes the scaled
mean free path (the Knudsen number).Since,we are modeling transport in
a crystal,as opposed to a vacuum,the velocity of a particle is given in terms
of its wave vector through of a dispersion relation of the form
vðkÞ ¼ r
k
"ðkÞ
where"ðkÞ denotes the band energy.In vacuum"ðkÞ ¼ constjkj
2
=2 holds and
the velocity vector and the wave vector are,up to a multiplicative constant,
identical.For the purposes of this paper,we will take the band energy
function to be spherically symmetric,so we will assume that"ðkÞ ¼"ðjkjÞ
holds.One of the most popular approximate choices for the band energy is
the Kane dispersion model of the form
"ðkÞ ¼
jkj
2
1 þ
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
1 þ jkj
2
p
,
but the approach outlined in this paper is certainly applicable to more
complicated band structures.The collision operator Q in Eq.(1a) models
the interaction of particles with a uniform phonon background,and can
+ [28.8.2002–7:05am] [431–452] [Page No.432] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
432 RINGHOFER
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
therefore be taken to be linear.In the scaled and dimensionless form (1),
the collision operator is of the form
Qð f Þðx,k,tÞ ¼
X
1

¼ 1
c


f ðx,k,tÞ
Z
R
3
k
dk
0
ð"ðkÞ "ðk
0
Þ þ

 
(

Z
R
3
k
dk
0
ð"ðk
0
Þ "ðkÞ þ
!Þ f ðk
0
Þ
 
)
:ð2Þ
So a particle can gain or loose an amount!of energy,with a probability
c
1
,c
1
,in the collision event.Because of the principle of detailed balance,
the constants c


have to be such that local Maxwellians of the form
f ðx,k,tÞ ¼ ðx,tÞe
"ðkÞ
are in the kernel of the collision operator.Thus the
constants c


has to satisfy the relation c


¼ c


e

!
.To express the
symmetry properties and integral invariants of the phonon collision opera-
tor it will be more convenient to rewrite Q as
Qð f Þ ¼
Z
R
3
k
dk
0
S","
0


fe
"
f
0
e
"
0
h i
,ð3Þ
where the scattering cross section Sð","
0
Þ is a symmetric function,so
Sð","
0
Þ ¼ Sð"
0
,"Þ holds.We use the short hand f
0
for f ðx,k
0
,tÞ in Eq.(3)
and from now on.S is then of the form
ðaÞ Sð","
0
Þ ¼
X
1

¼ 1
s


ð","
0
Þ ð" "
0
þ

ðbÞ s


ð","
0
Þ ¼ c


exp

! " "
0
2
 
:ð4Þ
The standard way to solve the Boltzmann Eq.(1) numerically is of
course the Monte Carlo (MC) method.This paper is concerned with an
alternative,deterministic approach,namely to expand the kinetic density
function f into a series of basis functions in the wave vector direction,
obtaining a system of hyperbolic conservation laws with source terms for
the coefficients.While such an approach cannot hope to reproduce MC
results in all physical detail,it has the advantage of being deterministic
and computationally considerably less expensive,particularly in fewer
space dimensions (d ¼1 or 2),since it has the ability to take advantage of
symmetries.Usually these approaches,and the closely related moment
closure methods,have to use sometimes quite crude approximations to
the collision operators to be effective (see Ref.[3]).Alternatively,they
+ [28.8.2002–7:05am] [431–452] [Page No.433] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
SEMICONDUCTOR BOLTZMANN EQUATION 433
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
have to use a difference method in the radial component of the wave vector
to deal with the functions in the integral kernel (see Refs.[2,5,12]).
This limits the use of these methods in more than one spatial dimension
because of the large number of necessary unknowns.In this paper,we present
a complete spectral approximation to the Boltzmann equation in the wave
vector direction using the complete phonon collision operator Q as given
in Eq.(3).The resulting system of hyperbolic conservation laws can be
extremely stiff in the spatial direction,due to very large potential gradients
and in the time direction due to small values of the Knudsen number .
In Sec.3 of this paper,we present several approaches to deal with these
stiffness problems.In Sec.4,we present a simple numerical test example
which serves to demonstrate that the resulting method is capable to simulate
transport quite far away from the equilibrium regime.
2.THE MOMENTUM DISCRETIZATION
The general idea of the approach outlined in this paper is to discretize
the Boltzmann Eq.(1) in the wave vector direction by a spectral Galerkin
method,thus obtaining a first order system of hyperbolic conservation laws
in position and time,which will contain source terms due to the acceleration
of particles by the electric field and the collision operator Q in Eq.(1).
This system will then be discretized by methods suitable for hyperbolic
conservation laws,which are described in the next section.An important
concept which is used in making the choice of basis functions is that of an
entropy,i.e.,a function hð f,x,kÞ which satisfies
ðaÞ
Z
R
3
k
dk
Z
R
d
x
dx hð f,x,kÞLð f Þ½  ¼ 0,
ðbÞ
Z
R
3
k
dk
Z
R
d
x
dx hð f,x,kÞQð f Þ½  0:ð5Þ
Astraight forward calculation yields,that for any function h satisfying
Eq.(5)
Z
R
d
x
dx
Z
R
3
k
dk½hð f,x,kÞ@
t
f  ¼
d
dt
Z
R
d
x
dx
Z
R
3
k
dkHð f,x,kÞ  0 ð6Þ
holds,where H denotes the antiderivative of h with respect to f.If the
function H,the entropy,is a convex function of the variable f this guaran-
tees that the Boltzmann Eq.(1) is well posed,and the basis functions for the
+ [28.8.2002–7:05am] [431–452] [Page No.434] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
434 RINGHOFER
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
wave vector discretization will be chosen such that the property (6) is pre-
served for the semi-discrete system.To this end we rewrite the Boltzmann
Eq.(1) in the variable gðx,k,tÞ ¼ hð f ðx,k,tÞ,x,kÞ,obtaining
@
t
ðg,x,kÞ þLððg,x,kÞÞ þ
1

Qððg,x,kÞÞ ¼ 0,ð7Þ
where  denotes the functional inverse of h;so g ¼ hð f,x,kÞ ()f ¼
ðg,x,kÞ holds.Equation (7) is now discretized by a standard Galerkin
procedure in the k-direction.So we express g as a linear combination of
basis functions in k,whose coefficients depend on x and t,and integrate
Eq.(7) against each basis function with respect to k.This gives the hyper-
bolic system of conservation laws for the expansion coefficients,which then
automatically satisfies the entropy estimate (6),and is thus well posed.
The physical entropy is given by logarithmic function (hð f,x,kÞ ¼ lnf,
Hð f,x,kÞ ¼ f ðln f 1Þ).However,since we are using only the linearized
collision operator,modeling the collision of electrons with phonons,
there is a much wider choice of possible entropy functions.As a matter of
fact,the equality (see Ref.[6])
Z
R
3
k
dk½hQð f Þ ¼
1
2
Z
R
3
k
dk
Z
R
3
k
dk
0
S fe
"
f
0
e
"
0

ðh h
0
Þ
h i
ð8Þ
implies that any function hð f,x,kÞ which is a monotonically increasing
function of fe
"ðkÞ
produces an entropy.The more complicated the chosen
entropy function,the more difficult it will be to compute the discretization
of the collision operator Q.Galerkin expansions which use the physical
entropy (the logarithm) usually have to replace the collision operator by
some BGK-type approximation to yield practical schemes (see Ref.[3]).
In the approach outlined in this paper we will instead choose the simplest
possible entropy,namely the linear function
hð f,x,kÞ ¼ e
"þqV
f,Hð f,x,kÞ ¼
1
2
e
"þqV
f
2
,ðg,x,kÞ ¼ e
" qV
g:
Note that this choice of entropy seems to implicitly imply a time
independent potential V.In the case,considered in this paper,when the
potential is computed self consistently via the Poisson equation,the above
choice of entropy can still be used.However the estimate (6),together with
the proof that the resulting entropy function H yields a convex functional,
becomes quite a bit more complicated.We refer the reader to Ref.[8] for the
details.
+ [28.8.2002–7:05am] [431–452] [Page No.435] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
SEMICONDUCTOR BOLTZMANN EQUATION 435
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
This leaves us with the choice of basis functions.The structure of the
collision operator (i.e.,the averaging over surfaces of equal energy) suggests
to rewrite the Boltzmann Eq.(1) in polar coordinates and to discretize
by spherical harmonic functions in the k-direction.The idea to expand
solutions of the Boltzmann equation into spherical harmonic functions goes
back a long time (see Ref.[1] for an overview).This set of basis functions
was first used numerically and in the context of semiconductors in Ref.[12]
and later on by Ref.[2],and more recently in Refs.[4,5].All of these
approaches use a difference approximation of the density function in the
radial direction.Because of the number of gridpoints necessary in the radial
direction to resolve the collision operator,they are,given current computer
resources,restricted to the case of one space dimension (d ¼1).In this paper,
we will instead use a spectral method in the radial direction as well,which
enables us to dramatically reduce the number of variables.According to the
entropy principles used,it is natural to use a polynomial basis in the radial
direction for the entropy variable g,and correspondingly,polynomials time
the equilibrium Maxwellian for the density function f.Thus we set
gðx,k,tÞ ¼
X
m2M
g
m
ðx,tÞ
m
ðr,,Þ,

m
ðr,,Þ ¼ P
m
1
ðrÞ
m
2
m
3
ð,Þ,m ¼ ðm
1
,m
2
,m
3
Þ
where the triple index mvaries in some finite index set M.r,, denote polar
coordinates and  is the spherical harmonics given by

m
2
m
3
ð,Þ ¼ 
m
3
m
2
ðcos ÞðsinÞ
m
3
expðim
3
Þ,k ¼
r cos 
r sin sin
r sin cos 
0
@
1
A
:
Here 
m
3
m
2
denote the associated Legendre polynomials and P
m
1
are some
suitable basis polynomials in the radial direction.Correspondingly the
Boltzmann density function f is approximated by
f ðx,k,tÞ ¼ ðg,x,kÞ ¼
X
m2M
f
m
ðx,tÞ
m
ðr,,Þe
"ðrÞ
,ð9Þ
where we have absorbed the factor e
qV
into the coefficients f
m
for simplicity.
The expansion (9) is inserted into the Boltzmann equation and integrated
against the basis functions for g of the form P
m
1
ðrÞ
m
2
m
3
ð,Þ with respect
to the wave vector k.Thus the whole procedure reduces,because of the linear
choice for h,to a Galerkin approximation for the Boltzmann density f with a
scalar product given by the weight function e
"ðjkjÞ
.Correspondingly,the
polynomials P
m
1
in Eq.(9) should be chosen orthogonal with respect to
+ [28.8.2002–7:05am] [431–452] [Page No.436] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
436 RINGHOFER
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
this weight function,so as to diagonalize the mass matrix.So the P
m
1
will
satisfy
Z
1
0
P
m
1
ðrÞP
n
1
ðrÞe
"ðrÞ
r
2
dr ¼
m
1
n
1
,P
0
ðrÞ ¼ const
With this choice the resulting system of conservation laws becomes
ðaÞ @
t
F þ
X
d

¼1
½A


@
x


F qð@
x


VÞB


F þ
1

CF ¼ 0,F ¼ ð f
111
,...Þ
ðbÞ A


ðm,nÞ ¼
Z
R
3
k
dk½
m
ð@
k


"Þe
"

n
,
ðcÞ B


ðm,nÞ ¼
Z
R
3
k
dk
m
@
k


½e
"

n

 
,
ðdÞ Cðm,nÞ ¼
Z
R
3
k
dk
m
Q½e
"

n

 
,
ðeÞ
m
¼ P
m
1
ðrÞ
m
2
m
3
ð,Þ ð10Þ
Once the matrices A


,B


,and C are computed,Eq.(10) constitutes a
system of hyperbolic conservation laws whose discretization in space and
time will be the subject of the next section.In the case of non-parabolic
bands,i.e.,"ðkÞ 6
¼ ðjkj
2
=2Þ these coefficient matrices will have to be computed
numerically.A considerable computational effort can be spent on this com-
putation since it only has to be performed once for a given band structure
and a given set of basis functions.However,the numerical evaluation of
the coefficient matrices should be done in such a way that the entropy
property (6) translates exactly to the semi-discrete system (10).We will
therefore briefly outline how to do this.First,integrating the identity
@
k


ðe
"

m

n
Þ ¼
m
ð@
k


"Þe
"

n
þ
m
@
k


½e
"

n
 þ
n
@
k


½e
"

m

with respect to the wave vector k implies that
A


¼ B


B
T


,
¼ 1,...,d ð11Þ
holds.This relation should hold exactly and therefore,we compute the
matrices B


by some numerical integration rule and define the matrices A


via Eq.(11).Second,using the identity (8),we write the collision matrix C as
Cðm,nÞ ¼
1
2
Z
R
3
k
dk
Z
R
3
k
dk
0
½Sð
m

0
m
Þð
n

0
n
Þ,ð12Þ
+ [28.8.2002–7:05am] [431–452] [Page No.437] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
SEMICONDUCTOR BOLTZMANN EQUATION 437
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
which implies that the matrix C is symmetric and positive semi-definite.
Moreover,the first row and column of C vanish identically since

000
¼ const holds.These two properties should also hold exactly for the
numerical approximation of the collision matrix C.Inserting the spherical
harmonic basis functions into Eq.(12) gives
Cðm,nÞ ¼
1
2
Z
R
3
k
dk
Z
R
3
k
dk
0
½Sð
m

0
m
Þð
n

0
n
Þ
¼
1
2
Z
1
0
dr
Z

0
d
Z


d
Z
1
0
dr
0
Z

0
d
0
Z


d
0
½r
2
sinðÞðr
0
Þ
2
sinð
0
ÞSð","
0
ÞðP
m
1

m
2
m
3
P
0
m
1

0
m
2
m
3
ÞðP
n
1

n
2
n
3
P
0
n
1

0
n
2
n
3
Þ:
Since the band energy function"depends only on the modulus r of the
wave vector,this sixfold integral can be simplified considerably by integrating
out the angular variables.A direct calculation gives
Cðm,nÞ ¼ 2
m
2
0

m
3
0

n
2
0

n
3
0
Z
1
0
dr
Z
1
0
dr
0
r
2
ðr
0
Þ
2
Sð","
0
Þ P
m
1
P
0
m
1


 P
n
1
P
0
n
1


þ4
m
2
n
2

m
3
n
3

m
2
0

m
3
0

n
2
0

n
3
0



Z
1
0
dr
Z
1
0
dr
0
r
2
ðr
0
Þ
2
Sð","
0
ÞP
m
1
P
n
1
 
,ð13Þ
where we have used the fact that the spherical harmonic functions 
m
2
m
3
are
mutually orthogonal and 
00
¼ 1=
ffiffiffiffiffiffi
4
p
holds.For the remaining integrals
in the radial direction one integral can be eliminated by using the delta
functions in the kernel Sð","
0
Þ.For this we need the inverse of the band
energy function,which we denote by ,so
u ¼"ðrÞ ()r ¼ ðuÞ
holds.For a general function Gðr,r
0
Þ
Z
1
0
dr
Z
1
0
dr
0
½Sð","
0
ÞGðr,r
0
Þ
¼
X
1

¼ 1
Z
1
0
dr c


e
"
Hð"þ
!ÞGðr,r
0


ðrÞÞ
d
du
ð"þ

 
ð14Þ
holds,where the variable r
0
is evaluated at r
0


ðrÞ ¼ ð"ðrÞ þ
!Þ and H
denotes the Heaviside function.The one dimensional integral on the right
hand side of Eq.(14) has to be evaluated numerically by some integration
formula.However,inserting Eq.(14) into Eq.(13) now yields a symmetric
+ [28.8.2002–7:05am] [431–452] [Page No.438] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
438 RINGHOFER
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
positive definite matrix C,independent of what integration formula is used
for this purpose.Because of the Kronecker symbols in Eq.(13),and since
the polynomial P
0
is a constant,the first row and the first column of the
collision matrix Cwill vanish as well.The semi-discrete version of the entropy
relations is now given through Eq.(11) and the fact that the matrix C is
positive semi-definite.Because of Eq.(11) we have
e
qV
F
T
½A


@
x


F qð@
x


V ÞB


F  ¼ @
x


½e
qV
F
T
B


F,
¼ 1,...,d
Thus,multiplying the semi-discrete Boltzmann Eq.(10) from the left
with e
qV
F
T
and integrating with respect to x gives
Z
R
d
x
dx e
qV
F
T
@
t
F
 
¼
1

2
Z
R
d
x
dx e
qV
F
T
CF
 
 0:ð15Þ
In the case of a time-independent potential,this immediately gives
the estimate that
R
R
d
x
dx½e
qV
jFj
2
 is non-increasing in time.In the case of a
self consistent potential V,the potential is computed from the Poisson
equation

x
V þq D
dop
ðxÞ
1

000
F
000
 
¼ 0,
1

000
F
000
¼ ðxÞ ¼
Z
f ðx,kÞ dk,ð16Þ
(remember that
000
is a constant,) and Eq.(15) only yields the estimate
@
t
Z
R
d
x
dx e
qV
jFj
2
 
 q
Z
R
d
x
dx ð@
t
VÞe
qV
jFj
2
 
which,using the Poisson equation again,can be used to show that some
appropriate convex functional of the solution F decays (see Ref.[8]).
3.THE SPACE–TIME DISCRETIZATION
We now turn to the discretization of the hyperbolic system of con-
servation laws (10) in space and time.In principle,any numerical method,
suitable for hyperbolic systems would be applicable.However,we will
make use of the special structure of the system which is due to the entropy
based construction in the previous section.Moreover,there are some
features of the system (10) which are special to the discretization of the
Boltzmann equation.Standard methods for hyperbolic conservation laws
+ [28.8.2002–7:05am] [431–452] [Page No.439] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
SEMICONDUCTOR BOLTZMANN EQUATION 439
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
would treat the field term and the associated matrices B


,the approxima-
tions of the wave vector derivative,as a source term,when in fact,it is a
transport term and part of the free streaming operator which as a whole is
skew self adjoint under the weight function e
qV
.Moreover,we will also
derive a spatial discretizaton,which is suitable for the direct solution of
the steady state problem,while in the usual methods for hyperbolic
conservation laws steady states are computed by letting the transient
problem converge to a steady state,which is a terribly expensive approach.
We will here only outline the basic idea of the approach and refer the
reader to the papers of Refs.[6,7] for the details.The first idea is based on
the observation that the free streaming operator L in Eq.(1) maps even
functions of k into odd functions and vice versa as long as the band energy
function"is even,which can always be assumed.So from now on,we
assume that"ðkÞ ¼"ð kÞ holds.Accordingly,we split the kinetic density
function f in its even and odd parts
f ¼ f
e
þf
o
,f
e
ðx, k,tÞ ¼ f
e
ðx,k,tÞ,f
o
ðx, k,tÞ ¼ f
o
ðx,k,tÞ,
and obtain the system of Boltzmann equations
@
t
f
e
þLð f
o
Þ þ
1

Q
e
ð f
e
,f
o
Þ ¼ 0,@
t
f
o
þLð f
e
Þ þ
1

Q
o
ð f
e
,f
o
Þ ¼ 0,
where Q
e
,Q
o
denote the even and odd parts of the collision operator,Q.
Moreover,a closer inspection of the collision operator Q in Eq.(2) yields
that,since the scattering cross sections are independent of the angular
variables,Q maps even functions into even functions of k and odd functions
into odd functions;so
Q
e
ð f
e
,f
o
Þ ¼ Qð f
e
Þ,Q
o
ð f
e
,f
o
Þ ¼ Q
o
ð f
o
Þ ¼ ð"Þ f
o
,
ð"Þ ¼ e
"
Z
Sð","
0
Þ dk
0
holds.If we split the basis function space of Maxwellians times
spherical harmonics into even and odd basis functions,and reorder the
coefficient vector F such that the first components correspond to even basis
functions and the last component correspond to odd basis functions,
we obtain that the matrices in the semi-discrete Boltzmann Eq.(10) are
of the form
F ¼
F
e
F
o
 
,A


¼
0 A
eo


A
oe


0
 
,B


¼
0 B
eo


B
oe


0
 
,C ¼
C
ee
0
0 C
oo
 
,
+ [28.8.2002–7:05am] [431–452] [Page No.440] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
440 RINGHOFER
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
and we obtain the system (10) in the form
ðaÞ @
t
F
e
þL
eo
F
o
þ
1

C
ee
F
e
¼ 0,
ðbÞ @
t
F
o
þL
oe
F
e
þ
1

C
oo
F
o
¼ 0,
ðcÞ L
eo
F
o
¼
X
d

¼1
A
eo


@
x


F
o
qð@
x


VÞB
eo


F
o
 
,
ðdÞ L
oe
F
e
¼
X
d

¼1
A
oe


@
x


F
e
qð@
x


VÞB
oe


F
e
 
ð17Þ
The symmetry property (11) becomes,expressed in this split form
A
eo


¼ ðA
oe


Þ
T
¼ B
eo


ðB
oe


Þ
T
ð18Þ
There are of course many approaches one can take for the spatial discretiza-
tion of the system (17),which will in general depend on the grid structure
used and the desired order of the space–time discretization.Rather than
going into the details of a specific discretization,we will outline the basic
underlying idea.Generally speaking,we would like the discretization to
exhibit the following features:
.
The spatial discretization should be suitable for a direct solution of
the steady state problem without having to iterate the transient
method to a steady state.
.
The discretization should reproduce the entropy estimate (15)
on the discrete level.
.
The discretization should yield the correct diffusion limit in the
Hilbert limit !0.
.
In regimes other than the Hilbert limit,where the behavior of the
solution is dominatedby transient waves,we would like to minimize
the amount of artificial diffusion necessary to solve the system.
Of course,the third and the fourth items are somewhat contradictory.
We will start with the spatial discretization.The basic idea is to use a
conservative difference discretization for the even moment Eq.(17a) and to
use a discretization for the odd moments which guarantees a discrete version
of the entropy estimate (15).To this end,it is essential to discretize the even
and odd coefficient vectors F
e
and F
o
on adjoint meshes,i.e.,we will use a
staggered grid approach.Difference operators approximating the partial
derivatives in x


direction will map one mesh into the other.One problem
arises in this approach,since acceleration terms B


F will connect these meshes
as well and therefore,some interpolation would be needed.The easiest way to
+ [28.8.2002–7:05am] [431–452] [Page No.441] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
SEMICONDUCTOR BOLTZMANN EQUATION 441
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
circumvent this problem by rewriting the operator L
eo
in Eq.(17a) as
L
eo
F
o
¼
X
d

¼1
@
x


A
eo


F
o
qVB
eo


F
o
ð Þ þqVB
eo


@
x


F
o
 
:ð19Þ
This way only spatial derivatives of the odd coefficient vector F
o
appear in the equation for F
e
,and,using a staggered grid,these spatial
derivatives are defined on the same mesh as F
e
.For the spatial discretization
of L
oe
in Eq.(17b),we use the fact that L is skew self adjoint under the L
2
scalar product using the weight function e
qV
.Therefore,L
oe
is the negative
adjoint of L
eo
under this weight function and is discretized as
L
oe
F
e
¼
X
d

¼1
ðA
eo


qVB
eo


Þ
T
e
qV
@
x


ðe
qV
F
e
Þ þðB
eo


Þ
T
e
qV
@
x


ðqVe
qV
F
e
Þ
 
:
ð20Þ
Using the identity (18),a direct calculation shows that this form of L
oe
is
equivalent to Eq.(17d).Note that,when written in this form the skew self
adjoint property
Z
R
d
x
dx e
qV
ðF
e
Þ
T
L
eo
F
o
þe
qV
ðF
o
Þ
T
L
oe
F
e
 
¼ 0
is immediate.It is in this form that the spatial discretization is carried out.
Given difference approximations D
eo


,D
oe


of the partial derivatives on the
even and odd meshes which satisfy a discrete integration by parts formula of
the form
I
e
½g
e
ðD
eo


g
o
Þ ¼ I
o
½g
o
ðD
oe


g
e
Þ ð21Þ
for any grid functions g
e
and g
o
defined on the corresponding meshes.Here,
I
e
,I
o
denote some suitable integration rules on the even and odd meshes.
The semi-discrete system (17) is then replaced by
ðaÞ @
t
F
e
þL
eo
F
o
þ
1

C
ee
F
e
¼ 0,
ðbÞ @
t
F
o
þL
oe
F
e
þ
1

C
oo
F
o
¼ 0,
ðcÞ L
eo
F
o
¼
X
d

¼1
½D
eo


ðA
eo


F
o
qVB
eo


F
o
Þ þqVB
eo


D
eo


F
o
,
ðdÞ L
oe
F
e
¼
X
d

¼1
h
A
eo


qVB
eo


ð Þ
T
e
qV
D
oe


e
qV
F
e


þðB
eo


Þ
T
e
qV
D
oe


qVe
qV
F
e


i
,ð22Þ
+ [28.8.2002–7:05am] [431–452] [Page No.442] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
442 RINGHOFER
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
where,for simplicity,we have used the same symbols L
eo
,L
oe
for the differ-
ence and differential operators.Because of reformulating the free streaming
operators L
eo
and L
oe
in the semi-discrete Boltzmann equation in the form
(19),(20) they now map the mesh for the even expansion coefficients into the
mesh for the odd coefficients and vice versa.However,the potential Vis now
needed on the even as well as on the odd meshes.More precisely,Eq.(22c,d)
should really read
ðaÞ L
eo
F
o
¼
X
d

¼1
D
eo


A
eo


F
o
qV
o
B
eo


F
o
ð Þ þqV
e
B
eo


D
eo


F
o
½ ,
ðbÞ L
oe
F
e
¼
X
d

¼1
½ðA
eo


qVB
eo


Þ
T
e
qV
o
D
oe


ðe
qV
e
F
e
Þ
þðB
eo


Þ
T
e
qV
o
D
oe


ðqV
e
e
qV
e
F
e
Þ,ð23Þ
where V
e
is given on the even mesh by the solution of the discrete Poisson
equation

X
d

¼1
D
eo


D
oe


V
e
þq D
dop

1

000
F
e
000
 
¼ 0:
Thus,we still need a suitable interpolation formula to compute V
o
in
Eq.(23) from V
e
.This formula is usually chosen in such a way that the
resulting scheme reduces to the well known Scharfetter-Gummel scheme for
the Drift-Diffusion equations in the Hilbert limit (see Refs.[6,7] for the
details of the interpolation formula and c.f.Ref.[11] for an overview of
the properties of the Scharfetter–Gummel scheme).
Note that,by virtue of construction,the entropy estimate
I
e
e
qV
ðF
e
Þ
T
@
t
F
e
 
þI
o
e
qV
ðF
o
Þ
T
@
t
F
o
 
¼
1

2
I
e
e
qV
ðF
e
Þ
T
C
ee
F
e
 
þI
o
e
qV
ðF
o
Þ
T
C
oo
F
o
  
 0
holds for the resulting systemof ordinary differential equations.The discret-
ization Eq.(22) is particularly convenient for the solution of the steady state
problem(setting @
t
¼ 0 in Eq.(22a,b)).In this case the odd coefficient vector
F
o
can be eliminated altogether by locally inverting C
oo
and inserting into
Eq.(22a).C
oo
is invertible since the kernel of the collision operator only exists
of even functions.So,in the steady state case the equation

2
L
eo
ðC
oo
Þ
1
L
oe
F
e
þC
ee
F
e
¼ 0 ð24Þ
+ [28.8.2002–7:05am] [431–452] [Page No.443] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
SEMICONDUCTOR BOLTZMANN EQUATION 443
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
is solved.By virtue of construction,the operator L
eo
ðC
oo
Þ
1
L
oe
F
e
is a self
adjoint positive semi-definite second order difference operator.If only one
expansion term (the charge density) for F
e
is used,Eq.(24) reduces to a
standard discretization for the resulting steady state drift-diffusion equations.
So,although the momentum as well as the space–time discretization have
been derived on the basis of entropy principles,which are inherently transient
concepts,the resulting structure of self adjoint andskewself adjoint operators
are essential for the numerical solution of the steady state problem as well.
It can be shown (see Ref.[7]) that the linearization of the resulting scheme in
the steady state case is stable precisely because of this structure.
When discretizing the conservation laws (17) in the temporal direction,
a distinction has to be made as to the purpose of solving the time dependent
problem.Either the time dependent problem is solved mainly to compute a
steady state.In this case the focus of the design will be on dealing with the
collision terms and with the resulting temporal stiffness close to the
Maxwellian regime efficiently.This is the idea pursued in Refs.[9,10].
The resulting methods will necessarily exhibit a large amount of diffusion
since they are designed to solve the system close to the diffusive regime.
We are,on the other hand interested in approximate solutions to the
Boltzmann equation quite far away from the Maxwellian regime,i.e.,for
values of the scaled Knudsen number  close to  ¼ 1.In this case the
behavior of the solution will be dominated by dispersive waves,due to the
acceleration induced by the electric field (see Ref.[6]).In order to compute
accurately the transient behavior of the system,it is necessary to resolve
these waves with as little artificial diffusion as possible.Splitting the system
into even and odd parts provides a natural way to do this,namely by
effectively discretizing a second order wave equation with central differences
in space and time.To compute a transient solution,we discretize Eq.(22) as
ðaÞ

t
½F
e
ðt þtÞ F
e
ðtÞ þL
eo
F
o
ðtÞ þ
1

C
ee
F
e
ðt þtÞ ¼ 0,
ðbÞ

t
½F
o
ðt þtÞ F
o
ðtÞ þL
oe
F
e
ðt þtÞ
þ
1

C
oo
F
o
ðt þtÞ ¼ 0:ð25Þ
The Scheme (25) is explicit in time (except,of course,for the solution
of the Poisson equation) in the following sense.Given the solution at time tF
e
at the next time step t þt is computed from Eq.(25a).Once F
e
ðt þtÞ is
known F
o
at the next time step can be computed from Eq.(25b).
The collision terms are always taken implicitly to avoid any resulting stiff-
ness problems for  << 1.A rather straight forward analysis (see Ref.[6])
+ [28.8.2002–7:05am] [431–452] [Page No.444] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
444 RINGHOFER
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
gives a CFL condition of the form t  Oðx=
max
Þ,where x denotes
the spatial mesh-size and 
max
denotes the maximal spectral radius of the
matrix A


in Eq.(10).A condition of this form has to be expected,since
the largest wave speeds in Eq.(17) will be of order Oð1=Þ.Moreover,
the Scheme (25) also reproduce the correct Hilbert limit for small Knudsen
numbers.For  << t,Eq.(25b) becomes approximately F
o
ðtÞ ¼
ðC
oo
Þ
1
L
oe
F
e
ðtÞ,and Eq.(25b) then becomes

t
½F
e
ðt þtÞ F
e
ðtÞ L
eo
ðC
oo
Þ
1
L
oe
F
e
ðtÞ þ
1

C
ee
F
e
ðt þtÞ ¼ 0,
which yields the right diffusive limit.More importantly,however,the
Scheme (25) minimizes diffusion.This can be seen by neglecting the collision
terms in Eq.(25).In the absence of collisions,differencing Eq.(25a) in time
and inserting from Eq.(25b) gives

2
t
2
½F
e
ðt þtÞ 2F
e
ðtÞ þF
e
ðt tÞ L
eo
L
oe
F
e
ðtÞ ¼ 0,
which is the centered difference in space and time approximation to the
second order wave equation

2
@
2
t
f
e
½ðr
k
"Þ  r
x
ðkÞ qðr
x
VÞ  r
k

2
f
e
¼ 0:
This scheme has the property that it exhibits no artificial diffusion,i.e.,
all the eigenvalues of the resulting linear operator lie on the unit circle as
long as the CFL condition is satisfied.It is usually considered as mildly
unstable,since round-off errors and other perturbations will grow linearly
in time.However,the damping terms introduced by the collision operators
usually take care of this linear growth in practice.Thus the only damping in
the Scheme (25) arises from the presence of the collision operators,which is
precisely the same behavior as for the analytical solution.This property of
minimizing the artificial diffusion is extremely important when computing
the transient behavior of the problem,since,c.f.‘‘rise times’’ (the time it
takes to get from one steady state to another approximate steady state) will
be extremely sensitive to the propagation of electron waves and artificial
diffusion will result in underestimating these rise times.
4.A NUMERICAL TEST EXAMPLE
In this section,we present as a numerical test example the simulation
of a standard n
þ
n n
þ
silicon diode in one spatial dimension and in
steady state.The purpose of this example is to demonstrate that the
+ [28.8.2002–7:05am] [431–452] [Page No.445] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
SEMICONDUCTOR BOLTZMANN EQUATION 445
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
presented method is capable to simulate situations which are quite far
from equilibrium.In one spatial dimension the choice of mesh difference
operators D
eo
1
,D
oe
1
in Eq.(21) is quite obvious.We define the grids for the
even and odd coefficient vectors F
e
,F
o
by
M
e
¼fx
0
<   <x
J
g,M
o
¼ y
j
:y
j
¼
1
2
ðx
j
þx
jþ1
Þ,j ¼0,...,J 1
 
and the difference operators by
ðD
oe
u
e
Þðy
j
Þ ¼
u
e
ðx
jþ1
Þ u
e
ðx
j
Þ
x
jþ1
x
j
,ðD
eo
u
o
Þðx
j
Þ ¼
u
o
ðy
j
Þ u
o
ðy
j 1
Þ
y
j
y
j 1
:
If we define the discrete integral operators I
e
,I
o
in the usual way D
eo
and D
oe
satisfy a discrete integration by parts formula of the form (21).
The n
þ
n n
þ
silicondiode has a channel lengthof 50 nm.This means
that the doping concentration D
dop
in Eq.(1c) is given by a step function of
the form
D
dop
ðxÞ
10
24
m
3
0 < x < 50nm
10
21
m
3
50 nm < x < 100 nm
10
24
m
3
100 nm < x < 150 nm
0
@
1
A
The results below were obtained by varying the applied bias (the
difference of the potential V between the right and left endpoint) from
0:1V to 0:5V.The phonon emission/absorption energy!in Eq.(4)
was chosen as 0.036 eV,and computations were performed using a parabolic
band structure ("ðkÞ ¼ ðjkj
2
=2Þ).The plots in Figs.1–8 were obtained using
a uniformmesh with 150 gridpoints in the x-direction (x ¼ 1 nm),and with
16 expansion terms in the energy direction and eight terms in the angular
direction.(For the one dimensional problemthe solution can be expected to
be cylindrically symmetric around the k
1
direction.) So a system of 128
conservation laws has been solved.It was verified that doubling the
number of expansion terms or the number of gridpoints does not signifi-
cantly change the obtained pictures.Figures 1–4 depict the potential as well
as the physical densities for electron,electron velocity and electron energy as
functions of the spatial variable x.Figure 1 shows the potential distribution
VðxÞ.Figure 2 shows the electron density h1iðxÞ ¼
R
f ðx,kÞ dk.Figure 3
shows the velocity distribution uðxÞ ¼ ðh@
k
1
"iðxÞ=h1iðxÞÞ with h@
k
1
"i ¼
R
@
k
1
"ðkÞ f ðx,kÞ dk.Figure 4 shows the corresponding energy densities
given by wðxÞ ¼ ðh"iðxÞ=h1iðxÞÞ.
Figures 5–8 depict the kinetic density f ðx,kÞ for the bias value of 0.4V
for various values of x.In one spatial dimension the kinetic density f will be
F1–F8
+ [28.8.2002–7:05am] [431–452] [Page No.446] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
446 RINGHOFER
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
Figure 2.Particle Density.
Figure 1.Potential.
+ [28.8.2002–7:05am] [431–452] [Page No.447] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
SEMICONDUCTOR BOLTZMANN EQUATION 447
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
Figure 3.Velocity.
Figure 4.Energy Density.
+ [28.8.2002–7:05am] [431–452] [Page No.448] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
448 RINGHOFER
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
Figure 5.r
2
*f:Bias ¼ 0.4 V,x¼30.0752nm.
Figure 6.r
2
*f:Bias ¼ 0.4V,x¼50 nm.
+ [28.8.2002–7:05am] [431–452] [Page No.449] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
SEMICONDUCTOR BOLTZMANN EQUATION 449
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
Figure 7.r
2
*f:Bias ¼ 0.4V,x¼100 nm.
Figure 8.r
2
*f:Bias ¼ 0.4V,x¼120.3008 nm.
+ [28.8.2002–7:05am] [431–452] [Page No.450] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
450 RINGHOFER
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
cylindrically symmetric around the k
1
axis.So Figs.5–8 show the function f
plotted over k
1
and
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
k
2
2
þk
2
3
p
for fixed values of x.
The purpose of these computations is to investigate howfar away from
the fluid dynamic regime the presented approach is applicable.(For such a
short channel and the given applied bias,we expect to see some distinctly
non-equilibrium phenomena).Figure 5 shows the distribution to the left of
the channel.It is essentially given by a forward and a backward traveling
Maxwellian of roughly the same amplitude.Figures 6 and 7 show the dis-
tribution at the beginning and the end of the channel.First,we notice that
the wave has developed a second peak in the forward as well as in the
backward traveling component.Moreover,it has become definitely asym-
metric in the k
1
direction at the end of the channel (in Fig.7).After the
electron has left the channel in Fig.8 remnants of the second peak are still
visible,but the solution has become symmetric in k
1
again.Figure 5 could
have been produced by a hydrodynamic model or even by drift-diffusion
equations,considering the low values of the group velocity in Fig.3 at this
point.The second peak in Fig.8 (after leaving the channel) could have been
produced by a SHE-model.However,the strongly asymmetric distribution
at the channel end (Fig.7) is a truly kinetic phenomenon and could not have
been produced by either of the above approximations.
5.CONCLUSIONS
We have presented a methodology to discretize the Boltzmann equation
for semiconductors based on entropy principles by a spectral approximation in
the wavevector directionandbydifferencemethods inthe space–time direction.
The use of entropy principles results in certain self adjoint properties of the
involved discrete operators which are essential for an efficient numerical
solution of the resulting discrete equations.One of the key ideas in the
presented approach was to split the kinetic density function in its even and
odd parts,and to effectively solve a systemof second order wave equations in
the transient case and a second order elliptic problemin the steady state case.
Although the wave vector discretization consists effectively of polynomial
corrections to Maxwellians the resulting discrete equations are capable of
handling situations quite far away fromthermodynamic equilibriumvery well.
ACKNOWLEDGMENT
This work was supported by NSF grants DMS 970-6792 and INT 960-
3253.(Submitted to TTSP,July 2001).
+ [28.8.2002–7:06am] [431–452] [Page No.451] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
SEMICONDUCTOR BOLTZMANN EQUATION 451
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
REFERENCES
1.Cercignani,C.The Boltzmann Equation and its Applications;Applied
Mathematical Sciences,Springer Verlag:1988;Vol.67.
2.Goldsman,N.;Henrickson,L.;Frey,J.A Physics Based Analytical–
Numerical Solution to the Boltzmann Equation for use in Semicon-
ductor Device Simulation.Solid State Electr.1991,34,389.
3.Levermore,D.Moment Closure Hierarchies for Kinetic Theories.
J.Stat.Phys.1996,83,1021–1065.
4.Majorana,A.;Pidatella,R.A Finite Difference Scheme for Solving the
Boltzmann–Poisson System for Semiconductor Devices,Manuscript
2000.
5.Majorana,A.Spherical Harmonics Type Expansion for the Boltzmann
Equation in Semiconductor Devices.Le Matematice LIII 1998,
331–344.
6.Ringhofer,C.Space–Time Discretization Methods for Series Expansion
Solutions of the Boltzmann Equation for Semiconductors.SIAM
J.Numer.Anal.2000,38,442–465.
7.Ringhofer,C.A Mixed Spectral-Difference Method for the Steady
State Boltzmann–Poisson System.Submitted,SIAM J.Numer.Anal.
2001.
8.Ringhofer,C.An Entropy-Based Finite Difference Method for the
Energy Transport System.Math.Mod.Meth.in Appl.Sci.2001,11,
769–795.
9.Schmeiser,C.;Zwirchmayr,A.Galerkin Methods for the Semiconduc-
tor Boltzmann Equation.Proc.ICIAM 95,Hamburg,1995.
10.Schmeiser,C.;Zwirchmayr,A.Convergence of Moment Methods for
the Semiconductor Boltzmann Equation.SIAM J.Num.Anal.1998,
36,74–88.
11.Selberherr,S.Analysis of Semiconductor Devices,2nd Ed.;Wiley:
New York,1981.
12.Ventura,D.;Gnudi,A.;Baccarani,G.;Odeh,F.Multidimensional
Spherical Harmonics Expansions for the Boltzmann Equation for
Transport in Semiconductors.Appl.Math.Let.1992,5,85–90.
Received July 4,2001
Revised May 19,2002
Accepted May 30,2002
+ [28.8.2002–7:06am] [431–452] [Page No.452] i:/Mdi/Tt/31(4-6)/120015508_TT_031_4-6_R1.3d Transport Theory and Statistical Physics (TT)
452 RINGHOFER
1
200
155
08
_
TT
_
03
1_4–
6
_
R
1.p
df
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924