J.
6 (2004)
J. Plasma
Plasma Fusion
FusionRes.
Res.SERIES,SERIES,Vol.
Vol.
6 (2004)17–22
000–000

Gyrokinetic Simulations of Tokamak Micro-Turbulence
Including Kinetic Electron Effects
IDOMURA Yasuhiro, TOKUDA Shinji and KISHIMOTO Yasuaki
Department of Fusion Plasma Research, Naka Fusion Research Establishment,
Japan Atomic Energy Research Institute, Ibaraki, 311-0193, Japan(Received: 9 December 2003 / Accepted: 6 July 2004)

Abstract
A gyrokinetic toroidal particle code for a 3-dimensional nonlinear turbulence simulation (GT3D) has beendeveloped to study the ion temperature gradient driven – trapped electron mode (ITG-TEM) turbulence in tokamak
plasmas. From linear zonal flow damping tests and nonlinear ITG simulations, it is shown that a new method basedon a canonical Maxwellian distribution is essential to simulate correct zonal flow dynamics in tokamaks. Recently,
GT3D has been extended including kinetic trapped electrons. A computational cost of ITG-TEM calculations are
drastically reduced by using a new bounce-averaged kinetic trapped electron model. A short wavelength unstableregion of the ITG-TEM is calculated using a gyrokinetic field solver with a Pade approximation. From preliminary
linear ITG-TEM calculations, the validity of these calculation models are confirmed.
Keywords:
gyrokinetic simulation, ion temperature gradient driven mode, trapped electron mode,
electron temperature gradient driven mode

1. Introduction

ITG-TEM calculations. A gyrokinetic field solver which cantreat a short wavelength unstable region is developed byapplying a Pade approximation to the ion polarization density.
Linear properties of the ITG-TEM is successfully recoveredby using these calculation models. Adding trapped electrons
not only increase the growth rate of the ITG mode, but also
produce another unstable electron mode, the TEM. Between
these two modes, the dominant mode changes depending on
kq and hi (= Ln/Lti), where kq is the poloidal wave number,and Ln and Lti are scale lengths of density and ion temperature
gradients, respectively.
We then discuss zonal flow dynamics in the ITG turbulence. Unlike a conventional method based on a local
Maxwellian distribution FLM (y, e, m), which is defined by aflux label y, the new FCM method avoids spurious linear
driving effect on axisymmetric modes. Through zonal flow
damping tests and ITG turbulence simulations, it is foundthat this spurious driving effect significantly affects zonal
flow dynamics and spurious zonal flow oscillations are
excited.
The reminder of the paper is organized as follows. In
Sec. 2, model equations and numerical methods are described.
In Sec. 3, preliminary linear calculations including kinetic
trapped electrons are shown and the validity of new models
are discussed. In Sec. 4, zonal flow dynamics in the ITG

A global gyrokinetic toroidal particle code for a 3dimensional nonlinear turbulence simulation (GT3D) [1] hasbeen developed to study the ion temperature gradient driven– trapped electron mode (ITG-TEM) turbulence in tokamak
plasmas. The code has been developed based on a finite
element particle-in-cell (PIC) method [2], except for several
extensions for a nonlinear simulation. Main features of GT3D
are summarized as follows. First, we have developed a new
method based on a canonical Maxwellian distribution FCM(Pj, e, m) [1], which is defined by using three constants ofmotion in the axisymmetric toroidal system, the canonicalangular momentum Pj, the energy e, and the magnetic moment m. Second, the conservation property of GT3D is greatly
improved using the optimized particle loading [3]. Third, we
use a quasi-ballooning representation [4] for nonlinear
perturbations, which enables linear and nonlinear global highm, n calculations. Fourth, an analysis of the ITG-TEM isenabled by using kinetic trapped electron models. Finally, thecode has been optimized for massively parallel scalar andvector machines, and it operates with 40 % and 25 % processing efficiency up to 512 processors respectively on the
JAERI Origin3800 system and on the Earth Simulator [5,6].
In this paper, first, we discuss code developments towards
the ITG-TEM calculations. A new bounce-averaged kinetic
trapped electron model enables order of magnitude low cost
Corresponding author’s e-mail: idomuray@fusion.naka.jaeri.go.jp
117©2004 by The Japan Society of Plasma
Science and Nuclear Fusion ResearchIdomura Y. et al., Gyrokinetic Simulations of Tokamak Micro-Turbulence Including Kinetic Electron Effects

