# Conservative stabilized Runge-Kutta methods for the Vlasov-Fokker-Planck equation

Nicolas Crouseilles, Ibrahim Almuslimani

► To cite this version:

Nicolas Crouseilles, Ibrahim Almuslimani. Conservative stabilized Runge-Kutta methods for the Vlasov-Fokker-Planck equation. Journal of Computational Physics, 2023, 488, pp.1-34. 10.1016/j.jcp.2023.112241. hal-03911417

HAL Id: hal-03911417

https://hal.science/hal-03911417

Submitted on 22 Dec 2022

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 Attribution 4.0 International License
Conservative stabilized Runge-Kutta methods for the
Vlasov-Fokker-Planck equation

Ibrahim Almuslimaniª and Nicolas Crouseillesᵇ

ªUniv Rennes, Inria (MINGuS) & IRMAR UMR 6625, France

ᵇUniv Rennes, Inria (MINGuS) & IRMAR UMR 6625 & ENS Rennes, France

Contents

<table><tr><td>1</td><td>Introduction</td><td>2</td></tr><tr><td>2</td><td>Approximation of the Fokker-Planck operator</td><td>4</td></tr><tr><td>2.1</td><td>Second order numerical scheme</td><td>4</td></tr><tr><td>2.2</td><td>Fourth order numerical scheme</td><td>6</td></tr><tr><td>3</td><td>Explicit stabilized Runge-Kutta methods</td><td>9</td></tr><tr><td>3.1</td><td>Solving the velocity part using RKC methods</td><td>10</td></tr><tr><td>3.2</td><td>Damping.</td><td>12</td></tr><tr><td>3.3</td><td>Adaptive time stepping strategy.</td><td>12</td></tr><tr><td>4</td><td>Coupling with the Vlasov equation</td><td>13</td></tr><tr><td>4.1</td><td>Classical splitting.</td><td>14</td></tr><tr><td>4.2</td><td>Non-classical splitting.</td><td>14</td></tr><tr><td>5</td><td>Numerical results</td><td>16</td></tr><tr><td>5.1</td><td>Homogeneous case</td><td>16</td></tr><tr><td>5.1.1</td><td>CFL condition and second order in time</td><td>16</td></tr><tr><td>5.1.2</td><td>Fourth order approximation</td><td>18</td></tr><tr><td>5.1.3</td><td>Adaptive time stepping and conservation.</td><td>18</td></tr><tr><td>5.2</td><td>Non-homogeneous case 1dx-1dv</td><td>21</td></tr><tr><td>5.3</td><td>Non-homogeneous case 1dx-2dv</td><td>23</td></tr><tr><td>6</td><td>Conclusion</td><td>27</td></tr></table>
Abstract

In this work, we aim at constructing numerical schemes, that are as efficient as possible in terms of cost and conservation of invariants, for the Vlasov–Fokker–Planck system coupled with Poisson or Ampère equation. Splitting methods are used where the linear terms in space are treated by spectral or semi-Lagrangian methods and the nonlinear diffusion in velocity in the collision operator is treated using a stabilized Runge–Kutta–Chebyshev (RKC) integrator, a powerful alternative of implicit schemes. The new schemes are shown to exactly preserve mass and momentum. The conservation of total energy is obtained using a suitable approximation of the electric field. An H-theorem is proved in the semi-discrete case, while the entropy decay is illustrated numerically for the fully discretized problem. Numerical experiments that include investigation of Landau damping phenomenon and bump-on-tail instability are performed to illustrate the efficiency of the new schemes.

# 1 Introduction

We consider the following model satisfied by $f(t, x, v)$ with the time $t \ge 0$, the space variable $x \in \mathbb{R}^{d_x}$ and the velocity variable $v \in \mathbb{R}^{d_v}$

$$ \partial_t f + v \cdot \nabla_x f + E \cdot \nabla_v f = \nabla_v \cdot ((v - u_f)f + T_f \nabla_v f) =: Q(f), \quad (1) $$

where the electric field $E(t, x) = -\nabla_x \phi(t, x)$ is obtained from the Poisson equation $\Delta_x \phi(t, x) = \int_{\mathbb{R}^{d_v}} f(t, x, v)dv - 1$, the mean velocity $u_f(t, x)$ and temperature $T_f(t, x)$ are defined by

$$ n_f u_f(t, x) = \int_{\mathbb{R}^{d_v}} v f(t, x, v) dv, $$

$$ n_f T_f(t, x) = \frac{1}{d} \int_{\mathbb{R}^{d_v}} |v - u_f(t, x)|^2 f(t, x, v) dv $$

$$ \text{with the density} \qquad n_f(t, x) = \int_{\mathbb{R}^{d_v}} f(t, x, v) dv. \qquad (2) $$

The goal of this work is to construct efficient numerical grid based schemes (ie using a grid of the phase space) for the numerical approximation of (1). Even if these methods are costly, their high accuracy enables a very good description of the phase space even in low density regions. Moreover, stability conditions relating space or velocity mesh to the time step usually lead to high CPU time, in particular when dissipative collision operators are considered. To overcome this drawback, we aim at using stabilized Runge-Kutta-Chebyshev (RKC) methods for the numerical simulations of collisional Vlasov equations (1).
RKC methods are explicit methods with an extended stability domain along the real negative axis which are very well suited for the time integration of problems whose eigenvalues lie in a narrow strip around the negative real axis, which is typically the case for parabolic problems. Indeed, when one is interested in high dimensional or/and nonlinear parabolic problems, implicit method require to solve large or/and nonlinear systems which may become very costly from a computational point of view. Inversely, explicit methods do not have to invert complex systems, but at the cost of a severe condition on the time step. Additionally, when the system comes from the discretization of a PDE which enjoys suitable properties, explicit methods may be preferred since the properties can be easily extended to the fully discrete case, and implicit methods may break the properties of the underlying numerical schemes. Hence, RKC methods turn out to be a good compromise since they are explicit and their domain of stability has a length proportional to $s^2$ with s is the number of stages, which allow for the use of large time steps (see for instance [23, 24, 22]).

Compared with high-order RK schemes, stabilized RKC schemes require few function calls, and have an extended stability region along the negative real axis, which covers the eigenvalues of the dissipative collision term more efficiently and allows larger time steps. Even if the error is of second order, the global error is relatively small when the full system is considered. In addition, to optimize time stepping, the second order RKC scheme is coupled with adaptive time stepping strategy: the time step is automatically chosen according to a local error estimate, and the number of stages is dynamically adjusted.

In the context of Vlasov equation, solving large systems is not acceptable due to the high dimensionality of the equations. Standard explicit methods, due to the fact that they require small time steps, also lead to a large number of iterations and then to a large cost. In plasma physics, even if binary collisions are not the main effect compared to the mean field interaction, they are not negligible and have been considered in the non homogeneous case using Eulerian methods in [8, 7, 9]. Collisions are currently taken into account for example in gyrokinetic simulations but in this context, a more simple choice to treat collisions consistently is to use the Fokker-Planck type operators which leads to a second order parabolic equation (see [2, 17, 6, 10]). As such, they have a moderate levels of stiffness, their spectrum belongs to the negative half space and is contained in a strip around the negative real axis.

For such collision operators, the use of implicit methods can be a strong obstacle due to cost and memory aspects but explicit methods may lead to very stringent stability constraints that also increase the computational cost. Hence, stabilized Runge-Kutta methods turn out to be an interesting alternative since in addition the conservation properties of the discretized collision operator are simply transferred to the fully discrete case. Let us mention that these methods have been already used in [8] for Vlasov-Poisson-Landau and in [2, 17, 6] for collisional gyrokinetic case. However in these works, the time integrator was not the main point and is therefore briefly discussed and one of the objectives of this work is to explain why these time integrators are a viable way for the weakly collisional plasma simulations.

The rest of the paper is organized as follows: first the velocity discretization of the Fokker-Planck operator is presented; second, the explicit stabilized Runge-Kutta methods are discussed
and in a third part, the coupling between the Vlasov and collision part is presented. Finally, a bunch of numerical results are presented in different configurations to highlight the advantages of the approach.

## 2 Approximation of the Fokker-Planck operator

Here, we recall the Fokker-Planck collision operator and its properties (conservation and H-theorem) at the continuous level. Then, we present the velocity discretization and the properties at the semi-discrete level (discrete in velocity and continuous in time). Extensions to higher order and higher dimensions are provided.

The collision operator $Q(f)$ can be rewritten in the standard form, log form or $L^2$ form

$$Q(f) = \nabla_v \cdot ((v - u_f)f + T_f \nabla_v f) = \nabla_v \cdot \left( T_f f \nabla_v \log \left( \frac{f}{M_f} \right) \right) = \nabla_v \cdot \left( T M_f \nabla_v \left( \frac{f}{M_f} \right) \right), \quad (3)$$

with $n_f = \int f dv, n_f u_f = \int v f dv, dn_f T_f = \int |v - u_f|^2 f dv$ and $M_f$ denotes the Maxwellian associated to $f$

$$M_f(v) = \frac{1}{(2\pi T_f)^{d/2}} \exp \left( - \frac{|v - u_f|^2}{2T_f} \right).$$

The operator $Q(f)$ enjoys the following properties

$$\int (1, v, |v|^2) Q(f) dv = 0, \quad \int \log(f) Q(f) dv \le 0, \quad \int Q(f) \frac{f}{M_f} dv \le 0,$$

which is important to preserve at the numerical level. However, since the Fokker-Planck operator has been widely studied, there exists some suitable discretizations in the literature (see [4]...) which ensure some of these properties. In particular, in [4], the authors made the link between the different formulations (3) of the Fokker-Planck operator (standard form, log form or $L^2$ form in (3)). According to the chosen form, the discretization will enjoy conservation properties or entropy decreasing property. In [4], the entropic averaging is proposed to make the link between the standard and log forms. However, using the log form may impose strong stability constraints on the numerical scheme to ensure the positivity of the solution. Thus, we follow the same path as in [4] to make the link between the standard form (well suited for the conservations) and the $L^2$ form (well adapted for the entropy).

### 2.1 Second order numerical scheme

To present the numerical scheme, we restrict ourselves to the case $d = 1$ with $v_j = j \Delta v$ ($j \in \mathbb{Z}$ and $\Delta v > 0$ the velocity mesh). We also introduce the discrete unknown $f_j \approx f(v_j)$. We look for a discrete function $\tilde{M}_{j+1/2}$ (which is supposed to approximate a Maxwellian) such that the
discretization of the standard form $\partial_v((v - u_f)f + T_f \partial_v f)$ and the discretization of the $L^2$ form $\partial_v(T_f M_f \partial_v (\frac{f}{M_f}))$ are equivalent. For the standard form, we consider the following discretization

