Implicit-Explicit Runge-Kutta schemes for the Boltzmann-Poisson system for semiconductors

Implicit-Explicit Runge-Kutta schemes for the Boltzmann-Poisson system for semiconductors

Giacomo Dimarco111 Institut de Mathématiques de Toulouse, Toulouse, France., Lorenzo Pareschi222Department of Mathematics and Computer Science, University of Ferrara, via Machiavelli 35, 44121 Ferrara, Italy., Vittorio Rispoli333Department of Mathematics and Computer Science, University of Ferrara, via Machiavelli 35, 44121 Ferrara, Italy.

In this paper we develop a class of Implicit-Explicit Runge-Kutta schemes for solving the multi-scale semiconductor Boltzmann equation. The relevant scale which characterizes this kind of problems is the diffusive scaling. This means that, in the limit of zero mean free path, the system is governed by a drift-diffusion equation. Our aim is to develop a method which accurately works for the different regimes encountered in general semiconductor simulations: the kinetic, the intermediate and the diffusive one. Moreover, we want to overcome the restrictive time step conditions of standard time integration techniques when applied to the solution of this kind of phenomena without any deterioration in the accuracy. As a result, we obtain high order time and space discretization schemes which do not suffer from the usual parabolic stiffness in the diffusive limit. We show different numerical results which permit to appreciate the performances of the proposed schemes.

Keywords:IMEX-RK methods, asymptotic preserving methods, semiconductor Boltzmann equation, drift-diffusion limit.

1 Introduction

The application of kinetic theory to the modeling of semiconductor devices simulations is a long dated but still very active research field because of the richness and diversity of observed phenomena [33, 35]. From the mathematical point of view, semiconductor devices can be accurately described by kinetic equations when the mean free path of the particles is large compared to a macroscopic characteristic length of the system ([2, 8, 9]). On the other hand, when the average time between particles’ collisions is small, the relevant scaling is the diffusive one and in such regimes these systems can be described by macroscopic drift-diffusion equations [39, 35, 40, 20]. Unfortunately, this passage from the microscopic to the macroscopic description leads to challenging numerical difficulties. Indeed, when the macroscopic reference length is several orders of magnitude larger than the mean free path, the kinetic equation contains stiff terms. This means that classical numerical methods need time step restrictions which make their use prohibitively expensive. In these cases, it becomes attractive to use domain decomposition strategies, which are able to solve the microscopic and the macroscopic models wherever it is necessary. These approaches have been largely studied for kinetic equations both for the diffusive [4, 12, 14, 29] and for the hydrodynamic scaling [13, 15, 16]. However, even if these methods are very efficient they are affected by some difficulties due to the fact that it is not always a simple task to define the different regions of the domain in which the use of a macroscopic model is fully justified. Thus in the recent past alternative strategies have been studied, the so-called Asymptotic Preserving (AP) schemes. They consist in solving the original kinetic model in the full domain avoiding the time step restriction caused by the presence of different stiff terms in the equations. The AP methods automatically transform the original problem in the numerical approximations of the relevant macroscopic model when the scaling parameter goes to zero [25, 21, 26, 27, 28, 30, 5, 32, 7, 17, 22, 23, 36, 37]. Recently, this approach has been considered in the framework of Implicit-Explicit (IMEX) Runge-Kutta schemes with the aim of deriving high order numerical methods which are accurate in all regimes [3, 11, 19, 18, 6, 38]. This means that such schemes are able to preserve the desired order of accuracy even in the limit when the scaling parameter tends to zero.

In this paper we develop high order schemes for the resolution of the Boltzmann-Poisson semiconductor model. In this kind of problems, since the characteristic speed of the hyperbolic part of the kinetic equation is of the order of (where is proportional to the mean free path), the CFL condition for an explicit approach would require . Of course, in the diffusive regime where , this would be too much restrictive since a parabolic condition would suffice to solve the limiting equation. This kind of difficulties has been studied in [26, 27, 28, 30, 36, 37] using different semi-implicit approaches. However, most of the previous literature on the subject, originates consistent low order explicit schemes for the limit model. Such explicit schemes clearly suffer from the usual stability restriction . Here, on the contrary to previous approaches, we are able to guarantee for all regimes a linear time step limitation , independently from , and high order accuracy in time and space. In order to accomplish this task, we first rewrite our system using the parity formalism then, following [6], we add and subtract the limiting diffusive flux to the convective kinetic flux. We then discretize the reformulated problem by Implicit-Explicit Runge Kutta method. As we will show this will allow to get Asymptotic Preserving high order schemes which uniformly work for the different values of and which automatically originate an IMEX Runge-Kutta method for the limiting convection-diffusion equation in which the diffusive term is discretized implicitly. Finally, in order to permit realistic simulations we consider the challenging case of complex interaction operators. These terms are handled by introducing a penalization technique which permit to avoid the inversion of such operators when their stiff character suggests an implicit treatment.

