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 Kulsrud3 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 ordering5 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 of8 

ω/Ω(kρi)e(ϕv·A)/TekρiO(ϵ),

where ω is the frequency of interest, Ω is the cyclotron frequency, ϕ and A are the perturbed electrostatic and vector potentials, respectively, k 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ρiO(1) and k/kO(ϵ). 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

Fαt+v·Fαx+qm[E+1cv×(B0+δB)]·Fαv=0,
(1)

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

E=ϕ(1/c)A/t,δB=×A,

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

L=12mv2qϕ+qcv·A,

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

Fαt+(v+qαAmαc)·Fαx+qm[(ϕ1cv·A)+1cv×B0]·Fα(v+qαA/mαc)=0,

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

Fαt+v·Fαx+qm[(ϕ1cv·A)1cAt+1cv×(B0+δB)]·Fα(v+qαA/mαc)=0,

where the subscripts and denote the direction parallel and perpendicular to B0, respectively, and the approximation of

v+qαAmαcv,

based on the ordering of

qαvA/cTα1
(2)

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· and v×(B0+δB) 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

vvb+cB0E×b,
(3)

where

E=(ϕv·A/c)(1/c)A/t,b=b̂0+δB/B0,b̂0=B0/B0,

and

δB=×A.

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, ϕ, and vector potentials, A and A, but without the toroidal effects, take the form of

Fαt+[vbcB0(ϕ¯v·A¯/c)×b̂0]·Fαxqm[(ϕ¯v·A¯/c)·b+1cA¯t]Fαv=0,
(4)

where

2ϕ+ωpi2Ωi22ϕ=4παqαFα¯dvdμ,
(5)
2A1vA22At2=4πcαqαvFα¯dvdμ,μBμ/Bconst.,
(6)

μv2/2,B=B0+δB, the second term on the LHS of Eq. (5) is valid only for kρi1 and the bar¯ quantities denote gyrophase averages. Note that the Darwin approximation ignores the transverse displacement current in Ampere's Law when (ω/kc)21, 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

vpL=(mic2/eB2)(ϕ/t),

and the transverse ion polarization drift of

vpT=(mic/eB2)(2A/2t),

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α, 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α in Eq. (4) is independent of the gyrophase angle associated with the rotation of v, based on the gyrokinetic ordering argument, the gyrophase-averaged quantity of v·A, using Stokes' Theorem, becomes12 

v·A¯=Ω2π02π0ρδBrdrdθ,
(7)

where ρ=v/Ω and Ω=eB0/mc. The energy conservation of the system, i.e., Eqs. (4)–(6), by assuming that ω2k2vA2, becomes

ddt(12v2+μ)(meFe+miFi)dvdμ+ωci2Ωi2|Φ|28π+|A|28πx=0,
(8)

where Φϕv·A¯/c, which includes both the electrostatic perturbation ϕ and the magnetic perturbation δB via Eq. (7), and x denotes spatial average. Here, the approximation of

v·vv·v+v·(v+2qαmαcA)

has been used to calculate the particle kinetic energy by ignoring higher order A terms. Thus, the energy conservation to the quadratic order in the perturbed potentials are independent of A, where δB×A for kk.

For comparison, using the guiding-center approximation of Eq. (3), we obtain, from Eq. (1), the governing drift kinetic equation as7 

Fαt+[vbcB0(ϕ+1cAt)×b̂0]·Fαxqm[ϕ·b+1cAt]Fαv=0.
(9)

Here, a special treatment of A/t needs to be developed with procedures similar to that for A/t.14 So, for numerical purpose, Eq. (4) is more preferable, since the quantity of A 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., vA/cϕ for Eq. (4), which is consistent with Eq. (2) and vAv/cϕ 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,

Fαt+dRdt·FαR+dvdtFαv=0,dRdt=vb*+v22Ωα0b̂0×lnB0cB0Φ¯×b̂0,
(10)

and

dvdt=v22b*·lnB0qαmα(b*·Φ¯+1cA¯t),

where

Fα=j=1Nαδ(RRαj)δ(μμαj)δ(vvαj),

Ωα0qαB0/mαc,b*b+(v/Ωα0)b̂0×(b̂0·)b̂0, and b=b̂0+×A¯/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 as7 

J(x)=J(x)+JM(x)+Jd(x)=αqαFα(R)(v+v+vd)δ(Rx+ρ)dRdvdμφ,

where

vd=v2Ωαb̂×(b̂·R)b̂+v22Ωαb̂×RlnB.

For kρi1, they can be written as7 

JM(x)=α×cb̂Bpα,

and

Jd=cBα[pα(×b̂)+pαb̂×(lnB)],

where

pα=mα(v2/2)Fα(x)dvdμ,

and

pα=mαv2Fα(x)dvdμ.

Now, the current density takes the form of

J=JM+Jd=cBα[b̂×pα+(pαpα)(×b̂)]cBαb̂×pα,
(11)

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α=pαpα can be used by assuming some kind of isotropization process or by assuming (×b̂) is small.

With the gyrokinetic Poisson's equation of

ωpi2Ωi22ϕ=4πρ,
(12)

by assuming that (k/k)2(ωpi/Ωi)21 in Eq. (5) and the parallel Ampere's law for the electrons of

2A=4πcJ,
(13)

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

J=cBb̂×p
(14)

