High-order Hamiltonian splitting for Vlasov-Poisson equations

We consider the Vlasov-Poisson equation in a Hamiltonian framework and derive new time splitting methods based on the decomposition of the Hamiltonian functional between the kinetic and electric energy. Assuming smoothness of the solutions, we study the order conditions of such methods. It appears that these conditions are of Runge-Kutta-Nystr{\"o}m type. In the one dimensional case, the order conditions can be further simplified, and efficient methods of order 6 with a reduced number of stages can be constructed. In the general case, high-order methods can also be constructed using explicit computations of commutators. Numerical results are performed and show the benefit of using high-order splitting schemes in that context. Complete and self-contained proofs of convergence results and rigorous error estimates are also given.


Introduction
Frequently, the Vlasov equation is solved numerically with particles methods. Even if they can reproduce realistic physical phenomena, they are well known to be noisy and slowly convergent when more particles are considered in the simulation. To remedy this, the so-called Eulerian methods (which use a grid of the phase space) have known an important expansion these last decades. Indeed, due to the increase of the machines performance, the simulation of charged particles by using Vlasov equation can be performed in realistic configurations. However, these simulations are still computationally very expensive in high dimensions and a lot has to be done at a more theoretical level to make simulations faster. For example, the use of high-order methods is classical when one speaks about space or velocity discretization. However, for the simulation of Vlasov-Poisson systems, the use of high-order methods in time is not well developed; generally, only the classical Strang splitting is used and analyzed; see however pioneering works of [23,19] following [24] or the recent work of [21] in the linear case. We mention also the work [11], which tells us that the increase of order of discretization in space should be followed with an increase of order in time.
On the other side, a literature exists around the construction of highorder methods for ODE (see [5,4,13,22]). The main goal of this work is to construct high-order splitting schemes for the nonlinear Vlasov-Poisson PDE system by the light of these recent references.
The starting point of our analysis relies on the fact that the Vlasov-Poisson equation is a Hamiltonian PDE for a Lie-Poisson bracket common to several nonlinear transport equations appearing in fluid dynamics, see for instance [15] and Section 2 below. Up to phase space discretization, the splitting schemes we construct preserve this structure and hence are geometric integrators in the sense of [13,14]. Moreover, each block is explicit in time, and can be used to derive high-order methods, taking into account some specific commutator relations.
We consider the following equation Vlasov-Poisson equation where f (t, x, v) depends on time t ≥ 0 and the phase space variables (x, v) ∈ T d ×R d , d = 1, 2, 3, and where for vectors (x 1 , . . . , x d ) ∈ T d and (y 1 , . . . , y d ) ∈ R d , we set x · y = x 1 y 1 + · · · + x d y d and |y| 2 := y · y. Here, T d denotes the d dimensional torus R d /(2πZ d ) which means that the domain considered is periodic in space. Note that formally, the same analysis is valid on more general domains, however, we will perform the analysis, in particular the convergence analysis in this simplified framework. Classically, the variable x corresponds to the spatial variable whereas v is the velocity variable. The electric potential φ(f ) solves the Poisson equation x i is the Laplace operator in the x variable acting of functions with zero average. The electric field depending only on x is defined as E(f ) = −∂ x φ(f ). The energy associated with equations (1.1)- (1.2) is The time discretization methods proposed in this paper are based on this decomposition of the energy. Indeed, the solution of the equations associated with T and U can be solved exactly (up to a phase space discretization, for example by interpolation in the framework of semi-Lagrangian methods). We denote by ϕ t T (f ) and ϕ t U (f ) the flows associated with T and U respectively (we postpone to Section 5 the precise definition of Hamiltonian flows). The first one corresponds to the equation for which the solution is written explicitly as For the flow ϕ t U , we have to solve the equation for which we verify that the solution is given by where E(f (0)) is the value of the electric field at time t = 0. Indeed, φ(f ) is constant along the solution of (1.4). Based on these explicit formulae, we will first consider numerical integrators of the form are real coefficients, and such that the numerical solution after time t = τ coincides with the exact solution up to terms of order τ p , i.e., for a given smooth function f , where ϕ τ H (f ) corresponds to the exact flow associated with (1.3). We will give a precise definition of smoothness in Section 5, and show that this condition ensures the convergence of order p of the numerical method.
As composition of exact flows of Hamiltonians T and U, such schemes are (infinite dimensional) Poisson integrators in the sense of [13,Chapter VII]. In particular they preserve the Casimir invariants for the structure for all times (e.g. all the L p norms of f ). Note that in this work, we do not address the delicate question of phase space approximation and focus only on time discretization effects (see [2,6,20]).
To analyze the order of the schemes (1.5), we will use the Hamiltonian structure of the flows. We will show that they can be expanded in suitable function spaces in terms of commutators, formally reducing the problem to classical settings based on the Baker-Campbell-Haussdorff (BCH) formula and the Lie calculus, see for instance [13,4]. A rigorous justification will be given in Section 5.
In the Vlasov-Poisson case, we will see that the functionals T and U in the decomposition (1.3) satisfy the following formal relation where [·, ·] is the Poisson bracket associated with the infinite dimensional Poisson structure (see Section 2). This property reduces the number of order conditions on the coefficients (a i ) s i=1 and (b i ) s+1 i=1 in formula (1.5). The situation is analogous to the case of splitting methods of Runge-Kutta-Nyström (RKN) type for ordinary differential equations (ODE) derived from a Hamiltonian function, see [13,4]. In dimension d = 1, the Vlasov-Poisson system even satisfies the stronger property where m is the total mass of f which is a Casimir invariant, preserved by the exact flow and the splitting methods (1.5). This means that we have naturally simpler algebraic order conditions than those of RKN type for the specific Vlasov-Poisson system in dimension 1. In any dimension, it also turns out that the exact flow of the Poisson bracket [[T , U], U] can be computed up to space discretization. We will retain these ideas to derive new high-order splitting integrators involving also the flow of this nested Poisson bracket with optimized coefficients and number of internal steps, in a similar way as in the ODE setting [5,4]. The paper is organized as follows: • In Section 2, we discuss the Hamiltonian Lie-Poisson structure of the Vlasov-Poisson equation, and give the expressions of some iterated Poisson brackets. They will form the cornerstone of our analysis. • In Section 3, we will first make the link between the standard Lie calculus and the Hamiltonian structure, and then derive high-order splitting methods based on the formula (1.5). We will then consider generalizations of these methods using explicitly calculable flows of iterated brackets. • In Section 4 we give numerical illustrations of the performances of the methods: we mainly exhibit the order of the method, but also address the question of Casimir invariant preservation (e.g. the L p norms), regarding the influence of phase space discretizations. • Finally, Section 5 is devoted to the mathematical analysis of the splitting method: we give convergence results in some function spaces.
To this aim, we first give a local existence result of the Vlasov-Poisson equation with precise estimates (following in essence [8]), then prove some stability estimates. The results presented in this section can be compared with the one in [10] for the Strang splitting, where however only compactly supported solutions are considered.

Hamiltonian structure
2.1. Poisson brackets. We define the microcanonical bracket {f, g} of two (sufficiently smooth) functions by the formula With this notation, we can write the Vlasov-Poisson equation as is the microcanonical Hamiltonian associated with f . Recall that for a given functional G(f ), its Fréchet derivative is the distribution δG δf (f ) evaluated at the point f , being defined by the formula for all smooth variation δf . In general, a Fréchet derivative is an operator acting on f , hence a rigorous writing of the previous formula necessitates a loss of derivative in f . We will discuss these issues in Section 5. Considering the two functionals T (f ) and U(f ) defined in (1.3), their Fréchet derivatives read explicitly where φ(f ) is given by (1.2). Due to the relation H = T + U, the Vlasov-Poisson equation can be written as which is a Hamiltonian equation for the Poisson structure associated with the following Poisson bracket: for two functionals H(f ) and G(f ), we set where the Fréchet functionals are evaluated in f . Note that the skewsymmetry is obtained using the relation We refer to [15] for discussions related to this structure. Note that to give a meaning to all the previous expressions, we usually have to assume smoothness for the function f and deal with loss of derivatives, see for instance (5.5) of Section 5. The Hamiltonian-Poisson structure defined above possesses Casimir invariants, meaning quantities preserved for every Hamiltonian system of the form (2.1), and not depending on the specific form of H. This is essentially a consequence of the fact that the nonlinear transport equation (2.1) involves divergence free vector fields. Let ψ : R → R be a smooth function, and consider the functional and hence the potential energy U can be written is the kernel of the inverse of the Laplace operator in the x variable.
The aim of this subsection is to prove the following result: For any smooth function f , the functionals T (f ) and U(f ), satisfy the following relation: is a constant of motion of (1.1), and where φ(f ) is defined in (1.2). In dimension d = 1, we have V(f ) = 0, and in any dimension d ≥ 1, the relation holds for all functions f .
Proof. Using (2.2), we calculate the following Let us calculate the Fréchet derivative of this functional. To this aim, we evaluate this functional at f + δf , where δf stands for a small perturbation f satisfying Hence, we get We see that, using an integration by parts in x, the third term can be written as − We deduce that . Using this relation, we calculate Now we see that the term involving the function Z(f )(x) vanishes, as the integral of ∂ v f (x, v) in v ∈ R d is equal to 0. We can thus write after an integration by parts In other words, we get (2.10) , and this yields (2.6). In dimension d = 1, we can further write that and conclude that V is identically equal to 0. In any dimension d, as the Fréchet derivatives of U and V depend only on x, we automatically obtain (2.8).

Flow of the Hamiltonian
where K satisfies Denoting by E x j , j = 1, · · · , d the components of the electric vector fields E, we get in the case d = 2 and in the case d = 1, where we used (2.9). Hence we have where E x j , j = 1, · · · , d, denote the components of the electric vector fields which can be solved, the right hand side being of zero average. Let us remark that the two last term are nothing but the determinant of the Jacobian matrix of the electric field E.
With the previous notations and Proposition 2.2, the equation associated with the Hamiltonian [[T , U], U] is given by Hence the flow associated with [[T , U], U] is explicit and given by This notation is justified by the fact that the equation (1.1) is equivalent to where G here are functionals acting on some function space. We will not discuss here the mathematical validity of such an equivalence, but taking for instance G(f ) as a norm of a function space (see Section 5), we can prove that the solution f (t) admits a formal expansion of the form By using similar expansions for the flows we have By using (3.1) and the Jacobi identity, we see that the following relation holds true. We deduce that the classical calculus of Lie derivatives also applies to our case.
Using this identification, (see also [4]) we can write formally the exact flows as ϕ t T =: e tT , ϕ t U =: e tU , and ϕ t H = e t(T +U ) , with the operators T and U satisfying the relation To derive splitting methods in dimension d ≥ 2, we will also consider numerical schemes containing blocks based on the exact computation of the flow associated with the Hamiltonian [[T, U ], U ]. In this section we will concentrate on the derivation of high-order splitting methods of the form (1.5) satisfying these formal relations.
Scheme (1.5) is at least of order 1 for the problem (1.1) if and only if the coefficients a i , b i satisfy the consistency condition We are mainly interested in symmetric compositions, that is, integrators In that case, they are of even order. In particular, a symmetric method verifying (3.3) is at least of order 2 [13]. Notice that the number of flows in the splitting method (1.5) or (3.4) is σ ≡ 2s + 1, but the last flow can be concatenated with the first one at the next step in the integration process, so that the number of flows ϕ τ U and ϕ τ T per step is precisely s. Restriction (1.6) imposes a set of constraints the coefficients a i , b i in the composition (3.4) have to satisfy. These are the so-called order conditions of the splitting method and a number of procedures can be applied to obtain them [13]. One of them consists in applying recursively the Baker-Campbell-Hausdorff (BCH) formula in the formal factorization (3.4). When this done, we can express ψ τ p as the formal exponential of only one operator (3.5) ψ τ p = e τ (T +U +R(τ )) , where The integrator is of order p if R(τ ) = O(τ p ), and thus the order conditions are p 21 = p 31 = p 32 = · · · = 0 up to the order considered. For a symmetric scheme one has R(−τ ) = R(τ ), so that R(τ ) only involves even powers of τ . In consequence, p 21 = p 41 = p 42 = · · · = p 2n,k = 0 automatically in (3.6) and we have only to impose p 31 = p 3,2 = · · · = p 2n+1,k = 0. The total number of order conditions can be determined by computing the dimension of the subspaces spanned by the k-nested commutators involving T and U for k = 3, 5, . . ., see [17].
For higher order terms in R(τ ). The number of order conditions is thus correspondingly reduced. In Table 1 we have collected this number for symmetric methods of order p = 2, 4, 6, 8 (line d > 1). Thus, a symmetric 6th-order scheme within this class requires solving 8 order conditions (the two consistency conditions (3.3) plus 2 conditions at order 4 plus 4 conditions at order 6), so that the scheme (3.4) requires at least 15 exponentials. In fact, it is a common practice to consider more exponentials than strictly necessary and use the free parameters introduced in that way to minimize error terms. In particular, in [5] a 6th-order splitting method involving 23 exponentials (11 stages) was designed which has been shown to be very efficient for a number of problems, including Vlasov-Poisson systems [12]. We have shown in the subsection 2. ψ τ p = e τ C 1 e a 1 τ T e τ C 2 · · · e τ C 2 e a 1 τ T e τ C 1 . In that case the order conditions to achieve order 6 are explicitly lead to a method of order four. With the inclusion of C i in the scheme, the number of exponentials can be significantly reduced (one has more parameters available to satisfy the order conditions): the minimum number of exponentials required by the symmetric method (3.7) to achieve order 6 is σ = 9 instead of 15 for scheme (3.4). There are several other systems where the evaluation of the flow associated with [[T, U ], U ] is not substantially more expensive in terms of computational cost than the evaluation of e τ U , and thus schemes of the form (3.7) have been widely analyzed and several efficient integrators can be found in the literature [4,18].
We have considered compositions of the form (3.7) with σ = 9, 11, and 13 exponentials. When σ = 9 there is only one real solution of equations (3.8). More efficient schemes are obtained by using more exponentials: the corresponding free parameters can be used to optimize the scheme (for instance, by annihilating higher order terms in R(τ ), reducing the norm of the main error terms, etc.). In Table 2 we collect the coefficients of the best methods we have found. The most efficient one (see Section 4) corresponds to σ = 13. In this case the two free parameters have been chosen to vanish the coefficient multiplying the commutator [

T, [T, [T, [T, [T, [T, U ]]]]]
] at order 7 and such that b 1 = b 2 . This procedure usually leads to very efficient schemes, as shown in [16,3]. The scheme reads (3.9) ψ τ 6 = e τ C 1 e a 1 τ T e τ C 2 e a 2 τ T e τ C 3 e a 3 τ T e τ C 4 e a 3 τ T e τ C 3 e a 2 τ T e τ C 2 e a 1 τ T e τ C 1 and its coefficients are collected in Table 3. Notice that all b i coefficients are positive and only one a i is negative. All methods from Table 2 will be tested and compared in Section 4.
In the one-dimensional case, d = 1, we have in addition [[T, U ], U ] = 2mU , so that the operators C i in scheme (3.9) are simply C i = (b i + 2c i mτ 2 )U . But in this case one can do even better since this feature leads to additional simplifications also at higher orders in τ . Specifically, a straightforward calculation shows that  Table 1 (d = 1). Although this reduction only manifests at orders higher than six, we can incorporate the flows of W 5,1 and W 7,1 into the composition, namely we can replace the e b i τ U in (3.4) by e τ D i , where In this way it is possible to reduce the number of exponentials in the composition and obtain more efficient integrators tailored for this special situation.
In the particular case of a 6th-order symmetric scheme it turns out that the d i and e i coefficients can be used to vanish some of the conditions at order seven, and thus reduce the overall error. The composition with and coefficients collected in Table 3 (d = 1) turns out to be particularly efficient, as shown in [1]. Here the parameters c i , d i and e i have been chosen to satisfy 4 out of 8 conditions at order 7.
The methods we have considered here are left-right symmetric compositions whose first flow corresponds to the functional U. It is clear, however, that similar compositions but now with the first flow corresponding to T can be considered. In that case, the schemes read (3.12) ψ τ p = e a 1 τ T e τ C 1 · · · e τ C 1 e a 1 τ T .  This corresponds to a different class of methods, in general with a different behavior and efficiency, since this problem possesses a very particular algebraic structure that is not preserved by interchanging T and U . We have also analyzed 6th-order schemes of this class, but we have not found better integrators than those collected in Tables 2 and 3 in our numerical experiments.

Numerical examples
This section is devoted to numerical illustrations of the previous splitting methods in the cases d = 1 and d = 2 (in (1.1)).
The splitting methods introduced above enable to reduce the numerical resolution of the Vlasov-Poisson problem (1.1) to one-dimensional linear transport problems of the form where z can denote the spatial direction x or the velocity direction v, a is a coefficient which does not depend on the advected direction z, and g denotes an initial condition given on a uniform grid of N points. Typically, a is the component of the vector v or of the electric field frozen at some grid point in the x-variable.
To deal with the one-dimensional advection equations, a semi-Lagrangian method is used (see [9,6,7]). Since the characteristics can be solved exactly in this case (a does not depend on z), the error produced by the scheme comes from the splitting (error in time) and from the interpolation step (error in x and v). Note that the interpolation is performed using highorder Lagrange polynomials (of order 17 in practice) so that the numerical solution of (4.1) writes where I is an interpolation operator (piecewise Lagrange interpolation in our case). We refer the reader to [2,9,6,7]) for more details. After each advection in the velocity direction (U part), the Poisson equation (1.2) is solved to update the electric potential φ. Note that in the case d = 2, the Hamiltonian splitting leads to 2-dimensional advections U and T . These subproblems are split again leading to one-dimensional advections; this does not introduce additional errors since it concerns linear advection for which this subsplitting is exact. The numerical resolution of the Poisson equations (1.2) and (2.12) to get φ and K is performed using a spectral method. Their derivatives are computed using high order finite differences.
We consider the following initial condition for (1.1) with d = 1 , v max = 8 and k = 0.5. In the case d = 2, the following initial condition for (1.1) is chosen We are interested in the total energy conservation H(f ) given by (1.3). Indeed, this quantity is theoretically preserved by (1.1) for all times, so it represents an interesting diagnostic. For a given time splitting, we introduce the discrete total energy H(f h )(t) (integrals in phase space are replaced by summations) where f h denotes the solution of the splitting scheme and we look at the following quantity (4.4) err H = max where t f > 0 is the final time of the simulation. We are also interested in the L 2 norm f h (t) L 2 of f h (which is also preserved with time) and we plot the quantity Different splitting will be studied regarding these quantities to compare their relative performances. First, we consider some splitting methods from the literature: the well-known 2nd-order Strang splitting (STRANG, σ = 3 flows per step size, even if we take σ = 2 in all the figures, since the last flow can be concatenated with the first flow at the next iteration), the so-called triple jump 4th-order composition [24] (3JUMP, σ = 7 flows) and the 6-th order splitting method proposed in [5] (06-23, σ = 23 flows). Then, the splitting methods introduced in this work are considered. When d = 1, the method of Table 3 (06-11, σ = 11 flows), and in the case d > 1 the schemes of Table 2: 06-9, 06-11 and 06-13, with σ = 9, 11 and 13 flows, respectively.
In the following figures, we choose a final time t f and the quantities (4.4) and Finally, we denote by N the number of points per direction we use to sample uniformly the phase space grid. In Figure 1, we first focus on the d = 1 case. We plot the quantity relative to the total energy err H defined in (4.4) for STRANG, 3JUMP, 06-23 and our 06-11 (see Table 3) using N = 256 points per direction and t f = 16. The expected orders of the different methods are recovered. However, even if 06-23 and 06-11 are both of 6th-order, 06-11 presents a better behavior since the total energy is better preserved up to two orders of magnitude than 06-23, with a comparable time CPU. Note that the 06-11 scheme has also been used with success in the one-dimensional context in [1]. In Figures 2, the same diagnostics as before is shown, but with two smaller values of N . For N = 64, we can also observe the plateau for small τ which reveals the phase space error. The level of this plateau can be decreased by increasing N .
On Figure 3, the time evolution of H(f h )(t) and f h (t) L 2 are displayed for the four splitting methods with different N at a almost constant CPU time: STRANG with τ = 1/8, 3JUMP with τ = 0.4, 06-23 with τ = 4/3 and 06-11 with τ = 0.64. It appears that the conservation of the total energy is very well preserved for 06-11. For the conservation of the L 2 norm, the benefit of high-order splitting is not so clear since all the curves are nearly superimposed. When N increases, we observe that the eruption time increases; the eruption time corresponds to the time at which the finest scale length of f reaches the phase space grid size. After this time, the error rapidly blows up (see [23]). Finally, on Figure 4, we display the whole phase space distribution function at time t f = 16 obtained with 06-11. We can observe the fine structures (filaments) which are typically developed in this nonlinear Landau test case. Next, we focus on the d = 2 case. In Figure 5 we plot the quantity relative to the total energy err H defined in (4.4) as a function of σ/τ for N = 64 and t f = 60. First, on the left part of Figure 5, methods of the literature are displayed: STRANG (order 2), 3JUMP (order 4), and 06-23 (order 6). This diagnostic enables us to recover the expected order of the different methods. Then, on the right part of Figure 5 we focus on our new methods 06-9, 06-11 and 06-13 (see Table 2). All the methods in this figure are of order 6, so that it enables to study the influence of the number of flows σ (the reference 06-23 method is also displayed) on the total energy conservation. Even if all the methods are of order 6, they have not the same precision. Indeed, adding some flows in the splitting method enables in our context to generate a more efficient method. Two explanations can be made: first, the coefficients can be chosen smaller and few of them are negative and second, the error term can be optimized (see Section 3). Finally, the 06-13 method appears to be the best method, reaching an error of about 10 −10 with τ ≈ 0.2. In Figure 6 we plot the quantity err L 2 defined in (4.5) as a function of σ/τ for the different methods in the case d = 2. This diagnostic enables to quantify the dissipation (small details of the solution are eliminated) of the numerical methods (see [7]). As mentioned in the case d = 1, the benefit of high-order time integrator is not very clear. We can also remark that the influence of the number of flows is not very significative (see right part of Figure 6);  indeed, when τ is small we observe that all the methods converges towards a constant (which is 0.04 for 3JUMP, 0.03 for 06-23 and a close value of 0.03 for the other methods). Here, 3JUMP and 06-9 are the worst methods; this may be linked with the fact that they contain large negative coefficients and then they present important zigzag (see [13]).

