Fusion Res.
Res. SERIES,SERIES, Vol.
Vol. 6 (2004) 000-000
3–7
J. Plasma Fusion

Drift Waves and Zonal Flows
BALESCU Radu
Association Euratom - Belgian State,
Physique Statistique et Plasmas, CP 231
Université Libre de Bruxelles,
B-1050, Brussels, Belgium(Received: 10 December 2003 / Accepted: 27 January 2004)

Abstract
Drift wave turbulence in a plasma is analyzed on the basis of the Wave Liouville Equation, describing theevolution of the distribution function of wave packets (quasiparticles) characterized by position x and wave vector k .
A closed kinetic equation is derived for the ensemble-averaged part of this function; it has the form of a nonMarkovian advection-diffusion equation describing coupled diffusion processes in x - and k -spaces. General forms ofthe diffusion coefficients are obtained in terms of Lagrangian velocity correlations. The latter are calculated in thedecorrelation trajectory approximation, duly accounting of the trapping of quasiparticles in the rugged electrostatic
potential. The analysis of individual decorrelation trajectories provides an illustration of the fragmentation of drift
wave structures in the radial direction and the generation of long-wavelength structures in the poloidal direction thatare identified as zonal flows.
Keywords:
plasma turbulence, zonal flow, wave kinetic equation, non-Markovian diffusion

1. Introduction

circumstances, the drift wave turbulence is able to generatespontaneously a large scale poloidal shear flow. This largescale turbulent poloidal flow is called a zonal flow. In order
to avoid confusion, let us stress that in the present work we
use a slightly more general definition of the term “zonal
flow” compared to other authors. We understand under thisterm a structure with poloidal wave vector k y much smaller
than the radial one k x , but not ecessarily zero. Its effect isagain, understandably, a tearing apart of the drift wave
structures and their fragmentation, strikingly seen in themassively parallel numerical simulations published in Ref.
[4]. The most elaborate explanation of this effect was
initiated by Diamond and his collaborators [5], and is still
being actively pursued (see e.g., Refs. [2,6,7]).
In the present paper we are mainly interested in themechanism of generation of the zonal flows. The “classical”
discussion of this problem is exposed in Ref. [2]. Its
argument can be summarized as follows. The fragmentation
produced by the complex zonal flows results in a diffusive
increase of the mean square of the radial wave-vector k x2 ,which is a measure of the inverse square of the correlation
length m -x 2 in the radial direction (i.e., the direction of thedensity gradient). As a result, the drift wave frequency~ d (k) = k y V) /( 1 + t 2s k 2 ) will decrease. (Here V) = t s c s /L n
is the electron diamagnetic velocity, t s the ion-acoustic
Larmor radius, c s the ion-acoustic velocity and L n the density

We consider a plasma in presence of a strong magnetic
field (such as a tokamak). Locally this magnetic field can beconsidered as constant, equal to B. A Cartesian reference
system is then defined by three unit vectors: e z directed along
B, e x in the direction of the density gradient, mimicking theradial direction in a tokamak (the temperature is taken
constant, for simplicity), and e y = e z # e x , representing thepoloidal direction (shearless slab geometry). If B issufficiently strong, the collisionless motion of the particles
can be approximated by the quasi-two-dimensional motion of
their guiding centres, characterized by their position vectorx = (x, y) . The plasma is in a turbulent state, produced by afluctuating electrostatic potential denoted by fz (x, t) , wherez is a dimensionless function. It is well known that in atokamak, under certain circumstances, a transport barrier
may appear. At its location the turbulence level and thetransport are significantly reduced.
An enormous amount of literature on transport barriers
has appeared in the past twenty years; the main results arereviewed in Refs. [1-3]. An important mechanism of
formation of such a barrier is produced by the action of ashear flow in the direction perpendicular to the magnetic field
and to the gradients (i.e., the poloidal direction in a tokamak).
This flow will carry along nonuniformly the drift wave
structures, thus deforming, and finally breaking them into
smaller substructures. It will be shown that under certain©2004 by The Japan Society of Plasma
Science and Nuclear Fusion Research

