An Asymptotic-Preserving Scheme for the Semiconductor Boltzmann Equation toward the Energy-Transport Limit

# An Asymptotic-Preserving Scheme for the Semiconductor Boltzmann Equation toward the Energy-Transport Limit

Jingwei Hu Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin, 201 East 24th St, Stop C0200, Austin, TX 78712 (hu@ices.utexas.edu).    Li Wang Department of Mathematics, University of California, Los Angeles, 520 Portola Plaza, Los Angeles, CA 90095 (liwang@math.ucla.edu).
###### Abstract

We design an asymptotic-preserving scheme for the semiconductor Boltzmann equation which leads to an energy-transport 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 even-odd decomposition to formulate the equation into a diffusive relaxation system. New difficulty arises in the two-scale 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, energy-transport system, asymptotic-preserving 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 drift-diffusion 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 non-trivial 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 asymptotic-preserving (AP) scheme , which was first introduced by S. Jin for transport equations in diffusive regimes . In the semiconductor framework, an initial effort toward the AP schemes was proposed in  for the linear Boltzmann equation with an anisotropic collision term, whose computation was further improved in . 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 . 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 drift-diffusion equation. Although this equation gives satisfactory simulation results for semiconductor devices on the micrometer scale, it is not able to capture the hot-electron effects in submicron devices . 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 electron-electron correlation as sub-dominant effects, one can pass on the asymptotic limit to obtain an energy-transport (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 . To design an AP scheme for such kinetic equation, we face two-fold challenge: 1. the convection terms are stiff; 2. two stiff collision terms live on different scales. The convection terms can be treated by an even-odd decomposition as in [23, 22]. For the collision terms, due to their complicated forms, we choose to penalize them with a suitable BGK operator . 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 electron-electron operator in next level takes into effect, the solution can be eventually driven to a fixed Fermi-Dirac 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):

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

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

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

4. Implicit terms can be implemented explicitly (free of Newton-type solvers).

An important ingredient in this AP scheme is the accurate numerical solvers for the collision operators. Since the electron-electron operator falls into a special case of the quantum Boltzmann operator, we adopt the fast spectral method developed in . 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 low-rank 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 two-scale 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 1-D 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 energy-transport limit

The Boltzmann transport equation that describes the evolution of electrons in the conduction band of a semiconductor reads [31, 8, 26]

 ∂tf+1ℏ∇kε(k)⋅∇xf+qℏ∇xV(x,t)⋅∇kf=Q(f),x∈Ω⊂Rd,  k∈B⊂Rd,  d=2,3, \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 energy-band diagram is given by

 ε(k)=ℏ22m∗|k|2,

where is the effective mass of electrons.

In principle, the electrostatic potential is produced self-consistently by the electron density with a fixed ion background of doping profile through the Poisson equation:

 ϵ0∇x(ϵr(x)∇xV(x,t))=q(ρ(x,t)−h(x)), \hb@xt@.01(2.2)

where

 ρ(x,t):=∫Rdf(x,k,t)g(2π)ddk

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:

 Q=Qimp+Qph+Qee,

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,

 Qimp(f)(k)=∫Rdϕ% imp(k,k′)δ(ε′−ε)(f′−f)dk′, Qph(f)(k)=∫Rdϕph% (k,k′){[(Nph+1)δ(ε−ε′+εph)+Nphδ(ε−ε′−εph)]f′(1−f) −[(Nph+1)δ(ε′−ε+εph)+Nphδ(ε′−ε−εph)]f(1−f′)}dk′, Qee(f)(k)=∫R3dϕ% ee(k,k1,k′,k′1)δ(ε′+ε′1−ε−ε1)δ(k′+k′1−k−k1) ×[f′f′1(1−f)(1−f1)−ff1(1−f′)(1−f′1)]dk1dk′dk′1,

where is the Dirac measure, , , , , , are short notations for , , , , , and respectively. is the phonon energy, and is the phonon occupation number:

 Nph=1eεphkBTL−1,

where is the Boltzmann constant and is the lattice temperature. The scattering matrices and are symmetric in and :

 ϕimp(k,k′)=ϕimp(k′,k),ϕph(k,k′)=ϕph(k′,k);

