Aspects of flow generation and saturation in drift-wave turbulence

The generation of large-scale flows by turbulence in plasmas and fluids has attracted a great deal of attention in the last years. The interplay between turbulence and sheared zonal flows has been identified as a possible candidate to explain the transition from L to H mode as well as the setup of transport barriers in general, while radial flows are made responsible for increased levels of transport. Using simplified equations for low-frequency plasma turbulence of the drift-wave type we consider some general aspects of flow generation and address the problem of saturation of drift-wave turbulence.


Introduction
Plasma turbulence, as a subset of general turbulence, is one of the most difficult and demanding fields in modern physics. For plasmas confined in magnetic fusion devices a complicated magnetic field geometry and an abundance of relevant time and spatial scales add to the complexity of the problem. Thus to obtain, under the given limitations of available computer power, solvable and tractable systems, first-principles modelling is most often sacrificed in favour of using experience from experiments. The models used to put numbers on quantities, such as the transport, thus incorporate much of the physics by using a parametrization for many of the fundamental, often nonlinear processes through the analysis of experiments.
On the other hand, the careful use-by means of the incorporated physics, geometry etc-of reduced first-principles models, often dubbed toy models, can raise understanding and physical insight. The latter though do not give conclusive answers to quantitative questions.
In this paper we shall focus on some qualitative aspects of drift-wave turbulence recoverable from simplified first-principles models. After a description of the linear and energetic properties of drift-waves and the basic dynamics of the turbulence generated, we shall in 28.2 section 3 focus on the topic of flow generation in terms of a competition between wave and turbulence effects. Section 4 will be devoted to the role of nonlinear vortical structures in the creation of zonal flows and section 5 will address the issue of saturation of drift-wave turbulence.