Corresponding author’s e-mail: rbalescu@tiscalinet.be13Balescu R., Drift Waves and Zonal Flows

gradient length; k y is the y-component of the wave vector k ,mimicking the poloidal component k i in a tokamak,k 2 = k x2 + k y2 .) It follows from the conservation of the action
density N (k) . e (k)/ ~ d (k) that the drift wave energy e (k)must also decrease, it being transferred to the large-scale
zonal flow energy. We are thus in presence of an instability
that produces an inverse cascade in k y -space. After its
saturation, a new state appears, in which the level of smallscale turbulence, hence the anomalous radial transport issignificantly reduced: a transport barrier has been created.
This process is described in Ref. [2] by a diffusion
coefficient in k -space, estimated in the quasilinear
approximation. The authors of Ref. [2] correctly note that
their quasilinear treatment of k -space diffusion is of limited
validity. When the turbulence level is not very weak, newfeatures become important. The rugged fluctuating potential
landscape (as it appears in the numerical simulations of Ref.
[4]) produces transient trapping of the particle trajectories in
its troughs or around its peaks, as well as formation ofcoherent structures [7]. The difficult nonlinear treatment ofthese processes has been approached in various ways in the
recent literature; it will not be further discussed here.
In the present paper we study the trapping effect on the
global transport. We base our treatment on the decorrelation
trajectory (DCT) method, introduced and developed inseveral works, starting in 1998 [8], and extended here inorder to include the random walk in k -space. The method
and the main results obtained sofar will be outlined here; thedetails of the calculations and additional references will befound in a recent publication [9].

2 t N = L W N / 7 H, N A .

The action density is a fluctuating quantity in the largescale domain; hence the Wave Liouville Equation is astochastic differential equation, whose associatedcharacteristic equations are typical Langevin equations.
The distribution function is separated into an averagepart and a fluctuating part: N (x, k, t) = n (x, k, t) + dN (x, k, t)
where n (x, k, t) = N (x, k, t) . Our first purpose is to derivea closed equation for the average distribution function. The
strategy is standard [9,10]: the equation for the fluctuation issolved by the method of characteristics, and the result is
substituted into the equation for the average, which becomes:

2 t n (x, k, t) + V g (k) $ 2 n (x, k, t)
2x
t
)
)
2
2
XX
= # dt1 ( $ L ^ t - t1h $+ 2 $ LXK ^ t - t1h $ 2 (4)
2x
2x 2x
2k
0
)
)
2
2
2
2
KX
KK
+ $ L ^ t - t1h $+$ L ^ t - t1h $ 2 n (x, k, t 1 ),2x 2 k2k2kwhere V g (k) = 2~ d (k)/ 2k is the (deterministic) groupvelocity of the drift waves. Equation (4) will be called the(true) Wave Kinetic Equation. It is a non-Markovian“advection-double-diffusion” equation for the average
distribution n (x, k, t) , describing two coupled diffusion
processes, respectively in x -space and in k -space, combinedwith propagation (advection) in x -space. Equation (4)contains the precise formulation of what Diamond and his
group call the “random walk in k -space”.
Equation (4) introduces four 2 # 2 Lagrangian
correlation tensors of the fluctuating x -velocities V and thek -velocities W / - 2 7 k $ V (x, t) A /2x . Typically:

2. The wave kinetic equation.
Of particular importance for the characterization of the
turbulent state is the action density N (x, k, t) , which is aconserved quantity. It is a quadratic functional of thepotential, averaged over small scales (with respect to t s ),and slowly varying in x -space. Its precise definition, in thecase of drift waves, was derived in Ref. [6], together with its
equation of evolution:

2 t N (x, k, t)
2H (x, k, t) 2N (x, k, t) 2H (x, k, t) 2N (x, k, t)
=
$
$
,
2k
2k
2x
2x

(3)

