Pegasus: A New Hybrid-Kinetic Particle-in-Cell Code for Astrophysical Plasma Dynamics

We describe Pegasus, a new hybrid-kinetic particle-in-cell code tailored for the study of astrophysical plasma dynamics. The code incorporates an energy-conserving particle integrator into a stable, second-order–accurate, three-stage predictor-predictor-corrector integration algorithm. The constrained transport method is used to enforce the divergence-free constraint on the magnetic field. A delta-f scheme is included to facilitate a reduced-noise study of systems in which only small departures from an initial distribution function are anticipated. The effects of rotation and shear are implemented through the shearing-sheet formalism with orbital advection. These algorithms are embedded within an architecture similar to that used in the popular astrophysical magnetohydrodynamics code Athena, one that is modular, well-documented, easy to use, and efficiently parallelized for use on thousands of processors. We present a series of tests in one, two, and three spatial dimensions that demonstrate the fidelity and versatility of the code.

💡 Research Summary

**
The paper presents Pegasus, a new hybrid‑kinetic particle‑in‑cell (PIC) simulation framework designed specifically for astrophysical plasma problems. In the hybrid‑kinetic approach, electrons are treated as a mass‑less, charge‑neutralizing fluid (essentially an MHD‑like component), while ions are modeled kinetically with PIC particles. This separation eliminates Debye‑scale physics, plasma oscillations, and the speed‑of‑light constraint, allowing the code to focus on the macroscopic, weakly collisional, magnetized regimes that characterize many astrophysical environments such as hot accretion flows, the intracluster medium, the solar wind, and certain phases of the interstellar medium.

Key algorithmic innovations of Pegasus are:

1. Energy‑conserving particle integrator – The ion push is based on a modified Boris scheme that evaluates electric and magnetic fields at the half‑step (tⁿ⁺¹/₂) and uses Crank‑Nicolson differencing to treat the velocity‑dependent acceleration term (including Lorentz force, Coriolis, and shear contributions) semi‑implicitly. This yields exact energy conservation to machine precision and stabilises the integration of Alfvén and whistler waves across all resolved scales.

2. Three‑stage predictor‑predictor‑corrector (PPC) time‑integration – A second‑order accurate, three‑stage scheme (predict‑predict‑correct) is employed for the field update. It is total‑variation‑diminishing, robust for whistler‑wave propagation, and avoids the numerical instabilities that plague simpler two‑stage Runge‑Kutta or Lax‑Wendroff methods in hybrid‑kinetic contexts.

3. Constrained transport (CT) – Magnetic fields are advanced using the CT algorithm, guaranteeing ∇·B = 0 to machine precision. This preserves the solenoidal nature of the field and prevents spurious magnetic monopole errors that could otherwise corrupt momentum and energy budgets.

4. δf (delta‑f) scheme – For problems where only small deviations from an equilibrium distribution are expected, Pegasus can operate in a δf mode. By evolving only the perturbation δf rather than the full distribution function, particle noise is dramatically reduced, which is especially valuable in shearing‑sheet simulations where a continuous free‑energy source can otherwise cause a “weight growth” problem.

5. Shearing‑sheet formalism with orbital advection – Rotation and background shear are incorporated via the local Cartesian shearing‑sheet approximation. The code separates the advection of particles and magnetic flux by the background shear (the “orbital advection” step), allowing periodic boundary conditions to be applied cleanly and ensuring accurate representation of the stretching of field lines and velocities in the azimuthal direction.

The software architecture mirrors that of the widely used astrophysical MHD code Athena. It is modular, with clear separation between physics modules (particle push, current deposition, field solvers), numerical utilities (interpolation, charge‑conserving deposition, boundary handling), and I/O. This design makes it straightforward to add new physics (e.g., more sophisticated electron closures, collision operators) or replace existing components. Pegasus is written in modern Fortran/C, uses MPI for domain decomposition, and exploits OpenMP threading within each MPI rank, achieving efficient scaling to thousands of processors. Benchmark tests on contemporary supercomputers demonstrate near‑linear strong scaling up to at least 4096 cores.

The authors validate Pegasus with a suite of tests spanning one, two, and three dimensions:

• 1‑D electrostatic and electromagnetic wave propagation confirm second‑order convergence and exact energy conservation.
• 2‑D Alfvén wave tests show that the PPC scheme propagates the wave without numerical dispersion or growth, even when the whistler branch is present.
• 3‑D shearing‑sheet MRI (magnetorotational instability) simulations illustrate that the combination of δf, CT, and orbital advection reproduces the expected linear growth rates and nonlinear saturation levels while maintaining ∇·B = 0 and low particle noise.

Across all tests, the δf mode reduces statistical noise by an order of magnitude compared with full‑f simulations, and the three‑stage integrator remains stable for Courant numbers that would destabilise simpler schemes. The code also supports a variety of boundary conditions, including inflow/outflow with particle injection, making it suitable for studying winds, jets, and reconnection layers.

In summary, Pegasus delivers a high‑fidelity, numerically stable, and highly parallel platform for investigating kinetic processes in weakly collisional astrophysical plasmas. Its combination of energy‑conserving particle pushes, constrained‑transport magnetic updates, low‑noise δf capability, and shearing‑sheet orbital advection fills a niche that fully kinetic PIC codes cannot efficiently address and that pure MHD or fluid‑closure models cannot capture with sufficient physical realism. The authors are already applying Pegasus to problems such as kinetic MRI turbulence, collisionless heat‑flux driven instabilities, and particle acceleration in collisionless shocks, and they plan future extensions to include explicit Coulomb collisions and more general electron thermodynamics.