Canonical Hamiltonian guiding center theory and classical intrinsic magnetic moment

Abstract

Guiding center theories are crucial in astrophysics, space plasmas, fusion research, and arc plasmas for addressing the multi-scale dynamics of magnetized plasmas. In this paper, we derive a new Lagrangian function of guiding center in 6D variables (\mathit{X},\dot{\mathit{X}}) different from the Littlejohn’s one by employing two different approaches. Based on the new Lagrangian function, we prove that the guiding center dynamics can be generally described as a constrained canonical Hamiltonian system with two constraints in six dimensional phase space. By explicitly expressing the Lagrangian multipliers, we can reformulate the constrained Hamiltonian system into equivalent Hamiltonian−Dirac equations in coordinates (\mathit{X},\mathit{p}). In these coordinates, the guiding center dynamics’ solution flow resides on a symplectic sub-manifold, ensuring the exact conservation of the symplectic structure. Thus, we identify the canonical coordinates of the guiding center. In this context, the guiding center behaves as a pseudo-particle with an intrinsic magnetic moment, effectively replacing charged particle dynamics over time scales longer than the gyro-period. The complete dynamical behaviors, including acceleration and force, of the guiding center pseudo-particle can be derived clearly and consistently from this theory. Additionally, this framework enables the systematic development of related theories, such as symplectic guiding center algorithms, canonical gyro-kinetic theory, and canonical particle-in-cell methods, enhancing the global accuracy of gyrokinetics and associated numerical techniques. The theory also sheds light on the origin of the intrinsic magnetic moment within the scope of classical mechanics.

Graphical abstract

Keywords

guiding center dynamics / canonical Hamiltonian system

1 Introduction

Plasmas at different characteristic scales are described by different theories, such as kinetic theory, gyrokinetics, drift kinetics, the two-fluid model, and magnetohydrodynamics (MHD). The completeness and self-consistency of these hierarchy theories for distinct scale levels are critical to handle complex plasma systems. The concept of guiding center and gyrokinetic theory has been widely applied to relieve the multi-scale problem arising in magnetized plasmas and has seen significant advancements and applications across various domains of physical phenomena above the gyro-motion scale in astrophysics [1], space physics [2], fusion research [2–8], arc plasmas [9], etc. One such application is in the analysis of guiding-center drift ions, a phenomenon that only occurs in the presence of strong magnetic fields [10].

The guiding-center model has also been applied in the first nonlinear gyrokinetic simulation of solar wind turbulence, providing a comprehensive resolution of the entire range of scales from the ion (proton) to the electron Larmor radius, with the correct mass ratio [1]. This work paved the way for understanding the complex nature of turbulence in space plasmas. Additionally, the model has been used to uncover a novel mechanism for spontaneous rotation in tokamaks, arising from the partially acoustic character of ITG turbulence, thus contributing to the study of plasma confinement and transport phenomena [3]. In gyrokinetic turbulence, particularly in the case of ion temperature gradient (ITG)-driven turbulence, the model has been instrumental in revealing the similarities between energy transfer in gyrokinetic turbulence and fluid turbulence, specifically through the study of the free energy cascade [11]. Further applications of the gyrokinetic equations have revealed how self-organisation of plasma microturbulence can lead to improved confinement states in tokamaks, highlighting the importance of the interface between plasma and material boundaries [5]. In this context, the propagation of turbulence activity beyond regions of convective drive was shown to play a critical role in the establishment of these improved confinement states.

In the field of laboratory plasma physics, the guiding-center model has been successfully employed to capture the motion of photoelectrons [9], providing insight into their behavior in various plasma environments. Similarly, gyrokinetic theory has been applied to study the physics of kinetic Alfvén waves, expanding our understanding of these waves and their role in plasma dynamics [4]. The model has also been utilized to simulate auroral electron acceleration, shedding light on the complex processes that drive these phenomena [2]. The model has also been used to investigate multi-physics problems in plasma dynamics, further demonstrating its versatility [6]. Recent experimental studies have provided conclusive evidence of a self-generated current in the KSTAR tokamak, a phenomenon not fully explained by existing theories, but successfully captured using the guiding-center model [7]. In addition, the gyrokinetic approach has been crucial in explaining the breakdown of the “pure Alfvénic state” in the short-wavelength kinetic regime [8]. Gyrokinetic theory has also been shown to rest on the intrinsic nature of particle dynamics as a boundary layer problem, offering a new perspective on plasma behavior at microscopic scales [12]. Lastly, the guiding-center dynamics in varying external Maxwell fields have been analyzed using a relativistically covariant action principle, extending the well-known Vandervoort expression for drift velocity to curved spacetime, further refining our understanding of particle motion in complex magnetic fields [13].

These applications highlight the guiding-center model’s critical role in advancing both theoretical and experimental plasma physics, providing a comprehensive framework for studying a wide range of phenomena across different plasma environments [12–26]. Guiding center serves as an approximate model above the gyro-motion timescale and spacescale. In applications and physical explanations, guiding centers are often conflated with the motions of charged particles. A good approximate model should retain as many physical properties of the original physical system as possible, such as its conservation laws, symmetries, and geometric structures. Therefore, it is essential to develop Lagrangian theories and canonical Hamiltonian theories for the guiding center system. Since the Lagrangian formulation of guiding center was first proposed by Littlejohn in 1983 [15], researchers have been looking for its canonical Hamiltonian formulation [14, 27–30]. As an independent hierarchical model, the final link of guiding center dynamics that has yet to be completed is its canonical Hamiltonian theory.

Hamiltonian mechanics was introduced by Sir Hamilton in 1833 as a reformulation of Lagrangian mechanics. It replaces generalized velocities used in Lagrangian mechanics with generalized momenta. Although both theories provide interpretations of classical mechanics and describe the same physical phenomena, Hamiltonian theory has a closer relationship with geometry, such as symplectic geometry and Poisson structures, and serves as a crucial link between classical and quantum mechanics. The Hamiltonian formulation of guiding center system forms a classic fundamental topic in the fields of plasma physics and astrophysics [14, 15, 17, 19, 31–33].

Meiss and Hazeltine discussed the existence of the canonical coordinates of guiding center systems, but their canonical scheme is not practical in real applications [34]. White, Zakharov, and Gao studied the canonical form of guiding center motion in magnetic fields with toroidal flux-surfaces in detail [35, 36]. On the other hand, the canonical Hamiltonian formulation can be asymptotically obtained from a noncanonical Hamiltonian form of the guiding center Euler-Lagrange equation through coordinate transformation based on the Darboux-Lie theorem [37]. The main idea of previous works is to find suitable coordinates to express the Hamiltonian function in a 4-dimensional space. These canonical coordinates thus depend on the specific expression of magnetic field. In real physical systems, not all magnetic fields have flux coordinates, and it is also challenging to articulate magnetic fields in terms of flux coordinates. In fact, in most cases we cannot know the specific formulation of the magnetic field in advance. A general and rigorous canonicalization of guiding center dynamics is therefore very important. A general canonical guiding center theory can provide clear physical descriptions of the guiding center system and further explain why the properties of magnetic fields are important to the guiding center dynamics. Without its canonical Hamiltonian formulation, the existing guiding center theory would not possess a canonical symplectic geometric structure. A famous way to obtain the canonical Hamiltonian formulation from the singular Lagrangian function is Dirac’s theory of constrained Hamiltonian system [38–40]. The Legendre transformation of the singular Lagrangian function deduces the dependent canonical momentum, which implies the corresponding constraints. Then in the variables of position and the canonical momentum, a constrained Hamiltonian system can be obtained.

