Proceedings of ITC/ISHW2007

Gyrokinetic Studies of Ion Temperature Gradient Turbulence andZonal Flows in Helical Systems
Hideo SUGAMA1,2) , Tomo-Hiko WATANABE1,2) , and Sergi FERRANDO-MARGALET1)
1)2)National Institute for Fusion Science, Toki 509-5292, Japan
Graduate School for Advanced Studies (Sokendai), Toki 509-5292, Japan

Gyrokinetic theory and simulation results are presented to investigate regulation of ion temperature gradient(ITG) turbulence due to E × B zonal ﬂows in helical systems. In order to examine eﬀects of changes in helicalmagnetic conﬁguration on anomalous transport and zonal ﬂows, magnetic ﬁeld parameters representing the standard and inward-shifted conﬁgurations of the Large Helical Device (LHD) [O. Motojima, N. Ohyabu, A. Komori,
et al., Nucl. Fusion 43, 1674 (2003)] are used. The linear gyrokinetic analyses show that the largest growth rate of
the linear ITG instability for the inward-shifted conﬁguration is higher than that in the standard one while, for theformer case, zonal ﬂows generated by given sources decay more slowly because of lower radial drift velocitiesof helical-ripple-trapped particles than for the latter as theoretically predicted. It is shown from the gyrokinetic
Vlasov simulation of the ITG turbulence that, in spite of the higher ITG-mode growth rate, the inward-shifted
plasma takes a smaller average value of the ion thermal diﬀusivity in the steady turbulent state with a higher
zonal-ﬂow level. These results imply that neoclassical optimization contributes to reduction of the anomaloustransport by enhancing the zonal-ﬂow level and give a physical explanation for the conﬁnement improvement
observed in the LHD experiments with the inward plasma shift. When equilibrium radial electric ﬁelds produce
poloidal E × B rotation of helically-trapped particles with reduced radial displacements, further enhancement of
zonal ﬂows and resultant transport reduction are theoretically expected.
Keywords: zonal ﬂow, ITG turbulence, gyrokinetic simulation, helical system, LHD

1 Introduction
In fusion science, numerous theoretical and experimental works have been done on zonal ﬂows which are now
well known to play a critical role in regulation of turbulent transport in plasmas [1, 2]. Therefore, in order to improve plasma conﬁnement in helical systems, where various geometrical conﬁgurations are explored [3, 4, 5, 6, 7],
it is very important to elucidate eﬀects of magnetic geometry on both microinstabilities and zonal ﬂows. This work
presents results from gyrokinetic theory and simulation toinvestigate regulation of ion temperature gradient (ITG)
turbulence due to E × B zonal ﬂows in helical systems.
It was shown in our previous papers [8, 9, 10, 11]
that, in helical systems, zonal ﬂows can be maintained
for a longer time by reducing the radial drift velocities ofparticles trapped in helical ripples. This implies a possibility that helical conﬁgurations optimized for reducing
the neoclassical transport can enhance zonal ﬂows and accordingly lower the turbulent transport as well because theneoclassical particle and heat ﬂuxes are also decreased byslowing down the radial drift of helical-ripple-trapped particles. In fact, it is observed in the Large Helical Device(LHD) [12] that not only neoclassical but also anomaloustransport is reduced by the inward plasma shift [13] which
decreases the radial particle drift but increases the unfavorable magnetic curvature to destabilize pressure-gradientauthor’s e-mail: sugama.hideo@LHD.nifs.ac.jp

driven instabilities such as ITG modes. This reduction of
anomalous transport by neoclassical optimization is a very
attractive property of helical systems to recent researcheson advanced concepts of helical devices [3, 4, 5, 6, 7]
It was shown by the ITG turbulence simulation in
our previous work [10, 11], in which model helical ﬁelds
for the standard and inward-shifted LHD conﬁgurations
were used in the gyrokinetic Vlasov (GKV) code [14], thatthe turbulent ion thermal transport in the inward-shiftedmodel, which has larger growth rates of the ITG stability,
was considerably regulated by the zonal ﬂows to a levelcomparable to the standard case although the thermal diffusivity χi for the inward-shifted case was slightly but stilllarger than for the standard case. However, in the recent
GKV simulation with more accurate conﬁguration models
installed [15], we ﬁnd that further stronger zonal-ﬂow generation occurs and makes χi smaller for the inward-shifted
conﬁguration [16].
Recently, Mynick and Boozer [17] predicted by using
the action-angle formalism that the collisionless residual
zonal-ﬂow level will be enhanced when the equilibrium
radial electric ﬁeld causes helical-ripple-trapped particles
to follow closed poloidal orbits with small radial displacements. Here, our gyrokinetic theory of zonal-ﬂow responseis also extended to analytically derive detailed expressions
for the eﬀects of the equilibrium electric ﬁeld on zonal
ﬂows in helical systems.

