# Deterministic particle method for Fokker-Planck equation with strong oscillations

Anaïs Crestetto, Nicolas Crouseilles, Damien Prel

► To cite this version:

Anaïs Crestetto, Nicolas Crouseilles, Damien Prel. Deterministic particle method for Fokker-Planck equation with strong oscillations. SMAI Journal of Computational Mathematics, 2023, 10, pp.141-173. 10.5802/smai-jcm.109 . hal-04343011

HAL Id: hal-04343011

https://hal.science/hal-04343011v1

Submitted on 13 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Deterministic particle method for Fokker-Planck equation with strong oscillations

Anaïs Crestetto*

Nicolas Crouseilles†

Damien Prel*

December 13, 2023

## Abstract

The aim of this paper is to investigate a deterministic particle method for a model containing a Fokker-Planck collision operator in velocity and strong oscillations (characterized by a small parameter $\epsilon$) induced by a space and velocity transport operator. First, we investigate the properties (collisional invariants and equilibrium) of the asymptotic model obtained when $\epsilon \to 0$. Second a numerical method is developed to approximate the solution of the multiscale Fokker-Planck model. To do so, a deterministic particle method (recently introduced for the Landau equation in [8]) 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 $\epsilon$ that enables to capture numerically the solution of the asymptotic model. Numerous numerical results are displayed, illustrating the behavior of the method.

**Keywords.** Vlasov equation, Fokker-Planck collision operator, highly oscillatory systems, multiscale numerical schemes, Particle method.

**MSC codes.** 65M25, 65M75, 35Q83.

## 1 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) \ge 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 induce 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 [4–6, 16]. In particular, in the limit $\epsilon \to 0$ (where $\epsilon$ 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 [4–6]), 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

<footnote>*Nantes Université, CNRS, Laboratoire de Mathématiques Jean Leray, LMJL, UMR 6629, F-44000 Nantes, France & Inria (MIN-GuS), France</footnote>

<footnote>†Univ Rennes, Inria (MINGuS) & IRMAR UMR 6625 & ENS Rennes, France</footnote>
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 [6, 16] 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 literature is not abundant ; let us quote however [1, 3, 4, 6, 26, 27]).

First, one has to study the averaged model of our beam-collision model. To do so, denoting by $\epsilon$ a parameter
characterizing the period of oscillations ($2\pi\epsilon$) induced by the transport part, a change of variables (usually called
filtering step) is performed to get the filtered system ; then, the (strong) limit $\epsilon \to 0$ can be performed leading to
the averaged model (see [6, 14, 15]). This asymptotic model corresponds to an averaging of the filtered collision
operator on one oscillation period. Using [6] 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 $\epsilon$ 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 [2,35]. 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 [7]. 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 [8, 28, 36]: 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 [20, 21].

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 [8, 28, 36], 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 [30] 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 $\epsilon$) can be constructed. The error of this scheme is
independent of the oscillation period $\epsilon$, and then the scheme is consistent with the asymptotic averaged model
obtained when $\epsilon$ 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 [16] 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 filtered-
averaged 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.

## 2 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 $\epsilon$ 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.

### 2.1 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) \ge 0$ the particles distribution function depending on time $t \in \mathbb{R}_+$, space $x \in \mathbb{R}$ and velocity $v \in \mathbb{R}$. We denote $z = (x, v) \in \mathbb{R}^2$ the phase space variable and introduce $\epsilon > 0$ which denotes a small parameter (typically a ratio between two characteristic lengths). The model we are interested in is the following

$$ \partial_t f - \frac{1}{\epsilon} J z \cdot \nabla_z f = \partial_v ((v-u)f + T \partial_v f) \quad (1) $$

where $J$ denotes the symplectic matrix

$$ J = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix} $$

and the macroscopic quantities $\rho$ (density), $u$ (mean velocity) and $T$ (temperature) are defined from the moments of $f$ by

$$ \rho = \int f \, dv, \quad \rho u = \int v f \, dv, \quad \rho \left( \frac{u^2}{2} + \frac{T}{2} \right) = \int \frac{v^2}{2} f \, dv. \quad (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] = \partial_v ((v-u)f + T \partial_v f) = \partial_v (T f \partial_v \log(f/\mathcal{M})), \quad (3) $$

where $\mathcal{M}$ is the Maxwellian associated to $f$ (through the moments (2))

$$ \mathcal{M}[f](t, x, v) = \frac{\rho}{\sqrt{2\pi T}} e^{-\frac{|v-u|^2}{2T}}. \quad (4) $$

It is interesting to know how the distribution function $f$ behaves as $\epsilon \to 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), \quad \text{where} \quad Z = e^{\frac{t}{\epsilon} J} z, \quad z = (x, v), \quad Z = (X, V), $$

which satisfies the filtered problem

$$ \partial_t F(t, Z) = Q^{\text{filt}}[F](t, t/\epsilon, Z), \quad (5) $$

where $Q^{\text{filt}}[F](t, s, Z)$ denotes the filtered collision operator (obtained by applying the change of variables to (3))

$$ Q^{\text{filt}}[F](t, s, Z) = \nabla_Z \cdot (T \mathcal{A}(s) \nabla_Z F) + \nabla_Z \cdot (\mathcal{A}(s)(Z - a(s)u)F), \quad (6) $$

with

$$ \mathcal{A}(s) = \begin{pmatrix} \sin^2(s) & -\sin(s) \cos(s) \\ -\sin(s) \cos(s) & \cos^2(s) \end{pmatrix}, \quad a(s) = (-\sin(s), \cos(s))^T, \quad b(s) = (\cos(s), \sin(s))^T, $$
such that $x = b(t/\epsilon) \cdot Z$ and $v = a(t/\epsilon) \cdot 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 $\bar{F}$ solution of the averaged model (see [14, 15, 19]) given by

$$ \partial_t \bar{F}(t, Z) = \langle Q^{\text{filt}} \rangle [\bar{F}](t, Z), \quad (7) $$

where $\langle Q^{\text{filt}} \rangle$ is the averaged of $Q^{\text{filt}}$ on the fast variable $s$, which is present in $u, T, A$ and $a$. It is defined by

$$ \langle Q^{\text{filt}} \rangle [\bar{F}](t, Z) = \frac{1}{2\pi} \int_{0}^{2\pi} Q^{\text{filt}}[\bar{F}](t, s, Z) \, ds. \quad (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^{\text{filt}}[F]$ given by (6) can also be expressed in the following divergence form

$$ Q^{\text{filt}}[F](t, s, Z) = \nabla_Z \cdot (T A(s) F \nabla_Z \log(F/M^{\text{filt}}[F](s))), $$

where $M^{\text{filt}}[F]$ denotes the Maxwellian sharing the same three first moments of $f$, but expressed in the filtered variables

$$ M^{\text{filt}}[F](t, s, Z) = \frac{\rho}{\sqrt{2\pi T}} \exp \left( -\frac{1}{2T} |a(s) \cdot Z - u|^2 \right). \quad (9) $$

In the definition of $Q^{\text{filt}}$ and $M^{\text{filt}}$, the macroscopic quantities $\rho, 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) \cdot Z$ (which corresponds to $x$ for $s = t/\epsilon$). For example, we have the following relations: $\rho(x) = \rho(b(s) \cdot Z)$ but also $Q^{\text{filt}}[F](t, t/\epsilon, X, V) = Q[f](t, x, v)$ and $M^{\text{filt}}[F](t, t/\epsilon, X, V) = M[f](t, x, v)$. We also define the filtered-averaged Maxwellian

