A Particle-in-cell Method for Plasmas with a Generalized Momentum Formulation, Part I: Model Formulation

TL;DR

This work introduces a particle-in-cell method for the relativistic Vlasov-Maxwell system based on a potential formulation under the Lorenz gauge, where particles evolve via a generalized-m momentum Hamiltonian. The field solver employs a MOL^T-inspired, unconditionally stable wave solver with a Green’s-function based integral approach, and spatial derivatives are computed analytically to match the field convergence rate. The resulting non-separable Hamiltonian PIC method (IAEM) demonstrates mesh-independent numerical heating, robustness at Debye-length under-resolution, and good agreement with standard PIC benchmarks across 1D and 2D problems, including plasma sheaths and relativistic beams. These results suggest practical efficiency gains and stability advantages for simulations in complex geometries, with future directions including enforcing the Lorenz gauge, higher-order time accuracy, and fully implicit formulations.

Abstract

This paper formulates a new particle-in-cell method for the Vlasov-Maxwell system. Under the Lorenz gauge condition, Maxwell's equations for the electromagnetic fields can be written as a collection of scalar and vector wave equations. The use of potentials for the fields motivates the adoption of a Hamiltonian formulation for particles that employs the generalized momentum. The resulting updates for particles require only knowledge of the fields and their spatial derivatives. An analytical method for constructing these spatial derivatives is presented that exploits the underlying integral solution used in the field solver for the wave equations. Moreover, these derivatives are shown to converge at the same rate as the fields in the both time and space. The field solver we consider in this work is first-order accurate in time and fifth-order accurate in space and belongs to a larger class of methods which are unconditionally stable, can address geometry, and leverage fast summation methods for efficiency. We demonstrate the method on several well-established benchmark problems, and the efficacy of the proposed formulation is demonstrated through a comparison with standard methods presented in the literature. The new method shows mesh-independent numerical heating properties even in cases where the plasma Debye length is close to the grid spacing. The use of high-order spatial approximations in the new method means that fewer grid points are required in order to achieve a fixed accuracy. Our results also suggest that the new method can be used with fewer simulation particles per cell compared to standard explicit methods, which permits further computational savings.

Figures (24)

• Figure 1: Space-time refinement of the solution and its spatial derivatives for the two-dimensional periodic example \ref{['subsubsec:field-solver periodic results']} obtained with the first-order BDF method.
• Figure 2: Space-time refinement of the solution and its spatial derivatives for the two-dimensional Dirichlet problem \ref{['subsubsec:field-solver dirichlet results']} obtained with the first-order BDF method.
• Figure 3: Trajectories for the single particle test obtained using the Boris method Boris1970, the AEM Gibbon2017Hamiltonian, and the IAEM. The particle rotates about a static magnetic field which points in the $z$-direction. Also shown is the time history of the Hamiltonian generated by each of the methods which is measured relative to the initial data. In particular, the AEM shows a growth in the overall energy causing the gyroradius to increase. This behavior is not observed in the improved method.
• Figure 4: Refinement study for the trajectory of a single particle obtained using the Boris method Boris1970, AEM Gibbon2017Hamiltonian, and the IAEM that uses a Taylor correction. Errors are measured in the $\ell_{\infty}$-norm against a reference solution obtained using the Boris method. Even though the AEM with the Taylor correction remains globally first-order accurate in time, its improvement over the AEM is quite apparent (roughly an order of magnitude).
• Figure 5: We plot the electrons in phase space obtained with the Poisson model for the two-stream instability problem at different times given in units of $\omega_{pe}^{-1}$. Results obtained using leapfrog time integration are shown in the left column, while the right column applies the AEM. The IAEM, which applies the Taylor correction, is not considered here because the contributions from the magnetic fields are ignored. The FFT is used to compute the scalar potential (as well as its gradient). Identical results are observed with both approaches.