$$Q_j = \frac{1}{\Delta v} [f_{j+1/2}(v_{j+1/2} - \tilde{u}) - f_{j-1/2}(v_{j-1/2} - \tilde{u})] + \frac{\tilde{T}}{\Delta v^2} (f_{j+1} - 2f_j + f_{j-1}), \quad (4)$$

with

$$f_{j+1/2} = \frac{f_{j+1} + f_j}{2}, \quad \text{and} \quad v_{j+1/2} = \frac{v_j + v_{j+1}}{2}, \qquad (5)$$

whereas for the $L^2$ form, we consider the following discretization

$$Q_j = \frac{\tilde{T}}{\Delta v^2} \left\{ \tilde{M}_{j+1/2} \left[ \left( \frac{f}{M} \right)_{j+1} - \left( \frac{f}{M} \right)_j \right] - \tilde{M}_{j-1/2} \left[ \left( \frac{f}{M} \right)_j - \left( \frac{f}{M} \right)_{j-1} \right] \right\}. \quad (6)$$

We need now to define the different quantities involved in the numerical schemes: first, for the macroscopic quantities $\tilde{u}$ and $\tilde{T}$, we choose (following [4])

$$n = \sum_j f_j \Delta v, \quad n \tilde{u} = \sum_j v_{j+1/2} f_{j+1/2} \Delta v \quad \text{and} \quad n \tilde{T} = \sum_j (v_{j+1/2} - \tilde{u})^2 f_{j+1/2} \Delta v. \quad (7)$$

and second the discrete Maxwellian $M_j$ is defined as

$$M_j = M(v_j) = \frac{n}{\sqrt{2\pi\tilde{T}}} \exp\left(-|v_j - \tilde{u}|^2/(2\tilde{T})\right). \quad (8)$$

As shown in [4] (the calculation are recalled in Appendix A)), the first formulation (4) enables to prove easily the conservations of mass, momentum and energy whereas the second one (6) enables to prove easily the entropy inequality. Hence, the goal is to find an expression of $\tilde{M}_{j+1/2}$ such that the two expressions are equivalent.

We introduce the flux

$$\mathcal{F}_{j+1/2} = f_{j+1/2}(v_{j+1/2} - \tilde{u}) + \frac{\tilde{T}}{\Delta v} (f_{j+1} - f_j), \quad (9)$$

so that (4) can be expressed as $Q_j = (\mathcal{F}_{j+1/2} - \mathcal{F}_{j-1/2})/\Delta v$. Similarly, defining

$$\mathcal{G}_{j+1/2} = \frac{\tilde{T}}{\Delta v} \tilde{M}_{j+1/2} \left[ \left( \frac{f}{M} \right)_{j+1} - \left( \frac{f}{M} \right)_j \right], \quad (10)$$

enables to express (6) as $Q_j = (\mathcal{G}_{j+1/2} - \mathcal{G}_{j-1/2})/\Delta v$. Hence, we identify $\mathcal{F}_{j+1/2}$ and $\mathcal{G}_{j+1/2}$ to define $\tilde{M}_{j+1/2}$

$$\mathcal{G}_{j+1/2} \equiv \frac{\tilde{T}}{\Delta v} \tilde{M}_{j+1/2} \left( \frac{f_{j+1} M_j - f_j M_{j+1}}{M_j M_{j+1}} \right) = f_{j+1/2} (v_{j+1/2} - \tilde{u}) + \frac{\tilde{T}}{\Delta v} (f_{j+1} - f_j) \equiv \mathcal{F}_{j+1/2}$$
which enables to define $\tilde{M}_{j+1/2}$ as

$$ \tilde{M}_{j+1/2} = \frac{\Delta v M_j M_{j+1} f_{j+1/2} (v_{j+1/2} - \tilde{u})}{\tilde{T} (f_{j+1} M_j - f_j M_{j+1})} + \frac{M_j M_{j+1} (f_{j+1} - f_j)}{f_{j+1} M_j - f_j M_{j+1}}. \quad (11) $$

Hence, the standard discretization (4) is equivalent to the discretization (6) where $\tilde{M}_{j+1/2}$ is defined by (11) which enables us to prove the following Proposition.

**Proposition 1.** The discretization (4)-(5)-(7) of the Fokker-Planck operator $Q$ (3) satisfies

$$ \sum_{j \in \mathbb{Z}} Q_j \begin{pmatrix} 1 \\ v_j \\ v_j^2/2 \end{pmatrix} \Delta v = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} \quad \text{and} \quad \sum_{j \in \mathbb{Z}} Q_j \frac{f_j}{M_j} \Delta v \le 0, $$

which means that the discrete entropy $\mathcal{E}(t) = \Delta v \sum_{j \in \mathbb{Z}} (f_j^2)/M_j$ satisfies $\frac{d}{dt} \mathcal{E}(t) \le 0$.

The proof is given in Appendix A.

## 2.2 Fourth order numerical scheme

Based on the discretization (4) of $Q$, the goal of this part is to derive a higher order velocity discretization for $Q$. First, integrating the standard form of $Q$ on $[v_{j-1/2}, v_{j+1/2}]$ leads to

$$ \begin{aligned} \bar{Q}_j &\equiv \frac{1}{\Delta v} \int_{v_{j-1/2}}^{v_{j+1/2}} Q(f)(v)dv \\ &= \frac{1}{\Delta v} \left( (v_{j+1/2} - u_f)f(v_{j+1/2}) - (v_{j-1/2} - u_f)f(v_{j-1/2}) + T_f(\partial_v f(v_{j+1/2}) - \partial_v f(v_{j-1/2})) \right), \end{aligned} $$

where the notations $u_f, T_f$ are given by (2). Using the following fourth order finite difference formula

$$ f(v_{j+1/2}) \approx \check{f}_{j+1/2} = -\frac{1}{16}f_{j-1} + \frac{9}{16}f_j + \frac{9}{16}f_{j+1} - \frac{1}{16}f_{j+2}, \quad (12) $$

$$ (\partial_v f)(v_{j+1/2}) \approx \frac{f_{j-1} - 27f_j + 27f_{j+1} - f_{j+2}}{24\Delta v}, $$

we get the following fourth order finite volume approximation $\bar{Q}_j^{[4]}$ of $\bar{Q}_j$

$$ \bar{Q}_j^{[4]} = \frac{(v_{j+1/2} - \check{u})\check{f}_{j+1/2} - (v_{j-1/2} - \check{u})\check{f}_{j-1/2}}{\Delta v} + \check{T} \left( \frac{-f_{j-2} + 28f_{j-1} - 54f_j + 28f_{j+1} - f_{j+2}}{24\Delta v^2} \right), \quad (13) $$

with $\check{f}_{j+1/2}$ given by (12) and the definitions of $\check{u}$ and $\check{T}$ have to be adapted in this case:

$$ n = \Delta v \sum_j f_j, \quad \check{u} = \Delta v \sum_j v_{j+1/2} \check{f}_{j+1/2}, \quad \text{and} \quad n\check{T} = \Delta v \sum_j (v_{j+1/2} - \check{u})^2 \check{f}_{j+1/2}. \quad (14) $$
To get a pointwize approximation (ie of $Q(f)(v_j)$), we will use the following relation between cell average and the pointwize value. For all smooth function $g$, we get

$$ \bar{g}_j = \frac{1}{\Delta v} \int_{v_{j-1/2}}^{v_{j+1/2}} g(v)dv = g(v_j) + \frac{\Delta v^2}{24} g''(v_j) + \mathcal{O}(\Delta v^4), \quad (15) $$

from which we deduce the following relation

$$ g(v_j) = \bar{g}_j - \frac{\Delta v^2 g(v_{j+1}) - 2g(v_j) + g(v_{j-1})}{24 \Delta v^2} + \mathcal{O}(\Delta v^4). $$

We then deduce a fourth order approximation $Q_j^{[4]}$ of $Q(f)(v_j)$ gathering the previous ingredients

$$ Q_j^{[4]} = \bar{Q}_j^{[4]} - \frac{1}{24} [Q_{j+1} - Q_j + Q_{j-1}] \quad \text{with } Q_j \text{ given by (4)-(5)-(7)}, \quad (16) $$

where $\bar{Q}_j^{[4]}$ is given by (13).

In the following Proposition, we gather the properties of this new discretization.

**Proposition 2.** The discretization (16)-(12)-(13)-(14) of the Fokker-Planck operator $Q$ (3) satisfies

$$ Q(f)(v_j) - Q_j^{[4]} = \mathcal{O}(\Delta v^4), $$

and

$$ \sum_j Q_j^{[4]} \begin{pmatrix} 1 \\ v_j \\ v_j^2/2 \end{pmatrix} \Delta v = \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix}. $$

*Proof.* First, we have from (15)

$$ \begin{aligned} \bar{Q}_j &= Q(f)(v_j) + \frac{\Delta v^2}{24}(\partial_v^2 Q(f))(v_j) + \mathcal{O}(\Delta v^4) \\ &= Q(f)(v_j) + \frac{1}{24}(Q(f)(v_{j+1}) - 2Q(f)(v_j) + Q(f)(v_{j-1})) + \mathcal{O}(\Delta v^4), \end{aligned} $$

and we also have $\bar{Q}_j = \bar{Q}_j^{[4]} + \mathcal{O}(\Delta v^4)$, from which we deduce

$$ Q(f)(v_j) = \bar{Q}_j^{[4]} - \frac{1}{24}(Q(f)(v_{j+1}) - 2Q(f)(v_j) + Q(f)(v_{j-1})) + \mathcal{O}(\Delta v^4). $$

It remains to check that $Q(f)(v_{j+1}) - 2Q(f)(v_j) + Q(f)(v_{j-1}) - (Q_{j+1} - Q_j + Q_{j-1}) = \mathcal{O}(\Delta v^4)$, where $Q_j$ is given by (4)-(5)-(7). First, we look at the consistency error $Q(f)(v_j) - Q_j$ by performing Taylor expansions of $f(v_{j\pm 1})$ in $Q_j$

$$ \begin{aligned} Q_j ={}& Q(f)(v_j) + (u_f - \tilde{u})(\partial_v f)(v_j) + (T_f - \tilde{T})(\partial_v^2 f)(v_j) \\ & + \Delta v^2 \left( \frac{1}{4}(\partial_v^2 f)(v_j) + \frac{1}{24}(v_j - \tilde{u})\partial_v^3 f(v_j) + \frac{\tilde{T}}{24}(\partial_v^4 f)(v_j) \right) + \mathcal{O}(\Delta v^4). \end{aligned} $$
Thus, we have

