J. Plasma Fusion Res. SERIES, Vol. 9 (2010)

Properties of freely decaying and driven turbulence of fusion plasmas using
gyrokinetic particle simulation
R. Ganesh∗
Institute for Plasma Research, Bhat Village, Gandhinagar - 382 428, India

W. W. LEE, S. ETHIER and J. MANICKAM
Princeton Plasma Physics Laboratory, Princeton, U.S(Received: 20 November 2009 / Accepted: 6 May 2010)
In a recent work, the eﬀects of nonlinear self-consistent proﬁle relaxation on the ITG turbulencein global delta-f particle simulations were studied [1] using Gyrokinetic Toroidal Code (GTC) [2]. In
contrast to several other gyrokinetic turbulence models, GTC shows that the turbulence is insensitive
to minute or drastic proﬁle changes and reaches the similar ﬂux levels irrespective of whether or not
it is a driven or freely decaying turbulence. In the present work, we address the spectral properties
of gyrokinetic turbulence and argue that correct inclusion of parallel velocity space nonlinearity
may bring forth certain important late-time features such as Geodesic Acoustic Modes like (GAM)
structures.
PACS numbers:
Keywords: gyrokinetics, plasma turbulence, simulation

I.

tial value problem leads to a new steady state with lower
ﬂux values. In this paper, we shall address the issues of
proﬁle relaxation, eﬀect of velocity space nonlinearity on
the spectral features of turbulence and its role in formation of Geodesic Acoustic Modes like (GAMs) structures.

INTRODUCTION

After the discovery of the elongated radial structuresassociated with the ITG modes (streamers) [3], and theestablishment of the importance of the E × B nonlinearity for the spatial trapping and de-trapping of the resonant particles [4, 5], powerful computing capabilities
have demonstrated to us the importance of the generation of zonal ﬂows [2], and the energy cascade throughmode-coupling processes in toroidal geometry [6]. However, it is not until recently, with the availability of enormous massively parallel computational power, that wehave the opportunity to examine the role played by the
global zonal ﬂow modes and the nonlinearity associatedwith the acceleration of the particles in velocity space dueto the perturbed ﬁeld for large tokamaks [7]. It is foundthat this particular nonlinearity, which has now drawn
widespread attention in the fusion community [8], plays
an important role in the generation of zonal ﬂows andthe ensuing nonlinear ion thermal transport, especially,
when we increase the simulation volume. In fact, it isfound that this nonlinearity is also responsible for thedamping of the zonal ﬂows. Without it, the simulation
would fail to converge. Another important question iswhether this resulting ion transport is a transient associated with a decay instability or it is indeed some type of
steady state phenomena prescribed by the initial proﬁles.
A recent study by taking into account the quasilinear proﬁle relaxation in the simulation has shown that the proﬁlerelaxation, which would aﬀect linear stability properties,
has negligible eﬀects on the nonlinear ion thermal transport [9]. In this study it is also shown that by removing
the wave-particle interaction eﬀects and evolving the latetime steady state density and temperature proﬁle as ini-II.

THE GOVERNING EQUATIONS

The governing equations of the simulation are basedon the gyrokinetic Vlasov-Poisson equations [11–13]. In
terms of the gyrokinetic units of Ωi (≡ eB/mi c), ρs (≡�
Te /Ti ρi ) for time and space, and eφ/Te for the perturbed potential, where ρi is the ion gyroradius, and using the Klimontovich-Dupree representation of
Fα =

N
�
j=1

δ(R − Rαj )δ(μ − μαj )δ(v� − v�αj ),(1)where F is the phase space distribution in the gyrocenter2
coordinates of (R ≡ x − ρ, μB , v� ), μ ≡ v⊥
/2 is the magnetic moment, α denotes species, x is the usual spatial
coordinate, ρ is the gyroradius, the equations of motions
are
dR
∂ φ̄
= v� b̂ + vd −
× b̂,
dt
∂R
dv�
= −b̂∗ ·
dt

μB ≡

(2)

�
� 2
∂ φ̄
v⊥ ∂
lnB +
,
2 ∂R
∂R

(3)

2
v⊥
∂
(1 − v� b̂ ·× b̂) ≈ const.,
2B
∂R

(4)
v2

