26 Guiding Center
Guiding-center theory provides the reduced dynamical equations for the motion of charged particles in slowly varying electromagnetic fields, when the fields have weak variations over a gyroradius in space and a gyroperiod in time.
The basic equation governing the motion of a charged particle with mass m and charge e is given by \[ \frac{\mathrm{d}\mathbf{u}}{\mathrm{d}t} = \frac{e}{mc}\left( c\mathbf{E} + \mathbf{v}\times\mathbf{B} \right) \] where \(\mathbf{u}=\gamma\mathbf{v}\) represents the particle’s four-velocity, \(\mathbf{v}\) is the particle velocity, c is the speed of light while \(\mathbf{E}\) and \(\mathbf{B}\) are the electric and magnetic field vectors, respectively.
In GCA, the orbital motion of the particle is detached from its instantaneous gyration center. This relaxes any constraints related to particle gyration, allowing systematic larger time steps to be taken and thus a considerable saving in computational time. The guiding center (GC) equations are obtained in the limit of slow-varying fields. This restricts its validity to situations in which the Larmor radius remains negligible with respect to the overall electromagnetic field scale and in which the cyclotron frequency is large enough so that adiabatic invariance holds. In this respect, the GCA is most suitable for particles with large charge-to-mass ratios such as electrons.
26.1 Classical Guiding Center Theory
Here are the essentials points in the derivation of the classical guiding center equations of motion (EOM). For a complete discussion refer to Vandervoort (1960).
26.1.1 Lorentz equation in tensor form
In order to derive the Guiding Center Approximation (GCA), we would use the tensor notations in relativity. The equation of motion for a charged particle in an EM field are written in a tensorial form as \[ \frac{\mathrm{d}x_\mu}{\mathrm{d}\tau^2} = F_{\mu\nu} \frac{\mathrm{d}x_\nu}{\mathrm{d}\tau} \tag{26.1}\] where \(x_\mu\) is the particle four-coordinate, \(F_{\mu\nu}\) is the EM tensor and \(\tau\) is the particle proper time.
The EM tensor is \[ F_{\mu\nu} = \frac{e}{mc} \begin{pmatrix} 0 & B_z & -B_y & E_x \\ -B_z & 0 & B_x & E_y \\ B_y & -B_x & 0 & E_z \\ E_x & E_y & E_z & 0 \end{pmatrix} \tag{26.2}\]
26.1.2 Fundamental assumptions
The GC formalism holds under two fundamental assumptions, namely:
- The gyration radius must remain small compared to the scale length upon which the electromagnetic field changes significantly.
- The particle undergoes many gyrations before the electromagnetic field changes appreciably (slowl-varying fields pproximation). If we let \(\rho\) be the particle gyroradius (apart from a factor √2), \(x_\mu\) the particle position, \(X_\mu\) the GC position, \(\tau\) the proper time and \(F_{\mu\nu}\) the EM field tensor, the two conditions stated above can be mathematically expressed as \[ \rho \left| \frac{\partial F_{\mu\nu}}{\partial x_\alpha} \right| \ll \left| F_{\mu\nu} \right| \tag{26.3}\] and \[ \frac{1}{\omega} \left| \frac{\partial F_{\mu\nu}}{\partial x_\alpha} \right| \left| \frac{\mathbf{d}X_\beta}{\mathrm{d}\tau} \right| \ll \left| F_{\mu\nu} \right| \tag{26.4}\] where \(\omega\) is the gyrofrequency. When the previous conditions hold, we can separate the particle trajectory into a gyration and a motion of the guiding center.
26.1.3 Derivations of GC equation
The general solution of Equation 26.1 is a linear combination of the two fundamental solutions related to the four EM tensor eigenvalues \(q=\pm i\omega\) and \(q=\pm\lambda\), where \[ \begin{aligned} \omega &= \frac{e}{mc}\sqrt{\frac{1}{2}(B^2 - E^2) + \frac{1}{2}\sqrt{(B^2 - E^2)^2 + 4(\mathbf{E}\cdot\mathbf{B})^2}} \\ \lambda &= \frac{e}{mc}\sqrt{-\frac{1}{2}(B^2 - E^2) + \frac{1}{2}\sqrt{(B^2 - E^2)^2 + 4(\mathbf{E}\cdot\mathbf{B})^2}} \end{aligned} \tag{26.5}\]
\(\omega\) is a generalized relativistic form of the Larmor frequency \(\omega_\xi = eB/(mc)\) in the case when \(E\neq 0\) (which can be verified by solving in the limit \(E\to 0\)). The general solution \(x_\mu\) of Equation 26.1 is \[ x_\mu = \xi_\mu\rho \cos(\omega\tau) - \eta_\mu \rho \sin(\omega\tau) + \alpha_\mu \nu \cosh(\lambda\tau) + \beta_\mu \nu \sinh(\lambda\tau) \tag{26.6}\]
Note that the first two terms are related to a periodic motion, that is the gyration around magnetic field lines with gyroradius \(\rho\). Here \(\rho, \nu\) are constants defined by the initial conditions and \(\xi_\mu,\eta_\mu,\alpha_\mu\) and \(\beta_\mu\) are four-versors normalized in the manner \[ \xi_\mu^2 = \eta_\mu^2 = 1,\quad \text{and}\quad \alpha_\mu^2 = -\beta_\mu^2 = 1 \tag{26.7}\]
Solutions for \(\rho = 0\) and \(\nu = 0\) must be valid separately, so by replacing \(x_\mu\) for these two particular cases inside Equation 26.1 we can see that \(\xi_\mu,\eta_\mu,\alpha_\mu\) and \(\beta_\mu\) form an orthogonal set of the Minkowsky space. The orthogonality relations between this set of four-vectors in particular states that the periodic motion (gyration) in the \((\xi_mu,\eta_\mu)\)-plane is perpendicular to the acceleration motion in the \((\alpha_\mu, \beta_\mu)\)-plane, and will be useful later. At this point, another useful relation is retrieved from the squared velocity invariance \[ \frac{\mathrm{d}x_\mu}{\mathrm{d}\tau}\frac{\mathrm{d}x_\mu}{\mathrm{d}\tau} = -c^2,\quad\text{which leads to}\quad \omega^2\rho^2 - \lambda^2\nu^2 = -c^2 \tag{26.8}\]
Since we are only interested in cases in which gyration appears, the only singular eigenvalue case we must consider is the one with \(\lambda=0, \omega\neq 0\), corresponding to \(\mathbf{E}\cdot\mathbf{B}=0\) and \(|\mathbf{E}|\to 0\). The associated solution is \[ x_\mu = \xi_\mu\rho \cos(\omega\tau) - \eta_\mu \rho \sin(\omega\tau) + U_\mu \tau \tag{26.9}\]
Here the gyration motion with Larmor frequency \(\omega_\xi = \omega\) is clearly shown, along with a uniform motion (drift) with velocity \(U_\mu = (\gamma c, \mathbf{U})\) in the plane perpendicular to the \((\mathbf{E}, \mathbf{B})\)-plane.1 In this particular case, the squared velocity invariance Equation 26.8 becomes \[ \omega^2\rho^2 + U_\mu^2 = -c^2,\quad\text{or}\quad U_\mu^2 = -c^2 - \omega^2\rho^2 < -c^2 \tag{26.10}\]
Equation 26.10 is of crucial importance in the GCA formalism, because the drift motion \(U_\mu\) is the mean motion separated from the gyration and corresponds to the GC motion itself, meaning that what we just found is a relation between the energy and momentum of a particle guiding center. It is very similar to the well-known formula of the squared four-velocity \({U_\mu^\prime}^2=-c^2\), the only difference being an additional term \(\omega^2\rho^2\). This suggests interpreting the motion of the GC as the one of a particle located in the center of gyration, performing no gyration and possessing a squared four velocity as stated in Equation 26.10, in agreement with part of the particle kinetic energy being stored in the gyration motion through the term \(\omega^2\rho^2\). By separating the spatial and temporal parts of \(U_\mu\) and rearranging some terms inside Equation 26.10, we can get a relation for the GC Lorentz factor \[ \gamma_\text{GC} = \sqrt{1+\frac{|\mathbf{U}|^2}{c^2} + \frac{\omega^2\rho^2}{c^2}} = \sqrt{1+\frac{\gamma^2|\mathbf{v}|^2}{c^2} + \frac{\omega^2\rho^2}{c^2}} \tag{26.11}\]
Although we only considered a limit case, the same equation can be proven to be generally valid in the GCA formalism.
We now want to isolate the gyration motion in a general case. This is done by separating the motion as follows: \[ x_\mu = \xi_\mu x + \eta_\mu y + x_\mu^\prime + X_\mu \tag{26.12}\] where \(x_\mu\) is the space-time position of the particle and \(X_\mu\) is the position of the guiding center. We associate two variables \(x\) and \(y\) to describe the periodic motion in the \((\xi, \eta)\)-plane, and the term \(x_\mu^\prime\) will account for the periodic motions which do not lie in the same plane. To more easily describe the gyration, we also choose a new pair of variables \[ \zeta = \frac{1}{\sqrt{2}}\left( x+ iy \right)\quad\text{and}\quad \zeta^\ast = \frac{1}{\sqrt{2}}\left( x - iy \right) \tag{26.13}\] so that \(\xi_\mu x + \eta_\mu y = \delta_\mu\zeta + \sigma_\mu\zeta^\ast\). The new variables will allow us to choose \(\zeta_\mu\) and \(\delta_\mu\) in a more convenient way and let some factor \(e^{i\phi}\) account for the initial conditions of the problem. We rewrite the particle coordinate using Equation 26.13 \[ x_\mu = \delta_\mu\zeta + \sigma_\mu\zeta^\ast + x_\mu^\prime + X_\mu \tag{26.14}\]
The next step is substituting this expression into Equation 26.1 and expand the second derivative, which leads to a lengthy expression containing eight unkowns: \(\zeta,\zeta^\ast\), two components of \(x_\mu\) (the condition with the \((\sigma, \delta)\)-plane sets the other two), and four components of \(X_\mu\).
The conditions for the GCA expressed in Equation 26.4 implicitly contain a parameter of smallness \(\epsilon =|u/\omega L|\ll 1\), where u is the typical particle four-velocity and L the scale length upon which the change in the EM fields \(F_{\mu\nu}\) is comparable to \(F_{\mu\nu}\) (slowly-varying fields approximation). This allows us to approximate \(F_{\mu\nu}(x_\mu)\) (particle coordinate) by using a Taylor expansion about \(X_\mu\) (the GC coordinate) \[ \begin{aligned} F_{\mu\nu}(x_\mu) =& F_{\mu\nu}(X_\mu) \\ &+ \frac{\partial F_{\mu\nu}}{\partial x_\xi}(\delta_\xi \zeta + \sigma_\xi \zeta^\ast) \\ &+ \frac{1}{2} \frac{\partial^2 F_{\mu\nu}}{\partial x_\xi\partial x_\pi}\left[ \delta_\xi\delta_\pi \zeta^2 + \sigma_\xi\sigma_\pi {\zeta^\ast}^2 + (\sigma_\xi \delta_\pi + \sigma_\pi \delta_\xi)|\zeta|^2 \right] \\ &+ \frac{\partial F_{\mu\nu}}{\partial x_\xi} x_\xi^\prime + ... \end{aligned} \]
We should now substitute \(F_{\mu\nu}(x_\mu)\) inside Equation 26.1 using the new coordinates in order to get a quite long equation containing all quantities evaluated at the GC position \(X_\mu\) up to 2nd−order, and group all 2nd−order terms on the right-end side as \(G_\mu\) \[ \begin{aligned} &\delta_\mu\frac{\mathrm{d}^2\zeta}{\mathrm{d}\tau^2} + \sigma_\mu\frac{\mathrm{d}^2\zeta^\ast}{\mathrm{d}\tau^2}+\left( i\omega\delta_\mu + 2\frac{\mathrm{d}\zeta}{\mathrm{d}\tau} \right)\frac{\mathrm{d}\zeta}{\mathrm{d}\tau} \\ &+\left( -i\omega\sigma_\mu + 2\frac{\mathrm{d}\sigma_\mu}{\mathrm{d}\tau} \right) \frac{\mathrm{d}\zeta^\ast}{\mathrm{d}\tau} + \frac{\mathrm{d}^2\delta_\mu}{\mathrm{d}\tau^2}\zeta + \frac{\mathrm{d}^2\sigma_\mu}{\mathrm{d}\tau^2}\zeta^\ast \\ &- \frac{\partial F_{\mu\nu}}{\partial x_\xi}\left( \delta_\xi \zeta + \sigma_\xi \zeta^\ast \right)\frac{\mathrm{d}X_\nu}{\mathrm{d}\tau} - F_{\mu\nu}\left( \frac{\mathrm{d}\delta_\nu}{\mathrm{d}\tau}\zeta + \frac{\mathrm{d}\sigma_\nu}{\mathrm{d}\tau}\zeta^\ast \right) \\ &- \frac{\partial F_{\mu\nu}}{\partial x_\xi}\left( \delta_\xi \zeta + \sigma_\xi \zeta^\ast \right) \left( \delta_\nu\frac{\mathrm{d}\zeta}{\mathrm{d}\tau} + \sigma_\nu\frac{\mathrm{d}\zeta^\ast}{\mathrm{d}\tau} \right) \\ &+ \left( \frac{\mathrm{d}^2 X_\nu}{\mathrm{d}\tau^2} - F_{\mu\nu}\frac{\mathrm{d}X_\mu}{\mathrm{d}\tau} \right) + \left( \frac{d^2 x_\mu^\prime}{\mathrm{d}\tau^2} - F_{\mu\nu}\frac{\mathrm{d}x_\nu^\prime}{\mathrm{d}\tau} \right) = G_\mu \end{aligned} \tag{26.15}\]
An expression for \(\zeta\) is needed to further proceed, but since the derivation is quite long we will not explicitly derive it here. The basic idea is to use the orthogonality equations for the new four-vectors \(\sigma_\mu\) and \(\delta_\mu\) and the antisymmetry property of the EM tensor inside Maxwell’s equations written in tensorial form, and then assume that a solution for \(X_\mu\) has already been found. Considered that \(x_\mu\) is a 1st−order term in the GCA, by ignoring 2nd−order quantities and writing terms independent from \(\zeta\) as a numerical value \(P\), then Equation 26.15 becomes \[ \frac{\mathrm{d}^2\zeta}{\mathrm{d}\tau^2} + i\Omega\frac{\mathrm{d}\zeta}{\mathrm{d}\tau} + \frac{i}{2}\frac{\mathrm{d}\Omega}{\mathrm{d}\tau}\zeta + \Omega a \zeta + P = 0 \tag{26.16}\] where \(\Omega a \zeta\) contains all the linear terms in \(\zeta\). Consequently, the corresponding solution \(\zeta\), by analogy to the Wentzel–Kramers–Brillouin (WKB) approximation in quantum mechanics, is \[ \zeta = \rho_0 \sqrt{\frac{\omega_0}{\Omega}}\exp\left[ -i\left( \Phi + \int_{\tau_0}^\tau \sigma_1 \mathrm{d}\tau \right) \right] \tag{26.17}\] to the 0th-order, being \[ \Omega = \omega - 2i\sigma_\mu \frac{\mathrm{d}\delta_\mu}{\mathrm{d}\tau},\quad \Phi = \int_{\tau_0}^{\tau}\Omega\mathrm{d}\tau - \phi,\quad\text{and}\quad \sigma = a - \frac{a^2}{\Omega} \tag{26.18}\] The subscripts 0, 1 indicate the associated order of approximation for quantities defined above.
Since now we know the expression for \(\zeta\), we can write the 1st−order EOM for \(X_\mu\). For this purpose, we require that in the left-end side of Equation 26.15 the terms containing the 1st− and 2nd−order derivatives of \(X_\mu\) are balanced by the remaining nonoscillatory terms, that is \[ \frac{\mathrm{d}^2 X_\mu}{\mathrm{d}\tau^2} - F_{\mu\nu}\frac{\mathrm{d}X_\nu}{\mathrm{d}\tau} = \frac{\partial F_{\mu\nu}}{\partial x_\xi}\left( \delta_\xi \sigma_\nu \frac{\mathrm{d}\zeta^\ast}{\mathrm{d}\tau} + \delta_\nu \sigma_\xi \frac{\mathrm{d}\zeta}{\mathrm{d}\tau} \right) \tag{26.19}\]
Having found the solution Equation 26.17 and making use once again of the orthogonality relations for \(\delta_\mu\) and \(\sigma_\mu\) and the tensorial Maxwell’s equations, we can write the equation of motion for the GC \[ \frac{\mathrm{d}^2 X_\mu}{\mathrm{d}\tau^2} - F_{\mu\nu}\frac{\mathrm{d}X_\nu}{\mathrm{d}\tau} + \rho_0^2 \omega_0\frac{\partial\omega}{\partial x_\mu} = 0 \tag{26.20}\] where \(\omega\) (Equation 26.5) is the generalized Larmor frequency, \(\rho_0^2 \omega_0 = \mu_0 = mc\gamma^2 v_\perp^2 / 2eB\) is the 0th-order approximation (apart from a constant factor) to the particle magnetic moment \(\mu\) and \(\omega_0 = eB/(mc)\) represents the lowest-order approximation to the Larmor frequency in the case \(E\ll B\).
More precisely, the magnetic moment should be regarded as constant only in the reference frame moving at \(\mathbf{v}_E\), because in that particular case \(\mathbf{E} = 0\). Indeed, a formal analysis (Northrop and Teller (1960)) shows that the corresponding constant of motion in a general reference frame is actually an asymptotic series expansion in a smallness parameter \(\epsilon=u/(\omega_0 L)\) in the form \(\mu = (e/c)\left[ \mu_0 + \epsilon \mu_1 + \epsilon^2 \mu_2 + ... \right]\) so that \(\mu_0\) is not constant and can still vary in compliance with the adiabatic theory. Nevertheless, extensive numerical testing confirms that the errors due to assuming \(\mathrm{d}\mu_0 / \mathrm{d}t \approx 0\) are at most of the same order of those introduced by the GC formalism, as a consequence of the slowly varying field condition which prevents sensible changes in the magnetic moment. Hence, we safely assume that \(\mu \approx (e/c)\mu_0\) is invariant.
26.1.4 First order solution for the GCA equations of motion
Although deceitful simple, the GC equation of motion Equation 26.20 is more conveniently cast in a form in which the GC velocity appears explicitly, under the nearly-crossed (EM) fields condition \[ \frac{\mathbf{E}_\parallel \cdot\mathbf{B}}{B^2 - E_\perp^2} \ll 1 \tag{26.21}\] where the subscripts \(\parallel\) and \(\perp\) indicate the parallel and perpendicular components of the electric field with respect to the magnetic field unit vector \(\hat{b}\). It can be proven that Equation 26.21 allows the EM tensor field \(F_{\mu\nu}\) to be split into a contribution \(F_{\mu\ni}^{(0)}\) constructed solely from \(\mathbf{E}_\perp\) and \(\mathbf{B}\), plus a correction term \(F_{\mu\nu}^{(1)}\) depending on \(\mathbf{E}_\parallel\). A formal analysis leads to the conclusion that the GC four-velocity \(U_\mu \equiv \mathrm{d}X_\mu/\mathrm{d}\tau\) can be decomposed into a 0th-order contribution \(\mathrm{U}^{(0)}\) (which contains the drift velocity \(\mathbf{v}_E = c\mathbf{E}_\perp \times\mathbf{B}/B^2\) as well as the velocity component parallel to the magnetic field line \(v_\parallel \hat{b}\)) plus 1st-order correction \(\mathbf{U}^{(1)}\), leading to \[ U_\mu = (\gamma c, \mathbf{U}) \simeq (\gamma c, \gamma v_\parallel \hat{b} + \gamma \mathbf{v}_E + \mathbf{U}^{(1)}) \tag{26.22}\] where \(\gamma\) is the particle Lorentz factor \(\gamma = (1- v^2/c^2)^{-1/2}\). An equation for \(\mathbf{U}^{(1)}\) can be derived from Equation 26.20 using a recursive approach as shown below, while from the spatial component one obtains an equation (to the same order) for the parallel component of the GC four-velocity \(\gamma v_\parallel\). These yield the (1st-order) GCA system of ordinary differential equations (ODEs) \[ \begin{aligned} \frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t} &= \mathbf{v}_E + v_\parallel \hat{b} + \frac{\gamma_E^2}{B}\hat{b}\times\left[ \frac{mc\gamma}{e}\left( v_\parallel\mathcal{L}(\hat{b}) + \mathcal{L}(\mathbf{v}_E) \right) + \mathcal{M}\left( \frac{B}{\gamma_E} \right) + \frac{v_\parallel E_\parallel}{c}\mathbf{v}_E \right] \\ \frac{\mathrm{d}(\gamma v_\parallel)}{\mathrm{d}t} &= \frac{e}{m}E_\parallel - \gamma\hat{b}\cdot\mathcal{L}(\mathbf{v}_E) - \frac{\mu}{\gamma m}\hat{b}\cdot\nabla\left(\frac{B}{\gamma_E}\right) \end{aligned} \tag{26.23}\] where \(\mathrm{d}\mathbf{X}/\mathrm{d}t = \mathbf{U}^{(1)}/\gamma\) represents the velocity of the guiding center, \(e/m\) is the particle charge-to-mass ratio and \(\mu_0 = c\mu/e\). The operators \(\mathcal{L}\) and \(\mathcal{M}\) are defined as \[ \begin{aligned} \mathcal{L}(x) &= \frac{\partial\mathbf{x}}{\partial t} + v_\parallel(\hat{b}\cdot\nabla)\mathbf{x} + (\mathbf{v}_E\cdot\nabla)\mathbf{x} \\ \mathcal{M}(x) &= \frac{\mu}{e\gamma}\left[ \frac{\mathbf{v}_E}{c}\frac{\partial x}{\partial t} + c\nabla x \right] \end{aligned} \tag{26.24}\] where \(\gamma_E = (1 - E_\perp^2/B^2)^{-1/2}\) is the Lorentz factor associated to the drift velocity \(\mathbf{v}_E\).
In the GCA Equation 26.23 information about motion perpendicular to the magnetic field is lost and the only component of the particle velocity that is actually evolved in time is \(u_\parallel\), because the magnetic moment \(\mu_0\approx c\mu/e\) is now considered a constant of motion. Each term on the right hand side of Equation 26.23 corresponds to a specific drift motion. Indeed, the first term represents the \(\mathbf{E}\times\mathbf{B}\) perpendicular drift, while the second one accounts for particle motion in the direction parallel to magnetic field. The first two terms in square bracket, \(\mathcal{L}(\hat{b})\) and \(\mathcal{L}(\mathbf{v}_E)\), describe the field curvature and polarization drifts2, respectively. The next term is the \(\nabla-B\) drift, while the last one represents a relativistic drift in the direction given by \(\hat{b}\times\mathbf{v}_E\).
Equation 26.23 lose their validity when either \(B\to 0\) (magnetic null) or \(E_\perp \ge B\), i.e. when the nearly-crossed EM fields condition, Equation 26.21, is violated. While the first condition can easily occur inside a current sheet, the second may occur in a non-ideal MHD regime only. Violation of either condition breaks down the GCA and can lead to severe numerical errors.
26.1.5 Recursive solution
An explicit solution for Equation 26.20 can be found under Equation 26.21, which introduces a new parameter of smallness \(\lambda / \omega \ll 1\).
26.1.6 Energy conservation
The time component of Equation 26.20 leads to an evolutionary ODE for the particle energy as a function of time, namely \[ \frac{\mathrm{d}\gamma c^2}{\mathrm{d}t} = \frac{e}{m}\frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t}\cdot\mathbf{E} + \frac{\mu}{m}\frac{\partial}{\partial t}\left( \frac{B}{\gamma_E} \right) \tag{26.25}\]
In the time-independent case, Equation 26.25 clearly shows that, when the GC velocity is perpendicular to the electric field, no acceleration occurs and the particle energy is conserved. For numerical purposes, however, Equation 26.25 is not solved for retrieving the Lorentz \(\gamma\)-factor, since large values of the right hand side could easily violate the condition \(\gamma\ge 1\) at the truncation level of the scheme, unless small time steps are taken. In practice, to an 1st-order approximation, the particle velocity can be decomposed into a parallel component \(v_\parallel\), a drift component \(v_E\) and a gyration component \(v_\perp\). This leads to the normalization condition \[ \gamma^2( c^2 - v_\parallel^2 - v_\perp^2 - v_E^2) = c^2 \tag{26.26}\]
From the definition of the magnetic moment and the Larmor frequency \(\mu_0 = mc\gamma^2 v_\perp^2/(2eB) = \gamma^2 v_\perp^2/(2\omega_0)\), we compute the particle Lorentz factor as \[ \gamma = \sqrt{\frac{c^2 + \gamma^2 v_\parallel^2 + 2\mu_0 \omega_0}{c^2 - v_E^2}} \tag{26.27}\]
Equation 26.27 replaces Equation 26.25 in common numerical solvers to get the Lorentz factor.
26.1.7 Numerical implementation
The GCA Equation 26.23 contains a set of 4 ODEs with unknowns \((\mathbf{X}, u_\parallel)\). The terms which require evaluation at grids are
- \(\mathbf{B}\)
- \(c\mathbf{E}\)
- \((\hat{b}\cdot\nabla)\hat{b}\)
- \((\hat{b}\cdot\nabla)\mathbf{v}_E\)
- \((\mathbf{v}_E\cdot\nabla)\mathbf{v}_E\)
- \((\mathbf{v}_E\cdot\nabla)\hat{b}\)
- \(\nabla(B/\gamma_E)\)
Mignone, Haudemand, and Puzzoni (2023) reports the implementation of GCA in the PLUTO MHD model.
26.1.8 Nonrelativistic case
In the nonrelativistic case, \(\gamma = 1\) so we don’t need to solve Equation 26.27. Equation 26.23 can be simplified into \[ \begin{aligned} \frac{\mathrm{d}\mathbf{X}}{\mathrm{d}t} &= \\ \frac{\mathrm{d}(v_\parallel)}{\mathrm{d}t} &= \end{aligned} \tag{26.28}\]
26.2 Hamiltonian Theory of Guiding Center Motion
Following the review by Cary and Brizard (2009).
Canonical and noncanonical Hamiltonian formulations of guiding-center motion offer improvements over non-Hamiltonian formulations: Hamiltonian formulations possess Noether’s theorem (hence invariants follow from symmetries), and they preserve the Poincaré invariants (so that spurious attractors are prevented from appearing in simulations of guiding-center dynamics). Hamiltonian guiding-center theory is guaranteed to have an energy conservation law for time-independent fields—something that is not true of non-Hamiltonian guiding-center theories. The use of the phase-space Lagrangian approach facilitates this development, as there is no need to transform a priori to canonical coordinates, such as flux coordinates, which have less physical meaning.
26.2.1 Noncanonical Hamiltonian guiding center theory
The guiding-center phase space consists of
the guiding-center position \(\mathbf{X}\), essentially the center of the helix;
the guiding-center parallel velocity variable \(u\equiv \hat{b}\cdot\dot{\mathbf{X}}\);
the (lowest-order) magnetic moment, \[ \mu \equiv \frac{m |\mathbf{w}|^2}{2B(\mathbf{X}, t)} \tag{26.29}\] where \(\mathbf{w} = \mathbf{v}_\perp - \mathbf{v}_E\) is the perpendicular velocity in the local frame moving with the \(\mathbf{E}\times\mathbf{B}\) drift velocity \(\mathbf{v}_E\equiv c\mathbf{E}\times\hat{b}/B\).
and the ignorable gyrophase \(\zeta\), which gives the location of the particle on the circle about the guiding center.
As there are still six variables parametrizing phase space, there is no loss of information in making the guiding-center transformation \((\mathbf{x},\mathbf{v}) \to (\mathbf{X},u,\mu,\zeta)\). For the sake of simplicity of notation, we occasionally use the gyroaction variable \(J \equiv mc/e\mu\) instead of the magnetic moment \(\mu\) whenever we need to refer to the action-angle coordinates \(J,\zeta\) associated with gyromotion.
The equations of motion for these variables are given by the guiding-center phase-space Lagrangian, \[ \mathcal{L}_\text{gc}(\mathbf{X},u,\mu,\zeta;t) = \left[ \frac{e}{c}\mathbf{A}(\mathbf{X},t) + mu\hat{b}(\mathbf{X},t) \right]\cdot\dot{\mathbf{X}} + J\dot{\zeta} - H_\text{gc} \tag{26.30}\] in which the guiding-center Hamiltonian is given by \[ H_\text{gc}(\mathbf{X},u,\mu;t) = \frac{m}{2}u^2 + \mu B(\mathbf{X},t) + e\phi(\mathbf{X},t) - \frac{m}{2}\left| \mathbf{v}_E \right|^2 \tag{26.31}\]
The arguments are shown here to emphasize that, for example, the magnetic-field strength \(B(\mathbf{X},t)\) is evaluated at the guiding-center position \(\mathbf{X}\), not at the particle position \(\mathbf{x}\). Here and throughout, the effects of a gravitational field may by found by adding the gravitational \(m\Phi_G\) to the electrostatic potential energy \(e\Phi\). In addition, we now drop the adjective phase space, as the guiding-center Lagrangian is always henceforth a phase-space Lagrangian.
The guiding-center Lagrangian Equation 26.30 comes not simply from gyrophase averaging, but from a transformation from the physical (particle) variables \((\mathbf{x},\mathbf{v})\) to the guiding-center variables \((\mathbf{X},u,\mu,\zeta)\). The relation between the physical location \(\mathbf{x}\) and the guiding-center position \(\mathbf{X}\) is \[ \mathbf{x} \equiv \mathbf{X} + \pmb{\rho} \tag{26.32}\] where \(\pmb{\rho}\) denotes the gyroradius displacement vector in the frame drifting with \(\mathbf{v}_E\). Here we simply note that the displacement vector \(\pmb{\rho}\equiv \tilde{\pmb{\rho}} + \bar{\pmb{\rho}}\) has a part (denoted by tilde) that is explicitly gyrophase dependent and a part (denoted by bar) that is gyrophase independent. In what follows, we show that the latter part \[ \bar{\pmb{\rho}} = \frac{\hat{b}}{\Omega}\times\mathbf{v}_E = \frac{c\mathbf{E}_\perp}{\Omega B} \tag{26.33}\] denotes the guiding-center polarization displacement (Kaufman, Phys. Fluids, 1986). …
Our guiding-center Lagrangian Equation 26.30 is obtained from an ordering in which the scalar potential \(\Phi\) is one order lower than the particle kinetic energy, unlike previous derivations of the Hamiltonian theory of guiding-center motion. In this ordering, the electric drift \(\mathbf{v}_E\) is of the same order as the perpendicular velocity \(\mathbf{w}\), as in some non-Hamiltonian calculations (Northrop (1963)). As we will see, this ordering allows us to obtain the polarization drift in the same order as the curvature and \(\nabla B\) drifts, although it appears differently in the theory. However, for consistency the parallel electric field \(E_\parallel = \mathbf{E}\cdot\hat{b}\) must be smaller by one order than the perpendicular field \(\mathbf{E}_\perp\) (Equation 26.21).
The variables \((J,\zeta)\) appear in canonical form in the symplectic part of the guiding-center Lagrangian Equation 26.30 as \(J\dot{\zeta}\) while the guiding-center Hamiltonian \(H_\text{gc}\) depends on \(J\) (or \(\mu\)) alone. The Hamilton equations for \((J,\zeta)\) are \[ \begin{aligned} \dot{J} &= -\frac{\partial H_\text{gc}}{\partial \zeta} \equiv 0 \\ \dot{\zeta} &= \frac{\partial H_\text{gc}}{\partial J} \equiv \Omega \end{aligned} \tag{26.34}\]
Equation 26.34 shows that the gyroaction (or magnetic moment) is conserved by the guiding-center equations of motion. This also follows from Noether’s theorem since the gyrophase \(\zeta\) is an ignorable coordinate, i.e., only its time derivative appears in the guiding-center Lagrangian Equation 26.30.
If one is concerned with only the motion of the guiding center and not the evolution of the gyrophase, the term linear in \(\dot{\zeta}\) can be dropped from the guiding-center Lagrangian, as it does not affect the equations of motion of the other variables, \(\mathbf{X}\) and \(u\). In the evolution equations for \(\mathbf{X}\) and \(u\), the adiabatic invariant \(\mu\) (or \(J\)) does appear but only as a guiding-center dynamical parameter.
Variation of the guiding-center Lagrangian Equation 26.30 with respect to the variable \(u\) gives the Euler-Lagrange equation \[ 0 = \frac{\partial\mathcal{L}}{\partial u} = m\hat{b}\cdot\dot{\mathbf{X}} - \frac{\partial H_\text{gc}}{\partial u} \] which yields \[ u \equiv \hat{b}\cdot\dot{\mathbf{X}} \tag{26.35}\]
Thus, the guiding-center Lagrangian Equation 26.30 dictates that \(u\) is the velocity of the guiding center in the direction of the magnetic field at the guiding center.
Last, we vary the Lagrangian Equation 26.30 with respect to the guiding-center position \(\mathbf{X}\), and obtain the Euler-Lagrange equation \[ \begin{aligned} m\dot{u}\hat{b} =& e\mathbf{E} - \mu\nabla B + \frac{m}{2}\nabla |\mathbf{v}_E|^2 - mu\frac{\partial\hat{b}}{\partial t} \\ & + \dot{\mathbf{X}}\times\left( \frac{e}{c}\mathbf{B}+mu\nabla\times\hat{b} \right) \\ \equiv& e\left( \mathbf{E}^\ast + \frac{1}{c}\dot{\mathbf{X}}\times\mathbf{B}^\ast \right) \end{aligned} \tag{26.36}\] where the effective EM fields \[ \mathbf{E}^\ast \equiv -\nabla\Phi^\ast - \frac{1}{c}\frac{\partial\mathbf{A}^\ast}{\partial t}\quad\text{and}\quad \mathbf{B}^\ast \equiv\nabla\times\mathbf{A}^\ast \tag{26.37}\] are defined in terms of the effective EM potentials \[ \begin{aligned} e\Phi^\ast &\equiv e\Phi + \mu B -\frac{m}{2}|\mathbf{v}_E|^2 \\ \mathbf{A}^\ast &\equiv \mathbf{A} + \frac{mc}{e}u\hat{b} \end{aligned} \tag{26.38}\]
The guiding-center canonical momentum is now simply expressed as \(e\mathbf{A}^\ast /c\) and the guiding-center Hamiltonian is \(e\Phi^\ast + mu^2/2\).
We obtain the rate of change of the variable \(u\) by taking the scalar product of Equation 26.36 with the effective magnetic field \(\mathbf{B}^\ast\) \[ \dot{u} = -\frac{\mathbf{B}^\ast}{m B_\parallel^\ast}\cdot\left( \nabla H_\text{gc} + \frac{e}{c}\frac{\partial\mathbf{A}^\ast}{\partial t} \right) \equiv \frac{e}{m}\frac{\mathbf{B}^\ast}{B_\parallel^\ast}\cdot\mathbf{E}^\ast \tag{26.39}\] with \(B_\parallel^\ast = \hat{b}\cdot\mathbf{B}^\ast\) the effective magnetic field in the parallel direction.
The time derivative Equation 26.39 contains terms that are higher order (in gyroradius) compared with the dominant terms, which are all that is usually kept. These higher-order terms, however, are needed for energy conservation.
The guiding-center velocity \(\dot{\mathbf{X}}\) comes from the vector product of Equation 26.36 with \(\hat{b}\) which, using Equation 26.35, yields \[ \begin{aligned} \dot{\mathbf{X}} &= \frac{\mathbf{B}^\ast}{m B_\parallel^\ast}\frac{\partial H_\text{gc}}{\partial u} + \frac{c\hat{b}}{eB_\parallel^\ast}\times\left( \nabla H_\text{gc} + \frac{e}{c}\frac{\partial\mathbf{A}^\ast}{\partial t} \right) \\ &= u\frac{\mathbf{B}^\ast}{B_\parallel^\ast} + \mathbf{E}^\ast\times\frac{c\hat{b}}{B_\parallel^\ast} \end{aligned} \tag{26.40}\]
If the effective fields Equation 26.37 were replaced by the standard fields \((\mathbf{E},\mathbf{B})\), Equation 26.39 and Equation 26.40 would be the equations of motion for a particle in straight, constant EM fields.
The guiding-center equations of motion Equation 26.39 and Equation 26.40 can be used to derive the guiding-center energy equation \[ \frac{\mathrm{d}H_\text{gc}}{\mathrm{d}t} = \frac{\partial H_\text{gc}}{\partial t} + \dot{\mathbf{X}}\cdot\nabla H_\text{gc} + \dot{u}\frac{\partial H_\text{gc}}{\partial u} = e\frac{\partial\Phi^\ast}{\partial t} -\frac{e}{c}\frac{\partial\mathbf{A}^\ast}{\partial t}\cdot\dot{\mathbf{X}} \tag{26.41}\] which implies that the guiding-center energy \(E_\text{gc}\equiv mu^2/2+e\Phi^\ast\) is a constant of the motion for time-independent fields.
An important remark must be made here concerning the polarization drift, which is absent from the guiding-center velocity Equation 26.40. This drift, however, is critical for obtaining the dielectric response of a low-frequency plasma. Instead, it appears in the transformation Equation 26.32 itself, i.e., the derivative of this relation gives \[ \dot{\mathbf{x}} = \dot{\mathbf{X}} + \dot{\tilde{\pmb{\rho}}} + \mathbf{v}_\text{pol} \tag{26.42}\] with \[ \mathbf{v}_\text{pol} \equiv \frac{\mathrm{d}\bar{\pmb{\rho}}}{\mathrm{d}t} \tag{26.43}\] representing the polarization drift, and where \(\dot{\tilde{\pmb{\rho}}}=\dot{\zeta}\partial \tilde{\pmb{\rho}}/\partial \zeta + ...\) consists of terms that oscillate on the gyroperiod time scale. The polarization drift is a pure derivative and, hence, can always be integrated. This implies, in particular, that the polarization drift cannot lead to diffusion even in a turbulent field. This is important, as the difference between the guiding center and the average location (found by dropping the second, oscillating term in Equation 26.32) is the polarization, the integral of the polarization drift. This difference must remain small or else the theory, which assumes that the particle remains close to \(\mathbf{X}\), would break down.
An alternate set of guiding-center equations of motion may be derived in which the polarization drift appears explicitly in the guiding-center velocity \(\dot{\mathbf{X}}\) by choosing \(\bar{\pmb{\rho}}\equiv 0\) instead of Equation 26.33.
For verification simply substitute \(x_\mu\) and \(U_\mu\) in the left and right-end side of Equation 26.1, respectively, to find that \(F_{\mu\nu} U_\nu = 0\), and rewrite the indexes of such equation explicitly.↩︎
Notice that the \(\mathcal{L}()\) operator is the Lagrangian derivative \(\mathrm{d}/\mathrm{d}t\).↩︎