Convergence estimates
This last section is devoted to the rigorous mathematical analysis of Vlasov-Poisson equations and their approximation by splitting methods of the form (1.5) satisfying the order conditions to ensure (1.6). For a given multi-index p = (p 1 , . . . , p d ) ∈ N d , we denote by ∂ p x the multiderivative ∂ p 1 As functional framework, we will consider the spaces H r ν equipped with the norms where, ∂ p x and ∂ q v denote the usual multi-derivative in the x and v variables. In such spaces -already considered in [8] -and for r and ν large enough, the Vlasov-Poisson equation is well-posed and satisfies stability estimates  ensuring the convergence of stable and consistent numerical methods, see Theorem 5.4 and Lemma 5.5 below for precise estimates.
Before giving a complete proof of these results, we will state some useful estimates. In the following, we will denote by L q x and L q for q = 2 and q = ∞ the standard L q spaces on T d and T d × R d respectively. Similarly, for r ≥ 0, H r x and H r denote the standard Sobolev spaces on T d and T d × R d respectively.
Lemma 5.1. Let ν > d/2. Then we have for p = (p 1 , . . . , p d ) ∈ N d Proof. For a given function g(x, v), we have , as soon as ν > d/2. Applying this formula, we first see with the definition (5.1) that We then obtain (5.2) by applying this formula to . After using the Gagliardo-Nirenberg inequality (easily obtained in Fourier) for ν > d/2, we then deduce the second equation (5.3) from (5.2).
We will first give a meaning to the expansion (3.2). To do this we will use the two following inequalities: for r ≥ 1 and ν > d/2, Indeed, the first is clear from the definition of the H r ν norm. To prove the second, we use the fact that for given (m, p, q) ∈ (N d ) 3 satisfying |p|+|q| ≤ r and |m| ≤ ν as in the definition of (5.1), we have where W r+1,∞ x denotes the standard Sobolev space in the x variable, controlling (r + 1) derivatives in L ∞ . Now, we conclude by using for all γ > d/2 (see (5.3) above) with d = 1, 2, 3 and γ < min(ν, 2). Using the Hamiltonian formalism ad Of course the same lemma holds for the exact flows ϕ t T and ϕ t U . Using these expansions and the identification between splitting method (1.5) and methods based on composition of exponentials (3.4) which is done using the relation (3.2), we can make precise the notion of order that we consider in this paper. The algebraic conditions analyzed in Section 3 and the previous estimates yields order p methods in the following sense: Definition 5.3. The method (1.5) is said to be of order p, if there exist ν 0 , ν p ≥ 0 and r 0 , r p ≥ 0 such that for all ν > ν 0 , r > r 0 and all bounded set B of H r+rp ν+νp , there exist τ 0 and C > 0 such that for all f ∈ B and all τ ≤ τ 0 , We will show in the next sections that the condition (5.6) implies the scheme is of order p in the sense that it approximates the solution in H r ν over a finite time interval [0, T ] with a precision O(τ p ), provided the initial data is in H such that for all t ∈ [0, T ], we have the estimate Moreover, for two initial conditions f 0 and g 0 satisfying the previous hypothesis, we have Equations (5.8) and (5.9) show that the flow is locally bounded in H r+2ν+1 ν and locally Lipschitz in H r ν for r large enough (ν being essentially d/2). Before proving the theorem, we will show a stability lemma that will be useful both for the local existence of solutions and the analysis of splitting methods.
Lemma 5.5. Let α, β ∈ R + , ν > d/2 and r > 3ν be given. There exists a constant C > 0 such that the following holds: Assume that g(t) ∈ H r ν and f (t) in H r ν are continuous function of the time, and let h(t) be a solution of the equation Then we have Now we consider the operators D = D m,p,q = v m ∂ p x ∂ q v for given muti-indices (m, p, q) ∈ N 3d such that |m| ≤ ν and |p| + |q| ≤ r. It is then clear that the second term in the right-hand side can be bounded by 2 h H r ν g H r ν , and we are led to prove that The operator L κ can be split into a linear combination of operators of the form . . , d. We compute that for any smooth function h where i is the multi-index with coefficients δ ij the Kronecker symbol, for j = 1, . . . , d (we make the convention that D m,p,q = 0 when p or q contains that for all f and g satisfying f H r+rp ν+νp and we deduce (5.17) from the estimate (5.3) for φ(f ). For the L 2 estimates, we have and we conclude in a similar way by using the fact that for α > d/2, for some constant C depending only on α and d. We then deduce the result by applying the estimate (5.2) on φ(f ).
With these estimates in hand, it is easy to show that an existence result like Theorem 5.4 holds for the flow ϕ t [[T ,U ],U ] with similar estimates. Convergence results for splitting methods (3.9) can then be easily shown as in the proof of Theorem 5.6.

Conclusion
In this work, new time splitting schemes are proposed for the Vlasov-Poisson system. They are based on the decomposition of the Hamiltonian H between the kinetic T and electric U part. In the one-dimensional case, the relation [[T , U], U] = 2m U enables to design very efficient (with optimized number of flows) high-order splitting using the modified potential approach. This can be generalized to arbitrary dimension, the price to pay being to compute the flow associated to the commutator [[T , U], U] which only depends on the spatial variables; in this case also, new high-order splitting are proposed which turns out to be very efficient compared to the existing splitting of the literature. Finally, a convergence result of such splitting methods applied to the Vlasov-Poisson system is obtained.