A Comparative Analysis of Relativistic Particle Pushers vis-à-vis Computation Time & Accuracy

Basic Information

• Paper ID: 2410.03352
• Title: A Comparative Analysis of Relativistic Particle Pushers vis-à-vis Computation Time & Accuracy
• Authors: Mohammad Yasir, Vikrant Saxena (Indian Institute of Technology Delhi)
• Classification: physics.plasm-ph (Plasma Physics)
• Publication Date: October 7, 2024
• Paper Link: https://arxiv.org/abs/2410.03352

Abstract

The performance of relativistic particle pushers has been a research focus in computational plasma physics, particularly in particle-in-cell (PIC) methods. Previous comparative studies have primarily targeted the ultrarelativistic regime. This paper employs the self-developed code PaTriC to investigate the core runtime of three particle pushers—Boris, Vay, and Higuera-Cary—across non-relativistic, relativistic, and ultrarelativistic particle regimes, subsequently comparing the accuracy and errors of the three integrators. The paper proposes a "fitness parameter" that serves as a comprehensive metric for determining which method is most suitable for specific simulation scenarios.

Research Background and Motivation

1. Problem Statement

Particle trajectory calculation is one of the most computationally intensive steps in relativistic plasma dynamics simulations. Multiple explicit integrators (Boris, Vay, Higuera-Cary, etc.) each have advantages and disadvantages, but systematic comparative studies across different relativistic regimes (non-relativistic, relativistic, ultrarelativistic) are lacking.

2. Problem Significance

• Computational Bottleneck: Plasma systems involve multiple temporal, energetic, and spatial scales (spatial scales spanning 10^8, temporal scales from femtoseconds to days), with limited computational resources
• Broad Applications: Particle trajectory calculations are not only used in plasma simulations but also applied to molecular dynamics, astrophysics, and other fields
• Speed-Accuracy Trade-off: Requires finding optimal balance between computational speed and simulation accuracy

3. Limitations of Existing Methods

• Research Bias: Previous comparative studies have primarily focused on the ultrarelativistic regime (e.g., Ripperda et al. 2017)
• Lack of Unified Standards: No unified evaluation metrics comprehensively considering computational cost and accuracy
• Unclear Application Scenarios: Researchers struggle to select the most appropriate integrator for specific problems

4. Research Motivation

Develop a comprehensive performance evaluation framework covering all relativistic regimes (LR, HR, UR) and propose quantitative metrics to help researchers select optimal integrators for specific applications.

Core Contributions

1. Systematic Comparative Study: First comprehensive comparison of Boris, Vay, and Higuera-Cary pushers across three regimes: non-relativistic (LR), relativistic (HR), and ultrarelativistic (UR)
2. Precise Computational Cost Analysis:
   • Pure magnetic field: Boris (24 FLOPs), Vay (41 FLOPs), HC (38 FLOPs)
   • Electromagnetic field: Boris (55 FLOPs), Vay (91 FLOPs), HC (88 FLOPs)
   • Measured runtime data (600,000 iterations)
3. Multi-dimensional Accuracy Assessment: Systematic evaluation of phase error, gyroradius error, relativistic factor γ error, and cross-field drift performance
4. Fitness Parameter Proposal: Innovatively proposes fitness parameter: f = (1/κ)e^(-ε), where κ is computational cost and ε is the logarithm of relative error, providing quantitative criteria for integrator selection
5. Custom Code Development: Constructed PaTriC (Particle Tracker in C++) for single-particle simulation and rigorous analysis

Methodology Details

Task Definition

Numerical integration of relativistic charged particle trajectories in electromagnetic fields:

• Input: Initial position x⃗, velocity u⃗ (u⃗ = γv⃗), electromagnetic fields E⃗ and B⃗, time step δt
• Output: Particle position and velocity at each time step
• Constraints: Maintain time reversibility, near-symplecticity, energy conservation (without electric field)

Leapfrog Integrator Framework

All three methods are based on the staggered leapfrog framework with equations of motion:

$\frac{d\vec{u}}{dt} = \frac{q}{m}(\vec{E} + \bar{\vec{v}} \times \vec{B})$

$\frac{d\vec{x}}{dt} = \bar{\vec{v}}$

where γ = 1/√(1-(v/c)²), with the key difference being different choices for average velocity v̄.