LKK
r | s (t - t 1 ) = W r 7 x, k, t A W s 8 x ^ t 1h, k ^ t 1h, t 1 B ,(5)and similar definitions for the three other tensors. Assuming
stationary turbulence, these functions only depend on the
difference of the two times. If Eq. (4) could be
Markovianized by neglecting the memory effect, i.e., settingn (x, k, t 1 ) . n (x, k, t) , it would become an “ordinary”
multidimensional advection-diffusion equation, with thefollowing coefficients, which are time-dependent, running
diffusion coeficients:

(1)

with:

t

H (x, k, t) = ~ d (k) + k $ V (x, t),

D rKK| s (t) = # dt 1 LKK
r | s (t 1 ),(2)0

where V = (fc/B)( e z # dz) is the (long-range part of the)electric drift velocity. The remarkable feature here is the
Hamiltonian structure of the evolution of a set of “quasiparticles”, i.e., wave packets, characterized by the
canonically conjugate variables ( x, k ). The action density
plays the role of the phase space distribution of the quasiparticles, and obeys the Liouville equation (1), which can be
written in terms of a Liouville operator LW , defined as instatistical mechanics by the Poisson bracket of H and N :(6)and similar definitions for the three other tensors. The
asymptotic values D rKK| s (3) / D rKK| s are the ordinary diffusion
coefficients. The Markovianization process is, however, notalways justified [10].

3. The DCT method for zonal flows
An essential role in the nature of strongly turbulent
phenomena is played by the trapping processes in the rugged
potential landscape. It is therefore natural to study these
phenomena by using the decorrelation trajectory (DCT)24

Balescu R., Drift Waves and Zonal Flows

and its first derivatives (through v ), but also its second
derivatives (through w ) enter the theory. The subensemble S
is therefore defined as the set of realizations in which thepotential z, the two components of the velocity v r , and thetwo components of the “k -velocity” w r have a given value attime 0:method, developed for “ordinary” drift wave turbulence inRef. [8]. This method requires an extension in the present
case.
We study the stochastic motion of a drift-wave
quasiparticle in presence of the fluctuating electrostatic
potential fz (x, t) . The fluctuations are defined statistically bya homogeneous and stationary Gaussian process, whose
Eulerian correlation introduces a correlation length m and acorrelation time x c , such that G z (0, 0)~ z (x, t) H . 0 whenexe & m and/or t & x c . We go over to dimensionless variables
defined as: x ( mx , k ( t -s 1 k , t = x c i . The Langevinequations of motion of the quasiparticle are written as:

dx ( i)= K d v g [k ( i)] + K v [x ( i), i],didk ( i)= Kw [x ( i), k ( i), i],dix (0) = 0,

k (0) = k 0 .

S: z (0, 0) = z 0, vr (0, 0) = vr0 , - k s0 d r vs (0, 0) = w r0 . (12)Let P0 (z 0 , v 0 , w 0 ) be the probability distribution (pdf)of these initial values, assumed to be Gaussian. Then thex - x component of the KK -Lagrangian correlation (10) is
provided (exactly) by the superposition:
0
0
0
0
0
0
0
LKK
x | x (i ) = # dz dv dw P0 (z , v , w ) w x(7)Sw x [x (i ), k (i ), i ] ,(8)where G ... HS denotes the average in the subensemble. Notethe important fact that in the Lagrangian correlation in S :G w x [0, k 0 , 0] w x [x (i), k (i), i] HS the first factor can be takenout of the average because of Eq. (12), hence the calculationof this quantity is reduced to the simpler calculation of theaverage Lagrangian velocity G w x [x (i), k (i), i] HS .
The pdf P0 as well as the Eulerian averages of thevelocities v and w in S can be calculated analytically as
appropriate conditional averages. Typically, the Eulerian
average k -velocity in the subensemble is of the form:G w (x, k, i) HS = w S (x, k) x (i) , with:(It is important to take an initial wave vector of non-zero
length: an initial | k 0 | = 0 would remain zero for all times.)
The three scaled velocities v g , v and w are related to thecorresponding dimensional quantities as follows:V g = (m/x c ) K d v g , V = (m/x c ) K v , W = (t s x c ) - 1 K w . We
introduced here the two basic dimensionless parameterscharacterizing the turbulence: the Kubo number K , related tothe intensity f of the fluctuations, and the diamagnetic Kubo
number, related to the gradient length through V) :K = x c cf ,m Bm

K d = x c V) .
mw S (x, k) = A z 0 + B $ v 0 + C $ w 0

(9)

(10)

We first note that all fluctuating quantities (including v
and w ) are derived from the potential z. As in all theories
based on Langevin equations, the primary Eulerian potential
autocorrelation has to be specified a priori. This quantity is
assumed to be of the same form as in Ref. [8]:z (0, 0) z (x, i ) = E (x) T (i ),(14)

where A, B, C are rather complicated functions of theEulerian correlation E and of its derivatives [9]; v S has asimilar form. Thus, w S and v S depend on the parametersz 0 , v 0 , w 0 defining the subensemble S , as well as on theinitial value of the wave vector k 0 .
The decorrelation trajectory (DCT) is defined as adeterministic trajectory [x S (i), k S (i)] in the subensemble S
of a fictitious quasiparticle whose motion is determined by
Eqs. (7), (8) in which v and w are replaced by v S and w S ,respectively [8]. Finally, the DCT approximation of the
Lagrangian correlations is obtained by replacing in Eq. (13)the subensemble Lagrangian average G w x [x (i), k (i), i] HSby the Eulerian average in S of the velocity, evaluated foreach time at the deterministic point [x S (i), k S (i)] :The transport problem requires the evaluation of the four
Lagrangian correlations entering Eq. (4). We concentratehere on the last one, which introduces all the new features;the three others are treated in a completely similar fashion. Its
corresponding dimensionless form, for a homogeneous andstationary turbulent state, is:0
LKKr | s ( i ) = w r (0, k , 0 ) w s [ x ( i ), k ( i ), i ]

(13)

0
0
0
0
0
0
0
LKK
x | x ( i ) = # d z dv dw P0 ( z , v , w ) w x(11)

w xS [x S (i ), k S (i )] x (i ).

where E (x) is a dimensionless function of the position, witha maximum at the origin, and T (i) a similar function of time.
All other Eulerian correlations of couples of fluctuating
quantities ( z, v, w ) are deduced from (11) by appropriate
differentiations [9].
The extension of the decorrelation trajectory (DCT)
method follows the line of Ref. [9]. The basic step is thedecomposition of the ensemble of realizations of the
turbulent ensemble into subensembles. This decomposition is
now different from Ref. [8], because not only the potential z(15)

The calculation of the Lagrangian correlations is thus
replaced by the calculation of Eulerian averages. The present
problem is, however, significantly more complicated in thepresent case than for the simple drift wave turbulence treatedin Ref. [8] because of the intimate coupling of the equations
for x S (i) and k S (i) . The Lagrangian correlations are now 5fold integrals, which makes their numerical evaluationdifficult and time-consuming. The final result depends on
more parameters, viz. the Kubo number K , the diamagnetic
3
5

Balescu R., Drift Waves and Zonal Flows

Kubo number K d , and also the initial wave vector k 0 . Some
important qualitative features can be obtained from ananalysis of the individual DCT trajectories.

direction. Note that the increase of k xS is monotonous, until itreaches a saturation value after i . i tr . The y- (“poloidal”)
component k yS quickly reaches a maximum: its monotonous
growth is then stopped and reversed after a certain time
( i . 1). It then changes sign, and finally ( i . i tr ) reaches anegative saturation value. This is the manifestation of thetrapping process in k -space, which only affects (in thepresent situation) the y-component. This process ends afteri . i tr ,
kS
after
which
remains
constant:S
3S
k (i) " k [and w (i) " 0] . It is important to note that (inthe present case) | k y3 |< k x3 ( k x3 = 1.54, k y3 = - 0.047 ). As aresult, in the asymptotic state ( i > i tr ), the average lengthscale of the turbulence in the radial direction is much smaller
than in the poloidal direction. This obviously explains thefragmentation process described qualitatively in the firstsection, and implies the generation of a (small-k y ) zonal flow
(in the subensembleS ).
To sum up: during the trapping time i tr , the wave vector
changes from its initial value to a constant asymptotic value
with | k y3 |< ek x3 e ; the fictitious quasiparticle is trapped duringthat time, and ends up moving with a new, deflected groupvelocity (refraction). In the final state, the length scale of the
turbulence in the y-(poloidal) direction is much greater thanin the x-(radial) direction.
This analysis shows that the decorrelation process is
now richer than in the case studied in Ref. [8]. In particular,
the position and the wave vector move in a strictly coupled
way and cannot be considered separately from each other.
There is weak trapping in x -space, but the trapping of the k ywave vector component is already quite significant at thisrelatively modest value of the Kubo number ( K = 0.5 ).
We consider next a situation of strong turbulence,
K = 10 , with the same subensemble parameters. The features
which were merely sketched in the case K = 0.5 are now
greatly enhanced. The (fictitious) quasiparticle is strongly
trapped in x -space: its motion is a “broken rotation” during
the time i tr (which has barely changed, i tr . 3 ). During thetrapping time it is more and more strongly deflected. Its
trajectory makes two turns, thus confining the quasiparticle toa finite region of space. But after i tr , the turbulent
component of the velocity vanishes, and the fictitious
quasiparticle moves away with the final group velocityv g (k 3 ) . Meanwhile, both the k xS - and the k yS -components of
the wave vector undergo a broken oscillation which stops at
i . i tr , when they reach the values: k x3 = 0.879 , k y3 = 0.965 .
We thus witness strong trapping in both k . In the present
subensemble we can no longer speak of zonal flow
generation: k y3 is now of order one, and moreover it
surpasses k x3 in absolute value. The strong turbulence
produces in this case a fragmentation both in the radial and in
the poloidal directions.
Clearly, the present discussion refers to a single
subensemble. Recently we obtained, in collaboration with I.
Petrisor, some preliminary results of a global nature (theseresults are not contained in Ref. [9]: they will be the object of
a forthcoming publication). Thus, the Lagrangian correlation

