ORB5: a global electromagnetic gyrokinetic code using the PIC approach in toroidal geometry

1. Introduction

• Microinstabilities are closely linked to tokamak confinement, with turbulence governing much of the heat and particle loss; energetic particles can destabilize Alfvén eigenmodes that deteriorate confinement.The paper situates its code update in the study of these performance-limiting phenomena.Source: §1 · p. 3
• Gyrokinetic theory reduces the phase space from 6D to 5D by removing fast gyro-angle dependence while retaining a kinetic description.Its use is motivated by the low frequency and small core amplitude of drift-wave microinstabilities and Alfvén eigenmodes relative to ion gyrofrequency and background quantities.Source: §1 · p. 3
• ORB5 solves the nonlinear global gyrokinetic Vlasov–Maxwell system with collisions and sources, and the paper reviews its physical and numerical developments.The code uses a Lagrangian PIC scheme; the paper also presents its parallel efficiency and illustrative simulations.Source: §1 · p. 3

2.1. Magnetic geometry, coordinate system, and normalization

• ORB5 uses ideal-MHD equilibria from CHEASE or an analytical equilibrium with circular concentric magnetic surfaces.The ideal-MHD equilibria are solutions of the Grad–Shafranov equation.Source: §2.1 · p. 3
• The code uses straight-field-line coordinates, labeling magnetic surfaces by s = ψ/ψedge and using ϕ for toroidal angle.The poloidal angle is defined using the safety-factor profile q(s) and geometric poloidal angle θ.Source: §2.1 · p. 4
• ORB5 normalizes quantities using ion mass and charge, magnetic field at the axis, and electron temperature at a reference surface.For multiple ion species, the user must select a reference ion species; derived units include inverse ion cyclotron frequency, ion sound speed, ion sound Larmor radius, and volume-averaged density.Source: §2.1 · p. 4

2.2. GK Equations for fields and particles

• The variational formulation permits approximations to be incorporated before deriving equations, yields conserved quantities for diagnostics, and provides weak forms suited to finite elements.The model includes first-order geometric corrections and electromagnetic corrections through second order under its stated ordering.Source: §2.2 · p. 4
• The ordering distinguishes gyrokinetic models with ϵ⊥ ∼ O(1) from drift-kinetic models with ϵ⊥ ≪ 1; both are implemented in ORB5.It treats background magnetic-geometry effects as one order smaller than relative electromagnetic fluctuations.Source: §2.2 · p. 4

2.2.2. Nonlinear gyrokinetic Vlasov equation

• The model linearizes the gyrokinetic Poisson and Ampère equations by using equilibrium distributions invariant under the unperturbed Hamiltonian, and excludes some nonlinear terms from particle characteristics.The neoclassical limit neglects electromagnetic fields and assumes small banana widths relative to characteristic system lengths.Source: §2.2 · p. 5, §2.2.2 · p. 7

