A&A 428, 1-19 (2004)
DOI: 10.1051/0004-6361:20034208

Covariant gyrokinetic description of relativistic plasmas

A. Beklemishev ¹ - M. Tessarotto ²

1 - Budker Institute of Nuclear Physics, 630090 Novosibirsk, Russia
2 - Department of Mathematical Sciences, University of Trieste, 34127 Trieste, Italy

Received 18 August 2003 / Accepted 20 July 2004

Abstract
A fundamental aspect of many plasma-related astrophysical problems is the kinetic description of magnetized relativistic plasmas in intense gravitational fields, such as in accretion disks around compact gravitating bodies. The goal of this paper is to formulate a gyrokinetic description for a Vlasov-Maxwell plasma within the framework of general relativity. A closed set of relativistic gyrokinetic equations, consisting of the collisionless gyrokinetic equation and corresponding expressions for the four-current density, is derived for an arbitrary four-dimensional coordinate system. General relativity effects are taken into account via the tetrad formalism. The guiding-center dynamics of charged particles and the gyrokinetic transformation are obtained accurate to the second order of the ratio of the Larmor radius to the nonuniformity scale length. The wave terms with arbitrary wavelength $(k\rho_{\rm L}\sim1)$are described in the second-order (nonlinear) approximation with respect to the amplitude of the wave. The same approximations are used in the derivation of the gyrophase-averaged Maxwell equations. The derivation is based on the perturbative Lagrangian approach with a fully relativistic, four-dimensional covariant formulation. Its results improve on existing limitations of the gyrokinetic theory.

Key words: plasmas - relativity - gravitation

1 Introduction

The gyrokinetic theory is an adequate, yet simplified, description of collisionless magnetized plasmas. It strives to utilize the conservation of the magnetic moment to reduce the number of effective degrees of freedom of charged particles. In this respect it is similar to the drift-kinetic approximation. However, it goes further than the drift approximation, since it allows short-wavelength (comparable to the Larmor radius, $\rho_{{\rm L } }$) perturbations of the background fields.

From the mathematical viewpoint, the dynamics of Hamiltonian systems can be expressed, in principle, in arbitrary hybrid (i.e., non-canonical) state variables. The search of phase-space transformations yielding "simpler'' equations of motion has long motivated theoretical research in physics and mathematical physics. Among such transformations, a particular case is provided by the gyrokinetic transformation, yielding hybrid variables defined in such a way that one of them, the so-called gyrophase angle, is ignorable. This transformation can only be introduced when particles move in a sufficiently strong magnetic field.

For a certain class of problems in plasma physics and astrophysics, the existing limitations of the standard gyrokinetic theory (Hahm et al. 1988; Brizard & Chan 1999; Cooper 1997; Littlejohn 1979,1984; Boozer 1996) make its use difficult or impossible. In particular, this involves the description of experiments in which the electric field may become comparable in strength to the magnetic field (so that the drift velocity becomes relativistic), and the study of relativistic plasma flows in gravitational fields, which are observed or assumed to exist in accretion disks and related plasma jets around neutron stars, black holes, and active galactic nuclei (Frank et al. 1992). The finite Larmor radius effects and the influence of short wavelength electromagnetic perturbations are also expected to play a fundamental role in particle dynamics. The present paper aims at extending the applicability range of gyrokinetic theory to include these features, and at providing an effective tool for use in plasma-physics-related problems in astrophysics.

The conventional non-relativistic gyrokinetic theory (Hahm et al. 1988; Littlejohn 1979) has been generalized to include only some relativistic effects, which are essential for the description of the confinement of fusion products in thermonuclear reactors (Brizard & Chan 1999; Cooper 1997; Littlejohn 1984; Boozer 1996), namely, the particle itself is considered relativistic, while its drift velocity is not. This deficiency has been pointed out by Pozzo & Tessarotto (1998), who attempted to write the theory in the reference frame moving with the fluid velocity, where all drifts should be first order in the expansion parameter, and the question about relativistic drifts should disappear. However, it proved to be more rewarding to create the covariant theory (Beklemishev & Tessarotto 1999), as it makes use of the term "relativistic'', i.e., with relation to any reference frame, in its full sense. This approach is inherently more general and symmetric than the non-relativistic treatment and its "relativistic'' generalizations. As a result, even the non-relativistic limit of the theory is found to have a broader applicability range than the standard derivation. In particular, covariant formulation allows self-consistent inclusion of the gravitational field, which makes the theory suitable for astrophysical applications.

In our previous paper (Beklemishev & Tessarotto 1999) the development of the covariant theory (derivation of the gyrokinetic transformation) was carried out through first order in the expansion parameter, and without the wave fields. It has been shown that the fundamental adiabatic invariant of a particle, i.e., its magnetic moment, as well as the ignorable gyrophase, can be defined as Lorentz invariant, independent of the reference frame. Also, the derivation itself becomes algebraically simple due to the internal symmetry of the electromagnetic field in the four-dimensional formulation. Covariant formulation allows the derived equations to be easily rendered for any coordinate system in four-dimensional Riemann space-time. It is important for astrophysics applications, as well as for problems where description in curvilinear magnetic coordinates is convenient. Although we do not use canonical variables, any coordinate system compatible with canonical variables can be imposed on the final result. Here it is necessary to note that the so-called relativistic drift Hamiltonian (Boozer 1996), and the gyrokinetic theories based on it (Cooper 1997), are limited by the applicability range of the drift equations given by Northrop (1963), i.e., assume non-relativistic drifts. The covariant drift equations require no such assumptions.

Covariant descriptions of guiding center dynamics have been published earlier, in the most complete form by Fradkin (1978). He first used the decomposition of motion of a charged particle into two quasi-independent motions in orthogonal invariant planes of the electromagnetic field tensor. Equations of motion were derived in the first order in inhomogeneity, without the wave fields. However, explicit use of the proper time for parameterization of the trajectory makes his expressions difficult to adapt to the gyrokinetic description, where parameterization should be applicable to the guiding center trajectory. Later, Beklemishev & Tessarotto (1999) showed that without such parameterization the choice of the decomposition planes becomes a matter of convenience rather than necessity.

The goal of the present work is two-fold. First, adopting the covariant gyrokinetic approach developed earlier (Beklemishev & Tessarotto 1999) we intend to describe particle dynamics in the presence of non-uniform and strong electromagnetic (EM) and gravitational fields, accurate to second order in the Larmor radius, and treating consistently the effect of short-wavelength EM perturbations. The basic features of the approach, allowing the relativistic drifts typical of many astrophysical plasmas, are the assumption of a non-vanishing parallel electric field and the description of particle dynamics in the field-related basis tetrad, which allows a simpler formulation of the relativistic equations of motion. Second, the perturbative theory is employed to derive the relativistic Vlasov-Maxwell equations expressed in gyrokinetic variables accurate through second order in  $\varepsilon$ (which is the ratio of the gyro-radius/gyrotime to the equilibrium gradient length/time) while the wave fields are taken into account in the second order in the amplitude. The gyrophase-averaged expression of the four-current density appearing in Maxwell equations is evaluated with the same accuracy.

