Boris Pusher

The native Boris pusher is a widely used algorithm for tracing charged particles in magnetic fields. It is second order, volume preserving, and exactly conserves energy in a static, uniform magnetic field, while being highly efficient for long-time tracing. A basic usage example is shown in Boris Method; this page covers the full family of Boris solvers available in TestParticle.jl.

1. Standard Boris Method

The Boris method advances the particle in a staggered (leapfrog) manner: the velocity is first rotated by the magnetic field, then the position is pushed using the updated velocity, and the process repeats. Because the magnetic rotation is applied as an exact rotation rather than an approximate update, energy is conserved to machine precision in a static, uniform field.

The method requires a fixed time step $\mathrm{\Delta }t$ much smaller than the gyroperiod $T=2\pi /\mathrm{\Omega }$ to accurately capture the orbital phase. The phase error scales as $\mathcal{O}(\mathrm{\Delta }{t}^2)$, so a finer step is needed when phase fidelity matters more than wall-clock time.

To reach extremely high fidelity over long time scales without shrinking $\mathrm{\Delta }t$ drastically, TestParticle.jl implements several advanced variants: Multistep Boris, Hyper Boris, and Adaptive Boris.

2. Multistep Boris ($n$-step Boris)

The Multistep Boris integrator artificially subdivides the standard update cycle into $n$ smaller sub-steps. The electric and magnetic fields are evaluated once at the full time step location, but the velocity rotation is broken into $n$ smaller partial rotations.

Let the normalized rotation vector for the standard Boris step be:

$$\mathbf{t}=\frac{q\mathrm{\Delta }t}{2m}\mathbf{B}$$

In the Multistep Boris scheme, the rotation vector is divided by $n$:

$${\mathbf{t}}_n=\frac{1}{n}\mathbf{t}$$

The rotation part of the Boris method is recursively applied $n$ times using ${\mathbf{t}}_n$. While reducing the numerical detuning and improving rotation accuracy, it does not intrinsically change the mathematical order of the phase error (it remains 2nd order) but substantially decreases the absolute error coefficient.

3. Hyper Boris (Gyrophase Correction)

The Hyper Boris integrator (Zenitani & Kato, 2025) achieves higher-order accuracy in tracking the precise gyro-orbital phase by introducing specialized correction factors to the electric and magnetic drift terms.

Given an order $N$ ($N\in \{4,6\}$), the normalized vectors are modified prior to the rotation:

$$\begin{array}{rl}{\mathbf{t}}_n^{\prime }& =f_N{\mathbf{t}}_n\\ {\mathbf{e}}_n^{\prime }& =f_N{\mathbf{e}}_n+c_N({\mathbf{e}}_n\cdot {\mathbf{t}}_n){\mathbf{t}}_n\end{array}$$

where ${\mathbf{e}}_n=\frac{q\mathrm{\Delta }t}{2mn}\mathbf{E}$.

The factors $f_N$ and $c_N$ are Taylor-expanded coefficients based on the magnitude $t_{mag}^2=|{\mathbf{t}}_n{|}^2$:

4th-Order Hyper Boris ($N=4$):

$$\begin{array}{rl}{f}_4& =1+\frac{1}{3}{t}_{mag}^2\\ c_4& =-\frac{1}{3}\end{array}$$

6th-Order Hyper Boris ($N=6$):

$$\begin{array}{rl}{f}_6& =1+\frac{1}{3}{t}_{mag}^2+\frac{2}{15}{t}_{mag}^4\\ c_6& =-\frac{1}{3}-\frac{2}{15}{t}_{mag}^2\end{array}$$

These analytically tuned correctors virtually eliminate phase disparity and energy fluctuation across standard $\mathrm{\Delta }t$ bounds.

4. Adaptive Boris Method

The AdaptiveBoris solver adjusts the time step $\mathrm{\Delta }t$ dynamically based on the local cyclotron frequency ${\mathrm{\Omega }}_c=|qB/m|$. This is particularly useful in systems with strong magnetic field gradients, such as magnetic mirrors or planetary magnetospheres, where the required resolution varies significantly along the particle's trajectory.

The time step is determined by:

$$\mathrm{\Delta }t=\eta T_c=\eta \frac{2\pi }{{\mathrm{\Omega }}_c}$$

