Arbitrary-Order Time-Accurate Semi-Lagrangian Spectral Approximations of the Vlasov-Poisson System

Arbitrary-Order Time-Accurate Semi-Lagrangian Spectral Approximations of the Vlasov-Poisson System

L. Fatone , D. Funaro , and G. Manzini Dipartimento di Matematica, Università degli Studi di Camerino, Italy; e-mail: lorella.fatone@unicam.it Dipartimento di Scienze Fisiche, Informatiche e Matematiche, Università degli Studi di Modena
e Reggio Emilia, Italy; e-mail: daniele.funaro@unimore.it
Istituto di Matematica Applicata e Tecnologie Informatiche, Consiglio Nazionale
delle Ricerche, via Ferrata 1, 27100 Pavia,
Group T-5, Applied Mathematics and Plasma Physics, Theoretical Division,
Los Alamos National Laboratory, Los Alamos, NM, USA; e-mail: gmanzini@lanl.gov
Abstract

The Vlasov-Poisson system, modeling the evolution of non-collisional plasmas in the electrostatic limit, is approximated by a Semi-Lagrangian technique. Spectral methods of periodic type are implemented through a collocation approach. Groups of particles are represented by the Fourier Lagrangian basis and evolve, for a single timestep, along an high-order accurate representation of the local characteristic lines. The time-advancing technique is based on Taylor developments that can be, in principle, of any order of accuracy, or by coupling the phase space discretization with high-order accurate Backward Differentiation Formulas (BDF) as in the method-of-lines framework. At each timestep, particle displacements are reinterpolated and expressed in the original basis to guarantee the order of accuracy in all the variables at relatively low costs. Thus, these techniques combine excellent features of spectral approximations with high-order time integration. Series of numerical experiments are performed in order to assess the real performance. In particular, comparisons with standard benchmarks are examined.

\@xsect

The Vlasov-Poisson system of equations describes the dynamics of a collisionless plasma of charged particles (electrons and ions), where the only relevant interaction is driven by a self-consistent electrostatic field [10]. Although the Vlasov-Poisson system is one of the simplest models that can be considered in plasma physics, its numerical treatment is quite challenging to the numerical modelers. In fact, each plasma species is described by a distribution function that is defined on a high-dimensional phase space. Since the beginning of numerical plasma simulations in the ’60s, a number of methods have been proposed to the scientific community and thoroughly investigated. We can roughly regroup them in a few big families: Particle-in-Cell (PIC) methods, Transform methods, Eulerian and Semi-Lagrangian methods.

The PIC method is very popular in the plasma physics community, as it is the most widely used method because of its robustness and relative simplicity [8]. There, the evolution of a plasma is described by the motion of a finite number of macro-particles in the physical space. These macro-particles are tracked along the characteristics of the Vlasov equation and their mutual interaction is driven by a nonlinearly coupled electric field, which solves the Poisson equation. The right-hand side of the Poisson equation depends on the charges carried by the macro-particles. The convergence of the PIC method for the Vlasov-Poisson system was proved in [21, 49, 50]. The PIC method has been successfully used to simulate the behavior of collisionless laboratory and space plasmas and provides excellent results for the modeling of large scale phenomena in one, two or three space dimensions [8]. Also, implicit and energy preserving PIC formulations that are suitable to long time integration problems are available from the most recent literature [11, 17, 18, 34, 35, 38, 46]. Nonetheless, PIC codes suffer from intrinsic drawbacks. As proved in [21], achieving high numerical resolution in multidimensional plasma physics simulations may require a huge number of particles, thus making such simulations infeasible even with the most powerful supercomputers currently available. Since only a relatively limited number of particles can be considered in practical calculations, the method is used in a suboptimal way and tends to be intrinsically noisy. Although research has been carried out to reduce PIC noise [39], the method remains effective mainly for problems with a low noise-to-signal ratio, and where the physics is not driven by fine phase space structures.

Based on the seminal paper [30], an alternative approach, called the Transform method, was developed at the end of the ’60s, which uses a spectral decomposition of the distribution function and leads to a truncated set of moment equations for the expansion coefficients [2]. To this end, Hermite basis functions are used for unbounded domains, Legendre basis functions for bounded domains, and Fourier basis functions for periodic domains, see, e.g., [36, 40, 33, 48, 47]. These techniques can outperform PIC [13, 14] in Vlasov-Poisson simulations. Moreover, they can be extended in an almost straightforward way to multidimensional simulations of more complex models, like Vlasov-Maxwell [23]. Convergence of various formulations of these methods was shown in [28, 37]. Transform methods offer a few indisputable advantages. First of all, they may be extremely accurate since they are based on a spectral approximations of the differential operators. Furthermore, physically meaningful discrete invariants (such as total number of particles, momentum and total energy) can be built directly from the expansion coefficients [42, 32]. The existence of such discrete invariants implies better stability properties in long-time integration problems. However, despite their good properties their implementation may be computational demanding. As a matter of fact, they suffer of the “curse of dimensionality” (i.e., a bad scaling of the computational complexity with the number of dimensions), when multidimensional basis functions are built by tensor product of one-dimensional ones.

