The Energy Conserving Particle-in-Cell Method

Stefano Markidis; Giovanni Lapenta

Source-linked AI summary

The Energy Conserving Particle-in-Cell Method

Stefano Markidis, Giovanni Lapenta

arXiv:1108.1959v1Aug 9, 2011physics.comp-ph

TL;DR

PIC methods can introduce spurious energy during plasma-instability simulations, motivating exact finite-time-step energy conservation. The paper uses an implicit midpoint, concurrently solved JFNK formulation and reports correct instability behavior with exact total-energy conservation, while identifying charge and particle-motion constraints.

Problem
PIC schemes commonly fail to conserve total energy exactly, allowing numerical heating or cooling to affect plasma-instability simulations.
Method
The method implicitly midpoint-differences particle and Maxwell equations and solves their coupled nonlinear system with a Jacobian-free Newton Krylov solver.
Results
The method eliminates finite-grid instability, correctly describes two-stream and Weibel instabilities, and conserves total energy exactly in the reported simulations.
Takeaways & Limitations
Computational time increases linearly with particle number and is relatively insensitive to grid-point number and time step.
Takeaways & Limitations
The method does not conserve momentum and requires limiting the maximum particle motion to avoid aliasing instability.

Abstract

from arXiv · show

A new Particle-in-Cell (PIC) method, that conserves energy exactly, is presented. The particle equations of motion and the Maxwell's equations are differenced implicitly in time by the midpoint rule and solved concurrently by a Jacobian-free Newton Krylov (JFNK) solver. Several tests show that the finite grid instability is eliminated in energy conserving PIC simulations, and the method correctly describes the two-stream and Weibel instabilities, conserving exactly the total energy. The computational time of the energy conserving PIC method increases linearly with the number of particles, and it is rather insensitive to the number of grid points and time step. The kinetic enslavement technique can be effectively used to reduce the problem matrix size and the number of JFNK solver iterations.

1. Introduction

The paper addresses spurious energy behavior in PIC simulations by introducing a fully implicit electromagnetic scheme that conserves total energy exactly at finite time steps.

Motivation: Explicit PIC schemes typically heat numerically, while implicit schemes typically cool numerically, compromising total-energy conservation during plasma-instability studies.Spurious numerical energy can feed instabilities and produce unphysical results.
Prior work: Earlier energy-conserving PIC schemes improved conservation at finite time steps but did not achieve exact energy conservation.The Lewis scheme conserves energy exactly only in the zero-time-step limit.
Contribution: The proposed scheme conserves total energy exactly at finite time steps, with precision determined only by the iterative solver’s error tolerance.It uses midpoint integration for particle and field equations together with suitable discrete spatial operators.
Method: The method solves coupled nonlinear particle-field equations concurrently using a Jacobian-free Newton Krylov solver.The formulation targets electromagnetic simulations rather than only electrostatic cases.
Paper scope: The paper presents the algorithm, conservation and stability properties, implementation, simulations, and performance results of the energy-conserving PIC method.Its scope includes finite-grid, two-stream, and Weibel-instability evaluations described elsewhere in the paper.

2. Algorithm

The algorithm advances particles and electromagnetic fields implicitly and concurrently, using midpoint time differencing, Yee-lattice spatial operators, and JFNK solution of coupled nonlinear equations.

Particle model: Computational particles represent electrons and ions through positions and velocities whose evolution follows the particle equation of motion.Each computational particle is associated with a species charge-to-mass ratio and interpolated electromagnetic fields.
Particle model: Cloud-in-Cell interpolation transfers particle and grid quantities using a piecewise-linear weighting function that vanishes beyond one cell spacing.The stated interpolation weight is 1 − |x_g − x_p|/Δx when |x_g − x_p| < Δx, and zero otherwise.
Time discretization: The implicit midpoint rule time-differences the particle and Maxwell equations, with barred variables denoting time averages.The nonlinear formulation is rewritten in terms of average particle velocities and new electromagnetic fields.
Spatial discretization: Yee’s lattice and centered second-order spatial differences preserve the relevant discrete vector identities used by the Maxwell solver.The discrete operators satisfy discrete divergence-curl and product identities.
Computational cycle: At each cycle, a JFNK solver concurrently solves the coupled particle and field equations, after which particle positions and velocities are updated.The cycle repeats after initialization, which sets particle quantities and electromagnetic fields self-consistently.
Divergence equations: The nonconservative current-density definition does not satisfy charge continuity, so Gauss’ law can accumulate numerical error.A pseudo-current correction diffuses this error away through a heat-equation-like evolution controlled by d.
Divergence equations: In two-stream simulations, pseudo-current corrections with d = 0.1c^2/ω_pe and 0.5c^2/ω_pe reduce Gauss-law error without major changes to instability onset or growth rate.Without correction, the numerical error builds up slowly; the corrected cases reduce the error associated with charge nonconservation.