The rest of the paper is organized as follows. First we present the kinetic semiconductor equation, its drift-diffusion limit and the reformulated parity system. The third section provides the basis of the numerical method describing in details the time discretization. The fourth section describes phase-space variables discretization, based on conservative finite difference schemes for the space variables and on a Gauss-Hermite approximation for the velocity variables. Numerical results for the proposed schemes are presented in the fifth section. A concluding section ends the paper.

2 The Boltzmann equation and its drift-diffusion limit

We consider the Boltzmann equation under the diffusive scaling which describes the time evolution of electrons inside semiconductor devices. Let be the density distribution function for particles at time , where position and velocity variables and are such that , with and or . Under these assumptions, the time evolution of the system is described by [29]


In this formula, is an integrable function of which models the generation and recombination process inside the semiconductor, is proportional to the mean free path and is the electric field which is self-consistently computed solving a suitable Poisson equation for the electric potential . Constants and are respectively the elementary charge and the effective mass of the electrons. The anisotropic collision term is defined by


where is the constant in time normalized Maxwellian at temperature

and is the anisotropic scattering kernel which is rotationally invariant and satisfies

for some given constant . We also assume that the collision frequency satisfies the following bound for some positive constant [30]


The collision operator is bounded and nonnegative on the Hilbert space, see [39, 30],
and it has a one-dimensional kernel spanned by . Moreover, if we let be the unique solutions to the problem

then, since is rotationally invariant, it follows that there exists a positive constant such that


Thus, defining the total mass as

one can show (using for example the Hilbert expansion with , see [32] for instance) that when , is approximated by with satisfying the drift-diffusion equation [29, 35]


In this equation, is the diffusion coefficient defined implicitly in terms of the cross section by (2.4), the constant is the so-called mobility given by the Einstein relation and is the integral of the generation recombination function

2.1 Even and odd parities

We now define the even and odd parities formalism which we will use in the development of the numerical schemes. To this aim we split equation (2.1) into two equations, one for and one for 


Next, we introduce the so called even parity and odd parity defined by


Now, adding and subtracting the two equations in (2.6) leads to


where is the collision frequency defined in (2.3) and where we used the property that

since is an odd function. We also assume is an even function, which is a physically consistent choice in semiconductor simulations. An important advantage of this formulation is that now only one time scale appears in our new system (2.9).

From the above formulation it is easy to see that in the limit , formally we get that the integral of the even parity with respect to the velocity variable is the solution to (2.5). In fact, the macroscopic variable can be expressed in terms of as


Taking the formal limit in system (2.9) we get


The first equation implies that being Ker Span. Now, replacing equation (2.12) and in the first equation of system (2.9) and integrating over the velocity space gives the drift diffusion equation (2.5). From such computation we recover, moreover, that is given by the following formula

In the following, in order to simplify notations, we will assume and , which implies .

3 Time discretization of the multi-scale system

We recall here the system (2.9) that we want to solve


and summarize the properties we demand to our time discretization scheme:

  1. The scheme has to be asymptotic preserving (AP). This ensures stability condition independently from . In other words, the scheme should satisfy the following definition

    Definition 3.1

    A consistent time discretization method for (3.1-3.2) of stepsize is asymptotic preserving (AP) if, independently of the initial data and of the stepsize , in the limit becomes a consistent time discretization for equation (2.5).

  2. The scheme has to be an high order asymptotically accurate (AA) method. This means we want our scheme to satisfy the following definition

    Definition 3.2

    A consistent time discretization method for (3.1-3.2)) of stepsize is asymptotically accurate (AA) if, independently of the initial data and of the stepsize , the order of accuracy in time is preserved for .

  3. The scheme should solve, in the limit , the drift-diffusion equation (2.5) with an implicit treatment of the diffusion term. This ensures a stability condition for the time step which is of the order of the space discretization: .

  4. We want to avoid the difficult inversion of complex collision operators which occurs when classical implicit solvers are used.

