Particle approximation of some Landau equations

We consider a class of nonlinear partial-differential equations, including the spatially homogeneous Fokker-Planck-Landau equation for Maxwell (or pseudo-Maxwell) molecules. Continuing the work of Fontbona-Gu\'erin-M\'el\'eard, we propose a probabilistic interpretation of such a P.D.E. in terms of a nonlinear stochastic differential equation driven by a standard Brownian motion. We derive a numerical scheme, based on a system of $n$ particles driven by $n$ Brownian motions, and study its rate of convergence. We finally deal with the possible extension of our numerical scheme to the case of the Landau equation for soft potentials, and give some numerical results.


Introduction and main results
1.1. The equation. Let S d be the set of symmetric d × d matrices with real entries, and S + d its subset of nonnegative matrices. For a : R d → S + d , we consider the partial differential equation where ∂ t = ∂ ∂t , ∂ i = ∂ ∂xi and where the unknown (f t ) t≥0 is a family of probability density functions (f t ) t≥0 on R d . The spatially homogeneous Landau (or Fokker-Planck-Landau) equation corresponds, in dimension d ≥ 2, to the case where for some κ : R + → R + , (2) a ij (z) = κ(|z| 2 )(|z| 2 δ ij − z i z j ).
Physically, one assumes that κ(r) = r γ/2 , for some γ ∈ [−3, 1]. One talks of soft potentials when γ < 0, Maxwell molecules when γ = 0, and hard potentials when γ > 0. We consider in this paper the case of Maxwell molecules, or of pseudo-Maxwell molecules, where κ is supposed to be smooth and bounded. This equation arises as a limit of the Boltzmann equation when all the collisions become grazing. We refer to Villani [10,11,12] and the many references therein for physical and mathematical details on this topic. See Cordier-Mancini [2] and Buet-Cordier-Filbet [1] for a review on deterministic numerical methods to solve (1).

