Proceedings of ITC/ISHW2007

Onset of Turbulence in a Drift Wave–Zonal Flow System
Ryusuke NUMATA1) , Rowena BALL1,2) , Robert L. DEWAR1) and Linda. STALS2)
1)

Department of Theoretical Physics, Research School of Physical Sciences & Engineering, The Australian National
University, Canberra, ACT 0200, Australia
2)Mathematical Sciences Institute, The Australian National University, Canberra, ACT 0200, Australia
Numerical analyses of bifurcation phenomena in the Hasegawa-Wakatani model are presented, that provide
new insights into the interactions between turbulence and zonal ﬂows in the tokamak plasma edge region. Thesimulation results show a regime where, after an initial transient, drift wave turbulence is suppressed throughzonal ﬂow generation. As a parameter controlling the strength of the turbulence is tuned, this zonal-ﬂowdominated state is rapidly destroyed and a turbulence-dominated state re-emerges. The transition is explainedin terms of the Kelvin-Helmholtz stability of zonal ﬂows. This is the ﬁrst observation of an upshift of turbulence
onset in the resistive drift wave system, which is analogous to the well-known Dimits shift in turbulence driven
by ion temperature gradients.
Keywords: bifurcation, transition, drift wave, turbulence, zonal ﬂow

1Introduction

Fusion plasmas and other turbulent ﬂows in quasi-twodimensional (2D) geometry can undergo spontaneous transitions to a turbulence-suppressed regime. Known as L–H(low-to-high conﬁnement) transitions in plasmas, they arestudied intensively because conﬁnement enhancement may
supervene through the concomitant abatement of anomalous or turbulent particle and heat ﬂuxes. In tokamakedge plasmas L–H transitions are associated with nonlinearly self-generated poloidal E × B shear or zonal ﬂows,
which absorb energy from drift waves and consume thesmall scale eddies that mediate turbulent transport. It is
now widely accepted that control of emergent zonal ﬂows
is crucial to achieving and sustaining improved conﬁnement [1].
In this paper we present the results of analytic and numerical investigations of transitions between turbulencedominated and zonal-ﬂow-dominated regimes, using theHasegawa–Wakatani (HW) model [2, 3] for electrostatic
resistive drift wave turbulence in 2D slab geometry. We
ﬁnd that bifurcations in the model correspond to the onset of drift wave driven turbulence, generation of zonal
ﬂows, and re-emergence of drift wave turbulence as thezonal ﬂows become unstable. This latter phenomenon is
analogous to the Dimits shift [4] described for turbulence
driven by ion temperature gradients (ITG).

2

Hasegawa–Wakatani model

The physical setting of the HW model may be consideredas the edge region of a tokamak plasma of nonuniform density n0 = n0 (x) and in a constant equilibrium magnetic
ﬁeld B = B0 ∇z. Following the drift wave ordering [5],author’s e-mail: ryusuke.numata@anu.edu.authe ion vorticity ζ ≡ ∇2 ϕ (ϕ is the electrostatic potential,∇2 = ∂2 /∂x2 + ∂2 /∂y2 is the 2D Laplacian) and the density
ﬂuctuations n are governed by the equations∂ζ + {ϕ, ζ} = α(ϕ − n) − D∇4 ζ,∂t∂
∂ϕn + {ϕ, n} = α(ϕ − n) − κ− D∇4 n,∂t
∂y(1)
(2)

