A recursive algorithm and a series expansion related to the homogeneous Boltzmann equation for hard potentials with angular cutoff

We consider the spatially homogeneous Boltzmann equation for hard potentials with angular cutoff. This equation has a unique conservative weak solution $(f_t)_{t\geq 0}$, once the initial condition $f_0$ with finite mass and energy is fixed. Taking advantage of the energy conservation, we propose a recursive algorithm that produces a $(0,\infty)\times\mathbb{R}^3$ random variable $(M_t,V_t)$ such that $E[M_t {\bf 1}_{\{V_t \in \cdot\}}]=f_t$. We also write down a series expansion of $f_t$. Although both the algorithm and the series expansion might be theoretically interesting in that they explicitly express $f_t$ in terms of $f_0$, we believe that the algorithm is not very efficient in practice and that the series expansion is rather intractable. This is a tedious extension to non-Maxwellian molecules of Wild's sum and of its interpretation by McKean.


Introduction
We consider a spatially homogeneous gas modeled by the Boltzmann equation: the density f t (v) ≥ 0 of particles with velocity v ∈ R 3 at time t ≥ 0 solves where The cross section B is a nonnegative function given by physics. We refer to Cercignani [4] and Villani [16] for very complete books on the subject. We are concerned here with hard potentials with angular cutoff: the cross section satisfies and some bounded measurable b : The important case where γ = 1 and b is constant corresponds to a gas of hard spheres. If γ = 0, the cross section is velocity independent and one talks about Maxwellian molecules with cutoff.
We classically assume without loss of generality that the initial mass R 3 f 0 (v)dv = 1 and we denote by e 0 = R 3 |v| 2 f 0 (v)dv > 0 the initial kinetic energy. It is then well-known, see Mischler-Wennberg [12], that (1) has a unique weak solution such that for all t ≥ 0, f t is a probability density on R 3 with energy R 3 |v| 2 f t (v)dv = e 0 . Some precise statements are recalled in the next section.
In the rest of this introduction, we informally recall how (1) can be solved, in the case of Maxwellian molecules, by using the Wild sum, we quickly explain its interpretation by McKean, and we write down a closely related recursive simulation algorithm. We also recall that Wild's sum can be used for theoretical and numerical analysis of Maxwellian molecules. Then we briefly recall how the Wild sum and the algorithm can be easily extended to the case of any bounded cross section, by introducing fictitious jumps. Finally, we quickly explain our strategy to deal with hard potentials with angular cutoff.
1.1. Wild's sum. Let us first mention that some introductions to Wild's sum and its probabilistic interpretation by McKean can be found in the book of Villani [16,Section 4.1] and in Carlen-Carvalho-Gabetta [1,2]. Wild [18] noted that for Maxwellian molecules, i.e. when γ = 0, so that the cross section B(v − v * , σ) = β v,v * (σ) does not depend on the relative velocity, (1) rewrites where, for f, g two probability densities on R 3 , dσ. It holds that Q(f, g) is also a probability density on R 3 , that can be interpreted as the law of V ′ = (V + V * + |V − V * |σ)/2, where V and V * are two independent R 3 -valued random variables with densities f and g and where σ is, conditionally on (V, V * ), a κ −1 β V,V * (σ)dσ-distributed S 2 -valued random variable. Wild [18] proved that given f 0 , the solution f t to (1) is given by (5) f t = e −κt
McKean [10,11] provided an interpretation of the Wild sum in terms of binary trees, see also Villani [16] and Carlen-Carvalho-Gabetta [1]. Let T be the set of all discrete finite rooted ordered binary trees. By ordered, we mean that each node of Υ ∈ T with two children has a left child and a right child. We denote by ℓ(Υ) the number of leaves of Υ ∈ T . If • is the trivial tree (the one with only one node: the root), we set is the subtree of Υ consisting of the left (resp. right) child of the root with its whole progeny. Then (5) can be rewritten as In words, (6) can be interpreted as follows. For each Υ ∈ T , the term e −κt (1 − e −κt ) |Υ|−1 is the probability that a typical particle has Υ as (ordered) collision tree, while Q Υ (f 0 ) is the density of its velocity knowing that it has Υ as (ordered) collision tree.
Finally, let us mention a natural algorithmic interpretation of (1) closely related to (6). The dynamical probabilistic interpretation of Maxwellian molecules, initiated by Tanaka [15], can be roughly summarized as follows. Consider a typical particle in the gas. Initially, its velocity V 0 is f 0 -distributed. Then, at rate κ, that is, after an Exp(κ)-distributed random time τ , it collides with another particle: its velocity V 0 is replaced by is the velocity of an independent particle undergoing the same process (stopped at time τ ) and σ is a κ −1 β V0,V * τ (σ)dσ-distributed S 2 -valued random variable. Then, at rate κ, it collides again, etc. This produces a stochastic process Consider now the following recursive algorithm.
Function velocity(t): .. Simulate a f 0 -distributed random variable v, set s = 0. .. While s < t do .. .. simulate an exponential random variable τ with parameter κ, Of course, each new random variable is simulated independently of the previous ones. In particular, line 7 of the algorithm, all the random variables required to produce v * =velocity(s) are independent of all that has already been simulated.
Comparing the above paragraph and the algorithm, it appears clearly that velocity(t) produces a f t -distributed random variable. We have never seen this fact written precisely as it is here, but it is more or less well-known folklore. In the present paper, we will prove such a fact, in a slightly more complicated situation.
In spirit, the algorithm produces a binary ordered tree: each time the recursive function calls itself, we add a branch (on the right). So it is closely related to (6) and, actually, one can get convinced that velocity(t) is precisely an algorithmic interpretation of (6). But entering into the details would lead us to tedious and technical explanations.
1.2. Utility of Wild's sum. The Wild sum has often been used for numerical computations: one simply cutoffs (5) at some well-chosen level and, possibly, adds a Gaussian distribution with adequate mean and covariance matrix to make it have the desired mass and energy. See Carlen-Salvarini [3] for a very precise study in this direction. And actually, Pareschi-Russo [13] also managed to use the Wild sum, among many other things, to solve numerically the inhomogeneous Boltzmann equation for non Maxwellian molecules.
A completely different approach is to use a large number N of times the perfect simulation algorithm previously described to produce some i.i.d. f t -distributed random variables V 1 t , . . . , V N t , and to approximate f t by N −1 N 1 δ V i t . We believe that this is not very efficient in practice, especially when compared to the use of a classical interacting particle system in the spirit of Kac [8], see e.g. [7]. The main reason is that the computational cost of the above perfect simulation algorithm increases exponentially with time, while the one of Kac's particle system increases linearly. So the cost to remove the bias is disproportionate. See [6] for such a study concerning the Smoluchowski equation, which has the same structure (at the rough level) as the Boltzmann equation.
The Wild sum has also been intensively used to study the rate of approach to equilibrium of Maxwellian molecules. This was initiated by McKean [10], with more recent studies by Carlen-Carvalho-Gabetta [1,2], themselves followed by Dolera-Gabetta-Regazzini [5] and many other authors.
with Φ bounded, e.g. by 1, we can introduce fictitious jumps to write (1) as } and something similar for v ′′ * . Hence all the previous study directly applies, but the resulting Wild sum does not seem to allow for a precise study the large time behavior of f t , because it leads to intractable computations.

