Plasma and Fusion Research: Regular Articles

Volume 7, 2403094 (2012)

Kinetic Simulations of Neoclassical and Anomalous Transport
Processes in Helical Systems∗)
Hideo SUGAMA1,2) , Tomohiko WATANABE1,2), Masanori NUNAMI1,2), Shinsuke SATAKE1,2) ,Seikichi MATSUOKA1) and Kenji TANAKA1)
1)

2)

National Institute for Fusion Science, Toki 509-5292, JapanThe Graduate University for Advanced Studies (SOKENDAI), Toki 509-5292, Japan(Received 10 January 2011 / Accepted 18 May 2012)

Drift kinetic and gyrokinetic theories and simulations are powerful means for quantitative predictions ofneoclassical and anomalous transport fluxes in helical systems such as the Large Helical Device (LHD). Theδ f Monte Carlo particle simulation code, FORTEC-3D, is used to predict radial profiles of the neoclassicalparticle and heat transport fluxes and the radial electric field in helical systems. The radial electric field profiles
in the LHD plasmas are calculated from the ambipolarity condition for the neoclassical particle fluxes obtainedby the global simulations using the FORTEC-3D code, in which eﬀects of ion or electron finite orbit widths
are included. Gyrokinetic Vlasov simulations using the GKV code verify the theoretical prediction that theneoclassical optimization of helical magnetic configuration enhances the zonal flow generation which leads to the
reduction of the turbulent heat diﬀusivity χi due to the ion temperature gradient (ITG) turbulence. Comparisons
between results for the high ion temperature LHD experiment and the gyrokinetic simulations using the GKV-X
code show that the χi profile and the poloidal wave number spectrum of the density fluctuation obtained fromthe simulations are in reasonable agreements with the experimental results. It is predicted theoretically andconfirmed by the linear GKV simulations that the E × B rotation due to the background radial electric field
Er can enhance the zonal-flow response to a given source. Thus, in helical systems, the turbulent transport islinked to the neoclassical transport through Er which is determined from the ambipolar condition for neoclassicalparticle fluxes and influences the zonal flow generation leading to reduction of the turbulent transport. In order
to investigate the Er eﬀect on the regulation of the turbulent transport by the zonal flow generation, the flux-tube
bundle model is proposed as a new method for multiscale gyrokinetic simulations.
c 2012 The Japan Society of Plasma Science and Nuclear Fusion ResearchKeywords: neoclassical transport, anomalous transport, ITG turbulence, zonal flow, gyrokinetic simulation,
helical system, LHD
DOI: 10.1585/pfr.7.2403094

1. Introduction
Quantitative predictions of transport fluxes of particles, momentum, and heat in magnetically confined plasmas are a critical issue for the design of fusion reactors.
In high-temperature fusion plasmas, where particle mean
free paths are much larger than system sizes, kinetic eﬀects
such as particle orbits, finite gyroradii, and Landau damping need to be taken into account for accurate analyses ofplasma transport processes. These eﬀects can be properly
described by kinetic model which treats particle distribution functions on phase space. However, kinetic model isgenerally more complicated to use than fluid model because the phase-space dimension is higher than the realspace dimension. Therefore, a large number of kinetictransport studies are done by with the help of large-scalecomputer simulation [1, 2].
In magnetically confined plasmas, there exist colliauthor’s e-mail: sugama.hideo@LHD.nifs.ac.jp∗) This article is based on the invited presentation at the 21st International
Toki Conference (ITC21).

sional and anomalous (or turbulent) transport processes. In
the presence of plasma density and temperature gradients,
the particle distribution function F deviates from the local
Maxwellian fM and is written as
F = fM + f1 + fˆ1 ,(1)where f1 and fˆ1 represent the quasisteady and fluctuatingparts of the perturbed distribution function [3]. Using the
gyroradius ρ and the equilibrium gradient scale length L,
the magnitudes of f1 and fˆ1 are scaled as
f1
fˆ1
ρ
∼
∼ .
fM
fM
L
Charged particles gyrate around the magnetic field and the
gyrophase angle can be chosen as one of the phase spacecoordinates. Then, the quasisteady perturbed distributionfunction f1 , which is associated with collisional (classical and neoclassical) transport [4], can be divided into the
gyrophase-averaged part f 1 and the gyrophase-dependent

2403094-1

c 2012 The Japan Society of Plasma
Science and Nuclear Fusion Research

Plasma and Fusion Research: Regular Articles

Volume 7, 2403094 (2012)

2. Drift Kinetic Simulation of Neoclassical Transport in Helical Systems

