Keywords: Schéma semi-Lagrangien, méthodes des volumes finis, méthodes conservatives, simulation numérique pour Vlasov, physique des plasmas

Méthodes semi-Lagrangiennes conservatives pour la résolution numérique de l'équation de Vlasov

Résumé : Des méthodes conservatives sont dveloppées pour la résolution numérique de l'équation de Vlasov dans le cadre de splitting directionnel unidimensionnel. Dans le cas de l'advection à coefficient constant, ces méthodes sont quivalentes aux méthodes semi-Lagrangiennes traditionnelles, mais pour les méthodes conservatives, il est possible d'ajouter des filtres qui garantissent la positivité de la solution. Dans le cas non-constant, les méthodes conservatives sont une alternative aux méthodes semi-Lagrangiennes classiques. Ces dernières peuvent prsenter des mauvaises conservations de la masse totale lorsqu'une technique de splitting est utilisée.

Introduction

To describe the dynamics of charged particles in a plasma or in a propagating beam, the Vlasov equation can be used to calculate the plasma response to the electromagnetic fields. The unknown f (t, x, v) which depends on the time t, the space x and the velocity v represents the distribution function of the studied particles. The coupling with the self-consistent electromagnetic fields is taken into account through the Maxwell or Poisson equation.

Due to its nonlinear structure, analytical solutions are available only in few academic cases, and numerical simulations have to be performed to study realistic physical phenomena. Nowadays, mostly two classes of methods are used to investigate the behaviour of the numerical solution to the Vlasov equation. On the one hand, Particle In Cell (PIC) methods, which are the most widely used, approach the plasma by macro-particles, the trajectories of which follow the characteristic curves of the Vlasov equation whereas the electromagnetic fields are computed by gathering the charge and current densities particles on a grid of the physical space (see [2]). On the other hand, Eulerian methods consist in discretizing the Vlasov equation on a grid of the phase space using classical numerical schemes such as finite volumes or finite elements methods for example (see [START_REF] Cheng | The integration of the Vlasov equation in configuration space[END_REF][START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF][START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF]).

Although PIC methods can theoretically and potentially resolve the whole 6 dimensional problem, it is well known that the inherent numerical noise makes difficult a precise description of low density regions, despite significant recent improvements. Hence, Eulerian methods offer a good alternative to overcome this lack of precision, even if problems of memory can arise when one deals with high dimensions. In particular, Vlasov codes seem to be appropriate to study nonlinear processes.

This last decade, gridded Vlasov solvers have been developed for 2D, 4D and even 5D phase space problems. Among them, the semi-Lagrangian method using a cubic spline interpolation (SPL) [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF] and the Positive Flux Conservative (PFC) method [START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF] have been implemented to deal with physical applications [START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF][START_REF] Ghizzo | A Hilbert-Vlasov code for the study of high-frequency plasma beatwave accelerator[END_REF][START_REF] Umeda | Comparison of numerical interpolation schemes for one-dimensional electrostatic Vlasov code[END_REF].

Recently, a parabolic spline method (PSM) has been introduced for transport equations arising in meteorology applications [START_REF] Zerroukat | A monotonic and positivedefinite filter for a Semi-Lagrangian Inherently Conserving and Efficient (SLICE) scheme[END_REF][START_REF] Zerroukat | The Parabolic Spline Method (PSM) for conservative transport problems[END_REF]. This mehod benefits from the best approximation property of the SPL method and from the conservation of mass and positivity (by applying a suitable filter) of the PFC method.

The aim of the present work is to study such a conservative method in the context of the Vlasov equation. Conservative methods present a lot of advantages. In addition to the inherent conservative property, slope limiters can be introduced in the reconstruction to ensure positivity; moreover, since they solve the conservative form of the equation, multi-dimensional problems can be solved by a splitting procedure so that the solution of the full problem is reduced to a succession of solution to only one-dimensional problems. Obviously, this property is of great interest from an implementation and algorithmic point of view.

We will focus here essentially on PSM, which has never been applied to our knowledge to Vlasov simulations. We will also introduce a new method based on a cubic splines approximation of the unknown; the characteristics curves are followed forwardly as in [START_REF] Crouseilles | A forward semi-Lagrangian scheme for the numerical solution of the Vlasov equation[END_REF], but the unknown is reconstructed in a conservative way using its values on the transported non-uniform mesh.

In our numerical experiments, we first consider the special case of directional splitting with constant advection (like the Vlasov-Poisson system). In that case, when no filter is applied, we prove that the advective scheme (e.g. SPL) and the conservative one (e.g. PSM) are equivalent. Note that in this setting, a mathematical proof of the convergence has been performed in [START_REF] Besse | Convergence of classes of high order semi-Lagrangian schemes for the Vlasov-Poisson system[END_REF]. We then discuss the choice of the filter in order to preserve the positivity. The filters introduced in [START_REF] Zerroukat | A monotonic and positivedefinite filter for a Semi-Lagrangian Inherently Conserving and Efficient (SLICE) scheme[END_REF][START_REF] Zerroukat | The Parabolic Spline Method (PSM) for conservative transport problems[END_REF] seem to modify too much the distribution function; we thus propose new filters which minimize in a certain sense the gap with the initial reconstruction while maintaining the positivity constraint. Numerical results are given for a classical plasma test case: the bump-on-tail instability.

We then focus on the case where we do not have a constant advection, as is the case for the guiding center model. In [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF], a 2D interpolation was proposed to approximate this model to ensure the conservation property (which is not true if a splitting procedure is used [START_REF] Huot | Instability of the time splitting scheme for the one-dimensional and relativistic Vlasov-Maxwell system[END_REF]). The time splitting in the conservative form has also been successfully tested for the PFC scheme [3], which is however more diffusive. Different works have been devoted to the study of the CIP method in its conservative form (see [START_REF] Nakamura | Exactly conservative semi-Lagrangian scheme for multi-dimensional hyperbolic equations with directional splitting technique[END_REF][START_REF] Nakamura | Cubic Interpolated Propagation Scheme for Solving the Hyper-Dimensional Vlasov-Poisson Equation in Phase Space[END_REF][START_REF] Tanaka | Constructing exactly conservative scheme in a non-conservative form[END_REF] and references therein). The time splitting in the advective form is often discarded since it can lead to bad mass conservation, especially for long time simulations ( [START_REF] Huot | Instability of the time splitting scheme for the one-dimensional and relativistic Vlasov-Maxwell system[END_REF][START_REF] Nakamura | Exactly conservative semi-Lagrangian scheme for multi-dimensional hyperbolic equations with directional splitting technique[END_REF]). We see here that with the conservative spline formulation (PSM), we can perform simulations with directional splitting not as diffusive as PFC, while maintaining the mass conservation. The time step is nevertheless limited by a CFL condition which seems luckily not to be so severe in our context. In particular, we obtain good numerical results for the classical test cases occuring in plasma physics.

The paper is organized as follows: first, semi-Lagrangian conservative methods are recalled and also introduced for one-dimensional general problems. Then, the constant advection case is investigated, and it is proved that, in this case, a conservative method and its advective counterparts are equivalent when no filter are considered; numerical results applied to the Vlasov-Poisson model are then discussed. Finally, we focus on the more interesting non-constant advection case for which numerical results illustrate the good behaviour of the new approaches.

Conservative methods

We are interested in the approximation of multi-dimensional transport equations of the form

∂g ∂t + ∇ x • (ag) = 0, x ∈ Ω ⊂ IR n , (1)
where the unknown g depends on time and on the multi-dimensional spatial direction x and a is a divergence free vector field ∇ x • a = 0 which can depend on time. The so-called conservative form (1) is then equivalent to the advective form ∂g ∂t

+ a • ∇ x g = 0, x ∈ Ω ⊂ IR n . (2)
In Vlasov type equations which enter in the class of equation of the form [START_REF] Besse | Convergence of classes of high order semi-Lagrangian schemes for the Vlasov-Poisson system[END_REF], Ω is a subset of the phase space which has up to 6 dimensions.

INRIA

Splitting the components of x into x 1 and x 2 , equation ( 1) can be written in the form ∂g ∂t

+ ∇ x1 • (a 1 g) + ∇ x2 • (a 2 g) = 0,
where a 1 and a 2 denote the component of the field a corresponding to x 1 and x 2 . It is well known (see [START_REF] Cheng | The integration of the Vlasov equation in configuration space[END_REF]) that a splitting procedure involves a successive solution of ∂g ∂t