where $\eta$ is a safety factor (typically between 0.01 and 0.1). This represents the ratio of the time step to the local gyroperiod.

Maintaining Time Reversibility and Energy Conservation

The standard Boris method is a second-order, volume-preserving integrator that exactly conserves energy in a static, uniform magnetic field. These properties are closely linked to its time-reversibility. However, naive adaptive time-stepping usually breaks this reversibility because it disrupts the staggered "leapfrog" synchronization between position and velocity.

To preserve the leapfrog structure and maintain energy conservation, TestParticle.jl employs a velocity resync procedure whenever the time step is updated. This technique involves moving the velocity back to the node (position location) and then repositioning it based on the new time step. Given a step from ${\mathbf{x}}_n$ to ${\mathbf{x}}_{n+1}$ with $\mathrm{\Delta }{t}_{old}$:

1. Advance position: The position is updated to ${\mathbf{x}}_{n+1}$ using the velocity ${\mathbf{v}}_{n+1/2}$ and $\mathrm{\Delta }{t}_{old}$.

2. Update time step: A new $\mathrm{\Delta }{t}_{new}$ is calculated based on the magnetic field at the new position ${\mathbf{x}}_{n+1}$.

3. Resynchronize velocity: The velocity ${\mathbf{v}}_{n+1/2}^{old}$ (centered at $t_{n+1}-\frac{1}{2}\mathrm{\Delta }{t}_{old}$) is moved to the node $t_{n+1}$ by a half-step Boris push, and then moved to a new half-step location $t_{n+1}-\frac{1}{2}\mathrm{\Delta }{t}_{new}$ by a half-step backward push.

This ensures that the velocity is always correctly centered relative to the current $\mathrm{\Delta }t$ before each step. This "re-centering" at the nodes allows the integrator to remain practically time-reversible and maintains excellent energy conservation even as the time step changes by orders of magnitude.

5. Using the Solvers

Recommended default

If you are unsure which member of the family to pick, start with MultistepBoris6(; n = 4), the sixth-order Hyper-Boris scheme with four sub-cycles. Over 300 gyroperiods at $\mathrm{\Delta }t=T_{\mathrm{gyro}}/20$ its gyrophase error stays at $3.7\times {10}^{-7}$ rad, where the standard Boris() has already saturated at $2.7$ rad, for about $2.8\times$ the cost of a standard step and with the energy error still at round-off. The extra sub-cycles are nearly free: the fields are evaluated only once per step, so raising $n$ tightens the absolute phase error without added field work. Reach for plain Boris() only when the cheapest possible step is the priority, and see AdaptiveBoris below when the local gyroperiod varies strongly along the trajectory.

Fixed-step Solvers

Fixed-step solvers are specified as the second argument to solve (after the problem) and require a dt keyword:

• Boris(): The standard Boris pusher.

• MultistepBoris2(n): Sub-cycling division count n.

• MultistepBoris4(n): 4th-order Hyper-Boris with sub-cycling n.

• MultistepBoris6(n): 6th-order Hyper-Boris with sub-cycling n.

# Standard Boris
sol = solve(prob, Boris(); dt)

# Multistep Boris (n=2)
sol = solve(prob, MultistepBoris2(n=2); dt)

# Hyper Boris (N=4, n=2)
sol = solve(prob, MultistepBoris4(n=2); dt)

Combining both $n>1$ and higher-order correction ($N>2$) ensures ultra-high stability tracking over drastically varying gradient fields.

One trajectory returns an ODESolution, the same type any other SciML solver returns.

Adaptive Boris

The adaptive solver adjusts the time step automatically based on the local gyroperiod.

# Adaptive Boris with safety factor 0.05 (20 steps per period)
alg = AdaptiveBoris(safety=0.05)
sol = solve(prob, alg)

Ensembles

Multiple particles are traced through a SciML EnsembleProblem, so the whole ensemble interface applies: EnsembleThreads(), EnsembleDistributed(), EnsembleSplitThreads(), trajectories, seed, batch_size, pmap_batch_size, and the output_func and reduction hooks.

eprob = EnsembleProblem(prob; prob_func, safetycopy = false)
sols = solve(eprob, Boris(), EnsembleThreads();
    dt, trajectories = 1000, seed = 1234)