2.2.1. Quasineutrality and Amp`ere equations

• The weak-form quasineutrality equation combines ion gyro-charge and polarization terms with electron drift-kinetic charge; drift-kinetic electrons contribute no linear polarization density.The Ampère equation is also derived from the variational principle for test functions.Source: §2.2.1 · p. 6

2.2.2. Nonlinear gyrokinetic Vlasov equation

• The Vlasov equation is reconstructed from linearized gyrocenter characteristics, with H = H0 + ϵδH1 defining the dynamics.The characteristics expose parallel, diamagnetic, magnetic-gradient, curvature, E × B, and A1∥-associated drifts; linear and neoclassical limits are also available.Source: §2.2.2 · pp. 6–7

2.3. Variants of the physical models

• Model variants generally change Hamiltonians H0, H1, and H2; the adiabatic-electron option additionally couples gyrokinetic equations to a fluid electron polarization density.The section also notes that the available approximations are not necessarily mutually exclusive.Source: §2.3 · p. 8

2.3.1. Long wavelength approximation

• The long-wavelength approximation replaces the nonlinear Hamiltonian with a second-order FLR expansion and changes only the polarization-charge term in quasineutrality.The Ampère equation and Vlasov equation remain unchanged by this approximation.Source: §2.3.1 · p. 8

2.3.3. Adiabatic electron model

• The adiabatic-electron model is electrostatic: electron response enters through a fluid contribution to the field term, while the Vlasov equation is for ions.Its quasineutrality equation modifies the gyro-charge term, and the ion dynamics can use full-FLR, Padé-approximated, or second-order long-wavelength Hamiltonians.Source: §2.3.3 · p. 9

2.3.4. Hybrid electron model

• The hybrid model treats trapped electrons kinetically; its update includes the passing-electron n = m = 0 density while keeping an adiabatic response for other passing components.The earlier model was intended to simulate linear electron modes such as TEM with a larger timestep than fully drift-kinetic electrons; in nonlinear regimes it did not ensure ambipolarity and could introduce spurious particle and momentum sources.Source: §2.3.4 · p. 10

2.3.2. Pad´e approximation

• The Padé approximation modifies the ion quasineutrality model, and the code multiplies the equation by an operator to avoid inverting it.The cited implementation note explains that the inverse of the discretized block-banded matrix is a full matrix.Source: §2.3.2 · p. 8, §2.3.5 · p. 11

2.3.5. Summary of the models

• ORB5 allows each ion polarization-density model to be combined with any electron model.The summary distinguishes ion gyrodensity, electron drift-kinetic density, and ion polarization-density contributions to quasineutrality.Source: §2.3.5 · p. 11

2.4. δf method and background distribution functions

• The δf method samples only the time-dependent perturbation, reducing statistical error relative to full-f sampling when |δf| ≪ f0.It separates the distribution into a known, time-independent f0 and a marker-sampled δf.Source: §2.4 · p. 11
• A local Maxwellian is not invariant under unperturbed dynamics, and retaining its derivative term can drive a spurious zonal flow during the linear phase.Neglecting that term, as done in many PIC codes, is inconsistent with the perturbative ordering used here.Source: §2.4 · p. 12
• A canonical Maxwellian can yield unrealistic parallel flows that prevent instability development, especially in small systems with large density and parallel-velocity gradients.Its effective density and temperature differ from the input profiles.Source: §2.4 · p. 12
• The corrected canonical Maxwellian minimizes the gap from the local Maxwellian while remaining an equilibrium distribution built from constants of motion.Its correction is zero for trapped particles and has opposite signs for forward and backward passing particles.Source: §2.4 · p. 12

2.5. Strong flows

• The strong-flow ordering permits uE/vth,i ∼ 1 and treatment of plasmas rotating toroidally near the Mach velocity.The polarization density uses a local Maxwellian, while the gyrokinetic Vlasov equation uses a canonical Maxwellian.Source: §2.5 · p. 12
• A background electrostatic potential modifies the unperturbed Hamiltonian and symplectic magnetic potential, while the field model remains electrostatic and the quasineutrality equation receives no strong-flow corrections.The gyrokinetic Vlasov equation is reconstructed from the modified characteristics.Source: §2.5 · p. 12
• For strong rotation, ORB5 uses FLOW to obtain toroidally rotating MHD equilibria consistent with collisionless kinetic equilibria in the large-system-size limit.The considered equilibria have temperature as a flux-surface function.Source: §2.5 · p. 13

2.5.1. Global gyrokinetic equilibria for rotating plasmas

• The strong-flow canonical momentum extends the conserved toroidal momentum using the toroidal component of the background E × B velocity.The conserved quantities include magnetic moment, unperturbed energy, passing-particle parallel-velocity sign, and toroidal canonical momentum.Source: §2.5.1 · p. 13
• The rotation-corrected canonical Maxwellian is chosen so flux-surface-averaged density remains close to n0 when rotation is introduced.In the local limit, this choice leads to an in-out density variation.Source: §2.5.1 · p. 13

2.6. Collisions

• Collisions are required for neoclassical transport and affect turbulence: they increase ITG turbulence through zonal-flow damping and reduce TEM turbulence through electron detrapping.Neoclassical transport is important for some particle classes, including heavy impurities.Source: §2.6 · p. 13
• ORB5 linearizes the Landau collision operator around a local Maxwellian and retains test-particle and background-reaction terms while neglecting the nonlinear perturbation term.FLR effects are neglected in the collisional dynamics.Source: §2.6 · p. 13
• The approximate background-reaction term avoids reconstructing δf and costly velocity-space integrals; combined with the test-particle term, it conserves mass, momentum, and energy and satisfies the H-theorem.The combined linear operator is self-adjoint and vanishes for shifted linearized Maxwellians.Source: §2.6 · p. 14
• The reduced electron-ion Lorentz operator models electron pitch-angle scattering and conserves mass and energy, but neglecting its background reaction means momentum conservation is not ensured.The operator assumes a large ion-to-electron mass ratio.Source: §2.6 · p. 14

2.7. Conservation laws and diagnostics

• ORB5 derives conserved quantities using Noether's method and implements them as diagnostics, including power-balance checks of simulation quality.The power balance follows from the energy conservation law associated with the model's Lagrangian.Source: §2.7 · p. 14
• The power-balance diagnostic separates particle and field energy contributions because particles advance continuously along characteristics while fields are evaluated on a fixed grid.The total electromagnetic energy is EEM = EF + Ekin.Source: §2.7 · p. 15
• The energy-balance expression is independent of the nonlinear Hamiltonian choice and remains valid for both long-wavelength and all-orders FLR polarization density.The equivalence follows because the equations of motion used to rewrite the energy come from the same field-particle Lagrangian.Source: §2.7 · p. 15
• The field energy has alternative expressions for full-FLR polarization, the long-wavelength approximation, and the Padé-approximated model.The expressions depend on the corresponding polarization model and nonlinear Hamiltonian.Source: §2.7 · p. 16
• For adiabatic electrons, their energy is counted as a field contribution, while the kinetic energy contains only the ion contribution.The adiabatic-electron energy expression uses a distinct variational formulation combining fluid and kinetic formalisms.Source: §2.7 · p. 16

3.1.2. Solving for the collisional dynamics

• Collisionless marker motion and weights use RK4, while collisions are applied afterward through a Langevin approach with randomized velocity-space kicks.The collision workflow converts weights to perturbations about a local Maxwellian, then restores their prior representation; the modified background-reaction operator conserves mass, momentum, and energy in each spatial bin.Source: §3 · p. 16, §3.1.2 · pp. 18–19

3. Numerical implementation

• ORB5 discretizes the gyrokinetic Vlasov-Maxwell system with a low-noise δf PIC method, operator splitting, RK4 collisionless integration, and Langevin collisions.The section also covers field discretization, noise control, sources, diagnostics, and parallelization.Source: §3 · p. 16

3.1. δf and equations of motion discretization

• The δf PIC scheme samples phase space with markers and evolves weights for a prescribed background and its perturbation.The standard δf method solves weight equations, whereas direct δf evaluates weights using invariance along nonlinear collisionless trajectories; collisionless runs need only the perturbed weight.Source: §3.1 · p. 17, §3.1.1 · pp. 17–18

3.1.3. Particle loading

• Markers are loaded using a Halton-Hammersley sequence, with configurable radial and velocity sampling and either white-noise or Fourier-mode weight initialization.White-noise initialization can delay physical-mode development as particle count increases; mode initialization is offered to accelerate it.Source: §3.1.3 · p. 19
• Near the magnetic axis, ORB5 changes coordinates to avoid a singularity; particles leaving the radial domain are reflected while conserving toroidal momentum, energy, and magnetic moment.Reflected particles receive zero weight to avoid unphysical perturbed-density accumulation at the edge.Source: §3.1.3 · p. 20

3.2. Quasineutrality and Amp`ere equations

