Mathematics | Discrete mathematics » Kinetic Ballooning Modes and Finite Beta Effects on Drift Type Modes

Source: http://www.doksinet Chapter 4 Kinetic Ballooning Modes and Finite E¤ects on Drift Type Modes 4.1 Introduction In Chapter 3, we have seen that how low frequency, electrostatic ion acoustic mode described by ! = kk cs is coupled to drift mode, ! = ! e . The drift mode can also couple to the electromagnetic Alfven mode ! = kk VA : The Alfven mode is one of the fundamental, low frequency electromagnetic modes. In a low plasma, the Alfven velocity far exceeds the ion acoustic velocity, VA cs : We have already seen in Chapter 2 that the ballooning instability caused by the combination of pressure gradient and (bad) magnetic curvature is essentially the destabilized Alfven mode. There, our treatment was hydrodynamic, with no account taken of the e¤ects of nite ion Larmor radius and kinetic resonances. In this Chapter, the gyrokinetic equation derived in Chapter 3 will be generalized to implement electromagnetic elds in order to assess kinetic modications to the MHD ballooning

mode and nite e¤ects on predominantly electrostatic modes, such as the drift and ion temperature gradient modes. 4.2 Kinetic Equations (Two-Potential Approximation) Electrostatic waves can be fully described in terms of the scalar potential alone as seen in Chapter 3. Electromagnetic modes, on the other hand, involve magnetic perturbations which can be deduced from the vector potential A: The electric and magnetic elds are then fully described by the two potentials, E= r B=r 1 1 @A ; c @t (4.1) A: (4.2) Source: http://www.doksinet In the Vlasov equation, e @f + v rf + @t m 1 E+ v c B @f = 0; @v (4.3) the perturbed magnetic eld B can be eliminated via Faraday’s law r E= 1 @B ; c @t (4.4) and considering the three independent components of the electric eld E is su¢ cient. In analysis using the two potentials and A, there are a total of four components. However, in either Coulomb gauge ( r A = 0) or Lorentz gauge r A+ 1 @ = 0; c2 @t one of the four

components can be eliminated, and we still have three independent components to work with. In analyzing low frequency (! i) electromagnetic waves in a plasma, it is often convenient to adopt the following three elds as the independent eld components, ; Ak and Bk ; (4.5) where Ak and Bk are the perturbed vector potential and magnetic eld both parallel to the unperturbed magnetic eld B. Ak produces a magnetic eld predominantly perpendicular to B, B? = r ? Ak ; (4.6) where r? is the perpendicular gradient operator. Ak can therefore be identied as the vector potential associated with the shear Alfven mode. The parallel magnetic eld Bk obviously describes the compressional Alfven wave (magnetosonic wave) which should be taken into account when plasma low is not negligible. In a plasma, the magnetosonic perturbation is ignorable and low frequency electromagnetic modes can be described in terms of the two potentials, and Ak with su¢ cient accuracy. The Alfven mode is

characterized by magnetic eld line bending which causes particle drift vk This is still an E B? : B (4.7) B drift with a motional electric eld given by E0 = 1 vk c B? ; (4.8) which yields c E0 B B? = vk : 2 B B 2 (4.9) Source: http://www.doksinet Combining with the electrostatic E B drift, we thus nd the total perturbed drift, vD = c r B? + vk : B2 B B (4.10) After linearization, the Vlasov equation becomes e df + vD rf0 + dt m 1 E+ v c B @f0 = 0; @v (4.11) where as before d @ = + (vk + VD ) r; dt @t with VD the magnetic drift velocity. If the unperturbed distribution f0 is Maxwellian, this reduces to df + vD rfM dt e E (vk + VD )fM = 0; T (4.12) where E 1 @Ak : c @t r (4.13) Following the same procedure as developed in Chapter 3, we obtain f= with e fM + T ! (4.14) = k? v? = : If needed, the magnetosonic perturbation can be readily implemented in a similar manner, f= e fM + T ! n ! ! b ! b D kk vk vk Ak J02 c where the Bessel functions J0;1 have the