Proceedings of ITC/ISHW2007

2 Basic Equations
The nonlinear gyrokinetic equation for the perturbed ion
gyrocenter distribution function fik⊥ with the wave number
vector k⊥ perpendicular to the magnetic ﬁeld B is writtenas!
∂+ vk b · ∇ + iωDi fk⊥
∂t
!
eφk

= −vk b · ∇ − iωDi + iω∗T i F0 J0 (k⊥ ρ) ⊥
Ti
c X
0
00
0
[b · (k⊥ × k⊥ )]J0 (k⊥ ρ)φk0⊥ fik00⊥ ,+B k0 +k00 =k
⊥

⊥

⊥

(1)
where F0 is the local ion equilibrium distribution function that takes the Maxwellian form, J0 (k⊥ ρ) is the zerothorder Bessel function, ρ = v⊥ /Ωi is the ion gyroradius,and Ωi = eB/(mi c) is the ion gyrofrequency. The two
frequencies ωDi and ω∗i are deﬁned by ωDi = k⊥ · vDi
and by ω∗T i = ω∗i [1 + ηi {mi v2 /2T i − 3/2}] respectively,
where vDi ≡ (c/eB)b × (e∇Φ + µ∇B + mi v2k b · ∇b)
is the ion gyrocenter drift velocity, ω∗i ≡ k⊥ · (b ×∇r)(cT i /eB)(d ln n0 /dr) is the ion diamagnetic frequency,and ηi ≡ Ln /LT i is the ratio of the density gradient scale
length Ln ≡ −1/(d ln n0 /dr) to the ion temperature gradient scale length LT i ≡ −1/(d ln T i /dr). In Eq. (1), fik⊥ isregarded as a function of the kinetic energy w ≡ 12 mi v2 , themagnetic moment µ ≡ mi v2⊥ /(2B), and the toroidal coordinates (r, θ, ζ), where r, θ, and ζ denote the ﬂux surface
label, the poloidal angle, and the toroidal angle, respectively. A closed system of equations to determine the perturbed ion distribution function fik⊥ and the electrostatic
potential φk⊥ in the ITG turbulence are given by Eq. (1)and the quasineutrality condition,Z
d3 v J0 fik⊥ −n0

eφk⊥

e
[1 − Γ0 (b)] = n0φk⊥ − hφk⊥ i ,(2)TiTewhere h· · ·i represents the ﬂux-surface average and Γ0 (b) isdeﬁned by Γ0 (b) ≡ I0 (b)e−b with the zeroth-order modiﬁed
Bessel function I0 (b) and b ≡ k⊥2 T i /(mi Ω2i ).
The magnetic ﬁeld is written as B = ∇ψ(r) × ∇(θ −ζ/q(r)), where 2πψ(r) is equal to the toroidal ﬂux within
the ﬂux surface labeled r and q(r) represents the safety factor. In the present work, the radial coordinate r is deﬁnedby ψ = B0 r2 /2. Following Shaing and Hokin [18], we here
consider helical systems with the magnetic ﬁeld strength B
written by
B/B0=

1 − 00 (r) − 10 (r) cos θ − L0 (r) cos(Lθ)
X−
h(n) (r) cos{(L + n)θ − Mζ}=

1 − 00 (r) − T (r, θ)|n|≤nmax− H (r, θ) cos{Lθ − Mζ + χH (θ)},(3)

Table 1 Parameters at the ﬂux surface r ' 0.6a in the standard
and inward-shifted conﬁgurations.

standard
inwardstandard
inward

q
1.9
1.7
s
-0.85
-0.96

r/R0
0.099
0.114
0
r00
/t
0.22
0.71

t
0.087
0.082
rt0 /t
1.02
1.00

h /t0.91
1.20
rh0 /t
1.96
2.44

− /t-0.28-0.74
r−0 /t
-0.63
-0.36

+ /t0-0.24
r+0 /t
0
-0.61

where
T (r, θ) =H (r, θ) =χH (r, θ) =C(r, θ) =10 (r) cos θ + L0 (r) cos(Lθ),
p
C 2 (r, θ) + D2 (r, θ),arctan[D(r, θ)/C(r, θ)],X
h(n) (r) cos(nθ),|n|≤nmaxD(r, θ) =Xh(n) (r) sin(nθ),(4)|n|≤nmaxand M (L) is the toroidal (main poloidal) period number ofthe helical ﬁeld. In the present work, we use L = 2 andM = 10 to consider the LHD conﬁgurations. Here, we
assume that L/(qM)  1. Multiple-helicity eﬀects can be
included in the function H (r, θ). Hereafter, we put L0 = 0,0nmax = 1, and 00 = 0 (but 00≡ d00 /dr , 0) in Eq. (3) atthe radial position r that we consider.
In order to model the standard and inward-shifted
LHD conﬁgurations, we use numerical values shown in Table I for the safety factor q, the magnetic shear parameter ŝ,the inverse aspect ration r/R0 , the Fourier components of
the ﬁeld strength (t ≡ 10 , h ≡ h(0) , − ≡ h(−1) , + ≡ h(+1) ),and their radial derivatives. The two sets of parametersin Table I for the standard and inward-shifted conﬁgurations are called S-B an I-B, respectively, in [15]. These
parameters correspond to the ﬂux surface at r ' 0.6a (a :the plasma surface label) and they are taken from the vacuum magnetic ﬁeld data, which describe the LHD conﬁgurations more accurately than the model ﬁeld used in ourprevious study [10, 11]. The use of vacuum ﬁeld data isjustiﬁed because low beta plasmas are considered here.

3 Linear Analyses of ITG Modes
and Zonal Flows
In this section, the linearized version of Eq. (1) and thequasineutrality condition given by Eq. (2) are numericallysolved by using the GKV code in order to obtain the linear
dispersion relation for the ITG instability and the zonalﬂow response to the initial perturbation in the standard andinward-shifted conﬁgurations with the parameters shownin Table I.

Proceedings of ITC/ISHW2007

3.1

non-adiabatic helically trapped particles. We have E = 0
and gi1 = 1 at t = 0 because helically trapped particles
give no shielding before they begin radial drift. On theother hand, E approaches a ﬁnite value and gi1 ' 0 fort  τc = 1/(kr vdr ) where τc represents the characteristic
time for the shielding due to helically-trapped particles tooccur. The response kernel KL (t) for the long-time behavior of the zonal-ﬂow potential takes the constant limiting
values,Linear ITG instability

Figure 1 shows real frequencies and growth rates of the linear ITG instability as a function of the normalized poloidal
wave number kθ ρti where ρti ≡ vti /Ωi is the ion thermal gyroradius. Here, ηi ≡ Ln /LT i = 3, Ln /R0 = 0.3,
T e /T i = 1, α ≡ ζ − qθ = 0, and the parameters in Table
I are used. The real frequencies and growth rates for the
inward-shifted conﬁguration take similar values to thosefor the standard conﬁguration. Compared with the resultsin our previous work [10, 11], where simple model ﬁeld
parameters are used, the diﬀerence in the growth rates between the standard and inward-shifted conﬁgurations are
reduced because of changes in values of q, ŝ and magnetic
curvature although the maximum growth rate for the lattercase is still larger than for the former case.

3.2

K< ≡ lim KL (t) =t/τc →+0

K>

Collisionless time dependence of the zonal-ﬂow potential,
which has the wave number vector k⊥ = kr ∇r perpendicular to the ﬂux surface, is analytically derived as [9]eφk⊥ (t)TiZ teφk⊥ (0)1K(t)+
dt0 K(t − t0 )Ti
n0 hk⊥2 ρ2ti i 0
(
E)−1
2D
0
1/2 
1 − gi1 (t − t , θ)(2H )× 1−
π*Z0d3 v e−ikr vdri (t−t ) Fi0 S ik⊥ (t0 )
×
κ2 <1
Z
n

o+
3
0
+
d v Fi0 S ik⊥ (t ) 1 + ikr ∆r − h∆r ipo ,
κ2 >1

(5)
where K(t) is deﬁned byK(t) = KGAM (t)[1 − KL (0)] + KL (t).

(6)

Here, KGAM (t) and KL (t) are written as
KGAM (t) = cos(ωG t) exp(γt),

(7)

and

(9)

and

Zonal-ﬂow response

=

1
,
1+G

≡

lim KL (t)hi= hk⊥2 ρ2ti i 1 − (2/π)h(2H )1/2 in× hk⊥2 ρ2ti i[1 − (3/π)h(2H )1/2 i + G]
o−1
+ (2/π)(1 + T i /T e )h(2H )1/2 i .
(10)t/τc →+∞In Eq. (10), the term proportional to T i /T e is derived fromtaking account of the radial drift of helical-ripple-trappedelectrons which cannot be described by the perturbed electron density model used in Eq. (2). Therefore, this term
should be neglected when using Eq. (10) for comparison
to numerical solutions of Eqs. (1) and (2).
Responses of the zonal-ﬂow potential to the initial perturbation K(t) = hφk⊥ (t)i/hφk⊥ (0)i obtained by the linear
gyrokinetic simulation for the standard and inward-shifted
conﬁgurations are shown in Fig. 2, where the initial condition for the perturbed ion gyrocenter distribution function
is given by fik⊥ (t = 0) = nk⊥ (t = 0) exp(−mi v2 /2T i ) withnk⊥ (t = 0) determined from φk⊥ (t = 0) through Eq. (2).
Figure 2 also shows KL (t) predicted by Eq. (8) for comparison to simulation results. We see that the change inthe zonal-ﬂow response K(t) = hφk⊥ (t)i/hφk⊥ (0)i between
the standard and inward-shifted conﬁgurations are well described by KL (t) except for the GAM oscillations. In order
to measure the change in the zonal-ﬂow response, we deﬁne the zonal-ﬂow decay time as
Z tf
τZF ≡
hφk⊥ (t)i/hφk⊥ (0)i,(11)
0KL (t) ≡1 − (2/π)h(2H )1/2 {1 − gi1 (t, θ)}i
,
1 + G + E(t)/ n0 hk⊥2 ρ2ti i(8)respectively. Detailed deﬁnitions of variables in Eqs. (5)–(8) are found in [9]. In Eq. (7), the real frequency anddamping rate of the geodesic acoustic mode (GAM) [19]are denoted by ωG and |γ| = −γ(> 0). Equation (6) represents that the GAM oscillations described by KGAM (t)are superimposed around the averaged zonal-ﬂow evolution expressed by KL (t). We note that K(0) = 1 andlimt→+∞ KGAM (t) = 0. In Eq. (8), G represents the ratio of
the neoclassical polarization due to toroidally trapped ionsto the classical polarization while E(t) and {1 − gi1 (t, θ)} areassociated with the shielding caused by the radial drift ofwhere t f represents the time at which the zonal ﬂow
reaches the ﬁnal residual level and t f = 24R0 /vti is usedin the present cases. Then, the simulation results giveτZF = 1.68 R0 /vti for the inward-shifted case which isabout 65% larger than τZF = 1.02 R0 /vti for the standard
case. The increase in τZF is attributed to the decrease in radial drift velocities of helical-ripple-trapped particles andthe resultant delay in their shielding of the zonal-ﬂow potential in the inward-shifted conﬁguration. The improvement of the zonal-ﬂow response due to the inward plasma
shift was also found in our previous work [10, 11] using simpler conﬁguration models although the degree of
the improvement is more evident in the present study (see[15]).

Proceedings of ITC/ISHW2007

Fig. 1 Real frequencies and growth rates of the linear ITG instability as a function of the normalized poloidal wave
number kθ ρti for the standard and inward-shifted conﬁgurations. Here, ηi ≡ Ln /LT i = 3, Ln /R0 = 0.3, T e /T i = 1,
α = 0, and the parameters in Table I are used.

