Variational formulation of particle algorithms for kinetic plasma simulations

Variational formulation of particle algorithms for kinetic plasma simulations

E. G. Evstatiev111Present address: FAR-TECH, Inc., 10350 Science Center Drive, Bldg. 14, Suite 150, San Diego, CA 92121. B. A. Shadwick Department of Physics and Astronomy, University of Nebraska-Lincoln, NE 68588-0111
July 5, 2019

Common time-explicit numerical methods for kinetic simulations of plasmas in the low-collisions limit fall into two classes of algorithms: momentum conserving [also known as Particle-In-Cell (PIC)] and energy conserving. Each has certain drawbacks. The PIC algorithm does not conserve total energy, which may lead to spurious numerical heating (grid heating). Its overall accuracy is at most second due to the nature of the force interpolation between grid and particle position. Energy-conserving algorithms do not exhibit grid heating, but because their formulation uses potentials, computationally undesirable matrix inversions may be necessary. In addition, compared to PIC algorithms for the same accuracy, these algorithms have higher numerical noise due to the restricted choice of particle shapes. Here we formulate time-explicit, finite-size particle algorithms using particular reductions of the particle distribution function. These reductions are used in two variational principles, a Lagrangian-based and a Hamiltonian-based in conjunction with a non-canonical Poisson bracket. The Lagrangian formulations here generalize previous such formulations. The Hamiltonian formulation is presented here for the first time. Many drawbacks of the two classes of particle methods are mitigated. For example, restrictions on particle shapes are relaxed in energy conserving algorithms, which allows to decrease the numerical noise in these methods. The Hamiltonian formulation of particle algorithms is done in terms of fields instead of potentials, thus avoiding solving Poisson’s equation. An algorithm that conserves both energy and momentum is presented. Other features of the algorithms include a natural way to perform coordinate transformations, the use of various time integrating methods, and the ability to increase the overall accuracy beyond second order, including all generalizations. For simplicity, we restrict our discussion to one-dimensional, non-relativistic, unmagnetized, electrostatic plasmas.

Numerical, Plasma, Kinetic, Variational, Energy Conserving, Momentum Conserving, Particle-In-Cell
journal: YJCPH\biboptions


1 Introduction

Particle-based simulation methods have a long and successful history in plasma physics. The original proposal by Hockney Hockney (1966), based on a method developed for fluid simulations Harlow (1964), used computational macro-particles moving in a self-consistent, mean field. Fields were approximated on a spatial grid and interpolated to the particle position to determine the Lorentz force. The essential physics was successfully captured but the simulations suffered from high numerical noise as -functions were used to represent the macro-particles. It was later realized Dawson (1983), that -function particles used in this way can lead to numerical instability Dawson (1960). A significant improvement was achieved by allowing the macro-particles to have a finite spatial extent Langdon and Birdsall (1970); these improved schemes originated the class of methods now known as Particle-In-Cell (PIC) algorithms. The PIC algorithm was first used to model plasmas with negligible collisions but later Coulomb collisions and atomic physics were included with techniques based on the Monte Carlo method Vahedi and Surendra (1995). PIC simulations are now widely used to study far ranging plasma systems including modeling laser-plasma interactions Faure et al. (2004); Geddes et al. (2008); Yin et al. (2009); Mori et al. (2011); Vay et al. (2011), Z-pinches Welch et al. (2011), astrophysical and magnetized plasmas Daughton and Karimabadi (2005); Lapenta et al. (2006); Brackbill and Lapenta (2008), plasma discharges and low temperature plasma processing Vahedi and Surendra (1995); Nanbu (2000), and numerous other applications.

PIC methods have subsequently undergone a significant development. Research on improvements in reliability and stability lead to the discovery of non-physical, purely numerical artifacts. The PIC algorithm does not conserve total energy exactly (even in the absence temporal discretization), which leads to a surprising phenomenon known as “grid heating” Langdon (1970); Okuda (1972). (See Ref. Cormier-Michel et al. (2008) for a overview of grid heating.) Grid heating is attributed to a kinetic instability where sub-grid structures are aliased to low frequency modes (due to finite grid resolution). Typically this instability saturates once the plasma has heated to the point where the Debye length is on the order of the grid spacing Langdon (1970). More recently, it has been shown that the choice of particle shape (the spline used for current and charge deposition and force interpolation) can lead to unphysical effects, especially when considering threshold phenomena such as self-trapping in a laser-plasma accelerator Cormier-Michel et al. (2008). A further limitation of the PIC algorithm is that its overall accuracy is at most second order in both space and time. This is due to the interpolation (typically splines) of quantities between the continuous particle position and the spatial grid.