An alternative to PIC and Transform methods is offered by the class of Eulerian and Semi-Lagrangian methods, which discretize the Vlasov equation on a grid of the phase space. Common approaches for the implementation are: Finite Volume Methods [25, 5], Discontinuous Galerkin [3, 4, 31], finite difference methods based on ENO and WENO polynomial reconstructions [20], or propagation of the solution along the characteristics in an operator splitting framework [1, 16, 19, 27, 26, 44, 22]. Semi-Lagrangian methods were first developed for meteorological applications in the early ’90s [6, 7, 45]. The aim was to take advantage of both Lagrangian and Eulerian approaches. Indeed, these methods allow for a relatively accurate description of the phase space using a fixed mesh and propagating the values of the distribution function along the characteristics curves forward or backward in time. High-dimensionality is typically addressed by a splitting operator strategy in order to advance the solution in time. Such a splitting makes it possible to approximate a multi-dimensional time-dependent problem by a sequence of one-dimensional problems. For the one-dimensional Vlasov-Poisson system, the splitting reformulates the Vlasov equation in two advection subproblems that advance the distribution function in space and velocity independently. High-order approximations are described in [41].

In this paper, we propose Semi-Lagrangian methods that provide the spectral accuracy of Transform methods. In particular, space and velocity representations are discretized using a spectral collocation approach and the approximation of the distribution function is advanced in time by following backward the characteristic curves. Furthermore, we do not resort to any time splitting of the Vlasov equation and the desired order of accuracy in time, e.g., or even higher, is attained by using well calibrated representations of the characteristic curves. The major advantage of our approach is to combine, in a simple and natural way, spectral accuracy with on purpose time discretization techniques, in principle of any order of convergence. The formulation of the method is the same for any space and velocity dimension, provided we adopt a multi-index notation. Finally, although we do not address these topics in the present paper, we note that an efficient implementation is possible by resorting to standard libraries such as the Discrete Fast Fourier Transform (DFT)] [12]. Moreover, we remark that unsplit algorithms, like the ones that we propose in this work, are more suited to task parallelization on multicore processors in comparison to split algorithms in the standard Semi-Lagrangian approaches.

The paper is organized as follows. In Section Arbitrary-Order Time-Accurate Semi-Lagrangian Spectral Approximations of the Vlasov-Poisson System, we present the continuous model. In Section Arbitrary-Order Time-Accurate Semi-Lagrangian Spectral Approximations of the Vlasov-Poisson System, we introduce the spectral approximation in the phase space. In Section Arbitrary-Order Time-Accurate Semi-Lagrangian Spectral Approximations of the Vlasov-Poisson System, we present a Semi-Lagrangian scheme based on a first-order accurate approximation of the characteristic curves, making use of a suitable Taylor expansion. In Section Arbitrary-Order Time-Accurate Semi-Lagrangian Spectral Approximations of the Vlasov-Poisson System, we derive more refined time discretization schemes, built in the framework of the method-of-lines, applying second-order and third-order multi-step Backward Differentiation Formula (BDF). To show the flexibility of our approach, we also present a single-step second-order approximation in time. In Section Arbitrary-Order Time-Accurate Semi-Lagrangian Spectral Approximations of the Vlasov-Poisson System, we investigate the conservation properties of the method and we show that the number of particles is always an exact invariant of the method, regardless of the order of the time discretization. Within a spectral accurate error, this is also true for momenta. Concerning the total energy, this is conserved up to an approximation error that depends on the accuracy of the time discretization. In Section Arbitrary-Order Time-Accurate Semi-Lagrangian Spectral Approximations of the Vlasov-Poisson System, we show the predicted convergence rate in time by using a manufactured solution. Furthermore, we assess the performance of the method on standard benchmark problems as the two stream instability and the Landau damping. In Section Arbitrary-Order Time-Accurate Semi-Lagrangian Spectral Approximations of the Vlasov-Poisson System, we present our final remarks and conclusions.

\@xsect\@xsect

The distribution functions , , solving the Vlasov-Poisson system describe the statistical evolution of a collection of collisionless charged particles of distinct species, subject to mutual electrostatic interactions [10]. From a physical viewpoint, each represents the probability of finding particles of species in an element of volume , at time and point in the phase space , where , . The 3D-3V Vlasov equation for the -th species with mass and electric charge reads as:

(1)

where represents the electric field. The initial condition for is given by a function , so that

(2)

The coupling with the self-consistent electric field is taken into account through the divergence equation:

(3)

where is the charge density of species . In (3) is the dielectric vacuum permittivity and , is the total charge density. We refer the reader interested in the theoretical analysis of the Vlasov-Poisson model and the properties of its solutions to [9, 29, 24].

\@xsect

To ease the presentation of the numerical scheme, we consider the 1D-1V Vlasov-Poisson formulation for the electron-ion coupled system. Consistently, we restrict the domain to and . Since positive ions (protons) are much heavier than electrons, we may assume that they do not move, so that their density distribution function is constant over . Without altering the generality of the exposition, we can set , , . By dropping out the label , we only have one distribution function for the electron species, so that the corresponding Vlasov equation and initial condition read as:

(4)
(5)

where the coupled electric field verifies the equation:

(6)

We recall that is the electron charge density defined by:

(7)

We assume the constraints (charge conservation):

(8)

where measures the extension of . By taking

(9)

equation (6) can be transformed into the Poisson equation for the potential field :

(10)

As far as boundary constraints in and are concerned, we will assume a periodic boundary condition for the Poisson equation and either periodic or homogeneous Dirichlet boundary conditions for the Vlasov equation.

In the continuum setting, the total number of plasma particles is preserved. Hence, from a straightforward calculation and using (8) it follows that:

(11)

Moreover, the distribution function solving the Vlasov-Poisson system satisfies the so-called -stability property for :

(12)

which holds for any In particular, we will be concerned with . In this case, (12) implies the -stability of the method [29] (sometimes called also energy stability in the literature).

Finally, we consider the total energy of the system defined by:

(13)

where the first term represents the kinetic energy and the second one the potential energy. The Vlasov-Poisson model is characterized by the exact conservation of the energy, i.e.:

(14)

If the electric field is smooth enough, for a “sufficiently small” , the local system of characteristics associated with (4) is given by the phase space curves solving

(15)

with the condition that when . Under suitable regularity assumptions, there exists a unique solution of the Vlasov-Poisson problem (4), (5), (6) and (7), see [29], which can formally be expressed by propagating the initial condition (5) along the characteristic curves that solve (15). Therefore, for every we have that

(16)

By using a first-order approximation of the characteristic curves given by:

(17)

the Vlasov equation is satisfied up to an error that decays as , for tending to . To achieve a higher order of convergence, we need a more accurate approximation of the characteristic curves, such as, for example, the one given by setting:

(18)

By direct substitution in (4), the Vlasov equation is satisfied at every point up to the quadratic remainder for tending to . Of course, (18) can be replaced by other more accurate expansions leading to a high-order remainder term proportional to for some integer . Without exhibiting the explicit formulas, which look rather involved, we point out this property as a possible extension for further generalizations.

In view of the expression above, it is also convenient to write the time derivative of the electric field by arguing as follows. We evaluate the time derivative of in (7) and use the Vlasov-Poisson equation:

(19)

where we observe that the integral of is zero for a periodic function or in presence of homogeneous Dirichlet conditions. Translated in terms of , the above equation implies the Ampère equation, which reads as:

(20)

(after an integration with respect to ). Finally, in order to preserve the conditions in (8), we must set in (20).

\@xsect

We propose a Semi-Lagrangian method to find numerical approximations to the self-consistent solutions of the 1D-1V Vlasov-Poisson problem defined by equations (4), (5), (6) and (7). The extension to higher-dimensional problems, e.g., the 3D-3V case, is straightforward and is discussed at the end of this section. Instead, in the subsequent sections, we will analyze suitable time discretization techniques. In view of imposing periodic boundary conditions, we start by considering the domain:

(21)

A function defined in is requested to be periodic in both and . This means that for any integer we must have:

(22)

and

(23)

where, as usual, the zero-th order derivative of the function (i.e., when ) is the given function itself.

Given two positive integers and , we consider the equispaced points in :

(24)

Hereafter, if not otherwise indicated, we will always use the indices and running from to to label the grid points along the -direction, and and running from to to label the grid points along the -direction.

Then, we introduce the Fourier Lagrangian basis functions for the and variables with respect to the nodes (24), that is:

(25)
(26)

It is known that

(27)

where is the usual Kronecker symbol.

Furthermore, we define the discrete spaces:

(28)

In this way, any function that belongs to can be decomposed as:

(29)

where the coefficients of the decomposition are given by:

(30)

For what follows, it will be useful to have the expression of the derivatives of the basis functions. For instance, one has:

(31)

and

(32)

More generally, will denote the -th derivative of evaluated at point , which is given by:

(33)

Analogously we can define:

(34)

