# Understanding nonlinear saturation in zonal-flow-dominated ion temperature gradient turbulence

## Abstract

We propose a quantitative model of ion temperature gradient driven turbulence in toroidal magnetized plasmas. In this model, the turbulence is regulated by zonal flows, i.e. mode saturation occurs by a zonal-flow-mediated energy cascade (“shearing”), and zonal flow amplitude is controlled by nonlinear decay. Our model is tested in detail against numerical simulations to confirm that both its assumptions and predictions are satisfied. Key results include (1) a sensitivity of the nonlinear zonal flow response to the energy content of the linear instability, (2) a persistence of zonal-flow-regulated saturation at high temperature gradients, (3) a physical explanation of the nonlinear saturation process in terms of secondary and tertiary instabilities, and (4) dependence of heat flux in terms of dimensionless parameters.

## 1 Introduction

Ion temperature gradient (ITG) driven turbulence is typically one of the largest causes of thermal losses in magnetized fusion plasmas. It is also an intensely studied type of plasma turbulence, for which a quantitative physical understanding is within grasp. Decades of contentious research seem to converge on a single fact: the interaction between plasma waves and spontaneously generated zonal flows (ZF) is the key problem to understanding the turbulence. This has been called the “drift wave-zonal flow” problem [1]. To clarify the key nonlinear physics, a comparison is often drawn between ITG and electron temperature gradient (ETG) turbulence, which are linearly isomorphic: whereas ETG turbulence displays anisotropic streamers [2], aligned radially with the temperature gradient, the ITG turbulence spectrum is characteristically smeared out among radial Fourier modes, resulting in a much more isotropic energy spectrum (and lower relative transport) – the reason for this difference is ZFs.

It can be argued that the difference between ITG and ETG turbulence is even more fundamental: while an ETG mode is thought to saturate by energy transfer directly into secondary modes [2], the energetic cost imposed on an ITG mode in the generation and sustainment of ZFs is small, and the real impact of the ZF is in its conservative “shearing” of the ITG mode; see [3]. The term “regulator” is thus applied to ZFs, evoking the analogy of a special device (e.g. a valve) that limits the operation of a machine.

Significant progress has been made recently in the understanding of ITG turbulence. A cascade model was proposed by [4] that eschews details of ITG mode saturation and ZF physics in favor of a simple set of scaling conjectures. Although excellent agreement with numerical simulations was found, the theoretical explanation seems incongruous with both previous [5] and subsequent [6] work, both of which support the paradigm of ZF regulation. In particular, it was demonstrated [6] that nonlinear energy transfer between ITG modes in the saturated steady state is dominated by wavenumber triads that are zonal-flow mediated (see Figure 11(b) of [6]); in other words, the energy transfer between ITG modes is accounted for mostly by ZF shearing. Adding further detail to the picture of ZF regulation, it was subsequently discovered [3] that the amplitude of ZFs in the turbulent state is controlled by energetic properties of the linear eigenmodes that drive the turbulence. In particular, it was shown that the generation of ZFs by the inverse cascade of electrostatic energy (i.e. its flow to large scales via nonlinear interaction) can be weakened or even reversed if the unstable modes possess sufficiently large free energy, a condition that is approached, e.g., at large temperature gradients (see Sec. 5).

Despite the existence of a large body of research, the challenge remains to understand in a detailed and quantitative manner the processes by which the ITG turbulent state forms. The present work addresses this challenge, focussing on the strongly nonlinear regime, well above marginal stability, where ZFs saturate by nonlinear decay (rather than by collisions or other linear damping mechanisms). We propose a phenomenological cascade model, based on the conservative nonlinear transfer of energy to small scales, where it is ultimately dissipated. The model is developed and tested as follows. By performing electrostatic gyrokinetic simulations (employing the GENE code) with parameters similar to the Cyclone Base Case (CBC), we demonstrate (1) that the rate of energy injection in ITG turbulence is determined by the linear growth rate and (2) a properly defined ZF shearing rate balances with this rate at energy-containing scales. We thus argue that saturation occurs by the nonlinear generation of ZFs and concurrent shearing of ITG modes. The key new ingredient in this picture is that the ITG mode intensity needed to produce sufficiently strong ZFs is sensitive to the energy content of the ITG modes; we explain this physically using secondary and tertiary instability theory. We propose a quantitative model that incorporates these ingredients and show that it agrees with the observed heat flux in the strongly-driven limit.

