Convergence of Moment Methods for Linear Kinetic Equations.
Numerical methods for linear kinetic equations based on moment expansions for a discretization in the velocity direction are examined. The moment equations are hyperbolic systems which can be shown to converge to the kinetic equation as the order of the expansion tends to infinity and to a drift-diffusion model as the Knudsen number tends to zero. A discretization of the moment equations with respect to time and space is presented, a stability result is proven, and some aspects of an implementation are discussed. In particular, an adaptive procedure is described where the order of the expansion is determined locally. Results of numerical experiments are presented.
Introduction
The numerical simulation of flows governed by macroscopic models in certain
subdomains and by kinetic models in others, has received significant attention
recently. Applications can be found in the fields of gas dynamics and
semiconductor device simulation. One obvious approach to these problems
is domain decomposition where finite difference or finite element methods
are used for the macroscopic equations and particle methods for the kinetic
equations. The coupling is performed by interface conditions. In the semi-conductor
field this approach has been carried out for a particular application
in [4]. Theoretical investigations can be found in [3]. The basic idea is that
the kinetic model should be solved in the whole domain, however, the computationally
much less expensive macroscopic model is sufficiently accurate
in regions where scattering processes dominate.
An alternative approach used in the semiconductor field is to extend the
basic macroscopic model, the drift-diffusion model, by adding equations for
higher order moments of the distribution function. The resulting models are
known as the hydrodynamic semiconductor equations. Recently, efforts have
been made to extend this idea in a systematic way. We mention expansions
of the distribution function in terms of spherical harmonics [12] and a second
approach [10] in the spirit of Grad's work [2] where an expansion procedure
is applied to a transformed kinetic equation with an adaptively determined
coordinate transformation.
Here we are concerned with moment expansions of arbitrary order and
with questions of stability and convergence. For the stability analysis it is of
great help to formulate the expansions within the framework of Levermore's
recent theory [6].
Our second emphasis is on the implementation of the moment methods.
In particular, we want to demonstrate that a transition between macroscopic
and kinetic regions can be performed by varying the order of the expansion.
In section 2 the kinetic model and its basic properties are presented. It
is a well known result that the drift-diffusion model is obtained in the limit
where the Knudsen number (mean free path/characteristic length) tends to
zero. The moment expansion procedure is presented in section 3. Section
4 contains stability and convergence results as the order of the expansion
tends to infinity. In section 5 it is shown that the moment equations also
reduce to the drift-diffusion model in the macroscopic limit. A semiimplicit
time discretization and a space discretization for the one-dimensional case
are presented in section 6. Finally, we discuss some further aspects of an
implementation and present results from our solutions of a model problem
in section 7.
2 The kinetic model
We consider the nondimensionalized transport equation
@t
with the Knudsen number (ratio of the mean free path to a characteristic
length) ", the electric potential \Phi(x; t) and a collision operator of the form
Z
R
denotes the normalized Maxwellian
and oe the collision cross section, assumed nonnegative and O(3)-
invariant, i.e., oe(Rv; Rv 0 orthogonal matrix R. The
prime in the above integral denotes evaluation at v 0 . We shall also need
boundedness away from zero of the cross section, i.e.,
as well as boundedness of the scattering rate, i.e.,
Z
The Boltzmann equation (2.1) models the transport of electrons in a semi-conductor
crystal under the physical assumptions of a parabolic energy band
and nondegeneracy. A further simplification is the positivity assumption
(2.2) meaning that Q is a regularized version of physically relevant models.
An important tool for the analysis of the Boltzmann equation is the following
symmetry property [9]:
Z
Q(f)g
dv
Z
f
!/
This suggests the introduction of the Hilbert space
with the corresponding inner product h\Delta; \Deltai H and norm k\Deltak H . Essentially all
the important properties of Q can be deduced from (2.4) [9]:
Lemma 2.1 Assume (2.2), (2.3). Then Q is a bounded operator on H. Also
. The orthogonal projection
of H onto kern Q is given by
R
and the following coercivity
estimate holds for any f 2 H:
The initial value problem for the Boltzmann equation is considered where
holds for (x; v; subject to the initial conditions
For the initial datum we require
The electric potential \Phi(x; t) is considered given and satisfying
Then, by standard transport theory, the problem (2.1), (2.6) has a unique
mild solution f 2 C([0; 1);H). An estimate for the solution is obtained by
taking the scalar product (in H) of the Boltzmann equation (2.1) with f ,
then multiplying by e \Gamma\Phi and integrating with respect to x:
dt
Z
e \Gamma\Phi kfk 2
Z
e \Gamma\Phi hQ(f); fi H dx
e \Gamma\Phi @ \Phi
@t
The semidefiniteness of Q and the Gronwall lemma imply
@ \Phi
@t
The significance of this estimate is in its uniformity with respect to ". In the
converges to %(x; t)M(v), where % solves
an initial value problem for the drift-diffusion equation [9]:
@%
@t
Z
f I (x; v)dv ; (2.11)
where the mobility is given by
Z
The fact that the drift-diffusion model is isotropic (i.e., the mobility is a
scalar) is a consequence of the O(3)-invariance of the collision cross section.
3 The Method
By Lemma 2.1 f 2 =(2M) is an entropy density for the Boltzmann equation
(2.1). Then, by Levermore's results [6], a moment method with favourable
properties is obtained in the form of a Galerkin method with trial functions
of the form ' k test functions m k (v).
Here we also assume m 0 such that the equilibrium distributions are
contained in the trial space, and that f' k g is a complete orthonormal system
in H. Thus, the orthogonality relations
Z
hold. Furthermore, we assume that the ' k are smooth and that they vanish
at infinity, i.e., ' k (v) ! 0 as jvj !1. An example for the choice of the m k
are Hermite polynomials (see [2]).
For we approximate the solution of (2.1), (2.6) by fK (x; v;
where the
coefficient vector ff satisfies hL(ff \Delta This system
of equations can be written in the form
@t
"A i
@ff
@ \Phi
where the summation convention has been used and the elements of the
matrices A i , B i and C are given by

@' l
The initial values of the components of ff are the Fourier coefficients of the
initial datum for f :
We collect some properties of the coefficient matrices:
Lemma 3.1 For 3, the matrix A i is symmetric, the first line of B i
vanishes and B i +B tr
holds. The matrix C is symmetric and has the
. ~
with ~
negative definite.
Proof: The symmetry of the A i is obvious. The elements of the first lines
of the B i are
@' l
@' l
Integration by parts gives
Z
@' l
Z
dv
The symmetry of C follows from (2.4). The elements of the first line are
Z
and, by symmetry, also the first column of C vanishes.
Given a vector we shall in the following use
the notation ~
. The orthogonality of the ' k implies
2.1 we have
~
ff tr ~
4 Stability and Convergence
By the symmetry of the A i , (3.1) is a symmetric hyperbolic system and
existence and uniqueness of a solution of (3.1), (3.2) follows from standard
results [7].
An equation analogous to (2.9) for the moment system is obtained by
taking the scalar product of (3.1) with e \Gamma\Phi ff and integrating with respect to
x:
dt
Z
Z
e \Gamma\Phi ff tr
e \Gamma\Phi @ \Phi
@t
Here we denote by j \Delta j the Euclidean norm of vectors. Observing that
I (x; \Delta)k H hold, we have analogously to (2.10)
@ \Phi
@t
Note that this estimate is uniform both in terms of " and K. We collect our
results:
Theorem 4.1 With the assumptions of section 1 the problem (3.1), (3.2)
has a unique solution ff
satisfies (4.2) for every t - 0.
For a convergence result as K ! 1 some smoothness of the solution of
the Boltzmann equation is needed. We therefore make additional assumptions
on the initial data,
and on the collision cross section,
Z
Z
Lemma 4.1 Assume (2.2), (2.3), (2.7), (2.8), (4.3), (4.4), and (4.5) hold.
Then for the solution of (2.1), (2.6) we have
Proof: We introduce the notation
. Formal differentiation
of the Boltzmann equation with respect to x i and v i gives
@t
@t
@
Z
@oe
Z
We consider (4.6), (4.7) as a system of 6 transport equations for
subject to the initial conditions
3: (4.9)
We intend to show that the right hand side of (4.6), (4.7) is the sum of an
element of (C([0; 1); H)) 6 and a linear, bounded (on H) operator applied
to F . For the right hand side of (4.6) this follows directly from Lemma 2.1
and (2.8).
Consider the operator ~
defined by ~
R
@oe
In
the same way as assumption (2.3) implies boundedness of Q, (4.5) implies
boundedness of ~
Q on H. Thus, by (2.10), ~
holds. The
same holds for the second term on the right hand side of (4.8) by the estimate
(Cauchy-Schwarz inequality)
Z
Z
'Z
dv
Z
Z
and assumption (4.4). The last term in (4.8) is taken care of by assumption
(2.3).
We conclude that a unique mild solution F 2 (C([0; 1); H)) 6 of problem
(4.6), (4.7), (4.9) exists. A standard argument shows that the components
of F are the weak derivatives of the solution of (2.1), (2.6).
We shall give a convergence result as K ! 1 for the case that the
functions m k are polynomials. Let
denote the one-dimensional Hermite polynomial of order k [1], normalized
such that
Then the functions
are a complete orthonormal system in H. In the following we shall need the
leading order coefficients:
We shall consider the moment method where f(x; v; t) is approximated by
m=0 ff klm (x; t)' klm (v).
Theorem 4.2 With the assumptions of Lemma 4.1 there exist sequences of
integers such that for the moment method with Hermite po-
lynomials, fKnLnMn converges to the solution f of (2.1), (2.6) in C([0; T
as n !1 for every T ? 0.
Proof: We introduce the projection
the system f' klm g is complete,
tends to zero in C([0; T 1. Writing the
error in the moment method as
r klm (x; t)' klm
it remains to prove that jrj tends to zero in C([0; T ]; L 2 (IR 3 )). We compute
This is an inhomogeneous version of (3.1) for r. We also have
Taking the scalar product of (4.11) with e \Gamma\Phi r and integrating with respect
to x gives (compare to (4.1))
dt
Z
@ \Phi
@t
Z
e \Gamma\Phi jrj 2 dx
Z
After a lengthy but straightforward computation using (4.10) we obtain
s
@f
s
@f
s
@f
For the last term we have
'Z
'Z
Z
Z
and the first factor on the right hand side tends to zero as minfK; L; Mg !
1.
By Lemma 4.1 we haveX
Z

@f
This implies that for a sequence K n !1
Z

@f
holds. Analogously, sequences L can be found such that
Z

@f
Z

@f
We conclude that the last term in (4.12) can be estimated by
Z
'Z
e \Gamma\Phi jrj 2 dx
with
lim
sup
Applying the Gronwall lemma to (4.12) gives the desired result.
5 The drift-diffusion limit
Before carrying out the macroscopic limit " ! 0 in the moment system we
rewrite (3.1) in terms of new variables introduced by
@%
@t
@
@t
@ \Phi
~
'i H
@%
@ \Phi
with
obtained from A i ,
by deleting the first rows and columns.
Theorem 5.1 With the assumptions of section 1, fK converges to %KM in
loc ((0; 1);H) weak* as " ! 0, where %K solves the drift-diffusion problem
(2.11) with the mobility tensor
'i tr
~
Proof: From the uniform estimate (4.2) weak* convergence of fK , and therefore
also of %, for a suitable subsequence is immediate. From (4.1) and the
coercivity of ~
C we obtain a uniform estimate for fi and, thus, also weak*
convergence for a further subsequence. The drift-diffusion problem for the
limit of % follows from going to the limit in (5.1) and (5.2) in the sense of
distributions.
From the negative definteness of ~
C and formula (5.3) it is immediate that
positive semidefinite. Since standard results for parabolic equations
imply uniqueness of the solution of the limiting problem, convergence is not
restricted to subsequences.
A few words have to be said about the mobility tensor. First, it is an
approximation for the exact, isotropic mobility from section 1, determined by
using the moment expansion method for computing approximations h iK 2
of the functions h
This has the solution
'(v) tr ~
Now we compute
Z
'i tr
~
The exact mobility is obtained whenever the h i are members of the orthogonal
system f' k g. A second comment is that the approximating mobility
tensor can itself be made isotropic by choosing the approximating subspaces
O(3)-invariant, a requirement also formulated in [6].
Lemma 5.1 The mobility tensor (5.3) converges to the mobility defined by
(2.12) as K !1.
Proof: Denoting by \Pi K the orthogonal projection of H onto span('
we write the error in h i as
The first term on the right hand side tends to zero in H because of h
and the completeness of the orthonormal system. We also have
From the coercivity of Q (or of ~
C) follows
implying convergence of h iK to h i in H by the continuity of Q. Convergence
of the mobility is an immediate consequence.
6 Time and Space Discretization
For a numerical implementation the problem (3.1), (3.2) has to be discreti-
zed. Finding an efficient discretization is nontrivial since the stiffness of the
problem goes to infinity both for " ! 0 and for K !1. Here we concentrate
on the stiffness caused by smallness of ".
First we want to make a general statement about the numerical solution of
problems parameter ". If the numerical method creates
a problem P ";N where N is a measure for the numerical effort, then an efficient
method should satisfy the requirement that the limits "
carried out in P ";N commute. Assuming that the method is convergent for
fixed means that P 0;N is a convergent
method for the limiting problem P 0 . This property is an indication that the
solution of the problem with small " can be approximated with a numerical
effort remaining bounded as " ! 0.
For stiff problems implicit time discretizations typically satisfy the above
requirement whereas explicit discretizations do not. Since particle methods
(Monte Carlo methods as well as deterministic particle methods) lead to
explicit discretizations they are inefficient for small Knudsen numbers.
The results of the preceding section show that the moment expansion
method satisfies the above requirement. However, because of the stiffness
this property can easily be destroyed by the discretization of the moment
system. In particular, discretizations explicit in time would suffer from the
same drawbacks as particle methods. The use of implicit discretizations
on the other hand seems to be out of the question because of the size of
the linear systems to be solved in this case. In the following we present a
semiimplicit discretization with respect to time having the property that it
reduces to a standard implicit discretization of the drift-diffusion equation
in the macroscopic limit. Afterwards a space discretization for the one-dimensional
case will be given.
We discretize the moment system in the form (5.1), (5.2) in time by evaluating
the terms on the right hand side at the new time step and the transport
terms on the left hand side of (5.2) at the old time step. Introducing the
time step \Deltat and t n\Deltat, we set \Phi n
and fi(x; t n ) by % n (x) and fi n (x), respectively. Then the semidiscretization
with respect to time is given by\Deltat
@
'i
\Deltat
@ \Phi n+1
~
'i
@% n+1
From the formal limit as " ! 0 the variable fi n+1 can be eliminated because
of the nonsingularity of ~
C. The result is the implicit Euler discretization of
the drift-diffusion equation (2.11) with the mobility (5.3).
The one-dimensional version of (6.1), (6.2) reads\Deltat
@
@x
'i
\Deltat
A
@x
@ \Phi n+1
@x
~
'i
@% n+1
@x
@x
For the space discretization we use upwinding (see, e.g., [5]) for the transport
terms in (6.4) and the Scharfetter-Gummel discretization (see, e.g., [8]) for
the drift-diffusion flux
@x
@x
We introduce x
respectively. The discrete equations read\Deltat
'i
\Deltat
'iJ n+1
with
(B(\Phi n+1
and
~
~
where \Phi n
z \Gamma1 is the Bernoulli function, and ~
~
A \Gamma is the standard decomposition of a symmetric matrix into matrices
with nonnegative and, respectively, negative eigenvalues. The Scharfetter-
Gummel discretization (6.8) of (6.5) is obtained by approximating J n+1 and
@ \Phi n+1
@x in (6.5) by constants and explicitly solving the resulting ODE for % n+1 .
It is an example of an exponentially fitted discretization and provides automatic
upwinding if the convection term in (6.5) dominates. The stability
result at the end of this section shows that upwinding and the Scharfetter-
Gummel discretization are reasonable choices. For the optimal resolution of
sharp fronts, however, more sophisticated discretizations (such as, e.g., ENO
methods [11]) might be preferable.
The system (6.6), (6.7) can be rewritten by eliminating fi n+1
'i
with
'i tr
\Deltat
and
\Deltat fi n
For every time step the following computations have to be carried out: For
j+1=2 at every point of the spatial grid a linear system with the symmetric,
positive definite coefficient matrix " 2
\Deltat
C has to be solved and an additional
system for the computation of -
-. Then the density is computed by solving a
linear system of the size of the spatial grid.
Finally, we check the property mentioned at the beginning of this section:
Obviously, as " ! 0, the above system reduces to a standard discretization of
the drift-diffusion equation commonly used in device simulation software. We
conclude this section with an unconditional stability result for the discrete
system for coarse spatial grids and time independent electric potential. We
shall need the norms
and an according definition of k% n k, where fi n
j+1=2 is replaced by % n
.
Theorem 6.1 Assume the potential is time independent and uniformly bo-
unded, i.e., \Phi n
. Then there are positive
constants - and c (independent of ", h and \Deltat) such that for "
sup
holds for the solution of (6.6)-(6.8).
Proof: We define an interpolated value of e \Gamma\Phi between grid points by
Now we multiply (6.6) by e \Gamma\Phi j % n+1
and (6.7) by (e \Gamma\Phi ) j+1=2 fi n+1
j+1=2 , add the
resulting equations and sum over j:\DeltatX
\DeltatX
After using the inequality u in the first two sums,
we multiply by \Deltat, sum from use the obvious estimate
and the negative definitness of ~
C:
e \GammaM k% N
e \GammaM oe
This implies the result of the theorem with
and an appropriately
chosen value of c.
7 Implementation and Numerical Results
Our main motivation for considering moment methods are problems with
both kinetic and macroscopic regions. Recently, strategies for the numerical
solution of such problems have been developped [3], [4] where interfaces between
macroscopic and kinetic regions are introduced. Then the discretized
drift-diffusion equation is solved in the macroscopic regions and a particle
method is used in the kinetic regions. A major issue in this approach is the
derivation of appropriate coupling conditions at the interfaces. The moment
expansion method considered in this work allows for a gradual transition
between macroscopic and kinetic regions by varying the order of the expansion
locally in space and time. In the implementation described below this
is done adaptively by a heuristic procedure. Theoretically motivated error
estimators still have to be developped.
Another implementation issue is the time step selection. It is subject
to accuracy as well as stability considerations. Since both the space and
time discretization are first order accurate, time steps of the same order of
magnitude as grid spacings seem to be reasonable from the accuracy point of
view, at least in the absence of appropriate estimators for the discretization
error. As far as stability is concerned, the result of the previous section shows
that time step limitations can only occur for O(") (or finer) grid spacings.
In this case, the discrete transport operator T h becomes important and we
expect a CFL condition (see [5]) for the time step. A quantitative criterion
for the time step selection is derived from the scalar model equation
@t
@x
with This is a hyperbolic equation with
a strong relaxation term on the right hand side and a caricature of (5.2).
Without loss of generality we assume a ? 0 and, thus, the discretized version
is
\Deltat
We perform a von Neumann stability analysis (see [5]). A simple computation
shows that u n
solves (7.2). The stability criterion j-j - 1 is satisfied if either
Case 1:
2"a
\Gammac and \Deltat -
or
Case 2: h -
2"a
Case 1 corresponds to a fine spatial grid O(")). The restriction for the
time step is a relaxed version of the CFL-condition for the case without the
relaxation term. Case 2 corresponds to the result of the previous section.
The stability analysis for the scalar equation (7.1) cannot be rigorously
applied to our discretization of the system (5.2) since in general the matrices
~
A and ~
C cannot be diagonalized simultaneously. Also the effect of the electric
field term of (5.2) is not present in (7.1). However, the above analysis is
used for the time step control in the implementation described below with
satisfactory results. The stability criterion is checked with a replaced by
the eigenvalue with largest modulus of ~
A and \Gammac replaced by the smallest
eigenvalue of \Gamma ~
C.
We solved Ringhofer's test problem [10], which is one-dimensional in space
with periodic boundary conditions. The collision operator is the relaxation
time approximation with oe(v; v 0 const. In this case the Boltzmann
equation has solutions of the form f(x
3 , where f solves a one-dimensional
(in both position and velocity) Boltzmann equation. We have
used the moment method with one-dimensional Hermite polynomials for the
and the Lax-Friedrichs discretization (see [5]) of the transport term
instead of upwinding (thus avoiding the decomposition ~
The time independent potential forms a well in the right part of the in-
terval, and the initial datum is a concentration of electrons in the left part
of the interval with a distribution in local thermal equilibrium. As
the solution converges to an equilibrium distribution with the electrons concentrated
in the potential well. The evolution of the density is shown in
Fig. 7.1. Fig. 7.2 contains a comparison of computed results at
the drift-diffusion model and the moment method with different orders. The
kinetic solutions differs strongly from the drift-diffusion solution. Reasonable
accuracy can already be achieved with order Finally, we have
implemented an adaptive method where the order of the moment expansion
is chosen individually for every grid point in space and time. Fig. 7.3 shows
the evolution of the orders. Note that the order is automatically reduced
to the minimal value (drift-diffusion model) as the solution converges to an
equilibrium distribution for t !1. This shows that an automatic transition
between macroscopic and kinetic regions is possible. In the presented example
the adaptive procedure reduces the computational effort by a factor of 4
as compared to the use of the maximal order throughout the computation.

Figure

7.1: Electron density vs position and time,
K=8
K=24
K=32
drift-diffusion
x

Figure

7.2: Electron density vs position at time
the drift-diffusion model and the moment expansion with different orders K.
order

Figure

7.3: Adaptive order of the moment expansion vs position and time.



--R


On the kinetic theory of rarefied gases

Simulation des Ladungstransportes in elektronischen Bau- elementen mit Hilfe der Monte-Carlo-Methode
Numerical Methods for Conservation Laws
Moment closure hierarchies for kinetic theories
Compressible Fluid Flow and Systems of Conservation Laws in Several Space Dimensions
The Stationary Semiconductor Device Equations
Diffusion approximation of the linear Boltzmann equation: Analysis of boundary layers
An expansion method for the Boltzmann transport equation using an adaptive Galerkin procedure
Efficient implementation of essentially non-oscillatory shock capturing schemes
Multidimensional spherical harmonics expansion of Boltzmann equation for transport in semi- conductors
--TR

--CTR
Matthias K. Gobbert , Timothy S. Cale, A kinetic transport and reaction model and simulator for rarefied gas flow in the transition regime, Journal of Computational Physics, v.213
moment expansion;semiconductors;kinetic equations