Analytical Calculation of the Orbital Spectrum of the Guiding Center Motion in Axisymmetric Magnetic Fields

Charged particle motion in axisymmetric toroidal magnetic fields is analyzed within the context of the canonical Hamiltonian Guiding Center theory. A canonical transformation to variables measuring the drift orbit deviation from a magnetic field line is introduced and an analytical transformation to Action-Angle variables is obtained, under a zero drift width approximation. The latter is used to provide compact formulas for the orbital spectrum of the drift motion, namely the bounce/transit frequencies as well as the bounce/transit averaged toroidal precession and gyration frequencies. These formulas are shown to have a remarkable agreement with numerically calculated full drift width frequencies and significant differences with standard analytical formulas based on a pendulum-like Hamiltonian description. The analytical knowledge of the orbital spectrum is crucial for the formulation of particle resonance conditions with symmetry breaking perturbations and the study of the resulting particle, energy and momentum transport.


I. INTRODUCTION
Charged particle dynamics in toroidal magnetic fields has been the key theoretical issue for the study of magnetically confined fusion plasmas for many decades. The understanding of the role of the magnetic field topology on single as well as on collective particle dynamics has been crucial for the design of fusion devices with magnetic fields having different types of symmetries [1]. These background magnetic fields also determine the interactions between particles and symmetry-breaking perturbations of a large range of spatial and temporal scales, mainly through resonance conditions, and the corresponding particle, energy and momentum transport.
The Guiding Center (GC) theory [2] has been used for a rigorous dynamical reduction by the systematic elimination of the rapidly varying gyro-angle variable related to the cyclotron motion around a magnetic field line. In the GC description, the magnetic moment is a constant of the motion and the particle dynamics are studied in terms of the drift motion of the center of the cyclotronic motion, leading to a gyro-kinetic Hamiltonian theory. Although the original derivation of the GC equations has been formulated in non-canonical variables, the utilization of magnetic coordinates has been shown to allow a Hamiltonian formulation in canonical variables [3][4][5].
The Hamiltonian formulation of the guiding center motion in canonical variables has several computational and conceptual advantages that are revealed under a transformation to Action-Angle variables [6,7]. Such a transformation is possible for the case of GC motion in an axisymmetric magnetic field where the Hamiltonian system is integrable. In the Action-Angle variable set, the topology of the motion is described by multi-dimensional tori and each orbit is labeled by a distinct set of the invariant values of the three Action variables [8]. This simple orbit parametrization allows for an orbit-based analysis of particle, energy and momentum transport which is particularly useful for the study of energetic particle dynamics in fusion plasmas in direct relation to velocity-space tomography techniques [9,10]. Moreover, all the orbital frequencies can be readily calculated in terms of the Action variables. The latter determine the resonance conditions between particles and symmetrybreaking perturbations [11,12] resulting in breaking of the Action invariance and diffusion in the Action space, that describes energy, momentum and radial particle transport [13][14][15][16].
The collective particle dynamics, under the presence of perturbations is characterized by the modification of the unperturbed particle distribution functions [17][18][19]. Another important feature of the Action-Angle description is that the different time scales of the motion are well separated in different degrees of freedom allowing for a systematic dynamical reduction to a hierarchy of evolution equations for the reduced distribution functions [15,20].
In this work, we utilize the canonical Hamiltonian formulation of the GC motion along with an appropriate canonical transformation in order to facilitate the calculation of Action-Angle variables for the case of an axisymmetric magnetic field. By doing so we fully exploit the advantages of the canonical Hamiltonian formulation in terms of the calculation of the orbital frequencies of all degrees of freedom corresponding to bounce/transit, bounceaveraged toroidal precession and gyration frequencies, as well as the calculation of the Action variables allowing for the dynamical reduction to a bounce-averaged system. Analytical results are obtained under a Zero Drift Width (ZDW) approximation providing compact formulas for the orbital frequencies that are in remarkable agreement with Full Drift Width (FDW) numerical calculations and have significant qualitative and quantitative differences with standard analytical formulas based on a pendulum-like GC Hamiltonian.
In Section 2, the canonical formulation of the GC theory is briefly presented for completeness. In Section 3, we introduce a canonical transformation to variables measuring Drift Orbit Deviation from a given field line. The transformation is general and applies to both axisymmetric and non-axisymmetric equilibria of arbitrary shape. We show that, for a Large Aspect Ratio (LAR) magnetic field equilibrium, the transformation takes a particularly simple form and we briefly comment on possibilities of treating equilibria with higher-order terms with respect to the inverse aspect ratio. The general form of the transformation to Action-Angle variables as well as the corresponding calculation of the orbital frequencies is also presented. In Section 4, we apply a ZDW approximation for the LAR equilibrium, according to which the GC orbit is considered to take place on a single flux surface, and we obtain simple analytical formulas for the Orbital Spectrum of the GC motion. We present a novel ZDW Hamiltonian retaining terms that are significant for particles with smaller pitch angles in comparison to standard pendulum-like Hamiltonians describing deeply trapped particles, and we compare the analytical results to FDW numerical calculations. The summary and conclusions are given in Section 5.

