Nicolas Crouseilles

Mohammed Lemou

S V Raghurama Rao

Ankit Ruhi

M Sekhar

Asymptotic Preserving scheme for a kinetic model describing incompressible fluids

The kinetic theory of fluid turbulence modeling developed by Degond and Lemou in [6] is considered for further study, analysis and simulation. Starting with the Boltzmann like equation representation for turbulence modeling, a relaxation type collision term is introduced for isotropic turbulence. In order to describe some important turbulence phenomenology, the relaxation time incorporates a dependency on the turbulent microscopic energy and this makes difficult the construction of efficient numerical methods. To investigate this problem, we focus here on a multi-dimensional prototype model and first propose an appropriate change of frame that makes the numerical study simpler. Then, a numerical strategy to tackle the stiff relaxation source term is introduced in the spirit of Asymptotic Preserving Schemes. Numerical tests are performed in a one-dimensional framework on the basis of the developed strategy to confirm its efficiency.

Introduction

Modeling and simulation of turbulence in fluid flows is a challenging task. Though the research in this topic has been pursued for nearly a century, the ideal model is still elusive. We refer to the books of Pope [START_REF] Pope | Turbulent Flows[END_REF] and Wilcox [START_REF] Wilcox | Turbulence Modeling for CFD[END_REF] for details of various turbulence models. The concepts from Kinetic Theory of Gases have been of significant use in some of the models for turbulence. Extending this approach further, Degond and Lemou [START_REF] Degond | Turbulence Models for Incompressible Fluids Derived from Kinetic Theory[END_REF] developed a model fully based on Boltzmann type equation and its expanded analogy for a description of turbulence. Chen et al. [START_REF] Chen | Extended Boltzmann Kinetic Equation For Turbulent Flows[END_REF] also followed a similar approach in developing a different model. This work is a first step, following the approach of [START_REF] Degond | Turbulence Models for Incompressible Fluids Derived from Kinetic Theory[END_REF] for study, analysis and simulation of turbulence from kinetic theory.

In [START_REF] Degond | Turbulence Models for Incompressible Fluids Derived from Kinetic Theory[END_REF] it is proposed to describe a turbulent incompressible fluid flow through a probability distribution function f (t, x, v) of fluid elements (or structures) depending not only on time and on the position x of these elements but also on their velocity. In this description, the fluid elements are assumed to interact in order to bring the fluid to an isotropic distribution in velocity (isotropic turbulence) and this is described by a relaxation collision operator. On the other hand, the pressure in the fluid acts in order to maintain the probability nature of the distribution function in such a way that the incompressibility (divergence-free) condition is satisfied. In the relaxation term (collision kernel), the relaxation time can be tuned in a such way that the model

∂f ∂t + v • ∇ x f -∇ x P • ∇ v f = 1 η Q u f (f ) + ηQ diss (f ), (1)
where η > 0 is a constant parameter. We shall now define the different terms in this equation. First the operator Q diss is a dissipative operator which allows to take into account the dissipation of turbulence in a fluid, that is its tendency to become laminar. We shall omit it here in our preliminary numerical study and refer to [START_REF] Degond | Turbulence Models for Incompressible Fluids Derived from Kinetic Theory[END_REF] for a detailed description and different possible choices of such operators. Secondly, the operator Q u f is the relaxation collision operator that brings the distribution function to an isotropic state and is given by

Q u f (f )(v) = 1 τ |v-u f | 2 2 [Πf (v) -f (v)] ,
where u f is defined from the following (nonlinear) equation

R d (v -u f ) τ |v-u f | 2 2 f (v)dv = 0, (2)
and the relaxation time τ is a function of the microscopic energy

|v-u f | 2 2
which will be made precise later on. The operator Π is the projection on to the space of isotropic functions around

u f Πf (v) = 1 |S d-1 | ω∈S d-1 f (u f + |v -u f |ω)dω,
where S d-1 is the sphere in R d (|S d-1 | denotes its measure). The distribution function f is constrained to satisfy

R d f (v)dv = 1, (3)
and this allows to see the pressure P in (1) as a Lagrange multiplier of this constraint. Note that this kinetic model was also used in [START_REF] Degond | Turbulence Models for Incompressible Fluids Derived from Kinetic Theory[END_REF] to derive a k-ε type turbulence model or Modified k-ε turbulence model (M-k-ε).

Defining the mean velocity ū = R d vf (v)dv, we remark that in general, we have u f = ū. However when the relaxation time τ is a (nonzero) constant function, then these two velocities coincide (see [START_REF] Bouchut | The BGK model for small Prandtl number in the Navier-Stokes approximation[END_REF]). Now we emphasize that the collision operator Q u f (f ) satisfies the important properties of conservation and leads to an entropy function (see [START_REF] Degond | Turbulence Models for Incompressible Fluids Derived from Kinetic Theory[END_REF] for details). In particular it preserves the mass, momentum and energy (see [START_REF] Degond | Turbulence Models for Incompressible Fluids Derived from Kinetic Theory[END_REF])

R d   1 v |v| 2   Q u f (f )dv = 0,
and it is worthwhile to see that the constraint (2) is necessary to ensure the momentum conservation property.

