Highorder Hamiltonian splitting for Vlasov–Poisson equations
Abstract.
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, highorder methods can also be constructed using explicit computations of commutators. Numerical results are performed and show the benefit of using highorder splitting schemes in that context. Complete and selfcontained proofs of convergence results and rigorous error estimates are also given.
Key words and phrases:
Vlasov–Poisson equations, Runge–Kutta–Nyström methods, Hamiltonian splitting, highorder splitting1991 Mathematics Subject Classification:
1. 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 socalled 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 highorder methods is classical when one speaks about space or velocity discretization. However, for the simulation of Vlasov–Poisson systems, the use of highorder 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 highorder 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 highorder methods, taking into account some specific commutator relations.
We consider the following equation Vlasov–Poisson equation
(1.1) 
where depends on time and the phase space variables , , and where for vectors and , we set and . Here, denotes the dimensional torus 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 corresponds to the spatial variable whereas is the velocity variable.
The electric potential solves the Poisson equation
(1.2) 
where is the Laplace operator in the variable acting of functions with zero average. The electric field depending only on is defined as . The energy associated with equations (1.1)(1.2) is
(1.3)  
The time discretization methods proposed in this paper are based on this decomposition of the energy. Indeed, the solution of the equations associated with and can be solved exactly (up to a phase space discretization, for example by interpolation in the framework of semiLagrangian methods). We denote by and the flows associated with and 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 , we have to solve the equation
(1.4) 
for which we verify that the solution is given by
where is the value of the electric field at time . Indeed, is constant along the solution of (1.4). Based on these explicit formulae, we will first consider numerical integrators of the form
(1.5) 
where and are real coefficients, and such that the numerical solution after time coincides with the exact solution up to terms of order , i.e., for a given smooth Êfunction ,
(1.6) 
where 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 of the numerical method.
As composition of exact flows of Hamiltonians and , 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 norms of ). 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 and in the decomposition (1.3) satisfy the following formal relation
(1.7) 
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 and 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 , the Vlasov–Poisson system even satisfies the stronger property
where is the total mass of 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 can be computed up to space discretization. We will retain these ideas to derive new highorder 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 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 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.
2. Hamiltonian structure
2.1. Poisson brackets
We define the microcanonical bracket of two (sufficiently smooth) functions by the formula
With this notation, we can write the Vlasov–Poisson equation as
(2.1) 
where
is the microcanonical Hamiltonian associated with . Recall that for a given functional , its Fréchet derivative is the distribution evaluated at the point , being defined by the formula
for all smooth variation . In general, a Fréchet derivative is an operator acting on , hence a rigorous writing of the previous formula necessitates a loss of derivative in . We will discuss these issues in Section 5.
Considering the two functionals and defined in (1.3), their Fréchet derivatives read explicitly
(2.2) 
where is given by (1.2). Due to the relation , the Vlasov–Poisson equation can be written as
(2.3) 
which is a Hamiltonian equation for the Poisson structure associated with the following Poisson bracket: for two functionals and , we set
(2.4) 
where the Fréchet functionals are evaluated in . Note that the skewsymmetry is obtained using the relation
for three functions of and the fact that the integral in of a Poisson bracket of two functions always vanishes. Moreover, this bracket satisfies the Jacobi identity
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 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 . This is essentially a consequence of the fact that the nonlinear transport equation (2.1) involves divergence free vector fields. Let be a smooth function, and consider the functional
(2.5) 
Its Fréchet derivative is and using the definition (2.4), we can observe that for all Hamiltonian functionals , we have
owing to the fact that for all functions and . Hence the functionals (2.5) are invariant under any dynamics of the form (2.3) (see (3.1) below). They are called Casimir invariants of the Poisson structure. Typical examples are given by the norms of the solution .
2.2. Relations between and
Let us remark that, using (1.2) we have
and hence the potential energy can be written
where
is the kernel of the inverse of the Laplace operator in the variable.
The aim of this subsection is to prove the following result:
Proposition 2.1.
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 , where stands for a small perturbation satisfying . First, we have
(2.9) 
Hence, we get
We see that, using an integration by parts in , the third term can be written as
We deduce that
Using this relation, we calculate
Now we see that the term involving the function vanishes, as the integral of in is equal to . We can thus write after an integration by parts
In other words, we get
(2.10) 
But we have with (2.7) and (1.2)
and this yields (2.6). In dimension , we can further write that
and conclude that is identically equal to . In any dimension , as the Fréchet derivatives of and depend only on , we automatically obtain (2.8).
2.3. Flow of the Hamiltonian
As mentioned above, the Fréchet derivative of the Hamiltonian only depends on . 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 , we can calculate its Fréchet derivative.
Proposition 2.2.
For any smooth function and using the notations introduced above, we have
where satisfies
(2.11) 
Denoting by , the components of the electric vector fields , we get in the case
(2.12) 
and in the case , .
Proof. Let us calculate
where we used (2.9). Hence we have
We deduce, that
(2.13) 
Let us consider its Laplacian of
where , , denote the components of the electric vector fields . Using that , with independent of , we obtain
(2.14) 
In dimension , 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 .
With the previous notations and Proposition 2.2, the equation associated with the Hamiltonian is given by
(2.15) 
Hence the flow associated with is explicit and given by
(2.16) 
because, depending only on and integrals of in , it is constant in the evolution of the flux associated with .
In dimension and , (and then ) can be easily computed in Fourier space by solving (2.11): in particular, the computational cost of a term of the form is essentially the same for (standard splitting) as for .
3. Derivation of highorder 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 equation (1.1) is equivalent to
(3.1) 
where here are functionals acting on some function space. We will not discuss here the mathematical validity of such an equivalence, but taking for instance as a norm of a function space (see Section 5), we can prove that the solution admits a formal expansion of the form
(3.2) 
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
with the operators and satisfying the relation in dimension , where is a constant, or the RKNtype relation . To derive splitting methods in dimension , we will also consider numerical schemes containing blocks based on the exact computation of the flow associated with the Hamiltonian . In this section we will concentrate on the derivation of highorder splitting methods of the form (1.5) satisfying these formal relations.
Scheme (1.5) is at least of order for the problem (1.1) if and only if the coefficients , satisfy the consistency condition
(3.3) 
We are mainly interested in symmetric compositions, that is, integrators such that , , so that
(3.4) 
In that case, they are of even order. In particular, a symmetric method verifying (3.3) is at least of order [13]. Notice that the number of flows in the splitting method (1.5) or (3.4) is , 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 and per step is precisely .
Restriction (1.6) imposes a set of constraints the coefficients , in the composition (3.4) have to satisfy. These are the socalled 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 as the formal exponential of only one operator
(3.5) 
where
(3.6) 
and are polynomials in the parameters , . Here we assume that the coefficients satisfy (3.3).
The integrator is of order if , and thus the order conditions are up to the order considered. For a symmetric scheme one has , so that only involves even powers of . In consequence, automatically in (3.6) and we have only to impose . The total number of order conditions can be determined by computing the dimension of the subspaces spanned by the nested commutators involving and for , see [17].
For the problem at hand identically, and this introduces additional simplifications due to the linear dependencies appearing at higher order terms in . The number of order conditions is thus correspondingly reduced. In Table 1 we have collected this number for symmetric methods of order (line ). Thus, a symmetric 6thorder scheme within this class requires solving order conditions (the two consistency conditions (3.3) plus conditions at order plus conditions at order ), so that the scheme (3.4) requires at least 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 thorder splitting method involving exponentials ( stages) was designed which has been shown to be very efficient for a number of problems, including Vlasov–Poisson systems [12].
Order  2  4  6  8 

2  4  8  18  
2  4  8  16 
We have shown in the subsection 2.3 that, besides the flow corresponding to , the flow associated to can also be explicitly computed in a similar way as . Moreover, since , both flows commute so that we can consider a composition (1.5) with the flow replaced by . Equivalently, in the composition (3.4) we replace by , where :
(3.7) 
In that case the order conditions to achieve order 6 are explicitly
(3.8)  
together with the consistency conditions (3.3). Here , , . The two equations in the first line of (3.8), together with (3.3), lead to a method of order four. With the inclusion of 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 is instead of for scheme (3.4). There are several other systems where the evaluation of the flow associated with is not substantially more expensive in terms of computational cost than the evaluation of , 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 , and exponentials. When 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 , 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 . In this case the two free parameters have been chosen to vanish the coefficient multiplying the commutator at order 7 and such that . This procedure usually leads to very efficient schemes, as shown in [16, 3]. The scheme reads
(3.9) 
and its coefficients are collected in Table 3. Notice that all coefficients are positive and only one is negative. All methods from Table 2 will be tested and compared in Section 4.









In the onedimensional case, , we have in addition , so that the operators in scheme (3.9) are simply . 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
(3.9) 