Pith. sign in

REVIEW 3 major objections 5 minor 1 cited by

An explicit, energy-conserving particle-in-cell scheme

T0 review · 3 major / 5 minor · reviewed 2026-08-12 · deepseek-v4-flash

Pith's one-line read An explicit particle-in-cell scheme conserves total energy exactly by rescaling each particle's final velocity with a closed-form factor, without implicit field or particle solves.

desk verdict A genuinely new local correction for explicit PIC energy conservation, with a real gap in the accuracy proof for near-zero-velocity particles that a referee should push on. read the letter →

arxiv 2411.09605 v1 pith:NW3CFWHH submitted 2024-11-14 math.NA cs.NAphysics.plasm-ph

classification math.NAcs.NAphysics.plasm-ph MSC 65M7565M7035Q8365M06 PACS 52.65.Rr
keywords particle-in-cellenergyconservationVlasovequationexplicittimediscretizationLagrangemultiplierplasmasimulationPSATDgridheating
verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Particle-in-cell (PIC) plasma simulations suffer from grid heating: standard explicit time stepping lets total energy drift secularly, and the usual cure is an implicit field or particle solve. This paper constructs an explicit PIC time discretization that conserves total energy exactly, up to round-off, by adding one per-particle correction to the final velocity update. The correction is the closed-form solution of a small optimization problem; it is local, cheap, and does not change the scheme's second-order temporal accuracy. If correct, the scheme gives explicit codes the energy stability of implicit ones, which would remove the main cause of the finite-grid instability while keeping explicit scalability.

What carries the argument

The machinery is a per-particle Lagrange multiplier calculation that produces a closed-form rescaling factor $\Gamma_p$. Starting from the velocity predictor $v_p^\dagger$, the scheme solves $\min_v \|v-v_p^\dagger\|^2$ subject to the scalar constraint that the kinetic-energy increment matches the work done by the interpolated field; the solution is $\Gamma_p = \sqrt{\|v_p^n\|^2 + 2 v_p^*\cdot(v_p^\dagger-v_p^n)}/\|v_p^\dagger\|$. The factor is applied only at the end of each push, so the scheme stays explicit and the correction is parallelizable. The other load-bearing piece is the discrete integration-by-parts identity (19) for the spatial discretization, which connects the field-energy change to the particle sum; the paper verifies it for Yee and pseudospectral spatial operators.

What would settle it

Run the corrected electrostatic scheme on a strongly driven or poorly resolved case and check both the fraction of particles with negative radicand in (41) and the resulting energy error; if that fraction is not negligible, the energy error should exceed the claimed $O(10^{-12})$ level. More directly, one can search over initial conditions and $\Delta t$ for a case where $\|v_p^n\|^2 + 2v_p^*\cdot(v_p^\dagger-v_p^n)<0$ for a significant fraction of particles, which would show the local constraint is not always realizable.

Watch

Extended reading notes

Core claim