part f˜1 :f1 = f 1 + f˜1 .
The gyrophase-dependent part f˜1 , which describes the particle gyro-motion, gives the diamagnetic plasma flows
and friction forces in the direction perpendicular to themagnetic field, from which classical transport fluxes ofparticles and heat are derived. The gyrophase-averaged
part f 1 represents the guiding-center distribution function,from which neoclassical transport fluxes are evaluated.
The guiding-center distribution function f 1 is governed by
the drift kinetic equation [5]. In addition to the abovementioned classical and neoclassical transport, anomaloustransport results from turbulent fluctuations driven by microinstabilities [6, 7]. Anomalous or turbulent transport
fluxes are derived from the fluctuating part fˆ1 of the perturbed distribution function, which is determined by solving the gyrokinetic equation [8–10].
Neoclassical and anomalous transport processes arestrongly influenced by toroidal magnetic configurations.
In helical systems such as stellarators and heliotrons [11],three-dimensional magnetic structures cause complex particle orbits which yield transport mechanisms diﬀerent
from those in tokamaks. In this paper, recent results fromkinetic simulation studies of neoclassical and anomaloustransport in helical systems are reported. As an interesting
feature of helical plasmas, the link between neoclassicaland turbulent transport processes arises through the background radial electric field. Because of the nonaxisymmetry, drift kinetic simulations can determine not only neoclassical particle and heat fluxes but also the backgroundradial electric field from the ambipolarity condition for theneoclassical particle fluxes. Then, the radial electric field
produces the E × B rotation of helical-ripple-trapped particles, which influences the zonal flow generation and accordingly the turbulent transport. In order to treat this effect of the E×B rotation on the turbulent transport, we propose the gyrokinetic simulation using the flux-tube bundlemodel in the present work.
The rest of this paper is organized as follows. Examples of drift kinetic simulation are presented in Sec. 2,
where radial electric field profiles in Large Helical Device(LHD) plasmas [12] are determined from the ambipolarity
condition for the neoclassical particle fluxes obtained byusing the global drift kinetic simulation code, FORTEC3D [13, 14]. Section 3 explains simulation studies using
gyrokinetic Vlasov simulation codes, GKV [15, 16] andGKV-X [17–19], to investigate ion temperature gradient(ITG) turbulence and zonal flows in LHD plasmas. In
Sec. 4, eﬀects of macroscopic radial electric fields on microscopic zonal flows are discussed and flux-tube bundlemodel is proposed for new multiscale gyrokinetic simulation to verify the enhancement of zonal flow generation
and resultant turbulent transport regulation predicted whenincreasing the radial electric field. Finally, conclusions are
given in Sec. 5.

Particles in helical systems are categorized into threeclasses which show diﬀerent types of orbits [11]. Twoclasses of them are passing particles and toroidally trapped
particles similar to those in tokamaks while a remaining
one consists of particles trapped locally within helical ripples. Drift motion of helical-ripple-trapped particles makes
dependence of neoclassical transport on the collision frequency ν and the radial electric field Er very diﬀerent fromthat in tokamaks. In the weakly collisional regime, theneoclassical diﬀusion coeﬃcient induced by helical-rippletrapped particles [20, 21] is qualitatively expressed by√ 2
νeﬀ
h vdr 2
νeﬀ + ω2E×B
 −1
νeﬀ(νeﬀ  |ωE×B |)∝νeﬀ /ω2E×B (νeﬀ  |ωE×B |) ,

D∼

(2)

√
where h , vdr , and νeﬀ ∼ ν/h represent the fraction, theaveraged radial drift velocity, and the eﬀective collision
frequency of helical-ripple-trapped particles, respectively,and ωE×B ∼ −cEr /rB denotes the frequency of the E × B
poloidal rotation.
It is well-known that, in tokamaks, the neoclassical transport fluxes of particles are intrinsically ambipolar; the ambipolarity condition is automatically satisfied
by the neoclassical particle fluxes for any radial electricfield Er [4, 5]. On the other hand, because of the abovementioned ripple transport, the ambipolarity of the neoclassical transport fluxes is nontrivial in helical systems
and the radial electric field Er can be determined fromthe neoclassical ambipolarity condition [11] which is written for a pure plasma consisting of electrons and a single
species of ions asΓencl (Er ) = Γincl (Er ).

(3)

Here, the functions Γencl and Γincl of Er represent the neoclassical particle fluxes of electrons and ions, respectively.
It was confirmed that the radial electric fields predictedfrom the neoclassical ambipolarity condition are in reasonable agreement with experimental results in helical devices [22, 23] even when the particle fluxes themselves are
dominated by anomalous transport. This fact is understandable by noting that the turbulent particle fluxes predicted from the gyrokinetic model are intrinsically ambipolar to the leading order in the gyroradius expansion [3] andtherefore they do not contribute to the ambipolarity condition for both helical systems and tokamaks.
A neoclassical transport simulation code, FORTEC3D, has been developed by Satake et al., which is applicable to general three-dimensional configurations [13, 14].
The FORTEC-3D uses the δ f Monte Carlo method andgives a global solution of the drift kinetic equation, whichis written as2403094-2

Plasma and Fusion Research: Regular Articles

Volume 7, 2403094 (2012)∂f1+ v b · ∇ f 1 + vd · ∇ f 1 + vd · ∇ fM = C( f 1 ), (4)∂twhere vd is the guiding-center drift velocity, C is
the linearized collision operator, f 1 represents the perturbed guiding-center distribution function on the fivedimensional phase space, and the phase space coordinates
are given by the guiding-center position (r, θ, ζ), the particle energy E ≡ mv2 /2 + eΦ, and the magnetic momentμ ≡ mv2⊥ /(2B). The third term vd · ∇ f 1 on the right-handside of Eq. (4) contains nonlocal eﬀects due to finite orbit
widths (FOWs). The radial electric field Er is also calculated by FORTEC-3D solving
⎛
⎞

⎜⎜
1
c2 ⎟⎟⎟ ∂Er2⎜⎜= e Γencl − Γincl ,|∇r| ⎝1 + 2 ⎟⎠
4π
∂t
vA

