Hasegawa–Wakatani and Modified Hasegawa–Wakatani Turbulence Induced by Ion-Temperature-Gradient Instabilities

We review some recent results that have been obtained in the investigation of zonal flow emergence, by means of a gyrokinetic trapped ion model, in the regime of ion temperature gradient instabilities for tokamak plasmas. We show that an analogous formulation of the zonal flow dynamics in terms of the Reynolds tensor applies in the fluid and kinetic regimes, where polarization effects play a major role. The kinetic regime leads to the emergence of a resonant mode at a frequency close to the drift frequency. With the objective of modeling both separate fluid and kinetic regimes of zonal flows, we used in this paper a methodology for deriving both Charney–Hasegawa–Mima (CHM) and Hasegawa–Wakatani models. This methodology is based on the trapped ion model and is analogous to the hierarchy leading from the Vlasov equation to the macroscopic fluid equations. The nature of zonal flows in the hierarchy of the Mima, Hasegawa and Wakatani models is investigated and discussed through comparisons with global kinetic simulations. Applications to the CHM equation are discussed, which applies to a broad variety of hydrodynamical systems, ranging from large-scale processes met in magnetically confined plasma to the so-called zonostrophy turbulence emerging in the case of small-scale forced, two-dimensional barotropic turbulence (Sukoriansky et al. Phys. Rev. Letters, 101, 178501, 2008).