+ ∇ x1 • (a 1 g) = 0, ∂g ∂t + ∇ x2 • (a 2 g) = 0, (3)
keeping high order accuracy in time for the whole equation [START_REF] Besse | Convergence of classes of high order semi-Lagrangian schemes for the Vlasov-Poisson system[END_REF]. However, the traditional semi-Lagrangian methods described in [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF] for example do not resolve the conservative form but the non-conservative form of the equations (2). Then, by solving only the advective form of (3), the corrective terms g∇ x1 • a 1 and g∇ x2 • a 2 are omitted and can lead to an important lack of accuracy in long time simulations (see [START_REF] Huot | Instability of the time splitting scheme for the one-dimensional and relativistic Vlasov-Maxwell system[END_REF][START_REF] Nakamura | Exactly conservative semi-Lagrangian scheme for multi-dimensional hyperbolic equations with directional splitting technique[END_REF]). An alternative way would be to solve the conservative form so that the solution of ( 1) can be performed by solving a succession of onedimensional problems. Hence in the sequel, we propose different conservative methods to solve one-dimensional problems; the methods are first presented in a general context but practical examples will be detailed in the next sections.

Conservative semi-Lagrangian methods for one-dimensional problems

The conservative methods enable to solve the conservative terms separately so that we restrict ourselves to the following one-dimensional problem

∂g ∂t + ∂(ag) ∂x = 0, x ∈ I ⊂ IR. ( 4
)
Let us consider a grid of the interval

I: x -1/2 < x 0 < x 1/2 < x 1 < ... < x N -1 < x N -1/2 < x N .
We consider the average quantity

ḡn i = 1 ∆x x i+1/2 x i-1/2 g(t n , x)dx, with ∆x = x i+1/2 -x i-1/2 , i = 0, .., N -1.
Thanks to the conservation of the volume, we can write the following equality

x n+1 i+1/2 x n+1 i-1/2 g(t n+1 , x)dx = x n i+1/2 x n i-1/2 g(t n , x)dx, (5)
where x n+1 i-1/2 and x n i-1/2 belong to the same characteristic curve defined by

dX(t) dt = a(t, X(t)), X(t n+1 ) = x n+1 i-1/2 , X(t n ) = x n i-1/2 , i = 0, .., N -1.
Assuming the values of g are known at time t n on each volume [x i-1/2 , x i+1/2 ], we can reconstruct the primitive function G(t n , x) of g(t n ) in I. Let us remark that the primitive function is only known on the mesh. To do this, we define G(t n ) as a cumulative function corresponding to ḡn :

G(t n , x i-1/2 ) = ∆x i-1 k=0 ḡn k , i = 0, .., N. (6)
Hence, using (5), we have

ḡn+1 i = x n+1 i+1/2 x n+1 i-1/2 g(t n+1 , x)dx = x n i+1/2 x n i-1/2 g(t n , x)dx = G(t n , x n i+1/2 ) -G(t n , x n i-1/2 ).
The only thing to do is to build a polynomial primitive function, assuming G(t n ) is known on the mesh (x i-1/2 ) i=0,..,N -1 . Several ways can be used to achieve a good approximation. Hence, as in the pointwize semi-Lagrangian method, the algorithm of conservative methods is composed of two main steps: the computation of the characteristic curves, and the reconstruction step.

Computation of the characteristic curves

In the semi-Lagrangian method, we have to compute the characteristics curves between two consecutive time steps. As remarked in [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF], second order accuracy can be reached using a two steps approach. Indeed, if we assume that g n-1 and g n are known, we can compute the advection term a(t n , X(t n )) which can depend on the solution g n . Let us remark that the advection term is assumed constant in time between t n-1 and t n+1 . Then, the equation to solve writes

dX(t) dt = a(t n , X(t)), (7)
with a final condition which coincides with the mesh point X i

X(t n+1 ) = X i = x 0 + i∆x. (8)
Here, [START_REF] Crouseilles | A forward semi-Lagrangian scheme for the numerical solution of the Vlasov equation[END_REF] is solved using a parabolic assumption (see [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF][START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF]). Let X i be the position of X(t n+1 ), then there exists a displacement α i ∈ IR such that

X n i = X i -α i , and X n-1 i = X i -2α i .
Hence, the displacement can be computed at second order by solving the following one-dimensional fixed-point

α i = ∆t a(t n , X i -α i ). (9)
In [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF], a Newton algorithm is used but every iterative methods can be employed. We also mention [START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF] in which a Taylor expansion of the right hand side of (9) is performed; this strategy is equivalent to a Newton algorithm in which two iterations are imposed. However, the drawback of these algorithms is that they require the evaluation of the Jacobian matrix of a(t n ). A fixed point algorithm can then be implemented. But, if we assume linear reconstruction of the advection term at points (X i -α i ) (as it is supposed in [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF][START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF]), an explicit algorithm can be used. The main steps of this new algorithm is detailed in the sequel.

INRIA

Starting from [START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF] and denoting by [x j , x j+1 ] the cell in which X i -α i falls, the linear reconstruction of a(t n ) writes

α i = ∆t [(1 -β)a(t n , x j ) + βa(t n , x j+1 )] , (10)
where β is such that

X i -α i = x j + β, x j = x 0 + j∆x, X i = x 0 + i∆x. (11)
Injecting the expression of α i into (10) leads

β [∆x + ∆t (a(t n , x j+1 ) -a(t n , x j ))] = (i -j)∆x -∆t a(t n , x j ), (12)
from which an expression of β can be deduced

β = [(i -j)∆x -∆t a(t n , x j )] / [∆x + ∆t(a(t n , x j+1 ) -a(t n , x j ))] . (13)
Now, it remains to determine the j index. To do that, it must be remarked that β given by ( 13) lives in the interval [0, 1]. Hence, from [START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF], we can deduce an expression for

x i = i∆x i∆x = j∆x + ∆t a(t n , x j ) + β [∆x + ∆t(a(t n , x j+1 ) -a(t n , x j ))] .
Using the fact that β ∈ [0, 1], and by remarking that [∆x + ∆t(a(t n , x j+1 ) -a(t n , x j ))] > 0 provided that ∆t is small enough, we deduce i∆x ∈ [M j , M j+1 ] , with M j = x j + ∆t a(t n , x j ).

Under the assumption that ∆t is small enough, the non-decreasing sequence (M j ) j=0,.,N -1 forms a non-uniform mesh from which it can be easily found the location of X i . The algorithm is then the following for each i = 0, .., N -1:

• determination of j such that X i ∈ [M j , M j+1 ]
• determination of β with (13)

• determination of α i with [START_REF] Ghizzo | A Hilbert-Vlasov code for the study of high-frequency plasma beatwave accelerator[END_REF] This algorithm has been proved to be faster than classical iterative based methods. Obviously, it leads to the same displacement α i .

Formulation by primitive

In this subsection, we present some reconstructions to update the unknown; once the feet of the characteristics is computed, we have to evaluate the integral of the function g(t n ). Two approaches to do that are presented in the sequel.

Lagrange reconstruction of the primitive. This paragraph presents a Lagrange reconstruction of the primitive, as done in [START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF] for example. The primitive is known on the mesh x i-1/2 , i = 1, ..., N , so that the cubic Lagrange polynomial can be written

G(t n , x) = i+2 j=i-1 G(t n , x j-1/2 )L j (x), ∀x ∈ [x i-1/2 , x i+1/2 ], (14)
where

L j (x) = j+2 k=j-1,k =j (x -x k-1/2 )/(x j-1/2 -x k-1/2 ).
This reconstruction corresponds to the PFC method introduced in [START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF] in which the slope limiters step is not performed. This approach and similar ones has been also introduced by [START_REF] Laprise | A class of semi-Lagrangian integrated-mass (SLIM) numerical transport algorithms[END_REF][START_REF] Colella | The piecewise parabolic method (PPM) for gas-dynamical simulations[END_REF]. See [START_REF] Lauritzen | An inherently mass-conservative semi-implicit semi-Lagrangian model[END_REF] for a more complete bibliography.

PSM: Cubic splines reconstruction of the primitive. We can also build a cubic spline for the interpolation conditions of which are given by ( 6). This can be written as

G(t n , x) = N -1 i=0 η i S x -x i-1/2 ∆x , (15)
where η i are the coefficients and S the cubic splines

6S(x) =    (2 -|x|) 3 if 1 ≤ |x| ≤ 2, 4 -6x 2 + 3|x| 3 if 0 ≤ |x| ≤ 1, 0 otherwise,
In addition to the interpolation conditions, we close the system satisfied by the cubic splines coefficients η i by considering the periodic case. But, even if ḡn is periodic (ḡ n 0 = ḡn N ), let us mention that it is not true for the cumulative function. However, we have the following relation on the primitive for

x > x N -1/2 G(t n , x) = x x -1/2 g(t n , x)dx = x N -1/2 x -1/2 g(t n , x)dx + x x N -1/2 g(t n , x)dx = M + x x N -1/2 g(t n , x)dx = M + G(t n , x),
where M denotes the total mass N -1 k=0 ḡn k . Then, using this last property, a periodic linear system can be recovered

Aη =             4 1 0 0 • • • 1 1 4 1 0 . . . 0 1 4 1 . . . . . . . . . . . . . . . . . . . . . 0 . . . 0 1 4 1 1 0 0 0 1 4                       η 0 η 1 . . . . . . η N -2 η N -1           = 6∆x           ḡ0 + M ḡ0 + ḡ1 . . . . . . N -2 i=0 ḡi N -1 i=0 ḡi -M          
The coefficients η i , i / ∈ [0, N -1] are deduced from the solutions of the previous linear system by

η -i = η -i+N -M, ∀i > 0, η i+N -1 = η i + M, ∀i ∈ [0, N -1].
This approach has been introduced in [START_REF] Zerroukat | The Parabolic Spline Method (PSM) for conservative transport problems[END_REF] as the Parabolic Spline Method. Their formulation refers to the reconstruction of the function g which is a C 1 piecewise parabolic function. The two formulations (by using primitive G or the function g) are completely equivalent, as already explained in [START_REF] Zerroukat | The Parabolic Spline Method (PSM) for conservative transport problems[END_REF].

INRIA

Forward update of the characteristics

Another strategy to update in a conservative way the unknowns consists in advancing in time the mesh points which are denoted by x n i-1/2 , i = 0, .., N . The ends of the characteristics starting at the mesh x n i-1/2 are called x n+1 i-1/2 , i = 0, .., N . These points form a non-uniform mesh of the domain in the general case. Then, due to the property of the conservation of the volumes [START_REF] Colella | The piecewise parabolic method (PPM) for gas-dynamical simulations[END_REF], it is possible to reconstruct a primitive function G(t n+1 ) of g(t n+1 ) with cubic spline polynomials in the spirit of paragraph 2.1.2. The main differences of the present case are twofold. On the one side, the mesh on which G(t n+1 ) is known and has to be reconstructed is not uniform, and on the other side, the characteristic curves are advanced in time.

In the rest of the section, the two steps of the method are detailed.

Computation of the characteristics curves

We have to compute the characteristics curves between two consecutive time steps. As previously, if we assume that g n-1 and g n are known, we can compute the advection term a(t n , X(t n )) which can depend on the solution g n . Let us remark that the advection term is supposed constant in time between t n-1 and t n+1 so that the equation to solve writes

dX(t) dt = a(t n , X(t)), (16)
with an initial condition which coicides with a mesh point X i in the considered case

X(t n-1 ) = X i = x 0 + i∆x. (17)
Even if the forward non-uniform approach can be expensive from a numerical point of view due to the solution of a new linear system at each time step, this algorithm has the opportunity to compute the characteristics curves in an explicit way. Indeed, contrary to the traditional backward method, the initial condition of ( 16) is given at the initial time. Hence, explicit algorithms can be implemented like Runge-Kutta ones. In our practical experiments, methods up to the fourth order Runge-Kutta algorithm have been implemented even if the formulation ( 9) is only second order accurate by essence. Let us recall the Runge-Kutta 2 method applied to our problem (we suppose that the advection field remains constant in time during the following steps, equal to its value at time t n )

k 1 = Ia(t n , X(t n-1 )), k 2 = Ia(t n , X(t n-1 ) + 2∆t k 1 ),
which leads to the following approximation of the characteristics at time t n+1

X(t n+1 ) = X(t n ) + 2∆t 2 (k 1 + k 2 ).
We denote by I an interpolation operator. In our experiments, a cubic spline polynomial has been used. Note that in the context of the time splitting studied in the guiding center model, we have also successfully applied Runge Kutta schemes for the backward method.

Interpolation using a non-uniform mesh

The so-reconstructed primitive G can then be interpolated on the uniform mesh x i-1/2 , i = 0, .., N to obtain the update unknown in the following way

G(t n+1 , x i+1/2 ) -G(t n+1 , x i-1/2 ) = 1 ∆x x i+1/2 x i-1/2 g n+1 (x)dx = ḡn+1 i .
As we said, the main step consists in the computation of the cubic splines coefficients on a non-uniform and periodic mesh. Cubic splines relations can be obtained by a repetition of convolution of the zero-th order splines (the box function) with itself (see [START_REF] De Boor | A practical guide to splines[END_REF]). Restricting to third order, we are faced to the following tridiagonal system to solve

Aη =             D 0 C 1 0 0 • • • B N -1 B 0 D 1 C 2 0 . . . 0 B 1 D 2 C 3 . . . . . . . . . . . . . . . . . . . . . 0 . . . 0 B N -3 D N -2 C N -2 C 1 0 0 0 B N -2 D N -1                       η 0 η 1 . . . . . . η N -2 η N -1           =           ḡ0 + B N -1 M ḡ0 + ḡ1 . . . . . . N -2 i=0 ḡi N -1 i=0 ḡi -C 0 M          
where the components of the matrix are defined for i = 0, .., N -1

B i = ∆x 2 i+1 /[(∆x i+1 + ∆x i + ∆x i-1 )(∆x i+1 + ∆x i )] D i = (∆x i-2 + ∆x i-1 )∆x i /[(∆x i + ∆x i-1 + ∆x i-2 )(∆x i + ∆x i-1 )] + (∆x i+1 + ∆x i )∆x i /[(∆x i+1 + ∆x i + ∆x i-1 )(∆x i + ∆x i-1 )] C i = ∆x 2 i-2 /[(∆x i + ∆x i-1 + ∆x i-2 )(∆x i-1 + ∆x i-2 )]
, and ∆x i is defined as follows: ∆x i = x n+1 i+1/2 -x n+1 i-1/2 , i = 0, .., N -1 and defined for all i by periodicity. As in subsection 2.1.2, the periodicity can be taken into account in the linear system, even if the primitive is not periodic.

Let us mention that a good behaviour of the present method requires a good approximation of the characteristics curves x n+1 i+1/2 at time t n+1 . Indeed, this strong dependence can be explained by the fact that they refer to the conditions of interpolation.

Slope limiters

In this subsection, we focus on the description of different filters which can be adapted to the previous reconstruction. Various approches have already been introduced (see [START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF][START_REF] Coquel | Second order entropy diminishing scheme for the Euler equations[END_REF][START_REF] Zerroukat | The Parabolic Spline Method (PSM) for conservative transport problems[END_REF]). We discuss here the use of specific filters for the PFC and PSM methods.

The PFC method In order to ensure the positivity, we are faced to the following problem : given ḡi-1 , ḡi , ḡi+1 , we have to find a reconstruction P = P (ḡ i-1 , ḡi , ḡi+1 ) which satisfies 1 ∆x The reconstruction on (x i-1/2 , x i+1/2 ) with the PFC method without the filter is the unique polynomial of degree ≤ 2 that satisfies 1 ∆x

x i+1/2 x i-1/2 P (x)dx = ḡi , P (x) ≥ 0, x ∈ (x i-1/2 , x i+1/2 ). (18)
x i-1/2 x i-3/2 P (x)dx = ḡi-1 , 1 ∆x x i+1/2 x i-1/2 P (x)dx = ḡi-1 , 1 ∆x x i+3/2
x i+1/2 P (x)dx = ḡi+1 .

By explicit computations, we have

P (x i-1/2 +ω∆x) = ḡi -(ḡ i -ḡ i-1 )((1-ω) 2 -1/3)/2+(ḡ i+1 -ḡ i )(ω 2 -1/3)/2, 0 ≤ ω ≤ 1.
We define ḡmin = min i=0,...,N -1

ḡi , ḡmax = max i=0,...,N -1 ḡi and suppose that ḡmin ≥ 0, which should be the case for our initial data.

We propose here to take the following filter, which is a little less restrictive version of the one used in [START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF] P

(x i-1/2 +ω∆x) = ḡi -ǫ -(ḡ i -ḡ i-1 )((1-ω) 2 -1/3)/2+ǫ + (ḡ i+1 -ḡ i )(ω 2 -1/3)/2, 0 ≤ ω ≤ 1, with (ǫ -, ǫ + ) =                                min(1, 4 ḡmax -ḡi ḡi -ḡi-1 ), min(1, 4 ḡi -ḡmin ḡi+1 -ḡi ) , if ḡi > ḡi-1 and ḡi+1 > ḡi , min(1, 3 ḡmax -ḡi ḡi -ḡi-1 ), min(1, 3 ḡmax -ḡi ḡi -ḡi+1 ) , if ḡi > ḡi-1 and ḡi+1 < ḡi , min(1, 4 ḡi -ḡmin ḡi-1 -ḡi ), min(1, 4 ḡmax -ḡi ḡi -ḡi+1 ) , if ḡi < ḡi-1 and ḡi+1 < ḡi , min(1, 3 ḡi -ḡmin ḡi-1 -ḡi ), min(1, 3 ḡi -ḡmin ḡi+1 -ḡi ) , if ḡi < ḡi-1 and ḡi+1 > ḡi . ( 19
)
The reconstruction then satisfies [START_REF] Nakamura | Exactly conservative semi-Lagrangian scheme for multi-dimensional hyperbolic equations with directional splitting technique[END_REF] and even ḡmin ≤ P (x) ≤ ḡmax , for x ∈ (x i-1/2 , x i+1/2 ). We will refer to this method as the PFC2 method (see Figs. The PSM method We denote by P i the reconstruction of PSM without the filter on the interval (x i-1/2 , x i+1/2 ). Explicit computations give for 0 < ω < 1,

P i (x i-1/2 +ω∆x) = g i-1/2 +ω(g i+1/2 -g i-1/2 )+(6ḡ i -3(g i-1/2 +g i+1/2 ))ω(1-ω), RR n°6856
where we have set g i+1/2 = P i (x i+1/2 ) = P i+1 (x i+1/2 ), for i = 0, . . . , N -1. In order to satisfy the positivity for P i (x), we proceed in two steps, like in [START_REF] Zerroukat | A monotonic and positivedefinite filter for a Semi-Lagrangian Inherently Conserving and Efficient (SLICE) scheme[END_REF][START_REF] Zerroukat | The Parabolic Spline Method (PSM) for conservative transport problems[END_REF], but we want to relax the filter used there which may be too severe in our context. 1. We take g new i+1/2 = max(ḡ min , min(ḡ max , g i+1/2 )), so that we have ḡmin ≤ g new i+1/2 ≤ ḡmax , ∀i. We then consider the first reconstruction for ω ∈ [0, 1]

R 1 (x i-1/2 +ω∆x) = g new i-1/2 +ω(g new i+1/2 -g new i-1/2 )+(6ḡ i -3(g new i-1/2 +g new i+1/2 ))ω(1-ω), 2. If R 1 (x i-1/2
+ ω∆x) has a strictly negative value for some ω 0 ∈ (0, 1), we consider for 0 ≤ β ≤ 1 the unique polynomial R 2,β of degree ≤ 2 which satisfies

x i+1/2 x i-1/2 R 2,β (x)dx = ḡi , R 2,β (x i-1/2 + β∆x) = 0, R ′ 2,β (x i-1/2 + β∆x) = 0.
(20) Clearly, this polynomial remains positive for x ∈ (x i-1/2 , x i+1/2 ). Now, among all β ∈ [0, 1] which satisfy [START_REF] Shoucri | Nonlinear evolution of the bump-on-tail instability[END_REF], we want to choose the closest polynomial R 2,β to R 1 in a certain sense. We choose here β * that minimizes a functional J(β). This functional can be chosen according to a property of interest which we want to respect as best as possible. For example, if we want to reconstruct a polynomial function R 2,β defined as in [START_REF] Shoucri | Nonlinear evolution of the bump-on-tail instability[END_REF], the boundary cell values of which are as close as possible to g new i±1/2 , we consider the following functional

J(β) = |g new i-1/2 -R 2,β (x i-1/2 )| + |g new i+1/2 -R 2,β (x i+1/2 )|. ( 21
)
The minimum of this functional is identified at some β ⋆ and thus we take R = R 2,β * . In the case where R 1 (x i-1/2 + ω∆x) has no strictly negative value, for ω 0 ∈ (0, 1), we do no other change, that is, we take R = R 1 . Hence the whole reconstruction satisfies the positivity of the polynomial reconstruction (see Figs.

2 in which for we compare the filter of [START_REF] Zerroukat | A monotonic and positivedefinite filter for a Semi-Lagrangian Inherently Conserving and Efficient (SLICE) scheme[END_REF] (PSM ZER), our optimal filter and the corresponding cubic splines). Note that the approach is general and we can use other functionals. We can for example optimize the following functional

J(β) = x i+1/2 x i-1/2 x n (R 1 (x) -R 2,β (x)) dx , n = 1, 2.
The functional reaches its minimum at β = β ⋆ . When n = 1, the R 2,β ⋆ reconstruction is as close as possible to the R 1 reconstruction regarding the first momentum which is preserved by the R 1 reconstruction. It is also possible to consider higher moments through the functional (n = 2 for example), but in these cases, the first momentum is not well preserved in our numerical tests, which pollutes the others diagnostics.

More generally, it is possible to dissociate the first reconstruction R 1 (which can be done with a Lagrange polynomial as for the PFC method) from the second one. We could also minimize in the L 2 norm (see e.g. [START_REF] Coquel | Second order entropy diminishing scheme for the Euler equations[END_REF] for an analog of [START_REF] Shoucri | A two-level implicit scheme for the numerical solution of the linearized vorticity equation[END_REF] with squares instead of absolute values). Generally, it is simpler to deal with squares than with absolute values for the search of a minimum. In that case, it turns out that it is the opposite case : for the squares, we have to search the roots of polynomials of degree 4 and only of degree 2 for absolute values.

The constant advection case

In this section, we are concerned with a constant advection field, so (4) can be rewritten equivalently in a conservative form

∂g ∂t + ∂(ag) ∂x = 0, x ∈ I ⊂ IR,
or in an advective form ∂g ∂t

+ a ∂g ∂x = 0, x ∈ I ⊂ IR.
Let us define for α ∈ R and N ∈ N * a transport operator T α,N : R N → R N . We will define several possibilities for the transport operator. If (g 0 , . . . , g N -1 ) is a discretization of a function g, then T α,N (g 0 , . . . , g N -1 ) should be a discretization of the shifted function x → g(x + α). In our case α denotes the displacement given by the solution of [START_REF] Crouseilles | A forward semi-Lagrangian scheme for the numerical solution of the Vlasov equation[END_REF] in which a is considered constant (in time and in space); we deduce that α = -∆t a.

The advective approach

We reconstruct a function ga which satisfies

ga (x i ) = g i , i = 0, . . . , N -1, (1)
and we take T α,N (g 0 , . . . , g N -1

) i = ga (x i + α). (2)
We suppose periodic boundary conditions, that is g i = g i+kN for k ∈ Z. We give here some standard interpolating functions g (see e.g. [START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF][START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF]).

Lagrange interpolation of order

d = 2p + 1 (Lag d ) : if x ∈ (x k , x k+1
), we take ga (x) = P (x), where P is the unique interpolating polynomial of degree

≤ d such that P (x i ) = g i , for i = k -p, . . . , k + p + 1.
Cubic spline interpolation (SPL) : we define the cubic spline S by

6S(x) =    (2 -|x|) 3 if 1 ≤ |x| ≤ 2, 4 -6x 2 + 3|x| 3 if 0 ≤ |x| ≤ 1, 0 otherwise,
and take ga (x) = N -1 i=0 η i S( x-xi ∆x ), where the coefficients η i , i = 0, . . . , N -1 are determined such that ga (x i ) = g i , for i = 0, . . . , N -1.

RR n°6856

The conservative approach

We reconstruct a function gc which satisfies 1 ∆x

x i+1/2 x i-1/2 gc (x)dx = ḡi , i = 0, . . . , N -1, (3)
and we then take

T α (ḡ 0 , . . . , ḡN-1 ) i = 1 ∆x x i+1/2 +α x i-1/2 +α gc (x)dx. ( 4
)
We have again several possibilities for the reconstruction of the function gc . We suppose here also periodic boundary conditions, that is ḡi = ḡi+kN for k ∈ Z.

Lagrange type interpolation of order

d = 2p (CLag d ) : if x ∈ (x k-1/2 , x k+1/2
), we take gc (x) = P (x), where P is the unique polynomial of degree ≤ d satisfying 1 ∆x

x i+1/2 x i-1/2 P (x)dx = ḡi , for i = k -p, . . . , k + p. (5)
Parabolic splines (CSPL 2 ) : it has already been previously described as the PSM method without limiters. We shortly recall one of its definition, in this context of notations. We define gc (x) = P i (x), for x ∈ (x i-1/2 , x i+1/2 ), where the polynomials P i , i = 0, . . . , N -1 are of degree ≤ 2 uniquely determined by the conditions 1 ∆x x i+1/2

x i-1/2 P i (x)dx = ḡi , i = 0, . . . , N -1, together with the continuity conditions

P i (x i+1/2 ) = P i+1 (x i+1/2 ), P ′ i (x i+1/2 ) = P ′ i+1 (x i+1/2 ), i = 0, . . . , N -1.
Note that CLag 2 is in fact the PFC method without limiters [START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF] and, as we have mentionned it just before, that CSPL 2 is PSM without limiters [START_REF] Zerroukat | A monotonic and positivedefinite filter for a Semi-Lagrangian Inherently Conserving and Efficient (SLICE) scheme[END_REF][START_REF] Zerroukat | The Parabolic Spline Method (PSM) for conservative transport problems[END_REF]. In the sequel we will make the link between the conservative and advective approaches and we will prove that CLag 2 is equivalent to LAG3 and that CSPL 2 is equivalent to SPL. We underline that this equivalence only holds for uniform meshes in a constant advection case.

Equivalence between conservative and advective approach

One interesting property is the following: the schemes deriving from CLag 2p or from Lag 2p+1 give rise to the same transport operator, and the same is true for CSPL 2 and SPL. We underline that this property only holds in this particular case of a constant 1D advection with a uniform mesh. On the other hand, we can make an equivalence between a conservative and an advective approach for a broad range of reconstructions. In the sequel, we will formalize this property and then check on examples for which we have the equivalence CLag 2p ∼ Lag 2p+1 and CSPL 2 ∼ SPL. We suppose for this that the transport operator has the following form :

T α,N (g 0 , . . . , g N -1 ) i = s ℓ=-s L ℓ (x i + α -x k )g k+ℓ , (6)

INRIA

where k is determined such that x i + α falls in [x k , x k+1 [ and s ∈ N * . The number s may depend on N . As examples, we will see that the reconstructions Lag 2p+1 and SPL fall in that category.

From advective to conservative. Now, if the transport operator comes from an advective approach, the reconstruction ga satisfies for 0 ≤ α < ∆x ga (

x i + α) = s ℓ=-s L ℓ (α)g i+ℓ
In particular, we have from the interpolation property (1), by taking α = 0,

g i = s ℓ=-s L ℓ (0)g i+ℓ , i = 0, . . . , N -1,
and thus we have

L 0 (0) = 1, L ℓ (0) = 0, if ℓ = 0. (7)
In order to make the link with the conservative approach, we want to define a reconstruction gc which satisfies (3) and which gives rise to the same transport operator. For this, we define for α ∈ [0, ∆x[ and i = 0, . . . , N -1,

G(x i+1/2 + α) = ∆x s ℓ=-s L ℓ (α) i+ℓ r=-s g r (8) and gc
(x) = G ′ (x), x 1/2 ≤ x < x N +1/2 , (9)
and we prolongate gc on R by periodicity. We then have for i = 1, . . . , N , by using (7)

1 ∆x

x i+1/2 x i-1/2 gc (x)dx = G(x i+1/2 ) -G(x i-1/2 ) ∆x = s ℓ=-s L ℓ (0) i+ℓ r=-s g r - s ℓ=-s L ℓ (0) i+ℓ-1 r=-s g r = i r=-s g r - i-1 r=-s g r = g i ,
and for i = 0, we get 1 ∆x

x 1/2 x -1/2 gc (x)dx = 1 ∆x x N +1/2 x N -1/2 gc (x)dx = g N = g 0 ,
and thus gc effectively satisfies (3). On the other hand, we have for α ∈ R and i = 0, . . . , N -1, 1 ∆x

x i+1/2 +α x i-1/2 +α gc (x)dx = G(x k+1/2 + x i-k + α) -G(x k-1/2 + x i-k + α) ∆x = s ℓ=-s L ℓ (x i-k +α) k+ℓ r=-s g r - s ℓ=-s L ℓ (x i-k +α) k-1+ℓ r=-s g r = s ℓ=-s L ℓ (x i-k +α)g k+ℓ ,
where k is determined such that x i + α falls in [x k , x k+1 [, and thus, we have T α,N (g 0 , . . . , g N -1 ) = ga (x i + α) = 1 ∆x

x i+1/2 +α x i-1/2 +α gc (x)dx,
which means that the conservative reconstruction gc derived from the advective reconstruction defines the same transport operator.

From conservative to advective. We suppose this time that the transport operator comes from a conservative approach. The reconstruction gc satisfies 1 ∆x

x i+1/2 +α x i-1/2 +α gc (x)dx = s ℓ=-s L ℓ (x i + α -x k )g k+ℓ , x i + α ∈ [x k , x k+1 [,
together with (3). In order to make the link with the advective approach, we simply define ga by ga (

x i + α) = s ℓ=-s L ℓ (x i + α -x k )g k+ℓ , x i + α ∈ [x k , x k+1 [, (10)
and it remains to prove that we have [START_REF] Crouseilles | A forward semi-Lagrangian scheme for the numerical solution of the Vlasov equation[END_REF], which comes straightforwardly from (3).

Equivalence for Lagrange reconstructions (CLag 2p ∼ Lag 2p+1 ).

From the advective approach, we have explicitely [START_REF] Cheng | The integration of the Vlasov equation in configuration space[END_REF], with

s = p + 1, L ℓ (α) = p+1 j=-p-1, j =ℓ α -j∆x (ℓ -j)∆x .
We then define gc by ( 8)- [START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF]. By putting P (x) = gc (x), for x ∈ (x k-1/2 , x k+1/2 ), we have to check that P is of degree ≤ p and satisfies [START_REF] Colella | The piecewise parabolic method (PPM) for gas-dynamical simulations[END_REF]. From ( 8)-( 9), we get

P (x) = ∆x p+1 ℓ=-p-1 L ′ ℓ (x -x k-1/2 ) k-1+ℓ r=-p-1 g r ,
and thus P is of degree ≤ p. We also have for ℓ = k -p, . . . , k + p 1 ∆x

x ℓ+1/2 x ℓ-1/2 P (x)dx = ∆x p+1 j=-p-1 L j (x ℓ+1/2 -x k-1/2 ) -L j (x ℓ-1/2 -x k-1/2 ) k-1+j r=-p-1 g r .
On the right hand side, we have

L j (x ℓ-1/2 -x k-1/2 ) = L j ((ℓ -k)∆x) = 0, if j = ℓ -k, since ℓ -k ∈ {-p -1, . . . , p + 1}, and L j (x ℓ-1/2 -x k-1/2 ) = 1, if j = ℓ -k. We similarly have L j (x ℓ+1/2 -x k-1/2 ) = L j ((ℓ + 1 -k)∆x) = 0, if j = ℓ+1-k, since ℓ+1-k ∈ {-p-1, . . . , p+1}, and L j (x ℓ+1-1/2 -x k-1/2 ) = 1, if j = ℓ + 1 -k. We thus get 1 ∆x x ℓ+1/2 x ℓ-1/2 P (x)dx = ∆x k-1+ℓ+1-k r=-p-1 g r -∆x k-1+ℓ-k r=-p-1 g r = g ℓ ,
which means that (5) is true.

Equivalence for spline reconstructions (CSPL 2 ∼ SPL).

From the advective approach, we have ga (

x i + α) = N -1 ℓ=0 η ℓ S(α/∆x + i -ℓ),

INRIA

and the η ℓ are determined by the conditions

η ℓ-1 + 4η ℓ + η ℓ+1 = 6g ℓ , α ℓ = α ℓ mod N , ℓ ∈ Z.
Following [START_REF] Kershaw | The explicit inverse of two commonly occuring matrices[END_REF], we can obtain explicitely for ℓ = 0, . . . , N -1

η ℓ = N -1 r=0 u N -|r-ℓ|-1 + u |r-ℓ|-1 2(1 -u N ) 6g r ,
by introducing

u N = (-2 + √ 3) N + (-2 - √ 3) N 2 , v r = (-2 + √ 3) r+1 + (-2 - √ 3) r+1 2 √ 3 , r = 0, . . . , N -1.
From the periodicity of the sequence (g i ), we get for ℓ ∈ Z

η ℓ = N/2-1 s=-N/2 v N -|s|-1 + v |s|-1 2(1 -u N ) 6g ℓ+s .
This finally gives ga (

x i + α) = N/2-1 j=-N/2 N -1 p=0 v N -|j|-1 + v |j|-1 2(1 -u N ) 6S(α/∆x + i -p)g p+j
An example where the advective approach is not conservative. Let N = 2p ∈ N * be an even number. We consider the reconstruction ga defined by

ga (x 2i + α∆x) = g 2i + (g 2i+1 -g 2i )α + g 2i+2 -2g 2i+1 + g 2i 2 α(α -1), 0 ≤ α < 2.
We then have for 0 ≤ α < 1

N -1 i=0 T α,N (g 0 , . . . , g N -1 ) i = p-1 i=0 (g 2i + (g 2i+1 -g 2i )α + g 2i+2 -2g 2i+1 + g 2i 2 α(α -1) + p-1 i=0 (g 2i + (g 2i+1 -g 2i )(α + 1) + g 2i+2 -2g 2i+1 + g 2i 2 (α + 1)α = (2 -α -(α + 1) + α(α -1) + (α + 1)α) p-1 i=0 g 2i + (α + (α + 1) -α(α -1) -(α + 1)α) p-1 i=0 g 2i+1 = (1 -2α + 2α 2 ) p-1 i=0 g 2i + (1 + 2α -2α 2 ) p-1 i=0 g 2i+1 ,
so that this quantity can depend on α which is not the case for a conservative approach, and thus we can not find gc satisfying (4) which yields to the same scheme.

Application to the Vlasov-Poisson model

As an application of the constant advection case, we are concerned with the Vlasov-Poisson model, the unknown of which f = f (t, x, v) is the electronic distribution function. It depends on the spatial variable x ∈ [0, L] where L > 0 is the size of the domain, the velocity direction v ∈ IR and the time t ≥ 0. The time evolution of this distribution function is given by the following phase space transport equation, the Vlasov equation

∂f ∂t + v∂ x f + E(t, x)∂ v f = 0, ( 11
)
with the initial condition

f (0, x, v) = f 0 (x, v).
The electric field E(t, x) is given by the coupling with the distribution function f through the Poisson equation

∂ x E(t, x) = ρ(t, x) -ρ i , L 0 E(t, x)dx = 0, (12)
where the electronic charge density ρ is given by ρ(t, x) = I R f (t, x, v)dv and ρ i denotes the ion density. In this work, we restrict ourselves to a uniform background of ions which leads to ρ i = 1 after a suitable choice of dimensionless parameters.

In view of finite volumes formulation, it will be convenient to re-write the Vlasov equation into a conservative form

∂f ∂t + ∂ x (vf ) + ∂ v (E(t, x)f ) = 0. ( 13
)
The Vlasov-Poisson model preserves some physical quantities with time which will be analysed and compared for the different numerical methods. First of all, the Vlasov-Poisson equation preserves the L p norms for p ≥ 0

d dt f (t) L p = 0. ( 14
)
The momentum d dt

L 0 I R vf (t, x, v)dxdv = 0,
and the total energy are also constant in time

d dt E(t) = d dt E k (t)+ d dt E e (t) = d dt L 0 I R f (t, x, v) |v| 2 2 dxdv+ 1 2 d dt L 0 I R |E(t, x)| 2 dxdy, (15)
where E e and E k denote the electric and kinetic energy respectively. From a numerical point of view, the good conservation of these different quantities is an important feature for Vlasov simulations.

The general algorithm

In this subsection, we review the main steps of a semi-Lagrangian method in the case of directional splitting with constant advection, which is applied for the discretization of the Vlasov-Poisson model. Discretization of the distribution function. The unknown quantities are then f n k,ℓ which are approximations of f (t n , x k , v ℓ ). We suppose periodic boundary conditions so that we only have to compute at each time t n

f n k,ℓ , for k = 0, . . . , N x -1, ℓ = 0, . . . , N v -1.
Transport operator. Let us define for α ∈ R and N ∈ N * a transport operator T α,N : R N → R N . We will define several possibilities for the transport operator.

If (f 0 , . . . , f N -1 ) is a discretization of a function f , then T α,N (f 0 , . . . , f N -1 )
should be a discretization of the shifted function x → f (x + α). We also denote by T x α : R Nx → R Nx (resp. T v α : R Nv → R Nv ) a transport operator which shifts along the x (resp. v) direction.

Algorithm. The time splitting algorithm then reads

Step 0. Initialization :

f k,ℓ = f 0 (x k , v ℓ ), k = 0, . . . , N x -1, ℓ = 0, . . . , N v -1.
Step 1. Half time step shift along the x-axis:

For each ℓ = 0, . . . , N v -1, (f k,ℓ ) Nx-1 k=0 → T x α ((f k,ℓ ) Nx-1 k=0 ) with α = -v ℓ ∆t/2.
Step 2. Computation of the charge density and the electric field by integrating [START_REF] Grandgirard | A drift-kinetic semi-Lagrangian 4D code for ion turbulence simulation[END_REF].

Step 3. Shift along the v-axis:

For each k = 0, . . . , N x -1, (f k,ℓ ) Nv-1 ℓ=0 → T v α ((f k,ℓ ) Nv-1 ℓ=0 ) with α = -E k ∆t.
For each ℓ = 0, . . . , N v -1, (f k,ℓ ) Nx-1 k=0 → T x α ((f k,ℓ ) Nx-1 k=0 ) with α = -v ℓ ∆t/2.
Step 5. n → n + 1 and loop to Step 2.

Note that if we make no diagnostic of the distribution function, we can simplify

Step 4.a-c into

Step 4. Shift along the x-axis:

For each ℓ = 0, . . . , N v -1, (f k,ℓ ) Nx-1 k=0 → T x α ((f k,ℓ ) Nx-1 k=0
) with α = -v ℓ ∆t. In the sequel, we present numerical results for the Vlasov-Poisson equation for which several choices of transport operators T α,N are performed.

Numerical results: Bump-on-tail test case

The present numerical schemes are validated on academic test cases introduced in [START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF][START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF][START_REF] Nakamura | Cubic Interpolated Propagation Scheme for Solving the Hyper-Dimensional Vlasov-Poisson Equation in Phase Space[END_REF]. In the present work, numerical results obtained by the methods of section 2 are applied on the bump-on-tail instability test case for which the initial condition writes (see [START_REF] Shoucri | Nonlinear evolution of the bump-on-tail instability[END_REF][START_REF] Nakamura | Cubic Interpolated Propagation Scheme for Solving the Hyper-Dimensional Vlasov-Poisson Equation in Phase Space[END_REF])

f 0 (x, v) = f (v)[1 + 0.04 cos(kx)], x ∈ [0, L x ], v ∈ [-v max , v max ],
with k = 0.5, v max = 9, L x = 20π; moreover, we have

f (v) = n p exp(-v 2 /2) + n b exp - |v -u| 2 2v 2 t RR n°6856
whose parameters are

n p = 9 10(2π) 1/2 , n b = 2 10(2π) 1/2 , u = 4.5, v t = 0.5.
The numerical parameters are N x = 128, N v = 128, ∆t = 0.1.

For this test case, we are interested in the time evolution of the total energy E together with the electric energy E e given by [START_REF] Kershaw | The explicit inverse of two commonly occuring matrices[END_REF]. We also look after the L p , p = 1, 2 norms of f (see ( 14)) which are conserved with time. The same is true for the total energy whereas the electric energy has an oscillatoring behaviour for large times (see [START_REF] Nakamura | Cubic Interpolated Propagation Scheme for Solving the Hyper-Dimensional Vlasov-Poisson Equation in Phase Space[END_REF]). Three vortices are then created in the phase space which are moving along the velocity v = v t (BGK equilibrium) and a lack of accuracy leads to a vorticies merging which has for consequence a loss of the oscillatoring behaviour of the electric energy.

Our goal is to compare the different methods we talked about using a Strang splitting (see the previous algorithm): the conservative methods PFC of [START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF], PFC2 (using the filter [START_REF] Nakamura | Cubic Interpolated Propagation Scheme for Solving the Hyper-Dimensional Vlasov-Poisson Equation in Phase Space[END_REF], PSM (using the filter ( 21)), PSM ZER of [START_REF] Zerroukat | A monotonic and positivedefinite filter for a Semi-Lagrangian Inherently Conserving and Efficient (SLICE) scheme[END_REF] and the advective methods LAG3 and SPL. Note that the Forward Update method introduced in section 2.2 is not shown in the present case; indeed the displacements are constant so that only uniform mesh are generated by the characteristics and the method is completely equivalent to the SPL or the PSM methods. Fig. 3 shows the time evolution of the electric energy for the different methods. The main features (see [START_REF] Shoucri | Nonlinear evolution of the bump-on-tail instability[END_REF][START_REF] Nakamura | Cubic Interpolated Propagation Scheme for Solving the Hyper-Dimensional Vlasov-Poisson Equation in Phase Space[END_REF]) of the expected behaviour is respected by all the methods: the electric energy presents a maximum at t ≈ 20 ω -1 p and then an slowing oscillatoring behaviour on which is superimposed the oscillation of the system at ω p . Nevertheless, two classes can be distinguished: splines and Lagrange based methods. For the methods based on Lagrange interpolation (LAG3, PFC), the oscillations of the electric energy due to the particles trapping are damped and the amplitude is decreasing for large time. It is not the case for the splines based methods (SPL, PSM) which keep the slow oscillatoring behaviour and a constant amplitude up to the end of the simulation. This can be explained by the fact that fine structures are developed in the vortices; they are quickly eliminated by the methods based on Lagrange interpolation whereas splines methods follow these thin details of the phase space solution for longer times. Hence, methods based on Lagrange interpolation seem to be not sufficiently accurate to describe this kind of phenomenon.

In Fig. 4, the L p , p = 1, 2 norms are plotted with respect to time for LAG3, PFC, PFC2, SPL and PSM. As expected, the conservative methods (PFC, PFC2, PSM) which can be coupled with a positive filter ensure the positivity of the solution so that the L 1 norm is preserved. It is not the case for the advective semi-Lagrangian methods LAG3 and SPL. In particular, LAG3 gives rise to bad results and the filters of the PFC methods has to act often to ensure positivity of the reconstruction. The SPL method also generates negative values, but in a reasonable proportion so that PSM filters do not intervene often. The two classes of methods are emphasized by the time evolution of the L 2 norm. Methods based on Lagrange interpolation (which are nearly superimposed) present a more diffusive behaviour compared to spline based methods. The smoothing effects of the methods based on Lagrange interpolation lead to a bad conservation of the L 2 norm since they do not follow the filamentation of the distribution function for long times. In Fig. 5, we distinguish the splines based methods (left figure) from the Lagrange ones (right figures). We add on the left figure the numerical results obtain by the PSM ZER method coupled with a monotone filter introduced in [START_REF] Zerroukat | A monotonic and positivedefinite filter for a Semi-Lagrangian Inherently Conserving and Efficient (SLICE) scheme[END_REF]. This filter is not very appropriate in our context and the optimal filter reconstructs a solution which is closer to the basic cubic splines solution.

Finally, in Fig. 6, the time evolution of the total energy is plotted for the different methods. This quantity is quite difficult to preserve at the discrete level (see [START_REF] Umeda | Comparison of numerical interpolation schemes for one-dimensional electrostatic Vlasov code[END_REF]4]). In particular, it is very difficult to ensure both positivity of the solution and conservation of the total energy for nonlinear tests. Numerical results are presented in Fig. 6; we first notice that no method preserves exactly this quantity. Moreover, we can observe the influence of the slope limiters on the results: indeed, the difference between PFC and Lagrange method is roughly the same as for the L 1 time evolution (see Fig. 4): PFC ensures positivity to LAG3 but it is to the detriment of the total energy conservation. This figure also emphasizes the good behaviour of the optimal filter applied to PSM. The conservation of the total energy is lower than 0.1%, which is similar to the SPL method.

As a conclusion, the PSM filter improves the numerical results of the original PSM method for every conserved quantity. The loss of total energy conservation is about a 0.1%, which remains very reasonable. Let us remark that a time discretization refinement has a great influence on the conservation of the total energy.

Non-constant case: the guiding-center model

In this work, we also deal with another type of Vlasov equation for which the advection term is not constant. The so-called guiding-center model enters in this category (see [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF]). This model, which has been derived to describe highly magnetized plasma in the transverse plane of a tokamak, considers the evolution of the particles density ρ(t, x, y) Figure 5: Time evolution of the L 2 norm for spline based methods PSM ZER (with monotone filter of [START_REF] Zerroukat | A monotonic and positivedefinite filter for a Semi-Lagrangian Inherently Conserving and Efficient (SLICE) scheme[END_REF]), PSM (optimal filter) and SPL (left) and methods based on Lagrange interpolation PFC, PFC2 (PFC with optimal filter) and LAG3 (right). N x = N v = 128, ∆t = 0.1 for the Bump-on-tail test.

∂ t ρ + E ⊥ • ∇ρ = 0, (1)
We denote by E ⊥ = (E y , -E x ). The specificity of (1) lies on the fact that one-dimensional splitting cannot (in principle) be applied (see [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF][START_REF] Huot | Instability of the time splitting scheme for the one-dimensional and relativistic Vlasov-Maxwell system[END_REF]) since the advection term E ⊥ depends on (x, y). Consequently, this model contains additional difficulties compared to the Vlasov-Poisson model and seems to be a good candidate to test numerical methods.

To that purpose, we briefly recall the conservation properties of (1) which should be preserved in the best manner by the numerical schemes. The guiding center model (1) preserves the total mass, the L 2 norm of the density (enstrophy) and the L 2 norm of the electric field (energy)

d dt ρ(t, x, y)dxdy = d dt ρ(t) L 2 = d dt E(t) L 2 = 0. (3)

The general algorithm

In this subsection, we review the main steps of a semi-Lagrangian method in the case of directional splitting which is applied for the discretization of the guiding-center-Poisson model. Grid notations. Let N x , N y ∈ N * , y max > 0, a time step ∆t > 0.

We define then classically as notations ∆x = L x /N x , ∆y = L y /N y x k = kL x /N x , y ℓ = ℓL y /N y for k = 0, . . . , N x , ℓ = 0, . . . , N y and t n = n∆t, n ∈ N. Discretization of the distribution function. The unknown quantities are then ρ n k,ℓ which are approximations of ρ(t n , x k , y ℓ ). We suppose periodic boundary conditions so that we only have to compute at each time t n ρ n k,ℓ , for k = 0, . . . , N x -1, ℓ = 0, . . . , N y -1.

Transport operator. Let us define for (α k ) ∈ R Nx+1 a transport operator T α : R Nx → R Nx . For the conservative approaches we detailed in section 2.1, this operator writes T α (ρ 0 , ρ1 , ..., ρNx-1 ) = 1 ∆x

x k+1/2 -α k+1/2 x k-1/2 -α k-1/2 ρ(x)dx k=0,..,N -1
.

The sequence α is determined following one of the algorithms detailed in sections 2.1.1 and 2.2.1.

Algorithm

Step 0. Initialization : ρ k,ℓ = ρ 0 (x k , y ℓ ), k = 0, . . . , N x -1, ℓ = 0, . . . , N y -1.

Step 1. Compute of the electric field (E 0

x , E 0 y ) by integrating (2).

Step 2. Compute ρ 1 k,ℓ using ρ 0 : Step 2.a. Half time step shift along the x-axis:

Compute the x-displacement for each ℓ α k = ∆t/4E 0 y (x k -α k , y ℓ ) For each ℓ = 0, . . . , N y -1, (ρ k,ℓ ) Nx-1 k=0 → T x α ((ρ k,ℓ ) Nx-1 k=0 ).
Step 2.b. Shift along the y-axis:

Compute the y-displacement for each k α ℓ = -∆t/2E 0 x (x k , y ℓ -α ℓ ) For each k = 0, . . . , N x -1, (ρ k,ℓ ) Ny-1 ℓ=0 → T y α ((ρ k,ℓ )
Ny-1 ℓ=0 ).

Step 2.c. Half time step shift along the x-axis: Compute the x-displacement for each k α k = ∆t/4E 0 y (x k -α k , y ℓ ) For each ℓ = 0, . . . , N y -1, (ρ k,ℓ ) Nx-1 k=0 → T x α ((ρ k,ℓ ) Nx-1 k=0 ).

Step 3. Compute the electric field (E 1

x , E 1 y ) by integrating (2).

k = ∆t/2E n y (x k -α k , y ℓ ) For each ℓ = 0, . . . , N y -1, (ρ k,ℓ ) Nx-1 k=0 → T x α ((ρ k,ℓ ) Nx-1 k=0 ). Step 4.b. Shift along the y-axis: Compute the y-displacement for each k α ℓ = -∆tE n x (x k , y ℓ -α ℓ ) For each k = 0, . . . , N x -1, (ρ k,ℓ ) Ny-1 ℓ=0 → T y α ((ρ k,ℓ )
Ny-1 ℓ=0 ). Different methods will be compared: PSM without filter, the forward Update Method (FUM), (detailed in subsection 2.2) and the traditional semi-Lagrangian method with cubic splines interpolation developed in [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF]. For this last approach, we look after the validity of the splitting. Indeed, as we already discussed in section 2, this method is not always conservative when the splitting procedure is performed (the method will be refered as SPL1D). However, if a full two-dimensional interpolation is performed, the method becomes conservative (the method will be refered as SPL2D). Hence, this point will be discussed.

Numerical results: Kelvin-Helmoltz instability test case

We consider the Kelvin-Helmoltz instability in the periodic-periodic case for which the growth rate of the instability can be computed a priori. This is of great importance to check at the qualitative point of view the accuracy of the code.

Following the computations of [START_REF] Shoucri | A two-level implicit scheme for the numerical solution of the linearized vorticity equation[END_REF], the linearization (1)-( 2) leads to the so-called stability Rayleigh equation. Considering as initial condition a periodic perturbation of the equilibrium solution to (1)-( 2), it is possible to start a Kelvin-Helmoltz instability. The difference between the Dirichlet-periodic case (which has been solved in [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF]) occurs in the neutrally stable solution which is INRIA equal to 1 in our case (instead of sin(y/2) in the Dirichlet-periodic case). Then, we can deduce the initial condition for (1)-(2), ρ(x, y, t = 0) = sin(y) + ε cos(kx),

where k = 2π/L x is the wave number associated to the length L x of the domain in the x-direction. The size of the domain in the y-direction is L y = 2π. Shoucri's analysis predicts an instability when k is chosen lower to 1. Otherwise, the initial perturbation remains unchanged, neither damped (since (1)-( 2) is only a fluid model, not a kinetic model), neither increased.

Various approaches can be employed to determine the instability growth rate of the chosen mode k. A finite difference numerical scheme has been applied to approximate the stability Rayleigh equation which leads to a eigenvalue problem. The results obtained by this way are very closed to those obtained by numerically solve the linearized problem as performed in [START_REF] Shoucri | A two-level implicit scheme for the numerical solution of the linearized vorticity equation[END_REF].

The numerical parameters are chosen as follows: k = 0.5, N x = N y = 128, and ∆t = 0.1.

Let us recall that periodic conditions are considered here; even if the present test bears similarities with the Dirichlet-periodic test presented in [START_REF] Sonnendrücker | Ghizzo The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF], the dynamics of the unknown is quite different in the present periodic-periodic context. For this test case, we are interested in the time evolution of the conserved quantities (3). We also look carefully at the conservation of the total mass in order to verify the difference between conservative and non-conservative methods. We focus our comparisons on splines based methods which appear to be more competitive for strongly nonlinear problems: SPL (SPL1D refers to nonconservative splitting and SPL2D refers to the full two-dimensional advection without splitting), FUM (Forward Update Method presented in subsection 2.2), and PSM (with splitting). As a diagnostics, it is also interesting to look after the 2D unknown to realize the fine structures developed along the simulation.

In Fig. 7, the time history of the total mass and the L 2 norm is plotted for the different methods. First, as discussed in [START_REF] Huot | Instability of the time splitting scheme for the one-dimensional and relativistic Vlasov-Maxwell system[END_REF][START_REF] Nakamura | Exactly conservative semi-Lagrangian scheme for multi-dimensional hyperbolic equations with directional splitting technique[END_REF][START_REF] Tanaka | Constructing exactly conservative scheme in a non-conservative form[END_REF], SPL1D does not preserve exactly the total mass whereas other methods do. This is expected since this approach solves the non-conservative form of the equation which is not appropriate with the splitting procedure. Then, we observe that the conservative methods present very similar behaviour compared to the method of reference SPL2D: they are conservative and the decay of the L 2 norm occurs at t ≈ 30 ω -1 p . This decay corresponds to the saturation of the instability. Very fine structures are created which can not be captured by the numerical schemes since their size becomes smaller than the grid size.

In Fig. 8, the logarithm of the first Fourier mode of the electric field E x is plotted as a function of time. The linear theory predicts an exponential growing, the rate of which can be computed a priori by solving an eigenvalue problem. This can be performed and the results can be compared to the numerical results. The numerical growth rate corresponds to the slope of the straight line which approximates the logarithm of the first Fourier mode of E x in the linear phase (between t ≈ 5 ω -1 p and t ≈ 10 ω -1 p ). Considering different values of the wave number k, it is possible to plot the quantity ω/k (where ω is the growth rate of the first Fourier mode of E x ) as a function of (k s -k) where k s = 1 in our case (k s = √ 3/2 in the Dirichlet-periodic case). This is performed in Fig. 9 (right); we can observe the very good agreement between the analytical and the numerical values. This kind of validation is of great importance since a quantitative comparison can be performed, at least on the linear phase. On Fig. 9 (left), the L 2 norm of the electric field is plotted as a function of time. This quantity is preserved with time by the continuous model. The conservative and splitting procedure based methods present a very good conservation of the energy whereas it is not the case of the non-conservative method SPL1D. The method SPL2D does not preserve very well the energy compared to FUM or PSM for example.

Finally, on Fig. 10 and Fig. 11 we plot the distribution function for the different methods at time t = 30 and t = 60ω -1 p . These results confirm the previous observations: first, the PFC method is diffusive (the thin structures are smoothed) and the SPL1D scheme leads to a bad behaviour since the main structures are not respected. In contrast, the FUM and PSM present a good behaviour, very similar to SPL2D.

As a conclusion, the conservative methods present a very good behaviour on this strongly nonlinear and large time case. The splitting procedure also enables to save memory since one-dimensional structures are often used (instead of twodimensional structures for the computation of the cubic spline coefficients in SPL2D for example). On the other side, we remarked that SPL2D seems to be able to support larger time steps, but we think that the implementation of high order numerical scheme in time could stabilize PSM or FUM when large time steps are used. This extension will be studied in a future work.

Conclusion

In this work, new conservative methods have been introduced and then compared to existing methods for equations occuring in plasma physics. Several properties make them very competitive. On the one side, their inherent conservation property enables the use of splitting procedure, which makes easier the implementation of multi-dimensional problems. On the other side, slope limiters can be introduced to ensure the positivity of the unknown. When they are compared to existing semi-Lagrangian methods, we first observed that for the guiding-center problem, as expected, the advective approaches lead to inaccurate results when splitting procedure is applied. This is not the case for conservative methods. Moreover, they are at least as accurate as the reference methods (SPL2D). For simpler cases like the Vlasov-Poisson model in which the advection term is constant, we proved that the advective methods and their conservative counterparts are equivalent. Obviously, this does not remain true in the non-constant advection case (like for the guiding-center model) which often occurs in particular in gyrokinetic models. The extension of the PSM method to such multi-dimensional system is currently investigated.

Moreover, the splitting procedure makes easier the use of high order numerical scheme in time for backward and forward approaches. The use of high order time splittings (see [START_REF] Yoshida | Construction of higher order symplectic integrators[END_REF]) is also under investigations.

Figure 1 :

 1 Figure 1: Comparison of the semi-Lagrangian LAG3, PFC and PFC2 reconsructions.

  1 for two examples of positive reconstructions and the corresponding polynomial Lagrange 3 (LAG3) reconstruction).

Figure 2 :

 2 Figure 2: Comparison of the semi-Lagrangian cubic splines SPL, PSM ZER of [26] and PSM reconstructions.

INRIA

  Step 4.a Half time step shift along the x-axis:For each ℓ = 0, . . . , N v -1, (f k,ℓ ) Nx-1 k=0 → T x α ((f k,ℓ ) Nx-1 k=0 ) with α = -v ℓ ∆t/2.Step 4.b We have f n k,ℓ = f k,ℓ , for k = 0, . . . , N x -1, ℓ = 0, . . . , N v -1. Step 4.c Half time step shift along the x-axis:

INRIA

Figure 3 :

 3 Figure 3: Time evolution of the electric energy for the different methods: semi-Lagrangian methods LAG3 and SPL (left) and conservative methods PFC and PSM (right). N x = N v = 128, ∆t = 0.1 for the Bump-on-tail test.

Figure 4 :

 4 Figure 4: Time evolution of the L 1 norm (left) and L 2 norm (right) for LAG3, PFC, PFC2, SPL and PSM. N x = N v = 128, ∆t = 0.1 for the Bump-on-tail test.

Figure 6 :

 6 Figure 6: Time evolution of the total energy for the different methods. N x = N v = 128, ∆t = 0.1 for the Bump-on-tail test.

Step 4 .

 4 Compute ρ n+1 k,ℓ using ρ n-1 , ρ n :Step 4.a. Half time step shift along the x-axis: Compute the x-displacement for each ℓ α

  Step 4.c. Half time step shift along the x-axis: Compute the x-displacement for each ℓ α k = ∆t/2E n y (x k -α k , y ℓ ) For each ℓ = 0, . . . , N y -1, (ρ k,ℓ ) Nx-1 k=0 → T x α ((ρ k,ℓ ) Nx-1 k=0 ). Step 5. Compute the the electric field (E n+1 x , E n+1 y ) by integrating (2). Step 6. n → n + 1 and loop to Step 4.

Figure 7 :

 7 Figure 7: Time evolution of the total mass and the enstrophy for SPL (with splitting (SPL1D) and without splitting (SPL2D)), FUM and PSM. N x = N y = 128, ∆t = 0.1 for the Kelvin-Helmoltz instability test.

INRIA

Figure 8 :

 8 Figure 8: Time evolution of the logarithm of the first Fourier mode for SPL (with splitting (SPL1D) and without splitting (SPL2D)), FUM and PSM. A zoom has been applied on the right figure. N x = N y = 128, ∆t = 0.1 for the Kelvin-Helmoltz instability test.

Figure 9 :

 9 Figure 9: Left figure: Time evolution of the energy for SPL (with splitting (SPL1D) and without splitting (SPL2D)), FUM and PSM. Right figure: normalized growth rate ω/k as a function of 1 -k. N x = N y = 128, ∆t = 0.1 for the Kelvin-Helmoltz instability test.

Figure 10 :

 10 Figure 10: Distribution function for the Kelvin-Helmoltz instability test at time t = 30ω -1 p . Respectively for top-left to bottom-right: PSM, FUM, SPL2D, PFC, SPL1D with N x = N y = 128, ∆t = 0.1 and SPL2D with N x = N y = 512, ∆t = 0.01.

INRIA

Figure 11 :

 11 Figure 11: Distribution function for the Kelvin-Helmoltz instability test at time t = 60ω -1 p . Respectively for top-left to bottom-right: PSM, FUM, SPL2D, PFC, SPL1D with N x = N y = 128, ∆t = 0.1 and SPL2D with N x = N y = 512, ∆t = 0.01.

RR n°6856

Crouseilles & Mehrenberger & Sonnendrücker