In contrast to non-covariant formulations, since the space and time variables are treated in a parallel fashion, the applicability limits of gyrokinetic theory are extended - the wave fields can be fast oscillating in time ( $\omega \rho_{{\rm L}}/c\sim1$), as well as in space ( $k\rho_{ {\rm L}}\sim1$), while the customary assumption $k_{\Vert}\rho_{{\rm L} }\ll1$ is not required. As a result, we find a contribution to the magnetic moment which is proportional to the usually neglected parallel component of the wave-field. An important feature of this work is its relative simplicity, which allows us to present results with the gyrophase-averaging performed explicitly. (The usual way of doing this is by leaving the averaging along the particle trajectory to the reader's computer; see, for example, Brizard & Chan 1999.) However, currently there is also a limitation linked to our approach, namely, in the case of a vanishing background parallel electric field, only certain polarizations of the wave field are allowed. Namely, the wave field should be describable by just two components of the vector potential, both orthogonal to the magnetic field (in the reference frame where the electric field vanishes), or have a large wavelength, to be described as part of the main field. We intend to eliminate this restriction in the future.

Another work dealing with the covariant representation of the gyrokinetic theory (for a flat space-time) is the Ph.D. Thesis of Boghosian (1987). It extends the covariant description of charged particle guiding center motion developed by Fradkin (1978). In line with non-covariant treatments Boghosian assumes a vanishing parallel component of the electric field and the usual gyrokinetic ordering for admissible wave forms, and thus has to modify Maxwell equations to maintain the assumed properties. His derivation is also conducted only to first-order both in nonuniformity and in the wave-amplitude.

The structure of the paper is as follows: The hybrid Hamilton variational principle is formulated in Sect. 2. In Sect. 3 we present the definition of the basis tetrad and its link to the Faraday tensor. In Sect. 4 the relativistic gyrokinetic transformation is derived to second order in  $\varepsilon$, and in the wave amplitude. The Euler equations for the gyrocenter trajectories are formulated in Sect. 5. Maxwell equations with the current and charge densities, expressed in terms of the distribution function as a function of the gyrokinetic variables, and integrated in the gyrophase, are presented in Sect. 6. In the Conclusion the results of the paper are summarized.

2 The variational principle

Construction of the gyrokinetic transformation for the one-particle Hamiltonian system requires a description of the particle dynamics in phase-space. It is possible to redefine the Hamilton variational principle in such a way that it becomes suitable for the phase-space description of motion in non-canonical variables (Cary & Littlejohn 1983). We refer to this as "the hybrid Hamilton variational principle''. In the framework of general relativity the four-dimensional formulation of the Hamilton variational principle is well known (Landau & Lifshits 1975). It makes it possible to describe the particle trajectory (or its world-line) as a four-dimensional extremal curve. In this section a hybrid form of this variational principle, suitable for the description of particle dynamics in the seven-dimensional phase-space (Beklemishev & Tessarotto 1999), is recalled.

The relevant tensor notations are standard. Thus, $g_{\mu\nu}$ denotes the metric tensor components, characterizing the coordinate system (and the underlying space-time structure). It also provides the connection between the co- and countervariant components of four-vectors (henceforth referred to as 4-vectors) $A_{\mu}=g_{\mu\nu}A^{\nu}$. The signature of  $g_{\mu\nu}$is assumed to be $\langle 1+,3-\rangle$, while the invariant interval ${\rm d}s$ is defined as

$${\rm d}s^{2}=g_{\mu\nu}{\rm d}x^{\mu}{\rm d}x^{\nu},$$ (1)

where the Greek indices are assumed to go through $\mu,\nu=0,...,3$.

The four-velocity $\tens{u}$ is then defined in terms of its countervariant components according to

$$u^{\mu}=\frac{{\rm d}x^{\mu}}{{\rm d}s}\cdot$$ (2)

The scalar product of any 4-vector $\tens{a}$ by itself, or its square length, $a^{2}=a_{\mu}a^{\mu},$ is not sign-definite. If its sign is negative, we shall call the 4-vector "space-like'', while if it is positive, the 4-vector is "time-like''. Unit 4-vectors are those with square lengths equal to $\pm1$. In these terms, because of (1) $, \tens{u}$ is the unit, time-like 4-vector

$$u_{\mu}u^{\mu}=1.$$ (3)

The motion or the world-line of a particle with rest-mass m_a and charge q_a in prescribed fields in phase space can be found from the Hamilton variational principle $\delta S=0$, as the extremal of the functional

$$S{=}\int {Q_{\mu }{\rm d}x^{\mu }=}\int (qA_{\mu }(x^{\nu })+u_{\mu }) {\rm d}{x^{\mu }},$$ (4)

where q=q_a/m_ac². Ideally, a variational principle in phase space should yield the standard relativistic equation of motion of the particle, plus the relationship between velocity and displacement, Eq. (2). As shown by Beklemishev & Tessarotto (1999), this is indeed the case for the functional (4), but only for variations of $u_{\mu }$ occurring on the hypersurface $u_{\mu }u^{\mu }=1$. In this sense the phase space is only seven-dimensional. The above variational principle is independent of the parameterization of the world lines.

Finally, we note that the above variational principle, as any Lagrangian variational principle, is gauge invariant, i.e., any change of the generating differential 1-form of the type

$$Q_{\mu}{\rm d}x^{\mu}\rightarrow Q_{\mu}{\rm d}x^{\mu}+{\rm d} F,$$ (5)

where F is an arbitrary smooth function of coordinates $x^{\mu}$ and zⁱ, does not change the extremal curve.

3 The Faraday tensor and the basis tetrad

The use of a tetrad of four-vectors is a well-known type of description of curved space-time (Landau & Lifshits 1975), which makes it possible to replace complex derivatives of the metric tensor by derivatives of the tetrad components. We adopt this approach below, while also finding it convenient to link our tetrad to properties of the electromagnetic field, or the Faraday tensor, $F_{\mu\nu}=\partial A_{\nu}/\partial x^{{\mu}}-\partial A_{\mu}/\partial x^{{\nu}}.$

We first introduce an orthogonal basis of unit 4-vectors $( \tens{\tau,l,l^{\prime },l^{\prime\prime}})$ so that the last three 4-vectors are space-like, and

$$\sqrt{-g}e_{\varsigma \lambda \mu \nu }\tau ^{\varsigma }l^{\lambda }l^{\prime \mu }l^{\prime \prime \nu }=1,$$ (6)

where $\sqrt{-g}e_{\varsigma \lambda \mu \nu }$ is the purely antisymmetric tensor. This last requirement is necessary for distinguishing between the left and right basis configurations. As long as the orientation of $( \tens{\tau,l,l^{\prime },l^{\prime\prime}})$ is arbitrary, it can be constructed with 6 free parameters (3 pure space rotations and 3 space-time rotations). The choice of orientation is assumed to depend only on $\tens{x}$.

A special choice of orientation links the basis $( \tens{\tau,l,l^{\prime },l^{\prime\prime}})$ to the electromagnetic field tensor, $F_{\mu \nu }$. We shall refer to this choice of basis as field-related and use it henceforth. With this choice the $(\tens{{l}^{\prime}{,l}^{\prime\prime}})$-plane coincides with the space-like invariant plane of the antisymmetric tensor $F_{\mu \nu }.$ In some sense this corresponds to the magnetic coordinates of non-relativistic treatments. Fradkin (1978) used these invariant planes of the electromagnetic field tensor for the decomposition of motion of a charged particle.

Any non-degenerate antisymmetric tensor of rank 4, $F_{\mu \nu }$ in particular, has two orthogonal invariant planes. In our case one of them is space-like, i.e., any vector belonging to it has negative length, while the other contains both time-like and space-like vectors. Each of the invariant planes has an associated eigenvalue. This means that if $(\tens{{l}^{\prime}{,l}^{\prime\prime}})$ is the first invariant plane, then $(\tens{l,\tau})$ is the other, and if H and E are the eigenvalues of $F_{\mu \nu }$, then

$$F_{\mu \nu }l^{\prime \prime \nu }=Hl_{\mu }^{\prime },\ \ F_{\mu \nu }l^{\prime \nu }=-Hl_{\mu }^{\prime \prime },$$ (7)

$$F_{\mu \nu }l^{\nu }=E\tau _{\mu },\ \ F_{\mu \nu }\tau ^{\nu }=El_{\mu }.$$ (8)

(Both equalities in the second line have the same sign due to the $\langle 1+,3-\rangle$signature of the metric tensor. The physical meaning is that H and E are the magnetic and electric field strengths in the reference frame where the electric and the magnetic fields are parallel. The signs are chosen in such a way that in this frame $\vec{E}=E\vec{l},\ \vec{B}=H\vec{l},$ if $\tens{l} =\left( 0,\vec{l}\right)$ is a pure-space vector. For this reason the sign of H is opposite to that in our previous paper (Beklemishev & Tessarotto 1999), consequently the sign of the magnetic moment is also the opposite.) With this choice one can assume $H\geq 0$, without loss of generality, while the sign of E is determined by the invariance of the scalar product $\vec{E}\cdot \vec{B}=EH.$

Due to these properties of the basis, the Faraday tensor can be fully expressed via its components as

$$F_{\mu \nu }=H\left( l_{\nu }^{\prime }l_{\mu }^{\prime \prim... ...t) +E\left( l_{\mu }\tau _{\nu }-l_{\nu }\tau _{\mu }\right) .$$ (9)

The bivector combinations entering the above expression

$$b_{\mu \nu }=l_{\nu }^{\prime }l_{\mu }^{\prime \prime }-l_{\... ... },\;\;c_{\mu \nu }=l_{\mu }\tau _{\nu }-l_{\nu }\tau _{\mu },$$ (10)

can be viewed as linear operators, which combine projection on an invariant plane with an orthogonal turn in it, i.e., each carries full information about an oriented plane. Then

$${F}=H{\rm b}+E{\rm c,}$$

while the pure projection operators $\mathcal{O}^{(a)},\mathcal{O}^{(b)}$ introduced by Fradkin (1978) can be expressed as

$$\mathcal{O}^{(a)}=-{\rm bb},\;\mathcal{O}^{(b)}={\rm cc}.$$

Since $\mathcal{O}^{(a)}+\mathcal{O}^{(b)}=\tens{I},$ and $\tens{bc}=0,$ it is easy to check that the inverse of the Faraday tensor, $D^{\nu \mu }=(F_{\mu \nu })^{-1}$ is

$${\rm D}=-\frac{1}{H}{\rm b}+\frac{1}{E}{\rm c}.$$ (11)

It exists, provided that $F_{\mu \nu }$ is not degenerate, i.e., $EH\neq 0$. This inverse of the Faraday tensor can also be expressed via its dual, namely,

$$\mathcal{F}^{\mu \nu }\equiv \frac{1}{2\sqrt{-g}}{\rm e}^{\mu \nu \lambda \zeta }F_{\lambda \zeta }.$$ (12)

Indeed,

$$D^{\nu \mu }=\frac{1}{EH}\mathcal{F}^{\nu \mu }.$$ (13)

4 The relativistic gyrokinetic transformation

The object of gyrokinetics is to introduce a new set of phase-space variables (called the gyrokinetic variables) such that the variable describing the rotation angle along the Larmor orbit (i.e., the gyrophase $\phi$) becomes ignorable. This happens, by definition, when the Lagrangian (or, more generally, the functional) is independent of $\phi$. Once an ignorable variable is found, the number of corresponding Euler equations is reduced by one, and the new variables allow simplified numerical calculations, as the motion is effectively integrated over the fast Larmor rotation. The one-to-one transformation from the original set of phase-space variables $(x^{\mu },u^{\nu })$ to the gyrokinetic variables is called the gyrokinetic transformation. In what follows to find these variables we use the Lagrangian perturbative approach, which is equivalent (in broad terms) to the Lie-transform method, though more direct.

Indeed, we can search for the necessary transformation, representing it as a Taylor-series expansion in small parameters (inhomogeneity and wave amplitude). Then, in each order we can require the transformed functional to be independent of the gyrophase (the ignorable variable), and thus find the transformation and the new phase-independent functional order by order. For this we should be sure that 1) the solution in the necessary order exists; and 2) that the procedure converges. The existence of the solution depends on the choice of the expansion form and is proven operationally, i.e., by direct construction. Then the true dynamical system becomes replaced by the new artificial one, described by the phase-independent Lagrangian of that order. There is also the problem of convergence, i.e., how accurate the gyrokinetic description is. This problem is common to all gyrokinetic (and drift) theories and has not yet been explored in detail. From the physical viewpoint, the artificial system has no cyclotron resonances, while the real one has them, but their width decreases with the small parameters. The gyrokinetic approach cannot work within this thin but possibly interlinking Arnold web of resonances, where the adiabatic invariants are not conserved. Outside the resonances one can claim that since the difference between the real and the artificial Lagrangians is limited and tends to zero within some suitable measure, the behaviour of the gyrokinetic system is similar to that of the real systems.

We begin the search for the gyrokinetic transformation by a suitable definition of the small parameters, i.e., we assume that the curvature radius of the space-time and the gradient lengths of the background electromagnetic fields are much larger than the Larmor radius characterizing the particle path. However, we allow for the existence of wave fields with sharp gradients [ $k\rho _{{\rm L}}\sim O(1),$ including $k_{\Vert }\rho _{ {\rm L}}\sim O(1)$], and rapidly varying in time [ $\omega \rho _{\rm L}/c\sim O(1)$], while such fields are assumed suitably smaller in strength than the background field. (We stress that, unlike in conventional formulations of gyrokinetic theory, this type of ordering is required in a covariant theory due to the reference-frame dependence of the ordering assumptions involving space and time scale lengths of the perturbations.) For this purpose we introduce the ordering scheme following the notation of Littlejohn (1984):

$${Q_{\mu }{\rm d}x^{\mu }}=\left\{u_{\mu }+q\left(\frac{1}{\va... ...on }A_{\mu }+\lambda a_{\mu }\right)\right\}{\rm d}{x^{\mu }},$$ (14)

where $\varepsilon$ and $\lambda$ are formal small parameters (they should be set to 1 in the final results) allowing distinction between the large-scale background field $A_{\mu },$ and the wave fields given by $a_{\mu }$. We search for the gyrokinetic transformation $(y^{i})\equiv (x^{\prime \alpha },\phi ,\hat{\mu},{u_{\parallel }})\leftrightarrow (x^{\alpha },u^{\beta })$ in the form of an expansion in powers of $\varepsilon$:

$$x^{\nu }=x^{\prime \nu }+\sum_{s=1}\varepsilon ^{s}r_{s}^{\nu }(y^{i}),$$ (15)

where $\phi$ is the ignorable phase variable (gyrophase), $\hat{\mu}$ and ${ \ \ u_{\parallel }}$ represent two more independent velocity characteristics (to be defined later), $x^{\prime \nu }$ is the 4-vector guiding center position, $\tens{{r}_{\rm s}}$ are arbitrary 4-vector functions of the new variables (yⁱ) to be determined. We require that $\tens{{r}_{\rm s}}$ are purely oscillatory in $\phi$, i.e., the $\phi$-averages of $\tens{{r}_{\rm s}}$ are zero, as part of the $x^{\prime \nu }$- definition. Note that the above descriptions of the new variables $(x^{\prime \alpha },\phi ,\hat{\mu},{\ u_{\parallel }})\$will acquire precise mathematical meaning only as a result of the search for the gyrokinetic transformation.

This search consists in applying the expansion (15) to the fundamental 1-form (14), and imposing the requirement that it is independent of $\phi$ in each order. Since the 1-form itself depends on $\lambda,$ expansion in this parameter will also appear.

Since the difference between $x^{\alpha}$ and $x^{\prime\alpha}$ is assumed small compared to the gradient lengths, the background field components (unlike the wave terms) can be expanded in the Taylor series

