Anaïs Crestetto

Nicolas Crouseilles

Damien Prel

Deterministic particle method for Fokker-Planck equation with strong oscillations

Keywords: Vlasov equation, Fokker-Planck collision operator, highly oscillatory systems, multiscale numerical schemes, Particle method. MSC codes. 65M25, 65M75, 35Q83

) is proposed for Fokker-Planck type operators. This particle method consists in reformulating the collision operator in an advective form and in regularizing the advection field in such a way that it conserves the geometric bracket structure. In the Fokker-Planck homogeneous case, the properties of the resulting method are analysed. In the non homogeneous case, the particle method is coupled with a uniformly accurate time discretization in ε that enables to capture numerically the solution of the asymptotic model. Numerous numerical results are displayed, illustrating the behavior of the method.

Introduction

When modelling charged particles systems, a kinetic description is often used to take into account non equilibrium phenomena. The unknown is a distribution function f (t, x, v) ≥ 0 depending on time t, space x and velocity v. The dynamics of f is driven by the so-called Vlasov models involving a (nonlinear) transport part in phase space (x, v) which may induced high oscillations due to the presence of a strong magnetic or electric field. Collisional effects can also be taken into account to describe binary interactions between charged particles through collisional operators acting on the velocity variable only. Collisional effects make the system relax towards a space dependent equilibrium (typically with a Maxwellian shape in velocity). One example is the Landau operator but due to its complexity (it is a nonlinear integro-differential operator), simpler operators, sharing similar properties, are preferred such as the Fokker-Planck or the Bhatnagar-Gross-Krook (BGK) operators.

The goal of this work is to investigate numerically the interplay between oscillation effects (arising from the transport part) and collisional effects, in the spirit of recent works [START_REF] Bostan | Finite Larmor radius approximation for collisional magnetic confinement. Part I: The linear Boltzmann equation[END_REF][START_REF] Bostan | Finite Larmor radius approximation for collisional magnetic confinement. Part II: The Fokker-Planck-Landau equation[END_REF][START_REF] Bostan | Finite larmor radius regime: Collisional setting and fluid models[END_REF][START_REF] Crestetto | Multiscale numerical schemes for the collisional vlasov equation in the finite larmor radius approximation regime[END_REF]. In particular, in the limit ε → 0 (where ε denotes the period of the oscillations), an asymptotic (or averaged) model can be derived, involving the collision operator. As a consequence, the equilibria of the averaged collision operator are modified, which leads to fluid models that are different from the usual ones. However, even if some explicit expressions of the averaged system can be found for some given collision operators (see [START_REF] Bostan | Finite Larmor radius approximation for collisional magnetic confinement. Part I: The linear Boltzmann equation[END_REF][START_REF] Bostan | Finite Larmor radius approximation for collisional magnetic confinement. Part II: The Fokker-Planck-Landau equation[END_REF][START_REF] Bostan | Finite larmor radius regime: Collisional setting and fluid models[END_REF]), these derivations require a huge amount of computations and leads to very complicated models. Moreover, in these works, the finite Larmor radius scaling was considered which involves a large number of phase space variables (2 in space and 2 in velocity) and thus making the numerical simulations very costly. In the present work, we rather consider a simplified model which contains the main ingredients, that is oscillations in phase space (1 dimension in space and 1 dimension in velocity) and collisions in the velocity direction. This enables us to better understand the interplay between oscillations and dissipation effects which drives the system to modified equilibria of the collisional operator and to simulate the multiscale models using dedicated numerical schemes.

The model we considered in this work shares some similarities with the highly oscillatory charged particles beam models used in [17, 22-25, 31, 33], but here collisional effects are taken into account through a nonlinear Fokker-Planck operator (or Lenard-Bernstein collision operator). Thus, this model may serve as a basis for the understanding of more complex models like the ones arising in strongly magnetized plasmas under the finite Larmor radius approximation [START_REF] Bostan | Finite larmor radius regime: Collisional setting and fluid models[END_REF][START_REF] Crestetto | Multiscale numerical schemes for the collisional vlasov equation in the finite larmor radius approximation regime[END_REF] coupled with the full Landau collision operator. Even if the situations in which oscillations and collisions interplay can be found in many applications (in plasma physics, it turns out that the litterature is not abondant ; let us quote however [START_REF] Ben Abdallah | Diffusion and guiding center approximation for particle transport in strong magnetic fields[END_REF][START_REF] Blanc | Asymptotic analysis of parabolic equations with stiff transport terms by a multi-scale approach[END_REF][START_REF] Bostan | Finite Larmor radius approximation for collisional magnetic confinement. Part I: The linear Boltzmann equation[END_REF][START_REF] Bostan | Finite larmor radius regime: Collisional setting and fluid models[END_REF][START_REF] Herda | Large-time behavior of solutions to vlasov-poisson-fokker-planck equations: From evanescent collisions to diffusive limit[END_REF][START_REF] Herda | Anisotropic boltzmann-gibbs dynamics of strongly magnetized vlasov-fokkerplanck equations[END_REF]).

First, one has to study the averaged model of our beam-collision model. To do so, denoting by ε a parameter characterizing the period of oscillations (2πε) induced by the transport part, a change of variables (usually called filtering step) is performed to get the filtered system ; then, the (strong) limit ε → 0 can be performed leading to the averaged model (see [START_REF] Bostan | Finite larmor radius regime: Collisional setting and fluid models[END_REF][START_REF] Chartier | High order averaging, formal series and numerical integrations i: B-series[END_REF][START_REF] Chartier | A formal series approach to averaging: exponentially small error estimates[END_REF]). This asymptotic model corresponds to an averaging of the filtered collision operator on one oscillation period. Using [START_REF] Bostan | Finite larmor radius regime: Collisional setting and fluid models[END_REF] in which a quite general framework is presented, we can identify the equilibria of the averaged collision kernel, without having to compute explicitly the averaged collision kernel. In particular, the equilibrium has Maxwellian shape not only on velocity (as in the usual case) but also on the space variable due to the oscillations occurring in phase space (x, v)

Second, our goal is to derive new numerical schemes to approximate the averaged model. But as its explicit expression is difficult, one may prefer to approximate the stiff filtered problem with a multiscale (uniformly accurate UA) numerical method [9-13, 17, 19]. Indeed, for such UA methods, the error does not depend on the value of the stiffness parameter ε and as such the schemes are asymptotic preserving. This UA approach is dedicated to the time approximation. Regarding the phase space approximation, which involves a transport part and a diffusion part, particle methods will be used in this work. Indeed, particle discretizations (such as Particle-In-Cell (PIC) method) are very popular in the computational plasma community since they are very efficient in high dimensions compared to Eulerian methods and they are able to capture filamentations with a reasonable computation cost. Moreover, they are easy to develop for transport equation [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF][START_REF] Sonnendrücker | Numerical Methods for the Vlasov-Maxwell equations[END_REF]. However, their extensions to the approximation of collisional operators (which involves diffusion operators) are more tricky. One way consists in solving a stochastic differential equation (SDE) associated to the diffusion term. However, these methods are low order in time and the conservation properties are not guaranteed. Recently, a deterministic particle method has been developed to discretize diffusion operators [START_REF] Carrillo | A blob method for diffusion[END_REF]. This method consists in interpreting the diffusion operator as a nonlinear advection operator and then in regularizing the advection field with a mollifier making the advection field well defined for particles unknown. A variant of this method has been developed for Landau operator [START_REF] Carrillo | A particle method for the homogeneous landau equation[END_REF][START_REF] Hirvijoki | Structure-preserving marker-particle discretizations of coulomb collisions for particle-in-cell codes[END_REF][START_REF] Zonta | Multispecies structure-preserving particle discretization of the landau collision operator[END_REF]: it consists in regularizing the entropy associated to the operator. By doing this, the geometric bracket structure of the equation is preserved, which directly implies the decreasing of the regularized entropy and the conservation of the invariants. Let us remark that similar ideas were also developed in [START_REF] Degond | The weighted particle method for convection-diffusion equations. I. The case of an isotropic viscosity[END_REF][START_REF] Degond | A deterministic approximation of diffusion equations using particles[END_REF].

In this work, due to its simpler form, the Fokker-Planck operator is considered. Since it also enjoys a geometric structure (bracket and entropy), we propose to adapt the formalism from [START_REF] Carrillo | A particle method for the homogeneous landau equation[END_REF][START_REF] Hirvijoki | Structure-preserving marker-particle discretizations of coulomb collisions for particle-in-cell codes[END_REF][START_REF] Zonta | Multispecies structure-preserving particle discretization of the landau collision operator[END_REF], which requires some specific developments due to the presence of macroscopic moments in the definition of the Fokker-Planck operator. We also mention a recent paper [START_REF] Jeyakumar | A structure-preserving particle discretisation for the lenard-bernstein operator[END_REF] in which the authors derive a different regularization technique that enables to preserve mass, momentum and energy. It turns out that the deterministic particle method is well adapted to our context since it enables to reproduce at the semi-discrete level (discretization in space and velocity) the main steps of the derivation of the averaged system. In particular, the deterministic particle method can naturally take into account the filtering step so that it can be combined with the two-scale framework developed in [9-13, 17, 19]. As a result, a uniformly accurate method (with respect to ε) can be constructed. The error of this scheme is independent of the oscillation period ε, and then the scheme is consistent with the asymptotic averaged model obtained when ε goes to 0.

This numerical framework enables us to numerically illustrate the modified equilibrium (called gyromaxwellian) of the averaged collision operator and also to investigate the difference between the so-called filtered Maxwellian and the gyromaxwellian. This difference has already been discussed in [START_REF] Crestetto | Multiscale numerical schemes for the collisional vlasov equation in the finite larmor radius approximation regime[END_REF] in the finite Larmor radius scaling. Some numerical tests are also performed for the Fokker-Planck operator (in the homogeneous case) to discuss the performances of the deterministic particle method in this context. And finally, some tests involving oscillations and collisional effects are performed.

The outline of the paper is the following. In Section 2, we present our model and study the associated gyrokinetic model. In particular, we focus on characterizing equilibriums and collisional invariants of the filteredaveraged collision operator. In Section 3, we present the particle discretization of the Fokker-Planck equation, and some properties associated. Section 4 is devoted to the presentation of the uniformly accurate scheme for the Fokker-Planck equation with oscillations. Numerical results illustrating the particle method for Fokker-Planck equation are presented in Section 5 whereas numerical results for full model and beam model are presented in Section 6. Some proofs are given in Appendices A and B.

Model and properties