## 2 Equations and definitions

We solve the electrostatic gyrokinetic equation for ions, which can be written (ignoring collisions) as

(1) |

where the coupling coefficient is , with , where is the flux surface label, and is the field line label. The ion gyrocenter distribution is , with denoting the coordinate along the magnetic field line, the wavenumber perpendicular to the magnetic field, an and the magnetic moment and energy respectively. The frequencies are defined , and , with , , , , and , and the ion charge in terms of the absolute electron charge . The equilibrium ion distribution function is , , where and are the bulk ion density and temperature. For the electron species we assume a modified Boltzmann response, so the quasi-neutrality constraint determining is

(2) |

where , and is the electron temperature. The non-zonal part of the potential is defined

(3) |

and the zonal part is defined by the usual flux surface average

(4) |

In the final line we have introduced the field line average . We will refer to and respectively as the zonal and non-zonal parts of . In our theoretical arguments, we invoke the strongly ballooning limit (see A), whereby the dynamics are nearly two-dimensional (i.e. do not vary significantly along ). In this case we may use the local approximation , and the local notation , where , , , , , , and is the radius of curvature for the magnetic field. Note that henceforth we will find it convenient to use (, ) in place of (, ) to conceal the presence of magnetic geometry and generally simplify the presentation.

Numerical simulations are done with flux tube domains, which are three-dimensional (non-uniform in ). These simulations are characterized by dimensionless parameters , , , and , where we have introduced a single scale ^{1}

To define a turbulent cascade, we require an energy quantity, which, following convention, we take to be the free energy

(5) |

where is the system volume. This quantity is conserved by all terms in the collisionless gyrokinetic equation (1) except the source term proportional to .

## 3 Turbulence phenomenology

In order to describe ITG turbulence in terms of an energy cascade, we will need to use a “phenomenological” description, by which we mean a shorthand notation for describing the cascade process, i.e. the injection and nonlinear turnover of energy, etc. Although it may be useful to use a loosely-defined phenomenology to develop initial understanding, it is clearly preferable to use precisely defined quantities, as we will do here. This sharpens the understanding and makes the theory more directly testable. For fluid turbulence, it was suggested in Chapter 6 of [8] that the velocity field at scale , denoted , should be defined as the root-mean-square-value of the field filtered around the wavenumber . That is, writing , and denoting the ensemble average (informally a time average) as , we can relate to the RMS Fourier component as follows:

(6) |

where , is the dimension of the space in which the cascade occurs, and means . Going between the first and second line we have assumed the turbulence is statistically homogeneous in space (Wiener–Khinchin theorem). Note that the first line can be considered as a definition of , but since this is a phenomenological quantity, it is only defined up to an overall multiplicative constant of order unity. Thus the factor of 2 in the range of is somewhat arbitrary. In deriving the final line we have assumed that the spectrum is such (e.g. a power law) that the effect of integrating over simply introduces a factor of . The final line demonstrates that the number of modes that participate in local energy cascade is effectively .

For present purposes, the electrostatic potential is the turbulent field, and we will need to divide it into zonal and non-zonal contributions, i.e. where recall is the flux-surface average of . ITG turbulence is not isotropic, but, if one excludes the zonal component, it does have a uniform distribution in for . Thus we define , the non-zonal electrostatic potential at scale , as

(7) |

where we have averaged over , and excluded the zonal component, i.e.

(8) |

where is the discrete delta function. Next, we define the zonal potential at scale as

(9) |

where