RIC EQUILIBRIUM
A general axisymmetric toroidal magnetic configuration consisting of nested toroidal flux surfaces can be represented in White-Boozer [5] coordinates as where ζ and θ are the toroidal and the poloidal angles. The toroidal flux ψ is related to the poloidal flux ψ p through the safety factor q(ψ) = dψ/dψ p . The functions g and I are related to the poloidal and toroidal currents and δ is related to the nonorthogonality of the coordinate system.
The guiding center motion of a charged particle is described by the Lagrangian where A and B are the vector potential and the magnetic field, v is the guiding center velocity, µ is the magnetic moment, ξ is the gyrophase, ρ is the velocity component parallel to the magnetic field, normalized to B, and is the Hamiltonian. The guiding center motion is given in normalized units where time is normalized to ω −1 0 , with ω 0 = eB 0 /m being the on-axis gyrofrequency, and distance is normalized to the major radius R, so that energy is normalized to mω 2 0 R 2 . According to the ordering of the guiding center approximation, the gyroradius is ρ = v/B << 1 and the magnetic moment µ = v 2 ⊥ /(2B) as well as the cross field drift are of order ρ 2 [5]. The three couples of canonically conjugate variables for this GC Hamiltonian are (µ, ξ), (P θ , θ) and (P ζ , ζ) [5] with providing the relation between the canonical momenta (P θ , P ζ ) and (ψ, ρ ) [5]. In terms of these canonical variables, Eq. (2) can be written as H(P θ , θ, P ζ , ζ, µ, ξ) = [P ζ + ψ p (P θ , P ζ )] 2 2g 2 (ψ(P θ , P ζ )) B 2 (ψ(P θ , P ζ ), θ) + µB(ψ(P θ , P ζ ), θ).
The gyroangle ξ does not appear in the guiding center Hamiltonian which is also independent of ζ due to axisymmetry of the magnetic field; therefore the corresponding canonical momenta, namely µ and P ζ , are constants of the motion and since the Hamiltonian does not depend explicitely on time (autonomous system), the system is integrable.
III. CANONICAL TRANSFORMATION TO DRIFT ORBIT DEVIATION VARI-

