The nonlocal radio-frequency response of a toroidal plasma

PhysicsLettersAl75(1993)45—52 North-Holland

PHYSICS LETTERS A

The nonlocal radio-frequency response of a toroidal plasma P.U. Lamalle Laboratoire dePhysique des Plasmas — Laboratorium voor Plasma!ysica, Association “Euratom—Etat Beige” — Associatie “Euratom—BeigischeStaaz”~ EcoleRoyale Militaire — Koninklijke Militaire School, B-1040 Brussels, Belgium Received 19 January 1993; accepted for publication 3 February 1993

Communicated by M. Porkolab

For the first time, a model of the nonlocal plasma radio-frequency response has been developed which, while using standard guiding centre variables, achieves the same high level of symmetry as Hamiltonian descriptions. The full richness of particle motion in tokamak geometry is taken into account. Various wave—particle phase decorrelation models incorporated in the theory enhance its physical relevance. Preliminary numerical results show important differences from simplerapproaches.

1. Introduction The theoretical description of radio-frequency linear wave propagation and absorption in tokamak plasmas relies on the Maxwell—Vlasov system [1] Vx (VxE)— (w/c)2E=iwpo(>J~+J5),

(1)

with the boundary condition Ex ii = 0 on perfectly conducting surfaces with normal vector n. The radiofrequency (rf) source angular frequency and current density are co andj~(r),respectively; E(r) is the rfelectric field (all rf quantities are represented by complex amplitudes). The index fi labels the different particle species (their charge is qp and their mass mp). jp(r)~.Cp.E=qpJfpvdv3

(2)

is the self-consistent rfcurrent density due to species ft. and ap the corresponding conductivity tensor operator. fp(r, v) is the perturbed distribution function, satisfying the linearized Vlasov equation; using the method of characteristics, f~can be expressed as fp(r, v)=

—

-~-

J

exp[—iw(t’—t)] [E(r’)+v’XB(r’)]•

~dt’,

(3)

where f~is the equilibrium distribution function and B the rf magnetic induction. The integral is evaluated along the unperturbed particle orbit passing through (r, v) at the instant t’=t. In view of solving eq. (1) numerically with the finite element method, we have based the present work on its Galerkinformulation, obtained by dot product with the complex conjugate of a vector field F(r) and integration by parts over the whole tokamak chamber ~

0375-9601/93/S 06.00 © 1993 Elsevier Science Publishers B.V. All rights reserved.

45

Volume 175, number 1

~$

PHYSICS LETTERS A

29 March 1993

~Ep=~iw,LoJF*.jSdr3,

(4)

2P~E]dr3—iwu

