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 diﬀerence discretizations

in space–time.The resulting class of approximate solutions

dissipate a certain type of entropy.

Key Words:Boltzmann equation;Galerkin methods;Finite

diﬀerences

AMS Classiﬁcation: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 ﬁeld qr

x

V which,using

a mean ﬁeld 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 diﬀerent 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 þ

ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

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 coeﬃcients.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 eﬀective (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 diﬀerence 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 stiﬀ 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

stiﬀness 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 ﬁrst order system of hyperbolic conservation laws

in position and time,which will contain source terms due to the acceleration

of particles by the electric ﬁeld 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 satisﬁes

ð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 coeﬃcients 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 coeﬃcients,which then

automatically satisﬁes 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 diﬃcult 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 ﬁrst 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 diﬀerence 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 ﬁnite 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 coeﬃcients 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 coeﬃcient matrices will have to be computed

numerically.A considerable computational eﬀort 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 coeﬃcient matrices should be done in such a way that the entropy

property (6) translates exactly to the semi-discrete system (10).We will

therefore brieﬂy 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 deﬁne 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-deﬁnite.

Moreover,the ﬁrst 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 simpliﬁed 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=

ﬃﬃﬃﬃﬃﬃ

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 deﬁnite 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 ﬁrst row and the ﬁrst 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-deﬁnite.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 ﬁeld 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 ﬁrst 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

coeﬃcient vector F such that the ﬁrst 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 speciﬁc 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 diﬀusion 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 artiﬁcial diﬀusion 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 diﬀerence 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 coeﬃcient vectors F

e

and F

o

on adjoint meshes,i.e.,we will use a

staggered grid approach.Diﬀerence 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 coeﬃcient vector F

o

appear in the equation for F

e

,and,using a staggered grid,these spatial

derivatives are deﬁned 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 diﬀerence 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

deﬁned 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 diﬀer-

ence and diﬀerential 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 coeﬃcients into the

mesh for the odd coeﬃcients 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-Diﬀusion 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 diﬀerential 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 coeﬃcient 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-deﬁnite second order diﬀerence 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-diﬀusion 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 stiﬀness close to the

Maxwellian regime eﬃciently.This is the idea pursued in Refs.[9,10].

The resulting methods will necessarily exhibit a large amount of diﬀusion

since they are designed to solve the system close to the diﬀusive 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 ﬁeld (see Ref.[6]).In order to compute

accurately the transient behavior of the system,it is necessary to resolve

these waves with as little artiﬁcial diﬀusion as possible.Splitting the system

into even and odd parts provides a natural way to do this,namely by

eﬀectively discretizing a second order wave equation with central diﬀerences

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 stiﬀ-

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 diﬀusive limit.More importantly,however,the

Scheme (25) minimizes diﬀusion.This can be seen by neglecting the collision

terms in Eq.(25).In the absence of collisions,diﬀerencing 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 diﬀerence 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 artiﬁcial diﬀusion,i.e.,

all the eigenvalues of the resulting linear operator lie on the unit circle as

long as the CFL condition is satisﬁed.It is usually considered as mildly

unstable,since round-oﬀ 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 artiﬁcial diﬀusion 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 artiﬁcial

diﬀusion 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 diﬀerence

operators D

eo

1

,D

oe

1

in Eq.(21) is quite obvious.We deﬁne the grids for the

even and odd coeﬃcient 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 diﬀerence 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 deﬁne 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

diﬀerence 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 veriﬁed that doubling the

number of expansion terms or the number of gridpoints does not signiﬁ-

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

ﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ

k

2

2

þk

2

3

p

for ﬁxed values of x.

The purpose of these computations is to investigate howfar away from

the ﬂuid 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 deﬁnitely 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-diﬀusion

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 directionandbydiﬀerencemethods 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 eﬃcient 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 eﬀectively 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 eﬀectively 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 Diﬀerence 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-Diﬀerence Method for the Steady

State Boltzmann–Poisson System.Submitted,SIAM J.Numer.Anal.

2001.

8.Ringhofer,C.An Entropy-Based Finite Diﬀerence 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

## Comments 0

Log in to post a comment