Linear gyrokinetic investigation of the geodesic acoustic modes in realistic tokamak configurations
Geodesic acoustic modes (GAMs) are studied by means of the gyrokinetic global particle-in-cell code ORB5. Linear electromagnetic simulations in the low- limit have been performed, in order to separate acoustic and Alfvénic time scales and obtain more accurate measurements. The dependence of the frequency and damping rate on several parameters such as the safety factor, the GAM radial wavenumber and the plasma elongation is studied. All simulations have been performed with kinetic electrons with realistic electron/ion mass ratio. Interpolating formulae for the GAM frequency and damping rate, based on the results of the gyrokinetic simulations, have been derived. Using these expressions, the influence of the temperature gradient on the damping rate is also investigated. Finally, the results are applied to the study of a real discharge of the ASDEX Upgrade tokamak.
The ion heat transport in the plasma core is governed by turbulence formed by a class of microinstabilities such as toroidal ion temperature gradient (ITG) driven modes . ITG turbulence is known to self-organize to form macroscopic structures . These structures take the form of a macroscopic radial electric field which depends only on the radial coordinate. poloidal flows associated with this electric field are referred to as zonal flows (ZFs) [3, 4, 5, 6].
The action of the toroidal magnetic field curvature on the ZF gives rise to oscillations of the radial electric field. These oscillations of the ZFs are called geodesic acoustic modes (GAMs) [7, 8]. The modes are observed predominantly in the edge region of the tokamak plasmas with characteristic frequency of the order of the sound frequency , where is the sound speed, is the major radius. One of the main linear damping mechanisms for the stationary ZF are collisional processes and for the GAM it is a collisionless wave-particle interaction, namely the Landau damping, and collisional damping at the very edge of the plasma, where equilibrium temperatures drastically decrease . A recent comparison of collisionless and collisional damping of GAMs, using existing analytical theories, for experimentally relevant plasmas was done in Ref. .
The importance of the ZF is that they can regulate the drift-wave (DW) turbulence . But it is still a question how the GAMs influence the ZF efficiency of the DW suppression [8, 12, 13]. On the other hand, the development of zonal structures can play a key role in the transition from the low to the high confinement regime (L-H transition) . In Ref.  interaction of the mean and oscillatory poloidal flows with the turbulence were experimentally observed. The turbulence suppression by the ZFs was observed in experiments described in Ref. . On the other hand, in Ref.  the role of the mean flow in the dynamic evolution towards the H-mode is emphasized. In Ref.  two predators - one prey system, including ZF, GAM and turbulence, was developed to study transitions between states with different combinations of the ZF and GAM.
In this paper, we investigate the GAM frequency and collisionless damping rate, carrying out linear collisionless simulations with kinetic electrons. The electromagnetic global gyrokinetic particle-in-cell code ORB5 is used [19, 20]. As it has been reported previously [21, 22], models, numerical or analytical, derived with adiabatic electrons, result in considerably smaller GAM damping rate in comparison to simulations performed with kinetic electrons. By adiabatic electron models, we mean here models treating the component of the electrons as adiabatic, and setting the zonal component of the electron density perturbation to zero. In simulations considered in this paper, electrondocumentclasss are treated drift-kinetically, and a realistic ion-electron mass ratio is used. Moreover, to study the influence of the plasma elongation on the GAM dynamics, magnetic equilibria with realistic plasma shapes are considered. To summarize the results obtained in different plasma regimes, interpolating formulae for the GAM frequency and damping rate, based on the gyrokinetic simulations with ORB5, are derived.
Due to the so-called phase mixing effect, the GAM damping rate is increased in the presence of a temperature gradient or the safety factor profile [8, 23, 24]. This effect arises when the damping rate of the wave depends on its wavenumber. In the case of the GAM the damping rate increases with the GAM radial wavenumber (more precisely, with the radial wavenumber of the radial electric field). Since the GAM frequency depends on the temperature and safety factor, the GAM oscillates with different frequencies at different radial points in presence of the temperature gradient or magnetic shear. Distorting the GAM radial structure and creating higher radial wavenumbers, this process can strongly increase the GAM damping rate [8, 25]. A section 4 of our paper is dedicated to the extension of previous works [22, 24], which were done treating the electrons as adiabatic, and in circular flux surfaces, to the inclusion of kinetic electrons and realistic tokamak configurations. Finally, the last section of this paper is dedicated to the investigation of a realistic discharge of ASDEX Upgrade, described in Ref. . In Appendix B we have shown a comparison between ORB5 and GENE for the case of non-flat temperature profile.
The gyrokinetic simulations presented in this work have been performed with the code ORB5 [19, 20]. ORB5 is a nonlinear gyrokinetic multi-species global particle-in-cell (PIC) code, which solves the Vlasov-Maxwell system in the electrostatic or electromagnetic limit, and has a capability of handling true MHD equilibrium for an axisymmetric toroidal plasma. The particle-in-cell method consists of coupling a particle-based algorithm for the Vlasov equation with a grid-based method for the computation of the self-consistent electromagnetic fields. Several physical models are available in ORB5, all of them derived from a systematic Hamiltonian theory [20, 27] to provide exact energy and momentum conservation. In this work, only one ion species (deuterium) has been considered while the electrons are assumed to be drift-kinetic. This corresponds to the following gyrokinetic total Lagrangian:
The velocity variables are the magnetic moment , the canonical parallel momentum and the gyroangle . The equilibrium magnetic field is , and are the mass and charge of the particle species and is the speed of light. The volume element of the velocity space is with , and ; denotes the volume element in physical space. Here is the distribution function for the species , while is the equilibrium time independent distribution function of the ions. In this system, only long wavelength electrostatic perturbation and magnetic perturbations perpendicular to the equilibrium magnetic field are considered. Note that no second order term in the fields is retained for the electrons, this is equivalent to neglect the electron polarization density in the Polarization equation (drift-kinetic approximation, see Ref.  for details). The first two terms in the total Lagrangian define the charged particles Lagrangian . The GK Hamiltonian in general depends on the electrostatic potentials and on the parallel component of the fluctuation magnetic potential . The third term in the total Lagrangian is the electromagnetic field Lagrangian, in which the electric field component has been neglected (quasi-neutrality approximation, see  for details). In this work we used the following Hamiltonian:
the gyroaveraging (Hermitian) operator , applied to an arbitrary function in configuration space, is defined by
where is the vector going from the guiding center position to the particle position.
In this work we have assumed for the electrons (drift-kinetic approximation).
The gyrokinetic equations for the particle distribution function and the GK
field equations can be derived from the GK Lagrangian using variational
principles. In summary, the GK model used in the following is:
gyrokinetic full-f Vlasov equation for the ions
drift-kinetic full-f Vlasov equation for the electrons:
having introduced the generalized potential
Linear polarization equation in the long wave-length limit (and drift-kinetic electrons):
Linear Ampère’s law:
where is the density associated with the equilibrium Maxwellian and . The skin depth is defined by , where , and it appears on the right-hand-side of the Ampére’s law because of the choice of the velocity space variables instead of the usual . The indexes and indicate ions and electrons respectively and is the ion charge while is the electron charge. Despite all the approximations made, this model is highly physically relevant and it can be used to describe not only the GAM and ZF dynamics, but also a large class of micro-instabilities excited by the density and temperature gradients, like ion temperature gradient (ITG) driven modes, trapped electron modes (TEM) or kinetic ballooning modes (KBM). It also contained the reduced MHD model as a subset (see, among other, Ref. ).
According to the PIC method the particle distribution function is discretized with macroparticles, known as markers. The motion of the markers is calculated using the equations of motions of the gyrokinetic model while the electromagnetic fields are evolved on a spatial grid using the two field equations. The charge and current density, that are necessary to solve the field equations, are calculated by projecting the marker weights on a spatial grid. After that, the fields are calculated using a finite elements method.
The code is based on a straight-field-line coordinate system . Here, radial coordinate is (where is the poloidal flux), is the straight-field-line coordinate (where and are the poloidal angle and the safety factor respectively) and is a toroidal angle. Two different kinds of magnetic equilibria are implemented: analytical equilibria with circular concentric magnetic surfaces and ideal MHD realistic equilibria. For the latter case, the ORB5 code is coupled with the CHEASE code , which solves the Grad-Shafranov equation with a fixed plasma boundary.
3 Frequency and Landau damping
3.1 Equilibrium and simulations parameters
Linear electromangetic gyrokinetic collisionless simulations with drift-kinetic electrons and a realistic electron - ion (deuterium) mass ratio have been performed. Electrostatic simulations with kinetic electrons are, in principle, faster than electromagnetic simulations, due to the smaller number of equations to be solved. Nevertheless, a high frequency oscillation, called the -mode , is observed to be often numerically unstable. To decrease the level of the high-frequency oscillations, electromagnetic simulations in the small- () limit have been performed instead of the electrostatic ones. MHD equilibria of the circular and elongated plasma have been calculated with an external code CHEASE . Simulations have been carried out with a flat density profile, which have been shown to not impact the GAM frequency and damping rate in linear simulations (see Appendix A). To focus on the Landau damping in the absence of the phase mixing effect, flat temperature profile has been considered in simulations used for the results of this section. Since the safety factor profiles have been taken from the CHEASE, there is a magnetic shear, that also causes the phase mixing, but its influence on the GAM damping rate is much smaller in comparison to the temperature gradient effect.
Plasma parameters have been taken close to the ASDEX Upgrade parameters near the plasma edge : the major radius m, the minor radius m (inverse aspect ratio is ), the magnetic field on the axis T. Since in the CHEASE code the plasma elongation is defined at the edge and changes gradually to the plasma center, the GAM frequency and damping rate have been measured at the same radial position to perform more accurate scan on the elongation. The temperature has been taken to be eV. It means that m/s, m and . Here, ( MHz) is the ion gyro-frequency, is an electron charge, and is a proton mass.
To weaken the constraint on the space step and to reduce the effect of the charge accumulation at the edge of the numerical work box, we have simulated only a ring from to in a poloidal cross section with the Dirichlet condition for the potential on inner boundaries ().
A typical simulation has the following parameters. Number of nodes in radial direction is taken to be , in toroidal directions and along the straight-field-line coordinate the number of nodes is . Time step is . The GAM damping rate and frequency have been calculated (see Appendix A) for different GAM radial wavenumbers (where m is the ion Larmor radius of the deuterium and is the thermal speed), the safety factor at and the plasma elongation at the edge. This is the regime where GAMs are typically observed in tokamak plasmas (see, for example, Ref. ). To simulate the GAM dynamics, the ORB5 simulations have been initialized by introducing an axisimmetric density perturbation designed to produce an initial electric potential field of the form , where (as in the so-called Rosenbluth-Hinton test ). All toroidal modes and poloidal modes have been filtered out. To study the GAM dynamics, the frequency and damping rate of the poloidally averaged radial electric field have been calculated.
3.2 Results of gyrokinetic simulations
In Fig. 1, a comparison of the GAM frequencies and damping rate obtained from numerical simulations with two analytical theories of Qiu 2009  and Gao 2010  is shown. A good agreement between numerical results and analytical predictions of the GAM frequency has been found. Nevertheless, the GAM damping rate, obtained from the theories, derived using adiabatic electrons, is smaller in comparison to numerical simulations with kinetic electrons, and the divergence increases for smaller values of the GAM radial wavenumber. Moreover, since the frequency stops increasing in the domain of higher wavenumbers (subplot of Fig. 1), a divergence between numerical results and analytical theories is observed. The same effect was observed in Ref. , where the GAMs were studied using drift reduced Braginskii equations.
The Gao 2010 theory describes the GAM dependence on the plasma elongation and it is in a good agreement with numerical results for the frequency. Although the Gao 2010 theory provides considerably smaller damping rate, it seems to give similar trend of the damping coefficient with the plasma elongation, i.e., the damping rate is weakened by the elongation. The Gao 2010 theory was derived in the large orbit drift width limit, where the dominant damping mechanism is the resonance (here, is a magnetic drift frequency, is a wave vector of the zonal potential in the radial direction and is a magnetic drift velocity) . As explained in Ref. , the GAM frequency decreases with the elongation less rapidly than the drift frequency. To satisfy the resonance particles have to have higher drift velocities, which involves fewer particles in the wave-particle interaction and, as a result, the GAM damping rate decreases.
3.3 Interpolating formulae
To provide a scaling of the GAM frequency and damping rate, corresponding interpolating expressions have been fitted to the results of the gyrokinetic simulations described in Sec. 3.2. For consistency, the regime has been chosen for the GAM wavenumbers in the range , safety factor and plasma elongation .
To derive an interpolating expression for the frequency several assumptions have been used. The experimentally obtained dependence  on the plasma elongation has been slightly modified to , where is an adjustable coefficient. The dependence on the safety factor has been taken in the form . In fact, the -dependence in a form of , that is given in Ref. , gives the same results. To describe how the frequency changes with the radial wavenumber, a polynomial has been taken. Moreover, to take into account the frequency saturation for higher wavenumbers we have introduced a function of the form . Here, , . The resulting frequency interpolating formula is the following one:
Among different tested functions, this form gives the best approximation to numerically simulated values of the GAM frequency, it has one of the smallest 95% confidential bounds and is not overfitted. The corresponding coefficients with their 95% confidential bounds (lower and upper bounds) are
For the damping rate we have derived the following expression (here, the damping rate is normalized to ):
with interpolating coefficients
Comparison between the results from the gyrokinetic simulations and the interpolation expression for the GAM damping rate is shown in Fig. 3 for some specific values of parameters taken as examples.
4 Phase mixing
To investigate the influence of the phase mixing on the GAM dynamics the same parameters as described in chapter 3.1 has been used, but the temperature gradient (the same for both the electrons and ions to have , eV) at a radial position has been introduced. The radial point has been chosen here to be in agreement with the section 3.1. Initial radial wavenumber of the radial electric field is . The safety factor is . We consider a temperature profile of the following form, similarly to Ref. :
where , . The temperature profiles and the corresponding temperature gradient profiles for different in a radial interval , are shown in Fig. 4. Dependence of the GAM half-decay time on the temperature gradient has been investigated in the domain .
A scan of gyrokinetic simulations with the temperature gradient has been performed, and the results are depicted in Fig. 5. In presence of a temperature gradient, the GAM is observed to oscillate with different frequencies at different radial points, that leads to the distortion of the initial GAM radial structure. Producing higher radial wavenumbers, this distortion amplifies the GAM damping. This combined effect, already investigated for a more simplified configuration in Ref. [23, 24], has been observed even more pronounced in the simulations described here. In fact, here the phase mixing effect is investigated using gyrokinetic simulations with kinetic electrons that significantly influences the GAM damping, and, as a consequence, the GAM half-decay time. For example, using the Sugama-Watanabe model, which is derived with adiabatic electrons, for the Landau damping and combining with phase mixing, we have obtained for the and for , that predicts much longer half-decay time of the GAM in comparison to the calculations based on the simulations with the kinetic electrons (compare with a Fig. 5).
In order to verify the results of gyrokinetic simulations, we have used a theoretical simplified model of the phase mixing, proposed in Ref. [8, 23, 24], where the linear growth in time of the radial wavenumber is considered.
In the phase mixing simulations a space point is considered with a certain temperature and temperature gradient . Initial radial electric field has the following radial structure:
with an initial amplitude and initial normalized radial wavenumber . The electric field is assumed to evolve in time at a point according to a simple rule
where is an amplitude of the electric field, that changes in time due to the damping, . The general form of the GAM frequency is
where describes frequency dependence on the radial wavenumber, the safety factor and the elongation. The safety factor profile is taken to be flat, and plasma with a circular cross-section is considered: , .
The damping rate is defined as
At the beginning of every time interval , new values of the damping rate and frequency are found with the scaling formulae given in Eq. (5) and (4), using a current value of the wavenumber . A new value of the electric field can be found, assuming that the damping rate is constant at the lapse of time :
After that, a new value of the wavenumber is calculated using the radial derivative of the frequency
With that, the wavenumber is assumed to change linearly in time as
where , . Another option, it is to estimate the time evolution of the radial wavenumber directly from numerical calculations in ORB5. Substituting new value of the normalized wavenumber into Eq. (5), we can find the damping rate at the next time point.
The results obtained with this reduced theoretical model are also shown in Fig. 5. As it can be seen in Fig. 5, the qualitative dependence of the half-decay time on the temperature gradient finds a good match of gyrokinetic simulations of ORB5 and analytical theory. The difference is due to the global dynamics of the ORB5 simulations, which is compared here with a theory where the phase mixing follows a local estimation given in Ref. .
5 Comparison with experimental data
The dispersion relations obtained in Sec. 3.3 as an interpolation of gyrokinetic simulations and given in Eqs. (4), (5) can be used to compare numerical estimations of the GAM behaviour to measurements of the GAM frequency, performed on ASDEX Upgrade tokamak  using Doppler reflectometry. More precisely, we consider the discharge AUG#20787 with the plasma elongation at the edge (and we assume that it is constant at the considered radial region ).
The GAM radial wavenumber is considered to be constant and is estimated to be from the experimental radial profile of the GAM amplitude (see Fig. 5f in Ref. ). Experimental safety factor and ion temperature profiles have been taken to estimate the GAM frequency and damping rate using the scaling formulae (4), (5) at different radial points . In Fig. 6 the GAM frequency profiles with corresponding theoretical prediction are depicted. A good general agreement is found in the central region of interest, where the GAM intensity, measured in the experiments, is peaked. On the other hand, the linear dispersion relation 4 can not explain neither the staircase nature (the plateaus) of the frequencies nor the GAM peak splitting that is observed experimentally at the radius positions or (although the presence of GAM eigenmodes has been suggested by simplified analytical models [37, 38], whose detailed analysis is out of the scope of this paper). For this reason, we can conjecture that the coherent phenomena at the basis of the formation of GAM extended eigenmode or frequency splitting must have a nonlinear origin.
For reference we have given here estimation of the GAM collisional damping rate using formulae, derived by Gao in Ref. . Introducing normalized ion collision rate , the collisional damping rate is calculated as
if , and as:
if . To find the ion collisional rate we have used classical expressions:
where , .
According to the Fig. 6, the collisional damping is found to be negligible in the radial domain where GAMs are experimentally measured, except in a very narrow region close to the separatrix, where it can be of the same order of magnitude as the Landau damping.
In tokamak plasmas, the drift-wave turbulence gives rise to the zonal flows that in their turn shear and distort convective and turbulent cells leading to the saturation of turbulence and, consequently, to a reduction of the radial heat transport. Action of the magnetic curvature results in the oscillatory zonal flows, so-called geodesic acoustic modes. The peculiarity of the GAM oscillations resides in the different shearing efficiency that the ZF have in relation to their oscillatory behavior. The nonlinear interactions between the GAM and the DW turbulence is defined in a high degree by the GAM damping rate. Lack of the experimental data of this characteristic of the GAM makes the results from linear gyrokinetic simulations particularly important for analytical and numerical investigation of the nonlinear GAM-DW systems.
In this work, linear gyrokinetic simulations have been performed with kinetic electrons to study the GAM dynamics. Numerical results have been compared to analytical theories, derived with adiabatic electrons. It has been shown that analytical theories, derived with adiabatic electrons, result in smaller values of the damping rate and for higher wavenumbers diverge from numerical calculations of the frequency. That is why, investigating the GAM dependence on the plasma safety factor, elongation and radial wavenumber, we have found approximating analytic expressions for the frequency and damping rate to predict the GAM behaviour in different plasma regimes. The derived expressions can be used to estimate the GAM linear characteristics used in analytical models of the nonlinear interactions between the GAM and the DW, such as different reduced models [39, 8]. Using these formulae, the phase mixing effect on the damping rate has also been calculated. Based on the gyrokinetic simulations with kinetic electrons, the results have shown smaller half-decay times of the GAM in comparison with the Ref. [23, 24].
The GAM is one of the special features of the I-mode and can be observed in the L-mode [15, 40, 41]. Comparison of the characteristic drive time of the GAM , which is given by the nonlinear coupling with the ion-temperature-mode (ITG)[8, 39], with the GAM half-decay time confirms the results of the Ref. [23, 24]. Indeed, we estimate the GAM drive time to be (where ) in the L-mode, in the I-mode and in the H-mode, according to Ref. . In this case, it can be seen from the Fig. 5 that the GAM half-decay time, which is defined by both the Landau damping and the phase mixing effect, is much higher than the drive-time in the L and I modes, , for all considered values of . This means that the energy transfer rate from the ITG turbulence to the GAM exceeds the Landau-phase mixing damping rate of the GAM. As a result, the GAM can be observed in the L and I modes, but not in the H-mode, where already for . This could be the explanation for the result that the GAM are not observed in the high-confinement mode, as proposed in Ref. [23, 24].
The approximating expressions have showed quite good agreement with experimental data, but estimations of the GAM frequency, obtained from linear gyrokinetic simulations, does not explain the staircase radial profile of the frequency and the GAM peak splitting. The frequency expression describes only the continuum or dispersive mode in contrast to eigenmode. The latter is characterized by the GAM mode frequencies which are predicted to remain constant over a large radial extent, but a significant radial overlap in the frequency radial profile can be observed, that lead to the GAM frequency peak splitting.
Part of this work was done while two of the authors, I. Novikau and A. Biancalani, were visiting LPP-Palaiseau (France), whose team is kindly acknowledged for the hospitality. The authors acknowledge discussions with F. Jenko, L. Villard, X. Garbet and T. Görler.
Appendix A Numerical convergence tests
To calculate the GAM damping rate and frequency, poloidally averaged radial electric field has been fitted, using the Levenberg–Marquardt algorithm , to a function of the form , where , are sought-for damping rate and frequency. Before the fitting it’s necessary to filter the radial electric field to get ride of the high-frequency Alfvén oscillations. In Fig. 7 the fitting is depicted for the case: , , . It is worth to mention that the choice of a time interval, where the fitting is performed, can influence the result damping rate, and it is not so crucial for the GAM frequency calculation. This ambiguity in the choice of the time interval can be explained by the fact that at the beginning of the simulations there are some transient processes that must be excluded from the damping rate measurements. Moreover, with the time, the global effects start to play a significant role, distorting initial radial structure of the radial electric field, that makes the damping rate to be variable in time.
The GAM dynamics (frequency and damping rate) doesn’t depend on the plasma density in linear calculations (see Fig. 8 and 9). But the frequency of the Alfvén waves decreases with the increase of the density, and for high values of the plasma density it becomes difficult to separate acoustic and Alfvénic time scales (see Fig. 9).
The transition from the simulations with the adiabatic electrons to the ones with the kinetic electrons applies additional restrictions on several numerical parameters such as the time step and the number of markers (see Fig. 10).
In projects with adiabatic electrons the normalized time step can be of the order of 20, but in case of the kinetic electrons it has to be significantly reduced till 2 because of the high parallel velocity of the passing electrons. Also electrostatic simulations of the kinetic electrons reveal high-frequency oscillations. These oscillations can lead to numerical instabilities in case of low number of markers. To reduce their level we have passed to electromagnetic simulations with small values of electron beta that gave us an opportunity to keep the number of markers on the level of .
The radial space step (or number of the points in the radial space grid) is determined by, among other parameters, the GAM wavenumber. To investigate the GAM dynamics with higher values of the radial wavenumber, we simulated a narrow poloidal ring near the edge instead of the full plasma cross-section to reduce the number of radial space points.
Appendix B Comparison with the code GENE
A complete cross-code verification between the gyrokinetic ORB5 and GENE [43, 44] codes has already been done in Ref.  on the linear collisionless dynamics of the GAMs with adiabatic and kinetic electrons in the specific case of flat temperature profiles. For completeness, in this paper a comparison between these different codes is shown including the additional phase mixing physical effect, which is driven by non-flat temperature profiles. The motivation behind this study is that although the linear physical models between ORB5 and GENE are equivalent , the numerical schemes are different. GENE is an Eulerian code, where the distribution function is not discretized with markers, but it is discretized on a 5D fixed grid in phase-space , where is the gyrocenter position, is the parallel velocity, and is the magnetic momentum.
The simulation plasma parameters have been taken as in Sec. 4 for both GENE and ORB5. A sinusoidal perturbation in the potential field is initialised, as defined in Sec. 3.1 and is let evolved in time. In GENE, the radial box size is . We have used 128 grid points in radial direction in order to have at least two points per ion Larmor radius. Along the field line 68 points have been used. In velocity space, 68 points and 128 equidistant symmetric grid points have been used for resolving respectively the and the space. The velocity space domain has been fixed to 3 and 9 times the thermal velocity, respectively in the and space. In order to avoid any recurrence problem, an hyperdiffusivity scheme has been used in the direction. In Fig. 11 peaks of the flux-surface averaged radial electric field, measured at the radial position for the ion-electron temperature gradient , are shown for both GENE and ORB5. Half-decay time, calculated in ORB5, is and for the case of GENE it is .
-  Horton, W. (1999) Rev. Mod. Phys, 71, 735.
-  Manz, P., Ramisch, M., and Stroth, U. (2009) Phys. Rev. Lett., 103, 165004.
-  Hasegawa, A., Maclennan, C. G., and Kodama, Y. (1979) Phys. Fluids, 22, 2122.
-  Rosenbluth, M. N. and Hinton, F. L. (1998) Phys. Rev. Lett., 80, 724.
-  Diamond, P. H., Itoh, S. I., Itoh, K., and Hahm, T. S. (2005) Plasma Phys. Controlled Fusion, 47, R35.
-  Gürcan, Ö. D. and Diamond, P. H. (2015) J. Phys. A: Math. Theor., 48, 293001.
-  Winsor, N., Johnson, J. L., and Dawson, J. M. (1968) Phys. Fluids, 11, 2448.
-  Zonca, F. and Chen, L. (2008) Europhys. Lett., 83, 35001.
-  Gao, Z. (2013) Phys. Plasmas, 20, 032501.
-  Simon, P., Conway, G. D., Stroth, U., Biancalani, A., Palermo, F., and the ASDEX Upgrade Team (2016) Plasma Phys. Control. Fusion, 58, 045029.
-  Biglari, H., Diamond, P. H., and Terry, P. W. (1990) Phys. Fluids B, 2, 1.
-  Miki, K. and Diamond, P. H. (2011) Nucl. Fusion, 51, 103003.
-  Scott, B. (2003) Phys. Lett., 320, 53.
-  P. Manz et al. (2012) Phys. Plasmas, 19, 072311.
-  Conway, G. D., Angioni, C., Ryter, F., Sauter, P., Vicente, J., and the ASDEX Upgrade Team (2011) Phys. Rev. Lett., 106, 065001.
-  Schmitz, L., Zeng, L., Rhodes, T. L., Hillesheim, J. C., Doyle, E. J., Groebner, R. J., Peebles, W. A., Burrell, K. H., and Wang, G. (2012) Phys. Rev. Lett., 108, 155002.
-  Cavedon, M., Pütterich, T., Viezzer, E., Birkenmeier, G., Happel, T., Laggner, F. M., Manz, P., Ryter, F., Stroth, U., and the ASDEX Upgrade Team (2017) Nucl. Fusion, 57, 014002.
-  Miki, K. and Diamond, P. H. (2011) Nucl. Fusion, 51, 103003.
-  Jolliet, S., Bottino, A., Angelino, P., Hatzky, R., Tran, T. M., Mcmillan, B. F., Sauter, O., Appert, K., Idomura, Y., and Villard, L. (2007) Comput. Phys., 177, 409.
-  Bottino, A. and Sonnendrücker, E. (2015) J. Plasma Phys., 81, 435810501.
-  Zhang, H. S. and Lin, Z. (2010) Phys. Plasmas, 17, 072502.
-  A. Biancalani et al. (2017) Phys. Plasmas, 24, 062512.
-  Palermo, F., Biancalani, A., Angioni, C., Zonca, F., and Bottino, A. (2016) Europhys. Lett., 115, 15001.
-  Biancalani, A., Palermo, F., Angioni, C., Bottino, A., and Zonca, F. (2016) Phys. Plasmas, 23, 112115.
-  Chen, L. and Zonca, F. (1995) Phys. Scr., T60, 81.
-  Conway, G. D., Tröster, C., Scott, B., Hallatschek, K., and the ASDEX Upgrade Team (2008) Plasma Phys. Control. Fusion, 50, 055009.
-  Tronko, N., Bottino, A., and Sonnendrücker, E. (2016) Phys. Plasmas, 23, 112115.
-  Sugama, H. (2000) Phys. Plasmas, 7, 466.
-  Miyato, M., Scott, B., and Strintzi, D. (2009) J. Phys. Soc. Japan, 78, 104501.
-  Lütjens, H., Bondeson, A., and Sauter, O. (1996) Comp. Phys. Comm., 97, 219.
-  Lee, W. W. (1987) J. Comput. Phys., 72, 243.
-  Gao, Z. (2010) Phys. Plasmas, 17, 092503.
-  Qiu, Z., Chen, L., and Zonca, F. (2009) Plasma Phys. Controlled Fusions, 51, 012001.
-  Rameswar Singh and Gürcan, Ö. D. (2017) Phys. Plasmas, 24, 022507.
-  Gao, Z. (2011) Plasma Sci. Technol., 13, 15.
-  Sugama, H. and Watanabe, T. H. (2008) J. Plasma Phys., 74, 139.
-  Sasaki, M., Itoh, K., Ejiri, A., and Takase, Y. (2008) Contrib. Plasma Phys., 48, 68.
-  Ilgisonis, V. I., Khalzov, I. V., Lakhin, V. P., Smolyakov, A. I., and Sorokina, E. A. (2014) Plasma Phys. Control. Fusion, 56, 035001.
-  Chen, L., Lin, Z., and White, R. (2000) Phys. Plasmas, 7, 3129.
-  P. Manz et al. (2015) Nucl. Fusion, 55, 083004.
-  G. Wang et al. (2013) Phys. Plasmas, 20, 092501.
-  Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. (2002) Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, .
-  Jenko, F., Dorland, W., Kotschenreuther, M., and Rogers, B. N. (2000) Phys. Plasmas, 7, 1904.
-  Görler, T., Lapillonne, X., Brunner, S., Dannert, T., Jenko, F., Merz, F., and Told, D. (2011) J. Comput. Phys., 230, 7053.
-  Tronko, N., Bottino, A., Görler, T., Sonnendrücker, E., Told, D., and Villard, L. (2017) Phys. Plasmas, 24, 056115.