where {a, b} ≡ (∂a/∂x)(∂b/∂y) − (∂a/∂y)(∂b/∂x) is thePoisson bracket, D is the dissipation coeﬃcient. The background density is assumed to have an unchanging exponential proﬁle: κ ≡ (∂/∂x) ln n0 . Electron parallel motion isdetermined by Ohm’s law with electron pressure pe = nT e ,jz = −enve,z = −)1 ∂ (Teϕ−ln n ,
η ∂ze

(3)

assuming electron temperature T e to be constant (isothermal electron ﬂuid). This relation gives the coupling between ζ and n through the adiabaticity operator α ≡−T e /(ηn0 ωci e2 )∂2 /∂z2 appearing in Eqs. (1) and (2). In our
2D setting α becomes a constant coeﬃcient when acting
on the drift wave components of ϕ and n by the replacement ∂/∂z → ikz , where 2π/kz = Lk  Ly is a lengthcharacteristic of the drift waves’ phase variation along theﬁeld lines. However, since zonal components of ﬂuctuations (ky = kz = 0 modes) do not contribute to the parallel
current, this resistive coupling term must be modiﬁed withcare [6]. Recalling that turbulence in the tokamak edge region, where there is strong magnetic shear, is consideredhere, ky = 0 should always coincide with kz = 0 because
any potential ﬂuctuation on the ﬂux surface is neutralizedby parallel electron motion. Let us deﬁne zonal and nonzonal components of a variable f as∫1
f dy, non-zonal: f˜ = f − h f i,zonal: h f i =Ly

Proceedings of ITC/ISHW2007
where Ly is the periodic length in y, and remove the contribution by the zonal components in the resistive coupling
term in Eqs. (1) and (2). Subtraction of the zonal components from the resistive coupling term α(ϕ − n) → α(ϕ̃ − ñ)
yields the modiﬁed HW (MHW) equations,(4)(5)

The variables in Eqs. (4) and (5) have been normalized byx/ρs → x, ωci t → t, eϕ/T e → ϕ, n1 /n0 → n,√
where ρs ≡ T e /mω−1
ci is the ion sound Larmor radius√(vsi ≡ T e /m is the ion sound velocity in the cold ionlimit), n1 is the ﬂuctuating part of the density.
Wakatani and Hasegawa found [3] that excitations ofwaves having kz that maximizes the linear growth rate (forgiven k x and ky ) are most likely to occur, since the plasma
can choose any parallel wavenumber (kz ). Using the parallel wave number of the maximum growth rate, α is given
by α = 4k2 ky κ/(1 + k2 )2 . This also gives α = 0 for thezonal mode.
The MHW model spans two limits with respect tothe adiabaticity parameter α. In the adiabatic limit
α → ∞ (collisionless plasma), the non-zonal component of electron density obeys the Boltzmann relationñ = n0 (x) exp(ϕ̃), and the equations are reduced to theHasegawa–Mima equation [5]. In the hydrodynamic limit
α → 0, the equations are decoupled. The vorticity is determined by the 2D Navier-Stokes equation, and the density
becomes a passive scalar. The advantage of our choice of
α as a free parameter is the capability for treating the limits
in a uniﬁed manner.
In the adiabatic, ideal limit (α → ∞, D → 0) theMHW system has two dynamical invariants, the energy E
and the potential enstrophy W,∫∫
11
E=(n2 + |∇ϕ|2 )dx, W =(n − ζ)2 dx, (6)22where dx = dxdy, which constrain the ﬂuid motion. Conservation laws are given bydE= Γn − Dα − DE ,dt

dW= Γn − DW ,dt

where ﬂuxes and dissipations are given by∫∂ϕ̃Γn = −κ ñ dx,∂y∫Dα = α (ñ − ϕ̃)2 dx,∫DE = D ((∇2 n)2 + |∇ζ|2 )dx,∫
DW = D (∇2 n − ∇2 ζ)2 dx.

(7)

3Simulation Results

The MHW equations are solved in a doubly periodic
square slab domain with box size L = 2π/∆k where thelowest wavenumber ∆k = 0.15 (L ∼ 42). The equations
are discretized on 256 × 256 grid points by the ﬁnite diﬀerence method. We examine the eﬀects of the parameters κand α on the nonlinearly saturated state, and ﬁx D = 10−4
throughout this paper.
We start simulations by imposing small amplitude random perturbations. The perturbations grow linearly in theinitial phase and generate drift waves, then the drift waves
undergo secondary instabilities which excite zonal ﬂows
until nonlinear saturation occurs. In the saturated state, weobserve that Γn ' Dα  DE , DW . The spatial structureof the saturated electrostatic potential is shown in Fig. 1.
We observe that zonally elongated structures of the electrostatic potential are generated in the MHW model because
the modiﬁcation removes the unphysical resistive dissipation of the zonal modes. The zonal ﬂows carry nearly all
the kinetic energy in the ﬁnal state — they have absorbed
nearly all the energy from the drift waves. The build-upof the zonal ﬂow and resulting transport suppression highlight the importance of the modiﬁcation of the model in thenonlinear regime [7].

y/ρs [azimuthal]∂ζ + {ϕ, ζ} = α(ϕ̃ − ñ) − D∇4 ζ,∂t∂
∂ϕn + {ϕ, n} = α(ϕ̃ − ñ) − κ− D∇4 n.
∂t∂yUnlike the Hasegawa–Mima model which is an energyconserving system, the MHW model has an energy source
Γn . Due to the parallel resistivity, ñ and ϕ̃ can ﬂuctuate outof phase which produces non-zero Γn . The system can absorb free energy contained in the background density proﬁle through the resistive drift wave instability.

20151050-5-10-15-20-20-15-10 -5 0 5 10 15 20x/ρs [radial]Fig. 1 Contour plot of ϕ in the saturated state. Zonally elongatedstructure of the electrostatic potential is clearly visible inthe MHW model.

Let us show how the parameters aﬀect the saturated
state in the MHW model. In Fig. 2, we plot∫ the ratio of thekinetic energy of the zonal ﬂow (F ≡ 1/2 (∂hϕi/∂x)2 dx)

Proceedings of ITC/ISHW2007∫Kinetic Energy [Zonal/Total]

to the total kinetic energy (E k ≡ 1/2 |∇ϕ|2 dx) against
α. It is clearly seen that there are two types of saturatedstates: one where zonal ﬂows prevail and the other dominated by isotropic turbulence. Conversion from one state
to the other occurs over a narrow range of the parameter
space. We can see that zonal ﬂows are generated in theadiabatic regime (α  1) while isotropic ﬂows are generated in the hydrodynamic regime (α  1). Transition to
the turbulent state also occurs if the drift wave instability isdriven strongly by increasing the density gradient parameter κ.

the eigenvalue problem is rather simple because it is notnecessary to consider the continuous spectra on the realω axis. If we ﬁnd the eigenvalue ω and the corresponding eigenfunction ϕ̂, the complex conjugate of ω is also an
eigenvalue and the corresponding eigenfunction is givenby the complex conjugate of ϕ̂. Thus, we can always restrict our quest for eigenvalues in the upper half plane of
the complex ω plane, and can neglect interactions between
the point spectra and the continuous spectrum.
For the case of ﬁnite α the complex conjugate of an
eigenvalue is not a eigenvalue, therefore we must solve
for negative ωi too. Moreover, there exist two continuous
spectra in this case:ω = ky V, ky V − iα where |V| ≤ V0 .

κ=1
κ=0.1
-3

10-2

10-1

11010(9)Both represent convective transport due to the background
ﬂow. One of them is damped by the resistivity. These continua may interact with the point spectrum. Thus the situation is much more complicated in the intermediate α casecompared with the adiabatic and hydrodynamic limits.

α4Stability of Zonal Flow

We examine the stability of the zonal ﬂows obtained
from the numerical simulations, and compare the stability
threshold and the transition point in this section. We consider the perturbation around the zonal ﬂow background.
The electrostatic potential and the density are decomposed
as ϕ = ϕ0 (x) + ϕ̂(x) exp i(ky y − ωt), and n = n̂(x) exp i(ky y −ωt) where dϕ0 /dx = V gives the background ﬂow in the ydirection. By linearizing the MHW equations, we obtainan eigenvalue equation containing the eﬀect of κ and α,[ 2]
ky V 00d2ϕ̂−k+yω − ky Vdx2()
ky κiα−1−ϕ̂ = 0. (8)ω − ky V + iα
ω − ky V
We neglect the viscosity. The generated zonal ﬂows in they direction are assumed to have a sinusoidal proﬁle, V =V0 sin(λx). The amplitude V0 and wavenumber λ = nλ π/Lare determined from the simulation results as V0 ∝ κ2 , andλ ≈ 0.3. We solve the eigenvalue equation by the standard
shooting method in the domain D = {x| − L/2 ≤ x ≤L/2}. The boundary is assumed to be rigid ϕ̂(±L/2) = 0for simplicity.
In two limits (α → 0 and α → ∞), some conditions
for stability are known (see [8] and references therein), and0.03
0.02
Growth Rate ωiFig. 2 Parameter dependence of the zonal kinetic energy normalized by the total kinetic energy. Transitions froma zonal-ﬂow-dominated state to a turbulence-dominatedstate occur.

0.01
0-0.01
-0.02
-0.03
01

2

3ky L

456Fig. 3 Growth rates for HW case. Solid curves show the grawth
rates for α = 0.0001, and thick lines near the real ω axis
denote two continuous spectra. Two branches from theα → 0 case (dotted line) are also shown for reference.

We ﬁrst show the eﬀect of α and neglect eﬀect of κ.
We consider nλ = 2 for simplicity. Figure 3 shows theimaginary parts of the eigenvalues for α = 0.0001. Twobranches from the α → 0 case (dotted line) are also shown
for reference, so that it is seen that ωi is slightly shifted
downwards for ﬁnite α. For increasing α, we observe thepositive eigenvalues disappearing at α ≈ 0.000417.
Next, we consider the eﬀect of κ in addition to α.
Since κ always appears in the form of κα and α is small inthe vicinity of the threshold, the eﬀect of κ is rather minor.
κ does not signiﬁcantly aﬀect the behavior of the eigenvalues except that κ controls the amplitude of ﬂow (V0 ∝ κ2 ).
We summarize the stability of zonal ﬂows by showing the bifurcation diagram in α-κ plane together with thenumerically obtained results. The only excitable mode thatcan be resolved in the numerical simulation is the ky = 0.15

Proceedings of ITC/ISHW2007
mode, which is the ﬁrst unstable mode of the primary instability (resistive drift wave instability). In Fig. 4, weshow the stability threshold of ky = 0.15 mode for the primary instability and the tertiary instability (KH instability).
Each mark in the ﬁgure denotes a numerically obtained saturated state: N, , • represent respectively the zonal-ﬂowdominated, transitional, and turbulence-dominated states.
In these states zonal ﬂows contain more than 90%, 2090%, and less than 20% of the total kinetic energy, respectively. The qualitative tendency of the thresholds in thebifurcation diagram shows agreement between the numerical simulations and the KH analysis, i.e. increasing α (κ)
is stabilizing (destabilizing). Zonal-ﬂow-dominated states
are observed in between the primary and the tertiary instability thresholds. The emergence of a turbulent state isshifted from the primary threshold to the tertiary threshold
due to the turbulence suppression eﬀect of the zonal ﬂow,
which is analogous to the Dimits shift observed in ITG turbulence.
The reasons for the quantitative discrepancy between
the boundary of the zonal and the turbulent states may be
because of the simpliﬁcation made in the KH analysis; thesimpliﬁed ﬂow proﬁle, the boundary condition and viscosity may also aﬀect the results.

-4Bifurcation Diagram (D=10 )210κ

10ZonalTransitional
TurbulentPrimary
Tertiary

1-1

10-3

10-2

10-1

10

1

10

2

10αFig. 4 Bifurcation diagram showing the correlation between thelinearized stability estimates described in the text and theregimes observed in our turbulence simulations.

5ConclusionIn summary, we have analyzed bifurcation phenomena in
two-dimensional resistive drift wave turbulence. First, wehave performed numerical simulations of the modiﬁed HWmodel to study bifurcation structures in a two-parameter(α-κ) space. We have shown that, in the MHW model,
zonal ﬂows are self-organized and suppress turbulence andturbulent transport over a range of parameters beyond thelinear stability threshold for resistive drift waves. By performing a systematic parameter survey, we have foundthat such zonal-ﬂow-dominated states suddenly disappear

as a threshold is crossed, being replaced by a turbulencedominated state.
The threshold of the onset of turbulence has been compared with the linear stability threshold of an assumed laminar zonal ﬂow proﬁle. Numerical analysis of the eigenvalue problem determining the stability of the assumed
zonal ﬂow proﬁle in the HW model conﬁrms the following
trend: κ determines the amplitude of the zonal ﬂows, thus,
large κ destabilizes the zonal ﬂows. On the other hand, theadiabatic response of parallel electrons given by α stabilizes them. The constructed bifurcation diagram in the α-κplane for the HW model conﬁrms the scenario of the onsetof turbulence in the drift wave/zonal ﬂow system being due
to the disruption of zonal ﬂows by KH instability.

Acknowledgments
The authors would like to acknowledge B. D. Scott for providing the simulation code used in this work. We thankP. N. Guzdar, W. Dorland, C. Tebaldi, and J. A. Krommes
for useful discussions. This work is supported by the Australian Research Council.
[1] P. H. Diamond, S.-I. Itoh, K. Itoh, and T. S. Hahm, Plasma
Phys. Control. Fusion 47, R35 (2005).
[2] A. Hasegawa and M. Wakatani, Phys. Rev. Lett. 50, 682(1983).
[3] M. Wakatani and A. Hasegawa, Phys. Fluids 27, 611 (1984).
[4] A. M. Dimits, G. Bateman, M. A. Beer, B. I. Cohen, W. Dorland et al., Phys. Plasmas 7, 969 (2000).
[5] A. Hasegawa and K. Mima, Phys. Rev. Lett. 39, 205 (1977).
[6] W. Dorland and G. W. Hammett, Phys. Fluids B 5, 812(1993).
[7] R. Numata, R. Ball, and R. L. Dewar, in Frontiers in Turbulence and Coherent Structures, eds. J. P. Denier and J. S.
Frederiksen (World Scientiﬁc, Singapore, 2007), vol. 6 ofWorld Scientiﬁc Lecture Notes in Complex Systems, pp. 431–442, Proceedings of the COSNet/CSIRO Workshop on Turbulence and Coherent Structures in Fluids, Plasmas andNonlinear Media, Canberra, Australia, 10-13 January 2006.
[8] R. Numata, R. Ball, and R. L. Dewar, accepted for publication in Phys. Plasmas (2007).