In this section, we consider a collisional kinetic model composed of two parts: a highly oscillatory (in phase space) transport term and a Fokker-Planck collision operator acting only in the v-direction. In Subsection 2.1, we present more precisely the model and derive the asymptotic model (i.e. when the period ε of the oscillations goes to zero). This asymptotic requires a change of variables called filtering step and averaging procedure which have to be applied to the Fokker-Planck operator. Thus, Subsection 2.2 is devoted to the study of some properties of the filtered-averaged Fokker-Planck operator.

Presentation of the Fokker-Planck equation with oscillations

We are interested in a kinetic model with a Fokker-Planck collision operator and strong oscillations in space and velocity. We set f (t, x, v) ≥ 0 the particles distribution function depending on time t ∈ R + , space x ∈ R and velocity v ∈ R. We denote z = (x, v) ∈ R 2 the phase space variable and introduce ε > 0 which denotes a small parameter (typically a ratio between two characteristic lengths). The model we are interested in is the following

∂ t f - 1 ε Jz • ∇ z f = ∂ v ((v -u)f + T ∂ v f ) (1)
where J denotes the symplectic matrix

J = 0 -1 1 0
and the macroscopic quantities ρ (density), u (mean velocity) and T (temperature) are defined from the moments of f by

ρ = f dv, ρu = vf dv, ρ u 2 2 + T 2 = v 2 2 f dv. (2)
The model ( 1) is equipped with a smooth initial condition f (t = 0, x, v) = f 0 (x, v) which is fastly decreasing at infinity.

In the following, we will use this notation for the Fokker-Planck collision operator

Q[f ] = ∂ v ((v -u)f + T ∂ v f ) = ∂ v (T f ∂ v log (f /M)) , (3)
where M is the Maxwellian associated to f (through the moments (2))

M[f ](t, x, v) = ρ √ 2πT e -|v-u| 2 2T . (4)
It is interesting to know how the distribution function f behaves as ε → 0, i.e. when the frequency of the oscillations becomes infinite. To study this limit, we consider the filtered distribution function

F (t, Z) = f (t, z), where Z = e t ε J z, z = (x, v), Z = (X, V ), which satisfies the filtered problem ∂ t F (t, Z) = Q filt [F ](t, t/ε, Z), (5)
where Q filt [F ](t, s, Z) denotes the filtered collision operator (obtained by applying the change of variables to (3))

Q filt [F ](t, s, Z) = ∇ Z • (T A(s)∇ Z F ) + ∇ Z • (A(s) (Z -a(s)u) F ) , (6)
with

A(s) = sin 2 (s) -sin(s) cos(s) -sin(s) cos(s) cos 2 (s) , a(s) = (-sin(s), cos(s)) T , b(s) = (cos(s), sin(s)) T , such that x = b(t/ε) • Z and v = a(t/ε) • Z.
This filtering step is important to identify the asymptotic limit. Indeed, under some assumptions, F solution of ( 5) is known to strongly converge to F solution of the averaged model (see [START_REF] Chartier | High order averaging, formal series and numerical integrations i: B-series[END_REF][START_REF] Chartier | A formal series approach to averaging: exponentially small error estimates[END_REF][START_REF] Crouseilles | Uniformly accurate Particle-In-Cell method for the long time two-dimensional Vlasov-Poisson equation with strong magnetic field[END_REF]) given by

∂ t F (t, Z) = Q filt [F ](t, Z), (7)
where Q filt is the averaged of Q filt on the fast variable s, which is present in u, T , A and a. It is defined by

Q filt [F ](t, Z) = 1 2π 2π 0 Q filt [F ](t, s, Z) ds. ( 8
)
Due to the nonlinear character of the collision operator, the filtered-averaged collision operator is not explicit. However, as we shall see in the next subsection, some properties can be derived.

Let us remark that the filtered collision operator Q filt [F ] given by ( 6) can also be expressed in the following divergence form

Q filt [F ](t, s, Z) = ∇ Z • T A(s)F ∇ Z log F/M filt [F ](s) ,
where M filt [F ] denotes the Maxwellian sharing the same three first moments of f , but expressed in the filtered variables

M filt [F ](t, s, Z) = ρ √ 2πT exp - 1 2T |a(s) • Z -u| 2 . ( 9
)
In the definition of Q filt and M filt , the macroscopic quantities ρ, u, T have to be computed from the moments (2) of the unfiltered function f associated to F , and then these x-dependent quantities are evaluated at b(s) • Z (which corresponds to x for s = t/ε). For example, we have the following relations:

ρ(x) = ρ(b(s) • Z) but also Q filt [F ](t, t/ε, X, V ) = Q[f ](t, x, v) and M filt [F ](t, t/ε, X, V ) = M[f ](t, x, v).
We also define the filtered-averaged Maxwellian

M filt [F ](t, Z) = 1 2π 2π 0 M filt [F ](t, s, Z) ds. ( 10
)
Remark 1. In some applications, one can assume u = 0 and T = 1 in (3) which corresponds to the interactions of the particles with a bulk with a given mean velocity and temperature. In this case, the Fokker-Planck operator is linear and the filtered-averaged operator can be explicitly computed as a 2-dimensional Fokker-Planck operator.

Indeed, in this case, we have

1 2π 2π 0 A(s) ds = 1 2 Id and then Q filt [F ] = 1 2 (T ∆ Z F + ∇ Z • (ZF )) .

Properties of the asymptotic model

We recall that, for a collision operator Q, a collisional invariant is a function c(x, v) such that for all distribution function f , Q[f ]c(x, v) dx dv = 0. Moreover, it is known that for the Fokker-Planck collision operator (3), the exponential of collisional invariants are the equilibriums. This comes from the entropic structure (see [START_REF] Bostan | Finite larmor radius regime: Collisional setting and fluid models[END_REF] for details). In the following, some properties of the filtered-averaged operator Q filt are given such as its collisional invariants and equilibrium. Several techniques to compute them have been introduced in [START_REF] Bostan | Finite larmor radius regime: Collisional setting and fluid models[END_REF] but we have to adapt this framework to our specific context.

We introduce the notation C s (Z) = C(e sJ Z) for a function C : R 2 → R 2 (similarly F s (Z) = F (e sJ Z) for F : R 2 → R).

Proposition 1. Let F : R 2 → R be a function, Q be the Fokker-Planck operator (3) and Q filt be the filtered Fokker-Planck operator [START_REF] Carrillo | A particle method for the homogeneous landau equation[END_REF]. Then we have the following assertions.

1. F is an equilibrium of Q filt if and only if F s is an equilibrium of Q, ∀s ∈ R. 2. C is a collisional invariant of Q filt if and only if C s is a collisional invariant of Q, ∀s ∈ R. 3. If C is a collisional invariant of Q filt , then V ∂ X C -X∂ V C is a collisional invariant of Q filt . 4. F is an equilibrium of Q filt if and only if log(F ) is a collisional invariant of Q filt .
Proof. The proof of the two first points is similar to [START_REF] Bostan | Finite larmor radius regime: Collisional setting and fluid models[END_REF] and will not be given. Let consider the third assumption. We first remark that

d ds C s (Z) = d ds C(e sJ Z) = Je sJ Z•(∇ Z C)(e sJ Z) = Je sJ Z•(∇ Z C) s (Z) = (JZ • ∇ Z C) s (Z) = (-V ∂ X C+X∂ V C) s (Z).
Assuming C is a collisional invariant of Q filt , then by the second assertion, it is equivalent to say C s is a collisional invariant of Q for all s ∈ R. Thus for all h > 0, C s+h -Cs h is also a collisional invariant by linearity and considering the limit h → 0 gives

d ds C s = (-V ∂ X C + X∂ V C) s is a collisional invariant of Q for all s ∈ R.
According to the second assertion, this is equivalent to say that (-V ∂ X C + X∂ V C) is a collisional invariant of Q filt . We prove now assumption 4. Let F be an equilibrium of Q filt . From assumption 1, it means that F s is an equilibrium of Q for all s ∈ R. It is equivalent to say that log(F s ) is a collision invariant of Q and, using assumption 2, it means that log(F ) is a collision invariant of Q filt .

The following theorem enables us to derive the invariants and the equilibrium of Q filt . Theorem 1. The collisional invariants of Q filt are linear combinations of 1, X, V, X 2 , XV and V 2 . Let the mass ϱ > 0, the mean velocity U ∈ R 2 and the symmetric positive-definite temperature tensor T ∈ R 2×2 , defined from the gyromoments of F :

ϱ = R 2 F (X, V ) dX dV, (11a) ϱU = ϱ U X U V = R 2 X V F (X, V ) dX dV, ( 11b
) ϱT = ϱ T XX T XV T XV T V V = R 2 X -U X V -U V ⊗ X -U X V -U V F (X, V ) dX dV. ( 11c
)
Then the only equilibrium sharing the same gyromoments as F is given by an anisotropic Maxwellian (called gyromaxwellian)

G[F ](Z) = ϱ det(2πT) e -1 2 (Z-U ) T T -1 (Z-U ) , Z = (X, V ) ∈ R 2 , (12)
that is,

R 2 C(X, V )(F (X, V ) -G[F ](X, V )) dX dV = 0, for C(X, V ) = (1, X, V, X 2 , XV, V 2 ) T .
Proof. First, we determine the collisional invariants. Let C(X, V ) be an invariant of Q filt . From the second assertion of Proposition 1, we know that C s is an invariant of Q for all s ∈ R. In particular for s = 0 we have that

C = C 0 is a collisional invariant of Q. So C(X, V ) = α(X)V 2 + β(X)V + γ(X), (13)
where α, β and γ are smooth functions. From the third assertion of Proposition 1, V ∂ X C -X∂ V C is also a collisional invariant of Q filt , and thus of Q. This means that

V ∂ X C(X, V ) -X∂ V C(X, V ) = α(X)V 2 + β(X)V + γ(X), (14)
where α, β and γ are three other smooth functions. Plugging [START_REF] Chartier | Derivative-free high-order uniformly accurate schemes for highly oscillatory systems[END_REF] in [START_REF] Chartier | High order averaging, formal series and numerical integrations i: B-series[END_REF], we get

α ′ (X)V 3 + (β ′ (X) -α(X))V 2 + (γ ′ (X) -2Xα(X) -β(X))V -( γ(X) + Xβ(X)) = 0.
The coefficients of this polynomial function in V should vanish, so that α ′ (X) = 0, (15a)

β ′ (X) = α(X), (15b)
γ ′ (X) = 2Xα(X) + β(X), (15c)
γ(X) = -Xβ(X). (15d)
From (15a), we deduce that α is constant, and so the coefficient in front of V 2 in (13) is constant. Since C represents any invariant of Q filt , all of its invariants will have a constant coefficient in front of