(5)

with the neoclassical particle fluxes obtained by substituting f 1 into

Γ

ncl

=

d3 v f 1 vd · ∇r .

(6)

Figure 1 shows profiles of the radial electric field in
the LHD configuration obtained from the global neoclassical transport simulations using FORTEC-3D and from theambipolarity condition using the local neoclassical transport model. In the simulations for the cases of Fig. 1, themagnetic field strength on the magnetic axis and the central densities and temperatures for electrons and ions are
given by B0 = 1.65 T, ne (0) = ni (0) = 2 × 1018 m−3 , andT e (0) = T i (0) = 1 keV. Under the conditions considered inFig. 1, the global simulations are done only for ions while
the the local neoclassical transport model is used for electrons, the FOW eﬀect of which is negligibly small. Details of density and temperature profiles used in the simulations are found in Ref. [14]. The major radius of thevacuum magnetic axis position is given by Rax = 3.7 m
and Rax = 3.6 m for the cases of top and bottom panelsof Fig. 1, respectively. Neoclassical transport is reduced in
the inward-shifted magnetic configuration with the vacuum
magnetic axis position Rax = 3.6 m than in the case with
Rax = 3.7 m because orbits of helical-ripple-trapped particles show larger radial displacements in the latter case. The
diﬀerence between the global simulation and local model
are not evident near the magnetic axis where helical ripples
are small. For the radial region r/a ≥ 0.4, the diﬀerence
becomes large because of the ion FOW eﬀect especially forthe case of Rax = 3.7 m, where the orbit widths of rippletrapped particles are wider than those for Rax = 3.6 m.
Thus, the FOW eﬀect tends to enhance the magnitude of
the radial electric field Er determined by the ambipolarity
condition for the neoclassical particle fluxes. However, it
is found that the diﬀerence between the resultant ambipolar neoclassical particle fluxes from the local model and the
FORTEC-3D is not so large as the change in Er (see Fig. 3
in Ref. [14]).

Fig. 1 Profiles of the radial electric field Er in the LHD configuration obtained from the global neoclassical transport
simulation using FORTEC-3D (red triangles) and fromthe ambipolarity condition using the local neoclassicaltransport model (open circles).

The FORTEC-3D code is also used by Matsuoka etal. [24,25] to investigate FOW eﬀects of electron drift motion on the radial electric field formation in Core ElectronRoot Confinement (CERC) plasmas [26]. The CERC plasmas obtained in several helical devices are characterized
by their high electron temperature (T e ) and steep T e gradient in the core region, and strong positive radial electricfield (Er > 0) called the electron root. Comparisons are
made between the Er profiles obtained from the FORTEC3D simulation, the local model, and the experiment for the
CERC plasma in LHD. It is confirmed that the diﬀerence
between the FORTEC-3D simulation and the local modeldue to the electron FOW eﬀect becomes large in the coreregion where the steep T e gradient appears. The Er profile obtained from the FORTEC-3D simulation shows a fairagreement with the experimental result except for the edge
region where stochastic magnetic fields neglected in thesimulation may cause additional electron transport.
As mentioned above, the FOW eﬀect generally enhances the negative Er for high-T i plasmas with T i ≥ T e
and it also reshapes the strong positive Er profile for CERC
plasmas with high T e although, in both cases, the resultant
ambipolar neoclassical particle fluxes are found to show
smaller changes due to the FOW eﬀect. We also note that,in general, the total transport is dominated by turbulent
transport. Thus, the FOW eﬀect is expected to influence
the plasma confinement not by changing the neoclassicaltransport fluxes but by altering Er which can regulate theturbulent transport as discussed in Sec. 4.

2403094-3

Plasma and Fusion Research: Regular Articles

Volume 7, 2403094 (2012)

3. Gyrokinetic Vlasov Simulationof ITG Turbulent Transport andZonal Flows in LHD Plasmas
The ITG turbulence is governed by the gyrokinetic
equation for the perturbed ion gyrocenter distribution function δ fi ,∂c+ (v b + vdi ) · ∇ + (b × ∇φ) · ∇ δ fi∂tBe∇φ
= (v∗i − vdi − v b) ·fiM + Ci (δ fi ),Ti(7)and the quasineutrality condition which is represented in
the perpendicular wave number vector (k⊥ ) space as eφ k⊥d3 v J0 (k⊥ v⊥ /Ωi )δ fik⊥ − 1 − Γ0 k⊥2 ρ2ti
Ti

e 
=
(8)
φ k⊥ − φ k⊥  ,Tewhere J0 is the zeroth-order Bessel function and Γ0 is defined by Γ0 (b) = I0 (b) exp(−b) with the zeroth-order modified Bessel function I0 . The gyrophase-averaged electrostatic potential φ is included in the nonlinear term on theleft-hand side of Eq. (7). The right- and left-hand sidesof Eq. (8) represent the perturbed ion and electron densities, respectively. Here, for the ITG mode, electrons are
assumed to make an adiabatic response to the electrostatic
potential. The fluctuation part fˆ1 of the particle distributionfunction in Eq. (1) is given in terms of the perturbed gyrocenter distribution function δ f and the potential fluctuation
φ asfˆ1k⊥ = δ f k⊥ exp(−ik⊥ · ρ)eφ k⊥ 
−1 − J0 (k⊥ ρ) exp(−ik⊥ · ρ) ,T