3. Energy Conservation

The proposed discretization conserves total energy exactly through midpoint time integration, current-density construction, and the resulting coupled particle–field equations. In an isolated system, the Poynting flux contributes zero to the total energy variation.

The discretized equations conserve the total energy exactly.The result follows from the numerical scheme presented for the coupled particle and field equations.
The field-energy variation contains particle work and a Poynting-flux term, whose contribution is zero in an isolated system.The particle-work term links the field-energy change to the particle dynamics.
Midpoint integration is used for particle and field equations to achieve energy conservation at finite time steps.
Energy conservation depends on the definition of the current density; charge-conserving PIC schemes using different current techniques do not necessarily conserve total energy.

4. Momentum Conservation

The method conserves energy but not momentum, because nonsmooth current deposition can produce spurious particle self-forces. This limitation can trigger phase-space instabilities and imposes a particle-motion constraint tied to grid spacing and time step.

The energy conserving PIC method does not conserve momentum because spurious particle self-forces arise from nonsmooth current deposition.The importance of these self-forces depends on interpolation functions, particles per cell, and grid spacing.
Two counter-streaming electron beams develop a phase-space instability associated with nonconservation of momentum.The simulation uses 10000 particles, velocities ±0.2c, 64 grid cells, and a time step of 0.5.
The aliasing instability disappears when the stated particle-motion condition is satisfied.
The modified Courant-Friedrichs condition restricts maximum particle motion to one and a half cells per time step.The constraint can still permit large time steps when grid spacing is chosen accordingly, and allows steps as large as those of implicit moment PIC under the stated comparison.

5. Numerical Stability

The numerical stability analysis finds real dispersion-relation roots for the energy-conserving scheme, indicating no exponential growth or damping for any time step in the analyzed cold-plasma case. Weibel simulations reproduce the instability while retaining radiation-field noise that numerical damping would remove.

The numerical stability analysis linearizes particle motion and derives a plasma dispersion relation from Fourier-transformed perturbations.The analysis assumes electric-field dependence proportional to exp(iωt).
For cold plasma, the dispersion-relation roots are always real, so the scheme is linearly unconditionally stable for any time step.The paper notes that grid-effect dispersion analysis was not carried out and was left for future work.
The implicit moment PIC comparison predicts growth for θ < 1/2, damping for θ > 1/2, and neither for θ = 1/2.For θ > 0.5, high-frequency unresolved waves are damped; θ = 0.5 corresponds to the energy-conserving scheme.
The Weibel instability appears in both simulations, while radiation-field noise remains at θ = 0.5 and is damped at θ = 0.6.The comparison is shown through the dispersion relation calculated from the magnetic-field component using Fourier transforms in space and time.

6. Implementation

The implementation uses a 1D3V energy-conserving PIC code that solves coupled particle and field equations implicitly with a Jacobian-free Newton–Krylov method. Particle updates, current deposition, and field equations are handled within the iterative solve.

Implementation: The implemented code is 1D3V with one-dimensional space and three-component particle velocities and electromagnetic fields.The paper states that extension to three dimensions is straightforward.
Implementation: At each computational cycle, a JFNK solver concurrently solves the implicitly midpoint-differenced particle and field equations.The equations are coupled through the current and are solved for average particle velocity and new electric and magnetic fields.
Newton–Krylov solver: The JFNK method solves G(x) = 0 through successive linear systems and updates the solution as x_i+1 = x_i + δx_i.The nonlinear iteration consists of solving a linear system followed by a solution update.
Newton–Krylov solver: GMRes provides the Krylov iterations, while the Jacobian is evaluated through a finite-difference Gateaux derivative rather than formed explicitly.This Jacobian-free construction avoids direct storage and calculation of the Jacobian matrix.
Implementation cycle: The residual evaluation sequentially estimates particle positions and current, interpolates fields, and computes the coupled particle-field residual.Particle equations, field equations, and current deposition are completed at each Krylov iteration.

7. Simulation Results

Tests in electrostatic and electromagnetic settings compare the energy-conserving PIC method with theory and explicit momentum-conserving PIC. The method preserves distributions and total energy while reproducing the tested instabilities.

