Tracing Particles under Vacuum Assumptions

12 Jun, 2025

I've been working on creating a GPU-parallelized particle tracer for energetic particles in stellarators. I thought it'd be informative (both for myself and for others) to publish a writeup of the equations I'm using in my solver.

Writeup

Consider a particle of mass $M$ and charge $q$ , moving in equilibrium field $𝐁_{0} = \nabla \times 𝐀_{0}$ in presence of shear-Alfven wave perturbation with magnetic field $δ 𝐁 = \nabla \times α 𝐁_{0}$ . The drift kinetics of a particle with charge $q$ is described with the expanded Lagrangian (Littlejohn)

L (𝐑, \frac{d 𝐑}{d t}, v_{| |}) = q (𝐀_{0} + α 𝐁_{0} + \frac{M v_{| |}}{q B_{0}} 𝐁_{0}) \cdot \frac{d 𝐑}{d t} - \frac{M v_{| |}^{2}}{2} - μ B_{0} - q δ Φ,

where $t$ is time, $v_{| |} = 𝐯 \cdot 𝐁_{0} / B_{0}$ is the component of particle's velocity $𝐯$ along the equilibrium magnetic field, $μ = M (v^{2} - v_{| |}^{2}) / (2 B_{0})$ is the magnetic moment, and $δ Φ$ is the scalar potential of the shear Alfven wave. Observe that scalar $δ Φ$ and vector $α 𝐁_{0}$ potentials are related for transverse wave, since

𝐁_{0} \cdot (\frac{\partial α 𝐁_{0}}{\partial t} + \nabla δ Φ) = 0 \Rightarrow \nabla_{| |} δ Φ = - B_{0} \frac{\partial α}{\partial t},

where $B_{0} \nabla_{| |} = 𝐁_{0} \cdot \nabla$ is derivative along the magnetic field, and $t$ is time.

ψ = \frac{1}{2 π} \int 𝐁_{0} \cdot d 𝐒

normalized by the total toroidal flux $ψ_{0}$ , $θ$ and $ζ$ are poloidal and toroidal Boozer angles. Co- and contravariant forms of the equilibrium $𝐁_{0}$ field are, respectively,

𝐁_{0} = ψ_{0} \nabla s \times \nabla θ - ι ψ_{0} \nabla s \times \nabla ζ = K \nabla s + I (s) \nabla θ + G (s) \nabla ζ .

where $ι$ is the rotational transform. In what follows, equilibrium current $𝐣_{0}$ is assumed to be straight in Boozer coordinates, $K = 0$ .

We know that we may choose a gauge where the vector potential satisfies:

A_{0} = ψ \nabla θ - ψ_{P} \nabla ζ

We also know that

\frac{d 𝐑}{d t} = (\frac{\partial R}{\partial ψ} \cdot \dot{ψ} + \frac{\partial R}{\partial θ} \cdot \dot{θ} + \frac{\partial R}{\partial ζ} \cdot \dot{ζ})

We assume (for unperturbed) that $α = 0$ . So the Lagrangian in Boozer coordinates $(ψ, θ, ζ, v_{‖})$ is:

L (ψ, θ, ζ, v_{‖}) = q (ψ \nabla θ - ψ_{P} \nabla ζ + \frac{M v_{‖}}{q B_{0}} (K \nabla ψ + I (ψ) \nabla θ + G (ψ) \nabla ζ)) (\frac{\partial R}{\partial ψ} \cdot \dot{ψ} + \frac{\partial R}{\partial θ} \cdot \dot{θ} + \frac{\partial R}{\partial ζ} \cdot \dot{ζ}) - \frac{M v_{| |}^{2}}{2} - μ B_{0}

Simplifying,

L (s, θ, ζ, v_{‖}) = q (\nabla θ (ψ + I (ψ) \frac{M v_{‖}}{q B_{0}}) + \nabla ζ (G (ψ) \frac{M v_{‖}}{q B_{0}} - ψ_{P}) + \nabla ψ (K \frac{M v_{‖}}{q B_{0}})) (\frac{\partial R}{\partial ψ} \cdot \dot{ψ} + \frac{\partial R}{\partial θ} \cdot \dot{θ} + \frac{\partial R}{\partial ζ} \cdot \dot{ζ}) - \frac{M v_{| |}^{2}}{2} - μ B_{0}

Using the dual property of the contravariant and covariant bases, we get

L (s, θ, ζ, v_{‖}) = q (\dot{θ} (ψ + I (ψ) \frac{M v_{‖}}{q B_{0}}) + \dot{ζ} (G (ψ) \frac{M v_{‖}}{q B_{0}} - ψ_{P}) + \dot{ψ} (K \frac{M v_{‖}}{q B_{0}})) - \frac{M v_{‖}^{2}}{2} - μ B_{0} .

Now, we find the canonical momenta of $s, θ, ζ, v_{‖}$ (assuming $K = 0$ ). Immediately, we know

p_{v ‖} = 0

since $L$ has no $\dot{v_{‖}}$ . We know that

