Numerical methods

The fields depend on three spatial coordinates and can be represented on a grid. The distribution function f(x, v, t) depends on six phase-space coordinates: a grid with 64 points per dimension would have 64⁶ ≈ 7·10¹⁰ cells for each species. Struphy therefore represents the fields by finite element coefficients and the distribution function by a finite set of weighted markers.

A kinetic time step then consists of three kinds of operation: updating the fields, pushing the markers, and transferring quantities between the two representations. Whether the coupled system conserves energy and charge depends largely on how this transfer is done.

The field equations rely on the identities ∇·(∇×A) = 0 and ∇×∇φ = 0. If a discretization satisfies them only up to truncation error, a magnetic field obtained as a discrete curl has a small non-zero divergence (a spurious magnetic charge), and a discrete gradient field has a small non-zero curl (a spurious electric field). Both can act on the plasma as non-physical forces.

Such errors accumulate over a run and can show up as artificial heating, spurious drifts or instabilities without a physical energy source. Common remedies, such as divergence cleaning, filtering or periodic projection, reduce the error after it has been introduced.

Struphy instead uses finite element exterior calculus (FEEC). Each field is assigned to one of four B-spline spaces forming a discrete de Rham complex: potentials, fields with a curl, fields with a divergence, and densities. The spaces are connected by discrete gradient, curl and divergence operators whose successive compositions vanish up to round-off. A magnetic field that is initially divergence-free therefore remains so, since its time evolution is given by a discrete curl.

FEEC basics shows the four spaces and which model variable sits in each.

A marker is a sample point in phase space that carries a weight; it does not represent a single physical particle. The distribution function is approximated by a sum over N weighted markers, and its moments (density, current) are obtained as Monte Carlo estimates. A simulation with 10⁵ markers thus samples a plasma with a much larger number of physical particles.

The number of markers determines the statistical error, which decreases as 1/√N for pseudo-random sampling. Markers can alternatively be drawn from low-discrepancy (Sobol) sequences, which can converge faster, or placed on a regular tesselation. The initial distribution is specified in physical terms, as a Maxwellian with given density and temperature, optionally with a drift or with different parallel and perpendicular temperatures.

For small perturbations of a known equilibrium f₀, the δf method reduces the noise: only the deviation f − f₀ is evolved, and f₀ serves as a control variate. The sampling noise then scales with the amplitude of the perturbation rather than with the full distribution, so that linear growth or damping rates can be measured with moderate marker numbers.

Kinetic backgrounds are the f₀ you subtract — the distributions themselves, with their velocity-space shapes.

Markers are pushed by the fields, and the fields are driven by currents and densities computed from the markers. If the interpolation of the fields to the marker positions and the deposition of marker quantities onto the grid are not adjoint to each other, the coupled system does not conserve energy.

Struphy uses the same spline basis for both operations. Pushing a marker evaluates the field splines at its position. The current or density that drives the fields is obtained by evaluating the same splines at all marker positions and summing over the markers. This is a Monte Carlo quadrature of the phase-space integrals in the weak form of the coupled equations, with the markers as quadrature points.

Since sampling noise is largest at short wavelengths, the deposited quantities can optionally be filtered before they enter the field solve, either with a local three-point stencil or by retaining selected toroidal mode numbers.

The 6 coupling propagators are the steps that do this, across 13 kinetic and hybrid models — see the propagator catalogue.

A time step is split into substeps (propagators), each corresponding to a distinct physical process: for example the rotation of velocities about the magnetic field, the streaming of markers along their velocity, the acceleration by the electric field, or the momentum exchange between an energetic species and a bulk fluid. Each substep is solved with a method suited to its structure. The rotation about the magnetic field, for instance, changes the direction but not the magnitude of the velocity and therefore leaves the kinetic energy unchanged for any step size.

Composing the substeps sequentially (Lie–Trotter) is first-order accurate in the step size. A symmetric composition (Strang), with half-steps, a full step and the half-steps in reverse order, is second-order accurate. For long runs, the energy error of a symmetric composition typically oscillates within a bounded range instead of drifting in one direction.

Time integration shows the splitting for any model, and the Runge–Kutta tableaux behind the explicit substeps.

Each model records its energies as scalar time series: field energy, kinetic energy, thermal energy and their sum. Conservation can be assessed by plotting the total energy, and the individual contributions show which of them changes.

In the electrostatic models, Gauss's law is imposed at t = 0 by solving a Poisson problem for the initial electric field; at later times it follows from the discrete scheme. The Vlasov–Maxwell model can record the residual as gauss_error, which should remain at round-off level. A growing gauss_error indicates that the discrete charge continuity equation is not satisfied.

Several models are additionally verified against known analytic results. The examples list the measured value next to the analytic one: a damping rate within about 1% of the exact −0.1533 for weak Landau damping, and a light-wave phase velocity within about 1% of c, both on grids coarse enough to run in seconds.

The examples run these end to end, each with its measured result next to the analytic one.