A Legendre-Fourier spectral method with exact conservation laws for the Vlasov-Poisson system

A Legendre-Fourier spectral method with exact conservation laws for the Vlasov-Poisson system

G. Manzini , G. L. Delzanno , J. Vencels , and S. Markidis T-5 Applied Mathematics and Plasma Physics Group, Los Alamos National Laboratory, Los Alamos, NM 87545, USA HPCViz Department, KTH Royal Institute of Technology, Stockholm, Sweden Istituto di Matematica Applicata e Tecnologie Informatiche, Consiglio Nazionale delle Ricerche (IMATI-CNR),
via Ferrata 1, I – 27100 Pavia, Italy,

We present the design and implementation of an -stable spectral method for the discretization of the Vlasov-Poisson model of a collisionless plasma in one space and velocity dimension. The velocity and space dependence of the Vlasov equation are resolved through a truncated spectral expansion based on Legendre and Fourier basis functions, respectively. The Poisson equation, which is coupled to the Vlasov equation, is also resolved through a Fourier expansion. The resulting system of ordinary differential equation is discretized by the implicit second-order accurate Crank-Nicolson time discretization. The non-linear dependence between the Vlasov and Poisson equations is iteratively solved at any time cycle by a Jacobian-Free Newton-Krylov method. In this work we analyze the structure of the main conservation laws of the resulting Legendre-Fourier model, e.g., mass, momentum, and energy, and prove that they are exactly satisfied in the semi-discrete and discrete setting. The -stability of the method is ensured by discretizing the boundary conditions of the distribution function at the boundaries of the velocity domain by a suitable penalty term. The impact of the penalty term on the conservation properties is investigated theoretically and numerically. An implementation of the penalty term that does not affect the conservation of mass, momentum and energy, is also proposed and studied. A collisional term is introduced in the discrete model to control the filamentation effect, but does not affect the conservation properties of the system. Numerical results on a set of standard test problems illustrate the performance of the method.

Vlasov-Poisson, Legendre-Fourier discretization, conservation laws stability

Collisionless magnetized plasmas are described by the kinetic (Vlasov-Maxwell) equations and are characterized by high dimensionality, anisotropy and a wide variety of spatial and temporal scales [15], thus requiring the use of sophisticated numerical techniques to capture accurately their rich non-linear behavior.

In general terms, there are three broad classes of methods devoted to the numerical solution of the kinetic equations. dimensional phase space via macro-particles that evolve according to Newton’s equations in the self-consistent electromagnetic field [3, 16]. PIC is the most widely used method in the plasma physics community because of its robustness and relative simplicity. The well known statistical noise associated with the macro-particles implies that PIC is really effective for problems where a low signal-to-noise ratio is acceptable. The Eulerian-Vlasov methods discretize the phase space with a six dimensional computational mesh [8, 30, 13]. As such they are immune to statistical noise but they require significant computational resources and this is perhaps why their application has been mostly limited to problems with reduced dimensionality. For reference, storing a field in double precision on a mesh with cells requires about terabytes of memory. A third class of methods, called transform methods, is spectral and is based on an expansion of the velocity part of the distribution function in basis functions (typically Fourier or Hermite), leading to a truncated set of moment equations for the expansion coefficients [2, 12, 18, 17, 29]. Similarly to Eulerian-Vlasov methods, transform methods might be resource-intensive if the convergence of the expansion series is slow.

In recent years there seems to be a renewed interest in Hermite-based spectral methods. Some reasons for this can be attributed to the advances in high performance computing and to the importance of simpler, reduced kinetic models in elucidating aspects of the complex dynamics of magnetized plasmas [27, 22]. Another reason is that (some form of) the Hermite basis can unify fluid (macroscopic) and kinetic (microscopic) behavior into one framework [5, 34, 11, 35]. Thus, it naturally enables the ’fluid/kinetic coupling’ that might be the (inevitable) solution to the multiscale problem of computational plasma physics and is a very active area of research [24, 10].

The Hermite basis is defined by the Hermite polynomials with a Maxwellian weight and is therefore closely linked to Maxwellian distribution functions. Two kinds of basis have been proposed in the literature (differing in regard to the details of the Maxwellian weight): symmetrically- and asymmetrically-weighted [17, 29]. The former features -stability but conservation laws for total mass, momentum and energy are achieved only in limited cases (i.e., they depend on the parity of the total number of Hermite modes, on the presence of a velocity shift in the Hermite basis, …). The latter features exact conservation laws in the discrete and the connection between the low-order moments and typical fluid moments, but -stability is not guaranteed [4, 29]. Earlier works pointed out that a proper choice of the velocity shift and the scaling of the Maxwellian weight (free parameters of the method) is important to improve the convergence properties of the series [17, 5]. Indeed, the optimization of the Hermite basis is a crucial aspect of the method, which however at this point does not yet have a definitive solution.