Discrete form: $\vec{u}_{i+1} = \vec{u}_i + (q/m)(\vec{E}_{i+1/2} + \bar{\vec{v}} \times \vec{B}_{i+1/2}) dt$

Common symbol definitions:

• f = qδt/(2m)
• ε⃗ = fE⃗
• β⃗ = fB⃗
• Γ(u⃗) = √(1 + u⃗·u⃗/c²)

1. Boris Method

Mathematical Expression: $\bar{\vec{v}} = (2\Gamma(\vec{u}_e))^{-1}(\vec{u}_{i+1} + \vec{u}_i)$

Update Steps:

u⃗_e = u⃗_i + ε⃗
τ⃗ = β⃗/Γ(u⃗_e)
s⃗ = 2τ⃗/(1 + τ⃗²)
u⃗_m = u⃗_e + (u⃗_e + (u⃗_e × τ⃗)) × s⃗
u⃗_{i+1} = u⃗_m + ε⃗

Computational Cost:

• Pure magnetic field: 24 FLOPs
• Electromagnetic field: 55 FLOPs

Characteristics:

• Classical method, widely used in EPOCH, SMILEI, and other software
• Near-symplectic integrator with excellent accuracy
• Lowest computational cost
• Uses simplified form (original form more accurate but 46% higher cost)

2. Vay Method

Design Motivation: Boris pusher exhibits spurious forces in non-trivial cases, leading to trajectory deviations in relativistic situations

Update Steps:

u⃗_{i+1/2} = u⃗_i + f(E⃗_{i+1/2} + v⃗_i × B⃗_{i+1/2})
u⃗_e = u⃗_{i+1/2} + ε⃗
u⃗_{i+1} = s(u⃗_e + (u⃗_e · t⃗)t⃗ + u⃗_e × t⃗)

where:

• u* = u⃗_e · β⃗/c
• σ = (Γ(u⃗_e))² - β⃗²
• γ_{i+1} = √((σ + √(σ + 4(β⃗² + u*²)))/2)
• t⃗ = β⃗/γ_{i+1}
• s = 1/(1 + t⃗²)

Computational Cost:

• Pure magnetic field: 41 FLOPs
• Electromagnetic field: 91 FLOPs

Characteristics:

• Prevents spurious forces, maintains Lorentz invariance
• Excellent cross-field drift performance
• Higher computational cost

3. Higuera-Cary (HC) Method

Design Motivation: Seeks an integrator that correctly captures cross-field drift while preserving phase space volume (Vay does not preserve volume)

Average Velocity Definition: $\bar{\vec{v}} = (1/2\bar{\gamma})^{-1}(\vec{u}_i + \vec{u}_{i+1})$

where: $\bar{\gamma} = \Gamma\left(\frac{\vec{u}_i + \vec{u}_{i+1}}{2}\right)$

Update Steps:

u⃗_e = u⃗_i + ε⃗
u⃗_m = s(u⃗_e + (u⃗_e · t⃗)t⃗ + u⃗_e × t⃗)
u⃗_{i+1} = u⃗_m + ε⃗ + (u⃗_m × t⃗)

Computational Cost:

• Pure magnetic field: 38 FLOPs
• Electromagnetic field: 88 FLOPs

Characteristics:

• Preserves phase space volume
• Excellent cross-field drift performance
• Smaller phase error
• Computational cost between Boris and Vay

Technical Innovations

1. Comprehensive Regime Coverage: Unlike previous studies concentrated on UR regime, this research covers γ ranging from 1.27 to 56.95
2. Actual Runtime Measurement: Beyond theoretical FLOPs calculation, actual runtime measured on Intel Xeon processor for 600,000 iterations
3. Multi-dimensional Error Analysis:
   • Phase error
   • Gyroradius error
   • Relativistic factor γ error
   • Cross-field drift slope error
4. Fitness Parameter Design:
   • Comprehensively considers computational cost and error accumulation
   • Uses error from 5 slowest process cycles for evaluation
   • Time step set to 1/50 of fastest process cycle
   • Exponential form enhances discrimination

Experimental Setup

Test Particle

Positron, with charge sign choice not affecting accuracy results

Relativistic Regime Classification

1. Low-speed Relativistic (LR):
   • γ = 1.27282
   • ω = 6.90×10^11 rad/s
   • Time step: dt ≈ T_c/18
2. High-speed Relativistic (HR):
   • γ = 3.589
   • T_c = 2.56545×10^-11 s
   • Time step: dt ≈ T_c/25
