Proceedings of ITC/ISHW2007

Zonal flows in 3D toroidal systems
H. E. MYNICK1) , and A. H. BOOZER2)
1)

2)

Plasma Physics Laboratory, Princeton University, Princeton, NJ, USA
Dept. of Applied Physics & Mathematics, Columbia University, New York City, NY, USA(Received 12 October, 2007)

A range of techniques for mitigating stellarator neoclassical (nc) transport has been developed, and attention
is now turning to also reducing the turbulent transport fluxes. As for tokamaks, zonal flows (ZFs) will be an
important tool in achieving this, and understanding the effect of machine geometry on these is important. This
paper discusses a theory of the shielding and time evolution of zonal flows in stellarators and tokamaks, which
attains greater generality and conciseness by use of the action-angle formalism. The theory supports the earlier
perspective that neoclassically-optimized devices should have less ZF damping, but it is pointed out that thefurther view, that this implies that optimized devices should therefore also have less turbulent transport, is overly
simplistic, neglecting the additional configuration dependence of the nonlinear source which drives the ZFs.
Keywords: stellarators, zonal flows, transport, turbulence, action-angle variables
DOI: 10.1585/pfr.1.001
Since the 1980s, a range of techniques for mitigatingstellarator neoclassical (nc) transport has been developed,
and attention is now turning to also reducing the turbulent
transport fluxes. As for tokamaks, an important tool forachieving this will be producing strong zonal flows (ZFs),
which can act to suppress the turbulence producing thetransport. It is thus of interest to understand how machinegeometry will affect the strength of these flows, and theireffectiveness in suppressing turbulent transport. This paper describes an analytic study of the shielding and time
evolution of ZFs in stellarators and tokamaks,[1] discusses
some of its implications, and its context in recent work on
ZFs.
ZFs are primarily poloidal E × B flows due to aradially-varying electrostatic potential φz (r, t) driven by thenonlinearity in the kinetic equation. For present purposes,
this is the Vlasov equation,(∂t + Ĥ0 )δ f (z, t) = −ĥ f0 − ĥδ f,(1)

with nonlinearity −ĥδ f ≡ −{δ f, h}. Here, b̂a ≡ {a, b} is
the Poisson bracket of a(z) with b(z) in the 6-dimensional
phase space z, and the Hamiltonian H(z, t) = H0 + h anddistribution function f (z, t) = f0 + δ f are divided into theirunperturbed (subscript 0) and perturbed portions, with f0
satisfying Ĥ0 f0 = 0. Here we consider electrostatic perturbations only, h(z, t) = eδφ(r, t), where r(z) is the particle
position. The dynamics of ZFs are determined by a selfconsistent loop between φz and the potential fluctuations
φk of the turbulence. Via the nonlinearity in Eq.(1), the φk
produce a source S ∼ |φk |2 for φz , which in turn affects
the growth and amplitudes of the φk .[2, 3, 4] The theory
developed in Ref. [1] follows earlier work[5, 6, 7, 8] inprimarily addressing the former of these 2 legs of the loop,
taking the nonlinearity −ĥδ f as a known source S f0 andauthor’s e-mail: mynick@pppl.gov