$$A_{\mu}=A_{\mu}^{\prime}+\varepsilon r_{1}^{\nu}\frac{\partia... ...\prime}}{\partial x^{\prime{\nu}}}\right) +o(\varepsilon^{2}),$$ (16)

where $A_{\mu}^{\prime}=A_{\mu}(\tens{x}^{\prime}).$ Furthermore, we assume that the gravitational field and the metric coefficients $g_{\mu\nu}$ have the same order of the gradient length, and thus can also be expanded around $\tens{x}^{\prime}.$ However, it is not necessary. It is sufficient to use the primed metric tensor for lifting and lowering indexes of primed tensor components. Finally, the ${\rm d}x^{\mu}$ vector components in the new variables (yⁱ) look like

$${\rm d}x^{\mu}=\frac{\partial x^{\mu}}{\partial y^{i}}{\rm d}... ...varepsilon }\frac{ \partial x^{\mu}}{\partial\phi}{\rm d}\phi.$$ (17)

Here $\left( {\rm d}x^{\mu}\right) _{\phi}$ means the differential of $x^{\mu}$ with respect to all other variables of the set $(x^{\prime \alpha },\phi ,\hat{\mu},{\ u_{\parallel }})\$ keeping $\phi$ constant. As usual, the factor  $1/\varepsilon$ is introduced to denote the fast dependence characterizing the Larmor rotation.

Equation (15) define only four out of seven transformation rules for the gyrokinetic transformation $(x^{\prime\alpha},\phi,\hat{\mu },{\ u_{\parallel}})\rightarrow(x^{\alpha},u^{\mu})$. The rest are chosen in the following way.

First, we express the four-velocity in the new variables as

$$u_{\mu }=w\left( l_{\mu }^{\prime }\cos \phi +l_{\mu }^{\prime \prime }\sin \phi \right) +\bar{u}_{\mu },$$ (18)

$$\bar{u}_{\mu }={u_{\parallel }}l_{\mu }+u_{o}\tau _{\mu },$$ (19)

which can be also regarded as the definition for the gyrophase $\phi$:it is defined as an angle in the velocity subspace, where we introduce the cylindrical coordinate system. This definition is covariant with respect to transformations of the space-time coordinate system, which may change the vector components, but not the vectors themselves. Note that the tetrad is linked to the position of the particle, x, rather than  ${x}^{\prime }$, so that its components should be transformed to  ${x}^{\prime }$ as well.

Furthermore, we assume that $w,{u_{\parallel }}$ and u_o are independent of $\phi$. The validity of this assumption is justified by the existence of the solution (at least for a non-degenerate Faraday tensor). Then, the $\phi$-independent part of the 4-velocity $\tens{{\bar{u}}}$ is not completely arbitrary, but satisfies certain restrictions following from the requirement $u_{\mu }u^{\mu }=1$ for all $\phi$:

$$u_{o}^{2}=1+w^{2}+{u_{\parallel }}^{2}.$$ (20)

Any two of three scalar functions w,u_o or ${ \ \ u_{\parallel }}$ can be considered independent characteristics of velocity, while the third can be expressed via Eq. (20).

Let us substitute expressions (15), (17) into (14) and drop terms of order  ${\varepsilon ^{2}}$ or smaller:

$${Q_{\mu }{\rm d}x^{\mu }}=\left( \frac{q}{\varepsilon }A_{\mu... ...}^{\mu }}{\partial \phi }\right) {\rm d}\phi +o(\varepsilon ).$$ (21)

Following Beklemishev & Tessarotto (1999) we construct the gauge-modified functional

$${\delta G}^{\prime } {Q_{\mu }{\rm d}x^{\mu }}-{\rm d}{}\left( \frac{1 }{2}q(A_{\mu }+... ... }-A_{\mu }^{\prime })({\rm d}{r_{1}^{\mu }+\varepsilon }{\rm d}{r_{2}^{\mu })}$$ =

$$-\frac{1}{2}q({r_{1}^{\mu }+\varepsilon r_{2}^{\mu }})({\rm d}A_{... ...rac{\partial r_{2}^{\mu }}{\partial \phi }\right) {\rm d}\phi +o(\varepsilon ),$$ (22)

so that after using Eq. (16) to represent $A_{\mu }$ and

$${\rm d}A_{\mu }=\left( \frac{\partial A_{\mu }^{\prime }}{\pa... ...x^{\prime {\varsigma }}}\right) {\rm d} \phi +o(\varepsilon ),$$ (23)

we get

$${\delta G}^{\prime }=\left( \frac{q}{\varepsilon }A_{\mu }^{\... ...ght) {\rm d}{r_{1}^{\mu }}+\Phi {\rm d} \phi +o(\varepsilon ),$$ (24)

where

$$\Phi =\varepsilon \left( r_{1}^{\nu }\frac{\partial u_{\mu }}... ...{\nu }^{\prime }}{\partial x^{\prime {\mu }}}\right] \right) ,$$ (25)

${\rm d}{r_{1}^{\mu }}$ is given by Eq. (17), and $F_{\mu \nu }$is the tensor of the electromagnetic field at $\tens{x}^{\prime }$:

$$F{_{\mu \nu }}=\frac{\partial A_{\nu }^{\prime }}{\partial x^... ...{\partial A_{\mu }^{\prime }}{\partial x^{\prime {\nu }}}\cdot$$ (26)

We also decomposed the tetrad vectors around ${x}^{\prime },$ $u_{\mu }\rightarrow u_{\mu }+\varepsilon r_{1}^{\nu }\partial u_{\mu }/\partial x^{\prime \nu },$ in Eqs. (24), (25), so that the definition of $u_{\mu }$is shifted to that point.

Now it is straightforward to eliminate from ${\delta G}^{\prime }$ terms oscillating in $\phi$ by properly defining displacements $\tens{r}_{{\rm s}}.$ We carry out the calculations for each order in $\varepsilon$ and $\lambda$.

4.1 Order ${\varepsilon }^{{1}}$

Initially we drop all terms in Eq. (24), which are of order  ${\varepsilon }^{{1}}$or higher, while retaining contributions of order $\lambda$. The conditions for $\phi$ to be ignorable for the gauge-modified functional ${\ \delta G}^{\prime \prime }={\delta G}^{\prime }-\varepsilon {\rm d}R$, where R is an arbitrary Gauge function look like

  

$$\tilde{u}_{\mu }+\lambda q\tilde{a}_{\mu }-qr_{1}^{\nu }F_{\mu \nu } =0;$$ (27)

$$\left\{ (u_{\mu }+\lambda qa_{\mu }-\frac{1}{2}qr_{1}^{\nu }F_{\m... ...mu }}{\partial \phi }\right\} ^{\sim } =\frac{ \partial R}{\partial \phi }\cdot$$ (28)

Here $\tilde{y}$ denotes the oscillating part of y, namely $\tilde{y}=y- \bar{y}$, where $\bar{y}=\langle y\rangle _{\phi }$ is the $\phi$-independent part of y. By a proper choice of the gauge-function R one can always satisfy Eq. (28). This means that

$$\tilde{u}_{\mu }+\lambda q\tilde{a}_{\mu }-qr_{1}^{\nu }F_{\mu \nu }=0$$ (29)

yields the only essential requirements. The physical meaning of Eq.( In the text 29 ) is the relationship between the first-order gyro-radius $\tens{{r}_{1}}$and the velocity $\tens{{\tilde{u}}}$ along the Larmor orbit. The time-like component describes changes in energy due to displacements of charge in the electric field.

If the above requirements (27) and (28) are satisfied and $\phi$ is ignorable, the phase-space variational principle in our approximation can be expressed as $\delta S^{\prime \prime }=0$ with

$$S^{\prime \prime }=\int \delta G^{\prime \prime }=\int \left\... ...}{\partial \phi } \right\rangle _{\phi }{\rm d}\phi \right\} .$$ (30)

Equation (29) can be formally solved for $r_{1}^{\nu }$ to yield

$$r_{1}^{\nu }=\frac{w}{qH}(l^{\prime \prime \nu }\cos \phi -l^{\prime \nu }\sin \phi )+\lambda D^{\nu \mu }\tilde{a}_{\mu },$$ (31)

where $D^{\nu \mu }=(F_{\mu \nu })^{-1}$ is given by Eq. (11), provided that $F_{\mu \nu }$ is not degenerate. As a result, the $\phi$-independent functional $S^{\prime \prime }$ becomes

$$S^{\prime \prime }=\int \left\{ \left( \frac{q}{\varepsilon }... ...\right) {\rm d}x^{\prime \mu }+\hat{\mu}{\rm d}\phi \right\} ,$$ (32)

where $\hat{\mu}$ is the wave-field-modified magnetic moment, accurate to order ${\varepsilon }^{{1}}$,

$$\hat{\mu}=\frac{w^{2}}{2qH}-\frac{\lambda w}{H}\left\langle (... ...ilde{a}_{\mu }}{\partial \phi }\tilde{a}_{\nu }\right\rangle ,$$ (33)

and $u_{o}=\sqrt{1+w^{2}+u{_{\parallel }}^{2}}.$

If the Faraday tensor is degenerate but the magnetic field is non-zero, i.e., E=0, $\ H\neq 0$, a solution of Eq. (29) exists for a special choice of the wave potential satisfying $c^{\nu \mu }\tilde{a}_{\mu }=0,$i.e., when the $\left( l\tau \right)$ projection of $\tilde{a}$ is zero. Then the above expressions are valid with the replacement $D^{\nu \mu }\rightarrow -b^{\nu \mu }/H$. For other choices of the wave potential there is no solution, meaning that we need to modify our representation of the gyrokinetic transformation, Eqs. (15) and (18). This case will be discussed elsewhere.

The evaluation of the gyrophase averages involving the wave field $a_{\mu }$is, of course, the difficult part. The reason for this is the necessity of transforming the given function of space-time coordinates $a_{\mu }(x^{\nu })$ into a function of new variables, while the transformation rule (15), (31) is itself dependent on $a_{\mu }$. A solution of the equation

$$a_{\mu }=a_{\mu }\left( x^{\prime \nu }+\varepsilon \left(\fr... ...hi \right)+\lambda D^{\nu \mu }\tilde{ a}_{\mu }\right)\right)$$ (34)

can only be obtained by expanding in powers of $\lambda$.

4.2 Order ${\varepsilon^{1}\protect\lambda^{0}}$

The solution of the problem in the limit $\lambda =0$ has already been published by Beklemishev & Tessarotto (1999). By varying the functional (32) with respect to $\phi$, we find ${\rm d}\hat{\mu}=0,$ i.e.,

$$\hat{\mu}\equiv \frac{w^{2}}{2qH},$$ (35)

the magnetic moment calculated in the reference frame where $\vec{E}$ and $\vec{B}$ are parallel is the first integral of motion in our approximation.

Now it is possible to rewrite the functional (32) using Eq. (19) and $2qH\hat{\mu}$ in place of w²:

