PaperlayerView exact revision on arXiv

Source-linked AI summary

An Energy-Conserving Unstaggered Electromagnetic-Potential Particle-in-Cell Method, Part II: Fully Relativistic Formulation with an Adapted Higuera--Cary Pusher

Andrew Christlieb, Luis Chacon, Sining Gong

arXiv:2609.36383v1math.NAastro-ph.HEcs.NA

TL;DR

Electromagnetic PIC needs compatible particle–mesh work exchange and field constraints, while preserving particle geometry. This paper adapts a generalized-momentum potential method to relativistic particles; its coupled scheme conserves energy to numerical tolerances, with particle-volume preservation proved for frozen orbit data.

  • Problem

    Electromagnetic PIC must coordinate charge conservation, field constraints, particle–mesh work exchange, and particle geometry for long-time reliability.

    Source: §1 · p. 2

  • Method

    The method combines a generalized-momentum Higuera–Cary particle update with Crank–Nicolson potential fields and orbit-consistent particle–mesh coupling.

    Source: §3.3 · p. 10, §6 · p. 24

  • Results

    Four projected runs conserve total energy within 6.18 × 10^-12 relative error; TSI diagnostics agree through T = 60, while Weibel nonlinear states differ at T = 100.

    Source: §6 · p. 24, §5.3 · p. 24

  • Takeaways & Limitations

    The converged coupled step conserves relativistic total energy up to solver tolerance, orbit-quadrature error, and roundoff, while volume preservation applies to the frozen-orbit map.

    Source: §6 · p. 24

    The volume-preservation theorem covers prescribed fields and frozen orbit data; a determinant theorem for the fully coupled implicit residual is left to future work.

    Source: §3.4 · p. 11

Abstract

from arXiv · show

We develop a fully relativistic extension of the energy-conserving, unstaggered electromagnetic-potential particle-in-cell method introduced in our earlier nonrelativistic work. As in the nonrelativistic generalized-momentum formulation, the fields are advanced by a Crank--Nicolson discretization of the Lorenz-gauge potential equations, the current is deposited from particle orbits, and charge is advanced through the discrete continuity equation. The new ingredient is a generalized-momentum adaptation of the Higuera--Cary relativistic particle pusher. The vector-potential force is decomposed into an orbit-discrete-gradient part, which enforces the same finite-difference chain rule used in our earlier nonrelativistic formulation, and a skew part, which generates a Higuera--Cary magnetic rotation in mechanical momentum. For prescribed fields and frozen orbit data, this generalized-momentum Higuera--Cary map is a unit-Jacobian conjugacy of the mechanical Higuera--Cary map and preserves particle phase volume in $(\xx,\PP)$. For the fully coupled implicit step, the same orbit data define the current deposit, mesh-to-particle gather, vector-potential chain rule, and relativistic kinetic-energy secant velocity. The resulting particle work cancels the Crank--Nicolson field-energy balance, yielding exact relativistic total-energy conservation up to nonlinear solver tolerance, orbit-quadrature error, and roundoff. A self-adjoint Fourier projection removes modes on even-grid Nyquist planes while preserving the particle--grid work identity. Numerical tests are organized around cold relativistic two-stream and Weibel/filamentation instabilities, with diagnostics for energy drift, Gauss's law, the Lorenz gauge, particle--grid work, and the orbit chain rule.

1. Introduction.

Reliable electromagnetic PIC requires compatible charge conservation, field constraints, particle–mesh work exchange, and particle geometry. This paper extends the prior generalized-momentum potential method to relativistic particles with an adapted Higuera–Cary push.

Motivation and contribution

  • Long-time PIC reliability depends on compatible charge conservation, electromagnetic constraints, particle–mesh work exchange, and particle geometry.The paper frames its goal alongside variational, compatible finite-element, and constrained-splitting approaches that preserve related structures.Source: §1 · pp. 1–2
  • The present construction retains Part I’s orbit-averaged maps and exact vector-potential chain rule, which supported energy conservation in the nonrelativistic method.Part I’s fully converged Crank–Nicolson fixed-point method is the foundation for this extension.Source: §1 · pp. 2–3
  • The method extends the generalized-momentum potential construction to relativistic particles by replacing its additive magnetic update with a Higuera–Cary rotation.The orbit-discrete-gradient contribution is retained in operator-applied-to-field form and decomposed into endpoint/chain-rule and skew magnetic parts.Source: §1 · p. 2
  • For the fully coupled implicit step, shared orbit data and a relativistic kinetic-energy secant velocity yield exact total-energy conservation.The coupled result relies on consistent orbit weights, the vector-potential chain rule, and particle–mesh work exchange.Source: §1 · p. 2
  • With prescribed fields and frozen orbit data, the generalized-momentum map has unit Jacobian in (x, P).This is a fixed-field particle-volume result, distinct from the fully coupled energy claim.Source: §1 · p. 2

