Enstrophy non-conservation and the forward cascade of energy in two-dimensional electrostatic magnetized plasma turbulence

A fluid system is derived to describe electrostatic magnetized plasma turbulence at scales somewhat larger than the Larmor radius of a given species. It is related to the Hasegawa- Mima equation, but does not conserve enstrophy, and, as a result, exhibits a forward cascade of energy, to small scales. The inertial-range energy spectrum is argued to be shallower than a -11/3 power law, as compared to the -5 law of the Hasegawa-Mima enstrophy cascade. This property, confirmed here by direct numerical simulations of the fluid system, may help explain the fluctuation spectrum observed in gyrokinetic simulations of streamer-dominated electron-temperature-gradient driven turbulence [Plunk et al., 2019], and also possibly some cases of ion-temperature-gradient driven turbulence where zonal flows are suppressed [Plunk et al., 2017].


Introduction
The turbulent cascade, a mechanism for the nonlinear transfer of energy across scales, is a key idea for understanding kinetic magnetized plasma turbulence. By considering simplified models, in uniform magnetic geometry, one can obtain a theoretical prediction for the spectrum of fluctuations, valid across an "inertial range" of scales, free from energy sources and sinks. Though such a theory is not able to fully describe the behavior of realistic turbulence, which hosts instabilities, damped modes, complicated magnetic geometries, etc., it nevertheless constitutes a quantitative prediction of nonlinear behavior of the underlying gyrokinetic equation, an equation which generally governs actual systems of practical interest. The existence of such theoretical test cases is valuable for validating the solution methods employed by gyrokinetic codes, and as a foundation for physical interpretation of the volumes of data they produce.
Here, a novel quasi-two-dimensional fluid system is derived from the electrostatic gyrokinetic system, to describe fluctuations that predominantly vary in the directions perpendicular to the mean magnetic field, i.e. in the "drift plane", at scales larger than the Larmor radius ρ, corresponding to a species of interest. The notion that quasi-two-dimensional behavior might underly magnetized plasma turbulence is intuitively justified by the fact that the magnetic guide field renders the dynamics inherently anisotropic. Furthermore, instabilities that drive electrostatic turbulence in fusion plasmas, e.g. the ion-temperature-gradient (ITG) and electron-temperature-gradient (ETG) modes, exhibit a kind of localization along the field line, accompanied by the domination of perpendicular dynamics over parallel dynamics. The fluid limit ρ is of particular importance, because the energy of the fluctuations is predominantly found at such scales -these are the scales of importance, most directly affecting the performance of fusion devices. Furthermore, the reduction of complexity afforded by fluid limits can reveal important features of the dynamics, that do not manifest in the analysis of the general gyrokinetic equations. Although similar systems as the one presented here have been proposed and studied in the past, most notably the Hasegawa-Mima (HM) equation, Hasegawa and Mima [1978], the present derivation takes special care in considering the consequences of the appearance of nonlinear finite-Larmor-radius (FLR) terms that appear in the dynamical equation for the electrostatic potential -i.e. the "vorticity" equation. Such terms introduce a closure problem in the fluid moment hierarchy, where lower moments are coupled to ever higher ones, generally without end. This motivates the cold ion limit that underlies the HM equation, which eliminates the inconvenient terms, but is however not generally appropriate for application to fusion plasmas. In the present work, it is noted that the presence of these terms introduces rapid dynamics, and a multiple-scales analysis is proposed in which the fluid moment hierarchy closes at the pressure moment, without using an ad hoc closure scheme, leading to a relatively simple system involving only two fields.
What is immediately apparent is that the presence of the additional field (the pressure perturbation) breaks the nonlinear conservation of enstrophy that is famously satisfied by the HM equation, and there causes an "inverse cascade" of energy to large scales. The new system, we argue, should exhibit distinct nonlinear behavior, including a shallower energy spectrum when the effect of the nonlinear FLR terms is sufficiently strong. Direct numerical simulation of the fluid model gives some confidence in these predictions. The results of this work may help to interpret observations of turbulence in parameter regimes where the dynamics tend toward the quasi-two-dimensional limit. We discuss possible examples, including cases explored in previous gyrokinetic turbulence simulations in tokamak and stellarator geometries.