$$ \langle M^{\text{filt}} \rangle [F](t, Z) = \frac{1}{2\pi} \int_{0}^{2\pi} M^{\text{filt}}[F](t, s, Z) \, ds. \quad (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 $\frac{1}{2\pi} \int_{0}^{2\pi} A(s) \, ds = \frac{1}{2} \text{Id}$ and then $\langle Q^{\text{filt}} \rangle [F] = \frac{1}{2} (T \Delta_Z F + \nabla_Z \cdot (ZF))$.

## 2.2 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$, $\iint 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 [6] for details). In the following, some properties of the filtered-averaged operator $\langle Q^{\text{filt}} \rangle$ are given such as its collisional invariants and equilibrium. Several techniques to compute them have been introduced in [6] 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: \mathbb{R}^2 \to \mathbb{R}^2$ (similarly $F_s(Z) = F(e^{sJ}Z)$ for $F: \mathbb{R}^2 \to \mathbb{R}$).

**Proposition 1.** Let $F: \mathbb{R}^2 \to \mathbb{R}$ be a function, $Q$ be the Fokker-Planck operator (3) and $\langle Q^{\text{filt}} \rangle$ be the filtered Fokker-Planck operator (8). Then we have the following assertions.

1. $F$ is an equilibrium of $\langle Q^{\text{filt}} \rangle$ if and only if $F_s$ is an equilibrium of $Q$, $\forall s \in \mathbb{R}$.

2. $C$ is a collisional invariant of $\langle Q^{\text{filt}} \rangle$ if and only if $C_s$ is a collisional invariant of $Q$, $\forall s \in \mathbb{R}$.

3. If $C$ is a collisional invariant of $\langle Q^{\text{filt}} \rangle$, then $V \partial_X C - X \partial_V C$ is a collisional invariant of $\langle Q^{\text{filt}} \rangle$.

4. $F$ is an equilibrium of $\langle Q^{\text{filt}} \rangle$ if and only if $\log(F)$ is a collisional invariant of $\langle Q^{\text{filt}} \rangle$.
*Proof.* The proof of the two first points is similar to [6] and will not be given. Let consider the third assumption. We first remark that

$$ \frac{d}{ds}C_s(Z) = \frac{d}{ds}C(e^{sJ}Z) = J e^{sJ}Z \cdot (\nabla_Z C)(e^{sJ}Z) = J e^{sJ}Z \cdot (\nabla_Z C)_s(Z) = (JZ \cdot \nabla_Z C)_s(Z) = (-V \partial_X C + X \partial_V C)_s(Z). $$

Assuming $C$ is a collisional invariant of $\langle Q^{\text{filt}} \rangle$, then by the second assertion, it is equivalent to say $C_s$ is a collisional invariant of $Q$ for all $s \in \mathbb{R}$. Thus for all $h > 0$, $\frac{C_s+h-C_s}{h}$ is also a collisional invariant by linearity and considering the limit $h \to 0$ gives $\frac{d}{ds}C_s = (-V \partial_X C + X \partial_V C)_s$ is a collisional invariant of $Q$ for all $s \in \mathbb{R}$. According to the second assertion, this is equivalent to say that $(-V \partial_X C + X \partial_V C)$ is a collisional invariant of $\langle Q^{\text{filt}} \rangle$. We prove now assumption 4. Let $F$ be an equilibrium of $\langle Q^{\text{filt}} \rangle$. From assumption 1, it means that $F_s$ is an equilibrium of $Q$ for all $s \in \mathbb{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 $\langle Q^{\text{filt}} \rangle$. $\square$

The following theorem enables us to derive the invariants and the equilibrium of $\langle Q^{\text{filt}} \rangle$.

**Theorem 1.** The collisional invariants of $\langle Q^{\text{filt}} \rangle$ are linear combinations of $1, X, V, X^2, XV$ and $V^2$. Let the mass $\rho > 0$, the mean velocity $\mathcal{U} \in \mathbb{R}^2$ and the symmetric positive-definite temperature tensor $\mathcal{T} \in \mathbb{R}^{2 \times 2}$, defined from the gyromoments of $F$:

$$ \varrho = \int_{\mathbb{R}^2} F(X, V) \, dX \, dV, \qquad (11a) $$

$$ \varrho \mathcal{U} = \varrho \begin{pmatrix} \mathcal{U}_X \\ \mathcal{U}_V \end{pmatrix} = \int_{\mathbb{R}^2} \begin{pmatrix} X \\ V \end{pmatrix} F(X, V) \, dX \, dV, \qquad (11b) $$

$$ \varrho \mathcal{T} = \varrho \begin{pmatrix} \mathcal{T}_{XX} & \mathcal{T}_{XV} \\ \mathcal{T}_{XV} & \mathcal{T}_{VV} \end{pmatrix} = \int_{\mathbb{R}^2} \begin{pmatrix} X - \mathcal{U}_X \\ V - \mathcal{U}_V \end{pmatrix} \otimes \begin{pmatrix} X - \mathcal{U}_X \\ V - \mathcal{U}_V \end{pmatrix} F(X, V) \, dX \, dV. \qquad (11c) $$

Then the only equilibrium sharing the same gyromoments as $F$ is given by an anisotropic Maxwellian (called gyromaxwellian)

$$ \mathcal{G}[F](Z) = \frac{\varrho}{\sqrt{\det(2\pi\mathcal{T})}} e^{-\frac{1}{2}(Z-\mathcal{U})^T \mathcal{T}^{-1} (Z-\mathcal{U})}, \quad Z = (X, V) \in \mathbb{R}^2, \qquad (12) $$

that is, $\int_{\mathbb{R}^2} C(X, V)(F(X, V) - \mathcal{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 $\langle Q^{\text{filt}} \rangle$. From the second assertion of Proposition 1, we know that $C_s$ is an invariant of $Q$ for all $s \in \mathbb{R}$. In particular for $s = 0$ we have that $C = C_0$ is a collisional invariant of $Q$. So

$$ C(X, V) = \alpha(X)V^2 + \beta(X)V + \gamma(X), \qquad (13) $$

where $\alpha, \beta$ and $\gamma$ are smooth functions. From the third assertion of Proposition 1, $V \partial_X C - X \partial_V C$ is also a collisional invariant of $\langle Q^{\text{filt}} \rangle$, and thus of $Q$. This means that

$$ V \partial_X C(X, V) - X \partial_V C(X, V) = \tilde{\alpha}(X)V^2 + \tilde{\beta}(X)V + \tilde{\gamma}(X), \qquad (14) $$

where $\tilde{\alpha}, \tilde{\beta}$ and $\tilde{\gamma}$ are three other smooth functions. Plugging (13) in (14), we get

$$ \alpha'(X)V^3 + (\beta'(X) - \tilde{\alpha}(X))V^2 + (\gamma'(X) - 2X\alpha(X) - \tilde{\beta}(X))V - (\tilde{\gamma}(X) + X\beta(X)) = 0. $$

The coefficients of this polynomial function in $V$ should vanish, so that

$$ \alpha'(X) = 0, \qquad (15a) $$

$$ \beta'(X) = \tilde{\alpha}(X), \qquad (15b) $$

$$ \gamma'(X) = 2X\alpha(X) + \tilde{\beta}(X), \qquad (15c) $$

$$ \tilde{\gamma}(X) = -X\beta(X). \qquad (15d) $$

From (15a), we deduce that $\alpha$ is constant, and so the coefficient in front of $V^2$ in (13) is constant. Since $C$ represents any invariant of $\langle Q^{\text{filt}} \rangle$, all of its invariants will have a constant coefficient in front of $V^2$. Since $V \partial_X C - X \partial_V C$ is an invariant of $\langle Q^{\text{filt}} \rangle$, $\tilde{\alpha}$ is constant. We can thus deduce from (15b) that $\beta(X)$ is a polynomial of degree 1 (and
thus the same for $\tilde{\beta}$). From (15c), we deduce in the same way that $\gamma$ (and $\tilde{\gamma}$) is a polynomial of degree 2. We have now the generic form of collisional invariants of $\langle Q^{\text{filt}} \rangle$ by replacing $\alpha, \beta$ and $\gamma$ in (13) 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 \mapsto 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 (12), it remains to show that $\mathcal{G}[F]$ and $F$ share the same gyromoments and $\mathcal{T}$ is positive definite. This can be done using classical computations and is postponed in Appendix A. $\square$

**Remark 2.** It is well known (see [6]) 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 $\mathcal{G}: F \to \mathcal{G}[F]$ defined by (11) and (12). We can do the same with the filtered-averaged Maxwellian $\langle \mathcal{M}^{\text{filt}} \rangle: F \to \langle \mathcal{M}^{\text{filt}} \rangle[F]$ defined by (10). It is worth to note that $\mathcal{G}[F]$ and $\langle \mathcal{M}^{\text{filt}} \rangle[F]$ are different except when $F$ is a gyromaxwellian. This is related to the non-commuting limits $t \to \infty$ and $\varepsilon \to 0$ in the filter model (5), this has been discussed in [16]. We detail this point further in Section 6.1.

It is worth to note that gyromoments (11) are also conserved when $F$ solves the filtered equation (5) for all $\varepsilon > 0$. Indeed, from Assertion 2 of Proposition 1, the collisional invariants of $\langle Q^{\text{filt}} \rangle$ are also collisional invariants of $Q$ when they are unfiltered form.

**Proposition 2.** $\forall \varepsilon \ge 0$, the gyromoments (11) are conserved $(\frac{d}{dt}(\varrho, \varrho\mathcal{U}, \varrho\mathcal{T}) = 0)$ and $\frac{d}{dt}\mathcal{G}[F] = 0$.

*Proof.* First, we have $\int Q[f]c(z) dz = 0$ for $c(z) = (1, x, v, x^2, xv, v^2)$ and we recall that $Z = e^{t/\varepsilon}z, \forall \varepsilon > 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/\varepsilon$). Then, we get $\int Q^{\text{filt}}[F]c(Z) dZ = \int Q[f]c(e^{t/\varepsilon}z) dz = 0$ by linearity. Integrate (5) against $c(Z)$ directly proves the conservation of gyromoments. As the gyromaxwellian $\mathcal{G}[F]$ depends on time only through gyromoments, we get $\frac{d}{dt}\mathcal{G}[F] = 0$. For $\varepsilon = 0$, $c(Z)$ is the vector of collisional invariants of $\langle Q^{\text{filt}} \rangle$ (from Theorem 1) so the conservation is obtained by definition. $\square$

# 3 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 \in \mathbb{R}^d$ where $d$ is the dimension, and all the integrals are considered over $\mathbb{R}^d$. As mentioned above, we consider particle methods in this work and will adapt the numerical method developed for the Landau operator in [8, 28, 36] 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.

## 3.1 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

$$ \partial_t f = Q(f) = T\Delta_v f + \nabla_v \cdot ((v-u)f), \quad f(t=0, v) = f_0(v), \quad v \in \mathbb{R}^d. \qquad (16) $$

As there is no transport in (16), the moments $\int 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 [28, 32]). In the Fokker-Planck case, the bracket is given by, for two functionals $\mathcal{F}$ and $\mathcal{G}$ depending on $f$,

$$ (\mathcal{F}, \mathcal{G})(f) = - \int f T \nabla_v \frac{\delta \mathcal{F}}{\delta f} \cdot \nabla_v \frac{\delta \mathcal{G}}{\delta f} \, dv, \qquad (17) $$
where $\frac{\delta \mathcal{F}}{\delta f}$ is the Fréchet derivative of $\mathcal{F}$. The driven entropy $S_M$ is the relative entropy to the Maxwellian $M$:

$$S_M[f] = \int f \log(f/M) \, dv.$$

The evolution of any functional $\mathcal{F}$ of $f$ is thus given by

$$\partial_t \mathcal{F}(f) = (\mathcal{F}, S_M)(f). \qquad (18)$$

For example, (16) can be recovered using the following functional $\mathcal{F}(f)(v) = \int \delta(v - v')f(v') \, dv'$ where $\delta(v)$ denotes the Dirac distribution. We also define the physical entropy

$$S[f] = \int f \log(f) \, dv. \qquad (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 $\rho$, $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

$$\frac{d}{dt} \int f \, dv = 0, \qquad \frac{d}{dt} \int v f \, dv = 0, \qquad \frac{d}{dt} \int \frac{|v|^2}{2} f \, dv = 0, \qquad \frac{d}{dt} S_M \le 0.$$

Moreover, the only equilibrium of the Fokker-Planck operator $Q$ is $f = M$.

*Proof.* The proof, using bracket formulation (18), is given in Appendix B.1. $\square$

We then consider (16) in the following advective form

$$\partial_t f = \nabla_v \cdot (f T \nabla_v \log(f/M)) = \nabla_v \cdot \left( f T \nabla_v \frac{\delta S_M}{\delta f} \right), \qquad (20)$$

where the advective field is the functional derivative of the relative entropy $T \frac{\delta S_M}{\delta f} = T \nabla_v \log(f/M)$.

## 3.2 Regularized Fokker-Planck operator

In view of using particle methods, (20) is not well adapted due to the presence of the logarithm. The idea is to regularize the distribution function using convolutions with a mollifier [7, 8, 28, 36] given by

$$\psi_{\mu}(v) = \frac{1}{(2\pi\mu)^{d/2}} e^{-\frac{|v|^2}{2\mu}}, \qquad (21)$$

where $\mu > 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^\mu$ given by

$$S_M^\mu[f] = \int f * \psi_\mu \log(f * \psi_\mu/M) \, dv, \qquad (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)

$$\partial_t f = \nabla_v \cdot \left( f T \nabla_v \frac{\delta S_M^\mu}{\delta f} \right) = \nabla_v \cdot \left[ f T \psi_\mu * \nabla_v \log(f * \psi_\mu/M) \right]. \qquad (23)$$

We can observe that the drift part is not affected by the regularization. Indeed, we have $\nabla_v \log(M) = -(v-u)/T$ so using properties of Gaussian, we have

$$-\psi_{\mu} * \nabla_{v} \log(M) = -\frac{1}{T} \left( \int w \psi_{\mu}(v-w) \, dw - u \int \psi_{\mu}(v-w) \, dw \right) = -\frac{v-u}{T}.$$
Consequently, (23) can be equivalently written as

$$ \partial_t f = (f, S_M^\mu) = \nabla_v \cdot \left[ f(T\psi_\mu * \nabla_v \log(f * \psi_\mu) + (v-u)) \right]. \quad (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 (17) is replaced by $T^M$ so that we recover (24) with $T = T^M$.

**Proposition 3.** (Properties of the regularized collision operator)

Let $\mathcal{M}$ be a Maxwellian characterized by $\rho$, $u$ and $T^M$ and let $f$ be the solution of (24) with mass $\rho$, mean velocity $u$ and temperature $T^f = \frac{1}{d\rho} \int |v-u|^2 f \, dv$. Initially, we assume $T^f = T^M$ at $t=0$. Then we have the following properties

$$ \frac{d}{dt} \int f \, dv = 0, \quad \frac{d}{dt} \int v f \, dv = 0, \quad \frac{d}{dt} S_M^\mu \le 0, \qquad (25) $$

and

$$ \frac{d}{dt} \int \frac{|v|^2}{2} f \, dv = d\rho (T^M - T^f) + \mu T^M D(f), \quad \text{with} \quad D(f) = \frac{1}{T^M} (S, S) (f * \psi_\mu). \qquad (26) $$

Moreover, the only equilibrium of the regularized Fokker-Planck operator is characterized by $f * \psi_\mu = \mathcal{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 = \frac{d\rho}{d\rho + \mu D(f)} T^f. \qquad (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 (19), that $\partial_t g = (g, S)(g)$. We deduce that the time evolution of the entropy $S[g]$ is then given by $\partial_t S[g] = (S, S)(g)$. To compute $D(f) = (S, S)(f * \psi_\mu)$, it is tempting to choose $g = f * \psi_\mu$, but $f$ is not the solution of the heat equation on $v \in \mathbb{R}^d$ (with vanishing condition at infinity). However, the following Proposition enlightens us.

**Proposition 4.** Let $g_t(\tau, v)$ be the solution of the heat equation $\partial_\tau g_t = \Delta_v g_t$ ($v \in \mathbb{R}^d$), with initial condition $g_t(\tau=0, v) = (f * \psi_\mu)(t, v)$. Then $D(f) = \frac{d}{d\tau} S[g_t]|_{\tau=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 $\partial_\tau g_t = (g_t, S)(g_t)$ and the evolution of $S[g_t(\tau, \cdot)]$ is given by $\frac{d}{d\tau} S[g_t(\tau, \cdot)] = (S, S)(g_t(\tau, \cdot))$. Thanks to the choice of the initial condition $g_t(\tau=0, v) = (f * \psi_\mu)(t, v)$, taking $\tau=0$ in the latter equation gives $\frac{d}{d\tau} S[g_t(\tau, \cdot)]|_{\tau=0} = (S, S)(g_t(\tau=0, \cdot)) = (S, S)((f * \psi_\mu)(t, \cdot)) = D(f)$. □

Thanks to Proposition 4, we can deduce an explicit expression for $D(f)$ defined in (26).

**Corollary 1.** $D(f)$ satisfies $D(f) = \int (f * \Delta_v \psi_\mu) \log(f * \psi_\mu) \, dv$ with $\Delta_v \psi_\mu = \frac{1}{\mu} \left( \frac{|v|^2}{\mu} - d \right) \psi_\mu$.

*Proof.* Let $g_t$ defined in Proposition 4, then we have

$$ \frac{d}{d\tau} S[g_t] = \int_{\mathbb{R}^d} \partial_\tau g_t (\log(g_t) + 1) \, dv = \int_{\mathbb{R}^d} \Delta_v g_t \log(g_t) \, dv. $$

Evaluating at $\tau=0$ and using the initial condition $g_t(\tau=0, v) = (f * \psi_\mu)(t, v)$ from Proposition 4 leads to

$$ D(f) = \frac{d}{d\tau} S[g_t]|_{\tau=0} = \int (f * \Delta_v \psi_\mu) \log(f * \psi_\mu) \, dv. $$

Using the expression (21) of $\psi_\mu$ directly gives $\Delta_v \psi_\mu$. □
## 3.3 Particle discretization

In this part, we present the particle method. Starting from (24), we approximate the continuous distribution function $f$ by a particle density given by

$$f_{N_p}(t, v) = \frac{\rho}{N_p} \sum_{p=1}^{N_p} \delta(v - v_p(t)), \qquad (28)$$

where $\rho = \int 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, \dots, N_p$,

$$\dot{v}_p(t) = -U_\mu[f_{N_p}](t, v_p(t)), \quad \text{with} \quad U_\mu[f_{N_p}](t, v_p) = T\psi_\mu * \nabla_v \log(f_{N_p} * \psi_\mu)(v_p) + (v_p - u). \qquad (29)$$

The particle density is solution in the sense of distribution of the regularized equation (24). In this work, we choose to take particles with same weights $\rho/N_p$ but the method also works with non constant weights. The computation of the regularized part of the advection field $U_\mu[f_{N_p}](t, v_q)$ can be done using a Gauss-Hermite quadrature. Indeed, using a change of variable, the first part can be written as

$$\begin{aligned} \psi_\mu * \nabla_v \log(f_{N_p} * \psi_\mu)(v_p) &= (\nabla_v \psi_\mu) * \log(f_{N_p} * \psi_\mu)(v_p) \\ &= \sqrt{\frac{2}{\mu\pi^d}} \int_{\mathbb{R}^d} e^{-|w|^2} w \log(f_{N_p} * \psi_\mu)(v_p + \sqrt{2\mu}w) dw = \int_{\mathbb{R}^d} e^{-|w|^2} g(w) dw, \end{aligned} \qquad (30)$$

where $g(w) = \sqrt{\frac{2}{\mu\pi^d}} w \log(f_{N_p} * \psi_\mu)(v_p + \sqrt{2\mu}w)$. The form of the integral motivates the use of Gauss-Hermite quadrature:

$$\int_{\mathbb{R}^d} e^{-|w|^2} g(w) dw \approx \sum_{i_1, i_2, \dots, i_d = 1}^{N_q} \omega_{i_1} \omega_{i_2} \dots \omega_{i_d} g(c_{i_1}, c_{i_2}, \dots, c_{i_d}),$$

where $N_q$ denotes the number of weights and points, $\omega_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_{N_p}$ and second point of Lemma 1 in Appendix B.2, we get

$$D(f_{N_p}) = \frac{\rho}{N_p} \sum_{p=1}^{N_p} \left( \Delta_v \psi_\mu * \log(f_{N_p} * \psi_\mu) \right) (v_p),$$

where the first convolution is approximated using Gauss-Hermite quadrature.

## 3.4 Time discretization

In this section, we investigate the impact of a forward Euler time discretization applied to (29). Let us mention that implicit integrators are used in [36] for Landau operator and [30] for Lenard-Bernstein operator where similar strategies are employed. Let consider $v_p^n \approx v_p(t^n)$ with $t^n = n\Delta t$ and $\Delta t > 0$ the time step, and $f_{N_p}^n$ given by (28) where $v_p(t)$ are replaced by $v_p^n$. Then we consider the following time discretization

$$\frac{v_p^{n+1} - v_p^n}{\Delta t} = -T^M \nabla_v \psi_\mu * \log(f_{N_p}^n * \psi_\mu)(v_p^n) - (v_p^n - u). \qquad (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 \in \mathbb{N}$, let define the discrete mass, momentum and energy as

$$\rho^n = \int_{\mathbb{R}^d} f_{N_p}^n(v) dv, \quad P^n = \int_{\mathbb{R}^d} f_{N_p}^n(v) v dv = \frac{\rho}{N_p} \sum_{p=1}^{N_p} v_p^n, \quad E^n = \int_{\mathbb{R}^d} f_{N_p}^n(v) \frac{|v|^2}{2} dv = \frac{\rho}{N_p} \sum_{p=1}^{N_p} \frac{|v_p^n|^2}{2},$$
and consider that particles evolution is given by numerical scheme (31). Suppose that initially, $\rho^0 = \rho$ and $P^0 = \rho u$. Then $\forall n \in \mathbb{N}$ we have

$$ \rho^n = \rho, \quad P^n = \rho u, \quad \frac{E^{n+1} - E^n}{\Delta t} = d\rho (T^M - T^{f,n}) + \mu(S, S)(f_{N_p}^n * \psi_\mu) - \frac{\Delta t}{2} T^M (S_M^\mu, S_M^\mu)(f_{N_p}^n), \quad (32) $$

where $T^{f,n}$ is defined by $E^n = \rho(|u|^2 + dT^{f,n})/2$, $T^M$ is the constant temperature of $M$ and the bracket, $S$ and $S_M^\mu$ are respectively defined by (17), (19) and (22).

*Proof.* It is clear that the mass $\rho$ is preserved. For momentum, apply $\rho \sum_{p=1}^{N_p}$ on (31) gives

$$ \begin{aligned} \frac{P^{n+1} - P^n}{\Delta t} &= -T^M \frac{\rho}{N_p} \sum_{p=1}^{N_p} \nabla_v \psi_\mu * \log(f_{N_p}^n * \psi_\mu)(v_p^n) - (P^n - \rho u) \\ &= -T^M \int_{\mathbb{R}^d} f_{N_p}^n \nabla_v \psi_\mu * \log(f_{N_p}^n * \psi_\mu)(v) \, dv - (P^n - \rho u) \\ &= -T^M \int_{\mathbb{R}^d} (f_{N_p}^n * \psi_\mu) \frac{\nabla_v f_{N_p}^n * \psi_\mu}{f_{N_p}^n * \psi_\mu} \, dv - (P^n - \rho u) = -(P^n - \rho u). \end{aligned} $$

Since $P^0 = \rho u$, we show by induction $P^n = \rho u$ for all $n \in \mathbb{N}$.
Regarding the energy, we multiply (31) by $(v_p^n + v_p^{n+1})/2$ and sum over $p$ to get

$$ \begin{aligned} \frac{E^{n+1} - E^n}{\Delta t} &= -T^M \frac{\rho}{2N_p} \sum_{p=1}^{N_p} (v_p^n + v_p^{n+1}) \cdot \nabla_v \psi_\mu * \log(f_{N_p}^n * \psi_\mu / M)(v_p^n) \\ &= -T^M \frac{\rho}{N_p} \sum_{p=1}^{N_p} v_p^n \cdot \nabla_v \psi_\mu * \log(f_{N_p}^n * \psi_\mu / M)(v_p^n) + \Delta t \frac{(T^M)^2}{2} \frac{\rho}{N_p} \sum_{p=1}^{N_p} \left| \nabla_v \psi_\mu * \log(f_{N_p}^n * \psi_\mu / M)(v_p^n) \right|^2 \\ &= -T^M \int_{\mathbb{R}^d} f_{N_p}^n v \cdot \nabla_v \psi_\mu * \log(f_{N_p}^n * \psi_\mu / M)(v) \, dv + \Delta t \frac{(T^M)^2}{2} \int_{\mathbb{R}^d} f_{N_p}^n \left| \nabla_v \psi_\mu * \log(f_{N_p}^n * \psi_\mu / M)(v) \right|^2 \, dv. \end{aligned} $$

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^\mu, S_M^\mu)(f_{N_p}^n)$. Eventually, we get (32). $\square$

In conclusion, the time discretization adds an extra term which also contributes to energy dissipation. As in the semi-continuous case, we will use (32) 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 = \frac{1}{T^M}(S, S)(f_{N_p}^n * \psi_\mu) \quad \text{and} \quad \tilde{D}^n = \frac{1}{T^M}(S_M^\mu, S_M^\mu)(f_{N_p}^n). \quad (33) $$

The division by $T^M$ makes $D^n$ independent of $T^M$ whereas $\tilde{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} = \frac{d\rho + \mu D^n - \sqrt{(d\rho + \mu D^n)^2 - 2d\rho \Delta t T^{f,n} \tilde{D}^n}}{\Delta t \tilde{D}^n}. \quad (34) $$

*Proof.* Using (33), (32) can be reformulated as

$$ \frac{E^{n+1} - E^n}{\Delta t} = -\frac{\Delta t}{2} \tilde{D}^n (T^M)^2 + (d\rho + \mu D^n) T^M - d\rho T^{f,n}, $$
and the energy is preserved if $T^M$ is defined as the root of the second order polynomial

$$T_{\pm}^{M} = \frac{d\rho + \mu D^{n} \pm \sqrt{(d\rho + \mu D^{n})^{2} - 2d\rho \Delta t T^{f,n} \tilde{D}^{n}}}{\Delta t \tilde{D}^{n}}$$

Notice that the discriminant is positive since $\tilde{D}^n \le 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). $\square$

**Remark 5.** It is clear from (32) that when $\Delta t \to 0$ or $\mu \to 0$, the energy preservation is recovered. We can also observe that when $\mu \to 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} + \mathcal{O}(\Delta t)$ since

$$\sqrt{(d\rho)^2 - 2d\rho \Delta t T^{f,n} \tilde{D}^n} = d\rho \sqrt{1 - 2 \frac{\Delta t T^{f,n} \tilde{D}^n}{d\rho}} = d\rho - \Delta t T^{f,n} \tilde{D}^n + \mathcal{O}(\Delta t^2),$$

so that

$$T_{-}^{M} = \frac{d\rho - d\rho + \Delta t T^{f,n} \tilde{D}^{n}}{\Delta t \tilde{D}^{n}} + \mathcal{O}(\Delta t) = T^{f,n} + \mathcal{O}(\Delta t).$$

**Remark 6.** Actually, the correction (34) is implicit since $\tilde{D}^n$, given by (33), also depends on $T^{M,n}$ through $\mathcal{M}$. We propose to approximate $\tilde{D}^n$ by $\hat{D}^n$ (which will be discussed below) in (34); 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 $\hat{D}^n$ is computed.

From bracket properties, we have $(S_M^\mu, S_M^\mu)(f_{N_p}(t^n)) = \frac{d}{dt}S_M^\mu[f_{N_p}]|_{t=t^n}$; then, we approximate the time derivative by $(S_M[f_{N_p}^{n+\delta t}] - S_M[f_{N_p}^n])/\delta t$ where $\delta t$ is a small time step and where $f_{N_p}^{n+\delta t}$ is obtained by advancing the particles from $t^n$ to $t^n + \delta 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 $\hat{D}^n = \frac{1}{T^{M,n-1}}(S_M[f_{N_p}^{n+\delta t}] - S_M[f_{N_p}^n])/\delta t$ as an approximation of $\tilde{D}^n$ in (34) so that $T^{M,n}$ can be computed explicitly.

Of course, solving the nonlinear relation (34) 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 $\tilde{D}$ becomes small). In this regime, we replace correction (34) by the semi-continuous one (27) since the extra term in (32) vanishes.

## 3.5 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 $\mathcal{O}(N_p^2)$ with $N_p$ the number of particles. More precisely, for each particle $p$, one has to compute $\nabla_v \psi_\mu * \log(f * \psi_\mu)(v_p)$: the first convolution is computed using a $d$-dimensional Gauss-Hermite quadrature with $N_q$ points and $f * \psi_\mu$ is computed with a sum over the $N_p$ particles. Finally, the complexity is $\mathcal{O}(N_q^d N_p^2)$ where $d$ is the phase space dimension. Some attempts to accelerate have been proposed using treecode in [8] or GPU programming in [36]. 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 * \psi_\mu$ in the logarithm convolution part. Indeed, let introduce $\tilde{N}_p < N_p$. We then split the set of $N_p$ particles into $M$ clusters of $\tilde{N}_p$ particles ($N_p = M\tilde{N}_p$). Let call $C_i$ the $i$-th cluster for $i = 1, \dots, M$ and the corresponding regularized unknown in each cluster writes

$$f_{\tilde{N}_p}^i * \psi_\mu(v) = \frac{\rho}{\tilde{N}_p} \sum_{q \in C_i} \psi_\mu(v - v_q).$$

Thanks to this definition, we replace (29) by the following time evolution of the particles velocities

$$\dot{v}_p(t) = -T \nabla_v \psi_\mu * \log(f_{\tilde{N}_p}^i * \psi_\mu)(v_p) - (v_p(t) - u(t)), \quad \forall p \in C_i.$$

For this equation, the complexity becomes $\mathcal{O}(N_q^d N_p \tilde{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 $\rho$ is just a constant number in the scheme. For momentum, denote $\rho u_i$ the momentum of cluster $C_i$ (momentum computed with $f_{\tilde{N}_p}^i$). As $u_i \neq u$ a priori, this momentum is not conserved and we have

$$ \partial_t (\rho u_i) = \rho(u_i - u). \qquad (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,

$$ \frac{1}{M} \sum_{i=1}^{M} u_i = \frac{1}{M} \sum_{i=1}^{M} \frac{1}{\tilde{N}_p} \sum_{q \in C_i} v_q = \frac{1}{\tilde{N}_p} \sum_{q=1}^{N_p} v_q = u. $$

Thus, we get

$$ \partial_t \rho u = \partial_t \frac{\rho}{\tilde{N}_p} \sum_{q=1}^{N_p} v_q = \frac{1}{M} \sum_{i=1}^{M} \partial_t \frac{\rho}{\tilde{N}_p} \sum_{q \in C_i} v_q = \frac{1}{M} \sum_{i=1}^{M} \rho(u_i - u) = \frac{1}{M} \sum_{i=1}^{M} \rho u_i - \rho 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.

# 4 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^{\text{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) \in \mathbb{R}^2$)

$$ Q^{\text{filt},\mu}[F](t, s, Z) = \nabla_Z \cdot (T A(s) F \nabla_Z \psi_{\mu} * \log(F * \psi_{\mu})(Z)) + \nabla_Z \cdot (A(s)(Z - a(s)u)F). \qquad (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 $\epsilon$, using techniques from [9, 17, 19].

## 4.1 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 $\psi_{\mu}: w \in \mathbb{R}^2 \to \psi_{\mu}(w) \in \mathbb{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 (36). This is explained in the Proposition below.

**Proposition 6.** Consider $f$ and $F$ two densities such that $F(Z) = f(e^{-\frac{t}{\epsilon}J}Z)$. Then we have the following equalities (with $z = (x, v) \in \mathbb{R}^2$)

$$ \psi_{\mu} * \log(f * \psi_{\mu})(z) = \psi_{\mu} * \log(F * \psi_{\mu})(Z), \qquad Q^{\mu}[f](t, z) = Q^{\text{filt},\mu}[F](t, t/\epsilon, Z). $$

*Proof.* The proof relies on the change of variable $\xi = e^{\frac{t}{\epsilon}J}\zeta$ which preserves the volume. Indeed, from the first convolution, we have

$$ \begin{aligned} f * \psi_{\mu}(w) &= \int f(w - \zeta) \psi_{\mu}(\zeta) \, d\zeta = \int f(e^{-\frac{t}{\epsilon}J} W - \zeta) \psi_{\mu}(\zeta) \, d\zeta \\ &= \int f(e^{-\frac{t}{\epsilon}J}(W - \zeta)) \psi_{\mu}(\zeta) \, d\zeta = \int F(W - \zeta) \psi_{\mu}(\zeta) \, d\zeta = F * \psi_{\mu}(W). \end{aligned} $$

Consequently, the double convolution becomes

$$ \begin{aligned} \psi_{\mu} * \log(f * \psi_{\mu})(z) &= \int \psi_{\mu}(z - w) \log(f * \psi_{\mu}(w)) \, dw = \int \psi_{\mu}(e^{-\frac{t}{\epsilon}J} Z - w) \log(f * \psi_{\mu}(w)) \, dw \\ &= \int \psi_{\mu}(Z - W) \log(F * \psi_{\mu}(W)) \, dW = \psi_{\mu} * \log(F * \psi_{\mu})(Z). \end{aligned} $$

The second equality is directly deduced from the first one. $\square$
Consider the particle density $f_{N_p}$ and the regularized particle density $F_{N_p}$ defined by

$$f_{N_p}(t, z) = \frac{\varrho}{N_p} \sum_{p=1}^{N_p} \delta(z - z_p(t)) \quad \text{and} \quad F_{N_p}(t, Z) = \frac{\varrho}{N_p} \sum_{p=1}^{N_p} \delta(Z - Z_p(t)), \qquad (37)$$

where the filtered particles are given by $Z_p(t) = e^{\frac{t}{\varepsilon}J}z_p(t)$ and $\varrho$ 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 (1) where $Q$ is replaced by its regularization version (24) gives

$$\dot{z}_p(t) = -\frac{1}{\varepsilon} J z_p(t) - U_\mu[f_{N_p}](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

$$\dot{Z}_p(t) = -U_\mu^{\text{filt}}[F_{N_p}](t, t/\varepsilon, Z_p(t)), \qquad (38)$$

where

$$U_{\mu}^{\text{filt}}[F_{N_p}](t, s, Z_p(t)) = T\mathcal{A}(s)\nabla_Z\psi_{\mu} * \log(F_{N_p} * \psi_{\mu})(Z_p) + \mathcal{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 (27) to preserve the quadratic gyromoments. Other gyromoments (i.e. those associated to $1, X, V$) are conserved up to the quadrature error.

## 4.2 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/\varepsilon$ is periodic of period $2\pi$, which enables to use the double scale strategy introduce in [9, 17, 19]. The main idea is to consider the fast scale $t/\varepsilon$ as an independent variable $s$ so that the original solution is recovered by an evaluation on the diagonal $s = t/\varepsilon$. By an abuse of notation, we still denote $Z_p(t, s)$ the double scale unknown which satisfies the following set of PDEs

$$\partial_t Z_p + \frac{1}{\varepsilon} \partial_s Z_p = -U_{\mu}^{\text{filt}}[F_{N_p}](t, s, Z_p).$$

We set $\mathcal{G}_p(t, s) = -U_{\mu}^{\text{filt}}[F_{N_p}](t, s, Z_p)$. Since $Z_p$ is periodic with respect to $s$, Fourier expansion is used in this direction: $Z_p(t, s) = \sum_{\ell} e^{i\ell s} \hat{Z}_{p,\ell}(t)$. We then get the following ODE on the Fourier modes $\hat{Z}_{p,\ell}$ ($\hat{\mathcal{G}}_{p,\ell}(t)$ denotes the Fourier modes associated to $\mathcal{G}_p(t, s)$)

$$\partial_t \hat{Z}_{p,\ell}(t) + \frac{i\ell}{\varepsilon} \hat{Z}_{p,\ell}(t) = \hat{\mathcal{G}}_{p,\ell}(t), \quad \ell \in \mathbb{Z}^*,$$

and $\partial_t \hat{Z}_{p,0}(t) = \hat{\mathcal{G}}_{p,0}(t)$ (which actually corresponds to the averaged model (7)).

Now we introduce the discrete unknown $\hat{Z}_{p,\ell}^n \approx \hat{Z}_{p,\ell}(t^n)$, $t^n = n\Delta t$, with $\Delta t > 0$ the time step. Applying the Duhamel formula combined with exponential integrators [29] finally gives a first order scheme

$$\hat{Z}_{p,\ell}^{n+1} = e^{-\frac{i\ell}{\varepsilon}\Delta t}\hat{Z}_{p,\ell}^n + \Delta t\varphi_1\left(-\frac{i\ell}{\varepsilon}\Delta t\right)\hat{\mathcal{G}}_{p,\ell}^n, \qquad (39)$$

whereas a second order scheme is given by

$$\hat{Z}_{p,\ell}^{n+1} = e^{-\frac{i\ell}{\varepsilon}\Delta t}\hat{Z}_{p,\ell}^n + \Delta t\varphi_1\left(-\frac{i\ell}{\varepsilon}\Delta t\right)\hat{\mathcal{G}}_{p,\ell}^n + \Delta t\varphi_2\left(-\frac{i\ell}{\varepsilon}\Delta t\right)\left(\hat{\mathcal{G}}_{p,\ell}^n - \hat{\mathcal{G}}_{p,\ell}^{n-1}\right), \qquad (40)$$

where the $\varphi_1, \varphi_2$ functions are given by

$$\varphi_1(z) = \frac{1}{z}(e^z - 1), \qquad \varphi_2(z) = \frac{1}{z^2}(e^z - 1 - z),$$

where $\hat{\mathcal{G}}_{p,\ell}^n$ is the Fourier transform of $-U_{\mu}^{\text{filt}}[F_{N_p}^n](t^n, s, Z_p^n)$. To compute $-U_{\mu}^{\text{filt}}[F_{N_p}^n](t^n, s, Z_p^n)$, the moments $u$ and $T$ need to be computed at time $t^n$. To do so, we unfilter $F_{N_p}^n$ and approximate them using second order spline interpolation on a sufficiently large space grid $[x_{\text{min}}, x_{\text{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_p^{n+1} = e^{-\frac{t^{n+1}}{\varepsilon}J}Z_p^{n+1}(s = t^{n+1}/\varepsilon)$, $Z_p^{n+1}(s)$ is reconstructed using a Fourier truncated series. For more details, we refer to [9, 16, 17, 19].
# 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 $\Delta t$ and the regularization parameter $\mu$. 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) = \frac{1}{2\sqrt{2\pi}} \left( e^{-\frac{(v-u_1)^2}{2}} + e^{-\frac{(v-u_2)^2}{2}} \right), \quad (41)$$

with $u_1 = -2.6$, $u_2 = 3.4$.

## 5.1 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: $\Delta t = 10^{-3}$, $\mu = 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_{N_p}(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_{N_p}(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_{N_p} * \psi_\mu$. 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_{N_p} * \psi$ has almost reached the Maxwellian equilibrium. To quantify this convergence, we plot in Figure 2 the difference in $L^1$, $L^2$ and $L^\infty$ norms between $f_{N_p} * \psi_\mu$ 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 $\psi_\mu$ 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$.

## 5.2 Noise of the method

In this part, we investigate the noise effect according to the numerical parameters $\Delta t$ and $\mu$. 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_{N_p} * \psi_\mu$ is plotted at $t=1.5$ with $N_p = 250$ ($\mu = 0.02$ and $N_q = 50$) and for different time steps: $\Delta t = 10^{-2}$, $1.2 \times 10^{-2}$, $1.5 \times 10^{-2}$, $2 \times 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.

## 5.3 Conserved properties

Here, we investigate the preservation of the moments according to the numerical parameters $\Delta t$ and $\mu$ (the other parameters are $N_p = 250$, $N_q = 50$ and the final time is $t=0.5$). We focus on the momentum $P = \rho u$ and the energy $E$ preservation. To do so, we give in Tables 1, 2 and 3 the following quantities

$$P_{\text{err}} = \max_{n=0,...,N} |P^n - P^0|,$$

$$E_{\text{err}} = \max_{n=0,...,N} |E^n - E^0|,$$

for $\Delta t = 10^{-3}$, $10^{-4}$ and $\mu = 10^{-1}$, $10^{-2}$ for the three different variants of the method presented in Section 3 (without correction, with correction (27) and with correction (34)). Since the initialisation procedure (rejection
<figure><img src="image_2.png" /><figcaption>Figure 1: $f_{N_p} * \psi_\mu$ 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$. $\Delta t = 10^{-3}, \mu = 0.02$, and $N_q = 50$.</figcaption></figure>

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_{N_p} * \psi_\mu - f_0\|_{L^2}^2$ on particles.

According to the variant, the amplitude of the energy error $E_{\text{err}}$ depends on $\Delta t$ and $\mu$. 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 $\tilde{D}$ (see Remark 6) in addition with Gauss-Hermite quadrature.
<figure><img src="image_3.png" /><figcaption>Figure 2: Time evolution of $||f_{N_p} * \psi_\mu - \mathcal{M}||$ for $N_p = 250$ (left) and $N_p = 10^3$ (right).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: $f_{N_p} * \psi_\mu$ for $v \in [6, 10]$: illustration of the bumps around isolated particles.</figcaption></figure>

<table><thead><tr><th>P<sub>err</sub></th><th>&mu; = 10<sup>-1</sup></th><th>&mu; = 10<sup>-2</sup></th></tr></thead><tbody><tr><td>&Delta;t = 10<sup>-3</sup></td><td>8.7714 &times; 10<sup>-9</sup></td><td>5.4729 &times; 10<sup>-10</sup></td></tr><tr><td>&Delta;t = 10<sup>-4</sup></td><td>8.7443 &times; 10<sup>-9</sup></td><td>5.3233 &times; 10<sup>-10</sup></td></tr></tbody></table>

<table><caption>Table 1: Conservation of momentum and energy without correction of $T^M$.</caption><thead><tr><th>E<sub>err</sub></th><th>&mu; = 10<sup>-1</sup></th><th>&mu; = 10<sup>-2</sup></th></tr></thead><tbody><tr><td>&Delta;t = 10<sup>-3</sup></td><td>9.2179 &times; 10<sup>-2</sup></td><td>4.3004 &times; 10<sup>-2</sup></td></tr><tr><td>&Delta;t = 10<sup>-4</sup></td><td>9.2882 &times; 10<sup>-2</sup></td><td>4.3799 &times; 10<sup>-2</sup></td></tr></tbody></table>

<table><thead><tr><th>P<sub>err</sub></th><th>&mu; = 10<sup>-1</sup></th><th>&mu; = 10<sup>-2</sup></th></tr></thead><tbody><tr><td>&Delta;t = 10<sup>-3</sup></td><td>8.4644 &times; 10<sup>-9</sup></td><td>5.5386 &times; 10<sup>-10</sup></td></tr><tr><td>&Delta;t = 10<sup>-4</sup></td><td>8.4531 &times; 10<sup>-9</sup></td><td>5.3865 &times; 10<sup>-10</sup></td></tr></tbody></table>

<table><caption>Table 2: Conservation of momentum and energy with the semi-discrete correction (27) of $T^M$.</caption><thead><tr><th>E<sub>err</sub></th><th>&mu; = 10<sup>-1</sup></th><th>&mu; = 10<sup>-2</sup></th></tr></thead><tbody><tr><td>&Delta;t = 10<sup>-3</sup></td><td>2.5359 &times; 10<sup>-3</sup></td><td>2.3464 &times; 10<sup>-3</sup></td></tr><tr><td>&Delta;t = 10<sup>-4</sup></td><td>2.5235 &times; 10<sup>-4</sup></td><td>2.3337 &times; 10<sup>-4</sup></td></tr></tbody></table>

<table><thead><tr><th>P<sub>err</sub></th><th>&mu; = 10<sup>-1</sup></th><th>&mu; = 10<sup>-2</sup></th></tr></thead><tbody><tr><td>&Delta;t = 10<sup>-3</sup></td><td>8.4758 &times; 10<sup>-9</sup></td><td>5.5325 &times; 10<sup>-10</sup></td></tr><tr><td>&Delta;t = 10<sup>-4</sup></td><td>8.4542 &times; 10<sup>-9</sup></td><td>5.3858 &times; 10<sup>-10</sup></td></tr></tbody></table>

<table><caption>Table 3: Conservation of momentum and energy with the correction (34) of $T^M$, using Remark 6.</caption><thead><tr><th>E<sub>err</sub></th><th>&mu; = 10<sup>-1</sup></th><th>&mu; = 10<sup>-2</sup></th></tr></thead><tbody><tr><td>&Delta;t = 10<sup>-3</sup></td><td>1.7000 &times; 10<sup>-7</sup></td><td>3.5123 &times; 10<sup>-7</sup></td></tr><tr><td>&Delta;t = 10<sup>-4</sup></td><td>4.9941 &times; 10<sup>-8</sup></td><td>1.7070 &times; 10<sup>-7</sup></td></tr></tbody></table>

## 5.4 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$ ($\mu = 0.04$, $\Delta 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 $\tilde{N}_p = 60$ particles
<figure><img src="image_5.png" /><figcaption>Figure 4: $f_{N_p} * \psi_\mu$ as a function of $v$ for different time steps: influence of the time step on the level of noise. $N_p = 250$, $\mu = 0.02$, and $N_q = 50$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: Time history of one particle $v_p$ for different time steps. $N_p = 250$, $\mu = 0.02$, and $N_q = 50$.</figcaption></figure>

(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 $\tilde{N}_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 $O(N_q^2 N_p \tilde{N}_p)$ instead of $O(N_q^2 N_p^2)$). 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.

# 6 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 $\langle M^{\text{filt}} \rangle$ and the
<figure><img src="image_7.png" /><figcaption>(c) Clustering method for $N_p = 1200$ split in 20 clusters of 60 (d) Random clustering method for $N_p = 1200$ split in 20 clusters of 60 particles.</figcaption></figure>

<caption>Figure 6: Comparison between the classical method and the clustering method.</caption>

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.

## 6.1 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 (12) and the filtered-averaged Maxwellian $\langle M^{\text{filt}} \rangle [F]$ defined by (10), both being associated to a density $F$. Moreover, we illustrate numerically the fact that $G[F]$ is a fixed point of the operator $\langle M^{\text{filt}} \rangle : F \mapsto \langle M^{\text{filt}} \rangle [F]$.

First, we recall that computing $G[F]$ from (12) requires to compute the gyromoments (11) of $F$. Second, computing $\langle M^{\text{filt}} \rangle [F]$ requires to compute $M^{\text{filt}}[F](s)$ for all $s \in [0, 2\pi]$ and then to take the average over $s$ as in (10). To compute $M^{\text{filt}}[F](s)$, we first unfilter $F$ by setting $f(z) = F(e^{sJz})$, so that we can compute the moments (2) of $f$, which determine the Maxwellian $M$ given by (4). Finally, $M^{\text{filt}}[F](s)$ given by (9) is nothing but the Maxwellian $M$ expressed in the filtered variables $Z = e^{-sJz}$.

Now, for all $m \in \mathbb{N}^*$, we recursively define $\langle M^{\text{filt}} \rangle^{\circ m} [F]$ by

$$ \langle M^{\text{filt}} \rangle^{\circ m} [F] = \langle M^{\text{filt}} \rangle^{\circ (m-1)} [\langle M^{\text{filt}} \rangle [F]]. $$

Below, we will illustrate the fact that $\langle M^{\text{filt}} \rangle^{\circ 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 $\langle M^{\text{filt}} \rangle^{\circ m} [F]$ for a given density $F(Z)$, $Z = (X, V) \in \mathbb{R}^2$.
We consider the following function (which is bi-Maxwellian in this example)

$$F(X, V) = \frac{1}{2\pi\sqrt{T_{X_1}T_{V_1}}}e^{-\frac{|X-u_{X_1}|^2}{2T_{X_1}}}e^{-\frac{|V-u_{V_1}|^2}{2T_{V_1}}} + \frac{1}{2\pi\sqrt{T_{X_2}T_{V_2}}}e^{-\frac{|X-u_{X_2}|^2}{2T_{X_2}}}e^{-\frac{|V-u_{V_2}|^2}{2T_{V_2}}}$$

with $u_{X_1} = -3, u_{V_1} = -2.5, u_{X_2} = 2.5, u_{V_2} = 3$ and $T_{X_1} = 1.75, T_{V_1} = 1.25, T_{X_2} = 1.5, T_{V_2} = 1.25$. On the one side, as discussed before, the gyromaxwellian $G[F]$ (given by (12)) requires the gyromoments of $F$ defined by $\int 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:

$$\langle \mathcal{M}^{\text{filt}} \rangle [F](X, V) = \frac{1}{2\pi} \int_{0}^{2\pi} \mathcal{M}^{\text{filt}}(s)[F](X, V) \, ds \approx \frac{1}{N_s} \sum_{j=0}^{N_s-1} \mathcal{M}^{\text{filt}}(s_j)[F](X, V) \Delta s, \quad (42)$$

where the integral is approximated using $N_s$ quadrature points $s_j = j\Delta s, j = 0, ..., N_s - 1$ and $\Delta s = 2\pi/N_s$. Now, it remains to compute the filtered-averaged Maxwellian $\mathcal{M}^{\text{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^{s_j J}z), z = (x, v)$ and then compute the moments $\int f_j(x, v)(1, v, v^2)^T \, dv$ of $f_j$. Once these moments are computed, the Maxwellian $\mathcal{M}_j(x, v)$ can be constructed for $j = 0, ..., N_s - 1$ following (4) and finally, we apply the change of variable

$$\mathcal{M}^{\text{filt}}(s_j)[F](Z) = \mathcal{M}_j(e^{-s_j J}Z), \quad j = 0, ..., N_s - 1.$$

We then compute the average using (42). We repeat this procedure iteratively to compute $\langle \mathcal{M}^{\text{filt}} \rangle^{\circ m}[F]$ for $m \ge 1$.

In Figure 7, we plot the level set of $F, \langle \mathcal{M}^{\text{filt}} \rangle^{\circ 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 $\langle \mathcal{M}^{\text{filt}} \rangle^{\circ m}[F]$ is far from the gyromaxwellian shape $G[F]$ for small $m$. However, when $m$ increases, the shape of $\langle \mathcal{M}^{\text{filt}} \rangle^{\circ 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 $\langle \mathcal{M}^{\text{filt}} \rangle^{\circ m}[F] - G[F]$ (using different discrete norms $L^1, L^2$ and $L^\infty$) 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 \mapsto \langle \mathcal{M}^{\text{filt}} \rangle^{\circ m}[F]$, that is $\langle \mathcal{M}^{\text{filt}} \rangle [G[F]] = G[F]$. In terms of operators, we have $\lim_{m \to \infty} \langle \mathcal{M}^{\text{filt}} \rangle^{\circ m} = G$ where $G$ denotes the gyromaxwellian operator $F \to G[F]$.

## 6.2 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 (7) towards the gyromaxwellian equilibrium $G[F]$. To do so, instead of solving the averaged model (7), we solve the filtered model (5) with a small $\varepsilon$ (with a bi-Maxwellian initial condition). Indeed, since the error between the averaged model (7) and the filtered model (5) is $O(\varepsilon)$ (and then negligible) and since our scheme is uniformly accurate in $\varepsilon$, the error is $O(\Delta t^2) + O(\varepsilon)$. The numerical parameters are chosen as follows: $\Delta t = 0.0133, \varepsilon = 10^{-10}, \mu = 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.

## 6.3 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 $\varepsilon$ 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 $\mu = 1$ and $N_q = 20$ whereas the space mesh (used to compute the moments) is defined by $x_{\text{min}} = -10, x_{\text{max}} = 10$ and $N_x = 100$. We consider the final time $t = 1$. We plot the order for different $\varepsilon = 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 $\varepsilon$.
<figure><img src="image_8.png" /><figcaption>Figure 7: Level set of the functions $F$ (top left), $\langle \mathcal{M}^{\text{filt}} \rangle^{\circ m}[F]$ for different $m$ (from top right to bottom left, $m = 1, 2, 3, 4, 8, 12$), and $\mathcal{G}[F]$ (bottom right).</figcaption></figure>

(second order is well recovered for small $\epsilon$). However, if we increase the number of Fourier modes to $N_s = 1024$, second order accuracy is observed uniformly in $\epsilon$ (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 [18]. 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 and 10d, we can observe second order accuracy for all $\epsilon$, even for $N_s = 64$.

## 6.4 A beam of particles with collisions

We consider a particle beam test as studied in [18,34], 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
<figure><img src="image_9.png" /><figcaption>Figure 8: $|| \langle \mathcal{M}^{\text{filt}} \rangle^{\text{om}}[F] - \mathcal{G}[F] ||$ as a function of $m$ (semi-log scale).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 9: Relaxation of the filtered unknown $F$ (red particles) towards its associated gyromaxwellian $\mathcal{G}[F]$.</figcaption></figure>

equation) in the transport part so that the model satisfied by $f(t, r, v)$ ($t \ge 0, r, v \in \mathbb{R}$) is

$$ \partial_t f - \frac{1}{\epsilon} Jz \cdot \nabla_z f + E_f \partial_v f = \nu Q[f], \quad (43) $$
<figure><img src="image_11.png" /><figcaption>Figure 10: Numerical order for the second order scheme.</figcaption></figure>

where $z = (r, v)$, $\nu \ge 0$ denotes the collision frequency. Moreover, the electric field $E_f$ solves a Poisson equation in polar coordinates

$$ \partial_r (r E_f) = r \int_{\mathbb{R}} f \, dv. $$

The radial variable $r > 0$ plays the same role as $x$. By symmetry, we can consider $r \in \mathbb{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 $\nu$ to illustrate the effect of the Fokker-Planck operator ($\nu = 0$ and $\nu = 0.1$) on the filamentation.

For this test, we consider the following initial condition

$$ f_0(x, v) = \frac{n_0}{\sqrt{2\pi v_{\text{th}}^2}} e^{-\frac{1}{2}(v/v_{\text{th}})^2} \mathbb{I}_{[-r_{\text{max}}, r_{\text{max}}]}(x), $$

where $r_{\text{max}} = 1.85$, $v_{\text{th}} = 0.1$, $n_0 = 12$ ($\rho = 2r_{\text{max}}n_0$) and $\mathbb{I}$ is the indicator function. The numerical parameters are chosen as follows: $N_p = 1600$, $\epsilon = 0.001$, $\mu = 0.005$, $N_q = 20$ and $\Delta 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 $\nu = 0$ and $\nu = 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 literature [18, 34]. Indeed, thin filaments are created which are well reproduced by the particle method. Second, when $\nu = 0.1$, the filaments
are dissipated and an equilibrium is reached. Again, the results obtained by the two methods are in very good agreement.

<figure><img src="image_12.png" /><figcaption>Figure 11: Snapshots of the solution of the beam problem (43). First line: particle method $N_p = 1600$, $\nu = 0$ (left) and $\nu = 0.1$ (right). Second line: reference (spectral) method $\nu = 0$ (left) and $\nu = 0.1$ (right).</figcaption></figure>

## 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 $\mathfrak{T}$ defined by (11c) is symmetric positive definite and that

$$ \int_{\mathbb{R}^2} C(X, V)(F(X, V) - \mathcal{G}[F](X, V)) \, dX \, dV = 0, \quad \text{with } C(X, V) = (1, X, V, X^2, XV, V^2)^T. \qquad (44) $$

We start by proving (44). To do so, we prove that the gyromoments of $\mathcal{G}[F]$ are equal to the gyromoments of $F$ denoted by $(\varrho, \varrho U, \varrho T)$ as in (11). First, denoting $Z = (X, V)$ and using the notations of Theorem 1, we have

$$ -\frac{1}{2} Z^T \mathfrak{T}^{-1} Z = \frac{-\mathfrak{T}_{XX}}{2 \det(\mathfrak{T})} \left( V - \frac{\mathfrak{T}_{XV}}{\mathfrak{T}_{XX}} X \right)^2 - \frac{1}{2 \mathfrak{T}_{XX}} X^2. \qquad (45) $$

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

$$ \int_{\mathbb{R}^2} \mathcal{G}[F](Z) \, dZ = \frac{\varrho}{\sqrt{\det(2\pi \mathfrak{T})}} \int_{\mathbb{R}^2} e^{-\frac{1}{2} Z^T \mathfrak{T}^{-1} Z} \, dZ = \varrho. $$
Regarding (11b), we have, still using the change of variable $Z \leftarrow (Z - U)$ and (45)

$$ \begin{aligned} \int_{\mathbb{R}^2} ZG[F](Z) \, dZ &= \frac{\varrho}{\sqrt{\det(2\pi\mathfrak{T})}} \int_{\mathbb{R}^2} Z \exp\left(-\frac{1}{2}Z^T\mathfrak{T}^{-1}Z\right) \, dZ + \varrho\mathcal{U} \\ &= \frac{\varrho}{\sqrt{\det(2\pi\mathfrak{T})}} \int_{\mathbb{R}} e^{\frac{-1}{2\pi\mathfrak{T}}X^2} \int_{\mathbb{R}} Ze^{\frac{-\mathfrak{T}X\mathfrak{X}}{2\det(\mathfrak{T})} \left(V - \frac{\mathfrak{T}X\mathfrak{V}}{\mathfrak{T}X\mathfrak{X}}X\right)^2} \, dV \, dX + \varrho\mathcal{U}. \end{aligned} $$

For $Z = X$, the integral above is zero by oddness in $X$. For $Z = V$, the integral above in $V$ is proportional to $(\mathfrak{T}_{XV}/\mathfrak{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

$$ \int_{\mathbb{R}^2} (Z - U) \otimes (Z - U)G[F](Z) \, dZ = \frac{\varrho}{\sqrt{\det(2\pi\mathfrak{T})}} \int_{\mathbb{R}^2} Z \otimes Z \exp\left(-\frac{1}{2}Z^T\mathfrak{T}^{-1}Z\right) \, dZ. \quad (46) $$

Let focus on the diagonal terms. For the $X^2$ term, we have

$$ \begin{aligned} \frac{\varrho}{\sqrt{\det(2\pi\mathfrak{T})}} \int_{\mathbb{R}^2} X^2 \exp\left(-\frac{1}{2}Z^T\mathfrak{T}^{-1}Z\right) \, dZ &= \frac{\varrho}{\sqrt{\det(2\pi\mathfrak{T})}} \int_{\mathbb{R}} X^2 e^{\frac{-1}{2\pi\mathfrak{T}}X^2} \int_{\mathbb{R}} e^{\frac{-\mathfrak{T}X\mathfrak{X}}{2\det(\mathfrak{T})} \left(V - \frac{\mathfrak{T}X\mathfrak{V}}{\mathfrak{T}X\mathfrak{X}}X\right)^2} \, dV \, dX \\ &= \frac{\varrho}{\sqrt{\mathfrak{T}_{XX}}} \int_{\mathbb{R}} X^2 e^{\frac{-1}{2\pi\mathfrak{T}}X^2} \, dX = \varrho\mathfrak{T}_{XX}, \end{aligned} \quad (47) $$

and symmetrically, for the $V^2$ term, we obtain $\int_{\mathbb{R}^2} V^2G[F](Z) \, dZ = \varrho\mathfrak{T}_{VV}$. For the extra-diagonal terms, we have

$$ \begin{aligned} \frac{\varrho}{\sqrt{\det(2\pi\mathfrak{T})}} \int_{\mathbb{R}^2} XV \exp\left(-\frac{1}{2}Z^T\mathfrak{T}^{-1}Z\right) \, dZ &= \frac{\varrho}{\sqrt{\det(2\pi\mathfrak{T})}} \int_{\mathbb{R}} X e^{\frac{-1}{2\pi\mathfrak{T}}X^2} \int_{\mathbb{R}} V e^{\frac{-\mathfrak{T}X\mathfrak{X}}{2\det(\mathfrak{T})} \left(V - \frac{\mathfrak{T}X\mathfrak{V}}{\mathfrak{T}X\mathfrak{X}}X\right)^2} \, dV \, dX \\ &= \varrho \frac{\mathfrak{T}_{XV}}{\mathfrak{T}_{XX}} \int_{\mathbb{R}} X^2 e^{\frac{-1}{2\pi\mathfrak{T}}X^2} \, dX = \varrho\mathfrak{T}_{XV}, \end{aligned} $$

which enables to get (11c). Let remark that using the above calculations, we also have $\iint_{\mathbb{R}^2} Z \otimes ZG[F](Z) \, dZ = \varrho(\mathcal{U} \otimes \mathcal{U} + \mathfrak{T})$.

We conclude the proof by showing that $\mathfrak{T}$ is positive definite. As $\mathfrak{T}$ is symmetric, it is diagonalizable in $\mathbb{R}$. Using Cauchy-Schwarz inequality, we get

$$ \left( \iint (X - \mathcal{U}_X)(V - \mathcal{U}_V)F \, dX \, dV \right)^2 < \iint (X - \mathcal{U}_X)^2 F \, dX \, dV \iint (V - \mathcal{U}_V)^2 F \, dX \, dV, $$

and so $\det(\mathfrak{T}) > 0$ (the equality case is not reached because $(X - \mathcal{U}_X)\sqrt{F}$ and $(V - \mathcal{U}_V)\sqrt{F}$ are linearly independant). Moreover, we have $\mathrm{tr}(\mathfrak{T}) \ge 0$. We deduce that the two eigenvalues are positive and thus, $\mathfrak{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 (17) and the reformulated equation (18) satisfied by any functional $\mathcal{F}$ of $f$, in particular the moments and the entropy. For the entropy inequality, as the bracket is negative, we deduce it immediately from $\frac{d}{dt}S_M = (S_M, S_M) \le 0$.

For the mass conservation, we have $\frac{\delta\rho}{\delta f} = 1$ and the mass is conserved from $(\rho, S_M) = 0$.

For momentum, we have $\frac{\delta(\rho u_i)}{\delta f} = v_i$. It follows that

$$ \begin{aligned} (\rho u_i, S_M) &= - \int_{\mathbb{R}^d} T f (\nabla_v v_i \cdot \nabla_v \log(f/M)) \, dv = - \int_{\mathbb{R}^d} T f \partial_{v_i} \log(f/M) \, dv \\ &= - \int_{\mathbb{R}^d} T \partial_{v_i} f + f(v_i - u_i) \, dv = - \int_{\mathbb{R}^d} f v_i \, dv + \int_{\mathbb{R}^d} f \, dv \, u_i = 0. \end{aligned} $$
For energy, from the relation $\frac{\delta(\rho|u|^2+d\rho T)}{\delta f} = |v|^2$, we get

$$ \begin{aligned} \frac{1}{2}(\rho(|u|^2 + dT), S_M) &= -\frac{1}{2} \int_{\mathbb{R}^d} T f \left( \nabla_v |v|^2 \cdot \nabla_v \log(f/M) \right) \, dv = -\sum_{i=1}^d \int_{\mathbb{R}^d} T v_i f \partial_{v_i} \log(f/M) \, dv \\ &= -\sum_{i=1}^d \int_{\mathbb{R}^d} T v_i \partial_{v_i} f + f v_i (v_i - u_i) \, dv = dT \rho - \rho (|u|^2 + dT) + \rho |u|^2 = 0. \end{aligned} $$

Let us focus on the equilibrium. It is clear that $M$ satisfies $Q[M] = 0$. Suppose now that $f$ satisfies $Q[f] = 0$. Thus we have $\int_{\mathbb{R}^d} Q[f] \log(f/M) \, dv = 0$; using (20) and an integration by part, we get

$$ 0 = \int_{\mathbb{R}^d} Q[f] \log(f/M) \, dv = - \int_{\mathbb{R}^d} T f \left| \nabla_v \frac{\delta S_M}{\delta f} \right|^2 \, dv = (S_M, S_M) = \frac{d}{dt} S_M. $$

It remains to prove that if $\frac{d}{dt}S_M[f] = 0$, then $f = M$. As $f > 0$ and $T > 0$, then $\left| \nabla_v \frac{\delta S_M}{\delta f} \right|^2$ has to be zero so that

$$ \nabla_v \frac{\delta S_M}{\delta f} = \frac{\nabla_v f}{f} + \frac{v-u}{T} = 0. $$

The solution are $f(v) = C e^{-\frac{|v-u|^2}{2T}}$ with $C \in \mathbb{R}$ determined by integrating on $v \in \mathbb{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.

**Lemma 1.** For all functions $f,g: \mathbb{R}^d \mapsto \mathbb{R}$, we have the following relations:

$$ \bullet \int f (\psi_{\mu} * g) \, dv = \int (f * \psi_{\mu}) g \, dv, $$

$$ \bullet \int f (\Delta_v \psi_{\mu} * g) \, dv = \int (f * \Delta_v \psi_{\mu}) g \, dv, $$

$$ \bullet \int f * \psi_{\mu} \, dv = \int f \, dv, $$

$$ \bullet \int (f * \psi_{\mu}) v \, dv = \int f v \, dv, $$

$$ \bullet \int (f * \psi_{\mu}) |v|^2 \, dv = \int f |v|^2 \, dv + d\mu \int f \, dv, $$

$$ \bullet \int v (f * g) \, dv = \int ((vf) * g) \, dv + \int (f * (vg)) \, dv. $$

**Proof of Theorem 3.**

The mass conservation in (25) is obtained similarly as for the original Fokker-Planck operator. For the momentum conservation (25), we replace $f$ by $\int f v_i \, dv$ ($i = 1, \dots, d$) in the the bracket form (24) to get

$$ \begin{aligned} \partial_t (\rho u_i) &= \left( \int_{\mathbb{R}^d} f v_i \, dv, S_M^\mu \right) = -T^M \int_{\mathbb{R}^d} f \partial_{v_i} \left( \psi_\mu * \log \left( f * \psi_\mu / M \right) \right) \, dv \\ &= - \int_{\mathbb{R}^d} T^M f \partial_{v_i} \left( \psi_\mu * \log \left( f * \psi_\mu \right) \right) \, dv + \int_{\mathbb{R}^d} f (v_i - u_i) \, dv = \int_{\mathbb{R}^d} T^M f \left( \psi_\mu * \frac{\partial_{v_i} f * \psi_\mu}{f * \psi_\mu} \right) \, dv \\ &= \int_{\mathbb{R}^d} T^M (f * \psi_\mu) \left( \frac{\partial_{v_i} f * \psi_\mu}{f * \psi_\mu} \right) \, dv = \int_{\mathbb{R}^d} T^M \partial_{v_i} f * \psi_\mu \, dv = 0. \end{aligned} $$
Considering now the energy relation, we proceed as follows (replacing $f$ by $\int f|v|^2/2$ $dv$ in (24)) to get

$$
\begin{align*}
\frac{1}{2} \partial_t (\rho(|u|^2 + dT^f)) &= \left( \int_{\mathbb{R}^d} f \frac{|v|^2}{2} \, dv, S_{\mathcal{M}}^{\mu} \right) = -T^{\mathcal{M}} \sum_{i=1}^{d} \int_{\mathbb{R}^d} f v_i \partial_{v_i} (\psi_{\mu} * \log(f * \psi_{\mu}/\mathcal{M})) \, dv \\
&= -T^{\mathcal{M}} \sum_{i=1}^{d} \int_{\mathbb{R}^d} f v_i \partial_{v_i} (\psi_{\mu} * \log(f * \psi_{\mu})) \, dv - \sum_{i=1}^{d} \int_{\mathbb{R}^d} f v_i (\psi_{\mu} * (v_i - u_i)) \, dv \\
&= -T^{\mathcal{M}} \sum_{i=1}^{d} \int_{\mathbb{R}^d} f v_i \partial_{v_i} (\psi_{\mu} * \log(f * \psi_{\mu})) \, dv - \sum_{i=1}^{d} \int_{\mathbb{R}^d} f v_i (v_i - u_i) \, dv. \tag{48}
\end{align*}
$$

As for the original Fokker-Planck operator, the second term is equal to $-dT^f \rho$ after an integration by part. Regarding the first term, we have

$$
\begin{align*}
\sum_{i=1}^{d} \int_{\mathbb{R}^d} f v_i \partial_{v_i} (\psi_\mu * \log(f * \psi_\mu)) \, dv &= \sum_{i=1}^{d} \int_{\mathbb{R}^d} f v_i \left( \psi_\mu * \frac{\partial_{v_i} f * \psi_\mu}{f * \psi_\mu} \right) \, dv = \sum_{i=1}^{d} \int_{\mathbb{R}^d} \left( (f v_i) * \psi_\mu \right) \frac{\partial_{v_i} f * \psi_\mu}{f * \psi_\mu} \, dv \\
&= \sum_{i=1}^{d} \int_{\mathbb{R}^d} v_i (f * \psi_\mu) \frac{\partial_{v_i} f * \psi_\mu}{f * \psi_\mu} \, dv - \sum_{i=1}^{d} \int_{\mathbb{R}^d} (f * (v_i \psi_\mu)) \frac{\partial_{v_i} f * \psi_\mu}{f * \psi_\mu} \, dv \\
&= \sum_{i=1}^{d} \int_{\mathbb{R}^d} v_i \partial_{v_i} f * \psi_\mu \, dv + \sum_{i=1}^{d} \int_{\mathbb{R}^d} (f * (\mu \partial_{v_i} \psi_\mu)) \frac{\partial_{v_i} f * \psi_\mu}{f * \psi_\mu} \, dv = -d\rho + \mu \int_{\mathbb{R}^d} \frac{|\nabla_v f * \psi_\mu|^2}{f * \psi_\mu} \, dv \\
&= -d\rho + \mu \int_{\mathbb{R}^d} (f * \psi_\mu) |\nabla_v \log(f * \psi_\mu)|^2 \, dv = -d\rho - \frac{\mu}{T^{\mathcal{M}}} (S, S) (f * \psi_\mu) = -d\rho - \mu D(f),
\end{align*}
$$

where we used the definition of the bracket (17), the entropy (19) and the definition of $D(f)$ in the last line. Gathering the two terms in (48) finally gives (26).

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

$$ \nabla_v \psi_\mu * \log(f * \psi_\mu / \mathcal{M}) = 0, $$

which is equivalent to

$$ \nabla_v \psi_\mu * \log(f * \psi_\mu) = -\frac{v-u}{T^{\mathcal{M}}}. $$

Integrating with respect to $v$ gives

$$ T^{\mathcal{M}} \psi_{\mu} * \log(f * \psi_{\mu}) = -\frac{|v|^2}{2} + u \cdot v + C, \quad (49) $$

which is more complicated to solve than in the original case. However, using Fourier techniques and following the lines in [8], it is possible to compute the solution $f$ of (49). Let denote $\mathfrak{F}[f]$ the Fourier transform in $v$ of $f$

$$ \mathfrak{F}[f](\xi) = \int_{\mathbb{R}^d} f(v)e^{-iv \cdot \xi} \, dv. $$

Using the properties of the Fourier transform, we can isolate $f * \psi_\mu$ in (49)

$$ f * \psi_{\mu} = \exp\left(\mathfrak{F}^{-1}\left[\mathfrak{F}\left[\frac{1}{T^{\mathcal{M}}}\left(-\frac{|v|^2}{2} + u \cdot v + C\right)\right] / \mathfrak{F}[\psi_{\mu}]\right]\right). $$

We compute now the right hand side step by step. First, applying $\mathfrak{F}$ to (49) gives

$$ T^{\mathcal{M}} \mathfrak{F}[\log(f * \psi_{\mu})] = \frac{1}{\mathfrak{F}[\psi_{\mu}]} \mathfrak{F}\left[-\frac{|v|^2}{2} + u \cdot v + C\right] = \frac{(2\pi\mu)^{d/2}}{\psi_{1/\mu}} \left(\frac{\Delta}{2} + iu \cdot \nabla + C\right) \delta_{\xi=0}. $$
Applying the inverse Fourier transform $\mathcal{F}^{-1}$ to this latter equation gives

$$ 
\begin{aligned}
T^{\mathcal{M}} \log(f * \psi_{\mu}) &= (2\pi\mu)^{d/2} \left(\frac{1}{2\pi}\right)^d \int_{\mathbb{R}^d} \left(\frac{2\pi}{\mu}\right)^{d/2} e^{\frac{\mu|\xi|^2}{2}} e^{iv \cdot \xi} \left(\frac{\Delta}{2} + iu \cdot \nabla + C\right) \delta_{\xi=0} \, d\xi \\
&= \int_{\mathbb{R}^d} \left(\frac{\Delta}{2} - iu \cdot \nabla + C\right) e^{\frac{\mu|\xi|^2}{2}} e^{iv \cdot \xi} \delta_{\xi=0} \, d\xi \\
&= \int_{\mathbb{R}^d} \left(\frac{1}{2}(d\mu + \mu^2|\xi|^2 - |v|^2 + 2\mu i\xi \cdot v) + u \cdot (v - i\mu\xi) + C\right) e^{\frac{\mu|\xi|^2}{2}} e^{iv \cdot \xi} \delta_{\xi=0} \, d\xi \\
&= \frac{1}{2}(d\mu - |v|^2) + u \cdot v + C = -\frac{1}{2}|v - u|^2 + \left(\frac{d\mu}{2} + C + \frac{|u|^2}{2}\right).
\end{aligned}
$$

Taking the exponential on both sides leads to $f * \psi_{\mu} = \tilde{C} \exp\left(-\frac{|v-u|^2}{2T^{\mathcal{M}}}\right)$ where $\tilde{C} = \exp\left(\frac{1}{T^{\mathcal{M}}}\left(\frac{d\mu}{2} + C + \frac{|u|^2}{2}\right)\right)$ is determined by integrating the equality

$$ \tilde{C}(2\pi T^{\mathcal{M}})^{d/2} = \int_{\mathbb{R}^d} \tilde{C} e^{-\frac{|v-u|^2}{2T^{\mathcal{M}}}} \, dv = \int_{\mathbb{R}^d} f * \psi_{\mu} \, dv = \int_{\mathbb{R}^d} f \, dv = \rho. $$

We finally find $f * \psi_{\mu} = \mathcal{M}$. If we continue computation to isolate $f$ using the same technic, we find that $f$ is a Maxwellian characterized by $\rho$, $u$ and $T^{\mathcal{M}} - \mu$.

## References

[1] N. Ben Abdallah and R. El Hajj. Diffusion and guiding center approximation for particle transport in strong magnetic fields. *Kinetic and Related Models*, 1(3):331–354, Sept. 2008.

[2] C. Birdsall and A. Langdon. *Plasma Physics via Computer Simulation*. Series in Plasma Physics and Fluid Dynamics. Taylor & Francis, 2004.

[3] T. Blanc, M. Bostan, and F. Boyer. Asymptotic analysis of parabolic equations with stiff transport terms by a multi-scale approach. *Discrete and Continuous Dynamical Systems - Series A*, 37(9):4637–4676, Sept. 2017.

[4] M. Bostan and C. Caldini Queiros. Finite Larmor radius approximation for collisional magnetic confinement. Part I: The linear Boltzmann equation. *Quarterly journal of pure and applied mathematics*, LXXII(2):323–345, 2014.

[5] M. Bostan and C. Caldini Queiros. Finite Larmor radius approximation for collisional magnetic confinement. Part II: The Fokker-Planck-Landau equation. *Quarterly journal of pure and applied mathematics*, LXXII(2):513–548, 2014.

[6] M. Bostan and A. Finot. Finite larmor radius regime: Collisional setting and fluid models. *Communications in Contemporary Mathematics*, 22(06):1950047, 2020.

[7] J. A. Carrillo, K. Craig, and F. S. Patacchini. A blob method for diffusion. *Calc. Var. Partial Differential Equations*, 58, 2019.

[8] J. A. Carrillo, J. Hu, L. Wang, and J. Wu. A particle method for the homogeneous landau equation. *Journal of Computational Physics*: X, 7:100066, 2020.

[9] P. Chartier, N. Crouseilles, M. Lemou, and F. Méhats. Uniformly accurate numerical schemes for highly oscillatory Klein-Gordon and nonlinear Schrödinger equations. *Numerische Mathematik*, 129(2):211–250, 2015.

[10] P. Chartier, N. Crouseilles, M. Lemou, F. Méhats, and X. Zhao. Uniformly accurate methods for Vlasov equations with non-homogeneous strong magnetic field. *Mathematics of Computation*, 88(320):2697–2736, 2019.
[11] P. Chartier, N. Crouseilles, M. Lemou, F. Méhats, and X. Zhao. Uniformly accurate methods for three dimensional vlasov equations under strong magnetic field with varying direction. *SIAM Journal on Scientific Computing*, 42(2):B520-B547, 2020.

[12] P. Chartier, N. Crouseilles, and X. Zhao. Numerical methods for the two-dimensional Vlasov-Poisson equation in the finite Larmor radius approximation regime. *Journal of Computational Physics*, 375:619-640, Dec. 2018.

[13] P. Chartier, M. Lemou, F. Méhats, and X. Zhao. Derivative-free high-order uniformly accurate schemes for highly oscillatory systems. *IMA Journal of Numerical Analysis*, 42(2):1623-1644, 05 2021.

[14] P. Chartier, A. Murua, and J. Sanz-Serna. High order averaging, formal series and numerical integrations i: B-series. *FOCM*, 10(6), 2010.

[15] P. Chartier, A. Murua, and J. Sanz-Serna. A formal series approach to averaging: exponentially small error estimates. *Discrete and continuous dynamical systems*, 32(9), 2012.

[16] A. Crestetto, N. Crouseilles, and D. Prel. Multiscale numerical schemes for the collisional vlasov equation in the finite larmor radius approximation regime. *Multiscale Modeling & Simulation*, 21(3):1210-1236, 2023.

[17] N. Crouseilles, M. Lemou, and F. Méhats. Asymptotic Preserving schemes for highly oscillatory Vlasov-Poisson equations. *Journal of Computational Physics*, 248:287-308, 2013.

[18] N. Crouseilles, M. Lemou, F. Méhats, and X. Zhao. Uniformly accurate forward semi-Lagrangian methods for highly oscillatory Vlasov-Poisson equations. *Multiscale Modeling and Simulation: A SIAM Interdisciplinary Journal*, 15(2):723-744, 2017.

[19] N. Crouseilles, M. Lemou, F. Méhats, and X. Zhao. Uniformly accurate Particle-In-Cell method for the long time two-dimensional Vlasov-Poisson equation with strong magnetic field. *Journal of Computational Physics*, 346:172-190, 2017.

[20] P. Degond and S. Mas-Gallic. The weighted particle method for convection-diffusion equations. I. The case of an isotropic viscosity. *Mathematics of Computation*, 93(345):485-507, 1989.

[21] P. Degond and F.-J. Mustieles. A deterministic approximation of diffusion equations using particles. *SIAM Journal on Scientific and Statistical Computing*, 11(2):293-310, 1990.

[22] P. Degond and P. Raviart. The paraxial approximation of the vlasov-maxwell equations. *Math. Models Methods Appl. Sci.*, 3:513-562, 1993.

[23] D. H. E. Dubin, J. A. Krommes, C. Oberman, and W. Lee. Nonlinear gyrokinetic equations. *Physics of Fluids*, 26(12):3524-3535, 1983.

[24] E. Frénod, F. Salvarani, and E. Sonnendrücker. Long time simulation of a beam in a periodic focusing channel via a two-scale pic-method. *Mathematical Models and Methods in Applied Sciences*, 19(02):175-197, feb 2009.

[25] E. Frénod and E. Sonnendrücker. Long time behavior of the two-dimensional vlasov equation with a strong external magnetic field. *Mathematical Models and Methods in Applied Sciences*, 10(4):539-553, 2000.

[26] M. Herda and M. Rodrigues. Large-time behavior of solutions to vlasov-poisson-fokker-planck equations: From evanescent collisions to diffusive limit. *Journal of Statistical Physics*, 170(5):895-931, 2018.

[27] M. Herda and M. Rodrigues. Anisotropic boltzmann-gibbs dynamics of strongly magnetized vlasov-fokker-planck equations. *Kinetic and Related Models*, 12(3):593-636, 2019.

[28] E. Hirvijoki. Structure-preserving marker-particle discretizations of coulomb collisions for particle-in-cell codes. *Plasma Physics and Controlled Fusion*, 63(4):044003, mar 2021.

[29] M. Hochbruck and A. Ostermann. Exponential integrators. *Acta Numer.*, 19:209-286, may 2010.

[30] S. Jeyakumar, M. Kraus, M. Hole, and D. Pfefferlé. A structure-preserving particle discretisation for the lenard-bernstein operator, arXiv, 2023.

[31] G. Knorr. Asymptotic state of the finite-larmor-radius guiding-centre plasma. *Journal of Plasma Physics*, 41(1):157-170, 1989.
[32] M. Krauss and E. Hirvijoki. Metriplectic integrators for the landau collision operator. *Physics of Plasmas*, 24(10):102311, mar 2017.

[33] W. Lee. Gyrokinetic particle simulation model. *Journal of Computational Physics*, 72(1):243-269, 1987.

[34] A. Mouton. Two-scale semi-lagrangian simulation of a charged particle beam in a periodic focusing channel. *Kinetic and Related Models*, 2(2):251-274, 2009.

[35] E. Sonnendrücker. *Numerical Methods for the Vlasov-Maxwell equations*. Lecture notes, 2015.

[36] F. Zonta, J. V. Pusztay, and E. Hirvijoki. Multispecies structure-preserving particle discretization of the landau collision operator. *Physics of Plasmas*, 29(12):123906, dec 2022.