4. Results and conclusions
We now consider the result of the numerical integrationof the decorrelation trajectories. For definiteness, we assume
the following form for the Eulerian potential autocorrelation (11):E (x, y) = exp (-x2+ y2), x (i ) = exp (- i )2(16)

We define a subensemble S by the following values ofthe parameters: z 0 = 2 , v x0 = v y0 = 0 , w x0 = w y0 = 1. We alsochoose the fixed value K d = 1. The choice v 0 = 0 implies thatthe (fictitious) quasiparticle starts at time zero with the initialgroup velocity v g (k 0 ) and ends after a sufficiently long time
( i & i tr ), when v nS = 0, w nS = 0 , with the final group velocityv g (k 3 ) . Note that, because of the factorization property
G w (0, 0) w (x, i) t HS = w 0 G w (x, i) HS (and similarly for v )[see Eq. (13)], the vanishing of v S , w S implies the vanishingof the Lagrangian velocity correlation in the subensemble S .
Thus the trapping time i tr is the time after which thefictitious quasiparticle is no longer correlated along its
trajectory with its initial value. The trapping time i tr isdetermined numerically from the shape of the trajectories. (Itshould be clear that i tr relates to a single DCT.) Because of
the limited amount of space available, we cannot reproduce
here the figures representing the characteristic features of the
DCT: they are found in Ref. [9]. We shall therefore
summarize the main facts in words.
The initial wave vector is chosen as: k x0 = 1, k y0 = 0 . As aresult, the (unperturbed) group velocity at the initial time is:v g (k 0 ) = (1/2) e y . Thus, in absence of turbulence (K = 0 ),the fictitious quasiparticle moves in a straight line in they -direction, with a group velocity that is constant, because
the wave vector k remains constant, Eq. (8).
We now consider a rather small value of the Kubo
number, K = 0.5 . Even in this relatively weak turbulence, thepicture departs radically from the unperturbed motion. The
beginning of a trapping process in x -space is evident: the(fictitious) quasiparticle starts with the initial group velocityv g (k 0 ) , but the turbulent velocity quickly overcomes thelatter and deflects the particle from its rectilinear motion. The
turbulent component v of the velocity v S vanishes after atime that will be called the trapping time i tr , of order of thecorrelation time; in the present case its value is seen to be:i tr . 3 , after which the quasiparticle moves again uniformly,
with its group velocity. Meanwhile, the latter has changed
because the wave vector has changed. Thus the quasiparticle
is deflected from its unperturbed motion in the y -direction
and moves now in an oblique direction. We thus witness herea refraction phenomenon due to the passage through the
turbulent medium. This refraction effect is also found (in adifferent context) in Ref. [7].
The wave vector starts from k 0 and initially increases inboth directions, the growth being largest in the x- (“radial”)
46

