Guiding-centre model

This page derives the four-dimensional canonical system that the symplectic integrators advance, its normalisation in the code, and the invariants and the trapping parameter reported in the outputs. Usage: Running SIMPLE.

Phase-space Lagrangian

A particle of mass \(m = n_d m_p\) and charge \(e = n_e e_0\) (n_d, n_e) moves in a static field \(\mathbf B = \nabla\times\mathbf A\) without electric field. In Gaussian units and arbitrary spatial coordinates \(x^k\), the guiding-centre phase-space Lagrangian to first order in \(\rho_L/L\) is (Littlejohn 1983)

(1)\[L = \Bigl(m v_\parallel h_k + \frac{e}{c} A_k\Bigr)\dot x^k + \frac{mc}{e}\,\mu\,\dot\gamma - H, \qquad H = \frac{m v_\parallel^2}{2} + \mu B ,\]

with covariant components \(h_k = B_k/B\), gyrophase \(\gamma\) and magnetic moment \(\mu = m v_\perp^2/(2B)\). \(\gamma\) is cyclic, so \(\mu\) is constant; \(H\) has no explicit time dependence, so the energy is constant. Position and \(v_\parallel\) remain: four phase-space variables.

Canonical form

If the coordinates \((r,\vartheta,\varphi)\) satisfy \(A_r = h_r = 0\) (construction in Coordinates and field representation), the \(\dot r\) term of (1) vanishes and

(2)\[L = p_\vartheta\dot\vartheta + p_\varphi\dot\varphi - H, \qquad p_\vartheta = m v_\parallel h_\vartheta + \frac{e}{c}A_\vartheta, \qquad p_\varphi = m v_\parallel h_\varphi + \frac{e}{c}A_\varphi ,\]

so \((\vartheta,p_\vartheta)\) and \((\varphi,p_\varphi)\) are canonical pairs. The integrators carry

\[z = (r,\vartheta,\varphi,p_\varphi)\]

and evaluate the field at \((r,\vartheta,\varphi)\). Everything else is a function of \(z\):

(3)\[v_\parallel(z) = \frac{p_\varphi - (e/c)A_\varphi}{m h_\varphi}, \qquad p_\vartheta(z) = m v_\parallel h_\vartheta + \frac{e}{c}A_\vartheta, \qquad H(z) = \frac{m v_\parallel^2}{2} + \mu B .\]

get_val in src/field_can.f90 implements (3); get_derivatives and get_derivatives2 add first and second \(z\)-derivatives by the chain rule. Since \(r\) replaces \(p_\vartheta\), Hamilton’s equations follow by the implicit-function theorem at fixed \(p_\vartheta(z)\):

(4)\[\dot\vartheta = \frac{\partial_r H}{\partial_r p_\vartheta}, \quad \dot p_\vartheta = -\partial_\vartheta H + \frac{\partial_r H}{\partial_r p_\vartheta}\,\partial_\vartheta p_\vartheta, \quad \dot\varphi = \partial_{p_\varphi} H - \frac{\partial_r H}{\partial_r p_\vartheta}\,\partial_{p_\varphi} p_\vartheta, \quad \dot p_\varphi = -\partial_\varphi H + \frac{\partial_r H}{\partial_r p_\vartheta}\,\partial_\varphi p_\vartheta .\]

With \(\partial_{p_\varphi} H = v_\parallel/h_\varphi\) and \(\partial_{p_\varphi} p_\vartheta = h_\vartheta/h_\varphi\) (normalised units) the third equation becomes \(\dot\varphi = (v_\parallel - h_\vartheta\,\partial_r H/\partial_r p_\vartheta)/h_\varphi\), as coded in sympl_euler1_advance_angles and gauss_canfield_residual. Albert, Kasilov and Kernbichler (2020, JCP) introduced this non-canonical quadrature. It requires \(\partial_r p_\vartheta \approx (e/c)\,\psi_\text{tor}'(r) \neq 0\), which fails at the axis for \(r=\rho\); for \(r=s\) the polar singularity at \(s=0\) remains (treatment in Integrators).

Normalisation

Inside the code all lengths are in cm, fields in G, and the reference speed is the birth speed

\[v_0 = \sqrt{2E_0/m}, \qquad E_0 = 3.5\,\text{MeV}/\texttt{facE\_al}, \qquad \rho_0 = \frac{m c v_0}{e}\ \ (\texttt{ro0},\ \text{G cm}),\]

set in src/params.f90 (params_init). Time is the path length \(\tau = v_0 t\). The public state of a marker is \(z_\text{out} = (r,\vartheta,\varphi,v/v_0,\lambda)\).

The symplectic integrators use a second set of normalised quantities, marked by a bar, in which the Hamiltonian has the form \(\bar H = \bar v_\parallel^2/2 + J_\perp B\):

Quantity

Definition

In the code

parallel velocity

\(\bar v_\parallel = \sqrt2\,v_\parallel/v_0 = \sqrt2\,(v/v_0)\lambda\)

f%vpar

perpendicular invariant

\(J_\perp = 2\mu/(m v_0^2) = (v/v_0)^2(1-\lambda^2)/B\)

f%mu (G-1)

Hamiltonian

\(\bar H = 2H/(m v_0^2) = (v/v_0)^2\)

f%H

gyroradius parameter

\(\bar\rho_0 = \rho_0/\sqrt2\)

f%ro0

canonical momenta

\(\bar p_k = p_k/(m v_0/\sqrt2) = \bar v_\parallel h_k + A_k/\bar\rho_0\)

si%z(4) \(=\bar p_\varphi\), f%pth

time

\(\bar t = \tau/\sqrt2\)

si%dt \(=\Delta\tau/\sqrt2\)

With these definitions Hamilton’s equations keep their form, \(d\bar q/d\bar t = \partial\bar H/\partial\bar p\), because \(\partial H/\partial p = (v_0/\sqrt2)\,\partial\bar H/\partial\bar p\). src/simple.f90 (encode_symplectic_state) is the single conversion from \(z_\text{out}\) to \((r,\vartheta,\varphi,\bar p_\varphi)\); the reverse map is \(v/v_0 = \sqrt{\bar H} = \sqrt{J_\perp B + \bar v_\parallel^2/2}\) and \(\lambda = \bar v_\parallel/(\sqrt2\,v/v_0)\). After a collision update_momentum re-encodes the state in the same way (Collisions).

Invariants and trapping parameter

For each marker the outputs report, at the start point (compute_pitch_angle_params in src/simple_main.f90),

\[J_\perp = \frac{(v/v_0)^2(1-\lambda^2)}{B}, \qquad \kappa = \Bigl((1-\lambda^2)\frac{B_\text{max}}{B} - 1\Bigr) \frac{B_\text{min}}{B_\text{max}-B_\text{min}},\]

where \(B_\text{min}\) and \(B_\text{max}\) are the extrema of \(B\) on the start surface (get_bminmax). \(\kappa < 0\) is equivalent to \(\lambda^2 > 1 - B/B_\text{max}\), i.e. a passing orbit; \(\kappa = 0\) is the trapped-passing boundary and \(\kappa = 1\) a deeply trapped particle at \(B_\text{min}\); contr_pp skips passing markers by \(\kappa\). Without collisions \(\bar H\), \(J_\perp\) and, in axisymmetry, \(\bar p_\varphi\) are invariants; src/orbit_invariants.f90 converts between states and invariants.

References