In this article, we apply Dirac’s theory of constrained Hamiltonian system to the new Lagrangian function of guiding center in 6D variables (\mathit{X},\dot{\mathit{X}}) obtained by employing two different approaches, and manipulate motion equations of guiding center in the 6-dimensional phase space. We rigorously prove that the guiding center dynamics can generally be described as a constrained canonical Hamiltonian system with two constraints, and that the solution flow of the guiding center lies on a canonical symplectic sub-manifold. This new canonicalization scheme is generalized and does not depend on the property of magnetic fields. The specific expression and information about magnetic fields are not required in advance. A similar application of Dirac’s theory of constraints to the Hamiltonian formulation of guiding center dynamics was also adopted in Refs. [41, 42]. Their treatments represent a direct implementation of Dirac’s framework of constrained Hamiltonian systems for Littlejohn’s Lagrangian function. The canonical form of guiding center dynamics with 6 variables presented in this paper can reduce to the original equations of motion given in Ref. [42], when restricted to the constraint manifold. In contrast to these works, our study applies Dirac’s theory to a new Lagrangian function. The advantage of our approach is that the resulting new canonical formulation of guiding center also enhances the physical comprehensions of the guiding center system. In the 6-dimensional canonical coordinates, a guiding center can thus be modeled as a pseudo-particle with its intrinsic magnetic moment, which properly replaces the charged particle dynamics on time scales larger than the gyro-period. The complete dynamical behaviors, such as the velocity and force, of the guiding center pseudo-particle can be clearly established from the canonical Hamiltonian, while two constrains play key roles. The velocity of guiding center exhibits clear origins. Its parallel velocity is constrained to the direction of the magnetic field by the two dynamical constraints. And its vertical velocity, which corresponds to the drift motion of guiding center, arises from the coupling of its intrinsic magnetic moment and the constraints in Lagrangian multipliers. As for the force acting on guiding center pseudo-particles, the two constraints supply an additional force to counteract the Larmor gyromotion. Furthermore, a series of related theories, such as symplectic numerical methods, the canonical gyro-kinetic theory, and canonical particle-in-cell algorithms can be systematically developed based on the canonical guiding center system. When extracting guiding center dynamics from charged particle dynamics in finer space-time scale, the intrinsic magnetic moment of guiding center pseudo-particle emerges naturally. This cross-scale treatment of guiding center system could shed light on the origin of the intrinsic magnetic moment of other more physically real particles, whereas in quantum physics magnetic moment and spin are considered fundamental intrinsic properties of particles and attributed to their internal symmetries.

This paper is organized as follows. In Section 2, we derive a new Lagrangian function in 6D variables (\mathit{X},\dot{\mathit{X}}), which is different from the Littlejohn’s one. Based on the new Lagrangian function, we rewrite the guiding center dynamics in 6-dimensional phase space, and obtain three equivalent canonical forms of guiding center dynamics in Section 3. In Section 4, beginning with the Hamilton−Dirac formulation of the guiding center, we derive a self-consistent framework that fully describes the guiding center’s behavior as a pseudo-particle possessing an intrinsic magnetic moment. Then to intuitively display and verify the canonical theory for guiding center pseudo-particle and its effectiveness in applications, we look into some physical systems from the numerical results obtained by the symplectic midpoint method in Section 5. Finally, in Section 6, we give a brief summary of this paper.

2 The derivation of new Lagrangian function in 6D variables

In this section, we derive the new Lagrangian function in 6D variables (\mathit{X},\dot{\mathit{X}}) following a path similar to Littlejohn’s one. The Lagrangian function of the charged particle can be written as

(1)L_{cp}=\left[{\displaystyle \frac{1}{\epsilon }}A(\mathit{x})+\mathit{v}\right]\dot{\mathit{x}}-\left[\mathrm{\Phi }(\mathit{x})+{\displaystyle \frac{1}{2}}{\mathit{v}}^2\right],

where \mathit{x} denotes the position of the charged particle, \mathit{A}(\mathit{x}) denotes the vector potential, B(\mathit{x})=||\mathit{B}(\mathit{x})||=||\mathrm{\nabla }\times \mathit{A}(\mathit{x})|| denotes the magnetic field strength, \mathit{b}(\mathit{x})=\frac{\mathit{B}(\mathit{x})}{B(\mathit{x})}=(b_1,b_2,b_3{)}^{\mathrm{T}} is the unit vector along the magnetic field, and \mathrm{\Phi }(\mathit{x}) denotes the scalar potential. Here, the mass and charge of the charged particle are normalized to unity. To derive the new Lagrangian function of guiding center, we continue to use the definition of \mathit{a} and \mathit{c} given in [15]

(2)\begin{array}{rl}\mathit{a}& =\cos \theta {\mathit{e}}_1-\sin \theta {\mathit{e}}_2,\\ \mathit{c}& =-\sin \theta {\mathit{e}}_1-\cos \theta {\mathit{e}}_2,\end{array}

where the instantaneous gyrophase \theta is defined implicitly. Let u={\mathit{v}}_{\parallel }, w={\mathit{v}}_{\perp }, and \mathit{X} denote the position of guiding center, we decompose the position of charged particle as \mathit{x}=\mathit{X}+\epsilon \frac{w\mathit{a}}{B}, which is the same as Littlejohn’s decomposition, and decompose the particle velocity to

(3)\begin{array}{rl}\mathit{v}& =\left(\mathit{b}\cdot \mathit{v}\right)\mathit{b}+w\mathit{c}\\ & =\left[\mathit{b}\dot{\mathit{X}}+\epsilon \mathit{b}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{w\mathit{a}}{B}}\right)\right]\mathit{b}+w\mathit{c}\\ & =\left(\mathit{b}\dot{\mathit{X}}\right)\mathit{b}+w\mathit{c}+O(\epsilon ),\end{array}

which is different from the velocity decomposition given by Littlejohn’s \mathit{v}=u\mathit{b}+w\mathit{c}. The three unit vector \mathit{a}, \mathit{b} and \mathit{c} satisfies \mathit{a}=\mathit{b}\times \mathit{c}. Then the Lagrangian function of charged particle can be written as

(4)\begin{array}{rl}{L}_{cp}=& \left[{\displaystyle \frac{1}{\epsilon }}A(\mathit{X})+\mathrm{\nabla }A\left({\displaystyle \frac{w\mathit{a}}{B}}\right)+(\mathit{b}\dot{\mathit{X}})\mathit{b}+w\mathit{c}+O(\epsilon )\right]\left[\dot{\mathit{X}}+\epsilon {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{w\mathit{a}}{B}}\right)\right]\\ & -\left[\mathrm{\Phi }(\mathit{X})+{\displaystyle \frac{1}{2}}(\mathit{b}\dot{\mathit{X}}{)}^2+{\displaystyle \frac{1}{2}}{w}^2+O(\epsilon )\right]\\ =& {\displaystyle \frac{1}{\epsilon }}A(\mathit{X})\dot{\mathit{X}}+\left[{\dot{\mathit{X}}}^T\mathrm{\nabla }A{\displaystyle \frac{w\mathit{a}}{B}}+(\mathit{b}\dot{\mathit{X}}{)}^2+w\mathit{c}\dot{\mathit{X}}+A(\mathit{X}){\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{w\mathit{a}}{B}}\right)\right]\\ & +\epsilon \left[(\dots )\dot{\mathit{X}}+\left(\mathrm{\nabla }A{\displaystyle \frac{w\mathit{a}}{B}}+(\mathit{b}\dot{\mathit{X}})\mathit{b}+w\mathit{c}\right)\cdot \left({\displaystyle \frac{\dot{w}\mathit{a}+w\mathit{c}\dot{\theta }}{B}}\right)\right]+O({\epsilon }^2)\\ & -\left[\mathrm{\Phi }(\mathit{X})+{\displaystyle \frac{1}{2}}(\mathit{b}\dot{\mathit{X}}{)}^2+{\displaystyle \frac{1}{2}}{w}^2+O(\epsilon )\right]\\ =& {\displaystyle \frac{1}{\epsilon }}A(\mathit{X})\dot{\mathit{X}}+\left[{\dot{\mathit{X}}}^T\mathrm{\nabla }A{\displaystyle \frac{w\mathit{a}}{B}}+(\mathit{b}\dot{\mathit{X}}{)}^2+w\mathit{c}\dot{\mathit{X}}+A(\mathit{X}){\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{w\mathit{a}}{B}}\right)\right]\\ & +\epsilon \left[(\dots )\dot{\mathit{X}}+\left({\displaystyle \frac{w^2}{B^2}}{\mathit{c}}^T\mathrm{\nabla }A\mathit{a}+{\displaystyle \frac{w^2}{B}}\right)\dot{\theta }+\left({\displaystyle \frac{w}{B^2}}{\mathit{a}}^T\mathrm{\nabla }A\mathit{a}\right)\dot{w}\right]+O({\epsilon }^2)\\ & -\left[\mathrm{\Phi }(\mathit{X})+{\displaystyle \frac{1}{2}}(\mathit{b}\dot{\mathit{X}}{)}^2+{\displaystyle \frac{1}{2}}{w}^2+O(\epsilon )\right].\end{array}

It is well known that by adding the term \frac{\mathrm{d}S}{\mathrm{d}t} with S being a scalar function to the Lagrangian function L_{cp}, the corresponding motion equation is invariant. Firstly, we choose the scalar function S_1 following Littlejohn’s way

S_1=-{\displaystyle \frac{w}{B}}\mathit{a}\cdot A(\mathit{X}),

where its derivative can be written as