V 2 . Since V ∂ X C -X∂ V C is an invariant of Q filt , α is constant.
We can thus deduce from (15b) that β(X) is a polynomial of degree 1 (and thus the same for β). From (15c), we deduce in the same way that γ (and γ) is a polynomial of degree 2. We have now the generic form of collisional invariants of Q filt by replacing α, β and γ in [START_REF] Chartier | Derivative-free high-order uniformly accurate schemes for highly oscillatory systems[END_REF] and one can conclude they are linear combinations of 1, X, V, X 2 , XV and V 2 . We verify that this is indeed collisional invariants by using Proposition 1: compose a linear combination of these invariants with the change of variable Z → e sJ Z for all s gives an invariant of Q.

From the fourth assertion of Proposition 1, the equilibrium are exponential of collisional invariants which are linear combinations of 1, X, V, X 2 , XV and V 2 . Since this is the case for the function [START_REF] Chartier | Numerical methods for the two-dimensional Vlasov-Poisson equation in the finite Larmor radius approximation regime[END_REF], it remains to show that G[F ] and F share the same gyromoments and T is positive definite. This can be done using classical computations and is postponed in Appendix A.

Remark 2. It is well known (see [START_REF] Bostan | Finite larmor radius regime: Collisional setting and fluid models[END_REF]) that strong oscillations bring dissipation also in space, so it is not surprising to get Maxwellian-type equilibrium which depends on both space and velocity.

Remark 3. We can define an operator G : F → G[F ] defined by [START_REF] Chartier | Uniformly accurate methods for three dimensional vlasov equations under strong magnetic field with varying direction[END_REF] and [START_REF] Chartier | Numerical methods for the two-dimensional Vlasov-Poisson equation in the finite Larmor radius approximation regime[END_REF]. We can do the same with the filtered-averaged Maxwellian M filt : F → M filt [F ] defined by [START_REF] Chartier | Uniformly accurate methods for Vlasov equations with non-homogeneous strong magnetic field[END_REF]. It is worth to note that G[F ] and M filt [F ] are different except when F is a gyromaxwellian. This is related to the non-commuting limits t → ∞ and ε → 0 in the filter model [START_REF] Bostan | Finite Larmor radius approximation for collisional magnetic confinement. Part II: The Fokker-Planck-Landau equation[END_REF], this has been discussed in [START_REF] Crestetto | Multiscale numerical schemes for the collisional vlasov equation in the finite larmor radius approximation regime[END_REF]. We detail this point further in Section 6.1.

It is worth to note that gyromoments [START_REF] Chartier | Uniformly accurate methods for three dimensional vlasov equations under strong magnetic field with varying direction[END_REF] are also conserved when F solves the filtered equation ( 5) for all ε > 0. Indeed, from Assertion 2 of Proposition 1, the collisional invariants of Q filt are also collisional invariants of Q when they are unfiltered form.

Proposition 2. ∀ε ≥ 0, the gyromoments [START_REF] Chartier | Uniformly accurate methods for three dimensional vlasov equations under strong magnetic field with varying direction[END_REF]

are conserved ( d dt (ϱ, ϱU, ϱT) = 0) and d dt G[F ] = 0. Proof. First, we have Q[f ]c(z) dz = 0 for c(z) = (1, x, v, x 2 , xv, v 2
) and we recall that Z = e t ε J z, ∀ε > 0 so 1, X, V, X 2 , XV, V 2 are linear combinations of 1, x, v, x 2 , xv, v 2 (with coefficients depending on t/ε). Then, we get

Q filt [F ]c(Z) dZ = Q[f ]c(e t ε J
z) dz = 0 by linearity. Integrate (5) against c(Z) directly proves the conservation of gyromoments. As the gyromaxwellian G[F ] depends on time only through gyromoments, we get d dt G[F ] = 0. For ε = 0, c(Z) is the vector of collisional invariants of Q filt (from Theorem 1) so the conservation is obtained by definition.

Particle discretization for the Fokker-Planck operator

In this section, for the sake of clarity, we first focus on the Fokker-Planck operator in the velocity direction and discuss its numerical approximation using particle methods which is well adapted for the coupling of the transport part that will be described in the next section. Without adding complexity, we consider v ∈ R d where d is the dimension, and all the integrals are considered over R d . As mentioned above, we consider particle methods in this work and will adapt the numerical method developed for the Landau operator in [START_REF] Carrillo | A particle method for the homogeneous landau equation[END_REF][START_REF] Hirvijoki | Structure-preserving marker-particle discretizations of coulomb collisions for particle-in-cell codes[END_REF][START_REF] Zonta | Multispecies structure-preserving particle discretization of the landau collision operator[END_REF] to the Fokker-Planck operator. To do so, some properties of the Fokker-Planck operator will be recalled in Subsection 3.1. In particular, the use of a metric bracket and an entropy enables us to define a regularized Fokker-Planck operator and to study its properties in Subsection 3.2. Then, this regularized operator will pave the way of deterministic particles method presented in Subsection 3.3 whereas Subsection 3.4 presents conservation properties at the fully discrete level (including time discretization), and finally Subsection 3.5 is devoted to the presentation of a variant of the method allowing parallelization.

The Fokker-Planck collision operator and its properties

Here, we focus on the collisional operator so we consider a density f (t, v) solution of

∂ t f = Q(f ) = T ∆ v f + ∇ v • (v -u)f , f (t = 0, v) = f 0 (v), v ∈ R d . ( 16
)
As there is no transport in ( 16), the moments f (t, v)(1, v, |v| 2 )dv are conserved with time, so are the mean velocity u and temperature T which are defined from the moments of the initial condition f 0 (v). In view of the deterministic particle method, we start with a reformulation of ( 16) using a metric bracket and a driven entropy (see also [START_REF] Hirvijoki | Structure-preserving marker-particle discretizations of coulomb collisions for particle-in-cell codes[END_REF][START_REF] Krauss | Metriplectic integrators for the landau collision operator[END_REF]). In the Fokker-Planck case, the bracket is given by, for two functionals F and G depending on f ,

(F, G) (f ) = -f T ∇ v δF δf • ∇ v δG δf dv, (17)
where δF δf is the Fréchet derivative of F. The driven entropy S M is the relative entropy to the Maxwellian M:

S M [f ] = f log(f /M) dv.
The evolution of any functional F of f is thus given by

∂ t F(f ) = (F, S M ) (f ). ( 18
)
For example, ( 16) can be recovered using the following functional F(f

)(v) = δ(v -v ′ )f (v ′ ) dv ′ where δ(v)
denotes the Dirac distribution. We also define the physical entropy

S[f ] = f log(f ) dv. ( 19
)
The following Theorem recalls some standard properties of the Fokker-Planck operator which will be proven using the metric bracket formalism.

Theorem 2. (Properties of the collision operator)

Let M be a Maxwellian characterized by ρ, u and T and let f be the solution of ( 16) which has the same mass, mean velocity and temperature as M. We have the following properties

d dt f dv = 0, d dt vf dv = 0, d dt |v| 2 2 f dv = 0, d dt S M ≤ 0.
Moreover, the only equilibrium of the Fokker-Planck operator

Q is f = M.
Proof. The proof, using bracket formulation [START_REF] Crouseilles | Uniformly accurate forward semi-Lagrangian methods for highly oscillatory Vlasov-Poisson equations[END_REF], is given in Appendix B.1.

We then consider [START_REF] Crestetto | Multiscale numerical schemes for the collisional vlasov equation in the finite larmor radius approximation regime[END_REF] in the following advective form

∂ t f = ∇ v • (f T ∇ v log(f /M)) = ∇ v • f T ∇ v δS M δf , (20)
where the advective field is the functional derivative of the relative entropy

T δS M δf = T ∇ v log(f /M).

Regularized Fokker-Planck operator

In view of using particle methods, [START_REF] Degond | The weighted particle method for convection-diffusion equations. I. The case of an isotropic viscosity[END_REF] is not well adapted due to the presence of the logarithm. The idea is to regularize the distribution function using convolutions with a mollifier [START_REF] Carrillo | A blob method for diffusion[END_REF][START_REF] Carrillo | A particle method for the homogeneous landau equation[END_REF][START_REF] Hirvijoki | Structure-preserving marker-particle discretizations of coulomb collisions for particle-in-cell codes[END_REF][START_REF] Zonta | Multispecies structure-preserving particle discretization of the landau collision operator[END_REF] given by

ψ µ (v) = 1 (2πµ) d/2 e -|v| 2 2µ , (21)
where µ > 0 is the regularization parameter. In order to keep the metric bracket structure (which implies the dissipation of relative entropy), we apply the regularization directly on the relative entropy. In other words, we replace S M by the regularized relative entropy S µ M given by

S µ M [f ] = f * ψ µ log(f * ψ µ /M) dv, (22)
which is now well defined for particle densities, as we shall see in the next subsection. The evolution equation associated to is given by (the notation f is kept for the unknown)

∂ t f = ∇ v • f T ∇ v δS µ M δf = ∇ v • f T ψ µ * ∇ v log(f * ψ µ /M) . (23)
We can observe that the drift part is not affected by the regularization. Indeed, we have ∇ v log(M) = -(v -u)/T so using properties of Gaussian, we have

-ψ µ * ∇ v log(M) = - 1 T w ψ µ (v -w) dw -u ψ µ (v -w) dw = - v -u T .
Consequently, ( 23) can be equivalently written as

∂ t f = (f, S µ M ) = ∇ v • f T ψ µ * ∇ v log(f * ψ µ ) + (v -u) . (24)
As we have modified the collision operator, we can naturally ask if the new operator shares the same properties as the original one. Even if mass and momentum are preserved, this is not the case for the energy which is dissipated, so is the temperature. As a consequence the temperature of the Maxwellian and the temperature of f will differ, this is why, in the following, we introduce different notations for these two quantities, namely T M and T f . Temperature T in the bracket [START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF] is replaced by T M so that we recover [START_REF] Frénod | Long time simulation of a beam in a periodic focusing channel via a two-scale pic-method[END_REF] with T = T M . Proposition 3. (Properties of the regularized collision operator) Let M be a Maxwellian characterized by ρ, u and T M and let f be the solution of [START_REF] Frénod | Long time simulation of a beam in a periodic focusing channel via a two-scale pic-method[END_REF] with mass ρ, mean velocity u and temperature T f = 1 dρ |v -u| 2 f dv. Initially, we assume T f = T M at t = 0. Then we have the following properties

d dt f dv = 0, d dt vf dv = 0, d dt S µ M ≤ 0, (25)
and d dt

|v| 2 2 f dv = dρ T M -T f + µT M D(f ), with D(f ) = 1 T M (S, S) (f * ψ µ ). ( 26
)
Moreover, the only equilibrium of the regularized Fokker-Planck operator is characterized by f * ψ µ = M.

Proof. The proof, using bracket formulation, is given in Appendix B.2.

During the evolution of ( 24), T M (t) = T M (0) whereas T f is not conserved. This is due to (26): the energy is not preserved by the regularized collision operator. However, we have an explicit expression of its evolution which motivates the following modification.

