J. Plasma
J.
Plasma Fusion
Fusion Res.
Res. SERIES,SERIES,Vol.
Vol.6 6(2004)(2004)23–27
000–000

Plasma Turbulent Transport and Fine-Scale Structures of Phase
Space Distribution Function
WATANABE Tomo-Hiko and SUGAMA Hideo
National Institute for Fusion Science / The Graduate University for Advanced Studies, Toki, 509-5292, Japan(Received: 12 December 2003 / Accepted: 7 April 2004)

Abstract
Recent results from kinetic simulations of steady and quasisteady states of the ion temperature gradient (ITG)
driven turbulence are reviewed in focus on the fine-scale structures of the phase space distribution function associatedwith the anomalous transport. Importance of proper treatment for the fine structures in numerical simulations of thecollisionless ITG turbulence is emphasized by a parameter survey for the velocity space resolution. Preliminary resultsof the flux tube simulation of the toroidal ITG mode are also given for the linear benchmark test and the collisionlessdamping in a tokamak configuration.
Keywords:
gyrokinetics, turbulence, transport, computer simulation

1. Introduction
the time-reversal symmetry causes the steady transport flux
with no dissipation mechanism. If the steady transport driven
by constant density and/or temperature gradients is observed
through coarse-graining of f with fine-scale structures, aquasisteady state of the collisionless turbulence should be
realized, where high-order velocity-space moments of fcontinue to grow while keeping low-order ones constant inaverage [3,4].
In comparison of collisionless turbulence simulation
results with actual plasma experiments, it is implicitly assumed that an asymptotic behavior of the transport coefficientin a low-collisionality limit agrees with the collisionless one
in the quasisteady state. The finite collisionality is indispensable to realizing the real statistically steady states for allorder-moments of f. The entropy variable associated withfluctuations is produced by the turbulent transport on a macro
velocity scale corresponding to the low-order moments. Then,
it is transfered in the phase space by the E ¥ B nonlinearity
and the phase mixing, and is, finally, damped by collisions
acting on a micro scale associated with the high-ordermoments.
In order to confirm the critical issues on collisionless andweakly-collisional turbulence given in the above twoparagraphs, we have investigated the slab ion temperature
gradient (ITG) driven turbulence by means of the Eulerian(Vlasov) kinetic simulations with high velocity-spaceresolution [5-7]. Results obtained in our previous works arereviewed in the next section with the aim of introducing the

The plasma turbulent transport has long been one of thecentral subjects in the magnetic confinement fusion research.
Since the collision frequency is much lower than characteristic ones of the turbulence in a core region of tokamaks withhigh temperature, in order to study the turbulent transport,
one needs to deal with a kinetic description of the magnetizedplasma such as the gyrokinetic equation. Numericalsimulations of the gyrokinetic plasma turbulence have rapidly
developed in the last decade [1], and revealed several
important aspects of the anomalous transport. For example,
turbulence suppression by the self-generated zonal flow has
been investigated by many authors [2], where interactions of
low-order velocity space moments of the one-body distribution function f and the electromagnetic field have been
discussed intensively. In the collisionless or weakly-collisional
plasma, however, a velocity space profile of f should havemore degrees of freedoms than those characterized by the
low-order moments such as the density, the flow velocity, thetemperature and so on. Even though a transport flux is definedby a correlation of the low-order moments and the electromagnetic field, fine-scale structures of f should be properly
treated in numerical simulations.
The one-body velocity distribution function f in acollisionless plasma turbulence is stretched and folded by the
Hamiltonian flow in the phase space while preserving its
amplitude. This is known as the phase mixing which leads tocontinuous generation of fine-scale fluctuations of f, and isclosely related to the paradoxical argument how thecollisionless turbulence described by the kinetic equation with

©2004 by The Japan Society of Plasma
Science and Nuclear Fusion Research

Corresponding author’s e-mail: watanabe.tomohiko@nifs.ac.jp
123Watanabe T. and Sugama H., Plasma Turbulent Transport and Fine-Scale Structures of Phase Space Distribution Function

