Spectral - Lagrangian methods for Collisional Models of Non - Equilibrium Statistical States

# Spectral - Lagrangian methods for Collisional Models of Non - Equilibrium Statistical States

Irene M. Gamba Dept. of Mathematics & Institute of Computational Engineering and Sciences, University of Texas Austin Sri Harsha Tharkabhushanam Institute of Computational Engineering and Sciences, University of Texas Austin
###### Abstract

We propose a new spectral Lagrangian based deterministic solver for the non-linear Boltzmann Transport Equation for Variable Hard Potential (VHP) collision kernels with conservative or non-conservative binary interactions. The method is based on symmetries of the Fourier transform of the collision integral, where the complexity in its computing is reduced to a separate integral over the unit sphere . In addition, the conservation of moments is enforced by Lagrangian constraints. The resulting scheme, implemented in free space is very versatile and adjusts in a very simple manner, to several cases that involve energy dissipation due to local micro-reversibility (inelastic interactions) or elastic model of slowing down process. Our simulations are benchmarked with the available exact self-similar solutions, exact moment equations and analytical estimates for homogeneous Boltzmann equation for both elastic and inelastic VHP interactions. Benchmarking of the simulations involves the selection of a time self-similar rescaling of the numerical distribution function which is performed using the continuous spectrum of the equation for Maxwell molecules as studied first in [13] and generalized to a wide range of related models in [12]. The method also produces accurate results in the case of inelastic diffusive Boltzmann equations for hard-spheres (inelastic collisions under thermal bath), where overpopulated non-Gaussian exponential tails have been conjectured in computations by stochastic methods in [49; 26; 46; 35] and rigourously proven in [34] and [15].

###### keywords:
Spectral Method, Boltzmann Transport Equation, Conservative/ Non-conservative deterministic Method, Lagrangian optimization, FFT

## 1 Introduction

In a microscopic description of a rarefied gas, all particles are assumed to be traveling in a straight line with a fixed velocity until they enter into a collision. In such dilute flows, binary collisions are often assumed to be the main mechanism of particle interactions. The statistical effect of such collisions can be modeled by collision terms of the Boltzmann or Enskog transport equation type, where the kinetic dynamics of the gas are subject to the molecular chaos assumption. The nature of these interactions could be elastic, inelastic or coalescing. They could either be isotropic or anisotropic, depending on their collision rates as a function of the scattering angle. In addition, collisions are described in terms of inter-particle potentials and the rate of collisions is usually modeled as product of power laws for the relative speed and the differential cross section, at the time of the interaction. When the rate of collisions is independent of the relative speed, the interaction is referred to as of Maxwell type. When it corresponds to relative speed to a positive power less than unity, they are referred to as Variable Hard Potentials (VHP) and when the rate of collisions is proportional to the relative speed, it is referred to as hard spheres.

The Boltzmann Transport Equation (an integro-differential transport equation) describes the evolution of a single point probability distribution function which is defined as the probability of finding a particle at position with velocity (kinetic) at time . The mathematical and computational difficulties associated to the Boltzmann equation are due to the non local - non linear nature of the collision operator, which is usually modeled as a multi linear integral form in -dimensional velocity space and unit sphere .

From the computational point of view, of the well-known and well-studied methods developed in order to solve this equation is an stochastic based method called ”Direct Simulation Monte-Carlo” (DSMC) developed initially by Bird [2] and Nanbu [48] and more recently by [54; 55]. This method is usually employed as an alternative to hydrodynamic solvers to model the evolution of moments or hydrodynamic quantities. In particular, this method have been shown to converge to the solution of the classical Boltzmann equation in the case of mono atomic rarefied gases [57]. One of the main drawbacks of such methods is the inherent statistical fluctuations in the numerical results, which becomes very expensive or unreliable in presence of non-stationary flows or non equilibrium statistical states, where more information is desired about the evolving probability distribution. Currently, there is extensive work from Rjasanow and Wagner [55] and references therein, to determine accurately the high-velocity tail behavior of the distribution functions from DSMC data. Implementations for micro irreversible interactions such as inelastic collisions have been carefully studied in [35].

In contrast, a deterministic method computes approximations of the probability distribution function using the Boltzmann equation, as well as approximations to the observables like density, momentum, energy, etc.,. There are currently two deterministic approaches to the computations of non-linear Boltzmann, one is the well known discrete velocity models and the second a spectral based method, both implemented for simulations of elastic interactions i.e. energy conservative evolution. Discrete velocity models were developed by Broadwell [20] and mathematically studied by Illner, Cabannes, Kawashima among many authors [41; 42; 21]. More recently these models have been studied for many other applications on kinetic elastic theory in [7; 24; 44; 59; 39]. These models have not adapted to inelastic collisional problems up to this point according to our best knowledge.