where , in (34) are obtained by replacing the nodes with the nodes in (31) and (32) and setting up the indices accordingly. As a special case we set: , . Moreover, it is easy to prove that there exists a constant , independent of , such that:

(35)

This estimate will be useful in the next section for studying the stability conditions in the time-marching schemes.

Furthermore, we remind that the following Gaussian quadrature formula:

(36)

which can be applied to any , is exact for every . For more details see [15, Section 2.1.2] and [43, Section 2.1.2].

In truth, given an integer , the derivative of order is trivially obtained by applying the first derivative matrix to the point-values of the -th derivative of a trigonometric polynomial. Such an operation can be performed by the fast Fourier transform algorithm, with an excellent cost reduction when the degree is relatively high and a power of 2, and very efficient implementations exist in freely available and commercial software libraries.

It is clear that, with little modifications, we can handle Lagrangian basis of nonperiodic type. Among these, the most representative ones are constructed on Legendre or Chebyshev algebraic polynomials, or Hermite functions (i.e., Hermite polynomials multiplied by a Gaussian function). In some preliminary tests, we observed that each one of these cases presents peculiar behavior in applications. A comparison between the different approaches would be too lengthy for the aims of the present paper. Therefore, we prefer to examine more deeply these extensions in a future analysis.

Now, consider the one-dimensional function . Given , by taking in formula (17), we define the new set of points where

(37)
(38)

where we recall that index is running through the range and index through the range . To evaluate a function at the new points through the coefficients in (30), we use the Taylor expansion. For a sufficiently smooth function , we have that

(39)

Applying (Arbitrary-Order Time-Accurate Semi-Lagrangian Spectral Approximations of the Vlasov-Poisson System) to , when , is defined in (37), we obtain:

(40)

Using (27), (33) and (34), we can rewrite (Arbitrary-Order Time-Accurate Semi-Lagrangian Spectral Approximations of the Vlasov-Poisson System) as:

(41)

Substituting (Arbitrary-Order Time-Accurate Semi-Lagrangian Spectral Approximations of the Vlasov-Poisson System) in (29), we obtain:

(42)

In compact form we can write:

(43)

where we set and . Finally, we truncate the summation with respect to at the integer to have a remainder term of order . The differentiation in the variables and can be computed exactly by multiplying the corresponding derivative matrices. Therefore, no approximation is introduced if we assume that the integer can range from 1 to infinity in (43).

\@xsect

The three-dimensional extension of (43) is straightforward by using the multi-index notation. To this end, we consider all indices in (43) as multi-indices of order three. More precisely, is the triplet of nonnegative integers and is the order of . The position vector is thus given by , and, a similar notation holds for the velocity position vector . A space vector subindexed by has to be interpreted as the grid point ; a velocity vector subindexed by has to be interpreted as the grid point . Consistently, we also have the double-subindexed vectors and . We use the standard notation for any given three-dimensional vector and multi-index , and we denote the partial derivatives of order of a generic function determined by the multi-index as:

A similar relation holds for the partial derivatives along . Finally, the three-dimensional basis functions are given by the tensor product of the one-dimensional basis functions:

Now, the three-dimensional version of equation (43) becomes:

(44)

where we set , , , ; the partial derivatives of the three-dimensional basis functions are given by

(45)

All considerations at the end of the previous section are still true here.

\@xsect

Given the time instants for any integer , we consider here the full approximation of the solution fields of the 1D-1V Vlasov-Poisson problem (4), (5), (6), (7):

(46)

where the function belongs to and the function belongs to . Taking into account (7), we define:

(47)

At any timestep , we evaluate in the following way:

(48)

where

(49)

In particular, at time , we use the initial condition for (see equation (5)) by setting

(50)

If we suppose that is given at step , we first define (take in (17)):

(51)

Since the solution of the Vlasov-Poisson system is expected to be constant along the characteristics, the most straightforward method is obtained by advancing the coefficients of as follows

(52)

where we used representation (48). This states that the value of , at the grid points and timestep , is assumed to be equal to the previous value at time , recovered by going backwards along the characteristics. Technically, in (51) we should use instead of , thus arriving at an implicit method. However, the distance between these two quantities is of the order of , so that the replacement has no practical effects on the accuracy of the first-order method. For higher order schemes, things must be treated more carefully.

Between each step and the successive one, we need to update the electric field. This can be done as suggested here below.

Let be fixed. Using the Gaussian quadrature formula (36) in (47) and (49) we write:

(53)

Indeed, it is possible to compute by using the Fourier series:

(54)

where the discrete Fourier coefficients and are determined, for , by the following formulas:

(55)

Actually, for strictly smaller than , the symbol “” can be replaced by the symbol “”.

Using equation (54) and equation (6) at , we conclude that: