The effort to obtain a set of MagnetoHydroDynamic (MHD) equations for a magnetized collisionless plasma was started nearly 60 years ago by Chew *et al*. [Proc. R. Soc. London, Ser. A **236**(1204), 112–118 (1956)]. Many attempts have been made ever since. Here, we will show the derivation of a set of these equations from the gyrokinetic perspective, which we call it gyrokinetic MHD, and it is different from the conventional ideal MHD. However, this new set of equations still has conservation properties and, in the absence of fluctuations, recovers the usual MHD equilibrium. Furthermore, the resulting equations allow for the plasma pressure balance to be further modified by finite-Larmor-radius effects in regions with steep pressure gradients. The present work is an outgrowth of the paper on “Alfven Waves in Gyrokinetic Plasmas” by Lee and Qin [Phys. Plasmas **10**, 3196 (2003)].

The search for a set of one-fluid hydromagnetic equations for collisionless plasmas from the Boltzmann equation was started sixty years ago by Chew, Goldberger, and Low (CGL).^{1} This attempt was made because the usual fluid equations were derived from collisional considerations.^{2} The CGL work was followed by Kulsrud^{3} as well as Frieman *et al*.^{4} However, the work on this interesting subject has not received much attention over the years, since the use of ideal MHD, based on usual MHD equations, even when applied to collisionless plasmas, has been proven empirically to be very useful. In the present paper, we will show the derivation of a set of collisionless gyrokinetic MHD equations by applying the gyrokinetic ordering^{5} on the Vlasov-Maxwell equations. The main ingredient of this connection is the contribution of the ion polarization drift to the quasineutrality condition,^{6} which we will explain. An initial attempt based on this methodology was made more than ten years ago.^{7}

Let us first re-visit the subject of the gyrokinetic approximation for the Vlasov-Maxwell equations based on the ordering of^{8}

where *ω* is the frequency of interest, Ω is the cyclotron frequency, $\varphi $ and **A** are the perturbed electrostatic and vector potentials, respectively, $k\u2225$ and *k*_{⊥} are the wave vectors parallel and perpendicular to the zeroth-order magnetic field, respectively, and *ϵ* is a smallness parameter. This ordering is less restrictive than the original ordering by Rutherford and Frieman,^{5} which is based on $k\u22a5\rho i\u223cO(1)$ and $k\u2225/k\u22a5\u223cO(\u03f5)$. This paper is closely related to that of Lee and Qin,^{7} but involves a new way of deriving the governing gyrokinetic equations, based on a set of Lagrange equations, as well as new physics insight in terms of gyrokinetic MHD equations arising from finite-Larmor-radius (FLR) effects.

The governing gyrokinetic Vlasov-Maxwell equations used in the present paper can be derived by first transforming the original Vlasov equation

to a new set of variables, where $F\alpha \u2261F\alpha (x,v,t)$ is the distribution function in six dimensional phase space, *α* denotes species

and $B0(x)$ is the equilibrium background magnetic field. Making use of the Lagrangian

as described, for example, by Corben and Stahle,^{9} we then obtain

where $\varphi $ and **A** are the perturbed scalar and vector potentials, respectively. Alternatively, by changing the phase variables, we can re-write the equation as

where the subscripts $\u2225$ and $\u22a5$ denote the direction parallel and perpendicular to **B**_{0}, respectively, and the approximation of

based on the ordering of

is used for the $dx/dt$ term in the above equation.

To derive the gyrokinetic Vlasov-Maxwell equations, one could follow the procedures used in Ref. 6 or, more formally, those of Ref. 10 based on the Lie transform and a non-canonical perturbation theory. For the present purpose, we derive them using the simplified method of Ref. 7 by applying the drift kinetic approximation for the velocities associated with the $v\xb7\u2207$ and $v\xd7(B0+\delta B\u22a5)$ terms in the above equation. The justification for the approximation is based on the fact that the formally derived gyrokinetic Vlasov equation, as given in Refs. 6 and 10, is actually a drift kinetic equation to the lowest order in *ϵ*. As such, one can use guiding center approximation for the particles as

where

and

This drift kinetic approximation for the velocities is consistent with the formulations given in Refs. 6 and 10. It can then be shown that the governing gyrokinetic equations based on the Darwin approximation,^{11} including both scalar potential, $\varphi $, and vector potentials, $A\u2225$ and $A\u22a5$, but without the toroidal effects, take the form of

where

$\mu \u2261v\u22a52/2,\u2009B=B0+\delta B$, the second term on the LHS of Eq. (5) is valid only for $k\u22a5\rho i\u226a1$ and the $bar\xaf$ quantities denote gyrophase averages. Note that the Darwin approximation ignores the transverse displacement current in Ampere's Law when $(\omega /kc)2\u226a1$, where *ω* is the frequency of interest. As such, the light waves can be eliminated from the simulation. The derivations of the gyrokinetic Poisson's equation, Eq. (5), and Ampere's law, Eq. (6), based on the longitudinal ion polarization drift of