∂
∂
b̂∗ = b̂ + v� b̂ × (b̂ · ∂R)b̂, vd = v�2 b̂ × (b̂ · ∂R)b̂ + 2⊥ b̂ ×∂
∂R lnB, where b̂ is the unit vector in the direction of∗ Electronic address: ganesh@ipr.res.in

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

R. Ganesh et al., Properties of Freely Decaying and Driven Turbulence of Fusion Plasmas Using Gyrokinetic Particle Simulation

the external magnetic ﬁeld, � and ⊥ are the directions
parallel and perpendicular to b̂, respectively. The transformation between the gyrocenter coordinates R and theusual particle coordinates
x associated with a gyroradius�ρ yields φ̄(R) = � φ(x)δ(x − R − ρ)dx�ϕ , where �· · · �ϕis the average over the gyro-angle ϕ. The corresponding
gyrokinetic Poisson’s equation [11] becomes(m, n) �= (0, 0), Eq. (5) can be solved using an iterative
scheme [17], whereas, for (m, n) = (0, 0), the gyrokinetic
2 2
ρi � 1 as
Poisson’s equation can be simpliﬁed for k⊥∇2⊥ φ = −n̄i + ne .

Numerically, the transformation between R and x canbe accomplished through a 4-point average process valid
for k⊥ ρi ≤ 2 [14]. The δf method, which is based on thescale separation between the background and perturbed
quantities [15] and is utilized for the simulations reported
in the present paper, is based on Eqs. (2), (3), (4) and
�
�
dw
∂ φ̄
Te
∂ φ̄
= −(1 − w) κ× b̂ · r̂ + (v� b̂ + vd ) ·
,
dt
∂R
Ti
∂R
(7)
where

III.

(8)

is the background inhomogeneity with κn ≡ 1/Ln andκT i ≡ 1/LT i representing the spatial gradients in density and temperature, respectively. [Note that the ion
acoustic speed cs is unity in these units.] The perturbed
distribution is deﬁned as
δfα =

N
�
j=1

THE SIMULATION CODE AND
PARAMETERS

The gyrokinetic toroidal code (GTC) [2, 16] is aparticle-in-cell (PIC) code in global geometry based on
the gyrokinetic Vlasov-Poisson equations, Eqs. (1)-(11),where the gyrophase averages associated with �· · · �ϕ werecarried out by representing a gyrokinetic particle as acharged ring [14]. An iterative scheme [17] was used tosolve the integral part of the gyrokinetic Poisson’s equation, Eq. (5), for (m, n) �= (0, 0) modes, and a direct
one-dimensional ODE solver in the radial direction wasused for (m, n) = (0, 0) mode for the zonal ﬂows, bothin conﬁguration space. As noted earlier, the electrons
are assumed to be adiabatic in the present model, andonly the ions are ”pushed” in the code. The code uses
the magnetic coordinates (ψ, θ, ζ) to describe a tokamakwith a circular cross section, where the concentric ﬂux
surfaces are represented by ψ(r) in the radial direction,and θ and ζ are the poloidal and toroidal angles, respectively. On the poloidal plane, a 2D mesh with a uniform but unstructured grid in ψ and θ is used with Δ(mesh size) ≈ ρi (ion gyroradius). The code also uses
a global ﬁeld-aligned mesh [18, 19], taking advantage ofthe basic property of microturbulence in tokamaks thatk� � k⊥ . As such, the resolution in the toroidal direction
is closely related to the resolution in the poloidal direction resulting in a reduction in number for the toroidal
grid [18, 20, 21]. A slight shift of the ﬁeld lines is thenwj δ(R − Rαj )δ(μ − μαj )δ(v� − v�αj ),where N is the total number of particle ions in the simulation, �Fα = Fα0 + δfα , Fα0 is the background Maxwellian
with Fα0 dx = 1,w ≡ δfα /Fα(11)

We should remark here that both Eqs. (5) and (11) are
two-dimensional equations perpendicular to the magnetic
ﬁeld.
The parallel nonlinearity is the last term on the righthand side of Eq. (3). This term has not been studied closely in the microturbulence community. However,
without this term, the energy conservation cannot be satisﬁed in the simulation [12, 13]. The only other nonlinearity in the simulation is the last term on the right handside of Eq. (2), the familiar E × B nonlinearity. It is the
key nonlinear term responsible for the nonlinear trappingof the resonant particles [5] and the generation of zonal
ﬂows, which, in turn, gives rise to the nonlinear saturation of ITG turbulence and it also regulates the nonlinear
state transport by breaking up the streamers. As shown
in Ref. [7], the nonlinear generation of zonal ﬂows is also
greatly inﬂuenced by the velocity space nonlinearity. The
terms associated the background inhomogeneity, κn andκT i in Eq. (8) are usually time independent. Recently,
simulations have been carried out by making these terms
time dependent to account for the quasilinear changes ofthe proﬁles [9], which we will describe.