[VXP~VXE—(co/c)

0~ p

where ~

(5)

JF*.jpdr3~qpJdr3Jdv3fpF*.v.

Requiring eq. (4) to hold for all continuous fields F satisfying the same metallic boundary conditions as E is equivalent to eq. (1) [21.In this formulation, the global quantity ~ (a bilinear form with respect to E and F) includes all information on the plasma if response: the case Fm F simply states Poynting’s theorem, and ~Ep= F~—iQp is the total complex power (active ix reactive) transferred from the fields to particle species p inside ~ choosing appropriate ö-functions for the components ofF yields backj(r) and the power deposition on the flux surfaces; finally, once discretized, ~#.EP is a direct input to finite element wave codes. Like in quasilinear diffusion theory, two types of formalisms have been used in the past to deal with the Maxwell—Vlasov system. Both assume the drift approximation to be valid, as this fundamental hypothesis yields integrability of particle motion [3}. They only differ by the set of variables describing this motion in the inhomogeneous equilibrium magnetic induction B0: either guiding centre (gc) variables historically the first formalism to be used, and the most appealing thanks to its physical transparency or the action—angle variables of the Hamiltonian formalism a powerful but much more abstract approach, first used in this context by Kaufman [4], and made accessible to numerical investigations by Gambier and Samain [5]. The gc formalism gave birth to various models of the conductivity tensor: a large amount of publications (see e.g. ref. [61) were first devoted to perpendicularly stratified geometries (i.e. inwhich B0~VI B0 I = 0); then Faulconer [71and Smithe et a!. [81 showed the important new effects arising from a parallel gradient in B0, using a “quasilocal” expansion of the integrand in eq. (3) near t’= t. Global wave numerical simulations in the actual tokamak geometry (where B0’V I B0 I is indeed nonzero) took place using some of these models locally in a simple, ad hoc fashion (see refs. [9,10] and a code of Lamalle to be published). At this very stage, one can say that the gc approach had lost its physical equivalence to the action—angle one: whereas the latter describes particle motion and wave propagation in a common geometry, the former has been applied solving the Maxwell system (eq. (1)) in a geometry after obtainingj~by eq. (2) in a different, simplified geometry. To our knowledge, nobody has stressed yet how dramatic the physical consequences of this hybrid use of the theory are: it is straightforward to show that the resulting total wave absorption by any species (P~=Re ‘IKEEP above) (i) is generally not linked to a wave—particle resonance condition as soon as the poloidal equilibrium magnetic field is taken into account (quite puzzling in a collisionless plasma!); (ii) cannot be demonstrated positive definite, even if the equilibrium is Maxwellian, in contradiction with the traditional pictures ofcollisionless damping and a number of works [11,121 (the Hamiltonian descriptions among others). Indeed the research exposed in this Letter was motivated by the numerical observation of a “Maxwellian plasma instability” in such a “hybrid” model. In the following sections, a general expression for the if response operator of a tokamak plasma (eq. (5)) is presented, based on the standard theory of particle motion in the drift approximation to first order in PLIL (PL is the Larmor radius and L a characteristic length of variation of B0). As the linearized Vlasov equation is solved in the actual curved magnetic equilibrium, the resultingfp (eq. (3)) is fully consistent with the geometry in which it has to be used. This bestows on the associated ~ remarkable symmetries, which allow global energy exchanges between waves and particles to take place only through resonant interactions. Before the present work, this fundamental property, valid for general equilibria, was the privilege of the Hamiltonian formalism; now both approaches have been reconciled. In the most general case where drift effects are retained, the if response is nonlocal over the whole plasma volume. It is valid for any tokamak poloidal cross section, includes nonuniform gc motion, banana (trapped) —

—

—

—

46

Volume 175, number 1

PHYSICS LETFERSA

29 March 1993

orbits with finite width, and tangent resonance interactions (i.e. Doppler-shifted wave—particle resonant layer tangent to the gc orbit). The if field structure is kept arbitrary (no geometric optics ansatz is made) and the equilibrium may be highly non-Maxwellian, an important characteristic of auxiliary heated plasmas. This Letter completes a preliminary study [13], which only dealt with Landau and fundamental cyclotron interactions, and extends it to all harmonics of the cyclotron frequency. The present work also permits a transparent analysis of different wave—particle phase decorrelation models, the role of collisions being crucial in reducing nonlocality (e.g. smearing out unrealistic high-order resonances between cyclotron and gc bounce motions). The simultaneous presence of these features in the theory and the minimal amount of assumptions required (basically the drift approximation and the ordering equation (6) below) also contribute to its originality. As a prelude to the full exploitation of the theory, a first numerical example typifies the approach and hopefully shows that more than an academic exercise has been performed.

2. The perturbed distribution function From now on, the species index

ft will be omitted whenever no confusion arises. We adopt the ordering (6)

27t/T~n,27t/1~QL,

W,Wc>>Wl,,>>

where w~,is the cyclotron frequency, w~, = 27/Tb a typical gc poloidal bounce frequency, Tb the poloidal bounce time, r~the collision time and TQL the characteristic evolution time of the equilibrium distribution. The last inequality ensures that j~is stationary on the time scales relevant to the present study. Another interesting form of eq. (3) for the study ofgeneral equilibria results from the Lagrangian potential of the if Lorentz force [141 and an integration by parts, (E+vXB)exp(—icot)as— ~4~(~dj~_V)[E.vexp(_iwt)]

leads to f(r, v)=

—

•E(r)+

—~--

5[~(~). ~-

+

~

.%7’]E(r’).~?exp[—ico(t’—t)] cIt.

(7)

In axisymmetric toroidal geometry and in the drift approximation, unperturbed particle trajectories are conveniently parametrized by the set {P,O,91,v,x,a,a}.

(8)

P= 9’— 2~tR2mçb/q is 2it/q times the gc toroidal canonical momentum, an approximate invariant of the motion in an axisymmetric tokamak; R is the distance to the torus axis and !Pthe poloidal flux function. Subscripts I and J~refer to B 0, and we use the local orthonormal basis (e1, e2, e1) defined by e1 IIV~P~ e2=e11 Xe1. 0 and ~‘ are respectively the poloidal and toroidal angles, a the gyrophase, a the sign of 2Ba/Bo v1, all at= time t; v is (p theis(in2p,Ba/mv2 the variant) particle velocity, and x the dimensionless adiabatic invariant x= (v1 /V) magnetic moment, B 0 I B0 I~Ba is a reference magnetic induction, e.g. B0 on the magnetic axis). The Larmor gyration around the gc is given by 0= a~+ v1 /.T, Ø( t’ = t) = a, where 1 /.~T=e1 (VX ~e1—Ve1 ~e2)is a geometrical coefficient depending on the local torsion and curvature of the magnetic field line [15,16]. The if field is decomposed into circularly polarized and parallel components relative to the local B0, —

—

E=E1e1 +E2e2 +E1e1 =E÷e÷ +E_e +E~e1, with e+ =(ei~ie2)/,~J~and E~= (E1 ±iE2)/%J~ (respectively left-, right-hand polarizations). Expanding in Taylor series about the gc, we have 47

Volume 175, number 1

PHYSICS LETTERS A

29 March 1993

~exp[—ip(Ø+~it)]M~[E(gc, t’)]

E(r’)=

Ordered in powers of the Larmor radius, the operator M~takes the form M~[...]

=

~

2kVkI Vk2

IPI+

k~0

k 2!

where V~=e~~V,k1 = k+ max (0, p) and k2 = k+ max (0, —p). (In the geometric optics approximation, this expression would reduce to a familiar Bessel series of the first kind [1,17].) Axisymmetry allows the independent study of toroidal modes (J~E÷,E, E1) = (J,; E~,,,E,,, E11,,) exp(in~),and we shall omit the mode index n wherever unambiguous. With the help of the previous_definitions, the integrand in eq. (7) can be parametrized with the set of variables of eq. (8) (e.g. v1 = V,.,/XB0/Ba, V11 = ov.,/i XB0/Ba, etc.). Before proceeding, eq. (7) merits an important remark; it was derived under the assumption of a collisionless plasma. Nevertheless, the range of integration therein extends to times t— t~>r~, rQ~during which collisions and heating cannot be neglected, as they ensure wave—particle phase decorrelation [8,18—20]. A more realistic description of the early times (t—t’--+oc) is obtained by introducing a suitable attenuation factor K(t—t’) in the integrand (,c(0)=l, ic(co) = 0). Smithe et al. [8] and Kasilov et al. [18] have obtained models of the form 5] (9) 1c(r)=exp[—(r/tL)’ which incorporate the effect of pitch angle scatteringby Coulomb collisions (the integer ö= 3 or 5, respectively; the decorrelation time tL < r~ 11).These models typically result in a much faster memory loss than Krook— Bhatnagar—Gross-like ones (KBG, ô=1, tL r~11).We assume tL to be a constant of the motion. In order to proceed, we still need to define a few quantities: the phase function yi~(t) cot+pØ( t) nço ( t), where p is the cyclotron harmonic number; the linear operators —

,

—

—

—

G±1~(E)=exp(_inco)v±M~[E±,,exp(inco)]/...J~, Go~(E)=exp(—inço)

V1M~[E,~exp(in~)]

in which V=v1e~+e1 X [v~(e11 ~V)e1~ +v~VB0/2B0]/co~ is the gc velocity; the equilibrium coefficients —

h

lofo+ — v ôv

2(

w~\Ofo+

v~\,~’°co) 8x

2mit.9f0 qw 19~•

(10)

Under the ordering of eq. (6), the equilibrium distribution only depends on the constants of the motion: fo = f0(P, v, x, a) for passing particles, and f0 =f~ (F, v, x) for trapped particles [21]. The coefficients h~are thus

also constants of the motion; an average over time r~,along the gc trajectory (i.e. over one poloidal bounce) will be indicated by <...>. As a consequence of periodicity of the gc motion in the poloidal plane, the orbit integral of eq. (7) finally reduces to a series of bounce harmonics (over index j), f(r, v)=f~(r,v)+ ~ h~A~~exp{i[nq’—p(a+ ~ (11) —

~ Here f and f~are evaluated at the particle position, whereas the exponential phase factor is evaluated at the gc. The other coefficients are constants of the motion; 48

Volume 175, number 1

PHYSICS LETTERS A

29 March 1993

: —o

.75

/-1/y

-‘

Fig. 1. The resonance curve 5(y) associated with the collision model (9) and 5=5. The dashed line shows the collisionless plasma case. Note that S( —y) = S~(y).

iq

(of0

Ofov

2xOfov

is only present when the equilibrium velocity distribution is anisotropic. The Fourier coefficients exp(_i11,j)) (12) 17N( t) ~v,,(t) t+JWbt have been introduced. (Note that in the limit of no and a new phase function toroidicity, w~~( 1) = <,ji,, > t, the bounce time harmonics reduce to poloidal angle ones, and the identification with the results of uniform plasma wave theory provides a useful check.) 4,, describes the wave—particle phase coherence over successive bounces, together with the smearing effect of collisions: in the KBG model, i”GL,,,..L(E)

—

generalizes the usual resonant denominator of homogeneous plasma theory; more generally, ( <~j’,,>—Jcob) tL], where the dimensionless resonance function S(y) = iK[ — 1Y/tL] /tL, K[s] being the Laplace transform of the attenuation ,c(r); S(y) is thus a characteristic of the collision model used, and eq. (9) allows its easy tabulation (see fig. 1 for d= 5). Equation (11) expresses the fully nonlocal if perturbation of the distribution function in tokamak geometry, including gc drift effects.

4,J= tLS[

3. The rf response operator The previous results directly yield the if current j and the conductivity tensor a; however, as discussed in the introduction, we prefer to concentrate on the global quantity ~ Changing from particle variables to gc 49

Volume 175, number 1

PHYSICS LETTERS A

29 March 1993

coordinates in eq. (5), the special form of the associated Jacobian suggests parametrizing the gc poloidal angle integrals as temporal ones [12]. Separating the contributions from passing and trapped particles, the integrals over the toroidal angle ço, gyrophase a and time are carried out (a step equivalent to angular integration in the canonical action—angle representation). The main result of this Letter follows, ~

~EP=~EPIoc+ ~=

(13)

JdVJdPJ~WFEP,

± I 0

Pi

0

where ~-J[N(r)r.E_r.($!(3fo+2xofo/ox)!dv3).E]dr3.

~EfiIoc=

(14)

This purely reactive term is associated with the anisotropic part off0 (N is the local density). 2~3t h * WEEP— ~ D ~ ~ ~J( ) .‘i ( ), ixq m a p,J=—~ ~‘

—

—

a=l forapassingparticle, = ~ for a trapped particle,

(15)

is a density of ~%V,E,in the constants of motion space (as above, summation over the toroidal mode index n is implicit). For trapped particles, ~ is independent of a. The most remarkable fact about this equation is that the operators acting on F and F are complex conjugate, whatever the equilibrium distribution function (the latter only intervenes through /z~).Radial drifts yield involved limits of integration X, P 1 and P2 in eq. (13), requiring a detailed study of the constants of motion space for every practical application (see the work of Rome and Peng [22]). If the drifts are neglected, X(P) =Ba/Bo( 9’=P, 0=0), P1 = 0 and P2= ~ where P= 9’e defines the plasma edge.

4. Discussion of the characteristics of the model The if power density transferred per (F, 2 3 7rqaVr~

~,

v, x, a) 2

ReWEEP= ~

orbit of species ft is (16)

(a is as above). In a collisionless plasma, this expression allows a power exchange only at the poles of 4,~,i.e. under the resonance condition co=p+ n <ç~>+fCOb, in agreement with the Hamiltonian theories [4,5,12], the test particle approach to quasilinear theory [20,231, and Puri’s recent work on toroidal Landau damping [24]. With a simple KBG collision model, Im4,~is always positive. A sufficient condition for positive power absorption density is then a uniform isotropic equilibriumfo=fo(V) with negative slope ensuring h~>~ 0— in agreement with ref. [12]. Equations (15), (16) thus cure the problems formerly faced by the gc variable approach and reconcile it with the action—angle formalism, as announced in the introduction. In more complicated situations where the equilibrium is anisotropic, includes a tail or is radially inhomogeneous, one cannot generally predict whether the plasma damps or amplifies the fields without solving eq. (4). One is led to the same conclusion when the collision models considered above are used with ö~ 3, as then —

50

Volume 175, number 1

PHYSICS LETTERS A

29 March 1993

Im4,~may be negative over parts of the j spectrum (fig. 1). This behaviour was already stressed by Auerbach in his study of the collisional damping of Langmuir waves [25].

5. The bounce integrals Equations (15), (16) require the evaluation ofthe Fourier transforms defined by eq. (12) over one gc bounce. Expanding the arbitrary if fields in poloidal harmonics exp(im0), this dependence can be included in the phase TI,,~.Under the basic ordering of eq. (1) for p~ 0, or in some cases of Cherenkov-like interaction for p= 0, the variations of fl~, are much larger than 27t, allowing an accurate evaluation ofA~1by the steepest descents method. Third order stationary points (sp) occur at tangent resonance (II,,,, =0: extrema in B0, vicinity ofbanana tip); they require a uniform asymptotic expansion covering cubic and quadratic phase time dependence as described by Dingle [26], (17) All rapid variations lie now in the phase factor. The sum is over all stationary phase points. A is a smooth function related to the Airy function and its derivative (Bécoulet et al. [19] have treated the limit of a circular large aspect ratio tokamak). Note that the factor A,,~in eq. (15) naturally describes the transition from correlated to uncorrelated sp (a delicate problem discussed by Catto and Myra [27]). In quasihomogeneous situations (for p,& 0, in a tiny region of radius p ~.Rocob/pwC around the magnetic axis; for p = 0, cases where 0, ~ const), or for very well trapped orbits, a Bessel expansion of the phase factor is more appropriate than eq. (17). Cattanei and Murphy have used that method [28] to study minority heating (p = ±1 terms) in large aspect ratio tokamaks, assuming uniform gc motion.

6. A first example Figure 2 illustrates some important effects included in the present theory: using typical parameters of the JET tokamak (major radius 3 m, Ba= 3.1 T, plasma current 4 MA), we study if heating of minority hydrogen ions at their fundamental cyclotron resonance (the resonant layer co= ~ is located around R = 3.4 m). We consider 10 keV ions with gc on a circular magnetic surface of minor radius 0.6 m (in this preliminary investigation, radial drifts and collisions are neglected), and a single toroidal (n = 20) and poloidal (m = 5) if field mode with4h~ left-hand polarization. The power (15) onlyindependent includes a p= 1 contribution: we norI E~I 2/2Ba to obtain simple units,density and a quantity of equilibrium distribution malize it by xv -

1

-

0

-

0

~

rt

~

-

I

I

120

140

‘

20

40

60

80

100

y,

Fig. 2. Normalized power absorption density versus outboard equatorial pitch angle. Comparison between the nonlocal model I (solid line) and a “quasiocal” one (dashed line). s indicates a 160 180 boundary between regions of passing and trapped motion. rt in-

(degree)

dicates the banana orbit grazing resonance with its tips.

51

Volume 175, number 1

PHYSICS LETTERS A

29 March 1993

function and if field amplitude. We show this quantity as a function of the particle outboard equatorial pitch angle y~(cos m = V11 /vJ o=o; mis used instead ofx and a to display both signs of v1 on a single figure). The dashed curve shows the corresponding quantity obtained with a “quasilocal” model in which constant V and B0 were used in evaluating the trajectory integrals of eq. (3). Both models agree on no absorption for Yo < 130 (these orbits do not cross the Doppler-shifted resonance layer), and behave alike for well passing orbits with y~—~ 180° (ordinary resonance crossing). m = 13°corresponds to a resonance layer tangent to the gc trajectory, giving a singular absorption in the local model, but a finite value in the nonlocal one where the resonance is crossed in a finite time. In the latter model, trapped particles correspond to
—

Acknowledgement It is a pleasure to thank Drs. D.W. Faulconer and R. Koch for their enthusiastic initiation to this domain of research and many constructive discussions, and Professor P.R. Vandenplas for his stimulating support. References [1] T.H. Stix, Wavesin plasmas (AlP, New York, 1992). [2] 0. Strang and G.J. Fix, An analysis of the finite element method, Seriesin automatic computation (Prentice-Hall, Englewood Cliffs, NJ, 1973). [3] A.J. Lichtenberg and M.A. Lieberman, Regular and stochastic motion (Springer, Berlin, 1983). [4] A.N. Kaufman, Phys. Fluids 15 (1972)1063; H.E. Mynick, J. Plasma Phys. 39 (1988) 303. [5] DJ. Gambier and A. Samain, Nucl. Fusion 25 (1985) 283. [6] M. Brambilla, Plasma Phys. Control. Fusion 31(1989) 723. [7] D.W. Faulconer, Plasma Phys. Control. Fusion 29 (1987) 433. [8] D. Smithe, P. Colestock, T. Kammash and R. Kashuba, Phys. Rev. Lett. 60 (1988) 801. [9] M. Brambilla and T. KrUcken, NucI. Fusion 28 (1988) 1813. [10] E.F. Jaeger and D.B. Batchelor, AlP Conf. Proc. 244 (1992) 159. [11] M. Brambilia and T. Krucken, Plasma Phys. Control. Fusion 30 (1988)1083. [12] D.N. Smithe, Plasma Phys. Control. Fusion 31(1989) 1105. [13] P.U. Lamaile, in: Europhysics Conference Abstracts 15 C (1991) IV-l93. [14] L.D. Landau and E.M. Lifshitz, Course oftheoretical physics, Vol.1 (Pergamon, Oxford, 1960); J.CJ. Nihoul, Cours moderne de mécanique rationnelle (Albin Michel, Paris, 1968). [15] A.I. Morozov and L.S. Solov’ev, in: Reviews ofplasmaphysics, Vol.2, ed. M.A. Leontovich (Consulants Bureau, New York, 1966). [16] R. Balescu, Transport processes in plasmas (North-Holland, Amsterdam, 1988) [17]R.Koch,Phys.Lett.A 157 (1991)399. [18] S.V. Kasiov, A.I. Pyatak and K.N. Stepanov, Nucl. Fusion 30 (1990) 2467. [191 A. Bécoulet, D.J. Gambier and A. Samain, Phys. Fluids B 3 (1991) 137. [201 D.W. Faulconer and M.P. Evrard, in: Europhysics Conference Abstracts 15 C (1991)111-297. [21] F.L. Hinton and R.D. Hazeltine, Rev. Mod. Phys. 48 (1976) 239. [22]J.A. Rome and Y.-K.M. Peng, Nucl. Fusion 19 (1979) 1193. [23] V.S. Belikov and Ya.I. Kolesnichenko, Plasma Phys. 24 (1982) 61. [24JS. Pun, in: Europhysics Conference Abstracts 16 E (1992) 13, 17. [25] S.P. Auerbach, Phys. Fluids 20 (1977) 1836. [261 R.B. Dingle, Asymptotic expansions: their derivation and interpretation (Academic Press, New York, 1973). [27] P.J. Catto and J.R. Myra, Phys. Fluids B 4 (1992) 187. [281G. Cattanei and A.B. Murphy, Nucl. Fusion 31(1991) 219.

52