Drift-wave dynamics
We shall discuss some aspects of the generation of flows invoking one of the simplest models available describing three-dimensional plasma turbulence at the edge of a magnetic confinement device, the Hasegawa-Wakatani equations (HWEs) [1]: with {·, ·} being the Poisson bracket and ω = ∇ 2 ⊥ ϕ. The usual scaling makes the system dimensionless, but we additionally renormalized the fluctuations with the size of the inverse background gradient scale length κ n = 1 Ln . This reduces the number of free parameters by one, so that the parallel collisionality ν is the only free parameter beside the fluctuation amplitudes and the system size. For the sizes we distinguish between the directions parallel to the static magnetic field B = B 0ẑ aligned with the z-axis and the directions perpendicular to it. The background density gradient is assumed to point in the x-direction. The plasma is strongly magnetized and the drift approximation is used for the perpendicular plasma-fluid motion. Note that the electron temperature is assumed to be constant; normalized density and pressure coincide. For further details the reader is referred to the original literature [1] or one of the review articles available [2]- [4].
The HWEs are formulated as a three-dimensional system. They include not only the dynamics of the fluctuations but also the dynamical variation of the density profile and the electrostatic potential profile. In this context, we distinguish between the background density, the variations of the profiles and finally the fluctuations of density and potential. The background density is inherently assumed to coincide with the equilibrium pressure profile, denoted by n 0 , something that is important to keep in mind when considering more advanced models, that include MHD equilibria with magnetic shear and curvature. Variations in the profiles we define as the flux surface average denoted by n and ϕ respectively. The flux surface average reduces in the slab geometry to an average over the zand y-coordinates in three dimensions and over the y-coordinate in two dimensions. The fluctuations are then defined to have zero mean bỹ n = n − n . Investigations of the three-dimensional system itself reveal interesting results [5], even in the absence of magnetic shear and curvature. Note that the HWEs are a fundamental building block for more complex models, such as found in [6].
One important aspect of turbulence in magnetized plasmas is that the convection, i.e. the dominant nonlinearity, is restricted to the plane perpendicular to the magnetic field. This is the reason for the turbulence to be regarded as mainly two dimensional in the sense that it shares most of its features with classical two-dimensional fluid turbulence. The dimension parallel to the magnetic field, however, allows for a delayed electron response via the resistivity ν. The resistivity and the parallel mode number thus determines the deviation from the Boltzmann distribution-ϕ ≈ n in normalized units. Without such a deviation the fluctuations have no access 28.3 to the free energy in the background density gradient, which they transform into fluctuations of the E × B velocity field and the density. They in turn enhance transport and ultimately flatten the density gradient. As the turbulence is two dimensional-within limitations only the sources/sinks are three dimensional-the use of a further reduced two-dimensional system is justified. The parallel dynamics is absorbed into a parameter, using a principal parallel scale-length δ = L 2 ν: This reduced two-dimensional model does not in principle allow us to track profile modifications as the parallel response of the electrons is fixed, and can be expressed in the form, motivating the non-adiabaticity parameter for δ, as for δ = 0 an adiabatic response is obtained, In the three-dimensional model, however, the parallel coupling between the density and vorticity equation is absent for modes which have ∇ = 0, breaking the relationship given by equation (5). Therefore variations of the density profile or developing flows in the two-dimensional models should be interpreted with much care. One frequently applied method to retain profile modifications in a two-dimensional model is to assume a separate evolution equation for the variational quantities [7], ruling out coupling via parallel electron motion. Averaging the HWEs results in the following equations for the profiles: demonstrating the relevance of the polarization drift nonlinearity {ω, ϕ} for flow generation and of the E × B nonlinearity {n, ϕ,} for density profile changes. However, the two-dimensional version of the HWEs already contains all basic aspects of nonlinear turbulent interaction. They are therefore suited to investigate initial flow generation, transport statistics and nonlinear saturation mechanisms, while we are unable to make decisive statements about for example the transport level or the final state of the generated flows.
Using the two-dimensional version of the HWEs we briefly discuss the dispersion relation, energy balance and overall route to saturated turbulence. For plane waves we obtain the dispersion relation: For δ → 0, that is an adiabatic electron response, we find the real part of the drift-wave dispersion ω = ky 1+k 2 . With this we obtain in the limit k → 0 the growth rate of the instability as γ = Im(ω) = δ kyk 2 1+k 2 . We see from this that the instability arises due to the finite resistivity the electrons experience in their parallel motion. Multiplying the vorticity equation of the HWEs by ϕ and the density equation by n, integrating over the whole space and adding the two contributions we obtain for the energy-like quantity E ∂ ∂t E = 1 2 This shows that the term that contains the resistivity really acts as a sink of energy, while the energy in the fluctuations is fed by the E × B density flux Γ = − n∂ y ϕ dV . Note that this density flux 28.4 is zero for an adiabatic electron response (δ = 0) as n = ϕ holds, consistent with a vanishing linear growth rate. The fluctuation energy is thus, as expected, taken out of the background density gradient, while both the nonlinearities {ϕ, ω} (originating from the divergence of the polarization drift) and {ϕ, n} conserve energy and only redistribute fluctuation energy between different scales. Analogously the evolution of the so-called generalized enstrophy W is obtained: The enstrophy is also conserved by the action of the nonlinearities. If the electrons do achieve a Boltzmann response to a potential perturbation, for example n = ϕ, for δ → 0, the flux Γ vanishes and the turbulence conserves both energy and generalized enstrophy. For regions of the wavenumber spectrum where instability is absent and there is no in-or out-flow of energy or enstrophy respectively-the so-called inertial ranges-a dual cascade of energy and enstrophy is found [8]. The energy will be transported under the influence of the nonlinearity to ever larger scales, while enstrophy is stretched out to ever smaller scales until it enters the region of wavenumber space where dissipation sets in (which is added to the numerical simulations in form of an artificial viscosity). This is the main difference between two-and three-dimensional turbulence, the latter having only energy as conserved quantity in the absence of dissipation. The general behaviour of the turbulence is then as follows: starting from noise the linear instability develops and waves form and grow, with a wavenumber of about unity in normalized quantities. The energy E grows until it reaches a level where the nonlinearities become important. The waves then break up into vortical structures undergoing merging processes until a saturation wavelength is reached. The turbulence then saturates with a finite influx of energy, for example a finite particle flux that is balanced by the viscosity taking energy out at very small scales. This overall behaviour is exemplified in figure 1 and might be compared with the dynamics developing in a simulation of a more complete three-dimensional model [6] including magnetic shear, curvature, coupling to sound waves and electromagnetic effects shown as a comparison in figure 2, taken from a simulation by the Risø TYR code [9].
Further simplifications consider the dynamics in the absence of instability, for example for n ≈ ϕ. This leaves us with the raw competition between nonlinearity and drift-wave motion given by the Hasegawa-Mima equation (HME) [10]. In the adiabatic limit we find δ → 0 from equations (1) and (2): giving us the possibility to investigate the action of the nonlinearity in the absence of instabilities. It is interesting to note that the HME is formally identical to the Charney equation [11] describing Rossby waves in geophysical fluid dynamics (see e.g. [12]). This allows us to reuse in the plasma context the extensive research performed on the generation of large-scale geophysical flows (see e.g. [13]). It should be noted, however, that for geophysical flows the Charney equation is also suitable to describe all later stages of flow formation. The above-mentioned problem, considering parallel electron response for profile variations in the two-dimensional version of the plasma equations, is clearly absent from the geophysical context.

