Conservative and non-conservative methods based on
Hermite Weighted Essentially Non Oscillatory
reconstruction for Vlasov equations.

Chang Yang, Francis Filbet

► To cite this version:

Chang Yang, Francis Filbet. Conservative and non-conservative methods based on Hermite Weighted Essentially Non Oscillatory reconstruction for Vlasov equations.. Journal of Computational Physics, 2014, 279, pp.18-36. hal-00918077

HAL Id: hal-00918077

https://hal.science/hal-00918077

Submitted on 13 Dec 2013

**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.

Public Domain
CONSERVATIVE AND NON-CONSERVATIVE METHODS BASED
ON HERMITE WEIGHTED ESSENTIALLY-NON-OSCILLATORY
RECONSTRUCTION FOR VLASOV EQUATIONS

CHANG YANG AND FRANCIS FILBET

ABSTRACT. We introduce a WENO reconstruction based on Hermite interpolation both for semi-Lagrangian and finite difference methods. This WENO reconstruction technique allows to control spurious oscillations. We develop third and fifth order methods and apply them to non-conservative semi-Lagrangian schemes and conservative finite difference methods. Our numerical results will be compared to the usual semi-Lagrangian method with cubic spline reconstruction and the classical fifth order WENO finite difference scheme. These reconstructions are observed to be less dissipative than the usual weighted essentially non-oscillatory procedure. We apply these methods to transport equations in the context of plasma physics and the numerical simulation of turbulence phenomena.

KEYWORDS. Finite difference method; semi-Lagrangian scheme; Hermite WENO reconstruction; Vlasov-Poisson model; Guiding-center model; Plasma physics.

CONTENTS

<table><tr><td>1. Introduction</td><td>1</td></tr><tr><td>2. The Vlasov equation and related models</td><td>3</td></tr><tr><td>3. Hermite WENO reconstruction for semi-Lagrangian methods</td><td>4</td></tr><tr><td>3.1. Third order Hermite WENO interpolation</td><td>5</td></tr><tr><td>3.2. Fifth order Hermite WENO interpolation</td><td>6</td></tr><tr><td>4. Hermite WENO reconstruction for conservative finite difference methods</td><td>7</td></tr><tr><td>5. Numerical simulation of Vlasov equation and related models</td><td>9</td></tr><tr><td>5.1. 1D transport equation</td><td>9</td></tr><tr><td>5.2. Simplified paraxial Vlasov-Poisson model</td><td>11</td></tr><tr><td>5.3. Guiding center model</td><td>14</td></tr><tr><td>6. Conclusion and perspective</td><td>17</td></tr><tr><td>Acknowledgment</td><td>17</td></tr><tr><td>References</td><td>18</td></tr></table>

1. INTRODUCTION

Turbulent magnetized plasmas are encountered in a wide variety of astrophysical situations like the solar corona, accretion disks, but also in magnetic fusion devices such as tokamaks. In practice, the study of such plasmas requires solving the Maxwell equations coupled to the computation of the plasma response. Different ways are possible to compute this response: the fluid or the kinetic description. Unfortunately the fluid approach seems to be insufficient when one wants to study the behavior of zonal flow, the interaction between waves and particles or the occurrence

<footnote>The authors are partially supported by the European Research Council ERC Starting Grant 2009, project 239983-NuSiKiMo.</footnote>
of turbulence in magnetized plasmas for example. Most of the time these plasmas
are weakly collisional, and then they require a kinetic description represented by
the Vlasov-Maxwell system. The numerical simulation of the full Vlasov equation
involves the discretization of the six-dimensional phase space (x, v) ∈ ℝ³×ℝ³, which
is still a challenging issue. In the context of strongly magnetized plasmas however,
the motion of the particles is particular since it is confined around the magnetic
field lines; the frequency of this cyclotron motion is faster than the frequencies of
interest. Therefore, the physical system can be reduced to four or five dimensions
by averaging over the gyroradius of charged particles (See for a review [1, 14]).

The development of accurate and stable numerical techniques for plasma turbu-
lence (4D drift kinetic, 5D gyrokinetic and 6D kinetic models) is one of our long
term objectives.

Actually there are already a large variety of numerical methods based on di-
rect numerical simulation techniques. The Vlasov equation is discretized in phase
space using either semi-Lagrangian [7, 8, 26, 27], finite element [17], finite difference
[6, 28, 10] or discontinuous Galerkin [4, 18] schemes. Most of these methods are
based on a time splitting discretization which is particularly efficient for classical
systems as Vlasov-Poisson or Vlasov-Maxwell systems. In that case, the characte-
ristic curves corresponding to the split operator are straight lines and are solved
exactly. Therefore, the numerical error is only due to the splitting in time and the
phase space discretization of the distribution function. Furthermore for such time
splitting schemes, the semi-Lagrangian methods on Cartesian grids coupled with
Lagrange, Hermite or cubic spline interpolation techniques are conservative [2, 8].
Hence, these methods are now currently used and have proved their efficiency for
various applications. In this context semi-Lagrangian methods are often observed to
be less dissipative than classical finite volume or finite difference schemes. However,
for more elaborated kinetic equations like the 4D drift kinetic [15] or 5D gyrokinetic
[16] equations, or even the two dimensional guiding center model [27], time split-
ting techniques cannot necessarily be applied. Thus characteristic curves are more
sophisticated and required a specific time discretization. For instance, in [15, 16]
several numerical solvers have been developed using an Eulerian formulation for
gyro-kinetic models. However, spurious oscillations often appear in the non-linear
phase when small structures occur and it is difficult to distinguish physical and
numerical oscillations. Moreover, for these models semi-Lagrangian methods are
no more conservative, hence the long time behavior of the numerical solution may
become unsuitable.

For this purpose, we want to develop a class of numerical methods based on
the Hermite interpolation which is known to be less dissipative than Lagrange
interpolation [8] together with a weighted essentially non-oscillatory (WENO) re-
construction applied to semi-Lagrangian and finite difference methods. Actually,
Hermite interpolation with WENO schemes were already studied in [23] in the con-
text of discontinuous Galerkin methods with slope limiters. A system of equations
for the unknown function and its first derivative is evolved in time and used in
the reconstruction. Moreover, a similar technique, called CIP (Cubic Interpolation
Propagation), has also been proposed for transport equations in plasma physics ap-
plications [22], but the computational cost is strongly increased since the unknown
and all the derivatives are advected in phase space. In [8], a semi-Lagrangian
method with Hermite interpolation has been proposed and shown to be efficient
and less dissipative than Lagrangian interpolation. In this latter case, the first
derivatives are approximated by a fourth order centered finite difference formula.
Here, we also apply a similar pseudo-Hermite reconstruction [8] and meanwhile introduce an appropriate WENO reconstruction to control spurious oscillation leading to nonlinear schemes. We develop third and fifth order methods and apply them to semi-Lagrangian (non-conservative schemes) and conservative finite difference methods. Our numerical results will be compared to the usual semi-Lagrangian method with cubic spline reconstruction [27] and the classical fifth order WENO finite difference scheme [20].

The paper is organized as follows: we first present the Vlasov equation and related models which will be investigated numerically. Then in Section 3, the semi-Lagrangian method is proposed with high order Hermite interpolation with a WENO reconstruction to control spurious oscillations. In Section 4, conservative finite difference schemes with Hermite WENO reconstructions are detailed. The Section 5 the one-dimensional free transport equation with oscillatory initial data is investigated to compare our schemes with classical ones (semi-Lagrangian with cubic spline interpolation and conservative finite difference schemes with WENO reconstruction). Then we perform numerical simulations on the simplified paraxial Vlasov-Poisson model and on the guiding center model for highly magnetized plasma in two dimension.

## 2. THE VLASOV EQUATION AND RELATED MODELS