• ORB5 represents fields with B-spline finite elements and solves the Galerkin quasineutrality and Ampère system in discrete Fourier space.The Fourier solve exploits axisymmetry and the concentration of relevant modes near field-aligned coefficients.Source: §3.2 · pp. 20–21

3.2.1. Gyroaveraging

• Gyroaveraging approximates the Larmor ring in the poloidal plane, using fixed or adaptive gyropoints to evaluate field quantities.Four points usually suffice for perturbations up to k⊥ρL ∼ 1; adaptive sampling smooths shorter-wavelength fluctuations.Source: §3.2.1 · p. 21

3.2.2. Fourier filter

• The field-aligned Fourier filter retains modes near m = −nq(s), removing nonphysical modes from deposited density and current to reduce noise and maximize timestep size.A typical filter width is Δm = 5.Source: §3.2.2 · p. 22

3.3. Noise control techniques

• Finite-marker sampling causes noise accumulation and weight growth, motivating controls that reduce noise without severe nonphysical artifacts.The section presents multiple noise-control methods to address these competing concerns.Source: §3.3 · p. 22

3.3.1. Krook operator

• The Krook operator weakly damps non-axisymmetric fluctuations while projecting out selected moments, including density, parallel velocity, zonal flows, and kinetic energy.When kinetic energy is not conserved, it can act as an auto-regulated heat source for temperature-gradient-driven simulations.Source: §3.3.1 · pp. 22–23

