About an H-theorem for systems with non-conservative interactions
We exhibit some arguments in favour of an H-theorem for a generalization of the Boltzmann equation including non-conservative interactions and a linear Fokker-Planck-like thermostatting term. Such a non-linear equation describing the evolution of the single particle probability of being in state at time , is a suitable model for granular gases and is indicated here as Boltzmann-Fokker-Planck (BFP) equation. The conjectured H-functional, which appears to be non-increasing, is with , in analogy with the H-functional of Markov processes. The extension to continuous states is straightforward. A simple proof can be given for the elastic BFP equation. A semi-analytical proof is also offered for the BFP equation for so-called inelastic Maxwell molecules. Other evidence is obtained by solving particular BFP cases through numerical integration or through “particle schemes” such as the Direct Simulation Monte Carlo.
The H-theorem is a consequence of the Boltzmann equation (whose validity, at least for diluted gases, is now rather well understood) and has a great relevance in statistical mechanics because it leads to two basic results Boltzmann (1872). First it provides a dynamical proof of the stationary Maxwell-Boltzmann distribution for the particle velocity and second it states that the approach of the probability distribution toward is monotonic. The latter result can be interpreted as an instance of the 2-nd law of the thermodynamics.
The literature about this subject is enormous and it is pretty impossible to enter into details. A discussion about the deep and intriguing issues concerning the reversibility and recurrence paradoxes, as well the rigorous derivation of the Boltzmann equation in the Grad-Boltzmann limit, can be found in Refs. Cohen (Editor) (1981); Cercignani (2006); Cercignani et al. (1994); Castiglione et al. (2008) .
After the seminal results by Boltzmann obtained employing his transport equation describing the behavior of diluted gases, similar H theorems had been obtained for systems described by Markov processes, that is systems ruled by master equations or Fokker-Planck equations van Kampen (1961). There is however an important difference with respect to the original Boltzmann H-theorem, where the functional depends only on and the knowledge of the asymptotic stationary probability is not required: for Markov processes the functional depends both on and the asymptotic stationary probability which must be determined, see next Section for details.
One of the basic feature of the Boltzmann equation is the presence of bilinear terms, describing the binary collisions. Although the Boltzmann equation had been originally derived for systems with conservative dynamics, it is not difficult, at least at formal level, to write down a similar evolution equation for the one-particle probability distribution for some dissipative systems (e.g. diluted granular gases). Due to the dissipative nature of granular gases, it is necessary to introduce an external mechanism pumping energy into the system, in order to have a statistical stationary state (mathematically a non trivial stationary probability distribution) . As a result the evolution equation for includes a linear term representing the coupling to the ”external bath” and a bilinear term accounting for binary collisions Puglisi et al. (1998).
In the classical derivation of Boltzmann a key feature used to prove the H-theorem is the presence of time reversal symmetry. The absence of such a property in dissipative systems is one of the technical difficulties to derive a more general H-theorem. Since the physical importance and the mathematical relevance the study of relaxation toward invariant probability, in terms of entropic functions, is an interesting issue which attracted the interest of many scientists Bena et al. (2006); Arnold et al. (2004); Toscani (1999); Garbazewski (2005).
Previous studies showed that both the stationary and dynamical statistical features are the combined results of bath and collisions Gradenigo et al. (2011). We will see that in granular systems H-theorems, if any, must be the outcomes of both the linear and bilinear part of the evolution equation. The control of such contributions is not easy: we are able to show an H-theorem in the particular limit of ”elastic granular” gases in an external bath. Such a system, although quite artificial, presents nontrivial dynamical features. In addition we give some semi-analytical treatment of a Maxwell model of granular gas, and detailed numerical computations supporting our idea.
The paper is organized as follows. Section 2 is a quick summary of known different H-theorems for the Boltzmann equation and Markov processes. In Section 3 we present a kinetic description of diluted granular gases and some results about the H-theorems for such a systems in particular limits. Section 4 treats numerical simulation of the granular gases in different regimes. In Section 5 the reader can find some conclusion. The Appendices are devoted to few technical details.
Ii Monotonic approach to invariant probability
For completeness we briefly review, in this Section, two known H-theorems.
ii.1 Boltzmann equation for elastic isolated gases
For the sake of notational simplicity we consider the case where the states of the system are discrete. Let call the probability of observing a state at time . Assuming the validity of the molecular chaos assumption its evolution is governed by the following non linear equation:
where denotes the transition rate of the ”collision” from the states and to the states and . By invariance under time inversion one has the following property
so that (1) takes the form
using such a structure it is easy to show the H-theorem, i.e. the function
is monotonically decreasing:
and reaches its minimum when the PdF corresponds to the Maxwell-Boltzmann distribution.
ii.2 Markov processes
An H-theorem also holds under rather general hypothesis van Kampen (1961) for processes governed by a Master equation of the form
Here is the transition rate from the state to the state .
In this case indicating with the invariant probabilities , which are the solutions of the following equation:
one has that the function
is non increasing, i.e.
and attains its minimum when . Let us note that is the conditional entropy, also known as Kullback-Leibler or relative entropy Cover and Thomas (2006).
In the case of continuous variables it is sufficient to replace the probability with the probability density (where is a state in a continuous space, e.g. velocity of a particle) and the sums with the integrals, the function becomes
ii.3 Some remarks
Let us notice that is an intrinsic property, that is switching to a new variable , if the transformation between and is invertible, is invariant, whereas is not Cover and Thomas (2006); Maynar and Trizac ().
Moreover, for certain classes of Markov processes it is possible to show Mackey and Tyran-Kamińska () a result stronger than the monotonic behavior of : there exists a constant such that
Finally, one may wonder if is non-increasing also in the case of the (elastic) Boltzmann case, Eq. (3). It is immediate to verify that it is true, indeed:
and being a linear combination of conserved quantities, the second term of the right hand side is constant, so that .
In Appendix A we report two derivations, necessary for the proof of result (14) below. Such derivations are nothing but the well known proofs of Eq. (5) and Eq. (7): the reader may verify that, by replacing with and with in Appendix A.2, one obtains the proof of Boltzmann H-theorem, and, replacing with in Appendix A.1, the H-theorem for Markov processes is proved.
Iii Granular gases with homogeneous energy injection
In the case of dilute granular gases interactions are dissipative (e.g. inelastic hard-core collisions) and the property in Eq. (2) does not hold. One of the consequences of energy dissipation is that the Boltzmann equation van Noije and Ernst (1998) without any external energy input has usually a trivial asymptotic state (e.g. the velocities of all particles vanish). Among the many models of energy injection Williams and MacKintosh (1996); van Noije et al. (1999); Gallas et al. (1992), experiments Gradenigo et al. (2011); Puglisi et al. (2012) have shown the relevance of a mechanism where all particles are coupled with a random energy reservoir Puglisi et al. (1998): such a simple mechanism well reproduces the effect of an interaction of all the particles with rough vibrating boundaries of the container.
In this model, under the hypothesis of molecular chaos and in the discrete representation introduced above, the probability obeys the following “Boltzmann-Fokker-Planck” (BFP) equation:
As anticipated, when collisions are inelastic the time inversion symmetry, Eq. (2), does not hold, i.e. . Notice that, in general, the existence of a solution satisfying detailed balance with respect to the Markov rates is not required to guarantee a stationary state of Eq. (9). However in the literature detailed balance with respect to has been often assumed Puglisi et al. (1998) and, for this reason, it is customary to consider such an energy injection mechanism equivalent to coupling the gas with a heat bath.
For continuous variables Eq. (9) is replaced by an evolution equation for the density :
where is a linear Fokker-Planck operator and is a bilinear integral operator (the inelastic Boltzmann-collision integral van Noije and Ernst (1998)). The first operator describes the interaction of the system with the heat-bath necessary to render the system stationary, while the second describes the collisions between the particles: the combination of the two operators produce a non trivial velocity distribution. Up to our knowledge, for Eq. (9) (or its counterpart with continuous velocities, Eq. (10)), no kind of “H-theorem” is known. The failure of the usual “H-theorem”, that for the functional, has been verified in Bena et al. (2006).
Notice that, assuming a relaxation toward equilibrium, i.e. , then it is easy to show that at large times is non-increasing. In fact, by writing with small, one has
We also mention that, for particular granular models, it is possible to prove that Eq. (10) in the elastic limit becomes equivalent to a Fokker-Planck equation: in that case the H-theorem for the functional is obviously verified Pareschi and Toscani (2006).
iii.1 An elastic granular gas
Let us analyze a particular limit of the previous model where the collisions are elastic , so that . The governing equation reads
and in addition we assume that the invariant probability is a stationary solution of both the linear master equation and the non-linear Boltzmann equation: i.e. for the satisfying Eq. (6), one has
if the collision is allowed.
It is important to realize that, even if the invariant probability is somehow trivial, the dynamics is not, since it depends upon the interplay between the bath and the collisions which may happen, for instance, on different timescales. The example discussed in Sec. IV.1, see Fig. 2, well illustrates this point, by showing the non-trivial dynamics of which is non-monotonic.
In such a system one shows that
In fact, we can write
Since we can rewrite as
It is now easy to show that both and are negative: it is enough to follow the standard proofs of the H-theorems for the Master equation and for the Boltzmann equation separately: those proofs are reported in the Appendix A for completeness.
iii.2 Granular Maxwell model
We discuss now an inelastic Maxwell model, a variation upon a theme, originally proposed by Ulam Ulam (1980); Ernst (1981) to study the approach to equilibrium, introduced by Ben Naim and Kaprivski, as a minimal kinetic model for granular gases Ben-Naim and Krapivsky (2000); Bobylev et al. (2000); Baldassarri et al. (2002) . The advantage of this model is that all moments can be explicitly computed Marini Bettolo Marconi and Puglisi (2002). In the 1d thermostatted version of the model, Eq. (10), takes the form
where is the restitution coefficient (when the collisions are elastic), the first term on the right hand describes the effect of the heat-bath at temperature and corresponds to , while the last term represents the non linear collisional term, which is the sum of a gain term and a loss term.
In general it is not possible to write explicitly , however, it is possible to obtain the evolution law of its moments defined as and prove that the high velocity tails of the PdF are gaussian. The evaluation of the -function requires the PdF, so that we used an approximation which correctly reproduces all moments up to a given order and displays the correct high velocity tails. To achieve that, we consider the following Sonine-Hermite representation of :
where is the non dimensional velocity squared and is the scaled PdF related to by the transformation
For the sake of simplicity we assumed that the distribution is an even function of the velocity so that all its odd moments vanish. The are the Sonine polynomials of order , given by the formula Liboff (2003):
having the property
and the are coefficients of the expansion related to the moments as shown in the appendix. They can be calculated from the averages
Since the evaluation of the moments of a given order requires only the knowledge of the moments of lower order one can proceed without excessive difficulty to any desired order. In practice, we carried on our calculation up to the eighth moment.
We can now evaluate the -function for the Maxwell model using the Sonine representation of by performing numerically the following integral:
In Figs. 1 we display the behavior of the and functions together with the evolution of the second moment of the PdF for two different initial conditions but having the same steady state. In the first case the second moment decreases towards its asymptotic value, while the function also decreases, whereas the function increases. In the second case instead the increases and both and decrease in time. The two examples show that it is necessary to consider always , while is not always monotonically decreasing.
A fast and simple way to prove numerically the hypothesis that the function is always decreasing is to consider its evolution for a short time interval starting from a distribution of the form (21) at the instant generated assuming initial values of the moments arbitrarily with the only constraint that this PdF is everywhere non negative. We then compute the evolution of the PdF over a small time interval, , using the governing equations for the moments and finally calculate the variation . For all possible choices of the initial values we have found that such a variation is which turns out to be always non positive. It is worth to comment that it is not necessary to follow the system evolution over a longer time interval, to verify the persistence of the sign of . In fact, any possible distribution which can be reached at time by the dynamics starting from the distribution , represents itself a good candidate as initial distribution, whose choice is arbitrary. We have sampled a large number of initial conditions of the form (21) randomly generated and verified that .
We conclude this section by saying that although we cannot prove analytically that for the inelastic Maxwell there exist an H-theorem, our numerical results provide a strong evidence that this is the case. The results of this section should be compared with those obtained for other granular systems, as discussed in the next Section.
Iv Numerical evidence for other examples of granular systems
iv.1 Granular gases with discrete states
We present here the simulation of a discrete BFP equation (9) for a choice of (conservative or non-conservative) collisions and with the presence of a thermal bath. In particular we have assumed that each state represents a possible value of the single particle energy , for instance where is some amount of energy. Collisions have been assumed to mix energy between colliding particles and, optionally, to dissipate a part of it, in order to reproduce the inelasticity in granular gases.
The collision model we adopted assigns the following collision rule to transform colliding energies into post-collision energies :
with being the amount of dissipated energy and the additional condition is enforced. Just for simplicity we assume that the probability of two particles of being chosen for a collision is independent of the relative velocity, as it occurs in so-called Maxwell models previously discussed Ernst (1981); Bobylev et al. (2000); Baldassarri et al. (2002). The single particle mean free time between collision is defined as . The stated collision model determines the rates in the equation (9). To avoid cumbersome expressions we do not reproduce here such rates.
The action of the heat bath is taken into account through the linear part of Eq. (9). In particular we have assumed that
where we have introduced , the characteristic bath time, and , the bath temperature. Rates satisfy detailed balance with respect to the “equilibrium” stationary distribution with a normalization constant. Such equilibrium stationary distribution is attained when collisions are switched off (i.e. ) as well as when they are elastic (i.e. ).
We have simulated Eq. (9) with the rates discussed above by means of a fourth order Runge-Kutta integration.
In Fig. 2 we show that when the collision are elastic () but a thermal bath acts on each particle, one has a monotonic decrease of but not of . Note that both mechanisms (elastic collisions and thermal bath), separately, guarantee the existence of an equilibrium steady state: however the thermal bath enforces a well defined temperature , while elastic collisions do not (the temperature, in the absence of the bath, would be chosen by initial conditions, i.e. initial average energy). When the two mechanisms act together, the final temperature is that of the bath, . Panel a) of Fig. 2 shows the evolution of (from initial conditions concentrated uniformly between and ) toward such an asymptotic equilibrium distribution. Correspondingly, panel b) shows the evolution of the average energy of the system. Panel c) shows the non-monotonic behavior of the usual Boltzmann function. Finally panel d) shows the fact, anticipated in Section III.1, that the evolution of is non-increasing.
In Fig. 3 we report the results for the case with dissipation (), in the presence of a thermal bath in order to guarantee the attainment of a steady state. In this case the knowledge of is not known a priori, therefore a first long run of the simulation is used to obtain it. A second run is then used to measure . We have repeated the simulation for three different initial conditions, as detailed in the figure caption. The evolution of the probability distribution for the particular case starting with concentrated between and is shown in panel a) of the Figure: it reaches an asymptotic distribution (the same for all initial conditions) different from the equilibrium one. In panel b) we observe the average energy which settles to a value smaller than the bath temperature because of the inelasticity of collisions. Panel c) and panel d) show that, while not always verifies the H-theorem, the validity of the H-theorem for is always verified.
iv.2 Granular gases with continuous states through the Direct Simulation Monte Carlo
Direct Simulation Monte Carlo (DSMC) Bird (1994); Cercignani et al. (1994) is usually considered an effective “solver” for Boltzmann equations and has been frequently used in the study of the kinetics of granular gases. It is a so-called “particles method”, since a finite number of particles is evolved stochastically: the statistics of those particles approximates, as the solution of the corresponding Boltzmann equation. For the purpose of the present paper, therefore, the study of the evolution of during the DSMC dynamics, with non-conservative interactions and with the presence of a heat bath, is meaningful as becomes larger and larger.
Here we use the DSMC algorithm discarding any spatial information, however the algorithm is often used with spatial coordinates, by dividing space in small cells. Our choice is equivalent to consider the space-homogeneous version of the BFP equation (10). In the algorithm time is advanced in time steps of length . At each time step two sub-steps are performed: 1) the heat bath step, where the velocity of each particle is advanced, from to , by the discretized solution of an Ornstein-Uhlenbeck stochastic process , i.e.
where is a random variable extracted from a normal distribution (different and different are all independent), is the bath temperature and is the typical interaction time of a particle with the bath; 2) collisions are performed by choosing random couples of particles and changing their velocities by the rule
and the rate of collision per particle is fixed at . The parameter represents the restitution coefficient, when collisions are elastic, otherwise they dissipate part of the kinetic energy. Note that in this version of the DSMC the random choice of particles is done uniformly: this is equivalent to solve a Boltzmann equation with a collision probability independent from the relative velocity of the colliding particle, as it happens in the Maxwell models discussed in Section III.2.
The results of the DSMC dynamics for a case with inelastic with heat bath are shown in Fig. 4. In panel a) we have reported the evolution of the probability distribution, which is started from a Gaussian distribution and . In panel b) the evolution of the average energy is shown, demonstrating that the initial temperature is forgot and the system attains a stationary state with an average kinetic energy (also called granular temperature) , because of inelastic collisions. In frame c) and d) we have reported the evolution of
where and are the empirical velocity probability and its long time limit respectively: the empirical probability at instant is obtained by choosing a velocity interval with and dividing it in sub-intervals, and taking to be the number of particles with velocity in the -th sub-interval. For and , we have where , so that
The results reported in frame C and D of Fig. 4 demonstrate the failure of the usual Boltzmann H-theorem (the one for ) together with the validity of .
It is interesting to notice that in the examples shown in Figures 2, 3 and 4, the functional has a behavior constituted by a first decrease followed by an increase. This seems a consequence of the relatively fast action of collisions superimposed to a slower action of the heat bath: collisions (even if inelastic) do not change dramatically the energy and therefore mainly contribute to “equilibrate” the initial distribution, so that satisfies the original H-theorem () at the beginning; when the heat bath action becomes dominant, the distribution is near to the equilibrium one, i.e. and the effect of the bath (a Ornstein-Uhlenbeck process which, being linear, conserves the Gaussian shape), is mainly a decrease of energy (initiated higher than ), implying .
iv.3 The connection between the space and the space
It is interesting to notice that the DSMC algorithm is a Markov process and therefore, as discussed in Section II.2, we already know that an H-theorem holds. Nevertheless the DSMC is a Markov process for an -dimensional vector, i.e. in the so-called “ space”, while the BFP Eq. (10) governs the evolution of i.e. the single particle velocity distribution, i.e. it lives in the “ space”.
By defining the probability density at time in the -dimensional space, we know that for the DSMC it obeys a Master Equation of the kind
and it is customary to assume that it reaches a stationary state . Therefore the function
is a non-increasing function of time, i.e. .
In the case of the DSMC discussed above, which includes the interaction with an external energy injection mechanism and pairwise collisions, one may separate
where is the operator representing the single particle Master equation, in our case that for the Ornstein-Uhlenbeck process, while
is the operator for the collision between the particles, where operates on a function of by mapping into which are the pre-collisional velocities, obtained by inverting Eq. (30). For non space-homogeneous systems and hard-core interactions, the collisional operator is different (for instance it enforces the condition of contact among particles, as well as the velocity dependence of the scattering probability, etc.), see for instance Brey et al. (2004).
It is interesting to notice that, by assuming a stronger version of Molecular Chaos, valid for all velocities (i.e. not only the pre-collisional ones)
it is immediate to get
We stress however that such a stronger form of Molecular Chaos is usually violated in non spatially homogeneous granular gases, as shown in recent experiments Gradenigo et al. (2011); Puglisi et al. (2012). The fact that it could work in the models discussed above is perhaps a consequence of spatial homogeneity together with the action of the external heat bath.
In the present paper we have offered several examples of Boltzmann-Fokker-Planck (BFP) models with conservative and non-conservative interactions, where an H-functional of the kind
appears to be non-increasing for the whole evolution from arbitrary initial conditions toward the asymptotic steady state . The only case where we are able to prove such a conjecture is the elastic BFP case, where the proof is a “superposition” of the proofs of the two different H-theorems for the elastic Boltzmann equation and for Markov processes: notwithstanding the simplicity of the proof and the triviality of the steady state, the elastic BFP model has a non-trivial dynamics where collisions and thermostat may act on different timescales, as demonstrated by the non-monotonous behavior of the Boltzmann functional. Establishing the monotonicity of or other entropic functionals for a rather general class of models is certainly a challenge for future research Arnold et al. (2004); Villani (2006). In a remark concluding the last section, we have recalled that the BFP equation may be obtained by marginalizing a master equation for the Markovian evolution of the many-particles vector in space: this remark, which is behind the operating principle of DSMC schemes, could be a possible starting point for further investigations on this complex and fascinating subject.
Acknowledgements.The authors acknowledge useful discussions with A. Baldassarri. They also acknowledge the support of the Italian MIUR under the PRIN 2009 grant n. 2009PYYZM5. A. P. acknowledges the support of the Italian MIUR under the grant FIRB-IDEAS n. RBID08Z9JE.
Appendix A H-theorem in the case of conservative interactions and under the action of a heat bath
a.1 Master equation contribution
Let us introduce the variables and . It is also useful to define the function , so that . It is easy to verify that
where we have exchanged the indexes in the second part of the sum, in order to collect . Then one notices that for a set of numbers the following relation holds
Now, since , it is immediate to see that
and therefore . The proof shown here is that found in van Kampen (2007).
a.2 Boltzmann equation contribution
Here it is useful to introduce the variable . The symmetry between collisions and implies the following identities , so that
and by means of the symmetry between and one gets
We notice that . Furthermore, since for every and , we finally get . The equality holds only when the are identically , i.e. . The above proof is the standard one contained in any textbook and is due to Ludwig Boltzmann Boltzmann (1872).
Appendix B Time evolution of the moments for the Maxwell model
We provide here some details regarding the calculations presented in the main text. Applying the Fourier transform to eq. (19), we obtain the following governing equation for characteristic function :
with . Substituting the representation in (47) and equating the coefficients of the same power of we derive a set of ordinary differential equations describing the evolution of the moments, whose integration is straightforward and leads to expressions of the type:
where the coefficients are
where the coefficients have the following expressions:
Notice that vanish in the elastic limit as the elastic collisions do not affect the PdF nor its moments. Since in eq. (48) the evolution of the moment of order n is coupled only to the evolution of moments of order smaller than n, the solution is very simple and can be achieved recursively.
In order to construct the PdF we have to use the Sonine-Hermite expansion and compute the coefficients which are proportional to the cumulants of the distribution. In terms of the moments we find:
Notice that being proportional to the cumulants vanish for the Gaussian distribution, but not for the granular gas.
- L. Boltzmann, Sitzungsber. Akad. Wiss. Wien Math.-Naturwiss. 66, 275 (1872).
- E. G. D. Cohen (Editor), Fundamental Problems in Statistical Mechanics V (Elsevier, 1981).
- C. Cercignani, The man who trusted atoms (Oxford University Press, 2006).
- C. Cercignani, R. Illner, and M. Pulvirenti, The Mathematical Theory of Dilute Gases (Springer-Verlag, New York, 1994).
- P. Castiglione, M. Falcioni, A. Lesne, and A. Vulpiani, Chaos and Coarse Graining in Statistical Mechanics (Cambridge University Press, 2008).
- N. G. van Kampen, Can. J. Phys. 39, 551 (1961).
- A. Puglisi, V. Loreto, U. M.