3.1 The Boscarino-Pareschi-Russo reformulation

We now focus on the first three points of the above list while we will discuss the last point at the end of the section. In order to construct time discretizations which satisfy the above requirements we reformulate system (3.1-3.2) following the approach introduced in [6]. We add to both sides of equation (3.1) a term that permits to choose which kind of time integrator (if explicit or implicit) we intend to use for the diffusive term in the drift-diffusion equation (2.5). This term reads

where is a positive function such that . The modified system reads


As we will see, the introduction of this new term allows to avoid the parabolic time step limitations for the limit drift diffusion equation (2.5). We will discuss later the different possible choices for .

In order to write the general time discretization formulation, we rewrite the previous system (3.3) in a more compact form as


where we defined

Remark: We point out that in this paper we consider the case in which does not depend on the space variable. More precisely, could depend on but not on and . Since in the general case instead depends on and , in the work we are assuming that in the space-time domain the range in which may vary is small enough to allow us to keep constant. We numerically observed that this condition is, indeed, true for a large set of values of . We remind to a future work the development of schemes which could take into account the possibility for to vary.

3.2 IMEX Runge-Kutta schemes

An IMEX Runge-Kutta scheme [3, 6] applied to the above system reads for the internal stages as


while the numerical solution is given by


In the above formulas, matrices , for and are matrices such that the resulting scheme is explicit in and implicit in , and . In general, an IMEX Runge-Kutta scheme, is characterized by the above defined two matrices and by the coefficient vectors , . Since computational efficiency is of paramount importance, in the sequel we restrict our analysis to diagonally implicit Runge-Kutta (DIRK) schemes for the source terms ( for ). The use of a DIRK scheme is enough to ensure that the two transport terms in the two equations of system (3.3) are always explicitly evaluated. In fact, observe that, when the integration for the odd parity is performed the values of are already available for the corresponding stage. This is due to the use of a partitioned Runge-Kutta approach for the time integration of our system. The type of schemes introduced can be represented with a compact notation by a double Butcher tableau,

where coefficients and are given by the usual relation and . IMEX schemes are a particular case of additive Runge-Kutta methods and so the order conditions can be derived as a generalization of the notion of Butcher tree; we refer to [24] for more details on the order conditions. Before stating the main properties concerning asymptotic preservation and asymptotic accuracy, we characterize the different IMEX schemes accordingly to the structure of the DIRK method. Following [6], we call an IMEX-RK method of type A if the matrix is invertible, or equivalently , while we call it of type CK (see [11]) if the matrix can be written as


with and the submatrix invertible, or equivalently , . We write also the matrix for the explicit Runge-Kutta method


where and . We now introduce two useful definitions to characterize the properties of the methods in the sequel. An IMEX-RK scheme is called implicitly stiffly accurate (ISA) if the corresponding DIRK method is stiffly accurate, namely . If in addition the explicit matrix satisfies , the IMEX scheme is said to be globally stiffly accurate (GSA) or simply stiffly accurate.

Note that for GSA schemes the numerical solution is the same as the last stage value, namely and . We recall that the order conditions for GSA type IMEX schemes are particularly restrictive since and . Another restrictive condition is, for high order methods, the request that the matrix is invertible for details see [19, 6]. All these requests make very difficult to derive IMEX GSA schemes of order higher than two.

The detailed analysis of the AP properties of the proposed IMEX schemes is reported in the Appendix. In that part we will give sufficient conditions which guarantee that the schemes are AP and AA. Here, we only recall the main results. Type IMEX schemes are Asymptotic Preserving and Asymptotically Accurate. If in addition they are GSA the distribution function is projected over the equilibrium at each time step. Two sufficient conditions for type IMEX schemes which guarantee the AP and AA properties is that they are and that the initial data are close to the equilibrium state (we say in this case that the initial data are consistent with the limit problem). Again in this case, we get also sufficient conditions to assure that the distribution function is projected over the equilibrium state at each time step.

