Spectral  Lagrangian methods for Collisional Models of Non  Equilibrium Statistical States
Abstract
We propose a new spectral Lagrangian based deterministic solver for the nonlinear Boltzmann Transport Equation for Variable Hard Potential (VHP) collision kernels with conservative or nonconservative 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 microreversibility (inelastic interactions) or elastic model of slowing down process. Our simulations are benchmarked with the available exact selfsimilar 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 selfsimilar 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 hardspheres (inelastic collisions under thermal bath), where overpopulated nonGaussian 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/ Nonconservative deterministic Method, Lagrangian optimization, FFT1 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 interparticle 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 integrodifferential 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 wellknown and wellstudied methods developed in order to solve this equation is an stochastic based method called ”Direct Simulation MonteCarlo” (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 nonstationary 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 highvelocity 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 nonlinear 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 nonlinear binary collisional integral, both for elastic or inelastic interactions. More recent implementations of spectral methods for the nonlinear 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 HardSphere interactions [18] using generalized Radon and Xray 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 nonlinear 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 nonlinear 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
2dimensions. 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
antidivergence terms due to dynamically (selfsimilar) 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 nonequispaced 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 spectralLagrangian 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 powerlike 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
(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
(2.2) 
where the parameter is the restitution coefficient
corresponding from sticky to elastic interactions, where .
We denote by and the precollision velocities
corresponding to and . In the case of microreversible
(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 postcollisional
specular reflection direction in the dimensional
sphere, referred as the Grad cutoff assumption, and that
is renormalized such that
(2.3)  
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:

(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)
monoatomic 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
(2.5)  
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
(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
(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
(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 NonEquilibrium Stationary States (NESS), our
computational approach we shall present several computational
simulations of nonconservative 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
(2.10) 
where is a constant. The second example relates to selfsimilar solutions of equation (2.9) for [45; 25], but dynamically rescaled by
(2.11) 
where
(2.12) 
Then, the equation for coincides (after omitting the tildes) with equation (2.9),for
(2.13) 
Of particular interest of dynamical timethermal 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
(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
(2.15) 
and satisfies the selfsimilar equation (2.9)
(2.16) 
We note that it has been shown that these dynamically selfsimilar
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
with a constant reference background or thermostat temperature (i.e. the average of and ). Define
(2.17) 
Then the corresponding evolution equation for is given by
(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
(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 nonlinear part of the collision integral, however we assume that the Grad cutoff assumption (2.3) is satisfied and that, in order to secure mass preservation, the corresponding differential cross section functions and , the nonlinear and linear collision kernels respectively, satisfy the renormalized condition
(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,
selfsimilar 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 nonlinear 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)
(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
(2.22)  
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
(2.23)  
where
(2.24)  
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 :
(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
(2.26)  
where . Let , For ,
Since the domain of computation is restricted to ,
(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
(3.1)  
(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 secondorder RungeKutta scheme or a Euler forward step method were used. Since a nondimensional 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 nondimensionalization of the Boltzmann Equation, such a reference quantity (time between collisions) comes up. With time discretizations taken as , the discrete version of the RungeKutta scheme is given by
The corresponding Forward Euler scheme with smaller time step is given by
(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
be the distribution vector at the computed time step and
be the corrected distribution vector with the required moments conserved. Let
and
be the vector of conserved quantities. Using the above vectors, the conservation can be written as a constrained optimization problem:
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
(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 ,
(3.6) 
And
(3.7) 
i.e. retrieves the constraints.
Solving for ,
(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
(3.9) 
Substituting into (3.6),
(3.10) 
Using equation for forward Euler scheme (3.4), the complete scheme is given by () :
(3.11) 
So,
with  identity matrix. Letting with  identity matrix, one obtains
(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 PseudoMaxwell
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 SelfSimilar 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 selfsimilar 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 SelfSimilar Solution for a nonnegative 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 selfsimilar 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 cutoff conditions (2.3). In particular, condition (2.20) is automatically satisfied.
In [14], Fourier transform of the isotropic selfsimilar solution associated to problem in (2.18) will take the form:
(4.1) 
where and and are related by
Note that corresponds to initial states with finite energy. It was shown in [14] for (i.e. cold thermostat), the Fourier transform of the selfsimilar, isotropic solutions of (2.18) is given by
(4.2) 
and its corresponding inverse Fourier transform, both for , and (as computed in [14]) is given by
(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 selfsimilar
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 selfsimilar solution for the
original equilibrium positive temperature (i.e. hot thermostat
case) for the linear collisional term, we denote, including time
dependence for convenience,
(4.4)  
Note that the solution constructed in (4.2) is actually . Then the selfsimilar solution for non zero background temperature, denoted by satisfies
(4.5)  
In particular, let then, taking the inverse Fourier Transform, we obtain the corresponding selfsimilar state, according to (2.15), in probability space