AN ASYMPTOTIC-PRESERVING SCHEME FOR THE SEMICONDUCTOR BOLTZMANN EQUATION WITH TWO-SCALE COLLISIONS: A SPLITTING APPROACH

. We present a new asymptotic-preserving scheme for the semiconductor Boltzmann equation with two-scale collisions — a leading-order elastic collision together with a lower-order interparticle collision. When the mean free path is small, numerically solving this equation is prohibitively expensive due to the stiﬀ collision terms. Furthermore, since the equilibrium solution is a (zero-momentum) Fermi-Dirac distribution resulting from joint action of both collisions, the simple BGK penalization designed for the one-scale collision [10] cannot capture the correct energy-transport limit. This problem was addressed in [13], where a thresholded BGK penalization was introduced. Here we propose an alternative based on a splitting approach . It has the advantage of treating the collisions at diﬀerent scales separately, hence is free of choosing threshold and easier to implement. Formal asymptotic analysis and numerical results validate the eﬃciency and accuracy of the proposed scheme.

1. Introduction.The semiconductor Boltzmann equation describes the transport of charge carriers (electrons or holes) in semiconductor devices [20,6,18].In this paper, we are interested in the following non-dimensionalized form [2,1,7]: where f (x, k, t) of position x, wave vector k and time t, is the electron distribution function in the conduction band of a semiconductor.The parameter α is the scaled mean-free path which can be either small (diffusive regime) or large (kinetic regime) depending on the device structure.ε(k) is the energy-band diagram given by (the parabolic band approximation is assumed for simplicity) V (x, t) is the electrostatic potential produced self-consistently through the Poisson equation: where D 0 is the square of scaled Debye length, ρ(x, t) = R d f (x, k, t) dk is the electron density, and h(x) is the fixed doping profile that takes into account the impurities due to acceptor and donor ions in the device.The right hand side of equation ( 1) models three different collision effects: 1. Q el is the elastic collision combining interactions between electrons and lattice impurities, and the elastic part of electron-phonon interactions; 2. Q ee is the interactions between electrons themselves; 3. Q inel ph is the remaining inelastic part of electron-phonon interactions.Specifically, Q el and Q ee are given as follows (the exact form of Q inel ph will not be needed in the following discussion and thus is omitted): where δ is the Dirac measure, ε , ε 1 , ε 1 , f , f 1 , f 1 are short notations for ε(k ), ε(k 1 ), ε(k 1 ), f (x, k , t), and etc, η is the typical distribution function scale characterizing the degree of degeneracy of the system.The scattering matrices Φ el and Φ ee are determined by the underlying interaction laws.They satisfy the symmetry conditions: To this date, there are numerous models describing the electron flow through a semiconductor device, ranging from macroscopic equations to kinetic equations, or even microscopic ones [20,6,18].The equation (1) we consider here is especially useful to capture the hot-electron effects in nanoscale devices [19].Since the delta function is involved in the collision operator, it is more realistic than many commonly used kinetic models that only deal with smoothed kernels [15,8,17,9].Furthermore, it includes the electron-electron interaction which is usually neglected under a low-density assumption [4,5] (not true in our case).Mathematically, the two-scale collisions with delta kernels bring many interesting features: the leadingorder elastic collision does not have a unique null space; only when the electronelectron collision at next level takes into effect, the solution can be driven to a fixed equilibrium -a (zero-momentum) Fermi-Dirac distribution.As a result, the macroscopic asymptotic limit when the mean free path goes to zero is a system of conservation laws for the electron mass and energy, the so-called energy-transport (ET) model [19].
Our goal in this paper is to design an efficient numerical scheme for the Boltzmann transport equation (1).The emphasis will be put in the situation that the parameter α takes a wide range of values: from α ∼ 1 (kinetic regime) to α 1 (diffusive regime).While individual solvers for both regimes are available (or possible), it is desirable to have a unified scheme working for different α, as in practice α may not be uniformly small or large in the entire domain of interest.To make it precise, we want a numerical scheme that is consistent to the kinetic equation ( 1), and when α approaches zero it automatically becomes a macroscopic solver for the limiting ET system, i.e., it is asymptotic-preserving (AP) [14].A natural thought is to use implicit rather than explicit schemes on stiff terms, but this is impractical since the collision operators are too complicated to invert (even forward computation is not an easy task).Inspired by [10], one can seek simple BGK-like operators to penalize these integral operators.However, as the collision possesses three different scales, a closer examination of the asymptotic behavior of the numerical solution reveals that the penalization should be performed wisely, otherwise it won't capture the correct asymptotic limit.The same problem was considered in our recent work [13], where we provided a solution by introducing a threshold on the penalization of the elastic collision.Here we propose an alternative based on a splitting approach.By carefully separating the different scales in equation (1), we are able to recover the correct ET limit at the discrete level.Moreover, compared to the previous idea, the new method is free of choosing threshold and hence is easier to implement.Finally, we comment that most of the AP schemes were developed for problems with two scales (the kinetic and the hydrodynamic ones, for example [14]).Here and in [13], the problem contains an additional, intermediate scale, thus offers new features and challenges for AP schemes.
The rest of the paper is organized as follows.In the next section we briefly review the properties of the collision operators and the energy-transport limit of equation (1).Section 3 describes our AP scheme: we first consider the spatially homogeneous case with an emphasis on the two-scale collision terms, and then include the spatial dependence to treat the full problem.In either case, the asymptotic property of the numerical solution is carefully analyzed.We present several numerical examples in Section 4 to illustrate the efficiency, accuracy, and AP properties of the new scheme.The paper is concluded in Section 5.
2. The energy-transport limit of the semiconductor Boltzmann equation.In this section, we give a brief review of the energy-transport limit of the semiconductor Boltzmann equation (1).The (formal) derivation is a combined procedure of the Hilbert expansion and moment method which mainly relies on the properties of the collision operators Q el and Q ee .We only list here the basic elements that are necessary to understand the numerical schemes.The complete treatment can be found in [7,13].Proposition 1.
1.For any "regular" test function g(ε(k)), 1. Conservation of mass and energy: is the zero momentum Fermi-Dirac distribution function.The macroscopic variables z and T are the fugacity and temperature.
Remark 1.The collision operator Q ee also conserves momentum: R d Q ee (f )k dk = 0.Only when fixed to ε-dependent functions, its null space is (5).
To derive the asymptotic limit of (1), one inserts the Hilbert expansion f = f 0 + αf 1 + α 2 f 2 + . . .and collects equal powers of α: Using the properties of Q el and Q ee , one can show that An important observation here: f 0 = M is not completely determined by ( 6) since the null space of Q el is not unique.It is the joint action of ( 6) and ( 7) that yields the above result.Now plugging f into (1), and taking the moments • (1, ε) T dk on both sides, to the leading order we have Theorem 2.1.In equation (1), when α → 0, the solution f formally tends to a Fermi-Dirac distribution (5), with the position and time dependent fugacity z(x, t) and temperature T (x, t) satisfying the so-called energy-transport (ET) model: where the density ρ and energy E are defined as the fluxes j ρ and j E are given by with the diffusion matrices and the energy relaxation operator W inel ph is Remark 2. Unlike the classical statistics, given macroscopic quantities ρ and E, finding the corresponding Fermi-Dirac distribution ( 5) is not trivial.In fact, ρ and E are related to z and T via ( [11]) where and Γ(ν) is the Gamma function.
3. Asymptotic-preserving schemes for the semiconductor Boltzmann equation based on operator splitting.Equation ( 1) contains three different scales, wherein explicit schemes become extremely expensive when α is small (convection part requires CFL condition ∆t = O(α∆x); collision part requires at least ∆t = O(α 2 )).To design a scheme whose stability is immune from the restriction induced by stiff terms, one usually expects a fully implicit scheme.However, as mentioned in the Introduction, this may result in a large algebraic system that is hard to invert.Our goal is to design an efficient numerical scheme that is uniformly stable regardless of the magnitude of α: only requires parabolic CFL condition ∆t = O(∆x 2 ), and the implicit terms can be treated in an explicit manner.While the stiff convection term has been successfully handled via an even-odd decomposition [15,16], the two-scale collisions is treated recently in [13] using a thresholded-penalization.There, a spatially dependent threshold is placed on the stiffer collision operator such that the evolution of the numerical solution resembles the Hilbert expansion at the continuous level.However, the choice of threshold depends on the accuracy of the elastic collision solver, so its optimal value is not easy to determine.Here we propose a different scheme based on a splitting approach (free of threshold).We will start with the spatially homogenous case whose splitting is somewhat straightforward.We then consider the inhomogeneous case, where the splitting needs to be done carefully in order to capture the correct ET limit.This is indeed inspired by [13].To facilitate the presentation, we always make the following assumptions without further notice: 1.The inelastic collision operator Q inel ph in ( 1) is assumed to be zero, since it is the weakest effect and its appearance won't bring extra difficulties to temporal and spatial discretizations.