One could of course envision a different spectral approach which considers a full polynomial expansion without any weight or free parameter. While any connection with Maxwellians is lost, such expansion could be of interest in presence of strong non-Maxwellian behavior and eliminates the optimization problem. The Legendre polynomials are a natural candidate in this case, because of their orthogonality properties. They are normally applied in some preferred coordinate system (for instance spherical geometry) to expose quantities like angles that are defined on a bounded domain. Indeed, Legendre expansions are very popular in neutron transport [28] and some application in kinetic plasma physics can be found for electron transport described by the Boltzmann equation [31]. Surprisingly, however, we have not found any example in the context of collisionless kinetic theory and in particular for the Vlasov-Poisson system.

The main contribution of the present paper is the formulation, development and successful testing of a spectral method for the one dimensional Vlasov-Poisson model of a plasma based on a Legendre polynomial expansion of the velocity part of the plasma distribution function. The expansion is applied directly in the velocity domain, which is assumed to be finite. It is shown that the Legendre expansion features many of the properties of the asymmetrically-weighted Hermite expansion: the structure of the equations is similar, the low-order moments correspond to the typical moments of a fluid, and conservation laws for the total mass, momentum and energy (in weak form, as defined in Sec. 4) can be proven. It also features properties of the symmetrically-weighted Hermite expansion: -stability is also achieved by introducing a penalty on the boundary conditions in weak form. This strategy is inspired by the Simultaneous Approximation Strategy (SAT) technique  [33, 20, 21, 32, 6, 26, 7, 25].

The paper is organized as follows. In Sec. A Legendre-Fourier spectral method with exact conservation laws for the Vlasov-Poisson system the Vlasov-Poisson equations for a plasma are introduced together with the spectral discretization: the velocity part of the distribution function is expanded in Legendre polynomials while the spatial part is expressed in terms of a Fourier series. The time discretization is handled via a second-order accurate Crank-Nicolson scheme. In Sec. A Legendre-Fourier spectral method with exact conservation laws for the Vlasov-Poisson system the SAT technique is used to enforce the -stability of the numerical scheme. In Sec. A Legendre-Fourier spectral method with exact conservation laws for the Vlasov-Poisson system conservation laws for the total mass, momentum and energy are derived theoretically. Numerical experiments on standard benchmark tests (i.e., Landau damping, two-stream instabilities and ion acoustic wave) are performed in Sec. A Legendre-Fourier spectral method with exact conservation laws for the Vlasov-Poisson system, proving numerically the stability of the method and the validity of the conservation laws. Conclusions are drawn in Sec. A Legendre-Fourier spectral method with exact conservation laws for the Vlasov-Poisson system.


We consider the Vlasov-Poisson model for a collisionless plasma of electrons (labeled “”) and singly charged ions (“”) evolving under the action of the self-consistent electric field . The behavior of each particle species with mass and charge is described at any time in the phase space domain by the distribution function , which evolves according to the Vlasov equation:


We assume the physical space to be periodic in , so that no boundary condition for is necessary at and , and that suitable boundary conditions, e.g., , are provided for at the velocity boundaries and for any time and any spatial position . We also assume that an initial solution is given at the initial time .

Remark 2.1

If the initial solution has a compact support in the phase space domain , then has also a compact support at any time . Moreover, the size of the support may increase in time in a controlled way, cf. [36, 14]. In such a case, it holds that until the size of the support equals the size of the velocity domain. This condition can be used to determine the final time at which a plasma simulation based on this numerical model is valid.

In the Vlasov-Poisson system, the electric field is the solution of the Poisson equation:


where is the dielectric constant and


is the total charge density of the plasma. By taking the time derivative of the Poisson equation and using the continuity equation

where is the total current density defined as


we obtain the Ampere equation


where is a suitable constant factor.

Remark 2.2

The Ampere equation can be used with the Vlasov equation instead of the Poisson equation to obtain the Vlasov-Ampere formulation. In the continuum setting, the two formulations are equivalent in the one-dimensional electrostatic case without any external electric field as the one considered in this work.


Consider the infinite set of Legendre polynomials , which are recursively defined for by [1, Chapters 8, 22]:


and normalized as follows


We remap the Legendre polynomials onto the velocity range through the linear transformation . Let be the inverse mapping from to . The -th Legendre polynomial is given by , where the scaling factor in front of is chosen to satisfy the orthogonality relation:


The first derivative of the Legendre polynomials is given by

where is a switch that takes value if is even, and if is odd. Using the chain rule and adjusting the normalization factor, we obtain the first derivative of the translated and rescaled Legendre polynomials:


The recursion relations that are used to expand the Vlasov equation on the Legendre basis are reported in appendix A for completeness.

Consider the spectral decomposition of the distribution function on the basis of Legendre polynomials given by


in the Vlasov equation (1). The boundary conditions are not exactly satisfied since they are imposed in weak form and the polynomials are not zero at the velocity boundaries.

A possible way to circumvent this issue is to consider the modified basis functions given by for . From the properties of the Legendre polynomials, it readily follows that for each and the expansion of on this set of functions will automatically satisfied the homogeneous conditions at the boundary of the velocity range. Nonetheless, we verified numerically in the first stages of this work that this approach may yield an unstable method and the numerical instability cannot be fixed as there is no mechanism that allows us to control the growth of the absolute value of the Legendre coefficients . Another possible choice is to consider for even and for odd . Although we have not implemented this second basis, a common characteristic of these choices is the loss of orthogonality, which we suspect may influence negatively the stability properties of the method. The alternative approach that we consider hereafter is to integrate by parts the velocity term in the Vlasov equations. This strategy allows us to set the boundary conditions in weak form, and, then, to introduce a penalty term to enforce the stability of the method through the boundary conditions (see Section A Legendre-Fourier spectral method with exact conservation laws for the Vlasov-Poisson system). To this end, we substitute (10) into (1), we multiply the resulting equation by and integrate over . Then, we use the recursion formulas (81a)-(81c) and the orthogonality relation (7) and we obtain the following system of partial differential equations for the Legendre coefficients :


where conventionally ,




is the boundary term resulting from an integration by parts of the integral term that involves the velocity derivative. The derivation of the coefficients and can be found in appendix A. If the distribution has compact support in , the homogeneous boundary conditions at and are imposed in weak form by assuming that in (13) is zero. However, since this term plays a major role in establishing the conservation laws and ensuring the stability of the discretization method, we will consider it in all the further developments and in the analysis of the next sections.

We truncate the spectral expansion of after the first Legendre modes by assuming that for and we approximate the distribution function by the finite summation:


The evolution of each coefficient with is still given by (11). To ease the notation, we will drop the subindex in by tacitly assuming that all the quantities containing are indeed numerical approximations dependent on the first modes of the truncated series.

Let be the vector that contains all the coefficients for , i.e., , and the vector containing the values of the Legendre shape functions evaluated at . It holds that . System (11) can be rewritten in the non-conservative vector form:




and . Since is a constant matrix, it follows that


Therefore, system (11) also admits the conservative form:


where is a real and symmetric matrix with real eigenvalues and eigenvectors.


We expand each Legendre coefficient on the first functions of the Fourier basis (for ) as follows


where each coefficient is a complex function of time . The Fourier basis functions satisfy the orthogonality relation


Substituting (20) in (11) and using (21), we derive the system for the coefficients , which reads as:


for and , and where denotes the convolution integral and denotes the mode of the Fourier expansion of the argument inside the square brackets. Explicit formulas for these quantities are given below. We also recall that if and are two given real functions of and and the coefficients of their Fourier expansion on the basis functions , then the -th Fourier mode of the convolution product is given by .

The Poisson equation for the electric field is similarly transformed by using (20) and the Fourier expansion of the electric field


into (2) to obtain


For the equation above becomes

which, according to the hypothesis of neutrality of the plasma, expresses the fact that the total charge in the system is zero. This implies that we can set the -th Fourier mode of the electric field to zero, i.e.,

For convenience of notation, we introduce the vector that contains the -th Fourier coefficients for all the Legendre modes , i.e., . System (22) can be rewritten in the vector form:




Note the vector expressions:

and for the -th Legendre components:

Consider the current density of species given by (4). We apply the Legendre decomposition (10), the Fourier decomposition (20) and we use (81b) to obtain the Legendre-Fourier representation of the total current density:


Taking the derivative in time of (24) and using (22) with , we obtain the Fourier representation of Ampere’s equation:


For and using definition (29) we reformulate Ampere’s equation as




For , the Fourier decomposition of Ampere’s equation (5) gives the consistency condition , the zero-th Fourier mode of the total current density .