The goal of this work is to provide numerical simulations to efficiently solve the model ( 1)-( 2)-(3) satisfied by f, u f and P. There are two main difficulties in constructing such efficient numerical schemes for this model. First, the determination of the velocity u f from formula [START_REF] Bouchut | The BGK model for small Prandtl number in the Navier-Stokes approximation[END_REF] requires the use of an iterative method like Newton-Raphson technique at each time step. This turns out to be numerically expensive and avoiding it would significantly accelerate the simulation, especially if one deals with multidimensional models. The second difficulty is with the stiffness of the collision operator. The usual way of solving stiffness problem is to use an implicit discretization for the collision operator, but this is not possible here since u f should be determined before. We note that advancing the mean velocity ū is possible but advancing in time the velocity u f is not possible since it is determined by [START_REF] Bouchut | The BGK model for small Prandtl number in the Navier-Stokes approximation[END_REF]. This means in particular that classical approaches based on Asymptotic Preserving schemes for a class of other stiff kinetic equations (see [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF][START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Dimarco | Asymptotic Preserving Implicit-Explicit Runge-Kutta Methods for Non-Linear Kinetic Equations[END_REF][START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF][START_REF] Klar | Numerical passage from radiative heat transfer to nonlinear diffusion models[END_REF][START_REF] Coron | Numerical passage from kinetic to fluid equations[END_REF][START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving compressible Navier-Stokes asymptotics[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF][START_REF] Pieraccini | Implicit-explicit schemes for BGK kinetic equations[END_REF]) do not work in this context.

The appellation Asymptotic Preserving has been introduced in [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF] for numerical schemes that are stable with respect to a small parameter (η in this work) and degenerate into consistent numerical schemes for the limit model when η goes to zero. This class of numerical schemes is particularly well adapted to our framework since we want to solve (1)-( 2)-( 3) for arbitrary small values of η, which may lead to a severe constraint on the time step (it has to be of order η for stability) when a standard explicit numerical scheme is employed.

The above mentioned difficulties in designing an AP scheme are overcome in this framework with an efficient strategy as described in the following steps. First we perform a suitable change of frame to make the full implicitization of the collision term possible. Usually, this implicitation is at the heart of designing AP schemes. Unfortunately this is not sufficient to get such AP schemes in our context, since an additional stiff transport term arises after the change of frame (as we shall see in the next section) and therefore this requires a specific treatment. In fact, this new stiff term originates with the constraint (2) in the new frame. The numerical scheme we propose is then based on two main steps. Firstly, we focus on the stiffness of the collision operator and make it fully implicit by following for example the strategies in [START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF]. As mentioned above, the obtained scheme still contains a stiff transport term. We then propose a suitable way to implicitize this new stiff term, by making a part of it implicit. We emphasize that the whole obtained implicit scheme is computationally explicit, which means that no additional computational step is needed to solve the implicit schemes compared to the resolution of an explicit one. In particular, this results in a first order numerical scheme which enjoys the Asymptotic Preserving property, and has the same computational cost as a standard explicit solver (it does not require any linear system to invert).

The rest of the paper is organized as follows. In section 2, we give a reformulation of the continuous kinetic model under study through a suitable change of frame. The asymptotic limit η → 0 is also studied. In section 3, a suitable time discretization is performed on the new formulation, and the spatial discretization is discussed. Finally, in section 6, some numerical results are presented in a one-dimensional setting to illustrate the behavior of the new scheme in various configurations.

A reformulation of the continuous kinetic model

In this section, the model is reformulated by using a suitable change of frame. In particular, the constraints (3) and (2) are reformulated in a more tractable way in view of the construction of efficient numerical schemes. The asymptotic limit η → 0 is then studied for this new formulation.

A change of frame

As underlined in the introduction above, the collision time τ allows to incorporate some important turbulence phenomenology (viscosity regime and turbulence regime following the Kolmogorov law). For instance it is shown that in the molecular viscosity regime (small turbulent energies ξ, with ξ = |v -u f | 2 /2), the relaxation time τ must satisfy τ (ξ) ∼ C/ξ whereas in the turbulent regime (large energies ξ), τ should satisfy τ (ξ) ∼ C/ √ ξ. Then a general power-law is considered for τ in [START_REF] Degond | Turbulence Models for Incompressible Fluids Derived from Kinetic Theory[END_REF],

τ (ξ) = ξ α with ξ = |v -u f | 2 /2, -1 2 ≥ α ≥ -1
, and extended k-ε models are derived from the corresponding kinetic equations. We will not detail this study here and refer to [START_REF] Degond | Turbulence Models for Incompressible Fluids Derived from Kinetic Theory[END_REF] for a thorough description of these aspects.

Here we will then consider τ as a general function of the particle kinetic energy

ξ = |v-u f | 2 /2, i.e., τ = τ |v-u f | 2 2
. One of the main difficulties at the numerical level is to determine u f which is the solution of the following nonlinear equation

R d v -u f τ ( |v-u f | 2 2 ) f dv = 0.
In the following we will use the notation u := u f . Hence, for a known function f , the computation of u needs typically a Newton type algorithm, which can be very expensive. Since we are working on a model where most of the parameters depend on (v -u), it makes sense to rewrite the equations after using the transformation v → (v -u) = v ′ (see [START_REF] Degond | On viscosity and termal conduction of fluids with multivalued internal energy[END_REF]). The greatest advantage of using this change of frame is that both τ and the projector Π can now be considered independent of time and space variables. In particular, the projector in the new frame will commute with the transport operator and the time derivative, a fact of great importance in the construction of efficient numerical schemes. Using this transformation, the constraints (3) and ( 2) can be rewritten as:

R d F dv = 1, (4)
R d v τ |v| 2 2 F dv = 0, (5)
where f (t, x, v) = F (t, x, v -u). Moreover, since R d vf dv = ū, we have

R d vF dv = ū -u. (6)
Let us now write the equation satisfied by F . The derivatives of F in terms of f can be written as

∂ t f = ∂ t F -∂ t u • ∂ v F, ∂ x f = ∂ x F -(∂ x u)∂ v F, ∂ v f = ∂ v F, with (∂ x u) i,j = ∂ xi u j .
Using the above relations we can rewrite (1) as follows

∂ t F + (v + u) • ∂ x F -E + (∂ x u) T v • ∂ v F = 1 η 1 τ |v| 2 2 (Π 0 F -F ) , (7)
where

E = ∂ t u + ∂ x P + (∂ x u) T u and Π 0 F = 1 |S d-1 | S d-1 F (|v|w)dw.
Hence, the unknowns (F, u, P, E) satisfy the model given by the equations

                   ∂ t F + (v + u) • ∂ x F -E + (∂ x u) T v • ∂ v F = 1 η 1 τ |v| 2 2 (Π 0 F -F ) , ∂ t u + (u • ∂ x )u = E -∂ x P, R d F dv = 1, R d v τ F dv = 0. ( 8
)
This model is equivalent to the original model ( 1)-( 3)-(2) satisfied by (f, u, P).

In the following, we rewrite the constraints R d F dv = 1 and R d v τ F dv = 0 in a more tractable way to construct numerical schemes.

Let us first consider the constraint R d F dv = 1. Integrating the first equation of ( 8) with respect to v leads to

∂ x • vF dv + u = 0, (9)
since -

[(∂ x u) T v] • ∂ v F dv = ∂ x • u.
Then, taking the divergence of the first moment of the first equation of ( 8) and summing the obtained equation to the divergence of the second equation of (8) leads to (using ( 9))

∂ x • v ⊗ (v + u)∂ x F dv -∂ x • v(E + (∂ x u) T v) • ∂ v F dv + ∂ x • [(u • ∂ x )u] = ∂ x • E -∆P. ( 10
)
The left hand side gives ∂ 2

x : (v + u) ⊗ (v + u)F dv + ∂ x • E -∂ x • [(∂ x • u)u]
which finally gives an equation for the pressure P

-∆P = ∂ 2 x : (v + u) ⊗ (v + u)F dv -∂ x • [(∂ x • u)u], (11)
where " : " denotes the contracted product of two tensors. We refer the reader to appendix 7.1 for the details of the computations.

Let us now consider the constraint R d v τ F dv = 0. Multiplying the first equation of (8) by v/τ and integrating with respect to v gives (using the constraint

R d v τ F dv = 0) v τ (v + u) • ∂ x F dv - v τ E • ∂ v F dv - v τ [(∂ x u) T v] • ∂ v F dv = - 1 η v τ 2 F dv.
This provides an explicit expression of

E v τ ⊗ ∂ v F dv E = 1 η v τ 2 F dv + ∂ x • 1 τ v ⊗ (v + u)F dv - v τ [(∂ x u) T v] • ∂ v F dv. ( 12
)
Let us remark that the

d × d matrix v τ ⊗ ∂ v F dv is invertible (see Appendix 7.2).
Finally, the system to solve is the following, satisfied by (F, u)

                     ∂ t F + (v + u) • ∂ x F -E + (∂ x u) T v • ∂ v F = 1 η 1 τ |v| 2 2 (Π 0 F -F ) , ∂ t u + (u • ∂ x )u = E -∂ x P, -∆P = ∂ 2 x : (v + u) ⊗ (v + u)F dv + ∂ x • [(∂ x • u)u], v τ ⊗ ∂ v F dv E = 1 η v τ 2 F dv + ∂ x • 1 τ v ⊗ vF dv - v τ [(∂ x u) T v] • ∂ v F dv. (13)
Remark 1. The case τ constant (equal to 1) leads to specific expression of E since the constraint (5) becomes R vF dv = 0 and as a consequence the stiff term in the last equation of ( 13) vanishes. Then, in this configuration, the equation for E reduces to

v ⊗ ∂ v F dv E = ∂ x • v ⊗ vF dv -v[(∂ x u) T v] • ∂ v F dv,
and the equation for F becomes

∂ t F + (v + u) • ∂ x F -(E + (∂ x u) T v) • ∂ v F = 1 η (Π 0 F -F ).
This equation has the same structure as the Vlasov-Poisson-BGK model in the hydrodynamic regime. Applying Π 0 to this equation enables to derive a non stiff macro equation. At the numerical level, IMEX methods proposed in [START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving compressible Navier-Stokes asymptotics[END_REF][START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF][START_REF] Dimarco | Asymptotic Preserving Implicit-Explicit Runge-Kutta Methods for Non-Linear Kinetic Equations[END_REF][START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF][START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF][START_REF] Pieraccini | Implicit-explicit schemes for BGK kinetic equations[END_REF] can be applied in a straightforward way.

Remark 2. The equation for F in (13) can be formulated in a conservative form

∂ t F + v • ∂ x F -E • ∂ v F -(∂ x u) T : ∂ v (vF ) + ∂ x • (uF ) = 1 η 1 τ |v| 2 2 (Π 0 F -F ) ,
from which the conservation of mass d dt f (t, x, v)dxdv = 0 is deduced easily.

Asymptotic models from the reformulation

The aim of this section is to derive the asymptotic model, i.e., the model obtained formally from (1) when η goes to zero, or equivalently, from [START_REF] Earl | Cosmic ray viscosity[END_REF]. It is of importance to understand the asymptotic limit at the continuous level when one wants to design Asymptotic Preserving schemes. To alleviate the notation, in the following we use brackets for the velocity integrals

• := R d • dv.
First of all, we rewrite (13) using the following notation

AF = (v + u) • ∂ x F -(∂ x u) T v • ∂ v F, (14)
so that the first equation of ( 13) reads

∂ t F + AF -E • ∂ v F = 1 ητ (Π 0 F -F ).
where E is given by

E = 1 η v τ ⊗ ∂ v F dv -1 v τ 2 F dv + v τ ⊗ ∂ v F dv -1 ∂ x • 1 τ v ⊗ vF dv - v τ [(∂ x u) T v] • ∂ v F dv .
This term makes an additional stiff term appear so that the operator

LF := M • ∂ v F + 1 τ (Π 0 F -F ) with M = v τ ⊗ ∂ v F dv -1 v τ 2
F dv needs to be studied. This is the purpose of the following proposition.

Proposition 2.1. The null space of L is the same as the null space of Q 0 , with Q 0 (F ) = 1 τ (Π 0 F -F ), i.e., the space of functions φ(|v|), where the function φ is arbitrary.

Proof. Considering the relation L(F )F := R d L(F )F dv = 0 gives 1 τ (Π 0 F -F )F = -M∂ v F F , (15)
with

M = v τ ⊗ ∂ v F dv -1 v τ 2 F dv . We now use the fact that 1 τ (Π 0 F -F )Π 0 F = 0. Indeed, since Π 0 is self-adjoint in L 2 and commute with isotropic function like 1/τ , we have, 1 τ (Π 0 F -F )Π 0 F = Π 0 1 τ (Π 0 F -F ) F = 1 τ Π 0 (Π 0 F -F )F = 0.
Then, (15) becomes

1 τ (Π 0 F -F )(I -Π 0 )F = -M • ∂ v F F = 0, so that 1 √ τ (Π 0 F -F ) 2 L 2 = 0.
The function τ being positive, we conclude that if L(F )F = 0, then F = Π 0 F , which means that F is isotropic, i.e., F depends on |v|.

We can now derive formally the asymptotic model from [START_REF] Earl | Cosmic ray viscosity[END_REF] satisfied by (G 0 , ū) which are the limits when η goes to zero of (F, u). This is presented in the following proposition.

Proposition 2.2. The asymptotic model obtained from (13) in the limit η goes to zero, satisfied by the limit (G 0 , ū) of (F, u), writes

             ∂ t G 0 + ū • ∂ x G 0 = 0, ∂ t ū + (ū • ∂ x )ū = E 0 -∂ x P 0 , -∆P 0 = 1 d ∆ R d |v| 2 G 0 dv + ∂ 2 x : (ū ⊗ ū), E 0 = - 1 d ∂ x R d |v| 2 G 0 dv. (16)
Proof. We first derive the asymptotic model from the first equation of [START_REF] Earl | Cosmic ray viscosity[END_REF]. To do that, we consider the Chapman-Enskog expansion of F :

F = G 0 + ηG 1 + O(η 2 )
, where G 0 = Π 0 F . Inserting this expansion in the first equation of ( 13) leads to

∂ t G 0 + AG 0 -E • ∂ v G 0 = - 1 τ G 1 + O(η). (17)
Applying Π 0 to this latter equation leads to the asymptotic equation for G 0

∂ t G 0 + Π 0 AG 0 -Π 0 (E • ∂ v G 0 ) = 0. ( 18
)
The term Π 0 (E

• ∂ v G 0 ) vanishes whereas the term Π 0 AG 0 becomes Π 0 AG 0 = Π 0 (v + u) • ∂ x G 0 -(∂ x u) T v • ∂ v G 0 = u • ∂ x G 0 -Π 0 ((∂ x u) T v • ∂ v G 0 ). Rewritting the isotropic function G 0 (v) := g(ξ) with ξ = |v| 2 /2, we have ∂ v G 0 = v∂ ξ g so that the last term becomes Π 0 ((∂ x u) T v • ∂ v G 0 ) = (∂ x u) T : Π 0 (v ⊗ v)∂ ξ g. Let us compute Π 0 (v ⊗ v) Π 0 (v i v j ) = 1 S d-1 S d-1 |v| 2 ω i ω j dω = 1 S d-1 |v| 2 δ i,j S d-1 ω 2 i dω = 1 dS d-1 |v| 2 δ i,j S d-1 d i=1 ω 2 i dω = 1 dS d-1 |v| 2 δ i,j S d-1 dω = 1 dS d-1 |v| 2 δ i,j S d-1 = 1 d |v| 2 δ i,j .
We then conclude (∂

x u) T Π 0 (v • ∂ v G 0 ) = 1 d |v| 2 ∂ x • u.
When η goes to zero, we have u converging towards ū since vF dv = ū -u with F which tends towards G 0 and vG 0 dv = 0; hence, the incompressibility condition ∂ x • ū = 0 makes this term null when η goes to zero. The asymptotic model of the first equation of ( 13) then writes

∂ t G 0 + ū • ∂ x G 0 = 0,
which is the first equation of the asymptotic model (16).

To get the limit of the second equation of ( 13), we need to derive the limit of E and ∂ x P when η goes to zero. The so-obtained model will be satisfied by the limit ū of u.

Let us first focus on the derivation of the limit of E. To do that, we apply (I -Π 0 ) to [START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF] and get (neglecting O(η) terms)

(I -Π 0 )AG 0 -E • ∂ v G 0 = - 1 τ G 1 .
We multiply by v and integrate with respect to v to get

v(I -Π 0 )AG 0 + E = 0. Indeed, we used (v/τ )G 1 = 0 since 0 = (v/τ )F = (v/τ )G 0 + η (v/τ )G 1 and (v/τ )G 0 = 0 (since G 0 is even)
, so that we deduce (v/τ )G 1 = 0. The expression of the limit E 0 of E when η goes to zero then becomes

E 0 = -v(I -Π 0 )AG 0 = -vAG 0 = -v (v + u) • ∂ x G 0 -(∂ x u) T v • ∂ v G 0 = -v(v + u) • ∂ x G 0 -v(∂ x u) T v • ∂ v G 0 ) = -∂ x • v ⊗ vG 0 = - 1 d ∂ x |v| 2 G 0 (19)
We then focus on the limit equation satisfied by the pressure P 0 , limit of P when η goes to zero. Considering the third equation of ( 13) when F is replaced by its limit G 0 and u by its limit ū, we have

-∆P 0 = ∂ 2 x : (v + ū) ⊗ (v + ū)G 0 + ∂ x • [(∂ x • ū)ū] = ∂ 2 x : (v + ū) ⊗ (v + ū)G 0 = ∂ 2 x : v ⊗ vG 0 + ∂ 2 x : (ū ⊗ ū) = 1 d ∆ |v| 2 G 0 + ∂ 2 x : (ū ⊗ ū). ( 20
)
The asymptotic of the second equation of (13) satisfied by ū which is the limit of u when η goes to zero, finally writes

∂ t ū + (ū • ∂ x )ū = E 0 -∂ x P 0 ,
where E 0 and P 0 satisfy [START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF] and [START_REF] Pieraccini | Implicit-explicit schemes for BGK kinetic equations[END_REF]. This corresponds to the second equation of ( 16), and with [START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF] and [START_REF] Pieraccini | Implicit-explicit schemes for BGK kinetic equations[END_REF] correspond to the third and fourth equation of ( 16). This concludes the proof.

Numerical schemes

In this section, we propose some numerical schemes for the kinetic models described previously. We divide our approach into the cases τ constant and τ non constant, the latter case being more complicated and needing different techniques. Our main goal is to provide semi-discretized (in time) numerical schemes which enjoy the Asymptotic Preserving (AP) property, that is: (i) they are stable with respect to the parameter η, (ii) they degenerate when η goes to zero towards numerical schemes which are consistent with the asymptotic model [START_REF] Klar | Numerical passage from radiative heat transfer to nonlinear diffusion models[END_REF]. Such schemes have been initiated in [START_REF] Jin | Efficient asymptotic-preserving (AP) schemes for some multiscale kinetic equations[END_REF] and have become, during the last decade, very popular to approximate stiff kinetic equations.

One important problem here to derive an AP scheme is to make the relaxation operator implicit since the determination of u requires the knowledge of f . This involves a nonlinear system to solve which would be computationally very costly and difficult to handle. As presented in the previous section, one issue is to consider the change of variable v ′ → v -u and the equation on F [START_REF] Degond | On viscosity and termal conduction of fluids with multivalued internal energy[END_REF]. However, an additional stiff term arises after this change of frame and the derivation of a numerical scheme where only the relaxation term 1 τ (Π 0 F -F ) is considered implicit does not provide the AP property of this scheme, as we shall see. We will also look at an exponential integrator following [START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF] but in this case also, this is not sufficient to get the AP property. Finally, we propose a new strategy which is a modification of the latter approaches, that ensures the numerical scheme to degenerate as η goes to zero towards a consistent numerical scheme of the asymptotic model [START_REF] Klar | Numerical passage from radiative heat transfer to nonlinear diffusion models[END_REF]. We emphasize that this numerical scheme has the same computational cost as a standard explicit solver.

In the following, we present two strategies from the literature (see [START_REF] Dimarco | Asymptotic Preserving Implicit-Explicit Runge-Kutta Methods for Non-Linear Kinetic Equations[END_REF][START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Pieraccini | Implicit-explicit schemes for BGK kinetic equations[END_REF]) and explain why they cannot be applied in our framework. We then propose our new scheme in two different versions. As we deal with a semi-discretization in time in this section, we introduce a time discretization t n = n∆t, n ∈ N, ∆t > 0 and we denote by F n an approximation of F at time t n .

Let us recall ( 13)

                       ∂ t F + AF -E • ∂ v F = 1 η 1 τ |v| 2 2 (Π 0 F -F ) , ∂ t u + (u • ∂ x )u = E -∂ x P, -∆P = ∂ 2 x : (v + u) ⊗ (v + u)F + ∂ x • [(∂ x • u)u], E = 1 η A(F ) -1 v τ 2 F + A(F ) -1 ∂ x • v τ ⊗ (v + u)F - v τ [(∂ x u) T v] • ∂ v F , (21)
where A is given by ( 14) and A(F ) is the d × d matrix defined by A(F ) = v τ ⊗ ∂ v F . We introduce the following notations for the time discretization of the term E

E n,n = 1 η A(F n ) -1 v τ 2 F n + E n 1 , (22)
and

E n,n+1 = 1 η A(F n ) -1 v τ 2 F n+1 + E n 1 , (23)
where E n 1 is given by

E n 1 = A(F n ) -1 ∂ x • v τ ⊗ (v + u n )F n - v τ [(∂ x u n ) T v] • ∂ v F n . ( 24
)
State of the art.

In this part, we look at two standard strategies to derive an AP numerical scheme for ( 21)

∂ t F + AF -E • ∂ v F = 1 τ η (Π 0 F -F ), ( 25
)
where A is given by ( 14).

As we shall see, the presence of the stiff term in the left hand side is the major obstacle in a direct application of standard strategies.

The first idea would be to consider implicitizing the relaxation term

F n+1 = F n -∆tAF n + ∆tE n,n • ∂ v F n + ∆t τ η (Π 0 F n+1 -F n+1 ).
To get an expression of Π 0 F n+1 , we apply Π 0 to this last equation

Π 0 F n+1 = Π 0 F n -∆tΠ 0 AF n + ∆tE n,n • Π 0 (∂ v F n ).
This macro equation still involves a stiff term in the term E n,n given by ( 22) which requires the time step ∆t to be lower than η. This is due to the fact that the other stiff term has not been considered implicit.

A second idea to design an AP scheme is to use an exponential integrator ( [START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Dimarco | Exponential Runge-Kutta methods for stiff kinetic equations[END_REF]). From (25), one has

∂ t (e t/(ητ ) F ) = -e t/(ητ ) AF + e t/(ητ ) η E • ∂ v F n + e t/(ητ ) τ η Π 0 F.
Integrating between t n and t n+1 leads to

F n+1 = e -∆t/(ητ ) F n - t n+1
t n e (t-t n+1 )/(ητ ) AF dt +

t n+1 t n e (t-t n+1 )/(ητ ) E • ∂ v F dt + t n+1 t n e (t-t n+1 )/(ητ ) τ η Π 0 F dt.
The key point of this strategy is the approximation of the integrals (see [START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF]). For the two first integrals, we choose the left rectangle quadrature whereas the last integral is computed exactly after considering Π 0 F implicit. This leads to

F n+1 = e -∆t/(ητ ) F n -∆te -∆t/(ητ ) AF n +∆te -∆t/(ητ ) E n,n • ∂ v F n + (1 -e -∆t/(ητ ) )Π 0 F n+1 .
We can remark that when η goes to zero, F n+1 goes to Π 0 F n+1 , which was not the case in the previous approach. One needs now to compute Π 0 F n+1 . To do that, we apply Π 0 to this last equation to get (recalling Π 0 commutes with isotropic functions)

Π 0 F n+1 = e -∆t/(ητ ) Π 0 F n -∆te -∆t/(ητ ) Π 0 (AF n ) +∆te -∆t/(ητ ) E n,n • Π 0 (∂ v F n ) + (1 -e -∆t/(ητ ) )Π 0 F n+1 . from which we get Π 0 F n+1 Π 0 F n+1 = Π 0 F n -∆tΠ 0 (AF n ) + ∆tE n,n • Π 0 (∂ v F n ), ( 26
)
where E n,n is given by [START_REF] Struchtrup | The BGK-model with velocity-dependent collision time[END_REF]. Here again, the macro equation ( 26) still contains a stiff term which cannot be considered implicit without requiring the inversion of a nonlocal operator. Let us remark that other choices of time integrals approximation do not lead to a semi-implicit Asymptotic Preserving scheme as well.

A slightly different strategy is then required since standard ones do not work. In the following, two schemes are proposed which enjoy the AP property. We point out that even if these new schemes involve an additional implicit term, namely E n,n+1 instead of E n,n , they have the same computational cost as an explicit numerical scheme.

New strategy: 1 (EXPO).

In the sequel, we present a modification of the previous exponential scheme. The numerical scheme for F writes

F n+1 = e -∆t/(ητ ) F n -∆te -∆t/(ητ ) AF n + ∆te -∆t/(ητ ) E n,n+1 • ∂ v F n + (1 -e -∆t/(ητ ) )Π 0 F n+1 . ( 27
)
The main problem now is to determine (v/τ 2 )F n+1 in E n,n+1 given by ( 23) without requiring the inversion of a linear system. Integrating this latter equation with respect to v, after multiplying by (v/τ 2 ), leads to

(v/τ 2 )F n+1 = (v/τ 2 )e -∆t/(ητ ) (F n -∆te -∆t/(ητ ) AF n ) + ∆t (v/τ 2 )e -∆t/(ητ ) E n,n+1 • ∂ v F n , (28) since (v/τ 2 )(1 -e -∆t/(ητ ) )Π 0 F n+1 = 0.
Let us focus on the last term to deduce (v/τ 2 )F n+1 in an explicit way

(v/τ 2 )e -∆t/(ητ ) E n,n+1 • ∂ v F n = 1 η (v/τ 2 )e -∆t/(ητ ) A(F n ) -1 v τ 2 F n+1 • ∂ v F n + (v/τ 2 )e -∆t/(ητ ) E n 1 • ∂ v F n = 1 η e -∆t/(ητ ) τ 2 (v ⊗ ∂ v F n )A(F n ) -1 v τ 2 F n+1 + (v/τ 2 )e -∆t/(ητ ) E n 1 • ∂ v F n , (29)
where E n 1 is given by (24). We can inject this term in (28) to get an expression of (v/τ 2 )F n+1

(v/τ 2 )F n+1 = I - ∆t η e -∆t/(ητ ) τ 2 v ⊗ ∂ v F n A(F n ) -1 -1 × (v/τ 2 )e -∆t/(ητ ) (F n -∆tAF n + ∆tE n 1 • ∂ v F n ) . = η ηI -∆t e -∆t/(ητ ) τ 2 v ⊗ ∂ v F n A(F n ) -1 -1 × (v/τ 2 )e -∆t/(ητ ) (F n -∆tAF n + ∆tE n 1 • ∂ v F n ) . (30)
Then, we have determined in an explicit way (v/τ 2 )F n+1 and we can remark that it is of order η. Then, the term E n,n+1 can be computed with an explicit complexity thanks to

E n,n+1 = 1 η A(F n ) -1 v τ 2 F n+1 + E n 1 , (31)
where E n 1 is given by (24).

The last thing to do is to express Π 0 F n+1 . To do that, we apply Π 0 to (27) (recalling that Π 0 commutes with isotropic functions)

Π 0 F n+1 = e -∆t/(ητ ) Π 0 F n -∆te -∆t/(ητ ) Π 0 (AF n ) +∆te -∆t/(ητ ) E n,n+1 • Π 0 (∂ v F n ) + (1 -e -∆t/(ητ ) )Π 0 F n+1 , from which we can compute Π 0 F n+1 Π 0 F n+1 = Π 0 F n -∆tΠ 0 (AF n ) + ∆tE n,n+1 • Π 0 (∂ v F n ). ( 32
)
This macro equation is the same as (26) except that the term (v/τ 2 )F n+1 in E n,n+1 is now implicit. Thanks to (30), this term is of order η making the last term of (32) non-stiff.

To advance the equation for u, the stiff term contained in E also has to be considered implicit. This can be done easily since v/τ 2 F n+1 has been computed previously and E n,n+1 is known. Hence, the numerical scheme for the equation on u writes

u n+1 = u n -∆t(u n • ∂ x )u n + ∆tE n,n+1 -∂ x P n , (33)
where E n,n+1 is computed from ( 23) and P n is computed from

-∆P n = ∂ 2 x : (v + u n ) ⊗ (v + u n )F n dv + ∂ x • [(∂ x • u n )u n ]. (34)
From (F n , u n ), the algorithm to compute (F n+1 , u n+1 can be summarized as follows

• compute P n with (34),

• compute (v/τ 2 )F n+1 with (30),

• compute E n,n+1 with (31),

• compute Π 0 F n+1 with (32),

• compute F n+1 with (27),

• compute u n+1 with (33).

Let us emphasize that the computational complexity of this algorithm is the same as the complexity of an explicit algorithm. The second step (computation of (v/τ 2 )F n+1 with (30)) only requires the inversion of a d × d matrix. The main properties of the so-obtained scheme (32)-( 30)-( 27)-( 33)-( 34)-(31) for [START_REF] Degond | On viscosity and termal conduction of fluids with multivalued internal energy[END_REF] are listed in the following proposition.

Proposition 3.1. The numerical scheme (32)-( 30)-( 33)-( 27)-( 34)-(31)

• is consistent and stable for (7);

• degenerates as η goes to zero towards a numerical scheme which is consistent with [START_REF] Klar | Numerical passage from radiative heat transfer to nonlinear diffusion models[END_REF].

Proof. We observe that formally, (27) degenerates when η goes to zero to F n+1 = Π 0 F n+1 . Thanks to (30), the term E is of order one since the stiffness has been stabilized by considering E n,n+1 for its approximation. Hence, applying Π 0 to (27) and considering the limit η goes to zero, the term E n,n+1 •Π 0 (∂ v F n ) goes to zero since F n converges towards Π 0 F n and Π 0 (∂ v Π 0 F n ) = 0. Finally (32) degenerates as η goes to zero to

Π 0 F n+1 = Π 0 F n -∆tΠ 0 (AΠ 0 F n ),
with AF n given by ( 14). The same computations as in the continuous case leads to

Π 0 F n+1 = Π 0 F n -∆tu n • ∂ x Π 0 F n ,
which is an explicit discretization of the asymptotic model for G 0 := Π 0 F (first equation of ( 16)).

Applying now (I -Π 0 ) to (27), and since

F n+1 = Π 0 F n+1 + O(η), one obtains 0 = -∆te -∆t/(ητ ) (I -Π 0 )A(Π 0 F n ) + ∆te -∆t/(ητ ) E n,n+1 • (I -Π 0 )∂ v (Π 0 F n ).
Multiplying by ve -∆t/(ητ ) /∆t and integrating with respect to v (following the computations of the continuous case) leads to

0 = -v(I -Π 0 )A(Π 0 F n ) + vE n,n+1 • ∂ v (Π 0 F n ) ,
or, after an integration by parts (using Π 0 F n = 1)

E n,n+1 = -v(I -Π 0 )A(Π 0 F n ) = - 1 d ∂ x |v| 2 Π 0 F n ,
which is the same expression as the limit equation satisfied by E (last equation of ( 16)). The equation for the pressure P does not present any difficulty since (34) degenerates into (u n converges towards ūn as η goes to zero)

-∆P n = ∂ 2 x : (v + ūn ) ⊗ (v + ūn )Π 0 F n dv + ∂ x • [(∂ x • ūn )ū n ].
Finally, the equation (33) for u becomes, with the previous notations

u n+1 = u n -∆t(u n • ∂ x )u n + ∆tE n,n+1 -∆t∂ x P n ,
which degenerates when η goes to zero into

ūn+1 = ūn -∆t(ū n • ∂ x )ū n - ∆t d ∂ x |v| 2 Π 0 F n -∆t∂ x P n ,
which is an explicit time discretization of the asymptotic equation on ū (second equation of ( 16)). Finally, we check that this numerical scheme enjoys the Asymptotic Preserving property.

New strategy: 2 (DIMP).

Here, a slight modification of the previous scheme is presented. Instead of using exponential integrators, we use here a simple implicit scheme, including the same treatment as before of the stiff transport term. The numerical scheme for F writes

F n+1 = F n -∆tAF n + ∆tE n,n+1 • ∂ v F n + ∆t ητ (Π 0 F n+1 -F n+1 ).
where A is given by ( 14) and E n,n+1 is given by ( 23), so that, with δ = 1/(1 + ∆t/(τ η)), one obtains

F n+1 = δ F n -∆tAF n + ∆tE n,n+1 • ∂ v F n + ∆t ητ Π 0 F n+1 . (35)
As previously, one needs to determine (v/τ 2 )F n+1 ; this is done as in the previous case by integrating this latter equation against (v/τ 2 )

(v/τ 2 )F n+1 = I - ∆t η δ τ 2 v ⊗ ∂ v F n A(F n ) -1 -1 × (δv/τ 2 )(F n -∆tAF n + ∆tE n 1 • ∂ v F n ) . ( 36
)
Then, the term E n,n+1 can be computed thanks to

E n,n+1 = 1 η A(F n ) -1 v τ 2 F n+1 + E n 1 , (37)
where E n 1 is given by ( 24). The last thing to do is to express Π 0 F n+1 . To to that, we apply Π 0 to (35) and get an expression of Π 0 F n+1

Π 0 F n+1 = Π 0 F n -∆tΠ 0 (AF n ) + ∆tE n,n+1 • Π 0 ∂ v F n . ( 38
)
The rest of the numerical scheme (computation of u n+1 , E n 1 and P n ) is the same as in the previous version. Hence, from (F n , u n ), the algorithm to compute (F n+1 , u n+1 can be written as follows

• compute P n with (34),

• compute (v/τ 2 )F n+1 with (36),

• compute E n,n+1 with (37),

• compute Π 0 F n+1 with (38),

• compute F n+1 with (35),

• compute u n+1 with (33).

As a conclusion, we derive the Asymptotic Preserving scheme (38)-( 36)-( 33)-( 35)-( 34)-( 37) for [START_REF] Degond | On viscosity and termal conduction of fluids with multivalued internal energy[END_REF], the properties of which are summarized in the following proposition.

Proposition 3.2. The numerical scheme (38)-( 36)-( 33)-( 35)-( 34)-(37)

• is consistent and stable for (7);

• degenerates as η goes to zero towards a numerical scheme which is consistent with [START_REF] Klar | Numerical passage from radiative heat transfer to nonlinear diffusion models[END_REF].

Proof. The proof is similar to the previous one and is left to the reader.

Spatial discretization

In this section, we focus on the one-dimensional case and we briefly detail the phase space discretization used to approximate the differential operators in (27) or in (35). Let us write the "EXPO" numerical scheme in the one-dimensional configuration (the DIMP version is similar)

F n+1 = e -∆t/(ητ ) F n -∆te -∆t/(ητ ) AF n + ∆te -∆t/(ητ ) E n,n+1 ∂ v F n + (1 -e -∆t/(ητ ) )Π 0 F n+1 , ( 39
) where AF n is given by

AF n = (v + u n )∂ x F n -∂ x u n v∂ v F n ,
and E n,n+1 is given by

E n,n+1 = 1 η (v/τ 2 )F n+1 (v/τ )∂ v F n + E n 1 ,
with E n 1 given by

E n 1 = 1 (v/τ )∂ v F n ∂ x (v 2 /τ )F n -∂ x u n (v 2 /τ )∂ v F n . ( 40
)
If moreover periodic boundary conditions are used, one can show that ū is a constant. Indeed, considering the first moment of (1) enables to write an equation for ū

∂ t ū + ∂ x v 2 f + ∂ x P,
and using the change of frame, we get

∂ t ū + ∂ x |v + u| 2 F + ∂ x P. ( 41
)
The incompressibility condition ∂ x ū = 0 leads to the equation satisfied by the pressure P (applying ∂ x to (41)) -∂ 2 x P = ∂ 2 x |v + u| 2 F . Periodic boundary conditions enable us to write -∂ x P = ∂ x |v + u| 2 F and we deduce from (41) ∂ t ū = 0. In addition to the incompressibility condition ∂ x ū = 0, we conclude that ū is constant. As a consequence, since vF = ū -u, u is then determined directly from this latter relation and the equation ( 33) is not needed in this simplified one-dimensional framework.

In the rest of this section, we describe the phase space discretization. We introduce a phase space uniform grid of size ∆x > 0 in the x-direction and ∆v > 0 in the v-direction: x i = x min +i∆x, i = 0, . . . , N x and v j = v min +j∆v, j = 0, . . . , N v , with N x , N v ∈ N and x min , v min ∈ R. In the sequel, we denote by F n i,j an approximation of F (t n , x i , v j ).

We consider upwind (forward/backward difference depending on wave-speed) schemes. Typically, to approximate [START_REF] Filbet | A class of asymptotic-preserving schemes for kinetic equations and related problems with stiff sources[END_REF]) and E n 1 given by (40), we use the following approximations

AF = (v + u)∂ x F -(∂ x u)v∂ v F (see
(AF n ) i,j ≈ (v j + u i ) + (D x -F n ) i,j + (v j + u i ) -(D x + F n ) i,j +(-v j (D x u) i ) + (D v -F n ) i,j + (-v j (D x u) i ) -(D v + F n ) i,j ,
and

(E n 1 ) i,j = (D x (v 2 /τ )F n ) i,j -(D x u n f ) i,j (v 2 /τ )(D v F n ) i,j (v/τ )(D v F n ) i,j ,
with the standard notation a ± = a±|a| 2 , where D x,v ± denotes a first order one-sided difference operator (for example (D

x + F n ) i,j = F n i+1,j -F n i,j ∆x
), D x denotes a centered finite difference operator

(D x F n ) i,j = F n i+1,j -F n i-1,j 2∆x
. All the integrals with respect to v are approximated by a rectangle quadrature.

More sophisticated spatial discretizations can be used, such as WENO schemes to achieve high orders, but our goal here is to validate the whole strategy.

An important point is to ensure the conservations of R F dv and R (v/τ )F dv at the discrete level. The previous spatial discretizations do not ensure these conservations. To overcome this problem, the scheme should be extended using some correction so that the above conservations are ensured irrespective of the discretization used. Let us mention that the above conservations can be ensured without an a posteriori correction by computing the discrete version of E and P by following the same strategy as in the continuous level. However, this is more restrictive and we choose the first approach which needs an a posteriori correction which we detail below.

To do that, we propose to a posteriori modify the obtained F n+1 i,j so that conservations are ensured. Say we use discretization B to the scheme to get F n+1 = BF n following the time and spatial discretizations described above. In general, even if F n satisfies the conservations j F n i,j ∆v = 1 and j (v j /τ j )F n i,j ∆v = 0, the obtained F n+1 does not preserve the properties which the corresponding equation in continuous level has. Then, we introduce an orthogonal projector P onto span{(1, (v/τ )) F n } in L 2 (dv/F n ), such that

P(F n+1 -F n ) = 0, which means R (F n+1 -F n ) 1 v/τ dv = 0.
We hence modify F n+1 as follows

F n+1 = F n+1 + P(F n -F n+1 ).
An example of such a projector can be given as

P(g) = g F n + v τ g v τ F n v 2 τ 2 F n
, from which we easily deduce P(g) = g and v τ P(g) = 0. At the discrete level, the integrals are replaced by the sums on the index j.

Introducing the defined projector P along with the scheme gives the benefit of using any higher order upwind discretization for space and velocity in the scheme. This way of correction does not depend on the scheme used, i.e., the scheme may be implicit or explicit.

Numerical results

This section is devoted to the validation of the numerical scheme presented above. Different schemes will be compared here in the one-dimensional configuration (d = 1). First, the scheme (32)-( 30)-(27) will be called "EXPO" and compared with explicit schemes such as Runge-Kutta schemes for [START_REF] Degond | On viscosity and termal conduction of fluids with multivalued internal energy[END_REF]. Let us remark that the "DIMP" scheme gives results which are very similar to those obtained by "EXPO" so that only results obtained by EXPO will be shown here. For these latter schemes, the time step ∆t will be chosen to satisfy ∆t ≤ η which is not the case for "EXPO" scheme.

The numerical tests will be of increasing difficulty: firstly, the homogeneous case with constant and non constant τ will be tested and lastly, the non homogeneous case with a constant and non constant τ will be also validated.

Homogeneous, constant τ case

To illustrate the effect of the collision operator, we first consider the simple case (τ = 1)

∂ t F = 1 η (Π 0 F -F ), v ∈ R. (42)
The following initial condition is chosen

F (t = 0, v)) = F 0 (v) such that vF dv = 0 and F dv = 1 F 0 (v) = (1 -α)M 1 (v) + αM 2 (v), with M i (v) = 1 √ 2π exp(-|v -u i | 2 /2),
where

u i , i = 1, 2 satisfies (1-α)u 1 +αu 2 = 0. The choice u 1 = -4, u 2 = 1 and α = -u 1 /(u 2 -u 1 )
is made in the following tests. The following numerical parameters are used:

N v = 256, ∆t = 0.1, v max = 16.
We are interested in entropies Θ(F ) , where Θ is a convex function. Indeed, for any convex function of F , Θ(F ) decreases with time in (42): multiplying (42) by Θ ′ (F ) and integrating with respect to v leads to

∂ t Θ(F ) = 1 η Θ ′ (F )(Π 0 (F ) -F ) = 1 η [Θ ′ (F ) -Θ ′ (Π 0 F )](Π 0 (F ) -F ) ≤ 0, since Θ ′ (Π 0 F ) is isotropic, we have Θ ′ (Π 0 F )](Π 0 (F ) -F ) = 0.
We are also interested in the time dependent quantity E(t) = F (t) -Π 0 F (t) L 2 for which we know an analytical solution. Indeed, since Π 0 F does not depend on time (applying Π 0 to (42) directly implies this), we get (F -

Π 0 F )(t, v) = e -t/(ητ ) (F -Π 0 F ), hence E(t) = e -t/(ητ ) E(0).
In Figure 1, the time evolution of E is plotted for different values of η and different schemes in time (explicit RK schemes and EXPO which refers to the scheme (32)-( 30)-( 27)). The curves "REF" corresponds to the analytical solution. We observe, in this simplified configuration, that the "EXPO" scheme captures the solution exactly whereas Runge-Kutta schemes produce an error of ∆t p where p is the order of the Runge-Kutta scheme used.

Figure 2 shows the time evolution of E for different values of η. Here also, the "EXPO" scheme captures the solution exactly, the time step being fixed as 0.1.

Figure 3 displays the time evolution of Θ(F ) for two choices of entropies Θ(F ) = F log F and Θ(F ) = F 4 , obtained using RK3 and EXPO schemes, η = 1 and η = 0.1. We observe that these two quantities decrease with time as expected, the rate being stronger as η is smaller. Here again, the results obtained with RK3 and EXPO are very close.

In Figure 4, we plot the velocity dependence of the distribution function for different times, for η = 1. We observe the relaxation the function towards a symmetric state (which corresponds to isotropic function in 1D), in v.

Homogeneous, non constant τ case

We now consider the homogeneous case with a non constant τ . Starting from (1) satisfied by f

∂ t f = 1 ητ |v-u| 2 2 (Π u f -f ),
with u = u(t) and together with the constraints f = 1 and (v/τ (|v -u| 2 /2)f = 0, the change of frame v → (v -u) leads to

∂ t F -∂ t u∂ v F = 1 ητ (Π 0 F -F ).
The first constraint F = 1 is automatically preserved by this model. The second constraint (v/τ )F = 0 should be reformulated. Taking (v/τ ) moment leads to The model to solve in this context is then

-∂ t u (v/τ )∂ v F = 1 η (v/τ 2 )(Π 0 F -F ) = - 1 η (v/τ 2 )F , which enables to determine ∂ t u ∂ t u = 1 η (v/τ 2 )F (v/τ )∂ v F .
∂ t F - 1 η (v/τ 2 )F (v/τ )∂ v F ∂ v F = 1 ητ (Π 0 F -F ), v ∈ R, (43)
originating from (1) in the homogeneous case. We detail the numerical scheme in this simplified configuration in Appendix 7. is chosen equal to 10 -2 and 10 -4 ; we observe that at the first iteration, the solution is very close to its projection (since E is nearly zero), even if the initial condition is not at equilibrium. This emphasizes the AP character of our numerical scheme.

Non homogeneous, constant τ case

In this subsection, we look at the non homogeneous case for constant τ . We consider the following initial condition

F 0 (x, v) = (1 -α(x))M 1 (x, v) + α(x)M 2 (x, v), x ∈ [0, 4π], v ∈ R, with M i (x, v) = 1 √ 2π exp(-|v -u i (x)| 2 /2). We choose u 2 (x) = 2, u 1 (x) = α(x)u 2 (x)/(α(x) -
(t) = (F -Π 0 F )(t) L 2 .
From top-left to bottom-right η = 1, 10 -1 , 10 -2 , 10 -5 . In Figure 9, the time history of E(t) = (F -Π 0 F )(t) L 2 (the L 2 norm is considered in x and v) is shown, for different values of η. The new scheme "EXPO" is compared to RK3 and RK1. The new scheme turns out to be competitive in the kinetic regime since the results are very close to those obtained with RK3 or RK1. When η is too small (η = 10 -5 ), the explicit scheme RK3 requires a very small time step so that only results obtained with EXPO scheme is shown. A very strong relaxation can be observed, confirming the fact that the solution is projected onto the set of isotropic functions at the first iteration.

Non homogeneous, non constant τ case

In this subsection, we look at the non homogeneous case for a non constant τ . We consider the following initial condition

F 0 (x, v) = (α(x)M 1 (x, v) + (1 -α(x))M 2 (x, v))(1 + |v|), x ∈ [0, 4π], v ∈ R, with M i (x, v) = 1 √ 2π exp(-|v -u i (x)| 2 /2). We choose u 2 (x) = -2, u 1 (x) = (1 + 0.3 cos(x/2))/2, and α(x) = (3u 2 (x) + u 2 (x) 3 )/(3u 2 (x) + u 2 (x) 3 -3u 1 (x) -u 1 (x) 3 ).
The numerical parameters are the following: ∆t = 0.2∆x/v max , N x = 512, N v = 512 and v max = 16.

(t) = (F -Π 0 F )(t) L 2 .
From top-left to bottom-right η = 1, 10 -1 , 10 -2 , 10 -3 . In Figure 10, the time history of E(t) = (F -Π 0 F )(t) L 2 (the L 2 norm is considered in x an v), for different values of η and different numerical schemes (Runge-Kutta 1 and 3, EXPO) is shown. For η = 10 -2 , comparison with explicit RK method needs ∆t = 0.05∆x/v max whereas the "EXPO" method still uses ∆t = 0.2∆x/v max . We observe that the three methods have a similar behavior which validates the EXPO method for η = 1, 10 -1 , 10 -2 . When η is smaller, explicit methods are too costly so that only results obtained by the EXPO method are presented. In particular, in Figure 11 when η = 10 -4 , one can see that the initial condition is immediately projected onto isotropic functions.

Conclusion

In this work, a new scheme is proposed for the numerical simulation of a kinetic theory based model of turbulence. The derivation of this scheme needed two important tools: (i) the change of frame enables to avoid the delicate issue of the computation of u; (ii) the improvement of standard AP schemes since the change of frame makes an additional stiff transport term appear.

Even if this scheme is implicit to avoid the severe constraint coming from the stiff relaxation source term and the additional stiff transport term, its numerical cost is the same as an explicit scheme; moreover, it enjoys the Asymptotic Preserving property, i.e., it degenerates when the stiff parameter goes to zero, to a numerical scheme consistent with the asymptotic model. The current work presents first results obtained with a prototype one-dimensional model. More general configurations can be tackled in future, in order to study turbulence from kinetic theory in a deeper way. In addition, the spatial discretization can be improved to capture physical turbulence phenomena; as mentioned above, high order WENO schemes can be quite easily adapted to our time discretization strategy.

Appendix

Computation of pressure

The equation for the computation of pressure P is presented here in detail. Starting from [START_REF] Degond | Existence of solutions of a kinetic equation modeling cometary flows[END_REF], we can obtain the following equality

v⊗(v+u)∂ x F dv-v(E+(∂ x u) T v)•∂ v F dv+[(u•∂ x )u] = ∂ x • (v + u) ⊗ (v + u)F dv +E-[(∂ x •u)u].
(44) Using the relations

-v(E + (∂ x u) T v) • ∂ v F dv = E + (∂ x u) T vF dv + ∂ x • u vF dv v ⊗ (v + u)∂ x F dv = ∂ x • v ⊗ (v + u)F dv -∂ x • u vF dv we can obtain v ⊗ (v + u)∂ x F dv -v(E + (∂ x u) T v) • ∂ v F dv + [(u • ∂ x )u] = ∂ x • v ⊗ (v + u)F dv -∂ x • u vF dv + E + (∂ x u) T vF dv + ∂ x • u vF dv + [(u • ∂ x )u] = ∂ x • (v + u) ⊗ (v + u)F dv -∂ x • u ⊗ (v + u)F dv + E + (∂ x u) T vF dv + [(u • ∂ x )u] = ∂ x • (v + u) ⊗ (v + u)F dv -∂ x • (u ⊗ u) -∂ x • u ⊗ vF dv + E + (∂ x u) T vF dv + [(u • ∂ x )u] = ∂ x • (v + u) ⊗ (v + u)F dv -∂ x • (u ⊗ u) + E + [(u • ∂ x )u] = ∂ x • (v + u) ⊗ (v + u)F dv -u(∂ x • u) + E.
Hence, from [START_REF] Degond | Existence of solutions of a kinetic equation modeling cometary flows[END_REF], we deduce by taking the divergence of (44) the equation for P -∆P = ∂ 2

x :

(v + u) ⊗ (v + u)F dv -∂ x • [(∂ x • u)u].

Computation of the matrix involved in the term E

Let us consider the matrix

v τ ⊗ ∂ v F dv v τ ⊗ ∂ v F dv = - 1 τ I -v ⊗ v τ ′ τ F dv = - F τ dv + v ⊗ v τ τ ′ F dv =: -S.
For an arbitrary vector ξ ∈ R d , let us consider the quantity (Sξ) • ξ

(Sξ) • ξ = F τ dv|ξ| 2 -(v • ξ) 2 τ ′ τ F dv = F τ |ξ| 2 -(v • ξ) 2 τ ′ dv.
For a decreasing function τ (τ has a power law shape τ (ξ) = ξ α , -1 ≤ α ≤ -1/2), this quantity is positive, which ensures the invertibility of the matrix v τ ⊗ ∂ v F dv.

Numerical scheme in the one-dimensional homogeneous case

Let us detail the steps of the algorithm in the simplified one-dimensional and homogeneous framework

∂ t F - 1 η (v/τ 2 )F (v/τ )∂ v F ∂ v F = 1 ητ (Π 0 F -F ), v ∈ R.
First, we rewrite this equation ∂ t (e t/(ητ ) F = e t/(ητ ) η (v/τ 2 )F (v/τ )∂ v F ∂ v F + e t/(ητ ) ητ Π 0 F, so that the equation to compute F n+1 writes (following (27))

F n+1 = e -∆t/(ητ ) F n + ∆t η e -∆t/(ητ ) (v/τ 2 )F n+1 (v/τ )∂ v F n ∂ v F n + (1 -e -∆t/(ητ ) )Π 0 F n+1 .

To compute (v/τ 2 )F n+1 , we use the following relation (see (30)), obtained by multiplying the last equation by (v/τ 2 ) and integrating with respect to v

(v/τ 2 )F n+1 = 1 - ∆t η (v/τ 2 )e -∆t/(ητ ) ∂ v F n (v/τ )∂ v F n -1
(v/τ 2 )e -∆t/(ητ ) F n = η (v/τ )∂ v F n η (v/τ )∂ v F n -∆t (v/τ 2 )e -∆t/(ητ ) ∂ v F n (v/τ 2 )e -∆t/(ητ ) F n .

The following equality enables us to compute Π 0 F n+1 (following (32))

Π 0 F n+1 = Π 0 F n + ∆t η (v/τ 2 )F n+1 (v/τ )∂ v F n Π 0 (∂ v F n ).

Figure 1 :Figure 2 :Figure 3 :Figure 4 :

 1234 Figure 1: τ = 1 constant and homogeneous case: time evolution ofE(t) = (F -Π 0 F )(t) L 2 .Comparison of the different schemes for η = 1.

3 .Figure 5 :Figure 6 :

 356 Figure 5: τ non constant and homogeneous case: time evolution of E. Comparison of the different schemes for η = 1.

  [START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving compressible Navier-Stokes asymptotics[END_REF] and α(x) = (1 + 0.3 cos(x/2))/2.The numerical parameters are the following: ∆t = 0.2∆x/v max , N x = 64, N v = 256 and v max = 16.

Figure 7 :

 7 Figure 7: τ non constant and homogeneous case: time evolution ofE(t) = (F -Π 0 F )(t) L 2 .RK3 (with ∆t = 0.01) and EXPO schemes for η = 0.5 (left) and η = 0.1 (right).

Figure 8 :Figure 9 :

 89 Figure 8: τ non constant and homogeneous case: time evolution of E. EXPO scheme for η = 10 -2 and η = 10 -4 .

Figure 10 :

 10 Figure 10: τ non constant and non homogeneous case: time evolution of E(t) = (F -Π 0 F )(t) L 2 . From top-left to bottom-right η = 1, 10 -1 , 10 -2 , 10 -3 .

Figure 11 :

 11 Figure 11: τ non constant and non homogeneous case: time evolution ofE(t) = (F -Π 0 F )(t) L 2 . η = 10 -4 .

Acknowledgment

The first four authors thank the Indo-French Centre for Applied Mathematics (IFCAM) for funding a collaborative research project which made this work possible.

with M i (v) = 1 √ 2π exp(-|v -u i | 2 /2), u i , i = 1, 2 satisfies u 1 = 1, u 2 = -2 and α = 7/9 (this value of α is chosen such that R (v/τ )F 0 (v)dv = 0). The τ function is τ (ξ) = (1 + √ ξ)/ξ which is physically relevant.

The following numerical parameters are used: N v = 256, ∆t = 0.1, v max = 16. In Figure 5, the time history of E(t) = F (t) -Π 0 F (t) L 2 is plotted for different numerical schemes (explicit RK schemes and EXPO) for η = 1. The same conclusions as before follow since EXPO and RK's are nearly superimposed. Note that the rate of relaxation is quite different from the previous case. Figure 6 displays the distribution function F (t, v) for different times t = 0, 1, 2, 5, 100, 1000 using EXPO scheme for η = 1. This illustrates the relaxation towards the isotropic state. The rate of isotropisation is not the same for all the velocities since we are dealing with a non constant relaxation time τ (see also [START_REF] Mieussens | Numerical comparison of BGK-models with proper Prandtl number[END_REF]). In Figure 7, E is plotted for different values of η using the EXPO scheme and RK3. For this latter scheme, the time step is fixed as ∆t = 0.01 for stability, whereas ∆t = 0.1 for EXPO. We observe for η = 0.5 that EXPO is very close to RK3, even if the time step is 10 times smaller for RK3. For η = 0.1, the two curves are not so close (Figure 8, right) and we also display the result obtained by EXPO with ∆t = 0.01; we then observe that the two curves are nearly the same, which confirms the convergence in time of EXPO. In Figure 8, results obtained with EXPO (∆t = 0.1) are shown, where the value of η


Nicolas Crouseilles

Pierre Navaro

Eric Sonnendrücker

Charge conserving grid based methods for the Vlasov-Maxwell equations

Keywords: Maxwell-Vlasov system, discrete continuity equation, Finite Volume method, Spectral method, semi-Lagrangian method

published or not. The documents may come

Charge-conserving grid based methods for the

Vlasov-Maxwell equations

Introduction

We consider the motion of particles in their self-consistent electromagnetic field, which can be described by the Vlasov-Maxwell equations. In this work, we restrict ourselves to the two-dimensional non relativistic model, which involves four phase space dimensions, namely x, y, v x , v y , but the ideas developed here can be extended in a straightforward manner to the 3D relativistic model. In our case the unknown quantities are the particle distribution function f (t, x, v), the electric field E(t, x) = (E x , E y , 0) and the magnetic field B(t, x) = (0, 0, B z ), where we denote by x = (x, y) and v = (v x , v y ). Then the Vlasov equation reads

∂f ∂t + v • ∇ x f + (E(t, x) + v × B(t, x)) • ∇ v f = 0, (1)
where ∇ x denotes the gradient in configuration space and ∇ v denotes the gradient in velocity space. The initial condition f (0, x, v) = f 0 (x, v) is given. The self-consistent electromagnetic field (E(t, x), B(t, x)) is computed thanks to Maxwell's equations

∂E ∂t -∇ x × B = -J, (2)
∂B ∂t + ∇ x × E = 0, (3)
∇ x • E = ρ, (4)
∇ x • B = 0. ( 5
)
The source of Maxwell's equations, namely the charge density ρ and current density J = (J x , J y ) are computed from the particle distribution f and the uniform neutralizing Maxwellian background particles defined by their density n b (x) = f 0 (x, v) dv thanks to

ρ(t, x) = R 2 f (t, x, v) dv -n b (x), J(t, x) = R 2 f (t, x, v)v dv. (6)
Note that integrating the Vlasov equation [START_REF] Barthelmé | Le problème de conservation de la charge dans le couplage des équations de Vlasov et de Maxwell[END_REF] with respect to the velocity variable v yields ∂ρ ∂t + ∇ x • J = 0, [START_REF] Crouseilles | Hamiltonian splitting for the Vlasov-Maxwell equations[END_REF] which is called the continuity equation and expresses the local conservation of charge. In the continuous setting, provided this continuity equation is satisfied and Gauss' law ( 4) is satisfied at time t = 0, if the electric and magnetic fields are computed using only Ampère's law (2) and Faraday's law [START_REF] Bouchut | On the discrete conservation of the Gauss-Poisson equation of plasma Physics[END_REF], Gauss' law is satisfied for all time. In general, numerical Vlasov-Maxwell solvers being PIC, Eulerian or semi-Lagrangian do not verify this and a correction scheme is generally used to enforce Gauss' law at each time step or from time to time. The aim of this paper is to propose a general procedure, valid for Vlasov solvers on a cartesian grid based on a splitting method between configuration and velocity space, that enables to establish a discrete continuity equation verified by the discrete charge and current densities computed from the Vlasov solver compatible with the Maxwell solver, so that only Ampère's and Faraday's law will need to be advanced by our Maxwell solver, Gauss' law being a consequence of those.

A number of work has been proposed to tackle this problem in PIC solvers (see [START_REF] Villasenor | Rigorous charge conservation for local electromagnetic field solvers[END_REF][START_REF] Chen | An energy-and charge-conserving, implicit, electrostatic particle-in-cell algorithm[END_REF][START_REF] Umeda | A new charge conservation method in electromagnetic particle-in-cell simulations[END_REF][START_REF] Marder | A method for incorporating Gauss's law into electromagnetic PIC codes[END_REF][START_REF] Zh | Exact charge conservation scheme for Particle-in-Cell simulation with an arbitrary form-factor[END_REF][START_REF] Barthelmé | Le problème de conservation de la charge dans le couplage des équations de Vlasov et de Maxwell[END_REF][START_REF] Barthelmé | Numerical charge conservation in particlein-cell codes, Numerical methods for hyperbolic and kinetic problems[END_REF]) but also in the Eulerian case. In this latter case, one can cite [START_REF] Crouseilles | A charge preserving scheme for the numerical resolution of the Vlasov-Ampère equations[END_REF] in which charge is preserved using a forward semi-Lagrangian Vlasov solver and [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D[END_REF] in which a conservative PPM method is used in such a way that the Gauss law is ensured. The present work combines ideas of [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D[END_REF] and [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for the Vlasov equation[END_REF]. Indeed, in [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for the Vlasov equation[END_REF], a general framework for conservative semi-Lagrangian methods is proposed so that [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D[END_REF] can be generalized to a class of Vlasov solvers. We also mention [START_REF] Crouseilles | Hamiltonian splitting for the Vlasov-Maxwell equations[END_REF] in which a Hamiltonian splitting is proposed for the Vlasov-Maxwell equations that enables to derive charge preserving high order in time methods.

The outline of the paper is as follows. First we will derive the discrete continuity equations associated to the Yee Maxwell solver and to a spectral solver, then we will introduced an algorithm for computing the discrete charge and current densities compatible with these discrete continuity equations for a conservative semi-Lagrangian algorithm and a spectral algorithm. Finally, after presenting the coupling algorithm, we will validate our method on some relevant test cases.

Discrete continuity equations

A discrete continuity equation is necessarily linked to the Maxwell solver, as the discrete curl and divergence operators need to be compatible. To this aim a first requirement of the Maxwell solver, is that the discrete divergence of the discrete curl vanishes. In this case a discrete continuity equation specific for each solver can be derived such that if this discrete continuity equation is satisfied and the Gauss law is satisfied at time step t n = n∆t it will also be satisfied at time step t n+1 .

The formulation of the discrete continuity equation for the Yee scheme, as well as for Finite Volume schemes for the Maxwell system on unstructured grids with the leap-frog and other time stepping schemes have already been obtained by Bouchut [START_REF] Bouchut | On the discrete conservation of the Gauss-Poisson equation of plasma Physics[END_REF].

We shall recall here the discrete continuity equation for the classical Yee solver and also introduce it for a spectral solver associated to a leap-frog method in time.

For the Yee Maxwell solver

Let us first consider the classical Yee solver on a staggered mesh. First, we introduce a mesh in the x and y directions: x i = x min + i∆x, ∆x = (x max -x min )/N x and y j = y min + j∆y, ∆y = (y max -y min )/N y . Then, we define the center of the cells x i+1/2 = x i + ∆x/2 and y j+1/2 = y j + ∆y/2. We also introduce the time discretization: t n = n∆t, n ∈ N, ∆t > 0, and t n+1/2 = t n + ∆t/2. The Yee scheme is satisfied by

E n x i+1/2,j ≃ E x (t n , x i+1/2 , y j ), E n y i,j+1/2 ≃ E y (t n , x i , y j+1/2 ), B n-1/2 z i+1/2,j+1/2 ≃ B z (t n-1/2 , x i+1/2 , y j+1/2 ),
for a given current J n+1/2

x i+1/2,j ≃ J x (t n+1/2 , x i+1/2 , y j ) and J n+1/2 y i,j+1/2 ≃ J y (t n+1/2 , x i , y j+1/2 ). Figure 1 displays the positons of the different components of the fields on a staggered cartesian grid in the Yee scheme. The Yee scheme approximating (2)-( 3) reads

E n+1 x i+1/2,j -E n x i+1/2,j ∆t = B n+1/2 z i+1/2,j+1/2 -B n+1/2 z i+1/2,j-1/2 ∆y -J n+1/2 x i+1/2,j , (8)
E n+1 y i,j+1/2 -E n y i,j+1/2 ∆t = - B n+1/2 z i+1/2,j+1/2 -B n+1/2 z i-1/2,j+1/2 ∆x -J n+1/2 y i,j+1/2 , (9)
B n+1/2 z i+1/2,j+1/2 -B n-1/2 z i+1/2,j+1/2 ∆t = E n x i+1/2,j+1 -E n x i+1/2,j ∆y - E n y i+1,j+1/2 -E n y i,j+1/2 ∆x . (10)
The associated discrete Gauss' law at time t n will then read

E n x i+1/2,j -E n x i-1/2,j ∆x + E n y i,j+1/2 -E n y i,j-1/2 ∆y = ρ n i,j . (11)
Now, taking the discrete divergence of Ampère's law ( 8)-( 9) yields 1 ∆t

E n+1 x i+1/2,j -E n+1 x i-1/2,j ∆x + E n+1 y i,j+1/2 -E n+1 y i,j-1/2 ∆y - E n x i+1/2,j -E n x i-1/2,j ∆x - E n y i,j+1/2 -E n y i,j-1/2 ∆y = - J n+1/2 x i+1/2,j -J n+1/2 x i-1/2,j ∆x + J n+1/2 y i,j+1/2 -J n+1/2 y i,j-1/2 ∆y . (12)
Let us now introduce, at time step t n the following discrete continuity

ρ n+1 i,j -ρ n i,j ∆t + J n+1/2 x i+1/2,j -J n+1/2 x i-1/2,j ∆x + J n+1/2 y i,j+1/2 -J n+1/2 y i,j-1/2 ∆y = 0. ( 13
)
If this is satisfied, it follows from (12) that 1 ∆t

E n+1 x i+1/2,j -E n+1 x i-1/2,j ∆x + E n+1 y i,j+1/2 -E n+1 y i,j-1/2 ∆y - E n x i+1/2,j -E n x i-1/2,j ∆x - E n y i,j+1/2 -E n y i,j-1/2 ∆y = ρ n+1 i,j -ρ n i,j ∆t .
So, if the discrete Gauss law ( 11) is satisfied at time t n , it follows that

E n+1 x i+1/2,j -E n+1 x i-1/2,j ∆x + E n+1 y i,j+1/2 -E n+1 y i,j-1/2 ∆y = ρ n+1 i,j ,
which is the discrete Gauss law at time t n+1 .

For a spectral Maxwell solver

We consider periodic boundary conditions in x and y so that discrete Fourier transform is adapted to the spatial discretization of the Vlasov and the Maxwell's equations. We denote by E n x,kx,ky the spatial discrete Fourier transform of E x (t n ) (similar notations are used for E y , B z , ρ, J x and J y ), with t n = n∆t, n ∈ N, ∆t > 0. Then the Fourier transform in (x, y) of the Faraday and Ampère's equations together with a leap-frog scheme in time lead to the following Maxwell solver:

E n+1 x,
B n+1/2 z,kx,ky -B n-1/2 z,kx,ky ∆t = ik y E n x,kx,ky -ik x E n y,kx,ky . (15)
The discrete spectral Gauss' law at time t n writes

ik x E n x,kx,ky + ik y E n y,kx,ky = ρ n kx,ky . (17)
Taking the discrete Ampère law in the spectral case, amounts to summing (14) multiplied by ik x and (15) multiplied by ik y . This leads to

1 ∆t (ik x E n+1 x,kx,ky +ik y E n+1 y,,kx,ky -ik x E n x,kx,ky -ik y E n y,kx,ky ) = -(ik x J n+1/2
x,kx,ky +ik y J n+1/2 y,kx,ky ).

Notice that here also, like in the Yee scheme the discrete curl of B z vanishes in the last equation, which is a necessary property for our Maxwell solvers if we want a discrete continuity equation. Since the Gauss law ( 17) is satisfied at time t n , we then have

ik x E n+1 x,kx,ky + ik y E n+1 y,kx,ky = ρ n kx,ky -∆t(ik x J n+1/2 x,kx,ky + ik y J n+1/2 y,kx,ky ). ( 18
)
The discrete version of the continuity equation ( 7) reads here

ρ n+1 kx,ky = ρ n kx,ky -∆t(ik x J n+1/2 x,kx,ky + ik y J n+1/2 y,kx,ky ). ( 19
)
If this is ensured, if follows from (18) that

ik x E n+1 x,kx,ky + ik y E n+1 y,kx,ky = ρ n+1 kx,ky ,
which is the Gauss law is satisfied at time t n+1 .

Computation of the discrete charge and current density from the Vlasov solver compatible with the discrete continuity equation

Having identified the discrete continuity equation for the Maxwell solver that we want to use, the next task is to find a way to compute the discrete charge and current densities from the Vlasov equation so that they satisfy this discrete continuity equation.

For Particle In Cell (PIC) Vlasov solvers associated to the Yee Maxwell solver there is a well-known procedure for doing that introduced by Villasenor and Buneman [START_REF] Villasenor | Rigorous charge conservation for local electromagnetic field solvers[END_REF] generalised to higher order deposition schemes in [START_REF] Barthelmé | Numerical charge conservation in particlein-cell codes, Numerical methods for hyperbolic and kinetic problems[END_REF][START_REF] Barthelmé | Le problème de conservation de la charge dans le couplage des équations de Vlasov et de Maxwell[END_REF]. Other options are also possible, like those introduced by Esirkepov [START_REF] Zh | Exact charge conservation scheme for Particle-in-Cell simulation with an arbitrary form-factor[END_REF], Umeda [START_REF] Umeda | A new charge conservation method in electromagnetic particle-in-cell simulations[END_REF] or Chen et al. [START_REF] Chen | An energy-and charge-conserving, implicit, electrostatic particle-in-cell algorithm[END_REF].

In the case of grid based Vlasov solvers, Sircombe and Arber [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D[END_REF] showed that this could be enforced for a split-eulerian Vlasov solver using the PPM method of Colella and Woodward [START_REF] Colella | The Piecewise Parabolic Method (PPM) for gas-dynamical simulations[END_REF] by computing J from the fluxes needed by the algorithm in the configuration space advections. Let us now show that this idea can be applied for general Finite Volume schemes, including as we will see conservative semi-Lagrangian schemes, following [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for the Vlasov equation[END_REF].

A splitting of the Vlasov equation between configuration space and velocity space advection will lead us to solve alternatively

∂f ∂t + v • ∇ x f = 0, (20)
and

∂f ∂t + (E(t, x) + v × B(t, x)) • ∇ v f = 0, (21)
The second equation does not modify ρ as can be seen integrating with respect to the velocity v. Hence it should not provide any direct contribution to J either for the discrete continuity equation to be satisfied. We shall only require that the discrete solver exactly conserves ρ at all spatial grid points, which is the case for all conservative solvers. So assuming a conservative solver is used for the velocity space advection, our problem is now to design a configuration space advection compatible with the discrete continuity equation given by the Maxwell solver. We first notice that the advection in configuration space (20) is for each given v via constant coefficient advection. Note that this would be true also in the relativistic case for each given p.

We shall further split the configuration space advection for fixed v into two 1D advections ( in x and y directions), which is simpler and less costly.

Note that this does not introduce an additional splitting error as constant coefficient advections commute.

Finally, we only need to consider constant coefficient 1D advections of the form ∂f ∂t

+ a ∂f ∂x = 0, (22)
with a being v x or v y at a given velocity grid point. We shall provide compatible schemes in the case of the Finite Volumes and spectral methods we consider.

Finite volume schemes

Let us define uniformly spaced control cells [x i-1/2 , x i+1/2 ] with x i+1/2x i-1/2 = ∆x. The unknown in a Finite Volume scheme will be the average value on a control cell that we shall denote by

f i = 1 ∆x x i+1/2 x i-1/2 f (x) dx.
Integrating (22) on a control cell and between time t n and t n+1 yields

f n+1 i = f n i - a ∆x t n+1 t n (f (t, x i+1/2 ) -f (t, x i-1/2 )) dt.
Let us denote by f n+1/2 i+1/2 a numerical approximation of 1 ∆t

t n+1 t n f (t, x i+1/2
) dt. Then our Finite Volume scheme becomes

f n+1 i = f n i - a∆t ∆x (f n+1/2 i+1/2 -f n+1/2 i-1/2 ). ( 23
)
Note that the actual computation of f n+1/2 i+1/2 makes up the specific Finite Volume scheme, this could be done with an upwind scheme or the PPM scheme or others. For our purposes, it will be sufficient to consider this generic form of the Finite Volume scheme.

Let us now come back to our split Vlasov solver, for which a Finite Volume scheme of the form (23) will be used in the x and y advection steps and any conservative scheme in the v advection step. We discretize the distribution function on a 4D grid with uniform steps in each direction, denoting by i the x index, j the y index, k the v x index and ℓ the v y index such that f n i,j,k,ℓ ≃ f (t n , x i , y j , v x k , v y ℓ ). Starting from f n i,j,k,ℓ at time step t n , the algorithm reads,

f n,1 i,j,k,ℓ = f n i,j,k,ℓ - v x k ∆t 2∆x (f n,1/2 i+1/2,j,k,ℓ -f n,1/2 i-1/2,j,k,ℓ ), f n,2 i,j,k,ℓ = f n,1 i,j,k,ℓ - v y ℓ ∆t 2∆y (f n,3/2 i,j+1/2,k,ℓ -f n,3/2 i,j-1/2,k,ℓ ), f n,3 i,j,k,ℓ ← f n,2 i,j,k,ℓ using a conservative advection in v space, f n,4 i,j,k,ℓ = f n,3 i,j,k,ℓ - v y ℓ ∆t 2∆y (f n,7/2 i,j+1/2,k,ℓ -f n,7/2 i,j-1/2,k,ℓ ),
f n+1 i,j,k,ℓ = f n,4 i,j,k,ℓ - v x k ∆t 2∆x (f n,9/2 i+1/2,j,k,ℓ -f n,9/2 i-1/2,j,k,ℓ ).
Now, the discrete charge density is linked to the discrete distribution function by

ρ n i,j = ∆v x ∆v y k ℓ f n i,j,k,ℓ ,
where ∆v x = (v x,max -v x,min )/N vx (resp. ∆v y = (v y,max -v y,min )/N vy ) denotes the velocity step in the v x (resp. v y ) direction. So, using the previous algorithm we can relate ρ n+1 i,j to ρ n i,j by summing the different lines with respect to k, ℓ. Then, we get

1 ∆v x ∆v y ρ n+1 i,j = k,ℓ f n+1 i,j,k,ℓ = k,ℓ f n,4 i,j,k,ℓ - ∆t 2∆x k,ℓ v x k (f n,9/2 i+1/2,j,k,ℓ -f n,9/2 i-1/2,j,k,ℓ ) = k,ℓ f n,3 i,j,k,ℓ - ∆t 2∆y k,ℓ v y ℓ (f n,7/2 i,j+1/2,k,ℓ -f n,7/2 i,j-1/2,k,ℓ ) - ∆t 2∆x k,ℓ v x k (f n,9/2 i+1/2,j,k,ℓ -f n,9/2 i-1/2,j,k,ℓ )
Then the conservativity of the advection in v space yields k,ℓ f n,3 i,j,k,ℓ = k,ℓ f n,2 i,j,k,ℓ . Thus, proceeding in the same manner with the first two steps of the algorithm we finally get

1 ∆v x ∆v y ρ n+1 i,j = k,ℓ f n i,j,k,ℓ - ∆t 2∆y k,ℓ v y ℓ (f n,7/2 i,j+1/2,k,ℓ + f n,3/2 i,j+1/2,k,ℓ -f n,7/2 i,j-1/2,k,ℓ -f n,3/2 i,j-1/2,k,ℓ ) - ∆t 2∆y k,ℓ v x k (f n,9/2 i+1/2,j,k,ℓ + f n,1/2 i+1/2,j,k,ℓ -f n,9/2 i-1/2,j,k,ℓ -f n,1/2 i-1/2,j,k,ℓ ). ( 24
)
Let us know denote by

J n+1/2 x i+1/2,j = ∆v x ∆v y k,ℓ v x k 1 2 (f n,9/2 i+1/2,j,k,ℓ + f n,1/2 i+1/2,j,k,ℓ ), ( 25
)
J n+1/2 y i,j+1/2 = ∆v x ∆v y k,ℓ v y ℓ 1 2 (f n,7/2 i,j+1/2,k,ℓ + f n,3/2 i,j+1/2,k,ℓ ). (26)
Then (24) becomes

ρ n+1 i,j -ρ n i,j ∆t + J n+1/2 x i+1/2,j -J n+1/2 x i-1/2,j ∆x + J n+1/2 y i,j+1/2 -J n+1/2 y i,j-1/2 ∆y = 0,
which is exactly the discrete continuity equation ( 13) needed by the Yee scheme. Hence expressions (25)-( 26) provide an expression of the discrete current density J consistent with a Finite Volume Vlasov solver and that satisfies a discrete continuity equation.

Finite Volume form of a semi-Lagrangian scheme

We have obtained in the previous section an explicit expression for the discrete current which satisfies a discrete continuity equation for any split Vlasov solver for which the x and y advection parts can be cast in the generic Finite Volume framework (23). We are now going to cast a conservative semi-Lagrangian scheme [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for the Vlasov equation[END_REF] for constant coefficient advection into this Finite Volume formalism. This will then enable us to construct charge conserving algorithms for a large class of split semi-Lagrangian schemes.

In both Finite Volume and conservative semi-Lagrangian schemes, the first step is to reconstruct a piecewise polynomial function on each cell. We shall call f R this reconstructed piecewise polynomial function. The reconstruction scheme does not matter for our purpose, it could be PPM, splines or something else, with limiters or not. The only property we shall need is that it is linked to the computed cell averages

f n i = 1 ∆x x i+1/2
x i-1/2 f (t n , x) dx, by

f n i = 1 ∆x x i+1/2 x i-1/2 f R (x) dx.
We shall also assume that the CFL condition a∆t ∆x ≤ 1 is verified. In a Finite Volume scheme for the 1D advection equation ( 22), we then compute

f n+1/2 i+1/2 = 1 ∆t t n+1 t n f (t, x i+1/2 ) dt = 1 ∆t t n+1 t n f (t n , x i+1/2 -a(t -t n )) dt ≃ 1 ∆t t n+1 t n f R (x i+1/2 -a(t -t n )) dt,
and by the change of variables

x = x i+1/2 -a(t -t n ) f n+1/2 i+1/2 = 1 a∆t x i+1/2 x i+1/2 -a∆t f R (x) dx.
Then using the formulation (23), we get the Finite Volume scheme

f n+1 i = f n i - 1 ∆x x i+1/2 x i+1/2 -a∆t f R (x) dx - x i-1/2 x i-1/2 -a∆t f R (x) dx . (27)
On the other hand, for a conservative semi-Lagrangian scheme, the distribution function is updated using the relation

f n+1 i = 1 ∆x X(x i-1/2 ) X(x i-1/2 ) f R (x) dx,
where X(x i+1/2 ) is the origin of the characteristic ending at x i+1/2 , that is in our case of constant advection at velocity a, X(x i+1/2 ) = x i+1/2 -a∆t.

Hence

f n+1 i = 1 ∆x x i+1/2 -a∆t x i-1/2 -a∆t f R (x) dx = 1 ∆x x i-1/2 x i-1/2 -a∆t f R (x) dx + x i+1/2 x i-1/2 f R (x) dx - x i+1/2 x i+1/2 -a∆t f R (x) dx = f n i + 1 ∆x x i-1/2 x i-1/2 -a∆t f R (x) dx - x i+1/2 x i+1/2 -a∆t f R (x) dx
which is the same expression as (27), so that both formalisms yield the same numerical scheme. Moreover, the finite volume flux f n+1/2 i+1/2 can be expressed for a semi-Lagrangian scheme with respect to the reconstructed function f R by

f n+1/2 i+1/2 = 1 a∆t x i+1/2 x i+1/2 -a∆t f R (x) dx.
In particular, if the reconstruction is performed using a primitive F R of f R , we have

f n+1/2 i+1/2 = 1 a∆t F R (x i+1/2 ) -F R (x i+1/2 -a∆t) .
It now remains to provide the fluxes for both half advections in x and y. For the half-advections in x, we have on the one hand, integrating on

[t n , t n+1/2 ] × [x i-1/2 , x i+1/2 ] f n,1 i = f n i - v x ∆t 2∆x t n+1/2 t n f (t, x i+1/2 ) -f (t, x i+1/2 ) dt = f n i - v x ∆t 2∆x f n,1/2 i+1/2 -f n,1/2 i-1/2 , ( 28
) with f n,1/2 i+1/2 = t n+1/2 t n f (t, x i+1/2
)dt. On the other hand

f n,1 i := 1 ∆x x i+1/2 x i-1/2 f n,1 (x)dx = 1 ∆x x i+1/2 -vx∆t/2 x i-1/2-vx∆t/2 f n (x)dx = f n i - 1 ∆x (φ n,1/2 i+1/2 -φ n,1/2 i-1/2 ), ( 29
)
where

φ n,1/2 i+1/2 = x i+1/2-vx∆t/2 x i+1/2
f n (x)dx. Hence, we can identify the fluxes in ( 28) and (29) to get

f n,1/2 i+1/2 = φ n,1/2 i+1/2 /(v x ∆t/2) if v x = 0, 0 if v x = 0.
In the same way for f n,9/2 i+1/2 , we have

f n,9/2 i+1/2 = φ n,9/2 i+1/2 /(v x ∆t/2) if v x = 0, 0 if v x = 0.
Then, it comes

φ n,1/2 i+1/2 = F R (x i+1/2 ) -F R (x i+1/2 -v x ∆t/2) with F R the primitive of f n , φ n,9/2 i+1/2 = F R (x i+1/2 ) -F R (x i+1/2 -v x ∆t/2) with F R the primitive of f n,4 .
For the half-advections in y, we proceed in the same way integrating on [t n , t n+1/2 ] × [y j-1/2 , y j+1/2 ] in order to obtain the finite volume formulation (as ( 28)) and identifying with the conservative semi-Lagrangian formulation (as (29)), we have

f n,3/2 j+1/2 = φ n,3/2 j+1/2 /(v y ∆t/2) if v y = 0, 0 if v y = 0. and f n,7/2 j+1/2 = φ n,7/2 j+1/2 /(v y ∆t/2) if v y = 0, 0 if v y = 0.
Then, we get similarly

φ n,3/2 j+1/2 = F R (y j+1/2 ) -F R (y j+1/2 -v y ∆t/2) with F R the primitive of f n,1 , φ n,7/2 j+1/2 = F R (y j+1/2 ) -F R (y j+1/2 -v y ∆t/2) with F R the primitive of f n,3 .

Spectral scheme

Denoting by f n (k x , k y , v x , v y ) the Fourier transform in space of f (t n ), the 1D advections in configuration space become for the spectral scheme

∂f ∂t + ik x v x f = 0, ( 30
)
for which the exact solution on a time step reads f (t+∆t) = f (t) exp(-ik x v x ∆t), and the same for the advection in y. So, using a spectral scheme in space coupled with a Strang splitting, we get the following algorithm

f n,1 kx,ky,k,ℓ = f n kx,ky,k,ℓ exp(-ik x v x ∆t/2), f n,2 kx,ky,k,ℓ = f n,1 kx,ky,k,ℓ exp(-ik y v y ∆t/2), f n,3
kx,ky,k,ℓ ← f n,2 kx,ky,k,ℓ using a conservative advection in v space, f n,4 kx,ky,k,ℓ = f n,3 kx,ky,k,ℓ exp(-ik y v y ∆t/2), f n+1 kx,ky,k,ℓ = f n,4 kx,ky,k,ℓ exp(-ik x v x ∆t/2).

As we did for the Finite Volume schemes, we need to extract the current from the space advections in a way that is compatible with the discrete continuity equation ( 19) for the Maxwell solver we want to couple to. To this aim we need to make a transport structure appear. Introducing f

n,1/2 kx,ky,k,ℓ f n,1/2 kx,ky,k,ℓ = f n kx,ky,k,ℓ [1 -exp(-ik x v x ∆t/2)] /(ik x ∆t/2) if k x = 0, 0 if k x = 0, (31)
the first advection in x writes

f n,1 kx,ky,k,ℓ = f n kx,ky,k,ℓ -ik x ∆t 2 f n,1/2
kx,ky,k,ℓ . Using the same notations, one has for the last x-advection

f n+1 kx,ky,k,ℓ = f n,4 kx,ky,k,ℓ -ik x ∆t 2 f n,9/2 kx,ky,k,ℓ ,
where f n,9/2 is deduced from (31) by replacing f n by f n,4 . Similarly, for the y-advections, it comes

f n,2 kx,ky,k,ℓ = f n,1 kx,ky,k,ℓ -ik y ∆t 2 f n,3/2 kx,ky,k,ℓ , f n,4 kx,ky,k,ℓ = f n,3 kx,ky,k,ℓ -ik y ∆t 2 f n,7/2 kx,ky,k,ℓ ,
where f n,3/2 (resp. f n,7/2 ) is given by the following expression with d = 0 (resp. d = 2)

f n,d+3/2 kx,ky,k,ℓ = f n,(d+1) kx,ky,k,ℓ [1 -exp(-ik y v y ∆t/2)] /(ik y ∆t/2) if k y = 0, 0 if k y = 0. (32)
Now, let us add up the different contributions to the current. First, the discrete charge density is linked to the discrete distribution function by

ρ n kx,ky = ∆v x ∆v y k ℓ f n kx,ky,k,ℓ .
So, using the previous algorithm we can relate ρ n+1 kx,ky to ρ n kx,ky by summing the different lines with respect to k, ℓ. Then, we get

1 ∆v x ∆v y ρ n+1 kx,ky = k,ℓ f n+1 kx,ky,k,ℓ = k,ℓ f n,4 kx,ky,k,ℓ -ik x ∆t 2 f n,9/2 kx,ky,k,ℓ = k,ℓ f n,3 kx,ky,k,ℓ -ik y ∆t 2 f n,7/2 kx,ky,k,ℓ -ik x ∆t 2 f n,9/2 kx,ky,k,ℓ .
Then the conservativity of the advection in v space yields k,ℓ f n,3 kx,ky,k,ℓ = k,ℓ f n,2 kx,ky,k,ℓ . Thus, proceeding in the same manner with the first two steps of the algorithm we finally get

1 ∆v x ∆v y ρ n+1 kx,ky = k,ℓ f n kx,ky,k,ℓ -ik y ∆t 2 k,ℓ (f n,7/2 kx,ky,k,ℓ + f n,3/2 kx,ky,k,ℓ ) -ik x ∆t 2 k,ℓ (f n,9/2 kx,ky,k,ℓ + f n,1/2 kx,ky,k,ℓ ). ( 33
)
Let us now denote by

J n+1/2 x,kx,ky = ∆v x ∆v y k,ℓ 1 2 (f n
,9/2 kx,ky,k,ℓ + f n,1/2 kx,ky,k,ℓ ), (34)
J n+1/2 y,kx,ky = ∆v x ∆v y k,ℓ 1 2 (f n
,7/2 kx,ky,k,ℓ + f n,3/2 kx,ky,k,ℓ ). (35)
Then (33) becomes

ρ n+1 kx,ky -ρ n kx,ky ∆t + ik x J n+1/2
x,kx,ky + ik y J n+1/2 kx,ky = 0, For more details, see [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for the Vlasov equation[END_REF]. Let us remark that high order upwind reconstruction LAG(2d + 1) and the extension to the two dimensional case are possible (see [START_REF] Crouseilles | Semi-Lagrangian simulations on polar grids: from diocotron instability to ITG turbulence[END_REF]).

Coupling the Vlasov solver with the Maxwell solver

The only point that is not straightforward in the coupling of our Vlasov solvers with their corresponding Maxwell solvers, is how to get the electric field E n+1/2 , needed for the velocity space advection before the full current has been computed. We use here for both our methods the strategy suggested in [START_REF] Sircombe | VALIS: A split-conservative scheme for the relativistic 2D[END_REF]. Ampère is advanced on ∆t/2 using the predict currents J n computed using

J n x,i+1/2,j = ∆v x ∆v y k,ℓ v x k 2 (f n i,j,k,ℓ + f n i+1,j,k,ℓ ), (36)
J n y,i,j+1/2 = ∆v x ∆v y k,ℓ v y k 2 (f n i,j,k,ℓ + f n i,j+1,k,ℓ ). (37)
This predicted electric field is then only used for the velocity advection and discarded. The electric field E n+1 will later be recomputed directly from E n using the current J n+1/2 verifying the discrete continuity equation. Another point comes from the fact that, before the advection in v, we have the electric field on staggered grid E n+1/2 x i+1/2,j , E n+1/2 y i,j+1/2 whereas we need it at the grid points (x i , y j ). A simple average enables to compute E n+1/2 x i,j , E n+1/2 y i,j . Finally, the initialization has to be done carefully. Indeed, one has to initialize f 0 such that its density ρ 0 satisfies the appropriate discrete Gauss law. In practice, for a given initial condition f 0 , we compute ρ 0 and solve the Poisson equation to get (E 0

x , E 0 y ) using any solvers. Then, ρ 0 is modified in such a way that it solves exactly the appropriate discrete Gauss law. Once ρ 0 is re-computed, f 0 is modified according to ρ 0 .

We then get in both cases the following algorithm where the half advections in x and y need to use the adequate formula for each solver introduced previously. ) → E n+1

Algorithm

) → E n+1/2 2. half-advection in x → f n,1 and J n,1/2 x (a) compute f n,1/2 and J n,1/2 x = k,ℓ v x k f n,1/2 ∆v x ∆v y (b) half-advection in x.

Numerical results

We focus on classical test cases which have the advantage to be compared easily to a Vlasov-Poisson solver. Indeed, we consider the following class of initial conditions

f 0 (x, v) = E(v)(1 + αg(x, y)),
where g(x, y) is a periodic function and E(v) a velocity function. Typically, we choose for the equilibrium

E(v) = 1 2π e -|v| 2 /2 (Landau) and E(v) = 1 2π e -|v| 2 /2 v 2 x (TSI).
For the spatial dependency, we consider g(x, y) = cos(k x x) ("1Dx" case) to recover the one dimensional case or g(x, y) = cos(k x x) cos(k y y) ("PROD" case). The parameters are chosen as follows: k x = k y = 0.5 = 2π/L x = 2π/L y where L x = [0, 4π] (resp. L y = [0, 4π]) denotes the spatial domain in the x direction (resp. y direction). The domain in velocity is v ∈ [-8, 8] 2 to ensure that the discrete integration in the v directions is sufficiently precise. We compare the present algorithm to a Vlasov-Poisson algorithm (VP) in which the electric field is computed just before the v advection. Spectral schemes. In Figures 2 and3, we plot the time history of the electric energy E x (t) 2 L 2 + E y (t) 2 L 2 for the reference method (VP) and the "new" method (VM), in the Landau case (α = 10 -3 ). Figure 4 shows the time evolution of the electric energy in the TSI case (α = 10 -3 ).

One can observe that the two methods are superimposed, illustrating the equivalence of the two algorithms. In particular, we checked that the electric field computed by VM at time t n+1 and the density ρ n+1 (computed from f n+1 ) solve the discrete Gauss law (which is spectral in this case) up to machine precision and up to the precision of the conservation of the v advection. Two other nice properties of the VM algorithm can be emphasized. First, the one-dimensional cases are preserved. Indeed, for the "1Dx" case, the electric field E y stays to zero for all times, the same being true in the symmetric case "1Dy". Second, if the magnetic field B z is zero at the initial time, this property remains true for all t > 0.

Conservative schemes. In this last part, numerical results obtained by the conservative methods are shown. The numerical parameters are chosen as previously. The properties described previously (symmetry between x and y and B z equal to zero) are also verified in the conservative case. We compare the new algorithm for Vlasov-Maxwell with the Vlasov-Poisson in which the Gauss law is solved using a spectral solver. Hence, the two solvers are not equivalent since the Vlasov-Maxwell solver ensures the following discrete Gauss law

E n x i+1/2,j -E n x i-1/2,j ∆x + E n y i,j+1/2 -E n y i,j-1/2 ∆y = ρ n i,j -1.
However, this difference is quite small and appears (see for example in Figure 6) when the amplitude of the electric energy is below a certain threshold. We can remark that the Vlasov-Maxwell solver benefits from the quality of the high order conservative methods used, which capture the long time behavior of the solution (see Figure 6) in a similar way as the spectral method. For a detailed comparison of these methods (regarding the conserved quantities), we refer to [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for the Vlasov equation[END_REF].

Conclusion

In this work, we have presented a class of Eulerian numerical schemes for the two dimensional Vlasov-Maxwell equations that ensure the discrete continuity equation, and as a consequence, the Gauss law. Moreover, since it is based on a directional splitting, this approach can benefit from an efficient 20 parallelization using domain decomposition and transposition techniques (see [START_REF] Coulaud | Parallelization of semi-Lagrangian Vlasov codes[END_REF]). Physically relevant numerical tests have to be performed; in particular, non periodic boundary conditions have to be adapted to the present approach. Also, higher order methods in time have to be studied to improve the accuarcy.

Figure 1 :

 1 Figure 1: Positions of the different components of the fields in the Yee scheme.

  3. half-advection iny → f n,2 et J n,3/2 y (a) compute f n,3/2 and J n,3/2 y = k,ℓ v y ℓ f n,3/2 ∆v x ∆v y (b) half-advection in y.

4 .

 4 advection in v with E n+1/2 and B n+1/2 z → f n,3 5. half-advection in y → f n,4 and J n,7/2 y (a) compute f n,7/2 and J n,7/2 y = k,ℓ v y ℓ f n,7/2 ∆v x ∆v y (b) half-advection in y. 6. half-advection in x → f n+1 and J n,9/2 x (a) compute f n,9/2 and J n,9/2 x = k,ℓ v x k f n,9/2 ∆v x ∆v y (b) half-advection in x.

  advance Ampère (14)-(15) on ∆t with E n , B

Figure 2 :

 2 Figure 2: "PROD" case. Time history of the electric energy (semi log-scale): comparison of Vlasov-Poisson (VP) and Vlasov-Maxwell (VM). 64 points per spatial direction and 128 points per velocity direction. ∆t = 0.01.

Figure 3 :

 3 Figure 3: "1Dx" case. Time history of the electric energy (semi log-scale): comparison of Vlasov-Poisson (VP) and Vlasov-Maxwell (VM). 64 points per spatial direction and 128 points per velocity direction. ∆t = 0.01.

Figure 4 :

 4 Figure 4: "TSI" case. Time history of the electric energy (semi log-scale): comparison of Vlasov-Poisson (VP) and Vlasov-Maxwell (VM). 64 points per spatial direction and 128 points per velocity direction. ∆t = 0.01.

Figure 5 :

 5 Figure 5: "PROD" case. Time history of the electric energy (semi log-scale): comparison of Vlasov-Poisson (VP) and Vlasov-Maxwell (VM). 64 points per spatial direction and 128 points per velocity direction. ∆t = 0.01. From left to right: SPL, LAG3, PPM2.

Figure 6 :

 6 Figure 6: "TSI" case. Time history of the electric energy (semi log-scale): comparison of Vlasov-Poisson (VP) and Vlasov-Maxwell (VM). 64 points per spatial direction and 128 points per velocity direction. ∆t = 0.01. From left to right: SPL, LAG3, PPM2.

  = k,ℓ v x k f n ∆v x ∆v y and J n y = k,ℓ v y ℓ f n ∆v x ∆v y .

	1. Maxwell prediction → B	n+1/2 z	and (E	n+1/2 x	, E	n+1/2 y	)
	(a) advance Faraday (16) on ∆t with B	n-1/2 z	and E n → B	n-1/2 z
	(b) compute J n x (c) advance Ampère (14)-(15) on ∆t/2 with E n , B	n+1/2 z	and (J n x , J n y
	. Starting from E n = (E n x , E n y ), B and f n+1 as	n-1/2 z	and f n , we compute E n+1 , B	n+1/2 z

Acknowledgements. Computer time grants are kindly provided by the Rechenzentrum Garching (RZG).

which is exactly the discrete continuity equation (19) needed by the spectral Maxwell scheme.

Hence expressions (34)-( 35) provide an expression of the discrete current density J consistent with a spectral Vlasov solver and that satisfies a discrete continuity equation.

Reconstruction

In this section, we recall the different reconstructions we used in this work. Let apply the conservative method to the following problem, using the notations introduced in the previous sections

Introducing the primitive function F (x) =

x x -1/2 f (y)dy, the conservative numerical scheme is

The interpolation conditions for the primitive are F (x i+1/2 ) = i j=0 f n j and we need to reconstruct the primitive F to evaluate at x i±1/2 -a∆t. Following [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for the Vlasov equation[END_REF], we use a Hermite reconstruction. If the feet of the characteristics satisfies

and the reconstruction writes for x ∈ [0, 1[

To approximate the derivative F ′ , we use the following formula leading to different numerical schemes

F ′ (x j+1/2 ) = 1 6 (5f j + 2f j+1 -f j-1 ), F ′ (x j-1/2 ) = 1 6 (5f j -f j+1 + 2f j-1 ) (LAG3).