(10) |

Note that the dimensional factor reflects the fact that shearing by ZFs is a one-dimensional process, i.e. the flow of energy is only in the direction.

If we now assume that the energy of turbulence peaks at a characteristic scale , then the total energy density obtained by integrating the spectrum can be expressed in terms of the fields at that scale. That is, following Equation (6), we can write , where denotes the system volume. In a similar manner, the heat flux is expressed phenomenologically as

(11) |

To obtain the final expression, we have assumed either that (a) the phase of and are uncorrelated, or (b) that and are not systematically in-phase to the degree that a significant cancellation occurs.

## 4 Energy injection and nonlinear transfer

Our model assumes that the turbulence has two key properties: (P1) the nonlinear turnover rate is determined by a local (in wavenumber space) ZF-mediated transfer of free energy and (P2) the injection of energy is largely determined by the linear instability.^{2}

Because large scales dominate the transport spectrum of ITG turbulence, let us consider the long-wavelength fluid limit of gyrokinetics, . At these scales, nonlinear phase-mixing [12] and other finite-Larmor-radius (FLR) effects are weak, and the density moment of the gyrokinetic equation (1) yields (see [3], Eqn. D.6)

(12) |

Here we recall that is the non-zonal part of (see Equation (3)), and represents linear source terms; the equation for the evolution of the zonal part is not shown here. The form of Equation (12) reflects the nonlocal response of the electrons, which travel rapidly along the magnetic field. It is well-known that this modification to the traditional Boltzmann electron response causes strong ZF generation. As Equation (12) reveals, another consequence of the electron response is that, for , the nonlinear evolution of the non-zonal potential is due entirely to conservative “shearing” by ZFs. Note however that the linear coupling () to moments other than density (e.g. parallel ion flow and perpendicular ion temperature fluctuations) underlies the mechanism of linear instability, and so one cannot know a priori whether the nonlinear interactions involving the non-zonal component of such moments might lead to mode saturation instead of the ZF shearing mechanism. Evidence from the literature helps rule out the possibility that these so-called “drift-wave/drift-wave” (DW/DW) interactions might be important: First, when ITG or ETG modes saturate without ZFs, the process generates fine structure parallel to the magnetic field, and a strongly anisotropic turbulence spectrum perpendicular to the magnetic field, as a consequence of the Cowley secondary instability [2, 13]. This, however, is not typical in ITG turbulence with CBC-like parameters, on which we focus here; an exception is found at low as discussed later. Second, and on a more empirical note, these DW/DW interactions have been investigated in numerical simulations and found not to contribute significantly to nonlinear transfer [5, 6], as previously discussed.

In short, we argue that Equation (12) captures the dominant nonlinear energy transfer mechanism. This mechanism can be thought of as a one-dimensional cascades in -space, i.e. energy is transported in -space under the action of zonal shearing, but not transferred significantly in -space (P1).

Let us describe this process phenomenologically. We first assume locality, i.e. that the nonlinear interactions tend to involve wavenumbers of a similar magnitude. Locality in gyrokinetic turbulence is well-established [14, 15, 16]. The turbulence can then be described in terms of the quantities and , as discussed in Section 3, bearing in mind their relationship to the RMS potential , henceforth abbreviated as simply .

Now we let us proceed to defining the nonlinear turnover rate . Equation (12) describes a one-dimensional shearing process. We model this as an energy cascade, where the nonlinear turnover (at the rate ) is due to the zonal velocity acting on non-zonal fluctuations. Using Equation (9) we can obtain the nonlinear turnover rate at perpendicular wavenumber (, ), i.e.

(13) |

where

(14) |

where is defined in Equation (10), is the wavenumber spacing, is the wavenumber of the ZF, and correspond to the non-zonal mode that is being sheared. The meaning of Equation (14) is that the energy in mode turns over at the rate due to shearing by ZFs having wavenumbers around . However, since the most unstable ITG modes satisfy , the further assumption of locality () is sufficient to express the required nonlinear turnover rate of such modes in terms of the single wavenumber :