4 Nonlinear Simulation of ITG
Turbulence and Zonal Flows
This section presents nonlinear simulation results of the
ITG turbulence and zonal ﬂows obtained by solving
Eqs. (1) and (2) with the GKV code (see also [16]). The
GKV code employs the toroidal ﬂux tube domain and wehere use the same local plasma parameters (ηi ≡ Ln /LT i =3, Ln /R0 = 0.3, T e /T i = 1, and α = 0) as in the linear calculations in Sec. 3.1. Figure 3 shows the turbulent ion thermal diﬀusivity χi as a function of time t obtained by the GKV simulation with the magnetic ﬁeld data
in Table I used for the standard and inward-shifted conﬁgurations. We see that, as expected from the results in
Sec. 3.1, χi grows faster for the inward-shifted conﬁguration in the early time stage (t < 40Ln /vti ) than for the
standard conﬁguration and the peak value χi ' 3.8ρ2ti vti /Lnfor the former case is about 50% larger than the peak value
χi ' 2.6ρ2ti vti /Ln for the latter case. However, in later time(t > 60Ln /vti ), the turbulent transport reaches statisticallysteady states and then the average ion thermal diﬀusivityχi ' 1.45ρ2ti vti /Ln for the inward-shifted case is about 20%
smaller than the average value χi ' 1.78ρ2ti vti /Ln for the
standard case. This reversal of the χi -value order resultsfrom a greater amount of zonal ﬂows generated by turbulence in the inward-shifted plasma as seen below.
The GKV simulation shows that radially-elongatededdy structures (streamers) are ﬁrst driven by the toroidal
ITG instability although they are destroyed into small eddies by the self-generated E × B zonal ﬂows in the later

Fig. 2 Time evolution of the zonal-ﬂow potential for the standard and inward-shifted conﬁgurations. Here, kr ρti = 0.1
is used. Solid curves are obtained from the linear gyrokinetic simulation. Dashed curves correspond to the longtime zonal-ﬂow response kernel KL (t) in Eq. (8) which
does not include the GAM oscillations.

steady turbulent state [16]. The time-averaged spectrumof the zonal-ﬂow potential φkx is plotted in Fig. 4. It
is found that the zonal-ﬂow amplitude for kr ρti ' 0.25
in the inward-shifted conﬁguration is about three times
larger than that in the standard conﬁguration. The stronger
zonal-ﬂow generation in the inward-shifted case is consistent with the larger zonal-ﬂow decay time as mentioned inSec. 3.2. A typical radial scale length of the zonal ﬂows
observed in the helical ITG simulations is shown to be
shorter than those found in the tokamak ITG case for the
Cyclone DIII-D base case parameters. Accordingly, thezonal-ﬂow potential spectrum in the low kr -region has relatively smaller amplitude than for the tokamak case. Thistendency is also expected from the kr -dependence of thezonal-ﬂow response expressed in Eqs. (8) and (10).

5 Eﬀects of Equilibrium Radial
Electric Fields on Zonal Flows
So far, we have neglected eﬀects of the equilibrium electrostatic potential Φ(r) which yields the radial electric ﬁeld
E = Er ∇r (Er = −dΦ/dr) and accordingly the E × B drift
velocity vE ≡ (c/B)Er ∇r × b in the direction tangential tothe ﬂux surface. Regarding the ITG modes, vE will justgive the Doppler shift k⊥ · vE to the real frequencies without changing the growth rates. For the zonal-ﬂow components with k⊥ = kr ∇r, at ﬁrst, the equilibrium electricﬁeld does not seem to inﬂuence the zonal-ﬂow response

Proceedings of ITC/ISHW2007

Fig. 3 Time evolution of the ion thermal diﬀusivity χi obtainedby the ITG turbulence simulations for the standard andinward-shifted conﬁgurations.