τ [φ(x) − φ̃(x)] = −n̄i (x) + ne (x),(5)
�where φ̃(x)≡<φ̄(R)Fi (R, μ, v� )δ(R − x +ρ)dRdμdv� >ϕ is the second average of the gyrophase angle ϕ on the LHS of the equation, whichtransforms the φ̄(R) back to the usual coordinates in x.
On the RHS, the densities are deﬁned as�n̄α (x) = � Fα (R)δ(R − x + ρ)dRdv� dμ�ϕ . (6)2κ = κn − (3/2 − v 2 /2vti
)κT i

2

(9)

and Fα ≡ δfα (wj = 1). For the adiabatic electron model
used here, only the evolution of the ion distribution function is followed in time together with a simpliﬁed electron
density response given by�φ(x),(m, n) �= (0, 0),
ne (x) = 1 +(10)0,(m, n) = (0, 0),
where (m, n) are the poloidal and toroidal mode numbers,
respectively. The approximation used for ne is adequate
for the studies described in the present paper and n̄i isgiven by Eq. (6), where n̄i0 = 1. For the modes with518

R. Ganesh et al., Properties of Freely Decaying and Driven Turbulence of Fusion Plasmas Using Gyrokinetic Particle Simulation

3

needed in order to match the mesh in the code [16, 18]. A
more complete form of Poisson’s equation in the presenceof adiabatic electrons can be written asτ [φ(x) − φ̃(x)] − [φ(x) − �φ(x)�] = −δn̄i (x),(12)

where �· · · � denotes ﬂux surface average and δni is theperturbed ion response. When long wavelength large
amplitude waves co-exist with short wavelength small
amplitude waves in a system such as ours, more accurate treatments are needed to solve this equation, which
has been implemented in our latest Gyrokinetic Tokamak
Simulation (GTS) code, which can use the experimental
measurements as input [22].
The Cyclone-based-parameters [23] have been used forthe simulation with 64 toroidal grid points for a/ρi = 125
on each poloidal plane, where a is the minor radius of thetokamak. Thus, the shortest wavelength modes that canbe resolved in the code are of the order of k⊥ ρi ≈ 2. The
other parameters are: n0 = 40 − 80 (number of particles per cell), R/LT = 6.9, R/a = 2.79, Ln /LT i = 3.13,
Ωi Δt = 7.6 and Te /Ti = 1. All these runs have been conducted by taking into account of all the physics describedin Eqs. (2) - (3), including the nonlinearly generated
zonal ﬂows and the velocity space nonlinearity. Typicallytwo kinds of proﬁles are used: one is an exponential pro6
ﬁle of the inhomogeneity is given by (1/L)e−[(r−rc)/rw ] ,where L represents either the temperature scale length
LT i or the density scale length Ln with rc /a = 0.5 andrw /a = 0.35, while the other is a cosh-proﬁle whose inhomogeneity is given by (1/L)cosh−[(r − rc )/rw ]2 Notethat while the exponential proﬁle is broad at the center and varies slowly, the cosh-proﬁle varies much more
rapidly so that proﬁle-eﬀects, if any, could come into play.
IV.

FIG. 1: Two diﬀerent density/temperature proﬁles (TopRow) are used for driving the ITG modes. Note that in the“Cosh-proﬁle”, proﬁle variation eﬀects, if any, would appear
whereas in the ”Exp-proﬁle”, the core is nearly ﬂat. Fivediﬀerent locations radial locations (1-5) where various quantities measured will be displayed. (Bottom Row) Fluxes at
5 locations shown as the function of time for upto 700 cs /atimes.