The scattering matrix Φ
Then the elastic collision (3) can be written as where In particular, for any odd function f (k), 3. The collision operators (3) and ( 4) can be written symbolically as where the forms of Q + el , Q ± ee should be clear from the definition.3.1.The spatially homogeneous case.In the spatially homogeneous case, equation (1) reduces to where f only depends on k and t.Following [13], we use a BGK-like penalization [10] to remove the stiffness on collision terms: where M is simply chosen as the Fermi-Dirac distribution (5) since it belongs to the intersection of N (Q el ) and N (Q ee ).The essence of penalization is to choose as small as possible so that it is non-stiff or less stiff.Viewing (12), we may choose This choice is sufficient to guarantee AP property and stability as illustrated by our later analysis and numerical results.Now to handle the two-scale collisions, a direct separation of scales in (13) suggests the following first-order splitting scheme: Note in the spatially homogeneous case ρ and E are conserved, so M = M (ε) is an absolute Maxwellian and can be obtained from the initial condition (first find ρ = R d f 0 dk, E = R d f 0 ε dk, and then invert (10) to get z and T so as to define M ).
3.1.1.Asymptotic properties of the numerical solution for the BGK model.In this subsection, we study the asymptotic behavior of the numerical solution to the scheme (15)(16).In the BGK model, where M ee is a general Fermi-Dirac distribution defined by1 , with z, u, T determined through the moments of f : Note that ρ, e are related to z, T via the same system (10) but with the energy E replaced by internal energy e (in fact, e and E have the relation e = E − 1 2 u 2 , thus in the zero momentum case, E = e and we have E rather than e in (10)).
Remark 3. The operator ( 17) is a simple, yet practical approximation to the complicated integral operator.The AP analysis for the full Boltzmann operator (4) is still lacking, but the analysis here sheds some light on the behavior of the solution.We still use (4) in numerical simulations.Now (15)(16) can be rewritten as where Notice that Q + el (f n ) only depends on ε, so when α is small, (18) drives f * toward a radially symmetric function as long as |γ 1 | < 1, i.e. β el > 1 2 λ el ; M * ee computed from f * then has a decreasing momentum and thus approaches M (M and M * ee share the same density and energy); and (19) . Therefore the loop containing these two stages resembles a Hilbert expansion within one time step.To be more precise, we need the following Lemma (its proof is given in the Appendix) characterizing the distance between M * ee and M .Lemma 3.1.Consider two Fermi-Dirac distribution functions M ee,1 and M ee,2 with the form and assume T 1,2 (x, t), z 1,2 (x, t) > 0 and Here Recall that the scheme (15-16) conserves mass and energy: If we take the moments •k dk on both sides of (18)(19), we have where the moments of Q + el and M vanish due to their radial symmetry.Therefore, thanks to Lemma 3.1.Notice that a combination of (18)(19) implies Applying the linear operator Q el on both sides of (22) yields where the first inequality is based on the facts and the second inequality is obtained using (21) since As a result of ( 21) and ( 23), for any m > 0, there exists an integer N such that for all n > N we have On the other hand, (18-19) can be written as and thus ), for n > N , because of (24).Since |γ 2 | < 1, therefore, no matter what the initial condition is, f n+1 will eventually be driven to the desired Fermi-Dirac distribution M , and the convergence rate depends on the magnitude of γ 2 .Remark 4. We would like to point out that in this splitting framework, the two stages in each time step alternate, that is, the solution is driven closer to its radially symmetric counterpart (with momentum decreased) in the first stage and then slightly to the Fermi-Dirac distribution (with momentum preserved) in the second stage, which makes the relaxation toward the final (zero-momentum) Fermi-Dirac distribution 'oscillatory'.This is different from the thresholded penalization scheme in [13] where the solution is driven all the way to a radially symmetric function at the beginning, and then toward the final Fermi-Dirac distribution after the onset of threshold.

3.2.
The spatially inhomogeneous case.We now include the spatial dependence to treat the full equation (1).Following [13], we reformulate it into a set of parity equations [15].Denote f + = f (x, k, t), f − = f (x, −k, t), then they solve we have where we have used the fact that j is an odd function in k, and thus Q el (j) = −λ el j.
For (25-26), the same penalization as in the homogenous case suggests Note here M = M (x, ε, t) is the local Fermi-Dirac distribution.The coefficients β el and β ee are chosen the same as in ( 14) except that β ee can also be made spatially dependent.
Rewrite the above equations into a diffusive relaxation system [16], we have where 0 ≤ θ(α) ≤ 1 α 2 is a control parameter chosen as θ(α) = min 1, 1 α 2 .Unlike the simple separation of O( 1α ) and O( 1 α 2 ) terms in the spatially homogenous case, we propose the following first-order splitting scheme for (27-28): • stage 1 where M * = M n since Q el conserves mass and energy.• stage 2 where M n+1 is computed via (10).Here ρ n+1 and E n+1 are first computed by taking the moments • (1, ε) T dk of (31), wherein the right hand side vanishes owing to the conservation of mass and energy.
3.2.1.Asymptotic properties of the numerical solution for the BGK model.Let's again assume that Q ee (f ) = M ee (f ) − f .A reformulation of (29-30) leads to Rearrange (31-32), we have Here M * ,± ee = M ee (f * ,± ) = M ee (r * ± αj * ), γ 1 , γ 2 are defined the same as in (20).First we point out that j n+1 (n ≥ 0) has magnitude O(1) since from (34), (36) Notice that j * is an odd function in k, thus M ee (r * ) and M ee (r * ± αj * ) share the same mass and energy.Since r * has zero momentum, we have M ee (r * ) = M * , and therefore Lemma 3.1 implies (37) Henceforth, when α 1 (so θ = 1), we have from (36) if all functions are smooth.Similar to the homogenous case, apply operator Thus there exists N such that for all n > N , Q el (r n ) = O(α).On the other hand, a reformulation of (33) and (35) yields Therefore, we have From the above discussion, we know that when α → 0, no matter what the initial condition is, immediately and after an initial transient time These are the desired AP property we want (compare with ( 8)): plugging them into (29), (31) and taking the moments • (1, ε) T dk, we get exactly a first-order time discretization for the ET system (9).

Numerical examples.
In this section, we present several numerical examples using our AP scheme (29-32).Here the space discretization follows the upwind scheme with flux limiter and wave vector discretization uses the spectral method.More details are refereed to [13].
In what follows, we always take k ∈ [−L k , L k ] 2 with L k = 9.2, and x ∈ [0, L x ] with L x = 1.N k is the number of points in each k direction, N x is the number of points in x direction.We assume the periodic boundary condition in x and choose L k large enough so f ≈ 0 at |k| = L k .The time step ∆t is chosen to only satisfy the parabolic CFL condition: ∆t = O(∆x 2 ) (independent of α).

AP property. Consider equation (1) with non-equilibrium initial data
The electric field ∂ x V is set to be one.
We check the asymptotic property by looking at the distance between r and M at each time step, i.e., errorAP n L ∞ = max x,k1,k2 The results are gathered in Figure 1, where we observe that, unlike the thresholded penalization which experiences a clear two-stage convergence, the solution by the splitting approach undergoes a faster convergence during the first few steps, and smoothly transits to the final equilibrium.This is because the two collisions, even though they are not in the same scaling, act together all along.
4.2.1-D n + -n-n + ballistic silicon diode.We next simulate a 1-D n + -n-n + ballistic silicon diode, which is a simple model of the channel of a MOS transistor.The initial data is taken to be For Poisson equation ( 2), we choose h(x) = ρ 0 (x) = f 0 dk, D 0 = 1/1000, with boundary condition V (0) = 0, V (L x ) = 1.The doping profile h(x) is shown in Figure 2. We consider two regimes: one is the kinetic regime with α = 1, where we compare our solution with the one obtained by the explicit scheme (forward Euler); the other is the diffusive regime with α = 1e − 3, where our solution is compared with that of the ET system using the kinetic solver [13].Good agreements are obtained in Figures 3, 4.Here the macroscopic quantities plotted are density ρ, energy E, temperature T , fugacity z, electric field E = −∂ x V , and mean velocity ū defined as ū = j ρ /ρ.This example is borrowed directly from [13] and the results are very similar to those computed using the thresholded penalization.
with α 0 = 1e − 3, thus it increases smoothly from α 0 to 1, and then suddenly drops back to α 0 .The plot of α is shown in Figure 5, which involves both kinetic and diffusive regimes.The initial condition is taken as in (38).We compare the solution using ∆x = 1/40 and ∆t = 0.2∆x 2 with a more refined solution with ∆x = 1/160 and ∆t = 0.05∆x 2 in Figures 6, 7 with η = 0.01 and 1, respectively.

Conclusion.
We designed an asymptotic preserving scheme for a Boltzmann-Poisson system that characterize the transport of charge carriers in the semiconductor.In a diffusive regime where the collisions are not in the same scales, the system approaches an energy-transport system.Besides the stiff convection terms, the two-scale collision operators pose new difficulties since the simple BGK penalization fails to capture the correct limit.Inspired by the two stages in the convergence to the equilibrium, we propose a splitting setting such that the relaxation of the numerical solution to the local equilibrium resembles the Hilbert expansion at the continuous level at each time step.Therefore, the numerical solution experiences an alternating path towards the equilibrium with one direction heading to the radially symmetric function and the other to the Fermi-Dirac distribution.The main advantage compared to the thresholded penalization in [13] is that it is free from the choice of threshold.We analyzed this asymptotic behavior using a simplified BGK  model.Several numerical results confirmed the asymptotic-preserving property for any non-equilibrium initial data, as well as the uniform stability of the scheme with respect to the mean free path, from kinetic regime to diffusive regime.