Spectral based models, which are the ones of our choice in this work, have been developed by Pareschi, Gabetta and Toscani [32], and later by Bobylev and Rjasanow [17] and Pareschi and Russo [52]. These methods are supported by the ground breaking work of Bobylev [4] using the Fourier Transformed Boltzmann Equation to analyze its solutions in the case of Maxwell type of interactions. After the introduction of the inelastic Boltzmann equation for Maxwell type interactions and the use of the Fourier transform for its analysis by Bobylev, Carrillo and one of the authors here [6], the spectral based approach is becoming the most suitable tool to deal with deterministic computations of kinetic models associated with Boltzmann non-linear binary collisional integral, both for elastic or inelastic interactions. More recent implementations of spectral methods for the non-linear Boltzmann are due to Bobylev and Rjasanow [19] who developed a method using the Fast Fourier Transform (FFT) for Maxwell type of interactions and then for Hard-Sphere interactions [18] using generalized Radon and X-ray transforms via FFT. Simultaneously, L. Pareschi and B. Perthame [51] developed similar scheme using FFT for Maxwell type of interactions. Later, I. Ibragimov and S. Rjasanow [40] developed a numerical method to solve the space homogeneous Boltzmann Equation on a uniform grid for a Variable Hard Potential interactions with elastic collisions. This particular work has been a great inspiration for the current work and was one of the first initiating steps in the direction of a new numerical method.

We mention that, most recently, Filbet and Russo [27][28] implemented a method to solve the space inhomogeneous Boltzmann equation using the previously developed spectral methods in [52; 51]. Afore mentioned work in developing deterministic solvers for non-linear BTE have been restricted to elastic, conservative interactions. Finally, Mouhout and Pareschi [47] are currently studying the approximation properties of the schemes. Part of the difficulties in their strategy arises from the constraint that the numerical solution has to satisfy conservation of the initial mass. To this end, the authors propose the use of a periodic representation of the distribution function to avoid aliasing. There is no conservation of momentum and energy in  [28], [27] and  [47]. Both methods ( [28],  [27],  [47]), which are developed in 2 and 3 dimensions, do not guarantee the positivity of the solution due to the fact that the truncation of the velocity domain combined with the Fourier method makes the distribution function negative at times. This last shortcoming of the spectral approach remains in our proposed technique; however we are able to handle conservation in a very natural way by means of Lagrange multipliers. We also want to credit an unpublished calculation of V. Panferov and S. Rjasanow  [50] who wrote a method to calculate the particle distribution function for inelastic collisions in the case of hard spheres, but there were no numerical results to corroborate the efficiency of the method. Our proposed approach is slightly different and it takes a less number of operations to compute the collision integral.

Our current approach, based on a modified version of the work in [17] and [40], works for elastic or inelastic collisions and energy dissipative non-linear Boltzmann type models for variable hard potentials. We do not use periodic representations for the distribution function. The only restriction of the current method is that it requires that the distribution function at any time step be Fourier transformable. The required conservation properties of the distribution function are enforced through a Lagrange multiplier constrained optimization problem with the desired conservation quantities set as the constraints. Such corrections to the distribution function to make it conservative are very small but crucial for the evolution of the probability distribution function according to the Boltzmann equation.

This Lagrange optimization problem gives the freedom of not conserving the energy, independent of the collision mechanism, as long momentum is conserved. Such a technique plays a major role as it gives the option of computing energy dissipative solutions by just eliminating one constraint in the corresponding optimization problem. The current method can be easily implemented in any dimension. A novel aspect of the presented approach here lays on a new method that uses the Fourier Transform as a tool to simplify the computation of the collision operator that works, both for elastic and inelastic collisions. It is based on an integral representation of the Fourier Transform of the collision kernel as used in [17]. If is the number of discretizations in one direction of the velocity domain in -dimensions, the total number of operations required to solve for the collision integral are of the order of . And this number of operations remains the same for elastic/ inelastic, isotropic/ anisotropic VHP type of interactions. However, when the differential cross section is independent of the scattering angle, the integral representation kernel is further reduced by an exact closed integrated form that is used to save in computational number of operations to . This reduction is possible when computing hard spheres in d+2 dimensions or Maxwell type models in 2-dimensions. Nevertheless, the method can be employed without much changes for the other case. In particular the method becomes , where , the number of each angular discretizations is expected to be much smaller than used for energy discretizations. Such reduction in number of operations was also reported in  [28] with number of operations, where the authors are assuming N to be the total number of discretizations in the -dimensional space (i.e. our and of order of unity).