of the zeroth-order proﬁles, given above, have then been
carried out. The comparisons of these two sets of results
should give us some indication of the eﬀects of proﬁlerelaxations on the simulations of ITG turbulence. As
shown in. Fig. 2, where the time evolution of the proﬁles are plotted, the proﬁle for the temperature gradienthas been modiﬁed by more than ±15% at the end of thesimulation, whereas the density gradient modiﬁcation issmall and mainly due to noise, as expected. Surprisingly,
as shown in Fig.3, no noticeable changes of the resulting
thermal ﬂuxes have been detected between the two cases.
These results raise the question concerning whether ITG
modes are really ”stiﬀ” and whether we have reachedsome kind of steady state. Details of the simulations canbe found in Ref. [9].
Further it is found that the system attains a quasisteady state in the absence of source or sink. It is expected that this is due to the nonlinear wave-particle
interaction. This idea can be tested in the following
manner. In GTC evolution where proﬁle changes areincluded, we start with R/LT � 6.9 as shown in Fig.4
for Cosh-proﬁle. Initial proﬁle is in continuous magenta
color. The ﬁnal state is shown in red. Clearly the turbulence has evolved at the cost of free energy in the gradients resulting in waves and interacting particles. The
corresponding ﬂuxes during the evolution is shown on theright panel. Note that the ﬂux at late times, has arrived
at a quasi-steady state. The background gradients havestopped evolving; this arrest comes about due to waveparticle interaction. If we now restart the simulations,after removing the wave-particle physics (by resetting

PROFILE RELAXATION AND ITG
TURBULENCE

To understand the issues of proﬁle relaxation in gyrokinetic PIC simulations using the δf method, we have performed simulations using the exponential/cosh proﬁles as
initial conditions for the background inhomogeneity andwe have also kept track of the quasilinear modiﬁcationsof the proﬁles during the course of the simulation, i.e.,κT i = −1 d(T0i + δT0i )
T0i
d�x

and
κn = −

1 d(n0i + δn0i )
,
n0i
d�x

where δT0i and δn0i are the ﬂux surface averaged evolution. However, there is no back reaction of the proﬁlechanges to the simulated turbulence. In Fig.1 (top row),the drives for both proﬁles are plotted and the corresponding ﬂuxes at 5 diﬀerent locations. Next, simulations that taking into account the proﬁle modiﬁcations

519

R. Ganesh et al., Properties of Freely Decaying and Driven Turbulence of Fusion Plasmas Using Gyrokinetic Particle Simulation

FIG. 3: Comparison of ﬂuxes at ﬁve diﬀerent locations indicated in Fig.1: Proﬁle changes shown in Fig.2 does not seem
to aﬀect the resultant ﬂux, thus casting doubts on issues of“stiﬀness” in ITG turbulence.

FIG. 2: Freely evolving turbulence : background proﬁle evolution of κT proﬁles is shown for Cosh-proﬁle. Density also
evolves, however, the relative changes are weaker compared totemperature gradient evolution. Note that as time evolve, inthe linear stage, substantial changes in κT are observed withless changes towards late times showing quasi-steady state.
Similar evolution is seen for Exp-proﬁle.

the particle weights w to zero) and with the ﬁnal background proﬁles (see Fig.4 lower left panel), we see thatthe new linear stage, followed by a new steady state witha new background proﬁle is seen. We have repeated thisprocess to show that the proﬁles alter the linear characteristics and wave-particle interaction leads to a newquasi-steady state. A similar physics is seen for Expproﬁle shown in Fig.5. Note that in the restart runs, thezonal ﬂow potential has not been included. Hence several
restarts may systematically bring down the ﬁnal R/LT
values. However, if one includes the zonal potential, thenthis process of instability followed by background proﬁlemodiﬁcation, generation of wave-particle interaction andsaturation should stop at a ﬁnite value of R/LT .

V.

4

FIG. 4: March towards marginality: For Cosh-proﬁle, in orderto determine the role of wave-particle interaction, after quasisteady state is reached, the particle weights w are reset to zeroand the system is restarted from the new steady state. This
successive relaxation evolves every time to a new quasi-steadystate with relative low ﬂux values.

v|| NONLINEARITY AND FORMATION OF
GAM-LIKE m = 1 STATE

Standard Geodesic Acoustic Modes [10] were ﬁrst predicted in 1968. These are electrostatic potential perturbations with mode number (m = 0, n = 0) and density
perturbations with (m = 1, n = 0). Generation, sustenance and damping of standard GAMs is subject of intense research [24, 25] as they are associated with zonal
ﬂows which regulate transport due to low frequency drift
wave turbulence. Standard GAMs have real frequency
� cs /R where cs is the ion sound speed and R is the ma-