because k⊥ · vE = 0. However, when treating helical conﬁgurations, we ﬁnd subtle points about the above argumentwith respect to the zonal-ﬂow response. In the previous
sections, we have used the ballooning representation andthe local ﬂux tube model, in which only the neighborhoodof a single ﬁeld line labeled by α ≡ ζ − q(r)θ is considered. For helical systems, the ﬁeld line label α explicitly
appears in the gyrokinetic equation in contrast to tokamak
cases although we have so far regarded α as a ﬁxed parameter based on the above-mentioned local model. But,
even if the zonal-ﬂow potential φ is independent of α, theexplicit appearance of α in the magnetic drift terms of the
gyrokinetic equation causes the perturbed gyrocenter distribution function δ f to depend on α. Therefore, in helicalconﬁgurations, we generally have vE · ∇δ f , 0 so that thezonal-ﬂow response can be aﬀected by the existence of theequilibrium electric ﬁeld.
Compared with passing and toroidally trapped particles, helical-ripple-trapped particles will have their orbitschanged more greatly by the equilibrium radial electricﬁeld Er . The radial displacements of helical-ripple-trapped
particles are signiﬁcantly reduced when the E × B drift
due to Er generates their rapid poloidal rotations as shownin Fig. 5. For such cases, neoclassical ripple transport isreduced and, in addition, higher zonal-ﬂow responses are
expected because the shielding of the zonal-ﬂow potential
by the helically-trapped particles is weakened. This scenario was ﬁrst presented by Mynick and Boozer [17], who
employed the action-angle formalism and pointed out theanalogy between the mechanisms of zonal-ﬂow shielding
and neoclassical transport.
Taking account of the dependence of the perturbed
distribution function on the ﬁeld line label α, our formulation of zonal-ﬂow response is extended to derive detailed
expressions for the Er eﬀects on the zonal-ﬂow response.

Fig. 4 Time-averaged radial-wave-number spectrum of thezonal-ﬂow potential obtained by the ITG turbulence simulations for the standard and inward-shifted conﬁgurations. The time average is taken over 60 ≤ vti t/Ln ≤ 250.

We now assume the bounce centers of helically-trapped
particles to draw poloidally-closed orbits with the poloidal
angular velocity ωθ ≡ −cEr /(rB0 ). Furthermore, considering the helical conﬁguration with the single-helicity component, for which H = h is independent of θ, we ﬁnd that,for t  1/ωθ , the shielding term E due to the helicallytrapped particles in Eq. (8) is replaced with EEr deﬁned by
EEr =

15
t vti(2h )1/2 (kr ρti )2
8π
rωθ

!2
1+

Te
Ti

!
(12)

Using Eq. (12), the collisionless long-time limit of thezonal-ﬂow response kernel, which represents the residual
zonal ﬂow level, is now given not by Eq. (10) but by
KEr =1
1 + G + EEr /(kr ρti )2


!2
!−1

15
T e 
1/2 t vti= 1 + G + (2h )1+
 (13)8πrωθ
Ti We see that, as Er increases, KEr increases and approaches
the same value 1/(1 + G) as in Eq. (9) because EEr is inversely proportional to the square of Er . It is noted thatEEr /kr2 given from Eq. (12) corresponds to the product of
the helically-trapped-particles’ fraction (∼ h1/2 ) and thesquare of the radial orbit width ∆E (∝ 1/ωθ ∝ 1/Er ) of
helically-trapped-particles’ poloidal rotation (see Fig. 5),
which agrees with Mynick and Boozer [17]. In helical conﬁgurations such as the inward-shifted LHD case, which are
optimized for reduction of neoclassical transport, the enhancement of zonal-ﬂow response due to Er is expected
to work more eﬀectively than in others because the neoclassical optimization reduces radial displacements ∆E ofhelically-trapped particles during their poloidal E × B rotation.

Proceedings of ITC/ISHW2007
Acknowledgments
This work is supported in part by the Japanese
Ministry of Education, Culture, Sports, Science andTechnology, Grant Nos. 16560727, 17360445, and17·05373, and in part by the NIFS Collaborative Research Program, NIFS06KDAD006, NIFS06KNXN060
and NIFS06KTAT038. Numerical simulations were carried out with the use of the Earth Simulator under the support from the Japan Agency for Marine-Earth Science andTechnology, and with the use of the Plasma Simulator andthe LHD numerical analysis system at the National Institute for Fusion Science.

Fig. 5 Orbit of bounce-center motion of helical-ripple-trapped
particles modiﬁed by the radial electric ﬁeld Er . Here,∆E represents the radial displacement of the orbit.