1.2.
Notation. Let P = P(R d ) be the set of probability measures on R d , and P k = {µ ∈ P, m k (µ) < ∞}, where m k (µ) = |x| k µ(dx). For x, y ∈ R d , we set |x| = ( d 1 x 2 i ) . Assume that |a(x)| + |b(x)| ≤ C(1 + |x| 2 ) (which is the case when a is defined by (2) with κ ∈ C 1 b ). A measurable family (P t ) t≥0 ⊂ P 2 is said to be a weak solution to (1) if for all t ≥ 0, sup [0,t] m 2 (P s ) < ∞ and for all ϕ ∈ C 2 b (R d ), where Lϕ(x, y) = 1 All the terms make sense due to our conditions on a, b, P t . See Villani [11] for a similar formulation.
1.3. Known results. To our knowledge, the first (and only) paper proving a rate of convergence for a numerical scheme to solve (1) is that of Fontbona-Guérin-Méléard [4]. Their method relies on a stochastic particle system. The aim of this paper is to go further in this direction. Let us thus recall briefly the method of [4], relying on the probabilistic interpretation of (1) developped by Funaki [6], Guérin [7]. Let σ : R d → S + d and b : R d → R d be Lipschitz continuous functions, and let P 0 Here independent of X 0 , with independent coordinates, each of which having covariance measure P t (dx)dt (see Walsh [14]). Existence and uniqueness in law for E 0 (P 0 , σ, b) have been proved in Guérin [7]. If furthermore σ(x)σ * (x) = a(x) and b i (x) = d j=1 ∂ j a ij (x), then (P t ) t≥0 is a weak solution to (1). The condition that σ and b are Lipschitz continuous is satisfied in the case of the Landau equation for Maxwell or pseudo-Maxwell molecules.
In [4], one considers an exchangeable stochastic particle system (X i,n t ) t≥0,i=1,...,n , satisfying a S.D.E. driven by n 2 Brownian motions. It is then shown that one may find a coupling between a solution (X 1 t ) t≥0 to E 0 (P 0 , σ, b) and such a particle system in such a way that under the condition that P 0 has a finite moment of order d+5. The proof relies on a clever coupling between the the white noise and n Brownian motions. In particular, one has to assume that P t has a density for all t > 0, in order to guarantee the uniqueness of some optimal couplings.
For each x ∈ R d , µ ∈ P 2 , a(x, µ) is a nonnegative symmetric matrix and thus admits an unique symmetric nonnegative square root a  Denote by W d the law of the d-dimensional Brownian motion, consider P 0 ∈ P 2 , and let (X 0 , B) ∼ P 0 ⊗ W d . We say that a R d -valued process (X t ) t≥0 solves E 1 (P 0 , a, b) (or E 1 (P 0 , a, b, X 0 , B) when This equation is nonlinear in the sense that its coefficients involve the law of the solution. Compared to (4), equation (5) is simpler, since it is driven by a finite-dimensional Brownian motion, and since the nonlinearity does not involve the driving process. However, one may check that at least formally, solutions to (4) and (5) have the same law. The link with (1) relies on a simple application of the Itô formula.
The natural linearization of (5) consists of considering n particles (X i,n t ) t≥0,i=1,...,n solving Here We thus use n Brownian motions. When linearizing (4), one needs to use n 2 Brownian motions, since the white noise is infinite dimensional. However, one may check that the solution to (6) and the particle system built in [4] have the same distribution (provided σσ * = a in [4, Equation (4)]).

1.5.
Main results. The main result of this paper is the following.
Theorem 3. Assume that b is Lipschitz continuous, that a is of class C 2 , with all its derivatives of order 2 bounded, and that P 0 ∈ P 2 . (i) There is strong existence and uniqueness for E 1 (P 0 , a, b): for any (X 0 , B) ∼ P 0 ⊗ W d , there is an unique solution (X t ) t≥0 to E 1 (P 0 , a, b, X 0 , B).
There is an unique solution (X i,n t ) t≥0,i=1,...,n to (6). Assume that P 0 ∈ P 4 , and consider the unique solution (X 1 t ) t≥0 to E 1 (P 0 , a, b, X 1 0 , B 1 ). There is a constant C T depending only on d, P 0 , a, b, T such that In the general case, we thus prove a rate of convergence in n −1/2 , which is faster than n −2/(d+4) . If we have some information on the nondegeneracy of a(x, P t ), then a 1 2 (x, µ) is smooth around µ ≃ P s , and we can get a better rate of convergence. Assume for example that a is uniformly elliptic (which is unfortunately not the case of (2), since a(x)x = 0 for all x ∈ R d ). Then sup x |a(x, P t ) −1 | ≤ sup y |a(y) −1 | < ∞, and we get a convergence rate in n −1 . In the case of the Landau equation for true Maxwell molecules, we obtain the following result.
Then a, b satisfy the assumptions of Theorem 3. Let P 0 ∈ P 4 , and adopt the notation of .
We finally consider the case of pseudo-Maxwell molecules.

Corollary 5. Consider the Landau equation for pseudo-Maxwell molecules, where a is given by
Assume that κ ′ has a bounded support. Then a, b satisfy the assumptions of Theorem 3-(ii). Assume furthermore that P 0 ∈ P 4 has a density with a finite entropy P 0 (x) log P 0 (x)dx < ∞, and that κ is bounded below by a positive constant. With the notation of Theorem 3, we have 1.6. Time discretization. To get a simulable particle system, it remains to discretize time in (6). Let N ≥ 1, and consider ρ N (s) = k≥0 k N 1 s∈[k/N,(k+1)/N ) . Consider the simulable particle system (X i,n,N t ) t≥0,i=1,...,n defined by ds.
Theorem 6. Assume that b is Lipschitz continuous, that a is of class C 2 , with all its derivatives of order 2 bounded, and that (6) and (X i,n,N t ) t≥0,i=1,...,n to (8). Then there is a constant C T depending only on d, P 0 , a, b, T such that 1.7. Conclusion. Choosing for example a, b, and P 0 as in Corollary 4-(ii) or as in Corollary 5, denoting by (P t ) t≥0 = (L(X 1 t )) t≥0 the weak solution to the corresponding Landau equation, we obtain for any ϕ ∈ C 1 b , by exchangeability, Thus if one simulates the discretized particle system (8), and if one computes 1 n n 1 ϕ(X i,n,N t ), we get an approximation of ϕ(x)P t (dx), with a reasonnable error.
1.8. Plan of the paper. In Section 2, we give the proofs of Theorems 3 and 6. Section 3 is devoted to the proofs of Corollaries 4 and 5. In Section 4, we briefly deal with the case of soft potentials, but our theoritical results do not extend well. Numerical results are given in Section 5. Finally an appendix lies at the end of the paper.

General proofs
In the whole section, we assume that P 0 ∈ P 2 , that a : R d → S + d is of class C 2 , with bounded derivatives of order two, and that b : R d → R d is Lipschitz continuous. We denote by C (resp. C T , C T,p ) a constant which depend only on a, b, d, P 0 (resp. additionally on T , on T, p) and whose value may change from line to line. [13] for many informations on the Wasserstein distance W 2 .

Preliminaries.
Our results are mainly based on the two following Lemmas.
Proof. Step 1. For µ ∈ P 2 fixed, we consider the map A : Step 2. We now fix x ∈ R d , and consider µ, ν ∈ P 2 . We introduce a couple (X, Y ) of random variables such that Step 3. The growth estimate (for a) follows from the Lipschitz estimate, since |a ). The growth estimate follows from the Lipschitz estimate, since |b(0, δ 0 )| 2 = |b(0)| 2 < ∞.