3.3.2. Coarse-graining

• Coarse-graining reduces weight-spreading and filamentation by dissipating fine-scale phase-space structure, unlike Krook damping, which can indiscriminately damp the distribution function.The authors report that Krook was unsuitable for neoclassical studies.Source: §3.3.2 · p. 23

3.3.3. Quad-tree particle-weight smoothing

• Quad-tree smoothing pairs nearby markers in velocity space and averages their weights using a distance-dependent weighting procedure.The operation conserves total weight, and typically one smoothing step is applied per timestep.Source: §3.3.3 · p. 24

3.3.4. Enhanced control variate

• The enhanced control variate mitigates the electromagnetic cancellation problem by discretizing skin terms and adiabatic currents with the same markers.The cancellation problem arises from different particle and field discretizations of large terms; production runs typically need fewer than 10 iterations.Source: §3.3.4 · pp. 25–26

3.4. Heating operators

• Heating operators control distribution-function moments rather than modeling detailed realistic heat-source physics.Thermal relaxation keeps the distribution near its initial profile without changing gyrocenter density or adding parallel momentum; fixed-power sources can represent core input and edge losses.Source: §3.4 · pp. 26–27

3.5. Parallelization

• ORB5 combines MPI domain cloning and toroidal domain decomposition, with OpenMP threading to limit clone-related communication and memory congestion.The code also has an OpenACC GPU port; clones require global reductions, while subdomains exchange guard-cell, grid, and particle data.Source: §3.5 · p. 27

4.1. Parallel scalability

• ORB5 retains more than 85% of ideal speed-up through 128 nodes, with a small superlinear gain from 2 to 16 nodes.At more than 256 nodes, MPI communications for parallel data transposes in field Fourier transforms mainly limit speed-up.Source: §4.1 · p. 28, Figure 2 · p. 28

4.2. Strong flows and toroidal rotation

• In the strong-flow simulation, temperature profiles relax as turbulence-driven transport begins, while the unconstrained parallel-flow profile relaxes slightly.The simulations include centrifugal and Coriolis drifts; with adiabatic electrons, little net momentum flux is observed, consistent with the expected nearly zero pinch-driven momentum flux.Source: §4.2 · p. 29, Figure 3 · p. 29, Figure 4 · p. 30

4.3. GK simulations of Alfvén modes in the presence of turbulence

• The simulations examine the nonlinear interaction of Alfvén modes and turbulence after prior linear benchmarking and nonlinear wave-particle studies.The case evolves three species and requires high spatial and temporal resolution, kinetic electrons, and many markers.Source: §4.3 · pp. 29, 31, 32
• Without energetic particles, heat transport has radial corrugations that are larger in the inner half of the domain, and corresponding long-lived structures appear in the E × B velocity profile.E × B velocity avalanches arise near s ~ 0.5 at the local GAM frequency and propagate outward at constant frequency; the nonzonal scalar potential grows linearly and then saturates.Source: §4.3 · p. 31, Figure 6 · p. 31, Figure 7 · p. 32
• The turbulent Alfvén-mode test finds an Alfvénic instability growing on top of turbulence when energetic particles are included.Zonal flows and geodesic acoustic modes mediate between small-scale turbulence and large-scale Alfvén modes.Source: §4.3 · pp. 29, 31, Figure 8 · p. 32

• ORB5 derives its equations from a variational Vlasov–Maxwell model that supports conservation properties for assessing simulation quality.Source: §5 · p. 32
• The code offers full-order, second-order FLR, or Padé polarization models, plus adiabatic, drift-kinetic, or hybrid electron responses.In the hybrid model, passing electrons are adiabatic and trapped electrons are drift-kinetic; the model has been corrected for its failure to respect ambipolarity.Source: §5 · p. 32
• A δf control variate, particle-weight constraints, and a Fourier filter reduce numerical noise, while hybrid parallelization scales to thousands of cores.Source: §5 · p. 33
• ORB5 agrees well with Lagrangian and Eulerian gyrokinetic codes and demonstrates simulations with strong flows, toroidal rotation, and shear Alfvén waves.Source: §5 · p. 33