The evolution of the density of particles $f(t, \mathbf{x}, \mathbf{v})$ in the phase space $(\mathbf{x}, \mathbf{v}) \in \mathbb{R}^d \times \mathbb{R}^d$, $d = 1, \dots, 3$, is given by the Vlasov equation,

$$ (1) \qquad \frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla_{\mathbf{x}} f + \mathbf{F}(t, \mathbf{x}, \mathbf{v}) \cdot \nabla_{\mathbf{v}} f = 0, $$

where the force field $\mathbf{F}(t, \mathbf{x}, \mathbf{v})$ is coupled with the distribution function $f$ giving a nonlinear system. We mention the well known Vlasov-Poisson (VP) and Vlasov-Maxwell (VM) models describing the evolution of particles under the effects of self-consistent electro-magnetic fields. We define the charge density $\rho(t, \mathbf{x})$ and current density $\mathbf{J}(t, \mathbf{x})$ by

$$ (2) \qquad \rho(t, \mathbf{x}) = q \int_{\mathbb{R}^d} f(t, \mathbf{x}, \mathbf{v}) d\mathbf{v}, \quad \mathbf{J}(t, \mathbf{x}) = q \int_{\mathbb{R}^d} \mathbf{v} f(t, \mathbf{x}, \mathbf{v}) d\mathbf{v}, $$

where $q$ is the single charge. The force field is given for the Vlasov-Poisson model by

$$ (3) \qquad \mathbf{F}(t, \mathbf{x}, \mathbf{v}) = \frac{q}{m} \mathbf{E}(t, \mathbf{x}), \quad \mathbf{E}(t, \mathbf{x}) = -\nabla_{\mathbf{x}}\phi(t, \mathbf{x}), \quad -\Delta_{\mathbf{x}}\phi = \frac{\rho}{\varepsilon_0}, $$

where $m$ represents the mass of one particle. For the Vlasov-Maxwell system, we have

$$ (4) \qquad \mathbf{F}(t, \mathbf{x}, \mathbf{v}) = \frac{q}{m} (\mathbf{E}(t, \mathbf{x}) + \mathbf{v} \wedge \mathbf{B}(t, \mathbf{x})), $$

and $\mathbf{E}, \mathbf{B}$ are solution of the Maxwell equations

$$ (5) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial \mathbf{E}}{\partial t} - c^2 \operatorname{curl} \mathbf{B} = - \frac{\mathbf{J}}{\varepsilon_0}, \\[1.5ex] \displaystyle \frac{\partial \mathbf{B}}{\partial t} + \operatorname{curl} \mathbf{E} = 0, \\[1.5ex] \displaystyle \operatorname{div} \mathbf{E} = \frac{\rho}{\varepsilon_0}, \quad \operatorname{div} \mathbf{B} = 0, \end{array} \right. $$

with the compatibility condition

$$ (6) \qquad \frac{\partial \rho}{\partial t} + \operatorname{div}_{\mathbf{x}} \mathbf{J} = 0, $$
which is verified by the Vlasov equation solution.

In the sequel we will also consider the so-called guiding center model [3], which has been derived to describe the evolution of the charge density in a highly magnetized plasma in the transverse plane of a tokamak. This model is described as follows