Indeed, we propose a technique to get the energy conservation. Instead of considering a constant (Maxwellian) temperature T M , we choose it in such a way that the time derivative of the energy in (26) equals 0. We thus take

T M = dρ dρ + µD(f ) T f . (27)
Remark 4. D(f ) is independent from T M and is negative from the negativity of the bracket.

The following Proposition and Corollary present a different way to express D(f ). This new formulation will be cheaper and easier to compute when considering particle discretization (detailed in Subsection 3.3). Suppose that g solves the heat equation, which means, using the bracket formulation and the physical entropy S given by [START_REF] Crouseilles | Uniformly accurate Particle-In-Cell method for the long time two-dimensional Vlasov-Poisson equation with strong magnetic field[END_REF], that ∂ t g = (g, S)(g). We deduce that the time evolution of the entropy S[g] is then given by ∂ t S[g] = (S, S)(g). To compute D(f ) = (S, S)(f * ψ µ ), it is tempting to choose g = f * ψ µ , but f is not the solution of the heat equation on v ∈ R d (with vanishing condition at infinity). However, the following Proposition enlightens us. Proposition 4. Let g t (τ, v) be the solution of the heat equation

∂ τ g t = ∆ v g t (v ∈ R d ), with initial condition g t (τ = 0, v) = (f * ψ µ )(t, v). Then D(f ) = d dτ S[g t ]| τ =0 .
Proof. Remark that since we have divided the bracket by T M , we consider T M = 1 in the proof without loss of generality. Since the function g t solves the heat equation, we have equivalently ∂ τ g t = (g t , S)(g t ) and the evolution of S[g t (τ, •)] is given by d dτ S[g t (τ, •)] = (S, S)(g t (τ, •)). Thanks to the choice of the initial condition

g t (τ = 0, v) = (f * ψ µ )(t, v), taking τ = 0 in the latter equation gives d dτ S[g t (τ, •)]| τ =0 = (S, S)(g t (τ = 0, •)) = (S, S)((f * ψ µ )(t, •)) = D(f ).
Thanks to Proposition 4, we can deduce an explicit expression for D(f ) defined in [START_REF] Herda | Large-time behavior of solutions to vlasov-poisson-fokker-planck equations: From evanescent collisions to diffusive limit[END_REF].

Corollary 1. D(f ) satisfies D(f ) = (f * ∆ v ψ µ ) log(f * ψ µ ) dv with ∆ v ψ µ = 1 µ |v| 2 µ -d ψ µ .
Proof. Let g t defined in Proposition 4, then we have

d dτ S[g t ] = R d ∂ τ g t (log(g t ) + 1) dv = R d ∆ v g t log(g t ) dv.
Evaluating at τ = 0 and using the initial condition g t (τ = 0, v) = (f * ψ µ )(t, v) from Proposition 4 leads to

D(f ) = d dτ S[g t ]| τ =0 = (f * ∆ v ψ µ ) log(f * ψ µ ) dv.
Using the expression (21) of ψ µ directly gives ∆ v ψ µ .

Particle discretization

In this part, we present the particle method. Starting from [START_REF] Frénod | Long time simulation of a beam in a periodic focusing channel via a two-scale pic-method[END_REF], we approximate the continuous distribution function f by a particle density given by

f Np (t, v) = ρ N p Np p=1 δ(v -v p (t)), (28)
where ρ = f 0 (v)dv is the total mass of the initial distribution function, N p is the number of particles and v p is the velocity of the p-th particle. Inserting (28) into (24) leads to the following set of ODEs for each p = 1, . . . , N p ,

vp (t) = -U µ [f Np ](t, v p (t)), with U µ [f Np ](t, v p ) = T ψ µ * ∇ v log f Np * ψ µ (v p ) + (v p -u). (29)
The particle density is solution in the sense of distribution of the regularized equation [START_REF] Frénod | Long time simulation of a beam in a periodic focusing channel via a two-scale pic-method[END_REF]. In this work, we choose to take particles with same weights ρ/N p but the method also works with non constant weights. The computation of the regularized part of the advection field U µ [f Np ](t, v q ) can be done using a Gauss-Hermite quadrature. Indeed, using a change of variable, the first part can be written as

ψ µ * ∇ v log f Np * ψ µ (v p ) = (∇ v ψ µ ) * log f Np * ψ µ (v p ) = 2 µπ d R d e -|w| 2 w log f Np * ψ µ (v p + 2µ w) dw = R d e -|w| 2 g(w) dw, ( 30
)
where

g(w) = 2 µπ d w log f Np * ψ µ (v p + √ 2µ w).
The form of the integral motivates the use of Gauss-Hermite quadrature:

R d e -|w| 2 g(w) dw ≈ Nq i1,i2,...,i d =1 ω i1 ω i2 . . . ω i d g(c i1 , c i2 , . . . , c i d ),
where N q denotes the number of weights and points, ω i and c i respectively denote the weights and points of the quadrature. Other quantities (like the entropies S and S M ) will be approximated using Gauss-Hermite quadrature using similar computations. Another example is the calculation of the correction ( 27): we need Corollary 1 and with f = f Np and second point of Lemma 1 in Appendix B.2, we get

D(f Np ) = ρ N p Np p=1 ∆ v ψ µ * log f Np * ψ µ (v p ),
where the first convolution is approximated using Gauss-Hermite quadrature.

Time discretization

In this section, we investigate the impact of a forward Euler time discretization applied to [START_REF] Hochbruck | Exponential integrators[END_REF]. Let us mention that implicit integrators are used in [START_REF] Zonta | Multispecies structure-preserving particle discretization of the landau collision operator[END_REF] for Landau operator and [START_REF] Jeyakumar | A structure-preserving particle discretisation for the lenard-bernstein operator[END_REF] for Lenard-Bernstein operator where similar strategies are employed. Let consider v n p ≈ v p (t n ) with t n = n∆t and ∆t > 0 the time step, and f n Np given by (28) where v p (t) are replaced by v n p . Then we consider the following time discretization

v n+1 p -v n p ∆t = -T M ∇ v ψ µ * log f n Np * ψ µ (v n p ) -(v n p -u). (31)
Theorem below shows conservation (or non-conservation) of discrete moments. We consider that the convolutions are computed exactly (and not approximated with Gauss-Hermite quadrature).

Theorem 3. For n ∈ N, let define the discrete mass, momentum and energy as

ρ n = R d f n Np (v) dv, P n = R d f n Np (v)v dv = ρ N p Np p=1 v n p , E n = R d f n Np (v) |v| 2 2 dv = ρ N p Np p=1 v n p 2 2 ,
and consider that particles evolution is given by numerical scheme [START_REF] Knorr | Asymptotic state of the finite-larmor-radius guiding-centre plasma[END_REF]. Suppose that initially, ρ 0 = ρ and P 0 = ρu. Then ∀n ∈ N we have

ρ n = ρ, P n = ρu, E n+1 -E n ∆t = dρ T M -T f,n + µ (S, S) (f n Np * ψ µ ) - ∆t 2 T M (S µ M , S µ M ) (f n Np ), ( 32
)
where T f,n is defined by E n = ρ |u| 2 + dT f,n /2, T M is the constant temperature of M and the bracket, S and S µ M are respectively defined by (17), ( 19) and [START_REF] Degond | The paraxial approximation of the vlasov-maxwell equations[END_REF].

Proof. It is clear that the mass ρ is preserved. For momentum, apply ρ Np Np p=1 on [START_REF] Knorr | Asymptotic state of the finite-larmor-radius guiding-centre plasma[END_REF] gives

P n+1 -P n ∆t = -T M ρ N p Np p=1 ∇ v ψ µ * log f n Np * ψ µ (v n p ) -(P n -ρu) = -T M R d f n Np ∇ v ψ µ * log f n Np * ψ µ (v) dv -(P n -ρu) = -T M R d f n Np * ψ µ ∇ v f n Np * ψ µ f n Np * ψ µ dv -(P n -ρu) = -(P n -ρu) .
Since P 0 = ρu, we show by induction P n = ρu for all n ∈ N.

Regarding the energy, we multiply ( 31) by (v n p + v n+1 p )/2 and sum over p to get

E n+1 -E n ∆t = -T M ρ 2N p Np p=1 v n p + v n+1 p • ∇ v ψ µ * log f n Np * ψ µ /M (v n p ) = -T M ρ N p Np p=1 v n p • ∇ v ψ µ * log f n Np * ψ µ /M (v n p ) + ∆t (T M ) 2 2 ρ N p Np p=1 ∇ v ψ µ * log f n Np * ψ µ /M (v n p ) 2 = -T M R d f n Np v • ∇ v ψ µ * log f n Np * ψ µ /M (v) dv + ∆t (T M ) 2 2 R d f n Np ∇ v ψ µ * log f n Np * ψ µ /M (v) 2 dv.
The first term is the same as in the regularized case (see Appendix B.2). For the second term, we recognize the bracket (S µ M , S µ M ) (f n Np ). Eventually, we get [START_REF] Krauss | Metriplectic integrators for the landau collision operator[END_REF].

In conclusion, the time discretization adds an extra term which also contributes to energy dissipation. As in the semi-continuous case, we will use [START_REF] Krauss | Metriplectic integrators for the landau collision operator[END_REF] in order to propose a modification that enables to conserve energy at the fully discrete level. This brings a new correction of T M which will be computed at each iteration (and thus denoted by T M,n ). We define

D n = 1 T M (S, S) (f n Np * ψ µ ) and D n = 1 T M (S µ M , S µ M ) (f n Np ). ( 33
)
The division by T M makes D n independent of T M whereas D n still depends on T M only through M.

Proposition 5. Scheme (31) conserved the energy (E n = E 0 with E n defined in Theorem 3) under the following correction of the Maxwellian temperature

T M,n = dρ + µD n -(dρ + µD n ) 2 -2dρ∆tT f,n D n ∆t D n . ( 34
)
Proof. Using (33), ( 32) can be reformulated as

E n+1 -E n ∆t = - ∆t 2 D n (T M ) 2 + (dρ + µD n ) T M -dρT f,n ,
and the energy is preserved if T M is defined as the root of the second order polynomial

T M ± = dρ + µD n ± (dρ + µD n ) 2 -2dρ∆tT f,n D n ∆t D n .
Notice that the discriminant is positive since D n ≤ 0 due to the property of the bracket. We select a positive solution T M -(the numerator and the denominator are negative in this case).

Remark 5. It is clear from (32) that when ∆t → 0 or µ → 0, the energy preservation is recovered. We can also observe that when µ → 0, the energy is preserved up to the time discretization error, and it is actually also true for the temperature correction. Indeed, from (34), we get

