On the velocity space discretization for the VlasovPoisson system: comparison between Hermite spectral and ParticleinCell methods.
Part 1: semiimplicit scheme
Abstract
We discuss a spectral method for the numerical solution of the VlasovPoisson system where the velocity space is decomposed by means of an Hermite basis. We describe a semiimplicit time discretization that extends the range of numerical stability relative to an explicit scheme. We also introduce and discuss the effects of an artificial collisional operator, which is necessary to take care of the velocity space filamentation problem, unavoidable in collisionless plasmas. The computational efficiency and the costeffectiveness of this method are compared to a ParticleinCell (PIC) method in the case of a twodimensional phase space. The following examples are discussed: Langmuir wave, Landau damping, ionacoustic wave, twostream instability, and plasma echo. The Hermite spectral method can achieve solutions that are several orders of magnitude more accurate (at a fraction of the cost) with respect to the PIC method.
keywords:
1 Introduction
The Vlasov equation represents the cornerstone for the kinetic modeling of collisionless plasmas. It describes the time evolution of the distribution function of a population of charged particles that respond selfconsistently to the effect of self and externally induced electromagnetic fields. The numerical solution of the Vlasov equation for collisionless plasmas represents a very active area of research. The main challenge resides in the fact that the distribution function lives in a six dimensional phase space. Unarguably, the most widely used technique to solve the Vlasov equation is the ParticleInCell (PIC) method (Birdsall and Langdon, 2005; Hockney and Eastwood, 2010). The main idea, originally developed in the 50’s, is to sample the distribution function in velocity space by means of a discrete number of superparticles. The electromagnetic field is represented on a grid in the computational domain, and the superparticles move through the grid according to the Lorentz force that is calculated by interpolation from the grid to the particles position. The particles interactions are therefore mediated by the grid, and in this way the number of operations per time step is reduced from (if the total force on a particle is calculated by summing the pair interaction with every other particle) to (with the number of superparticles). A comprehensive review of the PIC methodology can be found in (Verboncoeur, 2005).
The primary shortcoming of PIC simulations is related to the numerical noise: even starting from an equilibrium configuration, the discrete nature of the particles generates instantaneous (i.e. within the first time step) unphysical perturbations that produce a ’noise ground’ level in the fields, below which any physical signal is lost. The most obvious way to reduce the noise in PIC is by increasing the number of superparticles, i.e. refining the discretization in velocity space. The main problem is that while the simulation time scales roughly linearly with , the noise ground level decreases only as , as typical of Monte Carlo methods, implying that problems that require a high signaltonoise ratio require significant computational resources. Some examples that illustrate the impact of the particle noise on the efficiency of PIC codes are the recent studies in space plasmas on the acceleration of particles in the solar wind and the coupling between large scale turbulence and small scale (i.e. kinetic) dissipative effects (Valentini et al., 2008; Matteini et al., 2012; Haynes et al., 2013). PIC codes have demonstrated very poor performances for these problems. In twodimensional (2D) simulations, Camporeale and Burgess (2011) have shown that at least particles per cell (equivalent to in 3D) are needed to achieve a sufficiently large signaltonoise ratio. For comparison, one of the biggest and most advanced PIC simulations in the world (of a magnetic reconnection problem) was done with 250 particles per cell and required the use of the highest performance, petascale computer available at that time (Daughton et al., 2011; Bowers et al., 2009). We note that methods based on f or remapping of the distribution function have been developed in the PIC community to reduce particle noise. This has led to some very interesting work (see for instance (Chehab et al., 2005; Wang et al., 2011; Deng and Fu, 2014)), but has not yet resulted in a commonly accepted denoising technique adopted by the PIC community.
The nonoptimal scaling of the noise level with the number of particles in the PIC method has long been recognized in the computational plasma physics community (Byers, 1970; Denavit and Walsh, 1981; Sydora, 1999; Nevins et al., 2005), and it has perhaps been the main motivation for developing alternative ways of solving the Vlasov equation numerically. Nevertheless, despite its intrinsic noise, the PIC method still represents the preferred approach in the community, probably due to its simplicity of coding, the relative ease to parallelize it, and the recent impressive advances in available computing power. However, one should always keep in mind the resourceintensive nature of PIC and the related cost. As a figure of merit on the cost of large scale PIC simulations, a single run using the global gyrokinetic PIC code developed within the SciDAC Gyrokinetic Particle Simulation Center (GPSC) (Batchelor et al., 2007) for investigating the longtime evolution of turbulent transport in nuclear fusion devices was estimated, in 2008, to take approximately 25 million CPU hours, on 100,000+ cores (Tang, 2008). The monetary cost of such a simulation can be estimated at about $1,000,000 (Walker, 2009). In view of the upcoming exascale era and the complexity of multiscale plasma physics, it is therefore legitimate to question the computational efficiency of the PIC algorithm and wonder whether it is possible to devise algorithms that use the available computational resources more efficiently.
An alternative to the PIC method is represented by solving the Vlasov equation directly in phasespace (namely with a sixth dimensional computational grid in space and velocity),
by means of a socalled Eulerian Vlasov code.
This has been done, usually for problems with reduced dimensionality, with a variety of techniques such as highorder finite differences (Whealton et al., 1986; Guo and Qiu, 2012), finite elements (Besse and Sonnendrücker, 2003; Heath et al., 2012),
or finite volumes (Banks and Hittinger, 2010; Duclous et al., 2012).
The reader is referred to (Shoucri, 2008) for a recent review of different methods.
Notably, successful algorithms include the semiLagrangian schemes (Shoucri and Knorr, 1974; Cheng and Knorr, 1976; Sonnendrücker et al., 1999; Besse and Sonnendrücker, 2003; Umeda et al., 2006; Carrillo and Vecil, 2007; Crouseilles and Sonnendrücker, 2007; Imadera et al., 2009; Crouseilles et al., 2010; Qiu and Christlieb, 2010; Qiu and Shu, 2011; Rossmanith and Seal, 2011)
and positive flux conservative methods (Filbet et al., 2001).
Yet another class of techniques is represented by spectral methods where the velocity space is projected onto a complete orthogonal basis,
and therefore the phase space is effectively not gridded in the velocity coordinates (Engelmann et al., 1963; Armstrong et al., 1970; Denavit and Kruer, 1971; Klimas, 1983; Eliasson, 2003; Le Bourdiec et al., 2006).
Clearly, all these methods have advantages and shortcomings. Historically, the PIC method has been favored for largescale, multidimensional problems which are still
not doable with Vlasov solvers (mainly due to memory limitations) (Mangeney et al., 2002). On the other hand, Vlasov codes allow the study of fine scale details of the distribution functions that are typically inaccessible
in PIC codes, due to the above mentioned noise problem.
Vlasov codes, however, can suffer from serious numerical problems that can lead to a lack of discrete conservation properties, the occurrence of unphysical oscillations, and
the generation of nonpositive values in the distribution functions (Banks and Hittinger, 2010).
In this paper, we focus on a spectral method where the distribution function is projected onto an Hermite basis in velocity space.
This gives rise to a set of nonlinear, time and spacedependent PDEs for the coefficients of the expansion.
The use of Hermite functions to discretize the velocity space in the Vlasov equation dates back to the works of Grad (1949); Engelmann et al. (1963); Grant and Feix (1967), and Armstrong et al. (1970).
More recently this approach has been investigated by Holloway (1996); Schumer and Holloway (1998); Le Bourdiec et al. (2006) and, in the context of linearized equations, by Camporeale et al. (2006) and Siminos et al. (2011).
The expansion of the distribution function using an Hermite basis in velocity space can be appealing for several reasons. First, the Hermite functions are a complete orthogonal basis with respect to a Gaussian weight function. As such, they are the optimal basis to represent Maxwellianlike distributions. In fact, an exact Maxwellian in velocity can be represented by only one term of the Hermite basis. Second, as already noted by Grad (1949), the low order expansion coefficients are directly related to the low order moments of the distribution functions (i.e. density, mean velocity, energy, heat flux, etc.) that describe the plasma macroscopically and are usually the physical quantities of interest. As a consequence, the use of the Hermite basis allows a smooth transition from a fluid to a kinetic description by simply increasing the number of coefficients retained. This is an important and crucial feature of this method in approaching multiscale problems and in assessing the importance of kinetic effects, which is not available in PIC or Vlasov methods that are forced to treat the full distribution function.
This article describes a semiimplicit discretization in time of the truncated set of PDEs for the expansion coefficients.
A fullyimplicit time discretization is discussed in the companion paper (Part 2) (Camporeale et al., 2013).
The spatial dependence of the coefficients is represented by means of a discrete Fourier transform.
This paper builds largely upon the work of Schumer and Holloway (1998).
In (Schumer and Holloway, 1998) two different Hermite bases (so called ’asymmetrically’ and ’symmetrically’ weighted) and their properties were discussed, and a qualitative comparison between the FourierHermite (FH)
method and the PIC method was presented for simulations of Landau damping and bumpontail instability. Here we expand that work in three ways.
First, we propose a semiimplicit scheme that is numerically more stable than the explicit scheme presented in (Schumer and Holloway, 1998).
Second, we study the effects of an artificial collision operator, with the intent of controlling the filamentation process and suppressing/mitigating numerical instabilities.
Third, we quantitatively compare the FH method against a conventional explicit PIC code.
In order to assess which method must be preferred (and for which conditions), the key information is represented by the CPU time needed to obtain a solution with a certain accuracy
(this metric was not considered in Schumer and Holloway (1998)).
Hence, our comparisons are presented in terms of computational efficiency and efficacy. The former is essentially represented by the CPU time required to achieve a certain accuracy in the solution,
while the latter is a measure
of the costeffectiveness of the algorithm, i.e. how much an increased resolution and consequently a longer CPU time is paid off in terms of better accuracy (Lapenta and Chacón, 2006).
Although it is not surprising that a spectral method performs better than the superparticle discretization used in PIC, we are not aware of any other study that
presents such quantitative comparison in terms of computational efficiency and efficacy.
For the examples presented in this paper (involving mostly nearequilibrium Maxwellian plasmas in one dimension), the difference in performance between the FH method and the PIC is impressive, with the former
producing results that are several orders of magnitude more accurate for equal CPU time or, conversely, results with the same level of accuracy in a fraction of CPU time. Nevertheless, FH is still
in its ’infancy’ and the development of optimization techniques that can accelerate the convergence of the Hermite basis will be critical for multidimensional applications, particularly for plasmas
away from equilibrium.
The paper is organized as follows. Section 2 presents the mathematical foundation of the FH method, along with its conservation properties. Section 3 discusses the time discretization of the equations. Section 4 is devoted to a brief discussion of the filamentation problem in Vlasov codes, and how that translates to PIC codes. It also introduces the collisional term that we have implemented and used in some of the simulations. Section 5 presents the comparison between PIC and FH methods for five classical test cases: Langmuir wave, Landau damping, ionacoustic wave, twostream instability, and plasma echo. Finally, our conclusions are presented in Section 6.
2 Numerical Method: FourierHermite basis
We study the VlasovPoisson system in the 2dimensional phase space , where denotes position and velocity. The phase space is assumed to be periodic in . In order to describe the method we specialize to the case of a plasma consisting of an electron and a singly charged ion species. The quantities are normalized as follows. Time is normalized on the electron plasma frequency , velocities on the electron thermal velocity , lengths on the electron Debye length , the electric field on ( is the Boltzmann’s constant, is the electron temperature, is the electron mass and is the elementary charge), and densities on a reference density . The Vlasov equation for the species reads:
(1) 
Here is the particle distribution function, and are the charge and mass of the particles of species ( for electrons and ions, respectively) and is the electric field. The electric field is selfconsistently calculated through the Poisson equation:
(2) 
The discretization in velocity employs the asymmetricallyweighted Hermite basis:
(3)  
(4) 
where is the th Hermite polynomial. The distribution function is defined as:
(5) 
with , where and are two constant parameters for each species, and is the number of Hermite modes (equal for both species).
The following closure is used for both species: for .
The spatial dependence of the expansion coefficients will be treated later by means of a Fourier decomposition.
We note that the completeness of the Hermite basis does not depend on the choice of the free parameters and .
It is well known, however, that a proper choice of the rescaling coefficient can substantially accelerate the convergence of the series (Schumer and Holloway, 1998; Tang, 1993; Camporeale et al., 2006).
In this work, we allow even more
flexibility by introducing the free parameter , which is a shift in the Hermite function argument. The usefulness of such a parameter will be discussed in Section 5.4.
Formulas for calculating the optimal values for and were presented in Camporeale et al. (2006).
The Hermite basis has the following properties ( is the Kronecker delta):
(6)  
(7)  
(8) 
One can multiply Eq. (1) by , integrate in and, by using such properties, obtain:
(9) 
Similarly, the Poisson equation (2) becomes:
(10) 
One can notice the wellknown property that the electric field is related uniquely to the 0th order expansion coefficients , i.e., the density. By now treating and by means of a standard discrete Fourier decomposition (with the domain length, and modes):
(11)  
(12) 
and using the orthogonality of the Fourier basis, one can finally derive an infinite set of ordinary differential equations for the coefficients :
(14) 
where is the index associated with the Hermite (Fourier) basis. The expression for the electric field reads:
(15)  
(16) 
We note that the Fourier representation of Poisson equation (10) leaves the constant undefined, while imposing the constraint
(which physically means that the plasma is neutral).
However, the absence of an externally imposed electric field and the periodicity of the domain dictates that .
The solution of Equation (14) is the main objective of this paper.
An important point to realize is that the only nonlinearity of this equation is in the last term, which couples the electric field with the distribution function via a convolution.
Also, in the Hermite (velocity) space, the th coefficient is coupled only to , , to itself (if ), and to .
2.1 Conservation properties
As noted in Schumer and Holloway (1998), the decomposition in the asymmetric FourierHermite basis allows conservation of particle number and momentum exactly, irrespective of the chosen
time discretization scheme. On the other hand, exact energy conservation is not possible with the operator splitting scheme investigated in this paper.
It is however possible to recover exact energy conservation by employing a fullyimplicit time discretization, as discussed in the
companion paper (Part 2) (Camporeale et al., 2013).
Following (Schumer and Holloway, 1998), the total number of particles for each species is defined as
(17) 
By inspection of Eq. (14) one can see that , from which the conservation of mass follows.
The single species linear momentum is defined as:
(18) 
In order to prove the conservation of momentum, one can calculate the time derivative of the coefficients and :
(19)  
(20) 
One can use Eq. (15) to substitute in Eq. (20) and obtain:
(21) 
By taking into account that is the complex conjugate of (since the electric field in physical space is a real quantity),
it follows that and therefore the last term in parenthesis in Eq.(21) is zero. Hence =0.
3 Time discretization
In the literature it has been suggested that the VlasovPoisson system can be efficiently solved by means of an operator splitting procedure (Cheng and Knorr, 1976), i.e., by decoupling the advection and the acceleration terms, which are usually referred to as ’Xshift’ and ’Vshift’ operators, respectively. Therefore, we define the and the operators acting on as:
By following the secondorder splitting scheme introduced by Strang (1968), the Vlasov equation is solved by updating the distribution function from time to time
in the following three steps (subscript is omitted):