jor radius. They are either oscillatory or damped modes:
damping rate is typically small for large q values and reduces with collisionality. Thus it can be expected that ina collisionless plasma, standard GAM oscillations would
survive over ion-ion collision time and should be observable in late-time states of collisionless gyrokinetic simulations. It is also known that several GAM-like branches ofstandard GAM mode exists as has been reported recently
in [24, 25] as eigenmode problem. Standard GAMs have520

R. Ganesh et al., Properties of Freely Decaying and Driven Turbulence of Fusion Plasmas Using Gyrokinetic Particle Simulation

−5

2

−3

x 10

8

T=72.8 Phi

x 10

PNL=1 Exp−profile

0.02

T=109.3

1.5

6

0.015

1

4

0.01

0.5

2

0.005

0

0

0.5

1

1.5

0

−3

5

0.5

1

1.5

4

T=182.1

3

2.5

T=485.8

0.5

1
kθ ρ

1.5

2

0

0.5

1

1.5

x 10

T=728.6

1.5
1

1

1

0

2

2

2

T=145.7

−3

x 10

3

0

0

−3

x 10

4

0

0

5

0.5
0

0.5

1
kθ ρ

1.5

2

0

0

0.5

1
kθ ρ

1.5

2

FIG. 7: Evolution of Potential Fluctuation Spectrum withparallel velocity nonlinearity (PNL). As time progresses, thepower in the ﬂuctuations are seen to move towards larger
scales demonstrating at late times large scale ﬂuctuations are
dominant

FIG. 5: Same as Fig.4 but for Exp-proﬁlea/ρ=125, r0=0.5 a, detrended, centered0.06

a/ρ=125, r =0.5 a, detrended, centered
0.012

PNL=1

0.04

|φ(r0,ω)|2

φ(r0,t)

0
−0.02
−0.04

0.008
0.006
0.004
0.002

0

100

200

300 400t [a/cs]

500

0.06

600

0
−0.5

700

PNL=0

0.04

0
ω=2π/t|φ(r0,ω)|2

0
−0.02

0.5

PNL=0

0.015

0.02
φ(r0,t)

PNL=1

0.01

0.02

−0.06

0

0.01
0.005

−0.04
−0.06

0

100

200

300 400t [a/cs]

500

600

700

0
−0.5

ω=2π/t ,0
t normalized to a/cs

FIG. 8: Eﬀect of parallel velocity nonlinearity (PNL) on thefrequency spectrum of large scale GAM-like potential structure is shown. Top row shows electric potential time series
and its frequency spectrum measured at r = 0.5a with parallel nonlinearity on (PNL=1). On the bottom row is the casewith PNL=0. Frequency obtained for the case with PNL=1
appears to match with theoretical estimates.

FIG. 6: Mode structure of GAM-like m = 1 potential perturbation is shown when PNL=1. In our simulations, seed
for this global, low kr ρi mode forms just after the nonlinear
phase and sustains without damping throughout the late-time
evolution.

of this global GAM-like structure. Inclusion of this v|| nonlinearity not only produces a late time state as thatof an GAM-like structure but also results in a frequency
spectrum for our q−values close to the frequency of standard GAM modes. Note that the GTC simulations which
include v|| nonlinearity always produce global zonal ﬂow
structures with kr ρi � 1.
It is further demonstrated that the global GAM-likestructure is a result of an inverse cascade in poloidal mode
number kθ ρi . In Fig.7 we analyze the potential ﬂuctuations from linear stage in time to the fully nonlinear andradial structure akin to zonal ﬂows kr ρi � 1 where kr is
the radial mode number.
Here, undamped, global, GAM-like structures are reported with GTC simulations, which appear as the latetime states of gyrokinetic simulations. These GAM-like
potential perturbations have m = 1 and n = 0 and areglobal with kr ρi � 1, as shown in Fig.6. It is observedthat v|| -nonlinearity is fundamental to the development

521

0.5

R. Ganesh et al., Properties of Freely Decaying and Driven Turbulence of Fusion Plasmas Using Gyrokinetic Particle Simulation

6