$$ 
\begin{aligned}
Q_{j+1} - 2Q_j + Q_{j-1} &= Q(f)(v_{j+1}) - 2Q(f)(v_j) + Q(f)(v_{j-1}) \\
&\quad + (u_f - \tilde{u})\Delta v^2(\partial_v^3 f)(v_j) + (T_f - \tilde{T})\Delta v^2(\partial_v^4 f)(v_j) + \mathcal{O}(\Delta v^4) \\
&= \Delta v^2(\partial_v^2 Q(f))(v_j) + \mathcal{O}(\Delta v^4),
\end{aligned}
$$

since $u_f - \tilde{u} = T_f - \tilde{T} = \mathcal{O}(\Delta v^2)$ using second order quadrature. We deduce

$$ 
\begin{aligned}
Q(f)(v_j) - Q_j^{[4]} &= Q(f)(v_j) - \left[ \bar{Q}_j^{[4]} - \frac{1}{24}(Q_{j+1} - 2Q_j + Q_{j-1}) \right] \\
&= Q(f)(v_j) - \bar{Q}_j^{[4]} + \frac{\Delta v^2}{24}\partial_v^2 Q(v_j) + \mathcal{O}(\Delta v^4) \\
&= Q(f)(v_j) - \left[ \bar{Q}_j - \frac{\Delta v^2}{24}\partial_v^2 Q(v_j) \right] + \mathcal{O}(\Delta v^4) \\
&= \mathcal{O}(\Delta v^4).
\end{aligned}
$$

Next, we prove the conservation properties. We first prove that the conservation hold true for $\bar{Q}_j^{[4]}$. The flux formulation implies directly mass conservation. Next, we introduce the notation $(Df)_j = (-f_{j-2} + 28f_{j-1} - 54f_j + 28f_{j+1} - f_{j+2})/(24\Delta v^2)$ which is a fourth order approximation of $(\partial_v^2 f)(v_j)$. The momentum conservation follows from the property $\sum_j v_j (Df)_j \Delta v = 0$ which allows to follow the lines of the proof of Proposition 1. For the energy, we first check $\Delta v \sum_j (v_j^2/2)(Df)_j = n$ by two discrete integration by parts. Thus, we follow the lines of the proof of Proposition 1 to get $\sum_j \bar{Q}_j^{[4]}(1, v_j, v_j^2)\Delta v = (0, 0, 0)$.

Indeed, we have directly that $\sum_j Q_j^{[2]}(1, v_j, v_j^2)\Delta v = (0, 0, 0)$. It remains to check that $\sum_j [Q_{j+1} - 2Q_j + Q_{j-1}](1, v_j, v_j^2)\Delta v = (0, 0, 0)$. The first component is obvious using Proposition 1. For the momentum, we have

$$ 
\begin{aligned}
\sum_j Q_{j \pm 1} v_j \Delta v &= \sum_j Q_j v_{j \mp 1} \Delta v = \sum_j Q_j (v_j \mp \Delta v) \Delta v \\
&= \sum_j Q_j v_j \Delta v \mp \sum_j Q_j \Delta v^2 = 0,
\end{aligned}
$$

again from Proposition 1. For the energy, we proceed similarly to get $\sum_j Q_{j \pm 1} v_j^2 \Delta v = 0$ which enable to conclude the proof. $\square$

**Remark.** For the entropy, it is possible to follow the methodology used for the second order approach to construct (11). Indeed, we can derive the fluxes from (16)

$$ 
\begin{aligned}
\mathcal{F}_{j+1/2} &= (v_{j+1/2} - \tilde{u})\tilde{f}_{j+1/2} + \frac{\tilde{T}}{24\Delta v}(-f_{j+2} + 27f_{j+1} - 27f_j + f_{j-1}) \\
&\quad - \frac{1}{24}(\mathcal{F}_{j+3/2}^{[2]} - 2\mathcal{F}_{j+1/2}^{[2]} + \mathcal{F}_{j-1/2}^{[2]})
\end{aligned}
$$
with $\mathcal{F}_{j+1/2}^{[2]}$ the flux of the second order version given by (9), $n = \sum_j f_j \Delta v$ and

$$n\tilde{u} = \Delta v \sum_j v_{j+1/2} \tilde{f}_{j+1/2}, \quad n\tilde{T} = \Delta v \sum_j (v_{j+1/2} - \tilde{u})^2 \tilde{f}_{j+1/2},$$

$$n\tilde{u} = \Delta v \sum_j v_{j+1/2} f_{j+1/2}, \quad n\tilde{T} = \Delta v \sum_j (v_{j+1/2} - \tilde{u})^2 f_{j+1/2},$$

with $\tilde{f}_{j+1/2}$ is given by (12) and $f_{j+1/2} = (f_{j+1} + f_j)/2$. On the other side, we derive the fourth order fluxes for the $L^2$ form (with $g_j = (f/M)_j$)

$$\mathcal{G}_{j+1/2} = \frac{\tilde{T} \tilde{M}_{j+1/2}}{24 \Delta v} (g_{j-1} - 27g_j + 27g_{j+1} - g_{j+2}) - \frac{1}{24} (\mathcal{G}_{j+3/2}^{[2]} - 2\mathcal{G}_{j+1/2}^{[2]} + \mathcal{G}_{j-1/2}^{[2]})$$

with $\tilde{M}_{j+1/2}$ is to be computed and $\mathcal{G}_{j+1/2}^{[2]}$ the flux of the second order version given by (10). Identifying the fluxes $\mathcal{F}_{j+1/2}$ and $\mathcal{G}_{j+1/2}$ enables to derive $\tilde{M}_{j+1/2}$:

$$\tilde{M}_{j+1/2} = \frac{1}{\alpha_j} \left[ \mathcal{F}_{j+1/2} + \frac{1}{24} \left( \mathcal{G}_{j+3/2}^{[2]} - 2\mathcal{G}_{j+1/2}^{[2]} + \mathcal{G}_{j-1/2}^{[2]} \right) \right],$$

$$\text{with } \alpha_j = \frac{\tilde{T}}{24\Delta v} (g_{j-1} - 27g_j + 27g_{j+1} - g_{j+2}), \\ = \frac{\tilde{T}(f_{j-1}M_j M_{j+1} M_{j+2} - 27f_j M_{j-1} M_{j+1} M_{j+2} + 27f_{j+1} M_{j-1} M_j M_{j+2} - f_{j+2} M_{j-1} M_j M_{j+1})}{24M_{j-1}M_j M_{j+1}M_{j+2}\Delta v}$$

ensuring that the two forms are equivalent from a discrete point of view. Some calculations similar to the ones performed in the end of Appendix A enable to prove that $\tilde{M}_{j+1/2} = M(v_{j+1/2}) + \mathcal{O}(\Delta v^4)$.

### 3 Explicit stabilized Runge-Kutta methods

The space discretization of a diffusion dominated advection-diffusion PDE leads to a high dimensional system of ODEs of the form

$$\dot{y} = F(y), \quad y(0) = y_0, \tag{17}$$

where the eigenvalues of the Jacobian of the vector field $F$ are located in a narrow strip around the negative real axis. In order to study the stability of Runge-Kutta methods, it is common to consider the linear equation (see [22, Chap. IV.2])

$$\dot{y} = \lambda y, \quad y(0) = y_0,$$

where $\lambda \in \mathbb{C}$. Applying a RK method with time step $\Delta t$ to the above linear ODE leads to $y^n = R(\lambda \Delta t)^n y_0$ where $R(z)$ is called the stability function, and then the stability domain of the method is defined by

$$S = \{z \in \mathbb{C}; |R(z)| \le 1\}.$$
For stiff equations, the fact that the stability function of explicit methods is always a polynomial, which means bounded stability region, causes restrictions on the step size which turn out to be severe in the case of standard explicit schemes such as the forward Euler method. Implicit schemes are a good alternative for low dimensional problems, but their cost becomes prohibitive when solving high dimensional problems since they require the solution of high dimensional linear or/and nonlinear systems.

Explicit stabilized RKC methods (for Runge-Kutta-Chebyshev) have the advantages of both worlds. They do not require the solution of linear/nonlinear systems and, at the same time, they allow large step sizes thanks to their extended stability domains over the negative real axis. The literature of these methods is rich, and many schemes of this types where introduced to solve different types of stiff problems such as advection-diffusion-reaction equations [28, 30, 20], stochastic differential equations [29, 32, 33], and optimal control problems [31]. For a complete review on explicit stabilized RKC methods for stiff ODEs, we refer the reader to the review [23].

## 3.1 Solving the velocity part using RKC methods

Here, we detail how to solve numerically $\partial_t f = Q(f)$ or $\partial_t f + E\partial_v f = Q(f)$ depending on the followed splitting strategy (see Section 4). For the second equation, we combine the transport part $E\partial_v f$ with collisional part as

$$ \partial_t f = \tilde{Q}(f) := \partial_v((v - u - E)f + T\partial_v f). $$

we will present the algorithm for $\partial_t f = Q(f)$, which applies in the same way on $\partial_t f = Q(f)$. After a velocity discretization of $Q(f)$ following the methods presented in Section 2, it remains to integrate in time. Motivated by the fact that the eigenvalues of the Jacobian of $Q$ are localized in the negative half space (see Figure 1), we propose here to use a stabilized Runge-Kutta method which enables to satisfy several properties. Indeed, since it is an explicit method, it avoids a potentially costly implicit solver, moreover, the discrete properties of $Q$ can be transferred easily to the fully discrete case and of course it is not constrained by the severe parabolic CFL condition and enables to use large time steps.

Then, once discretized in velocity, we obtain

$$ \frac{df_j(t)}{dt} = Q_j, \quad \text{with } Q_j = Q(f_j(t)), \qquad (18) $$

which is amenable for stabilized RKC methods. The first order RKC method (RKC1) applied to
<figure><img src="image_2.png" /><figcaption>Figure 1: Eigenvalues of Q for $N_v = 512$, $T = 1.88$, $u = 0$, $\nu = 0.5$, and different values of $v_{max}$.</figcaption></figure>

(18) reads

$$ 
\begin{aligned}
K^{(0)} &= f_j^n \\
K^{(1)} &= K^{(0)} + \mu_1 \Delta t Q_j^n, \\
&\quad \text{for } \ell = 2, \dots, s \\
K^{(\ell)} &= \kappa_\ell K^{(\ell-2)} + \nu_\ell K^{(\ell-1)} + \mu_\ell \Delta t Q_j^{(\ell-1)}, \\
f_j^{n+1} &= K^{(s)},
\end{aligned} 
\quad (19) $$