3. Ultrarelativistic (UR):
   • γ = 56.95
   • T_c = 4.07×10^-10 s
   • Time step: dt ≈ T_c/20

Test Scenarios

1. Magnetic Gyration:

• Pure magnetic field: B⃗ = 5T (ẑ direction)
• Particle executes helical motion
• Evaluates trajectory accuracy, phase error, gyroradius error

2. Cross-Field Drift:

• Magnetic field: B⃗ = 10^-4 T (ẑ direction)
• Electric field: E⃗ = 100 V/m (x̂ direction)
• Drift velocity: |v⃗_D| = E/B = 10^6 m/s (negative ŷ direction)
• UR case analyzed in drift reference frame (κ = 1.002)

Computational Platform

• Processor: Intel Xeon W-1270 @ 3.40GHz
• Memory: 64GB
• Compiler: gcc 13.2.0
• Architecture: x86_64

Evaluation Metrics

1. Core Runtime: Actual CPU time for 600,000 iterations
2. Phase Error: Δφ = φ_simulated - φ_analytical (radians)
3. Gyroradius Relative Error: (r_simulated - r_analytical)/r_analytical
4. γ Relative Error: (γ_simulated - γ_analytical)/γ_analytical
5. Cross-field Drift Slope Error: Slope deviation of ⟨y⟩ over time

Implementation Details

• Fixed Time Step: Avoids variable step length breaking symplecticity
• Coarse Grid Testing: Uses T_c/18 to T_c/25 coarse grids to amplify errors
• Multiple Run Averaging: Each integrator run 10 times with averaging to eliminate noise
• Data Output Frequency: Record once every 25,000 iterations

Experimental Results

Main Results

1. Core Runtime Comparison

Pure Magnetic Field (600,000 iterations):

Electromagnetic Field (600,000 iterations):

Key Findings:

• Boris method fastest in all scenarios
• Vay and HC runtimes similar, consistent with theoretical expectations
• Electric field introduction affects different methods differently

2. Low-speed Relativistic (LR) Performance

Magnetic Gyration (γ=1.27, 8 cycles):

Phase Error:

• Boris: 0.489 rad
• Vay: 0.489 rad
• HC: 0.374 rad ✓ (optimal)

Gyroradius Error (theoretical value 173.74 µm):

• Boris: 0.057 µm (0.033%)
• Vay: 0.057 µm (0.033%) ✓ (optimal)
• HC: 0.45 µm (0.259%)

γ Relative Error:

• Boris: ~10^-9 ✓
• Vay: ~10^-8
• HC: ~10^-9 ✓

Cross-field Drift (7 cycles):

Phase Error:

• Boris: 0.24 rad
• Vay: 0.24 rad
• HC: 0.17 rad ✓ (optimal)

Drift Slope Error:

3. High-speed Relativistic (HR) Performance

Magnetic Gyration (γ=3.589, 6 cycles):

Phase Error:

• Boris: 0.1926 rad
• Vay: 0.1926 rad
• HC: 0.1583 rad ✓ (optimal)

Gyroradius Relative Error (theoretical value 0.4083 mm):

• Boris: 0.0325%
• Vay: 0.0325% ✓
• HC: 0.1141%

γ Relative Error: ~10^-6% (negligible for all methods)

Cross-field Drift (9 cycles):

Phase Error:

• Boris: 0.3264 rad
• Vay: 0.3263 rad
• HC: 0.2663 rad ✓ (optimal)

Drift Slope Error:

4. Ultrarelativistic (UR) Performance

Magnetic Gyration (γ=56.95, 9 cycles):

Phase Error:

• Boris: 0.4173 rad
• Vay: 0.4173 rad
• HC: 0.3213 rad ✓ (optimal, significant difference)

Gyroradius Error (theoretical value 7.4826 mm):

• Boris: 0.00235 mm ✓
• Vay: 0.00235 mm ✓
• HC: 0.015 mm

γ Error Behavior:

• Boris and Vay exhibit staggered pattern
• HC method shows more pronounced stepping
• Vay initial stability inferior

Cross-field Drift (drift reference frame, 8 cycles):

Phase Error:

• Boris: 0.263 rad
• Vay: 0.263 rad
• HC: 0.194 rad ✓ (optimal)

Gyroradius Relative Error (boosted frame):

• Boris: 0.937 ✓
• Vay: 0.939
• HC: 1.1002

γ Relative Error (boosted frame):