computing the resultant φz as a linear response problem,
k2 φz = 4πδρ xt /D, where δρ xt is the
R external charge-density
perturbation driven by S , δρ xt ∼ dtS (t), and D is the dielectric function. Regarding the latter leg, the effect of machine geometry on S is relatively unexplored to date, but is
an important issue in understanding the overall dynamics.
This is discussed further later in this paper.
In Ref. [1] the action-angle (aa) formalism[9] wasused to solve the kinetic equation without expansion ofthat equation in small parameters of radial excursions
and timescale, resulting in more general expressions forthe dielectric shielding, and extending results from earlier
work.[5, 6, 7, 8, 10] From these expressions, it was foundthat for each mechanism of collisional transport, there isa corresponding shielding mechanism, of closely related
form and scaling. Assuming the amplitude of the nonlinear source is unchanged, this correspondence supports
the suggestion raised in earlier work[11, 6, 12, 10] thatneoclassically-optimized stellarators will also have larger
ZFs, and consequently lower turbulent transport. On alonger, diffusive timescale, ZF evolution was shown to be
governed by a Langevin-like equation, with radial electricfield Er (t) moving diffusively about roots E a of the ambipolarity equation. The resultant probability distributionfunction is bounded, a balance between the turbulent fluctuations inducing diffusion and the neoclassical fluxes providing a restoring force toward E r = Ea . A fuller exposition of the analytic theory is given in Ref. [1]. Here werecap the elements of that theory, and discuss some of its
implications and related issues yet to be addressed.
Action-Angle formalism
In the aa formalism one parametrizes phase points z
with the 3 invariant actions J of the unperturbed motion
and their 3 conjugate angles θ, instead of the more directly
physical particle position r and momentum p. The unper-Proceedings of ITC/ISHW2007
turbed Hamiltonian is then independent of θ, H0 = H0 (J).
The key feature of aa variables is that they make the description of particle motion very simple. Hamilton’s equations are:= ∂J H = Ω(J) + ∂J h ' Ω(J),X
J̇ = −∂θ h = −ilhl (J, t) exp(il · θ),

(2)

θ̇

(3)

l

where ∂J (∂θ ) denotes a gradient in J (θ)-space, Ω(J) ≡∂J H0 , and l is the 3-component vector index, specifying
the harmonic of each component of θ in the Fourier deP
composition h(z) = l hl (J) exp(il · θ).
Following Refs. [5, 6, 7] in replacing the nonlinear
term −{δ f, h} with a specified source S (z, t) f0 , and Laplacetransforming in time and Fourier transforming in θ, one
readily obtains a solution for δ f , nonperturbative in the excursions made in a particle orbit,
δ fl (J, ω) = G0 [il · ∂J f0 hl (J, ω) +δ fl (J, t = 0) + S l (J, ω) f0 ],(4)with propagator G 0 ≡ (−iω + il · Ω + ν f )−1 , in which we
include an effective damping rate ν f , to later consider theeffect of collisions. Inserting Eq.(4) into the expression forthe
R charge density, (now showing species label s) δρ s (x) =dzρ(x|z)δ f s(z), where ρ(x|z) ≡ e s δ(x − r(z)) is the charge
density kernel and e s is the species charge, one obtains 3
contributions, labelled A, B, and C, corresponding to the 3terms on the right side of (4):Zδρ sA (x, ω) =dx0 K s (x, x0, ω)δφ(x0 , ω) (5)ZXδρ s,B+C (x, ω) = (2π)3 dJ
ρ∗l (x|J)G0 ×l[δ f sl (J, t = 0) + S sl (J, ω) f s0 ],(6)with response kernel[9]K s (x, x0, ω) = K sad (x, x0 ) + (2π)3ρ∗l (x|J)

Z

dJ

X
l

ω∂H0 f s0 + l · ∂J )H0 f s0
ρl (x0 |J).
l · Ω − ω − iν f

×(7)δρ sA , proportional to h or δφ, gives the self-consistent
response of the plasma, with response kernel K s . δρ sB, due
to the initial conditions of δ f , gives the transient ballistic
response, and the third term, δρ sC ,R is due to the nonlinear drive. K sad (x, x0) ≡ e s δ(x − x0 ) dzρ(x|z)∂H0 f s0 is the(generalized) adiabatic term, reducing to the usual adiabatic term when f0 is specialized to the local Maxwellianformn0
f M (J) ≡exp[−(H0 − eΦa )/T ),(8)(2πMT )3/2
where density n0 , ambipolar radial potential Φa , andtemperature T are functions of the drift-averaged minor
radius rd (J), and M is the particle mass.

