On the velocity space discretization for the Vlasov-Poisson system: comparison between Hermite spectral and Particle-in-Cell methods.
Part 1: semi-implicit scheme
We discuss a spectral method for the numerical solution of the Vlasov-Poisson system where the velocity space is decomposed by means of an Hermite basis. We describe a semi-implicit 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 cost-effectiveness of this method are compared to a Particle-in-Cell (PIC) method in the case of a two-dimensional phase space. The following examples are discussed: Langmuir wave, Landau damping, ion-acoustic wave, two-stream 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.
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 self-consistently 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 Particle-In-Cell (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 super-particles. The electromagnetic field is represented on a grid in the computational domain, and the super-particles 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 super-particles). 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 super-particles, 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 signal-to-noise 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 two-dimensional (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 signal-to-noise 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 re-mapping 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 non-optimal 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 resource-intensive 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 long-time 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 phase-space (namely with a sixth dimensional computational grid in space and velocity),
by means of a so-called Eulerian Vlasov code.
This has been done, usually for problems with reduced dimensionality, with a variety of techniques such as high-order 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 semi-Lagrangian 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 large-scale, 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 non-positive 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 space-dependent 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 Maxwellian-like 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 multi-scale 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 semi-implicit discretization in time of the truncated set of PDEs for the expansion coefficients.
A fully-implicit 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 Fourier-Hermite (FH) method and the PIC method was presented for simulations of Landau damping and bump-on-tail instability. Here we expand that work in three ways. First, we propose a semi-implicit 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 cost-effectiveness 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 super-particle 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 near-equilibrium 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, ion-acoustic wave, two-stream instability, and plasma echo. Finally, our conclusions are presented in Section 6.
2 Numerical Method: Fourier-Hermite basis
We study the Vlasov-Poisson system in the 2-dimensional 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:
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 self-consistently calculated through the Poisson equation:
The discretization in velocity employs the asymmetrically-weighted Hermite basis:
where is the -th Hermite polynomial. The distribution function is defined as:
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):
One can multiply Eq. (1) by , integrate in and, by using such properties, obtain:
Similarly, the Poisson equation (2) becomes:
One can notice the well-known property that the electric field is related uniquely to the 0-th 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):
and using the orthogonality of the Fourier basis, one can finally derive an infinite set of ordinary differential equations for the coefficients :
where is the index associated with the Hermite (Fourier) basis. The expression for the electric field reads:
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 non-linearity 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 Fourier-Hermite 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 fully-implicit 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
By inspection of Eq. (14) one can see that , from which the conservation of mass follows.
The single species linear momentum is defined as:
In order to prove the conservation of momentum, one can calculate the time derivative of the coefficients and :
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 Vlasov-Poisson 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 ’X-shift’ and ’V-shift’ operators, respectively. Therefore, we define the and the operators acting on as:
By following the second-order 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 fourth-order Runge-Kutta (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 semi-implicit scheme. The V shift is still treated with RK4, but the X-shift is treated implicitly.
We have tested both a first order implicit Euler scheme (BDF1) (Henrici, 1962) and a second-order Cranck-Nicolson (Crank and Nicolson, 1947)
scheme for the X-shift.
Note that the Strang splitting ensures the second order accuracy in time of the overall scheme.
In summary, the discretized form of the X-shift that advances the solution
for half time step is the following (where the superscript indicates the time step):
Of course, both the BDF1 and the CN schemes can be written in matrix form. For each Fourier component , we have:
where is the column-vector . The matrix is defined as
with for CN and for BDF1. For the BDF1 scheme the matrix is the identity matrix, and for the CN scheme is defined as:
Note that the matrices and commute in both cases (BDF1 and CN). Therefore the two consecutive X-shifts in steps 3) and 1): can be computed exactly as a single step:
We solve Eq. (34) by performing an LU decomposition of the matrix .
The V-shift is a standard fourth order Runge-Kutta integrator, reported here for completeness.
Note that the V-shift 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 semi-implicit scheme (X-shift treated with BDF1 or CN) with respect to the fully explicit scheme (X-shift 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 semi-implicit 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 well-studied example where filamentation occurs is linear Landau damping, i.e., the damping in time of an initial electric field perturbation due
to wave-particle 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 weakly-collisional 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:
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 Lenard-Bernstein collision operator (Lenard and Bernstein, 1958):
which, in the Fourier-Hermite space, reads:
An important point is that the Lenard-Bernstein operator transformed in the
Fourier-Hermite 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 fully-implicit time discretization is employed), by
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 zoom-in 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 Case-Van Kampen modes). Ng et al. (1999) have shown that a Lenard-Bernstein 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 super-particles 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 phase-space 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 well-known 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 signal-to-noise 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 signal-to-noise 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 super-particles 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.
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, ion-acoustic wave and two-stream 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 ’Cloud-in-Cell’ (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 super-particles, 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 cost-effectiveness 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 cost-effectiveness 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 two-stream instability tests the ions constitute a fixed background, while for the study of the ion-acoustic wave they evolve. The initial electric field is initialized as:
is the box length ( for all runs), and is the amplitude of the initial perturbation. In Fourier space such initialization corresponds to:
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 X-shift 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 (top-left panel) and FH codes (top-right 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 bottom-left 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 bottom-right 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 Ion-acoustic wave
The ion-acoustic wave is a typical example of a multi-scale 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 Two-stream 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 counter-streaming Maxwellians with equal temperature:
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 post-saturation phase not comparable between different runs (the smaller box in Figure 10 shows a zoom-in of the post-saturation 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 two-stream instability, only in the post-saturation 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 post-saturation 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 two-stream 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 one-dimensional 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 phase-mixing, the second order term loses its dependence at time , and therefore will not phase-mix 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 wave-wave interactions for small amplitude perturbations constitutes the building block of the weak-turbulence theory, important, for instance, for magnetized plasmas in the solar wind.
The focus of this work was to compare the discretization of the velocity space in the Vlasov-Poisson system based on the Hermite expansion with the use of super-particles 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 semi-implicit 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 cost-effectiveness 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, ion-acoustic wave, two-stream 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 two-stream 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 wave-wave 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 re-mapping 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 Fourier-Hermite 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, near-equilibrium 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 k-dependent 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 multi-scale plasma physics and the upcoming exascale computing era, it is obvious that numerical methods that can use efficiently the ever-increasing available computational power are critical. Despite the proof-of-principle 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.
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 DE-AC52-06NA25396.
- Armstrong et al.  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  J.W. Banks and J.A.F. Hittinger. A new class of nonlinear finite-volume methods for Vlasov simulation. Plasma Science, IEEE Transactions on, 38(9):2198–2207, 2010.
- Batchelor et al.  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  N. Besse and E. Sonnendrücker. Semi-Lagrangian schemes for the Vlasov equation on an unstructured mesh of phase space. Journal of Computational Physics, 191:341–376, November 2003. doi: 10.1016/S0021-9991(03)00318-8.
- Best  R.W.B. Best. Nonlinear plasma oscillations in terms of Van Kampen modes: The initial-value problem. Physica, 74(1):183–195, 1974.
- Birdsall and Langdon  C.K. Birdsall and A. B. Langdon. Plasma Physics Via Computer Simulaition. CRC Press, 2005.
- Bowers et al.  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  R. Bracewell. The Fourier Transform And Its Applications. McGraw-Hill Science/Engineering/Math, 1999.
- Brackbill  J.U. Brackbill. Virtual Plasma Wave Echoes. Physics of Fluids, 15:1358, 1972.
- Brambilla  M. Brambilla. Kinetic Theory of Plasma Waves: Homogeneous Plasmas, volume 96. Oxford University Press, 1998.
- Byers  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  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.  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.  E. Camporeale, G. L. Delzanno, B. K. Bergen, and J. D. Moulton. On the velociy space discretization for the Vlasov-Poisson system: comparison between Hermite spectral and Particle-in-Cell methods. Part 2: fully-implicit scheme. 2013.
- Canosa et al.  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  J. A. Carrillo and F. Vecil. Nonoscillatory Interpolation Methods Applied to Vlasov-Based Models. SIAM J. Sci. Comput., 29(3):1179–1206, May 2007. ISSN 1064-8275. doi: 10.1137/050644549. URL http://dx.doi.org/10.1137/050644549.
- Chehab et al.  JP Chehab, Albert Cohen, D Jennequin, JJ Nieto, Ch Roland, and JR Roche. An adaptive Particle-In-Cell method using multi-resolution analysis. Numerical methods for hyperbolic and kinetic problems, IRMA Lect. Math. Theor. Phys, 7:29–42, 2005.
- Chen et al.  G. Chen, L. Chacón, and D.C. Barnes. An energy-and charge-conserving, implicit, electrostatic particle-in-cell algorithm. Journal of Computational Physics, 230(18):7018–7036, 2011.
- Cheng and Knorr  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/0021-9991(76)90053-X.
- Crank and Nicolson  J. Crank and P. Nicolson. A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type. In Mathematical Proceedings of the Cambridge Philosophical Society, volume 43, pages 50–67. Cambridge Univ Press, 1947.
- Crouseilles and Sonnendrücker  G. Crouseilles, N.and Latu and E. Sonnendrücker. Hermite spline interpolation on patches for parallelly solving the Vlasov-Poisson equation. International Journal of Applied Mathematics and Computer Science, 17(3):335–349, 2007.
- Crouseilles et al.  N. Crouseilles, M. Mehrenberger, and E. Sonnendrücker. Conservative semi-Lagrangian schemes for Vlasov equations. Journal of Computational Physics, 229:1927–1953, March 2010. doi: 10.1016/j.jcp.2009.11.007.
- Daughton et al.  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  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  J. Denavit and J.M. Walsh. Nonrandom Initializations of Particle Codes. Comments on Plasma Phys. Controlled Fusion, 6:209–223, 1981.
- Deng and Fu  Wenjun Deng and Guo-Yong Fu. Optimization by marker removal for f particle simulations. Computer Physics Communications, 185(1):96–105, 2014.
- Duclous et al.  R. Duclous, B. Dubroca, F. Filbet, et al. Analysis of a high order finite volume scheme for the Vlasov-Poisson system. DCDS-S, 5:283–305, 2012.
- Eliasson  B. Eliasson. Numerical modelling of the two-dimensional Fourier transformed Vlasov-Maxwell system. Journal of Computational Physics, 190(2):501–522, 2003.
- Engelmann et al.  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  EG Evstatiev and BA Shadwick. Variational formulation of particle algorithms for kinetic plasma simulations. Journal of Computational Physics, 2013.
- Filbet et al.  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.  L Galeotti, F Califano, and F Pegoraro. ”Echography” of Vlasov codes. Physics Letters A, 355(4):381–385, 2006.
- Gould et al.  R_W Gould, TM O’Neil, and JH Malmberg. Plasma wave echo. Physical Review Letters, 19(5):219–222, 1967.
- Grad  H. Grad. On the kinetic theory of rarefied gases. Communications on pure and applied mathematics, 2(4):331–407, 1949.
- Grant and Feix  F. C. Grant and M. R. Feix. Fourier-Hermite 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  W. Guo and J.-M. Qiu. Hybrid semi-Lagrangian finite element-finite difference methods for the Vlasov equation. Journal of Computational Physics, 2012.
- Haynes et al.  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.  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  P. Henrici. Discrete variable methods in ordinary differential equations. New York: Wiley, 1962, 1, 1962.
- Hilscher et al.  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  R.W. Hockney and J.W. Eastwood. Computer simulation using particles. CRC Press, 2010.
- Holloway  J. P. Holloway. Spectral velocity discretizations for the Vlasov-Maxwell equations. Transport Theory and Statistical Physics, 25:1–32, January 1996. doi: 10.1080/00411459608204828.
- Hou et al.  YW Hou, ZW Ma, and MY Yu. The plasma wave echo revisited. Physics of Plasmas, 18:012108, 2011.
- Imadera et al.  K. Imadera, Y. Kishimoto, D. Saito, J. Li, and T. Utsumi. A numerical method for solving the Vlasov-Poisson 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  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  A. J. Klimas. A Splitting Algorithm for Vlasov Simulation with Filamentation Filtration. Journal of Computational Physics, 110:150–163, 1994.
- Klimas  Alexander J Klimas. A numerical method based on the Fourier-Fourier transform approach for modeling 1-D electron plasma evolution. Journal of Computational Physics, 50(2):270–306, 1983.
- Landau  L. Landau. On the vibrations of the electronic plasma. J. Phys.(UssR), 10(1):25–34, 1946.
- Lapenta and Chacón  G. Lapenta and L. Chacón. Cost-effectiveness of fully implicit moving mesh adaptation: A practical investigation in 1D. Journal of computational physics, 219(1):86–103, 2006.
- Le Bourdiec et al.  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  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  H Ralph Lewis. Energy-Conserving Numerical Approximations for Vlasov Plasmas. Journal of Computational Physics, 6:136, 1970.
- Malmberg et al.  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.  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.  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(1-4):373–396, 2012.
- Mouhot and Villani  C. Mouhot and C. Villani. On landau damping. Acta mathematica, 207(1):29–201, 2011.
- Nevins et al.  WM Nevins, GW Hammett, AM Dimits, W Dorland, and DE Shumaker. Discrete particle noise in particle-in-cell simulations of plasma microturbulence. Physics of plasmas, 12(12):122305–122305, 2005.
- Ng et al.  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  J. Parker and P. Dellar. Hermite expansions with hypercollisionality for velocity space degrees of freedom in ion-temperature-gradient driven instabilities. Bulletin of the American Physical Society, 57, 2012.
- Pezzi et al.  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  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  J.-M. Qiu and A. Christlieb. A conservative high order semi-Lagrangian 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  J.-M. Qiu and C.-W. Shu. Positivity preserving semi-Lagrangian discontinuous Galerkin formulation: Theoretical analysis and application to the Vlasov-Poisson system. Journal of Computational Physics, 230:8386–8409, September 2011. doi: 10.1016/j.jcp.2011.07.018.
- Rossmanith and Seal  J. A. Rossmanith and D. C. Seal. A positivity-preserving high-order semi-Lagrangian discontinuous Galerkin scheme for the Vlasov-Poisson equations. Journal of Computational Physics, 230:6203–6232, July 2011. doi: 10.1016/j.jcp.2011.04.018.
- Schumer and Holloway  J.W. Schumer and J.P. Holloway. Vlasov simulations using velocity-scaled Hermite representations. Journal of Computational Physics, 144(2):626–661, 1998.
- Shoucri  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  M. Shoucri and G. Knorr. Numerical Integration of the Vlasov Equation. Journal of Computational Physics, 14:84, January 1974. doi: 10.1016/0021-9991(74)90006-0.
- Siminos et al.  E. Siminos, D. Bénisti, and L. Gremillet. Stability of nonlinear Vlasov-Poisson equilibria through spectral deformation and Fourier-Hermite expansion. Phys. Rev. E., 83(5):056402, May 2011. doi: 10.1103/PhysRevE.83.056402.
- Sonnendrücker et al.  E. Sonnendrücker, J. Roche, P. Bertrand, and A. Ghizzo. The Semi-Lagrangian 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  G. Strang. On the construction and comparison of difference schemes. SIAM Journal on Numerical Analysis, 5(3):506–517, 1968.
- Sydora  RD Sydora. Low-noise electromagnetic and relativistic particle-in-cell plasma simulation models. Journal of computational and applied mathematics, 109(1):243–259, 1999.
- Tang  T. Tang. The Hermite spectral method for Gaussian-type functions. SIAM Journal on Scientific Computing, 14(3):594–606, 1993.
- Tang  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.  T. Umeda, M. Ashour-Abdalla, and D. Schriver. Comparison of numerical interpolation schemes for one-dimensional electrostatic Vlasov code. Journal of Plasma Physics, 72:1057, December 2006. doi: 10.1017/S0022377806005228.
- Valentini et al.  F. Valentini, P. Veltri, F. Califano, and A. Mangeney. Cross-scale effects in solar-wind turbulence. Physical review letters, 101(2):025006, 2008.
- Verboncoeur  J. P. Verboncoeur. Particle simulation of plasmas: review and advances. Plasma Physics and Controlled Fusion, 47(5A):A231, 2005.
- Walker  E. Walker. The real cost of a CPU hour. Computer, 42(4):35–41, 2009.
- Wang et al.  Bei Wang, GH Miller, and Phillip Colella. A Particle-In-Cell method with adaptive phase-space remapping for kinetic plasmas. SIAM Journal on Scientific Computing, 33(6):3509–3537, 2011.
- Whealton et al.  JH Whealton, RW McGaffey, and PS Meszaros. A finite difference 3-D Poisson-Vlasov algorithm for ions extracted from a plasma. Journal of Computational Physics, 63(1):20–32, 1986.