Introduction
Understanding and controlling turbulence is one of the key elements for a successful magnetic confinement of a tokamak plasma. Turbulence due to drift-wave is a ubiquitous feature of magnetically confined plasmas. In the studies of magnetized plasmas, theories and simulations have predicted that drift-wave turbulence or interchange turbulence [1], in the case of a toroidal geometry, can generate mesoscale structures, such as zonal flows (ZFs), sheared flows, but also streamers. These mesoscale structures, ZFs and streamers, have different behaviors and their impacts on turbulence-driven transport show strong contrasts. A spontaneous transition to a turbulence-suppressed regime is sometimes observed and is known as low-to-high (LH), confinement transition [2]. Streamers appear to be closely associated with avalanche-type transport events and contribute to enhancing the transport owing to their radial elongated structures. Motivated by the experimental discovery of the LH transition, experiment and theory in the last decade have focused on whether the turbulence associated with the H-states might be regulated by interactions with ZFs. The generation of ZFs and their feedback on turbulence itself by shearing turbulent eddies and/or drift wave packets have been studied in Refs. [3][4][5][6].
The theoretical understanding of ZFs and streamers has first come from studies on the drift-wave model but also in the gyrokinetic turbulence modeling [7], e.g., the ion temperature gradient (ITG) mode and the trapped ion modes (TIMs), thanks to the assumption of adiabatic invariance of the magnetic momentum and of the parallel kinetic energy of particles trapped in the so-called banana orbits [8,9], as well as gyro-fluid descriptions in a 2D or 3D space (see, e.g., Refs. [10,11]). ITG modes are tokamak instabilities of interchange type (i.e., akin to the classical Rayleigh-Taylor instability in fluids), which determine the growth of perturbations in the poloidal angle because of charge separation effects due to an unfavorable radial variation of the so-called ∇B-drift particle velocity (see Figure 1). Since the latter is proportional to the particle kinetic energy component perpendicular to the local magnetic field, it depends on the particle temperature according to v ∇B = T eB 3 B × ∇B, where T is the particle temperature in the poloidal plane, e is the particle charge and B is the amplitude of the magnetic field B. While in a tokamak the plasma temperature is radially decreasing outward from the core of the vessel, the magnetic field intensity B, because of stability design, decreases in the direction opposite to the tokamak major radius (the intensity of the dominant toroidal magnetic field component being approximatively inversely proportional to the tokamak major radius R 0 ). As a consequence, the "outer" part of the poloidal section of the plasma with respect to the tokamak rotation axis has the gradients of both temperature and magnetic induction intensity pointing toward the axis. A perturbation in the poloidal angular direction will therefore mix the inner, warmer plasma, having a higher ∇B-velocity, with the outer plasma, colder and then relatively slower. This difference in the angular velocity corresponds to a charge separation and to an angular electric field component, which will foster the mixing of the warm and cold portions of the plasma because of two oppositely oriented E × B-drift velocity components (v E×B = E×B B 2 ) along the radial direction. Left frame: sketch of a poloidal section of a tokamak, where the orientation of the gradients of the plasma temperature (yellow dashed arrows) and of the intensity of the magnetic field induction (white dashed arrows) is represented with respect to the curvature of the torus. In the right frame, a cartoon depicts the basic mechanism of the ion temperature gradient (ITG) instability in a portion of the plasma in the unstable region where ∇T and ∇B are equally oriented. A harmonic perturbation (white line) in the angular direction mixes warm plasma particles, which have a relatively larger v ∇B velocity (yellow curled arrows), with colder particles. The difference in velocity is associated to a charge separation and to an angular electric field perturbation (white arrows), oriented so that the corresponding v E×B velocity (black arrows) further increases the amplitude of the harmonic perturbation. The magnetic field toroidal component, which gives the dominant contribution to the E × B velocity, points in the toroidal direction out from the plane, as indicated at the bottom of the right frame.
A review of physics of ZFs can be found in Ref. [12]. The generation of ZFs was also predicted in the geostrophic vortex equation by Charney in [13] and by Hasegawa and Mima in Ref. [14] in the case of drift-wave turbulence. This model aims to capture the essential physics while involving the least number of scalar fields as possible, in order to get a deeper physical insight and in the attempt to unravel the complexity of the phenomena under investigation. In this context, great interest has been recently attracted by the two-field descriptions (involving ion density and vorticity) provided by the Hasegawa-Wakatani [15] (HW) and Modified Hasegawa-Wakatani [16,17] (MHW) models for drift-wave turbulence. They consist in a more complex formulation of the earlier and somehow simpler Charney-Hasegawa-Mima (CHM) model, successful in describing a wide variety of turbulent process, both in atmospheric fluids [13] and in plasmas [14]. On the other hand, barotropic instability is the atmosphere-ocean counterpart of plasma resistive drift-wave instability. In particular Charney in [18] showed that the CHM model may describe the regime of geostrophic turbulence pertinent to the largest planetary scales. Even in its simplified barotropic version (with infinite Rossby deformation radius), the commingling of strong nonlinearity, strong anisotropy and Rossby waves give rise to complicated dynamics (see for instance Refs. [19,20]).
ZFs are also met in natural phenomena and the physics of ZF formation was also studied within the geophysical fluid dynamics community. Recently, in CHM flows with a small-scale forcing, the simulations have shown that the inherent anisotropy inverse energy cascade can lead to a new regime of turbulence, referred as the zonostrophic regime [21], a particular regime of geostrophic turbulence as shown in [22,23]. This regime is characterized by an anisotropic spectrum, the formation of alternating zonal (east-west) jets and nonlinear (Rossby solitary) waves called zonons. These modes exhibit low frequency oscillations and display similarities with the resonant trapped-ion-mode (TIM) met in plasma physics. More recently, the data obtained from the Nasa spacecraft Cassini in Ref. [24] have shown that Jupiter's troposphere in sharp contrast to the Earth's atmosphere, conforms to the regime of zonostrophic turbulence.
While considerable progress has been achieved in the understanding of the zonal flow physics, many aspects of the ZF dynamics remain nevertheless poorly understand. The generation of ZFs and their feedback to turbulence and transport are essentially nonlinear processes. ZFs are non resonant and are generated through the Reynolds tensor in the drift-wave turbulence. A clear indication of the key role played by ZFs was the recent observation in Ref. [25] at the Experimental Advanced Superconducting Tokamak (EAST) or in the DIII-D tokamak in Ref. [26], associated with a low frequency signal at a few kilohertz, i.e., at a much lower frequency than the usual geodesic acoustic mode, the geodesic version of ion acoustic modes. Central to all these physical turbulent systems, be they in the core region of the tokamak plasma or in the shallow rotating atmosphere, is the generation of ZFs, exhibiting low frequency oscillations, which are believed to be responsible for suppressing small-scale turbulence and stabilizing the turbulence. This brings us to an important question about the physical origin of these low-frequency oscillations while the ZF frequency is usually attributed to be zero in the hydrodynamic approach. From the work described in Ref. [27], the answer appears to lie in the fact that the nature of the ZF is strongly modified by nonlinear effects induced by the resonant character of (kinetic) TIMs, at least for the case of turbulence met in tokamaks.
The kind of dynamics involved in all these systems has the same general structure, which explains the many qualitative similarities between the systems. In particular, in the case of the zonostrophic regime, the low frequency oscillations-referred as zonons-are interpreted as nonlinear perturbations which essentially depend on resonant triad interaction. This was theoretically predicted in Refs. [21,28] and can be recovered by a simple CHM model including a small-scale external driver. Even more surprisingly, the same general structure of the ZF, at least driven by the Reynold tensor, is found in the tokamak models of Hasegawa, Mima and Wakatani. As exemplified by the case of the study of zonostrophic turbulence, the possible occurrence of these nonlinear structures allows insight in the strong turbulence regime described by the CHM equation and in the role played by forcing term. However, although the reduced tokamak models of Mima, Hasegawa and Wakatani exhibit systematic deficiencies in representing the kinetic effects, since such effects are not included, accurate simulations of the ZF production require a correct representation of the resonant interaction energy transfer mechanism (in particular from small-scales).
In particular, the main interest in the HW and MHW models lays in their debated capability to account for a self-consistent generation of zonal flows by small-scale turbulence, and thus in the possibility to shed light in the spontaneous transition to a turbulence-suppressed regime. In this respect, these models display some differences.
While the feature of zonal-flow generation, driven by the Reynolds tensor, was already present in the fluid CHM framework, the CHM model seems to be inappropriate to describe the turbulence at small-scale without the introduction of an external forcing term. The HW model, which like the CHM model for plasmas evolves the ion vorticity (at the E × B-drift ordering) and the ion density fluctuations, was developed to describe anomalous edge transport induced by collisional drift-waves. It contains a term more than the CHM model. This term is due to the electron dynamics parallel to the magnetic field, as it is deduced from Ohm's law, and linearly couples the fluid stream function, i.e., the electrostatic potential, with the ion density. Such parallel electron dynamics contribution enters with a coefficient, which is usually named the "adiabaticity parameter", since it measures the electron adiabatic response along the magnetic field lines. The MHW differs from the HW model because the zonal components have been subtracted from such a linear coupling, and a further contribution proportional to the background inhomogenous density (assumed to be constant in time) times the drift velocity is included. These corrections to the HW model have been included in Ref. [16] because it was mentioned that, without these terms, the original HW model was otherwise incapable of describing the zonal flow generation. This assertion, mentioned in Ref. [17] to be at odds with the CHM limit of the HW model and with the fact that the CHM, was in turn capable of describing a self-consistent zonal flow generation, was pointed out by numerical simulations of the HW model to be true only for small values of the adiabaticity parameter [17], for which both the MHW and HW models converge toward the 2D Navier-Stokes equation.
Thus, we need a reference model for generation of ZFs and their impact on the turbulence suppression that allows to characterize the nature of ZFs generated in the different reduced models such as CHM, HW and MHW. The Hamiltonian trapped ion model is a very useful and instructive model, even if restricted to plasma physics. The results from the reduced Hasegawa, Mima and Wakatani models also have to be regarded with some caution because the nature of ZFs depends on the effects at play in each model. A particularly insightful way to utilize and compare these reduced models is to examine the interactions of the drift-wave (or interchange-type ) turbulence with ZFs, by starting from the global trapped ion model.
In this paper, we address these points by providing a rigorous derivation of a MHW-type model from a reduced kinetic model for trapped ions in banana orbits, apt to describe the production of ZF through a resonant energy transfer mechanism [1,8,9]. The trapped ion model can be seen as an extension of the HW model, and allows the description of TIMs, a prototype of collective modes of kinetic nature. This model allows a fine description of the interchange turbulence in which the resonant interaction between TIM and trapped ions takes place through the precession motion of the latter. In revising the hierarchy of CHM, HW and then MHW models, we will address how this hierarchy is linked to a closure condition obtained from the kinetic model, which allows us to recover the Reynolds tensor shear effects only in the HW approach, a mechanism that may impact the turbulence.
The MHW-type model differs however from the one proposed in Ref. [16], to which it will be compared together with the HW model of Ref. [17]. We then discuss the physical limits in which the kinetic equations reduce to the HW and CHM model, and we finally discuss some properties of the MHW set we have described. In particular, we present numerical results obtained from integration of the reduced kinetic, Hamiltonian gyrokinetic model for the trapped particles [9] based on an action-angle model first discussed in Ref. [8], which corresponds to a simplified version of the upgraded version called "TERESA", by referring to one of its numerical parallel hybrid (OpenMI-Message passing Interface) implementations [29] to describe the process of the ZF formation.