Our numerical study is performed for several examples of well establish behavior associated to solutions of energy dissipative space homogeneous collisional models under heating sources that secure existence of stationary states with positive and finite energy. We shall consider heating sources corresponding to randomly heated inelastic particles in a heat bath, with and without friction; elastic or inelastic collisional forms with anti-divergence terms due to dynamically (self-similar) energy scaled solutions [34; 15] and a particularly interesting example of inelastic collisions added to a slow down linear process that can be derived as a weakly coupled heavy/light binary mixture. On this particular case, when Maxwell type interactions are considered, it is shown that [13; 14; 12], on one hand dynamically energy scaled solutions exist, they have a close, explicit formula in Fourier space for a particular choice of parameters and their corresponding anti Fourier transform in probability space exhibit a singularity at the origin and power law high energy tails, while remaining integrable and with finite energy. On the other hand they are stable within a large class of initial states. We used this particular example to benchmark our computations by spectral methods by comparing the dynamically scaled computed solutions to the explicit one self similar one.

Convergence and error results of the Fourier Transform Lagrangian method, locally in time, are currently being developed [36], and it is expected that the proposed spectral approximation of the free space problem will have optimal algorithm complexity using the non-equispaced FFT as obtained by Greengard and Lin [38] for spectral approximation of the free space heat kernel.

Implementation of the space inhomogeneous case are also currently being considered. The spectral-Lagrangian scheme methodology proposed here can be extended to cases of Pareto tails, opinion dynamics and player games, where the evolution and asymptotic behavior of probabilities are studied in Fourier space as well. [53; 12].

The paper is organized as follows. In section 2, some preliminaries and description of the various approximated models associated with the elastic or inelastic Boltzmann equation are presented. In section 3, the actual numerical method is discussed with a small discussion on its discretization. In section 4, the special case of spatially homogeneous collisional model for a slow down process derived from a weakly coupled binary problem with isotropic elastic Maxwell type interactions is considered wherein an explicit solution is derived and shown to have power-like tails in some particular cases corresponding to a cold thermostat problem. Section 5 deals with the numerical results and examples. Finally in section 6, direction of future work is proposed along with a summary of the proposed numerical method.

## 2 Preliminaries

The initial value problem associated to space homogeneous Boltzmann Transport Equation modeling the statistical (kinetic) evolution of a single point probability distribution function for Variable Hard Potential (VHP) interactions is given by

 ∂∂tf(v,t) = Q(f,f)(v,t) = ∫w∈Rd,σ∈Sd−1[Jβf(′v,t)f(′w,t)−f(v,t)f(w,t)]B(|u|,μ)dσdw f(v,0) = f0(v), (2.1)

where the initial probability distribution is assumed integrable and is Jacobian of post with respect to pre collisional velocities which depend the local energy dissipation [22]. The problem may or may not have finite initial energy and the velocity interaction law, written in center of mass and relative velocity coordinates is

 u=v−w: the relative velocity% v′=v+β2(|u|σ−u), w′=w−β2(|u|σ−u) ,μ=cos(θ) = u⋅σ|u|: the cosine % of the scattering angle ,B(|u|,μ)=|u|λb(cosθ)with 0≤λ≤1,ωd−2∫π0b(cosθ)sind−2θdθ

where the parameter is the restitution coefficient corresponding from sticky to elastic interactions, where .

We denote by and the pre-collision velocities corresponding to and . In the case of micro-reversible (elastic) collisions one can replace and with and respectively in the integral part of (2.1). We assume the differential cross section function is integrable with respect to the post-collisional specular reflection direction in the dimensional sphere, referred as the Grad cut-off assumption, and that is renormalized such that

 ∫Sd−1b(u⋅σ|u|)dσ = ωd−2∫π0b(cosθ)sind−2θdθ (2.3) = ωd−2∫1−1b(μ)(1−μ2)(d−3)/2dμ =1,

where the constant is the measure of the dimensional sphere and the corresponding scattering angle is is defined by .