$$S^{\prime \prime }=\int \left\{ \left( qA_{\mu }^{\prime }+u{... ...\right) {\rm d}x^{\prime \mu }+\hat{\mu}{\rm d}\phi \right\} ,$$ (36)

with $u_{o}=\sqrt{1+2qH\hat{\mu}+u{_{\parallel }}^{2}}$.

Using the $\phi$-independent functional (36) in the variational principle, $\delta S^{\prime \prime }=0,$ defines the particle trajectory in terms of the new gyrokinetic variables $(x^{\prime \mu },\hat{\mu},{\ u_{\parallel }},\phi )$. This set is non-canonical, but further transformations of the variables (not involving $\phi$) also lead to $\phi$-independent functionals, and thus can be used to transform to canonical variables, if necessary.

4.3 Order ${\varepsilon ^{2}}$

We return to the approximation (24) with $\tens{{r}_{1}}$ determined by the lower-order requirement (29). This yields

$${\delta G}^{\prime }=\left( \frac{q}{\varepsilon }A_{\mu }^{\... ...^{\mu }})_\phi +(\Phi+\hat{\mu}) {\rm d}\phi +o(\varepsilon ).$$

Choosing the independent velocity variables to be $\phi ,w$ and ${ \ \ u_{\parallel }}$, and noting that $\tens{{r}_{1}}$ is independent of ${ u_{\parallel },}$ we can write

$$({\rm d}r_{1}^{\mu })_\phi=\frac{\partial r_{1}^{\mu }}{\part... ...}\varepsilon \frac{\partial r_{1}^{\mu }}{\partial w}{\rm d}w.$$ (37)

Inserting this into ${\delta G}^{\prime }$ leads to the following conditions for the gauge-modified functional, ${\delta G}^{\prime \prime }={\delta G} ^{\prime }-\varepsilon ^{2}{\rm d}R^\prime,$ with arbitrary $R^\prime$, to be independent of $\phi$:

  

$$\left\{ \left( \bar{u}_{\mu }+\lambda q\bar{a}_{\mu }+\frac{1}{2}... ...}_{\mu })\right) \frac{\partial r_{1}^{\mu }}{ \partial w}\right\} ^{\sim } =0,$$ (38)

$$\left\{ \left( qr_{2}^{\nu }F{_{\mu \nu }+}\frac{q}{2}r_{1}^{\nu ... ... \frac{\partial r_{1}^{\nu }}{ \partial x^{\prime {\mu }}}\right\} ^{\sim } =0,$$ (39)

$$\tilde{\Phi}=\frac{\partial R^\prime}{\partial \phi },$$ (40)

where $\Phi$ is given by Eq. (25). The last equation can always be satisfied by a suitable choice of $R^\prime$, Eqs. (39) define the second-order displacement $\tens{{r} _{2}{,}}$ while Eq. (38) can be regarded as a solvability condition. Using the solution for $\tens{{r}_{1}}$ it can also be rewritten as

$$\left\{ \left( \bar{u}_{\mu }+\lambda q\bar{a}_{\mu }+\frac{1... ...tial \tilde{a}_{\nu }}{\partial w}\right) \right\} ^{\sim }=0.$$ (41)

4.4 Order ${\varepsilon^{2}\protect\lambda^{0}}$

With $\lambda =0$, and in view of the antisymmetric properties of $D^{\mu \nu },$ Eq. (41) becomes

$$\left\{ \left( {u_{\parallel }}l_{\mu }+u_{o}\tau _{\mu }\rig... ...\phi +l_{\nu }^{\prime \prime }\sin \phi )\right\} ^{\sim }=0,$$

which is satisfied identically if the $(\tens{{l}^{\prime}{,l}^{\prime\prime}})$-plane is the invariant plane of $D^{\mu \nu },$ since it is orthogonal to $(\tens{l,\tau})$. Thus the field-related choice of the basis vectors (see Eqs. (7), (8)) is sufficient for the solvability of the second-order problem.

In the current approximation,

$$r_{1}^{\mu }=\frac{w}{qH}(l^{\prime \prime \mu }\cos \phi -l^{\prime \mu }\sin \phi ).$$

Then, Eq. (39) can be rewritten as

$$\left\{ q\left( r_{2}^{\nu }F{_{\mu \nu }+}\frac{1}{2}r_{1}^{... ...artial u_{\mu }}{\partial x^{\prime \nu }}\right\} ^{\sim }=0,$$ (42)

while the $\phi$-independent functional becomes

$${\delta G}^{\prime \prime }=\left( \frac{q}{\varepsilon }A_{\... ...\left\langle \Phi \right\rangle {\rm d} \phi +o(\varepsilon ).$$ (43)

The average $\left\langle \Phi \right\rangle$ is zero, and we get

$$\hat{\mu}=w^{2}/2qH$$

as in the lower order approximation, i.e., there are no second-order (in  $\varepsilon$) corrections to the magnetic moment. Thus, the second order modifications to the Lagrangian can be described as inhomogeneity contributions to the effective electromagnetic potential, while the form of the adiabatic invariant remains unchanged.

Using the explicit expression for $\tens{{r}_{1}}$ it is possible to evaluate averages, so that

$${\delta G}^{\prime \prime }=\left( \frac{q}{\varepsilon }A_{\... ...{\mu }\right) {\rm d} {x^{\prime \mu }}+\hat{\mu}{\rm d}\phi ,$$ (44)

where

$$\chi _{\mu }=l_{\nu }^{\prime }\frac{\partial l^{\prime \prim... ... }^{\prime \prime }}{ \partial {x^{\prime \nu }}}\right) \cdot$$ (45)

Describes the inhomogeneity of the electromagnetic field.

As one can see, it is not necessary to actually solve for $\tens{{r}_{2}}$to find the transformed Lagrangian. It is only necessary to be sure that the solution exists, which is true for the "field-related'' choice of the basis.

4.5 Order ${\varepsilon^{1}\protect\lambda^{1}}$

We proceed by considering the linear approximation in the wave-field amplitude without high order curvature effects. As noted above, the main problem here is the transformation of the highly-local wave field $a_{\mu }( \tens{x})$, where $\tens{x}$ is the particle position, to the guiding-center coordinates $(x^{\prime \alpha },\phi ,\hat{\mu},{\ u_{\parallel }})\$. The Taylor-expansion procedure used above for the transformation of the equilibrium field fails if the wavelength is sufficiently short. To avoid this problem, Eq. (34) is solved using Fourier analysis and expansion in powers of $\lambda$.

Assume that the wave field is given in terms of its Fourier components $a_{\mu}(\tens{k}),$ then Eq. (34) becomes

$$a_{\mu}(\tens{x})=\frac{1}{\left( 2\pi\right) ^{2}}\int{\rm d... ...\sin \phi)+\lambda D^{\nu\mu}\tilde{a}_{\mu}\right) \right\} ,$$ (46)

where only the second exponential factor depends on $\phi$.

Here the ordering of the wave-vector $\tens{k}$ is such that $\varepsilon k\sim O(1),$ as usual (Littlejohn 1984) , except that the time-like component of $\tens{k}$ is proportional to the wave frequency, and the ordering for it is equivalent to $\omega /\Omega _{\rm c}\sim O(1),$ where $\Omega _{\rm c}$ is the cyclotron frequency, and not $\omega /\Omega _{\rm c}\sim O(\varepsilon )$ as in conventional approaches; besides, the above ordering is also applicable to $k_{\Vert },$ i.e., $k_{\Vert }\rho _{ {\rm L}}\sim O(1)$ is assumed in place of the usual ordering $k_{\Vert }\rho _{{\rm L}}\sim O(\varepsilon )$. Therefore, our approach in this respect has a wider applicability range. This is inevitable, for the reference frame may move with relativistic velocity across the short-wavelength wave field, so that fast oscillations in time may result anyway. However, the particle, drifting across the wave field, will experience a quite definite oscillation frequency as well as the wavelength, which is taken into account by our definition of $\xi =k_{\bot }\rho _{{\rm L}}$, introduced below.

We can calculate the required averages in Eqs. (32) and (33) by using the expansion of Eq. (46) into a Fourier series in $\phi$. We use the following identity (Abramovits & Stegun 1972)

$$\exp \left[ {\rm i}\zeta \sin \psi \right] =\sum_{n=-\infty }... ... }J_{n}\left( \zeta \right) \exp \left( {\rm i}n\psi \right) ,$$ (47)

to get

$$a_{\mu }(\tens{x})=\frac{1}{\left( 2\pi \right) ^{2}}\int {\r... ...mbda \varepsilon k_{\nu }D^{\nu \mu }\tilde{a}_{\mu }\right) .$$ (48)

Here

$$\xi =k_{\bot }\rho _{{\rm L}}=(\varepsilon w/qH)\sqrt{\left( ... ...\nu }\right) ^{2}+\left( k_{\nu }l^{\prime \nu }\right) ^{2}},$$ (49)

where $k_{\bot }=\sqrt{\left( k_{\nu }l^{\prime \nu }\right) ^{2}+\left( k_{\nu }l^{\prime \prime \nu }\right) ^{2}}$ is the length of $\tens{k} _{\bot }$ , which is the projection of $\tens{k}$ on the $(\tens{{l}^{\prime}{,l}^{\prime\prime}})$-plane; $\phi _{0}=\arctan \left( k_{\nu }l^{\prime \prime \nu }/k_{\nu }l^{\prime \nu }\right)$ is the angle of $\tens{k} _{\bot }$ with respect to the $(\tens{{l}^{\prime}{,l}^{\prime\prime}})$-basis, $\rho _{{\rm L}}=\varepsilon w/qH$ is the Larmor radius.

To the order ${\varepsilon ^{1}\lambda ^{1}}$, $\lambda$ in the expression (48) can be set to 0, so that its zero-order solution is just

$$a_{\mu}^{(0)}(\tens{x})=\frac{1}{\left( 2\pi\right) ^{2}}\int... ...\left( \xi\right) \exp\left( {\rm i}n(\phi_{0} -\phi)\right) .$$ (50)

Note that neglecting the last exponent in Eq. (48) or solving by expansion in powers of $\lambda$ is possible only if the particle displacement from the equilibrium orbit due to the wave field is much smaller than the wavelength.

Now the averages, present in expressions (32), (33), become

 

$$\bar{a}_{\mu}^{(0)} =\frac{1}{4\pi^{2}}\int{\rm d}^{4}k\ a_{\mu }( \tens{k}){\rm e}^{{\rm i}\tens{kx}^{\prime}}J_{0}\left( \xi\right) ,$$ (51)

$$\left\langle (l^{\prime\prime\nu}\sin\phi+l^{\prime\nu}\cos\phi)\... ...}\left( \xi\right) (l^{\prime\nu}\sin\phi_{0}-l^{\prime\prime\nu}\cos\phi_{0}).$$ (52)

So that to the order ${\varepsilon ^{1}\lambda ^{1}}$

$$S^{\prime\prime}=\int\left\{ \left( \frac{q}{\varepsilon}A_{\... ...u}\right) {\rm d}x^{\prime\mu}+\hat{\mu}{\rm d} \phi\right\} ,$$ (53)

$$\hat{\mu}=\frac{w^{2}}{2qH}\left( 1-\frac{\lambda \varepsilon... ...b^{\mu \nu }\frac{J_{1}\left( \xi \right) }{\xi }\right) \cdot$$ (54)

Note that this definition of $\hat{\mu}$ is invariant under gauge transformations of $a_{\nu }$, and the adiabatic invariant can also be rewritten in terms of Fourier components of the electromagnetic field tensor of the wave $f_{\mu \nu }=\partial a_{\nu }/\partial x^{\mu }-\partial a_{\mu }/\partial x^{\nu }$:

$$\hat{\mu}=\frac{w^{2}}{2qH}\left( 1-\frac{\lambda \varepsilon... ...i \right) }{\xi }{\rm e}^{{\rm i}\tens{kx}^{\prime }}\right) ,$$ (55)

where $f_{\mu \nu }(\tens{k})l^{\prime \nu }l^{\prime \prime \mu }=f_{\mu \nu }b^{\mu \nu }/2$ is essentially the parallel component of the wave field. This effect is usually neglected in the standard gyrokinetic ordering, so that the first order contribution to the magnetic moment vanishes.