background of the simulation studies given in the latter partof this paper. Artificial effects of coarse numerical grid on
the transport coefficient are shown in Sec. 3, which
demonstrates importance of the velocity-space resolution insimulations of the kinetic plasma turbulence. Our recent
progress in the toroidal ITG simulation is also reported inSec. 4. A summary is given in Sec. 5.

term is calculated by the pseudo-spectral method with the
3/2-rule for de-aliasing.

2.2 Steady and quasisteady states of ITG
turbulence
From Eqs. (1) and (2), one finds a balance equation ofentropy variable defined by a functional, d S = Sk Ú dv|| | f˜k |2/
2FM, that is,

2. A review of the slab ITG simulation by
the Eulerian kinetic coded(dS + W ) = hi Qi + D,
dt

2.1 Simulation model
We consider the following simulation model which isdefined in a two-dimensional slab geometry with the translational symmetry in the z-direction. The periodic boundary
condition is used in the x and y directions. The uniform
magnetic field is set in the y-z plane such that B = B (ẑ + q ŷ)for q << 1. The governing equations are derived from a v^integral of the electrostatic gyrokinetic equations [8] byignoring the parallel nonlinear term and by assuming d fk (v|| ,v^ ) = f˜k (v|| ) FM (v^ ), where FM is the Maxwellian velocity
distribution. We also assume constant density and temperature
gradients of the background ions in the x-direction with scale–1–1
lengths of Ln  –d (ln n)/dx and LT  –d (ln Ti)/dx, where Ln,
LT >> L (L denotes the perpendicular system size). Therefore,
we arrive at the following equations given in the wave
number space k = (kx, ky), such that

where, Qi, W, and D mean the perpendicular ion heat flux,
the potential energy, and the collisional dissipation,
respectively. Their explicit definitions are as follows;2k-k 2 /2

Ú f˜ (v )dv - n˜ ,
k

||

||

W=

f k )v||2 f˜ - k / 2,

Â (2 - G ) | f | /2,
2

0

k

and

D=

Â Ú dv f˜ C ( f˜ ) / F
||

-k

i

M

k

k

= -n

ÂÚ dv |  f˜ + v f˜ | / F < 0.
2

||

v||

k

||

k

M

k

Since d S is rewritten as d S = SM – Sm within the second orderfor d f = f – FM (where SM = – Úd3v FM ln FM and Sm = –· Ú d3v
f ln f Ò represent macroscopic and microscopic entropies,
respectively), we call it the entropy variable, hereafter. Here,· ···Ò means the spatial average.
Existence of the quasisteady state in the collisionless ITG
driven turbulence is confirmed by our Eulerian kinetic
simulation with high resolution for the velocity space [5],where monotonic increase of d S balances with the transport
in a saturated turbulent flow, such that2

where the electric potential fk is related to Yk by Yk = e–k /2 fk
with k 2 = k x2 + k y2. Also, Q = q L n/ ri and hi = Ln/L T. The
background electron temperature Te = Ti and the adiabatic
electron response are also assumed. Here, we consider alimiting case with no zonal flow component of ky = 0 by fixing
f˜ky =0 = fky =0 = 0. The instability drive is contained in the firstgroup of terms on the right-hand side of Eq. (1). The last termon the right-hand side denotes the ion-ion collision term forwhich we employ the Lenard-Bernstein model collision
operator,

Ci ( f˜ k ) = n v|| [ v|| + v|| ] f˜ k (v|| ),-k 2 /2

y

k

(2)

e, k

||

k

and
0