Set-up of flows
The self-organization of flows from turbulence has, in the plasma context, gained much attention over the last couple of years. It has been widely accepted that not only is the turbulence responsible for the observed anomalous fluxes, but poloidal flows generated via the turbulence are accountable for some of the transport barriers or even the L-H transition observed in modern plasma fusion devices. While poloidal or zonal flows have k x finite and k y → 0, radial flows have k x → 0 and k y . The latter are also called streamers and made responsible for the very high levels of transport observed in some situations [14]. Here we shall, however, concentrate on zonal flow generation and in addition consider the dissipation-free situation only. Averaging the vorticity component of the HWEs over the poloidal and toroidal co-ordinates (respectively y and z), we find for the evolution of an averaged poloidal velocity variation V (x) This transport can be re-written as demonstrating that the divergence of the momentum flux, as given by the divergence of the Reynolds stress uv , is equivalent to the radial vorticity flux by waves or vortices. This result was found in the geophysical context by Taylor [15] and even earlier by Reynolds [16,17].
Recently the influence of the Reynolds stress was investigated numerically with an emphasis on experimental measurement methods [18]. Thus, if we interpret flow generation as a transport process, we readily anticipate that it needs an anisotropic driving. In the plasma context this anisotropy is provided by the background density gradient and in the geophysical fluid context by a gradient in the Coriolis force, the so-called β effect. Here it is appropriate to note that weak turbulence, for example resonant three-wave interaction, is unable to put energy into modes with k y = 0, as the interaction sums to zero by means of symmetry. As three-wave interaction is not an effective ingredient in flow generation, we have to look at the part of wavenumber space that allows for strongly nonlinear interactions. This will be the region in wavenumber space where the nonlinear terms dominate over the linear wavelike motion. Only in that domain can higher-order coupling, that is higher than three-wave interaction, between modes be effective. The inverse cascade will favour structure formation with low k y and finite k x , for example zonally extended structures.
The crossover from the linear to the nonlinear regime is estimated by comparing the relevant characteristic frequencies obtained from the flow generating nonlinearity ω turb = k 4 φ/(1 + k 2 ) and the linear term allowing wavemotion ω wave = k y /(1 + k 2 ). Using an average velocity U = kφ to describe the energy contained in the turbulence and approximating k y ≈ k we easily equate the isotropic, in the fluid context so-called Rhines length [19,20] as Note that within our normalization the inverse gradient scale length κ n has been absorbed into the definition of U . Reconsidering the anisotropy of the linear wave term we correspondingly arrive at One should note that the scale at which this transition occurs depends strongly on the energy contained in the fluctuations. Thus, in a driven system with growing energy we find that l Rhines is growing with time, while for decaying turbulence, smaller and smaller scales lose their capability of nonlinear interaction as the energy content decays. In figure 3 the dumbbell-shaped line marks the crossover from dominating nonlinear behaviour outside, while for smaller k, for example larger scales, linear, wavelike motion dominates [21]. For k-values inside the dumbbell linear decorrelation makes an effective nonlinear energy transfer into these modes unlikely; structures with a finite k x and k y → 0 are preferred by the inverse cascade that describes the general direction of energy transfer in two-dimensional turbulence. This corresponds to the nonlinear generation of zonal flows, while the generation of streamer-like structures, for example finite k y and k x → 0, is suppressed. This behaviour is seen in the decaying turbulence simulation depicted in figure 4. Starting with a ring in wavenumber space, the energy cascades inwards and then settles in a structure corresponding to that shown in figure 3. Naturally, this happens just at the crossover from the nonlinear to the linear type of dynamics. Initially the formation of round (as the nonlinearity is isotropic) localized vortical structures is seen. As the turbulence cascades towards the Rhines length l Rhines these structures become less and less pronounced and finally no longer play a role in the dynamics. However, in a driven situation these vortical structures have more importance and in the next section we shall focus on their dynamics and investigate whether and how they might contribute to the generation of flows.

