Nicolas Crouseilles

Xue Hong

Exponential DG methods for Vlasov equations

Keywords: exponential Lawson Runge-Kutta method, discontinuous Galerkin method, Vlasov-Poisson equation, Vlasov Maxwell equations, conservative property

In this work, an exponential Discontinuous Galerkin (DG) method is proposed to solve numerically Vlasov type equations. The DG method is used for space discretization which is combined exponential Lawson Runge-Kutta method for time discretization to get high order accuracy in time and space. In addition to get high order accuracy in time, the use of Lawson methods enables to overcome the stringent condition on the time step induced by the linear part of the system. Moreover, it can be proved that a discrete Poisson equation is preserved. Numerical results on Vlasov-Poisson and Vlasov Maxwell equations are presented to illustrate the good behavior of the exponential DG method.

Introduction

In this work, we are interested in the numerical simulation of Vlasov type equations using Eulerian based methods. Numerical approximation of Vlasov equations has been the subject of a lot works since these models are widely used to describe the dynamics of charged particles in a plasma through a distribution function f (t, x, v) with the time t ≥ 0, x the spatial variables and v the velocity variable. Hence, there have been a lot of numerical methods that have been proposed to numerically solve Vlasov equations. The so-called Particle-In-Cell(PIC) methods [START_REF] Verboncoeur | Particle simulation of plasmas: review and advances[END_REF][START_REF] Kraus | Gempic: geometric electromagnetic particle-in-cell methods[END_REF] in which the unknown is approximated by a sum of Dirac masses with a position and velocity (macro-particles) that solves a differential system. Even if these methods are efficient in high dimensions since only a spatial grid is required, they however suffer from some numerical noise which make them hardly get an accurate approximation. Indeed, the error slowly decreases when the number of macro-particles increases, which turns out to be a drawback in low density plasma region. On the other side, another family of methods have been developed which uses a grid of the phase space (x, v) like spectral methods [START_REF] Klimas | A splitting algorithm for Vlasov simulation with filamentation filtration[END_REF][START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF] or finite differences/volumes methods [START_REF] Banks | High-order accurate conservative finite difference methods for Vlasov equations in 2D+ 2V[END_REF][START_REF] Banks | A new class of nonlinear finite-volume methods for Vlasov simulation[END_REF]. These methods enable to get high order accurate approximation and as such, can capture fine physical phenomena like Landau damping or filamentation in Vlasov equations.

However, due to the phase-space grid, these methods are quite costly both in terms of memory and CPU point of view, in particular when high dimensions are considered. Moreover, their stability is controlled by the so-called CFL condition which imposes a constraint on the time step depending of the phase space mesh refinement, which makes them very costly in practice. To overcome this drawback, semi-Lagrangian have been developed [START_REF] Qiu | A conservative high order semi-Lagrangian WENO method for the Vlasov equation[END_REF][START_REF] Rossmanith | A positivity-preserving high-order semi-Lagrangian discontinuous Galerkin scheme for the Vlasov-Poisson equations[END_REF][START_REF] Einkemmer | A performance comparison of semi-Lagrangian discontinuous Galerkin and spline based Vlasov solvers in four dimensions[END_REF] or arbitrary Lagrangian-Eulerian methods [START_REF] Celia | An Eulerian-Lagrangian localized adjoint method for the advection-diffusion equation[END_REF][START_REF] Wang | A family of eulerian-lagrangian localized adjoint methods for multi-dimensional advection-reaction equations[END_REF][START_REF] Cai | An Eulerian-Lagrangian discontinuous Galerkin method for transport problems and its application to nonlinear dynamics[END_REF][START_REF] Hong | A generalized Eulerian-Lagrangian discontinuous Galerkin method for transport problems[END_REF], which allow extra large time stepping sizes with stability by tracking solutions along their characteristics. For high dimensional problems, since (high order) interpolation techniques are required which leads to huge computational cost. Moreover, reaching high order accuracy in time are quite complicated. One possibility is to use splitting methods which enable to deal with simple subproblems that can even be solved exactly. However, the number of stages required to get high order in time become prohibitive (see [START_REF] Crestetto | Comparison of high-order Eulerian methods for electron hybrid model[END_REF]).

For Vlasov type equations, the linear part induces the most stringent CFL condition since the electromagnetic fields (which induce the nonlinear part) are typically one order of magnitude smaller than the one of the linear advection part. Based on this observation, exponential time integrator have been proposed in which the linear part is solved exactly, and as such do not suffer from the stability condition induced by the linear part, whereas the nonlinear part is solved explicitely. These methods are very popular in a number of applications ( [START_REF] Hochbruck | Exponential integrators[END_REF] and references therein) and enables to derive easily high order methods in time since they are often based on a high order Runge-Kutta method. Regarding the use of such time integrators for Vlasov equations, we can quote [START_REF] Crouseilles | Exponential methods for solving hyperbolic problems with application to collisionless kinetic equations[END_REF][START_REF] Crestetto | Comparison of high-order Eulerian methods for electron hybrid model[END_REF][START_REF] Boutin | Modified lawson methods for Vlasov equations[END_REF], but these works are based on Fourier techniques in space to approximate the linear part, and despite its simplicity and its spectral accuracy, Fourier methods are quite limited in terms of applications (cartesian domains) and suffer from Gibbs phenomenon when non periodic boundary conditions are considered.

In the present work, we focus on exponential type method combined with Discontinuous-Galerkin (DG) method in space to approximate Vlasov type equations. The DG method is a class of finite element methods, in which the approximation space contains completely discontinuous, piecewise polynomials or other basis functions. High order accuracy can be obtained and complex geometries with boundary conditions can be handled (see DG review article [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convectiondominated problems[END_REF]) which is important for physically relevant problems (applications of the DG method to kinetic models also exist, e.g. [START_REF] Cockburn | Locally divergence-free discontinuous Galerkin methods for the Maxwell equations[END_REF][START_REF] Cheng | Discontinuous Galerkin solver for Boltzmann-Poisson transients[END_REF][START_REF] Ayuso De Dios | Discontinuous Galerkin methods for the one-dimensional Vlasov-Poisson system[END_REF][START_REF] De Dios | High order and energy preserving discontinuous Galerkin methods for the Vlasov-Poisson system[END_REF][START_REF] Cheng | Energy-conserving discontinuous Galerkin methods for the Vlasov-Ampere system[END_REF][START_REF] Cheng | Discontinuous Galerkin methods for the Vlasov-Maxwell equations[END_REF]); moreover, one element only communicates with its immediate neighbors which is very important for parallelization capability [START_REF] Biswas | Parallel, adaptive finite element methods for conservation laws[END_REF] but also, thanks to this local data structure, the matrices used in these methods are sparse which is an important point if one wants to combine DG methods in space with exponential methods in time.

Then, after a finite differences approximation of the Vlasov equation, a DG method is employed for the space approximation using central fluxes. Indeed, this choice is motivated by the fact that the DG matrix has a pure imaginary spectrum which is not the case when monotone fluxes are considered. Moreover, using central fluxes makes the treatment more easier compared to upwind flux. Indeed, the latter requires to split the flux of the linear part into two parts according to the sign of velocity v, which prevents a discrete Poisson equation from being satisfied because of a lack consistency in the nonlinear part. Once the semi-discrete equations is obtained, an ODE system has to be solved in time. To do so, exponential time integrators are used to overcome the stringent condition coming from the linear part, but the exponential of a large matrix has to be computed. Thanks to a one side coupling between the distribution function and the electromagnetic fields, it is possible to compute explicitely the exponential of the matrix and to derive an efficient numerical scheme which is high order in time, space and velocity, preserves the total mass and a discrete Poisson equation. Some numerical illustrates the good behavior of the method. In particular, the expected order are recovered on a two-dimensional linear advection and a good agreement is obtained when we compare the DG exponential method to the Fourier exponential method.

The rest of the paper is organized as follows. In Section 2, we present the exponential DG method for one-dimensional (1D) linear transport problems. Section 3 is dedicated to the construction of the exponential DG method for Vlasov equations includes Vlasov-Ampère (1dx-1dv) and Vlasov Maxwell (1dx-2dv) equations. In Section 4, we discuss some extensions of the exponential method to high dimensional Vlasov-Maxwell equations. In Section 5, the capability of the proposed exponential DG method is illustrated through several numerical tests. Finally, after some concluding remarks, several appendices details some specific aspects of the method.

2 Exponential DG method for 1D transport equation

We firstly consider the 1D transport equation:

u t + au x = 0, x ∈ [x a , x b ],
u(0, x) = u 0 (x).

(

) 1
For simplicity, we assume periodic boundary conditions, and the velocity field a is a constant.

Here we take a as 1 for simplicity. We perform a partition of the computational domain

x a = x 1 2 < x 3 2 < • • • < x N + 1 2 = x b as the mesh partition. Let I j = [x j-1 2 , x j+ 1 2
] denote an element of length

∆x j = x j+ 1 2 -x j-1 2
and define ∆x = max j ∆x j . For simplicity, we consider the uniform mesh in this paper with ∆x j = ∆x = (x b -x a )/N . We define the finite dimensional approximation space, V k h = {v h : v h | I j ∈ P k (I j )}, where P k (I j ) denotes the set of polynomials of degree at most k on I j . For any ψ ∈ V k h , we also denote the left limit of ψ at cell boundary as ψ -and the right limit as ψ + . Multiply (1) by the test function ψ ∈ V k h , integrate on cell I j and integrate by parts, we end up with the semi-discrete DG scheme: find

u h ∈ V k h such that ˆIj (∂ t u h ψ)dx = -F | x j+ 1 2 ψ -| x j+ 1 2 + F | x j-1 2 ψ + | x j-1 2 + ˆIj F ψ x dx, j = 1, ..., N, (2)
where F (u) . = u and F is chosen as either a central or upwind fluxes central flux:

F | x j± 1 2 = u -+ u + 2 | x j± 1 2 , upwind flux: F | x j± 1 2 = u -| x j± 1 2 , (3)
Next, we consider ξ m j (m = 0, 1, . . . , k) a basis of P k (I j ) and we choose a modal basis defined as ξ m j (x) = ((x -x j )/∆x) m so that we have the representation u h (t, x)| I j = k m=0 u m j (t)ξ m j (x), with u m j (t) the degree of freedom. The semi-discrete DG scheme can eventually be written as an ordinary differential equation (ODE) satisfied by the DG degrees of freedom u m j (t) for m = 0, . . . , k and j = 1, . . . , N . Introducing the vector u

(t) ∈ R (k+1)N u(t) = (u 0 1 , u 1 1 , . . . , u k 1 , u 0 2 , u 1 2 , . . . , u k 2 , . . . , u 0 N , u 1 N , . . . , u k N ) T (t), (4)
the semi-discretized problem simply becomes

du dt = Au, with A ∈ M (k+1)N,(k+1)N (R). (5)
The 'DG-matrix' A contains the DG approximation [START_REF] Banks | A new class of nonlinear finite-volume methods for Vlasov simulation[END_REF] for which the details are given in Appendix A. From this semi-discrete in space formulation, a Runge-Kutta discretization is classically used to get high order accuracy in time [START_REF] Cockburn | Runge-Kutta discontinuous Galerkin methods for convectiondominated problems[END_REF][START_REF] Zhang | Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws[END_REF][START_REF] Zhang | Stability analysis and a priori error estimates of the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws[END_REF]. But it seems also natural to use an exponential method for time discretization, which turns out to be exact in this simple linear transport case. Denoting u n ≈ u(t n ) the fully discretized unknown (with t n = n∆t, n ∈ N and ∆t > 0 the time step), the exponential-DG scheme thus writes

u n+1 = exp(A∆t)u n , ∀n ∈ N, (6)
with u 0 = u 0 (u 0 being the degrees of freedom of the initial condition u 0 in (1)). Note that from u n ∈ R (k+1)N (whose components are denoted by (u m j ) n , j = 1, . . . , N and m = 0, . . . , k, following (4)), it is possible to reconstruct a piecewise polynomial function

u n h ∈ V k h from u n h (x) = k m=0 (u m j ) n ξ m j (x), ∀x ∈ I j , j = 1, . . . , N. (7)
The properties of such an approximation obviously depends on the structure of the DGmatrix A which is discussed now. From the calculations (given in Appendix A), the DG-matrix A obtained with the central flux and periodic boundary conditions enjoys a circulant tri-diagonal block structure so that it can be written as:

A = 1 ∆x      C 1 C 2 0 . . . C 3 C 3 C 1 C 2 0 . . . 0 . . . . . . . . . C 2 C 2 0 . . . C 3 C 1      , with C j = M -1 D j (j = 1, 2, 3) (8)
where the matrix elements of the matrices M, D j ∈ M (k+1),(k+1) (R) are given by (for ℓ, m = 1, . . . , k + 1, see Appendix A for details)

M ℓ,m = (1/2) m+ℓ-1 m + ℓ -1 [1 -(-1) m+ℓ-1 ], (D 2 ) ℓ,m = (-1) m (1/2) m+ℓ-1 , (D 1 ) ℓ,m = (1/2) m+ℓ-2 ℓ -1 m + ℓ -2 - 1 2 [1 -(-1) m+ℓ-2 ], (D 3 ) ℓ,m = (-1) ℓ-1 (1/2) m+ℓ-1 , (9)
with (D 1 ) 1,1 = 0 by convention. The choice of central fluxes implies the matrix A is diagonalizable and the eigenvalues are pure imaginary. This has been checked numerically and some discussions are performed in the following remarks.

Remark 1. In [START_REF] Tee | Eigenvectors of block circulant and alternating circulant matrices[END_REF], the author proposes a way to deduce the eigenvalues of A ∈ M (k+1)N,(k+1)N (R) from the eigenvalues of some matrices of size (k + 1), which can be computed explicitly for small values of k (numerically for larger k). Considering (ρ j ) j=0,...,N -1 the N -th roots of the unity (ρ N j = 1 for j = 0, . . . , N -1), the (k + 1)N eigenvalues of A given by (8) can be deduced from the (k + 1) eigenvalues of C j = C 1 + ρ j C 2 + ρ N -1 j C 3 for j = 0, . . . , N -1, where C j ∈ M (k+1),(k+1) . Then, we checked numerically that the eigenvalues of C j are pure imaginary for all j = 0, . . . , N -1, and we deduce from [START_REF] Tee | Eigenvectors of block circulant and alternating circulant matrices[END_REF] that it is also true for the eigenvalues of A.

Remark 2. We explore another way to check the eigenvalues of A given by (8) are pure imaginary by using symbolic software. Denoting P A (λ) the characteristic polynomial of A, we made the following observations • odd case: (k+1)N = 2d+1. In this case, we have P A (λ) = λ d ℓ=0 a 2ℓ λ 2ℓ with a 2ℓ ∈ R and the roots can be written as 0, λ j , λj , j = 1, . . . , d , in particular 0 is a simple eigenvalue in this case. Since P A (-λ) = P A (λ), we deduce Re(λ j ) = 0.

• even case: (k + 1)N = 2d. In this case, we have P (λ) = λ 2 d-1 ℓ=0 a 2ℓ λ 2ℓ with a 2ℓ ∈ R and the roots can be written as 0, λ j , λj , j = 1, . . . , d -1 , in particular 0 is a double eigenvalue in this case. Since P A (λ) = P A (-λ), we deduce Re(λ j ) = 0.

We can now study the stability of the numerical scheme. To do so, we write the following proposition.

Proposition 1. Let us consider the matrix A ∈ M (k+1)N,(k+1)N (R) given by (8)-( 9) There exists C > 0 such that, for any time t and any k, N , we have ∥ exp(At)∥ ≤ C, with ∥ • ∥ an induced matrix norm.

Proof. First, we write A = ∆x -1 A 1 from (8) with A 1 independent of ∆x (with ∆x = (x bx a )/N ) and since A is diagonalizable from Remarks 1 & 2, there exist P invertible and D diagonal such that A = ∆x -1 P DP -1 . Let us remark from ( 8)-( 9) that A 1 does not depend on the space mesh ∆x (and then dos not depend on N ), so does the matrix P . Thus, there exists C > 0 (independent of N ) such that cond(P ) ≡ ∥P ∥∥P -1 ∥ ≤ C, where ∥ • ∥ denotes an induced matrix norm. Now, since the eigenvalues of A are pure imaginary for all k, N , we have D j,j = iλ j , λ j ∈ R, j = 1, 2, . . . , (k + 1)N . Finally, we get for all t

∥ exp(At)∥ = ∥P exp(∆x -1 D)P -1 ∥ ≤ cond(P ) ∥ exp(∆x -1 D)∥ = cond(P ) ≤ C.
We end this section by proving an error estimate for exponential DG method.

Proposition 2. Let u(t, x) the exact solution of 1D transport problem (1) with a smooth initial condition u 0 and let u n h ∈ V k h the numerical solution, n = 0, . . . , N (with N = T /∆t, T being the final time and ∆t the time step) obtained from (6)- [START_REF] Caliari | A µ-mode integrator for solving evolution equations in kronecker form[END_REF] where A is the DG-matrix given by (8)-( 9). Then we have the following L 2 -norm error estimate:

||u(t n , •) -u n h || L 2 ≤ C∆x k ,
Proof. First, we introduce u h (t, x) the exact solution of the semi-discrete DG scheme [START_REF] Boutin | Modified lawson methods for Vlasov equations[END_REF]. The classical DG projection analysis gives ||u(t n , •) -u h (t n , •)|| L 2 ≤ C∆x k , where C is a positive constant independent on ∆x (see the details in [START_REF] Richter | An optimal-order error estimate for the discontinuous Galerkin method[END_REF][START_REF] Lesaint | On a finite element method for solving the neutron transport equation[END_REF][START_REF] Johnson | An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF][START_REF] Peterson | A note on the convergence of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF][START_REF] Zhang | Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws[END_REF][START_REF] Zhang | Stability analysis and a priori error estimates of the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws[END_REF], but a simple proof of stability and error estimate is given in Appendix B). Since the exponential method exactly solves the semi-discrete DG scheme, we have ||u h (t n , •) -u n h || L 2 = 0. Finally we have