same argument 4.3 vk e Ak J02 ( ) fM ; c T ! ! b ! b D kk vk i o e v? J0 J1 A? fM ; c T (4.15) : Kinetic Ballooning Mode Equation The basic equations to govern low frequency modes are the charge neutrality condition ni = ne ; or Z fi dv = and parallel Ampere’s law, r2 Ak = 4 Jk = c 4 n0 e c 3 Z Z fe dv; vk (fi (4.16) fe )dv; (4.17) Source: http://www.doksinet where the perturbed ion distribution function fi can be found from Eq. (414), !+! b i J2 kk vk + ! b Di 0 e fM i + Ti ! fi = with ! b i (v 2 ) = cTi 1+ eB 2 M v2 2Ti i i 3 2 vk e Ak fM i ; c Ti (4.18) B0 ] k? ; (4.19) [r(ln n0 ) 1 2 v + vk2 (rB 2 ? cM eB 3 ! b Di (v) = k? v? B0 ) k? : (4.20) For the ballooning mode, we may assume that the mode frequency is much larger than the ion transit frequency j!j kk vT i : Then the ion density perturbation becomes electrostatic and the ion current parallel to the magnetic eld is ignorable, ni ( 1 + Ii ) e n0 ; Ti (4.21) where the function

Ii dened by Ii = Z !+! b i 2 J !+! b Di 0 k? v? fM i dv; (4.22) i involves ion kinetic resonance at ! + ! b Di (v) = 0: For the electrons, the nite Larmor radius e¤ect may be ignored. The perturbed electron distribution function can be written down analogously, e fM e Te fe = ! where ! b e v2 = cTe 1+ eB 2 ! b De (v) = cm eB 3 vk e Ak fM e c Te ! ! b e kk vk ! b De mv 2 2Te e 3 2 [r(ln n0 ) 1 2 v + vk2 (rB 2 ? B0 ] k? ; B0 ) k? : (4.23) (4.24) (4.25) In Eq. (423), e¤ects of trapped electrons are ignored They have relatively weak stabilizing in‡uence on the ballooning mode through a reduction in the electron parallel current. In the low frequency limit j!j ne = kk vT e ; the electron density perturbation can be approximated by * e n0 Te ! ! ! b e v2 !+! b De (v) kk vk ! e Ak ckk e n0 : Te vk Ak c + v e n0 Te (4.26) The parallel current is largely carried by the electrons, and can be evaluated from the 1st order moment of 4 Source:

http://www.doksinet the perturbed electron velocity distribution function, Jke = e Z vk fe dv n0 e2 (! kk Te (! !) + e ! e )(! ! De ) + ckk e ! e ! De Ak : (4.27) Substituting the ion and electron density perturbations into the charge neutrality condition, n e = ni ; and the parallel electron current into Ampere’s law, 4 Jk ; c r2? Ak = we obtain ( 1 + Ii ) e n0 = Ti ! ! e Ak ckk e n0 ; Te (4.28) and 4 n0 e2 (! c kk Te r2? Ak = e !) + (! ! e )(! These two equations form a closed set for the two unknowns, ! De ) + ckk e ! e ! De Ak : (4.29) and Ak . The parallel Ampere’s law can be rearranged as 2 kk k? kk Ak k2 = De ckk c2 (! e !) + (! ! e )(! ! De ) + ckk e ! e ! De Ak : (4.30) 2 where kDe = 4 n0 e2 =Te is the square of the electron Debye wavenumber. Eliminating Ak between Eqs (4.28) and (430) yields the following kinetic ballooning mode equation, 2 kk k? kk b + 2 kDe (! ! e )2 c2 1 + Ii (! ! De ) (! where b is a reduced scalar