simulations based on FCM and FLM is discussed. Finally, asummary is given in Sec. 5.

and growth rate spectra of the ITG-TEM. This means that
passing electrons respond almost adiabatically to the ITGTEM with w /k|| << vte, where w is the frequency of the ITGTEM. Since the high frequency wH mode [8] does not appear
in the latter model, the computational cost of the latter modelis order of magnitude lower than the former model. However,the computational cost of the drift-kinetic trapped electron
model is still high, because the Courant-Friedlichs-Lewy(CFL) condition is limited by the electron bounce frequency
wbe, and the ballistic noise is inevitable in the drift-kinetic
model. In order to resolve these problems, we have developed
a bounce-averaged kinetic electron model based on the
bounce-averaged kinetic equation [9]. In the bounce-averaged
kinetic formalism, the fast bounce motion is eliminated byapplying the Lie perturbation theory under the ordering, w /wbe ~ O(e). In the deeply trapped approximation, k 2 << 1,the bounce-averaged coordinates (bb, ab, Y, I, m, a) are givenas

2. Calculation model
2.1 Gyrokinetic equations
The present version of GT3D solves the electrostatic
gyrokinetic Vlasov-Maxwell system in a circular concentric
tokamak configuration. In the gyro-averaged coordinates, Z= (t ; R, v ||, m, a), the electrostatic gyrokinetic VlasovMaxwell system [7] is written as
dFs Fs dR
dv F
=
+
 — R Fs + || s = 0,t
dtdt
dt v||(1)dR
c= v||b +bdt
qs B||*()¥ qs — R f a + ms v||2 b  — R b + mB— R ln B ,

(2)

dv||
=
dt
-

B*
 qs — R f a + ms v||2 b  — R b + mB— R ln B , (3)ms B||*()

bb  y ,(6)a b  j - qq ,

(7)
(8)

dns ( x, t ) = Ú dfs ( R, v|| , m , t )d ([ R + r s ] - x ) Ds d 6 Z-

ns q sTs-— 2f = 4p

Â [1 - I (k r ) exp(-k r )]f exp(ik  x), (4)
0

2 2
^ ts

2 2
^ ts

k

(9)

k

Â q dn (x),sswhere q and j are the poloidal and toroidal angle, respectively, q is the safety factor, k is the pitch angle, R0 is the–
major radius, er = r/R0, and B is a flux surface averaged
magnetic field. In this coordinates, the bounce-averaged
Hamiltonian H0 and the corresponding kinetic equation are
obtained as,(5)s

where R is a position of the guiding center; v|| = v · b ; v^ =| v ¥ b |; m  msv^2 /2B ; Ws = qsB/msc; ms and qs are the massand charge of the particle species s, respectively; c is thevelocity of light; b = B/B is the unit vector in the direction of
the magnetic field; a is the gyro-phase angle; B*= B + (Bv|| /Ws)—R ¥ b; B*|| = b · B*; f and fk are the electrostatic potentials
in the configuration and wave number spaces, respectively;·fÒa is the gyro-averaged electrostatic potential; F0 and d fs
are the equilibrium and perturbed parts of the distributionfunction Fs, respectively; ns and dns are the equilibrium andperturbed parts of the particle density; Ts is the temperature;
Ds is the Jacobian of the gyro-averaged coordinates; In is then-th order modified Bessel function; k|| = k · b and k^ = |k ¥b | are the parapllel and perpendicular component of the wave
vector k and rts is the Larmor radius evaluated with thethermal velocity vts.