Specialization to Toroidal Geometry
The expressions given thus far are valid for any system where a complete set J of constants of the motion exists. We now specialize to toroidal geometries, including
tokamaks and stellarators. We represent position in termsof flux coordinates r = (ψ, θ, ζ), where 2πψ is the toroidal
flux within a flux surface, and θ and ζ are the poloidal andtoroidal azimuths. In terms of these, the magnetic field
may be written B = ∇ψ × ∇θ + ∇ζ × ∇ψ p = ∇ψ × ∇α p , with2πψ p the poloidal flux, Clebsch angle α p ≡ θ − ιζ, constantalong a field line, and ι ≡ q−1 ≡ dψ p /dψ the rotationaltransform. α p and momentum (e/c)ψ form a canonicallyconjugate pair for motion perpendicular to the field line.
It is also useful to define an average minor radius r(ψ) by¯ = 0) the average magneticψ ≡ B¯0 r2 /2, with B¯0 ≡ B(rfield strength on axis. We consider toroidal systems with
the nonaxisymmetric portion of magnetic field strength B
dominated by a single helical phase η0 ≡ n0 ζ − m0 θ,¯
B(x) = B(r)[1− t (r) cos θ − δh (x) cos η0 ],(9)

but with ripple strength δh (x) allowed to vary slowly overa flux surface, with flux-surface average h (r) ≡ hδh i.
In such configurations, collisionless particle motion
occurs on 3 disparate timescales, the gyro, transit/bounce,and perpendicular drift scales, denoted by subscripts g, b,
and d. The characteristic gyro, bounce, and drift frequencies satisfy Ωg  Ωb  Ωd , and the corresponding radial excursions particles make on those scales are the gyroradius ρg = v⊥ /Ωg , the radial bounce excursion/banana
width ρb ' vBt /Ωb , and the “superbanana width” excursion¯
ρd ' σh vBt /Ωd . Here, vBt = t µ B/(MΩg r) is the grad-B
drift amplitude, µ is the magnetic moment, and σh equals
1 for a ripple-trapped particle (in trapping state τ = h), and0 otherwise. Thus, passing or toroidally-trapped particles(τ = p and t, resp.) have ρd = 0.
The aa variables can be chosen so that motion on eachtimescale can be described by one of the 3 pairs (θi , Ji ). A
¯ withsuitable choice is θ = (θg , θb , ¯αp ), J = (Jg , Jb , (e/c)ψ),Jg ≡ (Mc/e)µ the gyroaction, θg the gyrophase, describing
the fastest time scale of the motion, Jb the bounce action,θb its conjugate bounce phase, ψ¯ the drift-orbit averaged
value of ψ, and its conjugate phase ᾱp , the orbit-averaged
Clebsch coordinate α p , describing the slow, drift timescale.
To make the periodicity of the drift angle 2π as for the other¯ we use the closely related
2 phases, instead of ( ᾱp , (e/c)ψ)¯
canonical pair (θd , Jd = (e/c)ψd ), with θd ≡ ᾱp /(1 − ιqmn0),ψ¯d ≡ ψ¯− ψ¯p qmn0 , where qmn0 ≡ m0 /n0 . For typical pa¯ Correrameters, ιqmn0  1, so that (θd ' ᾱp , ψ¯d ' ψ).
spondingly one has the characteristic frequencies of motion Ω ≡ (Ωg , Ωb , Ωd ), and vector index l ≡ (lg , lb , ld ).
We adopt an eikonal form for the structure of any
mode a,φa (x) = φ¯a (r) exp iηa (x),(10)
Rrwith wave phase ηa (x) ≡ [ dr0 kr (r0 ) + mθ + nζ], andslowly-varying envelope φ¯a (r). Then using Eqs. (5,6,7),

Proceedings of ITC/ISHW2007
and (10) in the Poisson equation, one obtains the radiallylocal response equation[1]:ei φ¯a (r) X −2 X ∗k2 D(k, ω)=
λ sihGla (J)×Tisli[δ f sl (t = 0)/ f s0 + S sl (ω)]
i(ω − l · Ω + iν f s )(11)Here, Hλ2si ≡R T i /(4πn s0e s ei ), k2 ≡ |k|2 , h. . .i ≡(2π)−2 dθdζ dp( f0 /n0 ) . . . is the flux surface andmomentum-space average over the unperturbed
distriH−3
bution function f0 , Gla (J) ≡ (2π)
dθ exp(−il ·θ) exp iδηa (z) the “orbit-averaging factor”, a concise expression for the interaction of mode a with particles withactions J, and δηa the portion of ηa oscillatory in θ (sohaving zero θ–average). Dielectric function D is given by
P
D(k, ω) ≡ 1 + s χ s (k, ω), with susceptibility χ s (k, ω) =(kλ s )−2 g s (k, ω)), and shielding functionfXω − ω∗si. (12)g s (k, ω) = 1 −h|Gla (J)|2ω − l · Ω + iν f slHere, ω∗f ≡ ω∗ [1 + η(u2 − 3)/2], with ω∗ ≡ −kα cT/(eBLn)
the diamagnetic drift frequency, η ≡ d ln T/d ln n, u ≡ v/v s
the particle velocity, normalized to the thermal speed v s ,kα ≡ ld /r, and L−1n ≡ −∂ ln n0 /∂r. The 1 in g s comes fromthe adiabatic term K sad in Eq.(7). The 2 terms on the rightside of Eq.(11) arise from δρ s,B+C . This response equation
is of the same form as that obtained in Refs. [5, 6], or of
any linear response calculation. The differences lie in theform of the dielectric D, and in the use of the aa form,
which facilitates dealing with the range of timescales andof orbit-averaging effects in complex geometries in a general manner.
To evaluate the G l , we describe the radial motion byr ' rd + δr(d) (θd ) + δr(b) (θb ) + δr(g) (θg ), where for eachi = g, b, d, we make a harmonic approximation of the motion in that phase, δr (i) (θi ) ' ρi cos θi . This is a very good
approximation for gyromotion, and a good approximation
for bounce motion not too near a trapping-state boundary.
For simplicity, we assume that superbanana (τ = h) particles do not detrap, but precess poloidally dominated by
E × B poloidal drift, Ωd ' ΩdE , which is roughly constant on a given orbit, while drifting radially as v Bt sin θ,
as usual. This yields radial motion of the given harmonic
form, with superbanana width ρd = σh vBt /ΩdE , as noted
above. For completeness, one may also include in thisdescription the finite banana widths ρbh of τ = h particles, which give rise to the helically–symmetric nc transport branch in straight stellarators[13], and a second typeof superbanana width, the finite radial excursions ρdt madeby τ = t particles on the drift timescale, which give rise to
the “banana-drift” transport branch.[14, 15, 16] One thenfinds
Gla (J) = Jlg (zg )Jlb (zb )Jld (zd )e−iξa ,(13)with zg,b,d = kr ρg,b,d , and ξa a phase factor, whose value isirrelevant, since G l enters only as |G l |2 .