and the transverse ion polarization drift of

respectively, can be found, for example, in Ref. 7. These additional terms associated with the density and current responses are the result of the gyrokinetic approximation for the distribution function, $F\alpha $, in Eq. (1), which breaks up into three different parts as given by Eqs. (4)–(6), respectively. Equation (4) is in agreement with the slab version of the equation used in Refs. 12 and 13. Assuming $F\alpha $ in Eq. (4) is independent of the gyrophase angle associated with the rotation of $v\u22a5$, based on the gyrokinetic ordering argument, the gyrophase-averaged quantity of $v\u22a5\xb7A\u22a5$, using Stokes' Theorem, becomes^{12}

where $\rho =v\u22a5/\Omega $ and $\Omega =eB0/mc$. The energy conservation of the system, i.e., Eqs. (4)–(6), by assuming that $\omega 2\u226ak2vA2$, becomes

where $\Phi \u2261\varphi \u2212v\u22a5\xb7A\u22a5\xaf/c$, which includes both the electrostatic perturbation $\varphi $ and the magnetic perturbation $\delta B\u2225$ via Eq. (7), and $\u27e8\cdots \u27e9x$ denotes spatial average. Here, the approximation of

has been used to calculate the particle kinetic energy by ignoring higher order $A\u22a5$ terms. Thus, the energy conservation to the quadratic order in the perturbed potentials are independent of $A\u22a5$, where $\delta B\u2225\u2248\u2207\u22a5\xd7A\u22a5$ for $k\u2225\u226ak\u22a5$.

For comparison, using the guiding-center approximation of Eq. (3), we obtain, from Eq. (1), the governing drift kinetic equation as^{7}

Here, a special treatment of $\u2202A\u22a5/\u2202t$ needs to be developed with procedures similar to that for $\u2202A\u2225/\u2202t$.^{14} So, for numerical purpose, Eq. (4) is more preferable, since the quantity of $A\u22a5$ in that equation can be easily calculated by Eq. (7). In the case when the physics of compressional Alfven waves is not important, the ordering to ignore these terms is also different for these two equations, i.e., $v\u22a5A\u22a5/c\u226a\varphi $ for Eq. (4), which is consistent with Eq. (2) and $vAv\u22a5/c\u226a\varphi $ for Eq. (9), which requires additional ordering.

Equations (4)–(8) are different from the pure electrostatic versions given by Refs. 6 and 7. These new equations include the all-important vector potential, **A**, in order for us to make connections to MHD equations. In fact, the present paper was inspired by the work of Porazik and Lin,^{12} who have used the equations, based on more rigorous Lie transformations, as given by Brizard and Hahm.^{15} The present paper can be viewed as a more intuitive way, similar to those of Refs. 6 and 7, to derive these equations with an approximate energy conservation of Eq. (8).

For the general toroidal geometry, the gyrokinetic Vlasov equation can be re-written as, e.g., Ref. 16,

and

where

$\Omega \alpha 0\u2261q\alpha B0/m\alpha c,\u2009b*\u2261b+(v\u2225/\Omega \alpha 0)b\u03020\xd7(b\u03020\xb7\u2207)b\u03020$, and $b=b\u03020+\u2207\xd7A\xaf/B0$. Again, the variables with subscript “0” represent equilibrium quantities.

Now let us look at the current density in original particle coordinates **x** calculated from the distribution function, *F*(**R**), in gyrocenter coordinates, which can be written as^{7}

where

For $k\u22a5\rho i\u226a1$, they can be written as^{7}

and

where

and

Now, the current density takes the form of

where the first expression resembles that of Chew, Goldberger, and Low.^{1} In order to recover the usual pressure balance equation for MHD, the approximation of $p\alpha =p\alpha \u2225\u2248p\alpha \u22a5$ can be used by assuming some kind of isotropization process or by assuming $(\u2207\xd7b\u0302)\u22a5$ is small.

With the gyrokinetic Poisson's equation of

by assuming that $(k\u22a5/k)2(\omega pi/\Omega i)2\u226b1$ in Eq. (5) and the parallel Ampere's law for the electrons of

by assuming that $(\omega /kvA)2\u226a1$ in Eq. (6), we can proceed to obtain a simple set of gyrokinetic MHD equations for collisionless plasmas. They are

for the current density associated for a given pressure profile, where *p* ≈ *p _{i}* and

as the vorticity equation, which can be obtained from Eqs. (10), (12), and (13). Together with parallel Ohm's law of the form

