pp. X–XX A NUMERICAL MODEL OF THE BOLTZMANN EQUATION RELATED TO THE DISCONTINUOUS GALERKIN METHOD

We propose a new deterministic numerical model, based on the discontinuous Galerkin 
method, for solving the nonlinear Boltzmann equation for rarefied gases. 
A set of partial differential equations is derived and analyzed. 
The new model guarantees the conservation of the mass, momentum and energy for homogeneous 
solutions. 
We avoid any stochastic procedures in the treatment of the collision operator of the 
Boltzmamn equation.

1. Introduction. The classical Boltzmann kinetic equation [9], [10] describes neutral particle transport phenomena. Today numerical solutions of the Boltzmann equation are requested to solve problems in different fields of real-world applications. There are two classes of computational methods, which are used to solve the kinetic equation. In the first of techniques, the well-known Direct Simulation Monte Carlo (DSMC) method, the molecular collisions are considered on a probabilistic rather than a deterministic basis. The literature of the applications of this method is very vast. In the second class, deterministic methods, the Boltzmann equation is discretized using a variety of methods and then solved directly or iteratively. The computational complexity is high due to the large number of independent variables. This heavy computational cost explains why kinetic equations are traditionally simulated by the Direct Simulation Monte Carlo methods. As examples of papers dealing with deterministic schemes for the Boltzmann equation, we cite some references [1], [6], [7] and the book [8], but there are many other interesting papers.
Recently, in the framework of the charge transport in semiconductor devices, deterministic solvers were considered in the literature. In this case the kinetic model, describing of electron flow in semiconductors, is given by the Boltzmann equation coupled to the Poisson equation. The methods applied in [2], [3], [4] provide accurate results which, in general, agree well with those obtained from DSMC simulations, sometimes at a comparable or even less computational time. In the last paper the authors employed the discontinuous Galerkin (DG) method. It is a finite element method using discontinuous piecewise polynomials as basis functions and numerical fluxes based on upwinding for stability. The method is widely used for solving partial differential equations, but it seems to be appropriate for solving also kinetic equations. The method has the advantage of flexibility for arbitrary unstructured meshes, with a compact stencil, and with the ability to easily accommodate arbitrary hp-adaptivity. These properties are important as for no regular domains as well for high dimensional problems. In this case, the total number of the unknowns may be reduced by using large cells where the solution is smooth and small cells only in the zone where the solution exhibits strong variations. For more details about DG scheme for convection dominated problems, we refer to the review paper [5].
According to the discontinuous Galerkin method the starting point for solving the Boltzmann equation is the weak formulation of the equation If we multiply both sides of the equation by a test function φ(x, ξ) and we integrate with respect to the coordinates x and the velocity ξ, then we obtain the equation where X denotes the x-domain. At this point, it is necessary an integration by parts to move the derivative with respect to x from the function f to the test function φ. This requires some information on boundary conditions and test functions. Choosing a suitable finite dimensional linear space V of functions of x and ξ alone, and assuming that f can be approximated by means of a linear combination of elements of this space, with coefficients depending only on the time t, then we can derive from Eq. (2) a closed set of ordinary differential equations, by using as many test functions φ(x, ξ) as the dimension of the space V. This is a very synthetic description of a standard application of the DG method; one needs many further details for concrete implementations.
In this paper, we are interested to show a partial application of the DG method, which, in general, will give a set of partial differential equations having some interesting properties. Now, we assume that the generic test function φ depends only on the velocity ξ. If we multiply both sides of the Boltzmann equation (1) by a test function φ(ξ) and we integrate with respect to the velocity ξ, then we obtain the equation (3) This recalls the classical moment method, but, in our framework, the test functions will have compact support. Now, choosing a suitable finite dimensional linear space of functions of ξ alone, and assuming that f can be approximated by a linear combination of elements of this space, with coefficients depending on the x and t, then we can derive from Eq. (3) a closed set of partial differential equations, using the right number of test functions φ(ξ). Some numerical scheme is required in order to solve this set of partial differential equations. Of course Eq. (3) reduces to an ordinary differential equation for spatial homogeneous solutions. There are some reasons to investigate on Eq. (3) instead of Eq. (2).
• The treatment of Eq. (2) requires the knowledge of the domain X and the boundary conditions on this set. This is not necessary for Eq. (3) since the domain of the velocities is the whole space R 3 ; so, we will give results, which are independent on the data concerning the domain X. • We will obtain explicit expressions of the right hand side of Eq. (3); so this term will not require a further analysis in solving the partial differential equations. We remark that this term is an integral of the collision operator, the most difficult part of the Boltzmann equation. • We will use the discontinuous Galerkin method to derive the set of partial differential equations, but Eqs. (3) can solved by means of different techniques. • The set of equations (3) will be ready for numerical investigations in the case of spatial homogeneous solutions.
The plan of the paper is the following. In Section 2, we introduce the Boltzmann equation and its weak form. Section 3 is devoted to the numerical model. A closed set of partial differential equations is obtained and some properties are shown. The numerical parameters of the partial differential equations are analyzed in Section 4. Section 5 treats conclusions and future work.
2. The Boltzmann equation. The purpose of this section is to briefly introduce the classical nonlinear Boltzmann equation for monatomic gases, to recall wellknown properties and to derive simple results. According to the standard notation, the Boltzmamn equation takes the form The one-particle distribution function f depends on time t, position x and velocity ξ. The kernel W of the collision operator is defined by The function K is related to the interaction law between colliding particles. The Dirac distributions guarantee momentum and energy conservation during the binary collisions. It is immediate to verify that the following symmetry properties hold. The collision operator of the Boltzmann equation (4) is usually written in different way, since Dirac distributions are employed and a five-fold integral is derived.
In this paper we find useful the form of collision term given in (4) or the partially reduced integral operator, where only the Dirac distribution describing momentum conservation is used to reduce the order of integration. The integration with respect to the variables ξ * and ξ * may be performed ab initio in the lost term of the collision operator, since the unknown f is not involved. To this aim, we must consider the total cross section It is related to the collision frequency. It is easy to prove (see Appendix) that the integral (7) is a function of ξ and ξ * only through |V| and it is equal to 2.1. Weak form of the collisional operator. Let φ : R 3 → R be a measurable function. We recall a well-known result. Assuming the existence of the integrals, the right hand side of Eq. (3) can be written as follows Here, as in the following, to simplify the notation, often we omit to write the variables t and x, explicitly. Now, if we define than the integral (9) becomes We note that the function K(φ; ξ, ξ * ) depends on (ξ, ξ * ) variables and also on the function φ, and plays the role of a kernel of the integral operator (11). Hence, if for some special test function φ we have a simple explicit expression of K(φ; ξ, ξ * ), then we can solve the six-fold integral (11) instead of the twelve-fold integral (9). There is a simple, but important case, where the function K is known explicitly. In fact, denoting by ψ one of the collision invariants 1, ξ, |ξ| 2 , we have From a numerical point of view, the use of (11) for solving the integral of Q(f, f ) φ seems better than to use (9). In fact, for any fixed test function φ, firstly we can solve (10) and then the integral in the RHS of Eq. (11); i.e. we have split the numerical problem.
3. The numerical model. The numerical treatment of a kinetic equation by means of finite differences or elements requires a bounded domain for the velocity space. We introduce a suitable characteristic function in the kernel of the collision operator, such that there exists a bounded domain D so that if f (0, x, ξ) = 0 for every ξ ∈ D and x ∈ X, then f (t, x, ξ) = 0 for every ξ ∈ D, x ∈ X and for all time t. Let E be a positive real number. We define the function χ E : (13) It is immediate to see that It is evident that the new modified kernel W E satisfies the same properties of symmetry of true kernel W . Moreover, W E guarantees that if the particles, before the impact, have velocities such that |ξ| 2 + |ξ * | 2 ≤ E, then the velocities after the impact satisfy the inequality |ξ | 2 + |ξ * | 2 ≤ E. If we expect a distribution function f centred in u 0 , then we can choose the function χ E (ξ − u 0 , ξ * − u 0 ) instead of χ E (ξ, ξ * ). Now, since we can change the origin of the velocity space, so that ξ → ξ − u 0 , and write the Boltzmann equation in the new velocity reference frame, then, without loss of generality, we avoid the use of this characteristic function. Finally, it is possible to use a smooth function instead of the characteristic function χ E to modify the kernel and to obtain the same previous conclusions, but this is useless in our numerical approach. We can choose D = ξ ∈ R 3 : |ξ| 2 ≤ E , and we look for solutions of the Boltzmann equation vanishing for velocities outside the set D. We consider N measurable sets C α (α = 1, 2, , .., N ) such that We denote by χ α the characteristic function on the set C α .
It is obvious that simple and explicit expressions of K are unrealistic in the whole space D × D and for any φ. Fortunately, the DG methods requires only a finite set of test functions. We consider the 5N test functions where ξ 1 , ξ 2 , ξ 3 are the three components of the vector ξ. Therefore we will have 5N partial differential equations (3).