Balescu R., Drift Waves and Zonal Flows

Acknowledgmentsof the radial k -velocities LKK
x | x (i) , Eq. (10) was calculated fortwo values of the Kubo number. For very weak turbulence,
K = 0.1, it is a function decreasing monotonously in timefrom a positive value towards zero. This is a typicalbehaviour in the quasilinear regime. The corresponding
running diffusion coefficient D xKK| x (i) , Eq. (6) is a positive
monotonously increasing function of time, reaching a finite
asymptotic value. For strong turbulence, K = 7 , the behavior
is quite different. The Lagrangian correlation decreases
rapidly from an initial positive value, becomes negative,
reaches a minimum at i . 0.4 , after which it performs abroad oscillation, and ends up with a long negative powerlaw tail reaching zero. The corresponding running diffusion
coefficient rises quickly to a maximum, and then decreases toa much smaller asymptotic value (remaining always positive,
as it should). These features are typical signatures of thetrapping effect.
In conclusion, we showed that the quasilinear
approximation, valid in the weak turbulence regime, leads toseriously overestimated values of the diffusion coefficients in
the regimes of large Kubo number. The decorrelation
trajectories indicate a fragmentation of the radial structures
and a spontaneous generation of zonal flows in particularsubensembles. A global study of the coupled advectiondiffusion process of drift waves requires a detailed scanningof the parameter space. This task will be the subject offorthcoming work.