ﬁnally at saturated stage. At the linear state, as canbe expected, the ﬂuctuation peaks at kθ ρi � 0.4 consistent with the linear stability models. As time evolves,due to quasi-two dimensional nature of the Tokamak
geometry k|| � k⊥ , an inverse cascade in power spectrum |φ|2 is seen clearly. For example, as time evolves,the spectral peaks at smaller and smaller values of kθ ρieventual reaching values where m � 1. Thus leading
to a global m = 1 GAM-like structure in the presenceof v|| nonlinearity. In Fig.8, we have obtained the frequency of ﬂuctuations at one point in (r, θ, φ) in thislate time state. Spectrum with and without v|| nonlinearity are clearly seen to be diﬀerent. Following Gao
et al [25], we consider for our√ values of�q, the expression ωGAM � √±0.2(1 + 0.4τ )( 2vti /qR) (1 − (1.2/q)2 )for ξ ≡ qRω/ 2vti � 1 and q > 1.2, ωGAM turns out
to be of the order ωGAM � 0.05cs /R for ξ = 0.14 and

R/a = 2.79, τ = 1. which is about the value obtained
including v|| nonlinearity.

[1] W. W. Lee, S. Ethier, R. Ganesh, R. Kolesnikov and W.
X. Wang, AIP Conf. Proc. 1069, 144 (2008).
[2] Z. Lin, T. S. Hahm, W. W. Lee, W. M. Tang and R.
White, Science, 281, 1835 (1998).
[3] S. E. Parker, W. W. Lee, and R. A. Santoro, Phys. Rev.
Lett. 71, 2042(1993).
[4] W. W. Lee and W. M. Tang, Phys. Fluids 31, 612 (1988).
[5] S. E. Parker et al., Phys. Plasmas 1, (5): 1461 (1994).
[6] W. X. Wang et al., Phys. Plasmas 14, 072306 (2007).
[7] W. W. Lee, S. Ethier, R. Kolesnikov, W. X. Wang, andW. M. Tang, ”Nonlinear Turbulent Transport in Magnetic Fusion Plasmas,” PPPL-4275, submitted to Computational Science and Discovery (2008).
[8] L. Villard et al., Plasma Phys. Contr. Fusion 46, B51(2004).
[9] R. Ganesh, W. W. Lee, R. Kolesnikov, S. Ethier, and J.
Manickam, ”Eﬀects on Proﬁle Relaxation on Ion Temperature Gradient Drift Instabilities,” invited talk, 20th International Conference on Numerical Simulations of Plasmas, Austin, TX (2007)[10] N. Winsor, J. L. Johnson, J. J. Dawson, Phys. Fluids 11,2248 (1968).
[11] W. W. Lee, Phys. Fluids 26, 556 (1983).
[12] D. E. Dubin, J. A. Krommes, C. Oberman and W. W.
Lee, Phys. Fluids 26, 3524 (1983).
[13] T. S. Hahm, Phys. Fluid 31, 2670 (1988).
[14] W. W. Lee, J. Comp. Phys. 72, 243 (1987).
[15] S. E. Parker, and W. W. Lee, Phys. Fluids B 5, 77 (1993).
[16] S. Ethier, W. M. Tang, and Z. Lin, J. Phys.: Conf. Series
16, 1 (2005).
[17] Z. Lin and W. W. Lee, Phys. Rev. E 52, 5646 (1995).
[18] Z. Lin, S. Ethier, T. S. Hahm, W M. Tang, Phys. Rev.
Lett. 88, 195004 (2002).
[19] A. M. Dimits, Phys. Rev. E 48, 4070 (1993).
[20] B.Scott, Phys. Plasmas. 5, 2334 (1998).
[21] W. X. Wang et al., Phys. Plasmas 13, 092505 (2006).
[22] W. X. Wang et al., Phys. Plasmas 14, 072306 (2007).
[23] A. M. Dimits, G. Batemen, M. A. Beer et al., Phys. Plasmas 7, 969 (2000).
[24] Z. Gao et al, Phys. Plasmas 13 100702 (2006)[25] Z. Gao et al, Phys. Plasmas 15 072511 (2008)

VI.

CONCLUSION

The eﬀects of nonlinear self-consistent proﬁle relaxation on the ITG turbulence in global delta-f particle
simulations using Gyrokinetic Toroidal Code (GTC) has
been reported. GTC shows that the turbulence is insensitive to minute or drastic proﬁle changes and reaches thesimilar ﬂux levels irrespective of whether or not it is adriven or freely decaying turbulence. Spatial and temporal spectral properties of gyrokinetic turbulence appear
to indicate that correct inclusion of parallel velocity space
nonlinearity may bring forth certain important late-timefeatures such as Geodesic Acoustic Modes (GAM).

522