4.6 Order ${\varepsilon ^{1}\protect\lambda ^{2}}$ (nonlinear corrections)

The second-order nonlinear corrections are responsible, in particular, for three-wave interactions. Calculating them for a general Fourier-decomposed wave field would allow subsequent construction of a weak-turbulence theory. For other effects, such as the estimate of the Miller force of the high-frequency field on a particle, just the eikonal representation of the wave is sufficient. Below we adopt the first, more complete description.

In this approximation we can use the $\phi$-independent variational principle (32) with properly defined parameters. Indeed, besides $\bar{a}_{\mu }^{(1)},$which is added to the effective vector potential, $\bar{a}_{\mu }=\bar{a}_{\mu }^{(0)}+\bar{a}_{\mu }^{(1)}$, there is one other significant contribution to the Lagrangian from the wave terms. It enters into the definition of the adiabatic invariant $\hat{\mu}$ via Eq. (33), so that $\hat{\mu}=\hat{\mu}^{(0)}+\hat{\mu}^{(1)}+\hat{\mu}^{(2)}$, where $\hat{\mu}^{(0)}+\hat{\mu}^{(1)}$is given by Eq. (54).

The second order contribution to the adiabatic invariant is given by

$$\hat{\mu}^{(2)}=-\frac{\lambda w}{H}\left\langle (l^{\prime \... ... }^{(0)}}{\partial \phi }\tilde{a}_{\nu }^{(0)}\right\rangle .$$ (56)

The second term in this expression can be found using the zero-order solution (50), while the first involves the first-order correction, $\tilde{a}_{\nu }^{(1)}(\tens{x}).$

Expanding the last exponent in Eq. (48) in powers of $\lambda$ we get

$$% a_{\mu}^{(1)}(\tens{x})={\rm i}\lambda\varepsilon\int{\rm d... ...eta}\tilde{a}_{\eta}( \tens{x})\exp\left(-{\rm i}n\phi\right),$$ (57)

where

$$a_{\mu n}(\tens{x}^{\prime},\tens{k})=\frac{1}{4\pi^{2}}\ a_{... ..._{n}\left( \xi\right)\exp\left({\rm i}n\phi_0(\tens{k})\right)$$ (58)

by definition. The average of Eq. (57) then looks like

$$\bar{a}_{\mu}^{(1)} ={\rm i}\lambda\varepsilon\int{\rm d}^{4}... ...},\tens{k})a_{\eta ,-n}( \tens{x}^{\prime},\tens{k}^{\prime}).$$ (59)

The sum of the Bessel-function products in the above expression can be simplified using the Graf addition theorem for Bessel functions (Abramovits & Stegun 1972). This theorem can be cast in the form

