Spectral - Lagrangian methods for Collisional Models of Non - Equilibrium Statistical States
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  and generalized to a wide range of related models in . 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  and .
keywords:Spectral Method, Boltzmann Transport Equation, Conservative/ Non-conservative deterministic Method, Lagrangian optimization, FFT
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  and Nanbu  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 . 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  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 .
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  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 , and later by Bobylev and Rjasanow  and Pareschi and Russo . These methods are supported by the ground breaking work of Bobylev  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 , 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  who developed a method using the Fast Fourier Transform (FFT) for Maxwell type of interactions and then for Hard-Sphere interactions  using generalized Radon and X-ray transforms via FFT. Simultaneously, L. Pareschi and B. Perthame  developed similar scheme using FFT for Maxwell type of interactions. Later, I. Ibragimov and S. Rjasanow  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 ,  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  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 ,  and . Both methods ( , , ), 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  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  and , 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 . 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  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 , 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  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.
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
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 . 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
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
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:
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
(), as well as all derivatives if initial so (see
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
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
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
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 ( for the elastic case, and  in the inelastic case) which, in the particular case of isotropic solutions depending only on , take the form
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
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  and  for
several examples for these type of models and additional references.
Following the work initiated in  and  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
Then, the equation for coincides (after omitting the tildes) with equation (2.9),for
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
and satisfies the self-similar equation (2.9)
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 .
The last source type we consider is given by a model, related to weakly coupled mixture modeling slowdown (cooling) process  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
Then the corresponding evolution equation for is given by
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
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
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)
with . In particular, taking , where is the Fourier variable, we get the Fourier Transform of the collision integral through its weak form as
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
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 :
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 . 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
where . Let , For ,
Since the domain of computation is restricted to ,
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
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
The corresponding Forward Euler scheme with smaller time step is given by
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
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
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 ,
i.e. retrieves the constraints.
Solving for ,
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
Substituting into (3.6),
Using equation for forward Euler scheme (3.4), the complete scheme is given by () :
with - identity matrix. Letting with - identity matrix, one obtains
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 ,  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  for a particular choice of zero background temperature (cold thermostat). For the sake of brevity, we refer to  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 , Fourier transform of the isotropic self-similar solution associated to problem in (2.18) will take the form:
where and and are related by
and its corresponding inverse Fourier transform, both for , and (as computed in ) is given by
Remark: It is interesting to observe that, as computed
originally in , 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
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,
Note that the solution constructed in (4.2) is actually . Then the self-similar solution for non zero background temperature, denoted by satisfies
In particular, let then, taking the inverse Fourier Transform, we obtain the corresponding self-similar state, according to (2.15), in probability space