The parameter regulates the collision frequency as a function of the relative speed . It accounts for inter particle potentials defining the collisional kernel and they are referred to as Variable Hard Potentials (VHP) whenever , Maxwell Molecules type interactions (MM) for and Hard Spheres (HS) for . The Variable Hard Potential collision kernel then takes the following general form:

•  B(|u|,μ)=Cλ(σ)|u|λ, (2.4)

with for Maxwell type of interactions; for Hard Spheres. In addition, if is independent of the scattering angle we call the interactions isotropic. Otherwise we referred to them as anisotropic Variable Hard Potential interactions.

For classical case of elastic collisions, it has been established that the Cauchy problem for the space homogeneous Boltzmann equation has a unique solution in the class of integrable functions with finite energy (i.e. ), it is regular if initially so, and converges in to the Maxwellian distribution associated to the -moments of the initial state . In addition, if the initial state has Maxwellian decay, this property will remain with a Maxwellian decay globally bounded in time ([33]), as well as all derivatives if initial so (see [1]).

Depending on their nature, collisions either conserve density, momentum and energy (elastic) or density and momentum (inelastic) or density (elastic - linear Boltzmann operator), depending on the number of collision invariants the operator has. In the case of the classical Boltzmann equation for rarefied (elastic) mono-atomic gases, the collision invariants are exactly , that is, according to the Boltzmann theorem, the number of polynomials in velocity space that generate , with . In particular, one obtains the following conserved quantities

 density   ρ(t) =∫v∈Rdf(v,t)dv momentum   m(t) = ∫v∈Rdvf(v,t)dv (2.5) energy   E(t) = 12ρ(t)∫v∈Rd|v|2f(v,t)dv.

Of significant interest from the statistical view point are the evolution of moments or observables, at all orders. They are defined by the dynamics of the corresponding time evolution equation for the velocity averages, given by

 (2.6)

where, the standard symmetric tensor product of with itself, times. Thus, according to (2), for the classical elastic Boltzmann equation, the first moments are conserved, meaning, for ; and . Higher order moments or observables of interest are

 Momentum Flow   M2(t) = ∫RdvvTf(v,t)dvEnergy Flow   r(t) = 12ρ(t)∫Rdv|v|2f(v,t)dv Bulk Velocity   V(t) = m(t)ρ(t) Internal Energy   E(t) = 12ρ(tr(M2)−ρ|V|2)Temperature   T(t) = 2E(t)kd (2.7)

with Boltzmann constant.

We finally point out that, in the case of Maxwell molecules (), it is possible to write recursion formulas for higher order moments of all orders ([5] for the elastic case, and [6] in the inelastic case) which, in the particular case of isotropic solutions depending only on , take the form

 mn(t) = ∫Rd|v|2nf(v,t)dv= e−λntmn(0)+n−1∑k=112(n+1)(2n+22k+1)Bβ(k,n−k)∫t0mk(τ)mn−k(τ)e−λn(t−τ)dτ;withλn=1−1n+1[β2n+n∑k=0(1−β)2k],Bβ(k,n−k) = β2k∫10sk(1−β(2−β)s)n−kds, (2.8)

for where and .

### 2.1 Boltzmann collisional models with heating sources

A collisional model associated to the space homogeneous Boltzmann transport equation (2.1) with grad cutoff assumption (2.2), can be modified in order to accommodate for an energy or ‘heat source’ like term , where is a differential or integral operator. In these cases, it is possible to obtain stationary states with finite energy as for the case of inelastic interactions. In such general framework, the corresponding initial value problem model is

 ∂∂tf(v,t) = ζ(t)Q(f,f)(v,t)+G(f(t,v)), f(v,0) = f0(v), (2.9)

where the collision operator is as in (2.1) and models a ‘heating source’ due to different phenomena. The term may represent a mean field approximation that allows from proper time rescaling. See [6] and [15] for several examples for these type of models and additional references.

Following the work initiated in [15] and [14] on Non-Equilibrium Stationary States (NESS), our computational approach we shall present several computational simulations of non-conservative models for either elastic or inelastic collisions associated to (2.9) of the Boltzmann equation with ‘heating’ sources. In all the cases we have addressed one can see that stationary states with finite energy are admissible, but they may not be Maxwellian distributions. Of this type of model we show computational output for three cases. First one is the pure diffusion thermal bath due to a randomly heated background [58; 49; 34], in which case

 G1(f)=μΔf, (2.10)