where,

$$ \mu_1 = \frac{\omega_1}{\omega_0}, \quad \mu_\ell = \frac{2\omega_1 T_{\ell-1}(\omega_0)}{T_\ell(\omega_0)}, \quad \nu_\ell = \frac{2\omega_0 T_{\ell-1}(\omega_0)}{T_\ell(\omega_0)}, \quad \kappa_\ell = 1 - \nu_\ell, \quad \ell = 2, \dots, s. \quad (20) $$

for $\omega_0 = 1 + \frac{\eta}{s^2}$, $\omega_1 = \frac{T_s(\omega_0)}{T_s(\omega_0)}$, and $\eta = 0.05$. The second order RKC method (RKC2) reads

$$ 
\begin{aligned}
K^{(0)} &= f_j^n \\
K^{(1)} &= K^{(0)} + \hat{\mu}_1 \Delta t Q_j^n, \\
&\quad \text{for } \ell = 2, \dots, s \\
K^{(\ell)} &= (1 - \hat{\kappa}_\ell - \hat{\nu}_\ell) K^{(0)} + \hat{\kappa}_\ell K^{(\ell-2)} + \hat{\nu}_\ell K^{(\ell-1)} + \hat{\mu}_\ell \Delta t Q_j^{(\ell-1)} - a_{\ell-1} \hat{\mu}_\ell \Delta t Q_j^n, \\
f_j^{n+1} &= K^{(s)},
\end{aligned} 
$$

where