Simulation scope: The tests cover finite-grid and two-stream instabilities electrostatically, followed by the Weibel instability electromagnetically.The stated goals are comparison with analytical theory and verification of exact energy conservation.
Solver tolerance: Smaller absolute and relative solver tolerances produce increased energy conservation in the two-stream simulation.The paper links tighter tolerances to improved energy conservation in the plotted histories.
7.1. Finite Grid Instability: 2% energy increase occurs in the explicit Maxwellian-plasma simulation, versus 10^-7% variation with energy-conserving PIC.The explicit simulation develops finite-grid-instability peaks, while the energy-conserving simulation retains the initial Maxwellian distribution.
7.2. Two-stream Instability: 0.35355 ω_pe is the predicted growth rate for the k = 1 ω_pe/c two-stream component, with simulation and linear theory in excellent agreement during the linear stage.Energy variation is approximately 5% for explicit PIC and 10^-4% for energy-conserving PIC.
7.3. Weibel Instability: 0.22ω_pe is the predicted growth rate for the k = 1ω_pe/c Weibel B_z component, with simulation and linear theory in good linear-regime agreement.The magnetic-field growth oscillates because of radiation-field noise; total energy varies within 10^-3%, while momentum oscillates between −0.02m_ec and 0.02m_ec.

8. Performance Results

Performance depends primarily on particle count, while grid size and time step have weaker effects. Kinetic enslavement reduces the large coupled problem associated with solving for every particle.

Computational scaling: A time step of 0.8ω_pe^-1 requires only 33% additional computational time compared with 0.1ω_pe^-1.The reported comparison comes from the Maxwellian-plasma performance study.
Computational scaling: Doubling the number of particles doubles computational time and reduces statistical noise, affecting Newton-step convergence.The performance study identifies particle count as the strongest computational dependence.
Solver tolerance: Decreasing JFNK error tolerance improves energy conservation but increases Newton iterations and computational time.The trade-off is illustrated through the solver-iteration and energy-conservation results.
Problem size: The method’s main disadvantage is a nonlinear system whose size grows with the number of computational particles.Typical simulations can involve millions of particles, producing matrices with ranks of order millions.
Kinetic enslavement: Kinetic enslavement restricts the problem unknowns to reduce the matrix size in large-particle simulations.The technique is evaluated by comparing Newton and average Krylov iterations with and without kinetic enslavement.

9. Conclusions

The energy conserving PIC method combines implicit coupled particle–field equations with JFNK solution to conserve total energy exactly. It reproduces key plasma-instability tests, with performance depending mainly on particle count and solver tolerance.

9. Conclusions: The method was tested against finite-grid, two-stream, and Weibel instabilities to assess correctness and exact total-energy conservation.The conclusions also report linear unconditional stability from implicit discretization.
9. Conclusions: The method is fully implicit and conserves total energy exactly through self-consistent particle and electromagnetic-field calculations in the JFNK solver.Its performance depends on the number of computational particles and solver error tolerance.
9. Conclusions: Computational time increases linearly with the number of particles and is relatively insensitive to the number of grid points and the time step.Lower solver error tolerance increases iterations and computational time.
9. Conclusions: Kinetic enslavement reduces the problem-matrix size and decreases the number of Newton iterations.The supplied figures compare Newton and average Krylov iterations with and without kinetic enslavement for the two-stream test.

Appendix A. ECpicES.m

Appendix A presents a skeleton electrostatic energy-conserving PIC implementation that initializes the plasma and field, advances the system with a Newton–Krylov solver, and records total energy.

Appendix A. ECpicES.m: The appendix provides a skeleton Matlab/Octave implementation for the electrostatic limit with electrons and motionless background ions.Additional solver files are required and are identified by name.
Appendix A. ECpicES.m: The code defines global simulation parameters, charges, initial particle positions and velocities, the initial electric field, and solver tolerances.The listed tolerance is tol = [1E−7, 1E−7].
Appendix A. ECpicES.m: The initial electric field is calculated using Gauss’ law and a Poisson solve before the computational cycle begins.The code also computes the initial field through a finite-difference expression involving the potential and grid spacing.
Appendix A. ECpicES.m: Each computational cycle invokes the Newton Krylov GMRes solver, updates particle positions and velocities, checks periodic boundaries, and computes and saves total energy.The average particle velocities and new electric field are obtained during the solver call.

Appendix B. residueEC.m

Appendix B defines the residual calculation used by the Newton–Krylov solver, combining particle average-velocity equations with electric-field equations and current-density evaluation.

Appendix B. residueEC.m: The residue function formulates the discretized equations of the energy-conserving PIC method for the Newton Krylov GMRes solver.The residual calculation includes particle and field components.
Appendix B. residueEC.m: Particle average-velocity equations and field equations are implemented as separate residual components.The appendix identifies the particle equations at line 22 and field equations at line 30.
Appendix B. residueEC.m: The residual calculation uses global domain, grid, time-step, particle, charge, background-density, and initial-state variables.These include L, dx, NG, DT, N, WP, QM, Q, rho back, x0, v0, and E0.
Appendix B. residueEC.m: The implementation evaluates particle positions at the n+1/2 time level, applies periodic-boundary checks, interpolates using neighboring grid indices, and calculates the average current J.The current is assembled from weighted particle contributions.