(15) |

To compare with this rate, we define the linear growth at wavenumber as

(16) |

Let us test how well the assumptions of the above model are satisfied by numerical simulations of the turbulence; Figure 1 summarizes our findings. For these plots the wavenumber is normalized to the peak of the heat transport spectrum (approximately the same as the peak of the free energy spectrum), denoted or alternatively . Note that if the nonlinear turnover of the turbulence increases sufficiently fast in , then the wavenumber will be set by the low- cutoff [17] of the (dominant) instability, rather than its peak; indeed is sometimes determined this way in quasilinear transport estimates [18]. The bottom panel of Figure 1 compares the maximum linear growth rate for ITG modes (i.e. for ) with the nonlinear turnover rate . The top panel compares the same growth rate with the rate of injection of free energy in the fully developed turbulence, . This quantity is formally defined by splitting the energy budget equation into linear and nonlinear contributions at each (see [19] for details): , where , and subscripts ‘L’ and ‘NL’ denote contributions from linear and nonlinear terms respectively. Then is

(17) |

The frequency is the correct quantity to compare with the nonlinear frequency since nonlinear turnover refers explicitly to the transfer of energy in -space. Our propositions (P1) and (P2) can now be succinctly stated as

(18) |

with values given by Equation (15), Equation (16) and Equation (17). However, these are not equalities, so their validity is tested not by how close the frequency ratios are to unity, but by whether or not they systematically change with the variation of parameters. Figure 1 demonstrates that the balance is satisfied uniformly well over a significant range of . However, there is a significant variation in ; this may reflect a different cascade process at small scales, as we will discuss later. The balance , is also satisfied uniformly well over .^{3}

## 5 Saturation of the turbulence

Having validated the assumptions of our model, let us now complete our description of the saturation process. Using Equation (15) and Equation (16), the balance implies a saturation rule for the zonal potential:

(19) |

where was defined in Equation (10). Note that this saturation rule differs from that of [5] in the mode-counting factor , and differs from that of [20] in that it does not assume timescale separation.

Equipped with Equation (19), we must now determine the relative amplitude of non-zonal fluctuations to complete the description of the saturated state. One might naively guess that the nonlinear rate associated with DW energy turnover, , should balance with the nonlinear rate of the zonal potential, (this follows from the vorticity equation for , i.e. the zonal part of Eqn. (D.18) of [3]), which would lead to the relative saturation amplitude . This turns out to be incorrect: it is missing a constant of proportionality, i.e. we should instead write , or in -space (see Equation (7) and Equation (9)) we have

(20) |

Equation (19) and Equation (20) raise a subtle conceptual issue: A zonal Fourier component only scales with the system size in the dimension, . This is in contrast to a Fourier component in isotropic turbulence, which must scale with the full system volume, i.e. also with and , because the number of modes per unit volume in -space increases with and so each individual Fourier mode becomes less important. Consequently, for isotropic turbulence, in the limit that any one of the dimensions of the system is large, any collection of Fourier components lying on a surface in -space must occupy a vanishing fraction of the -space volume and could be neglected without influencing the system dynamics. ZFs, however, occupy a single line in -space, and therefore if ZFs are to act as the regulator of the turbulence they must always be important and cannot diminish with increasing or . For this reason, one should not directly compare zonal and non-zonal Fourier modes as their relative amplitudes depends on system size, as specified in Equation (20).

Returning to the meaning of Equation (20), the factor quantifies a certain “non-universality”: in addition the scale () and rate (), it turns out that the saturation process also depends on an additional energetic property of the ITG mode itself, which causes the nonlinear response of the ZFs to weaken as increases.^{4}

