Nyström Type Exponential Integrators for Strongly Magnetized Charged Particle Dynamics

Solving for charged particle motion in electromagnetic fields (i.e. the particle pushing problem) is a computationally intensive component of particle-in-cell (PIC) methods for plasma physics simulations. This task is especially challenging when the plasma is strongly magnetized due numerical stiffness arising from the wide range of time scales between highly oscillatory gyromotion and long term macroscopic behavior. A promising approach to solve these problems is by a class of methods known as exponential integrators that can solve linear problems exactly and are A-stable. This work extends the standard exponential integration framework to derive Nyström-type exponential integrators that integrates the Newtonian equations of motion as a second-order differential equation directly. In particular, we derive second-order and third-order Nyström-type exponential integrators for strongly magnetized particle pushing problems. Numerical experiments show that the Nyström-type exponential integators exhibit significant improvement in computation speed over the standard exponential integrators.

💡 Research Summary

**
This paper addresses the computational bottleneck of particle pushing in particle‑in‑cell (PIC) simulations of strongly magnetized plasmas. In such regimes, charged particles undergo rapid gyromotion around magnetic field lines while simultaneously experiencing much slower drift motions. Conventional pushers such as the Boris and Buneman schemes are second‑order leap‑frog methods that require very small time steps when the cyclotron frequency is large, leading to prohibitive cost.

Recent work has shown that exponential integrators, which treat the linear part of the dynamics exactly via matrix exponentials or related φ‑functions, can be A‑stable and allow larger time steps. Existing exponential propagation iterative Runge‑Kutta (EPIRK) methods, however, are formulated for first‑order systems, doubling the dimensionality of the Jacobian and increasing the cost of evaluating matrix functions.

The authors propose a novel class of Nyström‑type exponential integrators that work directly on the second‑order Newtonian form of the particle equation, (\ddot{x}=q/m,(E+\dot{x}\times B)). By integrating position and velocity together in a Nyström framework, the Jacobian size is halved, and the required matrix functions involve only the original 6‑by‑6 system (for a single particle) rather than a 12‑by‑12 system. Two specific schemes are derived:

1. Second‑order Nyström exponential integrator (EP2‑Ny) – uses the φ₁ function to advance both position and velocity in a single step.
2. Third‑order Nyström exponential integrator (EPRK3‑Ny) – introduces an intermediate stage and a nonlinear residual, employing φ₁ and φ₃ to achieve third‑order accuracy.

A critical component of any exponential integrator is the efficient evaluation of φ‑functions. The paper adopts the Lagrange‑Sylvester interpolation technique: eigenvalues of the Jacobian are computed, the φ‑function values at those eigenvalues are interpolated by a polynomial of degree N‑1 (with N the matrix size), and the polynomial is evaluated at the matrix itself. When eigenvalues are repeated, derivative (osculating) conditions are added to maintain a full rank system. This approach is shown to be substantially faster than Krylov subspace methods for the low‑dimensional particle problem.

Numerical experiments cover a suite of benchmark problems: pure uniform magnetic field gyromotion, E×B drift in combined uniform fields, grad‑B drift in non‑uniform magnetic fields, and more complex non‑uniform configurations. Across all tests, the Nyström‑type integrators achieve the same accuracy as the reference EPIRK2 and EPIRK3 schemes while requiring 2–3 times less CPU time. The third‑order Nyström method, in particular, remains stable and retains its high order even when the time step is large enough that the cyclotron period is under‑resolved, demonstrating robustness for long‑time simulations. Energy conservation and phase accuracy are comparable to, or better than, the traditional pushers.

The authors conclude that Nyström‑type exponential integrators exploit the intrinsic second‑order structure of charged‑particle dynamics, reduce the dimensionality of matrix operations, and retain A‑stability. They are ready for direct incorporation into existing PIC codes and promise significant speedups for simulations of strongly magnetized plasmas. Future work is outlined, including extension to multi‑particle systems with collective fields, handling strongly nonlinear magnetic effects, and optimizing the Lagrange‑Sylvester implementation on GPUs and other high‑performance architectures.

In summary, the paper delivers a mathematically sound, computationally efficient, and practically applicable advancement in particle pushing algorithms, offering higher-order accuracy and substantial performance gains for strongly magnetized plasma simulations.