Basic equations.
In the framework of the discontinuous Galerkin method, a consistent approximation of the unknown is given using the same set of functions used as test functions. We denote by {ψ j : j = 0, 1.., 4} the ordered set 1, ξ 1 , ξ 2 , ξ 3 , |ξ| 2 . We now assume that The components of the five dimensional array η α (ξ) are functions, denoted by η α,i (ξ), which are linear combination of the collision invariants and such that for every i and j and for each cell C α . The vector functions g α (t, x) (α = 1, 2, .., N ) are the new unknowns instead of the distribution function f . It is immediate verify that where g α,j (t, x) are the components of g α (t, x). The physical meaning of the new unknowns is evident; for instance, g α,0 (t, x) is the density of particles at time t and position x having velocity belonging to the set C α . Inserting in Eq. (3) the approximation of the distribution function given by Eq. (16), where now φ(ξ) = ψ j (ξ) χ γ (ξ) (j = 0, 1, .., 4 and γ = 1, 2, ..., N ), we obtain the explicit set of partial differential equations with We note that the parameters V γ,ij and Λ αβγ,ijk are numerical constants, which depend only on the domain decomposition and the scattering kernel K. Proof. We define G is the part of the kernel K which derives from the gain term of the collision operator. Using the function G, we can write Now, we consider φ(ξ) = ψ(ξ) χ γ (ξ). If ξ and ξ * belong to D, then we have Now, we must take care that the assumption ξ and ξ * belonging to D is not sufficient to guarantee ξ and ξ * ∈ D; this holds, provided that we have the further hypothesis χ E (ξ, ξ * ) = 1. So, we have Hence This proves the theorem.
A simple consequence of this result is N γ=1 Λ αβγ,ijk = 0 for every α, β and i, j, k.