Here we show how the density values are obtained and that in the limit the schemes solve the drift-diffusion equation with an implicit treatment of the diffusion term. Observe in fact that in order to use schemes (3.5)-(3.8), we need to know the values of the density distribution and of its stages implicitly. These values are obtained by integrating with respect to velocity equation (3.5) and (3.7). We denote by , and for , the column vectors of the stages for and respectively and by the vector of the stages of the mass density ; it holds that . Moreover, we denote by the vector containing stage values for and similarly for . The equation for the internal stages in vector form reads

which we can also rewrite as


while the equation for the numerical solution is

In the above equation the integral of is implicit. However, it can be (explicitly) solved by inverting the matrix describing the discretized diffusion operator related to the term . This permits to know the density implicitly by the explicit knowledge of and . Finally, in the limit regime, as shown in the Appendix, the term is reduced to which means that we get as desired an implicit discretization of the diffusion term.

3.3 A linearization technique for the implicit collision term

In the numerical method described in the previous paragraph the collision operator has to be implicitly computed. Then, it is necessary that one is able to invert it. This is usually not the case since, in general, collisions are represented by nonlinear multidimensional operators which could be costly to compute or even more to invert. For this reason, we choose to penalize with a suitable operator which needs to be easier to invert and would not change the asymptotic behavior of the solution. The second condition is mathematically expressed by Ker Ker Span. This strategy has been proposed for the Boltzmann equation in [21] and subsequently studied in the context of IMEX schemes in [19].

In order to write the modified IMEX schemes, we add and subtract to the collision term an operator and then we combine the implicit and the explicit solvers as follows:

A possible choice for is a first order approximation of the original operator, obtained using an expansion of near the equilibrium distribution : .
Since it is not always possible or easy to compute analytically , one can choose to approximate it. A possibility is then


where is an upper bound of .

Regardless from the choice of , we apply the IMEX strategy to the penalized system in the following way


Observe that computing operator with the implicit solver stabilizes also the non-linear collision operator, without changing the asymptotic behavior of the solution. However, this stabilization is not straightforward, on the contrary in order to stabilize the reformulated system it is necessary that the coefficients of the scheme used for the time integration of the linearized collision operator dominate those used for the time integration of the original operator. Such technique allows us to treat very general collision operators. We prove in the appendix that both and type IMEX schemes if also Globally Stiffly Accurate are AP and AA. More in details, schemes needs the additional hypothesis of consistent initial data to assure that the asymptotic properties are satisfied. The request that the schemes are in this case becomes necessary for the stability in the limit of zero mean free path.

4 Phase-space discretization

We discuss in this section the discretizations of velocity and space variables. Concerning the velocity variable, because of the particular structure of the problem, it is convenient to approximate and using a Gauss-Hermite expansion. We decompose then the unknowns and as follows [29, 26]: and , with and . In this way it is possible to exploit the Gauss-Hermite approximations to efficiently and accurately compute the derivatives in and the collision operator (which is an integral in ). From (3.3) we have:


From equation (4.1) we obtain

from which


follows, with and . Similarly, from equation (4.2) we have


We conclude this section with an example for : in the particular case in which in (2.2), we get the co called relaxed time approximation (RTA). It is easy to see then that

and thus it holds that .

4.1 Velocity discretization

We describe the Gauss-Hermite approximation ([10, 30, 26]) in the monodimensional case. The multidimensional case is obtained applying the monodimensional rule dimension-by-dimension.

Let consider and , with

being the Hermite expansion. Here are the renormalized Hermite polynomials and coefficients and can be computed thanks to the inverse expansion (we refer to [30] for more details). The computation of the collision operator becomes

where are points and weights of the Gauss-Hermite quadrature rule. Finally, the derivatives with respect to , which are given by



Remark: Coefficients for any component of can be computed at the beginning of the simulation and stored in a matrix since they do not depend on functions and .

4.2 Space discretization

In this section we emphasize some requirements about the space discretization of the system. We want our scheme to work both in the kinetic regime (), in which the hyperbolic behavior is more relevant, and in the limit regime (), in which the system is characterized by diffusive behavior. Moreover, the characteristic speeds of the system (which are of the order of ) tend to infinity as and so shock capturing methods based on characteristics directions, such as, e.g., upwind methods, become useless. On the other hand, central differences schemes avoid excessive dissipation but, when is not small or when the limiting equations contain advection terms, may lead to unstable (or not accurate) discretizations.