T M = T f,n + O(∆t) since (dρ) 2 -2dρ∆tT f,n D n = dρ 1 -2 ∆tT f,n D n dρ = dρ -∆tT f,n D n + O(∆t 2 ),
so that

T M -= dρ -dρ + ∆tT f,n D n ∆t D n + O(∆t) = T f,n + O(∆t).
Remark 6. Actually, the correction (34) is implicit since D n , given by [START_REF] Lee | Gyrokinetic particle simulation model[END_REF], also depends on T M,n through M.

We propose to approximate D n by D n (which will be discussed below) in [START_REF] Mouton | Two-scale semi-lagrangian simulation of a charged particle beam in a periodic focusing channel[END_REF] ; even if the energy will not be preserved exactly with this approximation, the correction (34) will become explicit and will lead to an improved conservation. Now, let discuss how D n is computed.

From bracket properties, we have

(S µ M , S µ M ) (f Np (t n )) = d dt S µ M [f Np ] |t=t n ; then, we approximate the time derivative by (S M [f n+δt Np ] -S M [f n Np ]
)/δt where δt is a small time step and where f n+δt Np is obtained by advancing the particles from t n to t n + δt. This step requires the knowledge of T M (that we are looking for), hence we take T M = T M,n-1 during this step (and T M = T f at t = 0). Then, we set

D n = 1 T M,n-1 (S M [f n+δt Np ] -S M [f n Np ]
)/δt as an approximation of D n in (34) so that T M,n can be computed explicitly. Of course, solving the nonlinear relation [START_REF] Mouton | Two-scale semi-lagrangian simulation of a charged particle beam in a periodic focusing channel[END_REF] to get a more accurate approximation of T M,n is possible at the price of a fixed point procedure which turns out to be very costly. The previous procedure actually corresponds to one iteration of a fixed point procedure which already leads to a good energy preservation. Close to the equilibrium, definition of correction T M,n becomes numerically stiff (since D becomes small). In this regime, we replace correction (34) by the semi-continuous one [START_REF] Herda | Anisotropic boltzmann-gibbs dynamics of strongly magnetized vlasov-fokkerplanck equations[END_REF] since the extra term in (32) vanishes.

Clustering method

Even if the particles approach introduced in [8] is very attractive, a direct numerical implementation is quite expensive since its complexity is in O(N 2 p ) with N p the number of particles. More precisely, for each particule p, one has to compute ∇ v ψ µ * log (f * ψ µ ) (v p ): the first convolution is computed using a d-dimensional Gauss-Hermite quadrature with N q points and f * ψ µ is computed with a sum over the N p particles. Finally, the complexity is O(N d q N 2 p ) where d is the phase space dimension. Some attempts to accelerate have been proposed using treecode in [START_REF] Carrillo | A particle method for the homogeneous landau equation[END_REF] or GPU programming in [START_REF] Zonta | Multispecies structure-preserving particle discretization of the landau collision operator[END_REF]. Here, to reduce the complexity, we propose the following clusterization technique. The idea is to consider only a part of the set of particles to reconstruct f * ψ µ in the logarithm convolution part. Indeed, let introduce N p < N p . We then split the set of N p particles into M clusters of N p particles (N p = M N p ). Let call C i the i-th cluster for i = 1, . . . , M and the corresponding regularized unknown in each cluster writes

f i Np * ψ µ (v) = ρ N p q∈Ci ψ µ (v -v q ).
Thanks to this definition, we replace (29) by the following time evolution of the particles velocities

vp (t) = -T ∇ v ψ µ * log f i Np * ψ µ (v p ) -(v p (t) -u(t)) , ∀p ∈ C i .
For this equation, the complexity becomes O(N d q N p N p ). Moreover, each cluster can be treated separately which is well adapted for parallel calculation. This clustering approach does not destroy the conservation of mass and momentum. It is obvious for the mass since ρ is just a constant number in the scheme. For momentum, denote ρu i the momentum of cluster C i (momentum computed with f i Np ). As u i ̸ = u a priori, this momentum is not conserved and we have

∂ t (ρu i ) = ρ(u i -u). (35)
Equation ( 35) can be deduced from proof of Theorem 3 presented in Appendix B.2 by considering a density with mean velocity u i instead of u. Moreover we have,

1 M M i=1 u i = 1 M M i=1 1 N p q∈Ci v q = 1 N p Np q=1 v q = u.
Thus, we get

∂ t ρu = ∂ t ρ N p Np q=1 v q = 1 M M i=1 ∂ t ρ N p q∈Ci v q = 1 M M i=1 ρ(u i -u) = 1 M M i=1 ρu i -ρu = 0.
As for the classical method presented in previous subsection, the energy is not conserved but correction of T M similar as ( 27) can be derived to conserved energy.

Discretization for full model