Solve , evolving into by an half time step . Calculate the electric field from ;

Solve , evolving into by a full time step , using the electric field calculated in 1);

Solve , evolving into by an half time step .
Reference (Schumer and Holloway, 1998) has employed this scheme by using an explicit fourthorder RungeKutta (RK4) integration for both the X and V shifts.
They justify this choice by arguing that it offers a reasonable trade between accuracy and stability. However, we have verified that such a scheme has a very limited stability region,
even with relatively small time step.
In this paper, we improve on the stability of the explicit scheme by introducing a semiimplicit scheme. The V shift is still treated with RK4, but the Xshift is treated implicitly.
We have tested both a first order implicit Euler scheme (BDF1) (Henrici, 1962) and a secondorder CranckNicolson (Crank and Nicolson, 1947)
scheme for the Xshift.
Note that the Strang splitting ensures the second order accuracy in time of the overall scheme.
In summary, the discretized form of the Xshift that advances the solution
for half time step is the following (where the superscript indicates the time step):
Xshift (BDF1)
(22) 
Xshift (CN)
(23)  
(24) 
Of course, both the BDF1 and the CN schemes can be written in matrix form. For each Fourier component , we have:
(25) 
where is the columnvector . The matrix is defined as
(26)  
(27)  
(28) 
with for CN and for BDF1. For the BDF1 scheme the matrix is the identity matrix, and for the CN scheme is defined as:
(29)  
(30)  
(31) 
Note that the matrices and commute in both cases (BDF1 and CN). Therefore the two consecutive Xshifts in steps 3) and 1): can be computed exactly as a single step:
(33) 
that is
(34) 
We solve Eq. (34) by performing an LU decomposition of the matrix .
The Vshift is a standard fourth order RungeKutta integrator, reported here for completeness.
Vshift (RK4)
Note that the Vshift calculated via RK4 requires the evaluation of four convolutions for each Hermite mode .
Convolution is an expensive operation, which has a cost if performed explicitly. A cheaper alternative is to
exploit the property that the Fast Fourier Transform (FFT) of the convolution of two vectors is equal to the product of the FFT of the two vectors.
In this way, the cost is reduced to (see, e.g. (Bracewell, 1999)).
In practice, this consists in transforming the electric field and the coefficients from Fourier to physical space, and transforming their product back from physical to Fourier space at each time step.
In Figure 1, we present an empirical proof of the enhanced stability of the semiimplicit scheme (Xshift treated with BDF1 or CN) with respect to the fully explicit scheme
(Xshift treated with RK4). This example is for the study of Landau damping, discussed in more detail in the next Sections. The blue and red lines show the evolution of the electric field
using respectively a BDF1 and a CN scheme, with time step . The BDF1 shows its mildly dissipative character by having a stronger damping than CN.
The black circles are the solution for with CN. The simulation is still stable, although not very accurate.
On the other hand, when using RK4 even with smaller time steps ( for black line and for green line)
the simulation is unstable after only a few electron plasma periods (in all the Figures the time is in units of ).
For this particular case, the semiimplicit scheme presented here is at least a factor of 10 more stable.
4 Filamentation and artificial collisionality
The development of increasingly smaller phase space structures in a collisionless plasma is very well known in plasma physics and
typically referred to as the filamentation process.
A classical and wellstudied example where filamentation occurs is linear Landau damping, i.e., the damping in time of an initial electric field perturbation due
to waveparticle resonances (Landau, 1946).
In this case filamentation is due to the presence of the factor ( is the wavevector of the perturbation)
in the perturbed distribution function that creates oscillations with smaller
and smaller wavelengths in velocity space as time evolves.
Clearly, any discretization of the velocity space is associated with a minimum wavelength that can be resolved, and, therefore, any
numerical simulation of the filamentation process with fixed resolution is bound to fail after a certain time.
In the case of the FH expansion method, it is known that the time at which the simulation is unable to capture
filamentation due to lack of resolution in velocity space scales approximately as the square root of the number of Hermite modes, (Grant and Feix, 1967; Canosa et al., 1974; Schumer and Holloway, 1998). This time is known
as the recurrence time: the effect of the lack of resolution is to reconstruct a distribution function similar to the initial one, from which the Landau damping starts again in a recurrent way.
Figure 2 shows the recurrence effect for the linear Landau damping studied with the FH code (this case will be studied in more detail in Section 5.2).
Blue, black, and red lines show the amplitude of (the
Fourier component of the electric field perturbed at time ) for , respectively. One can appreciate that the time at which the simulation manifests the
recurrence phenomenon approximately doubles when is multiplied by a factor of 4, as obtained in Schumer and Holloway (1998).
Several fixes have been suggested in the literature in order to overcome the filamentation process in Vlasov codes. They involve
some form of filtering or smoothing of the high order moments of the distribution function (see, e.g. (Klimas, 1987, 1994)) or, equivalently, the introduction of a weaklycollisional operator (Grant and Feix, 1967).
In this work, we have tested the effect of a collisional operator. This is a purely numerical artifact, that does not represent physical collisions,
and must be used in a convergence sense, namely by tuning it such that the physical results of interest are unaffected, while the small filamentation structures are damped.
We have opted for a collision operator that acts
on the coefficient as:
(35) 
where is the collisional rate applied to the last Hermite coefficients . The collision operator of Eq. (35) can be constructed as a nonlinear combination of terms involving the LenardBernstein collision operator (Lenard and Bernstein, 1958):
(36) 
which, in the FourierHermite space, reads:
(37) 
An important point is that the LenardBernstein operator transformed in the
FourierHermite space acts on all the coefficients , including , while our operator is defined in such a way that it does not change the first three Hermite modes.
The reason for choosing the form in Eq. (35) is that it conserves charge and momentum (and energy if a fullyimplicit time discretization is employed), by
leaving unchanged.
The fact that this operator does not have a physical interpretation is not important, since our goal is to study collisionless plasmas.
Also, the regime of validity
of the simulation will be reduced to times for which the collisional rate is not dominant.
Of course, other forms of the collision operator might be employed: for instance (Parker and Dellar, 2012) proposes an iterative
version of Eq. (36). The crucial feature, however, is that the damping rate applied to the coefficient must increase (in absolute value) with increasing Hermite index .
The effect of the artificial collision operator is to damp the highest modes of the Hermite expansion. The convergence of the series implies that
for large enough . However, high modes can grow due to the filamentation process or even just due to
numerical errors. As we will show in Section 5.4, this can lead to numerical instabilities if the growth of the large modes is not suppressed artificially.
The effect of collisionality on the Landau damping study is presented in Figure 3. Here and three values of have been used: .
One can notice that the correct value of Landau damping (i.e. the one obtained before recurrence when )
is recovered for a long time, when and, therefore, in this case the use of collisionality is crucial to overcome
the recurrence due to filamentation. The small box shows a zoomin for time .
The plot in Figure 3 must be interpreted in light of the important result presented in Ng et al. (1999), and rediscussed in (Hilscher et al., 2013).
It is well known that Landau damping in a collisionless plasma is due to the effect of the destructive interference of a continuous spectrum of singular eigenmodes (the CaseVan Kampen modes).
Ng et al. (1999) have shown that a LenardBernstein collisional operator changes the spectrum of the linear Vlasov problem by replacing the singular continuous spectrum with
a set of proper discrete eigenmodes, and that the Landau damping is recovered as a discrete mode (along with other modes). In this context, Figure 3 clearly shows that the right
damping rate (consistent with Landau damping) can be recovered.
One has to keep in mind, however, that although the macroscopic nature of the plasma has been preserved, the
microscopic information associated with the high order moments of the distribution function is irreversibly modified by applying the collisional operator.
On the other hand, in most simulations of a kinetic collisionless plasma, the use of an artificial collisional operator is necessary because filamentation is an intrinsic feature.
Once again, the underlying assumption (to be verified through a convergence study) is that the use of artificial collisions will not affect the macroscopic evolution of the system.
Obviously, a PIC code is not immune to filamentation problems: the use of a discrete number of superparticles implies a finite
resolution in velocity space. However, the fact that any velocity value
can be assigned to a single particle and hence that the discretized phasespace is not gridded in any standard way makes it difficult to quantify
the relationship between the number of particles and the actual resolution in velocity space. On the other hand,
it is very wellknown that PIC simulations have difficulties in reproducing the fine details of a distribution function, such as high tails.
It is therefore legitimate to ask how the recurrence phenomenon shown in Figure 2 translates to PIC simulations.
In Figure 4, we show that although PIC simulations do not present recurrence similarly to Vlasov codes,
the lack of resolution in velocity space still manifests itself and produces inaccurate results.
Top, middle and bottom panels show PIC results (black line) for number of particles per cell , and , respectively. The red line
is a reference solution calculated with the FH code (with ).
For all three of these cases, the correct result is lost at some point, and the solution becomes essentially noise.
A useful way of understanding this deviation from the correct solution is to look at the signaltonoise ratio.
In Figure 5, we show the time evolution of the Landau damping test
for different values of the initial perturbation (see Section 5 for the discussion of the
initialization) from PIC simulations with 2,000,000 particles per cell.
The black line indicates the noise level, which has been calculated as the maximum value of in a case without initial perturbation ().
Blue and red lines denote simulations with and , respectively. Although the two simulations have the same number of particles per cell,
starting with a lower initial amplitude (blue line) (i.e. at a lower signaltonoise level), clearly impacts the result: the Landau damping is almost immediately lost for .
In other words, the noise ground level sets the amplitude of a signal at which the simulation loses any physical interpretation.
Of course, in the Hermite method, there is no noise ground level and the equivalent simulation (shown in Figure 2)
is independent of the value of the initial perturbation (given that the value is small enough to be in the linear regime).
A last remark is in order with regards to collisionality and PIC codes.
The discrete nature of the superparticles and the inevitable presence of numerical noise can in some sense be associated to the effect of a (weakly) collisional operator,
although its effect is not easily quantifiable.
The major disadvantage of PIC relative to the FH method is the bad scaling of the noise with respect to the number of particles, which, in certain applications,
might imply the need for significant computational resources to access a true collisionless regime.
Conversely, the collision operator does not affect the performance of the FH code and a convergence study on the collision rate can be conducted quickly.
5 Results
In this section we compare the computational performance of the FH spectral method with a standard explicit PIC code for the following four test cases: Langmuir wave, Landau damping, ionacoustic wave and twostream instability. These are the same test cases discussed in Reference (Chen et al., 2011) and routinely used to benchmark kinetic codes. Additionally, we present a result for an example of plasma echo. Each test tackles different aspects of kinetic plasma physics and the related difficulties encountered in the numerical simulations of collisionless kinetic plasmas.
In order to make the
comparison as fair as possible, the Poisson solver in the PIC code and the FH code are identical, i.e., the Poisson equation is solved by means of a Fourier decomposition in both cases.
The PIC code employs linear interpolation, usually referred to as ’CloudinCell’ (CIC) (Birdsall and Langdon, 2005).
Note that higher order interpolation schemes have been proposed, for instance, in (Lewis, 1970; Evstatiev and Shadwick, 2013).
Since the focus of this work is on the discretization in velocity space, i.e., the comparison
between the spectral Hermite basis and the use of superparticles, all simulations are made for the same choice of time step and grid size.
The comparison is characterized by the following three metrics: for each method we calculate
1) the error with respect to a ’reference’ highly accurate solution as a function of CPU time and velocity discretization (number of
Hermite modes for FH and number of particles per cell for PIC); 2) the error with respect to the ’previous’ less accurate solution;
3) the efficacy defined as the inverse of the product of CPU time and error.
The error used for all the runs is calculated as the norm of the difference between two solutions, averaged in time.
The error in 2) is what is actually used by a user who is performing a convergence study to decide when the solution is accurate enough.
The efficacy is a useful indicator of the costeffectiveness of an algorithm. It measures whether an additional cost in terms of CPU time is compensated by a gain in terms of accuracy.
Clearly an algorithm performs well if the efficacy increases notably with increasing CPU time. In this regard we can anticipate that the PIC algorithm, by construction, performs badly in terms of efficacy
since the error scales as , while the computing time scales roughly linearly with .
Therefore, the efficacy scales as the inverse of the square root of the CPU time,
that is, it actually decreases with increasing
CPU time. Hence, from a pure costeffectiveness point of view, it is never advantageous to increase the number of particles in a PIC code to reduce the error.
On the other hand, one is often forced to have a large number of particles such that
the physical signal is above the noise level (see e.g. (Camporeale and Burgess, 2011)).
For all the cases discussed below, we initialize the electrons (ions) with a Maxwellian distribution function with thermal velocity ().
For the Langmuir wave, the Landau damping and the twostream instability tests
the ions constitute a fixed background, while for the study of the ionacoustic wave they evolve.
The initial electric field is initialized as:
(38) 
is the box length ( for all runs), and is the amplitude of the initial perturbation. In Fourier space such initialization corresponds to:
(39) 
and the density is initialized consistently. For all runs the number of Fourier modes is equal to and the time step is . In all runs presented here we have used the CN scheme for the Xshift of the FH code. All the codes are written in MATLAB and run on an Intel Xeon 3.40 GHz Linux box.
5.1 Langmuir wave
The parameters are chosen as follows: , . In this case, the electrons are relatively cold and therefore are not subject to Landau damping. The initial perturbation causes the electrons to oscillate collectively around the fixed ion background. Figure 6 shows the error achieved with respect to the reference solution (red circles) and the previous solution (black circles) for PIC (topleft panel) and FH codes (topright panel), as a function of the resolution in velocity space. The theoretical scaling (black line) is clearly recovered for PIC. The difference in errors are huge, such that the most accurate solution obtained with PIC (102,400 particles per cell) is still less accurate than the solution obtained with only 4 Hermite modes. One can also notice that for the FH code the error with respect to the reference solution flattens for . The way this result is reflected in terms of CPU time is shown in the bottomleft panel of Figure 6. Here the CPU time is reported on the horizontal axis in seconds, and the error (with respect to the reference solution) is on the vertical axis. The black circles are PIC results, while the red circles are FH results. Strikingly, the most accurate solution for the FH code and the least accurate solution for the PIC code are obtained with roughly the same CPU time! The accuracy of these solutions differs by about 10 orders of magnitude. Finally, the efficacy is plotted on the bottomright panel of Figure 6 (black circles are PIC results, red circles are FH results). As anticipated the PIC efficacy decreases by increasing CPU time, and is consistent with the theoretical prediction. On the other hand, the FH efficacy increases by 7 orders of magnitude when the CPU time increases by less than a factor of 10.
5.2 Landau damping
Landau damping is a classical kinetic effect in warm plasmas. It is due to the effect of particles in resonance with an initial wave perturbation that results in the damping of the wave (Brambilla, 1998). The parameters are chosen as follows: , . Landau damping, as discussed in Section 4, is a particularly challenging problem for kinetic codes, because of the continuous filamentation in velocity space. For the FH code, this translates to the fact that more Hermite modes are needed as the simulation evolves, with respect to the Langmuir wave case. Figure 7 summarizes the results in the same format of Figure 6. The top left and right panels show the convergence study of error as a function of and for the PIC and FH codes, respectively. The bottom left panel shows the error as a function of CPU time, and the bottom right panel represents the efficacy as a function of CPU time. Similarly to the Langmuir wave case, the FH code greatly outperforms the PIC code: the efficacy for the FH method increases by about 6 orders of magnitude when the CPU time increases only by a factor of 6, in the interval from 10 to 60 seconds. Conversely, the difference in error between the two methods is as large as 8 orders of magnitude, for equal CPU time.
5.3 Ionacoustic wave
The ionacoustic wave is a typical example of a multiscale phenomenon in kinetic plasma physics, since it requires both electrons and ions dynamics.
Hence, the distribution function of both species is evolved. The ion acoustic wave is initialized by perturbing the ion population only.
The mass ratio between the two species is equal to 1836, and the ratio of electron and ion temperature is .
This results in a ratio of electron to ion thermal velocities of
. The initial perturbation is . Hence, the simulation is performed in a nonlinear regime.
For this case, we have used a collisionality with .
Figure 8 shows the time evolution of for PIC (black line, ) and FH (red line, , ). The results are in good agreement,
although PIC results are still very noisy even with such a large number of particles per cell.
Figure 9 shows the errors and the efficacy with the same format as before.
Once again, the scaling is approximately recovered for PIC, and the performance of the two methods differ
by orders of magnitude. For instance, to obtain an error of the order of the PIC code takes about 3 orders of magnitude
more CPU time than FH.
5.4 Twostream instability
The two stream instability is excited when the distribution function of a species is formed by two populations streaming in opposite directions with a large enough relative drift velocity. We initialize the electron distribution function as two counterstreaming Maxwellians with equal temperature:
(40) 
where is the drift velocity.
In principle, a distribution function of the form (40) would require a large number of Hermite polynomials in order to be described accurately
by the expansion in Eq. (5). However, a way to overcome this difficulty is to split the distribution function in (40)
into two distinct populations corresponding to each drifting Maxwellian,
and to solve Eq. (9) separately for each population. As already anticipated, including a free parameter
in the definition of the argument allows the description of a drifting Maxwellian by means of only one FH coefficient in Eq.(5) (by simply choosing ).
For this test case, we have chosen the following parameters: .
Figure 10 shows the evolution in time of obtained by PIC simulations with a different number of particles per cell (magenta=200; red=3,200; black=102,400; blue=409,600).
The black dashed line represents the theoretical growth rate of the instability, calculated by linear theory.
One can notice that the convergence of the PIC simulation is very slow, particularly due to the fact that the initial transient that gives rise to the instability
is very dependent on the number of particles, i.e. the noise level. Also, although each run saturates at approximately the same level, the saturation time is different, making the postsaturation
phase not comparable between different runs (the smaller box in Figure 10 shows a zoomin of the postsaturation dynamics, in linear scale).
For this reason, the errors used in this section have been calculated only in the linear regime region, i.e. for times smaller than 5.
Interestingly, capturing fine structures in the distribution function (in velocity space) becomes important, in the twostream instability, only in the postsaturation phase.
This can be clearly appreciated in Figure 11. Here, we plot the time evolution of , for FH simulations with different values of .
Even with a low number of Hermite modes (, red line), the linear stage of the instability is well captured. However, the simulation becomes numerically unstable
in the postsaturation phase. By increasing the number of Hermite modes (, blue line; , black line), the time at which the
simulation becomes unstable is postponed. This indicates that a large number of Hermite modes is needed because fine structures in the velocity space continue
to form after the instability has saturated, i.e. the plasma is developing filamentation. The solution indicated with circles is obtained with using a collisional operator with .
The numerical collisions do not change the result macroscopically, but allow longer simulations without the need of a prohibitively large number of Hermite modes.
The fact that the FH simulation becomes numerically unstable if collisions are not applied may be viewed as a shortcoming of the method.
On the contrary, we regard this feature as beneficial, especially when compared with the PIC code. In fact, the appearance of the numerical instability is a clear signal
of the lack of resolution in velocity space. The straightforward approach to overcoming this problem is by increasing and/or introducing artificial collisions. On the other hand,
in the PIC simulation the only way to understand that the results are not converged is to perform an expensive convergence study. In other words,
the PIC simulation will keep running, potentially for a long time, wasting computational resources, and producing results that have no physical meaning.
Figure 12 presents the errors and efficacy for the twostream instability case. The results are qualitatively similar to the previous cases, once again highlighting
the high performance of the FH method with respect to PIC.
5.5 Plasma echo
In contrast to the previous tests, the rationale for this case is not to compare results between PIC and FH methods but to provide an example
of a fundamental kinetic plasma phenomenon that is very hard to capture with PIC, but is still straightforward to simulate with the FH code.
The basic theory behind plasma echo has been studied in the seminal paper by Gould et al. (1967), and investigated experimentally by Malmberg et al. (1968).
Since then the theory has been revised and rediscussed by several authors (e.g. Brackbill, 1972; Best, 1974; Porkolab and Chang, 1978; Hou et al., 2011; Mouhot and Villani, 2011; Pezzi et al., 2013).
Plasma echo has been suggested as the ’archetype’ of numerical tests for onedimensional Vlasov codes in (Galeotti et al., 2006).
The plasma echo mechanism is briefly explained in (Porkolab and Chang, 1978) as the following.
An initial electric field perturbation with wavevector is initialized at time in a stable Maxwellian plasma.
This perturbation damps according to Landau damping and modulates the distribution function through a term of the form .
At time a second perturbation is initialized, with wavevector , which also damps according to
its linear damping rate. This second perturbation will act on the distribution function via a first order term of the form .
However, a second order contribution will arise from the interaction between the two perturbations, with the form .
While the first order terms will damp in time because of phasemixing, the second order term loses its dependence at time , and therefore will not phasemix around that time.
As a consequence, a new perturbation with wavevector will appear, peaking at time , and successively causing Landau damping.
The plasma echo is a nonlinear effect due to the interaction of two small amplitude damping waves. Interestingly, the time can be made much larger than and such that
the reformation of the electric field after the two parent perturbations have long been damped can effectively appear as an ’echo’.
The simulation is initialized with a perturbation in . At time , a second perturbation
is generated, again with amplitude but in .
The resulting echo appears at time for the Fourier mode .
Figure 13 shows the amplitudes of , , and respectively in the top, middle and bottom panels.
Black solid lines are for the FH code, and red lines are for PIC. The dashed black line indicates the theoretical
Landau damping of each mode.
Clearly, the FH results are in excellent agreement with the theoretical prediction discussed above, showing the peaking of
at time .
We note that we have used and in this example and some recurrence can be seen in the FH code for
The PIC simulations are run with . This is of course an extremely large number of particles
per cell, larger than what is routinely used. However, the PIC simulation completely misses the echo phenomenon, simply because a noise ground level comparable
with the echo signal is generated for since the beginning of the simulation (bottom panel), and both and have reached the noise ground level, by the time the echo should be generated.
It is important to emphasize that the plasma echo is not merely an academic exercise.
The physics of wavewave interactions for small amplitude perturbations constitutes the building block of the weakturbulence theory, important, for instance, for magnetized plasmas in the solar wind.
6 Conclusions
The focus of this work was to compare the discretization of the velocity space in the VlasovPoisson system based on the Hermite expansion with the use of superparticles adopted in the standard PIC method.
We have described the spectral method based on the FH expansion previously discussed in Schumer and Holloway (1998). We have introduced two important improvements with respect to the work in (Schumer and Holloway, 1998). First, we have employed a semiimplicit time discretization. Our scheme is still second order accurate in time, but we have shown that it is numerically more stable than the explicit scheme previously used. Second, we have discussed the use of an artificial collisional term, whose scope is to damp high order Hermite modes. The collision term has a double effect: it prevents the recurrence effect in perturbations which are damping, and avoids the growth of higher order modes that can lead to numerical instabilities.
As such, it is an essential ingredient of the Hermite spectral method as it takes care of the unavoidable filamentation problem typical of collisionless plasmas.
The comparisons between the FH method and PIC have been carried out from the perspective of computational efficiency and efficacy. The efficiency is given by the error achieved for a given CPU time, while the efficacy is a measure of the costeffectiveness of the algorithm and is defined as the inverse of the product of the error and the CPU time.
We have studied the following standard kinetic plasma test cases: Langmuir waves, Landau damping, ionacoustic wave, twostream instability. For each test case, we have evaluated the error of the solutions obtained with the two methods, and the CPU time for each run. In all these examples, the FH method obtains a much accurate solution than PIC, at a fraction of the computational cost.
For instance, for the twostream instability case, the FH method produces a result that is about 10 order of magnitudes more accurate in about 1/50 of the CPU time taken by PIC. We have also studied an example of the plasma echo phenomenon. For this case, we have shown that the standard PIC method cannot capture the nonlinear wavewave interactions underlying the echo physics, despite using an astonishing 2,000,000 particles per cell.
It has to be emphasized that whilst the PIC method has been investigated and developed for the last several decades, methods based on a spectral expansion of the distribution function (either Hermite or Fourier) are much less mature. For this reason, in our comparisons we have not taken into account any strategy that could reduce the PIC noise. Notably, methods based on f or remapping of the distribution function have been successful to reduce particle noise, and hence enhance the accuracy of PIC simulations (see, e.g. (Chehab et al., 2005; Wang et al., 2011; Deng and Fu, 2014)).
Of course, this aspect lies in the realm of algorithm optimization. For the sake of a fair comparison, however, we reckon that it is preferable to deal with the unoptimized method, as this exposes the real limitations of each algorithm. Furthermore, for what concerns PIC, there are many different ways to reduce the noise. This opens the possibility of trying and comparing several different algorithms, each with its own pros and cons, and it is likely that different test cases will be more or less favorable for a given optimization.
On the other hand, the optimization of the
FourierHermite method has practically never been investigated, although one could envision a strategy where the Hermite basis changes in time in order to use the least number of modes.
Another consideration is that the PIC method has, by construction, an unoptimal scaling of the error with the CPU time: its efficacy decreases with increasing CPU time.
Denoising techniques can be thought as an attempt to counteract such bad scaling. Of course, this might be beneficial for practical applications, but should not change the overall scaling of
the method with the number of particles. It is also not obvious whether such techniques will actually increase the PIC performance from an efficacy perspective.
As a final comment, we want to remark some of the limitations and challenges of the FH method relative to PIC. In general, while the benchmark tests studied here are standard in the PIC community, nearequilibrium Maxwellian plasmas favor FH over PIC (since the Hermite basis naturally links to Maxwellians). On the other hand, FH is not well suited for cold plasmas/beams, which are instead treated very effectively by PIC. A different spectral basis could be sought for these cases. Similarly, if the plasma develops very complicated structures in the distribution function, many Hermite modes might be necessary to achieve convergence. This can become a serious a problem in multidimensional applications, unless some form of basis optimization is successfully developed. This is likely the most critical need of FH for application to collisionless plasmas. We have also observed that numerical instabilities can arise in the FH method for long wavelengths and large initial perturbations. We have verified that a
possible fix/mitigation of this problem is to introduce a kdependent collision operator (a detailed analysis of this approach will be presented elsewhere).
On the contrary, the PIC algorithm is very robust.
In conclusion, considering the computational challenges associated with multiscale plasma physics and the upcoming exascale computing era, it is obvious that numerical methods that can use efficiently the everincreasing available computational power are critical. Despite the proofofprinciple nature of our study, our results indicate that a method based on a spectral expansion of the velocity part of the distribution function should be given serious consideration for this task. In addition to potential orders of magnitude higher efficacy (relative to PIC), we believe that another major strength of the Hermite spectral approach is its natural ability to bridge between regions where the plasma behaves like a
fluid (hence it is characterized by a low number of moments) and regions where the microscopic/kinetic physics is important (thus requiring a higher number of moments). Strategies that can accelerate the convergence of the Hermite basis will be critical for the success of the FH method in multidimensions, and we hope that this paper can stimulate some research in this area.
7 Acknowledgement
The authors would like to thank L. Chacon, W. Daughton, C. Huang, V. Roytershteyn, X. Tang for useful conversations. This work was funded by the Laboratory Directed Research and Development program (LDRD), U.S. Department of Energy Office of Science, Office of Fusion Energy Sciences, under the auspices of the National Nuclear Security Administration of the U.S. Department of Energy by Los Alamos National Laboratory, operated by Los Alamos National Security LLC under contract DEAC5206NA25396.
References
 Armstrong et al. [1970] T.P. Armstrong, R.C. Harding, G. Knorr, and D. Montgomery. Solution of Vlasov’s equation by transform methods. Methods Comput. Phys., 9:29, 1970.
 Banks and Hittinger [2010] J.W. Banks and J.A.F. Hittinger. A new class of nonlinear finitevolume methods for Vlasov simulation. Plasma Science, IEEE Transactions on, 38(9):2198–2207, 2010.
 Batchelor et al. [2007] D.A. Batchelor, M. Beck, A. Becoulet, R. V. Budny, C. S. Chang, P. H. Diamond, J. Q. Dong, G. Y. Fu, A. Fukuyama, T. S. Hahm, D. E. Keyes, Y. Kishimoto, S. Klasky, L. L. Lao, K. Li, Z. Lin, B. Ludaescher, J. Manickam, N. Nakajima, T. Ozeki, N. Podhorszki, W. M. Tang, M. A. Vouk, R. E. Waltz, S. J. Wang, H. R. Wilson, X. Q. Xu, M. Yagi, and F. Zonca. Simulation of Fusion Plasmas: Current Status and Future Direction. Plasma Science and Technology, 9(3):312, 2007. URL http://phoenix.ps.uci.edu/zlin/bib/batchelor07.pdf.
 Besse and Sonnendrücker [2003] N. Besse and E. Sonnendrücker. SemiLagrangian schemes for the Vlasov equation on an unstructured mesh of phase space. Journal of Computational Physics, 191:341–376, November 2003. doi: 10.1016/S00219991(03)003188.
 Best [1974] R.W.B. Best. Nonlinear plasma oscillations in terms of Van Kampen modes: The initialvalue problem. Physica, 74(1):183–195, 1974.
 Birdsall and Langdon [2005] C.K. Birdsall and A. B. Langdon. Plasma Physics Via Computer Simulaition. CRC Press, 2005.
 Bowers et al. [2009] K.J. Bowers, B.J. Albright, L. Yin, W. Daughton, V. Roytershteyn, B. Bergen, and T.J.T. Kwan. Advances in petascale kinetic plasma simulation with VPIC and Roadrunner. In Journal of Physics: Conference Series, volume 180, page 012055. IOP Publishing, 2009.
 Bracewell [1999] R. Bracewell. The Fourier Transform And Its Applications. McGrawHill Science/Engineering/Math, 1999.
 Brackbill [1972] J.U. Brackbill. Virtual Plasma Wave Echoes. Physics of Fluids, 15:1358, 1972.
 Brambilla [1998] M. Brambilla. Kinetic Theory of Plasma Waves: Homogeneous Plasmas, volume 96. Oxford University Press, 1998.
 Byers [1970] J.A. Byers. Noise suppression techniques in macroparticle models of collisionless plasmas. In Jay P. Boris and Ramy A. Shanny, editors, Fourth Conference on Numerical Simulation of Plasmas, 1970.
 Camporeale and Burgess [2011] E. Camporeale and D. Burgess. The dissipation of solar wind turbulent fluctuations at electron scales. The Astrophysical Journal, 730(2):114, 2011.
 Camporeale et al. [2006] E. Camporeale, G. L. Delzanno, G. Lapenta, and W. Daughton. New approach for the study of linear Vlasov stability of inhomogeneous systems. Physics of Plasmas, 13(9):092110, September 2006. doi: 10.1063/1.2345358.
 Camporeale et al. [2013] E. Camporeale, G. L. Delzanno, B. K. Bergen, and J. D. Moulton. On the velociy space discretization for the VlasovPoisson system: comparison between Hermite spectral and ParticleinCell methods. Part 2: fullyimplicit scheme. 2013.
 Canosa et al. [1974] J. Canosa, J. Gazdag, and J.E. Fromm. The recurrence of the initial state in the numerical solution of the Vlasov equation. Journal of Computational Physics, 15(1):34–45, 1974.
 Carrillo and Vecil [2007] J. A. Carrillo and F. Vecil. Nonoscillatory Interpolation Methods Applied to VlasovBased Models. SIAM J. Sci. Comput., 29(3):1179–1206, May 2007. ISSN 10648275. doi: 10.1137/050644549. URL http://dx.doi.org/10.1137/050644549.
 Chehab et al. [2005] JP Chehab, Albert Cohen, D Jennequin, JJ Nieto, Ch Roland, and JR Roche. An adaptive ParticleInCell method using multiresolution analysis. Numerical methods for hyperbolic and kinetic problems, IRMA Lect. Math. Theor. Phys, 7:29–42, 2005.
 Chen et al. [2011] G. Chen, L. Chacón, and D.C. Barnes. An energyand chargeconserving, implicit, electrostatic particleincell algorithm. Journal of Computational Physics, 230(18):7018–7036, 2011.
 Cheng and Knorr [1976] C. Z. Cheng and G. Knorr. The integration of the Vlasov equation in configuration space. Journal of Computational Physics, 22:330–351, November 1976. doi: 10.1016/00219991(76)90053X.
 Crank and Nicolson [1947] J. Crank and P. Nicolson. A practical method for numerical evaluation of solutions of partial differential equations of the heatconduction type. In Mathematical Proceedings of the Cambridge Philosophical Society, volume 43, pages 50–67. Cambridge Univ Press, 1947.
 Crouseilles and Sonnendrücker [2007] G. Crouseilles, N.and Latu and E. Sonnendrücker. Hermite spline interpolation on patches for parallelly solving the VlasovPoisson equation. International Journal of Applied Mathematics and Computer Science, 17(3):335–349, 2007.
 Crouseilles et al. [2010] N. Crouseilles, M. Mehrenberger, and E. Sonnendrücker. Conservative semiLagrangian schemes for Vlasov equations. Journal of Computational Physics, 229:1927–1953, March 2010. doi: 10.1016/j.jcp.2009.11.007.
 Daughton et al. [2011] W. Daughton, V. Roytershteyn, H. Karimabadi, L. Yin, B.J. Albright, B. Bergen, and K.J. Bowers. Role of electron physics in the development of turbulent magnetic reconnection in collisionless plasmas. Nature Physics, 7(7):539–542, 2011.
 Denavit and Kruer [1971] J. Denavit and W.L. Kruer. Comparison of Numerical Solutions of the Vlasov Equation with Particle Simulations of Collisionless Plasmas. Physics of Fluids, 14:1782, 1971.
 Denavit and Walsh [1981] J. Denavit and J.M. Walsh. Nonrandom Initializations of Particle Codes. Comments on Plasma Phys. Controlled Fusion, 6:209–223, 1981.
 Deng and Fu [2014] Wenjun Deng and GuoYong Fu. Optimization by marker removal for f particle simulations. Computer Physics Communications, 185(1):96–105, 2014.
 Duclous et al. [2012] R. Duclous, B. Dubroca, F. Filbet, et al. Analysis of a high order finite volume scheme for the VlasovPoisson system. DCDSS, 5:283–305, 2012.
 Eliasson [2003] B. Eliasson. Numerical modelling of the twodimensional Fourier transformed VlasovMaxwell system. Journal of Computational Physics, 190(2):501–522, 2003.
 Engelmann et al. [1963] F. Engelmann, M. Feix, E. Minardi, and J. Oxenius. Nonlinear Effects from Vlasov’s Equation. Physics of Fluids, 6:266–275, February 1963. doi: 10.1063/1.1706724.
 Evstatiev and Shadwick [2013] EG Evstatiev and BA Shadwick. Variational formulation of particle algorithms for kinetic plasma simulations. Journal of Computational Physics, 2013.
 Filbet et al. [2001] F. Filbet, E. Sonnendrücker, and P. Bertrand. Conservative Numerical Schemes for the Vlasov Equation. Journal of Computational Physics, 172:166–187, September 2001. doi: 10.1006/jcph.2001.6818.
 Galeotti et al. [2006] L Galeotti, F Califano, and F Pegoraro. ”Echography” of Vlasov codes. Physics Letters A, 355(4):381–385, 2006.
 Gould et al. [1967] R_W Gould, TM O’Neil, and JH Malmberg. Plasma wave echo. Physical Review Letters, 19(5):219–222, 1967.
 Grad [1949] H. Grad. On the kinetic theory of rarefied gases. Communications on pure and applied mathematics, 2(4):331–407, 1949.
 Grant and Feix [1967] F. C. Grant and M. R. Feix. FourierHermite Solutions of the Vlasov Equations in the Linearized Limit. Physics of Fluids, 10:696–702, April 1967. doi: 10.1063/1.1762177.
 Guo and Qiu [2012] W. Guo and J.M. Qiu. Hybrid semiLagrangian finite elementfinite difference methods for the Vlasov equation. Journal of Computational Physics, 2012.
 Haynes et al. [2013] C.T. Haynes, D. Burgess, and E. Camporeale. Reconnection and electron temperature anisotropy in electron scale plasma turbulence. arXiv preprint arXiv:1304.1444, 2013.
 Heath et al. [2012] RE Heath, IM Gamba, PJ Morrison, and C Michler. A discontinuous Galerkin method for the Vlasov–Poisson system. Journal of Computational Physics, 231(4):1140–1174, 2012.
 Henrici [1962] P. Henrici. Discrete variable methods in ordinary differential equations. New York: Wiley, 1962, 1, 1962.
 Hilscher et al. [2013] PP Hilscher, K Imadera, JQ Li, and Y Kishimoto. The effect of weak collisionality on damped modes and its contribution to linear mode coupling in gyrokinetic simulation. Physics of Plasmas, 20:082127, 2013.
 Hockney and Eastwood [2010] R.W. Hockney and J.W. Eastwood. Computer simulation using particles. CRC Press, 2010.
 Holloway [1996] J. P. Holloway. Spectral velocity discretizations for the VlasovMaxwell equations. Transport Theory and Statistical Physics, 25:1–32, January 1996. doi: 10.1080/00411459608204828.
 Hou et al. [2011] YW Hou, ZW Ma, and MY Yu. The plasma wave echo revisited. Physics of Plasmas, 18:012108, 2011.
 Imadera et al. [2009] K. Imadera, Y. Kishimoto, D. Saito, J. Li, and T. Utsumi. A numerical method for solving the VlasovPoisson equation based on the conservative IDO scheme. Journal of Computational Physics, 228:8919–8943, December 2009. doi: 10.1016/j.jcp.2009.09.008.
 Klimas [1987] A. J. Klimas. A method for overcoming the velocity space filamentation problem in collisionless plasma model solutions. Journal of Computational Physics, 68:202–226, January 1987.
 Klimas [1994] A. J. Klimas. A Splitting Algorithm for Vlasov Simulation with Filamentation Filtration. Journal of Computational Physics, 110:150–163, 1994.
 Klimas [1983] Alexander J Klimas. A numerical method based on the FourierFourier transform approach for modeling 1D electron plasma evolution. Journal of Computational Physics, 50(2):270–306, 1983.
 Landau [1946] L. Landau. On the vibrations of the electronic plasma. J. Phys.(UssR), 10(1):25–34, 1946.
 Lapenta and Chacón [2006] G. Lapenta and L. Chacón. Costeffectiveness of fully implicit moving mesh adaptation: A practical investigation in 1D. Journal of computational physics, 219(1):86–103, 2006.
 Le Bourdiec et al. [2006] S. Le Bourdiec, F. De Vuyst, and L. Jacquet. Numerical solution of the Vlasov–Poisson system using generalized Hermite functions. Computer physics communications, 175(8):528–544, 2006.
 Lenard and Bernstein [1958] A. Lenard and I.B. Bernstein. Plasma Oscillations with Diffusion in Velocity Space. Phys. Rev., 112:1456–1459, Dec 1958. doi: 10.1103/PhysRev.112.1456. URL http://link.aps.org/doi/10.1103/PhysRev.112.1456.
 Lewis [1970] H Ralph Lewis. EnergyConserving Numerical Approximations for Vlasov Plasmas. Journal of Computational Physics, 6:136, 1970.
 Malmberg et al. [1968] JH Malmberg, CB Wharton, RW Gould, and TM O’neil. Plasma wave echo experiment. Physical Review Letters, 20(3):95–97, 1968.
 Mangeney et al. [2002] A. Mangeney, F. Califano, C. Cavazzoni, and P. Travnicek. A numerical scheme for the integration of the Vlasov–Maxwell system of equations. Journal of Computational Physics, 179(2):495–538, 2002.
 Matteini et al. [2012] L. Matteini, P. Hellinger, S. Landi, P.M. Trávníček, and M. Velli. Ion kinetics in the solar wind: Coupling global expansion to local microphysics. Space science reviews, 172(14):373–396, 2012.
 Mouhot and Villani [2011] C. Mouhot and C. Villani. On landau damping. Acta mathematica, 207(1):29–201, 2011.
 Nevins et al. [2005] WM Nevins, GW Hammett, AM Dimits, W Dorland, and DE Shumaker. Discrete particle noise in particleincell simulations of plasma microturbulence. Physics of plasmas, 12(12):122305–122305, 2005.
 Ng et al. [1999] CS Ng, A Bhattacharjee, and F Skiff. Kinetic eigenmodes and discrete spectrum of plasma oscillations in a weakly collisional plasma. Physical review letters, 83(10):1974, 1999.
 Parker and Dellar [2012] J. Parker and P. Dellar. Hermite expansions with hypercollisionality for velocity space degrees of freedom in iontemperaturegradient driven instabilities. Bulletin of the American Physical Society, 57, 2012.
 Pezzi et al. [2013] O. Pezzi, F. Valentini, D. Perrone, and P. Veltri. Eulerian simulations of collisional effects on electrostatic plasma waves. arXiv preprint arXiv:1307.0683, 2013.
 Porkolab and Chang [1978] M Porkolab and RPH Chang. Nonlinear wave effects in laboratory plasmas: A comparison between theory and experiment. Reviews of Modern Physics, 50(4):745, 1978.
 Qiu and Christlieb [2010] J.M. Qiu and A. Christlieb. A conservative high order semiLagrangian WENO method for the Vlasov equation. Journal of Computational Physics, 229:1130–1149, February 2010. doi: 10.1016/j.jcp.2009.10.016.
 Qiu and Shu [2011] J.M. Qiu and C.W. Shu. Positivity preserving semiLagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the VlasovPoisson system. Journal of Computational Physics, 230:8386–8409, September 2011. doi: 10.1016/j.jcp.2011.07.018.
 Rossmanith and Seal [2011] J. A. Rossmanith and D. C. Seal. A positivitypreserving highorder semiLagrangian discontinuous Galerkin scheme for the VlasovPoisson equations. Journal of Computational Physics, 230:6203–6232, July 2011. doi: 10.1016/j.jcp.2011.04.018.
 Schumer and Holloway [1998] J.W. Schumer and J.P. Holloway. Vlasov simulations using velocityscaled Hermite representations. Journal of Computational Physics, 144(2):626–661, 1998.
 Shoucri [2008] M. Shoucri. Eulerian codes for the numerical solution of the Vlasov equation. Communications in Nonlinear Science and Numerical Simulations, 13:174–182, February 2008. doi: 10.1016/j.cnsns.2007.04.004.
 Shoucri and Knorr [1974] M. Shoucri and G. Knorr. Numerical Integration of the Vlasov Equation. Journal of Computational Physics, 14:84, January 1974. doi: 10.1016/00219991(74)900060.
 Siminos et al. [2011] E. Siminos, D. Bénisti, and L. Gremillet. Stability of nonlinear VlasovPoisson equilibria through spectral deformation and FourierHermite expansion. Phys. Rev. E., 83(5):056402, May 2011. doi: 10.1103/PhysRevE.83.056402.
 Sonnendrücker et al. [1999] E. Sonnendrücker, J. Roche, P. Bertrand, and A. Ghizzo. The SemiLagrangian Method for the Numerical Resolution of the Vlasov Equation. Journal of Computational Physics, 149:201–220, March 1999. doi: 10.1006/jcph.1998.6148.
 Strang [1968] G. Strang. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis, 5(3):506–517, 1968.
 Sydora [1999] RD Sydora. Lownoise electromagnetic and relativistic particleincell plasma simulation models. Journal of computational and applied mathematics, 109(1):243–259, 1999.
 Tang [1993] T. Tang. The Hermite spectral method for Gaussiantype functions. SIAM Journal on Scientific Computing, 14(3):594–606, 1993.
 Tang [2008] W.M. Tang. Scientific and computational challenges of the fusion simulation project (FSP). In Journal of Physics: Conference Series, volume 125, page 012047. IOP Publishing, 2008.
 Umeda et al. [2006] T. Umeda, M. AshourAbdalla, and D. Schriver. Comparison of numerical interpolation schemes for onedimensional electrostatic Vlasov code. Journal of Plasma Physics, 72:1057, December 2006. doi: 10.1017/S0022377806005228.
 Valentini et al. [2008] F. Valentini, P. Veltri, F. Califano, and A. Mangeney. Crossscale effects in solarwind turbulence. Physical review letters, 101(2):025006, 2008.
 Verboncoeur [2005] J. P. Verboncoeur. Particle simulation of plasmas: review and advances. Plasma Physics and Controlled Fusion, 47(5A):A231, 2005.
 Walker [2009] E. Walker. The real cost of a CPU hour. Computer, 42(4):35–41, 2009.
 Wang et al. [2011] Bei Wang, GH Miller, and Phillip Colella. A ParticleInCell method with adaptive phasespace remapping for kinetic plasmas. SIAM Journal on Scientific Computing, 33(6):3509–3537, 2011.
 Whealton et al. [1986] JH Whealton, RW McGaffey, and PS Meszaros. A finite difference 3D PoissonVlasov algorithm for ions extracted from a plasma. Journal of Computational Physics, 63(1):20–32, 1986.