Equations and definitions
We assume uniform magnetic geometry, where the magnetic guide field is constant and points in theẑ-direction. One species is assumed to be kinetic, with the other species satisfying a simple Boltzmann response model. We begin with a nondimensional form of the gyrokinetic system [Plunk et al., 2010], normalized relative to the kinetic species: v ⊥ /v th → v (with v th = » T /m, and T and m are the temperature and mass of the kinetic species) is the normalized perpendicular velocity and the normalized wavenumber is k ⊥ ρ → k where thermal Larmor radius of the kinetic species is ρ = v th /Ω c and Ω c = qB/m. The two-dimensional gyrokinetic equation is written as follows in terms of the perturbed gyrocenter distribution function g(R, v, t), where R =xX +ŷY is the gyrocenter position: where C[h] R is the collision operator (not treated here explicitly); the Poisson bracket is 2π 0 dϑA(R + ρ(ϑ)), where the Larmor radius vector is ρ(ϑ) =ẑ × v = v ⊥ (ŷ cos ϑ −x sin ϑ) and ϑ is the gyro-angle. (Note that the quantity inside of the collision operator is h = g + ϕ R F 0 . Note also that the spatial coordinate is R in the gyrokinetic equation and, formally, the spatial derivatives are to be interpreted in this variable, but for simplicity we avoid making the distinction explicit.) We mostly ignore the collision operator but note that some mechanism of coarse-graining will be necessary to get sensible solutions out of the equation. Quasi-neutrality yields the electrostatic potential ϕ(r, t), where r =xx +ŷy is the position-space coordinate: where the g is implicitly assumed to be integrated over parallel velocity so that 2π ∞ 0 vdv completes the integration over three-dimensional velocity space. The angle average is defined A(R) r = 1 2π 2π 0 dϑA(r − ρ(ϑ)), and the term τ ϕ is the adiabatic density response, and τ = T i /(ZT e ) for the case of ion scales and τ = ZT e /T i for the case of electron scales. For the ion case, this Boltzmann response might be considered reasonable if zonal flows are strongly suppressed. The operator Γ 0 φ = 2π ∞ 0 vdv F 0 (v) φ R r is more naturally expressed in Fourier space, assuming a Maxwellian background where I 0 is the zeroth-order modified Bessel function.

Fluid limit
We expand in the limit i.e. we assume that the scales of interest are larger than the Larmor radius of the species of interest. For electron scales, the limit is considered subsidiary to the adiabatic ion limit, so scales of the turbulence must remain much smaller than the ion Larmor scale, i.e., ρ e /ρ i k 1. Note that there may also be a minimum applicable k imposed by dynamics parallel to the magnetic field, but treating this explicitly is outside the scope of this work. We will only need the first two orders of the expansion in δ. The gyrokinetic equation, henceforth omitting explicit collisional effects, is written as and Eqn. 2 becomes We will denote v-moments of g as

Naive expansion
To give a taste for the issues that arise in the expansion, let us take an initial informal look at the moments of the gyrokinetic equation. We first examine the density moment. We include only the ostensibly dominant nonlinear terms. We stress that this equation is given only for illustrative purposes, and is not to be taken as a basis for the later derivations of the paper: Note that the term −∂ t ∇ 2 ϕ, which appears in the HM equation, should be neglected here because it is formally smaller than ∂ t ϕ by one power of the ordering parameter δ. Likewise, the term −∂ t ∇ 2 G 2 must be considered negligible if the ordering G 2 ∼ ϕ and ∂ t G 2 ∼ ∂ t ϕ holds. The situation is, however, a bit more subtle. The above equation couples to the v 2 moment of g, G 2 , and the equation for this and other such moments can be written, neglecting higher-order FLR terms, as What we now notice, examining these two equations, is that the density equation is driven by nonlinear terms that appear to be much smaller than those controlling the higher moments of the distribution function -that is, the equations for G n have dominant contributions from the E ×B nonlinearity, while the potential evolves under the influence of terms like the "polarization drift" nonlinearity, which is smaller by a factor of δ = k 2 . One possible resolution of this apparent imbalance is to consider G n itself to be large, as for instance in the non-resonant limit of the ITG/ETG mode, i.e. G n ∼ δ −1 ϕ. In this case, the polarization drift nonlinearity can be neglected, and we see the justification for retaining the additional time derivative term above, since ∂ t ϕ ∼ ∂ t ∇ 2 G 2 . This term can be evaluated from the Laplacian of Eqn. 8, yielding Eqns. 8-9 demonstrate a consistent fluid limit, but cannot describe ITG or ETG turbulence in the resonant limit, where ϕ ∼ G 2 . This corresponds the more physically reasonable scenario of a modest turbulence drive -i.e. not very far from the linear critical gradient, or considering the weakly unstable, large-scale modes that dominate the turbulence spectrum. To treat this limit properly, we must account for the fact that ϕ evolves much more slowly than G n . Physically, it can be argued that, in a turbulent state, Eqn. 8 will then describe rapid mixing of G n by E × B vortices, so that any initial variation along streamlines of constant ϕ will decay on a fast timescale (with the help of some explicit dissipation), leaving G n to be constant along those streamlines [Cowley, 2008]. To account for such processes more carefully, we abandon conventional perturbation theory in favor of the method of multiple scales (see, i.e. Bender and Orszag [1978]). We will henceforth disregard the equations presented here, in section 3.1, and proceed to derive equations that contain only terms justified by a set of explicitly stated ordering assumptions.

