Energetics of the interaction between electromagnetic ExB turbulence and zonal flows

A detailed analysis of the energy transfer system between ExB turbulence and zonal flows is given. Zonal flows, driven by the ExB Reynolds stress of the turbulence, are coupled to pressure disturbances with sinusoidal poloidal structure in toroidal geometry through the geodesic curvature. These pressure ‘sidebands’ are nonlinearly coupled not only back to the turbulence, but also to the global Alfvén oscillation whose rest state is the Pfirsch–Schlüter current in balance with the pressure gradient. The result is a statistical equilibration between turbulence, zonal flows and sidebands, and additionally the various poloidally asymmetric parallel dynamical subsystems. Computations in three-dimensional flux surface geometry show this geodesic transfer effect to be the principal mechanism which limits the growth of zonal flows in tokamak edge turbulence in its usual parameter regime, by means of both control tests and statistical analysis. As the transition to the magnetohydrodynamic (MHD) ballooning regime is reached, the Maxwell stress takes over as the main drive, forcing the Reynolds stress to become a sink.

The interaction between turbulence at medium to small scales and flows at larger scales is a topic of currently general interest to research in the dynamics of magnetized plasmas [1]- [9]. The usual situation of a gradient-driven turbulence in these plasmas is a competition between a dynamically incompressible ExB flow perpendicular to the magnetic field and a self-consistent response in the parallel direction [10]. This basic picture remains even in the core, where parallel dynamics render the electrons adiabatic [11], but the dynamics of the ion temperature continue to give rise to toroidal ion temperature gradient (ITG)-driven turbulence [12].
The large scale ExB flows with which the turbulent eddies interact are the zonal flows: the result of the flux surface ('zonal') average of the electrostatic potential. The interaction is largely through the Reynolds stress, with the energy transfer direction preferentially from the turbulence to the flows [2]. This gives the physical manifestation of the basic argument that an externally imposed ExB shear should lead to faster decorrelation of smaller-scale turbulent eddies and hence suppression of turbulence [1]. The correlation length, however, has no energetic relevance as the main point is the tendency of the eddies to pass energy to the background flow as they are tilted by the flow.
The Reynolds stress process is a variant of the more general inverse energy cascade tendency from smaller to larger scales in two-dimensional (2D), incompressible turbulence [13]. In detail, it proceeds through a statistical instability through which the zonal vorticity (ExB shear) due to the turbulence starts as random, but then self-amplifies, as localized regions of a given sign of vorticity intensify a Reynolds stress of the same sign in the same place and time. This instability, demonstrated in gyrokinetic computations [14], has been treated as a modulational instability [15], that is, a direct transfer to the flow, seeded by the flow itself, rather than a local cascade. In the case of imposed ExB shear, the turbulence suppression has been shown in drift-wave computations to proceed energetically [3]. Imposed ExB shear which is time-dependent was also studied as a model for the flows which occur self-consistently; the efficiency of turbulence suppression was found to degrade if the frequency of the flow shear was comparable to or greater than the frequencies of the turbulence [17]. Large scale gyrokinetic computations have shown the emergence of self-generated zonal flow shear layers to be very important in limiting the radial scale of the turbulence and consequently the resulting transport [5]. Their importance in fusion research lies in the fact that zonal ExB flow (electric field) shear is believed to underlie the transition and maintenance of the H-mode operation of tokamak confinement [18].
In toroidal geometry, the ExB flow is compressible, which is a weak effect on either magnetohydrodynamics (MHD) [19] or the turbulence [20] but a strong effect on the flows [21,22]. Due to the variation of a toroidal magnetic field with the major radius (B ∼ R −1 ), a zonal flow is faster on the outboard side of the torus and slower on the inside. The resulting compression gives a conservative energy transfer with the 'sideband' modes of the pressure (varying with the sine of the poloidal angle). The most basic manifestation of this coupling between flows and pressure is the geodesic acoustic oscillation [23], whose natural frequency is √ 2c s /R times numerical factors in general, where c s = √ T e /M i is the isothermal electron sound speed. It was recognized that this oscillation would be quite slow compared to the turbulence, which is broadband up to a maximum of c s /L ⊥ , given that the profile scale length L ⊥ is much smaller than R [8,9]. Under these conditions, the dynamics of the interaction between turbulence and zonal flows should change character.
Recently, it has been shown that the tendency of fully resolved drift-wave turbulence [24] to spontaneously generate zonal flows is very robust in slab geometry if the flow evolution is not somehow suppressed [25], much the same as in the core ITG case [4,5]. In highresolution computations in toroidal geometry, however, self-generation of a strong ExB shear layer was not found [26]. The toroidicity was found by a control test against sheared slab geometry, both in three dimensions, to be the agent preventing this scenario [27]. The cause was found to be the 'geodesic transfer effect', which is the combination of the conservative transfer of free energy from zonal flows to pressure sidebands, and the nonlinear damping of those sidebands [9]. It is important to note in the above context that these effects of geodesic curvature have not been part of many otherwise sophisticated computations to date, since the ExB flow is often taken as incompressible. These toroidal effects were present before [28], but their significance was only recognised later [8,9]. In the most sophisticated core turbulence computations to date, self-generated zonal flows are found to modulate the turbulence rather than suppressing it [29], and the geodesic oscillation per se was not found to dominate the zonal flow frequency spectrum [30]. Both results are consistent with the findings for edge turbulence in [9,26].
Herein, the energetics of this interaction is studied in detail. The main points of this study are (1) dynamics in the zonal flow/sideband system, (2) the energy theorem for the zonal flow/sideband system and (3) energetic interaction with the turbulence, focusing upon (4) the role of the geodesic transfer effect as the principal mechanism constraining zonal flow selfgeneration at edge turbulence parameters. The following sections present the basic model used in the analysis and computations, the general analysis of the interaction between zonal flows and sideband modes and the computational results. Both two 2D and 3D results are shown, making the contrast between an unopposed Reynolds stress flow drive in the 2D case and the statistical equilibration between zonal flow and pressure sideband free energy in the toroidal 3D case. The statistics of the zonal flow energetics is then presented.
The geodesic curvature transfer mechanism is found to be the principal agent limiting the self-generation of zonal flows by turbulence in any situation, where (1) the Landau damping is relatively weak, (2) the electron density ExB nonlinear advection is relatively strong and (3) plasma pressure is insufficient to reach the nonlinear ideal MHD regime; in other words, in any relevant regime of tokamak edge turbulence. The zonal flow energetics in this situation is characterized by a balance between ExB Reynolds stress acting as a drive and the geodesic transfer mechanism acting as a sink. Most of the zonal flow energy passes through pressure sidebands and then into parallel currents through the adiabatic response in the electron parallel dynamics. This is the global Alfvén oscillation whose rest state in a relaxation scenario is the Pfirsch-Schlüter current in balance with the pressure gradient. Dissipation of this oscillation by resistivity provides much of the ultimate sink.
As the ideal MHD regime is entered, more effects are introduced: the Maxwell stress, tending to cancel the Reynolds stress as more of the overall energy resides in shear Alfvén component of the system [42], and the increased free energy in the global Alfvén oscillation, tending to disrupt all the balances since the throughput of energy becomes larger than that in the Reynolds stress. This is experimentally unrealistic, but of physical interest nevertheless.
As the core regime is approached (passing electrons becoming adiabatic), the basic geodesic frequency becomes faster relative to the turbulence since (1) the turbulence is a weak turbulence controlled by a narrow set of modes at relatively long wavelength and (2) the scale ratio R/L ⊥ is not very large. Moreover, the dominant sideband-damping mechanism becomes ion Landau damping [31,32]. Initial computations with an electromagnetic gyrofluid model [26], modified to properly conserve free energy at all levels of finite gyroradius in the non-dissipative processes, indicate that the geodesic transfer effect continues to function even in this regime [33], although detailed analysis at the level presented herein is left for future work. Recent experimental work as well as two very recent theoretical efforts are briefly commented upon in the conclusions (section 7).

The drift Alfvén model DALF3
This study uses the same model as in [9], the isothermal version of the drift Alfvén model described in detail in [28]. Standard drift ordering is taken [34]. The geometry is that of a globally consistent flux tube, in which although fluxtube ordering is used [35], the parallel boundary condition is the same as that on the entire flux surface for every Fourier component or toroidal mode [36]. The shifted metric model is used to properly represent the slab-mode structure if this occurs [37]. Computations are set up as in [38]. The equations are normalized with standard gyro-Bohm forms: time to L ⊥ /c s , perpendicular spatial dimensions to ρ s and the parallel dimension to qR, where c 2 s = T e /M i and ρ 2 s = c 2 M i T e /e 2 B 2 following electron mobility and ion inertia, and L ⊥ is the perpendicular profile scale length and 2πqR is the parallel field line connection length. Relative amplitudes are scaled with an additional factor of δ = ρ s /L ⊥ , so that an amplitude ofφ = 1 in normalized units signifies eφ/T e = δ in physical units. The equations for the four dependent variables (vorticity˜ , electron pressurep e , parallel currentJ and parallel ion velocityũ ) are with consistency relations (Ampere's law and polarization) for the field potentials (electrostatic potentialφ and magnetic potentialÃ ), whereW =φ + τ ipe is the total ion flow stream function. Action of v E · ∇ or ∇ upon p e represents the gradient forcing terms and d/dt = ∂/∂t + v E · ∇ is the ExB advective derivative. The second term on the left-hand side of equation (1) results from the gyroviscous cancellation upon advection [39,40], which is also implicit in the appearance of d/dt in equation (4). The quantities τ i ,β,μ, C,ˆ and µ are constant parameters, discussed below. The geometry is described in terms of a field aligned flux tube geometry under the shifted metric treatment, with coordinates {x, y, s} representing the radial, perpendicular drift and parallel directions, noting that s is a projection of the parallel direction onto the poloidal angle [37]. The simple sheared flux surface model is used, neglecting all finite aspect ratio effects except for the existence of the curvature operator K. Hence, the divergence of the ExB velocity is represented by K(φ), while the ExB velocity appearing in the advection terms is divergence free. The normalized magnetic field strength is The ExB and parallel derivatives are for the ExB velocity and for the magnetic field. Action by v E · ∇ and ∇ in the {x, y} drift plane represent the ExB and magnetic flutter nonlinearities, respectively. The background gradient is in the negative x-direction, The manifestation of magnetic shear is a series of shifts in the y-coordinate while taking derivatives in the s-coordinate, where h s is the (equidistant) node spacing in s, and the shifts are given by whereŝ is the standard parameter representing the magnetic shear. This follows from the fact that the coordinate system is defined separately at each location in s, such that perpendicular operators are always evaluated with an orthogonal metric. In terms of a concentric circular model for a tokamak with minor radius and poloidal and toroidal angles {r, θ, ζ}, the coordinates are defined at a reference flux surface r = a as separately at each poloidal location θ = θ k , with shearŝ = d log q/d log r and with x and y k normalized to ρ s and s defined on [ − π, π]. The index k on the y-coordinate is understood throughout, though it is suppressed for clarity. The perpendicular Laplacian operator is given at s = θ k by The curvature operator, which will play a prominent role in this study, is given at s = θ k by defining {K x , K y } as a vector, where ω B is a parameter nominally equal to 2L ⊥ /R. The breakdown of this into the geodesic curvature and interchange forcing effects is different from a simple split into K x and K y , however, as discussed in the next section. The parameters controlling the parallel electron dynamics arê where V e is the electron thermal velocity given by V 2 e = T e /m e , ν e is the conventional collision frequency and the 0.51 comes from the parallel resistivity [41]. The sound waves are controlled byˆ = (c s /L ⊥ ) 2 (qR/c s ) 2 , just the square of the parallel/perpendicular scale ratio, and damped by parallel ion viscosity, µ , which in this model functions as a substitute for ion Landau damping.
Additionally, τ i = T i /T e is the background ion/electron temperature ratio,ŝ gives the magnetic shear and the curvature operator is scaled with ω B , which can be set independently (slab geometry is ω B = 0). The standard case is a set of nominal parameters corresponding to a typical tokamak edge for a neutral, single component deuterium plasma, witĥ roughly reflecting physical parameters: Cold ion cases are recovered by setting τ i = 0. We leave µ = 0 for most cases, but then use it to check for the role of parallel sound wave damping on the overall statistical saturation of the zonal flow/sideband system.

Zonal flows and geodesic coupling
In a simple 2D system, there are only the two coordinates {x, y} and hence Fourier components, or 'modes,' with wavenumbers {k x , k y }. The directions are distinguished by the gradient and the interchange forcing in the x-direction, both working through a finite ∂/∂y. The most important linear instabilities in either slab or interchange models with ω B = 0 or ω B > 0, respectively, will usually have k x k y , ideally forming radial plumes with k x = 0. Turbulence, by contrast, acts through the isotropic nonlinearities and hence tends to produce eddies with k x ∼ k y with both finite and of either sign. The zonal flow modes are the ones with k x = 0, with either sign, but with k y = 0. The Reynolds stress is given by (20) and is therefore mostly a result of the eddies. However, if the 2D spectrum in wavenumbers is statistically isotropic, the average of R E will vanish although there will generally be strong positive and negative samples. Dividingφ into the k y = 0 and k y = 0 parts, the energy transfer between them can be shown to be the zonal vorticity (two x-derivatives of the zonalφ) times the zonal Reynolds stress [2]. Each eddy component {k x , k y } contributes to the transfer through its contribution to the k y = 0 component of R E and the zonal vorticity. Hence, if a zonal flow is allowed to simply tilt the vortices according to the sense in which it is sheared, the resulting distribution of R E indicates that energy is transferred from eddies to the zonal flow, amplifying the latter. From a random initial phase, the part of the distribution in which zonal flows and the zonal Reynolds stress are aligned will self-amplify and the turbulence will tend to suppress itself by transferring its energy to the zonal flows. This result is well known and will be illustrated in section 6.1, below.
In three dimensions, the zonal flow mode is the flux surface average, with both k y and k wavenumbers zero. It is important to note that in a tokamak magnetic field, parallel is not toroidal and perpendicular is not poloidal. In terms of the poloidal and toroidal mode numbers {m, n}, the (unnormalized) drift direction and parallel wavenumbers in a simple torus are given by k y = nq/r and k = (m − nq)/qR. Indeed, the flux tube coordinate system represents the poloidal angle with the parallel coordinate. Nevertheless, once the y-direction is averaged, isolating the k y = 0 component, the field aligning and even the shifts in equations (12) and (13) do not enter, and an average over first y and then s is the same as a flux surface average. The new modes in the 3D model are those with k y = 0 but still k is finite, that is, n = 0 but m = 0. These are axisymmetric structures called 'sidebands', coupled to the zonal modes by the curvature operator K, which has the sinusoidal structure in equation (16).
The action of K upon the dynamical variables incorporates both geodesic and interchange curvature. Often these are referred to as simply the parts which act through finite ∂/∂x or ∂/∂y; in other words, K x or K y . The geodesic coupling effect, however, refers to the action of K solely upon the axisymmetric part (∂/∂y hence k y = 0). The k y = 0 modes are those which cannot carry transport, since for them the radial velocityṽ x E vanishes. Both the zonal mode and the sidebands are included in this set; they must be considered together because of the fact that the (sin s) structure of K x couples them, so that they conserve energy as a unit [43]. The split between geodesic and interchange curvature components for the turbulence (k y = 0) depends on the coordinate system (specifically, whether the shifted metric treatment is used), and so the entirety of the action of K upon the part with k y = 0 is referred to as the 'interchange forcing' effect. These modes collectively are the eddies which can carry transport. When geodesic coupling produces an energy transfer of a particular sign, it is then referred to as the geodesic transfer effect.
In a 3D sheared slab model, K vanishes but the sidebands are still coupled to the rest of the system by the nonlinearities, such that when in a local 2D drift plane the energy is transferred to the k y = 0 component, and the entire spectrum in k is excited; not all of the transfer goes directly into the zonal mode. Those portions with finite k feel the adiabatic response ∇ (p e −φ), so that much of their energy is dissipated. In a 3D model, this results in a somewhat weaker zonal flow drive than in 2D, such that it becomes very important to keep enough computational resolution in order to find the level of suppression of the turbulence-generally, the fact that the vorticity acts most strongly near ρ s scales becomes the most important reason that ρ s must be resolved using a proper computation [38]. With the Reynolds (and Maxwell) stress unopposed, saturation of the flow/turbulence system can only occur through nonlinear self-balance. The usual result for either 2D or 3D slab cases is a small Maxwell stress and a balance in the Reynolds stress taking opposite signs in the long and short wavelength regions in the spectrum, the same as in the case of externally applied flow shear [3].
In a toroidal 3D model, however, this changes. Although the low-frequency motion remains dynamically incompressible, quasi-static compression results from the ExB motion in the inhomogeneous magnetic field due mostly to the 1/R dependence of B. Energy which does transfer nonlinearly into the zonal mode is now transferred through the geodesic curvature to the sidebands; the zonal flow is the zonal mode of ∂φ/∂x and therefore is subject only to action by K x (here, it is useful to note that ∂/∂x commutes with the flux surface average as well as with poloidal dependences of quantities such as K x ). Physically, the zonal flow is a perpendicular rotation of the entire flux surface and a finite zonal vorticity indicates that this rotation is sheared. Due to the fact that the magnetic field is weaker on the outboard side of the torus, the flow is faster there and so there is a divergence on the top and bottom of the torus. This finite divergence of the ExB velocity excites a sinusoidal disturbance inp e by compressional work; hence, there is a coupling to the pressure which is ±1 in the poloidal mode number. The coupling is conservative, however, so that the first result is simply an oscillation between the zonalφ and the k qR = ±1 sidebands ofp e , in other words, the geodesic acoustic oscillation [23].
The significance this has for the interaction between the turbulence and zonal flows is that the oscillation as a conservative system forced by turbulence will tend towards a statistical equipartition of free energy between the two components of the oscillation. That part of the free energy residing in the sidebands is subject to dissipative Alfvén dynamics through the adiabatic response (∇ J and ∇ p e ), and 'diffusive mixing', which is a nonlinear interaction with the turbulence itself, tending to transfer energy from larger to smaller scales through the unsteady ExB advection v E · ∇p e . For the edge situation, L ⊥ qR, i.e. the coupling to MHD sound wave activity is relatively weak. A balance will result in which every piece of the dynamical system is statistically stationary. For the turbulence, it is struck between gradient drive and a combination of resistive Alfvénic dissipation and nonlinear transfer to subgrid scales as a sink. For the zonal flows, it is between the Reynolds stress as a drive and the geodesic transfer as a sink. For the pressure sidebands, it is between the geodesic transfer as a drive and a combination of the resistive Alfvén dynamics and diffusive mixing as a sink.
The only question is: how strong does the zonal flow amplitude and energy become before this balance is reached? Indeed, one could skeptically ask whether a balance is reached at all with a finite Reynolds stress drive. In other words, should the zonal flow statistical instability saturate with the large flow levels of the slab model, with low-k ⊥ Reynolds stress transferring from flows to turbulence and high-k ⊥ Reynolds stress transferring the other way, with direct cascade throughp e and inverse cascade throughφ making up the connections? Or should this geodesic transfer mechanism become so effective as to limit the zonal vorticity to levels too small to suppress the turbulence at all? We will find that the answer is closer to the second possibility, with a drastic difference between the slab and toroidal cases, and even between slab and toroidal cases modified such that the geodesic coupling effect is inserted into or removed from the model. To demonstrate this we need the computations, but before presenting them we will introduce the theory of all these coupling mechanisms in general, since there are many pathways among all the various degrees of freedom: one degree of freedom per mode per dependent variable, and more than one transfer pathway per degree of freedom.

Zonal flow and sideband dynamics
The theory of the zonal flow/sideband system follows from the dynamical equations (1-4). We will consider the degrees of freedom with k y = 0, both the zonal mode and the sidebands. There are six relevant degrees of freedom: the zonal averages of the two state variables and the sideband modes of all four variables. We define the electron parallel flow, the total perpendicular flow potential and vorticity, and the total pressure and electron force potential, and we note that in the free energy under drift ordering, the thermal energy is proportional to (1 + τ i ) times p 2 e . We consider the profile together with the disturbance and drop the 'tilde' symbol, so that This results in the Pfirsch-Schlüter equilibrium being treated self-consistently. The averaging conventions are angle brackets for the zonal average f = ds dy f (24) and braces for the radial domain average used in the energetics We note that · · · annihilates ∂/∂y and ∂/∂s and commutes with ∂/∂x. Particularly, only the geodesic piece survives the zonal average over a curvature term; moreover, only it survives the zonal average of any function of s times a curvature term.
The key assumption usually made in studies of rotation is that the sidebands are assumed to be in quasi-static equilibrium, and the zonal quantities are assumed to vary on the slow transport timescale, with the turbulence appearing as a set of forcing terms [31,32]. Herein, the computational results are used to diagnose what actually takes place, as the action of the turbulence is too complicated to be solved analytically. The phases of the four sideband modes are suggested by the action of K upon the two zonal modes. The sideband evolution equations are found by multiplying the state variable equations (1) and (2) by sin s and the flux variable equations (3) and (4) by cos s and then taking the zonal averages. The zonal evolution equations are found by taking the zonal averages of equations (1) and (2). In each case, it is important to keep the forcing terms provided by the turbulence. In the algebra, it is useful to recast the nonlinear terms as divergences, e.g. v E · ∇p → ∇ · (pv E ) or ∇ · (v E · ∇)∇ ⊥ W → (∇∇) : (v E ∇ ⊥ W). The components of the ExB velocity are given in equation (9), and the drift direction ion flow is given byũ y = ∂W/∂x. The components of the magnetic field are given in equation (10).
The sideband dynamical equations are If a quasi-static balance is assumed, then the first of these gives the Pfirsch-Schlüter current equilibrium and the second gives the associated electric field, assuming an MHD limit in which φ p e in the sidebands, i.e. the part with ∇ = 0. We will find, however, that the MHD limit is not a good assumption, and that the 'two fluid' terms-principally the adiabatic response terms, ∇ p e in Ohm's law and ∇ J in the pressure equation, but also the diamagnetic flow compression term, K(p e ) in the pressure equation-play a strong role even in the sideband balances.
The zonal evolution equations are If for a moment, we neglect the turbulence altogether and assume that the sidebands are in quasi-static balance, we would find that p sin s = −µ u cos s (35) for the sidebands, and hence for the zonal quantities. In other words, Pfirsch-Schlüter transport for the pressure and rotation damping is due to parallel viscosity. In the computations, we may therefore determine the role of dissipative flow damping effects by setting µ to finite values or to zero. In most cases we will have µ = 0, so that nonlinear processes which can disturb the sideband balances are necessary for the zonal vorticity (hence flow shear) to find statistical equilibrium in the presence of a robust Reynolds stress drive. A small residual u cos s always remains as a result of finite numerical dissipation, without which the computations cannot be managed. It is ensured, however, that this dissipation is <0.01 in normalized units, noting that L ⊥ /qR ∼ 0.01 is usual.

The geodesic acoustic oscillation
The geodesic acoustic oscillation is visible in equations (28) and (30), if we integrate one factor of ∂/∂x and drop the surface terms in equation (30) in the two partial time derivatives and the two curvature terms multiplied by ω B . Displaying just these terms, the dynamical system reads Starting with a zonal flow ũ y , the compression on top and bottom of the torus gives rise to the pressure sideband p sin s . The associated diamagnetic current itself has a finite quasi-static toroidal compression with a zonal component. This is balanced by a polarization current which keeps the total current divergence free. The MHD Lorentz force resulting from this polarization current acts to decelerate the zonal flow. It is a triviality that this restoration is due to the polarization current, since by definition the diamagnetic current balances the pressure gradient and therefore does not accelerate the plasma, and the parallel current does not contribute to the J × B Lorentz force. The natural frequency in normalized units is For edge turbulence this is rather slow compared to the native drift frequency c s /L ⊥ , since R/L ⊥ is of the order of 50 in the edge region. The initial reaction of the turbulence to a sudden change in parameters is typically faster than 0.01c s /L ⊥ (i.e. it is manifest within about 100 in normalized time units). This is robust enough for the turbulence to be sensitive to the geodesic acoustic oscillation, not only as a whole, but also during the peak and trough phases of a single oscillation. The turbulence is therefore able to interact nonlinearly with the oscillation. Ultimately, as we will find out, not only the the nonlinear terms in both equations (28) and (30) but also the two-fluid adiabatic response effects directly coupling p sin s to J cos s are very important as source and sink effects in the oscillation system, to the extent that the oscillation per se is not important to the turbulence and zonal flow interaction. The energy transfer channels, on the other hand, are.

Zonal flow and sideband energetics
Each degree of freedom in the dynamics has its own contribution to the overall free energy. The energetics for the turbulence in the drift Alfvén dynamical system is described in detail elsewhere [16,28,38]. In the present case, we have the two zonal modes, and p , evolving according to equations (30) and (31), and the four sideband modes, φ sin s and p sin s and J cos s and u cos s , evolving according to equations (26)- (29). The ExB energy is given by the domain average of −W × ∂ /∂t, the thermal free energy and sound wave energy are given similarly by (1 + τ i )p × ∂p/∂t andˆ u × ∂u /∂t, and the magnetic and electron parallel kinetic energies combine into the electron parallel energy, given by the domain average of J × (βA +μJ ). The various modes ({k y , k } pairs) contribute separately to these quadratic forms. Surface terms resulting from integration of ∂/∂x by parts are neglected.
The zonal flow energy evolves according to ∂ ∂t This shows two nonlinear drive/depletion terms, depending on their statistically average sign and a coupling term, which is also of indeterminate sign and could either act as a drive or a depletion. The first term is the Reynolds and Maxwell stress. For a shear Alfvén mode structure, these stresses cancel [42], although in general there is an imbalance even when the Maxwell stress is finite. The Reynolds stress mechanism is usually a drive, due to the statistical zonal flow instability. The drive is positive when the zonal vorticity is positively correlated with the zonal Reynolds stress ṽ x Eũ y . The last term in equation (40) quantifies the geodesic coupling. If the net drive due to the nonlinear terms is positive, energy can only leave the zonal flow by means of this effect. If the effect is to transfer energy out of the zonal flow in a net sense, it is referred to as geodesic transfer. The zonal flow energy can only be transferred to the pressure sideband.
The pressure sideband energy evolves according to Five terms are listed on the right-hand side. As in the case of equation (40), all have indeterminate sign. The first term gives the nonlinear coupling to the turbulence (k y = 0) modes, referred to as diffusive mixing if the transfer is preferentially out of the pressure sideband. The second term, arising from the toroidal compression of the ExB velocity, gives the geodesic coupling to the zonal flow, and since it appears in both equations (40) and (41) with opposite sign, we identify it as a conservative transfer effect between them. The third term gives the sound wave coupling to the parallel flow sideband. The last two terms give the adiabatic coupling to the parallel current sideband and the diamagnetic compressional coupling to the thermal background. The factor of cos 2s arises from sin 2 s and gives the coupling to the second sideband. These k qR = 2 modes are found in the computations to be quite negligible, however, and we will not consider them further. Energy leaving the pressure sideband can go basically anywhere else in the system. An important condition governing this is for the Pfirsch-Schlüter current balance given by equation (32) to be broken. If so, then a flow (vorticity) sideband is driven, as seen in equation (26). The flow energy sideband evolves according to The first term is the sideband of the Reynolds/Maxwell stress and the second and third terms give the MHD part of the global Alfvén oscillation, just the imbalance in the Pfirsch-Schlüter current correlated with the flow potential sideband W sin s . These give the Alfvénic coupling to the parallel current sideband and the diamagnetic compressional coupling to the thermal background. The Pfirsch-Schlüter balance can therefore only be broken in a statistically stationary state (with ∂/∂t annihilated by short-term temporal average) by the sideband component of the Reynolds/Maxwell stress. The electron parallel energy sideband representing J cos s evolves according to The first two terms are nonlinear ExB diffusion and magnetic flutter effects, usually found in the computations to be subdominant. The last three terms give the adiabatic and Alfvénic coupling to the pressure and vorticity sidebands (note the same terms with opposite sign in equations (41) and (42)) and the resistive dissipation, respectively. The coupling between J cos s and both W sin s and p sin s , noting that p − W is h e = p e − φ, completes the energetic description of the global Alfvén oscillation. Its rest state is the Pfirsch-Schlüter current balance, towards which it is damped by resistivity. The sideband Reynolds/Maxwell stress is required (equations (26) and (42) to disturb this balance. If this occurs, zonal flow energy may be ultimately damped by resistivity, after being transferred through the pressure sideband by geodesic curvature and to the current sideband by adiabatic coupling. It is most important to note that this can only occur in a two-fluid model, since this adiabatic coupling process is absent in an MHD model. In edge turbulence, the parallel ion flow (sound wave) sideband is found in the computations to be weak, since k c s is much slower than the turbulence and even (since q 3) slower than the geodesic acoustic oscillation. However, its energy equation is needed to close the system except for coupling to the turbulence. The sound wave sideband energy evolves according to Comparing with equation (41) we find the sound wave transfer to the pressure sideband (penultimate term on the right-hand side). The other three terms on the right-hand side are the damping effects: nonlinear ExB diffusion and magnetic flutter, and parallel viscous damping. In this model µ represents a substitute for ion Landau damping. There still remains the background thermal reservoir, given by the zonal pressure itself, whose gradient is the ultimate source for everything that transpires in this system. The evolution of the associated zonal pressure free energy is given by The first term is the direct drive of the turbulence, which self-adjusts to take energy out of the profile [44]. The other two terms are the geodesic coupling channels to the pressure and flow potential sidebands, respectively through toroidal compression of the diamagnetic and ExB flows.
The six equations (40)- (45) give the energetics of the complete dynamical system represented by zonal flows, the pressure and flow, the current and sound wave sidebands, and the thermal background, respectively. Since this is an electromagnetic system, the parallel electron dynamics has a magnetic inductive part (equation (3)). For even moderateβ within the range ofμ, the magnetic part is larger because of the extra factors of k 2 ⊥ in Ampere's law (equation (5)). Relaxation of the associated energy yields the Pfirsch-Schlüter current, which constantly adjusts itself on Alfvén timescales to a zonal pressure gradient whose changes are much slower. Nevertheless, the turbulence is able to pry this resistive 'Alfvén clamp' open sufficiently to cause some of the zonal flow/sideband energy to be dumped there. On parallel sound wave timescales, the parallel ion flow sideband and the pressure sideband would tend to keep each other in equilibrium with the zonal flow, if there were no turbulence.
The turbulence drives the zonal flow via Reynolds stress at a broad range of frequencies, keeping the zonal flow and pressure sideband at finite levels, with the (statistical) damping going through the turbulence, sound waves and the global Alfvén oscillation associated with the Pfirsch-Schlüter current. The principal energetic flow path is p →p →˜ → → p sin s → J cos s , through the gradient drive, adiabatic coupling, Reynolds stress, geodesic transfer and adiabatic transfer to the global Alfvén oscillation, which is resistively dissipated. Some of the energy in p sin s goes directly back to the turbulence via diffusive mixing of the pressure sideband and through sound waves to ion dissipation channels. The computations are required to determine the relative importance of all these processes.

Computational results
We will illustrate the functioning of all the elements in the zonal flow/sideband energy theorem with both 2D and 3D computations. The 2D case is useful as a pure control case illustrating the role of the Reynolds stress in the transfer of energy between the turbulence and the zonal flow mode, the role of the coupling mechanisms between the state variables in setting up the Reynolds stress/vorticity correlation, and the process of pure saturation in the Reynolds stress. The 3D model is used in both slab and toroidal forms; the slab version simply sets ω B → 0. The difference between the 3D slab and toroidal cases is found to be severe, but mostly due to the presence or absence of the geodesic coupling effect, which we define here as the k y = 0 component of the curvature operator. If this is put into the slab case or taken out of the toroidal case, the overall role of the zonal flows is found to be switched. The presence of the geodesic coupling effect is thereby found to inhibit the self-generation of large amplitude zonal flows. We note especially that this manipulation of the geodesic coupling effect is independent of the presence or absence of the interchange forcing effect, which is the action of the curvature on the k y = 0 component, that is, the turbulence itself. We therefore fully isolate the effect of geodesic coupling within the overall toroidal model.
It must be understood that these control tests do not involve dissipation or damping effects, but merely the subtraction of one or more energy transfer channels, as is clear upon examination of the equations describing the energetics, as they appear in section 5.
Further, we examine the statistics of the energy transfer processes, finding a displaced Gaussian for the probability distribution function (PDF) in most cases. The Reynolds stress drive and the geodesic transfer effect are found to have unambiguous directional tendencies, with the mean about two or three standard deviations away from zero in each case. Ultimately, the finding is that the principal saturation mechanism of self-generated zonal flows in toroidal geometry is the geodesic transfer effect. The ultimate sink of the pressure sideband energy is more subtle, with the indication being in favour of resistive damping of the global Alfvén dynamics to which the zonal flow/sideband system is coupled via the two-fluid adiabatic response.

2D results
The role of the various 3D phenomena are better understood when compared to a relatively simple 2D model, in which the basic precepts of the current zonal flow theoretical scenarios all function as predicted. We do this using a variant of the Wakatani-Hasegawa dissipative coupling model for collisional drift-wave turbulence [10]. The relevant curvature terms are added, noting that in a 2D model they represent purely interchange forcing. Other than this, the perpendicular dynamics of the state variables is the same as in the DALF3 model, facilitating separate assessment of the role of parallel dynamics and geodesic curvature in the latter. In the same language as for the DALF3 model, the equations are with the vorticity given bỹ with D a constant parameter representing the parallel electron dynamics, with p e under v E · ∇ representing the background gradient, (ω p is a switch, always set to unity except in appendix B), and with the curvature operator given solely by the interchange part, The terms involving D are done as in [46]; otherwise, numerical scheme is the same one as outlined for the 3D models in appendix A. Maintaining the accuracy, every term is evaluated explicitly, including the dissipative coupling, so that values of D are limited to 0.1 and smaller.
It is important here [47] to set D = 0 for the k y = 0 mode, since this allows for zonal flows (the use of D for k y = 0 follows from the non-vanishing of k for finite k y on closed flux surfaces [36]). We start with a random phase distribution forp e at RMS level of 3.0 (see [24] for details), withφ = 0. The grid is 64 × 256 nodes in {x, y} over a domain of L y = 4L x = 80π, and the boundaries are periodic. For cases with D = 0, we also startφ =p e with˜ set as in equation (48). The numerical hyperdiffusion coefficient (appendix A) was set to ν ⊥ = 0.01. All runs were taken for t = 1000, at which time snapshots were taken (unless otherwise indicated). Statistical quantities were averaged over 500 < t < 1000 with 50 samples. The timestep was 0.05 except for cases with D = 0.1, for which it was 0.005. The main point of this section is the functioning of the self-generation process for the zonal flows, specifically the the Reynolds stress/vorticity correlation. It is found to require coupling betweenp e andφ, as the control case with D = ω B = 0 found no zonal flow generation but merely the emergence of coherent vortices in the decaying 2D flow field [48], and the passive advection of the totalp e + p e . It is interesting to note that this numerical scheme is sufficiently well constructed to capture the coherent (isolated, spatially intermittent) vortex formation process even with a 64 × 64 grid resolution.
Cases with D = {0.00, 0.01, 0.1}, with ω B = {0, 0.05, 0.1, 0.2, 0.3} were examined. In all cases except D = 0 with ω B > 0, saturation at a finite turbulence and zonal flow amplitude was found. The suppression of the turbulence due to the flow was found to progress faster and remain deeper with the D = 0.1 cases than with D = 0.01. The saturation of the flows was always due to a balance in the Reynolds stress process between high-and low-k y spectral regions, with numerical dissipation for k y ρ s > 1 accounting for up to 30% of the flow sink. We define where · · · denotes the average over y. The flow generation process starts with R E > 0 as the initially random distributions of ˜ and ṽ x Eṽ y E statistically align. As the zonal vorticity amplitude increases, it begins to drive the low-k y component in the same 'Kelvin-Helmholtz assisted drift wave' process described in [3] for externally applied flow vorticity. At a particular amplitude of ˜ , the two spectral regions balance, tending towards R E = 0. This is the basic saturation process for the zonal flow vorticity in this geometry.
The Reynolds stress process itself requires coupling betweenφ andp e , as shown by the D = ω B = 0 case which did not result in significant zonal vorticity. The increased effect found with increasing D suggests a role for the adiabatic coupling: the zonal flow simply tiltsp e in the direction of the flow shear, with a sign sense such that ifφ were equal top e then the Reynolds stress energy transfer would be directed from the eddies to the zonal flow. In the absence of coupling, this has no effect onφ and hence no role in the Reynolds stress. But with finite D, the electron dynamics forcesφ to follow the structure ofp e , leaving the eddies in an energy losing relationship to the zonal flow. The interchange forcing has a similar but less direct effect, by each averaged over the spatial domain. All of A p , E e and A w representφ, but A p reflects mostly large scale flows and A w mostly the turbulence, with E e somewhat in between, due to the factors of ∇ ⊥ . The turbulence is very robust and rapidly approaches saturation, but just as quickly the A p signal rises and then the A w and Q e signals sharply drop and never recover. The fact that A p remains while A w is reduced shows thatφ concentrates towards the largest scales. This is the basic signature of spin-up-and-suppress in the time traces. The value for Q e was 0.565 ± 0.0744.      The spectra of the Reynolds stress energy transfer, given by R E (k y ), are shown for the case with D = 0.1 in figure 6, compared to the basic 3D slab case discussed in the next section. The plot is log-log, shown for both positive and negative values, and k y is in units of ρ −1 s . The result is that energy is transferred from eddies to flows at high k y , indicating suppression of turbulence, but from flows back to eddies at low k y , indicating energetic drive of large scale flow rolls. The zonal flow energy saturates when its shear becomes sufficient to create this long-wavelength drive Figure 6. Spectra of the Reynolds stress energy transfer (R E , labelled 'E') from eddies at each k y to the zonal flow, shown on log scales for both positive and negative values, for the nominal cases of the 2D and 3D slab models (the 3D cases are discussed in the next section). For the 3D case, the Maxwell stress transfer (R M , labelled 'M') is also shown. The saturation mechanism for zonal flow generation at finite adiabatic coupling levels is for turbulence suppression at higher k y to be balanced by turbulence drive at lower k y as the flow amplitude increases to statistical equilibrium.
region for the turbulence. The form of these curves is similar, indicating the same mechanism in both 2D and 3D slab cases. Previously, turbulence in the presence of prescribed ExB flow shear was found to exhibit this juxtaposition of drive and suppression regions if the flow shear was strong enough [3]. Both results are mutually consistent. The three dashed lines on the horizontal axis from left to right indicate the positions where 10, 20 and 50%, respectively, of the total positive part of the spectrum is accumulated. This demonstrates the importance of the k y ∼ 1 region to the positive part of R E , that is, the region in which the flows are driven by the eddies, even though most of the turbulence energy is closer to k y ∼ 0.1.
Of interest is the physical state in the phase during which the zonal flows are generated. This is very early, corresponding to about t < 100 in the time traces. A clear illustration of the emergence of the Reynolds stress flow drive is given by the profiles ˜ and ṽ   equal probability, their product R E is dominantly positive. This shows that the energy transfer proceeds from the eddies to the flows everywhere in the domain. The zonal flows emerge visibly somewhat later but are clearly visible by t = 70, as shown in figure 8. A time trace of R E confirms that it is robustly active throughout the phase during which the flows are generated and the turbulence is suppressed. At late times, R E is small, offset solely by the action of hyperdiffusion upon ˜ , and the zonal flow energy is saturated.
In the 2D model, then, the spin-up-and-suppress scenario works in detail as described by the theory [2], but only so long as one or both of the coupling processes (D or ω B ) are present. The physics of this situation is therefore quite different from purely hydrodynamic cases of spontaneous flow generation, so the universality of particular mechanisms as often cited [6] is questionable. The morphology and profile information give a clear picture of the appearance one finds when zonal flows dominate the turbulence although both are saturated at finite levels, and it is useful as a control when we turn to the 3D cases. The zonal flow shear (vorticity) amplitude saturates by the low-and high-k y spectral regions of the Reynolds stress/vorticity correlation coming into balance. We also find this balance in the 3D slab case. But in toroidal geometry, the geodesic transfer effect enters first, as we will presently show.

3D results-isolating geodesic coupling
Before the actual diagnosis of the physics as done for the 2D cases and as suggested by the zonal flow energy theorem in section 5, we give a more basic demonstration of its effect using the presence/absence test. The two basic cases are a toroidal and a slab model, given by the nominal parameter set in equation (18), with ω B = 0 for the slab case. A modified slab case is formed by putting the geodesic coupling effect into the slab model. A modified toroidal case is formed by taking this effect out of the toroidal model. We then have four cases, with and without the geodesic coupling effect and also with and without the interchange forcing present in both toroidal models. This test augments the one done previously with a different numerical scheme [9].
These four 3D cases are initialized equally: a random phase distribution forp e at RMS level of 3.0 in {x, y} and Gaussian envelope exp(−(8s/3π) 2 ) aligned to the magnetic field in s, screened from the x-boundaries by a relative decrement of exp(−(x ± L x /2) 2 /9), with the other dependent variables set to zero; the turbulence follows from the gradient driving and the adiabatic response [28]. The grid is 64 × 256 nodes in {x, y} over a domain of L y = 4L x = 80π, and 16 nodes in s over a domain of one connection length (−π < s < π). The boundaries are Dirichlet in x, periodic in y, and pseudo-periodic following global consistency and the shifted metric treatment in s as detailed in [37].
The time traces of the four cases at a single set of parameters, displayed in figure 9, shows the basic result of isolating these curvature effects. For the four cases, the time traces of A n , A p , A w and Q e are given, respectively representingp 2 e /2,φ 2 /2,˜ 2 /2 andp eṽ x E , averaged over all grid nodes in the 3D domain. Due to the adiabatic response,φ generally tracksp e for the k y = 0 modes, so large differences in A p and A n are a convenient way to track the zonal flows in these time traces. As in the 2D cases, A w most directly tracks the turbulence. The clearest finding from the appearance of these traces is that the two cases with the geodesic coupling effect present saturate well, while those which lack this effect continue to evolve at relatively late times (t > 2000). The numbers and diagnostics quoted below are averaged over the interval 2000 < t < 4000 for all cases except for those without geodesic curvature, for which the interval was 6000 < t < 8000. They are in gyro-Bohm units, so that A p = 2.0 would imply an RMS amplitude level of eφ/T e = ρ s /L ⊥ , and Q e = 1.0 would imply a transport diffusivity of χ GB = ρ 2 s c s /L ⊥ . In terms of this nominal parameter set, these are ρ s /L ⊥ = 0.0142 and χ GB = 0.45 m 2 s −1 , respectively.
In the basic toroidal case, the turbulence and the flow amplitude reach statistical saturation at a robust transport level, measured at A n = 10.3 ± 0.172, A w = 0.703 ± 0.052 and A p = 12.7 ± 3.12, while the transport stabilizes at Q e = 1.01 ± 0.117. These figures set the baseline against which the other cases are assessed. The time traces show that the native turbulence saturation time (neglecting sound wave and zonal flow energy) is relatively short, t < 500. In the basic slab case, the turbulence saturates in much the same way, but the flow amplitude continues to build and eventually the whole reaches statistical quasi-saturation (evolving with the flow amplitude) at a substantially lower transport level, measured at A n = 3.84 ± 0.17, A w = 0.495 ± 0.015 and Q e = 0.556 ± 0.021, with the flow amplitude at A p = 322 ± 5.9, saturating only well after t = 4000. The principal difference between the slab and toroidal cases is this behaviour of the flows.
The modified slab case finds a result with the same character as in the basic toroidal case: fast saturation, robust turbulence, moderate flow amplitude, robust transport level, measured at A n = 4.65 ± 0.34, A w = 0.447 ± 0.015, A p = 4.60 ± 0.43 and Q e = 0.528 ± 0.026, respectively. The change from the basic slab case to this one shows that the geodesic coupling introduced alone is decisively able to saturate the zonal flow level.
The modified toroidal case finds a result with the same character as in the basic slab case: slow saturation, weakened turbulence, strong flow amplitude, weakened transport level, measured at A n = 4.14 ± 0.21, A w = 0.515 ± 0.017, A p = 469 ± 13 and Q e = 0.574 ± 0.026, respectively. The change from the basic toroidal case to this one shows that geodesic coupling  figure 9. The two cases without geodesic curvature (slab and mod tor) find high-amplitude sinusoidal potential layers.
is not only able to but required to saturate the flow amplitude at the low levels A p < 100 seen generally in all the basic toroidal cases (at the variousβ and C discussed below).
The profiles at t = 4000 (at t = 8000 for the slab and modified toroidal cases) ofφ and˜ , averaged over all nodes in y and s, are given in figures 10 and 11, respectively. The vorticity profile directly measures the shear in the zonal flow. We find the results corresponding with the time traces above: the basic slab case and the modified toroidal cases are the ones with much stronger zonal flow level, with a single sinusoidal structure in φ(x) also seen in 3D periodic/unsheared slab cases [49]. Whether the end state of the zonal flow shows a single radial layer or more of them is found to be a delicate result of both boundary treatment and parameters. Particularly, more electrostatic (lowerβ) and more adiabatic cases (lowerβ and lower C) show narrower zonal flow layers and more of them.
The spectra in terms of k y of the four cases are shown in figure 12. They are all similar, broadband primarily within the range 0.1 < k y ρ s < 1.0, with the energy collecting towards the longer wavelengths, the transport intermediate and the zonal Reynolds stress drive at shorter wavelength. The entrophy spectrum is flat nearly to k y ρ s = 1. The fluctuation free energy is made up of four pieces, one for each dependent variable, following the ExB energy (|∇ ⊥ φ| 2 ), the thermal free energy (p 2 e ), the magnetic and electron kinetic energy (β|∇ ⊥Ã | 2 +μJ 2 ) and the sound wave kinetic energy (ˆ ũ 2 ), each divided by 2 and averaged over the grid nodes. The k y spectrum is the contribution due to each k y component to the total [28]. The transport is similar, Figure 11. Profiles (zonal averages, over y and s) of the ExB vorticity for the four cases in figure 9. The two cases without geodesic curvature (slab and mod tor) find high-amplitude sinusoidal shear layers. computed forp eṽ x E , and the zonal Reynolds stress drive (cf next section) is the zonal vorticity times the zonal Reynolds stress ṽ x Eṽ y E . The transport curves turn downward at slightly lower k y for the cases without geodesic coupling (strong flows), following the energetic suppression of the shortest wavelengths by the flow shear (cf [3]). Otherwise, the forms are the same. The predominance in each signal of a different region in the spectrum reflects the nonlinear nature of the turbulence.
The nominal case, with both interchange and geodesic curvature effects, shows a peak near k y ρ s = 0.1 superposed upon the rest of the broadband spectrum. This is a nonlinearly supported MHD mode, as diagnosed below. The region k y ρ s < 0.125 is found to account for 20% of the total transport, and either k y ρ s < 0.35 or 0.15 < k y ρ s < 0.55 for 50% of it.
It is important to note the strong contribution of the smaller scales to the zonal Reynolds stress drive, and to compare that to the spectrum ofp 2 e (which is very similar toφ 2 for k y ρ s < 1). This is a very nonlinear situation to which any discussion of a dominant mode or linear instability is poorly posed.
The parallel envelope structure (cf [38]) show the basic toroidal and slab forms, and in this case the presence or absence of geodesic coupling has no effect; only interchange forcing does. But the otherwise similarity of the basic toroidal and modified slab cases show that this 'ballooning' signal is not germane to the actual dynamics of the turbulence.  respectively. The stronger zonal flows reduce the scale of motion approximately isotropically, with the larger scale interchange driven range in the basic toroidal case seen in the spectra suppressed for the modified toroidal case and absent for the slab cases.
The morphology in the four cases is shown in figure 13. Only half of the y-domain is shown. The basic slab case shows dominance of the shear layer in the potential, withp e noticeably tilted into the y-direction. The basic toroidal case shows a faint shear, but its vorticity is not larger than that of the turbulence, which is why a strong amount of suppression does not occur. The cases show strong or weak shear layers according to whether geodesic coupling is absent or present, respectively. The shear levels of these weaker flows are comparable to the dynamical frequencies of the turbulence (about 0.1c s /L ⊥ ). The A p time traces (cf figure 9) for the weaker flows show them to be short-lived, comparable to the correlation time of the turbulence (about 6L ⊥ /c s ). These are the zonal flows which remain as part of the turbulence, leading in fact to moderate suppression but allowing it to remain at a robust amplitude. By contrast, the flow layers in the cases without geodesic coupling persist for the length of the run, and hence are more like Figure 13. Spatial morphology of the electrostatic potential and pressure disturbances for the four cases in figure 9. The two cases without geodesic curvature (slab and mod tor) find clear flow shear layers which cause the pressure disturbances to tilt with the zonal flow. The other two cases find nearly isotropic turbulence only marginally affected by the flows, and potential amplitudes comparable to the turbulence. mean flows even though they are self-generated. The basic toroidal case also shows the large scale interchange driven range in the superposition of large 2 < k y L y < 4 structures upon the turbulence. This is the spatial manifestation of the MHD component mentioned above.
The average ExB energy transfer spectra [38], in terms of k ⊥ defined as k 2 ⊥ = k 2 x + k 2 y (since the perpendicular metric is unit diagonal), are shown for the four cases in figure 14. These are the time-averaged contribution of each {k x , k y } component to the three transfer terms in the vorticity equation: nonlinear three wave transfer within the ExB energy through the eddy Reynolds stress, mainly linear transfer against the magnetic energy through the parallel currents, and strictly Figure 14. Energy transfer spectra for the four cases, showing average transfer into ExB eddies at the given k ⊥ due to polarization nonlinearity ('e'), parallel currents ('j') and interchange curvature ('k'). Interchange curvature is absent for slab cases, but the curves show how large it would be if it is switched on. The 'e' and 'j' curves show the basic drift-wave turbulence dynamics. For k ⊥ ρ s > 0.2 the polarization current is responsible for maintaining energy transfer into the ExB eddies through the parallel current, indicating the basic drift-wave mechanism maintaining the turbulence. The interchange forcing represented by the curvature enters into the toroidal case for k ⊥ ρ s ∼ 0.1.
linear transfer against the thermal free energy through the overall curvature term (interchange forcing). Only the contributions not involving the k y = 0 modes are shown. Note that for the two slab cases, the interchange forcing is zero, but the green lines show where the result would fall if the interchange curvature suddenly switched on. The modified slab case is important as a control case against the interchange forcing: the similarity of the other two forms shows that the main process maintaining flows and currents in the turbulence is the nonlinear slab drift wave physics, and that the presence or absence of interchange curvature forcing has no effect. So in general, in the basic toroidal cases, the interchange curvature is a perturbative effect adding to the general drift wave turbulence, producing ballooning but not affecting either dynamics or otherwise mode structure.
The frequency spectra of the zonal potential and zonal flow for the basic and modified toroidal cases is shown in figure 15. In both cases the zero frequency component is most prominent, but the presence of the geodesic coupling leads to activity up to about the nominal geodesic oscillation frequency of 0.035. It is important to note that no single coherent mode is dominant (and also that one might be misled in that direction by a linear-linear plot; cf [30]). The prominence of the frequency peak varies case by case (with parameters as discussed below), however, and a strong peak is not a necessary signal for the role of geodesic coupling. The zonal vorticity is forced upon by the turbulent Reynolds stress at all frequencies, and this is reflected in the noisy nature of the zonal potential frequency spectrum. With the zonal flow reflecting somewhat smaller scales due to the extra k ⊥ factor, the geodesic frequency peak is less pronounced. So although the geodesic transfer effect is clearly active energetically, this will not always result in a frequency peak visible by experimental diagnosis. The absence of such an observation does not therefore conclude the absence of the geodesic effect, even if the diagnostics should be able to discern between the signatures shown in figure 15. On the other hand, the activity of all frequencies except zero is sharply reduced in the absence of geodesic curvature.
Having performed presence and absence tests on both geodesic coupling and interchange forcing effects, we find that the main qualitative difference in drift Alfvén turbulence between slab and toroidal geometry is not interchange forcing, which produces linear instabilities, but geodesic coupling, whose role is to statistically saturate the zonal flows, hence act as geodesic transfer. Whether this effect actually enters with a pronounced geodesic acoustic oscillation is not central to the zonal flow result: the main point is that the geodesic coupling is effective in holding the zonal flow vorticity to levels not greater than that of the turbulence, and hence prevent the emergence of strong self-generated flow layers with enough vorticity to suppress the turbulence.

3D results-energetics and statistics
In this part we diagnose the energetics of the zonal flow/sideband system using the results of the analysis of section 5. Only the 'basic toroidal' cases are considered, i.e. those with all terms present in the equations. The emphasis is on the variation ofβ, i.e. how electromagnetic the adiabatic response of the parallel electron dynamics becomes [28]. Nominally, this sets in whenβ = 1, which is when the drift and Alfvén transit frequencies c s /L ⊥ and v A /qR are equal. However, the competition is between the electron inertia/resistivity and the inductive electric field in Ohm's law. This comparesβ toμ, that is, the traditional comparison between the electron dynamical beta (β e = c 2 s /v 2 A ) and the mass ratio (m e /M i ). But due to Ampere's law relatingÃ toJ , the mass ratio enters in combination with an extra factor of k 2 ⊥ ρ 2 s , so in the context of turbulence the electromagnetic induction begins forβ rather lower than either unity orμ. For linear waves there is a partial cancellation between ω and the diamagnetic frequency ω * (in normalized units the latter is just k y ), but for the turbulence ω has a broad distribution even for each individual wave, and so effects must be investigated using the computations. For the turbulence, the electromagnetic influence starts even at the lowestβ values which are affordable by an explicit computation. Asβ rises, a narrow range of modes with MHD character with k y ρ s ∼ 0.1 is superposed upon the drift wave turbulence with roughly 0.2 < k y ρ s < 1. The scale of motion of both components is independent ofβ, and below the transition there is little effect on either turbulence or transport. Forβω B > 1, the MHD component takes over and the drift wave component loses significance. In other words, this is not a linear boundary even though the transition is rather abrupt. The zonal flow energetics, however, starts changing somewhat before this boundary, as shown below.
In the diagnostics presented below, the zonal average · · · is over the grid nodes in y and s, and the radial average {· · ·} is over the grid nodes in x. All statistical results are averaged over the time interval 2000 < t < 4000.  (40). The time traces of R E and T G show a very robust signal with unequivocal sign, and are close to balance, while R M is small in both average and standard deviation. The signals are very unsteady, which is why it is necessary to carry these runs to t = 4000. The burst-like appearance is deceptive, however, as the probability distribution functions for these two quantities sampled in time, shown in figure 17, has a mainly Gaussian shape, although shifted about three standard deviations away from zero. Sampled in both x and t, however, both quantities show strong power-law tails to the positive side. These results confirm the basic energy transfer path for the zonal flows: they are driven by the Reynolds stress and depleted by the geodesic curvature effect.
The pressure sideband energy theorem is more complicated however. Not all of the energy goes back to the turbulence, nor is it dissipated by MHD sound waves. Recall that the zonal flow and sideband energy is energetically coupled to the global Alfvén oscillation. This goes into and out of the pressure sidebands through the 'geodesic diamagnetic coupling' and 'adiabatic coupling'given by the last two terms in equation (41), respectively; both are inherently non-MHD processes. It is necessary to measure all the terms in equation (41) for a complete accounting. We find that this global Alfvén dynamics has comparable energy transfer to the geodesic acoustic dynamics. Energy enters the pressure sidebands via geodesic transfer and geodesic diamagnetic coupling at about equal rates, and it exits mostly by adiabatic coupling into the magnetic sideband energy. In this case, the nonlinear decoupling gets randomized by all the Alfvén dynamicsit is robustly active but averages close to zero. The other terms, including the sound waves, second sideband (sin 2s) pieces and the magnetic flutter, are smaller. The transfer into and out of the pressure sideband energy is the result of all the terms in equation (41), as follows (all ×10 −3 ): diffusive mixing, out to turbulence: −0.67 ± 1.48; geodesic coupling, in from zonal flow energy: 3.00 ± 0.89; parallel MHD compression, out to sound wave sideband: −0.28 ± 0.08; adiabatic coupling, out to current sideband: −1.45 ± 0.45; diamagnetic compression, in from zonal pressure: 0.56 ± 0.28; and dissipation in edge diffusion layer: −0.70 ± 0.29. Transfer to the second sideband (cos 2s) via diamagnetic compression is small: 0.009 ± 0.013. The discrepancy from 100% is well within the error bars and consistent with saturation over the finite time of 4000L ⊥ /c s for the run. A test run with the edge diffusion applied to the zonal pressure only instead of all k y = 0 modes found that most of this 0.70 went into an increase of the adiabatic transfer toJ and about 20% into diffusive mixing back to the turbulence.
The next piece is the magnetic sideband energy, in equation (43). The nonlinear terms are found to be small and the only input is adiabatic transfer from the pressure sideband (+1.45 ± 0.451), with the energy exiting through Alfvén coupling to the flow sideband (−0.581 ± 0.394) and resistive dissipation (−0.866 ± 0.234). As expected, this component is in a nearly quasistatic balance, although the standard deviations are about one-third of the total throughput. The next piece is the flow sideband energy, in equation (42). In a slowly decaying equilibrium, this component is maintained by the Pfirsch-Schlüter current balance, driven by diamagnetic compression and damped by Alfvén transfer to the magnetic sideband. With the turbulence present, however, the global Alfvén dynamics reacts to it, so that the sign of respectively. Below these boundaries, the drift wave regime finds a balance between R E and T G , which extends into the resistive MHD regime albeit with larger error bars, following the more vigorous turbulence controlled by fewer k y modes at larger scale. Into the ideal MHD regime, both T G and R M change sign and R E is forced to be negative. The balance R E + R M = 0 which would hold for shear Alfvén waves is broken by T G , whose changes appear to control these transitions, although R M is the largest driving effect forβ 20.
these transfer effects is reversed: energy enters from the pressure sidebands through the Alfvén coupling (+0.581 ± 0.394), and leaves to the thermal background via diamagnetic compression (−0.185 ± 0.380), and to the turbulence through the Reynolds stress sideband (−0.263 ± 0.362). The Maxwell stress sideband is very small (−0.015 ± 0.055).
The sound wave sideband is maintained by the small fraction of the energy exiting the pressure sideband through parallel MHD compression (+0.280 ± 0.081), which is dissipated mostly linearly by the artificial dissipation given by ν , which in this case acts in place of µ which was set to zero. Diffusive mixing (−0.083 ± 0.023) accounts for the rest.
The zonal pressure drives the turbulence by diffusive mixing (−962 ± 106, much more robust than all the zonal flow/sideband dynamics), and drives the global Alfvén oscillation in a net sense (−0.561 ± 0.277 for diamagnetic compression to the pressure sidebands and +0.185 ± 0.308 in from the flow sidebands). The magnetic flutter transport is small and negative, with the transfer to the zonal pressure (+0.123 ± 0.079) much smaller than the ExB drive of the turbulence. It is interesting to note that the resistive dissipation which allows the magnetic sideband to take energy out of the flow/sideband system also results in the pressure and flow sidebands being driven by the zonal pressure-the sidebands are non-adiabatic. But the transfer is forced by the energy leaving the zonal flows for the pressure sidebands, which causes the adiabatic transfer to be so large, which reverses the sign of the pure Alfvén part of the overall global Alfvén oscillation.
It is very interesting to know how the main zonal flow drive and depletion rates scale with parameters, especiallyβ but also C. This scaling is shown in figure 18, for the zonal Reynolds stress drive R E , the geodesic transfer T G and the zonal Maxwell stress drive R M , the latter given byβ −1 { ˜ b xby }, as in equation (40). For all collisionality regimes, the basic roles remain unchanged: R E is positive, T G negative and R M forβ = 2 remains weak. For C > 10 the strength of these effects increases, mostly due to the increasing strength of the turbulence itself, but also due to the increasing impact of the longer wavelengths. The error bars are larger in this regime, following the increased intermittency when fewer k y modes are involved in the energetics. Note that the nominal resistive ballooning (MHD) boundary is Cω B = 1, which for this choice of ω B = 0.05 is C = 20.
The dependence uponβ is more complicated: asβ rises, the first thing to change is the geodesic transfer itself, approaching zero and then changing sign atβ between 7 and 10. The R M appears abruptly atβ = 10, but clearly after the observed changed in T G , suggesting the latter is driving the changes. A third regime is reached when forβ = 20 and above R M changes sign and becomes the principal zonal flow drive, while R E is forced to change sign and become the sink. This pattern of changes suggests R M is now driving the changes. Note that the nominal ideal ballooning (MHD) boundary isβω B = 1, which for this choice of ω B isβ = 20.
We next consider the Reynolds and Maxwell stress spectra, shown with respective labels 'E' and 'M' in figure 19 for the nominal case withβ = 2 and a 'high beta' case withβ = 10. The nonlinear interaction of the vorticity itself is limited to scales near ρ s by the adiabatic response [50,51], and even when the electrons are not adiabatic the parallel dynamics selectively constrains the longer wavelengths (the reason the drift wave nonlinear instability loses influence at the longest wavelengths). The Reynolds stress of the drift wave part of the spectrum (about k y ρ s > 0.2) is influenced by these effects. Measuring this spectrum by considering the contribution of each ±k y pair to R E , we find that it indeed peaks very near ρ s . Summing the Fourier components, we find that half of the Reynolds stress drive is done in the spectral range k y ρ s > 0.675, and 80% of it is done in the range k y ρ s > 0.45.
Analysis of simpler models suggests that for large enoughβ we should see R M arise and cancel R E [42]. This follows the fact that shear Alfvén waves transport no momentum (for them the cancellation is exact). Indeed forβ 10 we find a substantial fraction of the energy going into Alfvén waves, since the parallel structure ofφ (squared amplitude versus s) becomes flat while that ofp e becomes strongly ballooned-by contrast the Cω B > 1 regime shows a more conventional ballooning structure in both these state variables. The form of the spectra of R E and R M indicates this cancellation tendency even when R M R E , showing that it is a basic part of the drift wave mode structure. Asβ is increased to 10 the cancellation becomes affective. Forβ = 2 the total positive part of R E is 0.00234 and the total negative part of R M is 0.00015. Forβ = 10 the total positive part of R E is 0.00271 and the total negative part of R M is 0.00280.
But for the controlling effects of T G in the intermediate regime and the fact that R M rises past the level necessary to cancel R E , we have to look elsewhere. Quantitative analysis of the Pfirsch-Schlüter Alfvén system shows where: before the Maxwell stress cancellation becomes effective, the two main transfer channels describing the Pfirsch-Schlüter current balance in the flow sideband energy (equation 42) change sign. It indicates that with rising beta the turbulence begins putting more energy into the global Alfvén dynamics than the Pfirsch-Schlüter current drive does. Now the (measured) transfer path is eddy turbulence into zonal flows (Reynolds stress), into magnetic disturbances (Maxwell stress) as well as the pressure sideband (geodesic transfer), from the latter into the magnetic sideband (adiabatic coupling), from the latter into the flow sideband (Alfvénic coupling), and hence backing into the thermal background (sideband geodesic coupling), and finally from there into the pressure sideband (diamagnetic coupling). Ultimately, this causes the main geodesic transfer (the one affecting the zonal flow energy) to back up, eventually changing sign and causing the less intuitive effects in figure 18.  figure 6, to whose 2D and 3D slab cases these ones should be compared. There is a clear tendency for the Alfvén cancellation between R E and R M , with the relative amplitude varying withβ. At moderate β values, the Maxwell stress for both the slab and toroidal cases follows the Reynolds stress, although R M is smaller. In the toroidal cases, R E does not selfsaturate as in the slab cases. For moderateβ, the geodesic transfer balances a continued robust R E , and for largerβ the Maxwell stress enters quantitatively. At the entrance to the MHD regime (β = 20), the Reynolds/Maxwell stress cancellation is effective, but then the energetics is taken over by the Pfirsch-Schlüter current sideband system.  figure 18 that the strength of R E and T G mostly follows that of the turbulence itself.
Putting these diagnostics together, we find that even for nominal parameters well below the MHD beta limit the turbulence is robustly electromagnetic enough to knock the Pfirsch-Schlüter currents out of equilibrium so that their involvement in the energetics is comparable to the zonal flows and geodesic acoustic dynamics. Moreover, the drift wave component dominates the overall Reynolds stress zonal flow drive. One could not have obtained these results with an MHD model, nor with one which does not clearly resolve the drift scale.

Conclusions-drive and saturation of zonal flows
This study has focused upon drift Alfvén turbulence and zonal flows under tokamak edge conditions in which direct dissipation mechanisms in the perpendicular dynamics are very weak. In gyro-Bohm units, the collisional dissipation mechanisms in the electrons scaling with ρ 2 e ν e are of order 10 −3 and if warm ions are considered the dissipation in the ions scaling with ρ 2 i ν i are of order 10 −2 . These are comparable to or smaller than the numerical dissipation required to maintain statistical saturation in the turbulence with a perpendicular grid spacing comparable to the drift scale, ρ s . In these computations, these collisional dissipation mechanisms are therefore neglected. Moreover, parallel viscosity, which with the small values of L ⊥ /qR would be comparable to or smaller than the numerical parallel dissipation. Zonal flow damping is therefore forced to go either through the nonlinear cascade dynamics in the variables other than the vorticity (henceφ), or resistive dissipation in the Alfvén dynamics. The spurious damping problems brought up in the context of core turbulence [7,32] are therefore avoided. The principal findings of the study are as follows.

Drive of self-generated zonal flows
The zonal flow energy theorem is unambiguous as to the possible transfer channels by which energy can enter the zonal vorticity, hence flow shear. In both slab and toroidal models, the Reynolds stress drive, the correlation between zonal vorticity and zonal Reynolds stress, is found to be robust and responsible for the rise and maintenance of the zonal flows. Almost all of this drive is in the drift wave part of the spectrum, with the median k y ρ s approximately 0.7. This has been found to be a robust result, even given that the upwind method [26] used in [9] has been found to overestimate the Reynolds stress. Adiabatic coupling of the electron pressure and electrostatic potential through parallel currents is found to greatly strengthen the operation of the Reynolds stress mechanism: flow shear tiltsp e structures, the adiabatic response forcesφ to assume a similar orientation, resulting in the ExB eddies losing energy to the flow. Further details are in the statistics: with no zonal flow, the zonal Reynolds stress PDF is random, nearly Gaussian. Fluctuations in the zonal vorticity are also random. Where they correlate, the zonal flow energy is increased in the same locations. Where they anti-correlate, the zonal flow energy decreases. A positive correlation therefore develops naturally, much in the same way that a positive transport develops in the presence of a density gradient but in the absence of other forcing mechanisms [44]. The zonal flow energy is prevented from growing without bound, by the following saturation mechanisms.

Saturation of self-generated zonal flows
The zonal flow energy theorem finds that only two saturation mechanisms are possible, if the Reynolds stress is a drive: Maxwell stress or geodesic curvature. In slab geometry at the same moderate beta values (βω B < 1 for the toroidal cases), the Maxwell stress is relatively weak and the geodesic curvature is absent. The zonal flows can only saturate if the Reynolds stress drive self adjusts to zero. In slab geometry, that is what happens: the flow amplitude grows to sufficient amplitude that the flow shear drives the turbulence at longer wavelength (k y ρ s ∼ 0.1). This is the same result found earlier in simpler but still well-resolved drift wave models [3]. In toroidal geometry, however, the geodesic transfer mechanism always sets in before the zonal flows become so strong-the Reynolds stress drive remains robust, as it is measured to do in experiments [53]- [55]-but it is balanced by geodesic curvature and therefore the zonal flow energy and shear stop growing. In situations of experimental interest to edge turbulence (βω B < 1 and C > 1 and ω B < 0.1), this balance always establishes itself and hence the geodesic transfer mechanism is the principal mechanism of zonal flow saturation. The Maxwell stress remains interesting however: it always has the required spectral form to cancel the Reynolds stress even when it is too weak to do so. Forβ = 10, the Maxwell/Reynolds stress cancellation becomes effective, but this is also the effective beta limit. In actual experiments it remains possible that since the β limit is higher (due to finite aspect ratio and shaping [26]), there might be a window between the β values of stress cancellation and ideal ballooning onset, but this should be considered with a full gyrofluid model (required to resolve ρ i as well as ρ s ) and with a detailed tokamak geometry and is left to the future. Within the present model, the saturation of zonal flows by the geodesic transfer mechanism has also been found to be robust enough to prevail in both the earlier and present numerical schemes.

Zonal flow damping
The saturation mechanism found from the zonal flow energy theorem itself is the geodesic transfer, but for this to work the pressure sidebands must also be in balance, between the energy incoming from the zonal flows through geodesic transfer and some sink mechanism. In this tertiary process, the two numerical schemes have indeed differed. It was [9] found that nonlinear diffusive mixing was stronger. But the results herein have shown that more energy (about half the total) leaves the pressure sidebands via adiabatic coupling to the current sideband and is ultimately dissipated through resistivity, and the rest goes through diffusive mixing and a somewhat smaller part through sound waves which are ultimately also depleted through diffusive mixing. The part of the sideband energetics involved in the adiabatic coupling also involves drawing energy in from the background through diamagnetic compression. The turbulence/zonal flow/sideband system is therefore in energetic contact with the global geodesic Alfvén mode, which if allowed to relax undisturbed would result in the Pfirsch-Schlüter current. The sideband Ohm's law is indeed in close balance between the three linear terms: adiabatic and Alfvénic parallel dynamics and resistive dissipation. The vorticity sideband energy equation shows that the sideband Reynolds stress can upset the Pfirsch-Schlüter current balance sufficiently that the pressure sidebands can use the current sidebands as a sink. The Arakawa-Karniadakis scheme of [56] used herein is indeed superior, and where the results in [9] using the upwind scheme of [26] differ with the present ones, the latter are judged to be the more likely correct.

Other considerations
The drive and saturation mechanisms of the zonal flow energy itself are unambiguous, as found by their PDFs which show mean values of about 2-3 SD from zero. When the less realistic high β MHD regime is entered, the result is not as clear-first the overall drive is weakened, and second there is much more energy throughput in the global Alfvén dynamics than in the zonal flow/sideband system itself, so that the latter is no longer in control of the entire sideband dynamics. The average signs of both the Reynolds stress and the geodesic transfer are found to be reversed forβω B > 1, but the SD is now larger than the mean, and certainly no energetic transfer quantity outside the zonal flow energy equation has been found to have an unambiguous sign. One can take the position that the MHD equilibrium itself breaks down since the shear Alfvén dynamics is no longer able to maintain the sideband balances. Even Ohm's law is seriously disturbed in this regime. It may be physically interesting to pursue this further, but it is nevertheless an unrealistic regime.
It remains of interest what happens when the edge situation is in more of an ITG than a drift-wave regime, which occurs when η i = η e 2 according to the mode-structure signatures [38]. This is left to future work with the gyrofluid model, since the fluid model cannot treat it with this numerical scheme (the virulent cascade dynamics found when the ExB inverse cascade breaks down due to gyroviscosity forces the use of artificial dissipation values which are so high that the drift wave dynamics is no longer correctly reproduced). Some situations of experimental interest do however exhibit steeper density gradients, and even cold ions [57], though, so the present analysis is directly relevant there. It remains that a significant fraction of edge turbulence computational work is done with models even simpler than DALF3 and with strongly insufficient resolution, and perhaps with flawed numerical schemes. In this context, then, the considerations raised by the present work are very germane.
In terms of experimental interpretation, it is important to note that the circumstantial observation of the Reynolds stress process, experimentally [53][54][55] or even computationally [58], and even the finding as herein that it energetically drives zonal flow self-generation, do not imply that transport barriers observed in tokamaks to be coincidental with strong ExB flow shear [59] are actually caused by the Reynolds stress process. We have found herein that this process is balanced in toroidal geometry by the geodesic transfer effect, with the zonal vorticity held at moderate levels even in the face of robust Reynolds stress drive. The explanation for the actual formation of transport barriers involving ExB flow shear in tokamaks requires additional components, perhaps arising from the details of the equilibrium shear layer, as has been argued before [16,26].
Recent experimental evidence of geodesic acoustic oscillations [60]- [62] provide an interesting signpost that the geodesic compression of ExB flows is active, even if the oscillation itself (which is purely conservative) does not imply the presence or absence of the geodesic transfer effect on zonal flows.
A recent paper has looked at the zonal flow energetics in the core [63]. However, it does not contain enough toroidal modes in its calculations to allow the possibility of the Reynolds stress to have a different spectrum than the turbulence itself as it usually does. The low resolution in turn forced the use of a very unrealistic value of a/ρ s = 80, with an active region of only 0.6a = 48ρ s , that is, even less than the present edge turbulence computations. Finally, the use of background profile functions in the curvature terms caused the geodesic coupling process not to conserve energy (radially dependent profile functions do not commute with K x ). These processes are subtle and involve a high degree of scale separation. Computations which keep and resolve realistic scales are necessary. With present day resources, this also requires the use of field-aligned coordinates, as in these computations and those in [5].
Parallel to this paper is also an analysis of this problem by Naulin et al [64], which is in apparent disagreement with some of these results, especially at high β. That part of the problem is very delicate however, and the discussion in appendix A is germane. Specifically, the dissipation is too large at this resolution unless hyperviscosity is used, and the high-k ⊥ and -k dissipation should be applied to all quantities advected by the ExB velocity, not just the density/pressure and vorticity. However, the desired application is to plasmas which do not exhibit turbulence at the level of what could be called disruptions; that is, the regime of interest is below the β values at which large-scale MHD ballooning is excited, and in this regime the results are in agreement.

Appendix A. Notes on computational methods
Some notes concerning the computations in both the 2D and 3D models are necessary, as convergence of the zonal flow phenomena has been found much more subtle than that of the transport and turbulence mode structure itself. The reason is that while the drift wave nonlinear instability involves the entire range roughly 0.1 < k y ρ s < 1.0 [38], the zonal Reynolds stress is contained to a narrower region closer to k y ρ s = 1. Moreover, the zonal flow self-generation process through the Reynolds stress/vorticity correlation depends critically upon the underlying properties of the bracket structure describing the quadratic nonlinearities in the equations. We require a scheme which (1) is stable to hyberbolic advection phenomena, including waves, (2) confines the numerical dissipation to the smallest scales, and (3) leaves the numerical dissipation coefficients fixed when the parameters are changed.
The upwind method used in the previous computations with this model, specifically in [9] and outlined in [26], has been found generally problematic on this problem, mainly owing to its overestimation of the Reynolds stress but also in the dependence of the zonal flow transfer mechanisms upon theβ parameter, especially when entering the MHD regime (βω B > 1, the usual MHD ballooning criterion). An alternative method is available and is generally wellbehaved [56]. For  (and similarly for magnetic nonlinearities using −βÃ in place ofφ), the spatial discretization is one by Arakawa which preserves the symplectic structure [65]. Linear derivatives, ∂/∂s and K, are evaluated with centred differences. As in the upwind scheme, the combinationp e + p e is always discretized together. The time step is a third-order expansion derived by Karniadakis et al, using both the dependent variables and the right-hand sides of the equations, thus incorporating all mixed derivatives in the Taylor expansion to specified order [66]. This time step scheme is stable for waves, which removes the objection to the standard predictor-corrector methods (with derivatives done with either Fourier transforms or centred differences) which motivated the use of the upwind scheme in the first place.
In the equations, the only dissipation is resistivity (CJ ) and parallel viscosity (µ ). These are however insufficient to contain 3D turbulence, due to the vigorous cascade through v E · ∇p e [38]. Classical and neoclassical dissipation coefficients are too small to do this unless the physical dissipation range wherein ρ 2 e ν e k 2 ⊥ > ω is resolved. The range which must be resolved is down to a grid node spacing of h x = h y = ρ s , for which k ⊥ ρ s > 1. This is to provide a dissipation-free range all the way down to k ⊥ ρ s = 1. With the upwind scheme, the dissipation is part of the scheme and moreover increases for more vigorous turbulence. This work has found it necessary not only to keep the dissipation coefficients constant with parameters, since the nonlinear processes affecting the zonal flows always reach the high-k ⊥ range, but also to use a hyperdiffusion for the perpendicular coordinates, thereby narrowing the dissipation range necessary to provide saturation. Dissipation in the s coordinate is also necessary, as the direct cascade by v E · ∇p e also acts through k [67]. This is effected by adding the following dissipation operators to d/dt: where the ν ⊥ and ν terms are implemented in explicit form, but ordered in the scheme in the same way as implicit step parts are done in [66]. The values for every case to be computed are ν ⊥ = 0.01 and ν = 0.003. This means that the artificial dissipation for every degree of freedom for which k ⊥ ρ s < 0.9 is less than 10 −2 . Such a strictness in resolution is probably unprecedented in a set of 3D drift wave computations. It is very important to add the artificial dissipation to d/dt, which includes v E · ∇, but not ∂/∂t. This means the dissipation never acts directly upon the fieldsφ andÃ , but only upon the fluid variables,˜ ,J ,p e andũ . Furthermore, the edge dissipation incorporated to maintain the background gradient is not applied to˜ but only to the k y = 0 components ofp e , that is, not only upon the turbulence, but also not upon the zonal flows. The only boundary dissipation used is the simple damping of the k y = 0 component ofp e in an edge layer of Gaussian shape centred upon each boundary in x, with 1/e width of 1/10 the domain size, with damping coefficient 0.1. Aside from the terms in equation (A.2), no additional diffusive term is added for this profile maintenance.
It is important to note that these dissipation effects are to be understood as part of the numerical scheme, not the equations themselves. The problem being solved is the one with physical dissipation coefficients (such as ρ 2 e ν e as a diffusivity), but in the actual limit in which these are arbitrarily small [24], in the same spirit as high Reynolds number turbulence in hydrodynamics. It would be undesirable to give these coefficients vastly unrealistic values and then take them seriously as processes when in physical fact they are not active at all for k ⊥ ρ s ∼ 1 and k qR ∼ 1 with realistic parameters.
This method confirmed the basic result of [9], namely the geodesic transfer effect as the saturation process for zonal flow self-generation in experimentally relevant parameter regimes for tokamak edge turbulence, but the secondary details have changed, with nonlinear dissipation effects (diffusive mixing) not decisively dominant but concurrent with the processes in the parallel dynamics. For the pressure sidebands therefore, we find that resistive dissipation through the global Alfvén oscillation is an important part of the overall dynamics. The present result is judged to be the more reliable because of the superior properties of the numerical scheme.

Appendix B. Illustration of the Reynolds stress and geodesic transfer processes
This paper shows that the use of control cases is very helpful when diagnosing the various mechanisms entering in such complicated systems as 3D drift Alfvén turbulence-not only for the turbulence itself, as in [38], but also for its interaction with zonal flows. The process by which the Reynolds stress operates can only be elucidated by restricting to the simple 2D model, and switching the adiabatic coupling alternatively in and out of the computation, as in section 6.1. Similarly, identification of geodesic coupling as the process by which zonal flows are saturated is only decisive once the two pieces of the curvature operator (geodesic: acting upon k y = 0 modes, interchange: upon k y = 0 modes) are separated with the double control test as given in section 6.2.
What we examine in this appendix is animations of these control tests which make these processes clear at the intuitive level. We take the standard 2D model from section 6.1 with D = 0 and 0.1, and for each D with ω n = 0 and 1, with all cases run to t = 1000. This will serve to show the role of parallel dynamics in setting up the zonal flows, and the relative lack of a role for the background gradient. We then take the 0 < t < 1000 segments of the standard 3D model from section 6.2 with the parameter set in equation (18), with and without the geodesic coupling terms, watching the differences in the zonal flow morphology.
The first animation (D0w0, same format as in figures 2-4) shows the spatial morphology of the turbulence for D = 0 and ω n = 0. Here, we find as in previous work the emergence of coherent structure in the vorticity [48], and the dynamical alignment between vorticity and passively advected density disturbances [68]. Not only to provide a baseline for the physics, this serves to demonstrate the capability of the numerical scheme in what would have been thought of even in 1984 as a relatively low-resolution computation. The only reason the resultant mean flows are found in the x-direction is the domain aspect ratio of L y /L x = 4.
The second animation (D0w1) shows the result of the gradient by setting D = 0 and ω n = 1. Obviously the flows and vortices do not change, since there is still no back reaction from the density disturbances. The density however shows large scale disturbances being carved out of the gradient, with much of the small-scale manifestation being in the thinness of the sheets at the boundary of the disturbances.
The third animation (D003w0) shows the effect of the adiabatic coupling alone by setting D = 0.03 and ω n = 0. The zonal flows form most rapidly in this case. The fourth animation (D003w1) shows the combined result of the gradient and adiabatic coupling by setting D = 0.03 and ω n = 1. Now, the turbulence is continually forced, with energy flowing into the ExB eddies in the energy containing spectral range [24]. This actually competes with the Reynolds stress mechanism in the ExB energetics-this value of D is the borderline case. With too little adiabatic coupling there is no zonal flow formation and stronger turbulence, while with more coupling (D = 0.1, in section 6.1) there is much more zonal-flow effect. The dominant mechanism in the zonal-flow formation is clearly the adiabatic coupling: tilting of density disturbances, then parallel dynamics forcing the potential disturbances (hence the eddies) to have the same shape as the density disturbances, with the Reynolds stress according to this mechanism always having the sign required to transfer energy from eddies to zonal flows (the zonal vorticity and the Reynolds stress are positively correlated, cf equation (40)).
We now turn to the 3D cases. The fifth animation (nogd) shows the first 1000L ⊥ /c s of the case without the geodesic curvature effect, turning the K 'curvature terms' off for the k y = 0 component of all the variables (same format as in figure 13). With the parallel dynamics in the form of electromagnetic shear Alfvén waves and not merely dissipative coupling, the zonal flows are still found to form prominently.
The sixth animation (nominal) shows the nominal case with all terms having their expected form. Zonal flow layers are now observed to come and go. As found in section 6.2, this is due to the geodesic curvature, the only change between these last two animations. This test confirms that the presence of interchange curvature for both these cases is not decisive. It is important to note that the geodesic oscillation is forced at a broad range of frequencies, most of which are faster than the natural oscillation frequency of √ 2c s /R (corresponding to a period of about 200 in these units). This is why a coherent oscillation is not observed. Also, the slowness of the oscillation means that the turbulence mostly experiences the zonal flows as static, so that it is able to send energy into them. The energy then passes either nonlinearly back to the turbulence or resistively into the global shear Alfvén oscillation related to the Pfirsch-Schlüter current, as diagnosed in section 6.3.
The interplay between the larger eddies and the flows is shown by the 3D structure of the potential in the last two animations. The seventh animation (nominal3d) shows the nominal case, and the eighth (nogd3d) shows the case without the geodesic curvature effect. The latter shows much more dominance by the zonal component, while the nominal case shows the zonal component coming and going as part of the turbulence. Hence, zonal flows act as modulators of turbulence and not as a suppression mechanism [29].