In an attempt to correct for these deficiencies of the standard PIC algorithm, energy-conserving particle algorithms were devised Lewis (1970). While strict energy conservation eliminated the grid heating instability, these algorithms had their own drawbacks. They did not seem to have the same flexibility with respect to a choice of particle shapes as PIC algorithms, which in turn affected the level of numerical noise; i.e., for the same numerical accuracy, PIC algorithms had lower noise levels. Energy conserving algorithms also may require mass-matrix inversions, which is avoided by the field-based formulations of PIC. As a result, energy-conserving algorithms did not become as popular as PIC algorithms.

One major difference between PIC and energy-conserving algorithms lies in the way each is formulated. In PIC algorithms Hockney and Eastwood (1988); Birdsall and Langdon (1991), relations between the discretized electric (vector) potential, electric (magnetic) field, charge deposition, and current deposition were obtained by discretizing the corresponding continuous relations and equations. In this process, critical terms of the order of the accuracy of discretization are dropped (as is justified in an asymptotic procedure), but which led to the loss of energy conservation (and possibly violation of other conservation laws). In comparison, energy conserving algorithms were derived from a variational principle Lewis (1970); Eastwood (1991), using the fact that the Vlasov equation could be obtained from Low’s Lagrangian Low (1958). A number of benefits in using variational principle was pointed out: basic properties of the original system were retained in the reduced system; a natural way of making coordinate transformations was provided; and use of high accuracy space and time solvers was possible.

The goal of this paper is to generalize previous variational formulations and to offer a new formulation of energy conserving algorithms based on the Hamiltonian and the non-canonical Poisson bracket proposed by Morrison Morrison (1980); Weinstein and Morrison (1981); Morrison (1982). Our approach is to use particular reductions of the distribution function to a finite collection of terms as well as particular reductions of the continuous fields to finite number of degrees-of-freedom, either in the Lagrangian or in the Hamiltonian and Poisson bracket. As a result of our general method, we show how to avoid many of the previous drawbacks and deficiencies of both PIC and the energy-conserving methods.

In addition to the energy-conserving property of all algorithms in this work, we: (i) show that particle shapes in energy conserving algorithms can be chosen with more freedom instead of being a delta-function in space. In fact, the shape of an extended space particle has very few physical constraints; the shape may be symmetric about its centroid or may exploit the spatial symmetry of a particular physical problem, e.g., systems with azimuthal symmetry or symmetries in the gyrokinetic approximation Lee and Qin (2003); Lin et al. (2005); (ii) relax the method of time integration of the equations of motion (most often leapfrog previously) allowing for higher than second order accuracy. The choice of a time integrator becomes limited only by numerical stability. In this paper we emphasize the spatial discretization and leave time continuous, which allows for convenient formulations of time-explicit schemes. We prove conservation laws with continuous time and only in the last step do we choose a particular (explicit) time advancing scheme; (iii) show how to reduce continuous quantities with either grid-based reduction (using finite differences) or truncated bases; previous authors have only used truncated bases, which for the case of finite elements may necessitate mass matrix inversions. By using grid-based reduction, mass matrices do not appear (for an example, see B), which may be computationally advantageous; (iv) derive formulations in terms of fields that are based on the Hamiltonian and a non-canonical Poisson bracket. Previously only potential-based energy-conserving algorithms have been derived from a variational principle Lewis (1970); Eastwood (1991). Using field-based formulations eliminates the need for solving Poisson’s equation; and (v) derive a particle algorithm that conserves both total energy and total momentum. The arguments leading to this algorithm demonstrate the usefulness of the variational approach and exploit the relations between conservation laws and symmetries of the Lagrangian. Previous such models were derived by assuming a delta function for the particle shape Evstatiev et al. (2003, 2005). Throughout, we restrict our discussion to the case of a one-dimensional, nonrelativistic, unmagnetized, electrostatic plasma. This is done purely to streamline the discussion, elucidating the central ideas. The extension of these concepts to the fully relativistic, electromagnetic case, while involving numerous technical details, is largely straightforward and will be the subject of a future publication.

It has been emphasized before Lewis (1970); Eastwood (1991) that variational formulations naturally lead to higher overall accuracy particle algorithms (i.e., both in space and time). Here we show that this remains true with all of the above relaxed conditions. We show that force interpolation, field integrators, and time integrators may be chosen to increase the accuracy of a particle method beyond second order. In the course of all derivations, we point out where a certain property is being relaxed or is being lost.

