Numerical Methods

Efficient and accurate particle tracing requires a good understanding of how to set up initial conditions and choose the right numerical algorithms.

Setting Initial Conditions

The initial state of a particle is typically a 6-element vector or an SVector representing its position $\mathbf{x}$ and its velocity $\mathbf{v}$ (or a related quantity for relativistic cases).

Particle Types

TestParticle.jl provides several predefined species that can be passed to the prepare function using the species keyword argument:

• Proton: $m\approx 1.67\times {10}^{-27}$ kg, $q\approx 1.60\times {10}^{-19}$ C

• Electron: $m\approx 9.11\times {10}^{-31}$ kg, $q\approx -1.60\times {10}^{-19}$ C

• Ion: A generic ion type.

You can also define custom species using the Species(m, q) struct.

Non-relativistic Case

For non-relativistic simulations, the initial state is [x, y, z, vx, vy, vz].

x0 = [1.0, 0.0, 0.0]
v0 = [100.0, 0.0, 0.0]
stateinit = [x0..., v0...]

Relativistic Case

In relativistic simulations, the state is typically represented by its position $\mathbf{x}$ and the quantity $\gamma \mathbf{v}$, where $\gamma$ is the Lorentz factor. TestParticle.jl handles the conversion to velocity internally using the get_relativistic_v function.

x0 = [1.0, 0.0, 0.0]
v0 = [0.8c, 0.0, 0.0] # 80% of speed of light
# For relativistic equations, the state is (x, \gamma v)
gamma = 1 / sqrt(1 - (norm(v0)/c)^2)
stateinit = [x0..., (gamma .* v0)...]

Relativistic velocity

When using trace_relativistic!, the solver expects the last three elements of the state to be $\gamma \mathbf{v}$.

Choice of numerical algorithms

By default DifferentialEquations.jl applies Tsit5 to an ODE problem. If only OrdinaryDiffEq.jl is imported, then we need to explicitly select a solver.

However, not all solvers are guaranteed to work. For example, some cases with strong magnetic fields are tested to work only with fixed timestep algorithms like Euler and the Adams-Bashforth family.

Currently we recommend Vern9 as a starting point for adaptive timestepping, with additional fine tuning by changing reltol if needed.

Native Boris solver

You can also try out the native implementation of the Boris method in TestParticle.jl. The Boris pusher is specifically designed for particle tracing in magnetic fields and has the property of exactly conserving energy in a static, uniform magnetic field. It follows a similar interface for ease of adoption:

dt = 3e-11 # fixed time step
savestepinterval = 10
prob = TraceProblem(stateinit, tspan, param)
sol = TestParticle.solve(prob; dt, savestepinterval)[1]

Boundary check

Boundary check is performed in each iteration for the particle pusher via a specified isoutofdomain function if provided by the user. By default, no boundary check is performed, but the tracing may stop if NaN values are encountered when interpolating the EM fields or numerical instability occurs in the SciML solvers.

Boundary Check

When using SciML solvers (e.g. Vern9), it is recommended to use DiscreteCallback for checking boundary conditions instead of the isoutofdomain keyword argument for better performance.

isoutofdomain(u, p, t) = norm(u) < Rₑ
callback = DiscreteCallback(isoutofdomain, terminate!)
sol = solve(prob, Vern9(); callback)

However, for the native Boris pusher in TestParticle.jl, isoutofdomain is still the correct keyword argument to use.

Solution Interpolations

All the tracing solutions come with an interpolation method to make them continuous. Depending on the order of the scheme, different orders of interpolations are applied.

Note that if the scheme comes with a "lazy" interpolation (e.g. Vern family), you can either output the full solution (default) and interpolate later, or turn off the lazy interpolation and select part of the solution to save.

# Default: lazy interpolation enabled
sol = solve(prob, Vern9())

# Selective saving with lazy interpolation disabled
sol = solve(prob, Vern9(lazy=false), save_idxs=[1,2,3])

The interpolated solution would be wrong if lazy=true and save_idxs is not the full state vector [1,2,3,4,5,6].