Multiscale expansion
We introduce the fast and slow time variables t f , and t s , such that ∂ t f ∼ {ϕ, .} and ∂ ts ∼ {∇ 2 ϕ, .} ∼ δ∂ t f and expand the fields as We reiterate that the assumptions we have made are δ 1, the above multi-scale expansion, and the quasi-two-dimensional approximation, whereby variation in the direction along the magnetic field is neglected, and the non-kinetic species is assumed to follow a Boltzmann distribution, implying Eqn. 2; no further approximations will be made in this section. We proceed to examine the moments of gyrokinetic equation, order by order. The density moment at dominant order in δ is from which we formally establish that ϕ (0) depends only on the slow time variable. At next order in δ we obtain We introduce a time-average operator to extract the smooth-time behavior from this equation This time average extends over a period of time much longer than the short timescale (∆t −1 {ϕ, .}). Applying this average to Eqn. 13, we obtain The dominant-order equation for G n is from which, upon time averaging, we conclude that Ù G n is constant along closed streamlines of constant ϕ. Informally, we will say that Ù G n is a function of ϕ, although it can be multi-valued. For n = 2 we adopt the notation Note that, formally, we must exclude special points and lines where ∇ϕ = 0 (o-points, and the "separatrices" that include x-points) but these should occupy negligible volume in the x-y plane. At the next order, we will obtain the smooth evolution of G n , To this equation we apply two averages, the time average, and also an average along streamlines of constant ϕ. To define this average, we introduce a coordinate s which parameterizes these streamlines and satisfiesẑ × ∇ϕ · ∇s = 1 for convenience. Then we define The integral over s is either closed in the sense that the streamlines are closed, or effectively closed by periodic boundary conditions. The second term of Eqn. 18 is annihilated by the time average. Noting that {F (ϕ), A} = ∂ s (AF ), the third term is zero under the s-average, as is the last term after time average, since Ù G n+2 is a function of ϕ (0) by Eqn. 16. This is a crucial cancellation since the fluid moment hierarchy is consequently shown to be closed.
It is convenient to now introduce notation for the part of a field that varies on the fast timescale, i.e. the "fluctuating part", complementing the mean component defined by Eqn. 14: What results from the double average of Eqn. 18 can then be expressed To evaluate the nonlinear term of Eqn. 21, we must obtain dynamical equations for the fluctuating fieldsφ (1) andG (0) n . These come from Eqns. 13 and 16, respectively. The fluctuating part of Eqn. 16 is Subtracting Eqn. 15 from Eqn. 13, and using the Laplacian of Eqn. 22 to evaluate ∂ t f ∇ 2G (0) n , we find Finally, noting that¨∂ ts ϕ (0) ∂ s = 0 (from Eqn. 15), we obtain, from Eqn. 21, an expression determining the explicit time dependence of χ: where the partial time derivative is taken at constant ϕ (0) . To summarize, Eqns. 22 and 23 are the fast-time equations that determineφ (1) andG (0) n , which can be substituted into the slow-time equation 24 for χ, and coupled with the following equation (a repetition of Eqn. 15 written in terms of χ) to close the system: Noting that Û ϕ ≈ ϕ (0) andφ ≈φ (1) , we can, without introducing ambiguity, simply drop the superscripts in what follows.
The final system, Eqns. 22-25, has some noteworthy features. First, Eqn. 24 has the appearance of a heat transport equation, where the flux is carried by the rapidly varying pressure perturbationG 2 and the small amplitude fluctuating potentialφ. Note also how Eqns. 22 and 23 bear a strong resemblance to the fluid system given by Eqns. 8-9, where a similar ordering is satisfied, namelyφ G 2 .