||u(t n , •) -u n h || L 2 ≤ ||u(t n , •) -u h (t n , •)|| L 2 + ||u h (t n , •) -u n h || L 2 ≤ C∆x k .
Remark 3. It is worthy to be mentioned that the (k + 1)th order optimal convergence rate have been proved for DG with monotone flux (see [START_REF] Johnson | An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF][START_REF] Peterson | A note on the convergence of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF][START_REF] Zhang | Error estimates to smooth solutions of Runge-Kutta discontinuous Galerkin methods for scalar conservation laws[END_REF][START_REF] Zhang | Stability analysis and a priori error estimates of the third order explicit Runge-Kutta discontinuous Galerkin method for scalar conservation laws[END_REF]). However, kth order suboptimal convergence rate is proved for DG with central flux in [START_REF] Liu | Sub-optimal convergence of discontinuous Galerkin methods with central fluxes for linear hyperbolic equations with even degree polynomial approximations[END_REF] and a discussion is performed according to the oddness of k.

3 Semi-discretization of some Vlasov models with Discontinuous Galerkin method

In this section, we consider the numerical approximation of Vlasov-Maxwell equations using DG framework in space (as presented in the previous section) and finite differences in the velocity direction. The semi-discretization (in both space and velocity) is presented and we will see that the so-obtained ODE system is amenable to Lawson time integrators. We first present the methodology on the 1d x -1d v case on the Vlasov-Ampère system and then we consider the 1d x -2d v case on the Vlasov-Maxwell system.

Vlasov-Ampère equation

The equation we address is the following Vlasov-Ampère model satisfied by the distribution function f (t, x, v) ≥ 0 and the electric field

E(t, x) ∈ R with t ≥ 0, x ∈ [0, L] (L > 0) and v ∈ R,        ∂f ∂t + v ∂f ∂x + E ∂f ∂v = 0, ∂ t E = -ˆR vf dv + J, J = 1 L ˆL 0 ˆR vf dxdv, (10)
with the initial conditions (f 0 (x, v), E 0 (x)) such that the Poisson equation is satisfied initially ∂ x E 0 = ´R f 0 dv -ρ, with ρ = (1/L) ´L 0 ´R f 0 dvdx and periodic boundary conditions are imposed in space. The Vlasov-Ampère system is equivalent to the Vlasov-Poisson model where the electric field satisfies the Poisson equation ∂ x E = ´R f dv -ρ.

Semi-discretization

We shall use a DG method in the space direction x as presented in the previous section and we consider a truncated domain [-v max , v max ] in the velocity direction discretized by v j = -v max + j∆v, j = 0, . . . , N v , ∆v = 2v max /N v being the velocity mesh step. We firstly present DG discretization for E in the x-direction with

I i = [x i-1 2 , x i+ 1 2
], i = 1, ..., N x (N x being the number of cells):

E(t, x) ≈ E h (t, x) = Nx i=1 E h (t, x)| I i = Nx i=1 k m=0 E m i (t)ξ m (x).
Then, considering finite difference method for f in v direction, we consider the DG approximation in x-direction through

f (t, x, v j ) ≈ f h (t, x, v j ) = Nx i=1 k m=0 f m i (t, v j )ξ m (x), j = 0, ..., N v .
As in the previous section, we denote f j (t) ∈ R (k+1)Nx the vector of the DG coefficients f m i (t, v j ) of f h (t, x, v j ) using DG in space and evaluated at the velocity grid v j whereas E ∈ R (k+1)Nx denotes the vector of DG coefficients E m i of E h (t, x). For Vlasov-Ampère equation [START_REF] Cheng | Energy-conserving discontinuous Galerkin methods for the Vlasov-Ampere system[END_REF], we have the following DG scheme with the DG representation of E and f :

k m=0 ∂ t f m i (t, v j )ξ m , ξ ℓ I i -v j f m i (t, v j )ξ m , ∂ x ξ ℓ (x) I i + v j {f h (t, x, v j )}ξ ℓ i+1/2 i-1/2 + k m=0 k n=0 E n i ξ n (Df m i )(t, v j )ξ m , ξ ℓ I i = 0, (11)
where we used the central flux {f h (t, x, v j )}|

x i±1/2 = 1 2 (f h (t, x + i±1/2 , v j ) + f h (t, x - i±1/2 , v j )), ℓ = 0, 1, 2, ..., k, i = 1, 2, ...N x and Df (v j ) denotes a discrete approximation of (∂ v f )(v j ) (an example would be (Df )(v j ) = f (v j+1 )-f (v j-1 ) 2∆v
but any higher order finite difference approximation can be used). We denote

f j,i (t) = (f 0 i (t, v j ), f 1 i (t, v j ), ...., f k i (t, v j )) T , i = 1, 2, ..., N x , E i (t) = (E 0 i (t), E 1 i (t), ...., E k i (t)) T , i = 1, 2, ..., N x , and Df j,i (t) = ((Df 0 i )(t, v j ), (Df 1 i )(t, v j ), ...., (Df k i )(t, v j )) T , i = 1, 2, ..., N x .
We can rewrite the DG discretization as an ODE system of size (k + 1)N x for each j = 1, . . . , N v ∆x

     M M . . . M      d dt      f j,1 f j,2 . . . f j,Nx      -v j      D 1 D 2 . . . D 3 D 3 D 1 D 2 . . . . . . D 2 D 2 . . . D 3 D 1           f j,1 f j,2 . . . f j,Nx      +      B 1 B 2 . . . B Nx           (Df ) j,1 (Df ) j,2 . . . (Df ) j,Nx      = 0, (12)
where the matrices M, D 1 , D 2 , D 3 are the same as in the previous section (see also Appendix A) and B i (i = 1, . . . , N x ) are matrices of size k + 1 with elements (B i ) ℓ,m = ( k n=0 E n i ξ n ξ m , ξ ℓ ) I i . Introducing now the following vector containing the degrees of freedom of f h and E h

f j (t) = (f 0 1 (t, v j ), . . . , f k 1 (t, v j ), f 0 2 (t, v j ), . . . , f k 2 (t, v j ), . . . , f 0 Nx (t, v j ), . . . , f k Nx (t, v j )) T , E(t) = (E 0 1 (t), . . . , E k 1 (t), E 0 2 (t), . . . , E k 2 (t), . . . , E 0 Nx (t), . . . , E k Nx (t)) T ,
and Df j is defined similarly as f j , we can rewrite the DG scheme [START_REF] Cheng | Discontinuous Galerkin solver for Boltzmann-Poisson transients[END_REF] as

∂ t f j = v j Af j -Ẽ(Df ) j , (13)
A ∈ M (k+1)Nx,(k+1)Nx is the DG-matrix [START_REF] Califano | Spatial structure and time evolution of the weibel instability in collisionless inhomogeneous plasmas[END_REF] and Ẽ ∈ M (k+1)Nx,(k+1)Nx is a block diagonal matrix composed of N x block matrices of size (k + 1) × (k + 1) defined by (∆xM ) -1 B i , i = 1, . . . , N x .

Let us remark that we consider central finite differences method to approximate (Df ) j to avoid to discuss the sign of matrix Ẽ ∈ M (k+1)Nx,(k+1)Nx compared with upwind FD method. Let us now discuss the discretization of the Ampère equation. Since our goal is to find a consistent discretization that is compatible with a discrete Poisson equation, we first discuss how to solve the Poisson equation. Using the above discretization, we will use the DG matrix A which is an approximation of (-∂ x ). Thus, a direct approximation of the initial Poisson equation ∂ x E(0, x) = ´R f 0 (x, v)dv -ρ would be -AE 0 = j f 0 j ∆v -ρ (with E 0 and f 0 j the degrees of freedom of E(0, x) and f 0 (x, v j ). However, as mentioned in Remark 2, A ∈ M (k+1)Nx,(k+1)Nx is not invertible and we then introduce Π ∈ M (k+1)Nx,(k+1)Nx the projection onto the Ker(A) so that (A+Π) is invertible on R(A) with R(A) the range of A. Here we impose condition ΠE 0 = 0 to preserve the uniqueness of the solution E ∈ R (k+1)Nx , which is similarly as the constraint ´R E(x, v)dx = 0 for Poisson equation itself. We then consider the following discretized Poisson equation

-(A + Π)E 0 = (1 -Π) j f 0 j ∆v -ρ = (1 -Π) j f 0 j ∆v, (14)
with 1 the identity matrix of size (k + 1)N x and where in the last equality, we used the fact that constants belong to Ker(A) (see Appendix C for details).

We deduce the discretization of Ampère equation from the time derivative of the discretized Poisson equation inspired from [START_REF] Cockburn | Locally divergence-free discontinuous Galerkin methods for the Maxwell equations[END_REF], that is:

-(A + Π)E(t) = (1 -Π) j f j (t)∆v.
Considering the time derivative of the latter equation and using (13) leads to

-(A + Π)∂ t E(t) = (1 -Π)∂ t j f j (t)∆v = (1 -Π) j v j Af j (t)∆v = A(1 -Π) j v j f j (t)∆v = (A + Π)(1 -Π) j v j f j (t)∆v,
where we used j Df j = 0 and some relations between A and Π. Hence, we consider the following DG discretization of the Ampère equation

∂ t E(t) = -(1 -Π) j v j f j (t)∆v. (15)
Finally, gathering ( 13) and ( 15) enables to get the following semi-discretized scheme for the Vlasov-Ampère system

     ∂ t f j = v j Af j -Ẽ(Df ) j , ∂ t E = - j v j f j ∆v + j v j Πf j ∆v. ( 16
)
In view of the time discretization, we introduce the following vector of semi-discrete unknown NxNv . Then the previous system ( 16) can be rewritten as

U = ( ⃗ f , E) ∈ R (k+1)Nx(Nv+1) with ⃗ f = (f 1 , f 2 , . . . , f Nv ) ∈ R (k+1)
∂ t U = LU + N (U ), (17)
with L ∈ M (k+1)Nx(Nv+1),(k+1)Nx(Nv+1) given by

L =         v 1 A 0 Ñ , Ñ 0 Ñ , Ñ . . . 0 Ñ , Ñ 0 Ñ , Ñ 0 Ñ , Ñ v 2 A 0 Ñ , Ñ . . . 0 Ñ , Ñ . . . . . . . . . . . . . . . . . . 0 Ñ , Ñ . . . 0 Ñ , Ñ v Nv A 0 Ñ , Ñ -∆vv 1 (1 -Π) -∆vv 2 (1 -Π) . . . -∆vv Nv (1 -Π) 0 Ñ , Ñ         (18)
where we denote Ñ = (k + 1)N x and 1 = 1 Ñ , Ñ the identity matrix of size Ñ × Ñ . Finally, Nv+1) are given by