This effect was described by [3] as a cascade reversal, controlled by an energy ratio parameter , where for our purposes where . To obtain a quantitative estimate for , we will take a slightly different approach via the primary/secondary/tertiary instability framework conceived by [13, 21]. This approach allows one to probe the nonlinear physics by linearizing the dynamics about a nonlinearly-motivated initial condition. The secondary instability is calculated by linearizing the gyrokinetic equation about a large-amplitude “primary” ITG mode, and the tertiary mode is calculated by doing the same about a large-amplitude zonal mode. We argue that the relative amplitude of zonal and non-zonal fluctuations is set by balancing secondary and tertiary growth rates, i.e. and . The ratio of zonal and non-zonal amplitudes is set according to this balance, because if it deviates significantly, then one instability will overwhelm the other and guide the system back into a state of balance.

Although the secondary instability is usually expressed in terms of an isolated Fourier amplitude, we use the quantity to reflect the turbulent state. Thus, we take

(21) |

where is the wavenumber of the secondary mode and is the wavenumber of the primary mode; see Appendix D.5 of [3]. Defining we take

(22) |

where is the wavenumber of the tertiary [21].^{5}

(23) |

This expression depends on an unknown, the wavenumber of the tertiary. It was found by [21] that gives a maximal growth rate. However, there is a large range of scales accessible to the tertiary mode and the mode of peak growth rate is not necessarily the most important, so it is not clear how should be chosen.^{6}

(24) |

## 6 Strongly driven limit

Now we can consider the strongly driven (large ) limit, for which we can evaluate several quantities according to linear theory:

(25) | |||

(26) | |||

(27) |

The estimates 25 and 26 are determined by the strongly ballooning toroidal ITG mode (see i.e. [22] and A). The expression 27 can be obtained by balancing the sound wave transit rate with the growth rate 25, i.e. ; note that this balance is similar to what [4] call “critical balance,” and both yield the proportionality , though they disagree on the scaling in . The wavenumber corresponds to the transition from the toroidal branch to the subdominant slab branch, which is clearly observed in linear simulations. Naively, one might expect the linear estimate for , i.e. Equation (26), to be satisfied for the zonal component. However, although the -dependence is satisfied, the -dependence is not, at least for the range of parameters that we have investigated; instead we observe , which using Equation (24) yields

(28) |

The scalings above are supported by Figure 2. For the plotted quantity , we expect from Equation (20), Equation (27), and Equation (28) that . For , Equation (26) implies .

Now, using Equations 19, 20, 25, 26, 27, and 28, the heat flux may be calculated as (see Section 3), yielding

(29) |

which differs from [4] in its dependence; this difference is significant because it reflects completely different saturation physics than that assumed by [4]. Our simulations obey the scaling in Equation (29), as shown in Figure 3. Note that although nothing fundamental prohibits applying this to , the ITG mode becomes stabilized for the values of that we used. Note also that deviation from the theoretical line is expected at sufficiently low- because, given fixed , the weakening of the ZFs must ultimately cause nonlinear interactions among the ITG modes (DWs) to become important. Indeed, although all our other simulations show the signature of dominant zonal shearing, i.e. a flat energy distribution for , the spectrum of the low- outliers exhibit peaking at low- and associated streamer-like structures.

An important ingredient in the calculation of the overall saturation and heat flux is the wavenumber , which corresponds to the outer scale of the turbulence. As a measure for , the quantity is plotted in Figure 4 and compared with the theoretical prediction of Equation (27). This prediction differs from the outer-scale estimate of [4], that find , where . We suspect that these two definitions give different results because the spectrum of electrostatic energy tends to a constant at low and so the factor of preferentially weights large scales in the sum , i.e. the lowest available wavenumber dominates the computation of . Thus we conclude that is the correct measure of the dominant wavenumber contributing to heat flux.

## 7 Spectra

It is noteworthy that detailed knowledge of the energy spectrum is not important for determining the scaling of bulk properties of the turbulence (, , etc.) with respect to the system parameters. This is because the integral of a power-law spectrum will only introduce an overall dimensionless factor, corresponding to the steepness of the spectrum. Nevertheless, it is interesting that some universal behavior is exhibited, at least in the regime considered here.