potential dened by b = (1 + ! e) e ! e ! De b = 0; (4.31) Ii ) : (4.32) After ballooning transformation, Eq. (431) is converted into a di¤erential equation, d d [1 + (s where Ti = Te and sin )2 ] e = e + = have been assumed and the mode frequency is normalized by electron 4 n (1 + ) ( 1)[ 5 f ( )] + f ( ) ( 2 1)2 Ii ( ) d d = 0; (4.33) Source: http://www.doksinet diamagnetic drift frequency, ! ; ! D and (k? )2 = !=! e : In the MHD limit j!j 1; we readily recover the MHD ballooning mode equation analyzed in Chapter 2. It is noted that the safety factor q is absorbed in the ballooning parameter dened by = q2 R [ Ln i (1 + i) + (1 + e e )] ; and does not appear explicitly in the mode equation. This is because the ion acoustic transit e¤ect has been ignored by assuming ! kk cs : When this condition becomes marginal, an integral equation approach must be taken as will be discussed in Section 4.5 The ion integral, Ii ( ); involves double

integrations over v? and vk . Its explicit form is 4 Ii ( ) = p Z1 Z1 0 0 where x = v? =vT i ; y = vk =vT i (vT i = i + 1 + i (x2 + y 2 23 ) 2 hp 2 ( )x e J 0 1 2 + f ( ) 2 x + y2 p x2 y 2 xdxdy; (4.34) 2Ti =M ) are the normalized velocities, and ( ) is dened by, ( )=k p 1 + (s sin )2 : (4.35) The integral is to be evaluated numerically within a shooting code at every point in the ballooning space : Exact numerical evaluation of the integral is time consuming. An e¢ cient means to evaluate the velocity integral is the Gaussian-Hermite quadrature method which exploits the exponential function, e x2 y 2 , con- tained in the equilibrium Maxwellian velocity distribution function. The one-dimensional Gaussian-Hermite quadrature formula is, where fxi ; i = 1; Z 1 e x2 f (x)dx = 1 X Ai f (xi ); (4.36) i ; ng are the integral nodes and fAi ; i = 1; ; ng the weighting coe¢ cients. Both xi and Ai are independent of the integrand f (x). xi are the roots of the n th

