Charged Particles in 3D
The motion of charged particles in an electromagnetic field $(E,B)$ is governed by the Lorentz force,
\[\ddot{x} (t) = \frac{e}{m} \big[ E(x(t)) + \dot{x} (t) \times B(x(t)) \big] ,\]
where $m$ and $e$ denote the particle's mass and charge, respectively.
Canonical Formulation
The canonical form of the equations can be obtained from the Hamiltonian
\[H (x,p) = \frac{1}{2m} (p - e A(x))^2 + e \phi(x),\]
as
\[\begin{aligned} \dot{x} (t) &= \frac{\partial H}{\partial p} (x(t),p(t)) = \frac{1}{m} (p(t) - e A(x(t))), \\ \dot{p} (t) &= - \frac{\partial H}{\partial x} (x(t),p(t)) = \frac{e}{m} \nabla A(x(t)) \cdot (p(t) - e A(x(t))) - e \nabla \phi(x(t)) , \end{aligned}\]
where the fields $(E,B)$ are related to the potentials $(\phi, A)$ by
\[\begin{aligned} E (x) &= - \nabla \phi (x) , & B (x) &= \nabla \times A (x) . \end{aligned}\]
Noncanonical Formulation
The noncanonical form of the equations can be obtained from the phasespace Lagrangian
\[\begin{aligned} L (x,\dot{x},v,\dot{v}) &= (e A(x) + mv) \cdot \dot{x} - H(x,v) , & H(x,v) &= \frac{m}{2} v^2 + e \phi(x), \end{aligned}\]
as
\[\begin{aligned} \dot{x} (t) &= v (t) , & \dot{v} (t) &= \frac{e}{m} \big[ \nabla A (x(t)) \cdot \dot{x}(t) - \dot{A} (x(t)) - \nabla \phi(x(t)) \big] . \end{aligned}\]
Computing the time derivative of $A$ and using the relation between the potentials $(\phi, A)$ and the fields $(E,B)$, this can be rewritten as
\[\begin{aligned} \dot{x} (t) &= v (t) , & \dot{v} (t) &= \frac{e}{m} \big[ E(x(t)) + v (t) \times B(x(t)) \big] . \end{aligned}\]
This constitutes a noncanonical Hamiltonian system of the form
\[\dot{z} (t) = \Omega^{-T} (z(t)) \nabla H(z(t)) ,\]
with $z = (x,v)$ and the symplectic matrix $\Omega$ given by
\[\Omega = \frac{1}{m} \begin{pmatrix} \mathbb{0} & \mathbb{1} \\ - \mathbb{1} & e \hat{B} \\ \end{pmatrix} ,\]
and
\[\hat{B} = \begin{pmatrix} 0 & -B_3 & B_2 \\ B_3 & 0 & - B_1 \\ - B_2 & B_1 & 0 \\ \end{pmatrix} .\]
Modules
ChargedParticleDynamics.ChargedParticle3d — Module
The charged particle in a static electromagnetic field: the full Lorentz-force dynamics in three dimensions, with no averaging over the gyration.
Available in a canonical formulation on the phasespace $(x, p)$ and a noncanonical one on $(x, v)$, one module per equilibrium. See the Charged Particles in 3D page.
Each equilibrium is its own module, and each includes one of the two formulations — charged_particle_3d_canonical.jl on the phasespace $(x, p)$ or charged_particle_3d_noncanonical.jl on $(x, v)$. Each formulation's functions are therefore documented once, below, under TokamakSmallCylindrical and TokamakSmallNoncanonical respectively, and hold verbatim for every module sharing it. What differs between the modules is the chart, the equilibrium parameters and the initial conditions, which is what these docstrings record:
ChargedParticleDynamics.ChargedParticle3d.SingularField — Module
Charged Particle in a singular magnetic field of the form $B(x,y,z) = (x^2 + y^2)^{-3/2} e_z$.
ChargedParticleDynamics.ChargedParticle3d.SingularFieldCanonical — Module
Charged Particle in a singular magnetic field of the form $B(x,y,z) = (x^2 + y^2)^{-3/2} e_z$.
ChargedParticleDynamics.ChargedParticle3d.SymmetricField — Module
Charged Particle in an axisymmetric magnetic field of the form $B(x,y,z) = (1 + x^2 + y^2) e_z$.
ChargedParticleDynamics.ChargedParticle3d.ThetaPinchCanonical — Module
Charged Particle in an uniform magnetic field of the form $B(x,y,z) = B_0 e_z$.
ChargedParticleDynamics.ChargedParticle3d.ThetaPinchNoncanonical — Module
Charged particle in a uniform magnetic field of the form $B(x,y,z) = B_{0} e_{z}$, in the noncanonical formulation on the phasespace $z = (x, v)$.
The same equilibrium as ThetaPinchCanonical, in cartesian coordinates, so the metric is trivial and the Boris splitting of sodeproblem is a valid splitting of the model here.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallCartesian — Module
Small axisymmetric tokamak equilibrium, in cartesian coordinates $(x, y, z)$.
Major radius $R_{0} = 1$, magnetic field on axis $B_{0} = 1$ and safety factor $q_{0} = 2$.
The hard-coded μ and u of the initial_conditions_* below are inherited from earlier work on this package and their provenance is not recorded; they are known-good starting points for the charged particle model rather than values derived here.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallToroidal — Module
Small axisymmetric tokamak equilibrium, in toroidal coordinates $(r, \theta, \varphi)$.
Major radius $R_{0} = 1$, magnetic field on axis $B_{0} = 1$ and safety factor $q_{0} = 2$.
The hard-coded μ and u of the initial_conditions_* below are inherited from earlier work on this package and their provenance is not recorded; they are known-good starting points for the charged particle model rather than values derived here.
ChargedParticleDynamics.ChargedParticle3d.TokamakIterCylindrical — Module
ITER-size axisymmetric tokamak equilibrium, in cylindrical coordinates $(R, Z, \varphi)$.
Major radius $R_{0} = 6.2$, magnetic field on axis $B_{0} = 5.3$ and safety factor $q_{0} = \sqrt{2}$ — the ITER parameters of ElectromagneticFields.
The hard-coded μ and u of the initial_conditions_* below are inherited from earlier work on this package and their provenance is not recorded; they are known-good starting points for the charged particle model rather than values derived here.
ChargedParticleDynamics.ChargedParticle3d.SolovevIter — Module
Analytic ITER-like Solov'ev equilibrium, in cylindrical coordinates $(R, Z, \varphi)$.
The Solov'ev solution of the Grad-Shafranov equation fitted to ITER's shape, without an X-point. See SolovevIterXpoint for the diverted variant.
The hard-coded μ and u of the initial_conditions_* below are inherited from earlier work on this package and their provenance is not recorded; they are known-good starting points for the charged particle model rather than values derived here.
ChargedParticleDynamics.ChargedParticle3d.SolovevIterXpoint — Module
Analytic ITER-like Solov'ev equilibrium with an X-point, in cylindrical coordinates $(R, Z, \varphi)$.
The Solov'ev solution of the Grad-Shafranov equation fitted to ITER's shape, with the separatrix X-point included. See SolovevIter for the variant without one.
The hard-coded μ and u of the initial_conditions_* below are inherited from earlier work on this package and their provenance is not recorded; they are known-good starting points for the charged particle model rather than values derived here.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallCylindrical — Module
Small axisymmetric tokamak equilibrium, in cylindrical coordinates $(R, Z, \varphi)$.
Major radius $R_{0} = 1$, magnetic field on axis $B_{0} = 1$ and safety factor $q_{0} = 2$.
The hard-coded μ and u of the initial_conditions_* below are inherited from earlier work on this package and their provenance is not recorded; they are known-good starting points for the charged particle model rather than values derived here.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallCylindrical.compute_energy_error — Method
compute_energy_error(sol)ChargedParticleDynamics.ChargedParticle3d.TokamakSmallCylindrical.default_parameters — Method
The charged particle models are parameter-free — the electromagnetic field is injected as code rather than passed as parameters. The method exists so that every problem in this package can be constructed the same way.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallCylindrical.ω — Method
The symplectic two-form of the canonical formulation, $\Omega_{ij} = \partial \vartheta_{i} / \partial q^{j} - \partial \vartheta_{j} / \partial q^{i}$ with $\vartheta_{i} = g_{ii} v^{i} + A_{i}$, in the sign convention used throughout this package.
This is the $D \times D$ form on the configuration space that the LODE interface asks for, with $D = 3$ the length of q — not the $6 \times 6$ canonical form on $(q,p)$. The two are different objects, and GeometricIntegrators sizes the buffer it hands in from the problem dimension, so the canonical one does not fit.
It depends on v through the metric terms and reduces to the magnetic field, $\Omega_{ij} = \partial_{j} A_{i} - \partial_{i} A_{j}$, in cartesian coordinates.
Note that no GeometricIntegrators integrator currently evaluates the two-form of an LODE; it is carried by the problem for completeness and for use by projection methods.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallNoncanonical — Module
Small axisymmetric tokamak equilibrium, in toroidal coordinates $(r, \theta, \varphi)$, in the noncanonical formulation on the phasespace $z = (x, v)$.
Major radius $R_{0} = 1$, magnetic field on axis $B_{0} = 1$ and safety factor $q_{0} = 2$ — the same equilibrium as TokamakSmallToroidal.
The only curvilinear module of the four noncanonical ones, which makes it the one where the metric terms of the vector field matter. sodeproblem is therefore unavailable: the frozen-position kick is quadratic in $v$ in a curvilinear chart, so the Boris push is not its exact flow. Note also that the vector field carries a $1/g_{ii}$ and $g_{22} = r^{2}$, so it is singular on the magnetic axis.
The hard-coded initial velocity is inherited from earlier work on this package and its provenance is not recorded; it is a known-good starting point rather than a value derived here.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallNoncanonical.charged_particle_3d_sode_fv — Method
The velocity half of the splitting, at frozen position: $\dot{x} = 0$, $\dot{v} = E(x) + v \times B(x)$.
Writing $\hat{B}_{ij} = - \varepsilon_{ijk} B_{k}$, so that $\hat{B} v = v \times B$, the implicit midpoint rule applied to the frozen-position equation gives the update
\[\bigg( \mathbb{1} - \frac{h}{2} \hat{B} \bigg) v_{1} = \bigg( \mathbb{1} + \frac{h}{2} \hat{B} \bigg) v_{0} + h \, E(x) ,\]
which is the Boris push. For $E = 0$ it reduces to the Cayley transform of $\hat{B}$, an exact rotation, and then preserves $\vert v \vert$ exactly; with $E \neq 0$ the electric field does work on the particle and $\vert v \vert$ is not conserved.
Together with charged_particle_3d_sode_fx, which advances $\dot{x} = v$ at frozen velocity, this covers the whole of charged_particle_3d_v. The electric field used to be missing from both halves, so the splitting integrated the $E = 0$ system regardless of the equilibrium.
This is the cartesian Lorentz force. charged_particle_3d_v is no longer — it is now derived from $\Omega \dot{z} = - \nabla H$ and carries the metric — so the two agree only where the metric is trivial, which is exactly where sodeproblem will build a problem at all.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallNoncanonical.charged_particle_3d_sode_fx — Method
The drift half of the splitting, at frozen velocity: $\dot{x} = v$, $\dot{v} = 0$, whose flow is exact.
GeometricEquations calls the solution map of an SODE substep as q(q₁, t₁, q₀, t₀, params); the trailing params used to be missing from both maps here, so check_methods rejected the problem.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallNoncanonical.compute_energy_error — Method
compute_energy_error(sol)ChargedParticleDynamics.ChargedParticle3d.TokamakSmallNoncanonical.default_parameters — Method
The charged particle models are parameter-free — the electromagnetic field is injected as code rather than passed as parameters. The method exists so that every problem in this package can be constructed the same way.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallNoncanonical.dϑ — Method
The Jacobian of the one-form on the phasespace $z = (x, v)$, $(\mathrm{d}\vartheta)_{ij} = \partial \vartheta_{i} / \partial z^{j}$.
The velocity block is $\partial \vartheta_{i} / \partial v^{j} = g_{ij}$; it was previously set to zero, which is wrong in every coordinate system, cartesian included.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallNoncanonical.has_trivial_metric — Method
has_trivial_metric(t, q)Whether the coordinate system is metrically flat at q, i.e. $g_{ii} = 1$ and $\partial_{j} g_{ii} = 0$. True for the cartesian equilibria, false for the toroidal one.
This is the condition under which the Boris push of charged_particle_3d_sode_fv is the exact flow of its own substep, and hence the condition under which sodeproblem is a valid splitting of the model.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallNoncanonical.sodeproblem — Function
sodeproblem(qᵢ; kwargs...)The Boris splitting of the model into a drift at frozen velocity and a kick at frozen position.
Both the vector fields and their exact solution maps are supplied, index by index, which is the form SODE is happiest with: a composition method then uses the exact flows, while anything that wants to sub-integrate a substep has the vector field available.
The kick is solved by the Boris push, which is exact because $\dot{v} = E + v \times B$ is linear in $v$. In a curvilinear chart the kick picks up the metric-derivative terms of v₄–v₆ and becomes quadratic in $v$, so no Cayley transform solves it and the splitting has no exact flow. Rather than silently integrating the cartesian system in a chart where that is a different model, this constructor throws for any equilibrium whose metric is not trivial — in practice TokamakSmallNoncanonical, the one toroidal module of the four. Use odeproblem or iodeproblem there.
Passing nothing for the vector fields, as this used to, does not work: the typed constructor requires v::Tuple, and the fallback SODEProblem(v, args...) = SODEProblem(v, nothing, args...) then recurses on itself forever, so the constructor overflowed the stack rather than building.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallNoncanonical.v₁ — Method
The vector field, obtained from the noncanonical Hamiltonian form $\Omega \, \dot{z} = - \nabla H$ with the $\Omega$ of ω above. Writing that in blocks, with $G = \mathrm{diag}(g_{11}, g_{22}, g_{33})$ and $\Omega^{xx}_{ij} = \partial_{j} \vartheta_{i} - \partial_{i} \vartheta_{j}$,
\[\begin{pmatrix} \Omega^{xx} & G \\ -G & 0 \end{pmatrix} \begin{pmatrix} \dot{x} \\ \dot{v} \end{pmatrix} = - \begin{pmatrix} \nabla_{x} H \\ \nabla_{v} H \end{pmatrix} ,\]
whose lower block gives $\dot{x}^{i} = v^{i}$ and whose upper block gives
\[\dot{v}^{i} = \frac{1}{g_{ii}} \bigg( - \frac{\partial H}{\partial x^{i}} - \sum_{j} \Omega^{xx}_{ij} \, v^{j} \bigg) .\]
The contraction with $\Omega^{xx}$ is a cross product with the generalised magnetic field $\beta = \nabla \times \vartheta$ of β₁–β₃, so the result keeps the shape of the Lorentz force with $B$ replaced by $\beta$, $E$ by $-\nabla_{x} H$, and an overall factor $1/g_{ii}$.
This was previously the plain cartesian Lorentz force $\dot{v} = E + v \times B$, which is what this expression reduces to when $g_{ii} = 1$ and $\partial_{j} g_{ii} = 0$ — so the three cartesian equilibria are unaffected, while TokamakSmallNoncanonical, which is toroidal, was integrating a system inconsistent with its own one-form, Hamiltonian and two-form.
The $1 / g_{ii}$ is a new singularity that the cartesian expression did not have: in the toroidal chart $g_{22} = r^{2}$ vanishes on the magnetic axis, so the vector field blows up there. That is a property of the coordinate system rather than of the dynamics — the same orbit is perfectly regular in cartesian coordinates — but an orbit that approaches $r = 0$ will lose accuracy and eventually fail in the toroidal module.
ChargedParticleDynamics.ChargedParticle3d.TokamakSmallNoncanonical.ω — Method
The symplectic two-form of the noncanonical formulation, $\Omega_{ij} = \partial \vartheta_{i} / \partial z^{j} - \partial \vartheta_{j} / \partial z^{i}$ on the phasespace $z = (x, v)$, which is the convention shared with the 4D guiding centre and with GeometricProblems.
The velocity block is $\partial \vartheta_{i} / \partial v^{j} = g_{ij}$, the metric — not the identity, which is what it reduces to in cartesian coordinates only.