H0 ( b b , I , m ) = mB (1 - e r ) +1
qR0

e r mBI,
ms(10)

dFb Fb db b Fb da b Fb=++= 0, (11)dt
dt b bdt a b
t

fY
db b
=c
,
a b
dt

(12)

fY
da b
c H0
=-c
,
b b
dt
qs b b(13)

where ·fÒY is the bounce-averaged electrostatic potential.
Since the action integral I is an exact invariant in the bounceaveraged coordinates, the Hamiltonian H0 does not dependon the corresponding conjugate variable Y, and the fast
bounce motion and the associated ballistic noise are eliminated from the equation system. Accordingly, the CFL
condition of the bounce-averaged kinetic electron model isalmost comparable with that of ITG calculations, and thenumerical convergence against the number of bounce centers

2.2 Trapped electron models
In principle, kinetic electrons can be described by taking
a drift-kinetic limit of the gyrokinetic equation, Eqs. (1) – (3).
In the present study, we have developed two drift-kinetic
electron models. One is a full drift-kinetic electron model andthe other is a drift-kinetic trapped electron model with adiabatic passing electrons. From linear benchmark tests, we haveconfirmed that both models give the almost same frequency
2
18

Idomura Y. et al., Gyrokinetic Simulations of Tokamak Micro-Turbulence Including Kinetic Electron Effects

is drastically improved compared with the drift-kinetic electron model. In the electron (zero banana width) limit, thetrapped particle density is given as

È dR
ddf
dv F0 = -Í — R F0 + ||,dt
dt 1 v|| Î dt 1

where dR/dt|1 and dv|| /dt|1 show the perturbed part of nonlinear characteristics. Eq. (17) is valid when F 0 is anequilibrium solution, which satisfies dF0 /dt = 0. In theaxisymmetric toroidal system, an exact equilibrium solution
is given by a canonical Maxwellian distribution, FCM (Pj, e,m). However, the conventional d f method was developed
based on a local Maxwellian distribution, FLM (y, e, m). Since
y is not a constant of motion in the axisymmetric toroidal
system, FLM does not satisfy dF0/dt = 0. In the conventionald f method, Eq. (17) has been used by ignoring a variation of
FLM along the unperturbed characteristics, dR/dt |0 · —RFLM +dv|| /dt |0 v|| FLM. Obviously, this treatment violates the conservation property of the gyrokinetic equation. In a new d fmethod, we use FCM in calculating the r.h.s. of Eq. (17). In
the flux coordinates (y, q, j, v||, m, a), Eq. (17) is written as

dnb (y ,q , j ) = Ú dfb ( b b (y ), a b (y ,q , j ), I ¢, m ¢)d (Y ¢ - Y (y ,q , j , I ¢)) Db dY ¢dI ¢dm ¢da ¢. (14)The bounce-average in Eqs. (12), (13) and (14) is calculatedusing a finite sampling point average along the bounce orbit
whose turning points are determined by I or the pitch angle kof each particles. Actually, 20 ~ 80 sampling points arechosen along the field line by following Dq = v|| /qR0D t witha constant D t.
Another extension of the code for TEM calculations is atreatment of the ion polarization density. When we consider
the ITG turbulence with adiabatic electrons, where a linear
unstable region is limited for kqrti < 0.6 (see Fig. 3), a Taylor
expansion for the ion polarization density (the second termin Eq. (4)), 1 – I0 (b) e–b ~ b, is a relatively good approximation. Under a long wavelength approximation, k^2 r 2ti << 1, Eqs.
(4) and (5) yield the gyrokinetic Poisson equation for the ITG
mode,È dR Ê y C
ˆ F
y C
ddf
—q ˜ CM
—y +
= -Í
Á

q

dt
d
t
y
Ë
¯ y C1Îdv Ê y C FCM FCM ˆ + || Á
+
v|| ˜¯  ,dt 1 Ë v|| y C

