Abstract
Geodesic acoustic modes (GAMs) are oscillations of the electric field whose importance in tokamak plasmas is due to their role in the regulation of turbulence. The linear collisionless damping of GAMs is investigated here by means of analytical theory and numerical simulations with the global gyrokinetic particleincell code ORB5. The combined effect of the phase mixing and Landau damping is found to quickly redistribute the GAM energy in phasespace, due to the synergy of the finite orbit width of the passing ions and the cascade in wave number given by the phase mixing. When plasma parameters characteristic of realistic tokamak profiles are considered, the GAM decay time is found to be an order of magnitude lower than the decay due to the Landau damping alone, and in some cases of the same order of magnitude of the characteristic GAM drive time due to the nonlinear interaction with an ITG mode. In particular, the radial mode structure evolution in time is investigated here and reproduced quantitatively by means of a dedicated initial value code and diagnostics.
Decay of geodesic acoustic modes due to the combined action of phase mixing and Landau damping.
A. Biancalani, F. Palermo, C. Angioni, A. Bottino, F. Zonca.
MaxPlanck Institute for Plasma Physics, 85748 Garching, Germany
ENEA C. R. Frascati, Via E. Fermi 45, CP 6500044 Frascati, Italy
Inst. for Fusion Theory and Simulation, Zhejiang University, 310027 Hangzhou, P. R. China
contact of main author: www2.ipp.mpg.de/biancala
1 Introduction
Plasmas in magneticconfinementfusion plasmas are turbulent, due to the presence of strong spatial gradients. Turbulence is often observed accompanied by zonal (i.e. axisymmetric) radial electric fields, generating poloidal ExB flows in tokamaks. These zonal flows are usually constituted by two components, one nearly constant in time, dubbed zerofrequency zonal flows (ZFZFs) [1, 2, 3, 4, 5, 6, 7], and one oscillating with a frequency of the order of the sound frequency (, where is the sound speed and is the tokamak major radius), dubbed Geodesic Acoustic Modes (GAMs) [8, 9]. The importance of both components is due to their nonlinear interaction with turbulence, for they shear and distort convective and turbulent cells leading to the turbulence saturation and consequently to a reduction of the transport. The peculiarity of GAM oscillations resides in the different shearing efficiency that the zonal flows have in relation to their oscillatory behavior. In fact, while ZFZFs suppress turbulence efficiently, oscillations make the action of the zonal flow less effective [10, 11]. ZFs/GAMs are also known to transfer energy to the longwavelength components of the turbulence [12, 13].
Historically, the Landau damping has been recognized as one of the main damping mechanisms of GAMs. The Landau damping is due to the nonuniformity of the ion distribution function in velocityspace, and it is therefore an intrinsically kinetic effect. Its evaluation has been performed with increasing accuracy, where finiteorbitwidth (FOW) of the passing ions were initially neglected [14, 9], then included to the first order [15, 16, 17, 9], then included to higher orders [18, 19]. Finally, the effect of the flux surface shape was also included [20]. All analytical calculations of the Landau damping of GAMs have been performed so far by assuming a uniform plasma, i.e. by neglecting the effect of a gradient of the equilibrium density or temperature profiles.
More recently, the phase mixing, which was previously investigated for nonionized fluids [21], for Langmuir waves in plasmas [22], and for shearAlfvén waves in plasmas [23, 24, 25], has also been investigated for GAMs analytically [9, 26]. Contrarily to the Landau damping, the phase mixing can be studied also in fluid theory, whenever the local spectrum of oscillation of a wave  the continuous spectrum  varies in space (for the GAM in a tokamak, this happens due to the variation of the plasma temperature across the magnetic fluxsurfaces). In this case, the spatial shape of a perturbation will be deformed in time due to the effect of the different frequency of oscillation, occurring at different positions. When performing a Fourier transform in space, the effect of the phase mixing is that of shifting the peak of the perturbation in time, towards higher values of the wavenumber. Therefore, the amplitude of the perturbation in Fourier space calculated at the initial wavenumber, decays in time, giving rise to the continuum damping. Gyrokinetic simulations of GAMs with realistic profiles of ASDEX Upgrade had already shown a peculiar behaviour at the edge, which was not been possible to explain with a description based on the Landau damping alone [27]. In particular the damping rate at the edge, where the temperature profile is steep, was measured to be much stronger than that predicted by the theory of Landau damping. More recently, the investigation of this strong damping observed in gyrokinetic simulations has started, showing a new mechanism based upon the combined action of the phase mixing and the Landau damping, and the first results have been shown in Ref. [28].
In this work, we complete the level of understanding of this problem, by describing the details of the combined effect of the nonuniformities in velocity and in realspace, causing the cascade in the wavenumber and the effective damping at high wavenumbers, with a dedicated initialvalue code and diagnostics. In fact, as previously mentioned, the phase mixing, acting because of the nonuniformity in space, damps only each single mode in Fourier space by continuum damping, but in the whole real space the energy of the electric field is conserved, so the actual total damping is zero. Therefore, although our analysis is linear, we can state that the values of the damping rate calculated as pure Landau damping and as pure phase mixing do not just linearly sum up. In fact, the Landau damping of GAMs strongly increases with the radial wave number, due to the finite orbit width of the passing ions. Consequently the effect of the phase mixing is to greatly amplify the efficiency of the Landau damping. This general effect, recently proved to be crucial (see Ref. [28]), is explored here in specific cases chosen for their practical experimental relevance, and quantitatively proved to be indispensabile for an interpretation of realistic GAM dynamics. In particular, the phase mixing and Landau damping are first studied separately and then the resulting combined effect is shown to form a cascade of the energy in wavenumber (analogously to that generated in a turbulent system) with a strong absorption at the small scales [24, 25, 9], which is reproduced here quantitatively with a dedicated initial value code. The investigation of the Landau damping with analytical theory and numerical simulations with ORB5 was performed already in Ref. [40, 28], and recalled here for sake of completeness. The investigation of the phase mixing, on the other hand, is performed in this work with the help of an initial value code which solves the GAM dynamics by adopting the real part of the kinetic frequency, and neglects its imaginary part. This initial value code for the study of the phase mixing, is used to reproduce the evolution of the radial structure of the GAM perturbation, modified by the phase mixing alone.
In tokamaks, the combined effect of Landau damping and phase mixing can play a crucial role in the dynamics and in the decay of GAM at the tokamak edge, where the temperature gradients are very large. Depending on the confinement regime, low confinement (Lmode), improved (low) confinement (Imode) and high confinement (Hmode), different parameter values are achieved in the experiments, in particular changing the relative strength of temperature and density gradients at the edge of the plasma. The sequence of events in the transition phases between these regimes, particularly regarding the turbulence suppression, the development of the edge temperature and density gradients and the increase of the shear flow, is not completely established yet. GAMs are considered to be potential key players in the dynamics of the transition to I and to Hmodes [29], and their absence can enhance the effects of the shear flow on the turbulence [11]. GAMs are reported to be regularly observed in Lmode, and are also observed in Imode, whereas they are not observed in Hmode. It will be shown that the proposed mechanism of GAM decay is consistent with the observed existence or nonexistence of GAMs in the different confinement regimes [30].
Our investigation is carried out by means of analytical theory and numerical simulations. For the numerical simulations, we use the global particleincell code ORB5 [31], which now includes all extensions made in the NEMORB project [32, 33]. The ORB5 code uses a Lagrangian formulation based on the gyrokinetic VlasovMaxwell equations [34, 35, 36, 37, 38]. Due to the method of derivation of the GK VlasovMaxwell equations from a discretized Lagrangian, the symmetry properties of the starting Lagrangian are passed to the VlasovMaxwell equations, and the conservation theorems for the energy and momentum are automatically satisfied [39]. The code solves the fullf gyrokinetic Vlasov equation for ions and the driftkinetic equation for electrons. Only linear collisionless electrostatic simulations are considered in this paper, with electrons treated as adiabatic, and with magnetic flux surfaces with circular concentric poloidal sections.
This paper is structured as follows. In Sec. 2, the Landau damping in uniform plasmas is studied with ORB5 and compared with analytical theory. In Sec. 3, the phase mixing signatures are identified and compared with results of numerical simulations, obtained for realistic values of tokamak plasmas. Sec. 4 is devoted to a description of the combined action of the phase mixing and Landau damping. Finally, Sec. 5 discusses the application to predict the regimes where such damping overcomes the drive given by the turbulence, and the summary of the results.
2 Landau damping in uniform plasmas
In this Section, we recall the results about the Landau damping of GAMs in uniform plasmas, i.e. in the absence of gradients of the equilibrium profiles (i.e. the safety factor, density and temperature profiles are flat), shown in Ref. [28]. The GAM dynamics is independent on the radius to a good approximation, and the radial shape of the initial perturbation does not change in time, at least for the first few GAM oscillations (after which, global effects can occur). For our regimes of interest, a linear estimation of the frequency and damping rate is given by the analytical theory where firstorder ion finiteorbitwidth effects are retained [16, 17]. These simulations serve as a verification of ORB5, extending the work previously done in Ref. [40]. In particular, we study the dependence of the damping rate on the radial wavenumber , for different values of the safety factor .
We initialize a zonal perturbation of the radial component of the electric field, i.e. a perturbation independent of poloidal and toroidal angles, and depending only on the radius. In general in this paper, except where otherwise stated, we consider initial perturbations with a sinusoidal radial dependence, i.e. monochromatic in Fourier space. The perturbation evolves in time in a linear electrostatic collisionless simulation with adiabatic electrons. For each simulation, we investigate the evolution in time of the radial electric field, and in particular the values frequency and damping rate at each radial location, and the radial structure at each time. This type of numerical experiment is known as a RosenbluthHinton test [2].
Different simulations are investigated, with wavenumber chosen within the range 0.02 0.2, and safety factor within the range 1 q 3.5. The value of the plasma temperature is defined by . The ion Larmor radius is defined as , with and with being the ion/electron temperature calculated in the middle of the radial domain (in this Section, ). The magnetic field is calculated at the magnetic axis (r=0). We choose a tokamak with an inverse aspect ratio = 0.1. We consider analytical magnetic equilibria with circular flux surfaces. In this limit, the flux surface coordinate is a good approximation of the radial coordinate (normalized to the minor radius ).
The reference simulation has a spatial grid of () = () and a time step of 100 , with being the ion cyclotron frequency. All the simulations have been performed with ion markers. As a general remark, we find very good agreement between theory and simulations, consistently with Ref. [40]. In particular, our results with ORB5 are consistent with the analytical theory, showing that, at realistic values of q, the damping rate strongly increases with (see Ref. [16, 17, 28])).
3 Phase mixing in the presence of a radial nonuniformity
In this Section, we investigate the effect of a radial nonuniformity of the equilibrium, in generating the phase mixing. RonsebluthHinton tests similar to those described in the previous Section are performed, but with the main difference that the equilibrium radial profiles are not considered flat. This introduces a radial dependence in the GAM dynamics, for example of the GAM local frequency: . In this case, the GAM dynamics in its simplest form (no Landau damping considered here) can be described by the vorticity equation [9, 26]:
(1) 
This equation has been shown to describe the phase mixing [21, 22, 23, 24, 25], i.e. the cascade of energy from bigger to smaller spatial scales. We can recover the phase mixing by integrating Eq. 1 in space once, which yields the radial component of the electric field , and approximating for simplicity the frequency profile as linear in space: . Now it is sufficient to calculate the Fourier transform of the electric field in a radial region centered at , and with halfwidth , to obtain [28]:
(2) 
which can be written, in the limit of radially localized perturbation, as:
(3) 
As a first signature of the phase mixing, we observe the cascade of energy in Fourier space, in the form of a linear dependence of the wavenumber on time:
(4) 
As a second signature, we observe that the scalar potential dependence on time can be now estimated from , obtaining for its amplitude the characteristic decay which takes the name of continuum damping:
(5) 
Before proceeding further, we note that locally the phasemixing does not dissipate the energy of the GAM, being the local total energy (ExB + thermodynamic) equal to , with being the mass density, and being the ExB drift. The investigation of the global properties of GAMs, including the radial propagation of energy, are out of the scope of this paper, and will be discussed separately.
We now want to investigate the signatures of the phase mixing in the results of the numerical simulations in an equilibrium obtained from realistic values of parameters. To this purpose, we have selected a group of parameters defining the equilibrium of ASDEX Upgrade for the shot number 20787 [41], at the radial location where the GAM is detected (r 0.9 a), consistently with Ref. [27, 28]. These parameters will be used for the simulations described in the rest of the paper. The magnetic equilibrium is characterized by a magnetic field on axis of T, the major and minor radii of m and m, and a safety factor of . The plasma has a local deuterium and electron temperatures of eV. We want to investigate the results of several simulations, with different temperature profiles with temperature gradients in the range at r=0.5 (density and safety factor profiles are kept flat for simplicity). The GAM characteristic radial size is estimated as cm. Starting from these values we obtain the following quantities: an ion cyclotron frequency rad/s, an ion thermal velocity m/s and a Larmor radius m (with ). We initialize an electric field with sinusoidal profile with , equivalent to . We note that the GAM frequency depends strongly on the plasma temperature (for the order of magnitude we have , whereas does not depend on and depends weakly on in the considered regime of parameters), therefore we start by considering a simplified problem where the density and safety factor profiles are considered flat and the temperature profile has a radial dependence described as (see Fig. 1):
(6) 
where is the reference position in the simulation where the peak of the temperature gradient is located, and where the measurements are performed, is the value of the temperture gradient, and is the radial size of the temperature gradient peak (in this paper, electron and ion temperatures profiles are always considered equal).
The results of these simulations are investigated, with a particular attention in this paper to the evolution of the radial shape of the perturbation (see Fig. 2). In order to isolate the physics of the phase mixing, we repeat the simulation with a dedicated initial value code, which evolves the electric field at each position with its local frequency:
(7) 
where the frequency is given at each position by the analytical theory, and depends on the local temperature value. No Landau damping is considered in this code. This initial value code gives the evolution in time of the radial structure of the perturbation, due to the pure phasemixing. Effects of the Landau damping in modifying the radial structure are neglected in this code. The comparison of the results of ORB5 and the phasemixing code are shown in Fig. 2. The first feature of the radial structure evolution, observed in both codes, is the increase of the wavenumber in time. This can be observed to occur in both ORB5 and in the phasemixing code, with similar values. We conclude that the wavenumber calculated in gyrokinetic simulation is approximated to a good extent by the phasemixing code. As a second remark, we note that the results of ORB5 show a damping in time at r=0.5, whereas the phasemixing code does not show it, for construction. This damping will be investigated in the next Section. As a third remark, we note that the results of ORB5 show a slight displacement in time of the position of the peak near r=0.5. This secondorder global effect due to the radial GAM propagation, which is not present in the phasemixing code, will be investigated in a separate paper.
As time goes on, we measure the wavenumbers at the center of the radial domain, for both results of ORB5 and of the phasemixing code, and plot them in time. The result can be seen in Fig. 3, where the value of the analytical prediction given by the local approximation, Eq. 4, is also shown as a black line. One can see that the wavenumber measured in the numerical simulations is found to grow more slowly than the local analytical prediction on average. This is due to the local approximation, which does not consider the finite size of the temperature gradient peak. In fact the GAM perturbation experiences a lower temperature gradient away from the reference position r=0.5, and therefore the global effect is an average of the temperature gradient profile, rather than the value measured at r=0.5 only. In Fig. 3, the analytical prediction taking into account also this average to the first order is shown. One can see that the behaviour of the wavenumber measured with the numerical simulations is described for a longer time by the analytical prediction where the global effects are considered to the first order, and eventually higher order effects should be also taken into account. In conclusion, we observe that the phase mixing affects the numerical simulation and we can now extend this result to the combined effect of Landau damping and phase mixing.
4 Combined effect of phase mixing and Landau damping
We are now ready to investigate the combined effects of the phase mixing and Landau damping, which have been studied separately in the previous Sections. The same simulations as in Sec. 3 are considered, where the realistic conditions at the tokamak edge of AUG shot 20787 taken in Ref. [41], are reproduced at the center of the radial domain of the simulations. As shown in Fig. 3, the radial shape of the electric field perturbation evolves in time with higher and higher wavenumber, consistently with Ref. [9, 26]. In the same Figure, we can observe also that the value of the electric field at r=0.5 is lower and lower at each oscillation. In this Section, we want to investigate this damping with a proper theoretical model [28].
The values of q and allow us to treat the Landau damping with the analytical theory with firstorder FOW effects (i.e. with a quadratic dependence on ). We can then write the damping rate due to the combined action of phase mixing and Landau damping (PL) as a function of time:
(8) 
where f and g are the cofficients given in Ref. [16, 17]. Consequently, the evolution of the electric field envelop is described at each time by:
(9) 
Equation 9 is solved with an initialvalue code in order to calculate theoretically the evolution in time of the GAM electric field envelop. The frequency is also known from Ref. [16, 17], and does not change in time to a good approximation. The comparison of the analytical prediction with the result of ORB5 is shown in Fig. 4. We note that the GAM decay is described to a good extent by the PL model, whereas the Landau damping alone greatly overestimates the GAM amplitude at large times. This is due to the fact that the damping rate, which is the Landau damping corresponding to the initial wavenumber at t=0, grows quadratically in time due to the evolution of the wavenumber (see Fig. 4). We also note that some discrepancy is present, due to global effects, such as the initial formation of the GAM eigenmode during the first oscillation, and the radial displacement of the GAM perturbation at later times. The analysis of these corrections due to radial propagation will be done in a separate paper.
We now investigate the dependence of the GAM decay on the value of the temperature gradient, with both gyrokinetic simulations and the analytical theory of the PL mechanism. In particular, the evolution in time of the electric field is measured at r=0.5 in gyrokinetic simulations with different values of ranging from . The time of halfdecay is measured for each simulation and shown in Fig. 5, where it is compared with the analytical model of the PL mechanism. We observe that the decay time strongly decreases with the increase of the temperature gradient. The theory for the evolution of the envelop of the electric field, discussed in Sec. 4, describes well the strong damping with respect to the case = 0 in which only Landau damping acts. We note that for large values of , the halfdecay time in the simulations is slightly larger than the damping expected by the theory, which is due to the radial propagation not investigated in this paper.
5 Summary and conclusions.
Understanding the linear damping mechanisms of GAMs is a crucial step for reaching a complete understanding of their dynamics, and in particular of their nonlinear interaction with turbulence in tokamaks. Landau damping has been historically recognized as the main linear damping mechanism of GAMs, due to their low frequency (for instance with respect to other higher frequency, less damped modes, like incompressible shearAlfvén modes). Landau damping has been shown to depend on the radial shape of GAMs, in particular increasing with the local radial wavenumber. Moreover, the phase mixing has been predicted to affect GAMs, by means of analytical theory [9, 26]. The main signature of phase mixing has been shown to be a cascade of the GAM energy from lower to higher values of the wavenumber. In previous works [27, 28], the results of gyrokinetic simulations in the presence of a temperature gradient have been investigated, and in Ref. [28] the combined effect of the phase mixing and the Landau damping has been shown to generate a novel damping mechanism, responsible for a very efficient absorbption of the GAM energy by the bulk plasma [28]. This phase mixing / Landau damping (PL) mechanism, has been shown to increase up to an order of magnitude the effective damping of GAMs, for realistic tokamak conditions.
In this work, we have investigated the details of the radial structure evolution in time of the GAM electric field in the presence of a temperature gradient. The cascade of the wavenumber in time has been investigated with gyrokinetic simulations and with a specific initial value code dedicated to the investigation of the phasemixing. The radial structure of gyrokinetic simulations obtained with ORB5 has been shown to be approximated to a good extent by the one obtained with the initial value code. Second order effects, like the radial displacement of the GAM peak in time, have been neglected at this stage and have been left for a future work. Finally, consistently with Ref. [28], the combined effect of the phasemixing and Landau damping has been investigated, and shown to decrease the halfdecay time up to one order of magnitude, for realistic tokamak parameters, and using new simulations where the temperature gradient is peaked at the center of the radial domain.
From these results, we deduce that the PL mechanism can play an important role in the suppression of GAM oscillations in the regions characterized by a strongly nonuniform temperature profile such as in the H and I modes. We note that simulations performed with a density gradient different from zero have given results very close to the simulations performed with a flat density profile. This is in agreement with Eq. 8 that does not depends on the density gradient. Moreover, we note that in Eq. 8 depends also on q. However, the gradient of q has a weak influence on the phase mixing and consequently on . This aspect has been verified by numerical simulations showing that the main parameter in the PL mechanism is the temperature gradient.
In order to investigate and to quantify the importance of the PL damping mechanism on these regions, we compare it with the characteristic drive rate given by the nonlinear coupling with the ITG mode in Ref. [9], and discussed in Ref. [28]. The characteristic drive time can be calculated as , and should be compared with the PL decay time. If , this means that the PL damping rate exceeds the energy transfer rate from the ITG turbulence to the GAM, and as a consequence we may expect no observation of GAMs. The drive time is found to be for the Lmode, for the Imode, and for the Hmode, with being the sound time unit. Therefore, we observe that the PL model explains the observation of GAMs for L and I modes, but a desappearance is predicted for the Hmode, for (see Ref. [28] for the detailed calculation). We note that this analysis of orders of magnitude gives results which are consistent with experimental results that show the existence of GAM in Imode regime in spite of the strong temperature gradient comparable to that of Hmode [29, 30]. These results are also in agreement with the dynamics of GAM observed at the IH transition [29].
As a future work, a more detailed description of a particular experimental shot will be done by using the microturbulence features observed in the experiment of AUG and comparing with linear and nonlinear gyrokinetic simulations where global effects are also investigated.
Appendix A Numerical parameters and power balance.
The numerical simulations described in this paper, have been performed with the global gyrokinetic particleincell code ORB5, by solving electrostatic gyrokinetic equations for the ions and treating the electrons as adiabatic. Only linear collisionless simulations have been performed. A typical simulations has a spatial grid of = 1024 x 64 x 4 and a time step of 50 . The number of ion markers has been chosen as . A set of poloidal harmonics going from m=10 to m=10 has been selected by means a filter in Fourier space in the poloidal direction. Dirichlet boundary conditions on the potentials have been imposed at the outer boundary, and Neumann boundary conditions at the inner boundary.
In a collisionless gyrokinetic simulation, the GAM amplitude is damped in time due to the combined effect of phase mixing and Landau damping, where the phase mixing creates a cascade in the radial wavenumber and the Landau damping effectively damps the higher wave numbers [28]. Ultimately, the energy of the GAM goes into the bulk ion microscopic kinetic energy due to the Landau damping. In this Section, this energy channel among the global GAM mode and the microscopic ion energy is investigated.
In this model the rate of change of the particle kinetic energy must be equal to the power transfer from the particles to the field (power balance) (see Ref. [32]). The power balance not only gives an indication of the quality of the simulations, but also provides a measure of the instantaneous growth (decay) rate. Defining the rate of variation of kinetic energy and the power transfer, the relative error of the power balance is:
(10) 
where is the maximum value of in the time interval of interest. In Fig. 6, the relative error of two simulations is shown, namely one with flat temperature profile, and one with . For this specific runs, that falls within 3% at the time of half decay.
Acknowledgements
Interesting and useful discussions with G. Conway, E. Poli, P. Manz, Z. Qiu, B. Scott and ASDEX Upgrade team are kindly acknowledged. This work has been done in the framework of the European Enabling Research Project on “Verification and development of new algorithms for gyrokinetic codes“, WP15ER01/IPP01, and the European Enabling Research Project on “Microturbulence properties in the core of tokamak plasmas: close comparison between experimental observations and theoretical predictions”, WP15ER01/IPP02. Simulations were performed on the IFERCCSC Helios supercomputer within the framework of the VERIGYRO and the ORBFAST project.
References
 [1] A. Hasegawa, C.G. Maclennan, and Y. Kodama, Phys. Fluids 22, 2122 (1979)
 [2] M.N. Rosenbluth and F.L. Hinton, Phys. Rev. Lett. 80,4 724 (1998)
 [3] P.H. Diamond, S.I. Itoh, K. Itoh and T.S. Hahm, Plasma Phys. Controlled Fusion 47, R35 (2005)
 [4] F. Palermo, X. Garbet, A. Ghizzo, T. CartierMichaud, P. Ghendrih, V. Grandgirard and Y. Sarazin, Phys. Plasmas 22, 042304 (2015).
 [5] F. Palermo, X. Garbet, and A. Ghizzo, Europ. Phys. J. 22, 082304 (2015).
 [6] A. Ghizzo and F. Palermo, Phys. Plasmas 69, 1 (2015).
 [7] A. Ghizzo and F. Palermo, Phys. Plasmas 22, 082304 (2015).
 [8] N. Winsor, J.L. Johnson and J.M. Dawson, Phys. Fluids 11, 2448, (1968)
 [9] F. Zonca and L. Chen, Europhys. Lett. 83, 35001 (2008).
 [10] N. Miyato, Y. Kishimoto and J. Li, Phys. Plasmas 11, 5557 (2004).
 [11] P. Angelino, A. Bottino, R. Hatzky, S. Jolliet, O. Sauter, T.M. Tran and L. Villard, Plasma Phys. Controlled Fusion 48, 557 (2006).
 [12] B. Scott, Plasma Phys. Controlled Fusion 34, 12A, 1977 (1992)
 [13] B. Scott, New Journal of Phys. 7, 92 (2005)
 [14] F. Zonca, Liu Chen and R.A. Santoro Plasma Ph. Control. Fus. 38, 20112028 (1996)
 [15] F. Zonca, L. Chen, R.A. Santoro and J.Q. Dong Plasma Phys. Control. Fus. 40, 2009 (1998)
 [16] H. Sugama and T.H. Watanabe, J. Plasma Physics 72, 825 (2006)
 [17] H. Sugama and T.H. Watanabe, J. Plasma Physics 74, 139 (2008)
 [18] X. Q. Xu, Z. Xiong, Z. Gao, W. M. Nevins, and G. R. McKee, Phys. Rev. Lett. 100, 215001 (2008)
 [19] Z. Qiu, L. Chen, and F. Zonca, Plasma Phys. Controlled Fusion 51, 012001 (2009).
 [20] Zhe Gao, Phys. Plasmas 17, 092503 (2010)
 [21] K. M. Case, Phys. Fluids 3 2, 143 (1960)
 [22] Z. Sedláček, J. Plasma Phys. 5, 239 (1971)
 [23] H. Grad, Phys. Today 22, (12) 34 (1969)
 [24] A. Hasegawa and L. Chen, Phys. Rev. Lett. 32, 454 (1974)
 [25] L. Chen and A. Hasegawa, Phys. Fluids 17, 1399 (1974)
 [26] Z. Qiu, F. Zonca and L. Chen Plasma Sci. Technol. 13, 3 (2011).