U, N (U ) ∈ R (k+1)Nx(
U = (f 1 , f 2 , . . . , f Nv , E) T , N (U ) = -Ẽ(Df ) 1 , -Ẽ(Df ) 2 , . . . , -Ẽ(Df ) Nv , 0 T .

Time discretisation

The goal of this part is to present time discretization of ( 17) to get a fully discretized scheme of the Vlasov-Ampère system [START_REF] Cheng | Energy-conserving discontinuous Galerkin methods for the Vlasov-Ampere system[END_REF]. The form ( 17) is amenable to exponential scheme [START_REF] Hochbruck | Exponential integrators[END_REF][START_REF] Crouseilles | Exponential methods for solving hyperbolic problems with application to collisionless kinetic equations[END_REF] which is motivated by the fact that the linear part acts on a different scale compared to the nonlinear part in Vlasov type problems. Moreover, as discussed in Section 2, the linear part can be computed exactly thanks to the exponential. Among the exponential schemes, we shall use the Lawson class of methods for stability reasons [START_REF] Crouseilles | Exponential methods for solving hyperbolic problems with application to collisionless kinetic equations[END_REF].

Denoting U n = ( ⃗ f n , E n ) ≈ ( ⃗ f (t n ), E(t n )) = U (t n
) with t n = n∆t, ∀n ∈ N (∆t > 0 being the time step), the simplest (first order in time) Lawson method can be written as

U n+1 = exp(∆tL)U n + ∆t exp(∆tL)N (U n ). (19)
High order methods can be obtained from Runge-Kutta methods using the corresponding Butcher tableau [START_REF] Crouseilles | Exponential methods for solving hyperbolic problems with application to collisionless kinetic equations[END_REF].

The key point of exponential methods lies in the computation of exp(∆tL) with L given by [START_REF] Crouseilles | Exponential methods for solving hyperbolic problems with application to collisionless kinetic equations[END_REF]. To compute exp(At) ∀t > 0, one considers the linear part only ∂ t U = LU . First, we observe that the distribution function f j part is decoupled from the electric field part, so that it can be solved directly and we have, similarly as in Section 2

f j (t) = exp(v j At)f j (0). (20)
Now, let us consider the electric field equation

∂ t E = - j v j (1 -Π)f j (t)∆v,
which gives, after integrating it in time

E(t) = E(0) - j ˆt 0 v j (1 -Π)f j (s)∆vds.
Now, replacing f j (s) by exp(v j As)f j (0) from ( 20) enables to get an explicit expression of E(t). Indeed, using properties used to derive [START_REF] Cockburn | Locally divergence-free discontinuous Galerkin methods for the Maxwell equations[END_REF] we have

E(t) = E(0) - j ˆt 0 v j (1 -Π)f j (s)∆vds = E(0) - j ˆt 0 v j (1 -Π) exp(v j As)f j (0)ds∆v = E(0) - j ˆt 0 v j (1 -Π) ∞ k=0 (v j As) k k! Πf j (0) + (1 -Π)f j (0) ds∆v = E(0) - j v j (1 -Π) ˆt 0 Πf j (0) + ∞ k=0 (v j (A + Π)s) k k! (1 -Π)f j (0) ds∆v = E(0) - j v j (1 -Π) tΠf j (0) + ∞ k=0 (v j (A + Π)) k t k+1 (k + 1)! (1 -Π)f j (0) ∆v = E(0) - j (1 -Π) ∞ k=0 (A + Π) -1 (v j (A + Π)t) k+1 (k + 1)! (1 -Π)f j (0) ∆v = E(0) - j (1 -Π)(A + Π) -1 (exp(v j At) -1)(1 -Π)f j (0)∆v = E(0) - j (A + Π) -1 (exp(v j At) -1)(1 -Π)f j (0)∆v = E(0) + ∆v j (A + Π) -1 (1 -exp(v j At))f j (0).
Hence, denoting à = A + Π, we deduce from the above calculation the expression of exp(Lt):

        e v 1 At 0 Ñ , Ñ 0 Ñ , Ñ . . . 0 Ñ , Ñ 0 Ñ , Ñ 0 Ñ , Ñ e v 2 At 0 Ñ , Ñ . . . 0 Ñ , Ñ . . . . . . . . . . . . . . . . . . 0 Ñ , Ñ . . . 0 Ñ , Ñ e v Nv At 0 Ñ , Ñ ∆v Ã-1 (1 -e v 1 At ) ∆v Ã-1 (1 -e v 2 A∆t ) . . . ∆v Ã-1 (1 -e v Nv At ) 1         . (21)
Hence, the exponential-DG scheme for the Vlasov-Ampère equation corresponds to ( 19)- [START_REF] Filbet | Comparison of Eulerian Vlasov solvers[END_REF].

For this scheme, one can prove in the following proposition that a discrete Poisson equation is satisfied for each iteration.

Proposition 3. The exponential DG method (19)-( 21) (and its generalization to high order Lawson Runge-Kutta) satisfied by U n = ( ⃗ f n , E n ) preserves the following discretized Poisson equation:

(A + Π)E n = - j (1 -Π)f n j ∆v, ∀n ∈ N ⋆ ,
provided that it is satisfied at the initial time n = 0. Here, A is the DG matrix given by ( 8)-( 9), Π is the orthogonal projection onto Ker(A) and 1 is the identity matrix of size (k + 1)N x .

Proof. We present the proof for first order Lawson case (forward Euler), the proof can be generalized to arbitrary explicit Runge-Kutta scheme. First, we assume the Poisson equation

(A + Π)E 0 = -j (1 -Π)f 0
j ∆v holds at the initial time. Next, from the scheme [START_REF] De Dios | High order and energy preserving discontinuous Galerkin methods for the Vlasov-Poisson system[END_REF] with exp(Lt) given by ( 21), we have

f n+1 j = exp(v j A∆t)f n j -∆t exp(v j A∆t)( ẼDf ) n j , (22)
whereas for the E component, we have

E n+1 = E n + ∆v j (A + Π) -1 (1 -exp(v j A∆t))f n j -∆t∆v j (A + Π) -1 (1 -exp(v j A∆t))( ẼDf ) n j = E n + ∆v j (A + Π) -1 (1 -Π)(1 -exp(v j A∆t))f n j -∆t∆v j (A + Π) -1 (1 -Π)(1 -exp(v j A∆t))( ẼDf ) n j .
The last term can be split into two parts: the first one -∆t∆v(A + Π) -1 (1 -Π) j ( ẼDf ) n j vanishes thanks to the conservative properties of the discrete operator D whereas the second one, we use [START_REF] Hochbruck | Exponential integrators[END_REF] to get

∆t∆v(A + Π) -1 j (1 -Π) exp(v j A∆t)( ẼDf ) n j = ∆v(A + Π) -1 j (1 -Π)(exp(v j A∆t)f n j -f n+1 j ) = ∆v(A + Π) -1 j (1 -Π) exp(v j A∆t)f n j -(1 -Π)f n+1 j .
Finally, we then get

E n+1 = E n + ∆v j (A + Π) -1 (1 -Π)(1 -exp(v j A∆t))f n j +∆v(A + Π) -1 j (1 -Π) exp(v j A∆t)f n j -(1 -Π)f n+1 j = E n + ∆v j (A + Π) -1 ((1 -Π)f n j -(1 -Π)f n+1 j ). (23)
By induction, if the discrete Poisson equation is satisfied at iteration n, then it is satisfied at iteration n + 1 and the proof is complete.

Remark 4. In this remark, we discuss how to compute the matrix Π. As mentioned in Remark 2, when (k + 1)N x is odd (referred as odd case), 0 is a single eigenvalue and the associated eigenvector u 1 corresponds to the constant function in the DG space (see Appendix C). When (k + 1)N x is even, 0 has a double multiplicity and one has to find another eigenvector u 2 (see Appendix C). Once we get the eigenvectors u 1 , u 2 ∈ R (k+1)Nx of A associated to the eigenvalue 0, then Π is given by Πx = ⟨x,

u 1 ⟩u 1 + ⟨x, u 2 ⟩u 2 = (u T 1 ⊗ u 1 + u T 2 ⊗ u 2 )
x, for all x ∈ R (k+1)Nx , with the Kronecker product ⊗ and (v T ⊗ v) i,j = v i v j . Some examples are given in Appendix C.

Vlasov Maxwell equations 1dx-2dv

In this part, We consider the following Vlasov-Maxwell 1dx-2dv model satisfied by

f (t, x, v 1 , v 2 ), E 1 (t, x), E 2 (t, x), B(t, x), with t ≥ 0, x ∈ [0, L](L > 0) and (v 1 , v 2 ) ∈ R 2          ∂ t f + v 1 ∂ x f + E 1 ∂ v 1 f + E 2 ∂ v 2 f + B(v 2 ∂ v 1 f -v 1 ∂ v 2 f ) = 0, ∂ t B = -∂ x E 2 , ∂ t E 1 = - ˆR2 v 1 f dv 1 dv 2 + J1 , ∂ t E 2 = -∂ x B - ˆR2 v 2 f dv 1 dv 2 + J2 , (24)
with the initial conditions

(f 0 (x, v), E 0 1 (x), E 0 2 (x), B 0 (x)) such that the Poisson equation is satis- fied initially ∂ x E 0 1 (x) = ´R2 f 0 (x, v)dv 1 dv 2 -ρ, with ρ = 1 L ´L 0 ´R2 f 0 (x, v
)dxdv 1 dv 2 and periodic boundary conditions are imposed in space. Here Ji = 1 L ´L 0 ´R2 v i f dxdv 1 dv 2 ensures that the electric fields are zero average in space.

Semi-discretization

We follow the lines of the above subsection and use a DG method in the space direction x and we consider a grid in the velocity direction v ℓ,

j ℓ = -v ℓ,max + j ℓ ∆v ℓ , v ℓ ∈ [-v ℓ,max , v ℓ,max ], ℓ = 1, 2, ∆v ℓ = 2v ℓ,max /N v ℓ , N v ℓ ∈ N ⋆ .
The definitions of the different objects are a direct extension of the definitions introduced in the previous part. Indeed, we denote by f j 1 ,j 2 ∈ R (k+1)Nx the DG coefficient vector of

f j 1 ,j 2 (t, x) = f h (t, x, v j 1 , v j 2 ) ≈ f (t, x, v j 1 , v j 2 )
in space and evaluated at the velocity grid, and

E 1 , E 2 , B ∈ R (k+1)Nx are the DG coefficient vectors of (E 1,h , E 2,h , B h )(t, x) ≈ (E 1 , E 2 , B)(t, x). Moreover, ( F Df ) j 1 ,j 2 (with Fj 1 ,j 2 ∈ M (k+1)Nx,(k+1)Nx (R)) is obtained as previ- ously by a DG approximation of the nonlinear term (F • ∇ v f ), with F = (E 1 + Bv 2 , E 2 -Bv 1 ).
For fixed indices j 1 , j 2 , the derivation of the numerical scheme is very similar to the 1dx-1dv case and we then obtain the following semi-discretized (in space and velocity) scheme for

j ℓ = 1, . . . , N v ℓ (ℓ = 1, 2)          ∂ t f j 1 ,j 2 = v 1,j 1 Af j 1 ,j 2 -( F Df ) j 1 ,j 2 , ∂ t B = AE 2 , ∂ t E 1 = - j 1 ,j 2 v 1,j 1 (1 -Π)f j 1 ,j 2 ∆v 1 ∆v 2 , ∂ t E 2 = AB - j 1 ,j 2 v 2,j 2 (1 -Π)f j 1 ,j 2 ∆v 1 ∆v 2 , ( 25
)
where A is given by ( 8)-( 9) and 1 the identity matrix of size (k + 1)N x and Π the projection matrix onto Ker(A) introduced in the previous part. Let us denote

U (t) = ( ⃗ f (t), B(t), E 1 (t), E 2 (t)), where ⃗ f (t) ∈ R (k+1)NxNv 1 Nv 2 contains the DG coefficients f j 1 ,j 2 , B(t), E 1 (t), E 2 (t) ∈ R (k+1)Nx
denote the DG coefficients of the electromagnetic fields. Using the above notations and the following ones

⃗ v ℓ ∈ R Nv ℓ , (⃗ v ℓ ) j = v ℓ,j = -v ℓ,max + j∆v ℓ , j = 1, . . . , N v ℓ and ℓ = 1, 2, ⃗ f ⋆,j 2 ∈ R (k+1)NxNv 1 for j 2 = 1, . . . , N v 2 .
We also introduce the following compact notations for the size of the matrices: Ñ = (k + 1)N x and Ñ1 = (k + 1)N x N v 1 = Ñ N v 1 . The system (25) can be recast as

∂ t U = LU + N (U ), (26)
with

U = ⃗ f ⋆,1 , ⃗ f ⋆,2 , . . . , ⃗ f ⋆,Nv 2 , B, E 2 , E 1 T , ⃗ f ⋆,j 2 ∈ R Ñ1 , ∀j 2 = 1, . . . , N v 2 N (U ) = -( F D ⃗ f ) ⋆,1 , -( F D ⃗ f ) ⋆,2 , . . . , -( F D ⃗ f ) ⋆,Nv 2 , 0, 0, 0 T , (27)
L =              diag( ⃗ v 1 )⊗A 0 Ñ1 , Ñ1 0 Ñ1 , Ñ1 . . . 0 Ñ1 , Ñ1 0 Ñ1 , Ñ 0 Ñ1 , Ñ 0 Ñ1 , Ñ 0 Ñ1 , Ñ1 diag( ⃗ v 1 )⊗A 0 Ñ1 , Ñ1 . . . 0 Ñ1 , Ñ1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 0 Ñ1 , Ñ1 . . . 0 Ñ1 , Ñ1 diag( ⃗ v 1 )⊗A 0 Ñ1 , Ñ 0 Ñ1 , Ñ 0 Ñ1 , Ñ 0 Ñ , Ñ1 . . . 0 Ñ , Ñ1 0 Ñ , Ñ1 0 Ñ , Ñ A 0 Ñ , Ñ E 2,1 E 2,2 . . . E 2,Nv 2 A 0 Ñ , Ñ 0 Ñ , Ñ E 1 E 1 . . . E 1 0 Ñ , Ñ 0 Ñ , Ñ 0 Ñ , Ñ              (28)
where diag

( ⃗ v 1 ) ∈ M Nv 1 ,Nv 1 (R) denotes the diagonal matrix with ⃗ v 1 ∈ R Nv 1 on its diagonal, ⊗ denotes Kronecker product, diag( ⃗ v 1 )⊗A ∈ M Ñ1 , Ñ1
and 0 m,n is the zero matrix with m lines and n columns. Moreover, the matrices

E 1 , E 2,j 2 ∈ M Ñ , Ñ1 (R) for j 2 = 1, . . . , N v 2 are defined by E 1 = -∆v 1 ∆v 2 ⃗ v 1 ⊗ (1 -Π), E 2,j 2 = -∆v 1 ∆v 2 v 2,j 2 1 ⊗ (1 -Π) with 1 = (1, . . . , 1) ∈ R Nv 1 .
The size of the matrix L is ( Ñ1

N v 2 + 3 Ñ ) × ( Ñ1 N v 2 + 3 Ñ ) = (3 + N v 1 N v 2 )(k + 1)N x × (3 + N v 1 N v 2 )(k + 1)N x .
Even if it is a large matrix, one can see that L is sparse which will help to compute its exponential.

Time discretization

We now study the time discretization of ( 26) and as previously, we will design an exponential scheme. To do so, the discrete unknown

U n ≈ U (t n ) (t n = n∆t, ∆t > 0) is updated by U n+1 = exp(∆tL)U n + ∆t exp(∆tL)N (U n ), (29)
and one has to compute exp(∆tL). We give the following proposition to show the representation of exp(∆tL).

Proposition 4. The exponential exp(L∆t) of the matrix L given by ( 28) is given by

              e ∆tdiag( ⃗ v 1 )⊗A 0 Ñ1 , Ñ1 0 Ñ1 , Ñ1 . . . 0 Ñ1 , Ñ1 0 Ñ1 , Ñ 0 Ñ1 , Ñ 0 Ñ1 , Ñ 0 Ñ1 , Ñ1 e ∆tdiag( ⃗ v 1 )⊗A 0 Ñ1 , Ñ1 . . . 0 Ñ1
0 Ñ , Ñ e E 1 e E 1 . . . e E 1 0 Ñ , Ñ 0 Ñ , Ñ 1               where e B j 2 , e E 2,j 2 , e E 1 ∈ M Ñ , Ñ1 (R) for j 2 = 1, . . . , N v 2 are given by e B j 2 = ∆v 1 ∆v 2 v 2,j 2 (A + Π) -1 ⃗ α, e E 2,j 2 = ∆v 1 ∆v 2 v 2,j 2 (A + Π) -1 ⃗ β, e E 1 = ∆v 1 ∆v 2 (A + Π) -1 (1 ⊗ 1 -e ∆tdiag( ⃗ v 1 )⊗A )
where the matrices ⃗ α, ⃗ β are given by

⃗ α = [α 1 , α 2 , . . . , α Nv 1 ] ∈ M Ñ , Ñ1 (R) with α j 1 = -e A∆t 2(1 -v 1,j 1 ) - e -A∆t 2(1 + v 1,j 1 ) + e A∆tv 1,j 1 (1 -v 2 1,j 1 ) , ⃗ β = [β 1 , β 2 , . . . , β Nv 1 ] ∈ M Ñ , Ñ1 (R) with β j 1 = -e A∆t 2(1 -v 1,j 1 ) + e -A∆t 2(1 + v 1,j 1 ) + v 1,j 1 e A∆tv 1,j 1 (1 -v 2 1,j 1 )
.

Proof. First the ( Ñ1 N v 2 ) × ( Ñ1 N v 2 ) block is diagonal and the diagonal part is e ∆tdiag( ⃗ v 1 )⊗A
(N v 2 times). Second, the 3 × 3 right bottom block corresponds to the homogeneous Maxwell equations. Its exponential can be computed and is equal to

exp    0 Ñ , Ñ A∆t 0 Ñ , Ñ A∆t 0 Ñ , Ñ 0 Ñ , Ñ 0 Ñ , Ñ 0 Ñ , Ñ 0 Ñ , Ñ    =    exp(A∆t)+exp(-A∆t) 2 exp(A∆t)-exp(-A∆t) 2 0 Ñ , Ñ exp(A∆t)-exp(-A∆t) 2 exp(A∆t)+exp(-A∆t) 2 0 Ñ , Ñ 0 Ñ , Ñ 0 Ñ , Ñ 1    .
Finally, we compute the three last block lines of exp(L∆t).

Computation of e E 1 : solve E 1 First, we have for f j 1 ,j 2 (t):

f j 1 ,j 2 (t) = e v 1,j 1 A(t-t n ) f j 1 ,j 2 (t n ), which enables to compute E 1 (t n+1 ) E 1 (t n+1 ) = E 1 (t n ) - j 1 ,j 2 ˆtn+1 t n (1 -Π)e v 1,j 1 A(t-t n ) dtv 1,j 1 f j 1 ,j 2 (t n )∆v 1 ∆v 2 = E 1 (t n ) + ∆v 1 ∆v 2 (A + Π) -1 j 1 ,j 2 (1 -e v 1,j 1 A∆t )f j 1 ,j 2 (t n ),
from which we can thus deduce the last line of the exponential of the matrix.

Computation of e B j 2 , e E 2,j 2 : solve B, E 2 Next, we focus on the calculation of E 2 (t n+1 ) and B(t n+1 ) from known initial conditions E 2 (t n ) and B(t n ). Let us write down the equations for

(E 2 , B)(t) with t ∈ [t n , t n+1 ] d dt E 2 (t) = AB(t) - j 1 ,j 2 (1 -Π)e v 1,j 1 A(t-t n ) v 2,j 1 f j 1 ,j 2 (t n )∆v 1 ∆v 2 d dt B(t) = AE 2 (t)
which can be rewritten as dU dt = M U + R with U (t) = (E 2 (t), B(t)) and

M = 0 Ñ , Ñ A A 0 Ñ , Ñ , R(t) = R 1 (t) 0 Ñ ,1 , R 1 (t) = - j 1 ,j 2 (1 -Π)e v 1,j 1 A(t-t n ) v 2,j 1 f j 1 ,j 2 (t n )∆v 1 ∆v 2 (30)
Thus, one can write the variation of constant formula

U (t n+1 ) = e M ∆t U (t n ) + ˆtn+1 t n e -M (t-t n+1 ) R(t)dt. (31)
First, e M ∆t reads as, using its definition [START_REF] Peterson | A note on the convergence of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF] e M ∆t = 1 2

