High-order Hamiltonian splitting for the 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ö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.


Mathematics Subject Classification 37M15 Symplectic integrators · 65M25
Method of characteristics · 65M12 Stability and convergence of numerical methods · 65P10 Hamiltonian systems including symplectic integrators · 65Z05 Applications to physics · 82D10 Plasmas

Introduction
Frequently, the Vlasov equation is solved numerically with particles methods. Although they can reproduce realistic physical phenomena, they are well known to be noisy and converge slowly when more particles are considered in the simulation. To remedy this, so-called Eulerian methods (which use a grid of the phase space) have become increasingly popular in the last decades. Indeed, due to the increase in computer 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 [19,23] 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 high-order methods for ODEs (see [4,5,13,22]). The main goal of this work is to construct high-order splitting schemes for the nonlinear Vlasov-Poisson PDE system in 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 Sect. 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 Vlasov-Poisson equation in the following form 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 on functions with zero average. The electric field depending only on x and 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 the precise definition of Hamiltonian flows until Sect. 5). 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 chosen such that the numerical solution at 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 Sect. 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 the classical settings based on the Baker-Campbell-Haussdorff (BCH) formula and the Lie calculus, see for instance [4,13]. A rigorous justification will be given in Sect. 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 Sect. 2). This property reduces the number of order conditions on the 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 [4,13]. 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 [4,5]. The paper is organized as follows: • In Sect. 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 Sect. 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 Sect. 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, Sect. 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.

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 where 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 variations δ 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 Sect. 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 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 Sect. 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

Relations between T and U
Let us remark that, using (1.2) we have and hence the potential energy U can be written as 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: where is a constant of motion of (1.1), and 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 We can thus write after an integration by parts In other words, we get 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 [[T , U ], U ]
As mentioned above, the Fréchet derivative of the Hamiltonian [[T , U], U] only depends on x. Hence its exact flow can be calculated explicitly, making it possible to be included in the splitting methods blocks in any dimension. The situation is completely analogous to Hamiltonian systems in classical mechanics when the kinetic energy is quadratic in momenta, see e.g. [4,17,18]. From the expression (2.10) of the Poisson bracket [[T , U], U], we can calculate its Fréchet derivative.

Proposition 2 For any smooth function f and using the notations introduced above, we have
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 We deduce, that Now we consider the Laplacian applied to K where E x j , j = 1, · · · , d, denote the components of the electric vector fields E. Using (2.14) In dimension d = 2, we get 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.
In dimension d = 1, we have from (2.13) With the previous notations and Proposition 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 is essentially the same for γ = 0 (standard splitting) as for γ = 0.

Derivation of high-order methods
The splitting methods (1.5) are composition of exact flows of Hamiltonian equations of the form (2.1). To analyze their orders of approximation, we will use the algebraic structure of the Vlasov-Poisson equation. For a Hamiltonian equation of the form (2.3), let us define This notation is justified by the fact that the Eq. (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 Sect. 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 We are mainly interested in symmetric compositions, that is, integrators such that In that case, the resulting integrators 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 is done, we can express ψ τ p as the formal exponential of only one operator ψ τ p = e τ (T +U +R(τ )) , (3.5) where and p i j are polynomials in the parameters a i , b i . Here we assume that the coefficients satisfy (3.3).
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 τ . As a 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 the problem at hand [[[T, U ], U ], U ] = 0 identically, and this introduces additional simplifications due to the linear dependencies appearing at 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 to be of order 6. In fact, it is a common practice to consider more exponentials than strictly necessary and use the free parameters introduced 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 Sect.
In that case the order conditions to achieve order 6 are explicitly given by the following equations 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). Indeed, in this case, we have to solve five polynomial equations in five parameters. This has been done by applying Mathematica which produced 3 solutions as output, only one of them being real. 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 Sect. 4) corresponds to σ = 13. In this case the two free parameters have been chosen such that the coefficient multiplying the com-