For drift turbulence, which is driving the ZFs, one typically has k⊥d ρgi ∼ 0.3, and frequencies ωk ∼ ω∗ (k⊥d ). For
ZFs, one has much smaller kr and frequencies ωZ , downby an order ofpmagnitude, perhaps by the “mesoscale” ratio, krZ /k⊥d ∼ ρgi /a. Thus, for both species, one has theordering ωZ , Ωd  Ωb  Ωg , and zg < zb < 1. For the moment we leave the relative sizes of ωZ and Ωd unspecified.
Also, one may have zd & 1 for trapped particles, for ionsand also, notably, for electrons, as noted by [6]. Thus, as
opposed to tokamaks, in stellarators electrons can participate in orbit averaging, because their radial excursions on
the drift timescale can be comparable with those of ions.
Because ωZ  Ωb,g , the sum over l in Eq.(12) isdominated by the terms with lg,b = 0, an approximation
strengthened for zg,b  1, for which the factors Jl2g,b in |Gl |2
in Eq.(12) are negligible unless lg,b = 0. These reduce thePtriple sum there to a single sum ld . In that sum, if one
has ω  Ωd , then over the ld –range ∆ld ∼ zd over which
Jl2d in Eq.(12) is appreciable the integrand does not change
greatly, so that one can perform the summation, using theP
identity l Jl2 = 1, which eliminates the Jl2d factor, leaving
only the factor Jl2g Jl2b . In the other limit ω  Ωd , the sum
is dominated by the ld = 0 term, and the effect of Jl2d survives. Thus, for ω  Ωd , all of gyro-, bounce- and driftaveraging contribute. Neglecting ν f s , Eq.(12) becomesg s (k, ω) ' 1 − Λ0b (bg , bb ), (Ωd  ω),(14)