Time-implicit formulations of particle algorithms have significant attraction in simulating problems where long time evolution is necessary. Recently, authors have been successfully reformulated the PIC algorithm in terms of implicit time integration with the added benefit of energy conservation Chen et al. (2011); Markidis and Lapenta (2011). However, these formulations do not offer a general derivation and so is unclear how they can be extended. For example, these formulations use the Crank-Nicholson time integrator and charge-conserving particle shapes but it is an open question how to extend these methods to use more accurate time integrators.

While continuous-time formulations are the focus of this paper, we note that one may discretize the action principle in both space and time. In the context of particle methods for plasma simulations, this was first considered by Eastwood Eastwood (1991). While Eastwood’s method conserved energy exactly, it did so at the expense of an implicit time-advance. When the time-advance was altered to be fully explicit, exact energy conservation was lost. The notion of performing the temporal discretization in the action has been studied extensively by Marsden and co-workers (see for example Refs Wendlandt and Marsden (1997) and Marsden and West (2001)). There are a number of attractive features of this approach and it is a natural extension of the methods presented herein; this will be the subject of future work.

The paper is organized as follows. Section 2 is devoted to deriving algorithms based on a Lagrangian formulation of the Vlasov–Poisson system. This section presents the finite differencing formulation. In section 2.2 the particle models are derived from a Lagrangian formulation in terms of truncated bases. The important model which conserves both momentum and energy is derived there. In these derivations the fields are described in terms of electrostatic potential. Section 3 presents the Hamiltonian derivation of particle models using truncated basis and a reduction of the non-canonical Poisson bracket. The equations of motion are formulated in terms of electric field. Section 4 illustrates properties of the derived particle models with numerical examples. Conclusions are in section 5. A gives many examples of charge deposition rules, while B presents a hybrid cold fluid-kinetic particle model from a Lagrangian starting point. It demonstrates our general method with a different reduction of the particle distribution function. It also illustrates how the mass matrix and its inverse may be avoided by the use of grid-based reduction of the continuous quantities.

2 Lagrangian formulation

A plasma with negligible collisions is well described by a single-particle phase-space distribution function, , whose phase-space evolution is governed by the Vlasov equation Krall and Trivelpiece (1973)


where is the electric field, is the electric potential, and are the species mass and charge. For an initial phase-space distribution , the distribution at any later time is given by


where and are the macro-particle trajectories with initial conditions and : and . The particle trajectories correspond to characteristics of the Vlasov equation and (2) is simply the statement that the distribution function is constant on the characteristics. Vlasov dynamics can be obtained from the Lagrangian Low (1958); Galloway and Kim (1971); Ye and Morrison (1992)


where and are to be varied independently. (Here we consider a single-species plasma but the extension to multiple species is obvious.) In the usual way, demanding the action be stationary with respect to variations of the dynamical variables leads to the equations of motion. Variation with respect to particle positions gives


while variation with respect to the potential gives


We have used Galloway and Kim (1971)


in (5) and we have assumed either periodic boundary conditions or an infinite system to allow surface terms to be dropped. Note that (6) is a statement of particle number conservation and is equivalent to Gardner’s restacking theorem Gardner (1963).

The basic idea of Lagrangian macro-particle methods lies in the representation of the full distribution function as a sum of moving spatial volumes, , called macro-particles:


The choice of a delta function in velocity space is not essential but avoids the necessity to track stretching phase space volumes, which is why we adhere to it. In Eq. (7), are constant weights and the function is the fixed spatial extent of the computational particle (hereafter we use the terms particle, computational particle, and macro-particle interchangeably unless otherwise specified) and is normalized as


An additional simplification is made by assuming that all particles have the same shape. We note that the representation (7) is general and independent of whether both electric and magnetic fields are present in the system, i.e., it is valid for the general Vlasov–Maxwell system. For clarity of the presentation, in this paper we consider only electrostatic, non-relativistic models. Electromagnetic and relativistic models can be derived similarly to this presentation and will be presented in a future publication. We view (7) as a particular reduction of the particle distribution function. B gives another example of such a reduction.

Substituting our form of the distribution function, (7), into the Lagrangian and again using (6), we obtain a reduced Lagrangian




Although we have replaced a continuum of particles with labels and by macro-particles, we still have an infinite degree-of-freedom system due to the presence of the continuous field . The equations of motion are obtained from (9) by considering variations of the particle position and of the potential. For the particles, the usual Euler–Lagrange equation




Since the potential is a field, the Euler–Lagrange equation for the potential is


where denotes a functional derivative. Then