is symmetric pair-wisely for four variables:

 ϕee(k,k1,k′,k′1)=ϕee(k1,k,k′,k′1)=ϕee(k′,k′1,k,k1).

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:

 ρ0: typical density; ε0=qV0: typical kinetic energy of electrons, V0 is the applied bias; k0: typical norm of wave vector k such that ε(k0)=ε0/2; f0=(2π)dρ0gkd0=η: typical distribution function scale; t0: typical time scale; v0=ε0ℏk0: typical velocity scale; x0=t0v0: typical length scale; ϕimp,0,ϕph,0,ϕee,0: typical values of transition rates ϕimp(k,k′), ϕph(k,k′), ϕee(k,k1,k′,k′1),

and define also the dimensionless parameters:

 α2=εphε0, γ2=kBTLε0, νimp=ϕimp,0kd0t0ε0, νph=ϕph,0kd0t0ε0, νee=ϕee,0kd0t0ε0(2π)dρ0g.

After performing a change of variables as

 ~f=ff0,~t=tt0,~x=xx0,~k=kk0,~ε=εε0,~V=VV0,~ϕ∙=ϕ∙ϕ∙,0,

(LABEL:eqn:_semiB1) becomes

 ∂tf+∇kε⋅∇xf+∇xV⋅∇kf=Qimp(f)+Qph(f)+Qee% (f),

where we have dropped the tildes without ambiguity. The collision operators take the following dimensionless form

 Qimp(f)(k)=∫RdΦ% imp(k,k′)δ(ε′−ε)(f′−f)dk′, Qph(f)(k)=∫RdΦph% (k,k′){[(Nph+1)δ(ε−ε′+α2)+Nphδ(ε−ε′−α2)]f′(1−ηf) −[(Nph+1)δ(ε′−ε+α2)+Nphδ(ε′−ε−α2)]f(1−ηf′)}dk′, Qee(f)(k)=∫R3dΦ% ee(k,k1,k′,k′1)δ(ε′+ε′1−ε−ε1)δ(k′+k′1−k−k1) ×[f′f′1(1−ηf)(1−ηf1)−ff1(1−ηf′)(1−ηf′1)]dk1dk′dk′1, \hb@xt@.01(2.3)

where , , , and

 Nph=1eα2/γ2−1.

The energy-band diagram is now simply

 ε(k)=12|k|2. \hb@xt@.01(2.4)

The Poisson equation (LABEL:eqn:_Poisson1) becomes

 C0∇x(ϵr(x)∇xV(x,t))=ρ(x,t)−h(x), \hb@xt@.01(2.5)

where is the square of the scaled Debye length, and

 ρ(x,t)=∫Rdf(x,k,t)dk.

### 2.2 Elastic approximation of the electron-phonon 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.,

 α≪1.

 αγ∼O(1),

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 electron-phonon collision operator as

 Qph(f)(k)=∫Rd(2Nph+1)Φph(k,k′)δ(ε′−ε)(f′−f)dk′+α2Qinelph(f)(k),

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

 Q(f)=Qel(f)+Qee(f)+α2Qinelph(f),

with

 Qel(f)(k)=∫RdΦel(k,k′)δ(ε′−ε)(f′−f)dk′, \hb@xt@.01(2.6)

and

 Φel(k,k′)=Φimp(k,k′)+(2Nph+1)Φph(k,k′).

Since , it is reasonable to assume and are both of (refer to  for physical data). However, it is much delicate to estimate the electron-electron collision frequency as it depends on the distribution function itself. Following the discussion in , we assign this term . That is, the electron-electron 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

 ∂tf+∇kε⋅∇xf+∇xV⋅∇kf=Qel(f)+αQee(f)+α2Qinelph(f), \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 energy-transport limit

To derive a macroscopic model, we consider the time and length scales in a diffusive regime: , , then equation (LABEL:scaledB) rewrites as

 ∂tf+1α(∇kε⋅∇xf+∇xV⋅∇kf)=1α2Qel(f)+1αQee(f)+Q% inelph(f), \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 , 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 ,

 ∫RdQel(f)g(ε)dk=0.
2. is a self-adjoint, non-positive operator on .

3. The null space of is given by

 N(Qel)={f(ε(k)),∀f}.
4. The orthogonal complement of is

 N(Qel)⊥={g(k) | ∫Sεg(k)dNε(k)=0,∀ε∈R},

where the integral is defined through the “coarea formula” : for any “regular” functions and , it holds

 ∫Rdf(k)dk=∫R(∫Sεf(k)dNε(k))dε. \hb@xt@.01(2.9)

