N Besse

E Deriaz

É Madaule

Adaptive multiresolution semi-Lagrangian discontinuous Galerkin methods for the Vlasov equations

Keywords: Vlasov-Poisson, discontinuous Galerkin methods, multi-wavelets, multiresolution analysis, semi-Lagrangian methods, adaptive mesh renement

We develop adaptive numerical schemes for the Vlasov equation by combining discontinuous Galerkin discretisation, multiresolution analysis and semi-Lagrangian time integration. We implement a tree based structure in order to achieve adaptivity. Both multi-wavelets and discontinuous Galerkin rely on a local polynomial basis. The schemes are tested and validated using Vlasov-Poisson equations for plasma physics and astrophysics.

Introduction

Many numerical experiments are performed on the Vlasov-Poisson problem since it is a well known (but not fully understood) system from plasma physics and a major issue for future simulation of large scale plasmas. Our goal is to develop adaptive numerical schemes using discontinuous Galerkin discretisation combined with semi-Lagrangian description whose mesh renement is based on multi-wavelets. In [START_REF] Cohen | Full adaptive multiresolution nite volume scheme for conservation laws[END_REF], the authors start from the standard formulation of nite-volume schemes for conservation laws and transform them into adaptive nite-volume schemes by using multiscale representations of the data. In the same philosophy, starting from the standard formulation of the discontinuous Galerkin method for conservation laws the authors of Refs. [START_REF] Archibald | Adaptive discontinuous Galerkin methods in multiwavelets bases[END_REF][START_REF] Shelton | A multi-resolution discontinuous Galerkin method for unsteady compressible ows[END_REF][START_REF] Gerhard | Adaptive multiresolution discontinuous Galerkin Schemes for conservation laws: multi-dimensional case[END_REF][START_REF] Hovhannisyan | Adaptive multiresolution discontinuous Schemes for conservation laws[END_REF][START_REF] Gerhard | A High-Order Discontinuous Galerkin Discretization with Multiwavelet-Based Grid Adaptation for Compressible Flows[END_REF] use a multi-wavelet representation of the unknowns to put these schemes into an adaptive setting. Following the spirit of these works, and in order to solve the transport and continuity equations, we use a multi-wavelet representation of the unknowns to transform the semi-Lagrangian discontinuous Galerkin schemes of references [START_REF] Guo | A conservative semi-Lagrangian discontinuous Galerkin scheme on the cubed-sphere[END_REF][START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov Poisson system[END_REF] into adaptive multiresolution versions. The discontinuous Galerkin formulation enables high-order accuracy with local data for computation. It has recently been widely studied by Heath [START_REF] Heath | Analysis of the Discontinuous Galerkin Method Applied to Collisionless Plasma Physics[END_REF], Ayuso de Dios et al. 1 Laboratoire J.-L. Lagrange, UMR CNRS/OCA/UCA 7293, Université Côte d'Azur, Observatoire de la Côte d'Azur, Bd de l'Observatoire CS 34229, 06304 Nice Cedex 4 France 2 Institut Jean Lamour, UMR 7198 CNRS -Université de Lorraine, Faculté des Sciences et Technologies -Bâtiment 2ème cycle, F-54506 Vand÷uvre-lès-Nancy cedex [START_REF] Ayuso De Dios | Discontinuous Galerkin methods for the one-dimensional Vlasov-Poisson system[END_REF][START_REF] Ayuso De Dios | Discontinuous Galerkin methods for the multi-dimensional Vlasov-Poisson problem[END_REF], Cheng et al. [START_REF] Cheng | Study of conservation and recurrence of Runge-Kutta discontinuous Galerkin scheme for Vlasov-Poisson system[END_REF], Rossmanith et al. [START_REF] Rossmanith | A positivity-preserving high-order semi-Lagrangian discontinuous Galerkin scheme for the Vlasov-Poisson equations[END_REF], Restelli et al. [START_REF] Madaule | Energy conserving discontinuous Galerkin spectral element method for the Vlasov-Poisson system[END_REF], etc. Those are done with Eulerian description of the time resolution but Guo et al. [START_REF] Guo | A conservative semi-Lagrangian discontinuous Galerkin scheme on the cubed-sphere[END_REF] or Qiu and Shu [START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov Poisson system[END_REF], Crouseilles et al. [START_REF] Crouseille | Discontinuous Galerkin semi-Lagrangian method for Vlasov-Poisson[END_REF] or Bokanowski and Simarta [START_REF] Bokanowski | Simamarta Semi-Lagrangian discontinuous Galerkin schemes for some rst and second-order partial dierential equations[END_REF] performed semi-Lagrangian time resolution. We use multi-wavelets framework for the adaptive part (more precisely, for the multiscale representation of the distribution function of particles). Those have been heavily studied by Alpert et al. [START_REF] Alpert | Adaptive Solution of Partial Dierential Equations in Multiwavelet Bases[END_REF][START_REF] Alpert | A Class of Bases in L 2 for the Sparse Representation of Integral Operator[END_REF]. Some works merging multiscale resolution and discontinuous Galerkin methods have been described by Archibald et al. [START_REF] Archibald | Adaptive discontinuous Galerkin methods in multiwavelets bases[END_REF] for convection based problems, by Müller et al. for non-linear hyperbolic conservation laws in the nite volume framework [START_REF] Hovhannisyan | Adaptive multiresolution discontinuous Schemes for conservation laws[END_REF][START_REF] Gerhard | Adaptive multiresolution discontinuous Galerkin Schemes for conservation laws: multi-dimensional case[END_REF] or for compressible ows [START_REF] Gerhard | A High-Order Discontinuous Galerkin Discretization with Multiwavelet-Based Grid Adaptation for Compressible Flows[END_REF]. In the framework of relativistic Vlasov equation, Besse et al. [START_REF] Besse | A wavelet-MRA-based adaptive semi-Lagrangian method for the relativistic Vlasov-Maxwell system[END_REF] presented the advantage of using adaptive meshes. While they used wavelet representation, which requires large data stencil, multi-wavelet representation coupled to discontinuous Galerkin discretisation only requires local stencil. This favours the parallelisation but an ecient parallel version of our numerical schemes (presented hereafter) will be the matter of future work. Here we start with two semi-Lagrangian discontinuous Galerkin schemes which are presented in [START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov Poisson system[END_REF] and [START_REF] Guo | A conservative semi-Lagrangian discontinuous Galerkin scheme on the cubed-sphere[END_REF]. We call them the SLDG [START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov Poisson system[END_REF] and CDG [START_REF] Guo | A conservative semi-Lagrangian discontinuous Galerkin scheme on the cubed-sphere[END_REF] schemes. We then modify these schemes to obtain adaptive multiresolution schemes with an adaptive mesh coarsening and renement procedure. To achieve this aim, among other things detailed in Secs. 4, 5 and 6, the data are represented in a dierent nite-dimensional space which is constructed by using a multi-wavelet basis [START_REF] Alpert | Adaptive Solution of Partial Dierential Equations in Multiwavelet Bases[END_REF][START_REF] Alpert | A Class of Bases in L 2 for the Sparse Representation of Integral Operator[END_REF]. The multiscale representation, given by the multi-wavelet basis, allows to give a well-dened and natural link between local small scales of the distribution function (or data) and the local size of the adaptive mesh. Then multi-wavelet representation of data is the relevant tool to develop natural and ecient criteria for designing an adaptive mesh (renement versus coarsening) which follows the multiscale development of the distribution function. Therefore we obtain new numerical schemes that we call AMW-SLDG (Adaptive Multi-Wavelet SLDG) and AMW-CDG.

Let us note that the SLDG and CDG schemes are equivalent when we consider the one-dimensional advection or continuity equations on a one-dimensional grid with periodic boundary conditions. The proof is given in [START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov Poisson system[END_REF], where the authors assumed that all integrals involved in the schemes are exactly computed. Nevertheless it is important to start with two non-adaptive methods which can be proven equivalent to see if our adaptive procedure introduces a bias or some dierences between them. For the moment there is no mathematical proof that the AMW-SLDG and AMW-CDG schemes are equivalent. It is not only a difcult task but it remains a true open question. Actually we do not know if they are rigorously equivalent. What is certain is that the proof performed in [START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov Poisson system[END_REF] can not be extended for non-conformal adaptive n-dimensional mesh (with n ≥ 2), which is our case. In addition in our case, where the mesh is not uniform and non-conformal, integrals are not computed exactly. Nevertheless we may expect that the AMW-SLDG and AMW-CGD schemes are close to each other, and that our procedure does not break the "symmetry" between these two schemes. A way to conrm such intuition is to make systematic comparisons of the AMW-SLDG and AMW-CDG schemes.

The paper is organised as follows. Sec. 2 presents the Vlasov-Poisson equations and their conservation properties. Sec. 3 deals with several discontinuous Galerkin methods and their application in semi-Lagrangian framework. Sec. 4 describes the multi-wavelets used to get a multiscale reconstruction of the distribution function. In Sec. 5 we identify some numerical error sources which are responsible for a loss of some conservation laws. A brief overview of numerical implementation and a summary of the full algorithm are presented in Sec. [START_REF] Besse | Validity of quasilinear theory: refutations and new numerical conrmation[END_REF]. Some numerical results in plasma physics and astrophysics are presented in Sec. 7. Finally, the conclusion reminds the main points of our paper and gives future perspectives of work.

The Vlasov-Poisson equation

We considered the following Vlasov equation

∂ t f (x, v, t) + ∂ x (vf (x, v, t)) + ∂ v (E(x, t)f (x, v, t)) = 0, x ∈ Ω x , v ∈ R, f (x, v, 0) = f 0 (x, v), (1)
where f (x, v, t) is the phase-space distribution function and E(x, t) is the electric eld obtained from the Poisson equation

-∆Φ(x, t) = ρ(x, t) -1, x ∈ Ω x (2a) E(x, t) = -∂ x Φ(x, t), (2b)
with Φ(x, t) the electric potential and ρ(x, t) the charge density linked to the distribution function f according to

ρ(x, t) = R f (x, v, t) dv. ( 3
)
This equation is written in adimensionnal variables using a xed background of ions such that the global plasma is neutral and the integral of the distribution function on the complete phase-space is equal to 1. This system is closed using periodic boundary conditions and zero average for E(x, t). The set Ω x is a compact domain of physical space with periodic boundary conditions. It is impossible from a practical point of view to use an innite domain in v, except if the initial distribution function is compactly supported. For this reason most of the time we consider a domain Ω v = [-L, L] with L large enough to ensure Dirichlet boundary conditions. In our simulations, given our initial conditions, we know that the solution remains bounded in velocity space by a decreasing exponential (uniformly in x). If we take a domain large enough in the velocity direction, integrals are small enough to ensure that the error is bounded by a small exponential. Such estimation has been used in many semi-Lagrangian codes, such as the one from [START_REF] Besse | Validity of quasilinear theory: refutations and new numerical conrmation[END_REF]. We call

Ω = Ω x × Ω v .
Those equations represent the time evolution of a collisionless plasma distribution function and of its self-consistent electric eld. This system satises many conservation laws that allow to qualify performances of the numerical methods.

In order to preserve the maximum principle (in particular the positivity) for the distribution function, we could obviously use some limiters such as the positivity-preserving limiter of Qiu and Shu [START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov Poisson system[END_REF]. Nevertheless our aim is to evaluate the intrinsic properties of our schemes without adding an a posteriori procedure to enforce the positivity or the maximum principle. Moreover limiters might add a bit of diusion sometimes. Therefore in our numerical simulations no limiter is used.

Discontinuous Galerkin formulations

General formulation

We explain here the principle of the discontinuous Galerkin formulation on a one-dimensional conservation law. Cockburn and Shu presented a large analysis of this formulation, mostly coupled with Runge-Kutta time resolution in, e.g. [START_REF] Cockburn | RungeKutta Discontinuous Galerkin Methods for Convection-Dominated Problems[END_REF]. We perform the time resolution of the Vlasov-Poisson problem using Strang splitting. Indeed, here we use an (x,v)-area-preserving integrator in time obtained by decomposing the dynamics in integrable Hamiltonian steps. A natural and simple choice is the Strang time-splitting strategy also known as the centred leapfrog symplectic integrator. Roughly speaking, this intensively used scheme (see [START_REF] Cheng | The integration of the Vlasov equation in conguration space[END_REF][START_REF] Feix | A Universal Model: The Vlasov Equation[END_REF][START_REF] Besse | Convergence of classes of high-order semi-Lagrangian schemes for the Vlasov-Poisson system[END_REF][START_REF] Besse | Validity of quasilinear theory: refutations and new numerical conrmation[END_REF][START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov Poisson system[END_REF][START_REF] Colombi | Vlasov versus Nbody: the Hénon sphere[END_REF] and references therein) consists in splitting the full transport operator into two easily integrable Hamiltonian transport operators one in the physical space, the second one in the velocity and solving them successively in a right order to get high-order approximation in time of the complete transport operator. For the Strang splitting, the local (respectively global) error in time is of order three (respectively two) [START_REF] Besse | A wavelet-MRA-based adaptive semi-Lagrangian method for the relativistic Vlasov-Maxwell system[END_REF][START_REF] Einkemmer | Convergence Analysis of Strang Splitting for Vlasov-Type Equations[END_REF].

Let us consider the conservation equation

∂ t f (x, t) + ∂ x (a(x, t)f (x, t)) = 0. (9)
The discontinuous Galerkin formulation consists in solving a weak formulation of the problem on each cell of the mesh with a space of test functions. More precisely we rst multiply (9) by a test function g ∈ G with G the space of test functions which we do not specify yet. We then integrate it on each cell

I i = [x i-1 2 , x i+ 1 2
] of the domain partition to obtain the weak formulation

Ii ∂ t f (x, t)g(x) dx + Ii ∂ x (a(x, t)f (x, t))g(x) dx = 0. (10)
The next stage is to perform an integration by part to transfer the space derivative onto the test function.

Ii ∂ t f (x, t)g(x) dx - Ii a(x, t)f (x, t)∂ x g(x) dx + [a(x, t)f (x, t)g(x)] x i+ 1 2 x i-1 2 = 0.
(11) Here comes the choice of the test function space G and numerical uxes since the ux terms in [START_REF] Cheng | Numerical study of Vlasov-Poisson equations for innite homogeneous stellar systems[END_REF] is meaningless when using discontinuous reconstruction at the cell interfaces. This choice will determine the properties (stability, consistency, accuracy, conservation) of the numerical scheme. For each cell we consider {P k i } i , the set of polynomials of degree at most k on I i and null outside of I i . The test functions g are then taken as basis functions of P k i . Most of the time those are Legendre or Lagrange polynomials on I i , but depending on the situation, other polynomial basis may also be adequate. According to the philosophy of these methods, we do a L 2 -projection of f on the underlying polynomial space

f h (x, t) = i k l=0 f l,i (t)p l,i (x) (12)
where {p l,i } l=0..k,i is a basis of P k i . The major advantage of this formulation is that it is fully independent of the spatial order, which can be considered as a parameter. and According to the Taylor series representation of the distribution function, the error we commit with the polynomial projection is proportional to ∆x k+1 f (k+1) , where f (k+1) is the k + 1-th derivative of f . Let us now introduce the notation used in the next sections. We name 1 2 ] the intervals respectively along the x axis and along the v axis. The two-dimensional cell T i,j = I i × J j is the tensorial product of I i in direction x and J j in direction v as depicted on Fig. 1. The Gauss-Lobatto nodes are noted x ig (respectively v jg ) in direction x (respectively v) and the associated weight is w ig (respectively w jg ). Gauss weights are always expressed for a quadrature formula on the interval [-1, 1]. P k i is the space of polynomials of degree at most K on cell I i . In the two-dimensional case the polynomial basis also is 2D and is obtained as tensorial product of the 1D basis. We use the notation P k,i,j (x, v) = p k1,i (x)p k2,j (v) with k = (k 1 , k 2 ). For all our work we used Legendre polynomials of degree up to k as basis of P k ([-1, 1]) and a L 2 -normalisation. Superscript n is used to indicate time t n = n∆t and f n h (x) = f h (x, t n ). About characteristics, we call x the variable at time t n+1 and we call x * the foot of the characteristic at time t n .

I i = [x i-1 2 , x i+ 1 2 ] and J j = [v j-1 2 , v j+

The semi-Lagrangian discontinuous Galerkin formulation

This formulation is taken from [START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov Poisson system[END_REF]. We only adapt it to obtain a formulation in a two-dimension phase-space and abbreviate it SLDG. The main idea is to perform a time-integration of ( 11) and to transform it into a space-integration using characteristics. We start from [START_REF] Cheng | Numerical study of Vlasov-Poisson equations for innite homogeneous stellar systems[END_REF] and proceed to a time integration to obtain

Ii f n+1 (x)g dx = Ii f n (x)g(x) dx + t n+1 t n Ii af (x, t)g ′ (x) dx dt - t n+1 t n (af (x, t)g(x)) x - i+ 1 2 -(af (x, t)g(x)) x + i-1 2 dt, ∀g ∈ G. ( 13
)
t n+1 t n x i -1 / 2 x i -1 / 2 * Ω i-1 / 2 ⃗ n= ( n t n x ) = ( 0 1 ) ⃗ n= ( n t n x ) = ( -1 0 )
Figure 2: Characteristic scheme for SLDG method.

We respectively call the three integrals on the right hand side T 0 , T 1 and T 2 . The time integral in T 1 and T 2 are evaluated in a semi-Lagrangian fashion, We start from Fig. 2. On this gure, point x * i-1 2 is the foot of characteristic arriving at x i-1 2 at time t n+1 . Domain Ω i-1 2 is the region bounded by the three points

(x i-1 2 , t n+1 ), (x i-1 2 , t n ) and (x * i-1 2 , t n ).
We integrate (9) over the domain Ω i-1 2 and apply the divergence theorem:

0 = Ω i-1 2 ∂ t f (x, t) + ∂ x (a(x, t)f (x, t)) dx dt = ∂Ω i-1 2 f (x, t)n t + af (x, t)n x ds. (14)
Therefore,

x i-1 2 x * i-1 2 f (x, t n ) dx = t n+1 t n af (x i-1 2 , t) dt. ( 15
)
Hence it is possible to write

T 2 = g(x + i-1 2 ) x i-1 2 x * i-1 2 f n h (x) dx -g(x - i+ 1 2 ) x i+ 1 2 x * i+ 1 2 f n h (x) dx. (16)
By rst performing spatial integration using Gaussian quadrature rule on T 1 and following the same process, the semi-Lagrangian formulation gives

T 1 = t n+1 t n Ii af h (x, t)g ′ (x) dx dt = t n+1 t n ∆x i 2 ig af h (x ig , t)g ′ (x ig )w ig dt = ∆x i 2 ig w ig g ′ (x ig ) xi g x * ig f n h (x) dx. (17)
Foot of characteristics are obtained by solving backward the ODE

⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ dx dt = a(x, t), x(t n+1 ) = x i , (18)
between t n and t n+1 . As a consequence, [START_REF] Cheng | Numerical study of Vlasov-Poisson equations for innite homogeneous stellar systems[END_REF] becomes

Ii f n+1 h (x)g(x) dx = Ii f n h (x)g(x) dx + ∆x i 2 ig w ig g ′ (x ig ) xi g x * ig f n h (x) dx + g(x + i-1 2 ) x i-1 2 x * i-1 2 f n h (x) dx -g(x - i+ 1 2 )
x i+ 1 2

x * i+ 1 2

f n h (x) dx. (19)
One can clearly see that if g is chosen as a basis function of P k i , then the left hand side is simply f n+1 k,i and the rst line of the right hand side is f n k,i . The other terms are not simplied.

In the two-dimensional case, integration over cell T i,j requires an additional integration in the orthogonal direction. Polynomials are only dened on a reference element [-1, 1] and all variables and functions taken on this reference element are indicated with a tilde. The Gauss weights are not included in this statement since they are only used on this reference element. We dene k = (k 1 , k 2 ) such that k is the index for polynomials and coecients in twodimension and, k 1 is the corresponding index in direction x and k 2 the corresponding index in direction v. This leads to the integration in the x direction

f n+1 k,Ti,j = f n k,Ti,j + 2 ∆x i jg w jg ⎛ ⎝ ig w ig xi g x * ig f n h (x, v jg ) dx p′ k1 (x ig ) ⎞ ⎠ pk2 (ṽ jg ) + 2 ∆x i jg w jg ⎛ ⎝ x i-1 2 x * i-1 2 f n h (x, v jg ) dx pk1 (x + i-1 2 )- x i+ 1 2 x * i+ 1 2 f n h (x, v jg ) dx pk1 (x - i+ 1 2 ) ⎞ ⎠ pk2 (ṽ jg ), (20)
with

f n k,Ti,j = ∆x i ∆v j 4 f n k,Ti,j , (21)
Figure 3: Characteristic scheme for CDG. and in the v direction

f n+1 k,Ti,j = f n k,Ti,j + 2 ∆v j ig w ig ⎛ ⎝ jg w jg vj g v * jg f n h (x ig , v) dv p′ k2 (ṽ jg ) ⎞ ⎠ pk1 (x ig ) + 2 ∆v j ig w ig ⎛ ⎝ v j-1 2 v * j-1 2 f n h (x ig , v) dv pk2 (ṽ + j-1 2 )- v j+ 1 2 v * j+ 1 2 f n h (x ig , v) dv pk2 (ṽ - j+ 1 2 ) ⎞ ⎠ pk1 (x ig ). (22)

The characteristic-discontinuous-Galerkin formulation

We take the algorithm for this method from [START_REF] Guo | A conservative semi-Lagrangian discontinuous Galerkin scheme on the cubed-sphere[END_REF] and we call it characteristicdiscontinuous Galerkin and abbreviate it CDG. Nevertheless it has previously been called Lagrange-characteristic in [START_REF] Childs | Characteristic Galerkin Methods for Scalar Conservation Laws in One Dimension[END_REF] or Lagrange projection even sooner in [START_REF] Douglas | Numerical Methods for Convection-Dominated Diusion Problems Based on Combining the Method of Characteristics with Finite Element or Finite Dierence Procedures[END_REF] and [START_REF] Pironneau | On the transport-Diusion Algorithm and Its Application to the Navier-Stockes Equations[END_REF]. It also appeared in the framework of electron hole and drift diusion for semi conductors [START_REF] Douglas | Simulation of the transient behavior of a one dimensional semi-conductor device[END_REF][START_REF] Gamba | Simulation of the transient behavior of a one dimensional semiconductor device II[END_REF]. As for the SLDG formulation, we only adapt it to our two-dimensional problem. The starting point of the method is to determine [START_REF] Colombi | Vlasov versus Nbody: the Hénon sphere[END_REF], as presented on the left part of Fig. 3. The conservation equation [START_REF] Besse | Semi-Lagrangian schemes for the Vlasov equation of an unstructured mesh of phase space[END_REF] implies that

I * i = [x * i-1 2 , x * i+ 1 2 ] using
f n+1 h (x) = f n h (x * ) ∂x * ∂x . (23)
Multiplying [START_REF] Feix | A Universal Model: The Vlasov Equation[END_REF] by g(x), integrating the result in the variable x on cell I i , and using the change of variable from x to x * in the right hand side of the resulting equality, we then obtain

Ii f n+1 h (x)g(x) dx = I * i f n h (x * )g(x(x * )) dx * . ( 24
)
Note that in the right hand side, the integration variable is x * but test function g depends on variable x. These two variables are linked by [START_REF] Colombi | Vlasov versus Nbody: the Hénon sphere[END_REF]. Variable x * is at time t n while x is at time t n+1 . This second integral is sliced on each sub-interval I * i ∩ I l , and on every sub-interval the integral is computed using Gaussian quadrature rule so that [START_REF] Frénod | Long-time simulation of a highly oscillatory Vlasov equation with an exponential integrator[END_REF] becomes

I * i f n h (x * )g(x) dx * = l I * i ∩I l f n+1 h (x * )g(x(x * )) dx * = l ⎛ ⎝ ig w ig f n h (x * i,l,ig )g(x(x * i,l,ig )) Γ(I * i,l ) 2 ⎞ ⎠ , (25)
with Γ(I * i,l ) the length of I * i ∩ I l . Here we have

I * i = [x * i-1 2 , x * i+ 1 2 ], where x * i,l,ig
are foot of characteristics of Gauss points on I * i ∩ I l and x i,l,ig are the associated end of characteristics on I i . This process is illustrated on the right part of Fig. 3. Here again, taking g = p k,i , the left hand side equals to f n+1 k,i . Performing just the same as previously we integrate on the orthogonal direction over J j to obtain the formulation in direction x:

f n+1 k,Ti,j = jg w jg l ig w ig f n h (x * ig,l,jg , v jg )p k,l (x ig,l,jg , v jg ) Γ(I * i,l,jg ) ∆x i . (26)
Note that here, because the foot of characteristics depends on the orthogonal direction, Γ also depends on the orthogonal direction. We do the same in direction v:

f n+1 k,Ti,j = ig w ig l jg w jg f n h (x ig , v * jg,l,ig )p k,l (x ig , v jg,l,ig ) Γ(J * j,l,ig ) ∆v j . ( 27
)
and The following three sections describe how we modify the SLDG and CDG schemes to transform them into the adaptive AMW-SLDG and AMW-CDG schemes. These adaptive schemes rely on a multi-wavelet representation of the data (Sec. 4) and some special considerations concerning numerical evaluation of the integrals involved in the schemes (Sec. 5). These schemes also rely on a prediction procedure to achieve adaptive mesh renement (Sec. 6).

Multi-wavelets representation

The multi-wavelets are a tool to build a multiscale representation of any function belonging to a suitable functional space, typically L 2 . Their properties in one dimension are detailed in [START_REF] Alpert | Adaptive Solution of Partial Dierential Equations in Multiwavelet Bases[END_REF]. It then has been widely used, such as [START_REF] Hovhannisyan | Adaptive multiresolution discontinuous Schemes for conservation laws[END_REF] and [START_REF] Gerhard | Adaptive multiresolution discontinuous Galerkin Schemes for conservation laws: multi-dimensional case[END_REF]. An algorithm for the construction of one-dimensional multi-wavelets can be found in [START_REF] Alpert | A Class of Bases in L 2 for the Sparse Representation of Integral Operator[END_REF]. Let us remind the most important points from [START_REF] Alpert | Adaptive Solution of Partial Dierential Equations in Multiwavelet Bases[END_REF] about multi-wavelets.

Polynomial multi-wavelets form a multiscale basis of L 2 ([0, 1]). Let us consider a nested sequence of spaces

V k n = f [2 -n l ; 2 -n (l+1)] ∈ P k , l ∈ ⟦0 ; 2 n -1⟧ ,
where P k is the set of polynomials of degree at most k. The dimension of V k n is 2 n (k + 1) and

V k 0 ⊂ V k 1 ⊂ ⋯ ⊂ V k n ⊂ ⋯ ⊂ L 2 ([0, 1]), (28)
such that

clos L 2 ([0,1]) ∞ ⋃ n=0 V k n = L 2 ([0, 1]). (29)
Given an orthogonal basis {φ j } of V k 0 , it is possible to obtain an orthogonal basis {φ n j,l (x)} of V k n by dilatation and translation such that

φ n j,l (x) = 2 n 2 φ j (2 n x -l), j ∈ ⟦0 ; k⟧, l ∈ ⟦0 ; 2 n -1⟧, n ∈ N. (30)
It is always possible to nd a space W k n such that

V k n ⊕ W k n = V k n+1 , V k n W k n , (31)
with W k n of dimension 2 n (k + 1). As a consequence of (31), we have

V k n = V k 0 n-1 ⊕ η=0 W k η . (32)
Given an orthogonal basis {ψ j } of W k 0 , following [START_REF] Heath | Analysis of the Discontinuous Galerkin Method Applied to Collisionless Plasma Physics[END_REF], the piecewise polynomial basis functions of W k n are {ψ n j,l (x)} such that

ψ n j,l (x) = 2 n 2 ψ j (2 n x -l), j ∈ ⟦0 ; k⟧, l ∈ ⟦0 ; 2 n -1⟧. ( 33
)
Functions φ n j,l are called scaling functions while functions ψ n j,l are called multiwavelets functions. In our case, scaling functions are Legendre polynomial rescaled to the corresponding cell.

Once {φ j } and {ψ j } are known, it is then possible to build lters h 0 i,j , h 1 i,j , g 0 i,j and g 1 i,j according to [START_REF] Alpert | Adaptive Solution of Partial Dierential Equations in Multiwavelet Bases[END_REF] such that we have a restriction operator

φ i (x) = √ 2 k j=0 h (0) i,j φ j (2x) + h (1) i,j φ j (2x -1) , i ∈ ⟦0, k⟧, (34a)
ψ i (x) = √ 2 k j=0 g (0) i,j φ j (2x) + g (1) i,j φ j (2x -1) , i ∈ ⟦0, k⟧, (34b)
and an interpolation operator

φ i (2x) = 1 √ 2 k j=0 h (0) j,i φ j (x) + g (0) j,i ψ j (x) , i ∈ ⟦0, k⟧, (35a)
φ i (2x -1) = 1 √ 2 k j=0 h (1) j,i φ j (x) + g (1) j,i ψ j (x) , i ∈ ⟦0, k⟧. (35b)
Note that if you consider only polynomials of degree 0, you obtain Haar wavelets.

Starting from [START_REF] Heath | A discontinuous Galerkin method for the Vlasov-Poisson system[END_REF] it is possible to get the multiscale decomposition P n f of f at level m on scaling functions basis

P n f (x) = k j=0 s 0 j,0 φ j (x) + n-1 m=0 2 m -1 l=0 k j=0 d m j,l ψ m j,l (x), (36)
where coecients {s m j,l } and {d m j,l } are obtained from a L 2 -projection on the corresponding space,

s n j,l = 2 -n (l+1) 2 -n l f (x)φ n j,l (x) dx, (37a)
d n j,l = 2 -n (l+1) 2 -n l f (x)ψ n j,l (x) dx. (37b)
Filters can be built for any polynomial basis. Alpert et al. [START_REF] Alpert | Adaptive Solution of Partial Dierential Equations in Multiwavelet Bases[END_REF] present their construction for Legendre polynomial basis and Lagrange polynomial basis. In our code we use only Legendre polynomial basis.

According to (34) [START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov Poisson system[END_REF] it is possible to compute level m from level m + 1, using

s m i,l = k j=0 h (0) i,j s m+1 j,2l + h (1) i,j s m+1 j,2l+1 , (38a)
d m i,l = k j=0 g (0) i,j s m+1 j,2l + g (1) i,j s m+1 j,2l+1 , (38b)
and level m + 1 from level m, using

s m+1 i,2l = k j=0 h (0) j,i s m j,l + g (0) j,i d m j,l , (39a)
s m+1 i,2l+1 = k j=0 h (1)
j,i s m j,l + g (1) j,i d m j,l . (39b)
Formulae [START_REF] Shelton | A multi-resolution discontinuous Galerkin method for unsteady compressible ows[END_REF] (39) can be extended to two-dimensional case using tensorial products. Therefore space W k n becomes

W k n = {φ n i,l (x)ψ n j,l (y)}, {ψ n i,l (x)φ n j,l (y)}, {ψ n i,l (x)ψ n j,l (y)}, i, j ∈ ⟦0 ; k⟧, l ∈ ⟦0 ; 2 n -1⟧ . ( 40
)
In order to perform restriction or interpolation in the two dimensional case the restriction or interpolation operators are applied in the rst direction then in the second direction. This means that there are two equivalent ways to consider f h (x, v), the distribution function on a given cell, by using only coecients of levels m and m + 1. The rst way is to consider only the level m + 1 In this case, q is the index over the two dimensional basis functions and c is the index over the children cells. The second way is to consider the distribution function at level m plus the details

f h (x, v) = c=0,1,2,3 q φ m+1 q,[c] (x, v)s m+1 q,[c] . (41)
f h (x, v) = q φ m q (x, v)s m q + Θ=a,b,c k ψ m k,[Θ] (x, v)d m k,[Θ] (42)
where a, b and c respectively are the subspaces of W k n generated by {φ n i,l (x)ψ n j,l (y)}, {ψ n i,l (x)φ n j,l (y)} and {ψ n i,l (x)ψ n j,l (y)}. The thresholding operation consists in comparing the l 2 -norm of the details {d m ., [.] } to a given threshold and ignore it where it is smaller than this threshold, that is,

if ⎛ ⎝ Θ=a,b,c k d m k,[Θ] 2⎞ ⎠ 1 2 ≤ 0 , then ignore details of level m, (43)
and repeat this operation until all details kept are assessed signicant, proceeding from highest level to lowest level. This thresholding criterion is presented for l p -norm of details in [START_REF] Besse | A wavelet-MRA-based adaptive semi-Lagrangian method for the relativistic Vlasov-Maxwell system[END_REF]; here we use only the l 2 -norm. Parameter 0 and maximum level of renement are chosen empirically to have a good description of initial distribution function, with the possibility to follow details during simulation, but still big enough to have reasonable number of points and simulation time. We refer to [START_REF] Besse | A wavelet-MRA-based adaptive semi-Lagrangian method for the relativistic Vlasov-Maxwell system[END_REF] for a discussion concerning thresholding, adaptivity and optimality. and , Multi-wavelets are built and used with respect to the L 2 -norm, which means that part of the L 2 -norm depends on the multi-wavelet coecients and then the L 2 -norm is not conserved when coarsening. Nevertheless, contrary to some usual interpolation schemes, the interest and the power of the multi-wavelets representation is not only to control but also to minimise the loss of L 2 -norm by choosing locally in an optimal way the details which deserve to be conserved or discarded. Such property is conrmed by the numerical experiment of Sec. 7. Equations ( 41) and (42) also ensure that moments of the distribution function f are conserved during coarsening up to the degree of scaling functions, thanks to the orthogonality properties of the multi-wavelet functions ψ n j,l . This is conrmed with numerical experiment. We consider the distribution function For this distribution function, we consider uniform meshes with respectively 7, 6 and 5 level of renement and use polynomials of degree up to two for each direction. We rst express the distribution function on uniform mesh of level 7 and compute its rst momentum in velocity vf (x, v), that is

f (x, v) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 0 if v < 0, v 2 √ 2π (1 + 0.05 cos(0.5x)) exp -v 2 2 else, (x, v) ∈ [0, 4π] × [-6, 6]. (44) T 2 T 1 T 1 * x i -/ 3 2 x i -/ 1 2 x i + / 1 2 v j+ / 1 2 v j-/ 1 2
4π 0 6 -6 vf (x, v) dv dx. (45)
We then compute the distribution function on uniform mesh of level 6 following (38a) and compute the momentum. This means we only compute the scaling coecients s m i,l and eliminate the multi-wavelet coecients d m i,l . We do the same to get the distribution function on uniform mesh of level 5 and compute the momentum. The dierence of momentum computed is displayed in Table 1. This dierence is of the magnitude of machine error which validate our code. and

Error analysis

In this section we identify some sources of numerical error which are responsible for a loss of some conservation laws. Therefore, by explaining their origins and mechanisms we could nd a way to x them. In addition, event if the conservation laws are not exactly preserved a priori, such conservation laws are preserved asymptotically thanks to the convergence property of the AMW-SLDG and AMW-CDG schemes, which is numerically shown in Secs. 7.1 and 7.2.

Error linked to integration on transported cells

In this section we explain one of the limit and error source of our schemes. We present it on a uniform mesh since it is sucient to understand the source of error and it is much easier than on non-conformal mesh.

Let us start with the AMW-CDG method on uniform mesh and the equation

∂f ∂t (x, v) + g(v) ∂f ∂x (x, v) = 0, (46)
where g(v) is a polynomial of degree m and m smaller than or equal to K the degree of our Legendre polynomials. In our algorithm we then use K + 1 Gauss points. We consider time t n . On

T * 1 ∩ T 2 =[x i-1 2 , x * i-1 2 ] × [v j-1 2 , v j+ 1 2 ] = ⋃ v∈[v j-1 2 ,v j+ 1 2 ] [x i-1 2 , P(v)]. (47)
with

x * i-1 2 = P(v) = x i-1 2 + g(v)dt. (48)
Consequently, P ∈ P m . According to the method we then must solve

T * 1 ∩T2 f n h (x * , v)P k1,k2 (x(x * ), v) dx * dv = v j+ 1 2 v j-1 2 P(v) x i-1 2 f n h (x * , v)P k1,k2 (x(x * , v), v) dx * dv. (49)
According to the resolution of characteristics equations ( 18), one has

x(x * , v) = x * -g(v)dt. (50)
Here, x is a polynomial of degree one in the variable x * and a polynomial of degree m in the variable v.

We proceed to the following change of variable

x * = a -P(v) a -b x + a P(v) -b a -b , (51)
in order to integrate over rectangle

(x, v) ∈ [a, b] × [v j-1 2 , v j+ 1 2 ]
. Formula (47) then becomes

T * 1 ∩T2 f n h (x * , v)P k1,k2 (x(x * ), v) dx * dv = v2 v1 b a f n h a -P(v) a -b x + a P(v) -b a -b , v P k1,k2 x a -P(v) a -b x + a P(v) -b a -b , v a -P(v) a -b dx dv. ( 52
)
We must then search the maximum degree of

f n h a -P(v) a -b x + a P(v) -b a -b , v P k1,k2 x a -P(v) a -b x + a P(v) -b a -b , v a -P(v) a -b .
Let us write

f n h (x * , v)P k1,k2 (x(x * , v), v) =f n h (x * , v)P k1,k2 (x * -g(v)dt, v) =f n h a -P(v) a -b x + a P(v) -b a -b , v P k1,k2 x a -P(v) a -b x + a P(v) -b a -b -g(v)dt , v . (53)
The distribution function f n h is a polynomial of degree K in x and in v. Then, considering only monomials of highest degree, we obtain

(v m x) K v K f n h (v m x) K v K P K,K v m Jacobian of (51) = x 2K v 2K+2Km+m . ( 54
)
Using K + 1 Gauss points we can perform exact integration only for polynomials of degree up to 2K+1. Consequently, according to (54), integral (52) is computed exactly in direction x but not in direction v if m is larger than or equal to 1.

Then, as we implemented and use our AMW-CDG scheme, it does not conserve mass.

For the AMW-SLDG method, starting from (46) and term T 1 from (13), we get

T 1 = Jj ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ p k2 (v) Ii ⎛ ⎜ ⎜ ⎝ pk1 (x(x)) x x-g(v)dt x * f n h (ξ, v) dξ ⎞ ⎟ ⎟ ⎠ dx ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ dv = Ii pk1 (x(x)) Jj p k2 (v) x x-g(v)dt f n h (ξ, v) dξ dv dx. (55)
In the second line, for a given x, the problem is equivalent to the one presented above for the AMW-CDG scheme. This is also true for term T 2 in (13). Consequently, as for our AMW-CDG scheme, our implementation of the method AMW-SLDG does not conserve mass a priori. Nevertheless, in the numerical experiments of Sec. 7, we observe a rather good conservation of mass a posteriori. Such result is consistent with the fact that the schemes are convergent (see Secs. 7.1 and 7.2).

In both case (AMW-CDG and AMW-SLDG) we see that an exact integration would be a rst step to a priori mass conservation. An exact integration (considering a possibly piecewise polynomial P(v)) would require to compute the equation of the transported border. Since in our case, function g(v) is always known (it is either v er E n h (x) depending on the direction) and a polynomial, the transported border can always be computed. It is possible to integrate exactly polynomials on domains with polynomial boundaries but the computation cost increases. In addition, the fact that the boundary may be piecewise polynomial if g = E h also increases the complexity. However, this last point can be avoided if one consider level of renement for the electric eld E h coarse enough.

Error linked to the a priori lack of ux compensation

We perform most of our integration using quadrature formula with Gauss points. For a one-dimensional mesh (uniform or not) and for conformal ndimensional mesh, n ≥ 2, this formula is exact because the quadrature points are always adapted to each segment, as shown on Fig. 5a. In the case of nonconformal mesh, we use a quadrature formula which is exact for given polynomials. However, our distribution function is not a single polynomial on the domain, as depicted on Fig. 5b. It is piecewise polynomial. Therefore, integration performed with Gauss quadrature does not give the exact value. This particular case happen for the integration in direction v in equations ( 20) and [START_REF] Gamba | Simulation of the transient behavior of a one dimensional semiconductor device II[END_REF] and for the integration in direction x in equations ( 22) and [START_REF] Gerhard | A High-Order Discontinuous Galerkin Discretization with Multiwavelet-Based Grid Adaptation for Compressible Flows[END_REF].

In addition, in the case of non-conformal mesh, uxes are not computed at the same points for adjacent cells. This is illustrated on Fig. 6. The blue points are quadrature points for the big cell on the left. The red points are quadrature points for the small cell on the right. We clearly see that blue points do not match red points. Consequently, the distribution function, and thus uxes, computed at the blue points do not have the same values as the distribution function and uxes computed at the red points. Therefore there is no guaranty of a priori ux conservation. In a sense, to ensure a priori ux conservation, the computation of uxes must be independent of the scale of the interface. Therefore we can imagine a process which assigns to each interface a unique ux value from the ux values computed at the dierent scales considered at this interface.

Algorithm and numerical implementation

This section is devoted to the description of our algorithm and to a few implementation points. Our algorithm is as follows.

Step 1: initialisation. We project the distribution function f on the coarsest level of mesh (usually level 0) and compute the details at same level. If details are large enough, the cell is rened. We repeat this operation until description of f is precise enough according to (43). The thresholding used for renement is the same that will be used for coarsening on future steps.

Step 2: rst step of Strang splitting. The rst step of the Strang splitting is a half time step in direction x only (from t n to t n + dt 2). Each step of our Strang splitting can be divided into three stages.

(a) Prediction of the transport of the grey cell. The initial mesh is black. The considered cell is grey. Red and green arrows show the displacement of the corner of the cell and blue is the transported cell.

(b) Creation of the predicted mesh. In blue is the transported cell and in black are all the newly created cells. Prediction. We rst predict the mesh by using a time explicit Eulerian scheme to compute characteristics forward and to ensure that the cell containing end point of characteristics are at least one level ner than their departure point. This prediction will enable the mesh to follow the small structures and anticipate the creation of smaller ones. Fig. 7 is an illustration of this process. Fig. 7a represents the initial mesh M n . Fig. 7b represents the mesh obtained by transporting the initial mesh M n . To do so we use an explicit Euler scheme to transport the boundaries of each cell of Fig. 7a. We then make sure that all the area described by the transported cell has the expected level of renement. This process is illustrated in Fig. 8. Each corner (x, v) is transported to obtain a new corner (x + v dt, v). Then we rene the mesh such that the area described by the transported corners is at least of the same level of renement. Mesh of Fig. 7c is obtained by adding one level of renement to the mesh of Fig. 7b to anticipate creation of smaller structures during the computation of the time step. This choice of prediction is taken from [START_REF] Besse | A wavelet-MRA-based adaptive semi-Lagrangian method for the relativistic Vlasov-Maxwell system[END_REF]. While they add a criterion to chose whether to add or not one (or more than one) level of renement, we do not have this criterion and always make the mesh one level thinner. We refer to [START_REF] Besse | A wavelet-MRA-based adaptive semi-Lagrangian method for the relativistic Vlasov-Maxwell system[END_REF] for a discussion of accuracy and optimality during the prediction step. Mesh M n+1 of Fig. 7d is the merging of initial mesh M n (Fig. 7a) and the predicted mesh of Fig. 7c. The prediction of the mesh is done using an independent structure that contains the geometry but not the data. We manipulate two trees. The one on Figs. 7a and 7d contains the distribution function while the one presented on Figs. 7b and 7c contains only the geometry and its tree data. The mesh M n+1 (for computation) must be at least as locally rened as the initial mesh M n , otherwise signicant details of the distribution function f n h would be lost during the next computation step. The distribution function f n h known on the mesh M n (Fig. 7a) is projected on the mesh M n+1 (Fig. 7d) with all unknown thinner details set to zero. From [START_REF] Gerhard | Adaptive multiresolution discontinuous Galerkin Schemes for conservation laws: multi-dimensional case[END_REF] and thanks to the merging of the predicted mesh with the initial mesh, there is neither loss of information nor of moments conservation during the prediction step (for moments of order lower than the degree of scaling functions).

Computation. The chosen numerical scheme (SLDG or CDG) is applied on the nal mesh to update the distribution function. We solve

∂ t f + ∂ x (vf ) = 0 ( 56
)
between t n and t n+ 1 2 , as presented in Secs. 3.2 and 3.3.

Coarsening. We perform recursive coarsening operation according to (43) in order to keep only the nest signicant details.

Step 3: resolution of the Poisson equation. We solve the Poisson equation using the distribution function obtained after Step 2. The eld is always computed on the nest level of renement as its numerical cost is very small compared to the phase-space numerical cost.

Step 4: second step of Strang splitting. The second step of Strang splitting is a full time step in direction v only (from t n to t n+1 ). We solve

∂ t f + ∂ v (Ef ) = 0. ( 57
)
This step follows the scheme of Step 2 and uses data from Steps 2 and 3.

Step The mesh is organised as a tree structure, but unlike many tree structures we do not concentrate data on leaf cells. In our case, we remark that for computation of data of level m, considering or not cells of level higher than m where they were available, had absolutely no impact on the solution. This means that coarse cells are not aected by ne details during the time integration process. Therefore, storing the distribution function at every level according to (41) and using only the coarsest signicant level (m where available or the coarser level if level m is not available) reduces the number of computation needed. This operation requires to compute the distribution function on each cell of the tree at the end of steps 2, 4 and 5 (at the end of each step the distribution function is known only on leaves). The computational impact is signicant only when the mesh becomes highly heterogeneous in terms of local density of cell.

Here is another optimisation we use. Let p k be the k-th normalised Legendre polynomial. Then we have

1 -1 p k (x) dx = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ √ 2 2 if k = 0, 0 else. ( 58
)
This is used in the computation of (20) [START_REF] Einkemmer | Convergence Analysis of Strang Splitting for Vlasov-Type Equations[END_REF] to reduce the computational cost where the integration domain includes full cells with

Ii f (x) dx = Ii K k=0 f k p k (x) dx =f 0 Ii p 0 (x) dx =f 0 √ 2 2. ( 59
)
This simplication is particularly ecient when the time step is large compared to the size of cell.

7. Numerical results

Linear transport: rotation

We consider the rotation problem

∂ t f + v∂ x f -x∂ v f = 0, (x, v) ∈ [-10, 10] 2 , (60a) f (x, v, 0) = exp -(x -5) 2 -v 2 . ( 60b
)
The exact solution in time is

f (x, v, t) = exp -(x * cos(t) + v * sin(t) -5.) 2 -(v * cos(t) -x * sin(t)) 2 .
(61) we use ∆t = 10 -5 and 1, 000 time steps, which makes nal time T f = 10 -2 . Our goal here is to check the order of convergence to zero of the L 2 -error in phasespace. We perform simulations with polynomials of degree 2 and 3 as indicated in Table 2. The maximum level of renement is 8 but is not reached in these simulations. The threshold criteria are given in Table 2. For each threshold value, the average phase-space discretisation step ⟨h⟩ is the average length of cell side over phase-space and time, that is

⟨h⟩ = L x L v (N t + 1) ∑ Nt n=0 C n 1 2 , ( 62
)
with L x and L v respectively the length of the domain in direction x and v, N t the number of time steps, and C n the number of cells at time t n . For polynomials of degree k we expect the error to be of order k + 1. We observe in Table 2 that the rate of convergence is between k + 1 and k + 2 for k = 2, 3 and for both schemes.

Burgers equation

We consider the Burgers' equation

∂ t f (x, v, t) + ∂ x (f (x, v, t) 2 ) 2 = 0, (x, v) ∈ [0, 2π] × [-1, 1], (63a)
f (x, v, 0) = 1 + sin(x). (63b)
We perform simulations with ∆t = 10 -5 and 1, 000 time steps for a nal time T f = 10 -2 . The threshold criteria are given in Table 3. For each threshold value, average size ⟨h⟩ is the average length of cell side over space and time (62). For this test case, we do not enforce the maximum level of renement. These simulations are done either with polynomials of degree two or three, as indicated in captions of Table 3. For polynomials of degree k, the order of the error is expected to be k + 1. In the Table 3 we observe a convergence rate between k + 1 2 and k +1. Even if we do not always obtain the optimal convergence rate k +1, we must point out that, in nonlinear cases, convergence rate may crucially depend on the properties of the reconstruction operator (here multi-wavelets) and on the geometric features of the mesh (here a two-dimensional non-conformal adaptive mesh). Here, in order to understand better this behaviour, a careful and rigorous mathematical analysis of the schemes for the Burgers equations should be done. Nevertheless it is beyond the scope of this paper. Moreover the nature and behaviour of the Burgers nonlinearity is quite dierent from the quadratic nonlinearity appearing in Vlasov-Poisson equations. Therefore, in this paper, we do not pursue further with a ner analysis.

Swirling deformation

We consider solving

∂ t f -∂ x (g(x, v, t)f ) + ∂ v (g(v, x, t)f ) = 0, x ∈ [-π; π], v ∈ [-π; π], ( 64
) with g(x, v, t) = π cos 2 x 2 sin(v) cos πt T , (65)
and T = 1.5. This test case is the example 5.5 of [START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov Poisson system[END_REF].

The initial conditions are

f 0 (x, v) = ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 1 if x 2 + (v -π 2) 2 ≤ 0.30π and ( x ≥ π 20 or v ≥ 0.30π), 1 -x 2 + (v + π 2) if x 2 + (v + π 2) ≤ 0.30π, 1 4(1 + cos(π (x + π 2) 2 + v 2 ) if (x + π 2) 2 + v 2 ≤ 0.30π, 0 else.
The domain for this case is [-π, π] 2 . Final time is 1.5 and this simulation is done using 100 time steps for a time step ∆t = 0.015. The maximum level of renement is 8, the level 0 being the full domain, and the thresholding parameter is 0 = 0.01. Our polynomials are of degree up to 3 per direction.

The initial distribution function and the initial mesh are presented on Fig. 9. For convenience, oscillations at boundaries of the truncated cylinder have been truncated on this plot. On Fig. 9b, each cross indicates a cell corner. It is obvious that our mesh is non-conformal, and according to equations (20) ( 22) and ( 26)( 27), conform mesh is not required. Although the mesh may seem sparse, the distribution function on each cell is described using several polynomials. For polynomials of degree up to k, there are (k + 1) 2 polynomials per cell. Consequently, in this case, the distribution function is described using 16 polynomials per cell. Fig. 10 presents the distribution function and mesh of the swirling deformation during the deformation and at the nal time for the AMW-SLDG scheme. In both case the mesh follows very well details of the distribution function, especially the discontinuity of the cylinder. One can see that its sharp border matches with high mesh renement along time, which is the expected behaviour. The results of the two methods (AMW-SLDG and AMW-CDG) also look extremely similar during simulation and at nal time. Because for all our test cases distribution functions and meshes obtained with both methods always give results that can hardly be distinguished, we only present them for one method.

One can clearly see that the nal mesh is very dierent from the initial mesh, especially at the boundary of the truncated cone. Initially this part requires many points because the distribution function is not C 1 . However, the use of a polynomial basis smooths the distribution function on each cell. After just a few time steps, the distribution function acquires higher local regularity than it initially had.

Since this case projects a strongly discontinuous function on a polynomial space, the initial distribution function presents strong articial oscillations and it can hardly be considered good test case for conservation properties. Nevertheless, time evolution of numerical mass, L 1 -norm and L 2 -norm for this test case are presented on Fig. 11 and show good conservations.

Landau damping

Landau damping is a standard test case to validate codes for Vlasov-Poisson in two-dimensional phase-space for plasma physics. The initial distribution (a) Swirling at t=0.75. - function is a Gaussian function in velocity with a perturbation in space.

f 0 (x, v) = (1 -α cos(x 2)) √ 2π exp -v 2 2 , (x, v) ∈ [0, 4π] × [-L v , L v ], (66)
where L v is large enough to ensure Dirichlet boundary condition. Usually L v = 5 is sucient. Two main values are chosen for α. First, α = 0.01 corresponds to the linear Landau damping, also called weak Landau damping. With this small perturbation it is possible to get a linearised approximation on the electric energy along time, and especially the damping rate. Secondly, α = 0.5 corresponds to the non-linear Landau damping, also called strong Landau damping. Although there is no analytical solution for this case, it has been heavily simulated and there are many references to compare the damping rate in rst stage and the growth rate of second stage, such as [START_REF] Ayuso De Dios | Discontinuous Galerkin methods for the one-dimensional Vlasov-Poisson system[END_REF], [START_REF] Rossmanith | A positivity-preserving high-order semi-Lagrangian discontinuous Galerkin scheme for the Vlasov-Poisson equations[END_REF], [START_REF] Heath | A discontinuous Galerkin method for the Vlasov-Poisson system[END_REF] and many others. For both linear and non-linear cases, nal time is T = 100 with 10 3 time steps for a time step ∆t = 0.1. The maximum level of renement is 5 for the weak Landau damping test case and 8 for the strong Landau damping one. We use polynomials of degree up to 2 for each direction of space. Threshold is 0 = 10 -4 for weak Landau damping, while it is 0 = 10 -2 for strong Landau damping.

For the weak Landau damping, at nal time, the relative variation of mass and the relative variation of L 1 -norm are around 2.8 × 10 -7 , the relative loss of L 2 -norm is around 1.7 × 10 -5 , the relative variation of L ∞ -norm is almost -0.7%, the relative loss of total energy is slightly less than 0.4% and momentum remains between -1.5 × 10 -10 and 0. Fig. 12 depicts the evolution of electric eld for the AMW-CDG and AMW-SLDG methods (red and green curves) and the analytical damping (blue curve) for weak Landau damping. Here again the two methods give very similar results. (It is possible to see a dierence between the AMW-CDG and AMW-SLDG methods in the long time behaviour.) The expecting damping rate is also recovered until the details become ner than our thresholding criterion. On long time the electric eld increase again. This phenomenon is clearly due to the recurrence phenomena described in [START_REF] Cheng | The integration of the Vlasov equation in conguration space[END_REF] and [START_REF] Nakamura | Cubic interpolated propagation scheme for solving the hyper-dimensional VlasovPoisson equation in phase space[END_REF]. For uniform mesh, the recurrence time is T r = 2π (k∆v). If we consider the recurrence time computed using a uniform mesh with the maximum number of cells, we have ∆v = 0.3125, which gives us a recurrence time T r = 40.21. On Fig. 12, we observe two recurrence phenomena. The rst one is at time t = 38 and the second one is around time t = 76. The numerical recurrence time T r,num = 38 is in good agreement with the predicted one T r = 40.21, since the relative error is about 5.5% Fig. 13 shows the distribution function at various time for strong Landau damping while Fig. 14 presents the mesh. These are obtained with the AMW-CDG method. The distribution function can not be distinguished from the one obtained with the AMW-SLDG method and meshes are extremely similar. At time t=4 and t=8 it is clearly visible that the mesh is ner where details of the distribution function are more important. At large time, because ne details are present almost everywhere on the domain, the mesh is at its nest level almost everywhere for v ∈ [-4, 4] and looks like a static non-uniform mesh which would be thin for v ∈ [-4, 4] and coarse outside.

Fig. 15 indicates the variation of the main diagnostics with both methods. For this test case, the numerical mass is very well conserved but L 1 -norm varies more, which is conform to the loss of positivity. None of the methods preserve total energy but the damping rate and growing rate of electric energy match the literature values. Even if not perfect, momentum appears to be quite well conserved too. It will be easier to observe momentum time variation on bumpon-tail since initial momentum is then not zero. The large increase of total energy occurs during the rst plasma period (about 12 time steps). cells increases for a much longer time as the renement is local and small details spread among the domain.

Two stream instability

We consider the initial condition

f 0 (x, v) = v 2 √ 2π (1 + 0.05 cos(0.5x)) exp -v 2 2 , (x, v) ∈ [0, 4π] × [-6, 6]. ( 67
)
The nal time is still 100 with 10 3 time steps for a time steps ∆t = 0.1. Maximum level of renement still is 8 and our polynomial still are of degree up to 2. The threshold is 0 = 0.003. This initial condition and the domain for this test case are the same as in [START_REF] Heath | A discontinuous Galerkin method for the Vlasov-Poisson system[END_REF]. Fig. 17 illustrates the distribution function and the mesh at dierent times. Those gures are obtained using the AMW-SLDG scheme but distribution function obtained with the AMW-CDG scheme can not be distinguished and meshes are very similar. Here we can clearly see that the mesh follows details of the distribution function during the vortex creation and on longer time.

Fig. 18 presents evolution of numerical mass, L 1 and L 2 -norms, minimum of the distribution function, total energy and momentum. Positivity is not preserved, but numerical mass is well preserved since its relative variation is of order 10 -5 . Relative variation of energy is at most of 0.7%, and absolute time variation of momentum is less than 10 -5 (starting from null initial momentum). The minimum of the distribution function knows a peak to -0.0045 near t = 25 plasma period. Compared to the initial maximum of the distribution function (about 0.31), this is a bit more than 1.5%. Although it is not presented here, the maximum has a peak at the same time. This may be caused by Runge phenomena where we try to describe some variations in the phase-space that have high amplitude over short time. On Fig. 18, the relative variation of minimum is the variation of the minimum compared with the initial maximum of the distribution function.

Fig. 19 depicts the distribution function as a function of particular energy E(t) = v 2 2 -φ(t) at nal time with both schemes. We see that a stationary state of BGK type has not been reached at the end of our simulation (t = 100) since the distribution function is not a single-valued function of the energy E (especially for small values of E). We also observe that until time t = 100 numerical dissipation is very small (see the decrease of L 2 -norm on Fig. 18, less than 2%). This weak diusion in terms of L 2 -norm reects well the power of multi-wavelets representation in the Lebesgue space L 2 , because, contrary to some usual interpolations schemes, such representation controls and minimises the loss of L 2 -norm by choosing locally in an optimal way the details which deserve to be conserved on discarded. To see the eect of numerical dissipation we look at the behaviour of the solution on longer time. Using the same time step ∆t = 0.1 and a coarser description in space (polynomials of degree two, 8 levels of renement, 0 = 0.03) we perform some simulation up to t = 210. Fig. 20 shows the distribution function and its contour plot at nal time obtained with the AMW-SLDG scheme. The results obtained with the AMW-CDG scheme are extremely similar. On this gure we can see the impact of the coarsened mesh on the smoothed distribution function. The corresponding plots of the distribution function as a function of particular energy are on Fig. 21. At this point we clearly have a single-valued function of the energy E. Therefore it seems that at time t = 210, the solution has reached a stationary state of BGK type. Moreover, beyond time t = 100 and before time t = 210 we observe greater numerical diusion, almost 4%, as presented on Fig. 22. Our conjecture is that the numerical diusion might play the part of the physical mechanism which promotes the relaxation of the solution towards a stationary state of BGK type. Actually we can not say if the convergence towards such stationary state of BKG type comes from a true physical relaxation process or if it is a numerical diusion eect, because at this stage of the simulation, numerical diusion in term of L 2 -norm is not negligible.

Bump on tail

We consider the initial condition

f 0 (x, v) = (1 + 0.04 cos(0.5x)) 10 √ 2π 9 exp -v 2 2 + 2 exp -2(v -4.5) 2 , (x, v) ∈ [0, 20π] × [-9, 9]. ( 68
)
Final time is 300 with 3 × 10 3 time steps for a time step ∆t = 0.1. Maximum level of renement is 8 and polynomial degree is 2. Threshold is 0 = 3 × 10 -3 . Fig. 23 presents the distribution function and mesh obtained with the method AMW-SLDG. The method AMW-CDG gives extremely similar results.

Here we can see initial mesh slightly rened on the two beams. Along time it follows and matches vortices during their formation, existence, and during their merging.

Fig. 24 shows the number of cells and nest level of computation during simulation. During creation and existence of vortices the mesh is locally rened to its maximum and the number of cells increases.

The relative variation of mass, of L 1 and L 2 -norms and of total energy are presented on Fig. 25. The four quantities are very well conserved compared to some results presented in [START_REF] Heath | A discontinuous Galerkin method for the Vlasov-Poisson system[END_REF][START_REF] Besse | Semi-Lagrangian schemes for the Vlasov equation of an unstructured mesh of phase space[END_REF]. Relative variations of mass and of L 1 -norm show a dierence between the two schemes, especially during the merging of vortices.

Polar-like test cases

This part is devoted to experiment considering the following formulation of Vlasov equation from [START_REF] Frénod | Long Time Simulation of a Beam in a Periodic Focusing Channel via a Two-Scale PIC-Method[END_REF] and [START_REF] Frénod | Long-time simulation of a highly oscillatory Vlasov equation with an exponential integrator[END_REF]

∂ t f (r, v, t) + ∂ r v ε f (r, v, t) + ∂ v (E ε (r, t) + Ξ ε (r, t)) f (r, v, t) = 0, (69)
where the constant ε and elds E and Ξ will be dened for each test case. Under this formulation, the corresponding Poisson equation is

1 r ∂ r (rE ε (r, t)) = ρ(r, t). (70)
The boundary condition for equation ( 70) is E ε (0, t) = 0. This is necessary for a solution to exist due to factor 1 r in the equation.

Focusing beam

For this case ε = 0.01, E ε is obtained by solving (70), and Ξ ε (r, t) = -r ε. The initial distribution is

f 0 (r, v) = 1 √ 2πv th exp - v 2 2v 2 th 1 [-0.75,0.75] (r), (r, v) ∈ [-1, 1] × [-1.2, 1.2], (71)
with v th = 0.0727518214392 and

1 [a,b] (r) = ⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ 0 if r < a or r > b, 1 else.
The nal time is T = 100 with 10 3 time steps for a time step ∆t = 0.1. We use polynomials of degree up to 2. Maximum level of renement is 8 and threshold 0 = 3 × 10 -3 . Fig. 26 illustrates the distribution function and the mesh obtained with the method AMW-CDG. The method AMW-SLDG gives results that can not be distinguished. In this test case, we can see the formation of a central vortex with two coiling laments. Fig. 27 presents the relative variation of numerical mass, L 1 and L 2 -norms. Here again the variation of mass and of L 1 -norm indicate that the distribution function has some negative values.

36

Focusing channel

For this test case ε = 0.1, E ε is computed using (70) and

Ξ ε (r, t) = r -1 ε + cos 2 t ε .
The initial distribution function is

f 0 (r, v) = 3 4v th exp -v 2 2v 2 th 1 [-1.83271471003,1.83271471003] (r), (r, v) ∈ [-3, 3] 2 ,
(72) with v th = 0.0727518214392. This test case is detailed in [START_REF] Frénod | Long Time Simulation of a Beam in a Periodic Focusing Channel via a Two-Scale PIC-Method[END_REF].

The nal time is T = 100 with 5 × 10 3 time steps for a time step ∆t = 0.02. We use polynomials of degree up to 2. The maximum level of renement is 8 and the threshold is 0 = 10 -2 . Fig. 28 shows snapshots of the distribution function and the mesh obtained with the AMW-SLDG scheme. We obtain very similar result using the AMW-CDG scheme. The mesh accurately follows the lamentation of the distribution function.

Astrophysics test cases

In astrophysics, (1) becomes

∂ t f (x, v, t) + ∂ x (vf (x, v, t)) + ∂ v (E(x, t)f (x, v, t)) = 0, x ∈ Ω x , v ∈ R, f (x, v, 0) = f 0 (x, v), (73)
∆Φ(x, t) = 4π(ρ(x, t) -1), x ∈ Ω x (74a) E(x, t) = -∂ x Φ(x, t), (74b)
Here, the eld E is gravity eld. Note that we still add the constant -1 in Poisson equation for periodicity. The rst simulation of Vlasov-Poisson for gravitational case using discontinuous Galerkin method was done in [START_REF] Cheng | Numerical study of Vlasov-Poisson equations for innite homogeneous stellar systems[END_REF].

Cold layer

We consider initial condition

f 0 (x, v) = 1 0.15 √ 2π exp - (v -u(x)) 2 2 × 0.15 2 , (x, v) ∈ [0, 2π] × [-10, 10], (75)
with u(x) = 0.01 sin(x).

The nal time is 40 with 2 × 10 3 time steps for a time step ∆t = 0.02. We use polynomials of degree 2. The maximum level of renement is 8 and the threshold is 0 = 3 × 10 -3 . This test case is designed to always have cells on the nest level of renement. the AMW-CDG scheme. Finally, our mesh becomes coarse and dissipation of laments at time t = 14 and t = 20 appears. Here again the distribution functions from the two schemes (AMW-SLDG and AMW-CDG) can not be distinguished.

Periodicity is an artefact that we have introduced in our code and we can see its side eect at t = 6 on boundaries x = 0 and x = 2π, and later on the simulation. The next test case 7.8.2 will be adapted to avoid this numerical artefact. We see laments that go under x = 0 or above x = 2π and appear on the other side. This is a non-physical eect but a spurious result of periodicity. The lamentation and dissipation processes can be seen on Fig. 30 too. Filamentation occurs between times t = 2 and t = 5 and the number of cells increases a lot during this period. After time t = 12 dissipation becomes signicant, details are lost and number of cells starts to decrease. The solution is not physical any-more.

Gaussian initial condition

We consider the Gaussian initial condition from Sec. 2.1.1 from [START_REF] Colombi | Vlasov-Poisson in 1D: waterbags[END_REF]:

f 0 (x, v) = 4 exp - (x 2 + v 2 ) 0.08 , (x, v) ∈ [-2, 2] 2 (76)
with the proper normalisation to respect periodicity in Poisson equation. This test case is designed with open boundary conditions so we adapted our domain with their results to ensure an almost null distribution function everywhere on boundaries.

(a) Initial distribution function. (i) Distribution function at t=20. (g) Distribution function at t = 30. The nal time is T = 100 with 10 3 time steps for a time step ∆t = 0.1. We use polynomials of degree up to 3. The maximum level of renement is 9 and the threshold is 0 = 3 × 10 -3 . Fig. 31 displays the evolution of the distribution function and the mesh obtained with the method AMW-SLDG. Results from the method AMW-CDG are extremely similar. The simulation is performed on [-2, 2] 2 to enforce Dirichlet boundary conditions, but because nothing happen outside of [-1, 1], we only present a zoom on [-1, 1]. Because it is nearly impossible to distinguish laments on long time, Fig. 32 is a contour plot to highlight the very thin laments at nal time. We clearly see the very thin laments coiling without diusion, even thought the laments are extremely close from each others. In this test case the domain is large enough so that no numerical artefacts due to periodicity occur. Fig. 33 shows the increase of number of cell as thin laments appear and coil.

The main diagnostics are presented on Fig. 34. Fig 34b and34e show the correlation between the increase of the L 1 -norm and the fall of the minimum. The decrease of the L 2 -norm occurs at the same time. Step variation of mass is typical of amplied machine error variation.

Jeans instability

This test case is based on astrophysics formulation of the Vlasov equation (73) but uses the same normalisation of the Poisson equation as plasma physics [START_REF] Alpert | A Class of Bases in L 2 for the Sparse Representation of Integral Operator[END_REF]. The initial condition is

f (x, v, 0) = exp -v 2 2 √ 2π (1 -A cos(kx)) , (x, v) ∈ [0, 2π k] × [-V c , V c ]. (77)
We use the settings from [START_REF] Cheng | Numerical study of Vlasov-Poisson equations for innite homogeneous stellar systems[END_REF]. We consider A = 0.01, k = 0.8 and V c = 6. The nal time is T f = 100 with 1, 000 time steps, which gives ∆t = 0.1. We use polynomials of degree 2 in each direction. The maximum level of renement is 8 and the threshold parameter is 0 = 0.003. The distribution function and the mesh obtained with the AMW-CDG scheme are presented on Fig. 35. The results from the AMW-SLDG scheme are ex- tremely similar. Variations of mass, L 1 and L 2 -norms and momentum are displayed on Fig. 36. The two schemes (AMW-CDG and AMW-SLDG) give very similar results. Fig. 37 depicts the distribution function as a function of the particular energy E = v 2 2 -φ(x). This gure is very similar to Fig. 7a from [START_REF] Cheng | Numerical study of Vlasov-Poisson equations for innite homogeneous stellar systems[END_REF]. As in the two stream instability, Sec. 7.5, for small and moderate value of the energy E, the distribution function f is not a single-valued function of E. Then a stationary state of BGK type is not reached yet. We still observe that numerical diusion is very small until t = 100 (see Fig. 36c of L 2 -norm evolution, the loss of relative L 2 -norm is less than 0.5%). In order to see the eect of numerical diusion at T f = 100, we perform some simulations with the threshold parameter 0 = 0.01 and 7 levels of renement. The corresponding plots of the distribution function as a function of the particular energy are presented on Fig. 38. The corresponding distribution function at time T f = 100 is displayed on Fig. 39. We observe on Fig. 38 that the distribution function f is a single-valued function of the energy E. Then a stationary state of BGK type seems to be reached. Moreover for such test-case we observe a greater numerical diusion than in the rst parameter setting (i.e. 0 = 0.003, 8 levels of renement); let us see Fig. 40 where we observe a loss of relative L 2 -norm of 0.7%. Therefore, as for the two stream instability (see Sec. 7.5), we may conclude the same hypothesis about the role of numerical diusion on the behaviour of the solution.

Extension to the four-dimension problem

We extend the schemes to the four dimension problem. We still use directional splitting. Then the extension to four dimensions is just an addition of one dimension steps, as presented Secs. 3.2 and 3.3. Adaptivity is achieved by extended equations (40) to the four dimension problem.

We consider the equation The nal time of simulation is T = 2π with 200 time steps for a time step ∆t = π 100. We use polynomials of degree at most 3. The maximum level of renement is 6 and the threshold is 0 = 10 -10 . The maximum level of renement is always reached in this test case.

∂ t f + v x ∂ x f -v y ∂ y f -x∂ vx f + y∂ vy f = 0, (78)
The projections presented on Figs. 41 and 42 are obtained using the formulae P f (x, y) = with P f (x, y) the projection on plane (x, y). The initial projection on plane (v x , v y ) is null because the two Gaussian compensate each other. The nal time T = 2π is supposed to give the exact same distribution function than at initial time. However, it is clear with the projection on plane (v x , v y ) that we do not have an exact numerical rotation, but still very close to the exact solution. Fig. 42 is obtained with the AMW-CDG scheme. The AMW-SLDG scheme gives extremely similar results. The projection on plane (v x , v y ) enables to see the integral of the dierence between the numerical solution and the exact solution as the integral of the exact solution on plane (v x , v y ) is null. Fig. 43 shows the evolution of mass, L 1 , L 2 and L ∞ -norm of the two methods. The variation of total mass remains on range of slight amplication of machine error. Variation of L 1 and L 2 -norms are small, less than 2% at most for the L 1 -norm, and present some periodic peak. These are linked to the mesh and come from numerical artefacts. The total number of cells, displayed on Fig. 44, is very stable with a small diminution at the beginning. The total number of cells remains around 500, 000 cells. If one used a uniform mesh scaled on the nest mesh, this would requires 16, 777, 216 cells. This means that we use only 3% of the computation that would be required for a uniform-mesh simulation. Although the mass and norm variations are slightly better with the AMW-SLDG scheme, we see that it also uses a very few more cells. Such four-dimensional simulations show the potential of the AMW-SLDG and AMW-CDG schemes to be used and work well in four dimensions.

Concluding remarks

We have implemented and tested an adaptive version of two numerical schemes, namely the semi-Lagrangian discontinuous Galerkin scheme and the characteristics-discontinuous Galerkin scheme. Both schemes include a multiresolution discretisation using a multi-wavelet representation. Most of the time they give identical or extremely similar results. This observation is consistent with the results of Ref. [START_REF] Qiu | Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov Poisson system[END_REF] and conrm that our adaptive procedure does not introduce a bias or break the symmetry between these two schemes. It is an interesting and open problem to understand whether the AMW-SLDG and AMW-CDG schemes are equivalent or not. In the armative case, it is even ambitious but crucial to show it. Both methods have proven to react very well to multiscale decomposition and adaptive mesh. Our future works include a twodimensional relativistic Vlasov-Maxwell solver to compare our results with [START_REF] Besse | A wavelet-MRA-based adaptive semi-Lagrangian method for the relativistic Vlasov-Maxwell system[END_REF] or a possible extension to a four-dimensional Vlasov-Poisson solver. A challenging problem is to perform a rigorous mathematical analysis of the AMW-SLDG and AMW-CDG schemes to show their convergence and to obtain a priori error estimates. Finally an ecient parallel version of our algorithm will be a matter of future work.

Figure 1 :

 1 Figure 1: Mesh representation and notation.

Figure 4 :

 4 Figure 4: Transport of a 2D cell in one direction.

Fig. 4 , 1 and T 2

 412 domains T 1 and T 2 are bounded with black lines and domain T * 1 in red is obtained by transporting T 1 according to Sec. 3.3. Domains T 1 ∩ T * ∩ T * 1 are equivalent for the following analysis. Let us use only T 2 ∩ T * 1 . We start from

  (a) Quadrature points on interval [a, b]. Red stars are the quadrature points. These quadrature points give exact integration distribution function tested with polynomials of degree up to 3. (b) Quadrature points on interval [a, b]. Red stars are the quadrature points. Here polynomials are independently dened on [a, (a + b) 2] and on [(a + b) 2, b]. The quadrature points do not give exact integration.

Figure 5 :

 5 Figure 5: Quadrature points on interval [a, b] when [a, b] is single segment and when it is two segments.

Figure 6 :

 6 Figure 6: Points for ux computation. Blue stars correspond to the points where ux are computed on the big cell on the left. Red stars correspond to the points where ux are computed on the small cell on the right.

  (a) Initial mesh M n . (b) Predicted mesh. (c) Predicted mesh after renement. (d) Merged mesh M n+1 .

Figure 7 :

 7 Figure 7: Drawing of mesh prediction.

Figure 8 :

 8 Figure 8: Drawing of the transport prediction of the grey cell.

5 :

 5 third step of Strang splitting. The third and last step of Strang splitting is identical to Step 2 (rst step of Strang splitting) except that it uses data from Step 4 and its result is the distribution function at time t n+1 . The mesh is altered only at very specic steps. It is rened only during the renement stage of steps 2, 4 and 5, and coarsened only during the coarsening stage of these steps 2, 4 and 5. The mesh is xed for the computation and the evolution of the distribution function, and for solving the Poisson problem.

  (a) Initial distribution function.

Figure 9 :

 9 Figure 9: Initial distribution function and mesh for swirling deformation.

  Mesh at t=0.75. (c) Swirling at t=1.10.

  Mesh at t=1.10. (e) Swirling at t=1.50.

Figure 10 :

 10 Figure 10: Swirling deformation at various time for AMW-SLDG.

Figure 11 :

 11 Figure 11: Relative variations of numerical mass, L 1 and L 2 -norms for swirling deformation.

Figure 12 :

 12 Figure 12: Electric energy for weak Landau damping.

Fig. 16 displays

 16 the evolution of the number of cells and of the maximum level of renement. If the maximum level is reached very fast, the number of (a) Initial distribution function. (b) Distribution function at t=4. (c) Distribution function at t=8. (d) Distribution function at t=20. (e) Distribution function at t=50. (f) Distribution function at t=80.

Figure 13 :

 13 Figure 13: Distribution function for strong Landau damping for the AMW-CDG scheme.

Figure 14 :Figure 15 :

 1415 Figure 14: Mesh for strong Landau damping for the AMW-CDG scheme.

  Maximum level of renement.

Figure 16 :

 16 Figure 16: Number of cells and level of renement for strong Landau damping.

  Distribution function at t = 8 pp.

  Mesh at t = 8 pp.

  Distribution function at t = 20 pp.

  Mesh at t = 20 pp.

  Distribution function at t = 60 pp.

  Mesh at t = 60 pp.

Figure 17 :

 17 Figure 17: Distribution function and Mesh for two stream instability for the AMW-SLDG scheme.

  Mesh at t = 100 pp.

Figure 17 :

 17 Figure 17: Distribution function and Mesh for two stream instability for the AMW-SLDG scheme.

Figure 18 :

 18 Figure 18: Evolution of numerical mass, L 1 and L 2 -norms, minimum of distribution function, total energy and momentum for two stream instability.

Figure 19 :

 19 Figure 19: Distribution function as a function of particular energy at time T f = 100.

Figure 20 :

 20 Figure 20: Distribution function for two stream instability with the AMW-SLDG scheme at time t = 210.

Figure 21 :Figure 22 :

 2122 Figure 21: Distribution function as a function of particular energy at time T = 210.

  Distribution function at t = 60 pp.

  Mesh at t = 60 pp.

  Distribution function at t = 80 pp. Mesh at t = 80 pp.

  Distribution function at t = 200 pp.

  Mesh at t = 200 pp.

Figure 23 :

 23 Figure 23: Distribution function and Mesh for bump on tail with the AMW-SLDG scheme.

Figure 24 :

 24 Figure 24: Number of cells and maximum level for computation for bump on tail.

  Relative variation of total energy.

Figure 25 :

 25 Figure 25: Variation of mass, L 1 and L 2 -norms and total energy for Bump on tail.

Fig. 29 depicts

 29 the evolution of the distribution function and the mesh with the initial condition and the coiling of laments at t = 2, t = 4 and t = 6 with (a) Initial distribution function.

  Distribution function at t = 30.

  Mesh at t = 30. (e) Distribution function at t = 40. Mesh at t = 40.(g) Distribution function at t = 100. Mesh at t = 100.

Figure 26 :

 26 Figure 26: Distribution function and mesh for focusing beam with the AMW-CDG scheme.

Figure 27 :

 27 Figure 27: Relative variation of numerical mass, L 1 and L 2 -norms for focusing beam.

  Distribution function at t = 10.

  Mesh at t = 10. (e) Distribution function at t = 16.

  Mesh at t = 16.(g) Distribution function at t = 20.

  Mesh at t = 20.

Figure 28 :

 28 Figure 28: Distribution function and mesh for Focusing channel with AMW-SLDG.

  Distribution function at t=2.

  Mesh at t=4. (g) Distribution function at t=6.

Figure 29 :

 29 Figure 29: Distribution function and mesh for initial conditions (75) with the AMW-CDG scheme.

  Mesh at t=20.

Figure 29 :Figure 30 :

 2930 Figure 29: Distribution function and mesh for initial conditions (75) and during simulation with the AMW-CDG scheme.

  (a) Initial distribution function.

  Distribution function at t = 5.

  Mesh at t = 5.(e) Distribution function at t = 15. Mesh at t = 15.

  Mesh at t = 30.

Figure 31 :

 31 Figure 31: Distribution function and mesh for Gaussian initial distribution function and evolution in astrophysics framework with the AMW-SLDG scheme.

Figure 31 :

 31 Figure 31: Distribution function and mesh for Gaussian initial distribution function and evolution in astrophysics framework with the AMW-SLDG scheme.

  Distribution function on[-1, 1].

Figure 32 :

 32 Figure 32: Distribution function contours at t = 100 for Gaussian initial distribution function with the AMW-SLDG scheme.

Figure 33 :

 33 Figure 33: Number of cells and maximum level for Gaussian initial condition.

  Relative variation of L 1 -norm.

Figure 34 :

 34 Figure 34: Variation of the main diagnostics for Gaussian distribution function.

  Distribution function at t = 10.

  Mesh at t = 20.

  Distribution function at t = 20. Mesh at t = 20.

  Distribution function at t = 100.

  Mesh at t = 10.

Figure 35 :

 35 Figure 35: Distribution function and mesh for the Jeans instability with AMW-CDG scheme.

Figure 36 :

 36 Figure 36: Variation of the mains diagnostics for the Jeans instability.

Figure 37 :-

 37 Figure 37: Distribution function as a function of particular energy for the Jeans instability with few diusion at time T f = 100, with 0 = 0.003 and 8 levels of renement.

Figure 38 :

 38 Figure 38: Distribution function as a function of particular energy for the Jeans instability with higher diusion at time T f = 100 with 0 = 0.01 and 7 levels of renement.

  Distribution function at time T f = 100. Contour plot at time T f = 100.

Figure 39 :Figure 40 :

 3940 Figure 39: Final distribution function for Jeans instability with the AMW-CDG scheme, 0 = 0.01 and 7 levels of renement

with initial condition f 0 -v 2 y-

 02 (x, y, v x , v y ) = exp -(x + 1) 2 -(y -3) 2 -v 2 x exp -(x -1) 2 -(y + 3) 2 -v 2 x -v 2 y , (x, y, v x , v y ) ∈ [-10, 10] 4 . (79)Equation (78) depicts the combination of a rotation in plan (x, v x ) with angular momentum +1 with a rotation in plan (y, v y ) with angular momentum -1.

10 f

 10 (x, y, v x , v y ) dv x dv y Projection on plane(x, y).

  Projection on plane(x, vx).

  Projection on plane(y, vy).

  Projection on plane(vx, vy).

Figure 41 :

 41 Figure 41: Initial projections.

  Projection on plane(x, y).

  Projection on plane(x, vx).

  Projection on plane(y, vy).

  Projection on plane(vx, vy).

Figure 42 :

 42 Figure 42: Projections at nal time t = 2π with AMW-CDG.

Figure 43 :

 43 Figure 43: Variation of the main diagnostics for the four-dimensional test case.

Figure 44 :

 44 Figure 44: Number of cells in the four-dimensional test case.

Table 1 :

 1 Dierence of momentum between levels.

	levels	dierence
	7 with 6 9.05942 E -14
	6 with 5 7.10543 E -15
	7 with 5 9.76996 E -14

Acknowledgements

We are grateful to Stéphane Colombi and Thierry Sousbie for many fruitful discussions concerning numerical methods, astrophysics test cases and OpenMP handling. This work was supported by the VLASIX and EUROFUSION projects, respectively under grants No ANR-13-MONU-0003-01 and EURATOM-WP15-ENR-01/IPP-01. We are also grateful to referees whose constructive comments lead to some improvements of the paper.

(d) Error obtained with the AMW-CDG scheme for polynomials of degree 3. (d) Error obtained with the AMW-CDG scheme for polynomials of degree 3.

Table 3: Error for the Burgers' equation with AMW-SLDG and AMW-CDG schemes for polynomials of degree 2 and 3.