Structure-preserving collision integrator to complement the symplectic orbit push

[krystophny]

Context

The orbit push is symplectic (Albert, Kasilov & Kernbichler 2020), so the collisionless guiding-center invariants stay bounded over long integration. Collisions are optional (swcoll) and live in src/collis_alphas.f90 as a test-particle Coulomb operator against a fixed Maxwellian background:

• coleff/onseff give the Fokker-Planck coefficients: momentum-module diffusion dpp, pitch-angle diffusion dhh, and drag fpeff.
• stost applies one explicit Euler-Maruyama step. Pitch and momentum each take a diffusion kick sqrt(2 D dtau) * xi plus a drift * dtau, with reflecting boundaries at |lambda| = 1 and pmin.
• collide is operator-split from the orbit push, once per step (simple_main.f90).

So the reversible half is structure-preserving and the irreversible half is not. Euler-Maruyama is strong order 1/2, does not vanish on the Maxwellian exactly, and the boundary reflection at |lambda| = 1 distorts the pitch distribution. Over long slowing-down times this leaves a numerical drift in the test-particle energy and pitch statistics, which is the input the loss fraction is read from.

Proposal

Make the collision step structure-preserving as well: the discrete operator should vanish on the Maxwellian and relax to it at the analytic rate, with a discrete H-theorem (detailed balance). That is the metriplectic-consistency condition at the ensemble level, the dissipative counterpart of the symplectic orbit push. Two levels, by scope:

1. Drop-in, fixed background (matches SIMPLE today). Replace the linearized pitch kick with an exact norm-preserving rotation of the pitch variable (Boozer & Kuo-Petravic 1981), so |v| is conserved during pitch scattering and no boundary reflection is needed. Pair it with an energy step that preserves the Maxwellian moments (moment-preserving Monte-Carlo, e.g. the 2024 schemes below). Same cost class, better long-time statistics.

2. Self-consistent / multi-species (larger change). Once the background is no longer fixed, use the metriplectic particle discretization of the Landau operator (Kraus & Hirvijoki 2017 and successors): a discrete-gradient step that conserves mass, momentum, and energy and dissipates entropy monotonically. This realizes the full GENERIC structure (Poisson bracket plus metric bracket) shared with the symplectic push, via the reversible-irreversible split (Shang & Öttinger 2020).

What is realizable now: a single Langevin trajectory against a fixed thermostat has no closed energy to conserve, so the target for SIMPLE as it stands is option 1, preserving the equilibrium and its moments, not a per-trajectory energy-conserving step. Option 2 is the right object once collisions become self-consistent.

Plan

• Add exact pitch-angle scattering (norm-preserving rotation) behind a switch, alongside the current stost. Verify against the analytic relaxation of <lambda> and <lambda^2>.
• Add a Maxwellian-preserving energy step. Check that an initially Maxwellian test population stays Maxwellian to sampling noise, and a shifted population relaxes at the analytic rate.
• Add a discrete H-functional diagnostic. Confirm monotone relaxation and a collision step that vanishes on the Maxwellian.
• Benchmark loss fraction and slowing-down spectra against the current operator on a reference VMEC case. Quantify the statistical-fidelity gain at fixed marker count and step.
• (Stretch) Prototype the metriplectic Landau step for the self-consistent multi-species path.

Cost

The collision step is a small fraction of runtime; the field and spline evaluations in the orbit push dominate. The explicit Maxwellian-preserving variant (option 1) costs about the same, a few extra trig/exp per step, so the overhead is in the noise. An implicit discrete-gradient variant adds a small scalar solve inside the collision update, a few percent of total at most. The payoff is statistical accuracy and a correct equilibrium, not speed; it also opens the door to higher weak order at equal cost.

References

1. C. G. Albert, S. V. Kasilov, W. Kernbichler. Symplectic integration with non-canonical quadrature for guiding-center orbits in magnetic confinement devices. J. Comput. Phys. 403, 109065 (2020). arXiv:1903.06885.
2. C. G. Albert, S. V. Kasilov, W. Kernbichler. Accelerated methods for direct computation of fusion alpha particle losses within stellarator optimization. J. Plasma Phys. 86, 815860201 (2020).
3. A. N. Kaufman. Dissipative Hamiltonian systems: a unifying principle. Phys. Lett. A 100, 419 (1984).
4. P. J. Morrison. A paradigm for joined Hamiltonian and dissipative systems. Physica D 18, 410 (1986).
5. M. Grmela, H. C. Öttinger. Dynamics and thermodynamics of complex fluids: the GENERIC formalism. Phys. Rev. E 56, 6620 and 6633 (1997).
6. M. Kraus, E. Hirvijoki. Metriplectic integrators for the Landau collision operator. Phys. Plasmas 24, 102311 (2017).
7. E. Hirvijoki. Structure-preserving marker-particle discretizations of Coulomb collisions for particle-in-cell codes. Plasma Phys. Control. Fusion 63, 044003 (2021). arXiv:2012.07187.
8. J. A. Carrillo, J. Hu, L. Wang, J. Wu. A particle method for the homogeneous Landau equation. J. Comput. Phys. X 7, 100066 (2020).
9. Multispecies structure-preserving particle discretization of the Landau collision operator. Phys. Plasmas 29, 123906 (2022). arXiv:2212.00663.
10. Structure-preserving particle methods for the Landau collision operator using the metriplectic framework. arXiv:2404.18432 (2024).
11. Moment-preserving Monte-Carlo Coulomb collision method for particle codes. arXiv:2407.19151 (2024).
12. A multiscale hybrid Maxwellian-Monte-Carlo Coulomb collision algorithm for particle simulations. arXiv:2405.09573 (2024).
13. X. Shang, H. C. Öttinger. Structure-preserving integrators for dissipative systems based on reversible-irreversible splitting. Proc. R. Soc. A 476, 20190446 (2020).
14. A. H. Boozer, G. Kuo-Petravic. Monte Carlo evaluation of transport coefficients. Phys. Fluids 24, 851 (1981).
15. R. I. McLachlan, G. R. W. Quispel, N. Robidoux. Geometric integration using discrete gradients. Phil. Trans. R. Soc. A 357, 1021 (1999).