where is a constant. The second example relates to self-similar solutions of equation (2.9) for [45; 25], but dynamically rescaled by

 f(v,t)=1vd0(t)~f(~v(v,t),~t(t)),~v=vv0(t), (2.11)

where

 v0(t)=(a+ηt)−1, ~t(t)=1ηln(1+ηat),a,η>0. (2.12)

Then, the equation for coincides (after omitting the tildes) with equation (2.9),for

 G2(f)=−ηdiv(vf),η>0. (2.13)

Of particular interest of dynamical time-thermal speed rescaling is the case of collisional kernels corresponding to Maxwell type of interactions. Since the second moment of the collisional integral is a linear function of the energy, the energy evolves exponentially with a rate proportional to the energy production rate, that is

 ddtE(t) = λ0E(t),or equivalently E(t)=E(0)eλ0t, (2.14)

with the energy production rate. Therefore the corresponding rescaled variables and equations (2.11) and (2.9),(2.13) for the study of long time behavior of rescaled solutions are

 f(v,t)=E−d2(t)~f(vE12(t)) = (E(0)eλ0t)−d2~f(v(E(0)eλ0t)−12), (2.15)

and satisfies the self-similar equation (2.9)

 G2′(f)=−λ0xfx,where  x=vE−12(t)  is the self-similar variable. (2.16)

We note that it has been shown that these dynamically self-similar states are stable under very specific scaling for a large class of initial states [12].

The last source type we consider is given by a model, related to weakly coupled mixture modeling slowdown (cooling) process [14] given by an elastic model in the presence of a thermostat given by Maxwell type interactions of particles of mass having the Maxwellian distribution

 MT(v)=m(2πT)d/2e−m|v|22T ,

with a constant reference background or thermostat temperature (i.e. the average of and ). Define

 QL(f)(v,t)≐∫w∈Rd,σ∈Sd−1BL(|u|,μ)f(′v,t)MT(′w)−f(v,t)MT(w)] dσdw. (2.17)

Then the corresponding evolution equation for is given by

 ∂∂tf(v,t) = Q(f,f)(v,t)+ΘQL(f)(v,t) f(v,0) = f0(v). (2.18)

where , defined as in (2.1), is the classical collision integral for elastic interactions (i.e. ), so it conserves density, momentum and energy. The second integral term in (2.18) is a linear collision integral which conserves just the density (but not momentum or energy) since

 u =v−w % the relative velocityv′ =v+mm+1(|u|σ−u), w′=w−1m+1(|u|σ−u). (2.19)

The coupling constant depends on the initial density, the coupling constants and on . The collision kernel of the linear part may not be the same as the one for the non-linear part of the collision integral, however we assume that the Grad cut-off assumption (2.3) is satisfied and that, in order to secure mass preservation, the corresponding differential cross section functions and , the non-linear and linear collision kernels respectively, satisfy the renormalized condition

 ∫Sd−1bN(u⋅σ|u|)+ΘbL(u⋅σ|u|)dσ = 1+Θ. (2.20)

This last model describes the evolution of binary interactions of two sets of particles, heavy and light, in a weakly coupled limit, where the heavy particles have reached equilibrium. The heavy particles set constitutes the background or thermostat for the second set of particles. It is the light particle distribution that is modeled by (2.18). Indeed, corresponds to all the collisions which the light particles have with each other and the second linear integral term corresponds to collisions between light and heavy particles at equilibrium given by a classical distribution . In this binary -dimensional, mixture scenario, collisions are assumed to be isotropic, elastic and the interactions kernels of Maxwell type.

In the particular case of equal mass (i.e. ), the model is of particular interest for the development of numerical schemes and simulations benchmarks. Even though the local interactions are reversible (elastic), it does not conserve the total energy. In such a case, there exists a special set of explicit, in spectral space, self-similar solutions which are attractors for a large class of initial states. When considering the case of Maxwell type of interactions in three dimensions i.e. with a cooling background process corresponding to a time temperature transformation, such that as the models have self similar asymptotics [14; 12] for a large class of initial states. Such long time asymptotics corresponding to dynamically scaled solutions of (2.18), in the form of (2.16), yields interesting behavior in for large time, converging to states with power like decay tails in . In particular, such solution of (2.18) will lose moments as time grows, even if the initial state has all moments bounded. (see [14; 12] for the analytical proofs).

### 2.2 Collision Integral Representation