(9)where ρ ≡ b × v/Ω denotes the gyroradius vector. Simulation studies on ITG turbulence and zonal flows in helicalsystems are done by using the gyrokinetic Vlasov codes(GKV and GKV-X) which have been developed by Watanabe, Nunami et al. [15–19] The GKV and GKV-X codes
solve Eqs. (7) and (8) over a toroidal flux-tube domain [27]by using an Eulerian scheme with high phase-space resolution.
It is well-known that turbulent transport can be regulated by zonal flows [15, 28–36] which are generated byturbulent fluctuations themselves. Eﬃciency of the zonalflow generation due to the ITG turbulence is described by
the response function K(t). The function K(t) relates the
response of the zonal-flow potential φk (t) at time t to an
initially given source φk (0) with φk (t) = K(t)φk (0), where
k is the radial wave number. It is theoretically shown that
K(t) consists of the short- and long-time response partsdenoted by KGAM (t) and KL (t) [34]. The short-time response part KGAM (t) represents geodesic acoustic mode(GAM) oscillations [37] which vanishes due to the Landau damping in the long time limit t → +∞. The longtime response part KL (t), which shows slow temporal vari-ations, is more important for regulation of turbulent transport than KGAM (t). For helical systems such as the LHD,
the collisionless long-time response function KL (t) is written as [34]
1 − (2/π)1/2 (2H )1/2 {1 − gi1 (t, θ)}, (10)KL (t) =1 + G + E(t)/ n0 k⊥2 ρ2ti 
where definitions of H , gi1 (t, θ), G, and E(t) are foundin Ref. [34] and · · ·  represents the flux-surface average.
Equation (10) is derived for long radial wave numbers satisfying kρti  1, where ρti ≡ c(mi T i )1/2 /(eB) denotes thethermal ion gyroradius. In the denominator on the righthand side of Eq. (10), the geometrical factor G is associated with shielding of the zonal-flow potential by toroidally
trapped particles while E(t) represents the shielding effect due to helical-ripple trapped particles. The validity ofEq. (10) is verified by comparison to the linear gyrokinetic
simulation for the zonal-flow potential [34].
It is expected from Eq. (10) that zonal-flow generationis enhanced by neoclassical optimization which reduces radial drift of helical-ripple-trapped particles. This theoretical prediction is confirmed by using the GKV code for
gyrokinetic simulations of ITG turbulent transport in thehelical magnetic fields corresponding to the standard andneoclassically-optimized LHD configurations [15]. For
the neoclassically-optimized configuration with an inwardshifted LHD plasma, larger zonal-flow generation andsmaller turbulent ion heat transport are obtained by thesimulation than for the standard configuration. These theoretical and simulation results are consistent with the LHD
experiments [38] which show that not only neoclassical but
alsoanomalous transport is reduced in inward-shifted plasmas even though the inward shift increases unfavorable
magnetic curvature to destabilize pressure-gradient-driven
instabilitiessuch as the ITG mode. Direct comparisons
between gyrokinetic simulations and experiments can bemade by the GKV-X code [17–19] which uses detailed
geometric data of three-dimensional MHD equilibria obtained from the VMEC code [39]. The first validation testof the GKV-X code is done against the high-T i LHD experiment #88343 [40, 41]. Since a T i profile with a steep
gradient and a very flattened density profile are observed inthis experiment, trapped electron modes (TEMs) are stabilized and turbulent transport is anticipated to be dominantly driven by pure ITG modes which are verified to beconsiderably unstable in the radial region 0.3 < r/a < 0.9
by the linear GKV-X simulations [18].
Figure 2 shows turbulent potential structures aroundr/a = 0.65 obtained by the nonlinear GKV-X simulation
for the LHD experiment #88343. Three red open symbols in Fig. 3 show the turbulent ion heat diﬀusivity χi at
r/a = 0.46, 0.65, and 0.83 given by the GKV-X simulations of ITG turbulence under the conditions corresponding to the experiment. Error bars attached to these symbols
represent fluctuations of χi in the steady-state turbulence2403094-4

Plasma and Fusion Research: Regular Articles

Volume 7, 2403094 (2012)

Fig. 2 Turbulent potential structures around r/a = 0.65 obtainedby the nonlinear GKV-X simulation for in the LHD experiment #88343.

Fig. 3 Radial profiles of the ion heat diﬀusivities. The solid
curve represents the turbulent heat diﬀusivity obtained bysubtracting the neoclassical heat diﬀusivity from the totalheat diﬀusivity estimated by the LHD experiment #88343(dotted curve). The red symbols represent the turbulent
heat diﬀusivities at r/a = 0.46, 0.65, and 0.83 calculatedby the GKV-X simulations.

produced by the simulations. In Fig. 3, a solid curve represents the turbulent heat diﬀusivity obtained by subtracting
the neoclassical heat diﬀusivity from the total heat diﬀusivity estimated by the experiment (dotted curve). Error bars
for these curves are caused by inaccuracy of the observed
ion temperature gradient. We see that the GKV-X simulations reproduce the turbulent ion heat diﬀusivity χi givenby the experiment fairly well except that χi near the edge
(r/a = 0.83) is significantly underestimated by the simulation. Plasma beta values at r/a = 0.46, 0.65, and 0.83
are roughly given by β ∼ 0.4 %, 0.3 %, and 0.2 %, respectively. For these low beta values, electromagnetic microinstabilities such as kinetic ballooning modes (KBMs) are
predicted to be linearly stable while the ITG mode growth
rates are expected to be slightly reduced due to electromagnetic eﬀects according to Ref. [42]. Therefore, the present