g s (k, ω) ' 1 − Λ0d (bg , bb , bd ), (ω  Ωd ),where Λ0d (bg , bb , bd ) ≡ hJg2 Jb2 Jd2 i, Λ0b (bg , bb ) ≡2≡ J02 (zg,b,d ), bg ≡ kr2 ρ2gT ,Λ0d (bg , bb , bd = 0) ≡ hJg2 Jb2 i, Jg,b,dbb = bg q2 /(Ft t1/2 ), and bd ≡ kr2 ρ2dT , with ρgT ≡ vT /Ωg , vT
the species thermal velocity, ρdT ≡ ρd (v = vT ) ∝ v2T , andFt the fraction of toroidally-trapped particles.
The physics represented by Eqs.(14) is that if the theZF drive in a stellarator has a time variation slow compared with Ωd [cf. Eq.(14b)], τ = h particles have time topartially shield out φz by drifting along their collisionlesssuperbanana orbits, an averaging mechanism not available
to tokamaks. If the ZF drive varies rapidly compared with
Ωd [Eq.(14a)], this new mechanism for radial averaging islost.
The functions Λ0b and Λ0d succinctly describe theadditional contributions from finite bb , corresponding toshielding due to the “bounce” contribution gb to the shielding function computed in Refs. [5] and [6], and from finite bd , corresponding to a “drift” contribution gd to g, extending the result in Ref. [6]. Table 1 synopsizes someof the limiting cases covered by earlier work, extended inRef. [1], of which Eqs.(14) here are Eqs.(15) in [1] notedin the table. One notes that most of the entries are for collisionless theory, ν f = 0. In the tokamak limit (h , bd → 0),
Λ0d (bg , bb , bd ) → Λ0b (bg , bb ), so that Eq.(14a) again holds.
In the further cylindrical limit (t → 0), bb vanishes, andthe Λ’s in Eqs.(14) reduce to the more familiar Λ0 (bg ) ≡Λ0b (bg , bb = 0) ≡ hJg2 i = I0 (bg )e−bg , with I0 (b) the modi-Proceedings of ITC/ISHW2007
Parameter range
ν f = 0:tokamak limit (h , bd = 0), ω < Ωb :bb < 1bb arbitrarystellarators (h , bd > 0) :bd → ∞, ω < Ωdbd < 1, ω < Ωdbd < 1, Ωd < ω
bd arbitrary, ω < Ωdbd arbitrary, Ωd < ων f > 0:νh /Ωd > 1
νh /Ωd arbitrary

0

Ref.

[5]
[17]
[6, 7][1](16b)[1](16a)[1](15b)[1](15a)[8][1]

Table 1 Cases covered by the present theory

fied Bessel function of the first kind. For bg < 1, one hasΛ0 (bg ) ' 1 − bg , and thus g ' bg , the gyro- contribution
gg to g, corresponding to the classical (gyro) polarization
current J p,g .
Using the small-argument expansion J0 (z) ' 1−(z/2)2
for zd as well as zg,b in Λ0d , Eqs.(14) reduce to [1]g s (k, ω) ' bg + Ft cb bb , (Ωd  ω),(15)