Formally, we define the spectrum of zonal and non-zonal energies as and , respectively. The spectrum of zonal energy is predicted immediately by the saturation rule Equation (19), and the degree to which this is confirmed is established by the second panel of Figure 1. The dependence then implies that for at least some range.

The non-zonal spectrum at depends on the saturation process at those wavenumbers. As we have noted, the balance between injection, nonlinear turnover by zonal shearing, and linear growth is not uniformly satisfied in . Furthermore, there is evidence that zonal shearing only provides a small part of energy transfer at inertial-range (small) scales [23]. It is thus possible that an isotropic cascade (in ) overtakes the zonal-mediated transfer at these scales. However we argue that an isotropic cascade is not possible around because the associated nonlinear turnover rate is too weak to account for the observed energy injection; if, for example, this rate was computed according to the FLR corrections of Equation (12), one would obtain for . The spectral scaling [4] seems consistent with our numerical results in the large limit, although low- clearly has a steeper spectrum.

## 8 Discussion

We have presented a model of the saturation process of ITG turbulence, which agrees with numerical simulations in a detailed and quantitative way. We have made this description as explicit and physically transparent as possible, and a key part of our theoretical argument is a set of linearly calculable (primary, secondary, and tertiary) instabilities. We hope that these features will make it possible to generalize our results in the future.

Based on our findings, the physical picture of ITG turbulence is as follows. The saturation of the turbulent state is determined by a cascade process at an outer scale. The injection of free energy into the turbulence is determined by the linear growth rate of the ITG mode, and the flux of this energy is carried to high- by the shearing action of ZFs, which we have described as a one-dimensional cascade. The amplitude of the DWs relative to the ZFs is sensitive to the amount of free energy (relative to the electrostatic energy) that is entering the cascade, which is a feature of the ITG mode itself. This energy ratio is a fundamentally kinetic quantity, as it measures the amount of excitation of velocity-space structure in the distribution function. This property of the ITG mode can be quantified via the ratio of the temperature and potential fluctuations, and combined with the saturation rules, and outer scale estimate, yields the observed scaling of the ion heat flux .

Generally, our findings support the paradigm of ZF regulation, and we have shown that it applies even in the strongly-driven regime. Interestingly, we have seen that the ZFs are capable of regulating the turbulence, even as their amplitude (relative to DWs) is diminished. This fact should serve as a warning against using the relative amplitude of the ZFs as a measure of their strength. Instead the dynamical importance of ZFs is properly evaluated by comparing the rate of zonal shearing to the rate of energy injection by linear instability.

One issue that merits further investigation is the nonlinear decay process of the ZFs. We have cited the tertiary instability of [21] and given an argument for how it balances with ZF generation process (secondary instability), but a detailed description that accounts for a range of interacting scales, especially near marginal stability of the ITG mode, remains an open problem.

As a final point, we should note that our work is clearly only valid when the nonlinear decay of ZFs dominates over collisional decay. This is the regime of strong turbulence, corresponding to sufficiently large temperature gradient, i.e. those above the nonlinear critical gradient of Dimits [24]. It is not clear what regime is most relevant for present experiments, but future fusion devices with lower collisionality will exhibit lower transport in the weakly driven regime (in which ZF decay is collisional) and so the strong regime should be more accessible. Furthermore, if collisional decay and nonlinear decay must both be included to model a given experiment, it should prove useful to have a good understanding of the fully nonlinear regime.

## 9 Acknowledgements

We gratefully acknowledge Per Helander for helpful comments, and the Wolfgang Pauli Institute for hosting a series of workshops on Gyrokinetics. The research leading to these results has received funding from the European Research Council under the European UnionÕs Seventh Framework Programme (FP7/2007-2013)/ERC Grant Agreement No. 277870.

## Appendix A Properties of the strongly-driven toroidal ITG mode

The term “strongly ballooning” describes an asymptotic limit of the full three-dimensional ballooning mode calculation, whereby the two conditions are simultaneously satisfied: (1) the mode is strongly localized at the outboard midplane (location of bad curvature) and (2) the transit frequency (i.e. the rate associated with ion streaming a distance equal to the scale of the mode parallel to the magnetic field) is small. Because these two conditions may appear contradictory, it is useful to show explicitly that they can be satisfied simultaneously. A demonstration of this has been given in a recent work [22].

The strongly ballooning limit is essentially the strongly driven limit, which is also a non-resonant limit. The following local dispersion relation for the toroidal branch of the ITG mode is valid ():

(30) |

where . The distribution is

(31) |

where and . The strongly driven limit can be explicitly written as . Taking also , the solution can be written . To calculate the perturbed pressure (as it appears in the heat transport flux) we need only retain the dominant contribution, i.e. , so that

(32) | |||||

(33) |

Note that the temperature is out of phase with the potential, which means it contributes fully to transport.

## References

### Footnotes

- For a generic equilibrium magnetic field, arbitrarily many quantities may be introduced to parameterize its spatial variation; however, by limiting ourselves to the CBC model geometry, it is sufficient to use the single parameter , where is the safety factor.
- It has been noted by several authors that gyrokinetic turbulence remains “nearly linear,” as reflected by signatures of the linear eigenmodes in the fully saturated state of nonlinear simulations [9, 10, 11].
- In fact, in the purely linear system would be satisfied exactly (after a transient period) because the growth in energy would be due to the most unstable modes at each . Thus, the small value of the ratio at high signals nonlinear behavior, and may be due the nonlinear excitation of the damped spectrum of modes [11].
- Note that as increases, the ZFs weaken while the turbulence simultaneously goes to larger scales, causing the nonlinear terms due to the FLR (neglected in deriving Equation (12)) to also weaken. Thus, zonal shearing can remain the dominant nonlinear transfer mechanism.
- Note that our calculation neglects the effects that stabilize the tertiary at low- (such as finite ), so its validity is limited to values of significantly larger than the nonlinear critical gradient. However it seems generalizable to the weakly supercritical regime if more complete tertiary mode physics is included.
- Consider also that the nonlinear transfer of electrostatic energy shown in Figure 2(b) of [19] is nonlocal and involves a broad range of scales. This hints that the nonlinear zonal decay process may involve a range of scales.

### References