In order to overcome these well-known facts and to have the correct asymptotic behavior, we fix some general requirements for the space discretization:


correct diffusion limit: as we already observed in previous section, if we want a correct approximation in the limit case , we need that and to use the same space discretization for the transport terms in (4.3) and (4.4);


compact stencil: we want to use a scheme with a compact stencil in the diffusion limit . This property is satisfied if point i) is satisfied and we use a suitable discretization for the second order derivative that characterizes the diffusion limit;


shock capturing: the chosen scheme should be based on high order shock capturing fluxes for the convection part. This is necessary not only for large values of but also when we consider convection-diffusion type limit equations with small diffusion. The high order fluxes are then necessary for all space derivatives except for the second order term on the right hand side in (4.3);


avoid solving nonlinear algebraic equations: in order to have a more efficient method, we do not want to solve the nonlinear equation which comes from the implicit treatment of the space derivative in equation (4.4) for , (and in (3.14) for the odd parity ). To achieve this we have chosen a partitioned approach for the time integration, thanks to which the values of are already available from the previous solution of (4.3).

4.2.1 Modified fluxes

In order to satisfy the above requirements we choose to use a Lax-Friedrichs type flux with high order WENO reconstruction [41, 9, 8]. This gives us the ability to ensure accuracy and also to stabilize the solution in the presence of discontinuities or arising shocks.

Our strategy is the following: in the kinetic regime, where transport dominates the dynamics, we use the standard Lax-Friedrichs type flux with WENO reconstruction for the derivatives [6, 27, 36]. In this case, indeed, this is a proper strategy which allows us to obtain also high order accuracy. When considering the limiting regime instead, we cannot use the Lax-Friedrichs scheme as it is. As we will show later in this section, the numerical viscosity introduced by such scheme in this case is proportional to . Clearly, when the mean free path goes towards zero () such quantity is too large and causes loss of accuracy (see for instance [36, 37]). Thus, in such situation we decide to bound the numerical viscosity modifying the fluxes. This is possible because when becomes small the diffusive regime becomes dominant and thus stability is granted by the “physical” viscosity given by the system itself. We point out that, in this work, the stability of the proposed modified fluxes approach is supported by numerical evidence. Theoretical estimates for this technique and for this kind of scaling will be the subject of a future work.

We present here the modified fluxes approach using a model problem. Given the transport equation for the unknown

with a hyperbolic flux, we consider a complete discretization in conservation form


where and the parameter represents the numerical viscosity. The standard Lax-Friedrichs scheme requires the value . Such scheme is stable if the following two inequalities are satisfied


Formula (4.6) is the basis for the construction of a conservative numerical flux which can be implemented using ENO or WENO high order reconstructions [41].

To derive the modified fluxes, let consider now the prototype system


which shares the same structure of our original problem. In the limit , if , the above system leads to the drift-diffusion equation . We then write a semi-discrete approximation of (4.8) as

with numerical fluxes given by

where is the numerical viscosity as before. Rewriting the numerical fluxes as


we have


The stability conditions (4.7) in this case read:

As pointed out before, we want to avoid such restrictive time step. To this aim, close to the limit we modify the numerical viscosity setting


More in details, until we reach a regime in which the physical diffusion is not large enough to guarantee stability we need to satisfy (4.10). On the other hand, when physical diffusion becomes relevant, we can avoid such restrictions and we can choose the modified fluxes (4.11). The practical choice we did in our numerical tests is


To summarize, we write here the complete numerical discretization of system (3.13)-(3.14). We denote the values for and varying in the phase-space index set and at time . We denote in the same way the other variables appearing in the sequel. For stage vectors and we denote and for . Then we have


and and for the numerical solution. To make formulas more readable we defined some shorthands. Operator stands for the numerical discretization of the transport derivatives: , obtained as in (4.9) and it reads



The term is while and stand for the derivative with respect to and are obtained by means of (4.5).

The diffusion term in the r.h.s. of equation (4.13), i.e. the derivative , stands for the standard central second order finite difference technique, i.e.

when the second order time discretizations is used, and for the standard central fourth order finite difference technique when the third order time discretizations are used. Regarding the numerical viscosity in equations (4.13) and (4.14), if we are in the kinetic regime we choose the physical values given by (4.10) while in the limiting regime we consider the modified ones given by (4.11). The space discretization of the electric potential can be performed by standard methods [30, 40].