2. Relativistic model and orbit-averaged discretization.

The section sets out the nondimensional relativistic potential model and generalized-momentum particle variables, then defines compatible Crank–Nicolson field updates and orbit-averaged particle–mesh coupling. An orbit-discrete-gradient chain rule connects the particle update to the energy balance.

2. Relativistic model and orbit-averaged discretization.

  • The section adopts the nondimensional formulation introduced in for its relativistic model and orbit-averaged discretization.Source: §2 · p. 3

2.1. Nondimensional relativistic potential model.

  • The model recovers electromagnetic fields from the potentials, and its energy argument requires consistent coefficients across the potential update, midpoint Maxwell form, and field energy.The nondimensional speed of light is κ; σ1 and σ2 are source coefficients determined by reference scales.Source: §2.1 · p. 3

2.2. Relativistic generalized-momentum particles.

  • Relativistic particles use generalized momentum P = p + qA, whose formulation cancels the explicit time derivative of A from the particle force.The magnetic update is then advanced by an adapted Higuera–Cary rotation rather than an additive midpoint force.Source: §2.2 · pp. 3–4

2.3. Mesh-based field values and the Crank–Nicolson potential update.

  • The Crank–Nicolson field step deposits current from accepted orbits before advancing charge by discrete continuity, propagating the Lorenz gauge and Gauss’s law when their initial constraints hold.The field update uses the midpoint Lorenz gauge.Source: §2.3 · pp. 4–5

2.4. Symmetric orbit-averaged maps between particles and mesh.

  • Orbit-averaged scatter and gather use the same weights, making their adjointness the particle–mesh work equivalence used in the total-energy proof.The orbit velocity consistently defines the orbit, current, and work identity.Source: §2.4 · p. 5

2.5. Orbit-averaged discrete chain rule and discrete gradient.

  • The orbit-discrete-gradient matrix is the orbit-averaged spatial Jacobian of the time-interpolated particle interpolant, enforcing an exact finite-difference chain rule.The chain rule links generalized momentum to relativistic energy conservation, while the skew part supplies the effective magnetic generator; spline-orbit integrals are split at crossed knots for Gaussian quadrature.Source: §2.5 · p. 6

3. Generalized-momentum Higuera–Cary push and fixed-field volume preservation.

The section adapts the relativistic Higuera–Cary (HC) leapfrog rotation to generalized momentum using an orbit-generated skew matrix, while retaining the vector-potential chain rule. With prescribed fields and frozen orbit data, the resulting particle map preserves canonical particle volume.

3.1. HC as a leapfrog rotation map.

  • The HC step retains the kick–rotation–kick momentum geometry while replacing the final leapfrog line with the accepted orbit and its averaged velocity.The mechanical HC rotation is a centered magnetic update; electric half-impulses can change particle energy, whereas the magnetic step does no work.Source: §3.1 · p. 7
  • The HC magnetic substep is a centered rotation that preserves momentum magnitude and relativistic kinetic energy in a purely magnetic field.Its complete nonlinear map is volume preserving; that result depends on the incoming-momentum dependence of the HC Lorentz factor, not just the fixed-scalar Cayley matrix.Source: §3.1 · pp. 7–8, §3.2 · p. 9

3.3. Skew-matrix HC map used by the generalized-momentum method.

  • The generalized-momentum construction replaces the directly gathered magnetic field with a rotation generated by the orbit-gradient skew matrix.The discrete-gradient part participates in cancelling the explicit At contribution, while the skew part generates the HC magnetic rotation.Source: §3.1 · p. 8, §3.3 · p. 10