We gratefully acknowledge very fruitful discussions
with M. Vlad, I. Petrisor, M. Negrea, and J. Misguich.

References[1] P. Terry, Rev. Mod. Phys. 72, 109 (2000).
[2] P.H. Diamond, S. Champeaux, M. Malkov, A. Das, I.
Gruzinov, M.N. Rosenbluth, C. Holland, B. Wecht, A.
Smolyakov, F.L. Hinton, Z. Lin and T.S. Hahm, Nucl.
Fusion 41, 1067 (2001).
[3] A. Yoshizawa, S-I. Itoh and K. Itoh, Plasma and Fluid
Turbulence, Theory and Modelling (Inst. of Phys.,
Bristol, 2003).
[4] Z. Lin, T.S. Hahm, W.W. Lee, W.M. Tang and R.B.
White, Science 281, 1835 (1998).
[5] H. Biglari, P.H. Diamond and P. Terry, Phys. Fluids B
2, 1 (1990).
[6] A.I. Smolyakov, P.H. Diamond and M. Malkov, Phys.
Rev. Lett. 84, 491 (2000).
[7] P. Kaw, R. Singh and P.H. Diamond, Plasma Phys.
Control. Fusion 44, 51 (2002).
[8] M. Vlad, F. Spineanu, J.H. Misguich and R. Balescu,
Phys. Rev. E 58, 7359 (1998).
[9] R. Balescu, Phys. Rev. E 68, 046409 (2003).
[10] R. Balescu, Plasma Phys. Control. Fusion 42, B1 (2000).

57