• Boris: ~0.1% ✓
• HC: ~0.1% ✓
• Vay: unstable, inferior

Fitness Parameter Results

Definition: f = (1/κ)e^(-ε), where:

• κ: computational cost per iteration (FLOPs)
• ε: log₁₀(relative error over 5 cycles)

Calculation Results:

Interpretation:

• Pure Magnetic Field: Boris f parameter far exceeds others, optimal choice
• Electromagnetic Field: Boris still optimal but advantage reduced; Vay and HC comparable
• f parameter successfully synthesizes computational cost and accuracy, providing quantitative selection basis

Experimental Findings

1. Computational Efficiency:
   • Boris fastest in all scenarios
   • Vay and HC computational costs similar (theory and measurement consistent)
   • Electric field impact differs across methods
2. Accuracy Characteristics:
   • Phase Error: HC optimal across all regimes and scenarios
   • Gyroradius Error: Boris and Vay superior in most cases
   • Cross-field Drift: Vay and HC generally superior to Boris, Vay slightly better
   • Energy Conservation: All methods maintain good conservation in pure magnetic field
3. Relativistic Effects:
   • Larger γ makes relative phase error differences more pronounced
   • UR regime shows more obvious spurious force effects in Boris and Vay
   • HC's volume preservation advantage evident in HR and UR regimes
4. Error Evolution:
   • Phase error grows linearly
   • Gyroradius error exhibits periodic oscillation
   • γ error shows staggered pattern in UR regime

Related Work

Main Research Directions

1. Explicit Integrators:
   • Boris (1971): Classical method, widely applied
   • Vay (2008): Improved for relativistic cases
   • Higuera-Cary (2017): Preserves structural properties
2. Implicit Integrators:
   • Cohen et al. (1982): Larger steps, fewer errors
   • Pétri (2017): Fully implicit numerical integration
   • Zhang et al. (2024): Preserves volume, energy, and Lorentz invariance in strong magnetic fields
3. High-order Methods:
   • Winkel et al. (2015): High-order Boris integrator
   • Iterative methods achieving arbitrary order accuracy
4. Comparative Studies:
   • Ripperda et al. (2017): Comprehensive ultrarelativistic regime comparison
   • Vay (2008), Higuera-Cary (2017): Scenario-specific comparisons
   • This paper: First coverage of LR, HR, UR three regimes

Advantages of This Work

1. Coverage Range: Comprehensive analysis from γ=1.27 to γ=56.95
2. Measured Data: Actual runtime on real hardware, not just theoretical analysis
3. Unified Standard: Proposes fitness parameter f as selection criterion
4. Practical Orientation: Clear method recommendations for different scenarios

Conclusions and Discussion

Main Conclusions

1. Computational Efficiency Ranking: Boris > HC ≈ Vay (all scenarios)
2. Application Recommendations:
   • Pure Magnetic Gyration: Boris (fastest and sufficiently accurate)
   • Cross-field Drift: Vay (optimal accuracy) or HC (volume preservation)
   • Computationally Constrained Scenarios: Boris (efficiency priority)
   • Accuracy Priority Scenarios: HC (small phase error) or Vay (precise drift)
3. HC Method Evaluation:
   • Excellent phase error performance
   • Higher computational cost, other errors not necessarily superior
   • Volume preservation property valuable in specific applications
4. Fitness Parameter:
   • Successfully quantifies "cost-effectiveness"
   • Boris advantage clear in pure magnetic field (f=68.35 vs 19.91/43.17)
   • Boris still optimal in electromagnetic field but gap narrowed
5. No "Universal" Method: Different application scenarios require different choices; this research provides decision basis

Limitations

1. Test Range:
   • Single-particle dynamics only
   • Uniform fields, no field gradients or time-varying effects
   • Fixed time step, no adaptive step testing
2. Scenario Constraints:
   • No extreme nonlinear cases tested
   • No particle-particle interactions considered
   • Lack of systematic testing in actual PIC simulations
3. Fitness Parameter:
   • Energy error as sole metric potentially incomplete
   • Weight choices (1/κ and e^(-ε)) lack rigorous theoretical basis
   • Universal applicability across applications needs more verification
4. Implementation Dependence:
   • Runtime depends on specific implementation and hardware
   • Compiler optimization may alter relative performance
   • GPU accelerators not tested

Future Directions