Using this, from Eqs. (19) follows
These correspond to the first five moment equations. In fact, by (18) follows and, taking again into account Eq. (16), we have We must choose the decomposition of the set D in order to have the explicit values of the coefficients V γ,ij , but some of these are easily available in any case. For instance, for j = 0, we have where ρ * and u * are the numerical density and velocity, defined by Remark 2. The property (22) is sufficient to obtain the moment equations (25). 4. The numerical parameters. The partial differential equations (19) contain many numerical parameters, which must be evaluated for solving the equations numerically. The parameters V γ,ij are very simple to calculate, provided that the cell have a regular and simple form. The true difficulty is given by the coefficients Λ αβγ,ijk . Now, we need only to find K with high accuracy to guarantee the conservation equation (25). The six-fold integral (21) can be calculate with low precision and, nevertheless, Eqs. (25) hold. For instance, the simplest way to find the numerical values of the parameters (21) is the following. For each cell C α we choose one velocity, which is denoted by ξ α . Then where µ(C α ) is the measure of the cell C α . The approximation of Λ αβγ,ijk given by formula (26) is sufficient to give Eqs. (25) provided that K(ψ j χ γ ; ξ α , ξ β ) are evaluated exactly. Of course, we expect a low accuracy using formula (26); so, we can consider to most refined quadrature formulas. For instance, the application of Gaussian quadrature rule for six dimensional integrals furnishes again Eqs. (25), with the assumption that K(ψ j χ γ ; ξ, ξ * ) are evaluated in the nodes exactly. Some simple formulas can be derived and used for finding the numerical values of K(ψ j χ γ ; ξ, ξ * ). In general, since the function ν is usually given and does not require further studies, we must consider only the operator G. It is useful to define Taking into account the equations linking the velocity before and after a collision, it is immediate to verify the following identities.