which can be derived by using Eqs. (4) and (6) as shown earlier by Ref. 3, as well as the incompressible adiabatic equation of state of

implying that the energy and mass convect together, where

and *p _{e}* =

*n*for the electrons, we have a complete set of gyrokinetic MHD equations. Thus, with $J\u2225$ given by Eq. (13), Eq. (14)–(17) are the governing gyrokinetic MHD equations, which, for $E\u2225\u21920$ in Eq. (16), has the conservation property of

_{e}T_{e}for a collisionless plasma. It can be regarded as Poynting's Theorem for our system. Furthermore, when $\varphi \u21920$ in Eq. (15), we recover the MHD equilibrium as

Thus, we have finally accomplished what Strauss^{17,18} set out to do, without using the aspect ratio (*a*-minor radius/*R*-major radius) expansion for a tokamak. As a result, we recover the quasineutral condition as given by Eq. (19), which Strauss was not able to do. Nevertheless, it should be noted that the pioneering work by Strauss using the fluid approach mentioned here has inspired many researchers in this area for years including the present author. The difference here is that our approach is purely kinetic in nature similar to those of Chew-Goldberger-Low,^{1} Kulsrud,^{3} and Frieman-Davidson-Langdon.^{4} However, because of the use of gyrokinetic ordering in the present paper, the resulting equations are totally different from those previous attempts. Equation (15) now includes the correct $J\u22a5$ term and, in turn, gives us the MHD equilibrium of Eq. (19). Note that only a part of $J\u22a5$ in Eq. (19) was included in Eq. (32) of Lee and Qin.^{7}

Another interesting aspect of finite Larmor radius gyrokinetics is the existence of the equilibrium zonal flows associated with zeroth-order inhomogeneity. Let us elaborate. As first pointed out by Lee,^{6} the zeroth-order inhomogeneity also contributes to an extra ion particle density in addition to the ion gyrocenter density as given by Eq. (40) in that paper. A more complete expression is given by Eq. (17) in Ref. 19, as well as those in Ref. 20, and can be written as

where $pi\u2261niTi$ and $\rho i\u2261vti/\Omega i$ is the ion gyroradius and *v _{ti}* is the ion thermal velocity. From the gyrokinetic Poisson's equation, Eq. (12), it gives rise to an equilibrium

**E**×

**B**velocity of

where **x** is the direction of the zeroth-order inhomogeneity. The corresponding current is given by

Consequently, by taking into account the difference between the electrons and the ions for the **E** × **B** drift due to the finite Larmor radius effects, we obtain a new pressure balance equation modified by the current associated with the equilibrium zonal flows as

which can then be approximated as

by assuming that $\u2207\u22a52vE\xd7B\u22480$ and letting $pi\u2261p$. The correction term resembles the expansion of the ordinary Bessel function $J0(\rho i\u2207\u22a52pi/pi)$, usually associated with FLR corrections.^{6,21} Thus, Eq. (21) should be used instead of Eq. (14) in the regions with steep pressure gradient. The corresponding plasma pressure balance from Ampere's law in slab geometry now becomes

Thus, we have shown in this paper that the Darwin gyrokinetic model of Eqs. (4)–(6) can indeed be reduced to a set of MHD equations in the collisionless limit under certain conditions, a pursuit started sixty years ago using the Vlasov-Maxwell equations.^{1} The key to such a connection is the presence of the ion polarization density in the gyrokinetic Poisson's equation, Eqs. (5) and (12), which was first identified by Lee.^{6} Not surprisingly, this set of equations is different from the conventional MHD equations, notably the absence of the compressional Alfven waves, which can be ignored through the ordering argument with the present formulation. Nevertheless, it can still recover the conventional MHD equilibrium (Eq. (19)). Furthermore, the present paper also points out the corrections to these equations due to FLR effects in the collisionless limit, which have no relationship to Braginskii's equations.

In the future, it would be interesting to include the higher order fluid moments, such as heat fluxes, etc., as well as the compressional components of the Alfven waves in these gyrokinetic MHD equations. Studies including the use of Eq. (11), which differentiates between $p\u2225$ and $p\u22a5$, should also be very interesting. Moreover, the connection between these gyrokinetic equations and the MHD equilibria, as shown in the present paper, suggests that it is feasible to devise an iterative scheme between a gyrokinetic code and an MHD equilibrium code with the purpose of minimizing turbulence and anomalous transport in tokamaks based on an iterative procedure, which first “decouples” the transport problem from the equilibrium problem, and then “couples” them through global parameter exchanges.^{22}

The author wishes to thank Professor Russell Kulsrud of Princeton University for his interest in this work and his critical comments, also to Dr. Peter Porazik and Dr. Stuart Hudson of PPPL for useful discussions. This work was partially supported by U.S. DoE, Grant DE-AC02-09CH11466.