To control the filamentation effect, we modify system (25) by introducing the artificial collisional operator in the right-hand side [4]:


Consider the diagonal matrix whose -th diagonal entry is given by:


and where is an artificial diffusion coefficient whose value can be different from species to species. Then, the collisional term is given by . The effect of this operator is to damp the highest-modes of the Legendre expansion, thus reducing the filamentation and avoiding recurrence effects. This operator is designed to be zero for , in order not to have any influence on the conservation properties of the method.


Let be the time step, the time index, and each quantity superscripted by as taken at time , e.g., , , etc. We advance the Legendre-Fourier coefficients in time by the Crank-Nicolson time marching scheme [9]. Omitting the superscript “” in and to ease the notation, Vlasov equation (22) for each species and any Legendre-Fourier coefficient becomes:


Equation (35) provides an implicit and non-linear system for the Legendre-Fourier coefficients as each electric field mode for depends on the unknown coefficient that must be evaluated at the same time . In practice, we apply a Jacobian-free Newton-Krylov solver [19] to search for the minimizer of the residual given by (35).

Consider the difference of the Fourier representation of Poisson’s equation (24) at times and


By setting in (35), recalling that and noting that the collisional term does not give any contribution, we find that


Using (37) in (36) yields the discrete analog of Ampere’s equation that is consistent with the full Crank-Nicolson based discretization of the Vlasov-Poisson system:


where we have introduced the explicit symbol


to denote the boundary terms related to the behavior of all the distribution functions of the plasma species at the boundaries of the velocity domain. In Section A Legendre-Fourier spectral method with exact conservation laws for the Vlasov-Poisson system we make use of (38) and (39) to characterize the conservation of the total energy.


The distribution function solving the Vlasov equation satisfies the so-called -stability property for . To see this, just multiply equation (1) by and integrate over the phase space domain . Assuming that the velocity range is sufficiently large for having , a simple calculation shows that . This property is particularly useful for , which implies the stability of the method (sometimes called also “energy stability” in the literature). To derive a relation for the stability of the Legendre-Fourier method, we need the result stated by the following lemma. The proof of the lemma requires a few lengthy calculations and is reported, for the sake of completeness, in appendix C.

Lemma 3.1

Let be the vector containing the Legendre coefficients of the -th Fourier mode of the distribution function , the electric field and the matrix of coefficients defined in (17). Then, it holds that:


where denotes the zero-th Fourier mode of the argument inside the brackets, and denotes the conjugate transpose. All terms in (40) are real numbers.

The stability of the Legendre-Fourier method depends on the behavior of the distribution function at the boundaries and . This result is stated by the following theorem.

Theorem 3.1

The coefficients of the Legendre-Fourier decomposition have the property that:


Proof.  Multiply (33) from the left by , the conjugate transpose of , to obtain:

Add to this equation its conjugate transpose. Matrix is real and symmetric, is real and the spatial term cancels out from the equation. Summing over the Fourier index we end up with:


Since is a diagonal matrix with negative real entries and is a real quantity it holds that


The assertion of the theorem follows by applying the result of Lemma 3.1.     

If at the velocity boundaries (see Remark 2.1), then at any instant the time derivative in (42) is negative due to the collisional term and we have that

Note that in absence of the collisional term (take in (34)) the time derivative is exactly zero and is constant. We refer to this property as the stability because the orthogonality of the Legendre and Fourier basis functions implies that


(see appendix B), from which we immediately find the stability of the distribution function . However, and can be different than zero and in general they are non zero since the Legendre polynomials are globally defined on the whole domain and are non zero at the velocity boundaries. If the right-hand side of (41) becomes positive, the collisional term may be not enough to control the other term in the right-hand side of (41). Therefore, the method may become unstable and the time integration of is arrested.

According to [33] we can enforce the stability of the method by introducing the boundary conditions in weak form in the right-hand side of system (33) through the penalty coefficient . To this end, we modify system (33) as follows:


By suitably choosing the value of the penalty we minimize or set equal to zero the term in the right-hand side of (41) that may cause the numerical instability. This result is presented in the following theorem.

Theorem 3.2

The modified form (45) of the Legendre-Fourier method for solving the Vlasov-Poisson system is -stable for and any . The coefficients of the Legendre-Fourier decomposition have the property that:


Proof.  Repeating the proof of Theorem 3.1 yields:


Due to (40), the first term of the right-hand side of (47) is zero (when the coefficient that multiplies is non zero) by setting