$$% \sum_{n=-\infty }^{+\infty }J_{n}\left( \xi \right) J_{-l-n... ...( \xi ^{\prime \prime }\right) {\rm e}^{{\rm i}l(\pi -\chi )},$$ (60)

where l is arbitrary, while $\xi ,\xi ^{\prime },\xi ^{\prime \prime }$are proportional to the sides of a triangle, $\alpha$ being the angle between $\xi$ and $\xi ^{\prime },$ while $\chi$ is the angle opposite to $\xi$. These conditions are satisfied in our case (due to the interpretation of $\xi$ as proportional to the length of $\tens{k} _{\bot }$, (49)) with

$$\left\{ \begin{array}{lc} \xi ^{\prime }=\xi (\tens{k}^{\prim... ... _{0}^{\prime \prime }-\phi _{0}^{\prime }, \end{array}\right.$$ (61)

where $\phi_{0}^{\prime}=\phi_{0}(\tens{k}^{\prime})$, $\phi _{0}^{\prime \prime }=\phi _{0}(\tens{k+k}^{\prime }).$ The triangle is formed by four-vectors $\tens{k}_{\bot },\tens{k}_{\bot }^{\prime },\tens{k}_{\bot }+\tens{k}_{\bot }^{\prime }$ in the $(\tens{{l}^{\prime}{,l}^{\prime\prime}})$-plane.

Using the definition of $a_{\mu n},$ and taking into account the identity (60) with l=0, we reduce the first-order correction to $\bar {a} _{\mu}$ to the following expression

$$% \bar{a}_{\mu}^{(1)}=\frac{{\rm i}\lambda\varepsilon}{\left(... ...{0}\left( \xi\right) J_{0}\left( \xi^{\prime}\right) \right] .$$ (62)

Then, seeking the oscillating part of the first-order correction as a series

$$\tilde{a}_{\mu }^{(1)}(\tens{x})=\sum_{l\neq 0}a_{\mu l}^{(1)}(\tens{x} ^{\prime })\exp \left( {\rm i}l\phi \right) ,$$

and using Eq. (57), we find that the coefficients $a_{\mu l}^{(1)}$satisfy

$$% a_{\mu l}^{(1)}(\tens{x}^{\prime })={\rm i}\lambda \varepsi... ...{k})\ a_{\eta (-l-n)}( \tens{x}^{\prime },\tens{k}^{\prime }).$$ (63)

The sum under the integral can be found using the identity (60) as

$$\sum_{n+l\neq 0}a_{\mu n}(\tens{x}^{\prime },\tens{k})\ a_{\eta (-l-n)}( \tens{x}^{\prime },\tens{k}^{\prime }) \frac{a_{\mu }(\tens{k})a_{\eta }( \tens{k}^{\prime })}{\left( 2\... ...t) \exp \{{\rm i}n(\phi _{0}-\phi _{0}^{\prime })-{\rm i}l\phi _{0}^{\prime }\}$$ =

$$\frac{a_{\mu }(\tens{k})a_{\eta }(\tens{k}^{\prime })}{\left( 2\p... ...me }\right) J_{l}\left( \xi \right) {\rm e}^{{\rm i}l(\pi -\phi _{0})}\right] .$$ =

Fortunately, knowledge of the first harmonics $\tilde{a}_{\mu \pm 1}^{(1)}$is sufficient for the evaluation of $\hat{\mu}^{(2)}$. Indeed, the necessary average can be calculated as

$$\left\langle (l^{\prime \prime \mu }\sin \phi +l^{\prime \mu }\cos \phi ) \tilde{a}_{\mu }^{(1)}\right\rangle \frac{1}{2}\left[ l^{\prime \mu }\left( a_{\mu 1}^{(1)}+a_{\mu (-... ...}l^{\prime \prime \mu }\left( a_{\mu 1}^{(1)}-a_{\mu (-1)}^{(1)}\right) \right]$$ =

$$-\frac{\lambda \varepsilon }{\left( 2\pi \right) ^{4}}\iint {\rm ... ...}^{\prime \prime }-l^{\prime \prime \mu }\cos \phi _{0}^{\prime \prime }\right)$$ =

$$-J_{0}\left( \xi ^{\prime }\right) J_{1}\left( \xi \right) \left( l^{\prime \mu }\sin \phi _{0}-l^{\prime \prime \mu }\cos \phi _{0}\right) \big] .$$

After substituting the definitions of $\phi _{0},\phi _{0}^{\prime \prime }$ the above expression becomes

$$% \left\langle (l^{\prime \prime \mu }\sin \phi +l^{\prime \m... ...ight) \frac{J_{1}\left( \xi \right) }{\xi }k_{\zeta }\right] ,$$ (64)

where $b^{\mu \nu }=l^{\prime \prime \mu }l^{\prime \nu }-l^{\prime \mu }l^{\prime \prime \nu }.$

Proceeding now to calculation of the second term in Eq. (56) we note that

$$% \left\langle \frac{\partial \tilde{a}_{\mu }^{(0)}}{\partia... ...me }\right) {\rm e}^{{\rm i}n(\phi _{0}-\phi _{0}^{\prime })}.$$ (65)

Using again the addition theorem for Bessel functions (60) with l=0and Eq. (61) we get

$$% -\sum_{n\neq 0}{\rm i}nJ_{n}\left( \xi \right) J_{-n}\left(... ...}{\partial \alpha } J_{0}\left( \xi ^{\prime \prime }\right) ,$$ (66)

with

$$\xi ^{\prime \prime }=\sqrt{\xi ^{2}+\xi ^{\prime 2}-2\xi \xi ^{\prime }\cos \alpha }.$$ (67)

Then

$$\frac{\partial }{\partial \alpha }J_{0}\left( \xi ^{\prime \p... ...prime \zeta }-l^{\prime \zeta }l^{\prime \prime \eta }\right).$$ (68)

Finally, incorporating the above averages into Eq. (56), the second-order correction to the adiabatic invariant (magnetic moment) becomes

$$% \hat{\mu}^{(2)}=\frac{w^{2}}{qH^{2}}\frac{\lambda ^{2}\vare... ...{\rm i}(\tens{k+k}^{\prime }) \tens{x}^{\prime }}R^{\mu \nu },$$ (69)

where

$$R^{\mu \nu }=b^{\zeta \mu }\left[ \frac{J_{1}\left( \xi ^{\pr... ...{\prime \prime }}k_{\eta }k_{\zeta }^{\prime }b^{\zeta \eta }.$$ (70)

The above results, Eqs. (62), (69), look complicated because of the presence of general three-wave interactions. Indeed, they show that contributions to the particle parameters with the wave vector $\tens{k}_{\bot }^{\prime \prime }= \tens{k}_{\bot }+\tens{k}_{\bot }^{\prime }$ come from all pairs of electromagnetic perturbations with the same sum of the wave vector $\tens{k} _{\bot },\tens{k}_{\bot }^{\prime }$. If we are interested in the average Miller force only, it is sufficient to retain only the zero-sum pairs, $\tens{k}_{\bot }=-\tens{k}_{\bot }^{\prime }$, so that $\;\xi ^{\prime \prime }=0$, which simplifies the expressions significantly.

4.7 Gauge invariance

As one can see, the above results are not gauge invariant with respect to gauge transformations of the wave potential $a_{\mu }$. This problem can be resolved using gauge-invariance of the variational principle (4), (14). Indeed, any transformation of the potential

$$a_{\mu }\rightarrow a_{\mu }+\frac{\partial f(\tens{x})}{\partial {x^{\mu }}}$$ (71)

can be compensated by

$${Q_{\mu }{\rm d}x^{\mu }}\rightarrow {Q_{\mu }{\rm d}x^{\mu }-}{\rm \ d}f,$$ (72)

so that the functional remains the same. Conversely, any transformation of type (72) leads to an equivalent variational principle, while the difference can be hidden in $a_{\mu }$.

The above program of calculating the gyrokinetic transformation and the $\phi$-independent description of motion can be carried out for an arbitrary gauge function f. This means that expressions obtained in the above subsections should also be valid for

$$% a_{\mu }(\tens{k})\rightarrow a_{\mu }(\tens{k})+{\rm i}k_{\mu }f(\tens{k})$$ (73)

with an arbitrary function $f(\tens{k})$. Thus, any gauge transformations of the wave potential $a_{\mu }$ lead to different, but still valid expressions.

4.8 Summary

Finally, we summarize the results of this section by presenting the transformed variational principle valid through the second order in $\lambda$ and second order in $\varepsilon ,$ i.e., with terms of the order $\varepsilon \lambda ,\varepsilon \lambda ^{2}$ and ${\varepsilon ^{2}\lambda ^{0}}$ retained: $\delta S=0$ yields the particle phase-space trajectory with

$$S=\int \left( qA_{\mu }^{\prime }+u{_{\parallel }}l_{\mu }+\l... ...\mu }\right) {\rm d}{x^{\prime \mu } }+\hat{\mu}{\rm\ d}\phi ,$$ (74)

where $({x^{\prime \mu },}u{_{\parallel },}\hat{\mu},\phi )$ or $({x^{\prime \mu },}u{_{\parallel },}w,\phi )$ are the new gyrokinetic variables with

$$\hat{\mu}=\frac{w^{2}}{2qH}\left( 1-\frac{1}{2\pi ^{2}}\int {... ...{\rm e}^{{\rm i}\tens{kx}^{\prime }}\right) +\hat{\mu} ^{(2)},$$

where q=q_a/m_ac² is the signed charge-to-mass ratio,

$$% \bar{a}_{\mu }=\frac{1}{\left( 2\pi \right) ^{2}}\int {\rm ... ...{kx}^{\prime }}J_{0}\left( \xi \right) +\bar{a} _{\mu }^{(1)},$$ (75)

$\xi =k_{\bot }\rho _{{\rm L}}=\left( 2\hat{\mu}\left[ \left( k_{\nu }l^{\prime ... ...^{2}+\left( k_{\nu }l^{\prime \prime \nu }\right) ^{2}\right] /qH\right) ^{1/2}$, and the second-order corrections $\hat{\mu}^{(2)}$ and $\bar{a}_{\mu }^{(1)}$ are given by Eqs. (69) and (62), respectively. The orthogonal basis $( \tens{\tau,l,l^{\prime },l^{\prime\prime}})$ is chosen in such a way that $(\tens{{l}^{\prime}{,l}^{\prime\prime}})$ coincides with the space-like invariant plane of the antisymmetric tensor of the electromagnetic field $F{_{\mu \nu }}=\partial A_{\nu }^{\prime }/\partial x^{\prime {\mu }}-\partial A_{\mu }^{\prime }/\partial x^{\prime {\nu }}$, with H being the corresponding eigenvalue.

5 The kinetic equation and equations of motion

The single-particle distribution function can be written in general relativity either in the eight-dimensional phase space $\Phi (x^{\mu },u_{\nu }),$ $\ \mu ,\nu =0,...,3,$ or in the seven-dimensional phase space $f(x^{\mu },u_{\nu }),$ where only 3 components of the 4-velocity $u_{\nu }$ are independent, so that

$$\Phi (x^{\mu },u_{\nu })=f(x^{\mu },u_{\nu })\delta \left(\sqrt{u_{\zeta }u^{\zeta }}-1\right)\theta \left(u^{0}\right).$$

The $\delta$-function here reflects the fact that $u_{\zeta }u^{\zeta }=1$is the first integral of motion in the case of the eight-dimensional representation. However, the collisionless kinetic equation in both cases retains the same form, namely

$$u^{\mu }\frac{\partial f}{\partial x^{\mu }}+\left( \frac{{\r... ...u }}{{\rm d}s} \right) \frac{\partial f}{\partial u_{\nu }}=0,$$ (76)

although in the 7-dimensional case $\nu =1,2,3$ only, while u⁰ is the dependent variable. Here $\left( {\rm d}u_{\nu }/{\rm d}s\right)$ is a function of the independent variables $(x^{\mu },u_{\nu })$ found as the right-hand side of the single-particle dynamics equations. The kinetic equation can be multiplied by ${\rm d}s$ and thus represented in the parameterization-independent form as

$$\frac{\partial f}{\partial x^{\mu }}{\rm d}x^{\mu }+\frac{\partial f}{ \partial u_{\nu }}{\rm d}u_{\nu }=0,$$

where the differentials are tangent to the particle orbit.

Because of the general properties of variable transformations it is obvious that any non-degenerate transformation of the phase-space variables $(x^{\mu },u_{\nu })\rightarrow (y^{i})$ will lead to the same form of the kinetic equation

$$\frac{\partial f}{\partial y^{i}}{\rm d}y^{i}=0,$$

where the differentials are tangent to the particle orbit. In particular, this property is useful for the transformation to the gyrokinetic variables.

Let $(y^{i})\equiv (x^{\prime \alpha },\phi ,\hat{\mu},{u_{\parallel }})$, then the kinetic equation becomes

$$\frac{\partial f}{\partial x^{\prime \mu }}{\rm d}x^{\prime \... ... d}\hat{\mu}+\frac{\partial f}{\partial \phi } {\rm d}\phi =0.$$ (77)

By definition of the gyrokinetic variables the dynamic equations should be independent of $\phi$, i.e., expressions for $\left( {\rm d}x^{\prime \mu }/{\rm d}\phi \right) ,\left( {\rm d}u_{\parallel }/{\rm d}\phi \right) ,\left( {\rm d}\hat{\mu} /{\rm d}\phi \right)$ are independent of $\phi$, while $\partial f/\partial \phi$ is periodic in $\phi$. It follows that $\partial f/\partial \phi =0,$and, if $\hat{\mu}$ is the integral of motion, ${\rm d}\hat{\mu}=0$, we get the kinetic equation expressed in the gyrokinetic variables as

$$\frac{\partial f}{\partial x^{\prime \mu }}{\rm d}x^{\prime \... ...{ \partial f}{\partial u_{\parallel }}{\rm d}u_{\parallel }=0,$$ (78)

which we shall call the gyrokinetic equation. Here the coefficients ${\rm d}x^{\prime \mu }$ and ${\rm d}u_{\parallel }$ must be determined from the equations of motion in the gyrokinetic variables.

While the above derivation seems simple, there is a trick. Even if we prove that there is no $\phi$-dependence in all orders, this does not mean that there is no such dependence at all. In fact, this is what happens close to the cyclotron resonances - the exponentially small (and thus having zero expansion coefficients) resonant factor destroys the conservation of the magnetic moment, rendering the gyrokinetic theory invalid. It should always be kept in mind that the approximation works well for particles far from the low-order cyclotron resonances.

Another question is: why is it so important to have parameterization-independent equations of motion? Would it be possible to use, say, Fradkin's drift equations with the proper-time parameterization to construct the gyrokinetic equation in an alternative way? The answer follows from the transformation of the Liouville Eq. (76) with the transformation of the phase-space variables. Since its meaning is the conservation of world lines, the form is retained for canonical-like transformations (with a constant Jacobian). Otherwise, for an arbitrary set of non-canonical drift equations of motion, one has to prove that the Jacobian is at least independent of the gyrophase, in order to exclude this variable from the gyrokinetic equation. In other terms, kinetic Eq. (77) (with coefficients ${\rm d}x^{\prime \mu }/{\rm d}s,{\rm \ \ d}u_{\parallel }/{\rm d}s$independent of $\phi$) cannot be transformed into the gyrokinetic equation, unless ${\rm d} \phi /{\rm d}s$is also independent of $\phi$. This requirement is not fulfilled for Fradkin's equations, since the proper time s depends on the gyrophase with periodic corrections due to changes of velocity along the orbit.

5.1 The equations of motion

The equations of motion, or the relationships between differentials tangent to the particle orbit, can be obtained as Euler equations of the transformed variational principle (74). Let us find the first variation of Sassuming $(x^{\prime \alpha },\phi ,\hat{\mu},{\ u_{\parallel }})\$ to be independent:

$$\delta S \int \left( -qF_{\mu \nu }+d_{\nu \mu }\right) \delta x^{\prime \... ...}\tau _{\mu }}{u_{0}}\right) \delta u{\ _{\parallel }}{\rm d}{x^{\prime \mu } }$$ =

$$-\left( \tau _{\mu }\frac{\partial u_{0}}{\partial \hat{\mu}}+q\f... ...lta x^{\prime \mu }}+\delta \hat{\mu}{\rm d}\phi -{\rm d}\hat{\mu}\delta \phi .$$

Here we introduced a new tensor notation

$$d_{\nu \mu }\equiv {\rm Rot}_{\nu }\left[ u{_{\parallel }}l_{... ...mu }+q\bar{a}_{\mu }+\frac{1}{2}\hat{\mu}\chi _{\mu }\right] ,$$ (79)

where the operator ${\rm Rot}_{\nu }$ is defined by ${\rm Rot}_{\nu }f_{\mu }\equiv \partial f_{\mu }/\partial x^{\prime \nu }-\partial f_{\nu }/\partial x^{\prime \mu }.$

Extrema of the functional are achieved for world lines where $\delta S=0$for all variations of independent variables. This yields the Euler equations as

   

$${\rm d}\hat{\mu} =0,$$ (80)

$$\left( \tau _{\mu }\frac{\partial u_{0}}{\partial \hat{\mu}}+q\fr... ...t{\mu}}+\frac{\chi _{\mu }}{2}\right) {\rm d}{ x^{\prime \mu }+}{\rm d}\phi =0,$$ (81)

$$\left( l_{\mu }+\frac{u{_{\parallel }}\tau _{\mu }}{u_{0}}\right) {\rm d}{ x^{\prime \mu }} =0,$$ (82)

$$\left( -qF_{\mu \nu }+d_{\nu \mu }\right) {\rm d}{x^{\prime \mu }... ...\frac{u{_{\parallel }}\tau _{\nu }}{u_{0}}\right) {\rm d}u{\ _{\parallel }} =0.$$ (83)

The first equation here confirms that $\hat{\mu}$ is an integral of motion, the second equation determines the rate of rotation along the Larmor orbit and is not needed for our purposes here. By multiplying the 4-vector Eq. (83) by ${\rm d}{x^{\prime \nu }}$ and taking the sum in $\nu$ one can recover the third scalar Eq. (82) multiplied by ${\rm d}u_{\parallel }$, i.e., out of five equations for four unknown functions ${\rm d}{x^{\prime \mu }/}{\rm d}u{_{\parallel }}$only four are independent, as should be the case for solvability. This allows us to determine the non-trivial solutions (at least formally) and use them as coefficients of the gyrokinetic equation.

It is convenient to rewrite the 4-vector Eq. (83) through its projections on the orthogonal vectors of the basis $( \tens{\tau ,l,l^{\prime},l^{\prime\prime}}).$ Taking scalar products with $\left( \tens{l^{\prime},l^{\prime\prime}}\right)$ we get

  

$$\left( qHl_{\mu }^{\prime \prime }+l^{\prime \nu }d_{\nu \mu }\right) {\rm d}{x^{\prime \mu }} ={0,}$$ (84)

$$\left( -qHl_{\mu }^{\prime }+l^{\prime \prime \nu }d_{\nu \mu }\right) {\rm d}{x^{\prime \mu }} ={0.}$$ (85)

These two equations determine the drift velocity. Taking the scalar product with $\tens{l}$, we recover the following equation governing the acceleration parallel to the magnetic field

$${\rm d}u{_{\parallel }-}qE\tau _{\mu }{\rm d}{x^{\prime \mu }}+l^{\nu }d_{\nu \mu }{\rm d}{x^{\prime \mu }}=0.$$ (86)

Simplification occurred since the $(\tens{l,\tau})$ -plane for the field-related choice of the basis coincides with the second invariant plane of the electromagnetic field tensor, so that $F_{\mu \nu }l^{\nu }=E\tau _{\mu };\ \ F_{\mu \nu }\tau ^{\nu }=El_{\mu }$ (see Eqs. (7), (8)). Then, the last equation looks like

$$u{_{\parallel }}{\rm d}u{_{\parallel }+}u_{0}qEl_{\mu }{\rm d... ...mu }}-u_{0}\tau ^{\nu }d_{\nu \mu }{\rm d}{x^{\prime \mu }}=0,$$ (87)

which is the energy conservation law. However, as shown above, in its place one can use another, equivalent but much more compact Eq. (82). It describes the relationship between ${ \ \ u_{\parallel }}$ and the differentials of the guiding-center position

$$\left( u_{0}l_{\mu }+u{_{\parallel }}\tau _{\mu }\right) {\rm d}{ x^{\prime \mu }}=0.$$ (88)

It is possible to formally solve Eqs. (84)-(88) in terms of the basis vectors. To do this we introduce new notations,

$$v_{\mu }=u_{0}l_{\mu }+u{_{\parallel }}\tau _{\mu },$$ (89)

and $d_{pq}\equiv p^{\nu }d_{\nu \mu }q^{\mu },$ for arbitrary four-vectors $\tens{p,q}$. Note that $v_{\mu }\bar{u}^{\mu }=0,$ i.e., $\tens{v}$ is orthogonal to the particle world-velocity, and thus is orthogonal to the drift velocity too, as expressed by Eq. (88).

Searching for the solution in the form ${\rm d}{x^{\prime \mu }}=a\bar{u} ^{\mu }+bl^{\prime \mu }+cl^{\prime \prime \mu }$, and substituting it into Eqs. (84), (85), we get

$${\rm d}{x^{\prime \mu }=}\left[ \bar{u}^{\mu }+\frac{l^{\prim... ...-d_{l^{\prime }l^{\prime \prime }}}\right] {\rm d}s^{\prime },$$ (90)

$${\rm d}u{_{\parallel }=}\left[ \left( qE-d_{l\tau }\right) u_... ...-d_{l^{\prime }l^{\prime \prime }}}\right] {\rm d}s^{\prime }.$$ (91)

Note that here ${\rm d}s^{\prime }$ is just an arbitrary scalar function, which can be chosen as convenient. It means that the above equations determine just the direction of the world line of the gyrocenter in the phase space, but do not define its parameterization.

5.2 Symmetries and conservation laws

Consider a coordinate system chosen in such a way that the dependence of all functions on $x^{\prime \nu }$ (with some particular $\nu$) vanishes: $\partial /\partial {x^{\prime \nu }}=0$. Then

$$F_{\mu \nu }=\frac{\partial A_{\nu }^{\prime }}{\partial x^{\... ...u }+q\bar{a}_{\nu }+\frac{1}{2}\hat{\mu} \chi _{\nu }\right] ,$$

and the $\nu$-component of the Euler Eq. (83) becomes

$$\left( -q\frac{\partial A_{\nu }^{\prime }}{\partial x^{\prim... ...}+q\bar{a}_{\nu }+\frac{1}{2}\hat{\mu} \chi _{\nu }\right) =0.$$ (92)

Here the conservation of the magnetic moment, Eq. (80), and the absence of any dependence on $\phi$ are taken into account.

This means that

$$J=qA_{\nu}^{\prime}+u{_{\parallel}}l_{\nu}+u_{0}\tau_{\nu}+q\bar{a}_{\nu }+ \frac{1}{2}\hat{\mu}\chi_{\nu}$$ (93)

is the integral of motion corresponding to the symmetry with respect to $x^{\prime \nu }$.

6 Maxwell equations

In general, the kinetic description of plasmas involves a combination of the kinetic equation and Maxwell equations, which describe the evolution of the collective (Hartree-Fock) electromagnetic fields. The same is true for the gyrokinetic theory. However, since the gyrokinetic equation is written in specific gyrokinetic variables, and thus yields the distribution function in terms of these variables, one should have a special procedure of integration in velocity space to calculate the source terms of Maxwell equations. This procedure is presented below. Fortunately, the Jacobian of the gyrokinetic transformation in the present formulation is simple enough to allow integration in the gyrophase, so that expressions for the charge and current densities via the gyrokinetic distribution function are found explicitly.

The general form of Maxwell equations in the presence of an arbitrary gravitational field is well known (Landau & Lifshits 1975). The first pair of equations can be written as

$${\rm e}^{\varsigma \lambda \mu \nu }\frac{\partial F_{\mu \nu }}{\partial x^{\lambda }}=0,$$ (94)

and the second as

$$\frac{1}{\sqrt{-g}}\frac{\partial }{\partial x^{\nu }}\left( \sqrt{-g}F^{\mu \nu }\right) =\frac{4\pi }{c}j^{\mu },$$ (95)

where

$$j^{\mu }=c\sum_{\alpha }q_{\alpha }\int u^{\mu }f_{\alpha }\l... ...a \left( \sqrt{u^{\nu }u_{\nu }}-1\right) {\rm d} ^{4}\tens{u}$$ (96)

is the current density, expressed via the distribution function of particle species $\alpha$, and the signed particle charge  $q_{\alpha }$. The $\delta$-function under the integral makes it possible to do a partial integration, for example over u₀, and arrive at a more widely used form

$$\delta \left(\sqrt{u^{\nu }u_{\nu }}-1\right) {\rm d}^{4}\ten... ...ac{{\rm d}u_{1}{\rm d}u_{2}{\rm d}u_{3}}{u^{0}\sqrt{ -g}}\cdot$$

However, in the gyrokinetic transformation the four-velocity is expressed via Eqs. (18)-(20) as

$$u^{\mu }=w\left( l^{\prime \mu }\cos \phi +l^{\prime \prime \... ...\sin \phi \right) +{u_{\parallel }}l^{\mu }+u_{o}\tau ^{\mu },$$ (97)

so that ${\rm d}^{4}\tens{u}=w{\rm d}w{\rm d}\phi {\rm d} u_{\parallel }{\rm d}u_{o}$ (the sign is positive due to Eq. (6)), while the partial integration over ${\rm d}u_{o}$ leads to

$$\delta \left(\sqrt{u^{\nu }u_{\nu }}-1\right) {\rm d}^{4}\ten... ...row \frac{w{\rm d}w{\rm d}\phi {\rm d}u_{\parallel }}{u_{o}} ,$$ (98)

where $u_{o}=\sqrt{1+w^{2}+u{_{\parallel }}^{2}}.$ As a result, the expression for the components of the current density can be rewritten as

$$j^{\mu }=c\sum_{\alpha }q_{\alpha }\int \left( w\left( l^{\pr... ...\frac{w {\rm d}w{\rm d}\phi {\rm d}u_{\parallel }}{u_{o}}\cdot$$ (99)

In Eq. (99) the particle position $\tens{x},$ rather than its gyrocenter position $\tens{x}^{\prime },$ is kept constant while integrating over the particle velocity. This means, that if the distribution function $f_{\alpha }$ is given as a function of the gyrokinetic variables $f_{\alpha }=f_{\alpha }\left( x^{\prime \mu },w,u_{\parallel }\right)$, it is necessary to transform it to particle coordinates via the inverse transformation

$$x^{\mu }=x^{\prime \mu }+\sum_{i=1}\varepsilon ^{i}r_{i}^{\mu }$$ (100)

before integrating. Furthermore, here we assumed that the distribution function is defined in terms of the orbital velocity w rather than the magnetic moment $\hat{\mu}.$ The two representations are roughly equivalent, but integration in $\hat{\mu}$ would be much more difficult, since its relationship to w is defined via the electromagnetic field tensor at the gyrocenter position, which would also require transformation before integrating.

The gyrokinetic transformation has been found by expansion in orders of  $\varepsilon$ and $\lambda$, and this expansion has to be exploited here again:

$$j^{\mu }(\tens{x})=\sum_{i,j=0}\varepsilon ^{i}\lambda ^{j}j_{(ij)}^{\mu }.$$ (101)

The distribution function may have short-wavelength wave contributions, which can be separated from the slowly changing background by using the ordering scheme

$$% f_{\alpha }\left( x^{\prime \mu },w,u_{\parallel }\right) =... ...ft[ {\rm i}k_{\mu }x^{\prime \mu }\right] {\rm d}^{4}\tens{k}.$$ (102)

In these terms the zero-order current and charge density $j_{(00)}^{\mu }$is easy to find:

$$j_{(00)}^{\mu }=2\pi c\sum_{\alpha }q_{\alpha }\int \left( {u... ...arallel }\right) \frac{w{\rm d}w{\rm d}u_{\parallel }}{u_{o}},$$ (103)

while other orders require expansion.

6.1 Order ${\lambda ^{0}}$

The background distribution function $f_{\alpha }^{(0)}$ is slowly changing on the Larmor-radius/gyrotime-scale, and thus can be expanded into a Taylor series around the particle position. Note that the displacements $r_{i}^{\mu }$ are given by the gyrokinetic transformation as functions of $\left( x^{\prime \mu },w,u_{\parallel }\right)$ and thus should be expanded as well. Keeping terms of the order of  ${\varepsilon ^{2}}$ and larger, we find

$$f_{\alpha }^{(0)}\left( x^{\prime \mu }\right) =f_{\alpha }^{... ...al x^{\mu }\partial x^{\nu }}+o\left (\varepsilon ^{2}\right).$$ (104)

Here all functions on the right-hand side can be regarded as functions of the particle position $x^{\mu}$. Meanwhile, due to dependence of the tetrad in the velocity representation on $\tens{x}$, and due to the use of winstead of $\hat{\mu}$, all other functions under the integral in Eq. (99) do not require transformation.

The first order displacement of the gyrocenter is given by

$$r_{1}^{\mu }=\frac{w}{qH}(l^{\prime \prime \mu }\cos \phi -l^{\prime \mu }\sin \phi ).$$ (105)

Using it, we find the first-order current density as

$$j_{(10)}^{\mu }=-2\pi c\sum_{\alpha }q_{\alpha }\int \frac{w^... ...al x^{\nu }}\frac{w{\rm d}w {\rm d}u_{\parallel }}{u_{o}}\cdot$$ (106)

Its structure corresponds to the usual diamagnetic current.

Similarly, it is necessary to solve Eq. (42) for $r_{2}^{\mu }$ in order to calculate the second-order current density. The form of this equation (for $\lambda =0$) is such that $r_{2}^{\mu }$ may have components proportional to first and second powers of trigonometric functions of $\phi ,$ namely, $r_{2}^{\mu }=r_{21}^{\mu }+r_{22}^{\mu }.$ It is multiplied by the velocity, which has terms with the zero and first powers of the same, so that a range of first- to third-power terms are produced. However, of that range only the second-power terms can produce a non-zero contribution after integrating in $\phi$, so that in the following we ignore other terms. The same argument applies to the expansion (104) as a whole.

First we solve for $r_{21}^{\nu },$ which satisfies

$$qr_{21}^{\nu }F{_{\mu \nu }-}\bar{u}_{\nu }\frac{\partial r_{... ...\frac{\partial \bar{u}_{{\mu }}}{ \partial x^{\prime \nu }}=0,$$ (107)

then we find

$$\left\langle u^{\mu }r_{21}^{\lambda }\right\rangle =-\frac{w... ...frac{\partial \bar{u}_{{\kappa }}}{\partial x^{\nu }}\right) .$$ (108)

The next combination of terms is

$$\left\langle u^{\mu }\left( r_{22}^{\lambda }-\frac{1}{2}r_{1... ...tial x^{\nu }}\right\rangle \equiv \bar{u}^{\mu }X^{\lambda },$$ (109)

where $X^{\lambda }$ satisfies

$$qF{_{\mu \lambda }X^{\lambda }=}\frac{1}{2}\left\langle r_{1}... ...rtial r_{1}^{\nu }}{\partial x^{\prime {\mu }}}\right\rangle .$$ (110)

After substitutions we find

$${X^{\lambda }=}\frac{w^{2}}{4q^{2}H}D^{\lambda \mu }\left( l^... ...\prime \nu }{\rm Rot}_{\nu }l_{\mu }^{\prime \prime }\right) ,$$ (111)

where ${\rm Rot}_{\nu }f_{\mu }\equiv \partial f_{\mu }/\partial x^{\nu }-\partial f_{\nu }/\partial x^{\mu },$ as before.

The last remaining terms in Eq. (104) contribute to the current via combination

$$\left\langle r_{1}^{\nu }\frac{\partial }{\partial x^{\nu }}\... ...rtial f_{\alpha }^{(0)}}{\partial x^{{\mu }}}\right] \right) .$$ (112)

As a result, we can get the second-order current density as $j_{(20)}^{\mu }=j_{(20)\Vert }^{\mu }+j_{(20)\bot }^{\mu },$ where

$$j_{(20)\Vert }^{\mu }=2\pi c\sum_{\alpha }q_{\alpha }\int \fr... ...ambda }}}\right) \frac{w {\rm d}w{\rm d}u_{\parallel }}{u_{o}}$$ (113)

belongs to the $\left( {\rm l,\tau }\right)$-plane, while

$$j_{(20)\bot }^{\mu }=2\pi c\sum_{\alpha }q_{\alpha }\int \fra... ...l x^{{\lambda }}}\frac{w{\rm d}w{\rm d}u_{\parallel }}{ u_{o}}$$ (114)

has no $\left( {\rm l,\tau }\right)$ components.

The parallel component, $j_{(20)\Vert }^{\mu }$, represents the FLR corrections to the charge density and the parallel current density. It is rarely used, since it is two orders smaller than $j_{(00)}^{\mu }$ and has the same direction. The perpendicular component, $j_{(20)\bot }^{\mu }$, describes currents caused by the drifts in inhomogeneous fields.

6.2 Order ${\lambda ^{1}}$

First we calculate the transformation to the particle position variables

$$% f_{\alpha }^{(1)}\left( \tens{x},w,u_{\parallel },\phi \rig... ...u }{\sum_{s}}\varepsilon ^{s}r_{s}^{\mu }\right] {\rm d}^{4}k.$$ (115)

Retaining in the exponent terms of order $\varepsilon ^{0}$ only [while $\ k\varepsilon \sim O(1)$], we get

$$% f_{\alpha }^{(1)}\left( \tens{x}\right) =\frac{1}{4\pi ^{2}... ...mbda D^{\mu \nu }\tilde{a}_{\nu }\right) \right] {\rm d}^{4}k.$$ (116)

The last term in the exponent can be neglected in this order, and the rest can be expanded into the Fourier series using Eq. (47) as

$$f_{\alpha }^{(1)}\left( \tens{x}\right) =\frac{1}{4\pi ^{2}}\... ...( \xi \right) {\rm e}^{{\rm i}n(\phi -\phi _{0})}{\rm d}^{4}k,$$ (117)

where $\xi =k_{\bot }\rho ,\ \rho =w/qH,\ \ k_{\bot }=\sqrt{\left( k_{\nu }l^{\prime \nu }\right) ^{2}+\left( k_{\nu }l^{\prime \prime \nu }\right) ^{2}},$ as before, $\ \tan \phi _{0}=\left( k_{\mu }l^{\prime \prime \mu }\right) /\left( k_{\mu }l^{\prime \mu }\right)$. Then,

$$% j_{(01)}^{\mu }=\sum_{\alpha }\frac{cq_{\alpha }}{2\pi }\in... ...frac{w^{2}}{qH}\frac{ J_{1}\left( \xi \right) }{\xi }\right\}.$$ (118)