4.3 Boundary conditions

The treatment of boundary conditions for the Boltzmann-Poisson problem in the diffusive limit is, in the general case, a very hard task. It is necessary to tackle several difficulties, such as complex geometries of the boundaries and to take into account boundary layers. For a consistent treatment of boundary condition see, for instance, [30, 31, 32] and references therein. All these aspects are out of the scope of this article and we will only deal with assigned, constant in time boundary data, the so called maxwellian injection. We show now a possible strategy for a consistent numerical implementation of such conditions in the one-dimensional situation.
For a maxwellian injection is defined by

for , where and are assigned nonnegative functions proportional to the maxwellian distribution . We numerically approximate these conditions in two different ways, depending on the regime in which the system is.

In the kinetic regime, for we set , we extrapolate (the outgoing particles) from the values of and inside the domain and then we define and thanks to the parity formulas (2.7) and (2.8). At the right boundary a similar treatment is used.

In the diffusive regime instead, to get a boundary condition for and we use the relations (for positive only)


Then we consider equation (2.12), which gives a good approximation of when is small, i.e. , and applying it in (4.16) one gets

To approximate and we use one-sided finite difference discretizations of the desired order and to approximate we observe that from (4.16) it holds that

which in the end leads to (up to )


5 Numerical tests

In this section, we present several numerical results to test the performance of the proposed schemes. We show that our schemes are computationally very efficient, i.e. while at the same time they preserve high order of accuracy in all regimes. The setting is a monodimensional phase-space, i.e. .

In our computations we use two different scattering cross-sections: a simple isotropic case with a constant cross-section , this corresponds to the relaxation time approximation (RTA), for which the collision operator has the simple form

and a regularized anisotropic cross section [34] for electron-phonon interactions (EPI)

where is a smoothed delta function with a positive constant (we set ). When we consider the EPI model, we apply the penalization technique using operator given by (3.12) with (observe that this corresponds to nothing else but the RTA approximation).

Concerning the value of , we choose a simple form given by


As we already observed, in this work we assume has a constant for a given value of . We recall here that more accurate choices are possible for , i.e. , and we remind to a future work for a deeper analysis of this aspect.

In all our numerical tests the discretization in the velocity space is obtained using Gauss-Hermite quadrature points, with : there are 8 nodes for positive velocities and 8 for negative ones (we used scaled values in order to consider the range , with ). The influence of the number of quadrature points on the accuracy of the results is treated, for instance, in [30]. The IMEX schemes we used for our simulations are the second order IMEX ARS-(2,2,2) scheme [3] and the third order IMEX BPR-(3,5,3) scheme [6]. For the sake of completeness we report the Butcher tables of the schemes in the Appendix. We compare our results also with a simple first order IMEX scheme, obtained by combining the first order implicit and explicit Euler schemes. A reference solution is always reported for all the tests.

Test # 1

In this problem we have a potential well in the left half of the slab. We test the behavior of the scheme when the system is subject to a constant in time, electric field which varies along the -axis: it ranges from a minimum value of -10 to a maximum value of . We perform simulations in both kinetic and fluid regimes, using both scattering kernels.
In the kinetic regime, i.e. , we stop our simulations at time , with using grid points and with an initial distribution given by . At the boundaries we set the values and and we approximate them as described in previous section. The other parameters of the simulations are

with . Since we are in the kinetic regime, the time step is given by the hyperbolic condition


For this test we set the CFL constant to for all schemes.

Figure 1: Test # 1: comparison of the reference (line) and numerical ( () I order, () II order and () III order) mass distributions in the kinetic regime with potential well for the RTA model, . On the -axis the space variable, on the -axis the mass .

The computed solutions of the mass density are presented in figure 1. We show results for the RTA kernel since the obtained results for the EPI model are similar. We compare the reference density, obtained with a fourth order explicit RK scheme with third order WENO reconstruction using , with the first, second and third order IMEX approximations. As expected, the third order scheme gives a more accurate solution.

Next in figure 2 we report the results obtained in the fluid regime in the RTA case. We set in this case and stop the simulation at using grid points while the other parameters are the same as in the kinetic test case. Now the time step is