Vlasov versus reduced kinetic theories for helically symmetric equilibria 3 1 H. Tasso1, G. N. Throumoulopoulos2 0 2 1Max-Planck-Institut fu¨r Plasmaphysik, Euratom Association, n D-85748 Garching, Germany a J 2 University of Ioannina, Association Euratom-Hellenic Republic, 8 2 Section of Theoretical Physics, GR 451 10 Ioannina, Greece ] h p Abstract - m s Anew constant of motionforhelically symmetric equilibria inthe vicinity a l ofthemagneticaxisisobtainedintheframework ofVlasovtheory. Inview of p . this constant of motion the Vlasov theory is compared with drift kinetic and s c gyrokinetictheoriesnearaxis. Itturnsoutthatasinthecaseofaxisymmetric i s equilibria [H. Tasso and G. N. Throumoulopoulos, Phys. Plasmas 18, 064507 y h (2011)] the Vlasov current density thereon can differ appreciably from the p drift kinetic andgyrokinetic current densities. This indicates some limitation [ on the implications of reduced kinetic theories, in particular as concerns the 1 physics of energetic particles in the central region of magnetically confined v 7 plasmas. 1 4 6 . Accepted for publication in Physics of Plasmas 1 0 3 1 : v i X r a 1 I. Introduction Kineticequilibria mayprovide broaderandmoreprecise informationthan multifluidormagnetogydrodymanicequilibriaasthosegovernedbytheGrad- Shafranov equation. Since solving self consistently the kinetic equations is tough particularly in complicated geometries the majority of kinetic equilib- rium solutions are restricted to one dimensional configurations in plane ge- ometry, e.g. [1]-[8]. Of particular interest are equilibria with sheared toroidal and poloidal flows which play a role in the transition to improved confine- ment regimesintokamaksandstellarators, thoughunderstanding thephysics of this transition remains incomplete. Construction of kinetic equilibria is crucially related to the particle constants of motion which the distribution function depends on. In the framework of Maxwell-Vlasov theory only a couple of constants of motions are known for symmetric two-dimensional equilibria, i.e. the energy, H, and the momentum px3 conjugate to the ignor- able coordinate, x3, out of the four potential constants of motion. Therefore, for distribution functions of the form f(H,px3) only macroscopic flows and currents along the direction associated with x3 can be derived, e.g. toroidal flows for axisymmetric plasmas. The creation of poloidal flows requires addi- tional constant(s) of motion. This remains an open question despite the fact that in a previous study [9] we found locally in the vicinity of the magnetic axis the following new constant of motion: C = v +Iln|v |, where v is the z φ φ toroidal particle velocity, v the velocity component parallel to the axis of z symmetry and B = I/re the toroidal magnetic field near axis (r,φ,z are φ φ cylindrical coordinates). An additional open question remains a potential extension of the proof of non existence of magnetohydrodynamic axisym- metric equilibria with purely poloidal flows [10] to Vlasov equilibria. A third constant of motion was studied extensively in the Astrophysics literature [11]-[24]. In particular for axisymmetric astrophysical systems distribution functions depending on an approximate third constant of motion result in velocity ellipsoids with unequal axes in agreement with observational data for our Galaxy [14, 17, 24]. Owing to the contemporary and probably future limited computational efficiency for simulations in the framework of Vlasov theory approximate kinetic theories have been established in reduced phase spaces as the drift kinetic and gyrokinetic ones. In both theories the reduced phase space is five dimensional with three spatial components associated with the guiding center position, R(orthegyrocenter positionintheframeworkofgyrokinetic 2 theory), and a velocity component parallel to the magnetic field, v ; also, k the two components of the perpendicular particle velocity are approximated after a gyroangle averaging with the magnetic moment which is treated as an adiabatic invariant. A related underlying assumption for both reduced theories is that the ratio ǫ of the gyroperiod to the macroscopic time scale is small. In the drift kinetic theory ǫ is the same as the ratio of the gyroradius to macroscopic scale length while in the gyrokinetic theory small spatial variations are permitted. Because of the reduction of the phase space some information of the particle motion is missed. In this respect it may be noted that the reduced-phase-space kinetic theories are developed via expansions in ǫ, the convergence of which is not guaranteed. This gives rise to the question: is the missing information important? To address this question we compared axisymmetric Vlasov equilibria near the magnetic axis with drift kinetic and gyrokinetic ones [25]. It turned out that because of missing the aboveconstant ofmotion, C, inthelatter casethe onaxis current density can differ appreciably from the “actual” Vlasov current density. Also, unlikely Vlasov theory, the reduced kinetic theories can not distinguish a straight from a curved magnetic axis. The aim of the present contribution is twofold: first to extend the above local constant of motion, C, to the more general class of helical symmetric equilibria and second to compare the Vlasov helical symmetric equilibrium characteristics near the magnetic axis with respective drift kinetic and gy- rokinetic characteristics. The new local constant of motion is derived in Sec. II. Then, Vlasov theory is compared near axis with reduced kinetic theories in Sec. III. Sec. IV summarizes the conclusions. II. A third Vlasov constant of motion near magnetic axis We will employ the following form of Vlasov equation [26]: ∂f ∂f ∂f +v·∇f +e ·(E+v×B) +(e ·v×∇×v) = 0. (1) i i ∂t ∂v ∂v i i To derive (1) one uses general orthogonal coordinates (x1,x2,x3) with unit basis vectors e = ∇x /|∇x | (i = 1,2,3), expresses the “microscopic fluid” i i i velocity in the basis of e as v = v e and uses the “microscopic fluid” mo- i i i mentum equation, ∂v +v·∇v = E+v×B, ∂t 3 where E and B are the electric and magnetic fields consistent with Maxwell equations. The term “microscopic fluid” relates to the fact that the Vlasov equation is an approximation to the N-particles Liouville equation, which replaces the N particles by a continuum. This is sometimes termed “fluid approximation”. To avoid confusion, however, the usual term “particle” will be adopted in place of “microscopic fluid”. Also, for the sake of notation simplicityandwithoutlossofgeneralitywewillconsideronlyionsandemploy convenient units by setting m = q = c = 1 where m and q are the ion mass and charge and c is the velocity of light. In connection with the helical symmetry we specify the coordinates as (x1 = r,x2 = aφ+βz, x3 = −βφ+az), whereaandβ areparameters. Helical symmetry means that any quantity does not depend on x3. Translational symmetry and axisymmetry arethen recovered as particular cases for (a = 1, β = 0) and (a = 0, β = 1). The respective unit vectors are ae +βre are −βe e1 = er, e2 = (a2 +φ β2r2)1z/2, e3 = (β2 +za2r2)1φ/2. Note that these basis vectors are in general non orthogonal. A helical mag- netic axis however is located at a constant distance r from the z-axis. Since we will study the equilibrium near the magnetic axis we set r = 1 thereon 2 2 and make the choice a +β = 1. Consequently, the above basis vectors re- duce to (e1 = er, e2 = aeφ +βez, e3 = −βeφ +aez) and become orthogonal on axis, with e2 along and e3 perpendicular to the axis. Thus, (1) can be applied in the neighbourhood of axis. The magnetic field near axis can be written in the form I I ae2 −βe3 βe2 +ae3 B = reφ +B0ez = r a2 +β2 +B0 a2 +β2 . By the choice aB0 = βI we take on axis B = I/ae2 as it should be. Also, we examine the case of E = 0 on axis but ∇f 6= 0 thereon. The latter assumption may be regarded as extraordinary because ∇f is related to the density gradient on axis which for usual peaked density profiles is expected to vanish. The reason for assuming ∇f 6= 0 on axis is that for axisymmetric equilibria, if (1) is solved near axis by the method of characteristics under the assumption ∇f = 0 on axis the usual toroidal momentum constant of motion, p , is missed, while p is recovered if ∇f 6= 0 thereon [25]. φ φ 4 Using the above basis for the helically symmetric equilibria under consid- eration Eq. (1) near axis becomes ∂f ∂f ∂f ∂f ∂f v1 +v2 +w1 +w2 +w3 = 0, (2) ∂x1 ∂x2 ∂v1 ∂v2 ∂v3 where I I 2 w1 = (av2 −βv3) −v3 , w2 = av1(βv3 −av2), w3 = v1 +v1β(av2 −βv3). a a (3) The characteristics of (2) are given by the solutions of dx1 dx2 dv1 dv2 dv3 = = = = . (4) dv1 dv2 dw1 dw2 dw3 Integration of the last equation yields the new constant of motion: 2 2 4 2 2 3 3 β(β −a )(βv3 −av2)+aI ln|[a v2 −a β v2 +aβ v3 −β(I +a v3)]| C = v3 + (a2 −β2)2 for a 6= β, (5) 3 2 a (v2 −v3) C = 2v3 −v2 − for a = β. (6) I 2 The respective axisymmetric constant of motion, C = v + Iln|v |, is re- z φ covered from (5) for (a = 1,β = 0) and the translational symmetric con- stant of motion, C = v , for (a = 0,β = 1). Also, conservation of energy z 2 2 2 H = 1/2(v + v + v ) follows from the manipulation of equations (4) pre- 1 2 3 sented in Appendix. Distribution functions of the form f(H,C) could create a helical current on axis because of the dependence of C on v2 additionally to that in the logarithmic term. This is different from the axisymmetric case in which toroidal axisymmetric currents on axis are not possible at all for the same class of distribution functions because of the dependence of C on |v |. This difference may be expected because for translational symmetric φ equilibria which arerecoverable fromthe helically symmetric ones therespec- tive distribution functions f(H,C) = f(H,v ) can create currents on axis. z However, such a helical equilibrium is not physically acceptable because of the dependence of C on v3 (but not on v1) implying the creation of a flow perpendicular to the axis. Similar unphysical flows are created in the ax- isymmetric case because of the respective dependence of C on v (but not on z 5 vr = v1). As discussed in [9] possible reasons of this unphysical behaviour are: (i) the lack of a potential fourth constant of motion involving v1 in such a way that creation of regular poloidal flows be possible, (ii) non taking into account here the MHD property of coincidence of the magnetic surfaces with the flow surfaces, (iii) additional drifts because of the curvature and torsion of the magnetic field which are eliminated in translational symmetric geometry and (iv) potential damping of unphysical flows in the framework of a collisional kinetic theory. The above consideration shows that straight, circular and helical magnetic axes are well distinguished in the framework of Vlasov theory. Unlikely the axisymmetric case [25] the generalized momentum constant of motion, px3 = v3 + A3, where A(x1,x2) is the vector potential, can not follow from (4) though we have assumed ∇f 6= 0 on axis. This is obtained form the Lagrangian 1 2 2 2 L = (x˙1 ++x˙2 +x˙3)+Φ(x1,x2)+x˙ ·A(x1,x2) 2 even in the presence of an electric filed on axis associated with the electro- static potential Φ. Also, a generic non local derivation off px3 is provided in [27]. III. Comparison of Vlasov with drift kinetic and gyrokinetic theories. The near axis consideration of helically symmetric equilibria will be re- peated in the framework of the drift kinetic and gyrokinetic theories which were employed in Ref. [25] on an individual basis below. For convenience a brief review of these theories will also be given here because otherwise fre- quent reference to [25] would make reading of the present section tedious. Drift kinetic theory The drift kinetic theory established in Ref. [31] is based on the Little- john’s Lagrangian for the guiding center motion [32] extended to include the polarization drift in such a way that local conservation of energy is guaran- teed. Introducing the modified potentials, A⋆ and Φ⋆, this Lagrangian can be written in the concise form [33] L = A⋆ ·R˙ −Φ⋆, (7) 6 where A⋆ = A+v b+V , (8) k E 1 Φ⋆ = Φ+µB + (v2 +v2), (9) 2 k E E×B V = . (10) E B2 Here, Φ and A are the usual electromagnetic scalar and vector potentials and b = B/B. The drift kinetic equation for the guiding center distribution function f(R,v ,µ,t) (with µ˙ = 0) acquires the form k ∂f ∂f +V·∇f +v˙ = 0. (11) k ∂t ∂v k By means of the Euler-Lagrange equations following from (7) the guiding center velocity, V, and the “acceleration” parallel to the magnetic field, v˙ , k can be expressed as: B⋆ E⋆×b V = R˙ = v + , (12) kB⋆ B⋆ k k E⋆ ·B⋆ 1 v˙ = = V·E⋆, (13) k B⋆ v k k where ∂A⋆ ∂Φ⋆ E⋆ = − − , B⋆ = ∇×A⋆, (14) ∂t ∂R B⋆ = B⋆ ·b = B +v b·∇×b+b·∇×V . (15) k k E Explicit expressions for V and v˙ are given by Eqs. (3.24) and (4.17) of Ref. k [31]. Also, it is noted here that Eqs. (12) and (13) have similar structure as the respective gyrokinetic equations of Refs. [38] and [39] [Eqs. (5.39) and (5.41) therein]. Since B⋆ appears in the denominators of (12) and (13) a sin- k gularityoccursforB⋆ = 0. ForE = 0thissingularitycanbeexpressed bythe k critical parallel velocity v = −Ω/(b·∇×b), where Ω is the gyrofrequency. c Therefore, the theory is singular for large |v | at which V and v˙ diverge and k k consequently non-casual guiding center orbits occur and the guiding center conservation in phase space is violated [30]. It is the v -dependence of A⋆ k [Eq. (8)] that produces the singularity. In order to regularize the singularity vk in (8) can be replaced by an antisymmetric function g(z) with z = vk/v0, where v0 is some constant velocity [30, 31, 33]. The nonregularized theory 7 presented here for simplicity is obtained for g(z) = z. In the regularized theory g(z) ≈ z should still hold for small |z|. For large |z|, however, g must stay finite such that with v0 ≫ vthermal one has v0g(∞) < vc. A possible choice for g(z) is g(z) = tanhz. Since for a helically symmetric equilibrium the magnetic field on axis becomes purely helical and dependent only on r, one readily calculates near axis with E = 0 thereon: a a a ∇×b = ∇×e2 = rez, B⋆ = B+vkrez, Bk⋆ = B +vkrβ (16) dB 1 dB ∇B(r) = e , Φ⋆ = µB + v2, E⋆ = −µ e (17) dr r 2 k dr r and consequently 2 B v a µ dB V = vkB⋆ e2 + rBk ⋆ez − B⋆ dr e3, (18) k k k v˙ = 0. (19) k As expected on axis the guiding center velocity consists of a component parallel to B and the curvature and grad-B drifts. Also, the “acceleration” v˙ vanishes because there is no parallel force and the drift kinetic equation k (11) becomes ∂f V· = 0. (20) ∂R Becauseof(19)v isaconstant ofmotionandthereforedistributionfunctions k 2 of the form f(H,v ), where H = µB + 1/2v , either can or can’t produce k k helical currents by choosing f either odd or even function of v . This prop- k erty holds also for axisymmetric and translational symmetric equilibria [25] because (19) keeps valid irrespective of the kind of symmetry. Thus, un- like Vlasov theory the drift kinetic theory can not distinguish equilibria of straight, circular or helical magnetic axes. Also, according to the results of section III the respective obtainable near axis Vlasov class of distribution 1 functions f(H,C) can not produce physically acceptable currents on axis . Consequently, since near axis the overwhelming majority of the particles are 1Physically acceptable Vlasov currents can be produced by the class f(H,px3) which, however, corresponds to drift kinetic distribution functions of the form f(H,pd ), where x3 near axis pdx3 =A3+vkb3, as it follows from the Lagrangian(7). 8 passing, the parallel currents constructed in the framework of the drift ki- netic theory may differ from the “actual” ones. Also, unlike in the case of axisymmetric equilibria [25], the B⋆-singularity is present as indicated by Eq. k (16). Gyrokinetic theory We will use the gyrokinetic equations of Ref. [35] which have been em- ployed to a variety of applications (see for example Refs. [40, 41, 42]). Eq. (11) remains identical in form where f(R,v ,µ,t) is now the gyrocenter dis- k tribution function for ions. The gyrocenter velocity and “acceleration” are given by B0 V = vkb0 + B⋆ VE +V∇B +Vc , (21) 0k (cid:16) (cid:17) 1 v˙k = − V· ∇Φ+µ∇B0 . (22) v k (cid:16) (cid:17) Here, B0 is the equilibrium magnetic field, b0 = B0/B0, B0⋆k = (B0 +vk∇×b0)·b0, Φ stands for the perturbed gyroaveraged electrostatic potential, and the E×B-drift velocity V , the grad-B drift velocity V∇ , and the curvature E B drift velocity V are given by c ∇Φ×∇B0 V = − , (23) E 2 B 0 µ V∇ = B0×∇B0, (24) B ΩB0 µv2 B2 Vc = k2b0×∇ p0 + 0 . (25) ΩB0 2 ! Note that as in the case of drift kinetic theory a similar singularity occurs at B⋆ = 0. In numerical applications this singularity was “avoided” by 0k approximating B0⋆k = B0 (see for example Refs. [40, 41]). Consideration of the above equations for a helically symmetric equilibrium with E = 0 near axisyields relationssimilarto(18), (19)and(20). Therefore, theabovefound discrepancies of the drift kinetic theory with the Vlasov one persist in the 9 framework of the gyrokinetic theory. The structure of the reduced kinetic equations in conjunction with the symmetry of the equilibrium considered clearly indicate that this conclusion is independent of the particular drift kinetic or gyrokinetic equations. IV. Summary We have found a new constant of motion near the magnetic axis of heli- cally symmetric equilibria in the framework of Vlasov theory [Eqs. (5) and (6)]. On account of this constant of motion a comparison of the Vlasov equation with either the drift kinetic or the gyrokinetic equation indicates that the current densities near the magnetic axis in the former case may be different from the ones in the latter case. Also, unlike Vlasov, the reduced kinetic theories cannot distinguish equilibria with straight, circular or helical magnetic axes. This discrepancy is due to the loss of new Vlasov constant of motion in the reduced phase space thus indicating that this reduction, even if made rigorously so that local conservation laws and Liouvillean invariance of the volume element is guarantied, is associated with the loss of nontrivial physics. This could put certain limits on the conclusions from drift kinetic or gyrokinetic simulations for time scales large compared to the gyroperiod which is the case of an equilibrium. Though the derivations of both theories are formally correct the convergence of the ǫ-expansions is presumably not uniform for all times which could lead to the noticed discrepancies calculated in this study. In addition, a singularity which occurs in both drift kinetic and gyroki- netic theories for large parallel particle velocities, present for helically sym- metric equilibria, is usually eliminated in the literature by a rough approxi- mation. Alternative to the rather artificially imposed regularization reported in Sec. III, one might remove this singularity in the integrals associated with moments of the drift kinetic equations (see the generic relations (8.14)-(8.17) of [31] for the self consistent charge, current, energy and energy flux densi- ties) by Cauchy principle value integration provided that such an integration would be physically justifiable. This requires further investigation. Appendix: Energy conservation near axis From dx1 dv1 dx2 dv2 = and = v1 w1 v2 w2 10