g s (k, ω) ' bg + Ft cb bb + Fh cd bd , (ω  Ωd ),√where Fh = (2/π) 2h is the fraction of particles withτ =√ h, cd ' (15/2), and for a tokamak, one finds cb '10 2/(3π) ' 1.5, in approximate agreement with thevalue in Ref. [5].
One notes that the drift contribution gd = Fh cd bd 'Fh (kr ρd )2 in Eq.(15b) has a form analogous to the bounce
and gyro contributions, differing from the scaling gd ' Fhfound in Refs. [6, 7]. This is because in Refs. [6, 7],
the term Ωd ∂θd δ f was neglected in their kinetic equation,
thereby implicitly taking the limit bd → ∞ (see Table 1),
also of physical interest. Taking that limit in Eq.(14b) also
recovers that scaling.
As discussed in Ref. [1] and illustrated by Eqs.(15),there is a correspondence between the contributions g j tothe shielding function and the radial collisional (classical+nc) transport coefficients D j : the gyromotion producing the classical polarization term gg also gives rise to classical transport, the bounce motion producing gb gives riseto axisymmetric nc transport, and the drift motion yielding
gd also produces the “superbanana” branch of transport,
dominant in conventional stellarators. For each mechanism j, one may use the heuristic form D j ' F j ν f j (∆r j )2 ,with F j the fraction of particles participating in that mechanism, ∆r j the radial step in the random walk process, andν f j the effective stepping frequency in that random walk.
For example, for the 1/ν superbanana regime, one has
F j → Fh ' h1/2 , ∆r j → vBt /νh , and ν f j → νh ' ν/h .
And for the shielding contributions g j , one has the approximate form g j ' F j (kr ∆r j )2 , exemplified by Eqs.(15).

0Hence, g j /g j ' (D j /D j )(ν f j /ν f j0 ). Thus, for j → g andj0 → b, one expects the gyro- contribution gg in Eqs.(15) tobe smaller than the bounce contribution gb , because classical diffusion Dg is subdominant to banana diffusion Db .
Similarly, for j → b, j0 → d, one expects the bounce
contribution gb to dominate gd in Eq.(15b) as long as superbanana transport Dd is subdominant to Db . For NCSX,
for example, evaluations have shown[18] that, at a selfconsistent radial ambipolar field E a , NCSX should have
Dd down from Db by about an order of magnitude. Correspondingly, one expects that the ripple contribution gd to
ZF damping should be small compared with the tokamak
contribution gb .
It has been argued[11, 6, 12, 10] that neoclassicallyoptimized stellarators should also have lower turbulent
transport, due to less damping of ZFs. The basic idea ofmost nc optimization techniques has been to reduce ripple transport (Dd ) by reducing either F h , or by reducing
superbanana width ρd ' vBt /Ωd .[19] One notes from theresemblance of Dd to gd , characterized by the argument
1 2
2 2
2 hzd i ' F h kr ρdT , that the present theory supports this argument.
However, as noted in the introduction, the shieldingof ZFs from a given source, which most analytic ZF
studies to date (including the present one) address, isonly 1 of the 2 legs of the feedback loop in ZF dynamics.
That work demonstrates that neoclassically-optimized
stellarators will also tend to have lower damping of
ZFs. However, in general, different configurations will
have differing levels of drive for instabilities, and thus
differing strengths of the ZF source S . Thus, recent
gyrokinetic simulations[10] comparing configurations
modeling LHD in its (A)standard and (B)inward-shifted
operation, reported that the ZFs increased (by about50%) in going from A to B, as might be expected fromthe better nc-optimization of case B. However, because
case B also has much stronger (about 60%) ITG growth
rates, the turbulent flux in that case is increased (byabout 16%) over that in case A, in contrast with theexperimentally–observed reduction[11] in turbulent transport. Very recent further simulations[20] of this same
comparison, but with more realistic equilibrium profiles,
have brought the numerical transport trends more into
accord with experimental observations. Thus, it appears
that the now commonly-cited correlation between nc andturbulent optimization is too simple, and a more refined
understanding of this relationship must come from also
accounting for the effects of machine geometry on thesource strength.
Longer-time ZF evolution
As discussed in Ref. [1], the effects on ZFs from theg j (describing dielectric shielding) and the D j (describing
radial collisional fluxes) come together in the time evolution equation for the flux-surface averaged radial electricfield Er ≡ h∇r · Ei, obtained from the surface average ofProceedings of ITC/ISHW2007
Ampere’s law, plus an expression for the surface-averaged
radial current Jr ,∂t E rJr= −4πJr ,