Note that the factor appearing in (14) is the physical charge to mass ratio of the plasma species. It is not necessary to make the ad-hoc assumption that the macro-particle have the same charge to mass ratio as the plasma species, this is a consequence of the phase-space decomposition (7). Furthermore, the second form of the force in (14) may clearly be interpreted as the electric field averaged over the particle shape.

The substitution of (7) into (3) is equivalent to a choice of a trial function for that depends on a number of parameters, which in our case are the particle positions and velocities. The values of these parameters are obtained by solving the equations resulting from the variation (13). Other choices of trial functions may lead to models without particles at all Lewis et al. (1987); Shadwick et al. (2010).

A significant advantage of the variational formulation is the connection between symmetries and conservation laws as embodied in Noether’s theorem José and Saletan (1998). Our introduction of macro-particles through (7) neither results in explicit time-dependence in the Lagrangian nor breaks translational invariance of the Lagrangian, thus we should expect the equations of motion (14) and (16) to exactly conserve both energy and momentum. The total energy of the system is the sum of macro-particle kinetic energy and field energy,


Using the equations of motion, it is straightforward to see that is an invariant:


The total momentum of the system is simply


since the electrostatic field carries no momentum (the Poynting vector is zero in the electrostatic approximation). Now


where we have used (16).

At this point in our reduction, we have a finite number of macro particles representing the plasma but a continuous field for the potential. Here one must provide some approximate or exact solution of (16), which is then used to integrate the macro-particle equations of motion. One possibility would be to use methods based on evaluating the Green’s function, constructing as the superposition of the potentials due to each macro-particle. Even though the macro-particles interact via the mean field, computation of using Green’s functions scales as and is thus limited to relatively small systems. A more computationally advantageous alternative is to introduce a discrete representation for the potential. There are two general approaches: using a spatial grid, approximating the potential by its value at the grid point, or using a truncated set of (local or global) basis functions and representing the potential by its projection onto the basis.

The interaction term in the Lagrangian, (11), provides both the force in (14) as well as the charge density in (16). Of course, this will continue to be the case when the continuous potential is replaced by a discrete approximation. It will be necessary to approximate consistently with the choice of the discrete potential but this single approximation ultimately yields both the force term in the equation of motion and the charge density in a discrete analogue of Poisson’s equation. Thus we are guaranteed that these terms are consistently approximated.

2.1 Discretization using a spatial grid

We assume a fixed spatial grid with and grid spacing with being the numerical approximation to . We must now approximate two terms in (9), and . The interaction term requires knowledge of between the grid-points; so some manner of interpolation is required. Finite elements Becker et al. (1981) offer a consistent way to perform such interpolations to any accuracy. Let , be finite-element basis of some order. We interpolate between the grid points by


and thus (11) becomes




is the effective (projected) shape of the particle. Note the expression for can be computed analytically since the function is known. If are constructed from Lagrange polynomials, then = 1 and


This property means that the total charge deposited on the grid and at any instant of time is constant.

It remains to approximate in terms of . This can be approached in two ways which give roughly equivalent results. We can use (21) to write the integral in (12) as




we have


Alternatively, after integrating by parts in , we can approximate the integral as


While this appears to simply be using the trapezoidal rule to evaluate the integrand, with either periodic boundary conditions or an infinite domain, this approximation has spectral accuracy, that is, all modes supported by the grid are integrated exactly. We complete the approximation by choosing a finite-difference representation for the second derivative. Regardless of details of the finite-difference approximation, it can always be expressed as


for some integer . Thus we have


which has the same form as (27). Notice that while is always symmetric [cf. (26)], this need not be true for .

We now arrive at the finite degree-of-freedom Lagrangian


where we take either or . The dynamical equations are obtained by demanding that the action be stationary with respect to variations in both and . Taking these variations yields




The reason we choose to integrate by parts in (30) is now clear: by dong so, we are able to directly specify the difference method for the second derivative that appears in Poisson’s equation, (33).

Discretizing in (17) in the same manner as in the Lagrangian, we have


Using the equations of motion, we find


Introducing a spatial grid does not affect energy conservation. This is expected since the spatial discretization of does not introduce explicit time-dependence into the Lagrangian. An immediate advantage of the variational approach is that models derived in this way are automatically free of grid heating. A side effect of introducing a spatial grid is that it breaks the translation invariance of and consequently total momentum is no longer exactly conserved; see section 2.2 for a more complete discussion.

We conclude this section by providing a concrete example of this procedure to derive a model that is second-order accurate in . For simplicity, we consider the case of a charge-neutral electron plasma with an immobile ionic background and a spatially periodic domain. Since this system has no ion dynamics, we can forgo summing over species and simply introduce the ion density into the Lagrangian