6 Conclusions
In the present work, eﬀects of changes in helical magnetic
conﬁguration on anomalous transport and zonal ﬂows areinvestigated based on gyrokinetic theory and simulation of
ITG turbulence and zonal ﬂows. In order to represent a
speciﬁc ﬂux surface (r ' 0.6a) in the standard and inwardshifted LHD conﬁgurations, magnetic parameters such as
the Fourier components of the ﬁeld strength, their radialderivatives, the safety factor, and the magnetic shear areused, which describe the conﬁgurations more accurately
than our previous model parameters. We ﬁnd from the linear analyses that the largest growth rate of the linear ITG
instability for the inward-shifted conﬁguration is slightly
higher than that in the standard one while, for the former
case, the zonal-ﬂow response is more favorable to generation of low-frequency zonal ﬂows because of lower radial drift velocities of helical-ripple-trapped particles thanfor the latter as theoretically predicted. The nonlinear
gyrokinetic simulation shows that the turbulent ion thermal diﬀusivity χi for the inward-shifted plasma takes ahigher peak value in the early time stage but a lower average value in the later steady turbulent state with stronger
zonal-ﬂow generation. Thus, it is conﬁrmed that neoclassical optimization contributes to reduction of the anomalous transport by enhancing the zonal-ﬂow level. This
presents a physical mechanism to explain the conﬁnement
improvement observed in the LHD experiments with theinward plasma shift. Also, further enhancement of zonal
ﬂows and resultant transport reduction are theoretically expected when the equilibrium radial electric ﬁeld Er causes
poloidal E × B rotation of helically-trapped particles with
reduced radial displacements. Simulation studies on the Er
eﬀects require global treatment in the direction parallel tothe E × B drift velocity and remain as a future task.

[1] P. H. Diamond, S.-I. Itoh, K. Itoh, and T. S. Hahm, Plasma
Phys. Controlled Fusion 47, R35 (2005).
[2] K. Itoh, S.-I. Itoh, P. H. Diamond, et al., Phys. Plasmas 13,055502 (2006).
[3] M. Yokoyama, J. Plasma Fusion Res. 78, 205 (2002).
[4] H. E. Mynick, Phys. Plasmas 13, 058102 (2006).
[5] D. A. Spong, S. P. Hirshman, L. A. Berry, et al., Nucl. Fusion
41, 711 (2001).
[6] J. N. Talmadge, V. Sakaguchi, F. S. B. Anderson, et al., Phys.
Plasmas 9, 5165 (2001).
[7] G. Grieger, W. Lotz, P. Merkel, et al., Phys. Fluids B 4, 2081
(1992).
[8] H. Sugama and T.-H. Watanabe, Phys. Rev. Lett. 94, 115001(2005).
[9] H. Sugama and T.-H. Watanabe, Phys. Plasmas 13, 012501(2006); H. Sugama and T.-H. Watanabe, ibid. 14, 079902(2007).
[10] T.-H. Watanabe and H. Sugama, 21th IAEA Fusion Energy
Conference (Chengdu, China, 2006), EX/5-4.
[11] T.-H. Watanabe, H. Sugama, and S. Ferrando-Margalet,
Nucl. Fusion 47, 1383 (2007).
[12] O. Motojima, N. Ohyabu, A. Komori, et.al., Nucl. Fusion
43, 1674 (2003).
[13] H. Yamada, et al., Plasma Phys. Controlled Fusion 43, A55(2001).
[14] T.-H. Watanabe and H. Sugama, Nucl. Fusion 46, 24(2006).
[15] S. Ferrando-Margalet, H. Sugama, and T.-H. Watanabe,“Zonal ﬂows and ion temperature gradient instabilities in
multiple-helicity magnetic ﬁelds”, to be published in Phys.
Plasmas.
[16] T.-H. Watanabe, H. Sugama, and S. Ferrando-Margalet,
Joint Conference of 17th International Toki Conference on
Physics of Flows and Turbulence in Plasmas and 16th International Stellarator/Heliotron Workshop 2007 (Toki, Japan,2007), P2-06.
[17] H. E. Mynick and A. H. Boozer, Phys. Plasmas 14, 072507(2007).
[18] K. C. Shaing and S. A. Hokin, Phys. Fluids 26, 2136 (1983).
[19] N. Winsor, J. L. Johnson, and J. J. Dawson, Phys. Fluids 11,2248 (1968).