(16)

−1

= (4π) χ∂t Er + σ(Er − Ea ) + FS /B.

The first term in Jr , proportional to the time derivative of
Er , represents the polarization current J p , with χ containing the dielectric shielding contributions. The second term
represents the nonambipolar radial current due to nc transport, from a first-order expansion in E ≡ (E r − Ea ) =P−h∇r · ∇φz i of the nc radial current e s Γ s (Er ), where
Ea = −h∇r · ∇Φa i is the ambipolar value at which theion and electron particle fluxes are equal. F S is the force,
coming from the source S in Eq.(1), exerted by the turbulence within a magnetic surface normal to the magnetic
field, which acts as a source driving E r . Using Eq.(16b)in (16a) yields a Langevin-like equation, which in the ω
domain may be written−iωE(ω) + γE E(ω) = cS (ω),(17)

where γE (ω) ≡ 4πσ/D(ω), cS (ω) ≡ −4πFS /BD(ω), andD(ω) ≡ 1 + χ(ω) as before. γE , absent in a tokamak,provides the restoring force toward the point E r = Eaof ambipolarity. If D(ω) is ω-independent, then γ E is aswell, and in the time domain Eq.(17) reduces to a standard
Langevin equation for E,∂t E(t) + γE E(t) = cS (t).

(18)

The source cS that drives the ZFs is approximated as random. Thus, ensemble averaging (18), one has∂t hEi p = −γE hEi p ,(19)R
where h...i p ≡ dE...p(E, t) is the ensemble average withprobability distribution function (pdf) p(E, t). It satisfies∂t p = ∂E (DE ∂E p + γE E p),(20)
R∞
with DE ≡ 0 dτhcS (t)cS (t − τ)i p the diffusion coefficientin E-space. From this follows Eq.(19), and1∂t hE 2 i p = DE − γE hE 2 i p .
2(21)