1. Adaptive Time Stepping: Study symplectic-preserving adaptive algorithms and performance
2. Multi-particle Systems: Validate conclusions in actual PIC codes
3. More Integrators: Include implicit and high-order methods in comparison
4. Fitness Parameter Optimization:
   • Explore different weighting schemes
   • Multi-metric synthesis (beyond energy)
   • Application-specific customization
5. Extreme Condition Testing:
   • Strongly time-varying fields
   • Extreme relativistic regime (γ > 1000)
   • Strong nonlinear effects

In-depth Evaluation

Strengths

1. Strong Systematicity:
   • First coverage of complete LR, HR, UR regimes
   • Multi-dimensional error analysis (phase, radius, γ, drift)
   • Combines theoretical analysis with measured data
2. High Practical Value:
   • Provides clear application recommendations
   • Fitness parameter f enables quick decision-making
   • Measured runtime valuable for practical applications
3. Technical Rigor:
   • Detailed FLOPs calculation (Table 1)
   • Multiple run averaging eliminates noise
   • Coarse grid amplifies errors for analysis
4. Clear Presentation:
   • Complete mathematical derivations
   • Rich figures and tables (18 figures, 5 tables)
   • Clear logical structure
5. Open-source Potential: Self-developed PaTriC code can provide tools for community

Shortcomings

1. Theoretical Depth:
   • Lacks theoretical analysis of error propagation
   • Fitness parameter f mathematical foundation insufficient
   • No exploration of symplecticity-error relationship
2. Experimental Limitations:
   • Large gap between single-particle simulation and actual PIC
   • Uniform field assumption overly idealized
   • Limited test cases (2 scenarios per regime)
3. Comparison Gaps:
   • Original Boris form (tan form) not included
   • No implicit method comparison
   • Boris high-order variants not tested
4. Statistical Analysis:
   • Lacks statistical significance testing of errors
   • 10 runs may insufficient for complete performance characterization
   • Standard deviations or confidence intervals not reported
5. Practical Details:
   • No discussion of implementation in actual PIC codes
   • Lacks systematic parameter sweeps over field strength and particle energy
   • Multi-particle extrapolation (Section 3.3) oversimplified

Impact

1. Academic Contribution:
   • Fills systematic research gap in LR and HR regimes
   • Fitness parameter concept may be adopted by subsequent research
   • Provides quantitative criteria for integrator selection
2. Practical Value:
   • Plasma simulation researchers can directly apply recommendations
   • Guides PIC code development
   • Particularly valuable in computationally constrained scenarios
3. Reproducibility:
   • Detailed method description
   • Complete parameter specification
   • Code open-sourcing would greatly enhance reproducibility
4. Limitations:
   • Single-particle conclusions to multi-particle generalization needs verification
   • Hardware and implementation-specific dependencies
   • Fitness parameter universal applicability needs testing

Applicable Scenarios

Strongly Recommended:

• Plasma acceleration simulations (positron acceleration, etc.)
• Magnetic confinement fusion simulations
• Astrophysical relativistic jets
• Teaching and prototype development

Applicable but Requires Caution:

• Large-scale PIC simulations (needs further verification)
• Adaptive time step scenarios (requires modification)
• Strongly time-varying or non-uniform fields (extrapolation validity unknown)

Not Suitable:

• Scenarios requiring long-term energy conservation (consider implicit methods)
• Extreme relativistic regime (γ > 100, beyond test range)
• Scenarios with strict volume preservation requirements (HC potentially superior but insufficiently evaluated)

References

Key References

1. Boris (1971): Original Boris method proposal
2. Vay (2008): Physics of Plasmas 15:056701 - Vay integrator
3. Higuera & Cary (2017): ICOPS 2017 - HC integrator
4. Ripperda et al. (2017): ApJS 235 - Comprehensive ultrarelativistic regime comparison
5. Qin et al. (2013): Physics of Plasmas 20:084503 - Boris method advantages analysis
6. Zenitani & Umeda (2018): Physics of Plasmas 25(11) - Boris simplified form analysis

Overall Assessment: This is a high-quality, practically-oriented comparative research paper with strong systematicity and clear value to the plasma simulation community. The fitness parameter proposal demonstrates innovation, though theoretical foundations could be strengthened. Main shortcomings lie in the gap between single-particle simulation and practical applications, and lack of deeper theoretical analysis. Recommended future work includes validating conclusions in actual PIC codes and open-sourcing the PaTriC code to enhance impact.