Asymptotic-Induced Domain Decomposition Methods for Kinetic and Drift Diffusion Semiconductor Equations.
This paper deals with domain decomposition methods for kinetic and drift diffusion semiconductor equations. In particular accurate coupling conditions  at the interface between the kinetic and drift diffusion domain are given. The cases of slight and strong nonequilibrium situations at the interface are considered and numerical examples are shown.
Introduction
Semiconductor device simulations are usually done on the basis of drift diffusion
or hydrodynamic equations. However in view of the ongoing miniaturization
of semiconductor devices a more accurate modelling of the physics seems to be
necessary. This is achieved by kinetic semiconductor equations. We refer to the
books of Markowich et al. [16] and Selberherr [18] for a detailed description
of semiconductor equations and further references. However, it is in general
not necessary to model the whole region in the device by the computationally
expensive kinetic equation. Only in particular sensitive regions, where the
solution is far away from an equilibrium state, as for example in boundary
layers, kinetic equations have to be employed. In the remaining regions of
the device the usual drift diffusion equation may be valid, and will lead to
sufficiently accurate results.
Domain decompositon methods are a natural tool to design in this case accurate
numerical codes with reasonable computation times solving both types of
equations in their respective domains. This can be achieved similiar to the field
of gas dynamics, where domain decomposition methods have lead to the development
of hybrid codes for kinetic and hydrodynamic equations. Examples can
be found in [3, 15].
The general aim is to approximate the global kinetic solution by the solution
of the domain decomposition problem, computed with the hybrid code, as accurately
as possible, saving in this way a considerable amount of computing
time. Once the locations of the kinetic and drift diffusion domains in the device
are known, the main point of the problem is to obtain the correct coupling
conditions at the interface between the two regions.
In the present paper the physics of the semiconductor device is assumed to be
described by a kinetic transport equation. We restrict here to kinetic semiconductor
equations with a linear collision term and given electric field. However
the following analysis should be extendable to nonlinear equations and selfcon-
sistent treatments of the electric field.
In section 2 we state the physical problem and the equations describing it. In
section 3 we state coupling conditions for equilibrium situations at the inter-
face. Using the equality of fluxes, a condition that is due to Maxwell in the
gas dynamics case and to Marshak in the radiative transfer case, or simply the
equality of moments at the interface and an O(ffl 2 ) approximation of the distribution
function, where ffl is the mean free path, one obtains correct results.
These conditions are proven to yield O(ffl) approximations for a model situation.
In section 4 nonequilibrium situations at the interface are considered. Using
asymptotic analysis similiar to the usual boundary layer considerations as, e.g.,
in Bardos et al. [1] and Bensoussan et al. [2] for neutron transport, Cercignani
[4] and Sone and coworkers [19] in the gas dynamics case and Poupaud [17],
Golse/Poupaud [8] and Yamnahakki [20] for semiconductor equations, we develop
accurate coupling conditions in this case by an analysis of the interface
layer between the two domains. This leads to kinetic linear half space problem.
The essential point is then to find a fast approximate solution procedure of the
space problem. This will be discussed in the appendix. It is an extension of
the work done in Golse/Klar [7] and related to the so called variational meth-
ods, see Loyalka [14] and Golse [6] for gas dynamics and Larsen, Pomraning [13]
in the radiative transfer case. In section 5 we show some numerical results for
a one dimensional model problem. As one can see in the simulations the developed
coupling conditions lead to a very accurate approximation of the kinetic
solution by the solution of the domain decomposition problem.
2 The Basic Equations
We consider a domain D s in R 3 representing the semiconductor device geometry.
D s is divided into two not overlapping subdomains D 1 and D 2 ,
with smooth boundaries @D s and the interface I
The basic description of the physics in the device is given by the following linear
semiconductor Boltzmann equations in D
(v
are the distribution functions of the electrons in D i , x 2
e.g. Markowich et al. [16] for a discussion of these
equations.
is the electric field multiplied by the absolute value of the elementary
charge and divided by the effective mass.
Z
is the collision operator with
and
where ' is the temperature of the device. models the generation
recombination process. Defining
Z
and
Z
we get
The global kinetic solution is then defined by the solution of equations (1) in
with the initial and boundary conditions
and
where n(x) is the outer normal at @D s at the point x and
and
The interface condition is simply continuity of the distribution functions for
The mean free path is assumed to be small such that an approximation
of the kinetic equation by the drift diffusion equation is valid there, see,
e.g., Poupaud [17] for a mathematical treatment.
This means, f 2 (x; v; t) is approximated by ae(x; t)M(v), where ae(x; t) fulfills
with the diffusion coefficient
Z
where h i is the solution of
Z
and the mobility
We set in the following to simplify the notation.
The aim is now to approximate the global kinetic solution by the solution of
the following coupling problem: In D 1 the kinetic equation (1) is solved with
the above inital conditions in D 1 and the boundary condition at @D 1 " @D s . In
D 2 the drift diffusion equation (2) is solved with initial conditions
and boundary conditions
Equiping these equations with coupling conditions at the interface I will lead
to a properly stated problem. The solution depends however strongly on the
type of coupling conditions employed. In the next sections we discuss several
possibilities and investigate them.
Coupling Conditions for Near Equilibrium States
For equilibrium states at the interface coupling conditions can be straightforwardly
obtained. One starts by defining a 'macroscopic distribution function'
approximating the kinetic distribution function in D 2 to O(ffl 2 ). Using an approximation
to O(ffl) leads without overlapping of the domains to obviously
wrong results, see section 5, Figure 1.
Assuming, as in [1, 2, 17], that f 2 in D 2 is of the form
substituting it into the scaled equations (1) leads to
Q(f (0)
Q(f (1)
Therefore due to (3) and due to the properties of the collision kernel, see, e.g.,
Poupaud [17] we obtain
This means that the O(ffl) macroscopic distribution function would be given for
f macro (x; v;
(4) leads to
f (1)
Therefore we get for x 2 D 2
f macro (x; v;
as the O(ffl 2 ) macroscopic distribution function in D 2 . Using the latter one, one
obtains the following conditions:
The equality of moments:
Determine ae(-x; t); -
I by
Z
Z
f macro (-x; v; t)dv
Z
The ingoing function for the kinetic region is given by
for is the outer normal at @D 1 at - x.
The Maxwell-Marshak conditions:
Determine ae(-x; t); -
I by
Z
Z
Z
with n(-x). The ingoing function for the kinetic region is given as before.
The solution of equation (1) in D 1 and (2) in D 2 with the initial, boundary and
coupling conditions as above will be called the solution of the coupling problem.
It can be found numerically by a standard Schwarz type iteration procedure:
One determines the coupled solution by solving in turn the equations in the two
domains. In each step the above interface conditions are used. This is repeated
until a stationary state is obtained. See Klar [12] for an investigation of the
coupled solution and of the iteration procedure in a general case.
In the following we show for a model situation that these coupling conditions
lead to a correct approximation in case of near equilibrium situations.
We consider the stationary kinetic equationsffl i
(v
with ffl. The stationary drift
diffusion equation is for
with x 2 D 2 . We restrict to a slab x 2 D sub-domains
are B]. The boundary conditions are f(A;
ae field is of the form
Equation (5) in D 1 is coupled to equation (6) in D 2 by
one of the above coupling conditions. f 1 and ae denote the coupled solution in
the respective domains. The reference solution is the global stationary kinetic
solution, i.e. the solutions f of the kinetic equation (5) in D 1 and D 2 with
We denote the global solution by f(x; v; t); x 2 D s . One
can then prove:
Theorem 1
The solution of the coupling problem for the stationary equations (5) and (6)
with one of the above coupling conditions at the interface is an O(ffl) approximation
of the global stationary kinetic solution in the whole domain, i.e.
where C is some constant.
We need a trivial lemma:
Let ae(x) be the solution of the drift diffusion equation (6) in D
the Dirichlet boundary data ae(0) at Then the value
@ x ae(0) is given by
with
and
Proof of the Theorem:
We restrict to the Maxwell-Marshak conditions. They give together with the
Lemma Z
Z
Z
and for
From (7) we obtain
'Z
Substituting this into (8) a short calculation gives for the following
reflecting boundary condition for the kinetic part of the coupled solution
Z
with
Considering now the kinetic equation in D 1 with the fixed boundary conditions
ae A at A and the above reflecting boundary condition at 0, we observe that the
approximation of this kinetic problem for small ffl by a drift diffusion equation is
similar to a problem in neutron transport, that has been treated in Bensoussan
et al. [2], section 3.6.
A straightforward extension of the results obtained there leads to an O(ffl) approximation
of the kinetic solution f 1 in D 1 with the above reflecting boundary
conditions at 0 by a drift diffusion equation with the following boundary condition
at
with
Z
Z
dv
Z
Z
dv
Z
Z
Z
Z
Z
Z
dvH
Z
Z
Z
Z
This gives
as the boundary condition for a drift diffusion approximation of f 1 in D 1 .
Considering now the global kinetic solution in [A; B] we know that it is approximated
to O(ffl) by the global drift diffusion equation, see [17]. However,
the solution of the global drift diffusion equation in D s is in D 1 equal to the
solution of the drift diffusion equation in D 1 with the above boundary condition
at 0. This can be seen by an easy calculation.
Therefore finally we get an O(ffl) approximation of the global kinetic solution by
the coupled solution in D 1 . The corresponding property in D 2 is then obvious.
This result shows, that the equality of moments or fluxes is appropriate for near
equilibrium situations at the interface. If this is not the case, they might give
wrong results as shown in section 5.
In the next section we develop accurate however easy to handle coupling conditions
for this case.
Coupling Conditions for Nonequilibrium States
In contrast to section 3 non equilibrium situations might prevail here at the
interface between the two domains. We use here not only an O(ffl 2 ) approximation
in the bulk of the device, but also on the boundaries and interfaces to
develop the correct conditions. This includes an analysis of the interface layer
between the two domains. We concentrate in this section on the interface layer
neglecting the boundary layers.
Proceeding similarly as in the usual boundary layer treatment, see the references
in the introduction, one assumes the macroscopic distribution function in D 2
to be the same as in the last section plus an additional interface layer term
concentrated around the interface I with a thickness of the order of a mean free
are interface layer coordinates given by the
x be the point on the interface closest to x. The coordinate
along the normal to the interface is y 0 , i.e.
is the
outer normal at @D 1 at -
I. The interface itself is parametrized by ff i , i.e.
depend on the point x. See, e.g., [4],[19] for details
on boundary layer coordinates.
Requiring that the kinetic equation in D 2 is fulfilled by f macro to O(ffl 2 ), i.e.,
we obtain the following equations by substituting f macro from (9) into (10) and
comparing terms of the same order in ffl:
Up to order ffl \Gamma 1 we get the following halfspace problem at -
x
with n(-x). To order 1 we obtain
and another halfspace problem at the interface point -
x with y 2 [0; 1):
Remark that equation (12) is fulfilled since
To order ffl we must have
and a linear half space problem for -
W must be fulfilled. Equation (14) is
solvable, if the drift diffusion equation is fulfilled by ae(x; t).
Assuming continuity of the distribution functions at the interface to O(ffl 2 ), i.e.,
for
we obtain from (9) that the sum of the ingoing functions for the half space
problems for v has to fulfill
The splitting of this sum into -(0; ff
arbitrary. We
choose,
and
since the following becomes simpler for this choice. Moreover, the
asymptotic values of the halfspace problems must fulfill
since the interface layer term should vanish in the bulk.
We introduce two more notations:
be the solution of the half space problem (11) at - x with
we define the asymptotic value
which is independent of v and the Albedo operator
Let ~
now be the solution of the half space problem (11) at the
interface point - x with ingoing function
~
From (16) and (17) we have
From
This means, ae(-x; t); -
I is determined by the condition
Knowing this is an equation for ae(-x; t). Looking at equation (9) and
keeping equation (15) in mind one finds the ingoing function for the Boltzmann
region for -
\Delta(r x ae(-x;
The second equality is due to (16) and the definition of the Albedo operator.
Knowing ae this is an equation for f 1 giving f in terms of f It is
easily handeled in the framework of the Schwarz iteration procedure mentioned
in the last section.
It is now most important to have a fast however accurate enough solution procedure
to find the asymptotic state and the outgoing distribution of the half
space problem. We remark, that solving the whole half space problem, for ex-
ample, by a standard discretization procedure would need a lot of computing
time, in particular, since it has to be solved at each point of the interface. In
the appendix this will be achieved by a generalization of a procedure, we developed
in Golse/Klar [7], leading to easy to evaluate, however accurate coupling
conditions. These conditions yield a considerable improvement compared to
the ones obtained by equalizing moments or fluxes as in section 3.
We mention that in Golse [5] and Klar [11] coupling conditions for the coupling
of the Euler equations to the Boltzmann equation have been developed and
investigated that are also based on the analysis of a kinetic half space problem.
A different approach can be found in Illner/Neunzert [10].
Another remark concerns the situation, when the kinetic domain shrinks to 0
and only the drift diffusion equation is solved in the whole domain. In this
case the outgoing function is no longer needed and the condition on ae becomes
a standard slip boundary condition. An extensive treatment of such boundary
conditons for the gas dynamic case is given in the work of [18]. Here we
mention also the work of [20] on second order boundary conditions for the drift
diffusion equation. Using these boundary conditions one obtains a good quality
of approximation in the bulk of the device, but obviously no resolution of the
solution in the kinetic regions.
Finally, if higher order conditions are required, one may go to higher order in
the asymptotic expansion. However, this is very tedious and in most cases it
will not be necessary due to other limitations on the accuracy of the solution.
5 Numerical Results
In this section we investigate the coupling condition proposed in section 4 numerically
and compare them to the ones obtained by the equality of moments
or fluxes. The equality of moments or fluxes is straightforwardly implemented.
The coupling conditions based on the half space analysis are approximated by
the following:
+fflL a
- a (1; ff
where L a
x (k(\Delta)) is defined in (27) for arbitrary k and -
- a (1; v; t) is given by
formula (33).
The equation for the ingoing function for the Boltzmann region D 1 is explicitely
\Delta(r x ae(-x;
- a (0; ff
where A a
determined by (29) for arbitrary k and -
- a (0; ff
be approximated by 0.
In the following figures we calculated two examples. The solutions of the kinetic
and drift diffusion equations are determined by a scheme using a straightforward
explicit upwind discretization. To find the coupled solution we used, as
mentioned before, a Schwarz iteration procedure solving in turn the equations
in the kinetic and drift diffusion domain.
The kinetic equation is equation (1) in the slab D We have choosen
a relaxation model for the collision kernel,i.e. 1. The electric field E(x; t)
was choosen equal to The density of the
electrons is shown for a fixed time T such that a stationary state is obtained.
In

Figure

1 we consider an equilibrium situation. We take 3. The boundary
conditions are f 1 (0; v; The interface
is choosen at
We show the global kinetic solution, i.e. the reference solution, with mean
. Moreover the solution of the coupling problem is
shown, i.e., the solution of the kinetic equation with ffl 1 in D 1 and the solution
of the drift diffusion equation in D 2 . They are coupled together with the
Maxwell-Marshak condition using a first and second order approximation of the
distribution function in D 2 as explained in section 2. The boundary conditions
at for the drift diffusion equation are also found by the Maxwell Marshak
approximation at the boundary using an asymptotic expansion up to the
corresponding order.
Obviously the first order approximation is completely wrong as was to be ex-
pected. The other types of coupling conditions, i.e. the equality of moments
and the halfspace conditions, do in this case not differ very much from the
Maxwell Marshak conditions.11.41.82.2
x
KINETIC / DRIFT DIFFUSION
Kinetic
Second Order
First Order

Figure

1: Equilibrium situation - Maxwell-Marshak conditions
In

Figure

2 the same kinetic equation is shown but the boundary condition
at 3. The
interface between kinetic and aerodynamic domain is here located at
The mean free paths are ffl three types of coupling conditions
at the interface are shown. The boundary conditions at for the drift
diffusion equation are derived by the procedure corresponding to the one used
for the coupling condition. In this situation a nonequilibrium situation prevails
at the interface. The advantage of a more exact analysis at the interface is
clearly seen.11.41.82.2
x
KINETIC DRIFT DIFFUSION
Kinetic
Marshak
Halfspace
Moment

Figure

2: Nonequilibrium situation - Comparison of different coupling condition

In

Figure

3 the same as in Figure 2 is shown however we plotted only the region
around the interface.1.051.151.251.351.45
x
KINETIC DRIFT DIFFUSION
Kinetic
Marshak
Halfspace
Moment

Figure

3: Nonequilibrium situation - Interface region
6

Appendix

Approximate Solutions of Linear Half
Space Problems
In this section we derive approximations for the asymptotic values and outgoing
distributions of the halfspace problem (11) with arbitrary ingoing function and
of problem (13) with ingoing function 0.
Approximation of the solutions of (11):
The equation under consideration is
with
Z
In the following we will consider the above problem
at a fixed boundary point -
x with normal i.e. we consider
We proceed now with the extension of a scheme derived for simple model equations
in Golse/Klar [7].
Determination of the asymptotic value:
We consider a slightly more general equation introducing a velocity shift :
with e.g. Greenberg et al. [9] this equation has a unique
solution, if u is bigger than 0. If u is smaller than 0 it is necessary to fix the
to an arbitrary constant to obtain uniqueness of the solution. We assume u ? 0.
The adjoint equation to equation (20) is
Z
The reason for considering this equation will become clear in the following.
The condition on the flux is needed in this case to ensure uniqueness due to the
above considerations. We derive now an approximate series for /:
Shifting velocities using the definition of Q and
the fact that \Phi(\Gammav;
Z
Z
This equation will now be solved approximately by solving macroscopic equations
derived by a Chapman-Enskog type procedure: An approximate equation
for
Z
is obtained in the following way: Integrating (22) gives
@ y
Z
Z
A short calculation using the selfadjointness properties of the collision operator,
see [17], gives
Z
Z
Z
Z
Z
Approximating now / by its mean value, i.e. / - ae 1 M gives
Z
with D
R
R
Using (23) and (24) we get
This can be solved explicitely with two free parameters yielding
D u y +B:
The first approximation / 1 of the solution / of equation (22) is then obtained by
approximating /(w) in
R \Phi(v; w)/(w)dw by ae 1 M(w) and solving the resulting
equation
explicitely. A and B are determined by
Z
the closest analogue to
R
The / approximation can be iterated: From (22) and (25) we obtain the equation
for the remainder
\Gamma[
Z
Z
one can derive in a similar way as above an approximate
equation for ae 2 , namely
Z
Z
ae 2 is again uniquely determined up to 2 Parameters. The approximation / 2 of
given by the solution of
Z
The two parameters in the solution of (26) are determined by the condition
Z
Going on one obtains a series
Z
where
Z
and
Z
Assuming that this series converges, it can be easily seen by summing up the
's, that the series is equal to the solution / of the adjoint equation (22).
Hence one only has to transform v and u backwards, v ! \Gammav; u ! \Gammau to get
the desired approximation of (21), /(y; v).
Let ' be a solution of (20) and / one of (21) , then
Z
Z
Z
Z
Z
Z
Z
0:
Using this we get
Z
Z
and substituting gives
Z
Z
(v1 +u)?0
Or with
since '(1; v)M \Gamma1 (v) is independent of v
Z
(v1 +u)?0
where /(0; v) is approximated by the series
derived above. In this way one obtains an approximation of the asymptotic
value
Taking only the first step of the above procedure, i.e., and the limit for
u tending to 0 one obtains after some calculations the following approximation
for the asymptotic value L -
L a
R
R
Z
(v
R
R
dv
The first term in this formula is exactly the Maxwell-Marshak approximation.
This shows that actually the Maxwell-Marshak condition is a first approximation
of the asymptotic value of the solution of the halfspace problem.
Approximation of the outgoing distribution:
We consider equation (19) directly. Now the asymptotic value '(1; v) is known
from the above procedure. We derive for (19) in the same way as above the
macroscopic equation
The solution is ae some constant since ae 1 must be finite at infinity.
Then the first approximation ' 1 of the solution of (19) is obtained by solving
This gives
y
To obtain ' 1 (1; we have to put B equal to '(1; v)M \Gamma1 . Hence,
we get for
This can be approximated by the method derived above. To obtain in this
way the first nontrivial approximation of the outgoing distribution we have to
iterate the procedure once more: Considering the equation for
\Gamma[
Z
Z
one derives the second step gas dynamics equation:
Z
Z
As before approximated by the solution ' 2 of
Z
with Hence, solving this equation we obtain for
approximation for the outgoing distribution '(0; v) of (19).
Approximating the asymptotic value '(1; ff i ; v) of (18) by L a
x (k(\Delta)) yields in
this way the following approximation for the Albedo operator of equation (18):
A a
Z
x (k(\Delta))M (w)]dw:
Approximation of the higher order halfspace problem (13):
We consider finally the second relevant half space problem for -
at
the interface point -
namely
is the boundary layer
function defined in section 4 that solves (11) with ingoing function
rough approximation procedure for this equation will be given
in the following:
We start by approximating -(y; ff by the above procedure. Remarking
that -(1; ff is equal to 0 formula (28) gives
\Deltan
Putting this into equation (30) we obtain an approximate equation for -
With the same procedure as before we derive an approximate equation for
R
Z
with D as before, -
R
and the mean value
Z
Z
\Deltan
y g(-x; v; t)dv:
Solving (32) with ae(0; and the requirement of boundedness at infinity
gives ae(y; ff in particular the following approximation for -
- a (1; ff
Z
Z
(w
- a (0; ff is then the first approximation of -
Obviously here a more sophisticated approximation could be performed at the
expense of obtaining conditions that are more difficult to treat from a computational
point of view.

Acknowledgements

I am grateful to F. Golse, H. Neunzert and R. Wegener for helpful discussions.



--R

Diffusion Approximation and Computation of the Critical Size
Boundary Layers and Homogenization of Transport Processes
Numerical Coupling of Non-conservative or Kinetic Models with the Conservative Compressible Navier-Stokes Equations
The Boltzmann Equation and its Applications
Applications of the Boltzmann Equation Within the Context of Upper Atmosphere Vehicle Aerodynamics
Knudsen Layers from a Computational Viewpoint
A Numerical Method for Computing Asymptotic States and Outgoing Distributions for Kinetic Linear Half Space Problems


Linking Kinetic and Aerodynamic Descriptions

Convergence of Alternating
Initial and Boundary Conditions for Equilibrium Diffusion Theory
Approximate Method in the Kinetic Theory in Phys
Report for the Hermes Project DPH 6174/91

Diffusion Approximation of the Linear Semiconductor Equa- tion
Analysis and Simulation of Semiconductor Devices
Asymtotic theory of a Steady Flow of a Rarefied Gas Past Bodies for Small Knudsen Numbers
Second order boundary conditions of Drift Diffusion Equations of Semiconductors
--TR
kinetic semiconductor equations;asymptotic analysis;domain decomposition methods;drift-diffusion equation