r2
1
- — ^  2ti — ^f + 2 (f - f q )
l Di
l De
= 4pqi Ú dfi ( R, v|| , m , t )d ([ R + ri ] - x ) Di d 6 Z,

3. Linear ITG-TEM calculations
Simulations have been performed using Cyclone basecase parameters [13]: deuterium plasma, R0 = 1.3 m, a = 0.48
m ~ 150rti, B0 = 1.9 T, Te = Ti, R0 /Lti = 6.92, R0 /Ln = 2.22 at
a reference surface or r0 = 0.5a. The q profile is given as q =0.84 + 2.18(r/a)2 (q(r0) = 1.4, s(r0) = 0.776).
Figure 1 shows the hi dependence of the growth rate andfrequency. In the frequency spectrum, it is shown that the
ITG-TEM has two independent branches, the ITG mode andthe TEM, which propagate in the ion and electron
diamagnetic directions, respectively [14]. In the adiabatic
electron case, a critical hi exists near hi ~ 2. On the other
hand, in the kinetic electron case, a critical hi disappears, andthe ITG-TEM is destabilized by trapped electrons even at hi= 0 [15]. In a transition region near hi ~ 2.5, both the ITG
mode and the TEM coexist. In the figure, the bounceaveraged kinetic electron case shows reasonably good agreement with the drift-kinetic trapped electron case. Differencesof the real frequency (~ 20 %) and the growth rate (~ 10 %)in a low hi region, where the TEM becomes dominant, may
be due to a deeply trapped assumption k 2 << 1 or a large
aspect ratio assumption r/R0 << 1 used in the present bounceaverage formalism.
Figure 2 shows the kq spectrum of the ITG-TEM. As in
the adiabatic electron case (see Fig. 3), the growth rate
spectrum of the ITG mode peaks at kqrti ~ 0.3, but the growth

cp
rti2
2
2 — ^ f + (1 + — ^  r ti — ^ ) 2 (f - f q )
l Di
l De
È



ÍÎ



= 4p (1 + — ^  rti2 — ^ ) Íqi Ú dfid ([ R + ri ] - x ) Di d 6 Z + qednt  ,(18)where yC = –(c/qs)Pj (y, q, v||). Since the new method is constructed straightforwardly based on FCM, it keeps the conservation property of the system.

(15)

where ·fÒq is the flux surface averaged potential and lDs isthe Debye length. On the other hand, in analyzing the ITGTEM turbulence, which has a broad unstable spectrum up toa short wavelength region or kq rti 0.6, we use a Padeapproximation for the ion polarization density, 1 – I0 (b)e–b ~b/(1 + b) [10]. By using a Pade approximation, the gyrokinetic field equation for the ITG-TEM turbulence is obtained as,- —^ 

(17)(16)

where cp is a fraction of passing electrons and dnt is thetrapped particle density. In the code, the gyrokinetic field
equation is solved using a finite element method [11] with afixed boundary condition f = 0 at r = a and a natural
boundary condition at r = 0.

2.3 New d f method based on canonicalMaxwellian
In a nonlinear characteristic d f method [12], d f, which
is defined as a deviation from F0, is solved numerically using
a nonlinear characteristic method. In this method, an evolution
equation of d f along nonlinear characteristics or perturbed
particle orbits is described by3
19

Idomura Y. et al., Gyrokinetic Simulations of Tokamak Micro-Turbulence Including Kinetic Electron Effects

Fig. 1

The hi dependences of growth rate and frequency are plotted for the ITG-TEM (drift-kinetic trapped electrons:
open circle, bounce-averaged kinetic trapped electrons: open square) and the ITG mode (adiabatic electrons:
cross) in Cyclone base case plasmas kq rti = 0.275).

Fig. 2

The growth rate and frequency spectra of the ITG-TEM in Cyclone base case plasmas (R0 /Lti = 6.92, hi = 3.12).
Results obtained from gyrokinetic Poisson field solvers with a long wavelength approximation (Taylor) andwith a Pade approximation are compared.