One of the pivotal points in the derivation of the spectral numerical method for the computation of the non-linear Boltzmann equation lays in the representation of the collision integral in Fourier space by means of the weak form. Since for a suitably regular test function , the weak form of the collision integral takes the form (suppressing the time dependence in f)

 ∫v∈RdQ(f,f)ψ(v)dv=∫(w,v)∈Rd×Rd,σ∈Sd−1f(v)f(w)B(|u|,μ)[ψ(v′)−ψ(v)]dσdwdv, (2.21)

with . In particular, taking , where is the Fourier variable, we get the Fourier Transform of the collision integral through its weak form as

 ˆQ(ζ) = 1(√2π)d∫v∈RdQ(f,f)e−iζ⋅vdv (2.22) = ∫(w,v)∈Rd×Rd,σ∈Sd−1f(v)f(w)B(|u|,μ)(√2π)d[e−iζ⋅v′−e−iζ⋅v]dσdwdv.

We will use - the Fourier transform and for the classical inverse Fourier transform. Plugging in the definitions of collision kernel (which in the case of isotropic collisions would just be the Variable Hard Potential collision kernel) and

 ˆQ(ζ) = 1(√2π)d∫u∈RdGλ,β(u,ζ)∫v∈Rdf(v)f(v−u)e−iζ⋅vdvdu (2.23) = ∫u∈RdGλ,β(u,ζ)[f(v)f(v−u)]ˆ du,

where

 Gλ,β(u,ζ) = ∫σ∈Sd−1Cλ(σ)|u|λ[e−iβ2ζ⋅(|u|σ−u))−1]dσ (2.24) = |u|λ[eiβ2ζ⋅u∫σ∈Sd−1Cλ(σ)e−iβ2|u|ζ⋅σdσ−ω2].

Note that (2.24) is valid for both isotropic and anisotropic interactions. For the former type, a simplification ensues due to the fact the is independent of :

 Gλ,β(u,ζ)=Cλωd−2|u|λ[eiβ2ζ.usinc(β|u||ζ|2)−1]. (2.25)

Thus, it is seen that the dependence on i.e. the integration over the unit sphere is completely done independently and there is actually a closed form expression for this integration, given by (2.25) in the case of isotropic collisions. In the case of anisotropic collisions, the dependence of on is again isolated into a separate integral over the unit sphere as given in (2.24). The above expression can be transformed for elastic collisions into a form suggested by Rjasanow and Ibragimov in their paper [40]. The corresponding expression for anisotropic collisions is given by (2.24).

Further simplification of (2.23) is possible by observing that the Fourier transform inside the integral can be written in terms of the Fourier transform of since it can also be written as a convolution of the Fourier transforms. Let

 ˆQ(ζ) = ∫u∈RdGλ,β(u,ζ)[f(v)h(v)]ˆ du=∫u∈RdGλ,β(u,ζ)1(√2π)d(^f∗^h)(ζ)du (2.26) = ∫u∈RdGλ,β(u,ζ)1(√2π)d∫ξ∈Rd^f(ζ−ξ)^f(ξ)e−iξ⋅udξdu = 1(√2π)d∫ξ∈Rd^f(ζ−ξ)^f(ξ)^Gλ,β(ξ,ζ)dξ,

where . Let , For ,

 ^Gλ,β(ξ,ζ) = ∫r∫er2G(re,ζ)e−irξ⋅ededr = 16π2Cλ∫rrλ+2[sinc(2β|ζ|2)sinc(r|ξ−β2ζ|)−sinc(r|ξ|)]dr.

Since the domain of computation is restricted to ,

 ^Gλ,β(ξ,ζ)=16π2Cλ∫2√3L0rλ+2[sinc(2β|ζ|2)sinc(r|ξ−β2ζ|)−sinc(r|ξ|)]dr. (2.27)

A point worth noting is that the above formulation (2.26) results in number of operations, where is the number of discretizations in each velocity direction. Also, exploiting the symmetric nature in particular cases of the collision kernel one can reduce the number of operations to .

## 3 Numerical Method

### 3.1 Discretization of the Collision Integral

Coming to the discretization of the velocity space, it is assumed that the two interacting velocities and the corresponding relative velocity

 v,w,and w ∈[−L,L)d, (3.1) whileζ ∈[−Lζ,Lζ)d, (3.2)

where the velocity domain is chosen such that through an assumption that . For a sufficiently large , the computed distribution will not lose mass, since the initial momentum is conserved (there is no convection in space homogeneous problems), and is renormalized to zero mean velocity. We assume a uniform grid in the velocity space and in the fourier space with and as the respective grid element sizes. and are chosen such that , where = number of discretizations of and in each direction.