2.2.
Convergence proofs. We start this subsection with some moment estimates.
Using the Burkholder-Davies-Gundy inequality for the Brownian part, and the Hölder inequality for the drift part, we obtain, for all 0 ≤ t ≤ T , , whence the result by the Gronwall Lemma. Point (ii). Using the Cauchy-Scharz and Doob inequalities, we see that for 0 ≤ s ≤ t ≤ T , Proof of Theorem 3. We consider P 0 ∈ P 2 fixed.
Uniqueness. Assume that we have two solutions X, Y to E 1 (P 0 , a, b, X 0 , B), and set P t = L(X t ), Q t = L(Y t ). Using the Cauchy-Schwarz and Doob inequalities, we obtain, for 0 ≤ t ≤ T , We used Lemma 7 and the obvious inequality . The Gronwall Lemma allows us to conclude that X = Y . Existence. We consider the following Picard iteration: set X 0 t = X 0 , and define, for n ≥ 0, t ≥ 0, b(X n s , L(X n s ))ds.
We get as in (11), ds. Thus there classically exists (X t ) t≥0 such that lim n E[sup [0,T ] |X n t − X t | 2 ] = 0 for all T , which implies that lim n sup [0,T ] W 2 2 (L(X n t ), L(X t )) = 0. Passing to the limit in (12), we see that X solves E 1 (P 0 , a, b, X 0 , B).

Ellipticity estimates
We start with the Proof of Corollary 4. Recall here that a is given by (2) with κ ≡ 1 and b(z) = −(d − 1)z. Thus b is Lipschitz continuous, and the second derivatives of a are clearly bounded. We consider a weak solution (P t ) t≥0 to (1).