e A∆t + e -A∆t e A∆t -e -A∆t e A∆t -e -A∆t e A∆t + e -A∆t .

Second, one has to compute the integral term in (31)

ˆtn+1 t n e -M (t-t n+1 ) R(t)dt = 1 2   ´tn+1 t n e -A(t-t n+1 ) + e A(t-t n+1 ) R 1 (t)dt ´tn+1 t n e -A(t-t n+1 ) -e A(t-t n+1 ) R 1 (t)dt   = - 1 2   ´tn+1 t n j 1 ,j 2 e -A(t-t n+1 ) + e A(t-t n+1 ) e v 1,j 1 A(t-t n ) v 2,j 1 f n k,j 1 ,j 2 ∆v 1 ∆v 2 dt ´tn+1 t n j 1 ,j 2 e -A(t-t n+1 ) -e A(t-t n+1 ) e v 1,j 1 A(t-t n ) v 2,j 1 f n k,j 1 ,j 2 ∆v 1 ∆v 2 dt   = - j 1 ,j 2 I 1 v 2,j 1 f n k,j 1 ,j 2 ∆v 1 ∆v 2 j 1 ,j 2 I 2 v 2,j 1 f n k,j 1 ,j 2 ∆v 1 ∆v 2
where I 1 , I 2 are given by

I 1 = 1 2 ˆtn+1 t n e -A(t-t n+1 ) + e A(t-t n+1 ) e v 1,j 1 A(t-t n ) dt = (A + Π) -1 e A∆t 2(1 -v 1,j 1 ) - (A + Π) -1 e -A∆t 2(1 + v 1,j 1 ) - (A + Π) -1 v 1,j 1 e A∆tv 1,j 1 (1 -v 2 1,j 1 )
,

I 2 = 1 2 ˆtn+1 t n e -A(t-t n+1 ) -e A(t-t n+1 ) e v 1,j 1 A(t-t n ) dt = (A + Π) -1 e A∆t 2(1 -v 1,j 1 ) + (A + Π) -1 e -A∆t 2(1 + v 1,j 1 ) - (A + Π) -1 e A∆tv 1,j 1 (1 -v 2 1,j 1 )
.

Inserting these calculations in (31) leads to the following expression for E 2 (t n+1 ) and B(t n+1 )

E 2 (t n+1 ) = 1 2 (e A∆t + e -A∆t )E 2 (t n ) + 1 2 (e A∆t -e -A∆t )B(t n ) + ∆v 1 ∆v 2 (A + Π) -1 j 1 ,j 2 v 2,j 2 β j 1 f j 1 ,j 2 (t n ) B(t n+1 ) = 1 2 (e A∆t -e -A∆t )E 2 (t n ) + 1 2 (e A∆t + e -A∆t )B(t n ) + ∆v 1 ∆v 2 (A + Π) -1 j 1 ,j 2 v 2,j 2 α j 1 f j 1 ,j 2 (t n ) where ⃗ β = [β 1 , β 2 , . . . , β Nv 1 ] ∈ M Ñ × Ñ1 (R) and ⃗ α = [α 1 , α 2 , . . . , α Nv 1 ] ∈ M Ñ × Ñ1 (
R) are given by

β j 1 = -e A∆t 2(1 -v 1,j 1 ) + e -A∆t 2(1 + v 1,j 1 ) + v 1,j 1 e A∆tv 1,j 1 (1 -v 2 1,j 1 ) α j 1 = -e A∆t 2(1 -v 1,j 1 ) - e -A∆t 2(1 + v 1,j 1 ) + e A∆tv 1,j 1 (1 -v 2 1,j 1 )
, which conclude the proof.

Shared the same spirit with Vlasov-Ampère equation, we have the following discretized Poisson equation and error estimate for Vlasov-Maxwell equation with exponential Lawson RK DG FD discretization.

Proposition 5. The exponential DG method [START_REF] Muñoz-Matute | Exploiting the Kronecker product structure of φ-functions in exponential integrators[END_REF] where the exponential is given in Prop 4 (and its generalization to high order Lawson Runge-Kutta) satisfied by U n = ( ⃗ f , B, E 2 , E 1 ) n preserves the following Poisson equation

(A + Π)E n 1 = - j 1 ,j 2 (1 -Π)f n j 1 ,j 2 ∆v 1 ∆v 2 , ∀n ∈ N ⋆ ,
provided that it is satisfied at the initial time n = 0. Here, A is the DG matrix given by ( 8)-( 9), Π the orthogonal projection onto Ker(A) and 1 is the identity matrix of size Ñ = (k + 1)N x .

Proof. We present the proof for first order Lawson case (forward Euler), the proof can be generalized to arbitrary explicit Runge-Kutta scheme. First, we assume the Poisson equation

(A + Π)E 0 1 = -j 1 ,j 2 (1 -Π)f 0 j 1 ,j 2
∆v 1 ∆v 2 holds at the initial time. Next, from the scheme [START_REF] Muñoz-Matute | Exploiting the Kronecker product structure of φ-functions in exponential integrators[END_REF] and Prop 4, we have

f n+1 j 1 ,j 2 = e v 1,j 1 A∆t f n j 1 ,j 2 -∆te v 1,j 1 A∆t ( F Df ) n j 1 ,j 2 .
Regarding the E 1 component, we have

E n+1 1 = E n 1 + ∆v 1 ∆v 2 (A + Π) -1 j 1 ,j 2 (1 -e v 1,j 1 A∆t )f n j 1 ,j 2 -∆t∆v 1 ∆v 2 (A + Π) -1 j 1 ,j 2 (1 -e v 1,j 1 A∆t )( F Df ) n j 1 ,j 2 = E n 1 + ∆v 1 ∆v 2 (A + Π) -1 j 1 ,j 2 (1 -Π)(1 -e v 1,j 1 A∆t )f n j 1 ,j 2 -∆t∆v 1 ∆v 2 (A + Π) -1 j 1 ,j 2 (1 -Π)(1 -e v 1,j 1 A∆t )( F Df ) n j 1 ,j 2 = E n 1 + ∆v 1 ∆v 2 (A + Π) -1 j 1 ,j 2 (1 -Π)(1 -e v 1,j 1 A∆t )f n j 1 ,j 2 -∆t∆v 1 ∆v 2 (A + Π) -1 (1 -Π) j 1 ,j 2 ( F Df ) n j 1 ,j 2 + ∆v 1 ∆v 2 (A + Π) -1 j 1 ,j 2 (1 -Π)(e v 1,j 1 A∆t f n j 1 ,j 2 -f n+1 j 1 ,j 2 ) = E n 1 + ∆v 1 ∆v 2 (A + Π) -1 j 1 ,j 2 ((1 -Π)f n j 1 ,j 2 -(1 -Π)f n+1 j 1 ,j 2 ).
By induction, if the discrete Poisson equation is satisfied at iteration n, then it is satisfied at iteration n + 1 and the proof is complete.

Vlasov-Maxwell 2dx-2dv

We finally consider the 2dx-2dv Vlasov-Maxwell model satisfied by

f (t, x, y, v 1 , v 2 ), E 1 (t, x, y), E 2 (t, x, y), B(t, x, y), with t ≥ 0, (x, y) ∈ [0, L x ]×[0, L y ] and (v 1 , v 2 ) ∈ R 2                  ∂ t f + v 1 ∂ x f + v 2 ∂ y f + E 1 ∂ v 1 f + E 2 ∂ v 2 f + B(v 2 ∂ v 1 f -v 1 ∂ v 2 f ) = 0, ∂ t B = ∂ y E 1 -∂ x E 2 , ∂ t E 1 = ∂ y B - ˆR2 v 1 f dv 1 dv 2 + J1 , ∂ t E 2 = -∂ x B - ˆR2 v 2 f dv 1 dv 2 + J2 , ∂ x E 1 + ∂ y E 2 = ˆR2 f dv 1 dv 2 -ρ, ∂ x B + ∂ y B = 0, (32)
with initial conditions (f 0 (x, y, v 1 , v 2 ), E 0 1 (x, y), E 0 2 (x, y), B(x, y)) such that the Poisson equation is satisfied

∂ x E 0 1 + ∂ y E 0 2 = ´R2 f 0 dv 1 dv 2 -
ρ and periodic boundary conditions are considered in space. here, Ji = 1

LxLy ´Lx×Ly ´R2

v i f dxdydv 1 dv 2 , ρ = 1 LxLy ´Lx×Ly ´R2 f dxdydv 1 dv 2 , E = (E 1 , E 2 ), B = (B 1 , B 2 ), ∇ = (∂ x , ∂ y ).

Exponential DG discretization

Here we apply 2D DG method in (x, y) direction and consider the discretization on Cartesian meshes with a rectangular triangulation I i × I j . We define the space V h as the space of tensor product piece-wise polynomials of degree at most k in each variable on every element, i.e.

V k h = {v h : v h | I i ×I j ∈ Q k (I i × I j )}, where Q k (I i × I j )
is the space of tensor products of one dimensional polynomials of degree up to k.