rate itself is enhanced by trapped electrons which do notrespond to the ITG mode adiabatically. In a high kq region,
the ITG mode is stabilized by the finite Larmor radius (FLR)
effect, and the TEM becomes dominant. This high kq region
with kqrti > 0.6 can not be analyzed by using a gyrokinetic
Poisson equation with a long wavelength approximation,(kqrti)2 << 1. On the other hand, a gyrokinetic field solver witha Pade approximation successfully captures the high kq ITGTEM. It is noted that a linear ITG-TEM benchmark among
GT3D, GTC (a global gyrokinetic PIC code), and FULL (alocal gyrokinetic ballooning code) is now going on, and thesethree codes show reasonably good agreements for the linear
frequency and growth rate spectra [16].

 f a FCM
 f a FCMddf = q-q,
dt 0
j Pjte(19)

there is no linear driving effect on axisymmetric modes,
because the driving term associated with FCM /Pj (the pressure gradient in the Pj direction) disappears. It is noted thatin the zero orbit width limit, the Pj direction becomes theradial direction or the y direction. But, if we replace FCM(Pj, e, m) by FLM (y (Pj, q, e), e, m), we haveÈ  f a y
 f a y  FLMddf = Íq-qt e  ydt 0ÍÎ j Pj

4. Zonal flow dynamics in plasmas withcanonical Maxwellian distribution-q

In most of conventional gyrokinetic simulations, a local
Maxwellian distribution FLM (y, e, m) has been used as an
approximate equilibrium solution. Besides the problem related
to the conservation property of the d f scheme, a gyrokinetic
simulation based on F LM has the following problem.
According to the linear gyrokinetic equation in the canonicalcoordinates, (Pj, q, j, e, m, a) [17], f a FLM
t

e

.

(20)

In a FLM model, the driving term associated with FLM / y(the radial pressure gradient) can drive axisymmetric modes
through the second term in Eq. (20). For finite toroidal mode
number (n  0) modes, this spurious driving effect is just ahigher order minor correction compared with the dominant
driving effect in the first term, and the linear theory of micro4
20

Idomura Y. et al., Gyrokinetic Simulations of Tokamak Micro-Turbulence Including Kinetic Electron Effects

Fig. 3The eigenfrequency spectra of the ITG mode with adiabatic electrons are calculated for Cyclone base caseplasmas with FCM and FLM .

Fig. 5

Fig. 4

Zonal flow damping tests in Cyclone base case plasmas
with FCM and FLM . In the FCM case, the residual flow levelagrees well with the theory.

The fluctuation field energy in ITG turbulence simulations with FCM and FLM . Cyclone base case plasmas are
simulated with Np = 4 ¥ 107 marker particles. Eturbulence and En = 0, m = 0 show the turbulent and zonal flow components,
respectively (Etotal = Eturbulence + En = 0, m = 0 ). In the nonlinear phase, spurious zonal flow oscillations grow in the FLM
case.

instabilities has been developed successfully based on a FLM
model. However, for the axisymmetric perturbations, where
the dominant driving term disappears, the spurious driving
term may be more significant.
In order to examine these points, we have performed
linear calculations of the ITG mode and zonal flow damping
tests [18] in Cyclone base case plasmas with FCM and FLM. In
Fig. 3, both the FCM and FLM cases give similar growth rate
and frequency spectra, and both results agree well with theprevious linear benchmark calculations [13]. On the other
hand, in zonal flow damping tests (see Fig. 4), the FCM andFLM cases show quite different behaviors. In the FCM case,zonal flows are damped rapidly with m = 1 damping oscillations and the residual zonal flow level agrees well with thetheoretical prediction [18]. However, in the FLM case, spurious
zonal flow oscillations are excited. This result is inconsistent
with the linear gyrokinetic theory, which predicts no driving
effect on axisymmetric perturbations including zonal flows.
Figure 5 shows the time history of the fluctuation fieldenergy in ITG turbulence simulations. In the FLM case, zonal
flow oscillations grow after the nonlinear saturation. On theother hand, in the FCM case, the zonal flow energy keeps aquasi-steady state. These zonal flow oscillations are often
discussed as the geodesic acoustic mode (GAM). Since zonal
flows are linearly stable, some nonlinear turbulent drive isneeded to excite GAM. However, according to the saturation
amplitude and the turbulent spectra in the initial saturation
phase, such a nonlinear driving effect is expected to be almost
the same in both cases. Thus, it is considered that zonal flow
oscillations are excited by a spurious linear driving effect on
zonal flows in plasmas with FLM.