3.2. Mechanical Higuera–Cary rotation.

  • The mechanical HC magnetic map uses an incoming-momentum-dependent Lorentz factor, and its full nonlinear map—not merely its fixed-factor Cayley matrix—is volume preserving.The fixed-factor Cayley matrix is an orthogonal determinant-one rotation, but that fact alone does not establish the determinant of the complete map.Source: §3.2 · p. 9

3.4. No-explicit-At generalized-momentum HC step.

  • The canonical force avoids an explicit At gather by using discrete-gradient terms in its half transforms, which are algebraically equivalent to mechanical electric half-kicks.The vector-potential chain rule supports cancellation of the explicit temporal-potential contribution.Source: §3.4 · p. 11

3.5. Particle-volume theorem for prescribed fields.

  • With prescribed fields and frozen orbit data, the auxiliary GM–HC particle map preserves canonical particle volume.The theorem assumes the orbit data are fixed for the step; it does not establish a full-system volume theorem for the coupled implicit PIC residual.Source: §3.5 · p. 11, §3.4 · p. 11

4. Nondimensional relativistic energy balance for the orbit-averaged GM–HC–CN fixed point.

The relativistic GM–HC–CN step couples orbit-based particle work to a Crank–Nicolson field-energy balance. Under the stated discretization and convergence assumptions, particle and field work cancel to give exact total-energy conservation.

4.1. Nondimensional relativistic kinetic-energy secant velocity.

  • The relativistic kinetic-energy secant velocity supplies a discrete-gradient relation between momentum states for the kinetic-energy change.The orbit/current velocity is constructed from the two electric-work secant velocities.Source: §4.1 · p. 12

4.2. Mechanical split and magnetic zero work.

  • The HC substep supplies the magnetic rotation, whose skew generator does no work; the orbit velocity and weights must also define the orbit, current deposit, field gather, and particle work consistently.The same orbit data support the work cancellation between particles and fields.Source: §4.2 · p. 13, §4.3 · p. 14

4.3. CN field energy and total energy theorem.

  • Under the theorem’s periodic-boundary, summation-by-parts, gauge, coefficient-scaling, and converged-solver assumptions, the fully discrete relativistic total energy is conserved.In floating-point arithmetic, the observed defect is limited by nonlinear solver tolerance, orbit-quadrature error, chain-rule residual, and roundoff.Source: §4.3 · p. 14
  • The proof combines particle kinetic-energy change with the negative CN field-energy change, cancelling mesh work exactly.The field-energy balance uses the midpoint Lorenz gauge, compatible coefficient scaling, and periodic summation by parts.Source: §4.3 · p. 14

4.4. Diagnostics required by the proof.

  • The diagnostics monitor total relative energy error, Gauss and Lorenz-gauge residuals, the orbit chain-rule residual, and particle–grid work residual.The paper links a large chain-rule residual to checks of orbit splitting, quadrature, and the discrete-gradient construction.Source: §4.3 · p. 15, §4.4 · p. 15

4.5. Implicit energy-conserving PIC solver.

  • The self-adjoint Fourier projection removes even-grid Nyquist-plane modes while retaining the deposit–gather work identity, but does not make spline and spectral differentiation identical on all retained modes.The restriction removes a particular source of disagreement between mesh and particle magnetic representations.Source: §4.5 · p. 16
  • The solver eliminates fields and particle orbits from the outer nonlinear unknowns, which instead enforce consistency of the projected mesh current.A direct Fourier CN field solve is paired with inner particle-orbit iterations and an outer Anderson-accelerated current iteration.Source: §4.4 · p. 15, §4.5 · pp. 16–18

5. Numerical results.

Four production runs test cold relativistic two-stream and Weibel instabilities at smaller and larger fixed time steps. The diagnostics assess seeded-mode growth, discrete identities, time-step sensitivity, Nyquist projection, and spline-degree effects.

5.1. Cold relativistic two-stream instability with v0 = 0.9.

  • Fitted TSI growth rates are 0.09101901 and 0.09095955 for ∆t = 0.1 and 0.3, with relative errors 7.95×10−4 and 1.45×10−3 against the finite-grid reference.The modes subsequently enter nonlinear trapping, with similar first large maxima and T = 60 amplitudes differing by about 1.42%.Source: §5.1 · p. 20
  • Maximum relative total-energy errors are 3.53×10−12 and 6.17×10−12 for the two TSI runs, while Gauss, Lorenz, chain-rule, and particle–grid work diagnostics remain small.The longitudinal RMS velocities differ by approximately 0.050%; this moment comparison does not measure all fine-scale phase-space differences.Source: §5.1 · p. 20

