High-order Hamiltonian splitting for Vlasov–Poisson equations

High-order Hamiltonian splitting for Vlasov–Poisson equations

Fernando Casas IMAC and Departament de Matemàtiques, Universitat Jaume I, 12071-Castellón, Spain. Fernando.Casas@uji.es Nicolas Crouseilles INRIA-Rennes Bretagne Atlantique, IPSO Project and IRMAR (Université de Rennes I) nicolas.crouseilles@inria.fr Erwan Faou INRIA-Rennes Bretagne Atlantique, IPSO Project and IRMAR (Université de Rennes I) Erwan.Faou@inria.fr  and  Michel Mehrenberger IRMA (Université de Strasbourg and CNRS) and INRIA-Nancy-Grand Est, TONUS Project mehrenbe@math.unistra.fr

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.

Key words and phrases:
Vlasov–Poisson equations, Runge–Kutta–Nyström methods, Hamiltonian splitting, high-order splitting
1991 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 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 high-order 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 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


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


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 semi-Lagrangian 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


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


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 ,


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


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 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 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



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


where is given by (1.2). Due to the relation , 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 and , we set


where the Fréchet functionals are evaluated in . Note that the skew-symmetry 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


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


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.

For any smooth function , the functionals and , satisfy the following relation:




is a constant of motion of (1.1), and

where is defined in (1.2). In dimension , we have , and in any dimension , the relation


holds for all functions .

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


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


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


Denoting by , the components of the electric vector fields , we get in the case


and in the case , .

Proof. Let us calculate

where we used (2.9). Hence we have

We deduce, that


Let us consider its Laplacian of

where , , denote the components of the electric vector fields . Using that , with independent of , we obtain


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 .

In dimension , we have from (2.13)

as expected.  

With the previous notations and Proposition 2.2, the equation associated with the Hamiltonian is given by


Hence the flow associated with is explicit and given by


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 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 equation (1.1) is equivalent to


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


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 RKN-type 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 high-order 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


We are mainly interested in symmetric compositions, that is, integrators such that , , so that


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 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 as the formal exponential of only one operator




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 6th-order 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 th-order 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
Table 1. Numbers of independent order conditions to achieve order required by symmetric splitting methods when (), and when , with constant ().

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 :


In that case the order conditions to achieve order 6 are explicitly


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


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.

of the form (3.7)
with ,
of the form (3.7)
with ,
of the form (3.7)
with ,
Table 2. Coefficients for symmetric schemes of order 6 for the Vlasov–Poisson equation in the general case () for and .

In the one-dimensional 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