with being the given ion density. Linear finite elements yield second order accuracy interpolation Becker et al. (1981) (see Figure 1):


To determine we need to specify . Regardless of the choice of , the accuracy of the interpolation will be second order due to our basis choice. The choice of affects the quality of our approximation through the extent to which (7) is a good ansatz but has no influence on the formal order of the model. Arguably the simplest choice for is a top-hat cell-wide function:


While we have chosen to be exactly one grid-cell wide, this is by no means essential. The choice of support of is completely independent of the grid spacing; the particular choice in (39), allows us to make connection with the usual PIC particle shapes (see Figure 9 and Table 1).

Figure 1: Linear finite element basis functions. The basis functions are identified with the grid-point at which it takes on the value ; e.g., is the tent function with support .

We now use (23) to determine the grid charge deposition. With this , for any , there are only three values of for which . Take to be the grid point nearest and let . Clearly . It is straightforward to evaluate (23) to obtain


This is equivalent to the charge deposition obtained from the usual PIC quadratic particle shape Hockney and Eastwood (1988). It is possible to recover all of the usual smooth particle shapes. For example, taking


we obtain


which is equivalent to the usual quartic charge deposition rule. We take up the matter of particle shapes in some detail in A. In particular, we demonstrate a cubic spanning only three grid-points.

All that remains is to approximate (12) to second-order accuracy. Evaluating (26) for our linear basis, is straightforward. For any , we can see that has a non-zero overlap with only for and thus we have


Alternatively, we can use finite difference approximations and evaluate (12) using (30). Since we are considering a periodic domain, it is reasonable to choose a central difference approximation


which gives


Linear finite-elements give the same approximation to (12) as taking second-order central differences. This is essentially a coincidence; higher-order finite element bases (quadratic, cubic, etc.) do not yield expressions for that can be equated to conventional differencing schemes. Nothing, other than the symmetry of the problem, forces us to choose central differences; any second-order approximation would suffice. Note, using an expression for that is accurate beyond second-order will not increase the overall spatial order of the method unless a corresponding more accurate interpolation scheme is used to evaluate (23).

The macro-particle equation of motion, which is not alerted by the specifics of the spatial discretization, remains as in (32). For a uniform ion background, , we have


For completeness we restate Poisson’s equation including the ionic background and our approximation for


2.2 Discretization using a truncated basis and the question of momentum conservation

In the previous section the interaction and field parts of the Lagrangian (9) were reduced from an infinite to a finite degree of freedom quantities by discretizing on a grid. We noted that the introduction of a spatial grid breaks the translational invariance of , which leads to loss of momentum conservation in the reduced system. In this section we consider a reduction using a truncated global basis and investigate the question of momentum conservation in this case. We show that replacing the continuous potential by a finite collection of projections onto a truncated basis can result in a discrete system that retains translation invariance. (Of course, if the basis is not truncated, which is not useful from a computational perspective, then we would expect translation invariance to be maintained for any complete basis.)

Let , be the first elements of an orthonormal basis. We approximate the potential as




and is the dual to satisfying


for . The accuracy of this approximation will depend on the number of elements kept and the convergence properties of the basis. For a truncated system to possess translation invariance requires that the basis have certain properties. These properties are made evident by shifting the origin by an amount while leaving the physical system unchanged and requiring this transformation to be a symmetry of José and Saletan (1998). The kinetic term, , is obviously translation invariant under such a shift since is time independent. Consider . Let be the new coordinate, with . The particle coordinates and potential relative to are


The symmetry condition is then


Introducing our basis expansion into and suppressing prime symbols, we have




Combining these expressions and expanding to lowest order in , we have


Our symmetry condition requires that the term multiplying in (55) vanish. For the symmetry to exist independent of the particle shape and the details of the potential, this term must vanish for each . We are led to the condition


On a finite domain with periodic boundary conditions this condition is satisfied by the discrete Fourier basis; we are not aware of any other discrete basis that fulfills (56) on either a finite or infinite domain. Using the discrete Fourier basis, it is straightforward to show that is also translation invariant.

We now specialize our discussion to the case of a truncated Fourier basis. Let


where , and is the domain size. With this basis, the interaction term becomes






and we have used the relation . We also need to evaluate (12):


where, since is real, . Finally we arrive at the discrete form of the Lagrangian


To obtain the equations of motion, we require the action to be stationary with respect to variations of and (since and are not independent, we need only consider variations of ). The equation of motion are




Using (59) and (57) it is easy to show


allowing us to write the equation of motion as


The spatial charge density associated with a single particle is


and the corresponding total charge is