An AsymptoticPreserving Scheme for the Semiconductor Boltzmann Equation toward the EnergyTransport Limit
Abstract
We design an asymptoticpreserving scheme for the semiconductor Boltzmann equation which leads to an energytransport system for electron mass and internal energy as mean free path goes to zero. To overcome the stiffness induced by the convection terms, we adopt an evenodd decomposition to formulate the equation into a diffusive relaxation system. New difficulty arises in the twoscale stiff collision terms, whereas the simple BGK penalization does not work well to drive the solution to the correct limit. We propose a clever variant of it by introducing a threshold on the stiffer collision term such that the evolution of the solution resembles a Hilbert expansion at the continuous level. Formal asymptotic analysis and numerical results are presented to illustrate the efficiency and accuracy of the new scheme.
Key words. semiconductor Boltzmann equation, energytransport system, asymptoticpreserving scheme, fast spectral method.
AMS subject classifications. 82D37, 35Q20, 65N06, 65N12, 65N35.
1 Introduction
The semiconductor Boltzmann equation describes the transport of charge carriers in semiconductor devices. It is derived following a statistical approach by incorporating the quantum mechanical effects semiclassically, thus provides accurate description of the physics at the kinetic level [31, 8, 26]. A dimensionless form of this equation usually contains small parameters such as the mean free path or time. Besides the high dimensionality of the probability distribution function, the presence of these small parameters poses tremendous computational challenge since one has to numerically resolve the small scales.
To save the computational cost, in the past decades, various macroscopic models were derived from the Boltzmann equation based on assumptions of different dominating effects. One of the well accepted model, for example, is the driftdiffusion equation which consists of mass continuity equation for electrons (or holes) [33, 17]. Ideally, if the parameters in the kinetic equation are uniformly small in the entire domain of interest, then the macroscopic models suffice to describe the physical phenomena, and it is more efficient to just solve them [28, 9, 35]. In practical applications, however, the validity of these models may break down in part of the domain (parameters are not small anymore), and one is forced to solve the kinetic equation which contains mesoscopic information [5, 6]. A natural solution to this situation is a domain decomposition approach [4, 29], but finding the interface condition connecting the kinetic and macroscopic equations is a highly nontrivial task. Another line of research is to find a unified scheme for the kinetic equation such that when the small parameter goes to zero, it automatically becomes a macroscopic solver. This designing concept leads to the asymptoticpreserving (AP) scheme [21], which was first introduced by S. Jin for transport equations in diffusive regimes [20]. In the semiconductor framework, an initial effort toward the AP schemes was proposed in [22] for the linear Boltzmann equation with an anisotropic collision term, whose computation was further improved in [11]. Recently a higher order scheme with a less strict stability condition was constructed in the sense that the parabolic CFL constraint is relaxed to a hyperbolic one [13]. All these works consider a linear collision operator with smooth kernel which uniquely defines an equilibrium state. As a result, the corresponding macroscopic equation is in the form of a driftdiffusion equation. Although this equation gives satisfactory simulation results for semiconductor devices on the micrometer scale, it is not able to capture the hotelectron effects in submicron devices [27]. High field scaling deals with this problem to some extent, but it only works for the situation where the field effect is strong enough to balance the collision [24, 7].
In this work, we are interested in a more realistic semiconductor Boltzmann equation [1, 10]. By considering the elastic collision as dominant, and electronelectron correlation as subdominant effects, one can pass on the asymptotic limit to obtain an energytransport (ET) model. It consists of a system of conservation laws for mass and internal energy of charge carriers with fluxes computed through a constitutive relation [25]. To design an AP scheme for such kinetic equation, we face twofold challenge: 1. the convection terms are stiff; 2. two stiff collision terms live on different scales. The convection terms can be treated by an evenodd decomposition as in [23, 22]. For the collision terms, due to their complicated forms, we choose to penalize them with a suitable BGK operator [15]. However, unlike the usual collision operator with smoothed kernel, the leading elastic operator has non unique null space (the kernel is a Delta function). Only when the electronelectron operator in next level takes into effect, the solution can be eventually driven to a fixed FermiDirac distribution. A closer examination of the asymptotic behavior of the solution reveals that the penalization should be performed wisely, otherwise it won’t capture the correct limit. To this end, we propose a thresholded penalization scheme. Simply speaking, when the threshold is satisfied, we turn off the leading order mechanism and move to the next order, which in some sense resembles the Hilbert expansion in the continuous case. We will show that this new scheme, under certain assumptions, satisfies the following four properties ( denotes the small parameter; , , and are the time step and mesh size in spatial and wave vector (momentum) domain):

For fixed , it is consistent to the Boltzmann equation when .

For fixed , it becomes a discretization to the limiting ET system when .

It is uniformly stable for a wide range of , from to .

Implicit terms can be implemented explicitly (free of Newtontype solvers).
An important ingredient in this AP scheme is the accurate numerical solvers for the collision operators. Since the electronelectron operator falls into a special case of the quantum Boltzmann operator, we adopt the fast spectral method developed in [19]. For the elastic collision, it is desirable to evaluate it in the same spectral framework but the direct computation would be very expensive. We propose a new fast method by exploring the lowrank structure in the coefficient matrix.
The rest of the paper is organized as follows. In the next section we give a brief review of the scalings of the semiconductor Boltzmann equation and the derivation of the ET model through a systematic approach. Section 3 is devoted to a detailed description of our schemes. We will present it in a pedagogical way that the spatially homogeneous case is considered first with an emphasis on the twoscale stiff collision terms, and then embrace the spatial dependence to treat the full problem. In either case, the asymptotic property of the numerical solution is carefully analyzed. The spectral methods for computing the collision operators are gathered at the end. In Section 4 we give several numerical examples including the simulation of a 1D silicon diode to illustrate the efficiency, accuracy, and asymptotic properties of the scheme. Finally, the paper is concluded in Section 5.
2 The semiconductor Boltzmann equation and its energytransport limit
The Boltzmann transport equation that describes the evolution of electrons in the conduction band of a semiconductor reads [31, 8, 26]
\hb@xt@.01(2.1) 
where is the electron distribution function of position , wave vector , and time . is the reduced Planck constant, and is the positive elementary charge. The first Brillouin zone is the primitive cell in the reciprocal lattice of the crystal. For simplicity, we will restrict to the parabolic band approximation, where can be extended to the whole space , and the energyband diagram is given by
where is the effective mass of electrons.
In principle, the electrostatic potential is produced selfconsistently by the electron density with a fixed ion background of doping profile through the Poisson equation:
\hb@xt@.01(2.2) 
where
is the electron density (the spin degeneracy , with being the spin of electrons). and are the vacuum and the relative material permittivities. The doping profile takes into account the impurities due to acceptor and donor ions in the semiconductor device.
The collision operator explains three different effects:
where and account for the interactions between electrons and the lattice defects caused by ionized impurities and crystal vibrations (also called phonons); describes the correlations between electrons themselves. Specifically,
where is the Dirac measure, , , , , , are short notations for , , , , , and respectively. is the phonon energy, and is the phonon occupation number:
where is the Boltzmann constant and is the lattice temperature. The scattering matrices and are symmetric in and :
is symmetric pairwisely for four variables:
They are all determined by the underlying interaction laws.
2.1 Nondimensionalization of the Boltzmann equation
In order to nondimensionalize the Boltzmann equation (LABEL:eqn:_semiB1), we introduce the following typical values:
typical density;  
typical kinetic energy of electrons, is the applied bias;  
typical norm of wave vector such that ;  
typical distribution function scale;  
typical time scale;  
typical velocity scale;  
typical length scale;  
typical values of transition rates , , , 
and define also the dimensionless parameters:
After performing a change of variables as
(LABEL:eqn:_semiB1) becomes
where we have dropped the tildes without ambiguity. The collision operators take the following dimensionless form
\hb@xt@.01(2.3) 
where , , , and
The energyband diagram is now simply
\hb@xt@.01(2.4) 
The Poisson equation (LABEL:eqn:_Poisson1) becomes
\hb@xt@.01(2.5) 
where is the square of the scaled Debye length, and
2.2 Elastic approximation of the electronphonon interactions
We are interested in a high energy scale [2, 1], at which the relative energy gain or loss of electron energy during a phonon collision is very small, i.e.,
In addition, we assume that
which means that at the high energy scale, the phonon energy and the lattice thermal energy are considered as the same order of magnitude, and much smaller compared with the electron energy .
Treating as small parameter, one can expand the electronphonon collision operator as
where the first term is the elastic approximation and the second term is the inelastic correction. Therefore, the total collision operator can be recast as
with
\hb@xt@.01(2.6) 
and
Since , it is reasonable to assume and are both of (refer to [34] for physical data). However, it is much delicate to estimate the electronelectron collision frequency as it depends on the distribution function itself. Following the discussion in [10], we assign this term . That is, the electronelectron interactions are not as strong as elastic collisions, yet their density is not small enough to be safely neglected.
The final form of the scaled Boltzmann equation is thus
\hb@xt@.01(2.7) 
where , are given by (LABEL:Qel) and (LABEL:Qee) respectively (the exact form of will not be needed in the following discussion and is thereby omitted).
2.3 Diffusive regime and the energytransport limit
To derive a macroscopic model, we consider the time and length scales in a diffusive regime: , , then equation (LABEL:scaledB) rewrites as
\hb@xt@.01(2.8) 
which is the main kinetic equation we are going to study for the rest of the paper. This subsection is devoted to a formal derivation of the asymptotic limit of (LABEL:kinetic1) as . Our approach, following that of [10], is a combination of the Hilbert expansion and the moment method. To this end, we first list the required properties of the collision operators and . These will also be useful in designing numerical schemes.
Proposition 2.1
[1]