We follow the lines of the above subsections: we use a DG method in the 2D space direction (x, y) (with N x (resp. N y ) cells in direction x (resp. y) and a grid in the velocity direction v ℓ,j ℓ = v ℓ,min + j ℓ ∆v ℓ , ℓ = 1, 2, j ℓ = 1, . . . , N v . The 2D DG approximation for f is represented as (for j 1 , j 2 = 1, . . . , N v )

f (t, x, y, v j 1 , v j 2 ) ≈ f h (t, x, y, v j 1 , v j 2 ) = Nx i=1 Ny j=1 k m=0 k n=0 f mn ij (t, v j 1 , v j 2 )ξ m i (x)ξ n j (y).
For simplicity, we only present the 2D DG discretization for linear part of equation of f obtained by multiplying the Vlasov equation by ξ ℓ (x)ξ s (y) (for ℓ = 0, . . . , k and s = 0, . . . , k) on (x, y) ∈ I i × I j (for i = 1, 2, ...N x and j = 1, 2, ...N y ):

k m=0 k n=0 ∂ t f mn ij (t, v j 1 , v j 2 )(ξ m , ξ ℓ ) I i (ξ n , ξ s ) I j - k m=0 k n=0 v j 1 f mn ij (t, v j 1 , v j 2 )(ξ m , ∂ x ξ ℓ ) I i (ξ n , ξ s ) I j + v j 2 f mn ij (t, v j 1 , v j 2 )(ξ m , ξ ℓ ) I i (ξ n , ∂ y ξ s ) I j + k m=0 k n=0 (v j 1 {f h (t, x, y, v j 1 , v j 2 }ξ ℓ i+ 1 2 i-1 2 , ξ s ) I j + (v j 2 {f h (t, x, y, v j 1 , v j 2 }ξ s j+ 1 2 j-1 2 , ξ ℓ ) I i = 0, (33)
where we used the central fluxes in x and y, that is for the

x direction {f h (t, x, y, v j 1 , v j 2 )}| x i±1/2 = 1 2 (f h (t, x + i±1/2 , y, v j 1 , v j 2 ) + f h (t, x - i±1/2 , y, v j 1 , v j 2 )
). We consider f j 1 ,j 2 ∈ R (k+1) 2 NxNy the vector containing the degree of freedom f m,n i,j

f j 1 ,j 2 = [f 0,0 1,1 , f 1,0 1,1 , . . . , f k,0 1,1 , . . . , f k,0 Nx,1 , f 0,1 1,1 , . . . , f k,k Nx,Ny ] T j 1 ,j 2 (34)
and E 1 , E 2 , B ∈ R (k+1) 2 NxNy the vectors (defined as [START_REF] Tee | Eigenvectors of block circulant and alternating circulant matrices[END_REF]) containing the DG degree of freedom of the electromagnetic fields (E 1 , E 2 , B). Finally, we introduce ( F Df ) j 1 ,j 2 the DG approximation of the nonlinear term

(F • ∇ v f )(v j 1 , v j 2 ), with F = (E 1 + Bv 2 , E 2 -Bv 1 ) using similar
techniques as in the 1dx case. With these notations, we have the following semi-discretized scheme

                   ∂ t f j 1 ,j 2 = (1 y ⊗ v 1,j 1 A x )f j 1 ,j 2 + (v 2,j 2 A y ⊗ 1 x )f j 1 ,j 2 -( F Df ) j 1 ,j 2 ∂ t B = -(A y ⊗ 1 x )E 1 + (1 y ⊗ A x )E 2 , ∂ t E 1 = -(A y ⊗ 1 x )B - j 1 ,j 2 v 1,j 1 Pf j 1 ,j 2 ∆v 1 ∆v 2 , ∂ t E 2 = (1 y ⊗ A x )B - j 1 ,j 2 v 2,j 2 Pf j 1 ,j 2 ∆v 1 ∆v 2 , (35)
where A x ∈ M (k+1)Nx,(k+1)Nx (R) and A y ∈ M (k+1)Ny,(k+1)Ny (R) are the matrices coming from the DG semi-discretization in each space direction as before, ⊗ is the Kronecker product, P = (1 y -Π y ) ⊗ (1 x -Π x ) with 1 x (resp. 1 y ) the identity matrix of size (k + 1)N x (resp. (k + 1)N y ) and Π x (resp. Π y ) the projection onto Ker(A x ) (resp. Ker(A y )).

Before discussing the time discretization, we prove the following proposition.

Proposition 6. The semi-discretized system (35) satisfied by (f j 1 ,j 2 , B, E 1 , E 2 )(t) preserves the following discretized Poisson equation (provided that it is satisfied at time t = 0)

(1 y ⊗ A x )E 1 (t) + (A y ⊗ 1 x )E 2 (t) = - j 1 ,j 2 Pf j 1 ,j 2 (t)∆v 1 ∆v 2 , ∀t > 0.
Proof. Let us derive with respect to time the left hand side of the equality to get

∂ t ((1 y ⊗ A x )E 1 + (A y ⊗ 1 x )E 2 ) = (1 y ⊗ A x )∂ t E 1 + (A y ⊗ 1 x )∂ t E 2 = (1 y ⊗ A x )(-(A y ⊗ 1 x )B - j 1 ,j 2 v 1,j 1 Pf j 1 ,j 2 ∆v 1 ∆v 2 ) + (A y ⊗ 1 x )((1 y ⊗ A x )B - j 1 ,j 2 v 2,j 2 Pf j 1 ,j 2 ∆v 1 ∆v 2 ) = (A y ⊗ 1 x )(1 y ⊗ A x ) -(1 y ⊗ A x )(A y ⊗ 1 x ) B -(1 y ⊗ A x ) j 1 ,j 2 v 1,j 1 Pf j 1 ,j 2 ∆v 1 ∆v 2 ) -(A y ⊗ 1 x ) j 1 ,j 2 v 2,j 2 Pf j 1 ,j 2 ∆v 1 ∆v 2 ) = - j 1 ,j 2 ∂ t (Pf j 1 ,j 2 ) + P( F Df ) j 1 ,j 2 ∆v 1 ∆v 2 = -∂ t j 1 ,j 2 Pf j 1 ,j 2 ∆v 1 ∆v 2 ,
where we used the identities

(A y ⊗ 1 x )(1 y ⊗ A x ) = (A y 1 y ) ⊗ (1 x A x ) = A y ⊗ A x (1 y ⊗ A x )(A y ⊗ 1 x ) = (1 y A y ) ⊗ (A x 1 x ) = A y ⊗ A x
to pass from the third to the fourth equality. Integrating in time the obtained equality and assuming the discrete Poisson equation holds at time t = 0 leads to the result.

We end this part by giving some elements on the time discretization. First, in this case, it is difficult to compute the exponential of the linear part. However, we can consider the exponential of the f j 1 ,j 2 linear part (which corresponds to the (x, y) transport). Indeed, we observe from [START_REF] Verboncoeur | Particle simulation of plasmas: review and advances[END_REF] that this linear part writes

∂ t f j 1 ,j 2 = (1 y ⊗ v 1,j 1 A x ) + (v 2,j 2 A y ⊗ 1 x ) f j 1 ,j 2 = v 2,j 2 A y ⊕ v 1,j 1 A x f j 1 ,j 2 ,
where we used the definition of the Kronecker sum ⊕. The exact solution can be then written as

f j 1 ,j 2 (t) = exp (v 2,j 2 A y ⊕ v 1,j 1 A x ) t f j 1 ,j 2 (0) = exp v 2,j 2 tA y ⊕ v 1,j 1 tA x f j 1 ,j 2 (0).
It is well known that the exponential of a matrix with Kronecker sum structure is equal to the Kronecker product of the exponentials that is

f j 1 ,j 2 (t) = exp(v 2,j 2 tA y ) ⊗ exp(v 1,j 1 tA x )f j 1 ,j 2 (0),
which can be recast using the vec operation as

f j 1 ,j 2 (t) = exp(v 2,j 2 tA y )f j 1 ,j 2 (0) exp(∆v 1,j 1 t(A x ) T ), (36)
where f j 1 ,j 2 (s)=vec(f j 1 ,j 2 (s)) (for s = 0, t denotes the vectorization operation which takes the matrix f j 1 ,j 2 ∈ M (k+1)Nx,(k+1)Ny (R) as entry and gives the vector f j 1 ,j 2 ∈ R (k+1) 2 NxNy as a result. This means that the update of f requires matrix-vector products operations that only involves to assembly exponential of matrices A x and A y which are computed from the one-dimensional case (see ( 8)). Moreover, these matrix-vector products calculations can be performed in a very efficient way. This nice property has been exploited in the literature to design efficient routines for computing matrix exponentials [START_REF] Caliari | A µ-mode integrator for solving evolution equations in kronecker form[END_REF][START_REF] Croci | Exploiting Kronecker structure in exponential integrators: Fast approximation of the action of φ-functions of matrices via quadrature[END_REF][START_REF] Muñoz-Matute | Exploiting the Kronecker product structure of φ-functions in exponential integrators[END_REF].

Remark 5. The semi-discretized Vlasov-Maxwell system (35) can be degenerated to a semidiscretization of the Vlasov-Poisson system satisfied by (f, E 1 , E 2 ). In this case, the Lawson scheme only applies to the unknown f and then requires the calculation of exp(v 2,j 2 tA y ⊕v 1,j 1 tA x ) which can be performed efficiently thanks to [START_REF] Wang | A family of eulerian-lagrangian localized adjoint methods for multi-dimensional advection-reaction equations[END_REF]. The update of the electric field is performed using the Poisson equation thanks to the updated f .

Fourier based space discretization

In this part, we consider Fourier in space combined with finite differences in velocity to semidiscretize the Vlasov-Maxwell system (32) and we will see that in this case, it will be possible to compute explicitely the exponential of the linear part. Denoting fkx,ky,j 1 ,j 2 the Fourier coefficient of f in space and evaluated at the velocity grid introduced previously (k x , k y being the Fourier variables), Ê1,kx,ky , Ê2,kx,ky , Bkx,ky the Fourier coefficients of the electromagnetic fields (E 1 , E 2 , B), and introducing the force term F = (E 1 + Bv 2 , E 2 -Bv 1 ), we get the following semi-discretized scheme

                     ∂ t fkx,ky,j 1 ,j 2 + (v 1,j 1 ik x + v 2,j 2 ik y ) fkx,ky,j 1 ,j 2 + (FDf ) kx,ky,j 1 ,j 2 = 0, ∂ t Bkx,ky = ik y Ê1,kx,ky -ik x Ê2,kx,ky , ∂ t Ê1,kx,ky = ik y Bkx,ky - j 1 ,j 2 v 1,j 1 fkx,ky,j 1 ,j 2 ∆v 1 ∆v 2 + J1 , ∂ t Ê2,kx,ky = -ik x Bkx,ky - j 1 ,j 2 v 2,j 2 fkx,ky,j 1 ,j 2 ∆v 1 ∆v 2 + J2 , ( 37
)
with the initial conditions fkx,ky,j 1 ,j 2 (0), Bkx,ky (0), Ê1,kx,ky (0), Ê2,kx,ky (0) satisfying the Poisson equation ik x Ê1,kx,ky (0) + ik y Ê2,kx,ky (0) = j 1 ,j 2 fkx,ky,j 1 ,j 2 (0)∆v 1 ∆v 2 for (k x , k y ) ̸ = (0, 0).

For the semi-discretized system (37), we have a similar proposition as Prop [START_REF] Cai | An Eulerian-Lagrangian discontinuous Galerkin method for transport problems and its application to nonlinear dynamics[END_REF] in this Fourier case.

Proposition 7. The semi-discretized system (37) satisfied by ( fkx,ky,j 1 ,j 2 , Bkx,ky , Ê1,kx,ky , Ê2,kx,ky )(t) preserves the following discretized Poisson equation (provided that it is satisfied at time t = 0) ik x Ê1,kx,ky (t) + ik y Ê2,kx,ky (t) = j 1 ,j 2 fkx,ky,j 1 ,j 2 (t)∆v 1 ∆v 2 , (k x , k y ) ̸ = (0, 0), ∀t > 0, Proof. As in the proof of Prop ( 6), we take the derivative with respect to time of the left hand side to get ∂ t (ik x Ê1,kx,ky + ik y Ê2,kx,ky ) = ik x ik y Bkx,ky -

j 1 ,j 2 v 1,j 1 fkx,ky,j 1 ,j 2 ∆v 1 ∆v 2 + ik y -ik x Bkx,ky - j 1 ,j 2 v 2,j 2 fkx,ky,j 1 ,j 2 ∆v 1 ∆v 2 = j 1 ,j 2 ∆v 1 ∆v 2 ∂ t fkx,ky,j 1 ,j 2 + (FDf ) kx,ky,j 1 ,j 2 = ∂ t j 1 ,j 2 fkx,ky,j 1 ,j 2 ∆v 1 ∆v 2 ,
where we used the summation on j 1 , j 2 of the Vlasov equation together with the fact that D is a conservative finite difference operator. Finally, integrating the result in time between 0 and t and assuming the relation holds at t = 0 ends the proof.

To derive a fully discrete scheme, we introduce as previously the vector fkx,ky ∈ C Nv 1 Nv 2 and denote U (t) = ( f , B, Ê2 , Ê1 ) kx,ky (t) ∈ C Nv 1 Nv 2 +3 , then the system (37) can be rewritten as

∂ t U = LU + N (U ), ( 38
)
where L is given by

            -i(k x diag( ⃗ v 1 )+k y v 2,1 ) 0 Nv 1 ,Nv 1 . . . 0 Nv 1 ,Nv 1 0 Nv 1 ,1 . . . 0 Nv 1 ,1 0 Nv 1 ,Nv 1 -i(k x diag( ⃗ v 1 )+k y v 2,2
-i(k x diag( ⃗ v 1 )+k y v 2,Nv 2 ) 0 Nv 1 ,1 . . . 0 Nv 1 ,1 0 1,Nv 1 . . . 0 1,Nv 1 0 1,Nv 1 0 -ik x ik y -∆v 1 ∆v 2 v 2,1 1 -∆v 1 ∆v 2 v 2,2 1 . . . -∆v 1 ∆v 2 v 2,Nv y 1 -ik x 0 0 -∆v 1 ∆v 2 ⃗ v 1 -∆v 1 ∆v 2 ⃗ v 1 . . . -∆v 1 ∆v 2 ⃗ v 1 ik y 0 0             (39)
and fkx,ky = fkx,ky,1,1 , . . . , fkx,ky,Nv 1 ,1 , . . . , fkx,ky,Nv

1 ,Nv 2 T ∈ C Nv 1 Nv 2 , N (U ) = ( FDf ) kx,ky , 0, 0, 0 T ∈ C Nv 1 Nv 2 +3 ,
where we denote ⃗ v 1 the vector with components ( ⃗ v 1 ) j 1 = v 1,min + j 1 ∆v 1 and v 2,j 2 = v 2,min + j 2 ∆v 2 . In the same spirit as previously, diag( ⃗ v 1 ) denotes the diagonal matrix whose diagonal is composed of ⃗ v 1 . We now study the time discretization of (38) based on a Lawson scheme which requires to know exp(tL) with L given above. Similar (but more tedious) calculations to those performed in the proof of Prop 4 enable to get an explicit expression of exp(tL). To end this section, we prove that the following scheme, with U n ≈ U (t n ), t n = n∆t, ∆t > 0 and the notations introduced above

U n+1 = exp(∆tL)U n + ∆t exp(∆tL)N (U n ), (40)
that approximates the ODE (38) preserves a discrete Poisson equation. This is the object of the following proposition. f n kx,ky,j 1 ,j 2 ∆v 1 ∆v 2 , (k x , k y ) ̸ = (0, 0), ∀n ∈ N ⋆ , Proof. First, we need to know the shape of exp(∆tL). From the one-dimensional calculations and from [START_REF] Boutin | Modified lawson methods for Vlasov equations[END_REF], we have the representation of

e L∆t                e -i∆t(kx ⃗ v1+kyv2,1) 0 Nv 1 ,Nv 1 0 Nv 1 ,Nv 1 . . . 0 Nv 1 ,Nv 1 0 Nv 1 ,1 0 Nv 1 ,1 0 Nv 1 ,1 0 Nv 1 ,Nv 1 e -i∆t(kx ⃗ v1+kyv2,2) 0 Nv 1 ,Nv 1 . . . 0 Nv 1 ,
k 2 y cos(∆t|k|)+k 2 x |k| 2               
where we used the fact that the exponential of the homogeneous Maxwell part is

exp   t   0 -ik x ik y -ik x 0 0 ik y 0 0     =     cos(t|k|) -ikx sin(t|k|) |k| iky sin(t|k|) |k| -ikx sin(t|k|) |k| k 2 x cos(t|k|)+k 2 y |k| 2 kxky(1-cos(t|k|)) |k| 2 iky sin(t|k|) |k| kxky(1-cos(t|k|)) |k| 2 k 2 y cos(t|k|)+k 2 x |k| 2     ,
and the vectors ( e B, e E 2 , e E 1 ) kx,ky,⋆,j 2 ∈ C 3Nv 1 for all k x , k y and j 2 = 1, . . . , N v 2 will be given below. From the components on f , since we get a diagonal matrix, we have fkx,ky,j 1 ,j 2 (t) = e -i(t-t n )(kxv 1,j 1 +kyv 2,j 2 ) f n kx,ky,j 1 ,j 2

which can be inserted in the Maxwell part to compute the vectors ( e B, e E 2 , e E 1 ) kx,ky,⋆,j 2 . To do so, we consider the vector V (t) = ( B, Ê2 , Ê1 ) T (t) which solves the following ODE

∂ t V = J V +R with R(t) =    0 j 1 ,j 2 v 2,j 2 fkx,ky,j 1 ,j 2 (t)∆v 1 ∆v 2 j 1 ,j 2 v 1,j 1 fkx,ky,j 1 ,j 2 (t)∆v 1 ∆v 2    and J =   0 -ik x ik y -ik x 0 0 ik y 0 0   .
The variation of constant formula gives

V (t n+1 ) = e ∆tJ V (t n ) + ˆtn+1 t n e (t n+1 -t)J R(t)dt.
The calculations for the integral term involve the following integral term

I 1 ( f n ) = ˆtn+1 t n -ik x sin(|k|(t n+1 -t)) |k| R 2 (t) + ik y sin(|k|(t n+1 -t))) |k| R 3 (t) dt = j 1 ,j 2 (-ik x v 2,j 2 + ik y v 1,j 1 ) |k| C f n kx,ky,j 1 ,j 2 ∆v 1 ∆v 2 , = j 1 ,j 2
e B kx,ky,j 1 ,j 2 f n kx,ky,j 1 ,j 2 ,

I 2 ( f n ) = ˆtn+1 t n k 2 x cos(|k|(t n+1 -t)) + k 2 y |k| 2 R 2 (t) + k x k y (1 -cos(|k|(t n+1 -t))) |k| 2 R 3 (t) dt = j 1 ,j 2 v 2,j 2 k 2 x A + k 2 y B |k| 2 + v 1,j 1 k x k y (B -A) |k| 2 f n kx,ky,j 1 ,j 2 ∆v 1 ∆v 2 , = j 1 ,j 2
e E 2,kx,ky,j 1 ,j 2 f n kx,ky,j 1 ,j 2 ,

I 3 ( f n ) = ˆtn+1 t n k x k y (1 -cos(|k|(t n+1 -t))) |k| 2 R 2 (t) + k 2 y cos(|k|(t n+1 -t)) + k 2 x |k| 2 R 3 (t) dt = j 1 ,j 2 v 2,j 2 k x k y (B -A) |k| 2 + v 1,j 1 k 2 y A + k 2 x B |k| 2 f n kx,ky,j 1 ,j 2 ∆v 1 ∆v 2 , = j 1 ,j 2 e E 1,kx,ky,j 1 ,j 2 f n kx,ky,j 1 ,j 2 ,
where the time integrals are

A = ˆtn+1 t n cos(|k|(t n+1 -t))e -i(t-t n )(k•v j ) dt, B = ˆtn+1 t n e -i(t-t n )(k•v j ) dt = 1 ik • v j (1 -e -i∆tk•v j ), C = ˆtn+1 t n sin(|k|(t n+1 -t))e -i(t-t n )(k•v j ) dt, with k • v j = k x v 1,j 1 + k y v 2,j 2 .
To check the conservation of the Poisson equation, one focuses on the equations on Ê2 and Ê1 only. Thanks to the above calculations, we can write down the update of Ê2 , Ê1 using the first order Lawson scheme

Ên+1 2 = (e ∆tJ V n ) 2 + I 2 ( f n -∆t (FDf n )), Ên+1 1 = (e ∆tJ V n ) 3 + I 3 ( f n -∆t (FDf n )).
Thus, it remains to compute ik x Ên+1 1 + ik y Ên+1 2 using the last relations. First, one can check easily that

ik y (e ∆tJ V n ) 2 + ik x (e ∆tJ V n ) 3 = ik y -ik x sin(|k|∆t) |k| Bn + k 2 x cos(|k|∆t) + k 2 y |k| 2 Ên 2 + k x k y (1 -cos(|k|∆t)) |k| 2 Ên 1 + ik x ik y sin(|k|∆t) |k| Bn + k x k y (1 -cos(|k|∆t)) |k| 2 Ên 2 + k 2 y cos(|k|∆t) + k 2 x |k| 2 Ên 1 = ik y Ên 2 + ik x Ên 1 .
Then, we have

ik y Ên+1 2 + ik x Ên+1 1 = ik y Ên 2 + ik x Ên 1 + ik y I 2 ( f n -∆t (FDf n )) + ik x I 3 ( f n -∆t (FDf n )) = ik y Ên 2 + ik x Ên 1 + ik y j 1 ,j 2 v 2,j 2 ( f n -∆t (FDf n ))B∆v 2 ∆v 1 +ik x j 1 ,j 2 v 1,j 1 ( f n -∆t (FDf n ))B∆v 2 ∆v 1 = ik y Ên 2 + ik x Ên 1 + ∆v 2 ∆v 1 j 1 ,j 2 (ik • v j )( f n -∆t (FDf n )) 1 ik • v j (1 -e -i∆tk•v j ) = ik y Ên 2 + ik x Ên 1 + ∆v 2 ∆v 1 j 1 ,j 2 (1 -e -i∆tk•v j )f n -∆t j 1 ,j 2 (FDf n ) + ∆t j 1 ,j 2 e -i∆tk•v j (FDf n ) = ik y Ên 2 + ik x Ên 1 + ∆v 2 ∆v 1 j 1 ,j 2 (1 -e -i∆tk•v j ) f n + j 1 ,j 2 e -i∆tk•v j f n -f n+1 = ik y Ên 2 + ik x Ên 1 + ∆v 2 ∆v 1 j 1 ,j 2 ( f n -f n+1 ),
where we used the update for f :

f n+1 = e -i∆tk•v j f n -∆te -i∆tk•v j (FDf n
) and the conservation property of the discrete operator D. Then, if the Poisson equation is satisfied at iteration n, it is propagated to the next iteration, which concludes the proof.

Numerical experiments

In this section, we perform numerical experiments for linear transport problems and Vlasov equations. First, we study the different order of convergence on a linear problem and then, we present some numerical results of the exponential DG solutions for Vlasov equations in 1dx-1dv and 1dx-2dv cases.

2D linear passive-transport problems

We consider the following two-dimension linear transport equation

∂ t u + ∂ x u + ∂ v u = 0, (x, v) ∈ [0, 2π] 2 (42)
with the initial condition u(x, v, 0) = sin(x + v) and periodic boundary condition. The exact solution is u(x, v, t) = sin(x+v-2t) which enables us to check the different order of convergence. Indeed, for a Lawson scheme based on a underlying Runge-Kutta method RK(m, s) (order m, s stages, a DG space approximation with P k and a finite difference approximation in v of order 4 (which means D is chosen as a 4th order centered finite difference operator CD4), we expect the following estimate

∥u(t n ) -u n h ∥ L 2 ≤ C(∆x k + ∆v 4 + ∆t m )
Here we use the 3rd order Lawson-RK method for the time discretization, with a final time T = 1, and consider different parameters to test the convergence rates in t, x and v. We firstly take N v = 320, ∆t = 0.01, and consider different mesh size N x to check the convergence rate of DG in x direction for both central and upwind fluxes (in this linear case, upwind fluxes can be considereed easily). Table 1 shows the L ∞ and L 2 errors, the associated orders of convergence for DG-P k for k = 1, 2 in x-direction. The optimal convergence rate for DG is clearly obtained. In particular, the sub-optimal and optimal rates are observed according to the choice of the flux and to the oddness of k, as discussed in [START_REF] Liu | Sub-optimal convergence of discontinuous Galerkin methods with central fluxes for linear hyperbolic equations with even degree polynomial approximations[END_REF]. Then, we study the convergence in v-direction. We take k = 5, N x = 32, ∆t = 0.01, and consider different mesh size N v to check the convergence rate of the fourth order approximation of D (CD4). Table 2 shows the expected convergence (note that only central fluxes are considered in this case). Finally, to check time accuracy, we take k = 5, N x = 16, N v = 32 and ∆t = 0.0001 to compute a reference solution. Then we get the error table for different time step sizes ∆t. From Table 3, the expected 3rd order convergence is observed for all cases. 1: Linear transport equation: L ∞ and L 2 -norm space errors of the Lawson-DG scheme with P 1 , P 2 (CD4 in velocity (N v = 320) and RK [START_REF] Banks | High-order accurate conservative finite difference methods for Vlasov equations in 2D+ 2V[END_REF][START_REF] Banks | High-order accurate conservative finite difference methods for Vlasov equations in 2D+ 2V[END_REF] in time (∆t = 0.01)). 2: Linear transport equation: L ∞ and L 2 -norm velocity errors of the Lawson-DG scheme with CD4 (DG-P 5 in space (N x = 32) and RK(3,3) (∆t = 0.01)).

N x central flux upwind flux L ∞ -error order L 2 -error order L ∞ -
N v central flux L ∞ -

Vlasov-Ampère equation

We firstly consider the following initial condition for Landau damping

f 0 (x, v) = 1 √ 2π e -v 2 /2 (1 + α cos(kx)),
where we take x ∈ [0, 2π/k], k = 0.5, v ∈ [-9, 9] and α = 10 -3 . Here we still use DG method for space discretization, the finite difference method CD4 in v direction and the 3rd Lawson-RK method for time discretization (see Section 3.1). The numerical parameters are chosen as follows: ∆t = 0.1, N x = 31(P 2 ), N v = 121. In Figure 1, the time evolution of the electric energy ∥E(t)∥ L 2 is displayed in semi-log scale (with the corresponding damping rate in red) and the deviation of the total energy H(t) -H(0) with H(t) = ´v2 f (t, x, v)dxdv + ´E2 (t, x)dx. The expected behaviors (correct damping rate and good energy conservation) are recovered. We consider a second test called the two stream instability test with the initial condition

f 0 (x, v) = 1 √ 2π v 2 e -v 2 /2 (1 + α cos(kx)),
for which the same physical and numerical parameters as previously are kept except the final time which is T = 300. In Figures 2, we plot the time evolution of the electric energy in semi-log scale (and the corresponding instability rate in red) and the deviation of the total energy. For this test, a linear instability is first observed (up to t ≈ 30) during which a vortex in phase space is created (see 2), and it is followed by a nonlinear phase. These two behaviors are well reproduced by the scheme even if the mesh is quite coarse. In partiular, even if the vortex is well captured, we can observe spurious oscillations due to the use of central schemes. Note that the Poisson equation is satisfied in both cases up to machine accuracy.

Vlasov-Maxwell equations 1dx-2dv

We consider the Weibel instability [START_REF] Weibel | Spontaneously growing transverse waves in a plasma due to an anisotropic velocity distribution[END_REF] by consideering the Vlasov-Maxwell 1dx-2dv model studied in Section 3.2 with the initial distribution and fields are of the form

f (t = 0, x, v 1 , v 2 ) = 1 2πσ 1 σ 2 exp - 1 2 v 2 1 σ 2 1 + v 2 2 σ 2 2 (1 + α cos(k w x)), x ∈ [0, 2π/k), B(t = 0x) = β cos(kx), E 2 (t = 0, x) = 0,
and E 1 (x, t = 0) is imposed from the Poisson equation. We choose the parameters

σ 1 = 0.02/ √ 2, σ 2 = √ 12σ 1 , k = 1.
25, α = 0, β = -10 -4 for our test, which gives a growth rate of 0.02784 by solving the dispersion relation (see Weibel [START_REF] Weibel | Spontaneously growing transverse waves in a plasma due to an anisotropic velocity distribution[END_REF]). For the numerical simulations up to a final time T = 500, we still use DG method for space discretization in x-direction, finite difference method with CD4 in v direction and Lawson-RK(3,3) method for time discretization (Lawson-RK(3,3)-DG CD4) and consider ∆t = 0.5, N x = 21(P 2 ), N v 1 = N v 2 = 44. For comparison, we also consider Fourier method for space discretization in x-direction, finite difference with a third order upwind (UP3) in v direction and Lawson-RK [START_REF] Banks | High-order accurate conservative finite difference methods for Vlasov equations in 2D+ 2V[END_REF][START_REF] Banks | High-order accurate conservative finite difference methods for Vlasov equations in 2D+ 2V[END_REF] method in time (Lawson-RK(3,3)-Fourier UP3) with ∆t = 0.5,

N x = 64, N v 1 = N v 2 = 88.
In Figure 3, we show the time evolution of the electromagnetic energies ∥B(t)∥ L 2 , ∥E y (t)∥ L 2 , ∥E x (t)∥ L 2 (in semi-log scale) obtained by the two methods. First, we can observe that the theoretical growth rate is in very good agreement with the two numerical solution. Second, the two methods are very close up to time t ≈ 300 (which corresponds to the end of the linear phase) and slightly differs for larger times. We also show the evolution of the relative total energy |H(t) -H(0)|/H(0) in Figure 4 for Lawson-RK(3,3)-DG CD4 scheme without (left) and with (right) the energy correction step presented in [START_REF] Boutin | Modified lawson methods for Vlasov equations[END_REF] in the Vlasov-Ampère case. This projection approach enables to modify the unknown by a suitable coefficient which is of order the scheme so that the total energy is preserved almost up to machine error 10 -12 . The second test for Vlasov-Maxwell equation we considered is the streaming Weibel instability [START_REF] Califano | Spatial structure and time evolution of the weibel instability in collisionless inhomogeneous plasmas[END_REF][START_REF] Cheng | Discontinuous Galerkin methods for the Vlasov-Maxwell equations[END_REF] for which the initial condition is

f (t = 0, x, v 1 , v 2 ) = 1 2πσ 2 exp - v 2 1 2σ 2 δ exp - (v 2 -v 0,1 ) 2 2σ 2 + (1 -δ) exp - (v 2 -v 0,2 ) 2 2σ 2 , B(t = 0, x, t = 0) = β cos(kx), E 2 (t = 0, x) = 0,
and E 1 (t = 0, x) = 0 from the Poisson equation. We choose the parameters σ = 0.1/ √ 2, k = 0.2, β = -10 -3 , v 0,1 = 0.5, v 0,2 = -0.1 and δ = 1/6 for our test. The growth rate of E 2 is 0.03 [START_REF] Califano | Spatial structure and time evolution of the weibel instability in collisionless inhomogeneous plasmas[END_REF]. For the two schemes Lawson-RK(3,3)-DG CD4 and Lawson-RK(3,3)-Fourier UP3, we take the same parameters as in the previous case but ∆t = 0.1 for stability reasons. We show the result in Figure 5 in which the time evolution of the L 2 norm of the electromagnetic fields are displayed. First, we observe a good agreement with the theoretical growth rate for these two schemes and some deviations in the nonlinear phase can be observed. For this case, we also consider the correction on the total energy and plot the time history of the relative total energy in Figure 6 for Lawson-RK(3,3)-DG CD4 without (left) and with (right) energy correction step. The total energy can be well preserved for exponential DG approximations if it is stable, and the relative error is greatly improved with energy correction step without affecting the accuracy of the scheme. Note that the Poisson equation is satisfied in both cases up to machine accuracy.

Conclusion

In this work, we constructed and implemented a new exponential DG method for Vlasov equation, extending the previous works [START_REF] Crouseilles | Exponential methods for solving hyperbolic problems with application to collisionless kinetic equations[END_REF][START_REF] Crestetto | Comparison of high-order Eulerian methods for electron hybrid model[END_REF][START_REF] Boutin | Modified lawson methods for Vlasov equations[END_REF] on this topic where Fourier method in space were used. These methods allow to derive high order accuracy in time, space and velocity, still ensuring stability without the restrictive CFL type constraint coming from the linear part. Moreover, a discrete Poisson equation is satisfied and a projection technique enables to preserve the total energy. The extension to DG turns out to be an interesting alternative compared to previous approach based on Fourier which is restricted to cartesian domains with periodic boundary conditions. The approach only involves the calculation of exponential of DG-matrices of size (k + 1)N with k the DG degree and N the number of points in the space direction and we observe that thanks to the specific structure, this is also true in multi-dimensions.

One interesting extension is to perform an efficient implementation of the method in twodimension in space to exploit the Kronecker structure. We also plan to investigate the extension of this approach to problems involving boundary conditions, for which monotone fluxes are more appropriate but requires to study the stability of the scheme.

A Appendix A: The DG matrix construction

The goal of this appendix is to give some details on the construction of the DG-matrix ( 8) which is at the core of our scheme. To do so, we consider the 1dx transport equation ( 1) with coefficient a = 1, for which we wrote in (2) the semi-discrete DG scheme. To get the DG matrix, we consider ψ = ξ ℓ j and u

h (t, x) = k m=0 u m j (t)ξ m j (x) in (2) to get k m=0 (∂ t u m j (t)ξ m , ξ ℓ ) I j -(u m j (t)ξ m , ∂ x ξ ℓ ) I j + [{u h }ξ ℓ ] x j+1/2 x j-1/2 = 0, (43)
where ℓ = 0, . . . , k, j = 1, . . . , N and where the central flux

{u h }| x j±1/2 = (u + h (x j±1/2 ) + u - h (x j±1/2
))/2 is considered. Thus we obtain for the boundary term

[{u h }ξ ℓ ] x j+1/2 x j-1/2 = 1 2 k m=0 u m j (t)ξ m,-(x j+1/2 ) + u m j+1 (t)ξ m,+ (x j+1/2 ) ξ ℓ,-(x j+1/2 ) -u m j-1 (t)ξ m,-(x j-1/2 ) + u m j (t)ξ m,+ (x j-1/2 ) ξ l,+ (x j-1/2 ) = 1 2 k m=0 u m j (t)(1/2) m + u m j+1 (t)(-1/2) m (1/2) ℓ -u m j-1 (t)(1/2) m + u m j (t)(-1/2) m (-1/2) ℓ .
We denote u j (t) = (u 0 j (t), u 1 j (t), ...., u k j (t)) T , j = 1, 2, ..., N, then we can rewrite the DG discretization as

∆x      M M . . . M      ∂ t      u 1 u 2 . . . u N      -      D 1 D 2 . . . D 3 D 3 D 1 D 2 . . . . . . D 2 D 2 . . . D 3 D 1           u 1 u 2 . . . u N      = 0, (44)
where M, D i ∈ M (k+1),(k+1) (R), i = 1, 2, 3 are given by

M ℓ,m = (1/2) m+ℓ-1 m + ℓ -1 [1 -(-1) m+ℓ-1 ], (D 2 ) ℓ,m = (-1) m (1/2) m+ℓ-1 , (D 1 ) ℓ,m = (1/2) m+ℓ-2 ℓ -1 m + ℓ -2 - 1 2 [1 -(-1) m+ℓ-2 ], (D 3 ) ℓ,m = (-1) ℓ-1 (1/2) m+ℓ-1 ,
with (D 1 ) 1,1 = 0 by convention. Then we have ∆xM∂ t u -Du = 0, Lemma B.2. Let P ⊥ h q = q-P k q be the projection error. For any smooth function q(x), ∃ c > 0, such that

∥P ⊥ h q∥ D + h∥∂ x (P ⊥ h q)∥ D + h 1 2 ∥P ⊥ h q∥ ∞,D ⩽ ch k+1 |q| k+1,D , (47)
∥P ⊥ h q∥ Γ ⩽ ch k+ 1 2 ∥∂ k+1 x q∥, (48)
where the positive constant c is not dependent on h, solely depending on q, and D may be Ω or I j .

Furthermore, to avoid confusion with different constants, we denote a generic positive constant by C, which is independent of the numerical solution and the mesh size for our problem. But, the constant may dependent on the exact solution and may have a different value in each occurrence. Moreover, for problems considered in this paper, the exact solution is assumed to be smooth with periodic or compactly supported boundary condition. Therefore, the exact solution is always bounded.

We state the L 2 stability and L 2 -norm error estimate for the scheme and also give their proof.

Theorem B.1. For semi-discrete DG scheme (2) with central flux, we have the L 2 -stability:

d dt ||u h || 2 L 2 = 0.
Proof. Take the test function ψ = u h in the semi-discrete scheme (2), we have

1 2 d dt ˆIj u 2 h dx = -a{u h }| x j+ 1 2 u - h | x j+ 1 2 + a{u h }| x j-1 2 u + h | x j-1 2 + ˆIj au h (u h ) x dx, = -a{u h }| x j+ 1 2 u - h | x j+ 1 2 + a{u h }| x j-1 2 u + h | x j-1 2 + a 2 (u 2 h ) -| x j+ 1 2 - a 2 (u 2 h ) + | x j-1 2 . = - a 2 u + h u - h | x j+ 1 2 + a 2 u - h u + h | x j-1 2
.

(49) Sum over j of above equation, the L 2 -stability follows.

Theorem B.2. Let T > 0, u be the exact solution of problem (1), which is sufficiently smooth with bounded derivatives. Assume u h is the DG approximation of semi-discrete scheme (2) with the central flux and the approximation space V h is the space consisting of k-th piecewise polynomial. Then it holds that

∥u(T ) -u h (T )∥ L 2 ⩽ Ch k , ( 50
)
where C is a positive constant independent on ∆x.

Proof. Denote error as e u = u h -u. Notice that the scheme (2) is still satisfied with u h = u. So, we have the error equation ˆIj

(∂ t e u ψ)dx = -a{e u }| x j+ 1 2 ψ -| x j+ 1 2 + a{e u }| x j-1 2 ψ + | x j-1 2 + ˆIj ae u ψ x dx. ( 51
) Define e u = u h -u = (u h -P k u) -(u -P k u) = e u -P ⊥ k u. Then taking ψ = e u ,
+ a{P ⊥ k u}| x j+ 1 2 e u -| x j+ 1 2 -a{P ⊥ k u}| x j-1 2 e u + | x j-1 2 - ˆIj aP ⊥ k u( e u ) x dx. ( 52
) k = 2, u 2 = - 1 6 , 0, 1, 0, 0, -1 ∈R 6 , . . . , - 1 6
, 0, 1, 0, 0, -1

∈R 6 T ∈ R 3N . k = 3, u 2 = 0, - 3 20 , 0, 1 ∈R 4 , 0, - 3 20 , 0, 1 ∈R 4
, . . . , 0, -3 20 , 0, 1

∈R 4 T ∈ R 4N . k = 4, u 2 = - 3 280 , 0, 3 14 , 0,-1, 0, 0, - 3 14 , 0, 1 ∈R 10
, . . . , -3 280 , 0, 3 14 , 0, -1, 0, 0, -3 14 , 0, 1

∈R 10 T ∈ R 5N . k = 5, u 2 = 0, 5 336 , 0, - 5 18
, 0, 1

∈R 6 , 0, 5 336 , 0, - 5 18 , 0, 1 ∈R 6
, . . . , 0, 5 336 , 0, -

, 0, 1

∈R 6 T ∈ R 6N .
In the even case, when k is even, we observe a double pattern which is repeated N/2 (since when k is even, N is even to ensure (k + 1)N is even). Once we get (u 1 , u 2 ), the formula Πx = ⟨x, u 1 ⟩u 1 + ⟨x, u 2 ⟩u 2 enables to get Π.

D Appendix D: Lawson-Fourier method for Vlasov-Maxwell 1dx-2dv.

In this appendix, we extend the method presented in [START_REF] Boutin | Modified lawson methods for Vlasov equations[END_REF] to the Vlasov-Maxwell model in 1dx-2dv. This approach is compared to the Lawson-DG method in the numerical section 5.

Starting from the Vlasov-Maxwell 1dx-2dv model [START_REF] Johnson | An analysis of the discontinuous Galerkin method for a scalar hyperbolic equation[END_REF] satisfied by f (t, x, v 1 , v 2 ), E 1 (t, x), E 2 (t, x), B(t, x), with t ≥ 0, x ∈ [0, L] and (v 1 , v 2 ) ∈ R 2 , we shall use a Fourier method in the space direction x and we consider a grid in the velocity direction v ℓ,j ℓ = v ℓ,min + j ℓ ∆v ℓ , ℓ = 1, 2. Denoting fk,j 1 ,j 2 (t) the spatial Fourier coefficients of f (t, x, v j 1 , v j 2 ) and (E 1,k , E 2,k , B k )(t) the Fourier coefficients of (E 1 , E 2 , B)(t, x) then gives the following semidiscretized scheme for k = 1, . . . , N x (N x being the number of grid points in [0, L])

          
∂ t fk,j 1 ,j 2 + v 1,j 1 ik fk,j 1 ,j 2 + (FDf ) k,j 1 ,j 2 = 0, with F = (E 1 + Bv 2 , E 2 -Bv 1 ), ∂ t Bk = -ik Ê2,k , ∂ t Ê1,k =j 1 ,j 2 v 1,j 1 fk,j 1 ,j 2 ∆v 1 ∆v 2 , ∂ t Ê2,k = -ik Bkj 1 ,j 2 v 2,j 2 fk,j 1 ,j 2 ∆v 1 ∆v 2 , (55) with the initial condition fk,j 1 ,j 2 , Bk (0), Ê1,k (0), Ê2,k (0) satisfying the Poisson equation ik Ê1,k (0) = j 1 ,j 2 fk,j 1 ,j 2 (0)∆v 1 ∆v 2 for k ̸ = 0. Let us denote U = ( f , B, Ê2 , Ê1 ) ∈ M Nv 1 Nv 1 +3,Nv 1 Nv 1 +3 (C), then the previous system can be rewritten as . . .

∂ t U = LU + N (U ), (56)
0 1,Nv 1 0 1,Nv 1 0 -ik 0 -∆v 1 ∆v 2 ( ⃗ v 2 ) 1 1 -∆v 1 ∆v 2 ( ⃗ v 2 ) 2 1 . . . -∆v 1 ∆v 2 ( ⃗ v 2 ) Nvy 1 -ik 0 0 -∆v 1 ∆v 2 ⃗ v 1 -∆v 1 ∆v 2 ⃗ v 1 . . . -∆v 1 ∆v 2 ⃗ v 1 0 0 0              (57)
where we denote ⃗ v 1 the vector with components ( ⃗ v 1 ) j 1 = v 1,min + j 1 ∆v 1 and ⃗ v 2 the vector with components ( ⃗ v 2 ) j 2 = v 2,min + j 2 ∆v 2 . Moreover, diag( ⃗ v 1 ) denotes the diagonal matrix whose diagonal is composed of ⃗ v 1 , 1 ∈ R Nv 1 denotes the vector with components 1 and 0 A,B is a matrix with A lines and B columns with zeros. The size of the matrix L is (N v 1 N v 2 + 3) and in spite of its size, one can see that L is sparse.

A key point is to compute exp(L∆t) to design an exponential scheme approximating (56). Denoting U n ≈ U (t n ), t n = n∆t, ∆t > 0, the first order Lawson scheme is

U n+1 = exp(L∆t)U n + ∆t exp(L∆t)N (U n ).
The following proposition gives an explicit expression of exp(-∆tL).

Proposition 9. The exponential of the matrix L given by (57) is given by ( ⃗ v1) 0 Nv 1 ,Nv where the vectors ⃗ α, ⃗ β ∈ M 1,Nv 1 (C) whose components are given by α j 1 = -e -ik∆t 2(1 -v 1,j 1 ) -e ik∆t 2(1 + v 1,j 1 ) + e -ik∆tv 1,j 1 (1 -v 2 1,j 1 ) , β j 1 = -e -ik∆t 2(1 -v 1,j 1 ) + e ik∆t 2(1 + v 1,j 1 ) + v 1,j 1 e -ik∆tv 1,j 1 (1 -v 2 1,j 1 )

             e -ik∆tdiag
.

Proof. To compute exp(∆tL), we will solve exactly the linear part of the Vlasov-Maxwell system. First of all, we observe that the components associated to f is diagonal and can be solved independently so that the N v 1 × N v 2 top left block of exp(∆tL) is diagonal and is equal to e -ik∆t ⃗ v 1 (N v 2 times). Second, the 3 × 3 right bottom block corresponds to the homogeneous Maxwell equations and its exponential can be computed as It remains to compute the three last lines corresponding to the coupling between the Vlasov and Maxwell equations.

Computation of e E 1,k : solve Ê1,k First, we have for fk,j 1 ,j 2 (t) fk,j 1 ,j 2 (t) = e -ikv 1,j 1 (t-t n ) fk,j 1 ,j 2 (t n ), which enables to compute Ê1,k (t n+1 ) Ê1,k (t n+1 ) = Ê1,k (t n ) -

j 1 ,j 2 ˆtn+1
t n e -ikv 1,j 1 (t-t n ) dtv 1,j 1 f k,j 1 ,j 2 (t n )∆v 1 ∆v 2 = Ê1,k (t n ) + i∆v 1 ∆v 2 k j 1 ,j 2

(1 -e -ikv 1,j 1 ∆t ) fk,j 1 ,j 2 (t n ), from which we deduce the expression of e E 1,k .

Computation of e B k , e E 2,k : solve Bk , Ê2,k Next, we focus on the calculation of Ê2,k (t n+1 ) and Bk (t n+1 ). We write down the equations d dt Ê2,k (t) = -ik Bk (t)j 1 ,j 2 e -ikv 1,j 1 (t-t n ) v 2,j 1 f k,j 1 ,j 2 (t n )∆v 1 ∆v 2 d dt Bk (t) = -ik Ê2,k (t)

which can be rewritten as dU dt = M U + R with U (t) = ( Ê2,k (t), Bk (t)) and

M = 0 -ik -ik 0 and R = R 1 (t) 0 , with R 1 (t) = - j 1 ,j 2
e -ikv 1,j 1 (t-t n ) v 2,j 1 fk,j 1 ,j 2 (t n )∆v 1 ∆v 2 .

Thus, one can write the variation of constant formula

U (t n+1 ) = e M ∆t U (t n ) + ˆtn+1 t n
e -M (t-t n+1 ) R(t)dt.

First, e M ∆t reads as e M ∆t = exp 0 -ik∆t -ik∆t 0 = cos(kt) -i sin(kt) -i sin(kt) cos (kt) .

Second, one has to compute the integral term in (58) ˆtn+1

t n e -M (t-t n+1 ) R(t)dt = ´tn+1 t n cos(k(t -t n+1 ))R 1 (t)dt ´tn+1 t n i sin(k(t -t n+1 ))R 1 (t)dt = -j 1 ,j 2 I 1 v 2,j 1 f n k,j 1 ,j 2 ∆v 1 ∆v 2 j 1 ,j 2 I 2 v 2,j 1 f n k,j 1 ,j 2 ∆v 1 ∆v 2 (59)

where I 1 and I 2 are given by

I 1 = ˆtn+1 t n cos(k(t -t n+1
))e -ikv 1,j 1 (t-t n ) dt = ie -ik∆t 2k(1 -v 1,j 1 ) -ie ik∆t 2k(1 + v 1,j 1 ) -iv 1,j 1 e -ik∆tv 1,j 1 k(1 -v 2 1,j 1 )

,

I 2 = ˆtn+1 t n i sin(k(t -t n+1
))e -ikv 1,j 1 (t-t n ) dt = ie -ik∆t 2k(1 -v 1,j 1 ) + ie ik∆t 2k(1 + v 1,j 1 ) -ie -ik∆tv

Figure 1 :

 1 Figure 1: Vlasov-Ampère equation (Landau damping): time evolution of the electric energy in semi-log scale (left) and of the deviation of the total energy (right). Lawson-DG RK(3,3) and P 2 (∆t = 0.1, N x = 31(P 2 ), N v = 121).

  5.21E-06 2.99 2.31E-05 2.99 0.0250 6.51E-07 3.00 2.89E-06 3.00 0.0125 8.14E-08 3.00 3.61E-07 3.00 0.00625 1.02E-08 3.00 4.52E-08 3.00 Table3: Linear transport equation: L ∞ and L 2 -norm time errors of the Lawson-DG scheme with RK(3,3) (DG-P 5 in space (N x = 16) and CD4 (N v = 32)).

Figure 2 :

 2 Figure 2: Vlasov-Ampère equation (two stream instability): time evolution of the electric energy in semi-log scale (top left) and of the deviation of the total energy (top right), snapshot of f (t = 0) (bottom left) and snapshot of f (t = 300) (bottom right). Lawson-DG RK(3,3) and P 2 (∆t = 0.1, N x = 31(P 2 ), N v = 121).

Figure 3 :

 3 Figure 3: Vlasov-Maxwell equation (Weibel instability): time evolution of the electromagnetic energies in semi-log scale together with the analytic growth rate. Left: Lawson-RK(3,3)-DG CD4. Right: Lawson-RK(3,3)-Fourier UP3.

Figure 4 :

 4 Figure 4: Vlasov-Maxwell equation (Weibel instability): time evolution the relative total energy of Lawson-RK(3,3)-DG CD4 without (left) and with (right) energy correction step.

Figure 5 :

 5 Figure 5: Vlasov-Maxwell equation (streaming Weibel instability): time evolution of the electromagnetic energies in semi-log scale together with the analytic growth rate. Left: Lawson-RK(3,3)-DG CD4. Right: Lawson-RK(3,3)-Fourier UP3.

Figure 6 :

 6 Figure 6: Vlasov-Maxwell equation (streaming Weibel instability): time evolution the relative total energy of Lawson-RK(3,3)-DG CD4 without (left) and with (right) energy correction step.

withU 1 0

 1 = fk,⋆,1 , fk,⋆,2 , . . . , fk,⋆,Nv 2 , Bk , Ê2,k , Ê1,k T , fk,⋆,j 2 ∈ C Nv 1 , ∀k = 1, . . . , N x , ∀j 2 = 1, . . . , N v 2 , N (U ) = -(FDf ) k,⋆,1 , -(FDf ) k,⋆,2 , . . . , -(FDf ) k,⋆,Nv 2 , Nv 1 ,Nv 1 . . . 0 Nv 1 ,Nv 1 0 Nv 1 ,1 . . . 0 Nv 1 ,1

  Proposition 8. The Lawson scheme (40) satisfied by U n = ( f , B, Ê2 , Ê1 ) n kx,ky preserves the following Poisson equation (provided that it is satisfied at the initial time n = 0)

	ik x	Ên 1,kx,ky + ik y	Ên 2,kx,ky =
			j 1 ,j 2

  error order L 2 -error order

		10 1.62E-01	-	3.76E-01	-	4.78E-02	-	9.51E-02	-
	P 1	20 7.66E-02 1.08 1.86E-01 1.02 1.27E-02 1.91 2.42E-02 1.98 40 3.70E-02 1.05 9.24E-02 1.01 3.25E-03 1.97 6.07E-03 1.99
		80 1.82E-02 1.03 4.60E-02 1.00 8.23E-04 1.99 1.52E-03 2.00
		160 9.01E-03 1.01 2.30E-02 1.00 2.07E-04 1.99 3.79E-04 2.00
		10 2.29E-03	-	3.76E-03	-	2.52E-03	-	4.67E-03	-
	P 2	20 2.65E-04 3.11 4.53E-04 3.05 3.07E-04 3.04 5.83E-04 3.00 40 3.25E-05 3.02 5.63E-05 3.01 3.84E-05 3.00 7.29E-05 3.00
		80 4.05E-06 3.01 7.03E-06 3.00 4.79E-06 3.00 9.11E-06 3.00
		160 5.12E-07 2.98 8.99E-07 2.97 5.95E-07 3.00 1.16E-06 2.98
	Table

  we have ˆIj ∂

t ( e u ) e u dx = ˆIj ∂ t (P ⊥ k u) e u dx -a{ e u }| x j+ 1 2 e u -| x j+ 1 2 + a{ e u }| x j-1 2 e u + | x j-1 2 + ˆIj a e u ( e u ) x dx

0

  Nv 1 ,Nv 1 -ikdiag( ⃗ v 1 ) 0 Nv 1 ,Nv 1 . . . 0 Nv 1 ,Nv 1

					. . .	. . .	. . .
	. . . 0 Nv 1 ,Nv 1	. . . . . .	. . . 0 Nv 1 ,Nv 1	0 Nv 1 ,Nv 1 -ikdiag( ⃗ v 1 )	. . . 0 Nv 1 ,1	. . . . . .	. . . 0 Nv 1 ,1
	0 1,Nv 1

  1 0 Nv 1 ,Nv 1 . . . 0 Nv 1 ,Nv 1 0 Nv 1 ,1 0 Nv 1 ,1 0 Nv 1 ,1 0 Nv 1 ,Nv 1 e -ik∆tdiag( ⃗ v1) 0 Nv 1 ,Nv 1 . . . 0 Nv 1 ,Nv 1 where e B k,j 2 , e E 2,k,j 2 , e E 1,k ∈ M 1,Nv 1 (C) for j 2 = 1, . . . , N v 2 are given by

									
	. . . 0 Nv 1 ,Nv 1 e B k,1 e E 2,k,1	. . . . . . e B k,2 e E 2,k,2	. . . 0 Nv 1 ,Nv 1 . . . . . .	. . . . . . v1) 0 Nv 1 ,1 e -ik∆tdiag( ⃗ 0 Nv 1 ,Nv 1 e B k,Nv 2 cos(k∆t) e E 2,k,Nv 2 -i sin(k∆t) cos(k∆t) . . . . . . 0 Nv 1 ,1 -i sin(k∆t) 0 . . . . . . 0 Nv 1 ,1 0	           
	e E 1,k	e E 1,k	. . .	e E 1,k		0	0	1
	e B k,j 2 =	i∆v 1 ∆v 2 k	v 2,j 2 ⃗ α, e E 2,k,j 2 =	i∆v 1 ∆v 2 k	v 2,j 2	⃗ β, e E 1,k =	i∆v 1 ∆v 2 k	(1 -e -ik ⃗ v 1 ),

  1,j 1 k(1 -v 2 1,j 1 ).Inserting these calculations in (58) leads to the following expression for Ê2,k (t)Ê2,k (t n+1 ) = cos(k∆t) Ê2,k (t n ) -i sin(k∆t) Bk (t n ) + i∆v 1 ∆v 2 k j 1 ,j 2 v 2,j 2 β j 1 fk,j 1 ,j 2 (t n ), Bk (t n+1 ) = -i sin(k∆t)E 2,k (t n ) + cos(k∆t)B k (t n ) + i∆v 1 ∆v 2 k j 1 ,j 2 v 2,j 2 α j 1 fk,j 1 ,j 2 (t n ),where⃗ β = [β 1 , β 2 , . . . , β Nv 1 ] ∈ C Nv 1 and ⃗ α = [α 1 , α 2 , . . . , α Nv 1 ] ∈ C Nv 1 are given byWe conclude by writing the vectors e B k , e E 2,k corresponding to Bk and Ê2,k

	β j 1 = -	e -ik∆t 2(1 -v 1,j 1 )	+	e ik∆t 2(1 + v 1,j 1 )	+	v 1,j 1 e -ik∆tv 1,j 1 (1 -v 2 1,j 1 )	,
	α j 1 = -	e -ik∆t 2(1 -v 1,j 1 )	-	e ik∆t 2(1 + v 1,j 1 )	+	e -ik∆tv 1,j 1 (1 -v 2 1,j 1 )	.
	e B k,j 2 =	i∆v 1 ∆v 2 k	v 2,j 2 ⃗ α, e E 2,k,j 2 =	i∆v 1 ∆v 2 k	v 2,j 2	⃗ β.

where M is a block diagonal mass matrix of size N (k + 1), D is a block tridiagonal matrix of size N (k + 1), and u is the vector containing the degree of freedom u = (u 0 1 , u 1 1 , . . . , u k 1 , u 0 2 , u 1 2 , . . . , u k 2 , . . . , u 0 N , u 1 N , . . . , u k N ) T ∈ R (k+1)N .

Now we can rewrite the DG scheme as

where A is a block circulant matrix

or with the circblock notation

B Appendix B: Stability and error estimate for semi-discrete DG scheme.

In this appendix, we give some error estimate of the exponential-DG scheme for the one dimensional linear advection equation (1). To do so, we first define some notations about norms which will be used. For a given function x → v(x), we denote ∥v∥ j and ∥v∥ ∞,j as the L 2 -norm and L ∞ -norm of v on I j (j = 1, . . . , N ) respectively. Moreover,

where we express the value of v on the left and right limits of the grid point x j+ 1 2 with v -

respectively. Define the jump and the mean of v at x j-

) and {v} j-

B.1 Notations for projections and some properties of approximation space

The inverse properties of the finite space V h will be used.

where the positive constant C is independent of h and v.

Define the L 2 -projection P k of u into V h as follows:

The following lemma states the error of these projections [START_REF] Ciarlet | The finite element method for elliptic problems[END_REF].

By the the definition of the projections and some calculations, the right terms RHS of the error equation ( 52) become

Furthermore, from Holder's inequality, Lemma (B.1) and Lemma (B.2), we have

Thus by Gronwall's inequality, the conclusion in Theorem B.2 follows.

C Appendix C: Ker(A) and projection Π

Here, the projection matrix Π onto Ker(A) with A the DG-matrix ( 8) is discussed. A general expression (for any k, N ) turns out to be difficult and we compute Ker(A) (and the projection Π) for several practical cases. As mentioned in Remark 2, there are mainly two cases, according to the oddness of (k+1)N : (i) if (k+1)N is odd, 0 is a simple eigenvalue of A and Ker(A)=Span(u 1 ) where u 1 ∈ R (k+1)N correspond constants in the space P k ; (ii) if (k + 1)N is even, 0 is a double eigenvalue of A and Ker(A)=Span(u 1 , u 2 ) and u 2 ∈ R (k+1)N has to be determined. This second case recalls what happens for the second centered finite differences in which constant vector belongs to the kernel but also the sequence (-1) j .

• odd case: We can check that

is a eigenvector of A associated to the eigenvalue 0. By Πx = ⟨x, u 1 ⟩u 1 , we get the expression of the matrix Π Π = [u 1 , 0, . . . , 0, u 1 , 0, . . . , 0, . . . , u 1 , 0, . . . , 0 ∈ M (k+1)N,(k+1)N (R), with 0 ∈ R (k+1)N .

• even case: in addition to u 1 , we need to find a second eigenvector to construct Π. We give below the expression of u 2 for some k = 1 to k = 5 k = 0, u 2 = 0, 1, 0, 1, . . . , 0, 1

, . . . , 0, 1


A new MILP model for integrated lot-sizing and energy supply planning with onsite renewable energy

Ruiwen Liao 1 , Céline Gicquel 1

Université Paris-Saclay, LISN ruiwen.liao@lisn.fr, celine.gicquel@lri.fr

Mots-clés : energy-efficient manufacturing, production planning, proportional lot-sizing and scheduling problem, sequence-dependent setup costs and times, renewable energy

With the increasing concern on climate changes and global warming, industrial companies, one of the largest carbon emitters, are urged to reduce the carbon emissions of their manufacturing processes. The recent energy crisis also caused a sharp rise in the price of conventional energy resources and strong pressures on energy supply. One of the solutions to deal with these two challenges is to power, at least partially, manufacturing processes by electricity generated on site from renewable sources. However, the intermittence of the renewable energy sources as well as the time-of-use pricing policy frequently employed by energy providers make it necessary to accurately track the timing and amount of energy consumption, on-site generation and trading. This means that an integrated production planning and energy supply planning problem should be investigated.

The present work focuses on simultaneously planning the industrial operations and energy supply of the system depicted in Figure 1.

The production system comprises a single machine producing several types of products to satisfy an external customer demand. Backlogging is not allowed and demand must be met on time. In case finished products are produced in advance, they are kept in inventory, which is penalized by inventory holding costs. The machine may be either idle, carrying a changeover between two products, preserving its setup state for a given product or manufacturing a single type of product. Changing the setup state of the machine from one product to another one requires to carry out changeover operations, which incurs sequence-dependent changeover costs and times. The industrial production planning to be solved is thus a multi-item single-machine lot-sizing problem with sequence-dependent changeover costs and times.

Energy is consumed during changeover operations, manufacturing and preserving. The amount of energy consumed during a changeover is sequence-dependent, the one consumed during manufacturing is proportional to the number of items produced and the one consumed for preserving the machine setup status is proportional to the length of time spent in this state. No energy is consumed when the machine is idle. As for the energy supply system, it comprises three main elements : on-site power generation devices producing a time-varying amount of green electricity from renewable energy sources, an on-site energy storage system and the main electricity grid with which electricity may be traded at a time-varying price. Note that the one-site generated green electricity is assumed to be free and that trading energy with the grid incurs losses due to the electric transformer. This means that consuming the on-site generated green energy to power the manufacturing process will be in most cases more profitable than trading energy with the main grid.

We propose a new MILP model for this problem. It relies on a two-level time structure : a coarse time discretization to track the demand satisfaction (typically using periods corresponding to 8-hours shifts or days) and a fine time discretization to track the energy consumption,

). This fine time discretization is also used to plan the industrial production. The production planning problem is thus modeled as a Proportional Lot-sizing and Scheduling Problem (PLSP), a small-bucked lot-sizing model in which at most two different types of product may be produced in each micro-period. As our problem involves changeover operations that may overlap several planning micro-periods, we propose to extend the formulation recently introduced in [START_REF] Kaczmarczyk | Explicit modeling of multi-period setup times in proportional lot-sizing and scheduling problem with variable capacity[END_REF] to take into account sequence-dependent changeover costs and times and to integrate the energy supply planning into the model. Note that a similar integrated problem was studied in [START_REF] Gerhard Wichmann | Energyoriented lot-sizing and scheduling considering energy storages[END_REF]. The authors proposed to use three distinct time structures : one to track the demand satisfaction, one to build the production plan and one to track the energy supply. However, their model, based on an extension of the General Lot-sizing Problem (GLSP), leads to a large-size MILP formulation involving many big-M constraints and poses thus computational difficulties.

Numerical experiments were carried out on randomly generated instances to compare the quality of the energy supply and production plans obtained by our PLSP-based model with the one of the plans obtained with the GLSP-based model introduced in [START_REF] Gerhard Wichmann | Energyoriented lot-sizing and scheduling considering energy storages[END_REF]. Our preliminary computational results show that for small-size instances involving 3 products and 16 microperiods, the problem can be solved to optimality much faster by using the PLSP-based model than by using the GLSP-based model. Moreover, the cost of the production plan obtained with our model is the same as or even better than the one obtained in the GLSP-based model (since the GLSP-based model does not allow changeover overlapping two macro-periods). For medium and large instances involving up to 5 products and 96 micro-periods, the solver can hardly find a feasible solution when using the GLSP-based model, while it finds a relatively good feasible solution when using the PLSP-based model.