prob_func(prob, ctx) prepares each trajectory. It receives an EnsembleContext whose ctx.rng is derived from seed, which keeps a run reproducible. Because the prob_func of a TraceProblem is used to build the same ensemble, the shorter form is equivalent:

prob = TraceProblem(stateinit, tspan, param; prob_func)
sols = solve(prob, Boris(), EnsembleThreads();
    dt, trajectories = 1000)

sols.u is then a vector of ODESolutions, one per trajectory.

The ensemble algorithm is a required third argument in the shorter form, because solve(prob, alg) without one means a single trajectory. It cannot have a default either: SciML calls exactly that two-argument form once per trajectory it builds, so a default would ask every trajectory to build an ensemble of its own.

Bulk Initial States and Raw Output

On the kernel backend, solve(prob, Boris(), backend; ...) or solve(prob, Boris(), EnsembleKernel(backend); ...), two further keywords are accepted:

• u0: a (trajectories, 6) matrix of initial states, one row per particle. It is copied in a single transfer and replaces prob_func, which then never runs, so it works on a problem built without a prob_func; if both are given, u0 wins. The element type is converted to that of prob.u0, and a wrong number of rows is an ArgumentError.

• raw_output = true: return (; u, t, terminated) instead of an EnsembleSolution, where u is a (trajectories, 6, nout) array of the saved states, t the matching times, and terminated a flag per trajectory marking the ones that stopped early because they left the domain. It applies to every problem, whether the initial states come from prob_func or from u0, and it skips building one ODESolution per particle, which is what makes large ensembles affordable.

states = Matrix{Float64}(undef, 100_000, 6)   # one initial state per row
out = TestParticle.solve(prob, Boris(), CPU();
    dt, trajectories = 100_000, u0 = states, raw_output = true)

Both keywords belong to the kernel path only. EnsembleThreads(), EnsembleSerial() and the other SciML ensemble algorithms accept neither; passing u0 to them is an error. See GPU Ensemble Tracing for the backend setup.

Saving the Output

Every accepted step is saved by default. To choose the output times explicitly, pass saveat, either as a collection of times or as an interval:

# Save at the given times
sol = solve(prob, Boris(); dt, saveat = 0.0:5.0e-10:3.0e-8)

# Save every 5.0e-10 across the time span
sol = solve(prob, Boris(); dt, saveat = 5.0e-10)

saveat does not shorten any step, so the trajectory is identical to a run without it and only the reporting changes. Inside a step the state is interpolated linearly, which is the dense output these methods admit: a Boris step advances the position linearly with the half-step velocity ${\mathbf{v}}_{n+1/2}$, so the interpolant reproduces the stored positions exactly, whereas the velocity is only as accurate as the method itself. Requesting a value inside a step is therefore consistent with the solver's own order, not better or worse than the step values around it.

save_start and save_end (both true by default) add the ends of the time span to the requested times. save_fields = true and save_work = true keep appending their columns to every saved state.

Removed keyword: savestepinterval

savestepinterval = k, which saved every $k$-th step regardless of how the times were named, has been removed. Use saveat = k * dt to report the same times in a fixed-step run. Passing savestepinterval is now an ordinary unsupported-keyword error.

6. Where the solvers live

The Boris family is implemented in its own package, BorisPushers, which follows the SciML convention for an algorithm package: the methods are ordinary SciML algorithms, driven by the SciML loop, and they can be used on their own. TestParticle.jl depends on it and adds the parts that are specific to particle tracing: the TraceProblem container, the boundary callback, and the field and work columns.

Used directly, the solvers only need a parameter container that answers three questions, get_q2m, get_EField and get_BField. The default methods assume the layout (q2m, m, E, B, ...) that prepare produces, and any other type can support the solvers by adding methods to those three functions:

using BorisPushers, StaticArrays

param = (q2m, m, Efunc, Bfunc)   # or any type with the three accessors
prob = ODEProblem((u, p, t) -> nothing, SA[0.0, 0.0, 0.0, 1.0e5, 0.0, 0.0], tspan, param)
sol = solve(prob, AdaptiveBoris(safety = 0.1); dt)

An ODEProblem needs a right-hand side f(u, p, t), but a Boris method is a map from one state to the next, not a differential equation, so f is never called: at every step the method takes the charge-to-mass ratio and the fields from p, and (u, p, t) -> nothing is enough. Adaptivity is selected by the adaptive solve keyword, as for any other SciML solver, and the adaptive methods follow the local gyroperiod rather than an error estimate, since a Boris step forms none.

