pith. machine review for the scientific record. sign in

arxiv: 2604.20876 · v1 · submitted 2026-04-10 · ⚛️ physics.comp-ph

Recognition: unknown

High-Accuracy Numerical Solutions of Particle Motion in Static Magnetic Fields

Authors on Pith no claims yet

Pith reviewed 2026-05-10 15:34 UTC · model grok-4.3

classification ⚛️ physics.comp-ph
keywords Parker-Sochacki methodLorentz equationsnumerical integrationstatic magnetic fieldsenergy conservationRunge-Kutta methodsparticle trajectoriesplasma simulation
0
0 comments X

The pith

The Parker-Sochacki power series method conserves kinetic energy 4 to 13 orders of magnitude better than Runge-Kutta methods when integrating charged particle motion in static magnetic fields.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The paper tests the Parker-Sochacki method, which builds a time power series directly from the Lorentz equations, against fixed-step fourth-order Runge-Kutta, adaptive Dormand-Prince, and symplectic Gauss-Legendre Runge-Kutta on three static field geometries: uniform, hyperbolic tangent, and dipole. It reports that the power series approach improves kinetic energy conservation by four to thirteen orders of magnitude, runs faster than fourth-order Runge-Kutta at matched error tolerance, and is the only integrator that remains accurate and stable for both protons and electrons across all cases. In the dipole field the symplectic method fails on every electron trajectory while the power series method succeeds with lower error and, for protons, four to five times shorter runtime at fixed step size. These outcomes indicate that the method supplies a practical route to long, high-fidelity trajectory calculations without artificial energy drift.

Core claim

The Parker-Sochacki method constructs a power series solution tailored to the Lorentz force equations and advances the particle state by summing that series at each fixed time step. When applied to protons and electrons in uniform, tanh, and dipole magnetic fields, it produces kinetic energy errors four to thirteen orders of magnitude smaller than those of fourth-order Runge-Kutta or adaptive Dormand-Prince integrators and, for protons in the dipole field, four to five times lower error than the symplectic Gauss-Legendre Runge-Kutta method at identical fixed step size. The same power series integrator is the only method tested that remains stable for electron motion in the dipole field; the

What carries the argument

The Parker-Sochacki power series expansion, which recursively computes time derivatives from the Lorentz equations and sums the resulting series to advance the solution over each time step.

If this is right

  • At matched target kinetic energy error the Parker-Sochacki method completes integrations faster than fixed-step fourth-order Runge-Kutta.
  • The symplectic Gauss-Legendre Runge-Kutta method can exhibit secular growth in energy error during sufficiently long integrations in inhomogeneous fields.
  • Only the Parker-Sochacki method maintains both accuracy and stability for every tested combination of field geometry and particle species.
  • Fixed-time-step Parker-Sochacki runs achieve the reported accuracy levels without requiring adaptive step-size control.

Where Pith is reading between the lines

These are editorial extensions of the paper, not claims the author makes directly.

  • Because the method depends on rapid series convergence, it may be combined with occasional coefficient recalculation to handle mildly time-varying fields over short static intervals.
  • The observed stability advantage for light particles suggests the integrator could reduce the number of time steps needed in large-scale simulations of cosmic-ray or space-plasma trajectories.
  • Replacing standard Runge-Kutta codes with Parker-Sochacki in existing particle-tracing packages would lower the computational cost of reaching a target statistical precision in long-duration runs.

Load-bearing premise

The magnetic fields remain perfectly static and the power series converges rapidly enough for the chosen fixed time steps and field geometries tested.

What would settle it

Re-running the electron trajectories in the dipole field with the Parker-Sochacki method at the same fixed time step used for the failed Gauss-Legendre runs; if the energy error grows at a rate comparable to the Runge-Kutta methods or the integration becomes unstable, the reported superiority would be falsified.

Figures

Figures reproduced from arXiv: 2604.20876 by Heather Jiles, Robert Weigel.