### 3.2 Time Discretization

To compute the actual particle distribution function, one needs to use an approximation to the time derivative of . For this, a second-order Runge-Kutta scheme or a Euler forward step method were used. Since a non-dimensional Boltzmann equation is computed, for numerical computations the value of time step is chosen such that it corresponds in its dimensional form to times the time between consecutive collisions (which depends on the collision frequency). During the standard process of non-dimensionalization of the Boltzmann Equation, such a reference quantity (time between collisions) comes up. With time discretizations taken as , the discrete version of the Runge-Kutta scheme is given by

 f0(vj) = f0(vj) ~f(vj) = ftn(vj)+dt2Qλ,β[ftn(vj),ftn(vj)] ftn+1(vj) = ftn(vj)+dtQλ,β[~f(vj),~f(vj)] .

The corresponding Forward Euler scheme with smaller time step is given by

 ~f(vj)=ftn(vj)+dtQ(ftn,ftn). (3.4)

### 3.3 Conservation Properties - Lagrange Multipliers

Since the calculation of involves computing Fourier Transforms with respect to , we extensively use Fast Fourier Transform. Note that the total number of operations in computing the collision integral reduces to the order of for (2.23) and for (2.26). Observe that, choosing , the proposed scheme works for both elastic and inelastic collisions. As a note, the method proposed in the current work can also be extended to lower dimensions in velocity space.

In the current work, due to the discretizations and the use of Fourier Transform, the accuracy of the proposed method relies heavily on the size of the grid and the number of points taken in each velocity/ Fourier space directions. Because of this it is seen that the computed does not really conserve quantities it is supposed to i.e. for elastic collisions, for Linear Boltzmann Integral and for inelastic collisions. Even though the difference between the computed (discretized) collision integral and the continuous one is not great, it is nevertheless essential that this issue be resolved. To remedy this, a simple constrained Lagrange multiplier method is employed where the constraints are the required conservation properties. Let , the total number of discretizations of the velocity space. Assume that the classical Boltzmann collision operator is being computed. So and are conserved. Let be the integration weights where . Let

 ~f=(~f1~f2..~fM)T

be the distribution vector at the computed time step and

 f=(f1f2..fM)T

be the corrected distribution vector with the required moments conserved. Let

 C(d+2)×M=⎛⎜ ⎜⎝ωjviωj|vj|2ωj⎞⎟ ⎟⎠

and

 a(d+2)×1=(ρm1m2m3e)T