Here denotes the surface of constant energy , is the coarea measure, and is the energy density-of-states. Under the parabolic band approximation (LABEL:para), (LABEL:coarea) is just a spherical transformation:

 ∫Rdf(k)dk=∫∞0(∫Sd−1f(|k|σ)|k|d−2dσ)dε,

and .

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

6. For any , we have and .

###### Proposition 2.2



1. For any “regular” test function ,

 ∫Rd Qee(f)gdk=−14∫R4dΦee(k,k1,k′,k′1)δ(ε′+ε′1−ε−ε1)δ(k′+k′1−k−k1) ×[f′f′1(1−ηf)(1−ηf1)−ff1(1−ηf′)(1−ηf′1)](g′+g′1−g−g1)dkdk1dk′dk′1.

In particular, we have the conservation of mass and energy

 ∫RdQee(f)dk=∫RdQee(f)εdk=0.
2. H-theorem: let , then , and if ,

 ∫RdQee(f)H(f)dk=0⟺f=M(ε(k))⟺Qee(f)=0,

where

 M(ε(k))=1η1z−1eε(k)T+1 \hb@xt@.01(2.10)

is the Fermi-Dirac distribution function . 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 ,

 ¯¯¯¯¯¯¯¯¯Qee(F)(ε):=∫SεQee(F)(k)dNε(k).
###### Proposition 2.3



1. Conservation of mass and energy

 ∫R¯¯¯¯¯¯¯¯¯Qee(F)dε=∫R¯¯¯¯¯¯¯¯¯Qee(F)εdε=0.
2. H-theorem: let , then , and

 ∫R¯¯¯¯¯¯¯¯¯Qee(F)H(F)dε=0⟺F=M(ε)⟺¯¯¯¯¯¯¯¯¯Qee(F)=0.

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

 In equation (LABEL:kinetic1), when , the solution formally tends to a Fermi-Dirac distribution function (LABEL:FD), with the position and time dependent fugacity and temperature satisfying the so-called Energy-Transport (ET) model:

 ∂t(ρρE)+(∇x⋅jρ∇x⋅jE)−(0∇xV⋅jρ)=(0Winelph), \hb@xt@.01(2.11)

where the density and energy are defined as

 ρ(z,T)=∫Rdfdk=∫RdMdk,E(z,T)=1ρ∫Rdfεdk=1ρ∫RdMεdk; \hb@xt@.01(2.12)

the fluxes and are given by

 (jρ(z,T)jE(z,T))=−(D11D12D21D22)⎛⎜⎝∇xzz−∇xVT∇xTT2⎞⎟⎠ \hb@xt@.01(2.13)

with the diffusion matrices

 Dij=∫Rd∇kε⊗Q−1el(−∇kε)M(1−ηM)εi+j−2dk; \hb@xt@.01(2.14)

and the energy relaxation operator is

 Winelph(z,T)=∫RdQinelph(M)εdk. \hb@xt@.01(2.15)

Proof. Inserting the Hilbert expansion into equation (LABEL:kinetic1) and collecting equal powers of leads to

 O(α−2):Qel(f0)=0, \hb@xt@.01(2.16) O(α−1):Qel(f1)=∇kε⋅∇xf0+∇xV⋅∇kf0−Qee(f0). \hb@xt@.01(2.17)

From (LABEL:order-2) and Proposition LABEL:prop:_qel (3), we know there exists some function such that

 f0(x,k,t)=F(x,ε(k),t).

Plugging into (LABEL:order-1):

 Qel(f1)=∇kε⋅(∇xF+∇xV∂εF)−Qee(F). \hb@xt@.01(2.18)

The solvability condition for (Proposition LABEL:prop:_qel (4–5)) implies

 ∫Sε∇kε⋅(∇xF+∇xV∂εF)dNε(k)=∫SεQee(F)dNε(k).

Clearly the left hand side of above equation is equal to zero ( is odd in ), so

 ∫SεQee(F)dNε(k)=¯¯¯¯¯¯¯¯¯Qee(F)=0.

By Proposition LABEL:prop:_qeeb (2), we know that is a Fermi-Dirac distribution with position and time dependent and . Therefore, itself is equal to zero by Proposition LABEL:prop:_qee (2), and (LABEL:order-11) reduces to

 Qel(f1)=∇kε⋅(∇xM+∇xV∂εM).

Then using Proposition LABEL:prop:_qel (5–6), we have

 f1=−Q−1el(−∇kε)⋅(∇xM+∇xV∂εM). \hb@xt@.01(2.19)