$$ \hat{\mu}_\ell = \frac{2b_\ell \omega_2}{b_{\ell-1}}, \quad \hat{\nu}_\ell = \frac{2b_\ell \omega_0}{b_{\ell-1}}, \quad \hat{\kappa}_\ell = -\frac{b_\ell}{b_{\ell-2}}, \quad b_\ell = \frac{T_\ell''(\omega_0)}{(T_\ell'(\omega_0))^2}, \quad a_\ell = 1 - b_\ell T_\ell(\omega_0), \quad (21) $$
for $\omega_0 = 1 + \frac{\eta}{s^2}$, $\omega_2 = \frac{T'_s(\omega_0)}{T'_s(\omega_0)}$, and $\eta = 0.15$. In the above equations, $Q_j^{(\ell)} = Q(f_j^{(\ell)})$, $s \in \mathbb{N}\setminus\{0, 1\}$ is the number of stages of the method. The stability region of the RKC schemes contains a strip around the interval $[-C_\eta s^2, 0]$ where, typically, $C_\eta = 1.93$ for RKC1 and $C_\eta = 0.65$ for RKC2. The dependence on $\eta$ will be clarified in the next subsection. For a given step size $\Delta t$, the number of stages $s$ is chosen such that the product $\lambda_{\max}\Delta t$ lies in the stability region, where $\lambda_{\max}$ is the eigenvalue with the maximum real part in absolute value.

**Proposition 3.** Consider the homogeneous equation $\partial_t f = \nu Q(f)$ with $Q$ given by (3). Suppose that we discretize $Q(f)$ in velocity using (4)-(5)-(7) (second order version) or using (16)-(12)-(13)-(14) (fourth order version), and RKC in time (19). Let us denote by $f_j^n$ the approximation of $f(t^n, v_j)$, then the resulting numerical scheme preserves mass, momentum, and kinetic energy

$$ \text{For all } n \in \mathbb{N}, \quad \sum_j \begin{pmatrix} 1 \\ v_j \\ v_j^2/2 \end{pmatrix} f_j^n = \sum_j \begin{pmatrix} 1 \\ v_j \\ v_j^2/2 \end{pmatrix} f_j^0. $$

*Proof.* The proof follows easily by induction on the internal stages of the RKC methods. $\square$

## 3.2 Damping

The coefficients of the above methods depend on a parameter usually denoted $\eta$ and called damping parameter. In the absence of damping, ($\eta = 0$), the stability region has singular points that correspond to the value $\pm 1$ of the stability function. In order to make the schemes robust with respect to possible small perturbations of the eigenvalues, coming from advection terms for example, a positive value of $\eta$ is needed. For pure diffusion or diffusion dominated advection-diffusion operators, $\eta$ is usually fixed to 0.05 for RKC1 and 0.15 for RKC2 (hence the values 1.93 and 0.65 of $C_\eta$ respectively). However, increasing $\eta$ gives the stability region more width in the imaginary direction which gives the scheme the ability to handle larger advection terms with reasonable time step sizes. Figure 2 illustrates this idea. Note that this gain in the imaginary direction comes at the cost of reducing the constant $C_\eta$ and then a reduction in the length of the real negative interval included in the stability domain. A good estimation of this constant is $C_\eta = (1 + \omega_0)/\omega_1$.

## 3.3 Adaptive time stepping strategy

We will use the automatic time step selection procedure adapted to the second order RKC method from [24]. The estimator used to estimate the local error of the method is the following

$$ \text{Est}_{n+1} = \frac{1}{15}[12(f^n - f^{n+1}) + 6\Delta t(Q^n + Q^{n+1})]. $$

In order to reach a given accuracy, the adaptive time step selection is done in the spirit of the procedure explained in [24]. For reference, see also [21, P. 167]. After choosing the suitable time
<figure><img src="image_3.png" /><figcaption>Figure 2: Stability domains and stability functions of the Chebyshev method (RKC1) for $s = 7$ and different damping values $\eta = 0, 0.05, 3.98$.</figcaption></figure>

step, the number of stages is chosen using the formula

$$s := \left[ \sqrt{\frac{\Delta t \lambda_{\max} + 1.5}{C_{\eta}}} + 0.5 \right],$$

where the brackets mean rounding to the nearest integer. Note that the automatic time step control takes care as well of possible need to reduce the step size due to large advection during the integration process. This large effect of the advection term might appear when the maximum value of the velocity $v_{\max}$ is quite large.

# 4 Coupling with the Vlasov equation

To solve the full Vlasov-Poisson-Fokker-Planck system, we use splitting techniques between the different operators following two different strategies. The first strategy is classical is relatively cheaper, while the second can be smartly manipulated to exactly conserve the total energy.
## 4.1 Classical splitting

Even if a three steps splitting (between transport in space, velocity and collision) can be easily obtained, the velocity transport part can be merged with the collision part so that the resulting splitting called SL-RKC is

* Solve the transport in space $\partial_t f + v\partial_x f = 0$ using a semi-Lagrangian method or Fourier technique.

* Solve $\partial_x E = 1 - \int_{\mathbb{R}} f dv$ using Fourier transform.

* Solve $\partial_t f + E\partial_v f = Q(f)$ using a stabilized Runge-Kutta-Chebyshev method (RKC).

The first two parts are classical and the last part can be rewritten as

$$ \partial_t f = \partial_v((v - u_f - E)f + T_f \partial_v f) = \tilde{Q}(f), $$

which has the same structure as $Q(f)$ so that the numerical schemes presented in Section 2 can be applied.

## 4.2 Non-classical splitting

The second strategy is motivated by total energy (kinetic plus electric) conservation by using ideas from [25, 26] and applied to Vlasov-Ampère instead of Vlasov-Poisson. The algorithm called SL-RK2-RKC is

* Solve the transport in space $\partial_t f + v\partial_x f = 0$ using a semi-Lagrangian method or Fourier technique.

* Instead of solving a Poisson equation for the electric field $E$, solve the corresponding Ampère equation $\partial_t E = -\int_{\mathbb{R}d_v} vfdv$.

* Solve $\partial_t f + E\partial_v f = 0$ using a suitable second order Runge-Kutta method and suitable approximation of $E$ that guarantees the exact conservation of total energy. Details will be given later in the current section.

* Solve $\partial_t f = Q(f)$ using RKC.

For this non-classical splitting, we will benefit from a smart idea introduced in the papers [25, 26] of adapting the electric field approximation in order to construct a numerical scheme that preserves exactly the total energy of the system. Here, we will follow [26] and use Ampère equation for the electric field instead of Poisson equation which will be used only once at the beginning to get an approximation of the initial value of $E$. In what follows, we will set $d_x = d_v = 1$.

Let $D_x f$ be a (at least) consistent finite difference approximation of $\partial_x f$, and $D_v f$ be a (at least) second order finite differences approximation of $\partial_v f$ (centered or upwind). Let $f_{i,j}^n$ denote an
approximation of $f(t^n, x_i, v_j)$ with $t^n = n\Delta t$, $x_i = i\Delta x$, $v_j = j\Delta v$ with $\Delta t, \Delta x, \Delta v$ the time, space and velocity step meshes. We consider the following discretization of the VFP equation

$$ 
\begin{aligned}
f_{ij}^{(1)} &= f_{ij}^{n}, \\
f_{ij}^{(2)} &= f_{ij}^{n} + \frac{\Delta t}{2} F_{ij}^{(1)}, \\
f_{ij}^{n+1} &= f_{ij}^{n} + \Delta t F_{ij}^{(2)},
\end{aligned} 
\qquad (22) $$

where,

$$ F_{ij}^{(k)} = -v_j D_x f_{ij}^{(k)} - E_i^{(k)} D_v f_{ij}^{(k)} + Q_{ij}^{(k)}, \text{ with } E_i^{(1)} = E_i^n. $$

The time integrator used here is a 2-stage second order explicit RK method with $a_{21} = 1$ and weights $b_1 = 0, b_2 = 1$ (explicit midpoint rule). The aim is to find an adapted approximation of $E_i^{(2)}$ that guarantees exact conservation of the total energy

$$ \mathcal{E}_{tot}(t) = \mathcal{E}_{kin}(t) + \mathcal{E}_{elec}(t), $$

where

$$ \mathcal{E}_{kin}(t) = \frac{1}{2} \int_{\mathbb{R} \times [0, L]} |v|^2 f dxdv \quad \text{and} \quad \mathcal{E}_{elec}(t) = \frac{1}{2} \int_{0}^{L} E^2(t, x) dx. $$

In other words, we need to satisfy the equality $\frac{d\mathcal{E}_{tot}(t)}{dt} = 0$ at the numerical level. If we multiply (22) by $\frac{v_j^2}{2}\Delta x \Delta v$ and we sum over indices $i \in 1, \dots, N_x, j \in \mathbb{Z}$, we get

$$ \frac{\mathcal{E}_{kin}^{n+1} - \mathcal{E}_{kin}^{n}}{\Delta t} = \left( -\frac{1}{2} \sum_{i,j} v_j^3 D_x f_{ij}^{(2)} - \frac{1}{2} \sum_{i,j} E_i^{(2)} v_j^2 D_v f_{ij}^{(2)} + \frac{1}{2} \sum_{i,j} v_j^2 Q_{ij}^{(k)} \right) \Delta x \Delta v, $$

the first sum vanishes thanks to periodic boundary conditions in space, and the last sum also vanishes by Proposition 1. Now, using discrete integration by parts we get

$$ 
\begin{aligned}
\frac{\mathcal{E}_{kin}^{n+1} - \mathcal{E}_{kin}^{n}}{\Delta t} &= \Delta x \Delta v \sum_{i,j} E_i^{(2)} v_j f_{ij}^{(2)} \\
&= -\Delta x \sum_{i} E_i^{(2)} \frac{E_i^{n+1} - E_i^n}{\Delta t},
\end{aligned} $$

where $E_i^{n+1} = E_i^n - \Delta t \Delta v \sum_j v_j f_{ij}^{(2)}$ is an approximation of the electric field at time $t^{n+1}$ using the midpoint rule quadrature applied to Ampère equation. Thus, we consider the following approximation of the electric field

$$ E_i^{(2)} = \frac{E_i^{n+1} + E_i^n}{2}, \qquad (23) $$
which implies the conservation of the total energy

$$ \frac{\mathcal{E}_{kin}^{n+1} - \mathcal{E}_{kin}^{n}}{\Delta t} + \Delta x \sum_{i} \frac{(E^{n+1})^2 - (E^n)^2}{\Delta t} = 0. $$

The above scheme uses RK2 for the whole equation. The scheme that we would like to consider is the one described in the second strategy at the beginning of the current section (SL – RK2 – RKC2).

**Remark.** A high order upwind scheme is used to approximate the velocity derivative instead of a centered scheme as in [25] which improves the stability of the resulting scheme, but the strategy can be extended to the Vlasov-Maxwell case.

# 5 Numerical results

This section is dedicated to the presentation of some numerical tests in different configurations: homogeneous in space, one-dimension in space and velocity and one-dimension in space and two dimensions in velocity. Note that if we write RKC without precising the order, then we mean the second order method.

## 5.1 Homogeneous case

First, we consider the simple homogeneous problem in one dimension satisfied by $f(t, v)$

$$ \partial_t f = \nu \partial_v ((v-u)f + T \partial_v f), \qquad (24) $$

with the initial condition

$$ f(t=0, v) = M_{\alpha,0,T_c}(v) + M_{\frac{1-\alpha}{2},4,1}(v) + M_{\frac{1-\alpha}{2},-4,1}(v), \quad \text{with} \quad M_{\rho,u,T} = \frac{\rho}{\sqrt{2\pi T}} e^{-\frac{(v-u)^2}{2T}}. \qquad (25) $$

with $T_c = 0.2$ and for the current section, we choose the following initial data with $\alpha = 0.9$. We consider a velocity domain $[-v_{max}, v_{max}]$ with $v_j = -v_{max} + j\Delta v$, $\Delta v = 2v_{max}/N_v$. The collision operator is discretized by the scheme (4) with Dirichlet boundary condition at $v_{max}$.

### 5.1.1 CFL condition and second order in time

Let us discuss in practice the advantages of RKC methods applied to $\partial_t f = \nu Q(f)$, where $\nu = 0.1$, $v_{max} = 12$, $N_v = 256$, $s = 20$ stages. The preserved quantities $T = 1.8822$, $u = 0$ are calculated using the initial condition. For RKC1, the time step must be chosen such that

$$ \Delta t \le \frac{1.93s^2}{4\nu T} \Delta v^2 \approx 9.09, $$
<figure><img src="image_4.png" /><figcaption>Figure 3: Illustration of the CFL condition for RKC2: plots of $f(t_{max} = 300, v)$ with different time steps.</figcaption></figure>

and for RKC2, we get

$$ \Delta t \le \frac{0.65s^2}{4\nu T} \Delta v^2 \approx 3.06. $$

The CFL condition for standard explicit methods such as explicit Euler or RK2 reads

$$ \Delta t \le \frac{1}{2\nu T} \Delta v^2 \approx 0.023. $$

This means that for a final time $t_{max} = 100$, the cost of RKC1 is: $(t_{max}/\Delta t) \times s = (100/9.09) \times 20 \approx 221$ evaluations of $Q$, and for RKC2: $(100/3.06) \times 20 \approx 654$ evaluations of $Q$. On the other hand, the cost of explicit Euler method: $100/0.023 \approx 4348$ evaluations of $Q$ and that of RK2:$(100/0.023) \times 2 \approx 8696$ evaluations of $Q$. Note that here we are not taking into account the possible restriction coming from the advection in $v$, since for the above values the eigenvalues are all on the negative real axis (as in Figure 1a). To illustrate this for the RKC2 method, we plot in Figure 3 the solution for three different time steps at final time $t_{max} = 300$. We observe that the solution reaches the Maxwellian equilibrium for $\Delta t = 0.5, 3.05$ but for $\Delta t = 3.085$, the method becomes unstable, as predicted by the previous analysis.

Now, still with $N_v = 256$, $\nu = 0.1$, $v_{max} = 12$, and final time $t_{max} = 30$, we show in Figure 4 the second order in time of the RKC2 integrator. The values of $T$ and $u$ are calculated according to (2), and $s$ is calculated adaptively in order to guarantee stability according to strategies presented in 3.3. The used time steps are $\Delta t = 0.75/2^j$, $j = 1, \dots, 8$, and for the referenc solution we use $\Delta t = 3 \times 10^{-4}$.
<figure><img src="image_5.png" /><figcaption>Figure 4: $L^2$ error at the final time $t_{max} = 30$ versus time step of the RKC2 scheme.</figcaption></figure>

### 5.1.2 Fourth order approximation

We also consider the fourth order approximation (16) of the collision operator. First we illustrate the fourth order accuracy proven in Proposition 2. In Figure 5 (left), we plot the $L^2$ error

$$ \left( \sum_{j=0}^{N_v-1} |Q(f)(v_j) - Q_j|^2 \Delta v \right)^{1/2} $$

(for different $N_v = 2^k, k = 6, \dots, 13$ with $v_{max} = 14$) where $Q_j$ corresponds to the second or fourth order scheme, and $Q(f)(v_j)$ denotes the exact collision operator obtained for $f(v) = M_{0,1,4,1}(v) + M_{0,9,0,0,2}(v)$. The expected order of accuracy are recovered. Then, we consider the RKC method and we look at $f_j^\infty - M_{\rho,u,T}(v_j)$ where $(\rho, u, T) = (1, 0, 1.88)$ are the moments of the initial condition (which are preserved by the scheme) and $f_j^\infty$ is the numerical solution of (24) computed at a large time $t = 1000$ so that the entropy does not change and the equilibrium is reached. We consider the following numerical parameters $\Delta t = 2, s = 20, v_{max} = 14$ and two different values for $N_v$. The results in Figure 5 (right) show that the fourth order accuracy can be observed for the asymptotic numerical solution. Similar results are obtained using a standard RK2 time integrator, but using a much smaller time step.

### 5.1.3 Adaptive time stepping and conservation

In this section, we test the performance of the algorithm with time step control strategy presented in 3.3 by considering the homogeneous problem (24). The tolerance used throughout this section is $tol = 10^{-6}$.
<figure><img src="image_6.png" /><figcaption>Figure 5: Left: $L^2$ error on $Q$ as a function of $N_v$ ($N_v = 2^k, k = 6, \dots, 13$) for the second order scheme (4) and fourth order scheme (16). Right: $\max_j |f_j^\infty - M_{\rho,u,T}(v_j)|$ as a function of $N_v$ ($N_v = 100, 150, 200, 250, 300$) for $\Delta t = 2, t_f = 1000$ for the second order scheme (4) and fourth order scheme (16).</figcaption></figure>

<table><thead><tr><th></th><th>RK2</th><th>RKC2</th></tr></thead><tbody><tr><th>Q evals</th><td>8794</td><td>1013</td></tr><tr><th>CPU time</th><td>≈ 0.25s</td><td>≈ 0.05s</td></tr></tbody></table>

<table><thead><tr><th></th><th>RK2</th><th>RKC2</th></tr></thead><tbody><tr><th>Q evals</th><td>42752</td><td>1488</td></tr><tr><th>CPU time</th><td>≈ 1.6s</td><td>≈ 0.05s</td></tr></tbody></table>

<table><caption>Table 1: Cost comparison of the RKC2 and the RK2 methods for $t_{max} = 100$ (left), $t_{max} = 500$ (middle) and $t_{max} = 1000$ (right).</caption><thead><tr><th></th><th>RK2</th><th>RKC2</th></tr></thead><tbody><tr><th>Q evals</th><td>85200</td><td>1851</td></tr><tr><th>CPU time</th><td>≈ 5s</td><td>≈ 0.06s</td></tr></tbody></table>

## RKC2 versus RK2

We consider the same parameters as above ($N_v = 256, \nu = 0.1, v_{max} = 12$) and we compare the performance of the second order Runge's method (also called explicit midpoint rule or RK2) with Euler method for the error estimation and the RKC2 scheme with the error estimator explained in Section 3.3. Table 1 shows the advantage of RKC2 in terms of function evaluations and CPU time for different values of final time. Note that the number of rejected steps is negligible in all the cases for both methods.

## Conservation

Finally, we investigate the properties of the velocity discretization coupled with the time integrator. To do so, we consider RKC2 with the adaptive time stepping presented in 3.3 with the following parameters $N_v = 256, \nu = 0.1, v_{max} = 16, tol = 10^{-6}$ and $\Delta t_0 = 10^{-3}$ (note that $v_{max}$
<figure><img src="image_7.png" /><figcaption>Figure 6: History of the time steps and number of stages.</figcaption></figure>

is larger to neglect the effect of the truncation; otherwise, specific boundary conditions can be employed as in [4]). In Figure 7, the time history of the relative mass, momentum and energy is displayed, as well as the time history of the entropy $E(t) = \int f^2(t)/Mdv$ and of $\|f(t, \cdot) - M(\cdot)\|_{L^2}$ (with $M$ the constant Maxwellian associated to the initial condition 25). First, we observe the conservation of the three invariants up to machine accuracy. Second, the entropy enjoys a very rapid decay relaxation of the solution towards an equilibrium which corresponds to the Maxwellian $M$ which shares the three first moments of 25 as observed in the right figure 7. Even if the entropy decay can be proved in the semi-discrete case, it is not an easy task to prove it in the fully discrete case.

<figure><img src="image_8.png" /><figcaption>Figure 7: Left: Time history of the relative error of mass, momentum, and energy. Middle: Time evolution of $E(t) = \int \frac{f^2}{M} dv$. Right: Time evolution of $\|f(t, \cdot) - M\|_{L^2}$.</figcaption></figure>
## 5.2 Non-homogeneous case 1dx-1dv

Here, we solve the VPFP equation (1) for $d_x = d_v = 1$ with the two different methods proposed in Section 4, namely SL-RKC and SL-RK2-RKC. Different tests are proposed to illustrate the capability of the approach in the non-homogeneous case.

### Test 1

We start with the following initial data

$$f(0, x, v) = \frac{1}{\sqrt{2\pi}} e^{-\frac{v^2}{2}} (1 + 10^{-3} \cos(kx)), \quad x \in [0, 2\pi/k],$$

with $k = 0.5$. The results of this so-called (collisionless, $\nu = 0$) Landau damping test are displayed on Figures 8: the time evolution of the electric energy (left) and the total mass, momentum and total energy. First, the expected damping rate is observed on the electric energy and the three invariants are preserved up to machine accuracy. The method used here is SL-RK2-RKC with $\Delta t = 0.1$ is chosen quite small to ensure the stability of RK2 since in the collisionless regime the RKC is disabled and so its error estimator.

<figure><img src="image_9.png" /><figcaption>Figure 8: Test 1 (1dx-1dv): $N_x = 128$, $N_v = 256$, $v_{max} = 12$, $\nu = 0$.</figcaption></figure>

### Test 2

Then, we consider

$$f_h(v) = \frac{1}{\rho_h \sqrt{2\pi}} v^{\gamma} e^{-\frac{v^2}{2}} \quad \text{where} \quad \rho_h = \int_{v_{\min}}^{v_{\max}} f_h(v) dv, \quad f_c = \frac{1}{\sqrt{2\pi T_c}} e^{-\frac{v^2}{2T_c}},$$
so that the initial data is

$$f(0, x, v) = \alpha f_c(v) + (1 - \alpha) f_h(v)(1 + \beta \cos(kx)),$$

where $\gamma = 10, \beta = k = 0.5, T_c = 0.2$, and $\alpha = 0.9$. We consider here $\nu = 0.1$.

We compare the two approaches SL-RKC and SL-RK2-RKC. Both methods preserves mass and momentum but SL-RKC does not preserve the total energy (deviation of magnitude $10^{-4}$) whereas SL-RK2-RKC does. The time history of the mass, momentum and total energy are plotted in Figures 9.

<figure><img src="image_10.png" /><figcaption>Figure 9: Test 2 (1dx-1dv): Deviation of total energy for SL-RKC (left) and SL-RK2-RKC (right).</figcaption></figure>

## Test 3

We consider a test from [19] in which the initial condition is

$$f_0(x, v) = (M_{n_1,0,1} + M_{n_2/2,u,0.2} + M_{n_2/2,-u,0.2})(1 + \varepsilon \cos(kx)),$$

with $n_1 = 0.97, n_2 = 0.03, u = 4, x \in [0, L], L = 22, k = 2\pi/L, \varepsilon = 0.00056$. The numerical parameters are $N_x = 128, N_v = 256$ for the phase space discretization ($v_{max} = 14$) whereas $\Delta t = 1$ thanks to the use of SL-RK2-RKC.

In Figures 10, we display the time evolution of the electric energy (semi-log scale) for a collisionless case (left) and with $\nu = 5 \times 10^{-4}$ (right). First, an exponential growth is observed with rate $\gamma = 0.0746$ which is very close to the theoretical value proposed in [19]. Then, nonlinear effects saturate the exponential growth and the oscillations driving the particles trapping process
are observed. When $\nu > 0$, the evolution changes significantly: even if the exponential growth remains close to the collisionless case, the nonlinear phase is strongly affected by collisions since the saturation amplitude decreases as collisions smooth out the trapping effects and drives the system towards a Maxwellian shape. The influence of the collisions on the trapping process can be seen in Figure 11 where we plot the phase space distribution function at time $t_1 = 80$ and $t_2 = 220$ in the collisionless and collision case. The time $t_1$ corresponds to the end of the linear phase and $t_2$ corresponds to the nonlinear phase. In the collisionless case (figure 11, top row), the vortex around the phase velocity $v_\phi \approx 3.5$ can be observed whereas for $\nu > 0$ (figure 11, bottom row), the trapping structure is smoothed out. Finally in Figure 12, the time evolution of total mass, momentum and energy is plotted and we can see that in both collisionless and collision cases, these three invariants are preserved up to machine accuracy.

<figure><img src="image_11.png" /><figcaption>Figure 10: Test 3 (1dx-1dv): $N_x = 128$, $N_v = 256$, $v_{max} = 14$, $\Delta t = 1$. Left: $\nu = 0$, right: $\nu = 5 \times 10^{-4}$.</figcaption></figure>

### 5.3 Non-homogeneous case 1dx-2dv

Here, we solve the VPFP equation satisfied by $f(t, x, v)$ with $v = (v_x, v_y)$

$$ \partial_t f + v_x \partial_x f + E_x \partial_{v_x} f = \nabla_v \cdot ((v - u_f)f + T_f \nabla_v f), $$

with $n_f u_f = \int v f d v_x d v_y$ and $2n_f T_f = \int |v - u_f|^2 f d v_x d v_y$. In the following, we denote the Maxwellian by $M_{\rho, u_x, u_y, T} = \frac{\rho}{2\pi T} \exp(-[(v_x - u_x)^2 + (v_y - u_y)^2]/(2T))$. The numerical scheme (4) extends naturally to the multi-dimensional case and details are given in Appendix B. Let us remark that a three steps splitting is used.
<figure><img src="image_12.png" /><figcaption>Figure 11: Test 3 (1dx-1dv): $N_x = 128$, $N_v = 256$, $v_{max} = 14$, $\Delta t = 1$. First row: $\nu = 0$, $t = 80, 220$. Second row: $\nu = 5 \times 10^{-4}$, $t = 80, 220$.</figcaption></figure>

**Test 1** The initial condition is given by

$$f(t = 0, x, v_x, v_y) = M_{1,0,0,1}(1 + 10^{-4} \cos(kx)),$$

with $\alpha = 0.5$ whereas the spatial domain is $x \in [0, 2\pi/k]$ ($k = 0.3$) and the velocity domain is $v = (v_x, v_y) \in [-7, 7]^2$. We investigate here the weakly collisional regime ($\nu$ small) on the Landau damping phenomena (see [25, 27]).

For the weakly collisional Landau damping, we consider $N_x = 32$, $N_{v_x} = N_{v_y} = 64$. The linear theory predicts a damping in the collisionless case, ie the root of the dispersion relation is
<figure><img src="image_13.png" /><figcaption>Figure 12: Test 3 (1dx-1dv): $N_x = 128$, $N_v = 256$, $v_{max} = 14$, $\Delta t = 1$. Time evolution of the total energy. Left: $\nu = 0$, right: $\nu = 5 \times 10^{-4}$.</figcaption></figure>

$\omega = 1.1598 - 0.0126i$. In Figure 13, we plot the time history of the electric energy (in semi-log sale) for RKC2 and $\Delta t = 0.3, s = 6$. Let us remark that with this time step, RK2 is unstable since it requires $\Delta t \approx 10^{-2}$. For different values of the collision frequency $\nu = 0, 0.025, 0.05, 0.1$ we can observe that the electric energy is damped exponentially fast in time according to the collisionless Landau theory (the numerical damping $-0.012$ is in very good agreement with the linear theory $-0.0126$) in addition to the effect of the collisions which induces an additional damping rate (see [25, 27]).

**Test 2** Lastly, we consider the following initial condition

$$f(t = 0, x, v_x, v_y) = \left( \frac{(1 - \alpha)}{4} \sum_{i=1}^{4} M_{1, \pm 3, \pm 3, 1/2} + \alpha M_{1, 0, 0, 1} \right) (1 + 0.01 \cos(kx)),$$

with $\alpha = 0.5$ whereas the spatial domain is $x \in [0, 2\pi/k]$ ($k = 0.5$) and the velocity domain is $v = (v_x, v_y) \in [-18, 18]^2$. The velocity dependency of the initial condition is plotted in Figure 14 (left) for $x = \pi/k$. We choose the collisions frequency $\nu = 0.1$.

We compare different methods (applied on the collisional part) to emphasize the advantages of the stabilized Runge-Kutta methods. In Figure 15, we plot the time history of the electric energy in semi-log scale for different methods. A strange splitting between the transport in space, in velocity, and the collision operator is used, where transports in space and in velocity are treated using semi-Lagrangian methods.

First, in Figure 15a, we take $N_{v_x} = N_{v_y} = 96$. For the implicit method (implicit Euler only for the diffusion part of the collision operator) and the RKC methods we use $\Delta t = 0.5$, while RK2
<figure><img src="image_14.png" /><figcaption>Figure 13: Test 1 (1dx-2dv): $N_x = 32$, $N_{v_x} = N_{v_y} = 64$, $v_{max} = 7$, $\Delta t = 0.3$, $s = 3$. Influence of $v = 0, 0.025, 0.05, 0.1$ on the time evolution of the electric energy.</figcaption></figure>

is run with $\Delta t = 0.05$ for stability reasons. For RKC methods, we choose the smallest number of stages such that the scheme is stable ($s = 3$ for RKC1 and $s = 5$ for RKC2). For the explicit methods the decay stops before machine precision which means that the solution converges to a numerical equilibrium which is $O(\Delta v^2)$ far from the exact equilibrium as expected. Even if the implicit step might seem to be more optimized, the cost involved by this step turns out to be much more important compared to the explicit stabilized methods (factor 100 in our (non optimized) implicit solver). On the other side, using a standard Runge-Kutta method (RK2) needs the use of small time steps to ensure stability. Finally, RKC1 or RKC2 can run with time steps comparable to the ones used by implicit methods and avoid solving large linear systems.

In Figure 15b, we consider $N_{v_x} = N_{v_y} = 128$ and we use $\Delta t = 0.5$ for RKC methods ($s = 4$
for RKC1 and $s = 7$ for RKC2), while RK2 requires $\Delta t = 0.025$. We can see that, for explicit methods, the convergence to the Maxwellian state can be improved (to almost machine precision) by refining the velocity mesh. This shows again the high importance of stabilized methods since finer meshes increase the stiffness of the system which in turn increases the restriction on the time step for standard explicit methods, while stabilized methods can still run with the same large time steps but with just a few more internal stages.

In Figure 14 (right), the velocity dependence of $f$ at time $t = 50$ and $x = \pi/k$ is plotted which shows the relaxation of the non equilibrium initial condition to the Maxwellian state due to the effect of the collisions.

<figure><img src="image_15.png" /><figcaption>Figure 14: $N_x = 32$, $N_{v_x} = N_{v_y} = 96$, $v_{max} = 18$, $\nu = 0.1$. Initial data $f(t = 0, v_x, v_y, x = \pi/k)$ and final solution $f(t = 50, v_x, v_y, x = \pi/k)$.</figcaption></figure>

# 6 Conclusion

In this work, we used stabilized Runge-Kutta methods in the context of the numerical simulation of Vlasov-Fokker-Planck equations. These time integrators are particularly well adapted to the weakly collisional kinetic equations since they are explicit and allow the use of large time steps in contrast to classical implicit methods (which require to solve large linear systems) or explicit methods (which need to use small time steps for stability reasons). Regarding the Fokker-Planck operator, a new discretization is proposed based on [4] and a new time integrator is presented for the time integration of the Vlasov-Poisson-Fokker-Planck equation which enables to preserve the total energy, even when coupled with adaptive time step selection strategy. Several tests illustrate the very good behavior of the approach (both in terms of accuracy and in terms of computational cost) compared to standard approaches.
<figure><img src="image_16.png" /><figcaption>Figure 15: Test 2 (1dx-2dv): $N_x = 32$, $v_{max} = 18$. Comparison of the time evolution of the electric energy in semi-log scale for different time integrators: RK2 ($\Delta t$ = 0.05 (left) and $\Delta t$ = 0.025 (right)), implicit ($\Delta t$ = 0.5), RKC2 ($\Delta t$ = 0.5), RKC1 ($\Delta t$ = 0.5).</figcaption></figure>

Several extensions like the coupling with the Maxwell equations and taking into account the nonlinear Landau operator can be considered.

A Proof of Proposition 1

This part is dedicated to the proof of Proposition 1.

First, we prove the conservation property. Thanks to the fact that (4) can be written in a flux form (9), the mass conservation follows easily. For the conservation of momentum, using the definition (7) of $\tilde{u}$, we have

$$ 
\begin{aligned}
\Delta v \sum_j Q_j v_j &= \sum_j v_j \left[ f_{j+1/2}(v_{j+1/2} - \tilde{u}) - f_{j-1/2}(v_{j-1/2} - \tilde{u}) + \frac{\tilde{T}}{\Delta v} (f_{j+1} - 2f_j + f_{j-1}) \right] \\
&= \sum_j \left( v_{j+1/2} - \frac{\Delta v}{2} \right) f_{j+1/2} (v_{j+1/2} - \tilde{u}) - \sum_j \left( v_{j-1/2} + \frac{\Delta v}{2} \right) f_{j-1/2} (v_{j-1/2} - \tilde{u}) \\
&= \sum_j v_{j+1/2} f_{j+1/2} (v_{j+1/2} - \tilde{u}) - \sum_j v_{j-1/2} f_{j-1/2} (v_{j-1/2} - \tilde{u}) = 0,
\end{aligned}
$$

where we used the definition (7) of $\tilde{u}$ and discrete integration by parts.
For the conservation of energy, we insert $v_j^2 = (v_{j\pm 1/2} \mp \Delta v/2)^2$ to get

$$ \Delta v \sum_j Q_j \frac{v_j^2}{2} = \frac{1}{2} \sum_j \left( v_{j+1/2}^2 + \frac{\Delta v^2}{4} - v_{j+1/2} \Delta v \right) f_{j+1/2} \left( v_{j+1/2} - \tilde{u} \right) \\ - \frac{1}{2} \sum_j \left( v_{j-1/2}^2 + \frac{\Delta v^2}{4} + v_{j-1/2} \Delta v \right) f_{j-1/2} \left( v_{j-1/2} - \tilde{u} \right) + \tilde{T} \left( \sum_j f_j \Delta v \right). $$

Using the definition of $\tilde{u}$, we obtain

$$ \begin{aligned} \Delta v \sum_j Q_j \frac{v_j^2}{2} &= -\frac{\Delta v}{2} \sum_j (v_{j+1/2} - \tilde{u} + \tilde{u}) f_{j+1/2} (v_{j+1/2} - \tilde{u}) \\ &\quad -\frac{\Delta v}{2} \sum_j (v_{j-1/2} - \tilde{u} + \tilde{u}) f_{j-1/2} (v_{j-1/2} - \tilde{u}) + \tilde{T} n \\ &= -\frac{\Delta v}{2} \sum_j f_{j+1/2} (v_{j+1/2} - \tilde{u})^2 - \frac{\Delta v}{2} \sum_j f_{j-1/2} (v_{j-1/2} - \tilde{u})^2 + \tilde{T} n \\ &= -\Delta v \sum_j f_{j+1/2} (v_{j+1/2} - \tilde{u})^2 + n \tilde{T}, \end{aligned} $$

which is equal to zero thanks to the preservation of momentum and using the definition (5) of $\tilde{T}$.

Regarding the entropy inequality, thanks to the equivalence between (4) and (6), we consider (6) to get

$$ \begin{aligned} \sum_{j \in \mathbb{Z}} Q_j \left(\frac{f}{M}\right)_j \Delta v &= \sum_{j \in \mathbb{Z}} \frac{\tilde{T}}{\Delta v^2} \left\{ \tilde{M}_{j+1/2} \left[ \left(\frac{f}{M}\right)_{j+1} - \left(\frac{f}{M}\right)_j \right] \right\} \left(\frac{f}{M}\right)_j \Delta v \\ &\quad - \sum_{j \in \mathbb{Z}} \frac{\tilde{T}}{\Delta v^2} \left\{ \tilde{M}_{j-1/2} \left[ \left(\frac{f}{M}\right)_j - \left(\frac{f}{M}\right)_{j-1} \right] \right\} \left(\frac{f}{M}\right)_j \Delta v \\ &= \sum_{j \in \mathbb{Z}} \frac{\tilde{T}}{\Delta v^2} \left\{ \tilde{M}_{j+1/2} \left[ \left(\frac{f}{M}\right)_{j+1} - \left(\frac{f}{M}\right)_j \right] \right\} \left(\frac{f}{M}\right)_j \Delta v \\ &\quad - \sum_{j \in \mathbb{Z}} \frac{\tilde{T}}{\Delta v^2} \left\{ \tilde{M}_{j+1/2} \left[ \left(\frac{f}{M}\right)_{j+1} - \left(\frac{f}{M}\right)_j \right] \right\} \left(\frac{f}{M}\right)_{j+1} \Delta v \\ &= - \sum_{j \in \mathbb{Z}} \frac{\tilde{T}}{\Delta v^2} \left\{ \tilde{M}_{j+1/2} \left[ \left(\frac{f}{M}\right)_{j+1} - \left(\frac{f}{M}\right)_j \right]^2 \right\} \Delta v \le 0. \end{aligned} $$

To guarantee the last inequality, we need to ensure that $\tilde{M}_{j+1/2} > 0$. To do so, we prove in the sequel $\tilde{M}_{j+1/2} = M(v_{j+1/2}) + O(\Delta v^2)$ (with $M$ given by (8)) so that for $\Delta v$ small enough, we can conclude $\tilde{M}_{j+1/2} > 0$ since $M(v_{j+1/2}) > 0$.
To get the consistency of $\tilde{M}_{j+1/2}$, we consider (11) with $f_j$ replaced by $f(v_j)$ and perform Taylor expansions. By an abuse of notations, we will still denote $f(v_j)$ by $f_j$ (same for $M$) which allows us to perform Taylor expansions using $f'_j$ for the velocity derivative of $f$ at $v_j$.

Then, we perform Taylor expansions to reformulate $f_{j+1}M_j - f_jM_{j+1}$ as

$$f_{j+1}M_j - f_jM_{j+1} = \Delta v(f'_jM_{j+1/2} - f_jM'_{j+1/2} + \frac{\Delta v}{2}(f''_jM_{j+1/2} - f'_jM'_{j+1/2})) + \mathcal{O}(\Delta v^3).$$

Using $M_jM_{j+1} + M_{j+1/2}^2 + \mathcal{O}(\Delta v^2)$, the second term in (11) thus becomes

$$\frac{M_j M_{j+1} (f_{j+1} - f_j)}{f_{j+1} M_j - f_j M_{j+1}} = \frac{M_{j+1/2} (f'_j M_{j+1/2} + \frac{\Delta v}{2} f''_j M_{j+1/2}) + \mathcal{O}(\Delta v^3)}{(f'_j M_{j+1/2} - f_j M'_{j+1/2} + \frac{\Delta v}{2} (f''_j M_{j+1/2} - f'_j M'_{j+1/2})) + \mathcal{O}(\Delta v^3)}.$$

For the first term of (11), using $M'_{j+1/2} = -((v_{j+1/2} - \tilde{u})/\tilde{T})M_{j+1/2}$, we get

$$\frac{\Delta v M_j M_{j+1} f_{j+1/2} (v_{j+1/2} - \tilde{u})}{\tilde{T} f_{j+1} M_j - f_j M_{j+1}} = \frac{(M_{j+1/2}^2 + \mathcal{O}(\Delta v^2))(f_j + \frac{\Delta v}{2} f'_j + \mathcal{O}(\Delta v^2))(v_{j+1/2} - \tilde{u})}{\tilde{T} (f'_j M_{j+1/2} - f_j M'_{j+1/2} + \frac{\Delta v}{2} (f''_j M_{j+1/2} - f'_j M'_{j+1/2})) + \mathcal{O}(\Delta v^3)} \\ = \frac{M_{j+1/2} (-f_j M'_{j+1/2} - \frac{\Delta v}{2} f'_j M'_{j+1/2}) + \mathcal{O}(\Delta v^2)}{(f'_j M_{j+1/2} - f_j M'_{j+1/2} + \frac{\Delta v}{2} (f''_j M_{j+1/2} - f'_j M'_{j+1/2})) + \mathcal{O}(\Delta v^3)}.$$

Gathering now the two terms of (11) leads to

$$\tilde{M}_{j+1/2} = \frac{M_{j+1/2}(f'_j M_{j+1/2} + \frac{\Delta v}{2} f''_j M_{j+1/2} - f_j M'_{j+1/2} - \frac{\Delta v}{2} f'_j M'_{j+1/2}) + \mathcal{O}(\Delta v^2)}{(f'_j M_{j+1/2} - f_j M'_{j+1/2} + \frac{\Delta v}{2} (f''_j M_{j+1/2} - f'_j M'_{j+1/2})) + \mathcal{O}(\Delta v^3)} \\ = M_{j+1/2} + \mathcal{O}(\Delta v^2).$$

## B Extension of the Fokker-Planck discretization in the two-dimensional case

In this Appendix, we extend the numerical scheme presented in the one-dimensional case to the two-dimensional case. We denote $v = (v_x, v_y)$ the velocity variable so that the Fokker-Planck operator is

$$Q(f)(v) = \partial_{v_x}[(v_x - u_{f,x})f + T\partial_{v_x}f] + \partial_{v_y}[(v_y - u_{f,y})f + T_f\partial_{v_y}f],$$

where $n_f = \int f d\nu_x d\nu_y$, $n_f u_{f,x/y} = \int v_{x/y} f d\nu_x d\nu_y$ and $n_f T_f = \int |v - u|^2 f d\nu_x d\nu_y$ ($u = (u_x, u_y)$). The velocity space is discretized by $v_{x,i} = i \Delta v_x$, $v_{y,j} = j \Delta v_y$, $i, j \in \mathbb{Z}$ and the discretization of the
Fokker-Planck operator is

$$ 
\begin{aligned}
Q(f)(v_{x,i}, v_{y,j}) & \approx Q_{i,j} \\
& = \frac{1}{\Delta v_x} \left[ f_{i+1/2,j}(v_{x,i+1/2} - \tilde{u}_x) - f_{i-1/2,j}(v_{x,i-1/2} - \tilde{u}_x) \right] + \frac{\tilde{T}}{\Delta v_x^2} (f_{i+1,j} - 2f_{i,j} + f_{i-1,j}), \\
& + \frac{1}{\Delta v_y} \left[ f_{i,j+1/2}(v_{y,j+1/2} - \tilde{u}_y) - f_{i,j-1/2}(v_{y,j-1/2} - \tilde{u}_y) \right] + \frac{\tilde{T}}{\Delta v_y^2} (f_{i,j+1} - 2f_{i,j} + f_{i,j+1})
\end{aligned}
$$

with $f_{i+1/2,j} = (f_{i+1,j} + f_{i,j})/2$, $f_{i,j+1/2} = (f_{i,j+1} + f_{i,j})/2$ and

$$ 
\begin{aligned}
n &= \sum_{i,j} f_{i,j} \Delta v_x \Delta v_y, \\
n \tilde{u}_x &= \sum_{i,j} v_{x,i+1/2} f_{i+1/2,j} \Delta v_x \Delta v_y, \quad n \tilde{u}_y \sum_{i,j} v_{y,i+1/2} f_{i,j+1/2} \Delta v_x \Delta v_y, \\
2n \tilde{T} &= \sum_{i,j} (v_{x,i+1/2} - \tilde{u}_x)^2 f_{i+1/2,j} \Delta v_x \Delta v_y + \sum_{i,j} (v_{y,j+1/2} - \tilde{u}_y)^2 f_{i,j+1/2} \Delta v_x \Delta v_y,
\end{aligned}
\quad (26) $$

where $v_{x,i+1/2} = (v_{x,i} + v_{x,i+1})/2$, $v_{y,j+1/2} = (v_{y,j} + v_{y,j+1})/2$. Similar expressions are obtained for the fourth order discretization.

## References

[1] R. BEMEJO, J. CARPIO *An adaptive finite element semi-Lagrangian implicit-explicit Runge-Kutta-Chebyshev method for convection dominated reaction-diffusion problems*, Applied Numerical Mathematics 58 (2008), pp. 16-39.

[2] A. HAKIM, M. FRANCISQUEZ, J. JUNO, G. HAMMETT, *Conservative discontinuous Galerkin schemes for nonlinear Fokker-Planck collision operators* Journal of Plasma Physics, Volume 86, Issue 4, 2020.

[3] S. Y. CHO, S. BOSCARINO, G. RUSSO, S. B. YUN, *Conservative semi-Lagrangian schemes for kinetic equations. Part I: reconstruction*, J. Comput. Phys. 432, 110159, (2021).

[4] C. BUET, S. DELLACHERIE, R. SENTIS *Numerical solution of an ionic Fokker-Planck equation with electronic temperature*, SIAM Journal on Numerical Analysis, Vol. 39, No. 4 (2002), pp. 1219-1253.

[5] J. S. CHANG, G. COOPER, *A practical difference scheme for Fokker-Planck equations*, J. Comput. Phys., 6 (1970), pp. 1-16.
[6] P. CRANDALL, D. JAREMA, H. DOERK, Q. PAN, G. MERLO, T. GÖRLER, A. BANON NAVARRO, D. TOLD, M. MAURER, F. JENKO, *Multi-Species Collisions for Delta-F Gyrokinetic Simulations: Implementation and Validation with GENE*, Computer Physics Communications (2020).

[7] N. CROUSEILLES, F. FILBET, *Numerical approximation of collisional plasmas by high order methods*, J. Comput. Phys., 201, pp. 546-572, (2004).

[8] F. FILBET, L. PARESCHI, *A numerical method for the accurate solution of the Fokker-Planck-Landau equation in the nonhomogeneous case*, J. Comput. Phys., 340, pp. 470-497, (2017).

[9] C. ZHANG, I. GAMBA, *A conservative scheme for Vlasov-Poisson-Landau modelling collisional plasmas*, J. Comput. Phys., 179, pp. 1-26, (2002).

[10] PETER DONNEL, XAVIER GARBET, YANICK SARAZIN, VIRGINIE GRANDGIRARD, YUUICHI ASAHI, NICOLAS BOUZAT, ELISABETTA CASCHERA, GUILHEM DIF-PRADALIER, CHARLES EHRLACHER, PHILIPPE GHENDRIH, GUILLAUME LATU, CHANTAL PASSERON, *A multi-species collisional operator for full-F global gyrokinetics codes: Numerical aspects and validation with the GYSELA code*, Computer Phys. Comm., 234, pp. 1-13, (2019).

[11] J. DRONIOU, C. LE POTIER, *Construction and convergence study of schemes preserving the elliptic local maximum principle*, SIAM J. Numer. Anal. 49 (2), 459-490, 2011

[12] F. FILBET, C. NEGULESCU, C. YANG, *Numerical study of a nonlinear heat equation for plasma physics*, International Journal of Computer Mathematics vol. 89 issue 8 (2012), pp. 1060-1082.

[13] F. FILBET, C. NEGULESCU, *Fokker-Planck multi-species equations in the adiabatic asymptotics*, to appear J. Comput. Phys. (2023).

[14] J. HU, X. ZHANG, *Positivity-preserving and energy-dissipative finite difference schemes for the Fokker-Planck and Keller-Segel equations* https://arxiv.org/pdf/2103.16790.pdf

[15] J. HU, J.-G. LIU, Y. XIE, Z. ZHOU, *A structure preserving numerical scheme for Fokker-Planck equations of neuron networks: Numerical analysis and exploration*, J. Comput. Phys. 433 110195 (2021).

[16] H. LIU, H. YU, *An entropy satisfying conservative method for the Fokker-Planck equation equation of the finitely extensible nonlinear elastic dumbbell model*, SIAM Journal on Numerical Analysis, 2012, Vol. 50, No. 3 (2012), pp. 1207-1239

[17] Q. PAN, D. R. ERNST, P. CRANDALL, *First implementation of gyrokinetic exact linearized Landau collision operator and comparison with models*, Phys. Plasmas 27, 042307 (2020).
[18] S. TANAKA, K. YOSHIKAWA, T. MINOSHIMA, N. YOSHIDA, *Multidimensional Vlasov-Poisson simulations with high-order monotonicity and positivity-preserving schemes*, The Astrophysical Journal 849:76, 2017.

[19] O. PEZZI, F. VALENTINI, P. VELTRI, *Nonlinear regime of electrostatic waves propagation in presence of electron-electron collisions*, Phys. of Plasmas, 22 042112 (2015).

[20] C. ZBINDEN, *Partitioned Runge-Kutta Chebyshev methods for diffusion-advection-reaction problems*, SIAM J. Sci. Comput., Vol. 33, No. 4, pp. 1707-1725.

[21] E. HAIRER, S. NØORSET, G. WANNER *Solving ordinary differential equations I. Nonstiff problems*, Volume 8. Springer-Verlag, Series in Comput. Math., Berlin, 1993.

[22] E. HAIRER, G. WANNER *Solving ordinary differential equations II. Stiff and differential-algebraic problems*, Springer-Verlag, Berlin and Heidelberg, 1996.

[23] A. ABDULLE *Explicit Stabilized Runge-Kutta Methods*, Encyclopedia of Applied and Computational Mathematics, Springer Berlin Heidelberg, pages 460–468, 2015.

[24] B. SOMMEIJER, L. SHAMPINE, J. VERWER *RKC: an explicit solver for parabolic PDEs*, J. Comput. Appl. Math., 88:316–326, 1998.

[25] N. CROUSEILLES, F. FILBET *Numerical approximation of collisional plasmas by high order methods*, J. Comput. Phys., 201:546–572, 2004.

[26] R. B. HORNE, M. P. FREEMAN *A New Code for Electrostatic Simulation by Numerical Integration of the Vlasov and Ampère Equations Using MacCormack's Method*, J. Comput. Phys. 171:182–200, 2001.

[27] I TRISTANI, *Landau damping for the linearized Vlasov Poisson equation in a weakly collisional regime*, J. Stat. Phys. 169, 1 (2017), 107–125.

[28] I. ALMUSLIMANI *A fully adaptive explicit stabilized integrator for advection-diffusion-reaction problems*, To appear in BIT Numerical Mathematics, 2022.

[29] A. ABDULLE, I. ALMUSLIMANI, AND G. VILMART *Optimal explicit stabilized integrator of weak order 1 for stiff and ergodic stochastic differential equations*, SIAM/ASA J. Uncertain. Quantif., 6(2):937–964, 2018.

[30] A. ABDULLE AND G. VILMART *PIROCK: a swiss-knife partitioned implicit-explicit orthogonal Runge-Kutta Chebyshev integrator for stiff diffusion-advection-reaction problems with or without noise*, J. Comput. Phys., 242:869–888, 2013.

[31] I. ALMUSLIMANI AND G. VILMART *Explicit stabilized integrators for stiff optimal control problems*, SIAM J. Sci. Comput., 43(2):A721–A743, 2021.
[32] A. ABDULLE AND T. LI S-ROCK methods for stiff Ito SDEs, Commun. Math. Sci., 6(4):845-868, 2008.

[33] A. ABDULLE, G. VILMART, AND K. C. ZYGALAKIS Weak second order explicit stabilized methods for stiff stochastic differential equations, SIAM J. Sci. Comput., 35(4):A1792-A1814, 2013.