Note that the above expression describes the parallel current due to the wave, as well as the corresponding diamagnetic contribution.

6.3 Order ${\lambda ^{2}}$

The nonlinear corrections follow while the formerly neglected last term in the exponent of Eq. (116) is taken into account. The resulting first-order correction to $f_{\alpha }^{(1)}$ looks like

$$% \delta f_{\alpha }^{(1)}\left( \tens{x},w,u_{\parallel },\p... ... \xi \right) {\rm e}^{{\rm i}n(\phi -\phi _{0})} {\rm d}^{4}k,$$ (119)

where $\tilde{a}_{\mu }$ is taken at the particle position.

In principle, $a_{\mu }$ is a prescribed function of the particle position, but its oscillating part is not. It is necessary to transform the zero-order solution, Eq. (50),

$$% \tilde{a}_{\mu }(\tens{x})=a_{\mu }(\tens{x})-\bar{a}_{\mu ... ...{\rm e}^{{\rm i}\tens{k}^{\prime }\tens{x}^{\prime }}\right) ,$$ (120)

back to particle position coordinates via the inverse transformation (100) . Then

$$\tilde{a}_{\mu }(\tens{x})=\frac{1}{4\pi ^{2}}\int {\rm d}^{4... ...\right) {\rm e}^{{\rm i}m(\phi -\phi _{0}^{\prime })}\right] .$$ (121)

