A Vlasov Algorithm Derived from Phase Space Conservation
AA Vlasov Algorithm Derived from Phase Space Conservation
Jonathan P. Edelen and Stephen D. Webb a) RadiaSoft, LLC (Dated: 27 July 2020)
Existing approaches to solving the Vlasov equation treat the system as a partial dif-ferential equation on a phase space grid, and track in either an Eulerian, Lagrangian,or semi-Lagrangian picture. We present an alternative approach, which treats theVlasov equation as a conservative flow on phase space, and derives its equations ofmotion using particle-pushing algorithms akin to particle-in-cell methods. Depositionto the grid is determined from the convolution of local basis functions. This approachhas the benefit of allowing flexible definitions in the grid, which are decoupled fromhow the phase space flow evolves. We present numerical examples and comment onthe various properties of the algorithm. a) a r X i v : . [ phy s i c s . c o m p - ph ] J u l . INTRODUCTION The core problem of plasma physics is to solve for the evolution of the phase space densityunder external and self-consistent forces. Where analytic techniques, such as linearizedVlasov, fail, we are required to resort to numerical methods. The two dominant methods forcomputing the evolution of self-consistent plasma dynamics are the particle-in-cell approach– which follows the trajectories of macro-particles that sample the phase space – and thesolution of the Vlasov equation on a numerical grid – which advects the phase space densityon a mesh using various methods.Particle-in-cell methods trace macro-particles that sample the phase space distribu-tion as they move in self-consistent and external fields. This approach is quite mature,and has been used to solve a variety of self-consistent many-body problems ranging fromsemiconductor transport to large-scale astrophysics. This approach has the benefit of beingsubstantially faster than Vlasov solvers, but are artificially noisy due to the down-samplingof the number of particles. The noise is therefore amplified by a factor ∼ (cid:113) N real N macro. Vlasov solvers treat phase space density as a fluid, and advect the fluid based on theVlasov equation. These approaches treat the evolution of phase space by solving the Vlasovpartial differential equation. Early approaches, for example, used a local Chebyshev polyno-mial basis , or using splitting and spline or Fourier interpolation . Finite element approacheshave also been used . Modern approaches, such as flux balancing , use conservation ofphase space volume and local numerical integrals to implicitly back-track the characteristictrajectories of a local phase space volume. These characteristic methods rely on interpo-lation methods to compute the characteristics, which can be computationally expensive.Although the Vlasov equation is a partial differential equation describing the evolution of aphase space fluid, the evolution of the phase space density is taken not from the fluid equa-tion, but by tracing the characteristics of each volume of phase space, conceptually closerto particle-in-cell methods.We present here a hybrid method, which evolves the phase space density on a meshby tracing trajectories of those mesh points using particle-in-cell-like particle pushes. Westart by reviewing the formal solution of the Vlasov equation in terms of flows of pointsin phase space along time. We then translate the flow and phase space conservation intoa discrete mesh and make contact with particle-in-cell particle pushers, such as the Boris2usher or symplectic integrator methods. We conclude by demonstrating the method on twoexamples – trapped and streaming particles in the pendulum equation, and the high-gainone-dimensional free-electron laser model . II. THE VLASOV EQUATION AS FLOWS ON PHASE SPACE
The collisionless Vlasov equation can be formulated as a fluid equation on a 2 × D -dimensional phase space. The statement of the Vlasov equation as a partial differentialequation is ∂ t f + z (cid:48) · ∇ z f = 0 (1)where z = ( x , v ) are the configuration space coordinates. For Hamiltonian systems withHamiltonian H = p / V ( x ), for example, this becomes ∂ t f + p∂ x f + ( ∂ x V ) ∂ p f = 0 . (2)If z (cid:48) = v ( z , t ), then the solution for the individual trajectories can be given by a flow φ ,such that z ( t ) = φ t → t ◦ z ( t ), and the Vlasov equation can be solved by the method ofcharacteristics. Thus, the Vlasov equation is equivalent to the conservation of phase spaceover a flow φ that maps an initial coordinate z in. to z fin. related by z fin. = φ ◦ z in. . In thiscontext, neglecting collisions, the Vlasov equation is equivalent to the statement that f ( z in. ) = f ( z fin. ) (3)for the phase space density f .Particle pushing algorithms approximate the flow for a short time, the time step h ,using a discretization of the equations of motion. This could be by means of Runge-Kuttamethods, the Boris or Boris-related particle pushers, or with symplectic integrators forsingle-particle or self-consistent Hamiltonian systems . It is thus well-known how tocompute z fin. from z in. given the fields the particles move in. We therefore need only todetermine the final distribution in eqn. (3) in some computable representation.3 II. AN ALGORITHM BASED ON PHASE SPACE DENSITYCONSERVATION
The central concept of this algorithm is to take the relationship in eqn. (3) and develop adiscrete representation using local basis functions. There are many options for representingthe phase space density function, from Fourier series to splines. The derivation that followsapplies to an arbitrary choice of basis functions, although we will in the examples specializeto localized tent functions instead for the sake of concreteness.We discretize the phase space density as a set of overlapping basis functions: f ( x, p ) = (cid:88) j w j Λ( x − x j )Λ( p − p j ) . (4)from eqn. (3) f ( x in. , p in. ) = f ( x fin. , p fin. ) (5)Under this scheme, the individual ( x j , p j ) coordinates change by ( δx j , δp j ) in a single timestep. We can then write eqn. (3) in terms of the local basis functions (cid:88) j w ( fin. ) j Λ( x − x j )Λ( p − p j ) = (cid:88) k w ( in. ) k Λ( x − x k − δx k )Λ( p − p k − δp k ) . (6)By premultiplying by Λ( x − x j (cid:48) )Λ( p − p j (cid:48) ), we can reduce this to a matrix equation for apotentially overcomplete basis set: (cid:88) j Λ jj (cid:48) w ( fin. ) j = (cid:88) k ∆ kj (cid:48) w ( in. ) k (7)where Λ jj (cid:48) = (cid:90) dxdp Λ( x − x j )Λ( p − p j )Λ( x − x j (cid:48) )Λ( p − p j (cid:48) ) (8)and ∆ kj (cid:48) = (cid:90) dxdp Λ( x − x k − δx k )Λ( p − p k − δp k )Λ( x − x j (cid:48) )Λ( p − p j (cid:48) ) . (9)Note that if we choose an orthogonal basis set, then Λ will be diagonal. However, we arenot required to make this choice. The update sequence for each w ( fin. ) j is then given by thematrix equation w ( fin. ) j = (cid:88) j (cid:48) ,k Λ − jj (cid:48) ∆ kj (cid:48) w ( in. ) k (10)4hat depends on the matrix inverse of Λ – which, will also be diagonal if we use an orthogonalbasis – and ∆ , which contains the characteristic trajectory information over the time step.Global phase space volume conservation comes from the extent to which (cid:88) j w ( fin. ) j = ? (cid:88) j w ( in. ) j . (11)The global conservation can therefore be written with the condition (cid:88) j (cid:32)(cid:88) k Γ j,k w ( in. ) k (cid:33) = (cid:88) j w ( in. ) j (12)where Γ j,k = (cid:80) j (cid:48) Λ − jj (cid:48) ∆ kj (cid:48) . This is not satisfied for arbitrary choice of local basis functions,so the algorithm will not in general conserve the total phase space volume. The matrix Γ also has implications for positivity, and indeed there is no guarantee that the weights willremain positive for an arbitrary choice of basis. We will see an instance in our exampleswhere negativity in the weights appears at a discontinuity in the phase space density.The primary computational demand is in computing the matrix (cid:80) j (cid:48) Λ − jj (cid:48) ∆ kj (cid:48) , as both are N × N -dimensional matrices for N unique phase space volume cells. ∆ must be recomputedevery time step, while Λ − is a constant throughout the simulation. The convolving integralsthat define Λ and ∆ integrate over smooth functions the local phase space density, so thatfilamentation on scales smaller than the grid resolution are smoothed across the grid and donot present a computational challenge.Because this approach is built on established particle pushers, it can be extended to ar-bitrary dimensionality. The only dimension-dependent computations are the Λ and ∆ ma-trices, and those require integrating local regular basis shapes – these integrals are straight-forward to compute analytically. Furthermore, because the approach requires only a high-dimensional grid and is otherwise built on particle-in-cell techniques, the approach could bewell-suited as an addition to mature particle-in-cell codes. IV. NUMERICAL EXAMPLES
We demonstrate this algorithm on two example problems – the harmonic oscillator andpendulum problems, and the one-dimensional FEL equations. In each simulation we assumeperiodic boundaries in the coordinates and open boundaries in the momenta. As a metricfor performance and convergence, we will look at global phase space volume conservation5
IG. 1. Phase space evolution of a linear harmonic oscillator. and entropy. We expect to see some loss of phase space volume, as volume elements thatpass outside the range of the momentum grid will be lost.We do not make any comments on the timing performance of the algorithm, as theimplementation is more for proof-of-concept than to benchmark performance or, for example,scalability on many-core architectures. We do not that the basis matrix computations can bequite expensive – Λ − and ∆ are both M N × M N matrices, ∆ is a sparse matrix that mustbe recalculated every step, and Λ − is the inverse of a sparse matrix, which therefore maynot be itself sparse although it may be pre-computed once at the beginning of the simulation,although storing 8 × ( M N ) double-precision numbers in memory may be prohibitive forlarge grids. A. Harmonic Oscillators
In order to test that the algorithm is properly handling the phase-space deposition wetested both a linear harmonic oscillator and a nonlinear harmonic oscillator. For the linearoscillator we analytically advanced the phase space distribution to understand any inherentdiffusion in the algorithm. Figure 1 shows the evolution of a gaussian initial phase spacethrough one revolution around the origin. In this case the phase advance was π/ H = p / x ).Figure 2 shows the initial phase space and the phase space at four snapshots along thesimulation. The integration step was 2 π/
20 and the snapshots were taken at steps 10, 20,30, and 80 going from left to right. 6
IG. 2. Phase space evolution of the pendulum equation with a second order symplectic integrator.FIG. 3. Fractional change in the sum of the weights as a function of time for different gridresolutions for both the linear and nonlinear oscillator.
The colormap for this figure is red for positive values and blue for negative values. In thependulum example, negative weights arise and propagate from regions where filamentationis quite strong. We can see this effect occurring next to the hyperbolic fixed points in theplot of step 30 in fig. (2), which will have a sharp discontinuity in the phase space densityas any initial density beginning there will not move, while all the neighboring density valueswill tend to dilute.In these regions, the grid is under-resolving a phase space feature, creating local oscillatoryvalues in the phase space density. Figure 3 shows fractional change in the sum of the weightsas a function of time for both the linear and nonlinear oscillator. Remarkably, the total volume of phase space is well-conserved, even with the introduction of negative weights.The oscillatory behavior is similar in character to the Gibbs phenomenon in the Fourieranalysis of discontinuous functions, but which can also be seen in, for example, electrostaticPoisson solvers on distributions with hard edges.7
IG. 4. Fractional change in the system entropy as a function of time for different grid resolutionsfor both the linear and nonlinear oscillator.
These discontinuities in phase space at the fixed point can be well described by an in-creasing entropy of the system. Figure 4 shows the fractional change in entropy for boththe linear and nonlinear oscillators. For the linear oscillator as expected there is no changein the system’s entropy while for the nonlinear oscillator we are able to model the increasein entropy due to these hard edges in the phase space evolution.
B. FEL simulations
Next we simulated a FEL by numerically integrating the 1-D normalized FEL equationsusing a fourth-order Runge-Kutta algorithm. d Φ dz = δ + η ( z ) (13a) dηdz = − (cid:2)(cid:0) A ( z ) + iσ (cid:104) e − i Φ( z ) (cid:105) (cid:1) e i Φ( z ) + c.c. (cid:3) (13b) dAdz = (cid:104) e − i Φ( z ) (cid:105) (13c)Here Φ is the ponderomotive phase which is the position coordinate in phase space, η is the normalized energy, A is the self-consistent field magnitude, and z is our independentvariable. Using these equations and our Vlasov code we simulated a FEL starting up fromnumerical noise. To examine the effect of hard edge distributions and negative weighs wesimulated to distributions. One with a hard edge in energy, Figure 5, and one with a8 IG. 5. Phase space evolution for a self-consistent FEL, flat top initial distribution.FIG. 6. Phase space evolution for a self-consistent FEL, gaussian initial distribution.
Gaussian distribution in energy 6. Figures 5 and 6 show the initial phase-space distributionand the phase-space after 50, 100, and 150 steps. After 50 steps we can see the beginnings ofphase-space bunching which becomes more prominent after 100 steps. after 150 steps thereis clear filamentation in the phase-space caused by saturation of the FEL process.We can also see the same oscillatory behavior clearly in the evolution of a distributionthat has a hard edge in fig. (5), where the distribution is initially uniform for momentumbetween ± . IG. 7. Fractional change in the phase-space density as a function of time step for the two differentFEL simulations as we increase the mesh density.FIG. 8. Fractional change in the entropy as a function of time step for the two different FELsimulations as we increase the mesh density. boundaries in energy. However for the flat-top distribution we see a conserved phase-spacedensity while the entropy oscillates. This is attributed to the negative weights introducedby evolving a top-hat distribution on a finite grid.This Gibbs-like phenomenon resulting by depositing a hard edged distribution onto afinite grid is further demonstrated by examining a line out of the weights as a function ofenergy. Figure 9 shows the weights as a function of energy for a fixed position for threedifferent time steps. The negative weights are present shortly after the distribution startsevolving for the hard edged distribution while they only become apparent for the Gaussiandistribution when strong filamentation occurs in the simulations.10
IG. 9. Line-out plot of the weights from the FEL simulations at three different time-steps. Topis the for the flat-top distribution in energy and bottom is for the Gaussian distribution in energy.
V. DISCUSSION
We have presented a new approach to deriving algorithms for solving the Vlasov equationon a gridded representation of phase space, based on tracking the flows of the phase spacegrid points and projecting the shifted local basis functions back onto the original basis. Thisapproach uses particle tracking algorithms from particle-in-cell or other techniques to solvewhat is conventionally treated as a fluid equation.This approach is explicit in time, with the initial and final weights of a time step relatedby a matrix equation. Those matrices are determined by convolution integrals with basisfunctions for the phase space density, which are completely general in principle. This flexi-bility could allow for clever choices in basis functions for specific applications that captureessential physics intrinsically, allowing for a Vlasov algorithm that is particularly efficientfor certain types of problems.We demonstrated this algorithm in a handful of test cases – the linear oscillator andpendulum equations where the trajectories are single-particle, and for the one-dimensionalFEL equations, a plasma instability. For our implementation, we chose to use local tentfunctions. This specific implementation showed oscillatory behavior and negative-valued11eights emerging from under-resolved phase space features, such as discontinuities in thelocal density or phase space filaments that are on the order of a grid size. This feature wasnot present when the distribution is smooth. It is possible that smoother basis functionswill ameliorate this effect.A potential extension of this scheme is to provide a formal way to bridge the gap betweenmoment-based fluid equations and full Vlasov equations. For example, if the phase spacedensity is assumed to be locally thermal, so that the distribution is given by f ( v, x ) = N (cid:112) πσ ( x ) exp (cid:40) − ( v − v ( x )) σ ( x ) (cid:41) (14)then we can formally relate changes in σ to changes in v by carefully selecting local basisfunctions f = (cid:88) j w j (cid:112) πσ ( x j ) exp (cid:40) − ( v − v ( x j )) σ ( x j ) (cid:41) . (15)This may provide a formal method to stitch together various high-order fluid models acrossa simulation domain based on the local physics. VI. ACKNOWLEDGEMENTS
This work was supported in part by the Department of Energy Office of Science, Officeof Basic Energy Science under contract no. de-sc0017161.
REFERENCES C. K. Birdsall and A. B. Langdon,
Plasma Physics via Computer Simulation (McGraw-Hill, New York, 1985). R. W. Hockney and J. W. Eastwood,
Computer Simulation Using Particles (Taylor &Francis, 1989). M. Shoucri and G. Knorr, “Numerical integration of the Vlasov equation,” J. Comp. Phys. , 84–92 (1974). C. Z. Cheng and G. Knorr, “Integration of the Vlasov equation in configuration space,”J. Comp. Phys. (1976). S. Zaki, L. Gardner, and T. Boyd, “A finite element code for the simulation of one-dimensional Vlasov plasmas I. Theory,” J. Comp. Phys. , 184–199 (1988).12 S. Zaki, L. Gardner, and T. Boyd, “A finite element code for the simulation of one-dimensional Vlasov plasmas II. Applications,” J. Comp. Phys. , 200–208 (1988). F. Filbet, E. Sonnendr¨ucker, and P. Bertrand, “Conservative Numerical Schemes for theVlasov Equation,” J. Comp. Phys. , 166–187 (2001). M. Shoucri, “Eulerian codes for the numerical solution of the Vlasov equation,” Commun.Nonlinear Sci. Numer. Simul. , 174–182 (2008). R. Bonifacio, C. Pellegrini, and L. Narducci, “Collective instabilities and high-gain regimein free electron laser,” Opt. Commun. (1984). J. P. Boris, “Relativistic plasma simulation-optimization of a hybrid code,” in
Proceedingof the Fourth Conference on Numerical Simulations of Plasmas (1970). J.-L. Vay, “Simulation of beams or plasmas crossing at relativistic velocity,” Phys. Plasmas (2008). E. Chacon-Golcer and F. Neri, “A symplectic integrator with arbitrary vector and scalarpotentials,” Phys. Lett. A , 4661–4666 (2008). E. Forest, “Geometric integration for particle accelerators,” J. Phys. A: Math. Gen. ,5321–5377 (2006). H. Qin and X. Guan, “Variational Symplectic Integrator for Long-Time Simulations ofthe Guiding-Center Motion of Charged Particles in General Magnetic Fields,” Phys. Rev.Lett. (2008). R. D. Ruth, “A canonical integration technique,” IEEE Trans. Nucl. Sci.
NS-30 (1983). D. T. Abell, N. M. Cook, and S. D. Webb, “Symplectic modeling of beam loadingin electromagnetic cavities,” Phys. Rev. Acc. Beams (2017), 10.1103/PhysRevAccel-Beams.20.052002. E. Evstatiev and B. A. Shadwick, “Variational formulation of particle algorithms for kineticplasma simulations,” J. Comp. Phys. , 376–398 (2013). H. R. Lewis, A. Sykes, and J. A. Wesson, “A Comparison of Some Particle-in-Cell PlasmaSimulation Methods,” J. Comp. Phys. , 85–106 (1972). H. Ralph Lewis, “Energy-conserving numerical approximations for Vlasov plasmas,” J.Comp. Phys. , 136–141 (1970). H. Ralph Lewis, “Application of Hamilton’s Principle to the Numerical Analysis of VlasovPlasmas,” Methods. Comput. Phys. , 85–106 (1970).13 A. Stamm, B. Shadwick, and E. Evstatiev, “Variational Formulation of Macro-ParticleModels for Electromagnetic Plasma Simulations,” IEEE Trans. Plasma Science (2014). S. D. Webb, “A Spectral Canonical Electrostatic Algorithm,” Plasma Phys. ControlledFusion58