Joackim Bernier

Fernando Casas

Nicolas Crouseilles

SPLITTING METHODS FOR ROTATIONS: APPLICATION TO VLASOV EQUATIONS by

Keywords: 2000 Mathematics Subject Classification. - .splitting, rotation, Vlasov equations, high-order time integrators

In this work, a splitting strategy is introduced to approximate two-dimensional rotation motions. Unlike standard approaches based on directional splitting which usually lead to a wrong angular velocity and then to large error, the splitting studied here turns out to be exact in time. Combined with spectral methods, the so-obtained numerical method is able to capture the solution to the associated partial differential equation with a very high accuracy. A complete numerical analysis of this method is given in this work. Then, the method is used to design highly accurate time integrators for Vlasov type equations: the Vlasov-Maxwell system and the Vlasov-HMF model. Finally, several numerical illustrations and comparisons with methods from the literature are discussed.

Introduction

The main goal of this work is to introduce a splitting strategy to deal with rotations motions and to apply it to construct efficient high order time integrators for Vlasov type equations. The splitting is based on the fact that a rotation of angle θ can be decomposed into a product of three shear transformations (1.1)

1 -tan θ/2 0 1

1 0 sin θ 1 1 -tan θ/2 0 1 = cos θ -sin θ sin θ cos θ = e θJ ,
for θ = kπ, k ∈ Z and where J is the fundamental symplectic matrix (1.2) J = 0 -1 1 0 .

Note that this decomposition into shear matrices can be derived using formal computations and has been already introduced in the image processing community (see [START_REF] Paeth | A fast Algorithm for General Raster Rotation[END_REF][START_REF] Tanaka | A rotation method for raster image using skew transformation[END_REF][START_REF] Andres | The quasi-shear rotation[END_REF][START_REF] Welling | Rotation of 3D volumes by Fourierinterpolated shears Graphical Models[END_REF][START_REF] Chen | 3D volume rotation using shear transformation[END_REF]), in which several approaches have been developed to rotate an image on a computer screen. Moreover, this approach has also been used to design numerical methods for Gross-Pitaevskii equations (see [START_REF] Chin | Fourth-order algorithms for solving the imaginarytime Gross-Pitaevskii equation in a rotating anisotropic trap[END_REF] and [START_REF] Bader | Fourier methods for the perturbed harmonic oscillator in linear and nonlinear Schrödinger equations[END_REF], Lemma II.2) in which this underlying splitting is used to solve exactly the harmonic oscillator.

To make the link between (1.1) and the underlying partial differential equation, we introduce the following two-dimensional transport equation

(1.3) ∂ t u = Jx • ∇u, x ∈ R 2 ,
with the initial condition u(t = 0, x) = u in (x). The exact solution of (1.3) at time t writes u(t, x) = u in (e tJ x) which is nothing but the rotation of angle t of the initial condition u in . When the initial condition is not known analytically or when equation (1.3) is a part of a more complicated model, then one only has access to a discrete information of the initial condition and a numerical method is required to approximate (1.3). Our goal in this work is to introduce a directional splitting inspired by (1.1) which is exact with respect to the time variable.

Obviously, standard finite differences or finite volumes based methods can be used to approximate the spatial direction x and coupled to Runge-Kutta strategies in time. However, this leads to methods which usually suffer from strong CFL condition on the time step. Then, semi-Lagrangian methods are preferred, since they are free from stability condition still keeping Eulerian accuracy (see [START_REF] Sonnendrücker | The semi-Lagrangian method for the numerical resolution of the Vlasov equation[END_REF][START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF][START_REF] Zerroukat | The Parabolic Spline Method (PSM) for conservative transport problems[END_REF][START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF]). For (1.3), the feet of the characteristics can be computed exactly and a two-dimensional interpolation has to be performed to update the numerical unknown. However, high-dimensional interpolation is known to be non conservative and it is obviously more demanding in terms of complexity and time. Then, splitting methods are very competitive since they reduce the problem into very simple one-dimensional linear transport equations which can be solved efficiently with semi-Lagrangian methods (using high-order or even spectral interpolation). Moreover, in a splitting procedure, the variable that does not appear in the derivative is just a parameter so that a very simple parallelization can be performed by distributing the computation on the processors according to the values of this parameter.

For rotation dynamics however, the standard splitting strategy (like Strang or Lie splitting for example) can induce some error since it involves a wrong rotational velocity (see [START_REF] Boris | Relativistic plasma simulation-optimization of a hybrid code[END_REF]). Here, we propose a new splitting which enables to solve (1.3) exactly in time (like in [START_REF] Chin | Fourth-order algorithms for solving the imaginarytime Gross-Pitaevskii equation in a rotating anisotropic trap[END_REF][START_REF] Bader | Fourier methods for the perturbed harmonic oscillator in linear and nonlinear Schrödinger equations[END_REF]). Moreover, when this splitting is coupled with spectral methods (and under some assumptions detailed in the sequel), the so-obtained method is able to capture to a very high accuracy the exact solution (spectral accuracy in practice). A complete proof of convergence of the fully discretized numerical method is performed. We will see that this strategy and some simple extensions turn out to be very efficient compared to standard methods when applied to the following problems. First, it enables us to design high order (in time) methods for the Vlasov-Maxwell system. Second, when applied to the Vlasov-HMF model in the close-to-equilibrium regime, as in [START_REF] Horsin | Comportement en temps long d'équations de type Vlasov : études mathématiques et numériques[END_REF], this splitting turns out to be more accurate than the Strang one, at the same cost.