Nonlinear structures and flow generation
As was shown in the discussion of the Reynolds stress and vorticity transport the polarization drift nonlinearity {ω, ϕ} dominates the turbulence behaviour with respect to flow generation. This nonlinearity is isotropic and homogeneous; its action alone favours the formation of circular, nonlinear vortical structures. While these structures are solutions to the HME in the limit k 1, e.g. the two-dimensional Euler equation, they are not solutions to the HME itself and will radiate waves in that situation. However, vortical structures are a building block of the turbulence and therefore their behaviour can shed light on flow generation from a different kind of aspect. As a model nonlinear structure we consider a Gaussian pulse with an inverse length k ≈ 1 and an amplitude ϕ max = 0.1. Thus, we are in the linear regime where the nonlinear time t nl is shorter than the linear dispersion time of a wave package: or in other words the energy of the pulse results in a Rhines length longer than the scale size of the pulse, putting the pulse in the domain dominated by wavemotion. Figures 5(a) and (b) show that the pulse consequently disperses into waves. At t = 120 no sign of the initial structure is found in the contour plot of n. The corresponding density variation n(x) (figure 5(c)) stays constant in accordance with a vanishing radial E × B particle flux. The observed very small initial change in the profile, however, already reflects the influence of the nonlinearity. For a Gaussian pulse with amplitude ϕ max = 20, resulting in t nl /t disp = 20, the dynamics is dominated by the nonlinearity {ϕ, ∇ 2 ⊥ ϕ}. The corresponding Rhines length is-for this amount of energy-thus shorter than the pulse scale. Figures 6(a) and (b) show that the pulse survives for long times as a vortex structure. Note that a positive density vortex corresponds to an anti-cyclone as negative vorticity w = ∇ 2 φ is associated with it [22]. The vortex moves inwards in the radial and upwards in the poloidal direction. This is a well known behaviour for isolated anti-cyclonic vortices on the β-plane [23,24]. To keep their potential vorticity W = φ − ∇ 2 φ − x constant, they move towards lower x while φ − ∇ 2 φ decreases. While the structure propagates and radiates waves, the density profile changes significantly (figure 5(c)). Part of this change is directly connected to the radial movement of the coherent structure, whose outline is shaded in figure 5(c). Another part of the flow is due to large-amplitude drift-waves radiated from the structure. The change of the density profile demonstrates a non-vanishing flux, which originates, as will be seen below, from the polarization drift. Note, however, that this anti-cyclonic structure transports density into regions of higher background density. Symmetry of the HME shows that a cyclone will just move oppositely.
Let us now consider the transport associated with the polarization drift. Starting from the continuity equation we have v pol = −d t ∇ϕ. To lowest order in the drift-wave scaling holds. Thus the time evolution of the density profile is given by For the y-averaged flux we derive in the long-wavelength limit with d t · = ∂ t · +{ϕ, ·} In our normalized units n 0 ≈ L n /ρ s 1 holds, and consequently we neglected the radial dependence of the background density. Note that the first term, containing the background density, as well as the second one, vanishes when averaged radially. We thus only keep the last term on the rhs of equation (20). It can be split into two components using to lowest order the time evolution given by the HME: u and v being the E × B velocity components in x respectively y. Taking additionally the integral over x of (21), it takes the form of the Reynolds stress. Thus, here the buildup of a radial electric field, or zonal flow, is a consequence of the non-ambipolar transport connected to the polarization drift. 28.12