ABLES
The guiding center Hamiltonian describes all particle drifts due to the inhomogeneity of the magnetic field, causing the guiding center deviation from a field line. A canonical transformation with generating function [6] transforms to a new (barred) variable set, related to the original variables as This canonical tranformation is general and can be applied for either axisymmetric or non-axisymmetric equilibrium magnetic fields. The physical meaning of the new canonical variables becomes obvious for a Large Aspect Ratio (LAR) cylindrical equilibrium described by g = 1, I = 0, δ = 0 and B = 1 − r cos θ where r = √ 2ψ and the magnetic field is normalized to its on-axis value [3]. In this case, we have ∂ψ/∂P ζ = 0 and the original canonical momenta, as defined by Eq. (3), become P θ = ψ and P ζ = ρ − ψ p (P θ ), so that P ζ = ρ andζ = ζ. The new variables (P θ ,θ) provide the deviation of the guiding center orbit from a magnetic field line of reference, for which θ = ζ/q(ψ), intersecting the poloidal plane ζ = 0 at (ψ 0 , θ 0 ) = (P θ0 ,θ). The new canonical angleθ is directly related to the toroidal precession, defined asζ − qθ. In contrast to previous approaches [21] where an a-posteriori construction of a canonically conjugate pair of variables (ρ , s) for the LAR equilibrium case is necessary, within the context of our formulation, the canonical variables are a priori defined. Moreover, this also holds, not only for the LAR case but also, for an arbitrary equilibrium magnetic field.
For a magnetic field configuration that can be considered as a higher order perturbation of the LAR equilibrium, according to the standard tokamak ordering, we can write ψ = P θ + ǫ 2ψ (P θ , P ζ ), with ǫ being the inverse aspect ratio andψ(P θ , P ζ ) corresponding to corrections due to ellipticity and triangularity [5]. Therefore, for the new variableζ we have and to lowest orderζ = ζ, whereas the constant of the motion can be written as and to lowest orderP For the case of a smooth q profile we can neglect (locally) higher order derivatives of ψ p with respect to ψ = P θ , which is equivalent to neglecting the radial variation of the safety factor q (magnetic shear) within a guiding center drift orbit, and applying a Taylor expansion in with prime denoting differentiation of ψ p with respect to its argument. Since P ζ is constant, we can writeP with ρ 0 = P ζ + ψ p (P θ0 ) depending on both the invariant P ζ and the flux surface of reference related to P θ0 . This relation imposes a constraint on the two canonical momentaP ζ and P θ due to the axisymmetry of the magnetic field and the corresponding invariance of P ζ .
It is worth mentioning that for the case of drift orbits spanning a larger radial distance, we can also include the second order derivative ψ ′′ p (P θ0 ), related to the magnetic shear (q ′ ), that would result in a quadratic term with respect toP θ in (10) and (11). As a result the magnetic shear can modify the guiding center dynamics and orbital spectrum [22,23].
In the new canonical variables the Hamiltonian is given as It is clear that in the new variables there is no cyclic angle, so that neitherP θ norP ζ is a constant of the motion. However, the quantity P ζ is still a constant of the motion, due to axisymmetry, restrictingP θ andP ζ according to Eq. (11), so that integrability is preserved.
By substituting either of the canonical momenta as a function of the other from (11), the Hamiltonian system can be readily described in one degree of freedom either as H(P ζ ,ζ; ρ 0 , P θ0 ;μ) or as H(P θ ,θ; ρ 0 , P θ0 ;μ). In the former case, since the Hamiltonian does not depend onP θ , its canonically conjugate angleθ appears as an additive phase constant that can be omitted, whereas the same holds forζ in the latter case. This duality directly relates the radial and poloidal drift with the toroidal momentum and precession.
The Hamiltonian (12) provides a canonical description of the guiding center motion for a generic axisymmetric magnetic field equilibrium to the lowest order with respect to the standard tokamak ordering and drift orbit deviation from a magnetic field line. Higher order magnetic field equilibria, with respect to the inverse aspect ratio ǫ, can be considered by keeping higher order terms in (8). The one degree of freedom Hamiltonian H(P ζ ,ζ; ρ 0 , P θ0 ; µ) can be readily used to calculate the Action as well as the respective frequenciesω with E = H being the constant energy value of each orbit and the index (b, t) corresponding to the cases of trapped (bounce) or passing (transit) orbits. The mixed-variable generating provides the canonical transformation to Action-Angle variables also for the remaining canonical variable pairs (μ,ξ) → (J ξ ,ξ) It is worth emphasizing that in the other two pairs of canonical variables, the new momenta (Actions) are identical to the old momenta, i.e. J θ =P θ and J ξ =μ = µ but the new positions (angles) differ from the old ones, due to the last term of the generating function (15) [8]. The transformation to Action-Angle variables allows for the calculation of the orbital frequencies in the remaining two degrees of freedom. Therefore, The above equations show that the bounce/transit Actions J The frequencyω θ is directly related to the bounce/transit averaged toroidal drift as shown in the following expression where ∆θ T ζ =2π/ω ζ is the variation ofθ = θ−ζ/q in the time interval of a period T ζ = 2π/ω ζ .
Similarly, the frequencyω ξ corresponds to the bounce/transit averaged gyrofrequency [8] ω The Zero Drift Width (ZDW) approximation, under which the guiding center is considered as being fixed on a given flux surface P θ0 = ψ 0 , can be readily obtained by settinḡ P θ = 0. Moreover, for a relatively small drift widthP θ << P θ0 , ψ can be Taylor-expanded around P θ0 and withP θ substituted from (11) we have a single degree of freedom Hamiltonian for the canonical variables (P ζ ,ζ).