Concerning the Vlasov-Maxwell solvers, our goal was to improve the method introduced in [START_REF] Crouseilles | Hamiltonian splitting for the Vlasov-Maxwell equations[END_REF] in which a splitting into three parts has been proposed. Among these three parts, two were solved exactly in time whereas for the magnetic part, a standard directional splitting was performed. Here, the new method enables us to also solve this part exactly in time. This is then very helpful when designing high order splitting methods for the full Vlasov-Maxwell system. The resulting schemes are fourth order accurate in time and preserve the Gauss condition exactly. We also use the new splitting to approximate the solution of the Vlasov-HMF system, for which the close to equilibrium dynamics is driven by the linearized Hamiltonian part (see [START_REF] Horsin | Comportement en temps long d'équations de type Vlasov : études mathématiques et numériques[END_REF]). For such Hamiltonian, the new splitting has a good behavior (see [5]) and we compare its efficiency with the standard Strang splitting by studying perturbations of a non homogeneous equilibrium state.

The rest of the paper is organized as follows. First, the method is presented in the context of the numerical approximation of transport equation of the form (1.3) and a complete proof of convergence is performed with some numerical illustrations. Then, the Vlasov-Maxwell system is presented and we explain how the new method is used to design high-order Vlasov-Maxwell solvers. Finally, some numerical results are given to show the benefit of the new method in the Vlasov context.

Presentation of the method and its numerical analysis

In this section, we focus on the following two-dimensional equation (2.1)

∂ t u = Jx • ∇u, x = (x 1 , x 2 ) ∈ R 2 ,
supplemented with an initial condition u(t = 0, x) = u in (x). We intend to analyse the convergence of a splitting in time based numerical scheme coupled with a spectral method in space (i.e., in the x-direction). More precisely, we want to solve (2.1) on [t n , t n+1 ]; then we want to compute u n+1 (x), an approximation of u(t n+1 , x 1 , x 2 ), the solution at time t n+1 = t n + δ t (δ t > 0 being the time step and n ∈ N) of (2.1) with initial condition u in (x 1 , x 2 ) = u(t n , x 1 , x 2 ) at time t n = nδ t , n ∈ N. To do so, we propose a new splitting in which each step is a shear transformation.

Let us introduce some notations. For a given 2 × 2 matrix A, we denote by exp(δ t Ax • ∇)u n the solution at time t n+1 of (2.2)

∂ t u(t, x) = Ax • ∇u(t, x), x ∈ R 2 u in (x) = u n (x). ,
Then, inspired by (1.1), we search for a, b ∈ R so that the following relation holds true

(2.3) e -a 2 x 2 ∂x 1 e bx 1 ∂x 2 e -a 2 x 2 ∂x 1 u n = e δtJx•∇ u n ,
which can be written equivalently as

(2.4) e A 1 x•∇ e A 2 x•∇ e A 1 x•∇ u n = e δtJx•∇ u n , with (2.5
)

A 1 = 0 -a/2 0 0 , A 2 = 0 0 b 0 .
Using the method of characteristics, we have for (2.2)

e δtAx•∇ u n = u n • e δtA , δ t ≥ 0, so that (2.4) is nothing but u n (e A 1 e A 2 e A 1 x) = u n (e δtJ x).
Since A 1 and A 2 are nilpotent matrices, their exponential follows readily

e A 1 = 1 -a/2 0 1 , e A 2 = 1 0 b 1 ,
and it is clear from (1.1) that the choice a = 2 tan(δ t /2) and b = sin(δ t ) leads to an exact splitting in time, so that the scheme then writes as u n+1 (x) = u n (e A 1 e A 2 e A 1 x), with A 1 and A 2 given by (2.5). Let us remark that the usual Strang splitting corresponds to a = b = δ t .

In consequence, now we have to solve shear transformations which is nothing but one-dimensional linear advections. We consider here using a pseudospectral method. To do so, we discretize a square, of size R, centered in 0, (i.e.

-R 2 , R 2
2 ) with a regular grid with N ∈ N * points per direction. Its stepsize is h = R/N . We denote this grid by G 2 , with

(2.6) G = h - N -1 2 , N 2 .
Then, we define the discrete partial Fourier transforms

F 1 :    C G 2 → C G×G u → h g 1 ∈G u g 1 ,g 2 e -ig 1 ξ 1 and F 2 :    C G 2 → C G× G u → h g 2 ∈G u g 1 ,g 2 e -ig 2 ξ 2 , where G = η -N -1 2 , N 2
stands for the set of discrete frequencies with η = 2π/R. Now, we want to solve the continuous shear transformations (α ∈ R):

∂ t u = αx 2 ∂ x 1 u, ∂ t u = αx 1 ∂ x 2 u, (2.7)
which are the basic building blocks of the splitting presented above. These shear transformations are particularly simple to solve and we shall use a pseudo-spectral method. Then, for any parameter α ∈ R, we introduce two pseudo-spectral shear transformations, (2.8)

S α 1 :

C G 2 → C G 2 u → F -1 1 e iαξ 1 g 2 F 1 u and (2.9) S α 2 : C G 2 → C G 2 u → F -1 2 e iαξ 2 g 1 F 2 u Remark 2.1.
-If N is even, we have to pay attention to the mode N 2 associated to the frequency ηN 2 . Indeed, we can easily verify that

S α i R G 2 ⊂ R G 2 (for i = 1, 2) if and only if N is odd or α ∈ Z.
Finally, the numerical solution (u n ) n∈N of the numerical schemes we consider are defined by (for

δ t = kπ, k ∈ Z ) (2.10) u n = (L δt ) n u in |G 2 := (S δt 2 S -δt 1 ) n u in |G 2 , (Lie)
u n = (T δt ) n u in |G 2 := (S -δt/2 1 S δt 2 S -δt/2 1 ) n u in |G 2 , (Strang)
u n = (M δt ) n u in |G 2 := (S -tan(δt/2) 1 S sin(δt) 2 S -tan(δt/2) 1 ) n u in |G 2 , (New)
where u in |G 2 is the evaluation of the initial condition u in on the grid G 2 . The main goal of this section is now to perform a complete numerical analysis of these splittings defined in (2.10).

2.1. Numerical analysis. -We define some associated discrete Lebesgue norms. They are defined for u ∈ C G 2 by

u 2 L 2 (G 2 ) = h 2 g∈G 2 | u g | 2 and u L ∞ (G 2 ) = max g∈G 2 | u g |.
We also use the Schwartz space S (R 2 ) and introduce a scale of spaces, denoted (X s ) s≥0 , defined by

X s = u ∈ L 2 (R 2 ), u 2 X s := x 2s |u(x)| 2 dx + ξ 2s |F u(ξ)| 2 dξ < ∞ ,
where x = 1 + |x| 2 and F u denotes the Fourier transform of u.

2.1.1. Consistency. -First, we prove that the pseudo-spectral shear transformations (2.8) and (2.9) are consistent with the continuous ones (2.7). Let us remark that in addition to the analysis of the spectral consistency, we will also pay attention to the truncation R. The consistency error is of the pseudospectral shear transformations is stated in the following proposition for S α 1 but the result is also valid for S α 2 .

Proposition 2.2. -For all s > 1 and for all M > 0, there exists c > 0 such that for all u ∈ S (R 2 ), α ∈ (-M, M ), R > 0 and N ∈ N * we have

S α 1 u -v L 2 (G 2 ) ≤ c |α| R -s + h s √ h u X s+6 .
where

u = u |G 2 and v = v |G 2 with v(x) = u(x 1 + αx 2 , x 2 ).
Proof. -Applying the discrete Fourier-Plancherel isometry, we get

(2.11) h 2 g∈G 2 | (S α 1 u) g -v g | 2 = hη 2π (ξ 1 ,g 2 )∈ G×G | (F 1 S α 1 u) ξ 1 ,g 2 -(F 1 v) ξ 1 ,g 2 | 2 .
Thus, we are going to expand F 1 v and F 1 u with respect to u. More precisely, we apply the Poisson formula to get

F 1 u = h g 1 ∈hZ u(g 1 , g 2 )e -iξ 1 g 1 -h g 1 ∈G c u(g 1 , g 2 )e -iξ 1 g 1 = F 1 u(ξ 1 , g 2 ) + k∈Z * F 1 u(ξ 1 + 2kπ h , g 2 ) -h g 1 ∈G c u(g 1 , g 2 )e -iξ 1 g 1 ,
where

F 1 u(ξ 1 , x 2 ) = u(x)e -iξ 1 x 1 dx 1 is the continuous Fourier transform of u along the first direction and G c = hZ \ G. Consequently, since F 1 v(ξ 1 , x 2 ) =
e iαξ 1 x 2 F 1 u, we decompose the consistency error into three terms

(F 1 S α 1 u) ξ 1 ,g 2 -(F 1 v) ξ 1 ,g 2 = k∈Z * 1 -e iα 2kπ h g 2 e iαξ 1 g 2 F 1 u(ξ 1 + 2kπ h , g 2 ) (ε ξ ,g 2 ) + h g 1 ∈G c 1 -e iαξ 1 g 2 u(g 1 , g 2 )e -iξ 1 g 1 (ε ξ ,g 2 ) + h g 1 ∈G c [u(g 1 + αg 2 , g 2 ) -u(g 1 , g 2 )] e -iξ 1 g 1 . (ε ξ ,g 2 )
Now we bound each one of these three consistency errors.

Estimation of ε 1 : First, we have

| ε 1 ξ 1 ,g 2 | ≤ k∈Z * α 2kπ h g 2 F 1 u(ξ 1 + 2kπ h , g 2 ) = k∈Z * α 2kπ h g 2 ξ 1 + 2kπ h -s-1 F 1 (|∂ x 1 | s+1 u)(ξ 1 + 2kπ h , g 2 ) ≤ k∈Z * α 2kπ h g 2 (2|k| -1)π h -s-1 F 1 (|∂ x 1 | s+1 u)(ξ 1 + 2kπ h , g 2 ) ≤ 4|α|ζ(s + 1) 1 + g 2 2 π h s sup x 2 ∈R R (1 + x 2 2 )|(|∂ x 1 | s u)(x 1 , x 2 )|dx 1 ,
where ζ denotes the Riemann function. This estimate involves a norm of u that is neither usual nor isotropic. Furthermore, the estimates of ε 2 and ε 3 will lead to some other norms of this kind. Consequently, in order to get an estimate as readable as possible, we control these norms by the X s+6 norm. Such a control can be realized with classical techniques of pseudo-differential calculus. As these estimates are technical but not crucial here, we omit details (the interested reader could refer for example to [START_REF] Nicola | Global pseudo-differential calculus on Euclidean spaces, Pseudo-Differential Operators[END_REF]). Now, we observe that by monotonicity we have (2.12) h

g 2 ∈G\{0} 1 1 + g 2 2 ≤ R 1 1 + y 2 dy ≤ π.
Thus, since ε 1 ξ 1 ,0 = 0, there exists a constant c > 0, depending only on s, such that hη 2π

(ξ 1 ,g 2 )∈ G×G | ε 1 ξ 1 ,g 2 | 2 ≤ c|α| 2 R -1 h 2s (# G) u X s+6 ≤ c|α| 2 h 2s-1 u 2 X s+6 .
Estimation of ε 2 : First, naturally, we control ε 2 by

(2.13) | ε 2 ξ 1 ,g 2 | ≤ α|ξ 1 ||g 2 | h g 1 ∈G c u(g 1 , g 2 )e -iξ 1 g 1 .
In order to absorb the factor ξ 1 on the left, we perform a discrete integration by parts. So, we assume that ξ 1 = 0, we denote

ξ 1 = k 1 η and g 1 = g 1 = n 1 h, where k 1 ∈ -N -1 2 , N 2
and

n 1 ∈ Z \ -N -1 2 , N 2
. Then we introduce N + = 1 + N/2 and N -= -1 -(N -1)/2 . Consequently, we have

h g 1 ∈G c u(g 1 , g 2 )e -iξ 1 g 1 =h n 1 ≥N + u(g 1 , g 2 )e -2iπn 1 k 1 N + h n 1 ≤N - u(g 1 , g 2 )e -2iπn 1 k 1 N =h n 1 ≥N + u(g 1 , g 2 ) -u(g 1 + h, g 2 ) h h n 1 n=N + e -2iπnk 1 N (E + ) + h n 1 ≤N - u(g 1 , g 2 ) -u(g 1 -h, g 2 ) h h N - n=n 1 e -2iπnk 1 N . (E -)
To control (E + ), first we observe that since 0 ≤ |k 1 | ≤ N/2, we have

h n 1 n=N + e -2iπnk 1 N ≤ 2h |1 -e -2iπk 1 N | ≤ c hN 2π|k 1 | = c |ξ 1 | ,
where c is a universal constant. Then, by application of the mean value theorem, we get

|E + | ≤ c |ξ 1 |(1 + g 2 2 )
sup

x∈R 2 (1 + x 2 2 )|x 1 | s+1 |∂ x 1 u(x)| h n 1 ≥N + g s+1 1 ≤ c s |ξ 1 |(1 + g 2 2 ) u X s+6 R -s 1 N n 1 ≥0 N + + n 1 N s+1 ≤ c s |ξ 1 |(1 + g 2 2 ) u X s+6 R -s 1 N n 1 ≥0 1 2 + n 1 N s+1 (2.14)
where c s > 0 is a constant depending only on s. We recognize a Riemann sum, so we have

1 N n 1 ≥0 1 2 + n 1 N s+1 ----→ N →∞ 2 -s-2 ∞ 1 y -s-1 dy = 2 -s-2 s .
In particular, since this sequence converges, it is bounded by a constant depending only on s. Thus, we obtain the following bound for |E + |:

|E + | ≤ c s |ξ 1 |(1 + g 2 2 ) u X s+6 R -s ,
where c s is another constant depending only on s. Note that, by symmetry, the same control holds for E -. Finally, coming back to (2.13) and using (2.12), we have another constant, denoted c s , depending only on s such that hη 2π

(ξ 1 ,g 2 )∈ G×G | ε 2 ξ 1 ,g 2 | 2 ≤ c s |α| 2 R -1 R -2s (# G) u 2 X s+6 ≤ c s |α| 2 h -1 R -2s u 2 X s+6 .
(2.15)

Estimation of ε 3 : Let us introduce a useful technical lemma whose proof is postponed to the Appendix.

Lemma 2.3. -If y 1 , y 2 , y 3 , λ ∈ R are such that y 3 ∈ [y 1 ; y 1 + λy 2 ] then we have y 3 y 2 ≥ |y 1 | √ 1 + λ 2 .
Then applying the mean value theorem, for any g 1 , g 2 , α ∈ R, there exists m g 1 ,g 2 ,α in [g 1 ; g 1 + αg 2 ] such that

u(g 1 + αg 2 , g 2 ) -u(g 1 , g 2 ) = αg 2 ∂ x 1 u(m g 1 ,g 2 ,α , g 2 ).
Since |α| ≤ M , applying Lemma 2.3, we get

| ε 3 ξ 1 ,g 2 | ≤ h g 1 ∈G c |α||g 2 ||∂ x 1 u(m g 1 ,g 2 ,α , g 2 )| ≤ h g 1 ∈G c |α| 1 + g 2 2 m g 1 ,g 2 ,α g 2 -s-1 sup x∈R 2 (1 + x 2 2 )|x| s+1 |∂ x 1 u(x)| ≤ c s |α| 1 + g 2 2 R 2 √ 1 + M 2 -s-1 u X s+6   h g 1 ∈G c 2g 1 R -s-1   ,
where c s is a constant depending only on s.

Then, carrying out the same procedure as in (2.14), we get another constant c s depending only on s such that   h

g 1 ∈G c 2g 1 R -s-1   ≤ c s R.
Thus we have the estimate

| ε 3 ξ 1 ,g 2 | ≤ c s,M |α| 1 + g 2 2 R -s u X s+6 ,
where c s,M is a constant depending only on s and M . Consequently, we can repeat the same argument as (2.15) to get hη 2π

(ξ 1 ,g 2 )∈ G×G | ε 3 ξ 1 ,g 2 | 2 ≤ c s,M |α| 2 h -1 R -2s u 2 X s+6 ,
where c s,M is another constant depending only on s and M . We conclude by summing the different contributions of ε 1 , ε 2 and ε 3 .

2.1.2. Backward error analysis.

-We aim at describing the long time behavior of the splitting methods. So, we perform a general backward error analysis (1) for a large class of methods including Lie and Strang splittings but also the new splitting. Note that since we deal with a linear problem the expansions are convergent.

Proposition 2.4. -If a, b ∈ R satisfy ab < 2 then (2.16) e bx 1 ∂x 2 e -ax 2 ∂x 1 = e JL a,b x•∇ ,

and

(2.17) e -a 2 x 2 ∂x 1 e bx 1 ∂x 2 e -a 2 x 2 ∂x 1 = e JS a,b x•∇ , where

F (x) =      arcsin( √ x) √ x if 0 < x ≤ 1 asinh( √ -x) √ -x if x < 0 1 if x = 0.
Proof. -Considering the transport equation (2.2) which can be solved with the method of characteristics, we have

e tAx•∇ u 0 = u in • e tA .
Thus (2.16) is equivalent to exp 0 -a 0 0 exp 0 0 b 0 = e JL a,b , (1) The reader can refer to [START_REF] Hairer | Geometric numerical integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] for an overview on backward error analysis.

with J given by (1.2). These exponentials of matrices can be written as shear transforms. So (2.16) is equivalent to

(2.19) P a,b := 1 -a 0 1 1 0 b 1 = 1 -ab -a b 1 = e JL a,b .
Similarly, (2.17) is equivalent to

(2.20) 1 -a/2 0 1 1 0 b 1 1 -a/2 0 1 = e JS a,b .
First, we prove that if (2.19) holds with L a,b given by (2.18) then (2.19) also holds with S a,b given by (2.18). Indeed, observing that a Lie splitting is always conjugate to the Strang splitting we have

1 -a/2 0 1 1 0 b 1 1 -a/2 0 1 = 1 a/2 0 1 P a,b 1 -a/2 0 1 = 1 a/2 0 1 e JL a,b 1 -a/2 0 1 . But 1 -a/2 0 1 is symplectic, i.e. t 1 -a/2 0 1 J 1 -a/2 0 1 = J.
Thus, we have

1 a/2 0 1 e JL a,b 1 -a/2 0 1 = exp J t 1 -a/2 0 1 L a,b 1 -a/2 0 1 = e JS a,b ,
where S a,b is given by (2.18). So, now we aim at proving (2.19). The existence of such a L a,b is ensured by the following lemma (an elementary proof is given in the Appendix). But, by construction, e JL a,b = P a,b , so L a,b is an eigenvector associated with the eigenvalue 1 of the linear application

R a,b : S 2 (R) → S 2 (R) Q → t P a,b QP a,b .
By a straightforward calculation, we observe that

(2.21) Q a,b = b ab 2 ab 2 a satisfies R a,b (Q a,b ) = Q a,b .
Then we deduce of the following lemma (proven in Appendix) that if 0 < ab < 4 then there exists µ a,b ∈ R such that

(2.22) L a,b = µ a,b Q a,b . Lemma 2.6. -If 0 < ab < 4 the eigenspace of R a,b associated with the eigenvalue 1 is of dimension 1. Now, we just have to determine µ a,b . Since L a,b is symmetric, it is diago- nalizable in an orthonormal basis, i.e. ∃λ ∈ R 2 , ∃Ω ∈ O 2 (R), L a,b = Ω -1 λ 1 λ 2 Ω = Ω -1 DΩ.
So, since J and Ω commute, we have

P a,b = Ω -1 e JD Ω.
Since we assume that 0 < ab < 4, we deduce from (2.22) that L a,b is either positive or negative. In particular we have λ 1 λ 2 > 0. Thus we can define the symplectic matrix

K = 4 λ 1 /λ 2 4 λ 2 /λ 1 .
This matrix satisfies √ λ 1 λ 2 t KK = D and J t K = K -1 J . Thus, we have (2.23)

P a,b = (KΩ) -1 e √ λ 1 λ 2 J (KΩ) = (KΩ) -1 cos( √ λ 1 λ 2 ) -sin( √ λ 1 λ 2 ) sin( √ λ 1 λ 2 ) cos( √ λ 1 λ 2 ) (KΩ).
In particular, we have

Tr P a,b = 2 cos( λ 1 λ 2 ) = 2 cos( det L a,b ) = 2 cos(µ a,b det Q a,b ).
As a consequence, since det L a,b goes to zero when (a, b) goes to 0, we deduce of a straightforward calculation that if ab is small enough then

µ a,b = ±F (ab(1 -ab/4)).
Finally, we have to determine the sign of µ a,b . First, observe that by continuity, we have either µ a,b > 0 for all a, b small enough satisfying ab > 0 or µ a,b < 0 for all a, b small enough satisfying ab > 0. This second case is excluded because when a goes to zero we have

e -F (a 2 (1-a 2 /4))JQa,a = e -aJ+O(a 2 ) = P -a,-a + O(a 2 ) = P a,a + O(a 2 ).
To conclude, we have proved that if ab > 0 and (a, b) is small enough then

P a,b = e F (ab(1-ab/4))JQ a,b .
Furthermore, this relation is analytic with respect to a and b, so it can be extended to all a, b ∈ R such that ab < 2. Indeed, under this assumption we have ab(1 -ab/4) ∈ (-∞, 1) which is the domain of analyticity of F .

The classical splitting formulas of Lie and Strang correspond to the choice a = b = δ t in (2.16) and (2.17). However, as mentioned in the introduction, these choices are not necessarily the best. For the Strang like splittings, a straigthforward calculation proves that there exists an optimal choice for which the splitting is exact. This choice can be obtained by direct formal calculations by assuming a decomposition of the rotation matrix.

Lemma 2.7. -If δ t ∈ (-π, π) then we have S 2 tan(δt/2),sin(δt) = δ t 1 0 0 1 .
Note that due to the non-diagonal terms of L a,b , it is impossible to design an exact splitting based on the Lie splitting.

Convergence.

-We now consider the convergence of the pseudospectral splittings (2.10) to approximate our problem (2.1). Then, for a discrete initial condition u = u in |G 2 , the numerical solution at time t n = nδ t (n ∈ N) is given by n compositions of the operators defined in (2.10). For instance, for the standard Strang splitting, the numerical solution at time t n is (T δt ) n u. In the following theorem, we show that, up to a spectral spatial error, the dynamics generated by the Strang pseudo-spectral method T δt and by the Lie pseudo-spectral method L δt over very long times.

Theorem 2.8. -For all s > 0 there exists c > 0 such that for all N ∈ N * , all R > 0, all u ∈ S (R 2 ), all n ∈ N and all δ t ∈ [-1, 1] with t n = nδ t , we have

(L δt ) n u -e tnJS L δ t x•∇ u |G 2 L 2 (G 2 ) ≤ c t n R -s + h s √ h u X s+6 , and
(T δt ) n u -e tnJS T δ t x•∇ u |G 2 L 2 (G 2 ) ≤ c t n R -s + h s √ h u X s+6 , where u = u |G 2 , S L δt := L δt,δt /δ t = I 2 + O(δ t ) and S T δt := S δt,δt /δ t = I 2 + O(δ 2 t
). The definitions of S a,b and L a,b are given by (2.18) in Proposition 2.4, whereas L δt and T δt are given by (2.10).

Proof. -We focus only on proving the convergence estimate for the Lie splitting. The same proof could be applied to prove the estimate for the Strang splitting.

Let ε n ∈ L 2 (G 2 ) be the consistency error at time t n . It is defined by

ε n = L δt e tnJS L δ t x•∇ u |G 2 -e t n+1 JS L δ t x•∇ u |G 2 .
As usual, for linear schemes, the convergence error is given by

(L δt ) n u -e tnJS L δ t x•∇ u |G 2 = n-1 k=0 L n-1-k δt ε k .
Here, L δt is an isometry of L 2 (G 2 ), so we have

(L δt ) n u -e tnJS L δ t x•∇ u |G 2 L 2 (G 2 ) ≤ n-1 k=0 ε k L 2 (G 2 ) ≤ n sup k∈N ε k L 2 (G 2 ) .
Thus, we just have to bound ε k . Using formulas of Proposition 2.4, we decompose ε k into two consistency errors for the pseudo-spectral shear transformations

ε k = S δt 2 e -δtx 2 ∂x 1 e tnJS L δ t x•∇ u |G 2 -e δtx 1 ∂x 2 e -δtx 2 ∂x 1 e tnJS L δ t x•∇ u |G 2 + S δt 2 S -δt 1 e tnJS L δ t x•∇ u |G 2 -e -δtx 2 ∂x 1 e tnJS L δ t x•∇ u |G 2 .
Then applying Proposition 2.2, we get a constant c > 0, depending only on s > 0 such that

ε k L 2 (G 2 ) ≤ c |δ t | R -s + h s √ h s e tnJS L δ t x•∇ u X s+6 + e -δtx 2 ∂x 1 e tnJS L δ t x•∇ u X s+6 .
Now, we introduce a lemma to control these norms, whose proof is available in the Appendix.

Lemma 2.9. -For all κ > 0 and all s > 0 there exists a constant c > 0

such that if τ ∈ GL 2 (R) satisfies (2.24) ∀x ∈ R 2 , κ -1 |x| ≤ |τ (x)| ≤ κ|x| then for all u ∈ S (R 2 ) we have u • τ X s ≤ c u X s . We recall that if A ∈ M 2 (R) then e (Ax•∇) u = u in • e A .
∀δ t ∈ [-1, 1], ∀x ∈ R 2 , κ -1 |x| ≤ 1 δ t 0 1 x ≤ κ|x|.
For the other estimate, we observe that the associated quadratic form S L δt is a constant of the motion of exp(tJS L δt ): for all t ∈ R and all

δ t ∈ [-1, 1] we have (2.25) ∀x ∈ R 2 , t e tJS L δ t x S L δt e tJS L δ t x = t xS L δt x.
Furthermore, for all

δ t ∈ [-1, 1], we have det S L δt = δ -2 t arcsin 2 δ 2 t (1 -δ 2 t /4) > 0.
So, S L δt is either a positive or negative. Thus, since (S L δt ) 1,1 > 0, it is positive. Then, since δ t → S L δt is a continuous map, S L δt and (S L δt ) -1 are bounded uniformly with respect to δ t ∈ [-1, 1]. Consequently, there exists κ > 0 such that for all δ t ∈ [-1, 1] and all x ∈ R 2 we have

κ -1 t xS L δt x ≤ κ -1 |S L δt ||x| 2 ≤ |x| 2 ≤ κ|(S L δt ) -1 | -1 |x| 2 ≤ κ t xS L δt x.
Thus we deduce of (2.25) that for all t ∈ R, all δ t ∈ [-1, 1] and all x ∈ R 2 we have |e

tJS L δ t x| 2 ≤ κ t e tJS L δ t x S L δt e tJS L δ t x = κ t xS L δt x ≤ κ 2 |x| 2 . and |e tJS L δ t x| 2 ≥ κ -1 t e tJS L δ t x S L δt e tJS L δ t x = κ -1 t xS L δt x ≥ κ -2 |x| 2 .
As a corollary, we deduce the convergence error of these methods.

Corollary 2.10. -For all s > 0 and all h 0 > 0, there exists c > 0 such that for all N ∈ N * , all R > 0, all u ∈ S (R 2 ), all n ∈ N and all δ t ∈ [-1, 1] and h = R/N ≤ h 0 , denoting t n = nδ t , we have

(L δt ) n u -e tnJx•∇ u |G 2 L 2 (G 2 ) ≤ c t n R -s + h s √ h u X s+6 +c |e tnJ -e tnJS L δ t | u X 4 , and
(T δt ) n u -e tnJx•∇ u |G 2 L 2 (G 2 ) ≤ c t n R -s + h s √ h u X s+6 +c |e tnJ -e tnJS T δ t | u X 4 ,
where u = u |G 2 , L δt and T δt are given by (2.10)

Proof of Corollary 2.10. -We only focus on proving the convergence estimate for the Lie splitting, the case of the Strang splitting being similar. Applying Theorem 2.8 and the triangle inequality, we have

(L δt ) n u -e tnJx•∇ u |G 2 L 2 (G 2 ) ≤ c t n R -s + h s √ h u X s+6 + E u,δt,n,G , with E u,δt,n,G = e tnJS L δ t x•∇ u -e tnJx•∇ u |G 2 L 2 (G 2 )
. Consequently, to prove the corollary, we just have to bound E u,δt,n,G by

|e tnJ -e tnJS L δ t | u X 4 .
First, we introduce a technical lemma that will be proved in the Appendix.

Lemma 2.11. -There exists an universal constant c > 0 such that for all v ∈ H 2 (R 2 ), all R > 0 and all N ∈ N * we have

v |G 2 L 2 (G 2 ) ≤ u L 2 (R 2 ) + c h 2 ∆u L 2 (R 2 ) .
Since h ≤ h 0 , applying this lemma we get a constant c > 0, depending only on h 0 > 0, such that

E u,δt,n,G ≤ c (1 -∆) e tnJS L δ t x.∇ u -e tnJx.∇ u L 2 (R 2 ) .
Then applying the Fourier Plancherel isometry, we get

E u,δt,n,G ≤ c 2π (1 + |ξ| 2 ) F u • t e -tnJS L δ t -F u • t e -tnJ L 2 (R 2 ) .
Then introducing a Taylor remainder under its integral form, it comes

E u,δt,n,G ≤ (1 + |ξ| 2 ) 1 0 ∇ ξ F u(y α,ξ,n,δt ) • t e -tnJS L δ t -e -tnJ ξ dα L 2 (R 2 ) ≤ |e -tnJS L δ t -e -tnJ | max α∈(0,1) (1 + |ξ| 2 ) 3/2 ∇ ξ F u(y α,ξ,n,δt ) L 2 (R 2 ) (2.26)
where

y α,ξ,n,δt = t M α,n,δt ξ and M α,n,δt = I 2 -α e -tnJS L δ t -e -tnJ . Now, we distinguish two cases. If |e -tnJS L δ t -e -tnJ | ≤ 1/2 then we deduce that we have |M α,n,δt -I 2 | ≤ 1 2
. Consequently, we have

| det M α,n,δt | ≥ κ and |M -1 α,n,δt | ≤ 2,
where κ is an universal constant.

Thus, by performing a natural change of coordinates, we get

E u,δt,n,G ≤ |e -tnJS L δ t -e -tnJ | κ (1 + | t M α,n,δt -1 ξ| 2 ) 3/2 ∇ ξ F u L 2 (R 2 ) ≤ 8 |e -tnJS L δ t -e -tnJ | κ (1 + |ξ| 2 ) 3/2 ∇ ξ F u L 2 (R 2 ) ≤ c|e -tnJS L δ t -e -tnJ | u X 4
where c > 0 is an universal constant. Finally, we have to consider the case where |e

-tnJS L δ t -e -tnJ | ≥ 1/2.
Applying Lemma 2.11 and the Fourier-Plancherel isometry we get two constant c, κ > 0 depending only on h 0 such that

E u,δt,n,G ≤ u L 2 (G 2 ) + e tnJx•∇ u |G 2 L 2 (G 2 ) ≤ c (1 -∆)u L 2 (R) ≤ κ|e -tnJS L δ t -e -tnJ | u X 4 .
Next, we focus on the new splitting M δt . We provide a theorem showing that its dynamics corresponds, up to a spectral spatial error, to the rotation with the exact speed, for very long times.

Theorem 2.12. -For all s, ν > 0 there exists c > 0 such that for all N ∈ N * , all R > 0, all u ∈ S (R 2 ), all n ∈ N and all δ t ∈ R satisfying |δ t | < π -ν, denoting t n = nδ t , we have

(2.27) (M δt ) n u -e tnJx•∇ u |G 2 L 2 (G 2 ) ≤ c t n R -s + h s √ h u X s+6 ,
where u = u |G 2 , and M δt is given by (2.10).

Proof.

-By carrying out the same proof as in Theorem 2.8, we could easily prove that for all s, ν > 0 there exists c > 0 such that for all N ∈ N * , all R > 0, all u ∈ S (R 2 ), all n ∈ N and all δ t ∈ R satisfying |δ t | < π -ν, denoting t n = nδ t , we have

(M δt ) n u -e tnJS M δ t x•∇ u |G 2 L 2 (G 2 ) ≤ c t n R -s + h s √ h u X s+6
where u = u |G 2 and S M δt := S 2 tan(δt/2),sin(δt) /δ t where S a,b is given by (2.18). Thus, to conclude this proof, we just have to observe that by Lemma 2.7 we have S 2 tan(δt/2),sin(δt

) = δ t I 2 . Remark 2.13. -For all u ∈ L 2 (G 2 ) we have u L ∞ (G 2 ) ≤ h -1 u L 2 (G 2 )
, thus (2.27) gives a control of convergence error with the discrete L ∞ norm for very long times:

(M δt ) n u -e tnJx•∇ u |G 2 L ∞ (G 2 ) ≤ c t n R -s + h s h 3/2 u X s+6 .
2.2. Numerical illustrations. -In this subsection, we intend to illustrate the different results obtained previously, namely the spatial accuracy of pseudo-spectral method and the time accuracy of the time splitting.

Spatial accuracy

First, we present some numerical results to illustrate the estimates obtained in Proposition 2.2. To do so, we consider the following function

u(x) = exp - |x| 2 2 , x = (x 1 , x 2 ) ∈ R 2 ,
which is shifted by α = 0.01. We denote v |G 2 where v(x) = u(x 1 + αx 2 , x 2 ) the exact shifted solution, and we compute the (discrete) L 2 norm of the difference between S α 1 u |G 2 and v |G 2 . The spatial grid G 2 is defined by (2.6) where h = R/N , R = 15 and different values of N are considered to check the spatial accuracy. The results are displayed in Figures 1 and2. One can observe that for large h (or small N ), the term R -s is negligible and the term h s gives the exponential decreasing of the error which is the typical behavior of spectral methods. On the contrary, for very small values of h (or large values of N ), the term R -s /h -1/2 becomes prominent even if the error is quite small (around 10 -11 ).

Time accuracy

In this part, we give some numerical illustrations of the efficiency of the new splitting. To do so, we consider the following equation (2.28)

∂ t u = Jx • ∇ x u, x = (x 1 , x 2 ) ∈ R 2 ,
with the initial condition

u in (x) = 1 2πβ exp - (x 1 -0.3) 2 β + exp - (x 1 + 0.3) 2 β exp - x 2 2 β ,
with β = 0.01. The spatial truncated domain [-2, 2] 2 is discretized with the grid G 2 defined by (2.6) with R = 4 and a space step h = R/N , N = 243 = 3 5 . The time step is h ≈ 0.139 and the final time is T = 10 5 (the number of iterations is 71888). In Figures 3456, some results are displayed where we compare the exact solution, the solution given by (T δt ) n u in |G 2 (Strang splitting and spectral interpolation), the solution given (L δt ) n u in |G 2 (Lie splitting and spectral interpolation) and the solution given by the new method (M δt ) n u in (see (2.10)). First, in Figure 3, the three solution are plotted at the final time. We can observe that the exact solution and the solution obtained by the new method are very close whereas the solution obtained by the Strang splitting is not good due to the fact that the angular velocity of the Strang method is not exact. To precise these observations, we plot on Figure 4 (Figure 5 is a zoom) the relative L ∞ error of the new method, the Strang and the Lie methods. The error produced by the new method is close to 10 -13 which is the spectral error. On the contrary, the Strang and Lie methods periodically produce an error of order one. This is due to its wrong angular velocity: the solution moves away from the exact solution producing large error and at some times, the Strang method recover the exact solution so that the error becomes very small. These times can be computed from the above analysis. Indeed, from the rotation speed of the Strang splitting

ω δt = arcsin( √ δ 2 t (1-δ 2 t /4)) δt
, we deduce that the exact solution (which rotates with a speed ω ex = 1) and the numerical solution obtained by the Strang method will coincide every times T such that t n + ω ex T = t n + ω δt T [π] (the factor π (instead of a factor 2π) is due to our choice of a symmetric initial condition). Then, we have T = π/(ω δt -1) which gives with our choice of time step δ t ≈ 0.139, T ≈ 3888. We can observe a very good agreement on Figures 4 and5 and also on Figure 6 for which δ t = π/4 and then T ≈ 113. Finally, we study the performance of the new method. Indeed, we compare the new splitting and a direct two-dimensional solving of (2.28). The direct resolution is done by a semi-Lagrangian type strategy: at each time step, we first compute exactly the feet of the characteristics equations and we then use a two-dimensional spectral interpolation by means of the non uniform fast Fourier transform (the so-called nufft procedure introduced in [START_REF] Greengard | Accelerating the Nonuniform Fast Fourier Transform[END_REF]). We checked that this approach also leads to spectral accuracy, and we want here to compare the two spectral methods in terms of CPU time with respect to the total number of points N 2 (N being the number of points per direction). The results are displayed in Figure 7: the time execution (for 10 iterations) for both methods (new splitting and nufft) as a function of N 2 (for N = 2 5 , . . . , 2 11 ), in log -log scale. Even if both methods have the same complexity O(N 2 log(N )), the new approach clearly has a smaller constant (around ten times smaller). Moreover, in such a splitting procedure, a simple and efficient parallelization can be performed since the variable that does not appear in the derivative is just a parameter.

Application to the Vlasov-Maxwell equations

In this section, we intend to apply the above splitting to the context of the 1+1/2 Vlasov-Maxwell system. Indeed, the time discretization of this system is based on a time splitting, and one of the pieces (the so-called magnetic part) corresponds to a rotation in the velocity direction due to the presence of the self-consistent electromagnetic field. Then, instead of using a Strang splitting like in [START_REF] Crouseilles | Hamiltonian splitting for the Vlasov-Maxwell equations[END_REF], we shall use the exact splitting presented in the previous section, so that this magnetic part will be solved exactly in time and with a spectral accuracy in the velocity directions. This is very helpful to design high order methods for the full Vlasov-Maxwell system. After introducing the 1+1/2 Vlasov-Maxwell system we intend to solve, the splitting method introduced in [START_REF] Crouseilles | Hamiltonian splitting for the Vlasov-Maxwell equations[END_REF] is recalled and then high order methods dedicated to systems split into three parts are introduced. (E 1 , E 2 )(t, x 1 ) which are determined by solving the following system of evolution equations (3.1)

∂ t f + v 1 ∂ x 1 f + E • ∇ v f -BJv • ∇ v f = 0, ∂ t B = -∂ x 1 E 2 , ∂ t E 2 = -∂ x 1 B - R 2 v 2 f (t, x 1 , v)dv + J 2 (t), ∂ t E 1 = - R 2 v 1 f (t, x 1 , v)dv + J 1 (t), where v = (v 1 , v 2 ), J i (t) = 1/|L| L R 2 v i f (t, x 1 , v)dx 1 dv, i = 1, 2 (|L| denotes
the measure of L) and J denotes the symplectic matrix (1.2). This reduced system, which has been considered in several former studies (see [START_REF] Califano | Kinetic saturation of the Weibel instability in a collisionless plasma[END_REF][START_REF] Cheng | Discontinuous Galerkin methods for VlasovMaxwell equations[END_REF][START_REF] Crouseilles | Hamiltonian splitting for the Vlasov-Maxwell equations[END_REF]), has to be supplemented with the Gauss condition

(3.2) ∂ x 1 E 1 (t, x 1 ) = R 2 f (t, x 1 , v)dv -1, ∀t ≥ 0,
and with initial conditions

f (t = 0, x 1 , v) = f in (x 1 , v), E 2 (t = 0, x 1 ) = E in 2 (x 1 ) and B(t = 0, x 1 ) = B in (x 1 ). Notice that E in 1 (x 1
) is implied by the Gauss condition (3.2) at the initial time.

Splitting method.

-Here we propose to use the splitting method introduced in [START_REF] Crouseilles | Hamiltonian splitting for the Vlasov-Maxwell equations[END_REF], by reformulating the Vlasov-Maxwell system into

dF dt = H E (F ) + H f (F ) + H B (F ), F (0) = F in ,
where the fields H E (F ), H f (F ) and H B (F ) will be written below. We denote by F (δ t ) = (f, E 1 , E 2 , B)(δ t ) the solution of the Vlasov-Maxwell system (3.1). This solution can be formally written as

F (δ t ) = ϕ δt (F in ) := exp((H E + H f + H B )δ t )F in , where F in = (f in , E in 1 , E in 2 , B in
) denotes the initial condition. Now, we want to use a splitting method to approximate the system (3.1). To do so, we shall use the splitting introduced in [START_REF] Crouseilles | Hamiltonian splitting for the Vlasov-Maxwell equations[END_REF][START_REF] Crouseilles | An asymptotic preserving scheme for the relativistic Vlasov-Maxwell equations in the classical limit[END_REF] based on a decomposition into three parts corresponding respectively to the fields H E (F ), H f (F ) and H B (F ). Then, a first order Lie method based on this decomposition writes

χ δt (F in ) = ϕ δt (F in ) + O(δ t 2 ) with (3.3) χ δt = ϕ [H E ] δt • ϕ [H f ] δt • ϕ [H B ] δt
where ϕ

[H E ] δt , ϕ [H f ] δt , ϕ [H B ] δt
denotes the exact solutions corresponding to the fields H E , H f and H B . Using these notations, the adjoint [START_REF] Hairer | Geometric numerical integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF] of the Lie method χ t writes (3.4)

χ δt = ϕ [H B ] δt • ϕ [H f ] δt • ϕ [H E ] δt .
In the following we write down the equations associated to the fields H E , H f and H B :

ϕ [H E ] δt : ∂ t f + E • ∇ v f = 0, ∂ t E = 0, ∂ t B = -∂ x 1 E 2 , ϕ [H f ] δt : ∂ t f + v 1 ∂ x 1 f = 0, ∂ t E = - R 2 vf dv + J , ∂ t B = 0, ϕ [H B ] δt : ∂ t f -BJv • ∇ v f = 0, ∂ t E 1 = 0, ∂ t E 2 = -∂ x 1 B, ∂ t B = 0.
Then, as mentioned in [START_REF] Crouseilles | Hamiltonian splitting for the Vlasov-Maxwell equations[END_REF], ϕ

[H E ] δt
and ϕ

[H f ] δt can be computed exactly in time and efficiently in phase space using spectral methods. However, the computation of ϕ [H B ] δt was performed using a Strang splitting. Instead, we shall use the new splitting M δt introduced above to compute exactly in time ϕ

[H B ] δt
and efficiently in phase space using spectral methods. Let us remark that the application of the new splitting to the H B part requires a slight modification. Indeed, to solve ∂ t f -BJv • ∇ v f = 0 (with B constant in time during this part) on one time step δ t from an initial condition f in (defined on the velocity grid), we will use the new splitting with a modified time step Bδ t to capture the right rotation speed, i.e. (M Bδt ) with (M δt ) defined by (2.10).

Based on the fact that each step can be computed exactly in time, we now look for efficient integration methods for systems separable into three parts which enable us to design efficient high order methods in time. A simple and efficient way to achieve this goal is to consider compositions of a first-order method with its adjoint computed at fractional step sizes. This is the main subject of the next part.

Composition methods for systems separable into three parts.

-To simplify the presentation, we restrict ourselves to ordinary differential equations (ODEs). The so-obtained composition methods will then be used within the Vlasov-Maxwell framework.

Let us consider the following ODE

(3.5) dx dt (t) = u(x(t)), x(0) = x in ∈ R D ,
with D ∈ N , whose exact solution at time t = δ t will be denoted by x(δ t ) = ϕ δt (x in ). We are interested in problems where u in (3.5) can be split into three parts,

u(x) = u a (x) + u b (x) + u c (x)
in such a way that the exact flows ϕ

[a] δt , ϕ [b] δt , ϕ [c]
δt , corresponding to u a , u b , u c , respectively, can be computed exactly. One might consider then splitting methods of the form

(3.6) ϕ [a] asδt • ϕ [b] bsδt • ϕ [c] csδt • • • • • ϕ [a] a 1 δt • ϕ [b] b 1 δt • ϕ [c] c 1 δt
and fix the coefficients a i , b i , c i , i = 1, . . . , s so that it provides an approximation of order, say, p. It turns out, however, that the number of order conditions to be satisfied by these parameters grows very rapidly with the order. Thus, time-symmetric schemes of order p = 4 (resp. p = 6) require solving 11 (resp. 56) conditions. A more convenient way consists in considering compositions of χ δt and its adjoint χ δt , with (3.7)

χ δt = ϕ [a] δt • ϕ [b] δt • ϕ [c] δt and χ δt = ϕ [c] δt • ϕ [b] δt • ϕ [a]
δt . More specifically, we construct integrators within the family (3.8)

G 1 ≡ ψ δt = χ α 1 δt • χ α 2 δt • • • • • χ α 2s-1 δt • χ α 2s δt : s ≥ 1, (α j ) 1≤j≤2s ∈ R 2s ,
where χ δt and χ δt are given by (3.7), so that (3.9)

χ δt (x in ) = ϕ δt (x in ) + O(δ 2 t )
, and an analogous relation for χ δt . Composition integrators ψ δt ∈ G 1 are timesymmetric (self-adjoint) whenever they have left-right palindromic sequences of coefficients α i , i.e. if α 2s+1-i = α i , i = 1, . . . , s [START_REF] Hairer | Geometric numerical integration: Structure-Preserving Algorithms for Ordinary Differential Equations[END_REF].

Notice that one could achieve methods of order p within this family even if only first-order approximations to the flows ϕ

G 2 = {ψ δt = φ α 1 δt • • • • • φ αsδt : s ≥ 1, (α 1 , . . . , α s ) ∈ R s } ,
where φ δt : R D -→ R D is any second-order self-adjoint integrator. Notice that, if φ δt is chosen as

φ δt = χ δt/2 • χ δt/2 , then G 2 is contained in G 1 .
These integrators also enjoy the time-symmetric property if α s+1-i = α i , i = 1, . . . , s.

3.3.1. Analysis of the methods. -For the analysis, it is convenient to introduce the graded Lie algebra associated with the vector field defining the ODE (3.5) and its corresponding exact flow ϕ δt . As is well known, for each infinitely differentiable map g : R D -→ R, g(ϕ δt (x)) admits an expansion of the form

g(ϕ δt (x)) = e δtF [g](x) = g(x) + k≥1 δ t k k! F k [g](x), x ∈ R D ,
where F is the vector field associated with u,

(3.11) F = D i=1 u i (x) ∂ ∂x i .
Similarly, for the basic first-order method χ δt defined by (3.9), one has

g(χ δt (x)) = e Y δ t [g](x) with Y δt = k≥1 δ t k
Y k and for its adjoint χ δt also defined in (3.9), one has g(χ δt (x)) = e -Y -δ t [g](x). Then, one can formally compute the operator series associated to any integrator ψ δt ∈ G 1 defined by (3.8)

Ψ δt = exp(Y δtα 1 ) exp(-Y -δtα 2 ) • • • exp(Y δtα 2s-1 ) exp(-Y -δtα 2s ).
By repeated application of the Baker-Campbell-Hausdorff formula we can express formally Ψ δt as the exponential of an operator F δt ,

Ψ δt = e F δ t , with F δt = k≥1 δ t k F k , δ t k F k ∈ L k for each k ≥ 1 and L = k≥1 L k is
Y δtα i = δ t α i Y 1 + (δ t α i ) 2 Y 2 + (δ t α i ) 3 Y 3 + • • • -Y -δtα i = δ t α i Y 1 -(δ t α i ) 2 Y 2 + (δ t α i ) 3 Y 3 -• • • so that Ψ δt = exp δ t w 1 Y 1 + δ t 2 w 2 Y 2 + δ t 3 (w 3 Y 3 + w 12 [Y 1 , Y 2 ]) +δ t 4 (w 4 Y 4 + w 13 [Y 1 , Y 3 ] + w 112 [Y 1 , [Y 1 , Y 2 ]]) + O(δ t 5 ) ,
where w 1 , w 2 , . . . are polynomials in the coefficients α i . In particular [START_REF] Blanes | Splitting and composition methods in the numerical integration of differential equations[END_REF],

w 1 = 2s i=1 α i , w 2 = 2s i=1 (-1) i α 2 i , w 3 = 2s i=1 α 3 i , w 4 = 2s i=1 (-1) i α 4 i , w 12 = 1 2   2s-1 i=1 (-1) i+1 α 2 i 2s j=i+1 α j + 2s-1 i=1 α i 2s j=i+1 (-1) j α 2 j  
Then, a time-symmetric 4th-order method has to satisfy only consistency (w 1 = 1) and the order conditions at order three, that is w 3 = w 12 = 0. Let us remark that conditions at even order (w 2 = w 4 = 0) are automatically verified by symmetry. Notice, then, that the minimum number of maps to be considered in ψ δt ∈ G 1 is s = 3.

3.3.2. Methods of order 4. -It turns out, however, that methods involving the minimum number of maps (or stages) do not usually provide the best efficiency. In other words, considering additional stages (and thus some free parameters) leads to more efficient schemes, even when the computational cost per step is also higher. The difficulty then lies in the way the free parameters are fixed. In this respect, several objective functions have been considered in the literature. In particular we mention the following [START_REF] Blanes | Composition methods for differential equations with processing[END_REF] (let us recall that α = (α 1 , . . . , α 2s ) ∈ R 2s )

(3.12) E 1 (α) = 2s i=1 |α i | and E 2 (α) = 2s 2s i=1 α 5 i 1/4 .
The quantity E 2 is usually the dominant error term for a number of problems. The criterion we follow here will be to look for symmetric methods with small values of E 1 which, in addition, have also small values of E 2 . In the sequel, we consider composition methods in the class G 1 with s = 3, 4, 5, 6 (see (3.8)) which have been designed by optimizing both functions E 1 and E 2 .

Case s = 3. The integrator reads

(3.13) ψ [3] δt = χ α 1 δt • χ α 2 δt • χ α 3 δt • χ α 3 δt • χ α 2 δt • χ α 1 δt
and the unique (real) solution to the order conditions w 1 = 1, w 3 = w 12 = 0 is given by

α 1 = α 2 = 1 2(2 -2 1/3 ) , α 3 = 1 2 -2α 1 . If χ δt = ϕ [a] δt • ϕ [b] δt • ϕ [c]
δt , then it involves 13 maps (the minimum number).

For future reference, the values of the objective functions are

E 1 (α) = 4.40483, E 2 (α) = 4.55004. Remark 3.2.
-Notice that this corresponds to the familiar scheme of Yoshida [START_REF] Yoshida | Construction of higher order symplectic integrators[END_REF]

ψ δt = φ γδt/2 • φ βδt • φ γδt/2 in G 2 with γ = 1/(2 -2 1/3
). Moreover, this method is also recovered in [START_REF] Koseleff | Exhaustive search of symplectic integrators using computer algebra. Integration algorithms and classical mechanics[END_REF] when considering splitting methods of the form (3.6).

Case s = 4. The composition is

(3.14) ψ [4] δt = χ α 1 δt • χ α 2 δt • χ α 3 δt • χ α 4 δt • χ α 4 δt • χ α 3 δt • χ α 2 δt • χ α 1 δt
, involving 17 maps. Now we have a free parameter, which we take as α 1 . The minima of both E 1 and E 2 are achieved at approximately α 1 = 0.358, and so the coefficients are α 1 = 0.358

α 2 = -0.47710242361717810834 α 3 = 0.35230499471528197958 α 4 = 0.26679742890189612876 with E 1 (α) = 2.9084, E 2 (α) = 3.1527.
Case s = 5. Now the composition

(3.15) ψ [5] δt = χ α 1 δt •χ α 2 δt •χ α 3 δt •χ α 4 δt •χ α 5 δt •χ α 5 δt •χ α 4 δt •χ α 3 δt •χ α 2 δt
•χ α 1 δt involves 21 maps when applied to a system separable into three parts. By carrying out a similar analysis we conclude that the best solution according with the criterion adopted is achieved when

α 1 = α 2 = α 3 = α 4 = 1 2(4 -4 1/3 ) , α 5 = 1 2 -4α 1 , which give E 1 (α) = 2.3159, E 2 (α) = 2.6111. Remark 3.3.
-This method also belongs to G 2 since it can be written as

ψ δt = φ γδt • φ γδt • φ βδt • φ γδt • φ γδt belonging to G 2 with coefficients γ = 2α 1 , β = 2α 5 .
This method was originally presented in [START_REF] Suzuki | Fractal decomposition of exponential operators with applications to many-body theories and Monte Carlo simulations[END_REF].

Case s = 6. Analogously we have considered a composition involving three free parameters and 25 maps:

(3.16) ψ [6] δt = χ α 1 δt •χ α 2 δt •χ α 3 δt •χ α 4 δt •χ α 5 δt •χ α 6 δt •χ α 6 δt •χ α 5 δt •χ α 4 δt •χ α 3 δt •χ α 2 δt •χ α 1 δt .
A solution leading to small values of E 1 and E 2 is

α 1 = α 2 = 3 20 α 3 = 17 100 α 4 = -0.2628463256938681137 α 5 = 0.16217658484020533783 α 6 = 0.13066974085366277593 with E 1 (α) = 2.0513, E 2 (α) = 2.4078. Remark 3.4.
-Although the optimization criterion we have adopted here usually leads to good methods, one can find schemes in the literature with larger values of E 1 and E 2 which are very efficient in practice. Thus, in particular, we mention the fourth-order splitting method designed in [START_REF] Blanes | Practical symplectic partitioned Runge-Kutta and Runge-Kutta-Nyström methods[END_REF] which, once written as a method in G 1 , also involves s = 6 stages.

Numerical results

In this section, we show some numerical results to illustrate the efficiency and performance of the methods previously derived. We focus on Vlasov applications by considering the Vlasov-HMF and the Vlasov-Maxwell system.

4.1. Vlasov-Maxwell system.-The composition methods introduced in the previous sections can then be used to derive a global 4th order method for the Vlasov-Maxwell equation. As an example, the Yoshida (or triple-jump) method (s = 3) in the Vlasov-Maxwell context writes

ψ [3] δt = χ α 1 δt • χ α 2 δt • χ α 3 δt • χ α 3 δt • χ α 2 δt • χ α 1 δt , with α 1 = α 2 = 1 2(2-2 1/3 ) , α 3 = 1 2 -2α 1
, and where χ δt and χ δt are given by (3.3) and (3.4). Then, if we denote by F n an approximation at time t n = nδ t , n ∈ N of the Vlasov-Maxwell solution F (t n ), we have

F n = ψ [3] δt n (F in ),
and F n is a 4th order approximation of F (t n ). The other 4th order methods ψ δt which, with our notations writes

ψ [2] δt = χ δt/2 • χ δt/2 = ϕ [H E ] δt/2 • ϕ [H f ] δt/2 • ϕ [H B ] δt/2 • ϕ [H B ] δt/2 • ϕ [H f ] δt/2 • ϕ [H E ] δt/2 = ϕ [H E ] δt/2 • ϕ [H f ] δt/2 • ϕ [H B ] δt • ϕ [H f ] δt/2 • ϕ [H E ] δt/2 .
The Strang splitting for a decomposition into three parts involves 5 maps since, as usual, the first and the last maps can be concatenated.

We present some numerical results to illustrate the efficiency of the different methods. First of all, we used the methods ψ

[s]
δt , s = 2, 3, 4, 5, 6. In this context, one goal is to compare the new exact splitting for the rotation applied to the field H B to a standard Strang method. In the methods ψ given by (2.10). This means that in this method, χ δt is now replaced by χδt defined by χδt = ϕ

[H E ] δt • ϕ [H f ] δt • T [H B ] δt .
The global Strang splitting is then defined by ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt

ψ[2] δt = χδt/2 • χ δt/2 ,
and the definition of ψ[s] δt for s = 3, 4, 5, 6 follows directly. Let us remark that even if the magnetic part H B is not solved exactly in time, the global method ψ [s] δt still have the same order as ψ δt (i.e. of order 2 for s = 2 or of order 4 for s = 3, 4, 5, 6). We then want to investigate the impact of this approximation on the global error of the so-obtained splitting.

To do so, we consider the following initial condition for (3.1)

f in (x 1 , v 1 , v 2 ) = 1 πv th 2 √ T r e -(v 2 1 +v 2 2 /Tr)/v th (1 + α cos(kx 1 )),
and B in (x 1 ) = 10 + 3 cos(kx 1 ), E in 2 (x 1 ) = 0. We consider α = 10 -4 , k = 0.4, v th = 0.02, k = 0.4 and T r = 12. The phase space domain is (x 1 , v 1 , v 2 ) ∈ [0, 2π/k] × [-1, 1] 2 and the number of points is N x = 8 in space and N v = 513 per direction in velocity. The runs are performed up to a final time T = 2 and different values of the time step δ t are considered between 10 -3 to 0.4. The results are given in Figures 8 and9, where we have plotted the L ∞ error on the total energy with respect to δ t /M , where M is the number of maps. The total energy (which is conserved with time at the continuous level) is defined by (4.1)

H(t) = 1 2 L |E| 2 dx + 1 2 L |B| 2 dx + 1 2 L×R 2 |v| 2 f dvdx with L = [0, 2π/k],
H(t) -H(0) H(0) .
First, one can see that the order of convergence is well recovered for all the methods but some fourth order methods present some better efficiency. For instance, the two methods corresponding to s = 5 and s = 6 are clearly the best, and are much more efficient that the triple jump method (s = 3) or the Strang one (s = 2) even if they involve a larger number of maps. Second, we can observe than the error produced by the methods ψ δt (i.e. when the exact splitting is used for the part H B ) is smaller than the error performed by the methods ψ[s] δt (i.e. when a Strang splitting is used for the part H B ). Note that on Figures 8 and9, the lines indicating the order are kept fixed. Moreover, on Figure 9, in addition to the error produced by the methods ψ[s] δt , the error curves produced by the methods ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt and ψ [5] δt (labelled by s = 2 new (5) and s = 5 new (21) in the legend) have been displayed to ease the comparisons. Note that we have chosen to plot the ψ [5] δt method but very similar conclusions arise with the choice ψ [START_REF] Besse | Convergence of classes of high-order semi-Lagrangian schemes for the Vlasov-Poisson system[END_REF] δt since the two methods ψ [5] δt and ψ [START_REF] Besse | Convergence of classes of high-order semi-Lagrangian schemes for the Vlasov-Poisson system[END_REF] δt have a very close efficiency in our context. For the global Strang method the ratio between the error produced by ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt and ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt is about 2.5 whereas the ratio between the error produced by ψ [5] δt and ψ [5] δt is about 6 (the same ratio is observed between ψ [START_REF] Besse | Convergence of classes of high-order semi-Lagrangian schemes for the Vlasov-Poisson system[END_REF] δt and ψ [START_REF] Besse | Convergence of classes of high-order semi-Lagrangian schemes for the Vlasov-Poisson system[END_REF] δt ). Let us remark that, for a given method (i.e. a given s), the cost of a ψ[s] δt method is the same as a the one ψ

[s] δt method. We end this subsection by considering other splitting methods from the literature, namely the splitting methods of the form (3.6) from [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] which assume that each subpart is solved exactly, which is our case when the exact splitting is used for the magnetic part. The results are displayed in Figure 10 where we have tested second order methods (AK 3-2 and AK 5-2 involve 9 maps), a fourth order method (AK 11-4 involves 21 maps) and even a sixth order method (AY 15-6 involves 29 maps). We refer to [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] for more details on these methods. As previously we also added ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt (second order) and ψ [5] δt (fourth order) for comparison, whereas the slope 2 and 4 are the same as in Figures 8 and9. First, we observe that AK 3-2 is the best second order method. The third order PP method is not very attractive in this context compared to second order methods. Second, among the two fourth order methods (AK 11-4 and ψ [5] δt ), the method ψ [5] δt offers a better efficiency since the error is about 5 times smaller. Finally, the method AY 15-6 offers sixth order accuracy but this extra accuracy is only apparent for very small time steps. 4.2. Vlasov-Maxwell system: long time test. -We now present a test to highlight the fact that the new methods are able to capture the long time dynamics of the Vlasov-Maxwell solution. Then, we consider the same initial condition as in the previous test δt , s = 2, 3, 4, 5, 6 for the Vlasov-Maxwell system. The number of maps for each method is indicated into parenthesis. as B in (x 1 ) = 1 + 0.0001 sin(x 1 ). The number of points in space is N x = 32, whereas we took N v = 257 points per velocity direction.

f in (x 1 , v 1 , v 2 ) = 1 2πβ e -v 2 2 /β e -(v 1 -0.1) 2 /β + e -(v 1 +0.3) 2 /β , where (x 1 , v 1 , v 2 ) ∈ [0, 2π] × [-1, 1]
We compare the ψ [5] δt method (which is the best method according to the previous tests) with the two second order splittings ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt (referred as 's = 2 new') and ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt (referred as 'Strang'). Let us recall that these two second order splittings only differ in the solving of the magnetic part. The time step is chosen as δ t = 0.125 or δ t = 0.025 so that we will compare ψ [5] δt and ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt at a fixed computational cost and the final time is t f = 500.

In Figure 11 (left), we plot the time evolution of the relative total energy given by (4.2) for the two second order splittings with a small time step (δ t = 0.025). The cost of these two methods is the same, but we can see that the relative total energy is better preserved for the ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt method (about 2 × 10 -9 ) compared to the standard ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt method (about 2 × 10 -6 ). On the right part of Figure 11, we compare the methods for which the rotation is solved exactly, namely ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt (with δ t = 0.025 and 0.125) and ψ [5] δt (with δ t = 0.125). Let us mention that ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt for δ t = 0.025 and ψ [5] δt for δ t = 0.125 have the same number of stages so that the same computational cost. We can observe that the high δt , ψ [5] δt for the Vlasov-Maxwell system. The order lines 'slope 2' and 'slope 4' are the same as in Figure 8. The number of maps for each method is indicated into parenthesis. order method ψ [5] δt preserves very well the total energy (about 6 × 10 -12 ) which confirms the results obtained in the previous subsection.

Finally, the time history of first mode of E 1 is plotted in Figure 12, for ψ

[2] δt and ψ [5] δt method, with δ t = 0.125. We can observe that after a linear phase during which the amplitude of the mode grows exponentially, a saturation phase is well captured by the two methods, even if the saturation level is not the same (see the small figure with a zoom for t ∈ [400, 500]). Refining the time step by considering δ t = 0.025 enables the ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt method to recover the results obtained by ψ [5] δt .

4.3. Vlasov-HMF system.-Our goal is to solve numerically the Vlasov-HMF model satisfied by f (t, x, v), (x, v) ∈ L×R, with L = R/2πZ (see [START_REF] Barré | The Vlasov equation and the Hamiltonian Mean-Field model[END_REF][START_REF] Horsin | Comportement en temps long d'équations de type Vlasov : études mathématiques et numériques[END_REF]) δt , ψ [5] δt for the Vlasov-Maxwell system. The order lines 'slope 2' and 'slope 4' are the same as in Figure 8. The number of maps for each method is indicated into parenthesis. δt and ψ [START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] δt with δ t = 0.125. Right:

(4.3) ∂ t f + {f, H[f ]} = 0,
ψ [2]
δt with δ t = 0.125 and δ t = 0.025, ψ where {f, g}

= ∂ x f ∂ v g -∂ v f ∂ x g and H[f ] is given by H[f ] = v 2 2 -Φ[f ](x).
Finally, the potential is defined by We consider the following stationary solution (see [START_REF] Horsin | Comportement en temps long d'équations de type Vlasov : études mathématiques et numériques[END_REF] for more details)

(4.5) f eq (x, v) = γe -β v 2 2 -M 0 cos x with M 0 = L×R cos(y)f eq (y, u)dydu,
where γ, β, M 0 ∈ R will be explicitly given below. Following [START_REF] Horsin | Comportement en temps long d'équations de type Vlasov : études mathématiques et numériques[END_REF], the long time behavior of (4.3) is driven by the linearized Hamiltonian part, i.e. ∂ t f + {f, H[f eq ]} = 0, with H[f eq ] = v 2 2 -M 0 cos(x). We recognize the pendulum Hamiltonian for which a slight modification of the new splitting is able to capture the rotation phenomena with high accuracy compare to standard Strang splitting (see [5]). In this HMF context, the material introduced before has to be slightly modified.

First, let us introduce the discretization of the phase space L×[-v max , v max ], with v max > 0 a truncation of the velocity direction. We consider G x := h x 0, N x -1 the space grid (with h x = L/N x the stepsize and N x ∈ N * the number of points) and G v := h v -(N v -1)/2 , N v /2 the speed grid (with h v = 2v max /N v the stepsize and N v ∈ N * the number of points). We also introduce the set of discrete frequencies:

G x = η x -(N x -1)/2 , N x /2 and G v = η v -(N v -1)/2 , N v /2 with η x = 2π/L and η v = π/v max .
Then, we define the discrete partial Fourier transforms

F 1 :    C Gx×Gv → C Gx×Gv u → h x g 1 ∈Gx u g 1 ,g 2 e -ig 1 ξ 1 and F 2 :    C Gx×Gv → C Gx× Gv u → h v g 2 ∈Gv u g 1 ,g 2 e -ig 2 ξ 2 ,
whereas the shears are now defined by

S α 1 : C Gx×Gv → C Gx×Gv u → F -1 1 e iαξ 1 g 2 F 1 u and S α 2 : C Gx×Gv → C Gx×Gv u → F -1 2 e iαξ 2 E[u]g 1 F 2 u , where E[u] g 1 is deduced from the relation E[u](x) = -∂ x Φ[u](x) and (4.4) E[u] g 1 = sin(g 1 h x )h x h v (g 1 ,g 2 )∈Gx×Gv cos(g 1 h x ) u g 1 ,g 2 -cos(g 1 h x )h x h v (g 1 ,g 2 )∈Gx×Gv sin(g 1 h x ) u g 1 ,g 2 . (4.6)
Then, at time t n = nδ t , we denote by f n an approximation of the solution f (t n ) on the phase space grid computed by the Strang splitting T δt and the new splitting M δt which are defined by (4.7)

f n+1 = T δt f n := S -δt/2 1 S δt 2 S -δt/2 1 f n , (Strang) f n+1 = M δt f n := S -tc tan(δt/(2tc)) 1 S tc sin(δt/tc) 2 S -tc tan(δt/(2tc)) 1 f n , (New)
where f 0 := f in , and

t c = 1 √ M 0
is the characteristic time of the Vlasov-HMF model which has been introduced to capture the correct angular velocity. Let us remark that the electric field E[f ] has to be solved using (4.6) before the shear S α 2 in the splittings (4.7). To evaluate the performance of the new splitting compare to the Strang one, we consider an initial condition f in as a perturbation of the equilibrium solution (4.5) (with β = 10, M 0 = 0.9455421864232981 and α = 0.0001194365987897421) f in (x, v) = f eq (x, v)(1 + ε cos(x)), (x, v) ∈ [-π, π] × R, with ε = 10 -3 . We consider a truncated velocity domain of [-8, 8], the number of points in the spatial direction is N x = 128 whereas we considered N v = 256 points in the velocity direction, and the final time is T = 25. Note that the splitting can also be coupled to a semi-Lagrangian method; the shears S 1 and S 2 have to be modified accordingly (see [START_REF] Besse | Convergence of classes of high-order semi-Lagrangian schemes for the Vlasov-Poisson system[END_REF] for instance). We look at the L ∞ error between a reference distribution function (obtained with the new splitting with a small time step δ t = T /1000) and the one obtained by Strang or new splitting given by (4.7) (with t c = 1.0283940255) for different time steps δ t ∈ {T /50, T /100, T /150, T /200, T /250}. The results are displayed in Figure 13 in log-log scale. First we observe that, as expected, the two methods are second order accurate in time. But, one can remark that the error produced by the new splitting is much more smaller that the error produced by the Strang splitting, at the same cost (the number of maps is the same for the two methods).

Conclusion

In this work, we have studied a directional splitting which preserves exactly the rotations and apply to the PDE context. A careful numerical analysis of this splitting coupled with spectral interpolation techniques has been performed. These results are illustrated by some numerical experiments.

Then, this step serves as a building block of a splitting for the Vlasov-Maxwell system. Indeed, this system can be split into three parts which, thanks to this new splitting, can all be solved exactly. New high order composition methods are then designed to accurately and efficiently solve the full Vlasov-Maxwell system. Numerical results show the good behavior of these methods. Finally, a direct application of the new splitting for close to equilibrium simulations of the Vlasov-HMF model also shows very good results compared to the standard Strang splitting, with no additional cost.

The extension to the relativistic Vlasov-Maxwell equations in two or three dimensions in the velocity space are planned. The approach should be even more attractive in this context. To conclude, we just have to verify that P a,b is not the identity. First, we observe that P a,b is not a scalar matrix, so there exists x ∈ R 2 such that x is not an eigenvector of P a,b . Then, we consider a vector y ∈ R 2 \ {0} such that x and y are orthogonal. By construction, we have

|x| s (u•τ ) L 2 (R 2 ) = | det τ | -1 |τ (x)| s u L 2 (R 2 ) ≤ | det τ -1 ||τ | s u X s ≤ κ s+1 u X s .
6.5. Proof of Lemma 2.11. -First, we apply the Poisson formula and the discrete Fourier Plancherel isometry to get

v |G 2 L 2 (G 2 ) ≤ v |hZ 2 L 2 (hZ 2 ) = 1 2π k∈Z 2 F v(• + 2kπ h ) L 2 ((-π h , π h ) 2 ) .
Then we observe that if k ∈ Z 2 \ {0} and ξ ∈ (-π h , π h ) 2 then

|ξ + 2kπ h | ≥ π h (2|k| - √ 2).
Thus, we control v |G 2 L 2 (G 2 ) by

1 2π F v L 2 ((-π h , π h ) 2 ) + h 2 2π 3 k∈Z 2 \{0} 1 (2|k| - √ 2) 2 (|ξ| 2 F v)(• + 2kπ h ) L 2 ((-π h , π h ) 2 ) .
Finally, applying the Cauchy Schwarz inequality and the Chasles relation, we control the second term by

h 2 2π 3 |ξ| 2 F v L 2 (R 2 ) k∈Z 2 \{0} 1 (2|k| - √ 2) 4 .

( 2 . 18 ) 4 )

 2184 L a,b = µ a,b with µ a,b = F (ab(1 -ab/4)), where F is the continuous function on (-∞, 1] defined by

Lemma 2 . 5 .

 25 -If a, b are small enough, there exists a symmetric matrix L a,b such that L a,b goes to 0 as (a, b) goes to 0 and P a,b = exp(JL a,b ). Then, we have to determine a formula for L a,b . Since L a,b is a symmetric matrix, e JL a,b is a Hamiltonian flow at time 1. A fortiori, L a,b is a constant of motion. So we have t (e JL a,b )L a,b e JL a,b = L a,b .

|G 2

 2

Figure 1 .

 1 Figure 1. Spatial error (log -log scale) as a function of the number of points N between the exact shifted solution and the approximation computed using fast Fourier transform.

Figure 2 .

 2 Figure 2. Zoom of Figure 1.

Figure 3 .Figure 4 .

 34 Figure 3. Solution u(T, x) of (2.28). Left: Exact solution u(T, x). Middle: Numerical solution obtained by the new splitting. Right: Numerical solution obtained by the Strang splitting.

Figure 5 .Figure 6 .

 56 Figure 5. Time history of the relative errors (zoom of Figure 4 around T ≈ 3188 (left) and T ≈ 2 × 3188 (right).

3. 1 .Figure 7 .

 17 Figure 7. Time execution as a function of the total number of points (log -log scale). Blue: new method (new splitting and onedimensional fast Fourier transform). Red: exact computation of the feet of the characteristics and two-dimensional non uniform fast Fourier transform.

  δt are available, as long as one is able to construct the corresponding adjoint χ δt . Remark 3.1. -Another well-known class G 2 of integrators is formed by compositions (3.10)

  δt , s = 4, 5, 6 are defined by (3.14), (3.15) and (3.16) in Subsection 3.3.2. We also define the standard Strang splitting ψ[START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF]

  the Strang splitting T [H B ] δt

2 Figure 8 .

 28 Figure 8. Efficiency diagrams obtained by different composition methods ψ [s]

Figure 10 .

 10 Figure 10. Efficiency diagrams obtained by different methods from[START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF] and ψ[START_REF] Auzinger | Practical splitting methods for the adaptive integration of nonlinear evolution equations. Part I: Construction of optimized schemes and pairs of schemes[END_REF]

Figure 11 .

 11 Figure 11. Time history of the relative total energy. Left: ψ [2]

  δt with δ t = 0.125.

Figure 12 .

 12 Figure 12. Time history of the first mode of E 1 (semi-log scale): ψ [2]δt and ψ[5] δt with δ t = 0.125.

( 4 . 4 )

 44 Φ[f ](x) = cos x L×R cos(y)f (y, u)dydu + sin x L×R sin(y)f (y, u)dydu.

Figure 13 .

 13 Figure 13. Error as a function of the number of iterations for the HMF-Poisson system. Comparison of the Strang splitting ('strang') and the New splitting ('new').

6. 2 . 6 . 3 . 2 (

 2632 Proof of Lemma 2.5. -Notice that, if a and b are small enough, P a,b defined in (2.19 is close to the identity. Consequently, it admits a logarithm M a,b ∈ M 2 (R), defined by M a,b = n∈N (-1) n n + 1 (P a,b -I 2 ) n , and satisfying e M a,b = P a,b . A fortiori, we have exp(Tr M a,b ) = det P a,b = 1. Hence we have Tr M a,b = 0.Furthermore, the following application define an isomorphism of vector spaces (it is an injection between two spaces of dimension 3)S 2 (R) → sl 2 (R) L → JL , where sl 2 (R) = {M ∈ M 2 (R) | Tr M = 0}. As a consequence, there exists a symmetric matrix L a,b ∈ S 2 (R) such that M a,b = JL a,b . Proof of Lemma 2.6. -Since 0 < ab < 4, Q a,bis either positive or negative, and, as a consequence, the following Euclidean norm is well defined on S 2 (R)∀K ∈ S 2 (R), K 2 a,b := R t xKx) 2 e -| txQ a,b x| dx. Since det P a,b = 1, computing R -1 a,b K a,b , we deduce from a change of variables and from (2.21) that ∀K ∈ S 2 (R), R a,b K a,b = K a,b . This relation means that R a,b is an isometry for the Euclidean norm • a,b . A fortiori, we have det R a,b = ±1. But, since R 0,0 = I 2 and (a, b) → det R a,b is a continuous map, we deduce that det R a,b = 1. Consequently, R a,b is a rotation in a space of dimension 3. So, there are only two possibilities: either R a,b is the identity or the eigenspace of R a,b associated with the eigenvalue 1 is of dimension 1.

t

  yP a,b x = 0.Consequently, if K = y t y ∈ S 2 (R), we havet xR a,b (K)x = ( t yP a,b x) 2 = 0 = ( t yx) 2 = t xKx. Thus, we have R a,b (K) = K. 6.4. Proof of Lemma 2.9. -We have to bound |x| s (u • τ ) L 2 (R 2 ) and |ξ| s F (u • τ ) L 2 (R 2 ). However, a straightforward calculation shows thatF (u • τ ) = |det τ | -1 (F u) • t τ -1 ,and equation (2.24) is clearly equivalent to |τ | ≤ κ and |τ -1 | ≤ κ. Thus, since |τ | = | t τ |, if we get a bound on |x| s (u • τ ) L 2 (R 2 ) , uniform with respect to τ , we also a a bound on |ξ| s F (u • τ ) L 2 (R 2 ) uniform with respect to τ . Finally, to bound |x| s (u • τ ) L 2 (R 2 ) , we just have to apply a change of coordinates:

Acknowledgements. -FC has been supported by the Ministerio de Economía y Competitividad (Spain) through project MTM2016-77660-P (AEI/FEDER, UE), by the Mobility Grant PRX18/00145 and by the Centre Henri Lebesgue, whereas on sabbatical leave. He would like to express his gratitude to the Université de Rennes 1 for its hospitality. NC is supported by the French ANR project MOONRISE ANR-14-CE23-0007-01 and by the Enabling Research EUROFusion project MAGYK 2019-2020. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Appendix

In this Appendix, we gather the proofs of the different lemmas used in the proof of convergence of the pseudo-spectral splitting methods.

Else we have

Else if

This last quantity is a second order polynomial with respect to y 2 . Thus its infimum can be determined explicitly. More precisely, we have