aa r X i v : . [ phy s i c s . p l a s m - ph ] D ec Automation of The Guiding Center Expansion
J. W. Burby, J. Squire, and H. Qin
1, 2 Princeton Plasma Physics Laboratory, Princeton, New Jersey 08543,USA Dept. of Modern Physics, University of Science and Technology of China, Hefei,Anhui 230026, China (Dated: 12 September 2018)
We report on the use of the recently-developed Mathematica package
VEST (VectorEinstein Summation Tools) to automatically derive the guiding center transforma-tion. Our Mathematica code employs a recursive procedure to derive the transfor-mation order-by-order. This procedure has several novel features. (1) It is designedto allow the user to easily explore the guiding center transformation’s numerous non-unique forms or representations. (2) The procedure proceeds entirely in cartesianposition and velocity coordinates, thereby producing manifestly gyrogauge invariantresults; the commonly-used perpendicular unit vector fields e , e are never even in-troduced. (3) It is easy to apply in the derivation of higher-order contributions tothe guiding center transformation without fear of human error. Our code thereforestands as a useful tool for exploring subtle issues related to the physics of toroidalmomentum conservation in tokamaks.1 . INTRODUCTION The guiding center asymptotic expansion is both beautiful and revolting. Its beautystems from its simple physical underpinning; a strongly magnetized charged particle gyratesaround magnetic field lines much more rapidly than it drifts along or across them. Thissimplicity allows the approximation to be applied in a greater variety of settings than perhapsany other approximation scheme used in magnetized plasma physics. And in spite of theapproximation’s broad appicability, which might be expected to dilute its power, it affordssignificant practical benefits. Perhaps most notably, it enables gyrokinetic codes, such asthose discussed in Refs. 1 and 2, to work on the drift, rather than gyroperiod, time scale.The approximation begins to reveal its ugly side, however, when one endeavors to derivesuccessively higher-order contributions to the expansion . Aside from the usual prolifera-tion of terms common amongst higher-order perturbation expansions, the obstacles one en-counters include vector identities involving spatially varying unit vectors such as b = B / | B | and subtle issues related to gyrogauge invariance . Moreover, attempts to taylor the expan-sion to respect the Hamiltonian structure of the Lorentz force law encounter the so-called order-mixing issue, whereby different components of the coordinate transformation oneseeks appear at different orders in the transformed Lagrangian, thus complicating the pro-cedure used to find them.These abhorrent features can be frightening to the uninitiated. As a result, only a ded-icated minority have ever attempted delving into the calculation beyond the derivation ofdrifts proportional to first derivatives of the magnetic field. The reluctant majority, up untilfairly recently , could have justified their stance by proclaiming the higher-order correctionsto be practically unimportant, and therefore irrelevant. Recent advances, however, are mak-ing it more and more clear that at least corrections proportional to second derivatives ofthe magnetic field are important for resolving the physics of toroidal momentum conserva-tion in tokamaks . For this reason, certain largely unexplored aspects of these higher-ordercorrections now appear intriguing to study. In particular, the various representations of theguiding center expansion should be explored further.A representation of the guiding center expansion consists of a prescription for making allof the apparently arbitrary choices one must make in the process of deriving the expansion.Examples of different representations can be found In Ref. 4, where two representations are2resented, or in Littlejohn’s work in Refs. 9 and 10. There is nothing unphysical about thesedifferent representations - they merely arise from the fact that equations of motion whichare independent of gyrophase will remain so upon an arbitrary coordinate transformationthat commutes with the gyrosymmetry operation (see appendix B). Nevertheless, differentrepresentations lead to guiding center equations of motion with different numbers of terms.Thus, one could imagine optimizing the number of terms in the equations of motion overthe space of representations. It is also possible that different representations have differenttimes of validity. After all, Kruskal’s method , which provides the mathematical basis forthe guiding center expansion, can only guarantee equations of motion valid for times of order1 /ǫ , where ǫ is the ordering parameter ρ/L .In order to enable the study of these issues, a process which would surely involve de-riving the guiding center expansion in many different representations, we have developed,implemented, and verified an algorithm to automate the guiding center calculation usingthe newly-developed Mathematica package VEST (Vector Einstein Summation Tools) . Inparticular, we have slashed the time required to derive the expansion, and all but elim-inated the possible taint of human-made algebra errors in the derivation of higher-ordercontributions to the guiding center expansion.While other authors have presented algorithmic procedures for deriving the guiding centerexpansion in the past , the algorithm we present here is novel due to the combinationof the following.1) The algorithm has actually been implemented on a computer and used to derive theguiding center expansion in two different representations.2) Complicated, multi-term, vector identities are accounted for using the clever simplificationcapabilities of VEST .3) Issues related to gyrogauge invariance are completely avoided by working in cartesianposition and velocity coordinates. In particular, the only unit vector that plays a role is thephysical b = B / | B | .4) Gyroaverages and Fourier expansions in gyrophase are implemented in these coordinatesusing a coordinate-independent formulation of these operations.5) The approach manages to be manifestly Hamiltonian while addressing the order-mixingissue in a computationally attractive manner; for each m > n , the n ’th-order contributionto the perturbative coordinate transformation is determined without knowledge of any of3he details of the m ’th order contribution.6) The manner in which we address the order-mixing issue obviates the high degree offreedom in the form of the transformed Lagrangian.In what follows, we will describe our algorithm and report on the equations of motiongenerated in the two representations just alluded to. We will not evaluate these new repre-sentations in terms of their simplicity or time-validity properties; a properly thorough studyof these properties will appear in future work. We will begin with four sections describingwhat our algorithm is meant to do as well as our motivation for selecting an algorithm withthe novel features just described. In section II, we give a schematic overview of HamiltonianLie transform-based perturbation theory in order to remind the reader of the goal of theguiding center expansion. We then describe the motivation for selecting our algorithm viaa description of three difficulties we faced while developing it, and how we overcame them.In particular, sections III, IV, and V are devoted to discussing the difficulties presented bythe order-mixing issue; the desire for manifestly gyrogauge invariant results; and the taskof computing gyroaverages and gyroharmonics, respectively. With all of the motivationsin place, we present our algorithm in section VI. Finally, in section VII, we present theresults of automatically performing the guiding center expansion with our algorithm in twopreviously unstudied representations. II. A SCHEMATIC FOR HAMILTONIAN LIE TRANSFORMPERTURBATION THEORY
In this section, we will review the general structure and purpose of the guiding centerexpansion, and thereby indicate precisely what our algorithm is meant to do. We thendiscuss three key difficulties we faced while trying to develop the algorithm before actuallypresenting presenting it. The purpose of these first four sections is to provide a narrativeexplaining why the algorithm looks the way it does. Readers only interested in the algorithmitself can skip straight to section VI, but it may still be useful to skim these early sectionsin order to become familiar with our notation.We begin by recalling the coordinate-independent formulation of Hamiltonian dynamicalsystems . This formulation makes use of Cartan’s exterior calculus of differential forms; avery brief overview of the latter is provided in Appendix A. The phase space M is assumed4o be an even dimensional smooth manifold equipped with a symplectic two-form ω . Thedynamical equations are then specified by a function H : M → R known as the Hamiltonianfunction via Hamilton’s equations i X H ω = d H, (1)where X H is the vector field that specifies the time derivative of any particle’s phase spacelocation c ( t ) ∈ M , i.e. c ′ ( t ) = X H ( c ( t )). In any local coordinate system ( z i ) on M ,Hamilton’s equations become ˙ z i ω ij = ∂H∂z j , (2)where ˙ z i are the components of the vector field X H = ˙ z i ∂∂z i and ω ij = ω ( ∂∂z i , ∂∂z j ).In the guiding center problem, the phase space is the six-dimensional position-velocityspace, M = R × R , equipped with the symplectic form ω ǫ = − d ϑ ǫ , where the one-form ϑ ǫ is given in terms of the the magnetic vector potential A and the guiding center orderingparameter ǫ by ϑ ǫ = A · dx + ǫv · dx. (3)The equations of motion are then specified by the Hamiltonian function H = ǫ v · v . Ascan be readily verified, the vector field X H ( ǫ ) in the natural cartesian coordinates on M isgiven by ˙ v ( x, v ) = v × B ( x ) (4)˙ x ( x, v ) = ǫv. Strictly speaking, the placement of the ordering parameter ǫ = ρ/L in ϑ ǫ , and therefore itsplacement in the equations of motion, is only justified in appropriate dimensionless variables,as discussed in Ref. 3. However, we can regard Eq. (4) as a dimensional equation if we thinkof ǫ as a formal ordering parameter and if we normalize A by a particle’s charge-to-massratio so that B has the dimension of frequency.When ǫ = 0, which corresponds to the asymptotic limit where a particle undergoes gy-romotion with zero gyroradius and vanishingly slow drift, the equations of motion given inEq. (4) are gyrosymmetric (see section IV for the precise definition of gyrosymmetric ten-sors). Because the particle trajectories are periodic in this limit, Kruskal’s general theory M using a non-unique ǫ -dependent near-identity transformation T ǫ : M → M such that the transformed X H ( ǫ ), X H ( ǫ ′ ) ≡ T ǫ ∗ X H ( ǫ ), is gyrosymmetric to all orders in ǫ .The goal of the guiding center theory, and therefore the algorithm we will present later,is to find such a T ǫ . Because performing this task requires a degree of ingenuity, a numberof useful methods have been developed. Of particular relevance to the present work arethose methods that employ Lie transforms. In these cases, one posits that the desiredtransformation from the old phase space to the new, deformed phase space can be expressedin the form T ǫ = ... ◦ exp( − G n ( ǫ )) ◦ ... ◦ exp( − G ( ǫ )) , (5)where, for each n , G n ( ǫ ) : M → T M is a vector field that tends to zero as ǫ →
0, and does somore rapidly than does G m ( ǫ ) with m < n . The requirement that the transformed equationsof motion be gyrosymmetric then reduces to a sequence of requirements on the G n ( ǫ ). Thus,the Lie transform approach to finding T ǫ reduces to finding a sequence of G n ( ǫ ) that satisfythe latter requirements.One can derive these requirements in one of two ways. The direct method, which recentlymade an appearance in Ref. 21 (also see Ref. 22), consists of formally computing thetransformed X H ( ǫ ), X ′ H ( ǫ ) = T ǫ ∗ X H ( ǫ ), using Eq. (5) and then demanding that the resultbe gyrosymmetric to all orders. The Hamiltonian method, due to Littlejohn , consistsof calculating the transformed d ϑ ǫ and H , T ǫ ∗ d ϑ ǫ and T ǫ ∗ H , and then demanding thateach of be gyrosymmetric to all orders. These two methods are related by the general factthat if the symplectic form, ω , and Hamiltonian, H , appearing in Hamilton’s equations(1) admit a symmetry, then so does X H . Each method also involves making a number ofarbitrary decisions to completely determine the G n ( ǫ ); different choices lead to differentrepresentations.In principle, either approach can be automated on a computer. Indeed, historically thishas been one of the advertised “features” of the Lie transform approach to perturbationproblems in general. However, the Hamiltonian approach has the advantage of providingan attractive means for truncating the results of the expansion. Generally, given a near-identity transformation and a Hamiltonian system specified by a one-form and Hamiltonian,there are two ways to develop “finite” approximations, or truncations to the transformed6quations of motion. One approach is to directly truncate the transformed equations ofmotion at some order in the expansion parameter. The other approach is to truncate thetransformed one-form and Hamiltonian and use the vector field specified by the ensuingHamilton’s equations to approximate the transformed equations of motion. Either approachmay be used to generate arbitrarily accurate approximations to the transformed equationsof motion. However, the second approach always produces an approximation to the trans-formed equations of motion that is rigorously Hamiltonian. Because the original dynamicalsystem is Hamiltonian, the second “Hamiltonian” truncation scheme is theoretically prefer-able. Thus, an advantage offered by the Hamiltonian approach to deriving the guiding centertransformation is as follows. Because the Hamiltonian approach directly tracks the one-formand Hamiltonian through the near-identity transformation (5), the Hamiltonian truncationscheme is always immediately available; it is not necessary to compute the transformed one-form and Hamiltonian after finding the G n ( ǫ ) as it would be if the G n ( ǫ ) were calculatedusing the direct method. It is for this reason that we pursue the Hamiltonian approach inthe present work. III. DIFFICULTY 1: ORDER-MIXING
While developing the algorithm for automating the Hamiltonian Lie transform approachto finding T ǫ , we encountered three key difficulties. In this section and the two that followwe will describe each in turn, as well as the manner in which we overcame each difficulty.The first issue is rooted in the special form of ϑ ǫ given above. If one specifies the ǫ -dependence of the G n ( ǫ ) according to G n ( ǫ ) = ǫ n g n , then the one-form ϑ ǫ on the deformedphase space is given by ϑ ′ ǫ = A · dx + ǫv · dx + ǫL g ( A · dx ) + O ( ǫ ) , (6)where L g denotes the Lie derivative with respect to the vector field g (see appendix A).In order to find one of the transformations guaranteed by Kruskal’s theory, the combination v · dx + L g ( A · dx ) must be gyrosymmetric, modulo closed one-forms . If we write g = g x · ∂∂x + g v · ∂∂v , this condition can be satisfied by choosing g x = | B | v × b + αb , g v = Y ,where α and Y are arbitrary. However, as would become clear upon analyzing higher-ordercontributions to the transformed one-form, there are in fact constraints on α and Y , meaning7t least part of the Freedom in specifying g suggested by the first order change in the one-form is only apparent. This is a special case of the more general order-mixing issue; theconstraints on a given g n can only be deduced by considering multiple orders ϑ ′ ǫ , in particularorders whose form is effected by g m for m > n .While order-mixing does not prevent the success of the Hamiltonian method (see Ref.4 for one way of coping with it), from a computational point of view, it is bothersome. Itobfuscates the extent to which the choices one needs to make to find the various g n arecoupled across n -values. If the coupling were severe enough, then any algorithm one mightconstruct to automate these choices could be very complicated.In order to overcome this difficulty, we designed our algorithm to satisfy Resolution of Difficulty 1:
The rule for determining G n ( ǫ ) does not rely on any specificknowledge of any component of G m ( ǫ ) whenever m > n .In section VI, the precise manner in which the algorithm accomplishes this will becomeclear. IV. DIFFICULTY 2: MANIFEST GYROGAUGE INVARIANCE
In order to understand the second difficulty we faced in developing a good algorithm forautomating the guiding center calculation, one needs to understand the usual definition ofa gyrosymmetric tensor. This definition refers to a special type of coordinate system on M ,any instance of which we will call a fibered coordinate system . A fibered coordinate systemon M consists of an open subset U ⊂ M , 5 smooth functions ξ i : U → R , i = 1 , ...,
5, andone additional function θ : U → R mod 2 π satisfying:F1 : The six functions ξ i , i = 1 , ...,
5, and θ define a valid coordinate system on U F2: Holding the ξ i fixed, θ parameterizes, in a left-handed sense relative to b , the zero’thorder ( ǫ = 0) solutions to Eq. (4), which are called loops by Kruskal.The standard example of a family of fibered coordinate systems used in guiding center8 IG. 1. A typical arrangement of the perpendicular unit vectors e , e for a uniform magnetic fieldthat points out of the page. The two sets of arrows represent e and e . While in this case, e and e are not required to vary in space, for a more general sort of magnetic field, they would be.Reprinted from Phys. Plasmas 19, 052106 (2012). Copyright 2012 American Institute of Physics. theory is constructed as follows. First find a smooth unit vector field e perpendicular tothe magnetic field, e · b = 0. e ( x ) and e ( x ) = ( b × e )( x ) span the plane perpendicularto b ( x ) for each x in the domain of definition, D ⊂ R , of e , as depicted in Figure 1.As shown in Ref. 18, D cannot always be taken to be the entire 3-dimensional domainparticles move through. A fibered coordinate system can then be defined on the open subsetof phase space U ≡ { ( x, v ) ∈ M | x ∈ D and b ( x ) × v = 0 } . Labeling the ξ i according to( ξ , ξ , ξ ) = x, ξ = v ⊥ , ξ = v k , these functions are defined by the relation x = x ( x, v ) (7) v = v k ( x, v ) b ( x )+ v ⊥ ( x, v ) (cos( θ ( x, v )) e ( x ) − sin( θ ( x, v )) e ( x )) . A gyrosymmetric tensor can then be defined in terms of fibered coordinate systems asfollows. A tensor is gyrosymmetric if its components in an arbitrary fibered coordinatesystem do not depend on θ . Note that one doesn’t have to look at a tensor in every fibered9oordinate system to check this property; it is enough to check in a collection of fiberedcoordinate systems that cover M .This standard definition of gyrosymmetric tensors motivates the standard approach toderiving the constraints on the G n ( ǫ ). One first writes out the components of the transformed d ϑ ǫ and H in a family of fibered coordinate systems on M that cover M . Then one choosesthe local representatives of G n ( ǫ ) to eliminate the θ -dependence in these components in eachcoordinate system in the covering. For consistency , one also must demand that the localdefinitions of G n ( ǫ ) agree when changing from one fibered coordinate system in the coveringto another. This last consistency condition is one statement of the principle of gyrogaugeinvariance. Satisfying this consistency condition is equivalent to demanding that the localrepresentatives of G n ( ǫ ) have the gyrogauge invariant form first identified by Littlejohn inRef. 5.There is nothing conceptually wrong with this approach to finding the G n ( ǫ ), and it canbe made to work. However, there is a very practical problem with proceeding in preciselythis manner on a computer. In order to verify that a given expression for G n ( ǫ ) in a fiberedcoordinate system satisfies the principle of gyrogauge invariance, it is often necessary toaccount for non-trivial vector identities involving the perpendicular unit vectors e and b .For instance, the identity ∇ × (( ∇ e ) · e ) = 12 b (cid:18) Tr( ∇ b · ∇ b ) − ( ∇ · b ) (cid:19) (8)+ ( ∇ · b ) b · ∇ b − b · ∇ b · ∇ b, and even more complicated identities generated by taking derivatives of Eq. (8) must berecognized. Presently, there is no general method that would allow one to do this on acomputer in all cases one might encounter. Thus, one cannot guarantee that the G n ( ǫ )produced by a computer following the above procedure will manifestly exhibit gyrogaugeinvariance, i.e. it will not be obvious that G n ( ǫ ) is gyrogauge invariant, even if it actuallyis. In order to avoid this issue, we have chosen to avoid using fibered coordinate systemsaltogether. Resolution of Difficulty 2:
All tensors are expressed and manipulated in cartesian posi-tion and velocity coordinates. 10y proceeding in this manner, the results generated by our algorithm ( G n ( ǫ ), for ex-ample) will be expressed entirely in terms of v , | B | , b , and derivatives thereof, therebymaking our algorithm manifestly gyrogauge invariant; perpendicular unit vectors e , e andthe gyrophase coordinate θ are never even introduced. V. DIFFICULTY 3: FOURIER ANALYSIS WITHOUT INTRODUCINGADDITIONAL UNIT VECTORS
The third issue we wish to discuss arises as a result of our resolution of the difficultydiscussed in the previous section. Because part of the motivation for choosing to work inCartesian coordinates was to avoid introducing additional unit vectors, it would be a stepbackward if we had to introduce additional unit vectors in order to perform Fourier analysisin the gyrophase. Is there a method for computing the gyroaverages and gyroharmonics ofa tensor in Cartesian coordinates without introducing more unit vectors than the necessary b ? A conceptually appealing way to answer this question is to first derive some coordinate-independent properties of the gyrosymmetry that would allow one to answer this questionin an arbitrary coordinate system, and then specialize to cartesian coordinates. To ourknowledge, this interesting mathematical exercise has not been discussed elsewhere in theliterature, and so we will provide the details in the remainder of this section.First notice that in a fibered coordinate system ( ξ i , θ ) (this is shorthand for the sextuplet( ξ , ..., ξ , θ )), a function f : U → R is gyrosymmetric if and only if f ( ξ i , θ + ψ ) = f ( ξ i , θ ) (9)for all constants ψ . If, for each ψ ∈ R mod 2 π , we define the mapping Φ ψ : U → U usingthe formula Φ ψ ( ξ i , θ ) = ( ξ i , θ + ψ ) , (10)then the condition given in Eq. (9) can be re-expressed asΦ ∗ ψ f = f (11)11or each ψ ∈ R mod 2 π . Here Φ ∗ ψ denotes the pullback operator on functions, Φ ∗ ψ f = f ◦ Φ ψ .While the formula (9) only makes literal sense in a fibered coordinate system, the family ofmappings Φ ψ can actually be given a coordinate independent definition. Indeed, in cartesianposition and velocity coordinates we have Φ ψ ( x, v ) = (12)( x, v · b ( x ) b ( x ) + cos( ψ ) b ( x ) × ( v × b ( x )) + sin( ψ ) v × b ( x )) . Thus, gyrosymmetric functions f : M → R can be alternately characterized as those func-tions that satisfy the analogue of Eq. (11), Φ ∗ ψ f = f for each ψ ∈ R mod 2 π .What about more general tensor fields? Because the pullback operator of a mapping M → M is well defined on the entire tensor algebra, it is tempting to postulate that atensor field τ is gyrosymmetric if and only if Φ ∗ ψ τ = τ for all ψ ∈ R mod 2 π . This is indeedcorrect; it is a straightforward exercise to verify that this characterization is equivalent tothe usual one stated in the previous section.What is going on here? If we fix a ψ ∈ R mod 2 π , then the mapping Φ ψ : M → M canbe regarded as a global rearrangement, or relabeling, of points in M . If we regard Φ ψ aspointing from the “new arrangement” to the “old arrangement”, then Φ ∗ ψ τ is nothing morethan τ , regarded as a tensor in the old arrangement of M , expressed in the new arrangment.Thus, from this point of view, we see that gyrosymmetric tensors are precisely those tensorswhose form is invariant under any of the rearrangments in the family Φ ψ .With this coordinate-independent characterization of gyrosymmetric tensors in hand,we now seek a corresponding coordinate-independent version of Fourier analysis in the gy-rophase θ . The catch is that we do not desire to work with the gyrophase coordinate θ directly as the latter is only defined in fibered coordinate systems. Instead we will use theparameter ψ in the family of maps Φ ψ as a surrogate of sorts.Given an arbitrary tensor τ , set τ ψ = Φ ∗ ψ τ . τ ψ can be regarded as a periodic tensorfield-valued function of the single variable ψ with period 2 π . Therefore it admits a Fourierexpansion τ ψ = h τ i + ∞ X k =1 (Π k τ ) cos( ψ ) + ( ¯Π k τ ) sin( ψ ) , (13)12here the tensor fields h τ i , Π k τ , and ¯Π k f are given by h τ i = 12 π Z π (Φ ∗ ψ τ ) dψ (14)Π k τ = 1 π Z π (Φ ∗ ψ τ ) cos( kψ ) dψ ¯Π k τ = 1 π Z π (Φ ∗ ψ τ ) sin( kψ ) dψ. Note that Π k τ and ¯Π k τ are not gyrosymmetric tensors. Instead they satisfy the identitiesΦ ∗ ψ (Π k τ ) = cos( kψ )(Π k τ ) + sin( kψ )( ¯Π k τ ) (15)Φ ∗ ψ ( ¯Π k τ ) = − sin( kψ )(Π k τ ) + cos( kψ )( ¯Π k τ ) . However, as the notation suggests, h τ i is indeed gyrosymmetric.The validity of these formulae only relies on the fact that the mapping Φ ψ satisfies thegroup property Φ ψ + ψ = Φ ψ ◦ Φ ψ . Therefore, they may be applied to tensors on anymanifold equipped with such a family of mappings. These formulae also represent a shiftin perspective from the Fourier analysis one would usually employ in fibered coordinatesystems. To see this, consider a fibered coordinate system and the Fourier expansion of afunction f ( ξ i , θ ) = f ( ξ i ) + P k a k ( ξ i ) cos( kθ ) + b k ( ξ i ) sin( kθ ). It is not true that Π k f = a k ;instead (Π k f )( ξ i , θ ) = a k ( ξ i ) cos( kθ ) + b k ( ξ i ) sin( kθ ). Thus, the operators Π k and ¯Π k arenot merely calculating the usual Fourier coefficients a k ( ξ i ) , b k ( ξ i ). Moreover, Π k f is a gen-uine scalar whereas the usual Fourier coefficients a k ( ξ i ) , b k ( ξ i ) have non-trivial transforma-tion laws when passing from one fibered coordinate system to another (i.e. a change ofgyrogauge). Indeed, if ( ξ i , θ ) and ( ξ i , θ ′ ) are two fibered coordinate systems related by θ ′ = θ + φ ( ξ i ), then the usual Fourier coefficients in the primed coordinates a ′ k ( ξ i ) , b ′ k ( ξ i ) arerelated to the usual Fourier coefficients in the unprimed coordinates by a k ( ξ i ) = a ′ k ( ξ i ) cos( kφ ( ξ i )) + b ′ k ( ξ i ) sin( kφ ( ξ i )) (16) b k ( ξ i ) = b ′ k ( ξ i ) cos( kφ ( ξ i )) − a ′ k ( ξ i ) sin( kφ ( ξ i )) . (17)Therefore, the coordinate-independent Fourier analysis given by equations (13) and (14)calculate gyrogauge invariant combinations of the usual Fourier coefficients.The Fourier inversion formula, Eq. (14), together with the invariance properties given inEq. (15), is sufficient to solve all of the linear partial differential equations that one encounters13hile deriving expressions for the G n ( ǫ ) in any coordinate system. This is because: (a)all tensors encountered while deriving the guiding center expansion contain finitely manygyroharmonics; (b) the differential operator L ξ , where ξ = v × b · ∂∂v becomes an algebraicoperator on gyroharmonics, L ξ h τ i = 0 (18) L ξ (Π k τ ) = k ¯Π k τ (19) L ξ ( ¯Π k ) = − k Π k τ ; (20)and (c) L ξ is the only partial differential operator that ever needs to be inverted. Thus,we have effectively solved the problem of performing “Fourier analysis in θ ” in cartesianposition and velocity coordinates without ever referring to fibered coordinate systems. Wehave incorporated this solution into our algorithm as Resolution of Difficulty 3:
Gyroaverages and gyroharmonics are calculated in cartesianposition and velocity coordinates using Eqs. (14) and (15).
VI. THE ALGORITHM
As discussed in section II, the goal of the algorithm is to find a transformation T ǫ in theform given in Eq. 5 such that T ǫ ∗ d ϑ ǫ and T ǫ ∗ H are each gyrosymmetric to all orders in ǫ (section IV gives the general definition of a gyrosymmetric tensor). This T ǫ consists of aconcatenated sequence of transformations of the form exp( Y ). Thus, we are free to think of T ǫ as the result of many intermediate transformations, each closer to the identity transfor-mation than the last. Our algorithm proceeds by finding expressions for these intermediatetransformations (which amounts to specifying a G n ( ǫ )), one at a time, according to thefollowing recursive procedure.Suppose that some finite number of intermediate transformations have been performed.Let Θ ǫ and H ǫ denote the resulting one-form and Hamiltonian following this partial rear-14angement of M , and assume they have the form:Θ ǫ = ϑ + ǫϑ + ... + ǫ N ϑ N + ∞ X k =1 ǫ N + k α k (21) H ǫ = H + ... + ǫ N − H N − + ∞ X k =1 ǫ N − k h k , where N >
1, the ϑ j and H j are all gyrosymmetric, and the α j and h j are not necessarilyso. Suppose further that Ξ ǫ = ϑ + ... + ǫ N ϑ N satisfies the three propertiesND1 : d Ξ ǫ is a non-degenerate two-form.ND2: If β is an ǫ -independent one-form, then the vector field Y ( ǫ ) defined by i Y ( ǫ ) d Ξ ǫ = β (i.e. Y is the application of the Poisson tensor defined by Ξ ǫ to β ) is O ( ǫ − ).ND3: When β = − d H , the leading order behavior of Y ( ǫ ) is given by | B | ǫ ξ ≡ | B | ǫ v × b · ∂∂v .In this setting, which will serve as our inductive assumption, it is possible to find a trans-formation exp( − G ( ǫ )), for some small vector field G ( ǫ ), such that after this transformation,the one-form and the Hamiltonian have the same form as in Eq. (21), but with N replacedwith N + 1, i.e. the one-form and Hamiltonian are each gyrosymmetric to one higher orderthan previously. This also means that the transformed Ξ ǫ will automatically continue tosatisfy properties ND1-3. We will call a G ( ǫ ) that produces a transformation exp( G ( ǫ )) withthe latter two properties a recursive vector field.To see that one can in fact find many recursive vector fields under the inductive assump-tion, let G ( ǫ ) be a vector field that solves the algebraic equation (see appendix C for asolution method) i G ( ǫ ) d Ξ ǫ + ǫ N +1 α + ǫ N +1 d S = i h G ( ǫ ) i d Ξ ǫ + ǫ N +1 h α i , (22)where S is the unique function with h S i = 0 (see section V for the definition of the generaltensor gyroaverage operator hi ) that solves the partial differential equation (see appendix Dfor a solution method) h − | B | i ξ α − | B | i ξ d S = h h i − | B | i ξ h α i . (23)15ote that the oscillatory part of G ( ǫ ), ˜ G ( ǫ ) = G ( ǫ ) − h G ( ǫ ) i , is then uniquely determined,but the gyroaverage h G ( ǫ ) i is left completely free. Constrain the latter so that it satisfiesi h G ( ǫ ) i d Ξ ǫ = ǫ N +1 γ, (24)where γ is any ǫ -independent gyrosymmetric one-form. Note that h S i = 0 is not an arbitrarychoice, but is necessary in order for Eq. (22) to be self consistent. Indeed, upon gyroaveragingEq. (22), we see that d h S i = 0, which implies h S i is constant. This constant is clearlyinconsequential, and so we set it to zero.As is readily verified, such a G ( ǫ ) satisfies the following important properties.P1: G ( ǫ ) = O ( ǫ N − ) , but will be a rational function of ǫ .P2: Upon applying the transformation exp( − G ( ǫ )), the transformed one-form (moduloclosed one-forms) and Hamiltonian, Θ ′ ǫ and H ′ ǫ , becomeΘ ′ ǫ = ϑ + ... + ǫ N ϑ N (25)+ ǫ N +1 ( h α i + γ ) + O ( ǫ N +2 ) , and H ′ ǫ = H + ... + ǫ N − H N − (26)+ ǫ N − ( h h i + | B | i ξ γ ) + O ( ǫ N ) . Thus, the entire family of G ( ǫ ) just defined, a family which may be regarded as beingparameterized by the arbitrary gyrosymmetric one-form γ , consists of recursive vector fields.We refer the reader to Appendix E to clearly see the motivation for choosing Eqs. (22) and(23).With these recursive vector fields in hand, all that we must now show is that there is somebase case, consisting of a one-form and Hamiltonian in the form specified by Eq. (21), fromwhich our recursive algorithm can start. Unfortunately this base case clearly cannot be thenatural choice, Θ ǫ = A · dx + ǫv · dx and H ǫ = ǫ v · v , as this pair is not in the form specifiedby Eq. (21). However, this issue is easy to resolve. First, notice that if the transformed ǫ H is gyrosymmetric, then H will be too. Therefore, remove the ǫ from the Hamiltonian.Second, recognize that we are free to perform a preparatory transformation before finding16he G n ( ǫ ). In particular, we can apply a preparatory near-identity transformation of theform exp( − ǫG ) that transforms H ǫ = v · v and Θ ǫ = A · dx + ǫv · dx into the form specifiedby Eq. (21) with N = 2. For instance, with G = − v × b | B | · ∂∂x (27)+ (cid:18) ( v · b ) ∇ b · ( b × v ) | B | + ( v × b ) · ∇ b · vb | B | + ( v · b )( b · ∇ × b ) b × ( v × b ) | B | (cid:19) · ∂∂v then we arrive at the satisfactory starting pointΘ ǫ = A · dx + ǫ ( v · b ) b · dx (28)+ 12 | B | ǫ (cid:18) v × b · dv − ( v · b )[ ∇ b · ( v × b )] · dx (cid:19) + O ( ǫ ) H ǫ = 12 v · v + O ( ǫ ) , which can be verified by directly calculating Lie derivatives (using VEST , for instance).More generally, G must be chosen so that, after applying the preparatory transformationgenerated by G , the one-form has the form Θ ǫ = θ + ǫθ + ǫ θ + O ( ǫ ) for gyrosymmetricone-forms θ , θ , θ , and d Θ ǫ is non-degenerate. After applying such a preparatory transfor-mation, properties ND1-3 will automatically be satisfied, and the inductive procedure canbegin. We found the G just given, as well as the G given in the second example below bydirectly analyzing the transformed one-form and demanding that it satisfy these properties.On the other hand, it would be interesting to find the most general G that accomplishesthe preparatory transformation. We leave this to future studies.To summarize, our algorithm for finding the G n ( ǫ ) that generate T ǫ proceeds as follows. Because any real calculation can only calculate a finite number of the G n ( ǫ ), whenone stops calculating additional G n ( ǫ ), the one-form and Hamiltonian will be of the formspecified in Eq. (21) with N = N max . Therefore, specify the desired N max . Define three integers N , M , and l with initializations N = 2, M = N max −
2, and l = 1. 17 : Apply a preparatory transformation, such as that given in Eq. (27), so that the one-formand Hamiltonian have the form specified in Eq. (21). Only the first M non-gyrosymmetricterms must be calculated in each case. Find a recursive vector field G l ( ǫ ). Use Eq. (22) to find the unique oscillatory part˜ G l ( ǫ ), and specify the gyroaveraged part h G l ( ǫ ) i using an arbitrary gyrosymmetric one-form γ l according to Eq. (24). Store G l ( ǫ ). Set l = l + 1, N = N + 1, and M = M − Using the recursive vector field just derived to specify the transformation, expressthe new one-form and Hamiltonian in the form specified in Eq. (21). Only the first M non-gyrosymmetric terms must be calculated in each case. This amounts to applying exp(i G l ( ǫ ) d )and exp( L G l ( ǫ ) ) to the old one-form and Hamiltonian to generate the new α i and h i . If N = N max , stop. Else, return to step .Note that different representations of the guiding center expansion will be generated foreach choice of the sequence of gyrosymmetric one-forms γ l and the preparatory transforma-tion generated by G . In particular, if one does not attempt to constrain the form of thetransformed Hamiltonian, the O ( ǫ ) contribution to the transformed one-form can be anygyrosymmetric one-form whatsoever, including 0. Likewise, if one does not attempt to con-strain the form of the transformed one-form, then the O ( ǫ ) contribution to the transformedHamiltonian can be specified arbitrarily (at least away from those points in phase spacewhere v × b = 0 as i ξ γ = 0 at those points).Also note that the presence of a preparatory transformation implies that the completetransformation from the old phase space to the new, deformed phase space is given by T ǫ ◦ exp( − ǫG ), with T ǫ = ... ◦ exp( − G ( ǫ )) ◦ exp( − G ( ǫ )). In particular this motivatesour convention for indexing the G n ( ǫ ); by specifying the preparatory transformation asexp( − ǫG ), we obtain the appealing result that G n ( ǫ ) = O ( ǫ n ) in spite of the fact thatboth the preparatory transformation and the first transformation generated by the recursiveprocedure are O ( ǫ ). 18inally, note that now we can identity the precise manner in which this algorithm cir-cumvents the order-mixing issue. The solution is the combination of the preparatory trans-formation and the fact the the recursive procedure determines each non-preparatory G n ( ǫ )without any knowledge of the form of G m ( ǫ ) whenever m > n . The preparatory transforma-tion is necessary to ensure that the equation defining G ( ǫ ), Eq. (22), in the first recursivestep can be solved. In particular, it guarantees that d ( ϑ + ǫϑ + ǫ ϑ ) is invertible. Afterthis initial step is complete, the two-form inversion required to compute the higher G n ( ǫ )is always possible because the two form in question is always a small perturbation to atwo-form that is known to be invertible. Moreover, none of the quantities that appear in theequation defining G n ( ǫ ), Eq. (22), depend on knowledge of any of the G m ( ǫ ) with m > n .The key observation that lead to this approach to determining the G n ( ǫ ) was that it is notnecessary to assume the ǫ -dependence of each G n ( ǫ ) a priori . Indeed, the G n ( ǫ ) determinedby our scheme are ratios of polynomials in ǫ , whereas the usual approach assumes each G n ( ǫ )is a monomial ǫ . VII. TWO NEW REPRESENTATIONS
To illustrate the use of our algorithm, to suggest the relative ease of computing higher-order corrections to the guiding center expansion on a computer, and to emphasize thefact that there are still unexplored representations of the guiding center dynamics, we nowturn to presenting the results of using the algorithm to derive two previously unexploredrepresentations of the guiding center dynamics. Each representation we present here willchoose γ l + h α i = 0 in step , meaning each representation is closely related to the so-calledHamiltonian representation discussed in Ref. 4; this choice of the γ l leads to a transformedone-form that truncates at second order in ǫ . While this property does not completelycharacterize the Hamiltonian representation of Ref. 4, it is a characteristic thereof. The tworepresentations will differ from eachother in which preparatory transformation is used instep of the algorithm. Thus, these examples give a taste of the consequences of choosingdifferent preparatory transformations, but not of different schemes for choosing γ l . We willnot evaluate either representation in terms of simplicity or time-validity, but instead leavethis study to future work.Before we present these representations, it is important that we point out a subtle aspect19f our approach that follows from the fact that we work in cartesian position and velocitycoordinates. By employing these coordinates, all of our calculations are done in the full six-dimensional single-particle phase space. Thus, the transformed one-form and Hamiltonian,Θ and H , specified by our algorithm are quantities defined on a six-dimensional space. Thisimplies that the new dynamical vector field (which specifies the equations of motion) specifiesthe evolution of the six variables ( X, Y, Z, v x , v y , v z ) (note that x=(X,Y,Z) in our notation).In particular, µ , the magnetic moment is not a coordinate as it is in many treatmentsof guiding center theory. It is important to understand that in spite of this fact, thesetransformed six-dimensional equations nevertheless possess an exact conservation law for anytruncation of Θ and H . This follows from the fact that the full six-dimensional transformedequations of motion are Hamiltonian, meaning they satisfy Hamilton’s equations on thesix-dimensional phase space i X d Θ = − d H , where X = ( ˙ x, ˙ v ). Thus, the Hamiltonianversion of the Noether theorem implies that the function µ = i ξ Θ ( ξ = v × b · ∂∂v is theinfinitesimal generator of the gyrosymmetry) is constant along trajectories of the vectorfield X as a result of the invariance of Θ and H under the continuous family of symmetriesdefined by the gyrosymmetry Φ ψ . It is also important to realize that it is simple to obtainan expression of results produced by our algorithm in a coordinate system that uses theconserved quantity as a coordinate by using the six functions ( X, Y, Z, v k , µ, θ ), with θ agyrophase function defined relative to some perpendicular unit vectors, as coordinates. Inthis type of coordinate system, it will also be obvious how the one-form and Hamiltoniandescend to the reduced (four-dimensional) phase space parameterized by ( X, Y, Z, v k ) forfixed µ = i ξ Θ. To illustrate this fact, we will present the one-form in each representationboth in cartesian coordinates, as it is given directly by Mathematica, and in (
X, Y, Z, v k , µ, θ )coordinates, which we calculate by hand starting from the cartesian result. We will not dothe same for the Hamiltonian as it will be trivial to express the Cartesian Hamiltonian interms of ( X, Y, Z, v k ).For each representation, we will provide explicit expressions for the transformed one-formaccurate to all orders in ǫ . For the transformed Hamiltonians, we will provide H , H , and H . This information is enough to accurately express the full transformed equations ofmotion up to and including terms of order ǫ . Moreover, the reduced equations of motion,which describe the evolution of the guiding center position x and parallel velocity v k , canbe specified up to and including terms of order ǫ . For the sake of brevity, we will not20isplay the G n ( ǫ ). However, we stress that, when equipped with a copy of our code, findingthese vector fields that specify the transformation would be a simple task for any interestedreader. In fact, each of the results below takes about fifteen minutes to derive on a laptopcomputer equipped with VEST . We would also like to stress that all of the following resultshave been checked thoroughly. Using
VEST , we have explicitly expanded the pushforwardoperator T ǫ ∗ into a series of Lie derivatives along the G n ( ǫ ) and applied it to the zero’thorder Hamiltonian and one-form. The resulting expressions agree exactly with the resultsreported below up to the relevant order. They have also been checked indirectly by verifyingthat we reproduce the well-known first-order correction to the magnetic moment adiabaticinvariant (in untransformed variables), µ (see Ref. 28, for example): µ = 1 | B | (cid:18) v · ∇ b · ( v × b )( v · b ) (29) −
34 ( v × b ) · ∇ b · v ( v · b ) − b × κ · v ( v · b ) + 12 | B | ( v × b ) · ( v × b ) ∇| B | × b · v (cid:19) . Moreover, we have investigated µ using VEST . We have verified that the µ predicted bythe G n ( ǫ ) in each representation agree with one another. We have also directly verifiedthat the time derivative of µ + ǫµ + ǫ µ along the Lorentz force equations of motion˙ v = | B | v × b, ˙ x = ǫv is O ( ǫ ) (in general, the adiabatic invariant series must satisfy ddt ( P Nk =0 ǫ k µ k ) = O ( ǫ N +1 )). Finally, we have compared our expression for µ with thatgiven in Ref. 28, which is the only published µ applicable to general magnetic geometrywe are aware of. Interestingly, our expression is numerically different from the result re-ported in Ref. 28. However, this difference is most likely due to a likely typographical errorin Ref. 28, as explained in Appendix F. We also present our expression for µ in Appendix F. Example 1:
This representation will be defined by the use of the preparatory transformation alreadygiven in Eq. (27) and by always choosing γ l + h α i = 0. The consequences of these choicescome in the form of a transformed one-form equal to that given in Eq. (28) to all orders, thussimplifying the form of the transformed Poisson bracket. In fact, by expressing Eq. (28)in the usual sort of fibered coordinate system as well as cartesian position and velocity21oordinates, we see Θ ǫ = A · dx + ǫ ( v · b ) b · dx (30)+ 12 | B | ǫ (cid:18) v × b · dv − ( v · b )[ ∇ b · ( v × b )] · dx (cid:19) = A · dx + ǫv k b · dx + ǫ µ [ dθ − R · dx ]where µ = i ξ Θ = v ⊥ | B | and R = ( ∇ e ) · e , meaning the transformed Poisson bracket isexactly the same as that given in Ref. 5. Meanwhile the transformed Hamiltonian is givenby H = 12 v · v (31) H = 12 ( v · b ) µτH = µ (cid:18) ∇ · b ) + 116 κ · κ + 14 b · ∇ ( ∇ · b ) −
116 tr[ ∇ b · ∇ b + ∇ b · ( ∇ b ) T ] − ∇ ln | B | · ∇ ln | B | + 14 κ · ∇ ln | B | + 14 | B | ∇ | B | (cid:19) + µ ( v · b ) | B | (cid:18)
18 tr[3 ∇ b · ∇ b − ∇ b · ( ∇ b ) T ] + 18 ( ∇ · b ) + 12 b · ∇ ( ∇ · b ) + 138 κ · κ − κ · ∇ ln | B | (cid:19) − ( v · b ) | B | (cid:18) κ · κ (cid:19) where µ = i ξ Θ = ( v × b ) · ( v × b )2 | B | , τ = b · ∇ × b , and κ = b · ∇ b . Note that the H in this repre-sentation differs from the H in Brizard and Tronko’s Hamiltonian representation (also seeRef. 6), although it is similarly complicated. This difference is not just apparent; we haverigorously compared our H with Brizard and Tronko’s using VEST and found they arenot numerically equal. In order to recover Brizard and Tronko’s representation, it would benecessary to either: (a) choose h α i + γ l = d f l with appropriately chosen f l (we have chosen f l = 0 in this example); (b) alter the preparatory transformation Eq. (27); or (c) do both(a) and (b). Example 2: γ l + h α i = 0, but the preparatorytransformation exp( − ǫG ) will be specified by G = − v × b | B | · ∂∂x + 1 | B | (cid:18) ( v · b )( v × b ) · ∇ b (32) − v · b ) ∇ b · ( v × b ) + 14 v · [ ∇ b + ∇ b T ] · ( v × b ) b + 34 b · κ × v ( v · b ) b (cid:19) · ∂∂v This implies that the transformed one-form is given byΘ ǫ = A · dx + ǫ ( v · b ) b · dx (33)+ 12 | B | ǫ (cid:18) v × b · dv − ( v · b )[ ∇ b · ( v × b )] · dx − µ | B | τ b · dx (cid:19) = A · dx + ǫv k b · dx + ǫ µ [ dθ − ( R + 12 τ b ) · dx ] . Thus, the transformed Poisson bracket is the same as that given in Ref. 6. The transformedHamiltonian is given by H = 12 v · v (34) H = 0 (35) H = µ (cid:18) ∇ · b ) + 316 κ · κ + 14 b · ∇ ( ∇ · b )+ 116 tr[ ∇ b · ∇ b − ∇ b · ( ∇ b ) T ] − ∇ ln | B | · ∇ ln | B | + 14 κ · ∇ ln | B | + 14 | B | ∇ | B | (cid:19) + µ ( v · b ) | B | (cid:18)
18 tr[3 ∇ b · ∇ b − ∇ b · ( ∇ b ) T ] + 18 ( ∇ · b ) + 12 b · ∇ ( ∇ · b ) + 138 κ · κ − κ · ∇ ln | B | (cid:19) − ( v · b ) | B | (cid:18) κ · κ (cid:19) . Although the one-form, H , and H in this example are the same as those given by Parraand Calvo , H is in fact different. Using VEST , we have found that Parra and Calvo’s H has a numerically different value than the H in this example.23 III. CONCLUSION
We have reported, for the first time, on the automatic calculation of the guiding centerexpansion using a computer. In particular, we have implemented a novel Lie transform-basedalgorithm using the newly-developed Mathematica package
VEST and used it to derivetwo new representations of the guiding center equations of motion to the order relevantfor studying issues related to the physics of toroidal momentum conservation in tokamaks.By proceeding in this manner, we have avoided the pitfalls associated with hand-madealgebra errors and slashed the time required to perform the calculations from weeks tominutes. Readers interested in obtaining the Mathematica notebook we used to carry outour calculation can contact J. Squire via email at [email protected] are a number of opportunities for extending this work. Because our algorithmprovides the necessary tools to explore many representations of the guiding center expansion,it may be interesting to begin searching through different representations to find those withdesirable properties such as simple transformed equations of motion. Likewise, it may beinteresting to examine the time-validity of the transformed equations of motion in thesedifferent representations to see if some are worse than others. Kruskal’s theory guarantees1 /ǫ time-validity in all representations, but there may be representations that can do better.Yet another suitable application of our code would be pushing the calculation to higher orderthan previously calculated. For instance, it would be interesting to find µ and H . Finally,there should be no great difficulty in extending both our algorithm and our implementationin Mathematica using VEST to treat the gyrocenter transformation theory that forms thebackbone of modern gyrokinetic theory.
ACKNOWLEDGMENTS
The authors would like to express their gratitude to B. Faber for his help in editing thismanuscript. This work was supported by the U.S. Department of Energy under contractDE-AC02-09CH11466. 24 ppendix A: Elements of Exterior Calculus
In this appendix we will first list the basic identities commonly used when performingcalculations with the exterior calculus. Then we will give the component-form of the basicoperators d , L Y , i Y in the velocity phase space. For a much more thorough treatment ofthis topic, refer to Ref. 16.Let α k and β l be k - and l -forms on the manifold M , respectively. Let Y and Z be vectorfields on M . Then the following identities hold α k ∧ β l = ( − kl β l ∧ α k (A1)i Y ( α k ∧ β l ) = (i Y α k ) ∧ β l + ( − k α k ∧ (i Y β l ) (A2) d ( α k ∧ β l ) = ( d α k ) ∧ β l + ( − k α k ∧ ( d β l ) (A3) L Y ( α k ∧ β l ) = ( L Y α k ) ∧ β l + α k ∧ ( L Y β l ) (A4)i Y i Z = − i Z i Y (A5) L Y = d i Y + i Y d (A6) d L Y = L Y d (A7) dd = 0 . (A8)Let F : M → M and Φ : M → M be smooth mappings with a smooth inverses F − andΦ − . One example of this sort of mapping from the main text is Φ ψ , for fixed ψ , whoseinverse is Φ − ψ . The exterior calculus operations behave very well with respect to mappings.We will summarize this fact with a second list of identities. F ∗ Φ ∗ = (Φ ◦ F ) ∗ (A9) F ∗ ( α k ∧ β l ) = ( F ∗ α k ) ∧ ( F ∗ β l ) (A10) F ∗ (i Y α k ) = i F ∗ Y ( F ∗ α k ) (A11) F ∗ ( L Y α k ) = L F ∗ Y ( F ∗ α k ) (A12) F ∗ ( d α k ) = d ( F ∗ α k ) . (A13)When F = exp( Y ( ǫ )), with Y ( ǫ ) a vector field that tends to zero as ǫ →
0, we also have theasymptotic identities exp( − Y ( ǫ )) ∗ τ = exp( Y ( ǫ )) ∗ τ (A14)exp( Y ( ǫ )) ∗ τ = τ + L Y ( ǫ ) τ + 12 L Y ( ǫ ) T + ..., (A15)25here τ is an arbitrary tensor.The identities provided thus far, together with the fact that the wedge product is asso-ciative, are sufficient to verify all of the coordinate-independent manipulations of differentialforms in the main text. In order to perform exterior calculus using VEST it is also useful tohave component-based expressions for the operators d , i Y , and L Y . Actually, the relevantoperators for the sake of performing the guiding center calculation are d on functions, i Y onone-forms, and i Y d on both functions and one-forms.Let Y = Y xi ∂∂x i + Y vi ∂∂v i , where the indices are summed from i = 1 to i = 3 (although wedo so by habit, there is really no need to distinguish between covariant and contravariantindices in cartesian coordinates). Similarly, let α = α xi dx i + α vi dv i . Denote derivatives ofa scalar f with respect to the i ’th spatial argument and the i ’th velocity argument with f ,i and f ; j , respectively (note that ; does not denote a covariant derivative). Then we have d f = f ,i dx i + f ; i dv i (A16)i Y α = α xi Y xi + α vi Y vi (A17)i Y d f = f ,i Y xi + f ; i Y vi (A18)i Y d α = ( α xi,j − α xj,i ) Y xj dx i (A19)+ ( α xi ; j Y vj − α vj,i Y xj ) dx i + ( α vi,j Y xj − α xj ; i Y xj ) dv i + ( α vi ; j − α vj ; i ) Y vj dv i . Note that, by the identity given in Eq. (A6), the Lie derivative of a one-form, L Y α = d i Y α + i Y d α , can also be calculated using these component expressions. Appendix B: The Origin of Many Representations of The Guiding CenterExpansion
Suppose that a near-identity rearrangement of the phase space T ǫ : M → M is found thatrenders the transformed Lorentz vector field X ′ H = T ǫ ∗ X H gyrosymmetric. As explained byKruskal , at least one such transformation can be found using perturbation theory. In fact,as soon as one transformation in found, many more can be generated easily. This impliesthat there is much freedom in choosing T ǫ ; each choice leads to a different representation of26he guiding center equations of motion in the sense that X ′ H will be different in each case.In this appendix we will explain the origin of this freedom by completely characterizing it.First note that if F : M → M is a rearrangement of phase space (not necessarily near-identity) that commutes with the family of rearrangements Φ ψ that define the gyrosymmetry(see section V), i.e. F ◦ Φ ψ = Φ ψ ◦ F for each ψ ∈ R mod 2 π , then X ′′ H ≡ F ∗ X ′ H is alsogyrosymmetric. Indeed, Φ ∗ ψ X ′′ H = Φ − ψ ∗ F ∗ X ′ H (B1)= (Φ − ψ ◦ F ) ∗ X ′ H = ( F ◦ Φ − ψ ) ∗ X ′ H = F ∗ Φ ∗ ψ X ′ H = X ′′ H , where we have used the fact that X ′ H is gyrosymmetric and F ◦ Φ ψ = Φ ψ ◦ F . This tells usthat, given one of the near-identity rearrangements of phase space guaranteed by Kruskal,we can find many more by following the latter with any near-identity rearrangement of phasespace that commutes with Φ θ .In fact all of the rearrangements that fit into Kruskal’s theory can be found in thismanner. To be precise, suppose that R ǫ : M → M and Q ǫ : M → M are two near identityrearrangements that render X H gyrosymmetric, so that they fit into Kruskal’s theory. Then,by definition, Φ ∗ ψ ( R ǫ ∗ X H ) = R ǫ ∗ X H and Φ ∗ ψ ( Q ǫ ∗ X H ) = Q ǫ ∗ X H . Thus, both Q ǫ and R ǫ definesymmetry transformations on the original arrangement of phase space, ¯Φ Qψ = Q − ǫ ◦ Φ ψ ◦ Q ǫ and ¯Φ Rψ = R − ǫ ◦ Φ ψ ◦ R ǫ that leave X H invariant. Kruskal has proven in Ref. 11 that thesetwo symmetry transformations are in fact identical to all orders in ǫ ; ¯Φ ψ ≡ ¯Φ Qψ = ¯Φ Rψ . Itfollows then that the rearrangement F = Q ǫ ◦ R − ǫ commutes with Φ ψ . But this is preciselythe rearrangement that sends R ǫ ∗ X H into Q ǫ ∗ , which tells us that each representation ofthe guiding center equations of motion can be reached from a given one by applying anear-identity transformation that commutes with Φ θ .Note that if a rearrangement of the form exp( Y ), for some vector field Y , commutes withΦ θ , then it must be true that Y = h Y i . This is why one should expect complete freedom tochoose the h G n ( ǫ ) i for each G n ( ǫ ) appearing in the Lie transform ansatz given in Eq. (5).27 ppendix C: A General Formula For Inverting Exact Lagrange TensorsDefined On The Velocity Phase Space Step in our algorithm involves solving an algebraic equation of the form i X d Ξ = − β forthe vector field X given a non-degenerate two-form d Ξ and an arbitrary one-form β . In thisappendix we will present explicit expressions for the components of X = X xi ∂∂x i + X vi ∂∂v i interms of the components of Ξ = A i dx i + B i dv i and β = β xi dx i + β vi dv i .We proceed by making use of the linear isomorphism between the space of vector fieldsand the space of five-forms induced by the Liouville volume formΩ = 16 d Ξ ∧ d Ξ ∧ d Ξ . (C1)This isomorphism is given by X i X Ω. One can easily prove that this is an isomorphismusing the fact that the non-degeneracy of d Ξ implies that Ω is nowhere vanishing. Thereason this isomorphism is useful is that it is easier to find
X ≡ i X Ω than X . Indeed, uponwedge multiplying d Ξ ∧ d Ξ into both sides of the equation i X d Ξ = − β , we obtain X = − d Ξ ∧ d Ξ ∧ β. (C2)By explicitly calculating the right hand side of the last expression in components, and theninverting the isomorphism X i X Ω, we obtain X xn = − D (cid:18) ǫ ilk ǫ mjn β xm B k ; l ( B i,j − A j ; i ) (C3)+ ǫ lkm ǫ jin β vm A i,j B k ; l + 12 ǫ kim ǫ jln β vm ( B i,j − A j ; i )( B k,l − A l ; k ) (cid:19) X vn = 1 D (cid:18) ǫ jil ǫ kmn β vm A i ; j ( B k,l − A l ; k ) (C4)+ ǫ mji ǫ lkn β xm A i,j B k ; l + 12 ǫ mjl ǫ kin β xm ( B i,j − A j ; i )( B k,l − A l ; k ) (cid:19) , where the function D is defined by the relation Ω = D dx ∧ dx ∧ dx ∧ dv ∧ dv ∧ dv .28 ppendix D: Solving The Equation for S In order to complete step in our algorithm, the partial differential equation h − | B | i ξ α − | B | i ξ d S = h h i − | B | i ξ h α i (D1)must be solved for S given | B | , α , h , and the constraint h S i = 0. As is readily verified bygyroaveraging the equation, an equivalent condition on S is that it be chosen to eliminatethe non-zero gyroharmonics of the quantity h − | B | i ξ α − | B | i ξ d S .Let ν = h − | B | i ξ α . Using Eq. (14), we see that the Fourier expansion of ν ψ = Φ ∗ ψ ν isgiven by ν ψ = h ν i + ∞ X k =1 Π k ν cos( kψ ) + ¯Π k ν sin( kψ ) . (D2)Using the identities L ξ Π k S = k ¯Π k S and L ξ ¯Π k S = − k Π k S , we also see that the Fourierexpansion of Φ ∗ ψ | B | i ξ d S = | B | L ξ S ψ is given by | B | L ξ S ψ = ∞ X k =1 | B | k ¯Π k S cos( kψ ) − | B | k Π k S sin( kψ ) . (D3)Therefore, S must be given by S = ∞ X k =1 Π k S (D4)= ∞ X k =1 − ¯Π k ν | B | k . Appendix E: Motivation for Equations Defining The G n ( ǫ )In this appendix, we will motivate Eqs. (22) and (23).Suppose the one-form and Hamiltonian are expressed in the formΘ ǫ = ϑ + ǫϑ + ... + ǫ N ϑ N + ∞ X k =1 ǫ N + k α k (E1) H ǫ = H + ... + ǫ N − H N − + ∞ X k =1 ǫ N − k h k , with each θ i and H i gyrosymmetric and N >
1. Define Ξ = ϑ + ǫϑ + ... + ǫ N ϑ N . Considerchanging coordinates on the six-dimensional phase space using the mapping exp( − G ( ǫ ))29here G ( ǫ ) is some vector field that tends to zero as ǫ →
0. To be precise, considerexp( − G ( ǫ )) as the map that sends old coordinates into new coordinates. Note that we havenot specified the order of G ( ǫ ). Then the one-form and Hamiltonian expressed in the newcoordinates are Θ ′ ǫ = exp( − G ( ǫ )) ∗ Θ ǫ and H ′ ǫ = exp( − G ( ǫ )) ∗ H ǫ , respectively. Moreover,because G ( ǫ ) is small as ǫ →
0, we may expand these pushforward operators into a series ofLie derivatives, thereby obtaining the expressionsΘ ′ ǫ = Ξ + i G ( ǫ ) d Ξ + ǫ N +1 α + δ Θ (E2) H ′ ǫ = H + ... + ǫ N − H N − + i G ( ǫ ) d H + ǫ N − h + δ H , (E3)where δ Θ and δ H have not been displayed for simplicity, but are not necessarily higher orderin ǫ .These expressions for the transformed one-form and Hamiltonian are valid for any G ( ǫ )that tends to zero as ǫ →
0. Assuming the non-degeneracy conditions ND1-3 given in themain text, we will now use them to choose a specific G ( ǫ ) such that the transformed one-formand Hamiltonian are gyrosymmetric to one higher order than they were initially.First we demand that the one-form ν ≡ i G ( ǫ ) d Ξ + ǫ N +1 α should be gyrosymmetric up toan O ( ǫ N +1 ) exact differential. This implies ν + ǫ N +1 d S = h ν i , ori G ( ǫ ) d Ξ + ǫ N +1 α + ǫ N +1 d S = i h G ( ǫ ) i d Ξ + ǫ N +1 h α i . (E4)Note that this requirement forces the constraint h d S i = 0, as can be seen by gyroaveragingthe previous expression. Also note that this expression is identical to Eq. 22. Finally, notethat the gyroaverage of G ( ǫ ) completely disappears from the expression if we decompose G ( ǫ ) as h G ( ǫ ) i + ˜ G ( ǫ ). Indeed, we havei ˜ G ( ǫ ) d Ξ = − ǫ N +1 ( d ˜ S + ˜ α ) . (E5)Therefore, using the non-degeneracy conditions, we can impose the constraint i h G ( ǫ ) i d Ξ = ǫ N +1 γ , where γ is an arbitrary ǫ -independent gyrosymmetric one-form.By condition ND1 and ND2, for a given S , Eq. (E5) has a unique O ( ǫ N − ) solution for ˜ G ( ǫ )given by inverting the two-form d Ξ. Likewise, h G ( ǫ ) i = O ( ǫ N − ). Thus, G ( ǫ ) = O ( ǫ N − ).Using this result, it is straightforward to verify that, when G ( ǫ ) is chosen according toEq. (E4), δ Θ = O ( ǫ N +2 ). Therefore, Θ ′ ǫ is gyrosymmetric to one higher order in ǫ than it30as before applying the coordinate transformation. Even more, this conclusion still holdsfor any choice of the function S .This freedom in the selection of S can be used to make the transformed Hamiltonian H ′ ǫ gyrosymmetric to one higher order. To see this, first note that because G ( ǫ ) = O ( ǫ N − ), δ H = O ( ǫ N ). Thus, in order to make H ′ ǫ gyrosymmetric to one higher order in ǫ , all thatwe must do is ensure that the function f = i G ( ǫ ) d H + ǫ N − h is gyrosymmetric at O ( ǫ N − ).To see that this can be accomplished for a particular choice of S , define the vector field X H using the formula i X H d Ξ = − d H . This definition makes sense in light of condition ND1.In terms of X H , f can be written f = i X H i G ( ǫ ) d Ξ + ǫ N − h . Then, by Eq. (E4), we have f = ǫ N +1 i X H ( γ − ˜ α − d S ) + ǫ N − h . (E6)But by condition ND3, the leading-order contribution to X H is given by | B | ǫ ξ , which implies f = ǫ N − | B | i ξ ( γ − ˜ α − d S ) + ǫ N − h + O ( ǫ N ) . (E7)Thus, the condition that all gyroharmonics of f should be zero at O ( ǫ N − ) is | B | i ξ ( ˜ α + d S ) = ˜ h , (E8)which is precisely Eq. (23). Appendix F: µ for general magnetic geometry The expression for µ we have derived using VEST disagrees with the result given in Ref.28. However, upon close examination of the expression for µ given in Ref. 28, which isexpressed as a sum µ = | B | P n,m =0 v n k v m ⊥ a nm , we have identified a probable typographicalerror in the expression for a . Partially and temporarily adopting the notation used in Ref.28, the terms − ( n · ∇× b )( n · ∇× n − n · ∇× n ) and ( n · ∇× b )( n · ∇× n − n · ∇× n )were not combined in the obvious way. Given the simplicity of this simplification, it seemslikely that one or both of these terms was copied incorrectly.Define the vectors η = ( ∇ b ) · v and λ = v · ∇ b . Also define the scalars µ = ( v × b ) · ( v × b )2 | B | , v k = b · v , and γ = v · ∇ ln | B | . The expression for µ we have derived and verified using31 EST is µ = − µ | B | η · η − v k µ | B | λ · ∇ ln | B | (F1)+ 10 v k µ | B | ∇ ln | B | · η − v k µ | B | κ · η + 71 µ | B | ( κ · v ) − µ | B | γ + 71 v k µ | B | λ · κ − µ | B | λ · λ + 25 µ | B | λ · η + 5 v k | B | ( κ · v )( λ · v ) − v k | B | ( γ )( λ · v )+ 18 | B | ( λ · v ) + 217 v k µ | B | ( κ · v )( ∇ · b ) − v k µ | B | ( γ )( ∇ · b )+ 23 µ | B | ( λ · v )( ∇ · b ) − v k µ | B | ∇ ( b · v ) + 5 v k µ | B | v · ∇ ( ∇ · b )+ µ | B | vv : ∇∇| B | + 5 v k µ | B | bb : ∇∇ ( b · v )+ v k | B | vv : ∇∇ ( b · v ) − v k µ | B | b · ∇ ( ∇ · b ) − v k | B | η · η + 20 v k µ | B | κ · ∇ ln | B | − v k µ | B | κ · κ + 5 v k | B | ( κ · v ) − v k | B | ( κ · v )( γ ) − v k | B | λ · λ + 5 v k | B | λ · η − v k | B | ( λ · v )( ∇ · b ) − v k µ | B | ( ∇ · b ) + 55 v k µ | B | Tr( ∇ b · ( ∇ b ) T ) − v k µ | B | Tr( ∇ b · ∇ b ) + 5 v k | B | bv : ∇∇ ( b · v ) + 5 v k | B | λ · κ + 5 v k | B | ( κ · v )( ∇ · b ) + 5 v k | B | bb : ∇∇ ( b · v ) + 25 v k | B | κ · κ − µ | B | b · ∇ ( ∇ · b ) − µ | B | ∇ | B | + 15 µ | B | ∇ ln | B | · ∇ ln | B |− µ | B | κ · ∇ ln | B | − µ | B | κ · κ − µ | B | ( ∇ · b ) + 11 µ | B | Tr( ∇ b · ( ∇ b ) T ) − µ | B | Tr( ∇ b · ∇ b ) . In this expression, x and v are the untransformed position and velocity variables. In par-ticular, in the coordinate system used in this expression, the Lorentz force takes the usualform ˙ v = | B | v × b and ˙ x = ǫv . 32 EFERENCES F. Jenko, W. Dorland, M. Kotschenreuther, and B. N. Rogers, Phys. Plasmas , 1904(2000). G. Dif-Pradalier, P. H. Diamond, V. Grandgirard, Y. Sarazin, J. Abiteboul, X. Garbet,P. Ghendrih, G. Latu, A. Strugarek, S. Ku, , and C. S. Chang, Phys. Plasmas , 062309(2011). T. G. Northrop,
The Adiabatic Motion of Charged Particles , Interscience tracts on physicsand astronomy (Interscience Publishers, 1963). A. J. Brizard and N. Tronko, (2012), arXiv:1205.5772. R. G. Littlejohn, in
Fluids and Plasmas: Geometry and Dynamics , Contemporary math-ematics, Vol. 28, edited by J. E. Marsden (American Mathematical Society, 1984) pp.151–167. A. Brizard,
Nonlinear Gyrokinetic Tokamak Physics , PhD dissertation, Princeton Univer-sity, Department of Astrophysical Sciences (1990). J. Krommes, Ann. Rev. Fluid Mech. , 175 (2012). F. I. Parra and P. J. Catto, Phys. Plasmas , 056106 (2010). R. G. Littlejohn, Phys. Fluids , 1730 (1981). R. G. Littlejohn, J. Plasma Phys. , 111 (1983). M. Kruskal, J. Math. Phys. , 806 (1962). S. M. Omohundro,
Geometric Perturbation Theory in Physics (World Scientific PublishingCo. Pte. Ltd., Singapore, 1986). J. Squire, J. W. Burby, and H. Qin, “Vest: abstract vector calculus simplification inmathematica,” (2013), Comp. Phys. Comm. (submitted). M. Kruskal, “The gyration of a charged particle,” Project Matterhorn Report PM-S-33(NYO-7903) (Princeton University, 1958). J. Larsson, Physica Scripta , 342 (1986). R. Abraham and J. Marsden,
Foundations of Mechanics , AMS Chelsea publishing (AMSChelsea Pub./American Mathematical Society, 1978). J. M. Lee,
Introduction to Smooth Manifolds , Graduate Texts in Mathematics (Springer-Verlag New York, 2003). J. W. Burby and H. Qin, Phys. Plasmas , 052106 (2012).33 H. Qin, in
Fields Institute Communications , Vol. 46, edited by D. Levermore, T. Passot,C. Sulem, and P. Sulem (American Mathematical Society, 2005) pp. 171–192. For some quick intuition regarding this ansatz, recall that the time-advance map associatedto an arbitrary vector field Y : M → T M is given by exp( tY ). Thus, given x ∈ M , T ǫ ( x )is found by sequentially flowing along the vector fields G n ( ǫ ) for − x . L. De Guillebon and M. Vittot, (2013), arXiv:1211.5792. N. Bogoliubov,
Asymptotic Methods in the Theory of Non-Linear Oscillations , Interna-tional monographs on advanced mathematics and physics (Hindustan, 1961). R. G. Littlejohn, J. Math. Phys. , 742 (1982). Demanding the transformed d ϑ ǫ to be gyrosymmetric is equivalent to demanding that thetransformed ϑ ǫ be equal to the sum of a gyrosymmetric one-form and an arbitrary closedone-form. This follows from the fact that d ϑ ′ ǫ is left unchanged upon replacing ϑ ′ ǫ with ϑ ′ ǫ + α for an arbitrary closed one-form α , i.e. d α = 0. F stands for “fibered”. In spite of its simplicity and utility, this formula only seems to have been noticed recently.It can be found in the literature in Ref. 18. It was also independently discovered by ZhiYu , but never published. ND stands for “non-degenerate”. B. Weyssow and R. Balescu, J. Plasma Phys. , 449 (1986). F. Parra and I. Calvo, Plasma Phys. Control. Fusion , 045001 (2011).30