Going back to the original equation (LABEL:kinetic1), multiplying both sides by and integrating w.r.t. gives:

 =∫Rd[1α2Qel(f)+1αQee(f)+Qinelph% (f)](1ε)dk. \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

 ∂t(ρρE)+∫Rd1α(∇kε⋅∇xf+∇xV⋅∇kf)(1ε)dk=(0Winelph(f)), \hb@xt@.01(2.21)

where .

From the previous discussion, we know with being the Fermi-Dirac 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 ):

 ∂t(ρρE)+∫Rd(∇kε⋅∇xf1+∇xV⋅∇kf1)(1ε)dk=(0Winelph(M)). \hb@xt@.01(2.22)

Utilizing the special form of , one can rewrite as

 f1=−Q−1el(−∇kε)⋅(∇xzz−∇xVT+ε∇xTT2)M(1−ηM).

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  for a review and references therein). It can also be derived from the Boltzmann equation through the so-called SHE (Spherical Harmonics Expansion) model . In either case, the rigorous theory behind the formal limit is still an open question.

###### Remark 2.5

If the electron-electron 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 , but the diffusion coefficients are different. A rigorous result is available in this case .

###### Remark 2.6

If we consider even longer time scale such that the energy relaxation term equilibrates to the lattice temperature , and further assume that is the classical Maxwellian, then (LABEL:ET3) reduces to a single equation

 ∂tρ+∇x⋅[−D11(∇xρρ−∇xVTL)]=0.

This is the classical drift-diffusion model [36, 33].

###### Remark 2.7

Unlike the classical statistics, given macroscopic variables and (LABEL:rhoe), finding the corresponding Fermi-Dirac distribution (LABEL:FD) is not a trivial issue. Under the parabolic approximation (LABEL:para), and are related to and via 

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩ρ=(2πT)d2ηFd2(z),E=d2TFd2+1(z)Fd2(z), \hb@xt@.01(2.23)

where is the Fermi-Dirac function of order

 Fν(z)=1Γ(ν)∫∞0xν−1z−1ex+1dx,0

and is the Gamma function. For small (), the integrand in (LABEL:fermi) can be expanded in powers of :

 Fν(z)=∞∑n=1(−1)n+1znnν=z−z22ν+z33ν−….

Thus, when , behaves like itself and one recovers the classical limit.

## 3 Asymptotic-preserving (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:

1. 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.

2. The scattering matrices and are rotationally invariant:

 Φel(k,k′)=Φel(|k|,|k′|),Φee(k,k1,k′,k′1)=Φee(|k|,|k1|,|k′|,|k′1|).

Then it is not difficult to verify that (see Proposition LABEL:prop:_qel (4))

 Qel(f)(k)=λel(ε)([f](ε)−f(k)), \hb@xt@.01(3.1)

where

 λel(ε(k)):=∫RdΦel(k,k′)δ(ε′−ε)dk′=Φel% (|k|,|k|)N(ε),

and is the mean value of over sphere :

 [f](ε(k)):=1N(ε)∫Sεf(k)dNε(k).

In particular, for any odd function ,

 Qel(f)(k)=−λel(ε)f(k),Q−1el(f(k))=−1λel(ε)f(k).

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

 ∂tf=1α2Qel(f)+1αQee(f), \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 , i.e., penalize both and by their corresponding “BGK” operators:

 ∂tf=Qel(f)−βel(Mel−f)α2+βel(Mel−f)α2+Qee(f)−βee(Mee−f)α+βee(Mee−f)α. \hb@xt@.01(3.3)

Using the properties of and from the last section, can be naturally chosen as the Fermi-Dirac 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 non-stiff or less stiff, and , can be expressed symbolically as

 Qel(f)(k)=Q+el(f)(ε)−λel(ε)f(k);Qee(f)(k)=Q+ee(f)(k)−Q−ee(f)(k)f(k),

we hence choose

 βel≈maxελel(ε);βee≈maxkQ−ee(f)(k). \hb@xt@.01(3.4)

Other choices are also possible . Generally speaking, we only need the coefficient to be a rough estimate of the Frechet derivative of the collision operator around equilibrium.

Therefore, a first-order IMEX scheme for (LABEL:penalization) is written as

 fn+1−fnΔt= Qel(fn)−βel(M% el−fn)α2+βel(Mel−fn+1)α2 \hb@xt@.01(3.5) +Qee(fn)−βee(M−fn)α+βee(M−fn+1)α.

#### 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

 fn+