Hamiltonian Model for Ion Temperature Gradient (ITG) Turbulence
The model, we consider for a tokamak plasma, describes the dynamics of TIMs in which the gradient of temperature gives the source of free energy. TIMs play a major role in the range of frequencies that is well below the parallel transit frequency. TIMs are then obtained by averaging the particle dynamics over fast scales, that is, over the cyclotron and bounce motions in the toroidal geometry. This task is made easier in the framework of the Hamiltonian-Jacobi formalism using action-angle variables. Due to their curvature drift, the orbits of trapped ions in a tokamak display a "banana" shape centered on the low magnetic field side. The low-frequency response for TIM is obtained by making a phase-angle average over the cyclotron phase and the bounce motion. This is allowed by the invariance of the total energy E = 1 2 mv 2 G + µB (x G ) and of the so-called adiabatic invariant µ = mv 2 G⊥ 2B(x G ) , which is appropriate at the time scales of interest for this problem. Here, the label G is a conventional notation which refers to the guiding centre and x G refers to (r, θ) polar coordinates. In agreement with the experimental conditions, we consider low β pl values and a poloidal field B θ lower than the toroidal magnetic field B ϕ . The modulus of the magnetic field is given by where B 0 is the minimal value of the magnetic field amplitude B at θ = 0 (R = R 0 being then the major radius and r = r 0 the minor radius). In this configuration, the action-angle coordinates of the trapped particle distribution are the precession angle α and the poloidal flux ψ. The poloidal flux is related to the poloidal field by dψ = −B θ R 0 dr and B θ B = ε q(r) , where q is the safety factor and ε = r R 0 is the inverse aspect ratio. The precession angle is related to the toroidal and poloidal angle coordinates by α = ϕ − q 0 θ. The gyro-average operator J 0 is approximated by Pade's relation (see Refs. [30,31]): This operator introduces a "banana" scale δ b corresponding to the width of the particle's trajectory in the ψ direction, while the gyro-phase average on the Larmor radius ρ s acts along the direction of the precession angle α. Rather than working with the adiabatic invariant µ, it is interesting to introduce the pitch angle parameter κ defined by the relation κ 2 = sin 2 θ 0 2 = 1−λ 2ελ where λ = µB 0 E . The use of the pitch angle κ, rather than µ, will allow us to characterize ions in the form of trapped populations (κ < 1) or of passing particles (κ > 1). Following the work of Kadomtsev and Pogutse [32], the bounce and precession frequencies are respectively given by the following relations: is the magnetic shear. K (κ) and E (κ) are the complete elliptic integral of the first and second kind, respectively. Trapped ions in a banana orbit are described by two invariants, the energy E and the pitch-angle κ, and by the distribution function f = f E,κ (ψ, α, t) fulfills the Vlasov equation: with the advective-term depending on the gyro-average operator J 0 defined in Lable (1) and expressed by means of the Poisson brackets The right-hand-side (RHS) expresses numerical particle diffusion (see the end of the section). The polarization effects have been taken into account in the quasi-neutrality equation by the introduction of the polarization density (the Laplacian term). This is due to the difference between the real trapped particle density and the bounce-averaged banana centre density. Assuming an adiabatic response for electrons, the quasi-neutrality condition δn e = δn i , where δn e and δn i are the fluctuations of the particle densities, reads: Here, n 0 (ψ) is the equilibrium density, equal for both ions and electrons. The operator = , which describes polarization effects, introduces the spatial scales ρ s and δ b . The first term on the left-hand side of Equation (5) comes from the adiabatic condition of electrons while the second term is due to the polarization charge. The gyro-averaged ion density n i is obtained by replacing the distribution f κ,E with J 0 f κ,E in the expression of the ion density, n i , defined by C e and C i are constants accounting for f p , the fraction of trapped particles, and for the ratio τ = T i /T e of ion to electron temperatures, C e = (1 + τ) / f p and C i = C e f p /τ. Trapped ion turbulence develops on a length scale of the order of the banana width δ b and on a time scale determined by Here, δ b is constant since we have neglected the dependence in κ and is found to be close to ρ s q 0 √ ε . In Equation (4), the diffusion coefficient D = D (ψ) is introduced for numerical purposes. In practice, in the numerical code, this diffusion coefficient is set to zero everywhere, except in a small buffer region located on the boundaries of the numerical box at ψ = 0 and ψ = 1. This choice allows to control the level of turbulence on the box boundaries (the buffer regions occupies less than 10% of the radial domain, in which a small diffusion process is applied with a constant coefficient In what follows, we will focus on the domain of interest, namely the core turbulent region where D is zero. It is interesting to recall here some general properties of the linear analysis of TIMs. By linearizing Equations (4) and (5), for an equilibrium distribution F 0 (ψ) given by Equation (28) of Section 5, one obtains, by considering a standard potential perturbation mode in the form δφ n e i(nα−ωt) , the following relation: with the usual Landau prescription on the imaginary part of ω. Note that, in Equation (7), we have neglected the polarization term in first approximation. As first noted in Ref. [33], a convergence of the expansion for the pole's calculation in Label (7) can be obtained by considering a development around where we have taken the average over the κ and E variables. Here, ω e * < 0 is the electron diamagnetic frequency, k ψ is the wave vector (for the variable ψ) and ρ s is the ion Larmor radius. Thus, it is possible to excite an oscillating ZF component when the last term in (8) becomes of order of −n ω d κ . We then have In Equation (9), the resonant amplification is produced by a three-wave process in which an interchange mode (of frequency ω T = 5 2 n ω d (κ) κ and of toroidal number n) decays into a resonant collisionless TIM (of frequency ω R = 3 2 n ω d (κ) κ and of toroidal number n) and a ZF (with a zero toroidal number).

Recall of Previous HW, MHW-Type and CHM Models for Drift-Wave Turbulence
In the following, we are going to show the connection between the Hamiltonian model of Section 2 with HW-and MHW-type equations. Before deriving the specific form the latter takes in this framework, we recall the HW and MHW-type formulations which have been proposed in two previous works [16,17]. They consist of a two-field theory involving the ion vorticity ξ ≡ φ and the ion density fluctuation, which in our notation reads n i − 1, by assuming that J 0 → 1. In the coordinate system (α, ψ) and in the non-viscous regime we consider, the MHW formulation proposed by Numata et al. [16] reads where δφ = φ − φ α , δn i = n i − n i α and C ≡ T e /(n 0 ηΩ i e 2 )k 2 || is the adiabaticity parameter that depends on electron-ion collisions through the resistivity η, and on the toroidal wave-length component k || related to the poloidal scale length by 2π/k || ≡ L || L α , L α = 2π being the maximum angle α. The main difference with respect to the HW modelling (In a three-dimensional system, the distinction between the HW and MHW models disappears and the "adiabaticity parameter" is not a constant, but rather an operator proportional to ∇ 2 . However, since we are interested here in the radial contribution of the turbulence and not in the high-frequency contribution induced by high-frequency geodesic adiabatic mode (GAM) in zonal flow (and usually triggered by circulating particles), a two-dimensional approach is sufficient. Here, we have introduced the term "low-frequency" zonal flow (LFZF) to discriminate between the quasi-zero zonal flow (due to the Reynolds tensor) and the high-frequency GAM.) lays in having taken out the averages over α, so that the quantities δφ and δn i express the potential and ion density fluctuations, respectively. The parameter G is the normalized gradient of the density equilibrium profile, G ≡ (∂ ln n 0 /∂ψ). The particle flux in the ψ direction can be generated by out-of-phase fluctuations of φ and n i in α and depends on ψ itself, according to This acts as a source term for the conservation of the energy E of this MHW model (label (M) in the time derivatives below), and of the enstrophy W (M) , The unmodified HW set of equations, which has been considered by Pushkarev in Ref. [17], consists instead of the equations Differently from Equations (10) and (11), besides the absence of the last term on the right-hand-side (RHS) of the density equation, the coupling term does not depend on averages over the α-coordinate.
Labeling with " (U) ", the time derivatives of the energy and enstrophy in this unmodified HW model, they obey dE dt

Applications to CHM and HW Turbulence
At this point, we would like to come back to the observation of Numata et al. in Ref. [16], of fast emergence of ZF in the MHW model. It was claimed by these authors that the original HW model does not predict the formation of ZFs, which appears to be in contradiction with the CHM results of Connaughton et al. and reported in Ref. [34], considering the fact that the HW model includes the CHM model as a limiting case. We discuss here the importance of addressing the problem of drift-wave turbulence in the framework of the reduced Hamiltonian model for trapped particles, in order to gain a better understanding of the differences between the HW and MHW models. One of our objectives is to find a rigorous justification for the nature of ZF observed in the different models. This requires one to obtain a self-consistent dynamics equation. The problem is related to the choice of the closure condition whose origin comes from the nonlinearity and from the kinetic nature of the trapped ion model we consider. We choose to start with the CHM model.

The CHM Model
Assuming a homogeneous equilibrium density n 0 in ψ (a choice to simplify the presentation), it is possible to recover a CHM-type equation from Equation (A4) in the following form (see Appendix A for details): Note that now Equation (19) (20) which clearly indicates that the nature of the emerged ZF is driven only by the Reynolds tensor. The similarity between the reduced Hamiltonian model for trapped ions and the CHM equation had been already pointed out in Ref. [1].

The HW and MHW Model
The HW equation can be obtained with a similar analysis. By integrating Equation (A5) over α and by replacing ξ by −C i φ, we obtain From Equation (A4), the coupling function F(n i , φ, φ α ) in Equation (21) is straightforwardly identified as Using the quasi-neutrality Equation (5), a little algebra leads to Formula (A7) leads to Some remarks are due: (i) In the MHW, the mean coupling function F(n i , φ, φ α ) α is given by = CC i (δφ − δn i ) α , which is strictly zero. This indicates that only the Reynolds tensor plays a major role in the generation of the ZF since in that case the second term in the RHS of Equation (24) disappears and, as a consequence, both CHM and MHW models are equivalent as far as the emergence procedure of ZFs is concerned.
In the HW model, as shown in Equation (24), there are two basic effects in the dynamics: the first term at RHS shows the growth of ZF induced by the drift-wave turbulence (via the Reynolds tensor) while the second term indicates the possibility of the suppression of turbulence by a ZF shear mechanism, which takes place at a larger scale in the ψ space. Rather than proceeding scale by scale in a local manner, energy transfer is here induced by competing processes between the smallest and the largest scales. Thus, the strong growth of ZF is expected and usually observed in both CHM and MHW models, whereas the suppression of drift-wave turbulence induced by a ZF shear is only expected in the original HW model, provided that dn 0 dψ = 0.

The Kinetic Nature of ZF Driven by TIM
In the TIM model, it is also possible to connect the time variation of φ α of ZFs to the ion pressure fluctuations, here driven by an interchange-type turbulence, a physical mechanism which is absent in the HW model. Since the model has been already derived in a previous work [27], we just recall here the main features relevant to the ZF dynamics. Assuming that dn 0 dψ is zero, in order to simplify the treatment, we just consider an ITG instability driven only by an initial temperature gradient, which corresponds to suppressing the second term dn 0 dψ ψ ∂φ ∂ψ ∂φ ∂α α dψ in (24). We focus here on the effects induced by the ion pressure on the ZF dynamics. The evolution of the ZF is now described by the following equations: assuming that there is no dissipation (i.e. D → 0) and where In Equation (25), the first term at RHS denotes the Reynolds tensor, while the two following terms are a straightforward consequence of the kinetic effects. The third term at RHS of Equation (25) is linked to the heat flux Q α defined by Equation (26).
In a similar way, a mean ion pressure is also driven by nonlinear processes. Assuming that the gyro-average operator J 0 → 1 and D = 0, and that the second order effects (i.e., the second moment in energy of the distribution function f κ,E ) are negligible, the mean pressure obeys an equation expressing, on average over α, a Lagrangian advection by the flow φ (see Ref. [27] for more details): It is worth pointing out at this stage that choosing P i equal to zero in Equation (25) leads to the HW model expressed by (24) in the case of an homogeneous equilibrium density. Formula (25) makes evidence of a role of the Reynolds tensor, wider than it is usually evidenced in fluid HW-type models. In particular, this expression shows that the ZF is impeded by the resonant interaction with collisionless TIMs, leading to a low-frequency ZF mode with a frequency close to ω d0 .
Furthermore, when adiabatic conditions set in (or, equivalently, when C i → 0), from (5), we have

Numerical Simulations
In this section, we discuss the results of numerical simulations performed in order to elucidate some of the key features of the interchange turbulence in the presence of TIMs.
Our numerical scheme is based on the integration of Equations (4) and (5) in a self-consistent way. The numerical algorithm is based on a semi-lagrangian scheme detailed in Ref. [14], allowing the integration of the distribution function directly in phase space. It is worth noting that a semi-Lagrangian scheme was also applied in [35] to improve the efficiency of numerical models of the atmosphere. Such a numerical scheme allows an accurate integration of the Vlasov equation along its characteristics. The code employs a classic time splitting scheme to separate the treatment of the advection due to the drift-frequency ω d but uses a full 2D advection in space to treat the advection due to the electric potential. It must be pointed out that semi-Lagrangian Vlasov simulations are slowly introduced in place of the well-known Lagrangian particle-in-cell simulations for two main reasons: the lack of numerical noise (in the sense that, in the Eulerian approach, the graininess parameter g = 1 n 0 λ 3 D tends to zero and where λ D is the Debye length and n 0 the particle density) and the very good resolution of the distribution function in phase space provided the dimension of the momentum space is as lowest as possible (depending on the computational constraints).
We focus on the nonlinear generation of ZFs by varying the parameter C i related to polarization effects. Equations (4) and (5) have been solved using a backward semi-Lagrangian scheme [36]. Semi-Lagrangian Vlasov simulations have identified two different kinds of ZFs of somewhat different nature, the first kind of ZFs driven by the Reynolds tensor, as expected, in agreement with the results of the standard fluid HW model. The second type of ZFs is characterized by the resonant coupling with TIMs, a physical process that is not taken into account in the reduced HW model.
In numerical simulations, quantities are normalized as follows: the time is normalized to the inverse drift frequency ω −1 d0 , the poloidal flux ψ is given in ψ units (with ω d0 = q 0 T 0 eB 0 r 0 R 0 ), the energy E is normalized to T 0 . The electric potential is expressed in ω d0 ψ units and the constants C e and C i introduced in the quasi-neutrality Equation (5)  . The bounce and drift frequencies ω b and ω d depend explicitly on the pitch angle parameter κ (and of course on the energy E) and are given by Equations (2) and (3). The initial distribution function is given by: In Equation (28) the quantity φ 0 denotes the initial flow corresponding to the interchange case. Here, τ = ψ T 0 dT 0 dψ is the normalized ion temperature gradient. The starting point for an investigation of interchange turbulence is the reduced gyrokinetic Vlasov Equation (4) coupled with the quasi-neutrality equation condition (5) initiated in an equilibrium state with a perturbation term of type: The reason for introducing the sin (πψ) factor in Equation (29) on the right-hand side instead of the standard potential perturbation in α is that this function is the marginal solution obtained in the linear analysis of the TIMs. Here, we choose a perturbation on the potential φ pert = 10 −4 .
Without dissipation, three energetic subsystems interact to produce the complexity observed in the interchange-type turbulence: the kinetic energy of the plasma the energy of the zonal flow noted here by (31) and the potential fluctuation of turbulence The energy components defined in Equations (30)-(32) verify the conservation law Here, we have studied in detail the effects of the polarization term. A series of simulations has been performed for different values of the parameter C i , the other physical and numerical parameters being kept constant. In all simulations presented below, we have taken an ion temperature gradient of τ = 0.15, chosen above the threshold of the ITG instability given by τ s =  Figure 2, in correspondence with non-negligible values of the electric potential. We choose to start with a small value of C i = 0.25, i.e., well inside the Reynolds' tensor-dominated regime. Results are presented, in the top panel of Figure 2, the time evolution of the zonal flow energy E ZF (in red and in solid line) together with its turbulent energy counterpart E turb (in blue and in dotted line). As shown in Figure 2, the dynamics of the plasma are first governed by the growth of turbulence, as a result of the ITG instability. Here, it is the gradient in temperature that constitutes the energy source for the growth of TIMs, when the temperature gradient exceeds a critical value of τ s . The case of an initial density gradient and of a resulting shear flow generated by a Kelvin-Helmholtz instability was presented in Ref. [37]. Here, the instability linked to a density gradient is not excited in the linear phase of ITG, but can appear later when nonlinear effects take place leading to the generation of shear flows. Note that at the saturation of ZF, in top panel in Figure 2, the turbulent component of the energy has strongly decreased. In the bottom panel in Figure 2, the corresponding evolution of the zonal flow component φ α , taken at ψ = 1 2 (in ψ unit), is plotted in a logarithmic scale. Aside from the fast oscillation in the linear phase (linked to the linear TIM), no oscillatory behavior is observed in the saturation regime, a signature of the role played by the Reynolds tensor in the fluid regime. In the bottom panel, the time evolution of the squared mean potential φ α in a logarithmic scale. We observe that when the zonal flow saturates at high level, the turbulence is strongly reduced. Here, the polarization parameter C i is 0.25. In the bottom panel in Figure 2, the corresponding evolution of the zonal flow component φ α , taken at ψ = 1 2 (in ψ unit).
We will now address how ZF is modified in the kinetic regime by its coupling with TIMs leading to a low-frequency oscillatory behavior, a mechanism that is different from the standard Reynolds stress met in hydrodynamics, and which is not taken into account neither in the HW model, nor in the MHW model. A second simulation is performed with C i = 0.875. Results are shown in Figure 3 to be compared to those used in Figure 2. There is a qualitative change in the interchange turbulence, which now exhibits low-frequency oscillations just before the resonant amplification of the ZF component.
Zero-frequency ZF is quite non resonant and it is relatively easy for the Reynolds tensor to drive it up in a nonlinear way. While ZFs back-react upon turbulence by shearing, so weakening the source of their generation, it was shown in Ref. [37] that Kelvin-Helmholtz (KH) instability may act as a damping mechanism for ZFs. However, this latter mechanism becomes inefficient when the nature of ZF is modified leading to a time-varying ZF, which becomes now sensitive to resonant amplification through parametric three-wave scattering. Thus, a transition from the fluid regime to a kinetic one for ZF is expected when the polarization increases. Here, the LFZF generation is a nonlinear process in which shorter-scale fluctuations transfer their energy to larger-scale potential structures. In addition to the first phase of growth of TIM turbulence (10 ≤ tω d0 ≤ 20), the nonlinear interaction leads to the growth of an oscillating mode on a frequency close to the drift frequency ω d (κ). Such a behavior can be observed in the intermittent regime of turbulence and is usually associated with a turbulence burst. In order to illustrate the rich nonlinear dynamics that underlie the resonant amplification of ZF, we have performed a new numerical simulation over a long time, for the same value of the polarization parameter C i = 0.875, but by diminishing the initial perturbation level of the potential by a factor ten, in order to follow the time evolution of the system in the (strong) intermittent regime of turbulence. This time, the time step is smaller, tω d0 = 5 × 10 −4 the phase space sampling is N ψ = 256 (here, we have reduced the sampling in α to N α = 512) and we have used N κ N E = 16 × 128 values in pitch-angle and energy. Numerical results are presented in Figures 4-7. On the top panel in Figure 4, we have plotted the temporal evolution of the energies of the zonal flow component E ZF (in red and in solid line) and of the turbulent component E turb (in blue and in dotted line). The zonal flow energy displays three different stages during its evolution: a first phase characterized by a weak growth-rate for tω d0 < 18, followed by a first resonant peak at tω d0 18 and finally a second resonant process, now located at time tω d0 25, leading to a slow reduction of the turbulence after this time. On bottom panel in Figure 4, the temporal evolution of the mean potential is represented in a logarithmic scale: we see clearly the emergence of different behaviors during the time evolution. While the beginning of the simulation is characterized by high-frequency oscillations (due to TIMs), followed by the bursting behavior at a lower frequency, the oscillatory behavior disappears, at saturation, indicating that a fluid-type regime is now recovered and dominated by the Reynolds tensor. The zonal flow energy displays three different stages during its evolution: a first phase characterized by a weak growth-rate for time tω d0 < 18, followed by a first resonant peak at tω d0 18 and finally a second resonant process, now located at time tω d0 25, followed by a slow reduction of the turbulence. The plot of φ α was taken at ψ = 1 2 (in ψ unit). Figure 5 shows the spectrum in frequency of the electric potential for the considered numerical simulation (over the total time and at a fixed point ψ = 1 2 ): we observe clearly both components of ZF, the quasi-zero part driven by the Reynolds tensor and the oscillating component , at a frequency close to ω ZF, num 2.5 ω d0 , a somewhat higher value than the predicted value of ω ZF = ω d (κ) 2.11 ω d0 (using a magnetic shear of s = 2 ). It must be pointed out that the resonance takes place for the dominant streamer mode, i.e., n = 1. Indeed, both mode n = 1 and its harmonics n = 2 are excited, as can be seen in Figure 5, where we have plotted the electric potential. This feature resembles the excitation of the nonlinear zonon mode, usually triggered by the same zonal wave vector of Rossby-Haurwitz waves (RHWs). In Figure 5, we observe also the growth of the mode ω T 5 2 × 2.11 ω d0 = 5.275 ω d0 and of its harmonics 2ω T 10.55 ω d0 . In the linear regime, the frequency of TIM is reduced to a value of ω R . However, TIM can propagate in the direction of the precession motion of ions allowing strong resonance with precessing ions: therefore, a slow decrease of TIM frequency is also possible when nonlinear trapping effects are important. Such a process has been already predicted by Morales and O'Neil in [38]. Thus, due to strong nonlinear effects, the frequency of the kinetic mode ω R , initially close to the (linear) value 3 2 n ω d (κ) κ = 3 2 × 2.11 ω d0 3.1 ω d0 can slowly decrease to reach a value close to ω T 2 , which finally leads to an estimation, for the ZF frequency, of ω ZF ∼ ω T 2 ∼ 2.5 ω d0 . It is this way that harmonics of ZFs can be excited.  Figure 4, the spectrum in frequency of the electric potential is shown. We observe clearly both components of zonal flow (ZF), the quasi-zero part driven by the Reynolds tensor (fluid regime) and the low-frequency zonal flow (LFZF) component (in the kinetic regime), at a frequency close to ω ZF, num 2.5 ω d0 . Figure 6. Continuing the results of simulation shown in Figures 4 and 5, the electric potential is plotted at three different times when the resonance takes place, showing clearly that the dominant modes are n = 1 and n = 2 in the asymptotic regime. We have identified a nonlinear coupling of zonal flow (ZF) mediated by nonlinear collisionless trapped ion modes (TIMs), leading to a strong reduction of the interchange turbulence. Figure 6 shows the behavior of the electric potential at three different times when the resonance takes place, showing clearly that, in the asymptotic regime, the dominant modes are n = 1 and n = 2. We have therefore identified a nonlinear coupling of ZF mediated by nonlinear collisionless TIMs, leading to a strong reduction of the interchange turbulence.
In Figure 7, we have represented the mean ion pressure P i α , as a function of the poloidal flux ψ, and the corresponding electric potential at two times: the formation of a transport-barrier region (located around ψ ∼ 0.4), where the temperature gradient becomes negligible, is made evident. This also indicated that the ITG instabilities are now strongly reduced in this region. Here, it is the polarization source that triggers the barrier and reduces the level of turbulence (see Ref. [39] for more details). The last example, shown in Figures 8-10, describes the dynamics of the resonant mechanism in an accurate way, when we increase the polarization parameter to a value of C i = 1.25, i.e., well inside the kinetic regime of ZF. The numerical simulation was carried out using the same physical and numerical parameters of Figures 2 and 3. We now focus on the transition toward the kinetic regime of ZF, which is triggered by the resonance of TIM wave with precessional trapped ions. An intermittent regime is observed (as expected in the dynamics previously shown in Figure 4) and we focus here on the dynamics just before the occurrence of the first burst, when the nature of the ZF is modified. On the bottom panel in Figure 8, we have plotted the electric potential, in a logarithmic scale, as a function of time. The oscillating nature of the ZF is clearly visible, even at saturation. In this saturation phase, ZF manifests as linear TIMs with a frequency ω R = 3 2 n ω d (κ) κ 3.1 ω d0 and no frequency shift. The turbulent energy density peaks at the small n modes. These modes harbor the most energetic TIMs in the system. This mechanism is analyzed in more details through the phase space diagnostics shown in Figure 9 and the spectra in wave-numbers (i.e., in toroidal numbers) are plotted in Figure 10. Figure 9 illustrates the details of the electric potential, when the resonance takes place: we observe clearly the formation of streamers (i.e., nonlinear structures elongated along the ψ direction). The intermittence is characterized by the generation of a turbulence burst at time tω d0 25, as the result of the emergence of the resonant TIM with toroidal number n = 1 (note that the initial perturbed mode n = 5 has disappeared at that time). Diagrams in Figure 10 reveal the occurrence of (nonlinear) streamers with a broad spectrum (with a toroidal number of n ∼ [16][17][18][19][20][21][22][23][24], as can be seen in the plot shown at time tω d0 = 25.6. At time tω d0 = 26.8, an inverse cascade takes place leading to a dominant mode on the toroidal number n = 9, as indicated in the bottom panel in Figure 9. In the bottom panel, the time evolution of the squared mean potential φ α in a logarithmic scale. There is a qualitative change in the interchange turbulence, which exhibits a weaker simultaneous growth of E ZF and E turb until saturation, at time tω d0 25, and an oscillatory behavior of φ α at low frequency. Here, the physical parameter is C i = 1.25. In the saturation phase of the instability, zonal flow (ZF) features linear trapped ion modes (TIMs) (with a frequency ω R = 3 2 n ω d (κ) κ 3.1ω d0 with no frequency shift. The plot of φ α was taken at ψ = 1 2 (in ψ unit).

Figure 9.
Diagrams show that streamers are also occurring but on a higher mode (with a toroidal number of n = 20), as can be seen on the plot shown at time tω d0 = 25.6. Here, the polarization factor is C i = 1.25. Figure 10. Spectra in wave-numbers, i.e., in toroidal numbers. The spectra reveal the occurrence of (nonlinear) streamers with a broad spectrum (with a toroidal number of n∼16-24) as can be seen on the plot shown at time tω d0 = 25.6.

Connections with the Zonostrophic Turbulence in Geophysical Fluid Dynamics
It is well known that the CHM equation and the quasi-geostrophic equation have the same structure. This makes it possible to establish a direct link between drift-wave turbulence in magnetically confined plasmas and the quasi-geostrophic turbulence in geophysical fluid dynamics (GFD). Both systems are approximately two-dimensional respectively because of the strong guide field applied to a magnetized plasma and because of the fast planetary rotation and strong density stratification in the quasi-geostrophic limit. Curiously, the same symbol β is used, although with completely different meanings, to characterize both the strong-guide field limit in a plasma and the regimes of quasi-geostrophic turbulence. As already recalled in the introduction, in a plasma, β pl labels the ratio between the kinetic and magnetic pressure, and it is small in the strong guide field approximation. In a tokamak, for example, for which an approximated two-dimensional description for the poloidal dynamics is allowed by the strong toroidal (guide) magnetic field, typically β pl = 0.01. On the other hand, in geostrophic turbulence β denotes the ratio β = Ω R , where Ω and R are the angular velocity and the radius of a rotating fluid spherical shell. As mentioned earlier in the introduction indeed, computer simulations of 2D barotropic turbulence on a β-plane surface of a rotating sphere gave rise to a classification of flow regimes, which become possible in these geostrophic systems, depending on the value of β. One of these is the regime of zonostrophic turbulence, which is characterized by an anisotropic inverse cascade and by slowly varying structures of strong alternating zonal jets. The barotropic vorticity equation (on the surface of such rotating sphere) is then given by where ξ is the vorticity, φ is the stream function with ∇ 2 φ = ξ, f = 2Ωsinθ is the Coriolis parameter (the planetary vorticity), θ is the latitude and ϕ the longitude. In (34), the Jacobian J is defined by ν is the hyper-viscosity coefficient and p the power of the hyper-viscous operator (p = 4 is usually used) and λ is the linear friction coefficient that expresses the large-scale friction. In Ref. [22], the authors have considered, in the CHM model, a small-scale forcing term S that pumps energy into the system at a constant rate and that leads to a new regime of the geostrophic turbulence, the so-called zonostrophic regime. However, no justification on the origin of such source term in the forced CHM model was given. The development of this new regime of turbulence is characterized by the existence of a new class of nonlinear waves dubbled zonons.
Indeed, Equation (34) can be solved using the decomposition of the stream function in spherical harmonics Y m n (sinθ, ϕ), leading to a class of linear Rossby-Haurwitz waves (RHWs) of frequency given by Fixing the reference unit of length so that R = 1 eliminates the difference between the indices of the spherical harmonics and wave numbers. The decomposition takes the usual form where n and m are the meridional and zonal wave numbers, respectively (and N a truncation index). It must be pointed out that we have kept here the notation (n, m) for the standard toroidal and poloidal numbers, respectively, and that the roles of (n , m ) and (n, m) are inverted to tokamak's with respect to toroidal geometry since m = 0 refers to zonal jet while n = 0 represents ZF in magnetically confined plasmas. We recall indeed that, in the 2D description of drift-kinetic turbulence in a tokamak, the restriction to a planar geometry corresponds to the k ϕ = 0 condition, which holds for both zonal flows and streamers. In the 2D description of zonostrophic turbulence the restriction to a 2D geometry is given instead by the condition k r = 0. The correspondence between the (r, θ, ϕ) set of toroidal coordinates in drift-kinetic turbulence in tokamaks and the (r , θ , ϕ ) set of spherical coordinates in zonostrophic turbulence in the atmosphere, also with respect to the zonal and nonlinear structures encountered in the two frameworks, is summarized in the scheme of Figure 11. In the 2D description of drift-kinetic turbulence in a tokamak, the restriction to a planar geometry corresponds to the k ϕ = 0 condition, which holds for both zonal flows and streamers. In the 2D description of zonostrophic turbulence, the restriction to a 2D geometry is given instead by the condition k r = 0. The correspondence between the (r, θ, ϕ) set of toroidal coordinates in drift-kinetic turbulence in tokamaks and the (r , θ , ϕ ) set of spherical coordinates in zonostrophic turbulence in the atmosphere is shown here.
The situation observed in the zonostrophic regime resembles the situation met in previous sections where LFZF is generated (see Table 1 in which we have summarized the different waves implicated in the turbulent processes). We can say that, while ZFs are met in both plasmas and in the atmosphere, streamers in poloidal plasma turbulence are the correspective of zonons in zonostrophic turbulence. Our investigation reveals that, in tokamak turbulence, it is the small-scale wave-particle interactions (or nonlinear streamers) that determine the resonance conditions of the residual LFZF. Thus, a key physical point is that this low-frequency ZF can be excited by resonant triad. Moreover, in Ref. [28], the authors have proposed that RHWs can produce ZF through an energy transfer mechanism, triggered by resonant triads. Thus, the turbulence is here dominated by waves that are involved in triad interaction. What is more interesting is that, in the zonostrophy turbulence, the zonon frequency ω zonon (n , m ) is equal to the frequency of the most energetic RHWs (by fixing a value of n ) that excite zonons, as shown in Refs. [21][22][23][24]28,40,41]. In particular, ω zonon is proportional to m , a feature which exhibits a similarity to the property of linear resonant kinetic (TIM) modes ω R = 3 2 n ω d (κ) κ (note that streamers are the nonlinear version of TIM) . Thus, the relationship between both types of waves is quite straightforward. Concerning the LFZF generation, there exists also, in the zonostrophic regime, zonal zonons with a frequency ω zonon (n , 0) = 0 whereas ω RHW (n , 0) = 0 for all values of the meridional wave numbers. Table 1. Summary of different types of waves and coupling processes in both interchange-type turbulence and zonostrophic turbulence. Finally, we would like to comment on the relation between the CHM equation given in (19) and the trapped-ion model. Nonlinear small-scale effects are here introduced, in the right-hand-side member, which indeed disappears when taking the average over the α variable. This property shows that ZF is not affected by this term. It is not the case of the small-scale forced version of the CHM used in the study of the geostrophic turbulence: along with the linear RHWs, nonlinear small-scale structures (zonons) also emerge and an intense energy exchange between zonons and the ZFs can take place.

Conclusions
In reviewing some recent results about the emergence of zonal flows in the ITG-instability-driven turbulence in tokamak plasmas, we have compared a kinetic model for the description of TIMs with fluid CHM and HW-type models for turbulence. We have shown that the zonal flow generation in collisionless plasmas is controlled not only by the Reynolds tensor, usually recognized as the main actor in the fluid approach, but explicitly depends on the microscopic physics that govern the ion temperature gradient instability. We have found that the Hasegawa-Wakatani model can be recovered from the gyrokinetic trapped ion model when interchange-type turbulence is considered, and that the zonal flow dynamics is strongly impeded by polarization effects. It is the Laplacian operator acting on the electrostatic potential that determines the transition between a turbulent regime dominated by the Reynolds tensor and a new regime of kinetic nature, where streamers and shear-flow-driven vortices develop inside a low-frequency oscillating zonal flow. In this investigation, the use of a reduced HW model, or its modified version, has been particularly fruitful since, on the one hand, it contains the physics related to the Reynolds stress, a feature that is thought to play a fundamental role in the ZF emergence, and, on the other hand, it exhibits mathematical properties that allow us to make connections with the gyrokinetic trapped ion model in a rather compact and general form. The feature of ZF dynamics that we have discussed appear to be of quite general interest, as they are relevant to physical processes that take place in other natural systems, in particular in the zonostrophic regime of planetary atmospheric turbulence. This research will hopefully help to understand whether kinetic effects, occurring at small spatial scales, concur to the onset of zonotrophic turbulence regimes. as long as the quasi-neutrality condition (5) holds. The interest in this formulation is therefore in the identification of the coupling function F(n i , φ, φ α ) through the derivation of these equations from the Hamiltonian model discussed above, and in the comparison of Equations (A5) and (A6) with the MHW and HW models discussed in Refs. [16,17], which have been instead partially built on the basis of phenomenological assumptions rather than by analytical derivation. This comparison is discussed in Section 4.