for the current density associated for a given pressure profile, where p ≈ pi and

ddt2ϕ4πvA2c2·(J+J)=0
(15)

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

E1cAtb·ϕTee1pepex0,
(16)

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

dpdt=0,
(17)

implying that the energy and mass convect together, where

ddttcBϕ×b·,

and pe = neTe for the electrons, we have a complete set of gyrokinetic MHD equations. Thus, with J given by Eq. (13), Eq. (14)–(17) are the governing gyrokinetic MHD equations, which, for E0 in Eq. (16), has the conservation property of

t18π(|ϕ|2+vA2c2|A|2)dx=vA2c2E·Jdx,
(18)

for a collisionless plasma. It can be regarded as Poynting's Theorem for our system. Furthermore, when ϕ0 in Eq. (15), we recover the MHD equilibrium as

·(J+J)=0.
(19)

Thus, we have finally accomplished what Strauss17,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 term and, in turn, gives us the MHD equilibrium of Eq. (19). Note that only a part of J 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

ni|particleni=1+12ρi21pi2pi,
(20)

where piniTi and ρivti/Ωi is the ion gyroradius and vti is the ion thermal velocity. From the gyrokinetic Poisson's equation, Eq. (12), it gives rise to an equilibrium E × B velocity of

vE×B12b̂×pipicTieB,

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

JE×B(x)=αqαvE×B(R)Fα(R)δ(Rx+ρ)dRdμdvφ.

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

J=cBb̂×p+eniρi22[2vE×B+vE×Bpi2pi],

which can then be approximated as

JcBb̂×(p)[114ρi22pp],
(21)

by assuming that 2vE×B0 and letting pip. The correction term resembles the expansion of the ordinary Bessel function J0(ρi2pi/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

[B28π+p(114ρi22pp)]0.

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 and p, 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.

1.
G. F.
Chew
,
M. L.
Goldberger
, and
F. E.
Low
,
Proc. R. Soc. London, Ser. A
236
(
1204
),
112
118
(
1956
).
2.
S. I.
Braginskii
, “
Transport processes in a plasma
,” in
Review of Plasma Physics
, edited by
M. A.
Leontovich
(
Consultants Bureau
,
New York
,
1965
), Vol.
1
.
3.
R. M.
Kulsrud
,
Handbook of Plasma Physics
, edited by
M. N.
Rosenbluth
and
R. Z.
Sagdeev
(
North-Holland Publishing Co.
, 1983); Volume
1
:
Basic Plas Physics I
, edited by
A. A.
Galeev
and
R. N.
Sudan
(
North-Holland Publishing Co.
,
1983
).
4.
E.
Frieman
,
R.
Davidson
, and
B.
Langdon
,
Phys. Fluids
9
,
1475
(
1966
).
5.
P. H.
Rutherford
and
E. A.
Frieman
,
Phys, Fluids
11
,
569
(
1968
).
6.
W. W.
Lee
,
Phys. Fluids
26
,
556
(
1983
).
7.
W. W.
Lee
and
H.
Qin
,
Phys. Plasmas
10
,
3196
(
2003
).
8.
A. M.
Dimits
,
L. L.
LoDestro
, and
D. H. E.
Dubin
,
Phys. Fluids B
4
,
274
(
1992
).
9.
H. C.
Corben
and
P.
Stehle
,
Classical Mechanics
, 2nd ed. (
Wiley
,
1966
).
10.
D. E.
Dubin
,
J. A.
Krommes
,
C.
Oberman
, and
W. W.
Lee
,
Phys. Fluids
26
,
3524
(
1983
).
11.
C. W.
Nielson
and
H. R.
Lewis
, “
Particle-code models in non-radiative limit
,”
Methods in Computational Physics
(
Academic Press
,
1976
), Vol.
16
, p.
367
.
12.
P.
Porazik
and
Z.
Lin
,
Phys. Plasmas
18
,
072107
(
2011
).
13.
A.
Brizard
,
Phys. Fluids B
4
,
213
(
1992
).
14.
E.
Startsev
and
W. W.
Lee
,
Phys. Plasmas
21
,
022505
(
2014
).
15.
A. J.
Brizard
and
T. S.
Hahm
,
Rev. Mod. Phys.
79
,
421
(
2007
).
16.
W. W.
Lee
,
S.
Ethier
,
R.
Kolesnikov
,
W. X.
Wang
, and
W. M.
Tang
,
Comput. Sci. Discovery
1
,
015010
(
2008
).
17.
H. R.
Strauss
,
Phys. Fluids
19
,
134
(
1976
).
18.
H. R.
Strauss
,
Phys. Fluids
20
,
1354
(
1977
).
19.
W. W.
Lee
and
R. A.
Kolesnikov
,
Phys. Plasmas
16
,
044506
(
2009
).
20.
F. I.
Parra
and
P. J.
Catto
,
Plasma Phys. Controlled Fusion
50
,
065014
(
2008
).
21.
W. W.
Lee
,
J. Comput. Phys.
72
,
243
(
1987
).
22.
W. W.
Lee
,
E. A.
Startsev
,
S. R.
Hudson
,
W. X.
Wang
, and
S.
Ethier
, “
Multiphysics/multiscale coupling of microturbulence and MHD equilibria
,”
Bull. Am. Phys. Soc.
60
(
19
),
JP12 127
(
2015
).