order Hermite polynomial, Hn (x), and the coe¢ cients Ai are given by Ai = 2n 1 p 1)! (n 2n [Hn 2 1 (x)] : (4.37) Both xi and Ai are tabulated, or can be generated by the following recurrence relations, 8 > > H (x) > < n+1 H0 (x) > > > : H (x) 1 = 2xHn (x) = 1; = 2x: 2nHn 1 (x) ; : 6 (4.38) Source: http://www.doksinet The integral Ii ( ) in the ballooning mode equation involves two dimensional integrations of the form Z 1 Z 1 x2 y 2 f (x; y)e dxdy 0 0 XX i For the integral Ii ( ) it has been found that n Ai Aj f (xi ; yj ): (4.39) j 40 is su¢ cient to provide the required accuracy. An alternative is to nd an analytic approximation for the integral. Guided by the two ‡uid ion density perturbation in the long wavelength limit worked out in Chapter 3, it has been shown that the following expression provides a reasonable approximation for Ii ( ) which agrees qualitatively with the exact velocity integral, Ii ( ) = ! + 73 !

Di (! + ! i ) 2 ! + 53 ! Di i i ! Di ! i 10 2 9 ! Di ! + 53 ! Di ! Di ! + 53 ! Di 2 10 2 9 ! Di e e I0 ( ) [I0 ( ) I1 ( )] ; (4.40) 2 where ( ) = (k )2 [1 + (s sin ) ], and I0;1 are the zero-th and rst order modied Bessel function, respectively. The approximation in Eq (440) qualitatively agrees with the exact integral and may be useful for analytic stability analysis of low frequency electrostatic and electromagnetic modes. 4.4 Analysis of the Kinetic Ballooning Mode With the numerical techniques described in the preceding section, we now present the results of stability analysis of the kinetic ballooning mode. In Fig 41 (a) and (b), the growth rate and mode frequency normalized by the Alfven frequency , (k )2 = 0:01; i = e = 2; n =! A and ! r =! A (! A = VA =qR) are shown when s = 0:4; b0 = = Ln =R = 0:175: The growth rates predicted by the ideal MHD and two-‡uid analyses are also shown (dotted line). As far as the maximum growth rate is concerned, the

MHD and kinetic theories agree well. The critical for the onset of the ballooning mode from the kinetic theory is which is somewhat smaller than that obtained from the ideal MHD theory, c c = 0:35; = 0:39: The dashed line shows a second mode revealed by the kinetic analysis. The second stability regime predicted by the ideal MHD essentially disappears in kinetic analyses. The growth rate revealed from the kinetic analysis persists in the MHD second stability region. The kinetic ballooning mode in the MHD second stability regime requires a nite ion temperature gradient, i & 1, and is driven through the resonance contained in the non-adiabatic ion density perturbation, Ii ( ). At small shear, the critical for the (kinetic) ballooning mode becomes small and the instability becomes threshold-less and remains unstable at any : Fig. 42 shows the case s = 0:2 with other parameters unchanged from those in Fig. 41 7 Source: http://www.doksinet Figure 4-1: Growth rate =! A

and frequency ! r =! A of the kinetic ballooning mode vs. when s = 0:4; q = 2; b0 = 0:01; Ti = Te; Ln =R = 0:175; i = e = 2: The dotted line shows the growth rate of the ideal MHD ballooning mode and dashes lines show the second kinetc ballooning mode. 4.5 Finite E¤ects on the Drift Type Modes The toroidal drift mode and ion temperature gradient mode are predominantly electrostatic. In the analyses presented in Chapter 3, the magnetic perturbation Ak was ignored which, however, is justiable for discharges with negligible : The magnetic perturbation Ak is proportional to : This may be seen from the parallel Ampere’s law r2 A k = 4 Jk : c (4.41) Taking divergence of both sides and recalling the charge neutrality condition r J? + r Jk = 0; we obtain r r 2 Ak = 8 4 r J? ; c (4.42) Source: http://www.doksinet Figure 4-2: Same as Fig. 41 except s = 0:2: where the divergence of the cross eld current may be approximated by r J? For i n0 e2 Te [! + (1 + 2 i )! i ](k? s )

+ ! e (! De + ! Di ) ! : (4.43) well below the MHD ballooning limit, the ion polarization current is dominant and the parallel vector potential is related to the scalar potential through Ak !+! ckk ! pi i i 2 / : (4.44) Substituting this into the electron density perturbation, we nd ne = ! ! e Ak ckk e n0 = T 1 (! ! e )[! + (1 + (ckk )2 i )! i ] ! pi i 2 ! e n0 : Te (4.45) For modes having a frequency close to the electron diamagnetic frequency, electromagnetic corrections are small. For the long wavelength ion acoustic drift mode characterized by ! > ! e ; the electron density 9 Source: http://www.doksinet perturbation becomes smaller than that in the electrostatic limit, ne < e n0 ; Te (4.46) that is, electron response is less adiabatic which should be further destabilizing. In contrast, the toroidal mode is characterized by ! < 0 and j!j < (1 + i )! i : i Then, the electron response is more adiabatic in this case, ne > and

electromagnetic (nite ) corrections to the Finite e¤ects on the toroidal i i e n0 ; Te (4.47) mode are expected to be stabilizing. mode can still be analyzed by Eq. (433) as long as the ion transit frequency is negligible, j!j > kk vT i : In this limit, the ion dynamics remains electrostatic and electromagnetic e¤ects enter mainly through electron dynamics. Eq (433) describes both kinetic ballooning mode and mode corrected for nite i e¤ects. For the ballooning mode having a frequency j!j kk VA ; the condition is well satised relatively independent of the nite Larmor radius parameter k? : However, for the toroidal i mode, the eigenvalue ! scales with ! / k? and the condition j!j > kk vT i is satised only for comparatively short wavelengths. In contrast to the i mode, the ion acoustic drift mode is further destabilized by : Eq. (433) is, unfortunately, inapplicable to analyzing the ion acoustic mode because the mode frequency in this case is close to the ion

acoustic transit frequency kk cs kk vT i (Te Ti ) and the assumption ! down. A rigorous analysis on nite kk vT i breaks e¤ects on the ion acoustic mode requires integral equation formulation which has recently been developed in the the Plasma Physics Laboratory. Figure 4.3 shows stabilizing e¤ect of total conditions. In (a), Ln =R = 0:2; i = e (ballooning parameter) on the toroidal i mode in various = 2; in (b), Ln =R = 0:5 (nearly ‡at density prole), i = e = 4; and in (c), same condition as in (a) except trapped electrons are included (r=R = 0:2): Common parameters 2 are: (k ) = 0:1; s = 1; Ti = Te : Stabilization of the toroidal i mode occurs when the ballooning parameter exceeds a threshold which depends on discharge parameters. Trapped electrons have destabilizing e¤ect on the i mode. (This is in contrast to the case of the kinetic ballooning mode which tends to be stabilized by trapped electrons.) The critical e where " = r=R; e &

required for stabilization of the 2 1+ e p 3 (1 ") (1 + 2"n ) ( + 1) + = Te =Ti ; "n = Ln =R: 10 i 2 mode is approximately given by ; e (4.48) Source: http://www.doksinet Figure 4-3: Stabilization of the i mode by ; the ballooning parameter. Solid lines show growth rate =! e 2 and dashed lines mode frequency ! r =! e : s = 1; = Te =Ti = 1; b0 = (k ) = 0:1: (a) Ln =R = 0:2; i = e = 2; r=R = 0 (no trapped electrons). (b) Ln =R = 0:5; i = e = 4; r=R = 0: (c) Same as (a) except r=R = 0:2: 4.6 Stability Analysis based on Integral Equations When the mode frequency approaches either the ion transit frequency kk vT i or the electron transit frequency kk vT e ; the di¤erential formulation breaks down and formulation based on integral equations in the ballooning space must be used. ! kk vT i kk cs if Ti Te occurs in long wavelength regime (k ! kk vT e in short wavelength regime (k 2 i) 2 i) 1, and > 1: In this Section, derivation of integral equations

suitable for stability analysis of tokamaks will be outlined. Unfortunately, resultant integral equations can only be solved numerically. The perturbed distribution f derived in Section 4.1, f= e fM + T ! ! ! b ! b D kk vk vk e Ak J02 ( ) fM ; c T (4.49) is, strictly speaking, valid only if the operator nature of the parallel wavenumber kk is ignorable. If not, the 11 Source: http://www.doksinet nonadiabatic part h in f; e fM + h eik T f= = e fM + hJ0 (k? ) ; T (4.50) should be calculated by solving the gyro-kinetic equation for h; @"1 @t @h + VD0 rh + hVD1 i rfM + @t @fM = 0; @" (4.51) where VD0 is the zero-th order guiding center drift velocity, VD0 = vk + mc eB 3 1 2 v + vk2 B 2 ? rB; (4.52) VD1 is the perturbed guiding center drift due to wave motion, VD1 = c vk Ak B c c E? B =i 2 B2 B k; (4.53) and "1 is the energy perturbation, vk Ak : c "1 = e Note that the perturbed electric eld E? is E? h 1 (v B1 )? c i k? Ak + vk c vk Ak : c

= r? + = r? = ik? (4.54) i indicates gyro-phase averaging. For example, ik r 0e = ik rG 0e eik = ik rG J 0e (k ) ; where rG is the coordinates of the guiding center. With these preparations, Eq. (451) can now be cast into the form vk @h qR @ i[! ! D (v; )]h + i (! ! ) 12 vk e Ak J ( ) fM = 0; c T (4.55) Source: http://www.doksinet where ! D (v; ) mc 1 2 v + vk2 k (B rB) eB 3 2 ? mc 1 2 v + vk2 k [cos + (s eBR 2 ? = = ! =! 0 mv 2 2T 1+ 3 2 ; and k v? ( )= ! q 0 sin ) sin ] ; cT k (B eB 2 = r ln n0 ) = (4.56) cT k ; eBLn (4.57) 2 1 + (s sin ) : (4.58) Integration of Eq. (455) over can be performed separately for circulating (untrapped) particles and trapped particles. For circulating particles, we obtain vk > 0 : h+ C ( )= i vk < 0 : hC ( ) = Z i exp qR (! vk 0 i 1 Z 1 exp i where 0 = Z 0 where = 1 for vk > 0 and = 0 and 1 h ( 0) = A 0 00 !D d 00 ! ! 0 J0 ( 0 ) d J0 ( 0 ) d : 0 e fM ; T e fM

; T (4.59) (4.60) (4.61) 2; Z i exp i 0 0 ; (4.62) 0 1 for vk < 0; 0 and vk + Ak c ! ) qR ! vk For trapped particles with turning point angles hT ( ) = exp i ! ) qR (! vk 0 vk Ak c = = qR (! vk qR (! vk + A; e fM ; T (4.64) vk eAk fM ; c T (4.65) ! ) J0 ( 0 ) ! ) J0 ( 0 ) (4.63) is arbitrary angular location between the turning points, 1 < < 2: The boundary conditions are: h+ T ( 1 ) = hT ( 1 ) hT ( 1 ) ; (4.66) h+ T ( 2 ) = hT ( 2 ) hT ( 2 ) : (4.67) 13 Source: http://www.doksinet Since h+ T ( 1 ) = exp i h+ T ( 2) 0 i and hT ( 2 ) = exp i we nd 1 2 is arbitrary, we set it equal to 1: 1 Since 0 0 Z exp i 1 2 hT ( 1 ) = hT ( ) exp 2i exp i = Z 1 hT ( 1 ) + i 2 sin 21 Z i exp i 2 1 h exp i Z 2 exp i 0 + 1 0 d 0; (4.68) 0 (4.69) 1 Z 2 exp i 2 0 d 0; 1 0 2 0 + exp i 0 + 2 0 i d 0: (4.70) Then 2 1 h exp i 0 0 0 2 0 + exp i + 2 0 0 i d 0 d 0: (4.71) 1 Most

pressure gradient driven modes of interest fall in the frequency regime j!j ! be (the bounce frequency of trapped electrons ! be vT e =qR); and we may simplify hT ( ) for trapped electrons as hT ( ) where H 0 1 2 1 Z 2 1 n 0 +i h 2 0 0 H 2 1 i 0 A o d 0 +O ! ! be ; (4.72) is the Heaviside step function. Trapped ion modes pertain to the frequency regime well below the ion bounce frequency j!j ! bi and are ignored. Substitution of the perturbed distribution functions into the Poisson’s equation and parallel Ampere’s law yields the following integral equations, k 2 = 4 n0 X ej j 2 k? Ak ( ) = ej + Tj Z h i + gCj + gT+j J0 ( j ) dv Z h i 4 n0 X ej vk gCj + gT j J0 ( c j where gC(T )j = 1h + h 2 C(T )j 0 + hC(T )j 0 i : j ) dv; ; (4.73) (4.74) (4.75) An e¢ cient numerical algorithm to solve the set of integral equations has been developed by M. Elia (PhD thesis, 2000). Bench-marking studies with the previously found eigenvalues (! =

! r +i ), particularly those by Rewoldt and co-workers, have shown that eigenvalues in low beta regime agree. However, the growth rate of the ballooning mode found by Rewoldt seems overestimated. Fig 44 shows comparison among the eigenvalues 14 Source: http://www.doksinet from ve di¤erent methods, ideal MHD, Rewoldt’s comprehensive integral equation code, kinetic di¤erential equation method, semilocal kinetic code, and Elia-Hirose’s integral equation code (labelled “Nonlocal”) for discharge parameters pertinent to the Doublet III tokamak, r=R = 0:16; Ti = Te ; q = 1:14; e = 0:93; Ln =R = 0:07; when k i i = 0:87; = 0:28 (n = 20) : The growth rate found by Rewoldt is comparable with that of ideal MHD, while all other methods predict growth rates considerably smaller. Figure 4-4: Comparison among the eigenvalues found from 5 methods for discharge parameters pertinent to the Doublet III tokamak, r=R = 0:16; Ti = Te ; q = 1:14; i = 0:87; e = 0:93; Ln =R = 0:07; when k i

= 0:28 (n = 20) : 15