be the vector of conserved quantities. Using the above vectors, the conservation can be written as a constrained optimization problem:

 (∗){∥~f−f∥22→minCf=a;C∈Rd+2×M,f∈RM,a∈Rd+2.

To solve , one can employ the Lagrange multiplier method. Let be the Lagrange multiplier vector. Then the scalar objective function to be optimized is given by

 L(f,λ)=M∑j=1|~fj−fj|2+λT(Cf−a). (3.5)

Equation (3.5) can actually be solved explicitly for the corrected distribution value and the resulting equation of correction be implemented numerically in the code. Taking the derivative of with respect to and i.e. gradients of ,

 ∂L∂fj = 0;j=1,...,M ⇒ f = ~f+12CTλ. (3.6)

And

 ∂L∂λ1 = 0;i=1,...,d+2 ⇒ Cf = a (3.7)

i.e. retrieves the constraints.
Solving for ,

 CCTλ=2(a−C~f). (3.8)

Now is symmetric and because is the integration matrix, is positive definite. By linear algebra, the inverse of exists. In particular one can compute the value of by

 λ=2(CCT)−1(a−C~f). (3.9)

Substituting into (3.6),

 f=~f+CT(CCT)−1(a−C~f). (3.10)

Using equation for forward Euler scheme (3.4), the complete scheme is given by () :

 ~fj=fnj+dtQ(fnj,fnj)fn+1j=~fj+CT(CCT)−1(a−C~f). (3.11)

So,

 fn+1j = fnj+dtQ(fnj,fnj)+CT(CCT)−1(a−C~f) = fnj+dtQ(fnj,fnj)+CT(CCT)−1(a−a−dtCQ(fnj,fnj)) = fnj+dtQ(fnj,fnj)−dtCT(CCT)−1CQ(fnj,fnj) = fnj+dt[I−CT(CCT)−1C]Q(fnj,fnj),

with - identity matrix. Letting with - identity matrix, one obtains

 fn+1j=fnj+dtΛN(C)Q(fnj,fnj), (3.13)

where we expect the required observables are conserved and the solution approaches a stationary state, since .

Identity (3.13) summarizes the whole conservation process. As described previously, setting the conservation properties as constraints to a Lagrange multiplier optimization problem ensures that the required observables are conserved. Also, the optimization method can be extended to have the distribution function satisfy more (higher order) moments from (2.8). In this case, will include entries of from (5.1).

We point out that for the linear Boltzmann collision operator used in the mixture problem conserves density and not momentum(unless one computes isotropic solutions) and energy. For this problem, the constraint would just be the density equation. For inelastic collisions, density and momentum are conserved and for this case the constraint would be the energy and momentum equations. And for the elastic Boltzmann operator, all three quantities (density, momentum and energy) are conserved and thus they become the constraints for the optimization problem. The behavior of the conservation correction for Pseudo-Maxwell Potentials for Elastic collisions will be numerically studied in the numerical results section. This approach of using Lagrangian constraints in order to secure moment preservation differs from the one proposed in  [27][28] for spectral solvers.

## 4 Self-Similar asymptotics for a general elastic or inelastic BTE of Maxwell type or the cold thermostat problem - power law tails

As mentioned in introduction, a new interesting benchmark problem for our scheme is that of a dynamically scaled solutions or self-similar asymptotics. More precisely, we present simulation where the computed solution in properly scaled time approaches a self similar solution. This is of interest because of the power tail behavior i.e. higher order moments of the computed solution are bounded. For the completeness of this presentation, the analytical description of such asymptotics is given in the following two sub sections.

### 4.1 Self-Similar Solution for a non-negative Thermostat Temperature

We consider the Maxwell type equation from (2.18) related to a space homogeneous model for a weakly coupled mixture modeling slowdown process. The content of this section is dealt in detail in [14] for a particular choice of zero background temperature (cold thermostat). For the sake of brevity, we refer to [14] for details. However, a slightly more general form of the self-similar solution for non zero background temperature is derived here from the zero background temperature solution. Without loss of generality for our numerical test, we assume the differential cross sections for collision kernel of the linear and , the corresponding one for the nonlinear part, are the same, both denoted by , satisfying the Grad cut-off conditions (2.3). In particular, condition (2.20) is automatically satisfied.

In [14], Fourier transform of the isotropic self-similar solution associated to problem in (2.18) will take the form:

 ϕ(x,t) = ψ(xe−μt) = 1−a(xe−μt)p,   as   xe−μt→0,    withp≤1, (4.1)

where and and are related by

 μ = 23p2    and   Θ = (3p+1)(2−p)3p2.

Note that corresponds to initial states with finite energy. It was shown in [14] for (i.e. cold thermostat), the Fourier transform of the self-similar, isotropic solutions of (2.18) is given by

 ϕ(x,t) = 4π∫∞01(1+s2)2e−xe−2t3as2ds, (4.2)

and its corresponding inverse Fourier transform, both for , and (as computed in [14]) is given by

 fss0(|v|,t)=etF0(|v|et/3)withF0(|v|)=4π∫∞01(1+s2)2e−|v|2/2s2(2πs2)32ds. (4.3)

Remark: It is interesting to observe that, as computed originally in [9], for or in (4.1) yields , and one can construct explicit solutions to the elastic BTE with infinite initial energy. It is clear now that in order to have self-similar explicit solutions with finite energy one needs to have this weakly couple mixture model for slowdown processes, or bluntly speaking the linear collisional term added to the elastic energy conservative operator.

Finally, in order to recover the self-similar solution for the original equilibrium positive temperature (i.e. hot thermostat case) for the linear collisional term, we denote, including time dependence for convenience,

 ϕ0(x,t) = ϕ(x,t)Thermostat=0   and    ϕT(x,t) = ϕ(x,t)Thermostat=T (4.4) so that   ϕT(x,t) = ϕ0(x,t)e−Tx.

Note that the solution constructed in (4.2) is actually . Then the self-similar solution for non zero background temperature, denoted by satisfies

 ϕT(k,t) = 4π∫∞0e−|k|2e−2t/3as2/21(1+s2)2e−|k|2T/2ds (4.5) = 4π∫∞0e−|k|2[e−2t/3as2+T]/21(1+s2)2ds.

In particular, let then, taking the inverse Fourier Transform, we obtain the corresponding self-similar state, according to (2.15), in probability space

 fssT(|v|,t)=etFT(|v|e