Dimensionless Units and Normalization

This example shows how to trace charged particles in dimensionless units.

Fields from scale-free models such as MHD come in dimensionless units, so we can trace in those units and convert to any physical scale afterwards. Normalizing by the gyrofrequency also makes the grid spacing directly comparable to the gyroradius, which sets how many grid points a particle samples per gyration.

For the cosmic-ray-specific form of this normalization — the pitch angle, energy-conserving tracing, and the resolution argument — see the Cosmic Ray Propagation example.

Normalization

After normalization,     so that the gyroradius   . All the quantities are given in dimensionless units: if the magnetic field is homogeneous and the initial perpendicular velocity is 4, then the gyroradius is 4. To convert them to the original units,   and   .

The Lorentz equation in SI units is written as

It can be normalized to

with the following transformation

where all the coefficients with subscript 0 are expressed in SI units. All the variables with a prime are written in the dimensionless units.

The normalized equation is the general one with   , so tracing in dimensionless units only needs prepare(; E, B, q = 1, m = 1).

Dimensionless and Dimensional Tracing

We first solve the Lorentz equation in SI units, and then convert the quantities to normalized units and solve it again in dimensionless units.

using TestParticle, OrdinaryDiffEq, StaticArrays
using TestParticle: qᵢ, mᵢ, c
using CairoMakie

# Unit conversion factors between SI and dimensionless units
B_dim = 1.0e-8             # [T]
U_dim = c                # [m/s]
E_dim = U_dim * B_dim    # [V/m]
Ω_dim = abs(qᵢ) * B_dim / mᵢ # [1/s]
t_dim = 1 / Ω_dim        # [s]
l_dim = U_dim * t_dim    # [m]
# Electric field magnitude in SI units
Emag_dim = 1.0e-8           # [V/m]
### Solving in SI units
B_field(x) = SA[0, 0, B_dim]
E_field(x) = SA[Emag_dim, 0.0, 0.0]

# Initial conditions
x0_dim = [0.0, 0.0, 0.0] # [m]
v0_dim = [0.0, 0.01c, 0.0] # [m/s]
stateinit1 = [x0_dim..., v0_dim...]
tspan1 = (0, 2π * t_dim) # [s]

param1 = prepare(; E = E_field, B = B_field, species = Proton)
prob1 = ODEProblem(trace!, stateinit1, tspan1, param1)
sol1 = solve(prob1, Vern9(); reltol = 1.0e-4, abstol = 1.0e-6)

### Solving in dimensionless units
B_normalize(x) = SA[0, 0, B_dim / B_dim]
E_normalize(x) = SA[Emag_dim / E_dim, 0.0, 0.0]
# For full EM problems, the normalization of E and B should be done separately.
param2 = prepare(; E = E_normalize, B = B_normalize, m = 1, q = 1)
# Scale initial conditions by the conversion factors
x0_norm = x0_dim ./ l_dim
v0_norm = v0_dim ./ U_dim
tspan2 = (0, 2π)
stateinit2 = [x0_norm..., v0_norm...]

prob2 = ODEProblem(trace!, stateinit2, tspan2, param2)
sol2 = solve(prob2, Vern9(); reltol = 1.0e-4, abstol = 1.0e-6)

### Visualization
f = Figure(fontsize = 18)
ax = Axis(
    f[1, 1],
    xlabel = "x [km]",
    ylabel = "y [km]",
    aspect = DataAspect()
)

lines!(ax, sol1, idxs = (1, 2))
# Interpolate dimensionless solutions and map back to SI units
xp, yp = let trange = range(tspan2..., length = 40)
    sol2.(trange, idxs = 1) .* l_dim, sol2.(trange, idxs = 2) .* l_dim
end
lines!(ax, xp, yp, linestyle = :dashdot, linewidth = 5, color = Makie.wong_colors()[2])
invL = inv(1.0e3)
scale!(ax.scene, invL, invL)

We see that the results are almost identical, with only floating point numerical errors. Tracing in dimensionless units usually allows larger timesteps, which leads to faster computation.

Relativistic Tracing

In the relativistic case,

• Velocity is normalized by speed of light c,  ;

• Magnetic field is normalized by the characteristic magnetic field,  ;

• Electric field is normalized by   ,  ;

• Location is normalized by the  , where   , and

• Time is normalized by ,   .

Hence c = 1 in prepare; RelativisticOrbit makes the state (x, γv). At   the orbit closes after with gyroradius  , instead of after with radius 0.9 as in the non-relativistic case.

param = prepare(; B = xu -> SA[0.0, 0.0, 1.0], m = 1, q = 1, c = 1)
v = 0.9
γ = 1 / sqrt(1 - v^2)
stateinit = [0.0, 0.0, 0.0, γ * v, 0.0, 0.0]
tspan = (0.0, 2π * γ) # one gyroperiod
prob = TraceProblem{RelativisticOrbit}(stateinit, tspan, param)
sol = solve(prob, Vern9())

### Visualization
f = Figure(fontsize = 18)
ax = Axis(
    f[1, 1],
    title = "Relativistic particle trajectory",
    xlabel = "X",
    ylabel = "Y",
    aspect = DataAspect()
)

lines!(ax, sol, idxs = (1, 2))