$$ (7) \qquad \left\{ \begin{array}{l} \displaystyle \frac{\partial \rho}{\partial t} + \mathbf{U} \cdot \nabla \rho = 0, \\ -\Delta \phi = \rho, \end{array} \right. $$

where the velocity $\mathbf{U} = (-\partial_y\phi, \partial_x\phi)$ is divergence free.

Transport equations (1) or (7) can be recast into an advective form

$$ (8) \qquad \frac{\partial f}{\partial t} + \mathbf{A} \cdot \nabla f = 0, $$

where $\mathbf{A}: \mathbb{R}^{2d} \times \mathbb{R}^+ \to \mathbb{R}^{2d}$. Hence, classical backward semi-Lagrangian method can be applied to solve (8). Furthermore, under the assumption $\nabla \cdot \mathbf{A} = 0$, equations (1) or (7) can also be rewritten in a conservative form

$$ (9) \qquad \frac{\partial f}{\partial t} + \mathrm{div}(\mathbf{A}f) = 0. $$

for which a finite difference method can be used.

### 3. HERMITE WENO RECONSTRUCTION FOR SEMI-LAGRANGIAN METHODS

We introduce a high order Hermite interpolation coupled with a weight essentially non-oscillatory (HWENO) reconstruction for semi-Lagrangian methods. Actually, the semi-Lagrangian method becomes a classical method for the numerical solution of the Vlasov equation because of its high accuracy and its small dissipation [21, 27]. Moreover, it does not constrain any restriction on the time step size. Indeed, the key issue of the semi-Lagrangian method compared to classical Eulerian schemes is that it uses the characteristic curves corresponding to the transport equation to update the unknown from one time step to the next one. Let us recall the main feature of the backward semi-Lagrangian method. For a given $s \in \mathbb{R}^+$, the differential system

$$ \left\{ \begin{array}{l} \displaystyle \frac{d\mathbf{X}}{dt} = \mathbf{A}(t, \mathbf{X}), \\ \mathbf{X}(s) = \mathbf{x}, \end{array} \right. $$

is associated to the transport equation (8). We denote its solution by $\mathbf{X}(t; s, \mathbf{x})$. The backward semi-Lagrangian method is decomposed into two steps for computing the function $f^{n+1}$ at time $t_{n+1}$ from the function $f^n$ at time $t_n$:

(1) For each mesh point $\mathbf{x}_i$ of phase space, compute $\mathbf{X}(t_n; t_{n+1}, \mathbf{x}_i)$, the value of the characteristic at time $t_n$ who is equal to $\mathbf{x}_i$ at time $t_{n+1}$.

(2) As the function $f$ of transport equation verifies

$$ f^{n+1}(\mathbf{x}_i) = f^n(\mathbf{X}(t_n; t_{n+1}, \mathbf{x}_i)), $$

we obtain the value of $f^{n+1}(\mathbf{x}_i)$ by computing $f^n(\mathbf{X}(t_n; t_{n+1}, \mathbf{x}_i))$ by interpolation, since $\mathbf{X}(t_n; t_{n+1}, \mathbf{x}_i)$ is not usually a mesh point.

In practice, a cubic spline interpolation is often used [5, 13]. It gives very good results, but it has the drawback of being non local which causes a higher communication overhead on parallel computers. Moreover spurious oscillations may occur around discontinuities. On the other hand, the cubic Hermite interpolation is local, and has been shown in [8] to be less dissipative than Lagrange interpolation polynomial. However, it has still spurious oscillations for discontinuous solution.
Here, we develop a third and fifth order Hermite interpolation coupled with a weighted essentially non-oscillatory procedure, such that it is accurate for smooth solutions and it removes spurious oscillations around discontinuities or high frequencies which cannot be solved on a fixed mesh.

**3.1. Third order Hermite WENO interpolation.** Consider a uniform mesh $(x_i)_i$ of the computational domain and assume that the values of the distribution function $(f_i)_i$ and its derivative $(f'_i)_i$ are known at the grid points. The standard cubic Hermite polynomial $H_3(x)$ on the interval $I_i = [x_i, x_{i+1}]$ can be expressed as follows:

$$ (1) \qquad H_3(x) = f_i + \frac{f_{i+1} - f_i}{\Delta x}(x - x_i) + \frac{(f_{i+1} - f_i) - \Delta x f'_i}{\Delta x^2}(x - x_i)(x - x_{i+1}) \\ \qquad \qquad + \frac{\Delta x(f'_i + f'_{i+1}) - 2(f_{i+1} - f_i)}{\Delta x^3}(x - x_i)^2(x - x_{i+1}), $$

The polynomial $H_3(x)$ verifies:

$$ \begin{cases} H_3(x_i) = f_i, & H'_3(x_i) = f'_i, \\ H_3(x_{i+1}) = f_{i+1}, & H'_3(x_{i+1}) = f'_{i+1}. \end{cases} $$

Moreover, we define two quadratic polynomials on $I_i$ by

$$ \left\{ \begin{aligned} h_l(x) &= f_i + \frac{f_{i+1} - f_i}{\Delta x}(x - x_i) + \frac{(f_{i+1} - f_i) - \Delta x f'_i}{\Delta x^2}(x - x_i)(x - x_{i+1}), \\ h_r(x) &= f_i + \frac{f_{i+1} - f_i}{\Delta x}(x - x_i) + \frac{\Delta x f'_{i+1} - (f_{i+1} - f_i)}{\Delta x^2}(x - x_i)(x - x_{i+1}). \end{aligned} \right. $$

The polynomial $h_l$ verifies

$$ h_l(x_i) = f_i, \quad h_l(x_{i+1}) = f_{i+1}, \quad h'_l(x_i) = f'_i, $$

while $h_r$ verifies

$$ h_r(x_i) = f_i, \quad h_r(x_{i+1}) = f_{i+1}, \quad h'_r(x_{i+1}) = f'_{i+1}. $$

The idea of WENO reconstruction is now to apply the cubic polynomial $H_3$ when the function $f$ is smooth, otherwise, we use the less oscillatory second order polynomial between $h_l$ or $h_r$. Thus, let us write $H_3$ as follows

$$ H_3(x) = w_l(x) p_l(x) + w_r(x) p_r(x), $$

where $w_l$ and $w_r$ are WENO weights depending on $x$. When the function $f$ is smooth, we expect that

$$ w_l(x) \approx c_l(x) = \frac{x_{i+1} - x}{\Delta x} \quad \text{and} \quad w_r(x) \approx c_r(x) = 1 - c_l(x), $$

so that we recover the cubic Hermite polynomial. Otherwise, we expect that

$$ (2) \qquad w_l(x) \approx 1, w_r(x) \approx 0 \quad \text{or} \quad w_l(x) \approx 0, w_r(x) \approx 1 $$

according to the region where $f$ is less smooth. To determine these WENO weights, we follow the strategy given in [20] and first define smoothness indicators by integration of the first and second derivatives of $h_l$ and $h_r$ on the interval $I_i$:

$$ \left\{ \begin{aligned} \beta_l &= \int_{x_i}^{x_{i+1}} \Delta x (h'_l)^2 + \Delta x^3 (h''_l)^2 dx = (f_i - f_{i+1})^2 + \frac{13}{3} ((f_{i+1} - f_i) - \Delta x f'_i)^2, \\ \beta_r &= \int_{x_i}^{x_{i+1}} \Delta x (h'_r)^2 + \Delta x^3 (h''_r)^2 dx = (f_i - f_{i+1})^2 + \frac{13}{3} (\Delta x f'_{i+1} - (f_{i+1} - f_i))^2. \end{aligned} \right. $$

Then we set $w_l$ and $w_r$ as

$$ w_l(x) = \frac{\alpha_l(x)}{\alpha_l(x) + \alpha_r(x)} \quad \text{and} \quad w_r(x) = 1 - w_l(x), $$
where

$$ \alpha_l(x) = \frac{c_l(x)}{(\varepsilon + \beta_l)^2} \quad \text{and} \quad \alpha_r(x) = \frac{c_r(x)}{(\varepsilon + \beta_r)^2}. $$

where $\varepsilon = 10^{-6}$ to avoid the denominator to be zero.

Observe that when the function $f$ is smooth, the difference between $\beta_l$ and $\beta_r$ becomes small and the weights $w_l(x) \approx c_l(x)$ and $w_r(x) \approx c_r(x)$. Otherwise, when the smoothness indicator $\beta_s$, $s = l, r$ blows-up, then the parameter $\alpha_s$ and the weight $w_s$ goes to zero, which yields (2).

Finally, let us mention that here the value of the first derivative at the grid point $x_i$ is approximated by a fourth-order centered finite difference formula

$$ (3) \qquad f'_i = \frac{1}{12\Delta x} (8(f_{i+1} - f_{i-1}) - (f_{i+2} - f_{i-2})). $$

**3.2. Fifth order Hermite WENO interpolation.** We can extend previous method to a fifth order Hermite WENO (HWENO5) interpolation. In the same way, we first construct a fifth degree polynomial $H_5(x)$ on the interval $I_i$

$$ H_5(x_j) = f_j, j = i-1, i, i+1, i+2, \quad H'_5(x_{i-1}) = f'_{i-1}, H'_5(x_{i+2}) = f'_{i+2} $$

and then three third degree polynomials $h_l(x)$, $h_c(x)$, $h_r(x)$ verifying

$$ \left\{ \begin{array}{lll} h_l(x_j) = f_j, & j = i-1, i, i+1, & h'_l(x_{i-1}) = f'_{i-1}, \\ h_c(x_j) = f_j, & j = i-1, i, i+1, i+2, & \\ h_r(x_j) = f_j, & j = i, i+1, i+2, & h'_r(x_{i+2}) = f'_{i+2}, \end{array} \right. $$

where the first derivative $f'_i$ is given by a sixth order centered approximation

$$ f'_{i} = \frac{1}{60}((f_{i+3} - f_{i-3}) - 9(f_{i+2} - f_{i-2}) + 45(f_{i+1} - f_{i-1})). $$

Then the polynomial $H_5$ can be written as a convex combination

$$ H_5(x) = w_l(x)h_l(x) + w_c(x)h_c(x) + w_r(x)h_r(x), $$

where $w_l(x)$, $w_c(x)$, $w_r(x)$ are WENO weights depending on $x$. Similarly smoothness indicators are computed by integration of the first, second and third order derivatives of $h_l(x)$, $h_c(x)$, $h_r(x)$ on the interval $I_i$:

$$ \beta_j = \int_{x_i}^{x_{i+1}} \Delta x (h'_j)^2 + \Delta x^3 (h''_j)^2 + \Delta x^5 (h'''_j)^2 dx, \quad j = l, c, r. $$

Finally, the WENO weights are determined according to the smoothness indicators

$$ \left\{ \begin{array}{lll} w_l(x) = \frac{\alpha_l(x)}{\alpha_l(x) + \alpha_c(x) + \alpha_r(x)}, & \alpha_l(x) = \frac{c_l(x)}{(\varepsilon + \beta_l)^2}, & c_l(x) = \frac{(x - x_{i+2})^2}{9\Delta x^2}, \\ w_c(x) = \frac{\alpha_c(x)}{\alpha_l(x) + \alpha_c(x) + \alpha_r(x)}, & \alpha_c(x) = \frac{c_c(x)}{(\varepsilon + \beta_c)^2}, & c_c(x) = 1 - c_l(x) - c_r(x), \\ w_r(x) = \frac{\alpha_r(x)}{\alpha_l(x) + \alpha_c(x) + \alpha_r(x)}, & \alpha_r(x) = \frac{c_r(x)}{(\varepsilon + \beta_r)^2}, & c_r(x) = \frac{(x - x_{i-1})^2}{9\Delta x^2}. \end{array} \right. $$

This polynomial reconstruction allows to get fifth order accuracy for smooth stencil and the various stencils are expected to damp oscillations when filamentation of the distribution function occurs. Finally, let us observe that this technique can be easily extended to high space dimension on Cartesian grids.
## 4. HERMITE WENO RECONSTRUCTION FOR CONSERVATIVE FINITE DIFFERENCE METHODS

When the velocity $\mathbf{A}$ is not constant (8), the semi-Lagrangian method is not conservative even when $\mathrm{div}\mathbf{A} = 0$, hence mass is no longer conserved and the long time behavior of the numerical solution can be wrong even for small time steps. Therefore, high order conservative methods may be more appropriate even if they are restricted by a CFL condition. An alternative is to use the finite difference formulation in the conservative form and to use the semi-Lagrangian method for the flux computation [24, 25].

In this section, we extend Hermite WENO reconstruction for computing numerical flux of finite difference method. Suppose that $\{f_i\}_{1\le i\le N}$ is approximation of $f(x_i)$. We look for the flux $\{\hat{f}_{i+1/2}\}_{0\le i\le N}$ such that it approximates the derivative $f'(x)$ to k-th order accuracy:

$$ \frac{\hat{f}_{i+1/2} - \hat{f}_{i-1/2}}{\Delta x} = f'(x) + \mathcal{O}(\Delta x^k). $$

Let us define a function $g$ such that

$$ (1) \qquad f(x) = \frac{1}{\Delta x} \int_{x-\Delta x/2}^{x+\Delta x/2} g(s)ds, $$

then clearly

$$ f'(x) = \frac{1}{\Delta x} [g(x + \Delta x/2) - g(x - \Delta x/2)]. $$

Hence we only need

$$ \hat{f}_{i+1/2} \approx p(x_i + \Delta x/2). $$

Let us denote by $G$ one primitive of $g$

$$ G(x) = \int_{-\infty}^{x} g(s)ds, $$

then (1) implies

$$ G(x_{i+1/2}) = \sum_{j=-\infty}^{i} \int_{x_{j-1/2}}^{x_{j+1/2}} g(s)ds = \Delta x \sum_{j=-\infty}^{i} f_j := G_{i+1/2}. $$

Thus, given the point values $\{f_i\}_i$, the primitive function $G(x)$ is exactly known at $x = x_{i+1/2}$. We thus can approximate $G(x)$ by an interpolation method. Therefore,

$$ (2) \qquad g(x_{i+1/2}) = \left. \frac{dG}{dx} \right|_{x=x_{i+1/2}}. $$

Now let us interpolate the primitive function $G(x)$. Here we give the Hermite WENO scheme and outline the procedure of reconstruction only for the fifth order accuracy case.

The aim is to construct an approximation of the flux $f_{i+1/2}^{-}$ by the Hermite polynomial of degree five together with a WENO reconstruction from point values $\{f_i\}$:

(1) We construct the Hermite polynomial $H_5$ such that

$$ H_5(x_{j+1/2}) = G_{j+1/2}, j = -2, -1, 0, 1, \quad H'_5(x_{j+1/2}) = G'_{j+1/2}, j = -1, 0, $$
(2) We construct cubic reconstruction polynomials $H_l(x)$, $H_c(x)$, $H_r(x)$ such that:

$$ egin{cases} H_l(x_{j+1/2}) = G_{j+1/2}, & j = -2, -1, 0, \\ H_c(x_{j+1/2}) = G_{j+1/2}, & j = -2, -1, 0, 1, \\ H_r(x_{j+1/2}) = G_{j+1/2}, & j = -1, 0, 1, \end{cases} \quad \begin{aligned} H'_l(x_{i-1/2}) &= G'_{i-1/2}, \\ H'_c(x_{i-1/2}) &= G'_{i-1/2}, \\ H'_r(x_{i+1/2}) &= G'_{i+1/2}, \end{aligned} $$

where $G'_{i+1/2}$ is the sixth order centered approximation of first derivative. Let us denote by $h_l(x)$, $h_c(x)$, $h_r(x)$, $h_5(x)$ the first derivatives of $H_l(x)$, $H_c(x)$, $H_r(x)$, $H_5(x)$ respectively. By evaluating $h_l(x)$, $h_c(x)$, $h_r(x)$, $h_5(x)$ at $x = x_{i+1/2}$, we obtain

$$ h_5(x_{i+1/2}) = \frac{-8f_{i-1} + 19f_i + 19f_{i+1} + 3H'_{i-1/2} - 6H'_{i+1/2}}{27} $$

and

$$ h_l(x_{i+1/2}) = -2f_{i-1} + 2f_i + G'_{i-1/2}, $$

$$ h_c(x_{i+1/2}) = \frac{-f_{i-1} + 5f_i + 2f_{i+1}}{6}, $$

$$ h_r(x_{i+1/2}) = \frac{f_i + 5f_{i+1} - 2G'_{i+1/2}}{4}. $$

(3) We evaluate the smoothness indicators $\beta_l$, $\beta_c$, $\beta_r$, which measure the smoothness of $h_l(x)$, $h_c(x)$, $h_r(x)$ on the cell $I_i$.

$$ \begin{aligned} \beta_l &= \int_{x_i}^{x_{i+1}} \Delta x(h'_l(x))^2 + \Delta x^3(h''_l(x))^2 dx \\ &= \frac{1}{16} (835 f_{i-1}^2 + 139 f_i^2 + 300 (H'_{i-1/2})^2 - 674 f_{i-1} f_i - 996 f_{i-1} H'_{i-1/2} + 396 f_i H'_{i-1/2}), \end{aligned} $$

$$ \begin{aligned} \beta_c &= \int_{x_i}^{x_{i+1}} \Delta x(h'_c(x))^2 + \Delta x^3(h''_c(x))^2 dx \\ &= \frac{1}{12} (13 f_{i-1}^2 + 64 f_i^2 + 25 f_{i+1}^2 - 52 f_{i-1} f_i + 26 f_{i-1} f_{i+1} - 76 f_i f_{i+1}), \end{aligned} $$

$$ \begin{aligned} \beta_r &= \int_{x_i}^{x_{i+1}} \Delta x(h'_r(x))^2 + \Delta x^3(h''_r(x))^2 dx \\ &= \frac{1}{16} (55 f_i^2 + 367 f_{i+1}^2 + 156 (H'_{i+1/2})^2 - 266 f_i f_{i+1} + 156 f_i H'_{i+1/2} - 468 f_{i+1} H'_{i+1/2}). \end{aligned} $$

(4) We compute the non-linear weights based on the smoothness indicators

$$ \left\{ \begin{aligned} w_l &= \frac{\alpha_l}{\alpha_l + \alpha_c + \alpha_r}, & \alpha_l &= \frac{c_l}{(\varepsilon + \beta_l)^2}, \\ w_c &= \frac{\alpha_c}{\alpha_l + \alpha_c + \alpha_r}, & \alpha_c &= \frac{c_c}{(\varepsilon + \beta_c)^2}, \\ w_r &= \frac{\alpha_r}{\alpha_l + \alpha_c + \alpha_r}, & \alpha_r &= \frac{c_r}{(\varepsilon + \beta_r)^2}, \end{aligned} \right. $$
where the coefficients $c_l = 1/9$, $c_c = 4/9$, $c_r = 4/9$ are chosen to get fifth order accuracy for smooth solutions and the parameter $\epsilon = 10^{-6}$ avoids the blow-up of $\alpha_k$, $k = \{l, c, r\}$.

(5) The flux $f_{i+1/2}^{-}$ is then computed as

$$f_{i+1/2}^{-} = w_l h_l(x_{i+1/2}) + w_c h_c(x_{i+1/2}) + w_r h_r(x_{i+1/2}).$$

The reconstruction to $f_{i+1/2}^{+}$ is mirror symmetric with respect to $x_{i+1/2}$ of the above procedure.

## 5. NUMERICAL SIMULATION OF VLASOV EQUATION AND RELATED MODELS

We start with a very basic test on the one dimensional transport equation with constant velocity to check the order of accuracy and to compare the error amplitude of the various numerical schemes. Then we perform numerical simulations on the simplified paraxial Vlasov-Poisson model and on the guiding center model for highly magnetized plasma in the transverse plane of a tokamak.

In this section we will compare our Hermite WENO reconstruction with the usual semi-Lagrangian method with cubic spline interpolation without splitting [27], and with the classical fifth order finite difference technique [20] coupled with a fourth order Runge-Kutta scheme for the time discretization.

### 5.1. 1D transport equation.
We compare our Hermite WENO reconstruction with various classical methods for solving the free transport equation

$$(1) \qquad \frac{\partial f}{\partial t} + \frac{\partial f}{\partial x} = 0, \quad x \in [-1, 1], \quad t \ge 0,$$

with periodic boundary conditions.

Let us first consider a smooth solution, where the initial condition is chosen as

$$f(0, x) = \sin(\pi x), \quad x \in [-1, 1].$$

We present in Table 1, the numerical error for different methods. On the one hand for semi-Lagrangian methods, the Hermite WENO interpolation is compared with the cubic spline interpolation. The semi-Lagrangian method is unconditionally stable, we thus choose a CFL number larger than one, e.g. CFL = 2.5. We observe that the cubic spline and Hermite WENO reconstructions have both third order accuracy, and the numerical error has almost the same amplitude. The semi-Lagrangian method with a fifth order Hermite WENO reconstruction has fifth order accuracy, thus it is much more accurate than the previous third order methods.

On the other hand we focus on the finite difference method and compare the Hermite WENO reconstruction with the classical fifth order WENO reconstruction [20]. We observe that these two methods have fifth order accuracy, but the Hermite WENO interpolation method is much more accurate than the usual WENO method. Furthermore, for the same order of accuracy the semi-Lagrangian method is much more precise than the finite difference scheme, which is expected for linear problems since the error only comes from the polynomial interpolation.

We next consider a step function as follows

$$(2) \qquad f(0, x) = \begin{cases} 1, & \text{for } -1 \le x \le 0, \\ 0, & \text{otherwise.} \end{cases}$$

Comparisons of the methods are now summarized in Table 2. We first notice that all the methods can achieve order of accuracy of $\frac{p}{p+1}$, where p is degree of polynomial. On the one hand, it is clear that semi-Lagrangian methods are more precise than finite difference ones. More precisely for $n_x \le 10^3$, the cubic spline interpolation is more accurate than Hermite of degree three and five coupled
<table><caption>TABLE 1. 1D transport equation: Error in $L_1$-norm and order of convergence r for smooth solutions for semi-Lagrangian and finite difference methods. The final time is $T_{end} = 8$.</caption><thead><tr><th rowspan="2">n<sub>x</sub></th><th colspan="2">200</th><th colspan="2">400</th><th colspan="2">800</th></tr><tr><th>||·||<sub>1</sub></th><th>r</th><th>||·||<sub>1</sub></th><th>r</th><th>||·||<sub>1</sub></th><th>r</th></tr></thead><tbody><tr><td>Semi-Lagrangian cubic spline</td><td>1.03e-6</td><td>3.00</td><td>1.29e-7</td><td>3.00</td><td>1.61e-8</td><td>3.00</td></tr><tr><td>Semi-Lagrangian cubic HWENO</td><td>1.04e-6</td><td>3.03</td><td>1.29e-7</td><td>3.01</td><td>1.62e-8</td><td>3.00</td></tr><tr><td>Semi-Lagrangian HWENO 5th</td><td>9.28e-10</td><td>5.51</td><td>2.28e-11</td><td>5.35</td><td>8.63e-13</td><td>4.72</td></tr><tr><td>Finite difference WENO 5th</td><td>1.15e-7</td><td>4.99</td><td>3.60e-9</td><td>5.01</td><td>1.16e-10</td><td>4.95</td></tr><tr><td>Finite difference HWENO 5th</td><td>6.06e-8</td><td>4.99</td><td>1.92e-9</td><td>4.98</td><td>6.55e-11</td><td>4.87</td></tr></tbody></table>

with the WENO reconstruction. It illustrates perfectly the robustness of the semi-Lagrangian method with cubic spline interpolation. Nevertheless it also generates a lot of oscillations (see Table 2 (b)) which produce negative values of the distribution function. Furthermore the Hermite WENO5 reconstruction is less dissipative than usual WENO5 and it is more accurate. Both of them control well spurious oscillations (see Table 2 (b)).

<table><caption>(a) Error between exact solution and approximated solution</caption><thead><tr><th rowspan="2">n<sub>x</sub></th><th colspan="2">200</th><th colspan="2">400</th><th colspan="2">800</th></tr><tr><th>||·||<sub>1</sub></th><th>r</th><th>||·||<sub>1</sub></th><th>r</th><th>||·||<sub>1</sub></th><th>r</th></tr></thead><tbody><tr><td>Semi-Lagrangian cubic spline</td><td>2.47e-2</td><td>0.80</td><td>1.43e-2</td><td>0.79</td><td>8.52e-2</td><td>0.75</td></tr><tr><td>Semi-Lagrangian-HWENO 3rd</td><td>3.27e-2</td><td>0.78</td><td>1.89e-2</td><td>0.79</td><td>1.09e-2</td><td>0.80</td></tr><tr><td>Semi-Lagrangian-HWENO 5th</td><td>2.94e-2</td><td>0.84</td><td>1.63e-2</td><td>0.85</td><td>8.99e-3</td><td>0.86</td></tr><tr><td>Finite difference-WENO 5th</td><td>4.50e-2</td><td>0.83</td><td>2.53e-2</td><td>0.83</td><td>1.43e-2</td><td>0.82</td></tr><tr><td>Finite difference-HWENO 5th</td><td>4.07e-2</td><td>0.83</td><td>2.29e-2</td><td>0.83</td><td>1.29e-2</td><td>0.83</td></tr></tbody></table>

<table><caption>(b) Error of total variation</caption><thead><tr><th>n<sub>x</sub></th><th>200</th><th>400</th><th>800</th></tr></thead><tbody><tr><td>Semi-Lagrangian cubic spline</td><td>5.75e-1</td><td>5.12e-1</td><td>5.19e-1</td></tr><tr><td>Semi-Lagrangian-HWENO 3rd</td><td>5.71e-4</td><td>5.18e-4</td><td>4.37e-4</td></tr><tr><td>Semi-Lagrangian-HWENO 5th</td><td>1.09e-3</td><td>1.42e-3</td><td>1.46e-3</td></tr><tr><td>Finite difference-WENO 5th</td><td>9.54e-5</td><td>8.43e-5</td><td>6.58e-5</td></tr><tr><td>Finite difference-HWENO 5th</td><td>2.30e-3</td><td>2.47e-3</td><td>1.88e-3</td></tr></tbody></table>

<caption>TABLE 2. 1D transport equation: Comparison of different methods for the linear equation (1) with initial data (2). (a) Error in $L_1$ norm and r is the order of accuracy (b) Error on the total variation. The final time is $T_{end} = 8$.</caption>
We finally consider an oscillatory solution where the initial condition is given by [20],

$$ (3) \quad f(0, x) = \begin{cases} \frac{1}{6}[G(x, z - \delta) + G(x, z - \delta) + 4G(x, z)], & \text{for } -0.8 \le x \le -0.6, \\ 1, & \text{for } -0.4 \le x \le -0.2, \\ 1 - |10(x - 0.1)|, & \text{for } 0 \le x \le 0.2, \\ \frac{1}{6}[F(x, z - \delta) + F(x, z - \delta) + 4F(x, z)], & \text{for } 0.4 \le x \le 0.6, \\ 0, & \text{otherwise.} \end{cases} $$

where $G(x, z) = \exp(-\beta(x - z)^2)$, $F(x, a) = \{\max((1 - \alpha^2(x - a)^2)^{1/2}, 0)\}$ with $\alpha = 0.5$, $z = -0.7$, $\delta = 0.005$, $\alpha = 10$ and $\beta = (\log 2)/36\delta^2$.

We have similar observations of both regular solution and discontinuous solution cases in Figure 1. The usual semi-Lagrangian method with cubic spline interpolation approximates well exponential function, but involves too much oscillation in step function. The other methods with WENO or HWENO reconstruction avoid a lot spurious oscillations. Semi-Lagrangian methods is less dissipative than finite difference method, which can be seen in step function and peak function. Moreover, the finite difference method with fifth order Hermite WENO reconstruction is less dissipative than the one with usual fifth order WENO reconstruction.

### 5.2. Simplified paraxial Vlasov-Poisson model.
We apply the numerical methods presented in previous sections to the following Vlasov-Poisson system satisfied by $f(t, r, v)$, where $r \in \mathbb{R}$, $v \in \mathbb{R}$ [9, 3]

$$ (4) \quad \begin{cases} \frac{\partial f}{\partial t} + \frac{v}{\varepsilon} \frac{\partial f}{\partial r} + \left(E_f - \frac{r}{\varepsilon}\right) \frac{\partial f}{\partial v} = 0, \\ \frac{1}{r} \frac{\partial}{\partial r} (r E_f) = \int_{\mathbb{R}} f \, dv. \end{cases} $$

The electric field can be expressed explicitly as follows

$$ E_f(t, r) = \frac{1}{r} \int_0^r s \rho(t, s) \, ds, $$

where $\rho(t, r) = \int_{\mathbb{R}} f(t, r, v) \, dv$, hence we will compute $E_f$ by a simple numerical integration.

The initial condition is chosen as a Gaussian in velocity multiplied by a regularized step function in $r$:

$$ (5) \quad f_0(r, v) = \frac{4}{\sqrt{2\pi\alpha}} \chi(r) \exp(-\frac{v^2}{2\alpha}), $$

with $\chi(r) = \frac{1}{2}\operatorname{erf}(\frac{r+1.2}{0.3}) - \frac{1}{2}\operatorname{erf}(\frac{r-1.2}{0.3})$ and $\alpha = 0.2$. The Vlasov-Poisson system (4) conserves mass

$$ \frac{d}{dt} \int_{\mathbb{R}^2} f(t, r, v) dr dv = 0, $$

and also $L^p$ norm for $1 \le p < \infty$

$$ \frac{d}{dt} \|f(t, r, v)\|_{L^p(\mathbb{R}^2)} = 0. $$

Therefore, the evolution in time of these quantities will be observed for various numerical schemes. We will also investigate the time evolution of the kinetic energy of the Vlasov-Poisson system (4):

$$ (6) \quad \mathcal{E}(t) = \int_{\mathbb{R}^2} \frac{v^2}{2} f \, dr \, dv. $$
<figure><img src="image_1.png" /><figcaption>FIGURE 1. 1D transport equation: Plot solutions of the linear equation (1) with initial data (3). $n_x = 200$, $CFL= 2.5$ for semi-Lagrangian methods and $CFL= 0.85$ for finite difference methods. The final time is $T_{end} = 8$.</figcaption></figure>

A reference solution of kinetic energy is computed using a fifth order finite difference WENO method with very fine mesh ($n_x = 1025$, $\Delta t = 1/1600$).

In the following we take the parameter $\epsilon = 0.7$ and the computational domain is $(r, v) \in \Omega = [-4, 4]^2$.

Concerning the numerical resolution using semi-Lagrangian methods, we notice that we deliberately choose not to apply a time splitting in order to use this method in a general context. The characteristic curves corresponding to the Vlasov equation (4) cannot be solved explicitly. Then we apply a second order leap-frog scheme already developed in [27]. Finally, we interpolate the distribution function $f(r^n, v^n)$ by a tensor product for cubic spline or by a dimension by dimension Hermite WENO reconstruction.

In Figure 2 we compare the evolution of invariants (mass, $L^2$ norm, minimum $f$ the density) and the kinetic energy obtained from semi-Lagrangian and finite difference methods to the reference solution. We denote the linear phase for time interval $t \le 10$, and the nonlinear phase for time interval $t > 10$ where small filaments are generated.

As we can see from Figure 2, the kinetic energy obtained with the semi-Lagrangian method with cubic spline and Hermite WENO5 reconstruction is relatively close to the reference solution in the linear phase, but it diverges from the reference one
<figure><img src="image_2.png" /><figcaption>FIGURE 2. Simplified paraxial Vlasov-Poisson model : Comparison between semi-Lagrangian with cubic spline, Hermite WENO 5th and finite difference with Hermite WENO 5th methods for Beam test.</figcaption></figure>

in the nonlinear phase. Even with a fines mesh, we cannot improve the numerical results for large time.

Finally, we compare the distribution function $f$ obtained from semi-Lagrangian with cubic spline and Hermite WENO5 methods and finite difference with Hermite WENO5 reconstructions with a reference solution computed with a refined mesh ($\Delta t = 1/1600$ and $n_x = 1025$) in Figure 3. At time $t = 10$, the distribution function $f$ of semi-Lagrangian methods is very close to the reference solution, where two small filaments appear (see top of Figure 3). Then during the nonlinear phase, a large number of filaments are generated due to the non-linearity of the Vlasov-Poisson system and the distribution functions $f$ obtained with semi-Lagrangian and finite difference methods start to differ strongly at $t \ge 15$, which also correspond to the divergence of the kinetic energy of Figure 2. At time $t \ge 20$, the semi-Lagrangian methods generates a completely unstable beam who is not consistent with the results obtained for the reference solution.

These numerical simulations illustrate perfectly that semi-Lagrangian methods without splitting work well during the linear phase even with very large time step, but they do not seem very robust during the nonlinear phase when micro-structures appear (filamentation).
<figure><img src="image_3.png" /><figcaption>FIGURE 3. Simplified paraxial Vlasov-Poisson model : Distribution function for Beam test: (a) semi-Lagrangian with cubic spline; (b) semi-Lagrangian with Hermite WENO5; (c) finite difference with Hermite WENO5; (d) reference solution at time t = 10, 15 and 20. Mesh size is $n_x = 513$, $\Delta t = 1/800$.</figcaption></figure>

5.3. Guiding center model. We finally consider the guiding center model [3], which has been derived to describe highly magnetized plasma in the transverse plane

$$ (7) \qquad \begin{cases} \frac{\partial \rho}{\partial t} + \mathbf{U} \cdot \nabla \rho = 0, \\ -\Delta \phi = \rho. \end{cases} $$

where the velocity $\mathbf{U} = (-\partial_y\phi, \partial_x\phi)$. Here we consider the model in a disk domain

$$ D = \{(x, y) \in \mathbb{R}^2 : \sqrt{x^2 + y^2} \le R\} $$

and assume that the electric potential is vanishing at the boundary

$$ (8) \qquad \phi(x, y) = 0, \quad (x, y) \in \partial D. $$

Then if we ignore the effect of boundary conditions, the guiding center model verifies the following properties :

(1) Positivity of density $\rho$ : $0 \le \rho(t, x, y)$.

(2) Mass conservation : $\frac{d}{dt} \left( \int_D \rho dx dy \right) = 0$.

(3) $L^p$ norm conservation, for $1 \le p \le \infty$ :

$$ \frac{d}{dt} \|\rho\|_{L^p(D)} = 0. $$
(4) Energy conservation

$$ \frac{d}{dt} \left( \int_D |\nabla \phi|^2 dx dy \right) = 0. $$

To solve the system (7), we use a scheme based on Cartesian mesh and apply an Inverse Lax-Wendroff procedure to treat boundary conditions on kinetic equations [11, 12]. Actually, this system has been already solved in polar coordinates in [19]. But, the change of coordinate usually produce an artificial singularity at the origin, which needs a particular treatment. At contrast, with the Inverse Lax-Wendroff technique on Cartesian mesh, we do not have any singularity and it is not related to the numerical scheme since boundary effects and numerical schemes are treated independently. Furthermore, it is easy to adapt to other geometries [11, 12].

In this section, we only focus on discretization of boundary condition (8) of Poisson equation. The one for transport equation is trivial, since a homogeneous Dirichlet boundary condition will be used.

### 5.3.1. Discretization of Poisson equation.
A classical five points finite difference approximation is used to discretize the Poisson equation. However, to discretize the Laplacian operator $\Delta\phi$ near the physical boundary, we notice that some points of the usual five points finite difference formula can be located outside of interior domain. For instance, Figure 4 illustrates the discretization stencil for $\Delta\phi$ at the point $(x_i, y_j)$. The point $\mathbf{x}_g = (x_i, y_{j-1})$ is located outside of interior domain. Let us denote the approximation of $\phi$ at the point $\mathbf{x}_g$ by $\phi_{i,j-1}$. Thus $\phi_{i,j-1}$ should be extrapolated from the interior domain.

We extrapolate $\phi_{i,j-1}$ on the normal direction $\mathbf{n}$

$$ (9) \qquad \phi_{i,j-1} = \tilde{w}_p \phi(\mathbf{x}_p) + \tilde{w}_h \phi(\mathbf{x}_h) + \tilde{w}_{2h} \phi(\mathbf{x}_{2h}), $$

where $\mathbf{x}_p$ is the cross point of the normal $\mathbf{n}$ and the physical boundary $\Omega_{\mathbf{x}_\perp}$. The points $\mathbf{x}_h$ and $\mathbf{x}_{2h}$ are equal spacing on the normal $\mathbf{n}$, i.e. $h = |\mathbf{x}_p - \mathbf{x}_h| = |\mathbf{x}_h - \mathbf{x}_{2h}|$, with $h = \min(\Delta x, \Delta y)$, $\Delta x, \Delta y$ are the space steps in the directions $x$ and $y$ respectively. Moreover, $\tilde{w}_p, \tilde{w}_h, \tilde{w}_{2h}$ are the extrapolation weights depending on the position of $\mathbf{x}_g, \mathbf{x}_p, \mathbf{x}_h$ and $\mathbf{x}_{2h}$. In (9), $\phi(\mathbf{x}_p)$ is given by the boundary condition (8), whereas $\phi(\mathbf{x}_h), \phi(\mathbf{x}_{2h})$ should be determined by interpolation.

For this, we first construct an interpolation stencil $\mathcal{E}$, composed of grid points of $\Omega$. For instance, in Figure 4, the inward normal $\mathbf{n}$ intersects the grid lines $y = y_j, y_{j+1}, y_{j+2}$ at points $P_0^*, P_1^*, P_2^*$. Then we choose the three nearest points of the cross point $P_l^*$, $l = 0, 1, 2$, in each line, i.e. marked by a large circle. From these nine points, we construct a Lagrange polynomial $q_2(\mathbf{x}) \in \mathbb{Q}_2(\mathbb{R}^2)$ and evaluate the polynomial $q_2(\mathbf{x})$ at $\mathbf{x}_h$ and $\mathbf{x}_{2h}$, i.e.

$$ \phi(\mathbf{x}_h) = \sum_{\ell=0}^{8} w_{h,\ell} \phi(\mathbf{x}_\ell), $$

$$ \phi(\mathbf{x}_{2h}) = \sum_{\ell=0}^{8} w_{2h,\ell} \phi(\mathbf{x}_\ell), $$

with $\mathbf{x}_\ell \in \mathcal{E}$. Hence, we get $\phi_{i,j-1}$ which is approximated from the interior domain.

However, in some cases, we can not find a stencil of nine interior points. For instance, when the interior domain has small acute angle sharp, the normal $\mathbf{n}$ can not have three cross points $P_l^*$, $l = 0, 1, 2$ in interior domain, or we can not have three nearest points of the cross point $P_l^*$, $l = 0, 1, 2$, in each line. In this case, we alternatively use a first degree polynomial $q_1(\mathbf{x})$ with a four points stencil or even a zero degree polynomial $q_0(\mathbf{x})$ with an one point stencil. We can similarly construct the four points stencil or the one point stencil as above.
<figure><img src="image_4.png" /><figcaption>FIGURE 4. Spatially two-dimensional Cartesian mesh. ● is interior point, ■ is ghost point, □ is the point at the boundary, ○ is the point for extrapolation, the dashed line is the boundary.</figcaption></figure>

5.3.2. Numerical simulation of the diocotron instability. We now consider the diocotron instability for an annular electron layer. This plasma instability is created by two sheets of charge slipping past each other and is the analog of the Kelvin-Helmholtz instability in fluid mechanics. The initial data is given by

$$ \rho_0(\mathbf{x}_\perp) = \begin{cases} (1 + \varepsilon \cos(\ell\theta)) \exp(-4(r-6.5)^2), & \text{if } r^- \le \sqrt{x^2+y^2} \le r^+, \\ 0, & \text{otherwise,} \end{cases} $$

where $\varepsilon$ is a small parameter, $\theta = \operatorname{atan2}(y, x)$. In the following tests, we take $\varepsilon = 0.001$, $r^- = 5$, $r^+ = 8$, $\ell = 7$.

We have seen in the previous section that semi-Lagrangian methods may be not very appropriate. Indeed, the semi-Lagrangian method has some limitations during the nonlinear phase, when small filaments appear since a small time step must be used and the method is no more conservative. Therefore, we propose to apply a mixed method based on the Hermite interpolation with a WENO reconstruction: we use the semi-Lagrangian method for the linear phase with large time step; then we apply the conservative finite difference scheme for nonlinear phase with small time step respecting CFL condition. The criterion to pass from semi-Lagrangian to finite difference methods is as follows

$$ (10) \qquad \left| \int_{\mathbb{R}^2} [\rho_h(t_n) - \rho_h(t_{n-1})] dx \right| > h^3, $$

where $h$ is the smallest space step.

A comparison between the semi-Lagrangian with cubic spline method and the mixed method for the diocotron instability is presented in Figures 5 and 6. For a fair comparison an adaptive time step is also applied to the semi-Lagrangian
scheme. We first choose a CFL number $\lambda \approx 2$ during the linear phase, while we take $\lambda \approx 0.5$ for the nonlinear phase and the number of points in space is $n_x = n_y = 256$.

Using the semi-Lagrangian with cubic spline method, the relative error of mass oscillates a lot during the nonlinear phase, while relative error for the mixed method is more stable. We observe a very similar phenomenon for energy conservation. Unfortunately, relative error of $L^2$ norm for the mixed method is larger than that for the semi-Lagrangian with cubic spline method, but the price to pay is to generate a non negligible negative values. Indeed, the WENO reconstruction allows to control spurious oscillations generated from the discontinuous initial data and small structures.

Finally, the evolution of the density $\rho$ is presented in Figure 6. At first glance, we see the density of these two methods are very similar. At time $t = 40$, small filaments appear and then seven vortices are formed and move.

Looking more carefully, we observe that the numerical results obtained from the mixed method is a little bit more dissipative than the ones obtained from the semi-Lagrangian methods with cubic spline interpolation *i.e.* small structures of density are more thin. However, the semi-Lagrangian method is much more oscillatory than the mixed method, which can be observed from the minimum or maximum of density.

As a conclusion, although semi-Lagrangian method is less dissipative than the mixed HWENO method, it involves too much numerical instabilities in nonlinear phase. Therefore, the mixed method controlling spurious oscillations is more appropriate for long time simulation in plasma physics.

## 6. CONCLUSION AND PERSPECTIVE

In this paper, we have first developed a Hermite weighted essentially non-oscillatory reconstruction for semi-Lagrangian method and finite difference method respectively.

We illustrate that such a reconstruction is less dissipative than usual weighted essentially non-oscillatory reconstruction. Then we have compared our approach with the usual semi-Lagrangian method with cubic spline and finite difference WENO reconstruction. The semi-Lagrangian method is efficient and accurate for linear phase even with a large time step, however, it becomes less accurate for nonlinear phase and may lead to the wrong solution in some cases, for instance, the Beam test [3].

The finite difference method is stable under the classical CFL condition, but it is much more stable in nonlinear phase and it conserves mass. We thus apply a mixed method using the semi-Lagrangian method in linear phase and finite difference method during the nonlinear phase, called mixed HWENO5 method.

We finally apply the mixed HWENO5 method to the simulation of the diocotron instability and observe that although the mixed HWENO5 method is a little more dissipative than the semi-Lagrangian with cubic spline method, but it is much more stable during the nonlinear phase.

The next step is now to apply our mixed method to more realistic and high dimensional plasma turbulence simulations, for instance, 4D Drift-Kinetic simulation [15] or 5D Gyrokinetic simulation[14].

## ACKNOWLEDGMENT

The authors are partially supported by the European Research Council ERC Starting Grant 2009, project 239983-*NuSiKiMo*.
<figure><img src="image_5.png" /><figcaption>FIGURE 5. Guiding center model: *Comparison between semi-Lagrangian with cubic spline method and mixed semi-Lagrangian/finite difference Hermite WENO5 method*</figcaption></figure>

REFERENCES

[1] A. J. BRIZARD AND T.S HAHM, Foundations of nonlinear gyrokinetic theory, *Rev. Mod. Phys.* **79**, (2007) pp. 421-68.

[2] N. CROUSEILLES, M. MEHRENBERGER, E. SONNENDRUCKER, Conservative semi-Lagrangian schemes for Vlasov equations, *Journal of Computational Physics* **229** (2010), pp. 1927-1953.

[3] N. CROUSEILLES, M. LEMOU, F. MÉHATS, Asymptotic preserving schemes for highly oscillatory, *Journal of Computational Physics*, **248** (2013), pp. 287-308.

[4] B.A. DE DIOS, J.A. CARRILLO AND C.W. SHU, Discontinuous Galerkin methods for the multi-dimensional Vlasov-Poisson problem. *Math. Models Methods Appl. Sci.* **22** (2012), 1250042

[5] M.R. FEIX, P. BERTRAND, A. GHIZZO, Eulerian codes for the Vlasov equation, *Series on Advances in Mathematics for Applied Sciences, Kinetic Theory and Computing* **22**, (1994), pp. 45-81

[6] F. FILBET, Convergence of a Finite Volume Scheme for the One Dimensional Vlasov-Poisson System, *SIAM J. Numer. Analysis*, **39** (2001), pp. 1146-1169.

[7] F. FILBET, E. SONNENDRÜCKER, P. BERTRAND, Conservative numerical schemes for the Vlasov equation, *J. Comput. Phys.* **172** (2001), pp. 166-187

[8] F. FILBET, E. SONNENDRÜCKER, Comparison of Eulerian Vlasov solvers, *Compter Physics Communications*, **150** (2003), pp. 247-266.

[9] F. FILBET, E. SONNENDRÜCKER, Modeling and numerical simulation of space charge dominated beams in the paraxial approximation, *Mathematical Models and Methods in the Applied Sciences*, **16** (2006) pp. 763-791.

[10] F. FILBET, R. DUCLOUS AND B. DUBROCA. Analysis of a high order finite volume scheme for the Vlasov-Poisson system, *DCDS-S*, **5** (2012), pp. 283-305
<figure><img src="image_6.png" /><figcaption>FIGURE 6. Guiding center model : (a) semi-Lagrangian with cubic spline (b) mixed semi-Lagrangian/finite difference method with Hermite WENO5 at time $t = 40, 50$ and $60$.</figcaption></figure>
[11] F. FILBET AND C. YANG, An inverse Lax-Wendroff method for boundary conditions of Boltzmann equations, *Journal of Computational Physics*, **245**, (2013) pp. 43–61

[12] F. FILBET AND C. YANG, Numerical Simulations of Kinetic Models for Chemotaxis, submitted.

[13] A. GHIZZO, P. BERTRAND, M. SHOUCRI, T.W. JOHNSTON, E. FILJAKOW, M.R. FEIX, A Vlasov code for the numerical simulation of stimulated Raman scattering, *J. Comput. Phys.* **90** (1990) pp. 431–453.

[14] X. GARBET, Y. IDOMURA, L. VILLARD AND T.H. WATANABE, Gyrokinetic simulations of turbulent transport, *Nucl. Fusion* **50** (2010), 043002

[15] V. GRANDGIRARD, M. BRUNETTI, P. BERTRAND, N. BESSE, X. GARBET, P. GHENDRIH, G. MANFREDI, Y. SARAZIN, O. SAUTER, E. SONNENDRÜCKER, J. VACLAVIK, L. VILLARD, A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation, *Journal of Computational Physics*, **217** (2006), pp. 395–423.

[16] V. GRANDGIRARD et al. Computing ITG turbulence with a full-f semi-Lagrangian code Commun. Nonlinear Sci. Numer. Simul. **13** (2008) pp. 8187.

[17] W. GUO AND J.-M. QIU, Hybrid semi-Lagrangian finite element-finite difference methods for the Vlasov equation, *J. Comput. Phys.* **234** (2013), pp. 108-132.

[18] R. E. HEATH, I. M. GAMBA, P.J. MORRISON AND C. MICHLER, A discontinuous Galerkin method for the Vlasov-Poisson system *J. Comput. Phys.* **231** (2012), pp. 1140-1174.

[19] S. HIRSTOAGA, E. MADAULE, M. MEHRENBERGER AND J. PÉTRI, Semi-Lagrangian simulations of the diocotron instability, submitted

[20] G.-S. JIANG AND C.-W. SHU, Efficient implementation of weighted ENO schemes, *Journal of computational physics* **126**, (1996) pp. 202–228.

[21] CHENG, G. KNORR, *J. Comput. Phys.*, **22** (1976) pp. 330-348.

[22] T. NAKAMURA, T. YABE, Cubic interpolated propagation scheme for solving the hyper-dimensional Vlasov-Poisson equation in phase space. *Comput. Phys. Communications*, **120** (1999) pp. 122–154.

[23] J. QIU, CHI-WANG SHU, Hermite WENO schemes and their application as limiters for Runge-Kutta discontinuous Galerkin method: one-dimensional case, *Journal of Computational Physics*, **193**, (2013) pp. 115-135.

[24] J.M. Qiu and A. Christlieb, A conservative high order semi-Lagrangian WENO method for the Vlasov equation. *J. Comput. Phys.* **229** (2010), pp. 1130-1149.

[25] J.M. Qiu, C.-W. Shu, Conservative high order semi-Lagrangian finite difference WENO methods for advection in incompressible flow. *J. Comput. Phys.* **230** (2011), pp. 863-889.

[26] J. M. QIU AND C.-W. SHU, Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: theoretical analysis and application to the Vlasov-Poisson system. *J. Comput. Phys.* **230** (2011), pp. 8386-8409.

[27] E. SONNENDRÜCKER, J. ROCHE, The semi-Lagrangian method for the numerical resolution of the Vlasov equation, *J. Comput. Phys.*, **149** (1999) pp. 201–220.

[28] T. UMEDA, Y. NARIYUKI, D. KARIYA, A non-oscillatory and conservative semi-Lagrangian scheme with fourth-degree polynomial interpolation for solving the Vlasov equation *Comput. Phys. Commun.* **183** (2012), pp. 1094-1100.

FRANCIS FILBET

UNIVERSITÉ DE LYON,
UL1, INSAL, ECL, CNRS

UMR5208, INSTITUT CAMILLE JORDAN,
43 BOULEVARD 11 NOVEMBRE 1918,

F-69622 VILLEURBANNE CEDEX, FRANCE

CHANG YANG

UMR5208, INSTITUT CAMILLE JORDAN,
43 BOULEVARD 11 NOVEMBRE 1918,

VILLEURBANNE CEDEX, FRANCE

E-MAIL: filbet@math.univ-lyon1.fr

E-MAIL: yang@math.univ-lyon1.fr