mutator [T, [T, [T, [T, [T, [T, U ]]]
]]] at order 7 vanishes and such that b 1 = b 2 . This procedure usually leads to very efficient schemes, as shown in [3,16]. The scheme reads ψ τ 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 (3.9)   Table 2. 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 Sect. 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 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 ψ τ p = e a 1 τ T e τ C 1 · · · e τ C 1 e a 1 τ T . (3.12) 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 splitting methods introduced in the previous section in the cases d = 1 and d = 2 (in (1.1)). The splitting methods introduced above enable to reduce the numerical solution 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 in time.
To deal with the one-dimensional advection equations, a semi-Lagrangian method is used (see [6,7,9]). 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 high-order 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,6,7,9]) 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 solution 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 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 wellknown 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 (4.5) are plotted as a function of σ/τ , where σ is the number of flows of the considered method and τ is the time step used for the simulation. This choice ensures that all the diagnostics are obtained with a similar CPU cost. In the sequel, we consider 70 Case d = 1. We first focus on the d = 1 case. In Fig. 1, 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. Both methods require a similar amount of CPU time. Note that the 06-11 scheme has also been used with success in the one-dimensional context in [1]. In Fig. 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 Fig. 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 for a constant CPU time: STRANG with τ = 1/8, 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]).
On Fig. 4, we look at time evolution of err H for larger final time (t f = 64), considering N = 4096 (which ensures phase space convergence) and different time steps (τ are the same as before, but also τ/2 are considered). We observe the expected behavior for times smaller than 30. Then, for 06-11 and 06-23, err H increases, illustrating the fact that the fine structure is not well resolved. For STRANG and 3JUMP, err H still  remains at the same level. Indeed, since their time error is quite important, the lack of accuracy to capture fine structure does not affect it. As a conclusion, even for large times, the proposed splitting 06-11 presents the best behavior. Finally, on Fig. 5, 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.
Case d = 2. Next, we focus on the d = 2 case. In Fig. 6 we plot the quantity relative to the total energy err H defined in . This diagnostic enables us to recover the expected order of the different methods. Then, on the right part of Fig. 6 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 we are able to study the influence of the number of flows σ (the reference 06-23 method is also displayed) on err H . 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 Sect. 3). Finally, the 06-13 method appears to be the best method, reaching an error of about 10 −10 with τ ≈ 0.2. Concerning the L 2 norm conservation, the splitting preserves it exactly when x and v are supposed continuous. However, it is not true anymore when we are interested in the fully discretized case. We have remarked that the L 2 norm conservation mainly depends on the number of stages of the splitting method, but also on the L ∞ norm of the coefficients (see [13]).

Analysis of the Vlasov-Poisson equation
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 , . . . , As the 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 1 and Lemma 3 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 1 Let ν > d/2. Then we have for p
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 space) 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 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 Sect. 3 and the previous estimates yields order p methods in the following sense:

Definition 1
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 sets B of H r +r p ν+ν 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 r +r p ν+ν p .

Existence of solutions
The goal of this subsection it to prove the following result: 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 3
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 functions in time, and let h(t) be a solution of the equation (5.10) Then we have where [D, L κ ] = DL κ − L κ D is the commutator between the two operators. The first term in the previous equality can be written 2 Dh, L κ Dh L 2 = 1, L κ (Dh) 2 where L * κ is the L 2 adjoint of L κ , upon using the fact that Hamiltonian vector fields are divergence free. Hence we get 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 on the right-hand side can be bounded by 2 h H r ν g H r ν , and we are led to prove that | Dh, [D, L κ ]h L 2 | ≤ C(α + β f H r ν ) h and χ t

Splitting methods with iterated commutators
To end this section, we give arguments to show that Theorem 2 still holds for more general splitting methods defined by formulas (3.9) (3.11) using flows associated with iterated commutators. Let us first consider the one-dimensional case d = 1. In this case, the flow of [[T , U], U] = 2m U, as well as the flows of the high-order commutators W 5,1 , W 7,1 and W 7,2 (see (3.10)) are in fact flows of the Hamiltonian U scaled in time. Hence, all the previous convergence results extend straightforwardly to numerical schemes of the form (3.11).
In the higher dimensional cases d ≥ 2, the flow ϕ t [[T ,U ],U ] is given by the formulas (2.15), (2.16). We see that it has the same structure as the flow ϕ t U , but the potential φ( f ) is replaced by the potential K ( f ) given by the formula (2.14). This potential satisfies the following estimates (compare with Lemma 1)

Lemma 4 Let ν > d/2. Then we have for r
and Proof Using the definition (2.11) of K , we have We deduce that for r ≥ 2,  (3.9) can then be easily shown as in the proof of Theorem 2.

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.