Decaying turbulence
We will avoid the complications introduced by the instabilities that physically drive turbulence, and instead now consider decaying turbulence. (We note that a linear instability could be added to this fluid system using the non-resonant limit of the toroidal branch of the ITG or ETG mode, but this would require some care to maintain consistency with the ordering assumptions, as discussed in section 3.1.) Let us consider periodic boundary conditions, and include explicit dissipation using a fourth-order hyperviscosity term. Without drive terms, Eqn. 22 implies the rapid decay ofG n to zero, implying (∂ ts χ) ϕ = 0 (i.e. it only depends on the time via its dependence on ϕ). The variation of χ in ϕ (or more formally, its variation between distinct lines of constant ϕ) can be considered as an initial condition of our calculation. We need only then solve a single equation, which, neglecting superscripts for order and the subscripts of the time variable t s , becomes The electrostatic energy is conserved by the nonlinearity for arbitrary χ, which can be verified by multiplying the equation by ϕ and integrating over the x-y domain. Note that the resulting integral of the final nonlinear term of Eqn. 26 can be rewritten as − dxdy v E · ∇(ϕχ ∇ 2 ϕ), with v E =ẑ × ∇ϕ, which is zero using ∇ · v E = 0 and periodicity. Another quantity of interest is the enstrophy, which we will define here as The enstrophy balance equation is found by multiplying Eqn. 26 by −∇ 2 ϕ and integrating over x and y. Note that the presence of χ in the equation breaks enstrophy conservation if χ is a nonlinear function of ϕ. The nonlinear invariance of Z is associated with the inverse cascade of energy in Hasegawa-Mima turbulence. We thus expect to recover the spectra corresponding to the potential limit of the Hasegawa-Mima equation (i.e. where ∇ 2 ϕ ϕ; see Plunk et al. [2010]), if χ is small, and qualitatively different cascade when χ is sufficiently large.
The Hasegawa-Mima spectra can be derived in the rough "phenomenological" style, in terms of the fluctuation amplitude at scale , denoted ϕ , by assuming constancy of nonlinear flux of its nonlinear invariants (see e.g. Frisch [1995], Plunk et al. [2010]). For the forward cascade, i.e. at scales smaller than the scale of energy injection, the enstrophy flux, denoted ε Z , is assumed constant (independent of scale ), which is expressed as follows: with τ NL ( ) denoting the nonlinear turnover time. This leads to the scaling ϕ ∼ 2 ε 1/3 Z , implying a one-dimensional energy spectrum of E(k) ∼ k −5 . The constancy of the scale-by-scale flux of energy, expected for the inverse cascade at scales larger than the injection scale, is expressed as implying ϕ ∼ 4/3 ε 1/3 E and a spectrum E(k) ∼ k −11/3 . Because the additional nonlinear terms of Eqn. 26, henceforth called the "χ nonlinearity", formally break enstrophy conservation, we expect that if they are sufficiently strong, the inverse cascade should be eliminated and the forward cascade of Z replaced with a direct cascade of E. If this flux is carried by the HM nonlinearity, one might expect to observe the spectrum E(k) ∼ k −11/3 , as suggested by Plunk et al. [2019]. On the other hand, balancing the χ nonlinearity with the HM nonlinearity, scale-by-scale, implies the linear relation χ ∼ ϕ , i.e. χ ∝ ϕ, which would imply that enstrophy is actually a nonlinear invariant, preventing the forward cascade of E. For this reason, we may expect to observe an energy spectrum distinct from k −11/3 , whose steepness depends on the relationship between χ and ϕ , which itself depends on details of the turbulence.
Providing a definitive prediction of this relationship is beyond the scope of the present work, but a power law seems to be a reasonable possibility to explore, i.e. χ ∼ ϕ α . Note that any super-linear scaling α > 1 should lead to a spectrum shallower than k −11/3 , while a sub-linear scaling α < 1 would imply χ is not analytic in x and y. The fluctuating fieldsG 2 andφ could be especially active in regions of low E × B shear (see Eqn. 22), causing local extrema in the function χ, via Eqn. 24, so that a quadratic relationship prevails in such regions, χ ∼ ϕ 2 . Whether or not this seems plausible, assuming a simple nonlinear relationship will allow us to make the discussion now more concrete; qualitatively similar conclusions should apply for all α > 1. Let us consider the following form for χ: The nonlinear energy flux by the χ terms is then expressed as ε E ∼ λϕ 4 −4 , implying ϕ ∼ (ε E /λ) 1/4 , and the corresponding energy spectrum This spectrum should prevail in cases where the χ-nonlinearity dominates (e.g. large λ). At sufficiently low λ, one expects a return to the HM behavior, implying E(k) ∼ k −5 for the forward cascade. Some sort of hybrid behavior may also be possible, although the broad scale range needed for clear observation of this may be not be present for realistic conditions encountered in fusion plasmas. One might argue that, because the amplitude of fluctuations ϕ is generally expected to decrease as scales do, the cubic nonlinearity should be dominant at large scales, and subdominant at small scales. Thus, for sufficiently large λ, the energy cascade scaling ϕ ∼ (ε E /λ) 1/4 should hold from the injection scale, down to a transition scale, which can be found by balancing the HM nonlinearity with the χ-nonlinearity, i.e. ϕ 2 −4 ∼ λϕ 3 −4 . Defining the outer scale o as the scale of energy injection (or initial energy containing scale), and ϕ o = ϕ o , we can write the ϕ scaling as ϕ ∼ ( / o )ϕ o , so that the above balance occurs at the "transition" scale