Remark 3.
To evaluate the kernel K(ψ j χ γ ; ξ, ξ * ) we need to solve the two-dimensional integrals (33) and (34). This can be made with great accuracy without great difficulties. Moreover, we can use Eqs. (32) as a corrector. For instance, firstly we find the values of A γ (ξ, ξ * ) for some fixed ξ and ξ * and for all γ, using some quadrature rules. Hence, the first Eq. (32) can be used for checking the goodness of the results and to define new refined values, for every γ, by means of the simple formula The second equation (32) can be used in a similar way also for the function B γ (ξ, ξ * ).
4.1. The isotropic case. If the kernel K depends only on V and we look for solutions of the Boltzmann equation with the distribution function f = f (t, |ξ|), then we can introduce a simple decomposition of the velocity space, where Of course, 0 = r 1 < R 1 = r 2 < R 2 = .....r N < R N = √ E and the origin is included in the first cell. In this case the function A γ (ξ, ξ * ) and B γ (ξ, ξ * ) can be found explicitly (see Appendix) Iγ r dr otherwise and B γ (ξ, ξ * ) = 0 if |ξ + ξ * | = 0; otherwise, Here, I γ is the interval (depending also on ξ and ξ * ) or the empty set given by . This case can be an useful benchmark.

5.
Conclusions and future work. We have proposed a numerical model for the Boltzmann equation based on the framework of the discontinuous Galerkin method. We derive a set of partial differential equations, which contain the first five moment equations for the approximate density, momentum and energy. This fact guarantees the usual conservation laws. The coefficients of the equations are numerical constants, which depend only on the decomposition of the velocity space and the scattering kernel of the collision operator. A numerical application of the proposed model requires the computation of a great number of integrals, but this work is required only one time. Indeed, if we have N cells in the velocity space, then, in principle, the number of parameters is of order N 3 (many parameters are null), but the partial differential equations are 5N . We made some very preliminary tests, where we have used simple quadrature rules, and the results seem promising. 6. Appendix. We derive a few formulas, which are present in the literature, both to introduce the notation and to make the paper self-consistent. Using the properties of Dirac distributions, we show how to reduce the order of the integral Here F is a function such that the following operations are allowed. Introducing the change of variables ξ = ξ − A, ξ * = ξ * + B, the integral (35) becomes R 3 R 3 F (ξ, ξ * , ξ − A, ξ * + B) δ(A − B) δ(|ξ| 2 + |ξ * | 2 − |ξ − A| 2 − |ξ * + B| 2 ) dA dB and, after a simple integration, Here we have used the vector V, which is defined in (6). Now, performing the further change of variable C = 1 2 V − A, it is easy to verify that the integral (36) is equal to The last simple reduction of (35) is obtained performing the integration with respect to the velocity ξ * . The result is 6.1. The total cross section. A simple application of formula (36) allows to reduce the integral whereâ is the unit vector of A. If V = 0 the integral (39) is null, otherwise we introduce a reference frame in R 3 and spherical coordinates such that A = r (sin ϕ cos θ, sin ϕ sin θ, cos ϕ) and A · V = r |V| cos ϕ .
Hence the integral (39) reduces to 6.2. Reduction of functions related to G. We use the formula (37). Since taking into account the Dirac distribution, we can replace |C| 2 with 1 4 |V| 2 and we have Now, it is easy to verify that |V · n| 2 = 1 2 |V| 2 − V · C. It is possible to avoid the square root, by noting that 1 2 |V| 2 = 1 4 |V| 2 + |C| 2 . Therefore V · n = |C − 1 2 V|. We have