It remains to give the
Proof of Corollary 5. Recall here that a ij (x) = κ(|x| 2 )(|x| 2 δ ij − x i x j ) and that b(x) = −(d − 1)κ(|x| 2 )x, that κ is C 2 and that κ ′ has a bounded support, so that a has bounded derivatives of order 2, and b is Lipschitz continuous. We consider a weak solution (P t ) t≥0 to (1). As previously, we classically have m 2 (P t ) = m 2 (P 0 ). Furthermore, it is again classical and widely used that the entropy of P t is non-increasing, so that P t (x) log P t (x)dx ≤ P 0 (x) log P 0 (x)dx = C < ∞ for all times, see Villani [10,11,12]. If we prove that there is λ 0 > 0 such that for all t ≥ 0, x, y ∈ R d , (a(x, P t )y, y) ≥ λ 0 |y| 2 , then we deduce that |a(x, P t ) −1 | is uniformly bounded, so that the Corollary follows from (7). Observe that setting α ij (x) = |x| 2 δ ij − x i x j , we have (a(x, P t )y, y) ≥ λ 1 (α(x, P t )y, y), where λ 1 > 0 is a lowerbound of κ. But it is shown in Desvillettes-Villani [3,Proposition 4] that for E 0 ∈ R + , H 0 ∈ R + , there is a constant c E0,H0 > 0 such that for any probability density function f on R d such that m 2 (f ) ≤ E 0 and f (x) log f (x)dx ≤ H 0 , (α(x, f )y, y) ≥ c E0,H0 |y| 2 . Actually, they consider the case where α ij (x) = |x| γ (|x| 2 δ ij − x i x j ) for some γ > 0, but one can check that their proof works without modification when γ = 0. We finally obtain (a(x, P t )y, y) ≥ λ 1 c E0,H0 |y| 2 for all t ≥ 0, x, y ∈ R d , which concludes the proof.
On the other hand, we may apply the techniques introduced in [5] to estimate W 2 2 (P t , P ε t ), where (P t ) t≥0 = (L(X 1 t )) t≥0 is a weak solution to (1) with a and P 0 . We believe that, with a convenient coupling, it is possible to obtain something like sup [0,T ] One would thus get . This is of course an awfull rate of convergence. It does not seem reasonnable to handle a rigorous proof.
Simulation without cutoff. However, the particle system (8) is still well-defined and simulable for soft potentials (with γ ∈ [−3, 0]), at least if we replace 1 n k δ X k,n,N t by 1 n k =i δ X k,n,N t and if P 0 has a density. Based on the well-posedness result of [5], we hope that, at least when γ ∈ (−2, 0], one might obtain the same estimates as in Corollary 5 and Theorem 6 ( under additionnal conditions on P 0 ). The proof however seems to be quite difficult: we do not know how to get a sufficiently good estimate of quantities like |X i,n,N t − X j,n,N t | γ .

Numerics
Let us first observe that for the Landau equation 1 where a is given by (2) and b i = j ∂ j a ij , the simulable particle system (8) is conservative, in the sense that it preserves, in mean, momentum and kinetic energy: We consider here the Landau equation for soft potentials, for some γ ∈ [−3, 0], described in the previous section, in dimension d = 2. We use no cutoff procedure in the case γ < 0. We consider the initial condition P 0 with density P 0 (x 1 , x 2 ) = f (x 1 )g(x 2 ), where f is the Gaussian density with mean 0 and variance 0.1, while g(x) = (f (x − 1) + f (x + 1))/2. The momentum and energy of P 0 are given by (0, 0) and 1.02. Thus in large time, the solution P t should converge to the Gaussian distribution with mean (0, 0) and covariance matrix 0.51I 2 , see Villani [12]. We use the particle system (8) with n particles, and N steps per unit of time. Easy considerations show that the computation of (8) until time T is essentially proportionnal to T N n 2 , and should not depend too much on γ. However, it is consequently faster when γ = 0 for obvious computational reasons. Let us also remark that the law of (8) does not change when replacing a 1 2 by any σ such that σ(x, µ)σ(x, µ) * = a(x, µ). We thus use a Cholesky decomposition, which is numerically quite fast. Let us give an idea of the time needed to perform one time-step: with γ = 0, it takes around 7.10 −3 seconds (n = 500), 0.15 s (n = 2500), 3.5 s (n = 12500), and 13 s (n = 25000). The computations are around 10 times slower when γ < 0. Now we alway use n = 5000 particles, and N = 200 steps per unit of time. We draw, for different values of t and γ, the histogram (with 80 sticks) based on the second coordinates of (X i,n,N t ) i=1,...,n . The plain curve is the expected asymptotic Gaussian density, with mean 0 and variance 0.51. The convergence to equilibrium seems to be slower and slower as γ is more and more negative. For too small values of γ (say γ < −2.5), the numerical results are not so convincing. This is not surprising, since the coefficients are more and more singular as γ becomes smaller and smaller.

Appendix
The following Lemma can be found in Stroock  We also need the following estimates, which are probably standard. We now prove (ii). First observe that (A We conclude this annex with an elementary fact on the Wasserstein distance.