Thus the combination appearing in Eq. (119) can be rewritten as

$$% \tilde{a}_{\nu }\sum_{n=-\infty }^{+\infty }J_{n}\left( \xi... ... {\rm e}^{{\rm i}n(\phi -\phi _{0}^{\prime \prime })}\right] ,$$ (122)

where we have used the addition theorem for Bessel functions, Eq. (60), again, so that $\tens{k}^{\prime \prime }=\tens{k+k}^{\prime },$ etc.

Now we can calculate

$$% j_{(02)}^{\mu } -\sum_{\alpha }\frac{cq_{\alpha }}{\left( 2\pi \right) ^{3}} \int... ... \tens{k}\right) {\rm e}^{{\rm i}\left( \tens{k+{k}^{\prime}} \right) \tens{x}}$$ = (123)

$$\times \left[ \bar{u}^{\mu }\left( J_{0}(\xi )-J_{0}\left( \xi ^{... ...1}\left( \xi ^{\prime \prime }\right) }{\xi ^{\prime \prime }}\right) \right] ,$$ (124)

where $b^{\mu \nu }=l^{\prime \nu }l^{\prime \prime \mu }-l^{\prime \mu }l^{\prime \prime \nu }$ as before.

One can see that the nonlinear correction to the current density is due to the interaction of the perturbation of the distribution function with the perturbation of the particle motion, satisfying the three-wave interaction criteria.

7 Conclusion

A closed set of relativistic gyrokinetic equations, consisting of the collisionless gyrokinetic equation and the averaged Maxwell equations, is derived for an arbitrary four-dimensional coordinate system.

The basic features of the approach are as follows:

• the guiding-center dynamics of charged particles and the gyrokinetic transformation are obtained in the second order in terms of the ratio of the Larmor radius to the inhomogeneity scale;

• the wave terms with arbitrary $k\rho_{{\rm L}}$ are described in the second-order approximation in the amplitude of the wave;

• the present description is closed, in the sense that it includes the gyrophase-averaged expression for the four-current density, which is used with Maxwell equations;

• the perturbative Lagrangian approach with a fully relativistic, four-dimensional covariant formulation allows us to improve on the existing limitations of the gyrokinetic theory: in comparison to equations derived by other authors it is now possible to consider relativistic drift velocities, i.e., $E/B\sim O(1)$, and the wave field is no longer limited in frequency and wavelength, i.e., $\omega /\Omega _{\rm c}\sim O(1),$ $k_{\Vert }\rho _{ {\rm L}}\sim O(1)$, so that the class of admissible waves is broader than the usual drift-Alfvén perturbations and can include magneto-sonic waves, for example;

• the covariant description of the gyro-dynamics of charged particles in a strong magnetic field allows for the possibility of the gravitation influencing the space-time metrics, and thus is suitable for the investigation of astrophysical objects such as accretion disks and relativistic jets.

The applicability of the present gyrokinetic equations to modelling accretion disks can be assessed using the following scheme. If one is interested in quasi-equilibrium structures without short-wave perturbations, then the equations can be significantly simplified by throwing out the wave terms. The other requirements are as follows:

• The plasma should be magnetized, and the Larmor radius should be smaller than both characteristic scales: the Keplerian orbit radius, and the inverse gradient of the electromagnetic field across the accretion disc. The second scale is likely to be smaller for thin disks. Thus, significant FLR (finite-Larmor-radius) corrections are likely to arise when the disk thickness is comparable to the Larmor radius.

• Collisions, and the cyclotron/synchrotron radiation should not change the particle energy too much during the time scale of the processes under investigation. The second limitation is less severe if the plasma is optically thick, i.e., the cyclotron emission is collective.

• One can describe plasmas even within the last stable orbit, but no closer than one Larmor radius from the events horizon of the black hole.

For problems involving turbulence and dynamo effects, the Larmor radius should be compared with the typical wavelength of the perturbation. Then, besides the above two items, one should consider the following:

• If the wavelength is smaller than the Larmor radius, one can still throw out the gyrokinetic wave terms and use the zero-order theory. Otherwise, the FLR-corrections given by the wave terms are significant.

The important role of the FLR-corrections indicates when one should prefer the gyrokinetic theory to the (possible) covariant drift theory. However, in the opposite case these two should be equivalent (in results as well as in complexity). Non-relativistic gyrokinetic/drift equations can be applied if both the plasma temperature and the orbital velocity are non-relativistic, or, if one is working locally, in a moving reference frame. Finally, the generally preferred MHD (magnetohydrodynamic) approach is directly applicable for cold collisional plasmas of the outer reaches of the accretion disk, or just for fast MHD perturbations. The interesting subject of the covariant gyrofluid equations, which combine the relative simplicity of the MHD with explicit conservation of the magnetic moment in the collisionless plasmas, will be addressed in our subsequent papers.

Acknowledgements

This work has been conducted via the cooperation program between the Trieste University, Italy, and the Budker Institute of Nuclear Physics, Novosibirsk, Russia.

References

 
• Abramowitz, M., & Stegun, I. A. (ed.) 1972, Handbook of Mathematical Functions (New York: Dover) In the text  
• Beklemishev, A., & Tessarotto, M. 1999, Phys. Plasmas, 6, 4487 In the text  
• Boghosian, B. M. 1987, Covariant Lagrangian methods of relativistic plasma theory, University of California, Davis In the text  
• Boozer, A. H. 1996, Phys. Plasmas, 3, 3297 In the text  
• Brizard, A. J., & Chan, A. A. 1999, Phys. Plasmas, 6, 4548 In the text  
• Cary, J. R., & Littlejohn, R. G. 1983, Ann. Phys., 151, 1 In the text  
• Cooper, W. A. 1997, Plasma Phys. Control. Fusion, 39, 931  
• Fradkin, D. M. 1978, J. Physics A, 11, 1069 In the text  
• Frank, J., King, A., & Raine, D. 1992, Accretion Power in Astrophysics, Cambridge Astrophysics Series, 20, 2nd ed. (Cambridge University Press) In the text  
• Hahm, T. S., Lee, W. W., & Brizard, A. 1988, Phys. Fluids, 31, 1940 In the text  
• Landau, L. D., & Lifshits, E. M. 1975, The Classical Theory of Fields, 4th ed. (Oxford: Pergamon) In the text  
• Littlejohn, R. G. 1979, J. Math. Phys., 20, 2445 In the text  
• Littlejohn, R. G. 1984, Phys. Fluids, 27, 976 In the text  
• Northrop, T. G. 1963, The Adiabatic Motion of Charged Particles (New York: Interscience), 27 In the text  
• Pozzo, M., & Tessarotto, M. 1998, Phys. Plasmas, 5, 2232 In the text