One notes from this the balance between diffusion
and the restoration toward E r = Ea . In steady-state,Eqs.(20) and (21) yield p(E) = p0 exp(−γE E 2 /2DE ), andhE 2 i p = DE /γE . Since γE ∼ D−1 and DE ∼ D−2 , one
has hE 2 i p ∼ D−1 . Thus, assuming the turbulent forces
FS driving the ZFs are unaffected, the larger D implied
at low-ω by the drift-polarization shielding would reduce
γE , but reduce the diffusion DE even more, resulting in asmaller ZF amplitude hE 2 i1/2
p .
Discussion
In the tokamak limit σ, γE → 0, Eq.(18) or (21)predicts an unbounded diffusion. Then other restoringmechanisms, such as those given in the model Eq.(19) ofRef. [5], become important, and would provide an analogous bounded statistical evolution of E r , though for thatmodel equation the time-average value of E r would shift
from the stellarator value E a to 0.
The more robust
ambipolar field E a in a stellarator provides a turbulencesuppression mechanism additional to, enhancing or diminishing that of, the ZFs themselves. For example, internal transport barriers induced by jumps in E a (r) from theion to the electron root[21, 22, 23] (which enhances theshear in Er and resultant flow-shear) have been observedon W7AS[24], LHD[25], and on CHS[26]. Better understanding how the ambipolar and ZF-induced flow shear act
together presents an additional issue, and potential opportunity, for stellarators.
The theory from Ref. [1], as for earlier work[5,6, 7, 8] which it extends, treats one leg of the selfconsistent ZF loop, the time evolution of zonal flows,
given a specified turbulent source S . The relationship between the transport coefficients D j and shielding function
contributions g j established by that theory indicates thatneoclassically-optimized stellarators should have less ZF
damping, and thus supports the earlier view[11, 6, 12, 10]that neoclassically-optimized stellarators should also havelower turbulent transport. However, as noted here, the assumption that S remains fixed from one configuration toanother is unwarranted, and further study of the dependence of S on configuration is clearly indicated. In fact, therecent results[10, 20] from LHD simulations indicate that
this widely-held view is too simplistic, and that a fuller
perspective will include additional variables beyond only
the degree of a configuration’s neoclassical optimization.
Understanding the configuration-dependence of S isof course complicated, possibly the reason most analytic
work has focussed mainly on the first leg of the ZF loop.
Theoretical progress here can be aided greatly by analysis
of further numerical simulations designed to elucidate thisrelationship.
[1] H.E. Mynick, A.H. Boozer, Phys. Plasmas 14 072507(2007).
[2] P.H. Diamond, Y.-B. Kim, Phys.Fluids-B 3, 1626 (1991).
[3] H. Biglari, P.H. Diamond, P.W. Terry, Phys.Fluids-B 2, 1
(1990).
[4] L. Chen, Z. Lin, R.B. White, Phys. Plasmas 7, 3129 (2000).
[5] M.N. Rosenbluth, F.L. Hinton, Phys. Rev. Letters 80, 724(1998).
[6] H. Sugama, T.H. Watanabe, Phys. Rev. Letters 94, 115001(2005).
[7] H. Sugama, T.H. Watanabe, Phys. Plasmas 13, 012501(2006).
[8] K.C. Shaing, Phys. Plasmas 13, 052505 (2006).
[9] A.N. Kaufman, Phys. Fluids 15, 1063 (1972).
[10] T.-H. Watanabe, H. Sugama, S. Ferrando-Margalet, Nucl.
Fusion 47, 1383 (2007).
[11] H. Yamada, J.H. Harris, A.Dinklage, E. Ascasibar, F.Sano,
S.Okamura, U.Stroth, A.Kus, J.Talmadge, S.Murakami,
M.Yokoyama, C.Beidler, V.Tribaldos, K.Y.Watanabe, Proc.
31st EPS Conference on Plasma Physics, London, JuneProceedings of ITC/ISHW2007
29-July 2, 2004 (European Physical Society, Pepip-Lancy,2004), P5-9.
[12] K.C. Shaing, Physics of Plasmas 12, 082508 (2005).
[13] A. Pytte, A.H. Boozer, Phys. Fluids 24, 88 (1981).
[14] A.H. Boozer, Phys. Fluids 23, 2283 (1980).
[15] R. J. Goldston, R. B. White and A. H. Boozer, Phys.Rev.
Lett. 47, 647 (1981).
[16] R. Linsker, A.H. Boozer, Phys. Fluids 25, 143 (1982).
[17] Y. Xiao, P.J. Catto, Phys. Plasmas 13, 102311 (2006).
[18] D.R. Mikkelsen, H. Maassberg, M.C. Zarnstorff, C.D. Beidler, et al., Fusion Science and Technology 51 166-180(2007).
[19] H.E. Mynick, Phys. Plasmas 13, 058102 (2006).
[20] H. Sugama, T.-H. Watanabe, S. Ferrando-Margalet, (thismeeting, paper I-08).
[21] H.E. Mynick,W.N.G. Hitchon, Nucl. Fusion 23, 1053
(1983).
[22] K.C.Shaing, Phys. Fluids 27, 1567 (1984).
[23] D.E.Hastings, W.A. Houlberg, K.C.Shaing, Nucl. Fusion
25, 445 (1985).
[24] U Stroth, K Itoh, S-I Itoh, H. Hartfuss, H. Laqua, the ECRH
team, & the W7-AS team, Phys. Rev. Letters, 86, 5910
(2001).
[25] K. Ida, T. Shimozuma, H. Funaba, et al., Phys. Rev. Letters
91, 085003 (2003).
[26] T. Minami, A. Fujisawa, H. Iguchi, et al., Nucl. Fusion 44,
342 (2004).