The aim of this part is to deal with the discretization of the non homogeneous model (1) involving high oscillations in time and the collisional part. After a change of variable involving the main oscillation, we obtained the filtered formulation ( 5) where Q filt [F ] is given by (6). To apply the particle method on (6), we introduce similarly as in the previous section the regularized version of the filtered Fokker-Planck operator (with

Z = (X, V ) ∈ R 2 ) Q filt,µ [F ](t, s, Z) = ∇ Z • (T A(s)F ∇ Z ψ µ * log(F * ψ µ )(Z)) + ∇ Z • (A(s) (Z -a(s)u) F ) . (36)
The choice of particles approximation is motivated by the fact that the filtering operation commutes with the regularization step. This will be detailed in subsection 4.1. Then, Subsection 2 is devoted to the derivation of the time discretization which is uniformly accurate with respect to ε, using techniques from [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Crouseilles | Uniformly accurate Particle-In-Cell method for the long time two-dimensional Vlasov-Poisson equation with strong magnetic field[END_REF].

Particle approximation for the non homogeneous model

First, as the space and velocity variables are mixed in the filtered formulation, the particle unknown has also to be smoothed in space with a mollifier so that we consider in the sequel ψ µ : w ∈ R 2 → ψ µ (w) ∈ R given by ( 21) with d = 2. Secondly, we notice that the filtration preserves the regularization. In other words, the regularization and the filtration steps commute so that we can apply equivalently a particle method on the regularized version of (1) or on the filtered regularized formulation [START_REF] Zonta | Multispecies structure-preserving particle discretization of the landau collision operator[END_REF]. This is explained in the Proposition below.

Proposition 6. Consider f and F two densities such that F (Z) = f (e -t ε J Z). Then we have the following equalities (with z = (x, v) ∈ R 2 )

ψ µ * log (f * ψ µ ) (z) = ψ µ * log (F * ψ µ ) (Z), Q µ [f ](t, z) = Q filt,µ [F ](t, t/ε, Z).
Proof. The proof relies on the change of variable ξ = e t ε J ζ which preserves the volume. Indeed, from the first convolution, we have

f * ψ µ (w) = f (w -ζ)ψ µ (ζ) dζ = f (e -t ε J W -ζ)ψ µ (ζ) dζ = f (e -t ε J (W -ξ))ψ µ (ξ) dξ = F (W -ξ)ψ µ (ξ) dξ = F * ψ µ (W ).
Consequently, the double convolution becomes

ψ µ * log(f * ψ µ )(z) = ψ µ (z -w) log (f * ψ µ (w)) dw = ψ µ (e -t ε J Z -w) log (f * ψ µ (w)) dw = ψ µ (Z -W ) log (F * ψ µ (W )) dW = ψ µ * log(F * ψ µ )(Z).
The second equality is directly deduced from the first one.

Consider the particle density f Np and the regularized particle density F Np defined by

f Np (t, z) = ϱ N p Np p=1 δ(z -z p (t)) and F Np (t, Z) = ϱ N p Np p=1 δ(Z -Z p (t)), (37)
where the filtered particles are given by Z p (t) = e t ε J z p (t) and ϱ is the total mass of the densities. The ODE system satisfied by z p (t) = (x p (t), v p (t)) obtained by inserting (37) in [START_REF] Ben Abdallah | Diffusion and guiding center approximation for particle transport in strong magnetic fields[END_REF] where Q is replaced by its regularization version [START_REF] Frénod | Long time simulation of a beam in a periodic focusing channel via a two-scale pic-method[END_REF] gives

żp (t) = - 1 ε Jz p (t) -U µ [f Np ](t, v p (t)).
Since filtration and regularization steps commute, performing either the filtering procedure on this latter system or inserting (37) into (36) both gives

Żp (t) = -U filt µ [F Np ](t, t/ε, Z p (t)), (38)
where

U filt µ [F Np ](t, s, Z p (t)) = T A(s)∇ Z ψ µ * log F Np * ψ µ (Z p ) + A(s)(Z p -a(s)u). Remark 7.
From the relation (26) of Proposition 3, kinetic energy is dissipated locally in space which implies that quadratic gyromoments (i.e. those associated to X 2 , XV and V 2 ) are not conserved. It is possible to use the correction [START_REF] Herda | Anisotropic boltzmann-gibbs dynamics of strongly magnetized vlasov-fokkerplanck equations[END_REF] to preserve the quadratic gyromoments. Other gyromoments (i.e. those associated to 1, X, V ) are conserved up to the quadrature error.

Time integration

We are interested in this part by the time approximation of (38). To do so, we first observe that it is a multiscale problem where the fast variable t/ε is periodic of period 2π, which enables to use the double scale strategy introduce in [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Crouseilles | Uniformly accurate Particle-In-Cell method for the long time two-dimensional Vlasov-Poisson equation with strong magnetic field[END_REF]. The main idea is to consider the fast scale t/ε as an independent variable s so that the original solution is recovered by an evaluation on the diagonal s = t/ε. By an abuse of notation, we still denote Z p (t, s) the double scale unknown which satisfies the following set of PDEs

∂ t Z p + 1 ε ∂ s Z p = -U filt µ [F Np ](t, s, Z p ).
We set S p (t, s) = -U filt µ [F Np ](t, s, Z p ). Since Z p is periodic with respect to s, Fourier expansion is used in this direction: Z p (t, s) = ℓ e iℓs Z p,ℓ (t). We then get the following ODE on the Fourier modes Z p,ℓ ( S p,ℓ (t) denotes the Fourier modes associated to S p (t, s))

∂ t Z p,ℓ (t) + iℓ ε Z p,ℓ (t) = S p,ℓ (t), ℓ ∈ Z ⋆ ,
and ∂ t Z p,0 (t) = S p,0 (t) (which actually corresponds to the averaged model ( 7)). Now we introduce the discrete unknown Z n p,ℓ ≈ Z p,ℓ (t n ), t n = n∆t, with ∆t > 0 the time step. Applying the Duhamel formula combined with exponential integrators [START_REF] Hochbruck | Exponential integrators[END_REF] finally gives a first order scheme

Z n+1 p,ℓ = e -iℓ ε ∆t Z n p,ℓ + ∆tφ 1 - iℓ ε ∆t S n p,ℓ , (39)
whereas a second order scheme is given by

Z n+1 p,ℓ = e -iℓ ε ∆t Z n p,ℓ + ∆tφ 1 - iℓ ε ∆t S n p,ℓ + ∆tφ 2 - iℓ ε ∆t S n p,ℓ -S n-1 p,ℓ , (40)
where the φ 1 , φ 2 functions are given by

φ 1 (z) = 1 z (e z -1) , φ 2 (z) = 1 z 2 (e z -1 -z) ,
where S n p,l is the Fourier transform of -U filt µ [F n Np ](t n , s, Z n p ). To compute -U filt µ [F n Np ](t n , s, Z n p ), the moments u and T need to be computed at time t n . To do so, we unfilter F n Np and approximate them using second order spline interpolation on a sufficiently large space grid [x min , x max ] with N x points. Let remark that for the second order scheme (40), the first iteration can be done with the first order scheme (39). Finally, the original particles are reconstructed following z n+1 p = e -t n+1 ε J Z n+1 p (s = t n+1 /ε), Z n+1 p (s) is reconstructed using a Fourier truncated series. For more details, we refer to [START_REF] Chartier | Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations[END_REF][START_REF] Crestetto | Multiscale numerical schemes for the collisional vlasov equation in the finite larmor radius approximation regime[END_REF][START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Crouseilles | Uniformly accurate Particle-In-Cell method for the long time two-dimensional Vlasov-Poisson equation with strong magnetic field[END_REF].

5 Numerical results for the particle method applied to the Fokker-Planck operator

In this section, we focus on the homogeneous model to discuss the behavior and properties of the particle numerical method applied to the Fokker-Planck operator. Section 5.1 considers a relaxation test whereas Section 5.2 highlights the noise that the method can produce according to the time step ∆t and the regularization parameter µ. Section 5.3 is devoted to the conservation properties of the scheme and the temperature correction method. Finally, Section 5.4 illustrates the clustering method. In this section, we consider the following initial density

f 0 (v) = 1 2 √ 2π e -(v-u 1 ) 2 2 + e -(v-u 2 ) 2 2 , (41)
with u 1 = -2.6, u 2 = 3.4.

Relaxation to equilibrium test

In this part, we study a relaxation test starting initially from the bi-Maxwellian (41) which is known to converge towards its Maxwellian equilibrium. The following numerical parameters are chosen: ∆t = 10 -3 , µ = 0.02, the number of quadrature points is N q = 50 and the number of particles is N p = 250 or 10 3 . We initialize N p particles by a random draw of the initial condition f 0 , using a rejection sampling method. Let us remark that the so-obtained particle density f Np (t = 0) has not the same moments of f 0 (except mass). To determine the Maxwellian equilibrium, we compute the moments u and T of f Np (t = 0). Figure 1 shows the results of relaxation test for different times using N p = 250 (left column) and N p = 10 3 (right column). The red plot corresponds to the Maxwellian equilibrium M, the blue dots (on the horizontal axis) represents the particles v p whereas the blue line represents the regularized particles function f Np * ψ µ . At t = 0, the dashed blue line represents the initial condition f 0 . We can see at t = 0 the noise due to the initial sampling on the regularized particle function (which is less important when N p = 10 3 ) but after very short time, the solution becomes smoother. The plot corresponding to the final time t = 1 shows that f Np * ψ has almost reached the Maxwellian equilibrium. To quantify this convergence, we plot in Figure 2 the difference in L 1 , L 2 and L ∞ norms between f Np * ψ µ and M as a function of time in semi-log scale. We observe the (exponential) convergence for small time but a plateau is reached (the level of which is smaller using more particles). The explication is the following: in low density regions, the regularization of isolated particles makes appear the gaussian ψ µ centered at the (isolated) particle velocity. We thus observe a bump around the particle as illustrated in Figure 3. These bumps can not disappear, however this effect can be decreased using more particles N p = 10 3 .

Noise of the method

In this part, we investigate the noise effect according to the numerical parameters ∆t and µ. To do so, we consider the same test as in Subsection 5.1 (the initial condition is given by (41)). In Figure 4, the regularized solution f Np * ψ µ is plotted at t = 1.5 with N p = 250 (µ = 0.02 and N q = 50) and for different time steps: ∆t = 10 -2 , 1.2 × 10 -2 , 1.5 × 10 -2 , 2 × 10 -2 . One observes that when the time step is too large, the level of noise becomes large, which does not allow for a good representation of the solution. In particular, the presence of noise will deteriorate the quality of the Gauss-Hermite quadrature in the convolution terms, but also the conservation properties. Finally, in Figure 5, the time history of a given particles v p is plotted for different time steps. As discussed above, when the time step is too large, the particle does not converge to a right quantity and presents a chaotic behavior.

Conserved properties

Here, we investigate the preservation of the moments according to the numerical parameters ∆t and µ (the other parameters are N p = 250, N q = 50 and the final time is t = 0.5). We focus on the momentum P = ρu and the energy E preservation. To do so, we give in Tables 1, 2 and 3 the following quantities

P err = max n=0,...,N P n -P 0 , E err = max n=0,...,N E n -E 0 ,
for ∆t = 10 -3 , 10 -4 and µ = 10 -1 , 10 -2 for the three different variants of the method presented in Section 3 (without correction, with correction [START_REF] Herda | Anisotropic boltzmann-gibbs dynamics of strongly magnetized vlasov-fokkerplanck equations[END_REF] and with correction (34)). Since the initialisation procedure (rejection sampling) used to sample the particles according to the initial condition is an important source of noise in the method, we applied a minimization algorithm on the sampled particles before the time loop, in order to reduce the Gauss-Hermite quadrature error at the first time steps. One example to reduce the initial noise is to minimize the functional f Np * ψ µ -f 0 2 L 2 on particles. According to the variant, the amplitude of the energy error E err depends on ∆t and µ. When no correction for T M is applied (see Table 1), the momentum is preserved (up to some error due to the Gauss-Hermite quadrature) but not the energy, as expected from Theorem 3.3. In Table 2, the error from the correction proposed for the regularized Fokker-Planck operator is shown. As expected, the energy is not preserved, but one can see an important improvement compared to Table 1, and in particular the error is divided by 10 when the time step is divided by 10. Finally, in Table 3, the error from the correction proposed in (34) is displayed. The correction is able to preserve the error up to about 10 -7 (almost independently of the time step). The reason why machine accuracy is not reached relies on the approximation of D (see Remark 6) in addition with Gauss-Hermite quadrature. Table 3: Conservation of momentum and energy with the correction (34) of T M , using Remark 6.

Behavior of clustering method

The aim of this part is to give some preliminary results obtained with the clustering method presented in Section 5.4. On Figure 6, we plot the solution for t = 0.5 (µ = 0.04, ∆t = 10 -3 ) in different configurations: Figure 6a displays the solution of the classical method with N p = 60, whereas Figure 6b displays the classical method with N p = 1200, Figure 6c corresponds to the clustering method with N p = 1200 split in 20 clusters of Ñp = 60 particles (each cluster is kept unchanged during the simulation) and Figure 6d corresponds to the clustering method with N p = 1200 split in 20 clusters of Ñp = 60 particles (each cluster is randomly chosen at each iteration). As expected, the runtime observed for cases 6c and 6d is shorter than 6b even if there is the same number of particles (the computation of advection field is almost 20 times faster, as expected with complexity in

(a) ∆t = 10 -2 (b) ∆t = 1.2 × 10 -2 (c) ∆t = 1.5 × 10 -2 (d) ∆t = 2 × 10 -2
O(N 2 q N p N p ) instead of O(N 2 q N 2 p )).
Even if the quality of the solution is not as good as the one obtained with the classical method using N p = 1200 particles, in particular in the tail of the solution where there are few particles, the clustering method may be interesting for configurations where a lot of particles are required (for the transport part for example) and where the effect of the collision part is not too important ; this approach is able to strongly decrease the computational cost.

Numerical results for the Fokker-Planck equation with oscillations

In this section, we focus on the model (1) and on the numerical method presented in Section 4, which couples a Fokker-Planck operator in velocity and oscillations in phase space. The section is organized as follows: Subsection 6.1 investigates numerically the differences and links between the filtered-averaged Maxwellian M filt and the gyromaxwellian G. Subsection 6.2 presents some numerical results obtained with (40) to study the relaxation towards the gyromaxwellian equilibrium. In Subsection 6.3, some uniform accuracy results are provided and Subsection 6.4 presents a nonlinear charged particles beam test.

Difference between filtered-averaged Maxwellian and gyromaxwellian

The goal of this section is to illustrate numerically the differences and links between the gyromaxwellian G[F ] defined by [START_REF] Chartier | Numerical methods for the two-dimensional Vlasov-Poisson equation in the finite Larmor radius approximation regime[END_REF] and the filtered-averaged Maxwellian M filt [F ] defined by [START_REF] Chartier | Uniformly accurate methods for Vlasov equations with non-homogeneous strong magnetic field[END_REF], both being associated to a density F . Moreover, we illustrate numerically the fact that G[F ] is a fixed point of the operator

M filt : F → M filt [F ].
First, we recall that computing G[F ] from [START_REF] Chartier | Numerical methods for the two-dimensional Vlasov-Poisson equation in the finite Larmor radius approximation regime[END_REF] requires to compute the gyromoments (11) of F . Second, computing M filt [F ] requires to compute M filt [F ](s) for all s ∈ [0, 2π] and then to take the average over s as in [START_REF] Chartier | Uniformly accurate methods for Vlasov equations with non-homogeneous strong magnetic field[END_REF]. To compute M filt [F ](s), we first unfilter F by setting f (z) = F (e sJ z), so that we can compute the moments (2) of f , which determine the Maxwellian M given by (4). Finally, M filt [F ](s) given by ( 9) is nothing but the Maxwellian M expressed in the filtered variables Z = e -sJ z.

Now, for all m ∈ N ⋆ , we recursively define

M filt •m [F ] by M filt •m [F ] = M filt •(m-1) M filt [F ] .
Below, we will illustrate the fact that M filt •m [F ] converges towards G[F ] when m goes to infinity, and estimate the associated rate of convergence. We start by describing a procedure to compute G[F ] and

M filt •m [F ] for a given density F (Z), Z = (X, V ) ∈ R 2 .
We consider the following function (which is bi-Maxwellian in this example)

F (X, V ) = 1 2π T X1 T V1 e -| X-u X 1 | 2 2T X 1 e -| V -u V 1 | 2 2T V 1 + 1 2π T X2 T V2 e -| X-u X 2 | 2 2T X 2 e -| V -u V 2 | 2 2T V 2
, with u X1 = -3, u V1 = -2.5, u X2 = 2.5, u V2 = 3 and T X1 = 1.75, T V1 = 1.25, T X2 = 1.5, T V2 = 1.25. On the one side, as discussed before, the gyromaxwellian G[F ] (given by ( 12)) requires the gyromoments of F defined by F (X, V )C(X, V )dXdV with C(X, V ) = (1, X, V, X 2 , XV, V 2 ) T . On the other side, the filtered-averaged Maxwellian is defined and approximated as follows:

M filt [F ](X, V ) = 1 2π 2π 0 M filt (s)[F ](X, V ) ds ≈ 1 N s Ns-1 j=0 M filt (s j )[F ](X, V )∆s, ( 42
)
where the integral is approximated using N s quadrature points s j = j∆s, j = 0, ..., N s -1 and ∆s = 2π/N s . Now, it remains to compute the filtered-averaged Maxwellian M filt (s j )[F ](X, V ). To do so, we first apply the inverse change of variable on F for each s j (j = 0, . . . , N s -1) to get f j (z) = F (e sj J z), z = (x, v) and then compute the moments f j (x, v)(1, v, v 2 ) T dv of f j . Once these moments are computed, the Maxwellian M j (x, v) can be constructed for j = 0, . . . , N s -1 following (4) and finally, we apply the change of variable M filt (s j )[F ](Z) = M j (e -sj J Z), j = 0, . . . , N s -1.

We then compute the average using (42). We repeat this procedure iteratively to compute

M filt •m [F ] for m ≥ 1.
In Figure 7, we plot the level set of F, M filt •m [F ] (for m = 1, 2, 3, 4, 8, 12) and G[F ] in the filtered phase space (X, V ). We can see in particular that M filt •m [F ] is far from the gyromaxwellian shape G[F ] for small m. However, when m increases, the shape of M filt

Relaxation towards a gyromaxwellian

The goal of this part is to illustrate the capability of the second order scheme (40) to capture the relaxation of the solution F of the averaged model [START_REF] Carrillo | A blob method for diffusion[END_REF] towards the gyromaxwellian equilibrium G[F ]. To do so, instead of solving the averaged model [START_REF] Carrillo | A blob method for diffusion[END_REF], we solve the filtered model ( 5) with a small ε (with a bi-Maxwellian initial condition). Indeed, since the error between the averaged model ( 7) and the filtered model ( 5) is O(ε) (and then negligible) and since our scheme is uniformly accurate in ε, the error is O(∆t 2 ) + O(ε). The numerical parameters are chosen as follows: ∆t = 0.0133, ε = 10 -10 , µ = 0.075, N q = 20, N p = 200, N s = 32. On Figure 9, we plot the particles used to approximate the filtered unknown F for t = 0, t = 0.4933, t = 1 and t = 4. We also plot the level set of the gyromaxwellian G[F ] associated to the initial condition F (t = 0). We can observe that the particles enjoy some diffusion in both velocity and space directions as expected. Moreover, the particles spread out with time to reach the shape of the gyromaxwellian G[F ] for large time. More quantitative results would be required but there are quite difficult to obtain since a large number of particles would be required to decrease the noise. However, the method enables us to illustrate the relaxation phenomena even with a low number of particles.

Order of the method

In this part, we show that the scheme (40), solving equation (38), is of the second order and uniformly accurate in ε for sufficiently enough Fourier modes. To reduce the numerical cost we consider N p = 10 particles. The reference solution is computed using the same solver but with a smaller time step. We take µ = 1 and N q = 20 whereas the space mesh (used to compute the moments) is defined by x min = -10, x max = 10 and N x = 100. We consider the final time t = 1. We plot the order for different ε = 10 -1 , . . . , 10 -8 and for different N s on Figure 10a. We can observe that for N s = 64 Fourier modes, the error does not decrease with a second order for large ε (second order is well recovered for small ε). However, if we increase the number of Fourier modes to N s = 1024, second order accuracy is observed uniformly in ε (see Figure (10b)). These results may indicate the solution is not smooth enough, which may prevent the spectral convergence of the Fourier method used in the s direction. Similar behavior was observed in [START_REF] Crouseilles | Uniformly accurate forward semi-Lagrangian methods for highly oscillatory Vlasov-Poisson equations[END_REF]. This lack of smoothness can be imputed to the projection/interpolation during the computation of moment's step. It mixes a nonlinearity (when we compute u and T from moments) with low regularity function space (splines space). To illustrate this, we consider the linear Fokker-Planck operator with u = 0 and T = 1 so that we do not compute the moments anymore. On Figures 10c and10d, we can observe second order accuracy for all ε, even for N s = 64.

A beam of particles with collisions

We consider a particle beam test as studied in [START_REF] Crouseilles | Uniformly accurate forward semi-Lagrangian methods for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Mouton | Two-scale semi-lagrangian simulation of a charged particle beam in a periodic focusing channel[END_REF], but here, a collisional operator is added to investigate its interaction with the transport part. Here, we consider a self-consistent electric field E f (computed from a Poisson The radial variable r > 0 plays the same role as x. By symmetry, we can consider r ∈ R imposing that E f (-r) = -E f (r) and f (-r, -v) = f (r, v). The aim of this test is to investigate the effect of the collision operator on the thin structures (as filamentation) created by the nonlinear transport part (due to the presence of the self-consistent electric field E f ). From a numerical point of view, the scheme (40) can be easily extended to this nonlinear case and will be compared to a reference solution obtained with a direct method (spectral in r and v combined with a time splitting method and refined numerical parameters) to solve (43). We impose T = 1 and u = 0 in the Fokker-Planck operator in order to avoid the moment computation step and consider N s = 16. Finally, we will consider two values for ν to illustrate the effect of the Fokker-Planck operator (ν = 0 and ν = 0.1) on the filamentation. For this test, we consider the following initial condition

(t, r, v) (t ≥ 0, r, v ∈ R) is ∂ t f - 1 ε Jz • ∇ z f + E f ∂ v f = νQ[f ], (43)
f 0 (x, v) = n 0 2πv 2 th e -1 2 (v/v th ) 2 I [-rmax,rmax] (x),
where r max = 1.85, v th = 0.1, n 0 = 12 (ρ = 2r max n 0 ) and I is the indicator function. The numerical parameters are chosen as follows: N p = 1600, ε = 0.001, µ = 0.005, N q = 20 and ∆t = 0.02666. In Figure 11, the solutions f (t = 2, x, v) obtained by our particle scheme and by the spectral scheme are plotted with ν = 0 and ν = 0.1. First, in the collisionless regime, the results obtained by the particle method and the reference method are very close and also similar to the ones obtained in the litterature [START_REF] Crouseilles | Uniformly accurate forward semi-Lagrangian methods for highly oscillatory Vlasov-Poisson equations[END_REF][START_REF] Mouton | Two-scale semi-lagrangian simulation of a charged particle beam in a periodic focusing channel[END_REF]. Indeed, thin filaments are created which are well reproduced by the particle method. Second, when ν = 0.1, the filaments are dissipated and an equilibrium is reached. Again, the results obtained by the two methods are in very good agreement. Regarding (11b), we have, still using the change of variable Z ← (Z -U ) and (45)

R 2 ZG[F ](Z) dZ = ϱ det(2πT) R 2 Z exp - 1 2 Z T T -1 Z dZ + ϱU = ϱ det(2πT) R e -1 2T XX X 2 R Ze -T XX 2 det(T) V - T XV T XX X 2 dV dX + ϱU.
For Z = X, the integral above is zero by oddness in X. For Z = V , the integral above in V is proportional to (T XV /T XX )X which also gives an odd function in X, so that the integral is also zero in this case. Finally, for (11c), we first have

R 2 (Z -U ) ⊗ (Z -U )G[F ](Z) dZ = ϱ det(2πT) R 2 Z ⊗ Z exp - 1 2 Z T T -1 Z dZ. ( 46
)
Let focus on the diagonal terms. For the X 2 term, we have

ϱ det(2πT) R 2 X 2 exp - 1 2 Z T T -1 Z dZ = ϱ det(2πT) R X 2 e -1 2T XX X 2 R e -T XX 2 det(T) V - T XV T XX X 2 dV dX = ϱ √ T XX R X 2 e -1 2T XX X 2 dX = ϱT XX , (47)
and symmetrically, for the V 2 term, we obtain

R 2 V 2 G[F ](Z) dZ = ϱT V V .
For the extra-diagonal terms, we have ϱ

det(2πT) R 2 XV exp - 1 2 Z T T -1 Z dZ = ϱ det(2πT) R Xe -1 2T XX X 2 R V e -T XX 2 det(T) V - T XV T XX X 2 dV dX = ϱ T XV T XX R X 2 e -1 2T XX X 2 dX = ϱT XV ,
which enables to get (11c). Let remark that using the above calculations, we also have

R 2 Z ⊗ ZG[F ](Z) dZ = ϱ (U ⊗ U + T).
We conclude the proof by showing that T is positive definite. As T is symmetric, it is diagonalizable in R. Using Cauchy-Schwarz inequality, we get (X -U X )(V -U V )F dX dV 2 < (X -U X ) 2 F dX dV (V -U V ) 2 F dX dV, and so det(T) > 0 (the equality case is not reached because (X-U X ) √ F and (V -U V ) √ F are linearly independant). Moreover, we have tr(T) ≥ 0. We deduce that the two eigenvalues are positive and thus, T is symmetric positive definite.

B Conservation properties of Fokker-Planck and regularized Fokker-Planck equations

B.1 Proof of Theorem 2

We prove here Theorem 2. Even if the proof is classical, we use here the properties of bracket [START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF] and the reformulated equation ( 18) satisfied by any functional F of f , in particular the moments and the entropy. For the entropy inequality, as the bracket is negative, we deduce it immediately from d dt S M = (S M , S M ) ≤ 0. For the mass conservation, we have δρ δf = 1 and the mass is conserved from (ρ, S M ) = 0. For momentum, we have δ(ρui) δf = v i . It follows that has to be zero so that

(ρu i , S M ) = - R d T f (∇ v v i • ∇ v log(f /M)) dv = - R d T f ∂ vi log(f /M) dv = - R d T ∂ vi f + f (v i -u i ) dv = - R d f v i dv + R d f dv u i = 0.
R d T f ∇ v |v| 2 • ∇ v log(f /M) dv = - d i=1 R d T v i f ∂ vi log(f /M) dv = - d i=1 R d T v i ∂ vi f + f v i (v i -u i ) dv =
∇ v δS M δf = ∇ v f f + v -u T = 0.
The solution are f (v) = Ce -|v-u| 2

2T

with C ∈ R determined by integrating on v ∈ R d .

B.2 Proof of Theorem 3

We prove here Theorem 3 presented in Section 3.2, establishing some properties of the regularized Fokker-Planck collision operator. To prove the Theorem, we need the following technical and classical Lemma. Proof of Theorem 3. The mass conservation in [START_REF] Frénod | Long time behavior of the two-dimensional vlasov equation with a strong external magnetic field[END_REF] is obtained similarly as for the original Fokker-Planck operator. For the momentum conservation [START_REF] Frénod | Long time behavior of the two-dimensional vlasov equation with a strong external magnetic field[END_REF], we replace f by f v i dv (i = 1, . . . , d) in the the bracket form [START_REF] Frénod | Long time simulation of a beam in a periodic focusing channel via a two-scale pic-method[END_REF] to get

∂ t (ρu i ) = R d f v i dv, S µ M = -T M R d f ∂ vi (ψ µ * log (f * ψ µ /M)) dv = - R d T M f ∂ vi (ψ µ * log (f * ψ µ )) dv + R d f (v i -u i ) dv = R d T M f ψ µ * ∂ vi f * ψ µ f * ψ µ dv = R d T M (f * ψ µ ) ∂ vi f * ψ µ f * ψ µ dv = R d T M ∂ vi f * ψ µ dv = 0.
Considering now the energy relation, we proceed as follows (replacing f by f |v| 2 /2 dv in ( 24)) to get

1 2 ∂ t ρ(|u| 2 + dT f ) = R d f |v| 2 2 dv, S µ M = -T M d i=1 R d f v i ∂ vi (ψ µ * log(f * ψ µ /M)) dv = -T M d i=1 R d f v i ∂ vi (ψ µ * log(f * ψ µ )) dv - d i=1 R d f v i (ψ µ * (v i -u i )) dv = -T M d i=1 R d f v i ∂ vi (ψ µ * log(f * ψ µ )) dv - d i=1 R d f v i (v i -u i ) dv. ( 48
)
As for the original Fokker-Planck operator, the second term is equal to -dT f ρ after an integration by part.

Regarding the first term, we have

d i=1 R d f v i ∂ vi (ψ µ * log(f * ψ µ )) dv = d i=1 R d f v i ψ µ * ∂ vi f * ψ µ f * ψ µ dv = d i=1 R d ((f v i ) * ψ µ ) ∂ vi f * ψ µ f * ψ µ dv = d i=1 R d v i (f * ψ µ ) ∂ vi f * ψ µ f * ψ µ dv - d i=1 R d (f * (v i ψ µ )) ∂ vi f * ψ µ f * ψ µ dv = d i=1 R d v i ∂ vi f * ψ µ dv + d i=1 R d (f * (µ∂ vi ψ µ )) ∂ vi f * ψ µ f * ψ µ dv = -dρ + µ R d |∇ v f * ψ µ | 2 f * ψ µ dv = -dρ + µ R d (f * ψ µ ) |∇ v log (f * ψ µ )| 2 dv = -dρ - µ T M (S, S) (f * ψ µ ) = -dρ -µD(f ),
where we used the definition of the bracket [START_REF] Crouseilles | Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations[END_REF], the entropy [START_REF] Crouseilles | Uniformly accurate Particle-In-Cell method for the long time two-dimensional Vlasov-Poisson equation with strong magnetic field[END_REF] and the definition of D(f ) in the last line.

Gathering the two terms in (48) finally gives [START_REF] Herda | Large-time behavior of solutions to vlasov-poisson-fokker-planck equations: From evanescent collisions to diffusive limit[END_REF]. We end the proof of Theorem 3 by looking for the equilibrium of the regularized operator. Following the proof of the original Fokker-Planck operator, we have to solve ∇ v ψ µ * log(f * ψ µ /M) = 0, which is equivalent to

∇ v ψ µ * log(f * ψ µ ) = - v -u T M .
Integrating with respect to v gives

T M ψ µ * log(f * ψ µ ) = - |v| 2 2 + u • v + C, (49)
which is more complicated to solve than in the original case. However, using Fourier techniques and following the lines in [START_REF] Carrillo | A particle method for the homogeneous landau equation[END_REF], it is possible to compute the solution f of (49). Let denote F [f ] the Fourier transform in v of f

F [f ] (ξ) = R d f (v)e -iv•ξ dv.
Using the properties of the Fourier transform, we can isolate f * ψ µ in (49)

f * ψ µ = exp F -1 F 1 T M - |v| 2 2 + u • v + C /F [ψ µ ] .
We compute now the right hand side step by step. First, applying F to (49) gives

T M F [log(f * ψ µ )] = 1 F [ψ µ ] F - |v| 2 2 + u • v + C = (2πµ) d/2 ψ 1/µ ∆ 2 + iu • ∇ + C δ ξ=0 .

Figure 1 :

 1 Figure 1: f Np * ψ µ as a function of v for different time steps: relaxation of the bi-Maxwellian initial condition towards the Maxwellian. Left: N p = 250. Right: N p = 10 3 . From top to bottom, t = 0, t = 0.15, t = 0.35, t = 1. ∆t = 10 -3 , µ = 0.02, and N q = 50.

Figure 2 :

 2 Figure 2: Time evolution of ∥f Np * ψ µ -M∥ for N p = 250 (left) and N p = 10 3 (right).

Figure 3 :

 3 Figure 3: f Np * ψ µ for v ∈ [6, 10]: illustration of the bumps around isolated particles.

Table 1 : 4 Table 2 :

 142 Conservation of momentum and energy without correction of T M . P err µ = 10 -1 µ = 10 -2 ∆t = 10 -3 8.4644 × 10 -9 5.5386 × 10 -10 ∆t = 10 -4 8.4531 × 10 -9 5.3865 × 10 -10 E err µ = 10 -1 µ = 10 -2 ∆t = 10 -3 2.5359 × 10 -3 2.3464 × 10 -3 ∆t = 10 -4 2.5235 × 10 -4 2.3337 × 10 -Conservation of momentum and energy with the semi-discrete correction (27) of T M . P err µ = 10 -1 µ = 10 -2 ∆t = 10 -3 8.4758 × 10 -9 5.5325 × 10 -10 ∆t = 10 -4 8.4542 × 10 -9 5.3858 × 10 -10 E err µ = 10 -1 µ = 10 -2 ∆t = 10 -3 1.7000 × 10 -7 3.5123 × 10 -7 ∆t = 10 -4 4.9941 × 10 -8 1.7070 × 10 -7

Figure 4 :

 4 Figure 4: f Np * ψ µ as a function of v for different time steps: influence of the time step on the level of noise. N p = 250, µ = 0.02, and N q = 50.

  (a) Classical method for Np = 60. (b) Classical method for Np = 1200. (c) Clustering method for Np = 1200 split in 20 clusters of 60 particles. (d) Random clustering method for Np = 1200 split in 20 clusters of 60 particles.

Figure 6 :

 6 Figure 6: Comparaison between the classical method and the clustering method.

  •m [F ] becomes closer and closer to the one of G[F ]. This is confirmed by the error curves displayed on Figure 8 in which we plot M filt •m [F ] -G[F ] (using different discrete norms L 1 , L 2 and L ∞ ) as a function of m, in semi-log scale. The exponential rate of convergence (with a rate estimated at about 0.25, dashed line) clearly appears. An interpretation of these results would be that the gyromaxwellian G[F ] is a fixed point of the averaged filtered Maxwellian operator F → M filt •m [F ], that is M filt [G[F ]] = G[F ]. In terms of operators, we have lim m→∞ M filt •m = G where G denotes the gyromaxwellian operator F → G[F ].

7 :

 7 Level set of the functions F (top left), M filt •m [F ] for different m (from top right to bottom left, m = 1, 2, 3, 4, 8, 12), and G[F ] (bottom right).

Figure 8 :

 8 Figure 8: ∥ M filt •m [F ] -G[F ]∥ as a function of m (semi-log scale).

Figure 9 :

 9 Figure 9: Relaxation of the filtered unknown F (red particles) towards its associated gyromaxwellian G[F ].

  Ns = 64 and u = 0, T = 1 (d) Ns = 1024 and u = 0, T = 1

Figure 10 :

 10 Figure 10: Numerical order for the second order scheme.

Figure 11 :

 11 Figure 11: Snapshots of the solution of the beam problem (43). First line: particle method N p = 1600, ν = 0 (left) and ν = 0.1 (right). Second line: reference (spectral) method ν = 0 (left) and ν = 0.1 (right).

2 +

 2 For energy, from the relation δ(ρ|u| 2 +dρT ) dT ), S M ) = -1 2

dT ρ -ρ |u| 2 +

 2 dT + ρ|u| 2 = 0.Let us focus on the equilibrium. It is clear thatM satisfies Q[M] = 0. Suppose now that f satisfies Q[f ] = 0. Thus we have R d Q[f ] log(f /M) dv = 0 ; using[START_REF] Degond | The weighted particle method for convection-diffusion equations. I. The case of an isotropic viscosity[END_REF] and an integration by part, we get0 = R d Q[f ] log(f /M) dv = -R d T f ∇ v δS M δf 2 dv = (S M , S M ) = d dt S M .It remains to prove that if d dt S M [f ] = 0, then f = M. As f > 0 and T > 0, then ∇ v δS M δf 2

Lemma 1 .

 1 For all functions f, g : R d → R, we have the following relations:f (ψ µ * g) dv = (f * ψ µ ) g dv, f (∆ v ψ µ * g) dv = (f * ∆ v ψ µ ) g dv, f * ψ µ dv = f dv, (f * ψ µ ) v dv = f v dv, (f * ψ µ ) |v| 2 dv = f |v| 2 dv + dµ f dv, v (f * g) dv = ((vf ) * g) dv + (f * (vg)) dv.

Acknowledgments

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200 EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. This work has been also supported by ANR project MUFFIN (ANR-19-CE46-0004) and by Centre Henri Lebesgue, programme ANR-11-LABX-0020-01.

A End of the proof of Theorem 1

In this Appendix, we prove the last points of Theorem 1, which consists in establishing that T defined by (11c) is symmetric positive definite and that

We start by proving (44). To do so, we prove that the gyromoments of G[F ] are equal to the gyromoments of F denoted by (ϱ, ϱU, ϱT ) as in [START_REF] Chartier | Uniformly accurate methods for three dimensional vlasov equations under strong magnetic field with varying direction[END_REF]. First, denoting Z = (X, V ) and using the notations of Theorem 1, we have

To prove (11a), the change of variable Z ← (Z -U ) together with (45) give

Applying the inverse Fourier transform F -1 to this latter equation gives

Taking the exponential on both sides leads to f * ψ µ = C exp -(|v -u| 2 )/2T M where C = exp 1

is determined by integrating the equality

We finally find f * ψ µ = M. If we continue computation to isolate f using the same technic, we find that f is a Maxwellian characterized by ρ, u and T M -µ.