Direct numerical simulations
To explore the behavior of the model, Eqn. 26, and test the theoretical predictions, we perform direct numerical simulations, assuming the simple quadratic form of χ(ϕ) in Eqn. 31. This introduces a nonlinearity that is cubic in ϕ, which can be treated pseudo-spectrally using a padding factor of 2 for dealiasing; higher order nonlinearities require additional padding [Hossain et al., 1992]. The boundary conditions for the simulations are periodic in x and y, and τ = 1 for all simulations. Fig. 1 compares the simulation results with the theoretical scaling laws. All simulations are initialized with randomly phased fluctuation amplitude of ϕ ∼ 1 around k = 1, falling off exponentially at higher k. Note that although the model assumes k 1 there is no conflict in using k > 1 for the simulations, as scaling symmetries of the model allow the results to be reinterpreted for k 1. The spectrum found for the λ = 0 case is roughly consistent with the theoretical power law k −5 expected for the potential limit of the HM equation. We note that similar results (not shown here) are encountered for λ 0.1. At larger λ, the breaking of enstrophy conservation is indeed observed in the time trace of Z, as the energy fills in the spectrum at large k. For the case labeled λ → ∞ in Fig. 1, the spectrum seems consistent with the theoretical k −3 prediction at scales smaller than the injection scale. Note that this limit is obtained by actually setting λ = 1 and simply removing the HM nonlinearity (i.e. the first term of Eqn. 26) from the equation, as can be formally justified by rescaling Eqn. 26 in the limit λ → ∞. Similar behavior is observed for λ 1. Intermediate values of λ show intermediate behavior.
One example is shown in Fig. 2, which seems to show evidence of a transition scale between the two theoretical power laws, giving some support to the predictions of a hybrid scenario described theoretically in the previous section. A more extensive set of simulations would be needed to test the predictions in detail, for instance the dependence of the transition scale

Discussion
A novel fluid system has been derived to describe the behavior of certain classes of quasi-twodimensional electrostatic magnetized plasma turbulence. A possible application is to describe the energy cascade in cases of streamer-dominated ETG turbulence (note the spectrum, noted to be close to k −11/3 , in Fig. 5 of Plunk et al. [2019]), where the nonlinear stability of elongated turbulent eddies is believed to stem from the two-dimensional character of the dominant instabilities, e.g. the absence of sufficient variation of the mode structure in the direction along the magnetic field [Jenko and Dorland, 2002]. This turbulent state is, however, sensitive to magnetic geometry, and seems to vanish when, for instance, the global magnetic shear is varied in such a way as to induce stronger parallel electron flow to the ETG mode. The ensuing dynamics then depends on kinetic physics involving the parallel streaming term, absent from two-dimensional models. A second possible application of the present model might be to describe ITG turbulence in cases where the zonal flows are suppressed. One candidate is a case observed with simulations of the HSX stellarator having surprisingly steep fluctuation spectra [Plunk et al., 2017], found to be close to k −10/3 . Although the presented model has limited application, it fills a significant gap in present theories describing gyrokinetic turbulence cascades, as it accounts for the essential nonlinear terms that arise when the cold ion approximation is invalid. These terms, it is found, alter the conservative properties of the nonlinearity, with significant consequences on the cascade, so that, even in the two-dimensional limit, the inverse cascade of energy can be shut down. The numerical simulations confirm that the size of the pressure perturbation (χ) can control the cascade type, and HM-like behavior can be recovered if it is sufficiently small. This may underlie the slow secular growth of large-scale zonal flows [Guttenfelder and Candy, 2011] and other coherent structures [Nakata et al., 2010] in simulations of ETG turbulence, and the related appearance of a Dimits shift phenomenon in near-marginal cases [Colyer et al., 2017].