p_{ζ} = \frac{\partial L}{\partial \dot{ζ}} = q (G (ψ) \frac{M v_{‖}}{q B_{0}} - ψ_{P}) p_{θ} = \frac{\partial L}{\partial \dot{θ}} = q (ψ + I (ψ) \frac{M v_{‖}}{q B_{0}}) p_{ψ} = 0

We know that the energy is

E = \sum p \dot{q} - L,

but the first part of the expression cancels out so we get $E = \frac{M v_{‖}^{2}}{2} + μ B_{0} .$

The implicit EOMs are derived from the Euler-Lagrange equations. We derive the EOMs under a vacuum assumption (i.e. $I, K = 0$ ).

Our simplified Lagrangian is:

L = q (\dot{θ} ψ + \dot{ζ} (G (ψ) \frac{M v_{‖}}{q B_{0}} - ψ_{P})) - \frac{M v_{| |}^{2}}{2} - μ B_{0} .

and our new momenta for $θ$ is $q ψ$ .

We have the following EOM: From $ψ :$

\frac{\partial L}{\partial ψ} = \frac{d}{d t} \frac{\partial L}{\partial \dot{ψ}}

q \dot{θ} + (G^{'} (ψ) \frac{M v_{‖}}{B_{0}} - \frac{G (ψ) M v_{‖}}{B_{0}^{2}} \frac{\partial B_{0}}{\partial ψ} - q ι) \dot{ζ} - μ \frac{\partial B_{0}}{\partial ψ} = 0 .

From $θ :$

\frac{\partial L}{\partial θ} = \frac{d}{d t} \frac{\partial L}{\partial \dot{θ}}

q \dot{ψ} + \dot{ζ} \frac{M G (ψ) v_{‖}}{B_{0}^{2}} \frac{\partial B_{0}}{\partial θ} + μ \frac{\partial B_{0}}{\partial θ} = 0

From $ζ :$

\frac{\partial L}{\partial ζ} = \frac{d}{d t} \frac{\partial L}{\partial \dot{ζ}}

- \dot{ζ} G (ψ) \frac{M v_{‖}}{B_{0}^{2}} \frac{\partial B_{0}}{\partial ζ} - μ \frac{\partial B_{0}}{\partial ζ} = \frac{M}{B_{0}} (G^{'} (ψ) \dot{ψ} v_{‖} + G (ψ) {\dot{v}}_{‖}) - \frac{G (ψ) M v_{‖}}{B_{0}^{2}} (\frac{\partial B_{0}}{\partial ψ} \dot{ψ} + \frac{\partial B_{0}}{\partial θ} \dot{θ} + \frac{\partial B_{0}}{\partial ζ} \dot{ζ}) - q ι \dot{ψ}

From $v_{‖} :$

\frac{\partial L}{\partial v_{‖}} = \frac{d}{d t} \frac{\partial L}{\partial {\dot{v}}_{‖}} = 0

\frac{\partial L}{\partial v_{‖}} = G (ψ) \dot{ζ} \frac{M}{B_{0}} - M v_{‖} = 0 .

\dot{ζ} = \frac{B_{0} v_{‖}}{G}

Substituting this into our equation from $θ$ , we find that

q \dot{ψ} + \frac{M v_{‖}^{2}}{B_{0}} \frac{\partial B_{0}}{\partial θ} + μ \frac{\partial B_{0}}{\partial θ} = 0,

\dot{ψ} = - \frac{1}{q} \frac{\partial B_{0}}{\partial θ} (\frac{M v_{‖}^{2}}{B_{0}} + μ) .

Substituting $\dot{ζ}$ into our equation from $ψ$ , we find that

q \dot{θ} + \frac{M G^{'} v_{‖}^{2}}{G} - \frac{M v_{‖}^{2}}{B_{0}} \frac{\partial B_{0}}{\partial ψ} - \frac{q ι B_{0} v_{‖}}{G} - μ \frac{\partial B_{0}}{\partial ψ} = 0,

\dot{θ} = \frac{1}{q} (\frac{\partial B_{0}}{\partial ψ} (μ + \frac{M v_{‖}^{2}}{B_{0}}) + \frac{q ι B_{0} v_{‖}}{G} - \frac{M G^{'} v_{‖}^{2}}{G})

Finally, we solve for $v_{‖} .$

Solving, for ${\dot{v}}_{‖}$ , we have

{\dot{v}}_{∥} = - \frac{G^{'}}{G} v_{∥} \dot{ψ} + \frac{v_{∥}}{B_{0}} (\frac{\partial B_{0}}{\partial ψ} \dot{ψ} + \frac{\partial B_{0}}{\partial θ} \dot{θ}) + \frac{B_{0}}{M G} (q ι \dot{ψ} - μ \frac{\partial B_{0}}{\partial ζ}) .

We know that

\dot{ψ} = - \frac{1}{q} \frac{\partial B_{0}}{\partial θ} (\frac{M v_{‖}^{2}}{B_{0}} + μ)

\dot{θ} = \frac{1}{q} (\frac{\partial B_{0}}{\partial ψ} (μ + \frac{M v_{‖}^{2}}{B_{0}}) + \frac{q ι B_{0} v_{‖}}{G} - \frac{M G^{'} v_{‖}^{2}}{G}) .

These two have a common term of

\frac{\partial B_{0}}{\partial θ} (\frac{M v_{‖}^{2}}{B_{0}} + μ)

that gets cancelled when we sum them in our expression:

\frac{\partial B_{0}}{\partial ψ} \dot{ψ} + \frac{\partial B_{0}}{\partial θ} \dot{θ} = \frac{\partial B_{0}}{\partial θ} \frac{1}{G} (ι B_{0} v_{‖} - \frac{M G^{'} v_{‖}^{2}}{q}) .

From here, it is relatively easy to substitute $\dot{ψ}$ into the rest of our equation:

{\dot{v}}_{∥} = - \frac{G^{'}}{G} v_{∥} \dot{ψ} + \frac{v_{∥}}{B_{0}} (\frac{\partial B_{0}}{\partial ψ} \dot{ψ} + \frac{\partial B_{0}}{\partial θ} \dot{θ}) + \frac{B_{0}}{M G} (q ι \dot{ψ} - μ \frac{\partial B_{0}}{\partial ζ})

= \underset{T_{1}}{\underset{⏟}{- \frac{G^{'}}{G} v_{∥} (- \frac{1}{q} \frac{\partial B_{0}}{\partial θ} (\frac{M v_{∥}^{2}}{B_{0}} + μ))}} + \underset{T_{2}}{\underset{⏟}{\frac{v_{∥}}{B_{0}} \frac{1}{G} \frac{\partial B_{0}}{\partial θ} (v_{∥} B_{0} ι - \frac{M v_{∥}^{2}}{q} G^{'})}} + \underset{T_{3}}{\underset{⏟}{\frac{B_{0}}{M G} (- ι \frac{\partial B_{0}}{\partial θ} (\frac{M v_{∥}^{2}}{B_{0}} + μ) - μ \frac{\partial B_{0}}{\partial ζ})}}

= \underset{T_{1}^{'}}{\underset{⏟}{\frac{v_{∥}}{q G} \frac{\partial B_{0}}{\partial θ} (\frac{M v_{∥}^{2}}{B_{0}} + μ) G^{'}}} + \underset{T_{2}^{'}}{\underset{⏟}{\frac{1}{G} \frac{\partial B_{0}}{\partial θ} (v_{∥}^{2} ι - \frac{M v_{∥}^{3}}{q B_{0}} G^{'})}} + \underset{T_{3}^{'}}{\underset{⏟}{- \frac{v_{∥}^{2}}{G} ι \frac{\partial B_{0}}{\partial θ} - \frac{μ B_{0}}{M G} (ι \frac{\partial B_{0}}{\partial θ} + \frac{\partial B_{0}}{\partial ζ})}}

notice T_{1}^{'} : \frac{M v_{∥}^{3}}{q G B_{0}} \frac{\partial B_{0}}{\partial θ} G^{'} cancels with - \frac{M v_{∥}^{3}}{q G B_{0}} \frac{\partial B_{0}}{\partial θ} G^{'} from T_{2}^{'}, \frac{1}{G} \frac{\partial B_{0}}{\partial θ} v_{∥}^{2} ι cancels with - \frac{v_{∥}^{2}}{G} ι \frac{\partial B_{0}}{\partial θ} from T_{3}^{'} .

⟹ {\dot{v}}_{∥} = \frac{μ v_{∥} G^{'}}{q G} \frac{\partial B_{0}}{\partial θ} - \frac{μ B_{0}}{M G} (ι \frac{\partial B_{0}}{\partial θ} + \frac{\partial B_{0}}{\partial ζ}) .

Note that, if $B_{0}$ is invariant with respect to $ζ$ , then we have

\frac{\partial L}{\partial ζ} = - \dot{ζ} G (ψ) \frac{M v_{‖}}{B_{0}^{2}} \frac{\partial B_{0}}{\partial ζ} - μ \frac{\partial B_{0}}{\partial ζ} = 0 = \frac{d}{d t} p_{ζ} .

Note also that if $B_{0}$ depends only on $ψ$ and $χ = θ - N ζ$ then we have a new Lagrangian:

L^{'} (ψ, χ, ζ, v_{‖}) = L (ψ, χ - N ζ, ζ, v_{‖}) .

(of course, the time derivatives of each variable are also in the Lagrangian).

We can now treat $θ$ as a function of $χ, ζ$ . In this new field, we have $L^{'}$ independent of $ζ$ and will have

p_{ζ} = \frac{\partial L^{'}}{\partial \dot{ζ}} .

We know that

\dot{θ} = \dot{χ} + N \dot{ζ} .

Substituting into our old Lagrangian:

L = q (\dot{θ} ψ + \dot{ζ} (G (ψ) \frac{M v_{‖}}{q B_{0}} - ψ_{P})) - \frac{M v_{| |}^{2}}{2} - μ B_{0}

yields

p_{ζ} = q N ψ - q ψ_{P} + G \frac{M v_{‖}}{B_{0}} .

as a conserved quantity.