- P H Diamond, S-I Itoh, K Itoh, and T S Hahm. Zonal flows in plasmaâa review. Plasma Phys. Control. Fusion, 47(5):R35, 2005.
- F. Jenko, W. Dorland, M. Kotschenreuther, and B. N. Rogers. Electron temperature gradient driven turbulence. Physics of Plasmas (1994-present), 7(5):1904–1910, 2000.
- G G Plunk, T Tatsuno, and W Dorland. Considering fluctuation energy as a measure of gyrokinetic turbulence. New J. Phys., 14(10):103030, 2012.
- M. Barnes, F. I. Parra, and A. A. Schekochihin. Critically balanced ion temperature gradient turbulence in fusion plasmas. Phys. Rev. Lett., 107:115003, Sep 2011.
- R. E. Waltz and C. Holland. Numerical experiments on the drift wave–zonal flow paradigm for nonlinear saturation. Phys. Plasmas, 15(12):122503, 2008.
- M. Nakata, T.-H. Watanabe, and H. Sugama. Nonlinear entropy transfer via zonal flows in gyrokinetic plasma turbulence. Phys. Plasmas, 19(2):022303, 2012.
- W. K. Hagan and E. A. Frieman. Nonlinear gyrokinetic theory, the direct interaction approximation, and anomalous thermal transport in tokamaks. Phys. Fluids, 29(11):3635–3638, 1986.
- U. Frisch. Turbulence: The Legacy of A. N. Kolmogorov. Cambridge University Press, 1995.
- Tilman Dannert and Frank Jenko. Gyrokinetic simulation of collisionless trapped-electron mode turbulence. Phys. Plasmas, 12(7):–, 2005.
- T. Görler and F. Jenko. Multiscale features of density and frequency spectra from nonlinear gyrokinetics. Phys. Plasmas, 15(10):–, 2008.
- D. R. Hatch, P. W. Terry, F. Jenko, F. Merz, and W. M. Nevins. Saturation of gyrokinetic turbulence through damped eigenmodes. Phys. Rev. Lett., 106:115003, Mar 2011.
- W. Dorland and G. W. Hammett. Gyrofluid turbulence models with kinetic effects. Phys. Fluids B, 5(3):812–835, 1993.
- S. C. Cowley, R. M. Kulsrud, and R. Sudan. Considerations of ion-temperature-gradient-driven turbulence. Phys. Fluids B, 3(10):2767, 1991.
- T. Tatsuno, M. Barnes, S. C. Cowley, W. Dorland, G. G. Howes, R. Numata, G. G. Plunk, and A. A. Schekochihin. Gyrokinetic simulation of entropy cascade in two-dimensional electrostatic turbulence. J. Plasma Fusion Res. Ser., 9:509, 2010.
- A. Bañón Navarro, P. Morel, M. Albrecht-Marc, D. Carati, F. Merz, T. Görler, and F. Jenko. Free energy cascade in gyrokinetic turbulence. Phys. Rev. Lett., 106:055001, Jan 2011.
- Bogdan Teaca, Alejandro Bañón Navarro, Frank Jenko, Stephan Brunner, and Laurent Villard. Locality and universality in gyrokinetic turbulence. Phys. Rev. Lett., 109:235003, Dec 2012.
- F. Jenko, W. Dorland, B. Scott, and D. Strintzi. Simulation and theory of temperature gradient driven turbulence. In O. Sauter J.W. Connor and E. Sindoni, editors, Theory of Fusion Plasmas, page 157. Societa Italiana di Fisica, Bologna, 2002.
- J. Citrin, C. Bourdelle, P. Cottier, D. F. Escande, Ö. D. Gürcan, D. R. Hatch, G. M. D. Hogeweij, F. Jenko, and M. J. Pueschel. Quasilinear transport modelling at low magnetic shear. Physics of Plasmas, 19(6):062305, 2012.
- A. Bañón Navarro, P. Morel, M. Albrecht-Marc, D. Carati, F. Merz, T. Görler, and F. Jenko. Free energy balance in gyrokinetic turbulence. Phys. Plasmas, 18(9):092303, 2011.
- T. S. Hahm, M. A. Beer, Z. Lin, G. W. Hammett, W. W. Lee, and W. M. Tang. Shearing rate of time-dependent e x b flow. Physics of Plasmas, 6(3):922–926, 1999.
- B. N. Rogers, W. Dorland, and M. Kotschenreuther. Generation and stability of zonal flows in ion-temperature-gradient mode turbulence. Phys. Rev. Lett., 85(25):5336–5339, Dec 2000.
- G. G. Plunk, P. Helander, P. Xanthopoulos, and J. W. Connor. Collisionless microinstabilities in stellarators. iii. the ion-temperature-gradient mode. Phys. Plasmas, 21(3):–, 2014.
- Bogdan Teaca, Alejandro Bañón Navarro, and Frank Jenko. The energetic coupling of scales in gyrokinetic plasma turbulence. Phys. Plasmas, 21(7):–, 2014.
- A. M. Dimits, G. Bateman, M. A. Beer, B. I. Cohen, W. Dorland, G. W. Hammett, C. Kim, J. E. Kinsey, M. Kotschenreuther, A. H. Kritz, L. L. Lao, J. Mandrekas, W. M. Nevins, S. E. Parker, A. J. Redd, D. E. Shumaker, R. Sydora, and J. Weiland. Comparisons and physics basis of tokamak transport models and turbulence simulations. Phys. Plasmas, 7(3):969, 2000.