5. Summary
Through zonal flow damping tests and ITG simulations,a spurious driving effect of FLM on zonal flows are identified.
This result shows an importance of a new FCM model in gyrokinetic turbulence simulations where zonal flows play asignificant role. For the ITG-TEM calculations, GT3D isextended including kinetic trapped electrons. A computational
cost of ITG-TEM calculations is drastically reduced by a
newly developed bounce-averaged kinetic trapped electron
model. A short wavelength unstable region of the ITG-TEM
5
21

Idomura Y. et al., Gyrokinetic Simulations of Tokamak Micro-Turbulence Including Kinetic Electron Effects[4] R. Gruber, S. Semenzato, F. Troyon et al., Comput. Phys.
Commun. 24, 363 (1981).
[5] Y. Idomura, M. Adachi, K. Gorai et al., J. Plasma Fusion
Res. 79, 172 (2003).
[6] Y. Kishimoto, Annual Report of the Earth Simulator
Center April 2002 – March 2003, 201 (2003).
[7] T.S. Hahm, Phys. Fluids 31, 2670 (1988).
[8] W.W. Lee, J. Comput. Phys. 72, 243 (1987).
[9] B.H. Fong and T.S. Hahm, Phys. Plasmas 6, 188 (1999).
[10] G. Hammett, W. Dorland and F.W. Perkins, Phys. Fluids
B 4, 2052 (1992).
[11] A. Bottino, T.M. Tran, O. Sauter et al., Theory of Fusion
Plasmas, Int. Workshop, (2001).
[12] S.E. Parker and W.W. Lee, Phys. Fluids B 5, 77 (1993).
[13] A.M. Dimits, G. Bateman, M.A. Beer et al., Phys.
Plasmas 7, 969 (2000).
[14] J.Q. Dong, S.M. Mahajan and W. Horton, Phys. Plasmas
4, 755 (1997).
[15] G. Rewoldt and W.M. Tang, Phys. Fluids B 2, 318(1990).
[16] G. Rewoldt, Z. Lin and Y. Idomura, Bull. Am. Phys. Soc.
48, 219 (2003).
[17] S. Wang, Phys. Rev. E 64, 056404 (2001).
[18] M.N. Rosenbluth and F.L. Hinton, Phys. Rev. Lett. 80,
724 (1998).

is successfully calculated using a gyrokinetic field solver witha Pade approximation. The validity of these models are confirmed through the preliminary linear ITG-TEM calculations.
The results show that trapped electrons drastically change ionscale micro-instabilities. The ITG-TEM turbulence simulation
is essential for a comprehensive study of the anomalous ion/
electron heat transport and the anomalous particle transport.
In a future work, nonlinear ITG-TEM simulations will be
addressed.

Acknowledgments
One of the authors (Y.I.) would like to acknowledge G.
Rewoldt, Z. Lin, A. Bottino, and T.S. Hahm for their useful
discussions on ITG-TEM calculations. The authors thank M.
Kikuchi, H. Ninomiya, and W.M. Tang for their support. The
simulations were performed on the JAERI Origin3800
system.

References[1] Y. Idomura, S. Tokuda and Y. Kishimoto, Nucl. Fusion
43, 234 (2003).
[2] M. Fivaz, S. Brunner, G. Ridder et al., Comput. Phys.
Commun. 111, 27 (1998).
[3] R. Hatzky, T.M. Tran, A. Könies et al., Phys. Plasmas 9,898 (2002).

6
22