ÂÚ dv (-ik e

Qi =

(1)

[1 - G (k )]f = e(4)d(dS ) / dt ª hi Qi .

(5)

It means that d f never reaches to the steady state in thecollisionless case, even if a steady transport flux defined by
the low-order moments is observed. Continuous generationof fine-scale structures of d f in the velocity space, due to thephase mixing caused by the parallel advection term on theleft-hand-side of Eq. (1), is responsible for the increase ofd S.
In the weakly-collisional case, not only low-order but
also all-order moments of the distribution function can bestatistically steady. This is because the fine-scale structuresof d f in the velocity space are smoothed out by the collision
operator with the second-order derivative. It is also numerically confirmed that d (d S)/dt ª 0 and(3)with the collision frequency n normalized by vti/Ln. Morecomplete descriptions of the simulation model as well as thenormalization are given in Refs. [5,6].
Equations (1) and (2) are numerically time-integrated bymeans of the nondissipative integrator for the collisionlesscase [9-11] or the fourth-order Runge-Kutta-Gill (RKG)
method for the weakly collisional ones. The E ¥ B advection

hi Qi ª - D > 0

2
24

(6)

Watanabe T. and Sugama H., Plasma Turbulent Transport and Fine-Scale Structures of Phase Space Distribution Functionin the real statistically steady state of the slab ITG turbulence[6].

simulations of collisionless plasma turbulence to keep sufficient resolution to the fine-scale structures of d f generated
by the phase mixing. We have observed artificially enhanced
transport flux in cases with a coarser grid in the velocity space
as shown below. For comparison, we have carried out thecollisionless slab ITG turbulence simulations for Dv / vti = 5/
16, 5/32, 5/64, and 5/256 with the same physical and othernumerical parameters as given in Ref. [5]. Time-evolutionsof the transport coefficient ci are plotted in Fig. 1 (top andmiddle) for different D v. As typically seen in cases for Dv /vti = 5/16 and 5/32, ci gradually increases in the turbulence
after t ~ 300 or 400. The higher values of ci averaged from t= 800 to 1000 are observed for larger D v in the bottom of
Fig. 1. For D v / vti = 5/16 and 5/32, ci is twice larger thanthat reported in Ref. [5] where the sufficient resolution was
kept during the simulation (Dv / vti = 5/4096). As D v becomes
smaller, the time-averaged ci approaches to its proper value(shown by the horizontal dashed line) obtained in thequasisteady state of the collisionless turbulence. This isbecause the finest scale of d f in the velocity space reaches to
the grid scale at earlier time for larger D v. Thus, the result
from the simulation with a coarser grid suffers from the larger
aliasing error. The obtained result demonstrates theimportance of keeping enough velocity-space resolution in thecollisionless ITG turbulence simulation. In order to avoid thealiasing error, one should stop running the collisionlesssimulation before the finest scale-length of d f becomescomparable to the grid size. Then, if the collisionlessturbulence had not reached to the quasisteady state yet, it isnecessary to employ the finer resolution.

2.3 Collision frequency dependence
We have also investigated collision frequency dependence of the ion thermal transport coefficient ci defined as ci
 Qi / hi. A parameter survey for a wide range of n shows alogarithmic dependence of the transport coefficient on
relatively large values of n while ci approaches a value in thecollisionless case for sufficiently low collision frequency. Thisconclusion agrees with a concept that the quasisteady state is
an idealization of the real steady state in the weak-collisionality limit [12].
A spectral analysis of d f in the velocity space has also
been made by means of the Hermite-polynomial expansion[6,7]. The obtained spectrum has almost the same values at n
= 1, 2, and 3 even for different values of n (where n denotes
the order of the Hermite-polynomials), which agrees with cibeing independent of n in the low collisionality cases. This isalso consistent with the conjecture by Krommes and Hu that‘flux determines dissipation’ [3]. The entropy variable
produced in the macro velocity scale on the low-n side istransfered toward the micro scale with high-n values and thelarge wave numbers (k) by the phase mixing and the E ¥ Bnonlinearity. Then, it is dissipated by the collision term on
the high-n side.

2.4 Methodological issues in numerical simulations of the collisionless plasma turbulence
In order to simulate the collisionless plasma dynamics
with the time-reversibility, we have developed an Eulerian
simulation method based on the symplectic integrator [9-11].
The nondissipative simulation scheme for solving Eq. (1) isdesigned so as to preserve the time-reversibility and thephase-space integral of f 2 that is an invariant of the Vlasov
and the collisionless drift (gyro-) kinetic equations. The
numerical accuracy in dealing with the time-reversible plasma
dynamics including the phase mixing process has beenverified by the plasma echo simulation [11] as well as the 3mode ITG problem [9,10]. The sufficient velocity-spaceresolution to the fine-scale structures of d f is indispensable
to keep soundness of the numerical simulation, which isdemonstrated in the next section by comparing the resultswith different grid spacings in the velocity space.
The kinetic-fluid closure model is based on existence ofthe quasisteady state of the collisionless plasma turbulence.
The kinetic simulation result for the collisionless ITG turbulence agrees well with the collisionless fluid one with thenondissipative closure model (NCM) [12,13]. By employing
the NCM in the fluid equations, the time-reversibility of
unstable linear kinetic eigenmodes can be preserved, whichis essential to fluid simulations not only of the 3-mode ITG
problem [12] but also of the turbulent transport [13].

4. Development of a toroidal ITG simulation code for a flux tube geometry
We are developing the gyrokinetic-Vlasov simulation
code for the toroidal ITG mode in a flux tube geometry for a
tokamak configuration. The coordinate system and theboundary condition are based on the work by Beer [14],where concentric circular magnetic surfaces with a large
aspect ratio are assumed. We consider the gyrokinetic equation for ions [15,16], the adiabatic electron response, and thequasi-neutrality. As the first trial, we deal with passing ionsonly, while neglecting the mirror force term. Spatial
coordinates in radial (x) and field-line-label (y) directions are
discretized by the Fourier expansion, while the parallel (z)derivatives are approximated by the fifth-order upwind finite
difference. The parallel velocity (v ||) and the magnetic
moment (m) are chosen for the velocity space coordinates
which are discretized by grid points. The RKG method is used
for the time-integration. The simulation code is welloptimized in order to achieve high efficiency for vector andparallel operations.
The linear growth rate of the toroidal ITG modes for the
Cyclone DIII-D base case parameters [1] are shown in Fig.
2, where the solid and dashed lines represent the real andimaginary parts (denoted by wr and g ) of the eigenfrequency
obtained by the linear gyrokinetic code [17], respectively.

3. Effect of velocity space resolution
As summarized in the above, it is important in kinetic
3
25

Watanabe T. and Sugama H., Plasma Turbulent Transport and Fine-Scale Structures of Phase Space Distribution Function

Fig. 1

Fig. 2

Real frequency (wr) and linear growth rate (g ) of thetoroidal ITG modes for the Cyclone DIII-D base caseobtained by the gyrokinetic-Vlasov simulation code.
Solid and dashes lines indicate wr and g obtained by thelinear gyrokinetic code [17].

Fig. 3Collisionless damping of the density perturbation in atokamak configuration in absence of the electric field(solid). Dashed line represents the asymptotic behavior(µ t –2) predicted in Ref. [17].

[17], since not only the parallel advection term but also thetoroidal magnetic drift terms contribute to generation of finescale structures of the distribution function in the phase space.
The collisionless damping process can be successfully
reproduced by our simulation as shown in Fig. 3, where afiner numerical grid for the (v||, m)-space is employed, suchas (±192, 64) grid points, in order to continue the run up to t
= 100 Ln/vti. In lack of the resolution, otherwise, ñ unphysically grows at earlier time. The result demonstrates that, also
in a tokamak configuration, treatment of the fine-scale
structures of the distribution function is one of the key issues
for simulating the collisionless damping.

(top and middle) Time-evolutions of the ion thermaltransport coefficient (ci) obtained by the collisionless
ITG turbulence simulations for Dv / vti = 5/16, 5/32, 5/64,
and 5/256. (bottom) Time-averaged ci from t = 800 to1000 for different Dv in comparison to the case with thehigh velocity-space resolution of Dv / vti = 5/4096 (horizontal dashed line).

Here, we have employed (84, ±20, ±32, ±64, 32) modes/gridpoints in the five-dimensional (kx, ky, z, v||, m)-space, where kx
and ky denote the wave numbers in the x- and y-directions,
respectively. Solid squares and open circles indicate wr and ggiven by the gyrokinetic-Vlasov simulation results which
agree well with the linear code prediction.
In the absence of the electric field, the initial density
perturbation ñ with the ballooning type mode structure isdamped due to the phase mixing associated with the toroidal
particle drift. Its asymptotic behavior is proportional to t –2

5. Conclusions
In this paper, we have reviewed our recent progress inkinetic simulations of the steady and quasisteady states of theslab ITG turbulent transport. Even though the transport flux
is directly related to the correlation between the low-ordermoments and the electric field, high-order moments
associated with the fine-scale structures of the distribution
426

Watanabe T. and Sugama H., Plasma Turbulent Transport and Fine-Scale Structures of Phase Space Distribution Function[6] T.-H. Watanabe and H. Sugama, Phys. Plasmas 11, 1476
(2004).
[7] H. Sugama and T.-H. Watanabe, Spectrum of the Velocity
Distribution Function in the Slab Ion Temperature Gradient Driven Turbulence, presented in this conference.
[8] D.H.E. Dubin, J.A. Krommes, C. Oberman and W.W.
Lee, Phys. Fluids 26, 3524 (1983).
[9] T.-H. Watanabe, H. Sugama and T. Sato, Phys. Plasmas
7, 984 (2000).
[10] T.-H. Watanabe, H. Sugama and T. Sato, J. Phys. Soc.
Jpn. 70, 3565 (2001).
[11] T.-H. Watanabe and H. Sugama, to appear in Transp.
Theo. Stat. Phys. (2005).
[12] H. Sugama, T.-H. Watanabe and W. Horton, Phys.
Plasmas 8, 2617 (2001).
[13] H. Sugama, T.-H. Watanabe and W. Horton, Phys.
Plasmas 10, 726 (2003).
[14] M.A. Beer, S.C. Cowley and G.W. Hammett, Phys.
Plasmas 2, 2687 (1995); M.A. Beer, Ph.D thesis,
Princeton University, 1994.
[15] E.A. Frieman and L. Chen, Phys. Fluids 25, 502 (1982).
[16] T.S. Hahm, Phys. Fluids 31, 2670 (1988).
[17] H. Sugama, Phys. Plasmas 6, 3527 (1999).

function should also be taken into account accurately. This isbecause they are related with each other through the transferof the entropy variable in the phase space. In numericalsimulations of the kinetic plasma turbulence, thus, it isessential to keep enough resolution for the fine structures
which are spontaneously generated by the phase mixing in
the velocity space. Lack of the resolution may lead to artificially enhanced transport flux through the aliasing error as
demonstrated in this paper. In addition, taking into accountof the time-reversibility of the collisionless kinetic equationis also a key issue in numerical studies of the collisionlessturbulence, which was a fundamental motivation in development of the nondissipative closure model [12,13]. In a toroidal
configuration, the phase mixing by the toroidal magnetic drift
and the parallel advection terms generates the fine-scale
structures. A toroidal version of our Eulerian kinetic simulation code is successfully applied to the collisionless damping
in a tokamak configuration. Simulations of the toroidal ITG
turbulent transport is currently in progress and will bereported elsewhere.

References[1] A.M. Dimits, G. Bateman, M.A. Beer et al., Phys.
Plasmas 7, 969 (2000).
[2] P.W. Terry, Rev. Mod. Phys. 72, 109 (2000).
[3] J.A. Krommes and G. Hu, Phys. Plasmas 1, 3211
(1994).
[4] H. Sugama, M. Okamoto, W. Horton and M. Wakatani,
Phys. Plasmas 3, 2379 (1996).
[5] T.-H. Watanabe and H. Sugama, Phys. Plasmas 9, 3659
(2002).

Acknowledgment
This work is supported by the Japanese Ministry of
Education, Culture, Sports, Science, and Technology, Grant
Nos. 12680497 and 14780387.

527