GKV-X simulation results based on the electrostatic model
are not considered to be drastically changed by electromagnetic eﬀects although elaborate studies using the electromagnetic model remain as future tasks.
For the LHD experiment #88343, poloidal wave number spectra of the turbulent density fluctuation obtained
from the phase contrast imaging (PCI) measurement [43]and from the simulation are compared with each other [19].
The both spectra show similar shapes for high poloidal
wave numbers kθ > 0.5ρ−1
ti . The peak position of the spectrum obtained from the experiment is found at a higher
poloidal wave number than that obtained from the simulation. However, this quantitative diﬀerence appearing in
the low poloidal wave number spectra is not conclusive
when we take account of a large ambiguity of the experimental results due to the cutoﬀ at kθ ∼ 0.4ρ−1
ti and coarse
resolution in the PCI measurement for low poloidal wave
numbers. The above-mentioned results showing reasonable agreements between the LHD experiment and simulations encourage us to further gyrokinetic simulation studies
for anomalous transport in helical systems.
The GKV-X simulations for the LHD experiment#88343 also show that zonal flows are a critical factor todetermine the turbulent ion heat diﬀusivity χi . It is shown
from the nonlinear simulation results [19] that χi is approximately proportional to T/Z 1/2 where T and Z 1/2 represent the square of the turbulent potential amplitude andthe mean square root of the zonal-flow potential, respectively, in the saturated turbulence state. As found by M.
Nunami et al. [19], this relation of χi to (T, Z) reflects thefact that, even if the turbulence amplitude T has the samevalue, the larger zonal flow amplitude Z causes the spectral
transfer of turbulent fluctuations into higher radial wave
number regions, in which ITG modes are less unstable andturbulent transport occurs ineﬃciently. We also find that
T and Z have significant correlations with the linear ITG
growth rate and the linear zonal-flow response. Therefore,
it is expected that modeling of χi can be done by using results from linear analyses of the unstable modes and zonal
flows.

4. Radial Electric Field Eﬀects on
Zonal Flows
In the gyrokinetic turbulence simulations explainedin the previous section, a local flux-tube domain is usedbased on the idea of the ballooning representation [44], inwhich only the neighborhood of a single field line labeledby α ≡ θ − ζ/q(r) is considered. For helical systems, thefield line label α explicitly appears in the gyrokinetic equation in contrast to tokamak cases although α is regarded asa fixed parameter in the local flux-tube model. However,even if the zonal-flow 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, when tak-2403094-5

Plasma and Fusion Research: Regular Articles

Volume 7, 2403094 (2012)

ing account of the equilibrium or macroscopic radial electric field Er produced by the ambipolarity condition for
neoclassical radial particle fluxes (see Sec. 2), the E × B
rotation term, ωE×B ∂δ f /∂α, enters the zonal-flow component of the gyrokinetic equation and influences the zonalflow response in helical configurations. It is theoretically
predicted that E × B rotation induced by Er enhances thezonal-flow potential amplitude by reducing radial displacements of helical-ripple-trapped particles [45–48]. Whenthe Er eﬀect is included, the collisionless long-time response of the zonal-flow potential φk (t) to a initially given
source φk (0) is given for kρti  1 by [46]φk (t) =φk (0),−2
1 + Gp + Gt + Mp (Ght + Gh )(1 + T e /T i )(11)where geometrical factors Gp , Gt , Ght , and Gh represent
shielding eﬀects of neoclassical polarization due to diﬀerent types of particle orbits existing in helical systems andtheir definitions are given in Ref. [46]. The radial electric field Er enters the denominator in the right-hand sideof Eq. (11) through the poloidal Mach number defined by
Mp ≡ |(cEr B0 )/(rvti /Rq)| where q denotes the safety factor.
We see from Eq. (11) that zonal-flow generation can be enhanced when the geometrical factors Ght and Gh associatedwith orbits of helical-ripple-trapped particles are reducedby neoclassical optimization and/or when Mp is increased
by increasing Er or using ions with a heavier mass. Sincea higher zonal-flow response is predicted for a heavier ion
mass under the identical conditions on the magnitude of
Er and the magnetic geometry, the turbulent transport is
expected to show a more favorable ion-mass dependence
than the conventional gyro-Bohm scaling [46]. In order
to investigate the Er eﬀect, the simulation domain needs
to be extended in the α-direction from the flux tube to theshell region covering the whole magnetic flux surface. Linear gyrokinetic simulations of the zonal-flow response aredone by the extended GKV code [49], from which it is confirmed that the linear zonal-flow response is enhanced byincreasing Er as predicted by Eq. (11).
Because of helical ripples, gyrokinetic simulation
generally requires higher resolution of the phase space forhelical systems than for tokamaks. Especially, nonlinear
gyrokinetic simulation for the shell region over the whole
flux surface in helical systems needs much more computational memory and time than for a single flux-tube domain. Here, we propose a flux-tube bundle model as a new
method for multiscale gyrokinetic simulation which treats
the macroscopic E × B rotation and the microscopic zonal
flows with a smaller computational burden than direct simulation for the whole shell region. Based on the scale separation concept, we use the macroscopic (or large scale
length) coordinates (r, α, z) and the microscopic (or small
scale length) coordinates (x, y). Here, r, α ≡ θ − ζ/q(r),and z ≡ θ represent the radial coordinate, the field-line