The central claim is that exact energy conservation is not lost when the energy-conserving Eulerian update is translated to particles, provided the final velocity is rescaled as $v_p^{n+1} = \Gamma_p v_p^\dagger$, where $v_p^\dagger$ is the standard second-order predictor and $\Gamma_p$ is the analytic solution (35) of a per-particle constraint. The constraint forces $v_p^*\cdot(v_p^{n+1}-v_p^n) = \tfrac12(\|v_p^{n+1}\|^2-\|v_p^n\|^2)$, which makes the particle kinetic-energy change exactly equal to the negative field-energy change provided the spatial discretization obeys the integration-by-parts identities (19). The authors prove $\Gamma_p=1+O(\Delta t^3)$ for the basic scheme and $1+O(\Delta t^4)$ for the corrected variant (44), so accuracy is preserved; the same correction works with the Amp\`ere-form electrostatic solve and with electromagnetic field solves based on Crank-Nicolson, leapfrog/Yee, or PSATD.

Load-bearing premise

Exact conservation depends on the square root in $\Gamma_p$ being a real number for every particle at every time step, since an imaginary value must be replaced by $\Gamma_p=1$, breaking conservation; the paper gives empirical evidence that imaginary values are extremely rare but no a priori bound.

Editorial extensions

If this is right

  • Existing explicit PIC codes can gain machine-precision energy conservation by appending a per-particle rescaling step and, for the corrected variant, one extra current deposition, with no implicit solve.
  • The secular energy growth known as grid heating should be eliminated for problems where the radicand stays positive, addressing the finite-grid instability that has limited explicit PIC for decades.
  • The scheme retains second-order accuracy because the correction factor differs from unity by $O(\Delta t^3)$ or $O(\Delta t^4)$, so conservation does not force a smaller time step.
  • Because the correction works with leapfrog/Yee and PSATD field solvers, existing electromagnetic PIC codes with good light-wave dispersion can adopt it without changing their field solve.
  • The magnetic field does zero work by construction, as in the Boris push, which protects long-time accuracy in magnetized problems.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

  • The per-particle constraint suggests that a CFL-like sufficient condition could be derived to guarantee the radicand is positive, turning the observed rarity of imaginary $\Gamma_p$ into a provable property for a restricted class of problems.
  • The same rescaling idea should extend to relativistic PIC by replacing the kinetic-energy norm with the relativistic kinetic energy; the paper lists relativity as future work.
  • Because the correction is local, it could combine with adaptive subcycling of individual particles or multi-rate field updates without a global synchronization point, which global Lagrange-multiplier schemes cannot offer.
  • A direct head-to-head cost-accuracy comparison against implicit PIC at the same energy error would clarify where the explicit scheme wins; the paper demonstrates conservation but does not quantify the runtime advantage.
Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Referee Report

3 major / 5 minor

Summary. The paper proposes an explicit, energy-conserving particle-in-cell scheme for the Vlasov equation. Starting from a two-stage explicit Eulerian discretization, the authors translate the method to PIC and show that the direct translation fails to conserve energy because the half-step velocity v* is not the exact average of v^n and v^{n+1}. They then enforce a scalar energy condition (32) through a per-particle Lagrange-multiplier correction, obtaining the closed-form velocity update v^{n+1} = Gamma_p v†_p with Gamma_p given by (35) and (41). The scheme is claimed to be second-order accurate because Gamma_p = 1 + O(Delta t^3) for the basic scheme (27)/(31) and Gamma_p = 1 + O(Delta t^4) for the corrected scheme (44). The construction is extended to the electromagnetic case with a Crank-Nicolson field solve, a leapfrog/Yee option, and a PSATD option, with energy-conservation proofs given for each. Charge conservation is discussed theoretically via sub-stepped shape functions. Numerical tests cover linear Landau damping, the two-stream instability, and the Weibel instability.

Significance. If the central claims hold, this is a valuable contribution: the correction is local and analytic, requires no fitted parameters, preserves the Boris-style zero-work property of the magnetic field, and is compatible with widely used spatial discretizations including Yee and PSATD. The energy-conservation proofs are clean and the numerical experiments are informative. The main weakness is that the exact-conservation and order-of-accuracy statements are made under an unstated regularity condition on the particle velocities, and the paper's own numerical results show that the practical scheme can deviate from exact conservation when the correction factor is imaginary. The core idea is sound, but the headline claims need to be qualified or strengthened before publication.

major comments (3)
  1. [Sec. 3.2, Eqs. (36) and (41)] The estimate Gamma_p = 1 + O(Delta t^3) is obtained by dividing by ||v†_p||^2, but no lower bound on ||v†_p|| is provided or established. For a particle with v^n_p = 0 in an electrostatic field, v†_p = Delta t E^{n+1/2}_p = O(Delta t), so the correction term in (36) becomes ((E^{n,*}_p - E^{n+1/2}_p) · E^{n+1/2}_p)/||E^{n+1/2}_p||^2, which is O(1), not O(Delta t^3). Moreover, the radicand in (41) is then (E^{n,*}_p · E^{n+1/2}_p)/||E^{n+1/2}_p||^2, which is negative whenever the electric field reverses between the stage and half-step; such reversals are generic in oscillatory or two-stream dynamics. Thus the realness of Gamma_p and the stated order of accuracy are not guaranteed for low-velocity particles without an additional hypothesis such as a uniform lower bound on ||v†_p||.
  2. [Sec. 3.3 and Appendix A, Eq. (44)] The improved estimate Gamma_p = 1 + O(Delta t^4) for the corrected scheme repeats the same division by ||v†_p||^2. The boundedness assumptions in Appendix A concern grid current and charge densities and their derivatives; they do not control the pointwise magnitude of v†_p for a single particle. Therefore the accuracy statement for ESEC2 is not uniform over phase space and can fail for particles whose velocity is O(Delta t) at the update, exactly the regime in which the radicand in (41) can become negative. The paper should either prove a uniform lower bound, restrict the theorem to a class of initial data and time steps where it holds, or modify the update so that the correction remains well defined and energy-preserving for all particles.
  3. [Abstract and Sec. 3.2, after Eq. (41)] The assertion of exact energy conservation is conditional on Gamma_p being real for every particle at every step. When Gamma_p is imaginary, the implemented fallback Gamma_p = 1 violates the constraint (32), and the energy proof (33) no longer applies. The paper provides only empirical evidence that imaginary values are rare; the two-stream test in Sec. 5.2 (Fig. 6) already shows a non-zero energy error for ESEC2 caused by exactly this effect. To support the headline claim, the paper needs either a quantitative a priori bound on the frequency and energy impact of imaginary radicands, or a fallback update that preserves energy even when (41) has no real root.
minor comments (5)
  1. [Sec. 5.1, Fig. 4] The quantity plotted as 'maximum fractional energy error' should be defined precisely; in particular, state whether it is max_n delta_n over the full simulation and how the time horizon affects the reported values.
  2. [Sec. 3.5] The charge-conservation argument is presented as a compatibility result and is not implemented in the numerical tests; the text should state this limitation more prominently so that readers do not infer that the tested schemes also conserve charge exactly.
  3. [Throughout] There are several typographical errors: 'electostatic' in Section 6, 'unforunate' and 'defintion' in Section 3.2, 'born out' should be 'borne out' in Section 3.3, and 'istropic' in reference [18].
  4. [Sec. 2.2, Eq. (12) and later] The notation introduces distinct shape functions S^h_rho, S^h_j, and S^h_E, but most subsequent formulas use a single S^h. The paper should state explicitly which shape functions are assumed identical in the main derivation and which are allowed to differ.
  5. [Sec. 3.2, after Eq. (37)] The Taylor expansion of the shape functions assumes that x*_p and x^{n+1/2}_p lie in the same cell. This is reasonable in the small-time-step limit but should be stated as a hypothesis of the accuracy analysis, and the relation to the cell-subdivision procedure used for charge conservation should be clarified.

Circularity Check

0 steps flagged · score 1.0 of 10

No significant circularity identified: energy conservation is enforced by construction, and the cited spatial-discretization identities receive independent proof in Appendix B.

full rationale

The paper's central energy-conservation claim is a direct consequence of its construction rather than a hidden reuse of its inputs. The scheme imposes the scalar constraint (32), v*_p·(v†_p - v^n_p) = (||v^{n+1}_p||^2 - ||v^n_p||^2)/2, and then derives the closed-form rescaling (35) as the analytic solution of the associated constrained minimization problem. The proof in (33) verifies that this constraint makes the particle kinetic-energy change telescope against the field energy change via the discrete integration-by-parts identities (19). No parameter is fitted to data, and no predicted quantity is an input in disguise. The only load-bearing external input is the discrete integration-by-parts identity for pseudospectral spatial discretizations, attributed to the authors' earlier work [23]; however, the paper supplies an independent proof in Appendix B using antisymmetry of the pseudospectral differentiation matrix, so this is not an unexplained self-citation. Charge conservation similarly relies on shape-function pairs established in previous works [8,15,23], which are cited as standalone results rather than derived from the present scheme. The admitted possibility of imaginary Gamma_p values in (41) is a robustness limitation, not a circular step, because the paper states it explicitly and monitors it empirically. Overall, the derivation is self-contained and the circularity score is minimal.

Assumptions & free parameters 0 free parameters · 6 assumptions · 0 invented entities

The method introduces no fitted parameters and no new physics. Its conservation properties rest on known discrete vector-calculus identities and on an empirical assumption about the rarity of imaginary correction factors.

assumptions (6)
  • domain assumption The spatial discretization satisfies the discrete integration-by-parts identities (19) for arbitrary grid quantities.
    These identities are the backbone of the energy conservation proof in (29); they are stated as assumptions in Section 2.2.3 and known to hold for Yee lattices and pseudospectral methods with binomial filtering.
  • domain assumption The electromagnetic spatial discretization satisfies the vector identity (50).
    Required for the EM energy conservation algebra in (49); proven for pseudospectral discretization in Appendix B and stated to hold for the Yee lattice.
  • domain assumption The same shape function S_h is used for current deposition and electric field interpolation.
    The equality in (29) between the particle-field coupling and the grid current requires S^j = S^E; the scheme definitions in (28) and (47) use a single S_h.
  • domain assumption For sufficiently small Δt, a particle's sub-step path lies within a single cell so the shape function is smooth on the segment connecting x^*_p and x^{n+1/2}_p.
    Used in the Taylor expansion analyses of ϵ2 in (39) and in Appendix A to derive the Γ accuracy; fails for particles crossing cell boundaries at finite Δt.
  • domain assumption Charge-conserving shape function pairs satisfying (54) exist and can be constructed as in prior work.
    The charge conservation claim in Section 3.5 relies on the existence of pairs (S^ρ, S^j) that satisfy the local quadratic condition; the paper cites [8,15,23] for explicit constructions and does not test them.
  • standard math The pseudospectral differentiation matrix is real and antisymmetric.
    Used in Appendix B to prove (50); cited from [30].

how reviews work

0 comments
Cite this review

Pith. "Pith review of An explicit, energy-conserving particle-in-cell scheme." pith.science (2026). https://pith.science/paper/NW3CFWHH

@misc{pith2026241109605,
  author       = {Pith},
  title        = {Pith review of: An explicit, energy-conserving particle-in-cell scheme},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/NW3CFWHH}},
  note         = {Machine review of arXiv:2411.09605}
}
abstract

We present an explicit temporal discretization of particle-in-cell schemes for the Vlasov equation that results in exact energy conservation when combined with an appropriate spatial discretization. The scheme is inspired by a simple, second-order explicit scheme that conserves energy exactly in the Eulerian context. We show that direct translation to particle-in-cell does not result in strict conservation, but derive a simple correction based on an analytically solvable optimization problem that recovers conservation. While this optimization problem is not guaranteed to have a real solution for every particle, we provide a correction that makes imaginary values extremely rare and still admits $\mathcal{O}(10^{-12})$ fractional errors in energy for practical simulation parameters. We present the scheme in both electrostatic -- where we use the Amp\`{e}re formulation -- and electromagnetic contexts. With an electromagnetic field solve, the field update is most naturally linearly implicit, but the more computationally intensive particle update remains fully explicit. We also show how the scheme can be extended to use the fully explicit leapfrog and pseudospectral analytic time-domain (PSATD) field solvers. The scheme is tested on standard kinetic plasma problems, confirming its conservation properties.

Figures

Figures reproduced from arXiv: 2411.09605 by the authors.

Figure 1
Figure 1. Evolution of electrostatic potential energy for linear Landau damping test case. [PITH_FULL_IMAGE:figures/full_fig_p029_1.png] view at source ↗
Figure 2
Figure 2. Fractional change in total energy over time for the three tested schemes. [PITH_FULL_IMAGE:figures/full_fig_p030_2.png] view at source ↗
Figure 3
Figure 3. Number of problematic particles with imaginary correction factors Γ [PITH_FULL_IMAGE:figures/full_fig_p031_3.png] view at source ↗
Figures from the paper (5 more)
Figure 4
Figure 4. Figure 4: Maximum fractional energy errors for the linear Landau damping test as a [PITH_FULL_IMAGE:figures/full_fig_p032_4.png]
Figure 5
Figure 5. Figure 5: Electrostatic potential energy for the two-stream instability test case, showing [PITH_FULL_IMAGE:figures/full_fig_p033_5.png]
Figure 6
Figure 6. Figure 6: Left: Fractional energy errors for the two-stream instability case as function of time for the three schemes tested. Right: Number of problematic particles at each time-step for ESEC1 and ESEC2. of existing explicit schemes with the energy conservation enjoyed by recen…
Figure 7
Figure 7. Figure 7: Left: Total potential energy as a function of time for both the pseudo-spectral analytic time-domain (PSATD) and Crank-Nicolson (CN) temporal discretizations of the field solve. Right: Breakdown of various sources of potential energy, again as a function of time. Appen…
Figure 8
Figure 8. Figure 8: Left: Fractional energy change as a function of time for both PSATD and CN schemes applied to the Weibel instability. Right: Problematic particles with imaginary Γ n p at each time-step in the Weibel instability problem. As with other examples, such particles are extre…

Discussion (0). Continue with ORCID to comment.

Forward citations

Cited by 1 Pith paper

Reviewed papers in the Pith corpus that reference this work. Sorted by Pith novelty score. Full citation record

  1. Intrinsic local Gauss's law preserving PIC method: A self-consistent field-particle update scheme for plasma simulations

    physics.plasm-ph 2025-06 conditional novelty 6.0 of 10

    A locally updated particle-in-cell scheme that maintains discrete Gauss's law exactly and avoids Poisson or Ampère solves is derived and tested on Landau damping, two-stream instability, and diocotron benchmarks.

Reference graph

Works this paper leans on

30 extracted references · 14 canonical work pages · cited by 1 Pith paper

  1. [1]

    Ku, C.-S

    S. Ku, C.-S. Chang, P. H. Diamond, Full-f gyrokinetic particle simula- tion of centrally heated global itg turbulence from magnetic axis to edge pedestal top in a realistic tokamak geometry, Nuclear Fusion 49 (11) (2009) 115021. 38

  2. [2]

    Fi´ uza, R

    F. Fi´ uza, R. A. Fonseca, L. O. Silva, J. Tonge, J. May, W. B. Mori, Three-dimensional simulations of laser–plasma interactions at ultrahigh intensities, IEEE Transactions on Plasma Science 39 (11) (2011) 2618– 2619

  3. [3]

    Y. Chen, S. E. Parker, Electromagnetic gyrokinetic δf particle-in-cell turbulence simulation with realistic equilibrium profiles and geometry, Journal of Computational Physics 220 (2) (2007) 839–855

  4. [4]

    J.-L. Vay, A. e. Huebl, A. Almgren, L. Amorim, J. Bell, L. Fedeli, L. Ge, K. Gott, D. Grote, M. Hogan, et al., Modeling of a chain of three plasma accelerator stages with the warpx electromagnetic pic code on gpus, Physics of Plasmas 28 (2) (2021)

  5. [5]

    C. K. Birdsall, A. B. Langdon, Plasma physics via computer simulation, CRC press, 2018

  6. [6]

    Friedman, A

    A. Friedman, A. Langdon, B. Cohen, Direct method for implicit particle-in-cell simulation, Comments Plasma Phys. Controlled Fu- sion;(United Kingdom) 6 (6) (1981)

  7. [7]

    A. B. Langdon, B. I. Cohen, A. Friedman, Direct implicit large time- step particle simulation of plasmas, Journal of Computational Physics 51 (1) (1983) 107–138

  8. [8]

    G. Chen, L. Chac´ on, D. C. Barnes, An energy-and charge-conserving, implicit, electrostatic particle-in-cell algorithm, Journal of Computa- tional Physics 230 (18) (2011) 7018–7036

Show all 30 references
  1. [9]

    G. Chen, L. Chacon, A multi-dimensional, energy-and charge- conserving, nonlinearly implicit, electromagnetic vlasov–darwin particle-in-cell algorithm, Computer Physics Communications 197 (2015) 73–87

  2. [10]

    Chac´ on, G

    L. Chac´ on, G. Chen, A curvilinear, fully implicit, conservative electro- magnetic pic algorithm in multiple dimensions, Journal of computa- tional physics 316 (2016) 578–597

  3. [11]

    D. C. Barnes, L. Chac´ on, Finite spatial-grid effects in energy-conserving particle-in-cell algorithms, Computer Physics Communications 258 (2021) 107560

  4. [12]

    Markidis, G

    S. Markidis, G. Lapenta, The energy conserving particle-in-cell method, Journal of Computational Physics 230 (18) (2011) 7037–7052. 39

  5. [13]

    Lapenta, Exactly energy conserving semi-implicit particle in cell formulation, Journal of Computational Physics 334 (2017) 349–366

    G. Lapenta, Exactly energy conserving semi-implicit particle in cell formulation, Journal of Computational Physics 334 (2017) 349–366

  6. [14]

    Bacchini, J

    F. Bacchini, J. Amaya, G. Lapenta, The relativistic implicit particle- in-cell method, in: Journal of Physics: Conference Series, Vol. 1225, IOP Publishing, 2019, p. 012011

  7. [15]

    G. Chen, L. Chacon, L. Yin, B. J. Albright, D. J. Stark, R. F. Bird, A semi-implicit, energy-and charge-conserving particle-in-cell algorithm for the relativistic vlasov-maxwell equations, Journal of Computational Physics 407 (2020) 109228

  8. [16]

    Gonoskov, Explicit energy-conserving modification of relativistic pic method, Journal of Computational Physics 502 (2024) 112820

    A. Gonoskov, Explicit energy-conserving modification of relativistic pic method, Journal of Computational Physics 502 (2024) 112820

  9. [17]

    L. Ji, Z. Yang, Z. Li, D. Wu, S. Jin, Z. Xu, An asymptotic-preserving and energy-conserving particle-in-cell method for vlasov–maxwell equa- tions, Journal of Mathematical Physics 64 (6) (2023)

  10. [18]

    Yee, Numerical solution of initial boundary value problems involving Maxwell’s equations in istropic media, IEEE Trans

    K. Yee, Numerical solution of initial boundary value problems involving Maxwell’s equations in istropic media, IEEE Trans. Antennas Propag. 14 (1966) 302–307

  11. [19]

    J.-L. Vay, I. Haber, B. B. Godfrey, A domain decomposition method for pseudo-spectral electromagnetic simulations of plasmas, Journal of Computational Physics 243 (2013) 260–268

  12. [20]

    Leh´ e, J.-L

    R. Leh´ e, J.-L. Vay, et al., Review of spectral maxwell solvers for electro- magnetic particle-in-cell: Algorithms and advantages, in: Proceedings of the 13th International Computational Accelerator Physics Confer- ence, Key West, FL, USA, 2018, pp. 20–24

  13. [21]

    H. Qin, S. Zhang, J. Xiao, J. Liu, Y. Sun, W. M. Tang, Why is boris algorithm so good?, Physics of Plasmas 20 (8) (2013)

  14. [22]

    Hairer, C

    E. Hairer, C. Lubich, Energy behaviour of the Boris method for charged- particle dynamics, BIT Numer. Math. 58 (2018) 969–979

  15. [23]

    L. F. Ricketson, G. Chen, A pseudospectral implicit particle-in-cell method with exact energy and charge conservation, Computer Physics Communications (2023) 108811

  16. [24]

    B. J. Sturdevant, L. Chac´ on, Eliminating finite-grid instabilities in gy- rokinetic particle-in-cell simulations, Journal of Computational Physics 464 (2022) 111330. 40

  17. [25]

    T. Z. Esirkepov, Exact charge conservation scheme for particle-in-cell simulation with an arbitrary form-factor, Computer Physics Commu- nications 135 (2) (2001) 144–153

  18. [26]

    L. F. Ricketson, A. J. Cerfon, Sparse grid techniques for particle-in-cell schemes, Plasma Physics and Controlled Fusion 59 (2) (2016) 024002

  19. [27]

    N. A. Krall, A. W. Trivelpiece, R. A. Gross, Principles of plasma physics, American Journal of Physics 41 (12) (1973) 1380–1381

  20. [28]

    T. H. Stix, Waves in plasmas, Springer Science & Business Media, 1992

  21. [29]

    Cheng, A

    Y. Cheng, A. J. Christlieb, X. Zhong, Energy-conserving discontinuous galerkin methods for the vlasov–maxwell system, Journal of Computa- tional Physics 279 (2014) 145–173

  22. [30]

    L. N. Trefethen, Spectral methods in MATLAB, SIAM, 2000. 41

Pith tools

Reviewed August 12, 2026 · model on record in the stance chip above.