(5)\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{S}_1}{\mathrm{d}t}}=& -A(\mathit{X}){\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{w\mathit{a}}{B}}\right)-{\left({\displaystyle \frac{w\mathit{a}}{B}}\right)}^{\mathrm{T}}\mathrm{\nabla }A\dot{\mathit{X}}\\ =& -A(\mathit{X}){\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{w\mathit{a}}{B}}\right)-{\dot{\mathit{X}}}^{\mathrm{T}}\mathrm{\nabla }A{\displaystyle \frac{w\mathit{a}}{B}}+{\dot{\mathit{X}}}^{\mathrm{T}}\mathrm{\nabla }A{\displaystyle \frac{w\mathit{a}}{B}}\\ & -{\left({\displaystyle \frac{w\mathit{a}}{B}}\right)}^{\mathrm{T}}\mathrm{\nabla }A\dot{\mathit{X}}\\ =& -A(\mathit{X}){\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{w\mathit{a}}{B}}\right)-{\dot{\mathit{X}}}^{\mathrm{T}}\mathrm{\nabla }A{\displaystyle \frac{w\mathit{a}}{B}}+{\displaystyle \frac{w\mathit{a}}{B}}(\mathrm{\nabla }{A}^{\mathrm{T}}-\mathrm{\nabla }A)\dot{\mathit{X}}\\ =& -A(\mathit{X}){\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{w\mathit{a}}{B}}\right)-{\dot{\mathit{X}}}^{\mathrm{T}}\mathrm{\nabla }A{\displaystyle \frac{w\mathit{a}}{B}}+{\displaystyle \frac{w\mathit{a}}{B}}\dot{\mathit{X}}\times \mathit{B}\\ =& -A(\mathit{X}){\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left({\displaystyle \frac{w\mathit{a}}{B}}\right)-{\dot{\mathit{X}}}^{\mathrm{T}}\mathrm{\nabla }A{\displaystyle \frac{w\mathit{a}}{B}}-w\mathit{c}\dot{\mathit{X}}.\end{array}

Then it can cancel the terms of O(1) in the original Lagrangian function in Eq. (4), and the remaining term of O(1) is only (\mathit{b}\dot{\mathit{X}}{)}^2. By choosing another scalar function S_2 given in [15]

S_2=-\epsilon {\displaystyle \frac{w^2}{2B^2}}{\mathit{a}}^{\mathrm{T}}\mathrm{\nabla }A\mathit{a},

with its derivative being

(6)\begin{array}{rl}{\displaystyle \frac{\mathrm{d}{S}_2}{\mathrm{d}t}}=& -\epsilon {\displaystyle \frac{w\dot{w}}{B^2}}{\mathit{a}}^{\mathrm{T}}\mathrm{\nabla }A\mathit{a}-\epsilon {\displaystyle \frac{w^2}{2B^2}}({\mathit{c}}^{\mathrm{T}}\mathrm{\nabla }A\mathit{a}+{\mathit{a}}^{\mathrm{T}}\mathrm{\nabla }A\mathit{c})\dot{\theta }+O(\epsilon )\\ =& -\epsilon {\displaystyle \frac{w\dot{w}}{B^2}}{\mathit{a}}^{\mathrm{T}}\mathrm{\nabla }A\mathit{a}-\epsilon {\displaystyle \frac{w^2}{B^2}}{\mathit{c}}^{\mathrm{T}}\mathrm{\nabla }A\mathit{a}\dot{\theta }+\epsilon {\displaystyle \frac{w^2}{2B^2}}({\mathit{c}}^{\mathrm{T}}\mathrm{\nabla }A\mathit{a}\\ & -{\mathit{a}}^{\mathrm{T}}\mathrm{\nabla }A\mathit{c})\dot{\theta }+O(\epsilon )\\ =& -\epsilon {\displaystyle \frac{w\dot{w}}{B^2}}{\mathit{a}}^{\mathrm{T}}\mathrm{\nabla }A\mathit{a}-\epsilon {\displaystyle \frac{w^2}{B^2}}{\mathit{c}}^{\mathrm{T}}\mathrm{\nabla }A\mathit{a}\dot{\theta }+\epsilon {\displaystyle \frac{w^2}{2B^2}}\mathit{a}(\mathit{c}\times \mathit{b})B\dot{\theta }\\ & +O(\epsilon )\\ =& -\epsilon {\displaystyle \frac{w\dot{w}}{B^2}}{\mathit{a}}^{\mathrm{T}}\mathrm{\nabla }A\mathit{a}-\epsilon {\displaystyle \frac{w^2}{B^2}}{\mathit{c}}^{\mathrm{T}}\mathrm{\nabla }A\mathit{a}\dot{\theta }-\epsilon {\displaystyle \frac{w^2}{2B}}\dot{\theta }+O(\epsilon ),\end{array}

we can rewrite the Lagrangian as

(7)\begin{array}{rl}{L}_{cp}^{{}^{\prime }}=L_{cp}+{\displaystyle \frac{\mathrm{d}{S}_1}{\mathrm{d}t}}+{\displaystyle \frac{\mathrm{d}{S}_2}{\mathrm{d}t}}=& {\displaystyle \frac{1}{\epsilon }}A(\mathit{X})\dot{\mathit{X}}+(\mathit{b}\dot{\mathit{X}}{)}^2+O(\epsilon )\dot{\mathit{X}}\\ & +\left[\epsilon {\displaystyle \frac{w^2}{2B}}+O({\epsilon }^2)\right]\dot{\theta }+O({\epsilon }^2)\\ & -\left[\mathrm{\Phi }(\mathit{X})+{\displaystyle \frac{1}{2}}(\mathit{b}\dot{\mathit{X}}{)}^2+{\displaystyle \frac{1}{2}}{w}^2+O(\epsilon )\right].\end{array}

In the above equation, the Lagrangian is independent of \theta, and \frac{\mathrm{\partial }{L}_{cp}^{{}^{\prime }}}{\mathrm{\partial }\dot{\theta }}=\epsilon \frac{w^2}{2B}=\epsilon \mu is an adiabatic invariant. By replacing w^2 by 2\mu B, canceling the high order terms of \epsilon, and setting \epsilon =1, we obtain a new Lagrangian function of guiding center with 6 variables (\mathit{X},\dot{\mathit{X}})

(8)\begin{array}{rl}L(\mathit{X},\dot{\mathit{X}})=& {\displaystyle \frac{1}{2}}{\dot{\mathit{X}}}^{T}M(\mathit{X})\dot{\mathit{X}}+\mathit{A}(\mathit{X})\dot{\mathit{X}}\\ & -\left[\mu B(\mathit{X})+\mathrm{\Phi }(\mathit{X})\right],\end{array}

where (\mathit{b}\dot{\mathit{X}}{)}^2={\dot{\mathit{X}}}^{T}M(\mathit{X})\dot{\mathit{X}} and the matrix M(\mathit{X}) is

M(\mathit{X})=\left(\begin{array}{ccc}{b}_1^2& b_1{b}_2& b_1{b}_3\\ b_1{b}_2& b_2^2& b_2{b}_3\\ b_1{b}_3& b_2{b}_3& b_3^2\end{array}\right).

Supposing \mathit{B}(\mathit{X})\ne 0 for magnetized plasmas, the matrix M(\mathit{q}) is singular with a constant rank \mathrm{R}\mathrm{a}\mathrm{n}\mathrm{k}(M(\mathit{X}))=1, and satisfies M^2(\mathit{X})=M(\mathit{X}) and M^T(\mathit{X})=M(\mathit{X}).

Actually, the new Lagragian function Eq. (8) can also be obtained directly from the Littlejohn’s one in Ref. [15]. The Littlejohn’s guiding center Lagrangian in the 8D variables (\mathit{X},\dot{\mathit{X}},u,\dot{u}) is

(9)\begin{array}{rl}L(\mathit{X},\dot{\mathit{X}},u,\dot{u})=& [\mathit{A}(\mathit{X})+u\mathit{b}(\mathit{X})]\cdot \dot{\mathit{X}}\\ & -\left[{\displaystyle \frac{1}{2}}{u}^2+\mu B(\mathit{X})+\mathrm{\Phi }(\mathit{X})\right].\end{array}

Here u has a ambiguous physical meaning, because it appears in guiding center Lagrangian but is usually interpreted as the parallel velocity of the original charged particle with respect to the magnetic field. Similarly, \mu denotes the magnetic moment of the original particle. The Euler−Lagrange equation concludes \mathit{b}\dot{\mathit{X}}=u. To clarify the physical model and explore the mathematical structure, we replace the variable u in Eq. (9) by \mathit{b}\dot{\mathit{X}}. Then the Lagrangian of guiding center can be rewritten in a purely 6D guiding center variables (\mathit{X},\dot{\mathit{X}}) in form of Eq. (8).

For the new Lagrangian in Eq. (8), it is still degenerate. Compared with Littlejohn’s Lagrangian function, it depends on fewer variables and is more similar to the Lagrangian of charged particle. Because \mu is an adiabatic invariant and comes from gyro-symmetry, we then define it as the intrinsic magnetic moment of the guiding center pseudo-particle, and its meaning is no longer related to the original particle.

3 Derivation of three canonical guiding center forms in 6-dimensional phase space

From the new Lagrangian function of guiding center in Eq. (8), the Hamiltonian formulation of guiding center can be obtained through Legendre transformation. If the Lagrangian is non-degenerate, i.e., \frac{{\mathrm{\partial }}^{2}L}{\mathrm{\partial }\dot{\mathit{X}}\mathrm{\partial }\dot{\mathit{X}}} is nonsingular, a canonical Hamiltonian system can be deduced using the standard Legendre map Leg(\mathit{X},\dot{\mathit{X}})=(\mathit{X},\mathit{p}), where the canonical momentum is defined by \mathit{p}=\frac{\mathrm{\partial }L}{\mathrm{\partial }\dot{\mathit{X}}}=(p_1,p_2,p_3{)}^{\mathrm{T}}. The corresponding Hamiltonian function is defined as H(\mathit{X},\mathit{p})={\mathit{p}}^{\mathrm{T}}\dot{\mathit{X}}-L(\mathit{X},\dot{\mathit{X}}(\mathit{X},\mathit{p})). For degenerate Lagrangian, generalized Hamiltonian theory has also been well developed by Dirac [38] and A. Van der Schaft [43]. The Dirac’s theory of constrained Hamiltonian system has been well developed [38–40]. When we apply the generalized Hamiltonian theory to the guiding center Lagrangian in Eq. (8), the generalized Legendre transformation leads to

(10)\mathit{p}={\displaystyle \frac{\mathrm{\partial }L}{\mathrm{\partial }\dot{\mathit{X}}}}=M(\mathit{X})\dot{\mathit{X}}+\mathit{A}(\mathit{X})\iff \begin{array}{rl}{p}_1& =b_1\mathit{b}\cdot \dot{\mathit{X}}+A_1(\mathit{X}),\\ p_2& =b_2\mathit{b}\cdot \dot{\mathit{X}}+A_2(\mathit{X}),\\ p_3& =b_3\mathit{b}\cdot \dot{\mathit{X}}+A_3(\mathit{X}).\end{array}

It outputs 3 dependent generalized momentum variables p_1,p_2 and p_3, which are connected by 3-1=2 primary constraints [43]. From Eq. (10), we conclude that \frac{p_i-A_i(\mathit{X})}{b_i(\mathit{q})} are equal for i=1,2,3, i.e., \mathit{b}\cdot \dot{\mathit{X}}=\frac{p_1-A_1(\mathit{X})}{b_1(\mathit{X})}=\frac{p_2-A_2(\mathit{X})}{b_2(\mathit{X})}=\frac{p_3-A_3(\mathit{X})}{b_3(\mathit{X})}. The two primary constraints are then expressed explicitly as

(11)\begin{array}{rl}{g}^1(\mathit{X},\mathit{p})& =b_1(\mathit{X})\left(p_2-A_2(\mathit{X})\right)-b_2(\mathit{X})\left(p_1-A_1(\mathit{X})\right)=0,\\ g^2(\mathit{X},\mathit{p})& =b_1(\mathit{X})\left(p_3-A_3(\mathit{X})\right)-b_3(\mathit{X})\left(p_1-A_1(\mathit{X})\right)=0.\end{array}

These two primary constraints force \mathit{p}-\mathit{A}(\mathit{X}) to follow the direction of magnetic field. They are steady motion constraints that depend on both position \mathit{X} and canonical momentum \mathit{p}, different from other more regular constrained Hamiltonian systems. The corresponding Hamiltonian of this non-degenerate system is

(12)\begin{array}{rl}H(\mathit{X},\mathit{p})& ={\mathit{p}}^{\mathrm{T}}\dot{\mathit{X}}-L(\mathit{X},\dot{\mathit{X}}(\mathit{X},\mathit{p}))\\ & ={\displaystyle \frac{1}{2}}{\dot{\mathit{X}}}^{\mathrm{T}}M(\mathit{X})\dot{\mathit{X}}+\mu B(\mathit{X})+\mathrm{\Phi }(\mathit{X})\\ & ={\displaystyle \frac{1}{2}}(\mathit{b}\cdot \dot{\mathit{X}}{)}^2+\mu B(\mathit{X})+\mathrm{\Phi }(\mathit{X}).\end{array}

We further express (\mathit{b}\cdot \dot{\mathit{X}}{)}^2 in Eq. (12) using coordinates (\mathit{X},\mathit{p}). Because \mathit{b}(\mathit{X}) is a unit vector, i.e., \sum _{i=1}^3{b}_i^2=1, and b_i(\mathit{b}\cdot \dot{\mathit{X}})=p_i-A_i(\mathit{X}) in Eq. (10), we can write (\mathit{b}\cdot \dot{\mathit{X}}{)}^2 in the form of (\mathit{b}\cdot \dot{\mathit{X}}{)}^2=\sum _{i=1}^3{b}_i^2(\mathit{b}\cdot \dot{\mathit{X}}{)}^2=\sum _{i=1}^3(p_i-A_i(\mathit{X}){)}^2. The Hamiltonian is finally written as

(13)\begin{array}{rl}H(\mathit{X},\mathit{p})& ={\displaystyle \frac{1}{2}}\sum _{i=1}^3(p_i-A_i(\mathit{X}){)}^2+\mathrm{\Phi }(\mathit{X})+\mu B(\mathit{X}).\end{array}

It consists of two parts corresponding to the motion of guiding center and the intrinsic magnetic moment, respectively. The canonical Hamiltonian in Eq. (13) together with the two primary constraints in Eq. (11) determines the complete guiding-center dynamics. Calculated according to [43], the dynamical equations of guiding-center as a constrained Hamiltonian system with two constraints can be obtained

(14)\{\begin{array}{l}\dot{\mathit{y}}=J{\mathrm{\nabla }}_{\mathit{y}}H(\mathit{y})+\sum _{k=1}^2{\lambda }^{k}J{\mathrm{\nabla }}_{\mathit{y}}{g}^k(\mathit{y}),(\mathit{y}=({\mathit{X}}^{\mathrm{T}},{\mathit{p}}^{\mathrm{T}}{)}^{\mathrm{T}}\in R^6),\\ 0=g^k(\mathit{y}),(k=1,2),\end{array}

where J=\left(\begin{array}{cc}0& I_3\\ -I_3& 0\end{array}\right) is the standard symplectic matrix, and {\lambda }^k(k=1,2) denotes two Lagrangian multipliers. According to Eq. (14), the motion of the guiding center can be regarded as a constrained pseudo-particle with its intrinsic magnetic moment in 6-dimensional phase space. The primary constraints g^k(k=1,2) must be preserved along the solution and satisfy \dot{g^k}={\left({\mathrm{\nabla }}_{\mathit{y}}{g}^k\right)}^T\dot{\mathit{y}}=0(k=1,2), which are also known as the Casimir functions and explicitly written as

(15)\begin{array}{rl}& {\left(g_{\mathit{X}}^1\right)}^{\mathrm{T}}{H}_{\mathit{p}}-{\left(g_{\mathit{p}}^1\right)}^{\mathrm{T}}{H}_{\mathit{X}}+{\lambda }^2\left\{g^1,g^2\right\}=0,\\ & {\left(g_{\mathit{X}}^2\right)}^{\mathrm{T}}{H}_{\mathit{p}}-{\left(g_{\mathit{p}}^2\right)}^{\mathrm{T}}{H}_{\mathit{X}}+{\lambda }^1\left\{g^2,g^1\right\}=0,\end{array}

where \{\cdot ,\cdot \} is the Poisson bracket defined by \{f,l\}=\sum _{i=1}^3\left(\frac{\mathrm{\partial }f}{\mathrm{\partial }{\mathit{X}}_i}\frac{\mathrm{\partial }l}{\mathrm{\partial }{\mathit{p}}_i}-\frac{\mathrm{\partial }f}{\mathrm{\partial }{\mathit{p}}_i}\frac{\mathrm{\partial }l}{\mathrm{\partial }{\mathit{X}}_i}\right). Without loss of generality, \left\{g^1,g^2\right\}\ne 0 holds. The Lagrangian multipliers can then be calculated from Eq. (15) as

(16)\begin{array}{rl}{\lambda }^1(\mathit{X},\mathit{p})=& {\displaystyle \frac{\left\{g^2,H\right\}}{\left\{g^1,g^2\right\}}}\\ =& {\displaystyle \frac{(p_3-A_3)\mathrm{\nabla }{b}_1\cdot (\mathit{p}-\mathit{A})-(p_1-A_1)\mathrm{\nabla }{b}_3\cdot (\mathit{p}-\mathit{A})+\mu (\mathit{b}\times \mathrm{\nabla }B{)}_2+(\mathit{b}\times \mathrm{\nabla }\mathrm{\Phi }{)}_2}{b_1[(\mathit{p}-\mathit{A})\cdot \mathrm{\nabla }\times \mathit{b}+\mathit{b}\cdot (\mathrm{\nabla }\times \mathit{A})]}},\\ {\lambda }^2(\mathit{X},\mathit{p})=& -{\displaystyle \frac{\left\{g^1,H\right\}}{\left\{g^1,g^2\right\}}}\\ =& {\displaystyle \frac{-(p_2-A_2)\mathrm{\nabla }{b}_1\cdot (\mathit{p}-\mathit{A})+(p_1-A_1)\mathrm{\nabla }{b}_2\cdot (\mathit{p}-\mathit{A})+\mu (\mathit{b}\times \mathrm{\nabla }B{)}_3+(\mathit{b}\times \mathrm{\nabla }\mathrm{\Phi }{)}_3}{b_1[(\mathit{p}-\mathit{A})\cdot \mathrm{\nabla }\times \mathit{b}+\mathit{b}\cdot (\mathrm{\nabla }\times \mathit{A})]}}.\end{array}

Substituting the expressions of {\lambda }^k(\mathit{y}) into Eq. (14), we can rewrite it as

(17)\{\begin{array}{l}\dot{\mathit{y}}=J{\mathrm{\nabla }}_{\mathit{y}}H(\mathit{y})+\sum _{k=1}^2{\lambda }^k(\mathit{y})J{\mathrm{\nabla }}_{\mathit{y}}{g}^k(\mathit{y}),\\ 0=g^k(\mathit{y}),(k=1,2).\end{array}

With the explicit expression of Lagrangian multipliers above, we can prove that the two constraints are preserving along the solution flow of Eq. (17),

\begin{array}{rl}\frac{\mathrm{d}{g}^i}{\mathrm{d}t}& ={\left(\mathrm{\nabla }{g}^i\right)}^{\mathrm{T}}\dot{\mathit{y}}\\ & ={\left(\mathrm{\nabla }{g}^i\right)}^{\mathrm{T}}\left[J\mathrm{\nabla }H+\sum _{k=1}^2{\lambda }^k(\mathit{y})J\mathrm{\nabla }{g}^k(\mathit{y})\right]\\ & ={\left(\mathrm{\nabla }{g}^i\right)}^{\mathrm{T}}J\mathrm{\nabla }H+\sum _{k=1,k\ne i}^2{\lambda }^k(\mathit{y}){\left(\mathrm{\nabla }{g}^i\right)}^{\mathrm{T}}J\mathrm{\nabla }{g}^k(\mathit{y})\\ & =0,i=1,2,\end{array}

where the last equal sign holds because the definition of {\lambda }^k(k=1,2). Then the global constraints can be replaced by the constraints at the initial point, and we obtain the equivalent guiding center dynamics

(18)\{\begin{array}{c}\dot{\mathit{X}}={\displaystyle \frac{\mathrm{\partial }H}{\mathrm{\partial }\mathit{p}}}+\sum _{k=1}^2{\lambda }^k(\mathit{X},\mathit{p}){\mathrm{\nabla }}_p{g}^k(\mathit{X},\mathit{p}),\\ \dot{\mathit{p}}=-{\displaystyle \frac{\mathrm{\partial }H}{\mathrm{\partial }\mathit{X}}}-\sum _{k=1}^2{\lambda }^k(\mathit{X},\mathit{p}){\mathrm{\nabla }}_X{g}^k(\mathit{X},\mathit{p}),\\ g^k({\mathit{y}}_0)=0,k=1,2.\end{array}

The above equation is also called Hamilton−Dirac equation based on the Dirac’s theory of constrained Hamiltonian system [38–42]. For the guiding-center dynamics written in the form (18), its solution lies on the submanifold M=\{(\mathit{y})\in R^6|g^k(\mathit{y})=0,k=1,2\}. Since the solution flow {\phi }_t:M\to M of (14) is a transformation on M, its derivative \frac{\mathrm{\partial }{\phi }_t({\mathit{y}}_0)}{{\mathit{y}}_0} is a mapping between the corresponding tangent spaces. We can prove that the the solution flow is symplectic, i.e., its Jacobi matrix satisfies

{\left({\displaystyle \frac{\mathrm{\partial }{\phi }_t({\mathit{y}}_0)}{{\mathit{y}}_0}}\right)}^{\mathrm{T}}J\left({\displaystyle \frac{\mathrm{\partial }{\phi }_t({\mathit{y}}_0)}{{\mathit{y}}_0}}\right)=J,\mathrm{f}\mathrm{o}\mathrm{r}{\mathit{y}}_0\in M.

Based on the symplecticity of the solution flow, we know that in the coordinates (\mathit{X},\mathit{p}), the solution flow of the guiding center dynamics in Eq. (18) is equipped with a canonical symplectic structure, while the 4-dimensional guiding-center dynamics has a non-canonical symplectic structure. Moreover, we can also rewrite Eq. (18) as a standard canonical Hamiltonian system form by applying the preserving of two constraints

(19)\{\begin{array}{c}\dot{\mathit{X}}={\displaystyle \frac{\mathrm{\partial }\tilde{H}}{\mathrm{\partial }\mathit{p}}},\\ \dot{\mathit{p}}=-{\displaystyle \frac{\mathrm{\partial }\tilde{H}}{\mathrm{\partial }\mathit{X}}},\\ g^k({\mathit{y}}_0)=0,k=1,2,\end{array}

where

(20)\tilde{H}(\mathit{X},\mathit{p})=H(\mathit{X},\mathit{p})+{\displaystyle \sum _{k=1}^2{\lambda }^k(\mathit{X},\mathit{p})g^k(\mathit{X},\mathit{p}).}

In this sense, we call the coordinates \mathit{X} and \mathit{p} are the canonical coordinates of guiding center. The canonical coordinates of guiding center can be obtained simultaneously in the 6-dimensional phase space independent on the specific magnetic field. The guiding center motion has 4 local degrees of freedom, since it’s on a sub-manifold in 6-dimensional space with two constraints. The properties of magnetic field, such as the existence of magnetic surface, explicitly influence the guiding center dynamics through the two constraints. Compared to the noncanonical Hamiltonian expression of the guiding center dynamics in 4-dimensional space, the constrained Hamiltonian expression in Eq. (14) has significant advantages in finding the canonical coordinates.

In this section, we obtain three equivalent canonical guiding center forms, including canonical Hamiltonian system with two constraints in Eq. (14), Hamilton−Dirac equation in Eq. (18) and canonical Hamiltonian system without constraints in Eq. (19).

4 Dynamics of the canonical guiding center pseudo-particle

The three different forms of guiding center dynamics in Eq. (14), Eq. (18) and Eq. (19) are all equivalent. However, in practical applications and simulations, we usually utilize the form of Hamilton−Dirac equation of Eq. (18), because only the constraints are required at the initial value compared to Eq. (14) and it has less terms than Eq. (19). For the guiding center dynamics Eq. (18), we can substitute the expressions of {\lambda }^1 and {\lambda }^2 in Eq. (16) and use the conservation of two constraints. Then the guiding center dynamical equations Eq. (18) can be expressed explicitly as

(21)\{\begin{array}{l}\dot{\mathit{X}}=\mathit{p}-\mathit{A}(\mathit{X})+{\displaystyle \frac{\mathit{\xi }+\mu \mathit{b}\times \mathrm{\nabla }B+\mathit{b}\times \mathrm{\nabla }\mathrm{\Phi }}{(\mathit{p}-\mathit{A})\cdot \mathrm{\nabla }\times \mathit{b}+B}},\\ \dot{\mathit{p}}=-\mu \mathrm{\nabla }B-\mathrm{\nabla }\mathrm{\Phi }+({\displaystyle \frac{\mathrm{\partial }\mathit{A}}{\mathrm{\partial }\mathit{X}}}{)}^T\dot{\mathit{X}}\\ +({\displaystyle \frac{\mathrm{\partial }\mathit{b}}{\mathrm{\partial }\mathit{X}}}{)}^{\mathrm{T}}{\displaystyle \frac{(p_1-A_1)/b_1\mathit{\xi }+(\mathit{p}-\mathit{A})\times (\mu \mathrm{\nabla }B+\mathrm{\nabla }\mathrm{\Phi })}{(\mathit{p}-\mathit{A})\cdot \mathrm{\nabla }\times \mathit{b}+B}},\\ g^k({\mathit{X}}_0,{\mathit{p}}_0)=0,k=1,2,\end{array}

where

\mathit{\xi }=\left[\begin{array}{l}({\mathit{p}}_2-{\mathit{A}}_2)\mathrm{\nabla }{\mathit{b}}_3\cdot (\mathit{p}-\mathit{A})-({\mathit{p}}_3-{\mathit{A}}_3)\mathrm{\nabla }{\mathit{b}}_2\cdot (\mathit{p}-\mathit{A})\\ ({\mathit{p}}_3-{\mathit{A}}_3)\mathrm{\nabla }{\mathit{b}}_1\cdot (\mathit{p}-\mathit{A})-({\mathit{p}}_1-{\mathit{A}}_1)\mathrm{\nabla }{\mathit{b}}_3\cdot (\mathit{p}-\mathit{A})\\ ({\mathit{p}}_1-{\mathit{A}}_1)\mathrm{\nabla }{\mathit{b}}_2\cdot (\mathit{p}-\mathit{A})-({\mathit{p}}_2-{\mathit{A}}_2)\mathrm{\nabla }{\mathit{b}}_1\cdot (\mathit{p}-\mathit{A})\end{array}\right].

This formulation based on canonical coordinates (\mathit{X},\mathit{p}) offer a complete and self-sufficient description of guiding center dynamics, which is distinct from the Littlejohn’s interpretation as the approximation of charged particle dynamics.

Based on Eq. (21), the guiding center can be regarded as a pseudo-particle in 6-dimensional phase space. The first equation of Eq. (21) determines the velocity of pseudo-particle, depicted in Fig.1. The parallel velocity V_{||}=p-A(\mathit{X}) is constrained in the direction of magnetic field by two primary constraints, where p is now the parallel canonical momentum of the guiding center. The velocity component perpendicular to the magnetic field takes a more complex form as

V_{\perp }={\displaystyle \sum _{k=1}^2{\lambda }^k(\mathit{X},\mathit{p}){\mathrm{\nabla }}_p{g}^k(\mathit{X},\mathit{p})={\displaystyle \frac{\mathit{\xi }+\mu \mathit{b}\times \mathrm{\nabla }B+\mathit{b}\times \mathrm{\nabla }\mathrm{\Phi }}{(\mathit{p}-\mathit{A})\cdot \mathrm{\nabla }\times \mathit{b}+B}}.}

Through the Lagrangian multipliers {\lambda }^k(\mathit{X},\mathit{p}), the intrinsic magnetic moment and two primary constraints interfere the vertical velocity and cause the drift motion of guiding center. In the perpendicular velocity, \frac{\mathit{\xi }}{(\mathit{p}-\mathit{A})\cdot \mathrm{\nabla }\times \mathit{b}+B} corresponds to the curvature drift, \frac{\mu \mathit{b}\times \mathrm{\nabla }B}{(\mathit{p}-\mathit{A})\cdot \mathrm{\nabla }\times \mathit{b}+B} corresponds to the gradient drift, and \frac{\mathit{b}\times \mathrm{\nabla }\mathrm{\Phi }}{(\mathit{p}-\mathit{A})\cdot \mathrm{\nabla }\times \mathit{b}+B} corresponds to the electric drift. The second equation of Eq. (21) establishes the force acting on the pseudo-particle

(22)\begin{array}{rl}\dot{\mathit{V}}=& \dot{\mathit{p}}-({\displaystyle \frac{\mathrm{\partial }A}{\mathrm{\partial }\mathit{X}}})\dot{\mathit{X}}+\dot{V_{\perp }}=-\mu \mathrm{\nabla }B(\mathit{X})-\mathrm{\nabla }\mathrm{\Phi }(\mathit{X})\\ & +\dot{\mathit{X}}\times \mathit{B}(\mathit{X})+[-{\lambda }_1(p_2-A_2)\mathrm{\nabla }{b}_1+{\lambda }_1(p_1-A_1)\mathrm{\nabla }{b}_2\\ & -{\lambda }_2(p_3-A_3)\mathrm{\nabla }{b}_1+{\lambda }_2(p_1-A_1)\mathrm{\nabla }{b}_3+\dot{V_{\perp }}],\end{array}

where, -\mu \mathrm{\nabla }B is the force related to the intrinsic magnetic moment, -\mathrm{\nabla }\mathrm{\Phi } denotes the electrostatic force, \dot{\mathit{X}}\times \mathit{B}(\mathit{X}) corresponds the magnetic force. The last term is the force caused by the two constraints, which exactly counteracts the Larmor gyromotion, as shown in Fig.3. The third equation in Eq. (21) determines the initial values of the guiding center, where the initial value of the position {\mathit{X}}_0 and the initial value of canonical momentum {\mathit{p}}_0 have to satisfy the primary constraints Eq. (10). The initial canonical momentum of the guiding center pseudo-particle is defined as {\mathit{p}}_0=b({\mathit{X}}_0)V_{||0}+A({\mathit{X}}_0), where V_{||0} is its initial parallel velocity.

So far, we have provided a general canonicalization scheme for the guiding center system without specific information of magnetic field. The guiding center canonical coordinates (\mathit{X},\mathit{p}) has been explicitly written out in a 6-dimensional phase space. The guiding center canonical Hamiltonian system is completely described by Eq. (13) together with the two primary constraints in Eq. (11). The guiding center is then modeled as a pseudo-particle independent of the original charged particle with an intrinsic magnetic moment, and its motion is restrained by two primary constraints on a submanifold. Its complete dynamics is described by the dynamical equation in Eq. (21). In canonical coordinates, the symplectic structure of the guiding center system is clear, and symplectic algorithms for guiding center are ready to be constructed. The canonical gyrokinetic equation for guiding center distribution function F(\mathit{X},\mathit{p},t) without considering collisional effects can also be easily obtained, which is

(23){\displaystyle \frac{\mathrm{\partial }F}{\mathrm{\partial }t}}+\dot{\mathit{X}}\cdot {\displaystyle \frac{\mathrm{\partial }F}{\mathrm{\partial }\mathit{X}}}+\dot{\mathit{p}}\cdot {\displaystyle \frac{\mathrm{\partial }F}{\mathrm{\partial }\mathit{p}}}=0,

where \dot{\mathit{X}} and \dot{\mathit{p}} satisfy Eq. (21). Compared with previously proposed versions, this canonical gyrokinetic equation governs the collective evolution of guiding-center pseudo-particles, which has a clear physical picture and is equivalent to the evolution of charged particles on larger time scales.

5 Validation of the canonical guiding center theory through numerical experiments

To intuitively display and verify the canonical theory for guiding center pseudo-particle and its effectiveness in applications, we look into some physical systems based on it. In these numerical experiments, we apply the symplectic mid-point method as the numerical algorithm to the system Eq. (21), and obtain the corresponding numerical results.

5.1 Example 1

Firstly, we numerically calculate the motion of a guiding center pseudo-particle and corresponding charged particle for reference, using the mid-point rule, in a given magnetic field

\mathit{B}(\mathit{X})=(y,-x,100),

with the corresponding electromagnetic potentials being \mathit{A}(\mathit{X})=(-50y,50x,\frac{x^2+y^2}{2}) and \mathrm{\Phi }(\mathit{X})=\frac{1}{2}{x}^2+\frac{1}{2}{y}^2+\frac{1}{2}{z}^2. We set the initial values of guiding center to be {\mathit{X}}_0=(0.301,0.207,-1.4) and {\mathit{p}}_0=(-10.35,15.05,0.2650), and utilize midpoint rule to the canonical form of guiding center in Eq. (21). We also apply the mid-point rule to the dynamics of charged particle at the initial values of {\mathit{x}}_0=(0.3,0.2,-1.4) and {\mathit{p}}_0=(-10.7,15.08,0.265). Fig.1 compares the orbits of the guiding center governed by Eq. (21) and its original charged particle, respectively. The trajectory of the guiding center perfectly depicts the motion of the charged particle on time scales larger than the gyroperiod. In Fig.1, it is demonstrated that \mathit{p}-\mathit{A}(\mathit{X}) evolves exactly as the parallel velocity of guiding center pseudo-particle. The Lagrangian multipliers and the constraints contribute to the perpendicular velocity of pseudo-particle \sum _{k=1}^2{\lambda }^k{\mathrm{\nabla }}_p{g}^k and cause the drift motion.

Moreover, we also apply a non-symplectic method of RK4 to the guiding center dynamics for comparison. Here, we choose the time step h=0.02, and give the comparison of the relative errors of energy obtained by midpoint rule and non-symplectic RK4 after N=6\times {10}^5 steps in Fig.2. As shown in Fig.2, the relative energy error of the midpoint rule can be confined to a very small magnitude, whereas the relative energy error of RK4 grows progressively as the number of simulated steps accumulates. It implies the advantages of symplectic methods when applied directly to the canonical form of guiding center in Eq. (21) in term of the conservation of energy over long-term simulation.

5.2 Example 2

Secondly, we calculate a guiding center trajectory in a magnetic mirror field, which is widely used to confine charged particles or plasmas for a variety of purposes [44, 45]. The magnetic mirror field in the numerical cases is set to be

(24)\begin{array}{rl}\mathit{B}(\mathit{X})=& {\displaystyle \frac{{\mu }_0{p}_m}{4\pi }}\{\left({\displaystyle \frac{2(x-a{)}^2-y^2-z^2}{[(x-a{)}^2+y^2+z^2{]}^{5/2}}}+{\displaystyle \frac{2(x+a{)}^2-y^2-z^2}{[(x+a{)}^2+y^2+z^2{]}^{5/2}}}\right)\mathit{i}\\ & +\left({\displaystyle \frac{3(x-a)y}{[(x-a{)}^2+y^2+z^2{]}^{5/2}}}+{\displaystyle \frac{3(x+a)y}{[(x+a{)}^2+y^2+z^2{]}^{5/2}}}\right)\mathit{j}\\ & +\left({\displaystyle \frac{3(x-a)z}{[(x-a{)}^2+y^2+z^2{]}^{5/2}}}+{\displaystyle \frac{3(x+a)z}{[(x+a{)}^2+y^2+z^2{]}^{5/2}}}\right)\mathit{k}\},\end{array}

where {\mu }_0 is the vacuum magnetic permeability, p_m=I\cdot S={10}^6\mathrm{A}\cdot 0.04\pi {\mathrm{m}}^2, and a=1\mathrm{m}. The magnetic axis of the magnetic mirror field is in the x direction. In this numerical example, we choose the initial position of guiding center to be {\mathit{X}}_0=(-0.4,0.2,0) \mathrm{m}, and the initial parallel velocity is V_{||0}={10}^6 \mathrm{m}/\mathrm{s}. The magnetic moment is set to be \mu =4.4458\times {10}^{-18}, and the time step of simulations is set to \mathrm{\Delta }t=5.882\times {10}^{-12} \mathrm{s}. In this condition, the guiding center is bounced back by the magnetic mirror, with its reflection point denoted by a black solid circle.

In Fig.3, we investigate the contributions of the intrinsic magnetic moment and two primary constraints in the the guiding center dynamics by depriving them in Eq. (18) respectively. If we drop both the IMM term in the canonical guiding center Hamiltonian by setting \mu =0 and the two constraints by setting {\lambda }^{1,2}(\mathit{X},\mathit{p})=0 in Eq. (18), the obtained Hamiltonian function is the same as that of the charged particle dynamics. However, because the initial values are chosen to satisfy the condition of guiding center, the corresponding trajectory particularly gives the absence of visible gyromotion, which is depicted by the orange curve in Fig.3. It implies that without the force from the intrinsic magnetic moment and two constraints, the particle gradually deviates the guiding center orbits and finally runs out of the mirror field. If we only drop the two primary constraints in the the guiding center dynamics by setting {\lambda }^{1,2}(\mathit{X},\mathit{p})=0 in Eq. (18), the corresponding system becomes a dynamics with the complete Hamiltonian but without constraints, and its trajectory approaches the guiding center orbit, but fast Larmor gyromotion still exists, see the light blue hollow circles in Fig.3. The complete canonical guiding center dynamics with constraints successfully eliminates the small timescale effects from the dynamics, see the red curve in Fig.3.

To compare the corresponding numerical results of midpoint rule and RK4, we also apply non-symplectic RK4 method to the guiding center dynamics of Eq. (21) in the magnetic mirror using dimensionless parameters. Here, we choose the parameter to be \mu =0.05, initial values of guiding center to be {\mathit{X}}_0=(-0.4,0.2,0) and V_{||0}=0.5. The time step is set to be h=0.02, and after N=5\times {10}^4 steps simulation, the comparison of the relative errors of energy obtained by midpoint rule and non-symplectic RK4 are given in Fig.4. As shown in Fig.4, the midpoint rule exhibits bounded relative energy error with minimal magnitude, while RK4 demonstrates cumulative error of energy over prolonged simulation steps, even though RK4 method is of order 4. This highlights the superiority of symplectic methods when applied directly to the canonical guiding center equations (Eq. (21)) for long-term energy conservation.

5.3 Example 3

Thirdly, to verify the canonical gyrokinetic equation in Eq. (23), we calculate the evolution of distribution function in the magnetic mirror, using both the guiding center model and charged particle model. In the magnetic mirror fields given in Eq. (24), and a radial electric field \mathit{E}=\left(0,-200y/\sqrt{y^2+z^2},-200z/\sqrt{y^2+z^2}\right), we sprinkle a total number of {10}^5 samples of both guiding centers and charged particles and follow their trajectories to solve the corresponding Vlasov equations, respectively. The calculated distribution function evolutions for both guiding centers and charged particles are shown in Fig.5.

The initial spatial distribution of guiding centers and charged particles are the same, see Fig.5(a). The initial positions for all samples in the x direction are x=-0.4\mathrm{m}. In the y−z plane, all samples are uniformly distributed in a Tai Chi region, see Fig.6(b) at the beginning. Here, the Tai Chi diagram is formed according to three circles with different radii, i.e., R_1, R_2, and R_3, where R_1 = 0.5 m, R_2 = R_1/2, R_3 = R_2/6, see Fig.6(a). The initial parallel velocity of all the samples is set to 10\mathrm{m}/\mathrm{s}, and the initial perpendicular velocity of charge particles is Maxwellian distributed with the expectation being zero and variance {\sigma }^2=kT/{\mathrm{m}}_e, where kT=4\mathrm{e}\mathrm{V}. The initial values of the adiabatic invariant of guiding center pseudo-particles can be obtained correspondingly. For charged particles, the initial gyro-phase is randomly distributed from 0 to 2\pi. The initial position of each sample in the magnetic mirror determines its constant flux surface, and the color bar for s in Fig.5 of the paper represents the corresponding normalized magnetic flux for each sample.

Due to the limitation of gyro-period, the time step for simulating charged particles is \mathrm{\Delta }t=5.882{\mathrm{e}}^{-13}\mathrm{s}. And the guiding center pseudo-particles are simulated using the time step \mathrm{\Delta }t=5.882{\mathrm{e}}^{-12}\mathrm{s}, which is 10 times larger. We plot the spatial distribution as well as corresponding normalized magnetic flux value at three moments t = 2 μs, 100 μs, 210 μs, respectively, in Fig.5(b) and (c). As shown in Fig.5, the spatial evolution of the two models are obviously consistent. When transferring canonical guiding center distributions to corresponding particle representations, the number distribution with respect to pitch angle, parallel velocity, vertical velocity, and the position in the x direction calculated by two different models are all compared and show strongly consistent, as depicted in Fig.7. In this case, the guiding center model describes the same physical distribution evolution as the charged particle model with only one-tenth of the time steps and therefore much less computational cost.

6 Conclusion

This paper introduces a comprehensive Hamiltonian-Dirac description of the guiding center dynamics based on a new Lagrangian function of guiding center in 6D variables. Guiding center system can be generally expressed using canonical coordinates in a 6D phase space. In this form, the guiding center is interpreted as a pseudo-particle with intrinsic magnetic moment and constraints independent of the particle motion. The velocity and force of the pseudo-particle are determined by canonical dynamical equations of the pseudo-particle, where two primary constraints are key to eliminate fast-scale gyromotion. Geometric structures, symplectic algorithms, and canonical gyrokinetic theories are readily developed based on the generalized canonicalization of guiding center systems. In the future work, we will continue developing the theory of canonical gyro-kinetics and the canonical gyrokinetic PIC algorithms based on the canonical guiding center theory. We will study the geometric structures and structure-preserving algorithms of guiding center system and their applications to real plasma problems. Furthermore, this canonicalization scheme endows the pseudo particle with an intrinsic magnetic moment by restraining small scale dynamics, which is also an interesting topic and worth exploring.

References

[1] G. G. Howes, J. M. TenBarge, W. Dorland, E. Quataert, A. A. Schekochihin, R. Numata, and T. Tatsuno, Gyrokinetic simulations of solar wind turbulence from ion to electron scales, Phys. Rev. Lett. 107(3), 035004 (2011)

[2] J. W. Schroeder, G. Howes, C. Kletzing, F. Skiff, T. Carter, S. Vincena, and S. Dorfman, Laboratory measurements of the physics of auroral electron acceleration by Alfvén waves, Nat. Commun. 12(1), 3103 (2021)

[3] L. Wang and P. Diamond, Gyrokinetic theory of turbulent acceleration of parallel rotation in Tokamak plasmas, Phys. Rev. Lett. 110(26), 265006 (2013)

[4] L. Chen, F. Zonca, and Y. Lin, Physics of kinetic Alfvén waves: A gyrokinetic theory approach, Rev. Mod. Plasma Phys. 5(1), 1 (2021)

[5] G. Dif-Pradalier, P. Ghendrih, Y. Sarazin, E. Caschera, F. Clairet, Y. Camenen, P. Donnel, X. Garbet, V. Grandgirard, Y. Munschy, L. Vermare, and F. Widmer, Transport barrier onset and edge turbulence shortfall in fusion plasmas, Commun. Phys. 5(1), 229 (2022)

[6] S. Maeyama, T. H. Watanabe, M. Nakata, M. Nunami, Y. Asahi, and A. Ishizawa, Multi-scale turbulence simulation suggesting improvement of electron heated plasma confinement, Nat. Commun. 13(1), 3166 (2022)

[7] Y. Na, J. Seo, Y. Lee, G. Choi, M. Park, S. Park, S. Yi, W. Wang, M. Yoo, M. Cha, B. Kim, Y. H. Lee, H. Han, B. Kim, C. Lee, S. K. Kim, S. M. Yang, C. S. Byun, H. S. Kim, J. Ko, W. Lee, and T. S. Hahm, Observation of a new type of self-generated current in magnetized plasmas, Nat. Commun. 13(1), 6477 (2022)

[8] Z. Qiu, L. Chen, and F. Zonca, Gyrokinetic theory of toroidal Alfvén eigenmode saturation via nonlinear wave–wave coupling, Rev. Mod. Plasma Phys. 7(1), 28 (2023)

[9] J. Dubois, S. Berman, C. Chandre, and T. Uzer, Capturing photoelectron motion with guiding centers, Phys. Rev. Lett. 121(11), 113202 (2018)

[10] D. H. Dubin, Electronic and positronic guiding-center drift ions, Phys. Rev. Lett. 92(19), 195002 (2004)

[11] A. Bañón Navarro, P. Morel, M. Albrecht-Marc, D. Carati, F. Merz, T. Görler, and F. Jenko, Free Energy Cascade in Gyrokinetic Turbulence, Phys. Rev. Lett. 106(5), 055001 (2011)

[12] C. Wang, M. Gong, X. Zhao, Q. W. Nan, X. Y. Yu, Y. Cheng, V. Kimberg, X. J. Liu, O. Vendrell, K. Ueda, and S. B. Zhang, Rebuilding the vibrational wavepacket in TRAS using attosecond X-ray pulses, Commun. Phys. 7, 1 (2024)

[13] D. T. Son and M. Stephanov, Relativistic guiding-center motion: Action principle, kinetic theory and hydrodynamics, Phys. Rev. Lett. 133(14), 145201 (2024)

[14] R. G. Littlejohn, A guiding center Hamiltonian: A new approach, J. Math. Phys. 20(12), 2445 (1979)

[15] R. G. Littlejohn, Variational principles of guiding centre motion, J. Plasma Phys. 29(1), 111 (1983)

[16] A. Hasegawa and Y. Kodama, Guiding-center soliton, Phys. Rev. Lett. 66(2), 161 (1991)

[17] H. Qin, W. Tang, W. Lee, and G. Rewoldt, Gyrokinetic perpendicular dynamics, Phys. Plasmas 6(5), 1575 (1999)

[18] H. Qin, W. Tang, and W. Lee, Gyrocenter-gauge kinetic theory, Phys. Plasmas 7(11), 4433 (2000)

[19] H. Sugama, Gyrokinetic field theory, Phys. Plasmas 7(2), 466 (2000)

[20] A. J. Brizard and T. S. Hahm, Foundations of nonlinear gyrokinetic theory, Rev. Mod. Phys. 79(2), 421 (2007)

[21] J. C. W. Song, G. Refael, and P. A. Lee, Linear magnetoresistance in metals: Guiding center diffusion in a smooth random potential, Phys. Rev. B 92(18), 180204 (2015)

[22] G. G. Howes, J. M. TenBarge, W. Dorland, E. Quataert, A. A. Schekochihin, R. Numata, and T. Tatsuno, Gyrokinetic simulations of solar wind turbulence from ion to electron scales, Phys. Rev. Lett. 107(3), 035004 (2011)

[23] C. S. Chang, S. Ku, G. R. Tynan, R. Hager, R. M. Churchill, I. Cziegler, M. Greenwald, A. E. Hubbard, and J. W. Hughes, Fast low-to-high confinement mode bifurcation dynamics in a tokamak edge plasma gyrokinetic simulation, Phys. Rev. Lett. 118(17), 175001 (2017)

[24] J. W. Burby, Some mathematics for quasi-symmetry, J. Math. Phys. 61, 093503 (2020)

[25] T. Rubin, J. Rax, and N. Fisch, Guiding centre motion for particles in a ponderomotive magnetostatic end plug, J. Plasma Phys. 89(6), 905890615 (2023)

[26] A. J. Brizard, Particle and guiding-center orbits in crossed electric and magnetic fields, Phys. Plasmas 30(4), 042113 (2023)

[27] R. G. Littlejohn, Hamiltonian formulation of guiding center motion, Phys. Fluids 24(9), 1730 (1981)

[28] R. G. Littlejohn, Hamiltonian perturbation theory in noncanonical coordinates, J. Math. Phys. 23(5), 742 (1982)

[29] N. Tronko and A. J. Brizard, Lagrangian and Hamiltonian constraints for guiding-center Hamiltonian theories, Phys. Plasmas 22(11), 112507 (2015)

[30] A. Neishtadt,A. Artemyev, Hamiltonian for guiding center motion: Symplectic structure approach, arXiv: 1905.12316 (2019)

[31] J. R. Cary and A. J. Brizard, Hamiltonian theory of guiding-center motion, Rev. Mod. Phys. 81(2), 693 (2009)

[32] S. Middelkamp, P. Torres, P. Kevrekidis, D. Frantzeskakis, R. Carretero-González, P. Schmelcher, D. Freilich, and D. Hall, Guiding-center dynamics of vortex dipoles in Bose-Einstein condensates, Phys. Rev. A 84(1), 011605 (2011)

[33] J. Xiao and H. Qin, Slow manifolds of classical Pauli particle enable structure-preserving geometric algorithms for guiding center dynamics, Comput. Phys. Commun. 265, 107981 (2021)

[34] J. Meiss and R. Hazeltine, Canonical coordinates for guiding center particles, Phys. Fluids B 2(11), 2563 (1990)

[35] R. Whiteand L. E. Zakharov, Hamiltonian guiding center equations in toroidal magnetic configurations, Phys. Plasmas 10(3), 573 (2003)

[36] H. Gao, J. Yu, J. Nie, J. Yang, G. Wang, S. Yin, W. Liu, D. Wang, and P. Duan, Hamiltonian canonical guiding center equations constructed under a nested magnetic flux coordinates system, Phys. Scr. 86(6), 065502 (2012)

[37] R. Zhang, J. Liu, Y. Tang, H. Qin, J. Xiao, and B. Zhu, Canonicalization and symplectic simulation of the gyrocenter dynamics in time-independent magnetic fields, Phys. Plasmas 21(3), 032504 (2014)

[38] P. A. M. Dirac, Generalized Hamiltonian dynamics, Canadian Journal of Mathematics 2, 129 (1950)

[39] P. A. M. Dirac, Generalized Hamiltonian dynamics, Proc. R. Soc. Lond. A 246(1246), 326 (1958)

[40] K. Sundermeyer, Constrained Dynamics, Lecture Notes in Physics, Vol. 169, Springer-Verlag, New York, 1982

[41] D. Pfirsch and P. J. Morrison, The energy-momentum tensor for the linearized Maxwell–Vlasov and kinetic guiding center theories, Phys. Fluids B 3(2), 271 (1991)

[42] D. A. Kaltsas, Constrained Hamiltonian systems and physics-informed neural networks: Hamilton-Dirac neural networks, Phys. Rev. E 111(2), 025301 (2025)

[43] A. J. Schaft, Equations of motion for Hamiltonian systems with constraints, J. Phys. Math. Gen. 20(11), 3271 (1987)

[44] G. Gibson, W. C. Jordan, and E. J. Lauer, Containment of Positrons in a Mirror Machine, Phys. Rev. Lett. 5(4), 141 (1960)

[45] E. Anderson, C. Baker, W. Bertsche, N. M. Bhatt, G. Bonomi, . Observation of the effect of gravity on the motion of antimatter, Nature 621(7980), 716 (2023)

RIGHTS & PERMISSIONS

Higher Education Press