Saturation
Let us now turn to the question of saturation of the instabilities. This is important as for a situation with instability we have ever more energy cascading to lower wavenumbers. This results in l Rhines → ∞. Nothing sets a typical size for the turbulence observed; moreover, we expect a catastrophic growth of the fluctuations resulting in a violation of the drift-wave scaling order used to derive the model. If we look at the corresponding geophysical situation we note that saturation of Rossby waves takes place by two means: either via a relaxation of the driving background gradient or due to the presence of Eckman damping [25], that is due to the friction with a bottom layer. These mechanisms are absent from the plasma model, where we assume that the background gradient is kept up by heating and a large-scale Ekmantype damping is absent. A sort of quasilinear reduction of the background gradient might however take place for flows bounded in the direction of the gradient, so that the flux is not able to leave the system [26]. We thus have to look for another mechanism that can lead to saturation in drift-wave turbulence and sets the characteristic spatial scales observed. We therefore describe the turbulent energy exchange by statistical means via the eddy-damped quasi-normal Markovian (EDQNM) method [27]. We start from the reduced HWEs close to the adiabatic regime (δ < 1) [28,29]: which allows us to investigate the influence of the so-called E × B nonlinearity {n, φ} in a simple way. Note that this nonlinearity vanished in the adiabatic n = φ case, but is relevant whenever there is an instability due to a deviation from the Boltzmann response of the electrons. After Fourier transform we write with the nonlinear term and coupling coefficients and Within the EDQNM theory we are mainly interested in expressing the action of the nonlinearities as a sort of correction to the linear growth rates and frequencies. Details of the calculation are given in the appendix. We identify three contributions to the nonlinear dispersion relation: the polarization drift nonlinearity which both add to a growth rate, but not to a frequency shift; finally we have a cross term which leads to a nonlinear frequency shift, To qualitatively evaluate the nonlinear growth rate (27) and (28) we use the linear times for τ ww and approximate |φ k | 2 by the two-dimensional Kolmogorov spectrum [30], that we obtain in the energy inertial range as and for the enstrophy inertial range. Figure 7 shows the nonlinear growth rates γ P ol and γ E×B in arbitrary units. Both nonlinearities conserve energy and just transport energy between scales. The polarization drift nonlinearity acts in the sense of an inverse cascade, leading to a growth rate of modes with small k and taking out energy from the high-k modes. The E × B just works the other way round, but is usually much weaker and only for k 1 can it dominate the polarization drift nonlinearity. Figure 8 shows γ ef f = γ P ol + γ E×B (δ) for different values of δ. Note that k krit , at which the E × B nonlinearity stops the inverse energy cascade, is found at k krit ≈ 1 for δ of order one. Also shown in figure 8 for two values of δ is the value of k av at which the inverse cascade is stopped, as obtained from a numerical solution of the HWEs. Note that for small values of δ close to adiabaticity k krit decreases dramatically, so that only for very large amplitudes or large scales can saturation be achieved.
The nonlinear frequency shift ω cross is depicted in figure 9. We note that in the HWEs this frequency shift influences the linear growth rate even in the lowest order Iterating with the linear growth rate we obtain Thus the nonlinear frequency shift can especially for k < 1 change the instability into a damping, while for other k-values it might enhance the linear instability. For k 2 k y > 1 the growth rate is enhanced for rather small k x . This is an indication for the possible occurrence of structures strongly elongated in x-streamers-as the numerical simulations show for certain parameter regimes. These would be a reflection of the nonlinear dynamics on the growth rate, rather than self-organized nonlinear structures of the turbulence.

Conclusion
Using strongly simplified, but nonlinear models of drift-wave turbulence, we have demonstrated mechanisms of flow generation and saturation. It was demonstrated that the nonlinear behaviour is only weakly connected to features deferred from linear wave properties, for example the typical length scales of the fluctuations. Much can be learned about the turbulent system by using rather simple arguments, even before detailed analysis starts. The basic behaviour observed within the used toy models can in great part be recovered as being embedded in more complex models. It proves useful to classify features of the more complex models. Last but not least we can infer methods and ideas to influence and control turbulence in experiment [31], by testing and learning from simplified models.
Discussions on the topic of vortex motion with J Juul Rasmussen and A H Nielsen are gratefully acknowledged.

Appendix. EDQNM theory
To evaluate the nonlinear dispersion relation in the EDQNM approximation we write in a shorthand notation [32] L(1 2)φ(2) for the evolution equation. Here and are defined utilizing the summation convention.
(A. 22) ω has to be determined iteratively. To lowest order one recovers the real part of the drift-wave dispersion relation ω k = ky 1+k 2 . Note that ω is on the fast timescale τ , while growth rates belong to the slow timescale T .