[27]
A. Biancalani, A. Bottino, Ph.W. Lauber, 40th European Physical
Society Conference on Plasma Physics, Espoo, Finland (2013), P1.157
http://ocs.ciemat.es/EPS2013PAP/pdf/P1.157.pdf  [28] F. Palermo, A. Biancalani, C. Angioni, F. Zonca and A. Bottino, Europhys. Letters in press (2016)
 [29] I. Cziegler, P.H. Diamond, N. Fedorczak, P. Manz, G. R. Tynan, M. Xu, R.M. Churchill2, A. E. Hubbard, B. Lipschultz, J.M. Sierchio, J.L. Terry and C. Theiler, Phys. Plasmas 20, 055904 (2013).
 [30] G. D. Conway, C. Angioni, F. Ryter, P. Sauter, and J. Vicente Phys. Rev. Lett. 106, 065001 (2011).
 [31] S. Jolliet, A. Bottino, P. Angelino, R. Hatzky, T.M. Tran, B.F. Mcmillan, O. Sauter, K. Appert, Y. Idomura and L. Villard Comput. Phys 177, 409 (2007)
 [32] A. Bottino, E. Sonnendrücker J. Plasma Phys. 81, 435810501 (2015).
 [33] A. Biancalani, A. Bottino, A. Mishchenko, B. McMillan, and L. Villard, “Adapting a gyrokinetic global code for Alfvén modes”, to be submitted to Comput. Phys. Commun.
 [34] T. S. Hahm, W. W. Lee and A. Brizard, Phys. Fluids 31, 1940 (1988)
 [35] H. Sugama, Phys. Plasmas 7, 466 (2000)
 [36] A. J. Brizard and T. S. Hahm, Rev. of Modern Phys. 79 421 (2007)
 [37] A. Bottino, T. Vernay, B. Scott, S. Brunner, R. Hatzky, S. Jolliet, B.F. McMillan, T.M. Tran and L. Villard, Plasma Phys. Controlled Fusion 53, 124027 (2011)
 [38] B. D. Scott, and J. Smirnov Phys. Plasmas 17, 112302 (2010)
 [39] B. D. Scott Phys. Lett. A 320, 53 (2003).
 [40] A. Biancalani, A. Bottino, Ph. Lauber and D. Zarzoso, Nucl. Fusion 54 104004 (2014)
 [41] G. Conway, C. Tröster, B. Scott, K. Hallatschek and the ASDEX Upgrade Team, Plasma Phys. Control. Fusion, 50 055009 (2008).
 [42] L. Chen, Z. Lin and R. White, Phys. Plasmas, 7 3129 (2000).
 [43] G. D. Conway, private communication, (2016).
 [44] A. Bottino, A. G. Peeters, O. Sauter, J. Vaclavik, L. Villard, and ASDEX Upgrade, Phys. Plasmas 11, 198 (2004)