1.4.
Hard potentials with angular cutoff. Of course, the angular cutoff (that is, we assume that κ < ∞) is crucial to hope for a perfect simulation algorithm and for a series expansion in the spirit of Wild's sum. Indeed, κ = ∞ implies that a particle is subjected to infinitely many collisions on each finite time interval. So our goal is to extend, at the price of many complications, the algorithm and series expansion to hard potentials with cutoff. Since the cross section is unbounded in the relative velocity variable, some work is needed.
We work with weak forms of PDEs for simplicity. First, it is classical, see e.g. [16,Section 2 As already mentioned, one also has R 3 |v| 2 f t (dv) = e 0 , so that g t (dv) = (1 + e 0 ) −1 (1 + |v| 2 )f t (dv) belongs to P(R 3 ) for all t ≥ 0. A simple computation shows that, for all reasonable φ ∈ C b (R 3 ), This equation enjoys the pleasant property that the maximum rate of collision, given by Λ 0 (v) = κ sup v * ∈R 3 (1 + |v * | 2 ) −1 (1 + e 0 )|v − v * | γ , is finite. Hence, up to some fictitious jumps, one is able to predict when a particle will collide from its sole velocity, knowing nothing of the environment represented by g t . The function Λ 0 is not bounded as a function of v, since it resembles 1+|v| γ , but, as we will see, this it does not matter too much. On the contrary, the presence of (1 + |v| 2 ) −1 (1 + |v ′ | 2 ) in front of the gain term is problematic. It means that, in some sense, particles are not all taken into account equally. To overcome this problem, we consider an equation with an additional weight variable m ∈ (0, ∞). So we search for an equation, resembling more the Kolmogorov forward equation of a nonlinear Markov process, of which the solution (G t ) t≥0 ⊂ P((0, ∞) × R 3 ) would be such that ∞ 0 mG t (dm, dv) = g t (dv) for all times. Then, one would recover the solution to (1) as f t (dv) = (1 + e 0 )(1 + |v| 2 ) −1 ∞ 0 mG t (dm, dv). All this is doable and was our initial strategy. However, we then found a more direct way to proceed: taking advantage of the energy conservation, it is possible to build an equation of which the solution Here Λ is any (explicit) function dominating Λ 0 , the additional variable a is here to allow for fictitious jumps, and the post-collisional characteristics (m ′′ , v ′′ ), depending on m, v, m * , v * , σ, a are precisely defined in the next section. The perfect simulation algorithm for such an equation is almost as simple as the one previously described, except that the rate of collision Λ(v) now depends on the state of the particle. On the contrary, this state-dependent rate complicates subsequently the series expansion because the time and phase variables do not separate anymore.
1.5. Plan of the paper. In the next section, we expose our main results: we introduce an equation with an additional variable m, state that this equation has a unique solution F t that, once integrated in m, produces the solution f t to (1). We then propose an algorithm that perfectly simulates an F t -distributed random variable, we write down a series expansion for F t in the spirit of (6) and discuss briefly the relevance of our results. The proofs are then handled: the algorithm is studied in Section 3, the series expansion established in Section 4, well-posedness of the equation solved by F t is checked in Section 5, and the link between F t and f t shown in Section 6.

Weak solutions.
We use a classical definition of weak solutions, see e.g. [16, Section 2.3].
Definition 1. Assume (3) and recall (4). where Everything is well-defined in (7) by boundedness of mass and energy and since For any given f 0 ∈ P(R 3 ) such that R 3 |v| 2 f 0 (dv) < ∞, the existence of unique weak solution starting from f 0 is now well-known. See Mischler-Wennberg [12] when f 0 has a density and Lu-Mouhot [9] for the general case. Let us also mention that the conservation assumption is important in Definition 1, since Wennberg [17] proved that there also solutions with increasing energy.

An equation with an additional variable.
We fix e 0 > 0 and define, for v ∈ R 3 , We also introduce E = (0, ∞) × R 3 and, for y = (m, v) and y * = (m * , v * ) in E and z ∈ H, with a small abuse of notation. We finally consider, for y = (m, v) and y * = (m * , v * ) in E, which is a measure on H with total mass ν y,y * (H) = κ, see (3).
All is well-defined in (8) thanks to the conditions on F and since |BΦ(y, y * )| ≤ 2κ||Φ|| ∞ Λ(y) ≤ C Φ (1 + |v| γ ) (with the notation y = (m, v)). As already mentioned in the introduction, the important point is that the function Λ does not depend on y * . Hence a particle, when in the state y, jumps at rate κΛ(y), independently of everything else.
We will also verify the following estimate.
Finally, the link with the Boltzmann equation is as follows.
and if the quantity e 0 used to define the coefficients of (A) is precisely e 0 = R 3 |v| 2 f 0 (dv), then (f t ) t≥0 is the unique weak solution to (1) starting from f 0 .

2.3.
A perfect simulation algorithm. We consider the following procedure. Algorithm 6. Fix e 0 > 0 and F 0 ∈ P(E). For any t ≥ 0 we define the following recursive function, of which the result is some E × N-valued random variable.
Of course, each time a new random variable is simulated, we implicitly assume that it is independent of everything that has already been simulated. In particular, line 7 of the procedure, all the random variables used to produce (y * , n * ) =(value(s),counter(s)) are independent of all the random variables already simulated. By construction, counter(t) is precisely the number of times the recursive function calls itself.
Algorithm 6 a.s. stops and thus produces a couple A series expansion. We next write down a series expansion of F t , the solution to (A), in the spirit of Wild's sum (6). Unfortunately, the expressions are more complicated, because the time (t ≥ 0) and phase (y ∈ E) variables do not separate. This is due to the fact that the jump rate Λ depends on the state of the particle.
Observe that (8) may be written, at least formally, We finally consider the set T of all finite binary (discrete) ordered trees: such a tree is constituted of a finite number of nodes, including the root, each of these nodes having either 0 or two children (ordered, in the sense that a node having two children has a left child and a right child). We denote by • ∈ T the trivial tree, composed of the root as only node. Proposition 8. Assume (3). Let F 0 ∈ P(E) such that E |v| γ F 0 (dm, dv) < ∞. The unique solution (F t ) t≥0 to (A) starting from F 0 is given by where Υ ℓ (resp. Υ r ) is the subtree of Υ consisting of the left (resp. right) child of the root with its whole progeny.
We will prove this formula by a purely analytic method. We do not want to discuss precisely its connection with Algorithm 6, but let us mention that in spirit, the algorithm produces a (random) ordered tree Υ t of interactions together with the value of Y t , and that Γ t (J Υ (F 0 )) can be interpreted as the probability distribution of Y t restricted to the event that Υ t = Υ.
(a) Gathering Propositions 7 and 5, we find that Algorithm 6 used with e 0 = R 3 |v| 2 f 0 (dv) and with F 0 produces a random variable

Indeed, we know from Proposition 7 that
2.6. Discussion. It might be possible to prove Proposition 5 assuming only that F 0 ∈ P(E) (1) is known to be well-posed as soon as the initial energy is finite, see [12,9]. However, it would clearly be more difficult and our condition is rather harmless.
Observe that (A) is well-posed under the condition that F 0 ∈ P(E) satisfies E |v| γ F 0 (dm, dv), which does not at all imply that e 0 = E m|v| 2 F 0 (dm, dv) < ∞. But, recalling that the e 0 has to be finite for the coefficients of (A) to be well-defined, this is not very interesting.
The series expansion of Proposition 8 is of course much more complicated than the original Wild sum, since (a) we had to add the variable m, (b) we had to introduce fictitious jumps, (c) time and space do not separate. So it is not clear whether the formula can be used theoretically or numerically. However, it provides an explicit formula expressing f t as a (tedious) function of f 0 . Algorithm 6 is extremely simple. Using it a large number of times, which produces some i.i.d.
For a central limit theorem to hold true, one needs E[M 2 t ] = E m 2 F t (dm, dv) to be finite. We do not know if this holds true, although we have some serious doubts. Hence the convergence of this Monte-Carlo approximation may be much slower that N −1/2 . The main interest of Algorithm 6 is thus theoretical.
Step 1. We now consider the following procedure. It is an abstract procedure, because it assumes that for each t ≥ 0, one can simulate a random variable with law G t and because the instructions are repeated ad infinitum if the cemetery point is not attained.
Observe that in the last line, we may have s < ∞ either because after a finite number of steps, the simulation of (y * , n * ) with law G s has produced (△, ∞), or because we did repeat the loop ad infinitum, but the increasing process N became infinite in finite time.
At the end, this produces a process (Y t , N t ) t≥0 and one easily gets convinced that for each t ≥ 0, (Y t , N t ) is G t -distributed. Indeed, if one extracts from the above procedure only what is required to produce (Y t , N t ) (for some fixed t), one precisely re-obtains Algorithm 6 if (Y t , N t ) = (△, ∞) (and in this case Algorithm 6 stops), while (Y t , N t ) = (△, ∞) implies that Algorithm 6 never stops.
, all y ∈ E and all n ∈ N.
Step 2. Here we handle a preliminary computation: for all y, y * ∈ E, we have Writing y = (m, v), y * = (m * , v * ) and recalling Subsection 2.2, A(y, y * ) equals Step 3. We now prove that (Y t , N t ) t≥0 actually does not explode nor reach the cemetery point, For A ∈ N * , we introduce ζ A = inf{t ≥ 0 : N t ≥ A}. The process (Y t , N t ) t≥0 does not explode nor reach the cemetery point during [0, ζ A ), so that we can write, with Ψ(y, n) = n ∧ A, (recall that N 0 = 0 and that t → N t is a.s. non-decreasing), whence, by the Gronwall lemma, We next choose Ψ(y, n) = Λ(y)1I {n<A} and write, as previously, But L s Ψ(y, n) equals Gathering (10) and (11) and letting A increase to infinity, we first conclude that E[N t ] < ∞ for all t ≥ 0. In particular, N t < ∞ a.s. for all t ≥ 0, and the process (Y t , N t ) t≥0 does a.s. not explode and never reach (△, ∞). (11). Finally, we easily conclude from (10) Step 4. By Step 3, we know that G t (which is the law of (Y t , N t )) is actually supported by E × N for all t ≥ 0. Hence Algorithm 6 a.s. stops. The process (Y t , N t ) t≥0 is thus an inhomogeneous Markov with generatorL t defined, for Ψ ∈ C b (E × N), bỹ L t Ψ(y, n) = Λ(y) H E×N Ψ(h(y, y * , z), n + n * + 1) − Ψ(y, n) G t (dy * , dn * )ν y,y * (dz) and we thus have Let now F t ∈ P(E) be the law of Y t (so F t is the first marginal of G t ). It starts from F 0 and solves (A).
is locally bounded by Step 3 and for all Φ ∈ C b (E), applying the above equation with Ψ(y, n) = Φ(y), we findL t Ψ(y, as desired. Finally, we have already seen in We have proved Proposition 7, as well as the existence part of Proposition 3.
Step 3. We fix k ∈ N * and denote by T k ⊂ T the finite set of all ordered binary trees with at most k nodes. We introduce F k t = Υ∈T k Γ t (J Υ ). By Step 2, we know that t → E Λ(y)F k t (dy) is locally bounded. We claim that for all Φ ∈ C b (E), all t ≥ 0, and then that, for Υ ∈ T k \ {•}, Since J Υ (ds, dy) = Q(Γ s (J Υ ℓ ), Γ s (J Υr ))(dy)ds by definition, the result follows.
Step 6. By Step 5, the series of nonnegative measures F t = Υ∈T Γ t (J Υ ) converges, satisfies F t (E) ≤ 1, and we know that t → E Λ(y)F t (dy) is locally bounded. Passing to the limit in the time-integrated version of (13), we find that for all Φ ∈ C b (E), all t ≥ 0, To justify the limiting procedure, it suffices to use that t → E Λ(y)F t (dy) is locally bounded, as well as . We used that the map Υ → (Υ ℓ , Υ r ) is bijective from T \ {•} into T × T , as well as the bilinearity of Q. By the same way, (16) rewrites as see (15). To conclude that (F t ) t≥0 solves (A), it only remains to verify that t → E |v| γ F t (dm, dv) is locally bounded, which follows from the fact that E |v| γ F t (dm, dv) ≤ E Λ(y)F t (dy), and that F t (E) = 1 for all t ≥ 0. Applying the previous equation with Φ = 1 (for which BΦ = 0), we find that F t (E) = 1 + t 0 (F s (E) − 1)α s ds, where α s = κ E Λ(y)F s (dy) is locally bounded. Hence F t (E) = 1 for all t ≥ 0 by the Gronwall lemma. The proof of Proposition 8 is complete.

Well-posedness of (A)
We have already checked (twice) the existence part of Proposition 3. We now turn to uniqueness. Let us consider two solutions F and G to (A) with F 0 = G 0 . By assumption, we know that where for µ a finite signed measure on E, |µ| = µ + + µ − with the usual definitions of µ + and µ − . We fix Φ ∈ C b (E) such that ||Φ|| ∞ ≤ 1 and we use Definition 2 to write where . Using only that |BΦ(y, y * )| ≤ 2κ||Φ|| ∞ Λ(y) ≤ 2κΛ(y), we get We next recall that ||Φ|| ∞ ≤ 1 and that E G t (dy * ) = 1 and we write is a nonnegative measure bounded by 2(F t + G t )(dy), we may write, for any M ≥ 1, Integrating in time (recall that u Φ 0 = 0) and taking the supremum over Φ ∈ C b (E) such that ||Φ|| ∞ ≤ 1, we find u t e κMt ≤ t 0 [(2κα s + κM )u s + 2κǫ M s ]e κMs ds. Recall the following generalized Gronwall lemma: if we have three locally bounded nonnegative Applying this result with v t = u t e κMt , g t = 2κǫ M t e κMt and h t = 2κα t + κM , we get Recalling that α is locally bounded and that t 0 ǫ M s ds tends to 0 as M → ∞, we conclude that u t = 0, which was our goal. The proof of Proposition 3 is now complete.

Relation between (A) and the Boltzmann equation
It remains to prove Proposition 5. In the whole section, we consider the solution F to (A) starting from some F 0 ∈ P(E) such that E [|v| γ + m(1 + |v| 2+2γ )]F 0 (dy) < ∞. We define the nonnegative measure f t on R 3 by f t (A) = E m1I {v∈A} F t (dy) for all A ∈ B(R 3 ) and we assume that f 0 ∈ P(R 3 ) and that R 3 |v| 2 f 0 (dv) = e 0 , where e 0 was used in Subsection 2.2 to build the coefficients of (A). We want to prove that f = (f t ) t≥0 is a weak solution to (1).
The main difficulty is to establish properly the following estimate, of which the proof is postponed at the end of the section.
Next, we handle a few preliminary computations.
It only remains to prove Lemma 9.