How Struphy computes

The fields live in three dimensions, so you can put them on a grid. The distribution function f(x, v, t) lives in six, and you cannot. A modest 64 points per dimension is 69 billion cells in 6D — before you have resolved anything interesting, and for each species separately.

So Struphy grids the fields and samples the distribution. That split is why a kinetic plasma code has the structure it has: a field solver, a particle pusher, and — the part that decides whether the answer is physics — the machinery that lets the two exchange energy without either one inventing it.

∇·B = 0 is not a boundary condition you impose, it is a property of the magnetic field. If your discrete curl and divergence only satisfy ∇·(∇×A) = 0 to truncation error, then your simulated B acquires a small monopole charge — and a monopole exerts a force. The same goes for ∇×∇φ = 0 on the electrostatic side: a curl that does not vanish on a gradient is a spurious electric field that does work on your plasma.

These errors do not average out. They accumulate with the run, and they look like physics: anomalous heating, drifts that should not be there, instabilities with no free energy behind them. The usual remedy is to keep cleaning up after them — divergence cleaning, filtering, periodic re-projection.

Struphy instead puts each field in the spline space where its identity holds exactly: potentials in one space, fields with a meaningful curl in the next, fields with a meaningful divergence in the next, densities in the last, connected by discrete gradient, curl and divergence operators whose compositions vanish to machine precision. An initially divergence-free B stays divergence-free for the whole run because the discrete curl cannot produce divergence — not because anything corrects it.

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

The single most useful thing to know about the particle side: a marker is not an ion. It is a sample point of the distribution function, carrying a weight that says how much phase-space volume it represents. A run with 10⁵ markers is not modelling 10⁵ ions — it is doing a Monte-Carlo integral of a plasma containing 10²⁰ of them.

That framing explains the knobs. The number of markers sets the statistical error of that integral, which falls like 1/√N — so halving the noise costs four times the markers. Struphy can draw them pseudo-randomly, or from low-discrepancy (Sobol) sequences that cover phase space more evenly and beat 1/√N in practice, or place them on a regular tesselation. The starting distribution is a physical one: a Maxwellian at a given density and temperature, optionally drifting or anisotropic in the parallel and perpendicular directions.

And it explains delta-f. If the mode you are chasing is a 0.1% perturbation on a Maxwellian, sampling the whole distribution wastes nearly all your statistics resolving the equilibrium you already know analytically. Instead, subtract it: evolve only the departure f − f₀, with the known background as a control variate. The sampling noise then scales with the size of the perturbation rather than the size of the plasma, which is what makes a linear growth rate measurable at all.

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

Particles feel fields, and fields feel particles. Numerically these are opposite directions, and the temptation is to treat them as separate conveniences: interpolate the field to the particle however is cheapest, smear the particle's charge onto the grid however is cheapest. That is where energy conservation goes to die — if the two operations are not adjoint, the coupled system quietly gains or loses energy every step.

Struphy uses the same spline basis in both directions. Pushing a marker evaluates the field's splines at the marker's position. Building the current or density that drives the fields evaluates the same splines and sums over markers — which is precisely a Monte-Carlo quadrature of the phase-space integrals that appear in the coupled equations, with the markers as quadrature points. Deposition is not an approximation bolted onto the field solve; it is the field solve's integral, evaluated with the samples you have.

Because sampling noise is loudest at the shortest wavelengths — where it is least physical and most destabilizing — the deposited quantities can optionally be filtered before they reach the field solve, either with a local three-point stencil or by keeping 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 not solved in one piece. It is split into substeps, and Struphy splits along physical lines: one step rotates velocities about the magnetic field (the gyration), one streams markers along their velocity, one accelerates them in the electric field, one exchanges momentum between an energetic species and a bulk fluid. Each of those is a process you can name, and each is solved in a way that respects its own structure — the magnetic rotation, for instance, changes a particle's direction but not its speed, so it cannot heat the plasma no matter the step size.

Composing them in sequence is first-order accurate; composing them symmetrically — half-steps out, a full step, half-steps back — is second-order. The practical consequence is the one that matters for long runs: the energy error oscillates within a bounded envelope instead of drifting in one direction, so a simulation run ten times longer does not have ten times the energy error.

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

None of the above is a promise you have to take on faith. Every model tracks its energies as scalar diagnostics over the run — field energy, kinetic energy, thermal energy, and their total — so a structure-preserving claim is a curve you can plot. If the total drifts, something is wrong, and the decomposition tells you which channel is leaking.

The electrostatic models go further: Gauss's law is imposed once at t = 0 by solving a Poisson problem for the initial field, and after that it is a consequence of the discrete scheme rather than something re-imposed. The Vlasov–Maxwell model can track the residual as gauss_error, which should stay at round-off. A growing gauss_error means the charge continuity your scheme depends on is not holding.

And the physics itself is the final check. Several models carry verification tests that reproduce a known answer, and the examples put the measured number next to the analytic one: Landau damping comes out at -0.1268 against an exact -0.1533, and a light wave propagates at 0.9832 c against an exact 1 c — on a grid coarse enough to run in seconds.

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