5.2. Cold relativistic Weibel/filamentation instability.

  • Fitted Weibel growth rates are 5.10911239 × 10−1 and 5.11347526 × 10−1 for WBI-0.04 and WBI-0.2, with relative errors 1.07 × 10−3 and 1.93 × 10−3.Both follow the finite-grid prediction during the fitted interval and show similar initial saturation followed by bounded oscillations, though later nonlinear evolution differs.Source: §5.2 · p. 21
  • For WBI-0.04 and WBI-0.2, maximum absolute relative total-energy changes are 3.74 × 10−14 and 4.13 × 10−13, and maximum chain-rule RMS residuals are 6.68 × 10−16 and 5.57 × 10−16.Maximum Gauss and Lorenz RMS residuals and deposit–gather work residuals also remain small in both runs.Source: §5.2 · p. 21
  • Within each instability, the larger time step requires more outer iterations on average, and all runs remain below the outer limit of 80 and inner limit of 64.The particle populations and physical problems differ, so iteration counts are not direct measures of relative computational cost across the two problems.Source: §5.2 · p. 21

5.3. Cross-test solver behavior and summary.

  • At T = 60, charge-weighted transverse RMS is 0.366930 unprojected and 2.28724 × 10−12 projected in the longitudinal TSI setup; this comparison does not establish how unprojected motion develops.Projection removes Nyquist components from the coupled field–particle update rather than directly resetting transverse velocities.Source: §5.3 · pp. 21–22
  • Finite-offset regularity diagnostics give maximum two-sided differences of 2.83 and 4.45 for r = 1, versus 5.96 × 10−8 and 8.81 × 10−8 for r = 2.These diagnostics are consistent with discontinuous first derivatives for r = 1 and continuous first derivatives for r = 2; convergence-triggered refinement occurs only for degree one, not as an accuracy-controlled subcycling study.Source: §5.3 · pp. 23–24
  • TSI shows close agreement in reported longitudinal diagnostics through T = 60, while Weibel conserves energy and resolves early growth but differs appreciably in its nonlinear state at T = 100.Neither conservation nor Nyquist-feedback suppression alone establishes temporal accuracy.Source: §5.3 · p. 24

6. Conclusions.

The relativistic generalized-momentum Higuera–Cary push splits the vector-potential force into an orbit-chain-rule contribution and a skew magnetic rotation. Its frozen-orbit map preserves canonical phase-space volume, while the converged coupled step conserves total relativistic energy within solver, quadrature, and roundoff errors.

6. Conclusions.

  • The generalized-momentum Higuera–Cary push splits the vector-potential force into an orbit-discrete-gradient chain-rule part and a skew part advanced by magnetic rotation.The force avoids an explicit finite-difference approximation of ∂tA.Source: §6 · p. 24
  • The converged coupled step conserves total relativistic energy up to nonlinear solver tolerance, orbit-quadrature error, and roundoff.The four main projected runs remain within 6.18 × 10−12 relative energy error over their reported intervals, without particle time subcycling.Source: §6 · p. 24
  • For prescribed fields and frozen orbit data, the particle map preserves canonical phase-space volume; this result is restricted to the frozen-orbit map.In the coupled step, a shared orbit construction makes particle work cancel the Crank–Nicolson field-energy balance.Source: §6 · p. 24
  • The self-adjoint Nyquist projection preserves particle–grid work and suppresses substantial transverse motion in the matched two-stream comparison.The projected runs’ early fitted growth rates agree with finite-grid predictions within 0.2%.Source: §6 · pp. 24–25
  • The TSI time-step pair agrees closely in reported longitudinal diagnostics through T = 60, whereas the Weibel pair differs in late-time magnetic-mode histories.The auxiliary spline comparison supports quadratic particle-shape robustness; convergence-triggered particle refinement is needed only for the degree-one calculation. The results distinguish conservation and nonlinear convergence from temporal accuracy, and particle enslavement provides a basis for future large-scale parallel implementations with accuracy-controlled subcycling and physics-based preconditioning.Source: §6 · p. 25
Loading 2609.36383v1…