Figure 1
Figure 1. Figure 1: Comparison of PS, RK4, and RK45 methods for a 100 eV electron in a uniform 0.01 T magnetic field using 64-bit [PITH_FULL_IMAGE:figures/full_fig_p010_1.png] view at source ↗
Figure 2
Figure 2. Figure 2: Comparison of PS, RK4, and RK45 methods for a 10 keV electron in a hyperbolic tangent magnetic field ( [PITH_FULL_IMAGE:figures/full_fig_p014_2.png] view at source ↗
Figure 3
Figure 3. Figure 3: Comparison of PS, RK4, and RK45 methods for a 10 keV electron in a hyperbolic tangent magnetic field ( [PITH_FULL_IMAGE:figures/full_fig_p015_3.png] view at source ↗
Figure 4
Figure 4. Figure 4: Comparison of PS, RK4, and RK45 methods for a 100 keV proton in a hyperbolic tangent magnetic field ( [PITH_FULL_IMAGE:figures/full_fig_p016_4.png] view at source ↗
Figure 5
Figure 5. Figure 5: Comparison of PS, RK4, and RK45 methods for a 10 keV electron in a hyperbolic tangent magnetic field ( [PITH_FULL_IMAGE:figures/full_fig_p017_5.png] view at source ↗
Figure 6
Figure 6. Figure 6: Comparison of the PS, RK45, and RKG methods for a 100 keV proton in Earth’s dipole magnetic field, initialized [PITH_FULL_IMAGE:figures/full_fig_p022_6.png] view at source ↗
Figure 7
Figure 7. Figure 7: From the same run as Figure 6, comparison of the PS, RK45, and RKG methods for a 100 keV proton in Earth’s [PITH_FULL_IMAGE:figures/full_fig_p023_7.png] view at source ↗
Figure 8
Figure 8. Figure 8: Comparison of the PS, RK4, RK45, and RKG methods for protons in Earth’s dipole magnetic field, [PITH_FULL_IMAGE:figures/full_fig_p024_8.png] view at source ↗
Figure 9
Figure 9. Figure 9: Comparison of bounce and drift periods for protons with equatorial pitch angle [PITH_FULL_IMAGE:figures/full_fig_p026_9.png] view at source ↗
Figure 10
Figure 10. Figure 10: Comparison of the PS, RK4, and RK45 methods for a 100 MeV electron in Earth’s dipole magnetic field, initialized [PITH_FULL_IMAGE:figures/full_fig_p027_10.png] view at source ↗
Figure 11
Figure 11. Figure 11: From the same run as Figure 10, comparison of the PS, RK4, and RK45 methods for a 100 MeV electron in Earth’s [PITH_FULL_IMAGE:figures/full_fig_p028_11.png] view at source ↗
Figure 12
Figure 12. Figure 12: Comparison of PS, RK4, RK45, and RKG methods for an electron in Earth’s dipole magnetic [PITH_FULL_IMAGE:figures/full_fig_p029_12.png] view at source ↗
Figure 13
Figure 13. Figure 13: Comparison of bounce and drift periods for electrons with equatorial pitch angle [PITH_FULL_IMAGE:figures/full_fig_p030_13.png] view at source ↗
Figure 14
Figure 14. Figure 14: Error-matched comparison between the PS method and the fixed-time-step RK4 method for a 100 keV proton in [PITH_FULL_IMAGE:figures/full_fig_p032_14.png] view at source ↗
read the original abstract

The Parker-Sochacki (PS) method is investigated as an alternative to Runge-Kutta (RK) methods for solving the Lorentz equations of motion for a charged particle in a static magnetic field. Traditional methods, including fixed-time-step fourth-order RK, adaptive Dormand-Prince RK, and Gauss-Legendre Runge-Kutta (RKG), advance the solution by sampling derivative estimates at selected points to approximate the solution over a time increment. In contrast, the PS method uses a power series expansion in time that is specific to the system of equations, which is a fundamentally different approach. We assess the accuracy and long-term stability of the RK, RKG, and PS methods for three static magnetic fields: uniform, hyperbolic tangent, and dipole, with the RKG method included only for the dipole problem. The PS method results in a 4 to 13 orders-of-magnitude improvement in kinetic energy conservation over the RK methods. When the methods are compared at matched target kinetic energy error, the PS method was substantially faster than RK4, the method with the shortest runtime under identical fixed-time-step conditions. For the dipole field problem, the PS method had the lowest kinetic energy error and had runtimes 4 to 5 times shorter than RKG when using the same fixed time step for proton runs. The PS method was the only method in this study to maintain accuracy and stability for all problems for both protons and electrons; the RKG method failed on all electron runs in the dipole problem. We further show that, over sufficiently long integrations in inhomogeneous magnetic fields, the symplectic RKG may exhibit secular growth in energy error. Overall, these results indicate that the PS method provides a computationally efficient and highly accurate alternative to the symplectic RKG and standard RK methods.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

2 major / 1 minor

Summary. The manuscript presents a comparison of the Parker-Sochacki (PS) power series method with Runge-Kutta (RK4, adaptive Dormand-Prince) and Gauss-Legendre Runge-Kutta (RKG) methods for integrating the Lorentz equations of motion for charged particles in static magnetic fields. Using test cases of uniform, hyperbolic tangent, and dipole magnetic fields, the authors report that the PS method achieves 4 to 13 orders of magnitude better conservation of kinetic energy, is faster at matched error levels, and is the only method that remains stable and accurate for both protons and electrons across all problems, while RKG fails for electrons in the dipole field. It also notes potential secular energy error growth in long RKG integrations.

Significance. If the numerical results hold under scrutiny, the PS method could represent a valuable tool for high-accuracy trajectory integrations in magnetized plasmas, offering better long-term conservation without adaptive stepping. The direct head-to-head comparisons on fixed problems provide useful benchmarks for the field. However, the reliance on fixed time steps in inhomogeneous fields where frequencies vary spatially requires careful validation to ensure the improvements are not due to unequal effective resolution.

major comments (2)
  1. [Numerical results for dipole field] The claim that PS is the only stable method for electrons rests on fixed Δt comparisons, but the local cyclotron frequency varies by orders of magnitude in the dipole (B ∝ 1/r³); without reporting Δt relative to the minimum gyro-period or local truncation error estimates, the stability advantage may not generalize or may reflect insufficient resolution for RK/RKG rather than inherent superiority of PS.
  2. [Methods and implementation] Details on the PS series truncation (order selection, remainder bound) and the precise protocol for measuring kinetic energy error (integration duration, sampling, definition of 'error') are missing, making it difficult to assess the reported 4-13 orders-of-magnitude gains or to reproduce the runtime comparisons.
minor comments (1)
  1. [Abstract] The abstract mentions 'chosen fixed time steps' but could specify the criteria used for selecting those steps across methods and particles.

Simulated Author's Rebuttal

2 responses · 0 unresolved

We thank the referee for their careful and constructive review of our manuscript. We address each major comment below and have revised the paper to incorporate the requested details and clarifications.

read point-by-point responses
  1. Referee: The claim that PS is the only stable method for electrons rests on fixed Δt comparisons, but the local cyclotron frequency varies by orders of magnitude in the dipole (B ∝ 1/r³); without reporting Δt relative to the minimum gyro-period or local truncation error estimates, the stability advantage may not generalize or may reflect insufficient resolution for RK/RKG rather than inherent superiority of PS.

    Authors: We agree that reporting the time step relative to the local and minimum gyro-period is essential for fair assessment in inhomogeneous fields. In our tests, Δt was chosen to resolve the highest cyclotron frequency encountered in each domain (for electrons in the dipole, Δt ≈ T_gyro,min / 100 near the origin). We have added a dedicated subsection in the revised manuscript that tabulates Δt / T_gyro,min for every particle and field, together with a priori local truncation-error estimates for each integrator. These additions confirm that all methods operated at comparable effective resolution, and the observed stability advantage of PS is not an artifact of under-resolution. We have also expanded the discussion of generalization limits. revision: yes

  2. Referee: Details on the PS series truncation (order selection, remainder bound) and the precise protocol for measuring kinetic energy error (integration duration, sampling, definition of 'error') are missing, making it difficult to assess the reported 4-13 orders-of-magnitude gains or to reproduce the runtime comparisons.

    Authors: We thank the referee for highlighting these omissions. The revised Methods section now specifies: (i) the PS truncation rule, in which the power-series order is increased until the remainder bound (computed from the next term via the Cauchy estimate) drops below 10^{-15}; (ii) the exact kinetic-energy error protocol—integrations of 10^4 gyroperiods, error sampled at every step, and error defined as max_t |K(t) − K(0)| / K(0). These additions enable full reproduction of the accuracy and runtime results. revision: yes

Circularity Check

0 steps flagged

No circularity: performance claims rest on direct numerical benchmarks

full rationale

The paper reports measured outcomes from fixed-time-step integrations of the Lorentz equations in three static magnetic fields (uniform, tanh, dipole) for both protons and electrons. The 4–13 order-of-magnitude kinetic-energy conservation advantage, runtime comparisons at matched error, and stability statements are presented as empirical results from those runs; no equation, parameter fit, or self-citation is shown to reduce the claimed superiority to a tautology or to the input data by construction. The PS power-series description is a standard, non-circular algorithmic choice for the ODE system and does not feed back into the performance metrics. The derivation chain is therefore self-contained against external test problems.

Axiom & Free-Parameter Ledger

0 free parameters · 1 axioms · 0 invented entities

The work is a numerical methods comparison and introduces no new physical parameters, axioms, or entities.

axioms (1)
  • domain assumption The magnetic fields are time-independent and the solution remains sufficiently smooth for the chosen integrators.
    Required for both the Lorentz equations and the power-series construction to remain valid.

pith-pipeline@v0.9.0 · 5622 in / 1166 out tokens · 41604 ms · 2026-05-10T15:34:21.002739+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

35 extracted references · 21 canonical work pages

  1. [1]

    Boyce and R

    W. Boyce and R. DiPrima.Elementary Differential Equations and Boundary Value Problems. 10th. Wiley, 2017

  2. [2]

    Chapra and R

    S. Chapra and R. Canale.Numerical Methods for Engineers. McGraw-Hill, 2015

  3. [3]

    J. C. Butcher.Numerical Methods for Ordinary Differential Equations. 3rd. Wiley, 2016

  4. [4]

    A Comparison of Explicit Runge-Kutta Methods

    S. Walters, R. Turner, and L. Forbes. “A Comparison of Explicit Runge-Kutta Methods”. In:The ANZIAM Journal64 (2022), pp. 227–249

  5. [5]

    Iserles.A First Course in the Numerical Analysis of Differential Equations

    A. Iserles.A First Course in the Numerical Analysis of Differential Equations. Cambridge University Press, 1996

  6. [6]

    Implementing the Picard Iteration

    G. Parker and S. Sochacki. “Implementing the Picard Iteration”. In:Neural, Parallel, and Scientific Computations4 (1996), pp. 97–112

  7. [7]

    An Effective Approach for Solving Nonlinear Fractional Initial Value Problems: The Fractional Legendre-PicardIterationMethod

    S. Ansari and M. H. Akrami. “An Effective Approach for Solving Nonlinear Fractional Initial Value Problems: The Fractional Legendre-PicardIterationMethod”.In:JournalofMathematicalExtension18.3(2024),pp.1–29.doi:10.30495/JME.2024. 2941

  8. [8]

    A comparison of HPM, NDHPM, Picard and Picard–Pade´methods for solving Michaelis–Menten equation

    H. Vazquez-Leal et al. “A comparison of HPM, NDHPM, Picard and Picard–Pade´methods for solving Michaelis–Menten equation”. In:Journal of King Saud University27 (2015), pp. 7–14.doi:10.1016/j.jksus.2014.11.001

  9. [9]

    J. W. Rudmin.Application of the Parker-Sochacki Method to Celestial Mechanics. Tech. rep. Physics Dept., James Madison University, 1998

  10. [10]

    AnAdaptive,HighlyAccurateandEfficient,Parker-SochackiAlgorithmforNumericalSolutionsto Initial Value Ordinary Differential Equation Systems

    J.GuentherandM.Wolf.“AnAdaptive,HighlyAccurateandEfficient,Parker-SochackiAlgorithmforNumericalSolutionsto Initial Value Ordinary Differential Equation Systems”. In:SIAM Undergraduate Research Online12 (2019).issn: 2327-7807. doi:10.1137/19s019115

  11. [11]

    Adiabatic Charged-Particle Motion

    T. G. Northrop. “Adiabatic Charged-Particle Motion”. In:Reviews of Geophysics1.3 (1963), pp. 283–304.doi:10.1029/ RG001i003p00283

  12. [12]

    InsolubilityofTrappedParticleMotioninaMagneticDipoleField

    A.DragtandJ.Finn.“InsolubilityofTrappedParticleMotioninaMagneticDipoleField”.In:JournalofGeophysicalResearch 81.13 (1976), pp. 2327–2347.doi:10.1029/JA081i013p02327

  13. [13]

    On a plasma sheath separating regions of oppositely directed magnetic field

    E. G. Harris. “On a plasma sheath separating regions of oppositely directed magnetic field”. In:Il Nuovo Cimento23.1 (Jan. 1962), pp. 115–121.issn: 1827-6121.doi:10.1007/bf02733547

  14. [14]

    Knuth.The Art of Computer Programming, Vol

    D. Knuth.The Art of Computer Programming, Vol. 1: Fundamental Algorithms (3rd ed.)See Chapter 1, Section 1.2.9. Addison-Wesley, 1997

  15. [15]

    Spiking neural network simulation: numerical integration with the Parker–Sochacki method

    R. Stewart and W. Bair. “Spiking neural network simulation: numerical integration with the Parker–Sochacki method”. In: Journal of Computational Neuroscience27.1 (2009), pp. 115–133.doi:10.1007/s10827-008-0121-5

  16. [16]

    IEEE Standard for Floating-Point Arithmetic.IEEE Std 754-2019 (Revision of IEEE 754-2008)(2019), 1–84

    “IEEEStandardforFloating-PointArithmetic”.In:IEEEStd754(2019),pp.1–84.doi:10.1109/IEEESTD.2019.8766229

  17. [17]

    Machine behaviour

    C.R.Harrisetal.“ArrayprogrammingwithNumPy”.In:Nature585.7825(Sept.2020),pp.357–362.doi:10.1038/s41586- 020-2649-2

  18. [18]

    K., Pitrou, A., & Seibert, S

    S.K.Lam,A.Pitrou,andS.Seibert.“Numba:ALLVM-BasedPythonJITCompiler”.In:ProceedingsoftheSecondWorkshop on the LLVM Compiler Infrastructure in HPC. 2015, pp. 1–6.doi:10.1145/2833157.2833162

  19. [19]

    1980 , issn =

    J.R. Dormand and P.J. Prince. “A family of embedded Runge-Kutta formulae”. In:Journal of Computational and Applied Mathematics6.1 (Mar. 1980), pp. 19–26.issn: 0377-0427.doi:10.1016/0771-050x(80)90013-3

  20. [20]

    E., et al

    P. Virtanen et al. “SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python”. In:Nature Methods17 (2020), pp. 261–272.doi:10.1038/s41592-019-0686-2

  21. [21]

    F. F. Chen.Introduction to Plasma Physics and Controlled Fusion (3rd ed., Vol. 1). Springer, 2016, pp. 19–21

  22. [22]

    Birdsall and A

    C. Birdsall and A. Langdon.Plasma Physics via Computer Simulation. See Chapter 4, Section 4.2. Institute of Physics Publishing, 1991

  23. [23]

    Kivelson and C

    M. Kivelson and C. Russell.Introduction to Space Physics. Cambridge University Press, 1995, pp. 250–252

  24. [24]

    Imprints of Quasi-Adiabatic Ion Dynamics on the Current Sheet Structures Observed in the Martian Magnetotail by MAVEN

    E. E. Grigorenko et al. “Imprints of Quasi-Adiabatic Ion Dynamics on the Current Sheet Structures Observed in the Martian Magnetotail by MAVEN”. In:Journal of Geophysical Research: Space Physics122.10 (2017), pp. 10, 176–10, 193.doi: 10.1002/2017JA024216

  25. [25]

    Statistical Survey of Thin Current Sheets in Earth’s Magnetotail: MMS Observations

    Z. Zhang et al. “Statistical Survey of Thin Current Sheets in Earth’s Magnetotail: MMS Observations”. In:Journal of Geophysical Research: Space Physics129.5 (2024).doi:10.1029/2024JA032575

  26. [26]

    Stiffness and Nonstiff Differential Equation Solvers, II: Detecting Stiffness with Runge–Kutta Methods

    L. F. Shampine. “Stiffness and Nonstiff Differential Equation Solvers, II: Detecting Stiffness with Runge–Kutta Methods”. In: ACM Transactions on Mathematical Software3.1 (1977), pp. 44–53.doi:10.1145/355719.355723

  27. [27]

    Walt.Introduction to Geomagnetically Trapped Radiation

    M. Walt.Introduction to Geomagnetically Trapped Radiation. Cambridge University Press, 1994. 34

  28. [28]

    Symplecticintegration:Anewapproachtotracingchargedparticlemotioninthegeomagneticfield

    H.YugoandT.Iyemori.“Symplecticintegration:Anewapproachtotracingchargedparticlemotioninthegeomagneticfield”. In:Journal of Geophysical Research106.A11 (2001), pp. 26, 075–26, 079.doi:10.1029/2000JA000424

  29. [29]

    Symplectic tracing of high-energy chraged particles in the inner magnetosphere

    H. Yugo and T. Iyemori. “Symplectic tracing of high-energy chraged particles in the inner magnetosphere”. In:Journal of Geophysical Research112.A06243 (2007).doi:10.1029/2006JA011601

  30. [30]

    WORLD SCIENTIFIC, 1986.isbn: 9789814503150.doi: 10.1142/0287.url:http://dx.doi.org/10.1142/0287

    S M Omohundro.Geometric Perturbation Theory in Physics. WORLD SCIENTIFIC, 1986.isbn: 9789814503150.doi: 10.1142/0287.url:http://dx.doi.org/10.1142/0287

  31. [31]

    J. G. Roederer.Dynamics of Geomagnetically Trapped Radiation. Physics and Chemistry in Space. Berlin, Heidelberg: Springer-Verlag, 1970.isbn: 978-3-642-49300-3

  32. [32]

    Pitchangledistributionofmagnetospherictrappedparticles:Atest-particlesimulation

    P.K.Soni,B.Kakad,andA.Kakad.“Pitchangledistributionofmagnetospherictrappedparticles:Atest-particlesimulation”. In:Advances in Space Research68.8 (2021), pp. 3381–3390.issn: 0273-1177.doi:10.1016/j.asr.2021.06.004

  33. [33]

    On the validity of the guiding-center approximation in a magnetic dipole field

    A. J. Brizard and D. G. Markowski. “On the validity of the guiding-center approximation in a magnetic dipole field”. In: arXiv:2111.05353 [physics.plasm-ph](2021).url:https://arxiv.org/abs/2111.05353

  34. [34]

    ModelingtheEffectsofDriftOrbitBifurcationonRadiationBeltElectrons

    J.Huang,W.Tu,andW.W.Eshetu.“ModelingtheEffectsofDriftOrbitBifurcationonRadiationBeltElectrons”.In:Journal of Geophysical Research: Space Physics127.11 (2022).issn: 2169-9402.doi:10.1029/2022ja030827

  35. [35]

    Advanced numerical techniques for time integration of relativistic equations of motion for charged particles

    T. Umeda and T. Ozaki. “Advanced numerical techniques for time integration of relativistic equations of motion for charged particles”. In:Earth, Planets and Space75.1 (2023).issn: 1880-5981.doi:10.1186/s40623-023-01902-8. 35 A Supplemental Tables Table A.1: Comparison of mean and maximum relative kinetic energy error and runtime for particle simulations i...