What TestParticle.jl still owns is the layer around the loop, in src/boris/boris_solve.jl: the conversion of a TraceProblem into an ODEProblem, the gyroperiod-based initial step of the adaptive methods, the isoutside callback, and the appended field and work columns.

To maximize performance, fixed-step Boris solves on standard domains automatically take a specialized fast path inside solve!. This eliminates generic SciML integrator loop bookkeeping and preallocates trajectory arrays, achieving raw native performance. Whenever adaptive stepping (AdaptiveBoris), boundary checking (isoutside), intermediate save targets (saveat), or custom SciML callbacks are present, solve! routes transparently through the full SciML integrator loop.

7. GPU Portability with KernelAbstractions

The Boris algorithms are fully portable across CPU and GPU architectures. SciML integrator loops typically rely on mutable integrator caches that cannot be executed within device kernels. To achieve backend-agnostic execution without code duplication, BorisPushers provides a stateless step API that separates the physics and stepping math from the integrator bookkeeping:

• advance_boris(v_half, r, dt, t, p, alg): advances the staggered position r and half-step velocity v_half by a time step dt.

• update_velocity_half(v, r, dt, t, p, alg): shifts the node velocity back by dt / 2 to initialize the leapfrog staggering.

• update_velocity_node(v_half, r, dt, t, p, alg): reconstructs the synchronous velocity at the node for saved states.

Because all algorithm types in AbstractBoris are isbits structs, they can be passed directly into GPU kernels.

Switching Between CPU and Device Backends

TestParticle.jl routes execution based on the chosen solver interface and backend:

• Host CPU solves: Standard calls to solve(prob, alg; ...) or solve(prob, alg, EnsembleThreads(); ...) route through BorisPushers. This path handles single-particle trajectories, adaptive stepping (AdaptiveBoris), and traces requiring SciML callbacks.

• Device ensemble solves: Passing a KernelAbstractions backend to solve(prob, alg, backend; trajectories, ...) or using the SciML ensemble interface solve(prob, alg, EnsembleKernel(backend); trajectories, ...) routes execution through the native KernelAbstractions driver. The device kernels call the exact same stateless advance_boris and update_velocity_half functions from BorisPushers, ensuring identical numerical behavior between CPU and GPU.

Boundary Checks on Device

The isoutside keyword is honored on the device path as well. A step whose end state leaves the domain is rolled back and the trace stops at the last state inside it, the same convention the native CPU solvers follow, and the solution is reported with ReturnCode.Terminated. Since t is shared by the whole ensemble, raw_output = true cannot shorten it per trajectory: a terminated trace instead repeats its last state in the slots it did not reach, and the returned terminated flags mark which trajectories stopped early.

SciML callbacks cannot run inside a device kernel, so they remain a CPU-only feature: solve(prob, alg) and the host ensemble algorithms route through the integrator loop, which handles them.

The check is free when the state is reported at every step, because the node velocity it tests is then built anyway. With sparse output it costs one extra node reconstruction, hence one extra field evaluation, per step: measured at about 8% of a 100,000-particle ensemble over 20 steps of an interpolated field. The check itself is skipped entirely when no isoutside is given, so a solve that does not use it pays nothing for the possibility.

On a device the predicate is passed to the kernel as an argument, so it has to be isbits, like every other kernel argument. A check written over the state alone, or over numbers it closes over, is fine; one that closes over an interpolated field or an array held on the device is not, and belongs on the host path instead.

Algorithm and Field Portability

• Supported Solvers: All fixed-step solvers—including standard Boris() and the higher-order Hyper-Boris variants MultistepBoris{N}—run natively on device backends (such as CUDA, AMDGPU, oneAPI, Metal, or multi-threaded CPU).

• Adaptive Solvers: AdaptiveBoris dynamically adjusts the time step per particle based on the local gyrofrequency. Because non-uniform per-particle time stepping diverges across SIMD threads on device hardware, adaptive methods run exclusively on the host CPU through the SciML loop.

• Fields: Both analytic and discrete grid interpolations are supported. Interpolated fields are automatically transferred to device memory via Adapt.jl.

For detailed usage, backend setup, and large-scale ensemble benchmarks, see GPU Ensemble Tracing.