Fig. 4 The flux-tube bundle model. A bundle of flux tubes used
for simulation domains are shown in red color.

label, and the poloidal angle, respectively, which are included in the equilibrium (or background) variables; forexample, the equilibrium magnetic field strength is written
as B = B(r, α, z). The (x, y) coordinates represent the sameas (r, α) although (x, y) are used to describe microscopic
scale variation of turbulent variables on the plane perpendicular to the magnetic field. For example, the fluctuating
potential is written as φ = φ(x, y, z; r, α) where (x, y) and(r, α) are treated as independent pairs of coordinates to separately represent microscopic and macroscopic variations.
Figure 4 shows a bundle of flux tubes distributed overa given flux surface labeled by the macroscopic radial coordinate r. Each flux tube is specified by assigning certain constant values to r and α. We note that z ≡ θ isregarded as the coordinate along the magnetic field line
because the direction parallel to the magnetic field is given
by changing z with all other coordinates fixed. When microscopic fluctuations are considered for each flux tube,
unstable modes such as the ITG instability need to have finite wave numbers ky  0 in the y-direction. For ky = 0,
fluctuations give linearly stable modes such as zonal flows.
The gyrokinetic equation for the zonal (ky = 0) component of the perturbed ion gyrocenter distribution function
δ fi (x, z, E, μ; r, α) is written as

v ∂
∂
∂
+
+ vdr
∂t qR ∂z∂x∂− Ci δ fi (x, z, E, μ; r, α)
+ωE×B
∂α


v ∂
e
∂
+ vdr
φ(x, z; r, α)= − fiM
TiqR ∂z∂x+N(φ, δ fi ),(12)

where N(φ, δ fi ) denotes the ky = 0 component of thenonlinear E × B convection term that is regarded as asource of zonal flows. The eﬀect of the equilibrium ra-

2403094-6

Plasma and Fusion Research: Regular Articles

Volume 7, 2403094 (2012)dial electric field Er on the zonal flow generation appears
from the E × B rotation term ωE×B ∂δ fi /∂α on the lefthand side of Eq. (12). This term causes the interaction
between the zonal (ky = 0) modes distributed over thebundle of flux tubes with diﬀerent values of α and thus
influences zonal flow generation and turbulent transport.
When diﬀerent flux-tube regions are coupled through thebackground E × B rotation in the flux-tube bundle model,
the perturbed distribution function is considered to show
poloidally-global structures reflecting poloidal excursionsof helical-ripple-trapped particles which cannot be grasped
by conventional single-flux-tube simulations. We expectfrom this coupling that the turbulent transport processes arenecessarily linked to the neoclassical transport because theambipolar neoclassical particle fluxes determine the E × B
rotation which can, in turn, enhance zonal-flow generation
and accordingly regulate the turbulent transport.
Without the E × B rotation, turbulent fluctuations ineach flux-tube domain are determined independently of
those in other flux tubes in the same way as in the case of asingle flux tube simulation. In axisymmetric systems suchas tokamaks where equilibrium variables do not depend onα, δ fi is set to be independent of α and accordingly the
E × B rotation term has no eﬀect even if ωE×B (∝ Er )  0.
Gyrokinetic turbulence simulation using the flux-tube bundle model described above is planned to confirm the Er effect on the zonal flows and the ITG turbulent transport, andits results will be reported elsewhere.
In the LHD experiment #88343, for which the GKVX simulations are done as described in Sec. 3, the data of
the background radial electric field are obtained only near
the peripheral region (r/a ≥ 0.8) (see Fig. 1 in Ref. [18]).
The maximum E × B drift velocity is given from these databy (vE×B )max  1.4 × 103 m/s at r/a = 0.83, which corresponds to the poloidal Mach number Mp  0.028. For
this small Mach number, the enhancement of the residual
zonal flow is anticipated to be weak according to Eq. (11)and Ref. [49]. Thus, as far as the radial electric field iswithin this level of Mp , we can justify the GKV-X fluxtube simulations in Sec. 3. However, if Mp becomes larger
by increasing the radial electric field and/or by using ionswith a heavier mass, larger zonal-flow generation and resultant turbulent transport regulation are expected.
It was shown by Anderson and Kishimoto [50] andby Uzawa et al. [51] that sheared mean flows can weaken
the generation rate of zonal flows by aﬀecting the modulational process. Such an eﬀect of the sheared mean flow on
the zonal flow generation is not included in Eqs. (11) and(12) where the background E × B flow shear is neglected
as a smaller term of higher order in ρ/L. The modulational
instability analysis of zonal flows including eﬀects of thesheared mean flow in helical systems still remains as oneof interesting and complicated subjects.