For a LAR equilibrium the Full Drift Width (FDW) Hamiltonian is written as
where we have substitutedP ζ = ρ and dropped bar in ζ for simplicity. It is worth emphasizing that obtaining the Hamiltonian (23) is enabled by the utilization of the normalized parallel guiding-center velocity ρ (normalized to the magnetic field), instead of the regular parallel guiding-center velocity used in previous works [25]. Moreover, according to Eq.

IV. ZERO DRIFT WIDTH APPROXIMATION
We consider the Zero Drift Width approximation under which the guiding center motion is considered fixed on a given flux surface ψ = ψ 0 and the frequencies of θ and ζ are simply related as ω ζ = q(ψ 0 )ω θ . The Hamiltonian is written as or equivalently Further approximations can be based on the relative magnitude of the three terms in the curly brackets of the above equation. For ρ 2 /µ << 1, corresponding to particles with large pitch angles α = tan −1 (v ⊥ /|v |), the Hamiltonian is reduced to which resembles the Hamiltonian of a pendulum, with the trapped and passing orbits corresponding to a libration and rotation type of motion [21,25,26]. It is worth mentioning that, under the ZDW approximation, both H ZDW and H ′ ZDW scale with µ when ρ 2 is also divided by µ, which clearly is not the case for the FDW Hamiltonian [Eq.(23)].
The above approximations have significant quantitative and qualitative differences: They provide different values for the total energy of a specific orbit and also result in different separatrices betweeen trapped and passing orbits in the phase space, as shown in Fig. 1.
Therefore, a particle that is described as being trapped according to H ′ ZDW can be actually passing according to H ZDW and vice versa. Both ZDW Hamiltonians describe guiding center orbits that are symmetric with respect to ρ = 0 whereas this is not the case with the FDW Hamiltonian, according to which positive (co-passing) and negative (counter-passing) orbits are not symmetric. These differences are not uniform across the phase space of the system and depend strongly on the pitch angle and the flux surface of reference ψ 0 .
The ZDW Hamiltonians are very useful for obtaining analytical forms for the frequencies of the guiding center motion and their dependence on the constants of the motion, parametrizing each orbit. These frequencies determine the Orbital Spectrum (OS) of different particle species, including bulk and energetic particles, as well as their resonance conditions with any type of non-axisymmetric perturbations. The latter are crucial for the energy, momentum and particle transport in a toroidal magnetic field configuration.
For both H ZDW and H ′ ZDW bounce and transit motion is characterized by 0 ≤ k < 1 and 1 < k, respectively, with k is the trapping parameter, defined as and The bounce/transit frequencies and Actions corresponding to H ′ ZDW are given as follows (Appendix). and where K and E are the complete elliptic integrals of the first and second kind [27], and is the frequency of the deeply trapped bounce motion, corresponding to k = 0. These are the standard analytical expressions for the frequencies and Actions obtained from the pendulum-like Hamiltonian H ′ ZDW [5,25,26]. The respective frequencies and Actions for H ZDW are (Appendix) and with Π being the complete elliptic integral of the third kind [27], and For a given set of values of the magnetic moment (µ) and the flux surface (ψ 0 ), the value of the trapping parameter (k) is determined by the total energy (E) of the particle.
The dependence of the bounce and transit frequencies on the energy E(k) according to Hamiltonian in the original canonical variable set, given in Eq. (4), are also shown. Note that under the ZDW approximationθ =ζ/q =ζ/q = ω b,t /q =ω ζ /q, and we have taken q = 1 for simplicity. The analytical formula for ω b (35) shows a remarkable agreement with the numerically calculated frequencies based on the FDW Hamiltonian and significantly deviate from the analytical formula for ω ′ b (30) with the deviation being more pronounced for larger values of ψ 0 , as shown in Fig. 2. For the transit frequencies, shown in Fig. 3, the FDW Hamiltonian describes asymmetric orbits and consequently different frequencies for co-passing (ρ > 0) and counter-passing (ρ < 0) orbits. The analytical formula for ω t (36) corresponds to an intermediate frequency with respect to the two branches of the numerically calculated FDW frequencies, whereas the analytical formula for ω ′ t (31) tends to follow one of the two branches.
The two analytical formulas for the bounce freqencies have also a significant qualitative difference. In fact, in contrast to ω ′ b , ω b is non-monotonic with respect to k (energy) and, in accordance to FDW numerical calculations, predicts a local maximum of the frequency of the trapped orbits with respect to the energy, indicating the existence of two trapped orbits with different energy but equal frequency, as shown in Fig. 2(b). The location of this local maximum in the space (E(k), ψ 0 ) is depicted in Fig. 3. The nonmonotonic dependence of the frequency on the particle energy has important implications for particle and momentum transport under the presence of non-axisymmetric perturbations that break the integrability of the system. In such cases, the locations of the phase space regions where orbits are significantly modified due to perturbations are determined by resonance conditions. The nonmonotonicity suggests that the same resonance can modify two regions of the phase space leading to extended transport under resonance overlap conditions. Moreover, the existence of energy values where the derivative of the frequency with respect to energy is zero can be related to an intrinsic degeneracy condition which is crucial for the effect of perturbations [7].
The bounce-averaged toroidal precession and gyration frequencies can be analytically caclulated according to Eqs. (19) and (20)  The knowledge of the Orbital Frequencies is crucial for determining the resonance conditions under particle interaction with non-axisymmetric perturbations that affect energy, momentum and particle transport in toroidal plasma configurations. Moreover, the transformation to Action-Angle variables allows for the application of standard canonical perturbation methods as well as the systematic dynamical reduction and the formulation of a bounce gyrkinetic description.
The equation of motion dζ/dt = ∂H/∂ρ can be written as where χ = ζ/(2q(ψ 0 )). The integration of the above relation over a closed orbit provides the period of the bounce and the transit motion as with the respective frequency given as ω = 2π/T .
where we have takenζ = 0 when χ = χ min . Note that this relation defines the transformation from ζ = 2q(ψ 0 )χ toζ in an implicit form. For the case of the Hamiltonian H ′ ZDW corresponding to σ = 0, the relation can be inverted with the use of Jacobi elliptic functions, whereas for the case of H ZDW corresponding to σ = 1, the inversion of the relation, to the best of our knowledge, cannot be expressed in a convenient form [27].