For any “regular” test function ,
In particular, for any ,

is a selfadjoint, nonpositive operator on .

The null space of is given by

The orthogonal complement of is
where the integral is defined through the “coarea formula” [14]: for any “regular” functions and , it holds
\hb@xt@.01(2.9) Here denotes the surface of constant energy , is the coarea measure, and is the energy densityofstates. Under the parabolic band approximation (LABEL:para), (LABEL:coarea) is just a spherical transformation:
and .

The range of is . The operator is invertible as an operator from onto . Its inverse is denoted by .

For any , we have and .
Proposition 2.2
[2]

For any “regular” test function ,
In particular, we have the conservation of mass and energy

Htheorem: let , then , and if ,
where
\hb@xt@.01(2.10) is the FermiDirac distribution function [32]. The variables and are the fugacity and (electron) temperature. Alternatively, can be defined in terms of the chemical potential and .
We also need the properties of the operator , an energy space counterpart of : for any ,
Proposition 2.3
[1]

Conservation of mass and energy

Htheorem: let , then , and
Now that the mathematical preliminaries are set up, we are ready to derive the macroscopic limit. The main result is summarized in the following theorem.
Theorem 2.4
[10] In equation (LABEL:kinetic1), when , the solution formally tends to a FermiDirac distribution function (LABEL:FD), with the position and time dependent fugacity and temperature satisfying the socalled EnergyTransport (ET) model:
\hb@xt@.01(2.11) 
where the density and energy are defined as
\hb@xt@.01(2.12) 
the fluxes and are given by
\hb@xt@.01(2.13) 
with the diffusion matrices
\hb@xt@.01(2.14) 
and the energy relaxation operator is
\hb@xt@.01(2.15) 
Proof. Inserting the Hilbert expansion into equation (LABEL:kinetic1) and collecting equal powers of leads to
\hb@xt@.01(2.16)  
\hb@xt@.01(2.17) 
From (LABEL:order2) and Proposition LABEL:prop:_qel (3), we know there exists some function such that
Plugging into (LABEL:order1):
\hb@xt@.01(2.18) 
The solvability condition for (Proposition LABEL:prop:_qel (4–5)) implies
Clearly the left hand side of above equation is equal to zero ( is odd in ), so
By Proposition LABEL:prop:_qeeb (2), we know that is a FermiDirac distribution with position and time dependent and . Therefore, itself is equal to zero by Proposition LABEL:prop:_qee (2), and (LABEL:order11) reduces to
Then using Proposition LABEL:prop:_qel (5–6), we have
\hb@xt@.01(2.19) 
Going back to the original equation (LABEL:kinetic1), multiplying both sides by and integrating w.r.t. gives:
\hb@xt@.01(2.20) 
Terms involving and vanish due to Propositions LABEL:prop:_qel (1) and LABEL:prop:_qee (1). Recall that is the difference between and an elastic operator, and both of them can be easily seen to conserve mass, so . Thus (LABEL:ET1) simplifies to
\hb@xt@.01(2.21) 
where .
From the previous discussion, we know with being the FermiDirac distribution (LABEL:FD), and given by (LABEL:f1). Plugging into (LABEL:ET2), to the leading order we have ( term drops out since is an even function in ):
\hb@xt@.01(2.22) 
Utilizing the special form of , one can rewrite as
Then a simple manipulation of (LABEL:ET21) yields the ET system (LABEL:ET3).
The ET model (LABEL:ET3) is widely used in practical and industrial applications (see [27] for a review and references therein). It can also be derived from the Boltzmann equation through the socalled SHE (Spherical Harmonics Expansion) model [16]. In either case, the rigorous theory behind the formal limit is still an open question.
Remark 2.5
If the electronelectron interaction is assumed as one of the dominant terms in (LABEL:kinetic1), i.e., the same order as elastic collision (which could be the physically relevant situation for very dense electrons), one can still derive the ET model (LABEL:ET3) via a similar procedure [2], but the diffusion coefficients are different. A rigorous result is available in this case [3].
Remark 2.6
Remark 2.7
Unlike the classical statistics, given macroscopic variables and (LABEL:rhoe), finding the corresponding FermiDirac distribution (LABEL:FD) is not a trivial issue. Under the parabolic approximation (LABEL:para), and are related to and via [18]
\hb@xt@.01(2.23) 
where is the FermiDirac function of order
\hb@xt@.01(2.24) 
and is the Gamma function. For small (), the integrand in (LABEL:fermi) can be expanded in powers of :
Thus, when , behaves like itself and one recovers the classical limit.
3 Asymptoticpreserving (AP) schemes for the semiconductor Boltzmann equation
Equation (LABEL:kinetic1) contains three different scales. To overcome the stiffness induced by the and terms, a fully implicit scheme would be desirable. However, neither the collision operators nor the convection terms are easy to solve implicitly. Our goal is to design an appropriate numerical scheme that is uniformly stable in both kinetic and diffusive regimes, i.e., works for all values of ranging from to , while the implicit terms can be treated explicitly.
We will first consider a spatially homogeneous case with an emphasis on the collision operators, and then include the spatial dependence to treat the convection terms. To facilitate the presentation, we always make the following assumptions without further notice:

The inelastic collision operator in (LABEL:kinetic1) is assumed to be zero, since it is the weakest effect and its appearance won’t bring extra difficulties to numerical schemes.

The scattering matrices and are rotationally invariant:
Then it is not difficult to verify that (see Proposition LABEL:prop:_qel (4))
\hb@xt@.01(3.1) 
where
and is the mean value of over sphere :
In particular, for any odd function ,
This observation is crucial in designing our AP schemes.
3.1 The spatially homogeneous case
In the spatially homogeneous case, equation (LABEL:kinetic1) reduces to
\hb@xt@.01(3.2) 
where only depends on and . An explicit discretization of (LABEL:homo), e.g., the forward Euler scheme, suffers from severe stability constraints: has to be smaller than . Implicit schemes do not have such a restriction, but require some sort of iteration solvers for and which can be quite complicated.
To tackle these two stiff terms, we adopt the penalization idea in [15], i.e., penalize both and by their corresponding “BGK” operators:
\hb@xt@.01(3.3) 
Using the properties of and from the last section, can be naturally chosen as the FermiDirac distribution (in the homogeneous case and are conserved, so is an absolute Maxwellian and can be obtained from the initial condition), whereas can in principle be any function of such that and share the same density and energy (the choice of is not essential as we shall see, and we will get back to this when we consider the spatially inhomogeneous case). As the goal of penalization is to make the residue as small as possible so that it is nonstiff or less stiff, and , can be expressed symbolically as
we hence choose
\hb@xt@.01(3.4) 
Other choices are also possible [15]. Generally speaking, we only need the coefficient to be a rough estimate of the Frechet derivative of the collision operator around equilibrium.
Therefore, a firstorder IMEX scheme for (LABEL:penalization) is written as
\hb@xt@.01(3.5)  
3.1.1 Asymptotic properties of the numerical solution
To better understand the asymptotic behavior of the numerical solution, in this subsection we assume that . Then scheme (LABEL:homoscheme) becomes