5. Conclusions
In this paper, recent results from kinetic studies ontransport processes in helical systems such as the Large
Helical Device (LHD) are reported. These results showthat drift kinetic and gyrokinetic theories and simulations
are useful tools for quantitative predictions of neoclassicaland turbulent transport fluxes with accurate treatment ofcollective particles’ motion in complex three-dimensional
magnetic field structures.
In helical systems, ripple-trapped particles give adominant contribution to neoclassical transport although
the ripple transport can be significantly reduced with increasing the radial electric field Er . In contrast to axisymmetric systems such as tokamaks, Er profiles in helicalsystems can be determined by the condition of ambipolar neoclassical radial particle fluxes because of rippletrapped particles’ drift motion causing diﬀerent dependence of neoclassical transport on the collisionality ν and
Er from that in the axisymmetric case. Using the δ f Monte
Carlo particle simulation code, FORTEC-3D, which solves
the drift kinetic equation in the global toroidal region, radial profiles of the neoclassical particle and heat transport
fluxes and the radial electric field in helical systems canbe calculated. The FORTEC-3D is used to investigate effects of ion finite orbit widths (FOWs) on the neoclassicaltransport and the Er profile in the LHD configurations forthe case of equal ion and electron temperatures, where theion FOW eﬀect are shown to be significant under the conditions corresponding to large radial excursions of rippletrapped particles. Neoclassical transport simulations using
the FORTEC-3D code are also done for Core ElectronRoot Confinement (CERC) plasmas with high electron
temperature T e and steep T e gradient in the core region.
Electron FOW eﬀects and a reasonable agreement between
the Er profile obtained from the simulation and the LHD
experiment are confirmed for the core region in the CERC
plasma.
ITG turbulent transport and zonal flows in helical systems are studied by gyrokinetic Vlasov simulation codes,
GKV and GKV-X, which employ a flux tube domain. The
GKV simulations show that the turbulent heat diﬀusivity
χi is reduced by enhanced zonal flow generation for theneoclassically-optimized magnetic configuration as theoretically predicted. The radial profile of χi and the poloidal
wave number spectrum of the density fluctuation obtained
from the high ion temperature LHD experiment #88343
are compared with the results from the GKV-X simulations
which can include precise geometrical data of MHD equilibria corresponding to experiments. The comparisons ofthe χi profiles and the fluctuation spectra show fairly good
agreements between the simulations and the experimental
results, which imply that we are on the right track towardfurther studies of anomalous transport in the LHD experiments and future helical reactors based on the gyrokinetic
simulations. The GKV-X simulation results also show that
χi is approximately proportional to the ratio between the2403094-7

Plasma and Fusion Research: Regular Articles

Volume 7, 2403094 (2012)square of the turbulent potential amplitude and the meansquare root of the zonal-flow potential in the saturated turbulence state.
It is theoretically predicted that, in helical systems, themacroscopic or background radial electric field Er can enhance generation of zonal flows leading to the further reduction of turbulent transport and the favorable isotope effect of the ion mass. The enhanced zonal-flow response
due to the increased Er is verified by the linear GKV simulations, for which the simulation domain is extended fromthe flux tube to the shell region covering the whole magnetic flux surface. In order to confirm the Er eﬀect on thezonal flows and the ITG turbulent transport, the flux-tube
bundle model is proposed as a new method for multiscale
gyrokinetic turbulence simulation.
As future tasks, we plan to extend the gyrokinetic simulation codes for including the flux-tube bundle model andother physical processes such as turbulent electron transport and electromagnetic fluctuations.

Acknowledgments
The authors thank Dr. M. Yokoyama and Dr. S.
Toda for fruitful discussions on neoclassical transportin helical systems, and also thank the LHD experiment group for providing experimental data anduseful discussions on them. This work is supported
in part by the Japanese Ministry of Education, Culture, Sports, Science and Technology (Grant Nos.
21560861 and 22760660) and in part by the NIFS
Collaborative Research Program (NIFS10KNTT003,
NIFS10KNST006, NIFS11KNST017, NIFS11KNXN229,
and NIFS11KNST014). Numerical simulations werecarried out with the use of the Plasma Simulator and the
LHD Numerical Analysis System at the National Institutefor Fusion Science.
[1] Y. Idomura, T.-H. Watanabe and H. Sugama, Comptes Rendus Physique 7, 650 (2006).
[2] X. Garbet, Y. Idomura, L. Villard and T.-H. Watanabe, Nucl.
Fusion 50, 043002 (2010).
[3] H. Sugama, M. Okamoto, W. Horton and M. Wakatani, Phys.
Plasmas 3, 2379 (1996).
[4] P. Helander and D.J. Sigmar, Collisional Transport in Magnetized Plasmas (Cambridge University Press, Cambridge,2002).
[5] S.P. Hirshman and D.J. Sigmar, Nucl. Fusion 21, 1079
(1981).
[6] W. Horton, Rev. Mod. Phys. 71, 735 (1999).
[7] W.M. Tang, Nucl. Fusion 18, 1089 (1978).
[8] A.J. Brizard and T.S. Hahm, Rev. Mod. Phys. 79, 421 (2007).
[9] E.A. Frieman and L. Chen, Phys. Fluids 25, 502 (1982).
[10] H. Sugama, Phys. Plasmas 7, 466 (2000).
[11] M. Wakatani, Stellarator and Heliotron Devices (Oxford
University Press, Oxford, 1998).
[12] A. Komori, H. Yamada, S. Imagawa et al., Fusion Sci. Technol. 58, 1 (2010).
[13] S. Satake, H. Sugama and T.-H. Watanabe, Nucl. Fusion 47,
1258 (2007).

[14] S. Satake, M. Okamoto, N. Nakajima, H. Sugama and M.
Yokoyama, Plasma Fusion Res. 1, 002 (2006).
[15] T.-H. Watanabe, H. Sugama and S. Ferrando-Margalet,
Phys. Rev. Lett. 100, 195002 (2008).
[16] T.-H. Watanabe, H. Sugama and S. Ferrando-Margalet,
Nucl. Fusion 47, 1383 (2007).
[17] M. Nunami, T.-H. Watanabe and H. Sugama, Plasma Fusion
Res. 5, 016 (2010).
[18] M. Nunami, T.-H. Watanabe, H. Sugama and K. Tanaka,
Plasma Fusion Res. 6, 1403001 (2011).
[19] M. Nunami, T.-H. Watanabe, H. Sugama and K. Tanaka,
Phys. Plasmas 19, 042504 (2012).
[20] K.C. Shaing, Phys. Fluids 27, 1567 (1984).
[21] S. Toda and K. Itoh, J. Plasma Fusion Res. 78, 582 (2002).
[22] M. Yokoyama, A. Wakasa, S. Murakami et al., Fusion Sci.
Technol. 58, 269 (2010).
[23] H. Maaßberg, C.D. Beidler, U. Gasparino et al., Phys. Plasmas 7, 295 (2000).
[24] S. Matsuoka, S. Satake, M. Yokoyama, A. Wakasa and S.
Murakami, Phys. Plasmas 18, 032511 (2011).
[25] S. Matsuoka, S. Satake, M. Yokoyama and A. Wakasa,
Plasma Fusion Res. 6, 123016 (2011).
[26] M. Yokoyama, H. Maaßberg, C.D. Beidler et al., Nucl. Fusion 47, 1213 (2007).
[27] M.A. Beer, S.C. Cowley and G.W. Hammet, Phys. Plasmas
2, 2687 (1995).
[28] A. Hasegawa and M. Wakatani, Phys. Rev. Lett. 59, 1581
(1987); H. Sugama, M. Wakatani and A. Hasegawa, Phys.
Fluids 31, 1601 (1988).
[29] Z. Lin, T.S. Hahm, W.W. Lee, W.M. Tang and R.B. White,
Science 281, 1835 (1998).
[30] P.H. Diamond, S.-I. Itoh, K. Itoh and T.S. Hahm, Plasma
Phys. Control. Fusion 47, R35 (2005).
[31] K. Itoh, S.-I. Itoh, P.H. Diamond et al., Phys. Plasmas 13,055502 (2006).
[32] A. Fujisawa, K. Itoh, H. Iguchi et al., Phys. Rev. Lett. 93,
165002 (2004).
[33] M.N. Rosenbluth and F.L. Hinton, Phys. Rev. Lett. 80, 724(1998).
[34] H. Sugama and T.-H. Watanabe, Phys. Plasmas 13, 012501(2006).
[35] P. Helander, A. Mishchenko, R. Kleiber and P. Xanthopoulos, Plasma Phys. Control. Fusion 53, 054006 (2011).
[36] P. Xanthopoulos, A. Mischchenko, P. Helander, H. Sugama
and T.-H. Watanabe, Phys. Rev. Lett. 107, 245002 (2011).
[37] N. Winsor, J.L. Johnson and J.J. Dawson, Phys. Fluids 11,2248 (1968).
[38] H. Yamada, A. Komori, N. Ohyabu et al., Plasma Phys.
Control. Fusion 43, A55 (2001).
[39] S.P. Hirshman and O. Betancourt, J. Comput. Phys. 96, 99(1991).
[40] K. Ida, M. Yoshinuma, M. Osakabe et al., Phys. Plasma 16,056111 (2009).
[41] K. Tanaka, C.A. Michael, L.N. Vyacheslavov et al., Plasma
Fusion Res. 5, S2053 (2010).
[42] H. Sugama and T.-H. Watanabe, Phys. Plasmas 11, 3068
(2004).
[43] K. Tanaka, C.A. Michael, L.N. Vyacheslavov et al., Rev.
Sci. Instrum. 79, 10E702 (2008).
[44] R.D. Hazeltine and J.D. Meiss, Plasma Confinement(Addison-Wesley, Redwood City, 1992) p. 298.
[45] H. Sugama, T.-H. Watanabe and S. Ferrando-Margalet,
Plasma Fusion Res. 3, 41 (2008).
[46] H. Sugama and T.-H. Watanabe, Phys. Plasmas 16, 056101

2403094-8

Plasma and Fusion Research: Regular Articles

Volume 7, 2403094 (2012)(2009).
[47] H. Sugama and T.-H. Watanabe, Control. Plasma Phys. 50,
571 (2010).
[48] H.E. Mynick and A.H. Boozer, Phys. Plasmas 14, 072507(2007).
[49] T.-H. Watanabe and H. Sugama, Nucl. Fusion 51, 123003(2011).
[50] J. Anderson and Y. Kishimoto, Phys. Plasmas 13, 102304(2006).
[51] K. Uzawa, Y. Kishimoto and J.Q. Li, J. Phys. Soc. Jpn. 77,
034501 (2008).

2403094-9