The nonlinear dynamics of energetic-particle (EP) driven geodesic acoustic modes (EGAM) is investigated here. A numerical analysis with the global gyrokinetic particle-in-cell code ORB5 is performed, and the results are interpreted with the analytical theory, in close comparison with the theory of the beam-plasma instability. Only axisymmetric modes are considered, with a nonlinear dynamics determined by wave-particle interaction. Quadratic scalings of the saturated electric field with respect to the linear growth rate are found for the case of interest. The EP bounce frequency is calculated as a function of the EGAM frequency, and shown not to depend on the value of the bulk temperature. Near the saturation, we observe a transition from adiabatic to non-adiabatic dynamics, i.e., the frequency chirping rate becomes comparable to the resonant EP bounce frequency. The numerical analysis is performed here with electrostatic simulations with circular flux surfaces, and kinetic effects of the electrons are neglected.
Saturation of energetic-particle-driven geodesic acoustic modes due to wave-particle nonlinearity
A. Biancalani, I. Chavdarovski, Z. Qiu, A. Bottino, D. Del Sarto, A. Ghizzo, Ö Gürcan, P. Morel, I. Novikau.
Max-Planck Institute for Plasma Physics, 85748 Garching, Germany
Institute for Fusion Theory and Simulation and Department of Physics, Zhejiang University, 310027 Hangzhou, People’s Republic of China
Institut Jean Lamour- UMR 7198, University of Lorraine, BP 239 F-54506 Vandoeuvre les Nancy, France
LPP, CNRS, Ècole polytechnique, UPMC Univ Paris 06, Univ. Paris-Sud, Observatoire de Paris, Université Paris-Saclay, Sorbonne Universités, PSL Research University, 91128 Palaiseau, France
contact of main author: www2.ipp.mpg.de/~biancala
Two main issues related to magnetic confinement fusion are the turbulent transport, and the dynamics of energetic particles (EP), produced by fusion reactions, or injected for heating purposes. Zonal (i.e. axisymmetric) electric fields are observed to interact with turbulence in tokamaks, in the form of zero-frequency zonal flows (ZF) [1, 2, 3] and finite-frequency geodesic acoustic modes (GAM) [4, 5]. Geodesic acoustic modes, due to their finite frequency, can also interact with EP via inverse Landau damping, leading to EP-driven GAM (EGAM) given that EP drive is sufficient to overcome the threshold condition induced by GAM Landau damping and continuum damping [6, 7, 8, 9, 10, 11, 12]. Understanding the dynamics of EGAMs is crucial due to their interaction with turbulence, which can modify the turbulent transport [13, 14]. Moreover, a strong nonlinear interaction of EGAM with EP is observed in tokamaks [15, 16], leading potentially to a strong redistribution of the EP in phase space.
In order to predict the importance of the interaction of EGAMs with turbulence and with EP, it is important to understand their nonlinear saturation mechanisms. One of the reasons of saturation for EGAMs is the wave-particle nonlinear interaction. In this case, the saturation occurs due to the EP nonlinear redistribution in phase space, and the consequent decrease of the energy exchange between the EP and the EGAM. Another possible reason of saturation is the wave-wave coupling. This can occur for an EGAM interacting with another EGAM and generating zonal side-bands, or for an EGAM with non-zonal modes, e.g., turbulence . In this paper, we focus on the wave-particle nonlinear interactions. For an investigation of the EGAM self-coupling, see Ref. .
The nonlinear saturation of EGAMs is investigated here by means of electrostatic simulations with the gyrokinetic particle-in-cell code ORB5 [18, 19, 20]. ORB5 has been succesfully verified against analytical theory and benchmarked against other gyrokinetic codes, for the linear dynamics of GAMs , and EGAMs [22, 23]. A detailed comparison with the beam-plasma instability (BPI) [24, 25] in a 1-dimensional uniform plasma is also done, following the scheme anticipated in Ref. [8, 26]. Similarly to the BPI, the saturation level of EGAM is shown to scale quadratically with the linear growth rate. A similar investigation was also previously done for Alfvén modes (see, for example, Ref. ). The EGAM frequency is shown to evolve in time when approaching the saturation, like the BPI (see, for example, Ref. [28, 29]). The chirping rate is observed to get of the order of magnitude of the squared of the EP bounce frequency near the saturation. This denotes a transition from adiabatic to non-adiabatic regimes.
The paper is organized as follows. The adopted gyrokinetic model is described in Sec. 2, and the equilibrium and case definition in Sec. 3. The linear dynamics is described in Sec. 4. The saturation levels are investigated in Sec. 5. The regimes of different adiabaticity are investigated by means of the analysis of the frequency, in Sec. 6. Finally, a summary of the results is given in Sec. 7.
2 The model
The main damping mechanism of GAM and EGAM is the Landau damping, which makes the use of a kinetic model necessary. In this paper we use the global gyrokinetic particle-in-cell code ORB5 . ORB5 was originally developed for electrostatic ITG turbulence studies, and recently extended to its multi-species electromagnetic version in the framework of the NEMORB project [19, 20]. In this paper, collisionless electrostatic simulations are considered.
The model equations of the electrostatic version of ORB5 is made by the trajectories of the markers, and by the gyrokinetic Poisson law for the scalar potential . These equations are derived in a Lagrangian formulation , based on the gyrokinetic Vlasov-Maxwell equations of Sugama, Brizard and Hahm [32, 33]. The equations for the marker trajectories (in the electrostatic version of the code) are :
The coordinates used for the phase space are , i.e. respectively the gyrocenter position, canonical parallel momentum and magnetic momentum (with and being the mass and charge of the species). and are respectively the parallel and perpendicular component of the particle velocity. The gyroaverage operator is labeled here by the tilde symbol . , where and are the equilibrium magnetic field and magnetic unitary vector.
Kinetic effects of the electrons are neglected. This is done by calculating the electron gyrocenter density directly from the value of the scalar potential as :
We focus on the dynamics of zonal perturbations, by filtering out all non-zonal components (this is to avoid interactions of zonal/non-zonal modes). Wave-wave coupling is neglected, by evolving the bulk-ion and electron markers along unperturbed trajectories. This means that, in Eqs. 1, 2, 3 for the bulk ions, the last terms, proportional to the EGAM electric field, are dropped. The nonlinear wave-particle dynamics is studied by evolving the EP markers along the trajectories which include perturbed terms associated with the EGAM electric field. This means that the EP markers are evolved with Eqs. 1, 2, 3 where the terms proportional to the EGAM electric field are retained.
The gyrokinetic Poisson equation is :
with being here the total plasma mass density (approximated as the ion mass density). The summation over the species is performed for the bulk ions and for the EP, whereas the electron contribution is given by . Here is the gyrocenter perturbed distribution function, with and being the total and equilibrium (i.e. independent of time, assumed here to be a Maxwellian) gyrocenter distribution functions. The integrals are over the phase space volume, with being the velocity-space infinitesimal volume element.
3 Equilibrium and simulation parameters
The tokamak magnetic equilibrium is defined by a major and minor radii of m and m, a magnetic field on axis of T, a flat safety factor radial profile, with , and circular flux surfaces (with no Grad-Shafranov shift). Flat temperature and density profiles are considered at the equilibrium. The bulk plasma temperature is defined by , with , with being the sound speed and the cyclotron frequency. Three increasing values of bulk plasma temperature are considered, for investigating the dependence of our results on the Landau damping, corresponding to: , , and ( for all cases described in this paper).
All parameters defined so far, are adopted by ORB5 with no consideration of the mass of the bulk ion species. The choice of the bulk ion mass is done only during the post-processing of the results of ORB5, if we need to know the values of equilibrium or perturbed quantities in non-normalized units. In particular, in the case of a hydrogen plasma, we get a value of the ion cyclotron frequency of . With this choice of bulk specie, we can calculate the plasma temperature and sound frequency. The three values of plasma temperature are 515 eV, 2060 eV, and 8240 eV. The sound frequency is defined as (with , which for reads ). We obtain the following three values of the sound velocity: m/s, m/s, m/s. These correspond to the following three values of the sound frequency: rad/s, rad/s, and rad/s.
The energetic particle distribution function is a double bump-on-tail, with two bumps at (see Fig. 1), like in Ref. . In this paper, is chosen. In order to initialize a distribution function which is function of constants of motion only, the modified variable is used instead of (similarly to Ref. [10, 22, 23]). Neumann and Dirichlet boundary conditions are imposed to the scalar potential, respectively at the inner and outer boundaries, and .
4 Linear dynamics
The linear dynamics of EGAMs in an equilibrium similar to the one considered in Sec. 3, has been investigated in Ref.  for a bulk plasma temperature given by , by means of ORB5 simulations and analytical theory. A scan on the EP concentration is performed here, similarly to what was done in Ref. , but with three different values of bulk plasma temperature, corresponding to three different values of , as described in Sec. 3. The dependence of the linear dynamics (frequency and growth rate) on the EP concentration is shown in Fig. 2. Both the frequency and the growth rates are observed to follow the qualitative scalings as described in Ref.  and . Note, in particular, that the growth rate does not grow linearly with .
The dependence of the frequency on the EP concentration is not observed to change with . Regarding the growth rate, no change is observed when going from to , meaning that the measured growth rate is basically given by the value of the drive, and the Landau damping here is negligible for the chosen values of EP concentration. On the other hand, when further increasing the temperature, and going to , a smaller value of growth rate is measured, meaning a higher Landau damping. The transit resonance velocity of the EP can be calculated by knowing the EGAM frequency of a specific simulation. Considering a case with as an example, the frequency is measured as: . For comparison, the GAM frequency for these parameters is . Then, the transit resonance velocity in the linear phase is calculated as , with being the EGAM linear frequency.
5 Scaling of the saturated amplitudes
In this Section, we focus on the value of the saturated electric field , and we investigate its dependence on the value of the linear growth rate and of the damping. The corresponding scaling of sheds light on the mechanism which is responsible for the saturation. A one-to-one comparison with the saturation mechanism of the BPI is also described.
The amplitude of the EGAM at saturation has been measured in different simulations performed with ORB5, for different values of the bulk temperature, and different values of the EP concentrations. As an example, the radial structure of a nonlinear simulation with , is depicted in Fig. 3-a. No sensible change in the radial wave-number is observed when going from the linear phase, to the saturation, and after the saturation. This confirms that in this particular configuration, where all equilibrium radial profiles are flat, EGAM can be treated as a 1-dimensional problem where the radial direction does not play an important role.
When varying only the bulk temperature, both the linear frequency and growth rate are observed to scale with the sound frequency, which is a good normalization frequency (consistently with Fig. 2). The saturation level increases with the linear growth rate, similarly to other instabilities like the BPI in a uniform system and the Alfvén instabilities in tokamaks. This is depicted in Fig. 3-b, where non-normalized units are used (in particular, the ion cyclotron frequency is selected as a time unit not depending on the temperature).
The scalings with the energetic particle concentration are also investigated. The results are shown in Fig. 4a. We obtain that, in the considered regime, the saturated level scales as the quadratic power of the linear growth rate. This quadratic scaling is typical for marginally stable bump-on-tail instabilities, as derived by O’Neil [24, 25].
We can consider the problem to be similar to a monochromatic beam-plasma system, in which particles are moving in the potential well of the perturbed electric field. Depending on the energy, some particles are trapped inside the well and execute bounce motion with frequency in the frame moving with the wave phase velocity. The resonant particles exchange energy with the mode, causing the amplitude to grow and the particles to redistribute in phase space, flattening the velocity distribution in the vicinity of the resonant parallel velocity . The drive is due to the positive slope of the particle distribution in the velocity space at the resonant parallel velocity, which acts as an inverse Landau damping. In the initial stage , the mode grows exponentially with a linear growth rate, making more and more particles to become trapped in the phase space. After some significant particle velocity redistribution, the power exchange between the wave and the particles is balanced, causing the wave amplitude to saturate.
Since the initial perturbation is negligible, the saturation level is determined by the exchange of energy between the mode and a band of resonant particles . The chirping of the mode seen in Fig. 5 is a strongly non-linear effect that occurs when the amplitude is large enough to have the trapped (and more generally resonant) particles with drastically change the dynamics of the mode through the modification of the distribution function and non-perturbative fast particle response. Namely, the mode dynamics is determined by all resonant particles that exhibit a continuous oscillation, trapping and detrapping in the potential well of the mode, thus contributing (non-perturbatively) to the non-adiabatic behavior observed in the simulation.
The quadratic scaling of the saturation level with the damping obtained with our simulations is similar to the saturation of the BPI, where it occurs due to wave-particle trapping . For the BPI, the original reference of M. B. Levin (1972) gives a value of at saturation  (for comparison, note that more recent numerical calculations find [34, 35]). For EGAMs, the bounce frequency is given by :
with being the mass of the energetic particle specie, considered equal to the bulk ion mass in this paper, the velocity matching the resonance condition, and the magnetic curvature drift. Therefore we have:
We emphasize that the value of depends on . This is a main difference with respect to the BPI, where there is only one value , with being the plasma frequency.
The dependence of the maximum electric field on the linear growth rate can be measured with the results of the numerical simulations. For the simulations shown in Fig. 4, we find:
The values of are found to depend on the bulk temperature. For the three chosen increasing values of , i.e. 0.0039, 0.0078 and 0.0156, we have respectively V/m, V/m, and V/m. Finally the relationship between the EP bounce frequency and the linear growth rate is obtained:
where is calculated as , which yields:
Note that here, depends on the EGAM frequency, which changes with the intensity of the drive, given here by the EP concentration (see Fig. 2). As a comparison, note that in the problem of the BPI, solved in Ref. [24, 25, 30], on the other hand, the mode frequency is assumed to be constant and equal to the frequency measured in the absence of EP. On the other hand, the values of are found not to depend on the damping, which changes with the three different bulk temperatures: an interpolation can be drawn for all considered simulations, and shown to depend on only (see Fig. 4b). In this sense, the formula given in Eq. 10 is universal for the chosen regime, because it does not depend on the bulk plasma temperature.
The considered equilibrium has been chosen in order to excite EGAMs out of a GAM, and not out of a Landau pole, as described in Ref. . This choice of the mode is done, in order to have a one-to-one correspondence with the BPI, where the mode which is excited by the energetic electron beam is the Langmuir wave which is an eigenmode of the system in the absence of EP. Following this consideration, we can consider the interpolation of the results shown in Fig. 4, and take the extrapolation to , which is the limit assumed in the resolution of the BPI. In this case, the extrapolation gives a unique value, which defines the EGAM instability, i.e. . This is to be compared with the value of obtained for the BPI [30, 34, 35], i.e. (originally estimated as by Levin).
In this section, we show the results of the measurement of the time evolution of the EGAM frequency. In Sec. 5, we have shown that a quadratic scaling of the saturated electric field on the linear growth rate is found. This quadratic scaling, has a one-to-one correspondence on the Langmuir wave problem investigated by O’Neil, where the saturation occurs due to wave-particle trapping. The wave-particle trapping mechanism, is usually referred to as adiabatic, meaning that a slowly increasing potential well traps more and more energetic particles. In this adiabatic regime, the mode frequency varies very slowly with respect to the bounce frequency. On the other hand, in the EGAM case considered here, we show that the saturation is not strictly adiabatic, but a transition between adiabatic and non-adiabatic regime occurs at the time of the saturation.
In this section, we take again the EGAM case with , , as a typical case, and we investigate the variation of the frequency in time in comparison with the bounce frequency. For the measurement of the frequency, we use the radial zonal electric field measured at s=0.5. The measurement of the frequency is performed in two different ways: a) as an average of the period between several EGAM oscillation peaks, as shown in Fig. 5-a; b) with a short-time Fourier transform (STFT), as shown in Fig. 5-b.
With the first technique, namely measuring the frequency by inversion of the period between neighbouring peaks, an upward chirping is observed in the nonlinear phase, of the order of 10% of the linear frequency. This means that the resonance condition changes in time, with resonance velocity slightly increasing at the time of the saturation or in the later phase (see Fig. 5-a). A dominantly upward chirping of EGAMs was previously observed and documented in Ref. [36, 11].
The second technique, consists in measuring the frequency with a short-time Fourier transform (STFT) on a Hamming time-window. With this technique, the error bar in frequency is large (due to the small number of oscillations in the nonlinear regime around the saturation), namely of the order of 10-20% (see Fig. 5-b). With such a big error-bar, no clear upward chirping is observed. Near the time of the saturation, i.e. , only one mode is observed. This is the condition of application of the direct technique of measurement of the frequency described above, where the frequency can be measured as the inverse of the period among peaks.
As mentioned above, the EP bounce frequency depends on the mode amplitude. Its value for the considered simulation is shown in Fig. 6-a. When the EGAM frequency evolves slowly in time with respect to the inverse of the bounce frequency, then the EP can bounce back and forth many times, and this is called adiabatic dynamics. The adiabaticity parameter, defined as , measures the level of adiabaticity of the dynamics. The time evolution of the adiabaticiy parameter for the considered simulation is depicted in Fig. 6-b (with being the instantaneous mode frequency). A transition from adiabatic to non-adiabatic dynamics near the saturation is observed. In particular, near the saturation, the EP do not have the time to perform many bounces during the nonlinear modification of the wave. From this respect, the EGAM dynamics and the BPI are not in analogy. The details of the EGAM saturation mechanism will be investigated with a power-balance diagnostics in phase space, and discussed in a separate paper.
7 Summary and discussion
The importance of understanding the nonlinear dynamics of EGAM, i.e. energetic particle (EP) driven geodesic acoustic modes (GAM), is due to their possible role in interacting with turbulence and with the EP population present in tokamak plasmas as a result of nuclear reactions and/or heating mechanisms. In particular, the level of the nonlinear saturation of EGAM is directly related to their efficiency in regulating turbulence or redistributing EP in phase space.
In this paper, we have investigated the problem of the nonlinear saturation of EGAM with the gyrokinetic particle-in-cell code ORB5, focusing on the wave-particle nonlinearity. Electrostatic collisionless simulations have been considered, with circular flux surfaces magnetic equilibria, and neglecting the kinetic effects of electrons.
The level of the saturated electric field has been shown to scale quadratically with the linear growth rate, similarly to the beam-plasma instability (BPI) in 1D uniform plasmas, and to Alfvén Eigenmodes (AE) near marginal stability. We note that, in the case of beta-induced AE, due to the finite radial mode structure compared to resonant particle radial excursion, deviation from the is observed due to the competition of resonance detuning and radial decoupling [37, 38]. A similar deviation is anticipated as one further increases the drive of EPs, from the correspondence of EP pitch angle scattering by frequency chirping EGAMs and the radial wave-particle pumping by frequency chirping finite-n BAEs. This will be the content of our next publication.
We have also investigated the relationship of the bounce frequency of the EP in the field of the EGAM, with the linear growth rate, finding a linear dependence, similarly to the BPI. The linear coefficient, in the case of the EGAM problem, is proportional to the square root of the EGAM frequency, and is not a constant like in the beam-plamsa instability. In the limit of , which corresponds to the BPI problem when the mode frequency does not move sensibly from the value of the plasma frequency , a unique value of the linear coefficient can be estimated: , to be compared with the factor of the BPI. This value of defines the EGAM in the regime of interest. The investigation of the dependence of on equilibrium parameters like the safety factor , the magnetic flux surface elongation , the EP distribution function, and the effect of kinetic electrons is left to a future work.
Finally, we have investigated the temporal variation of the EGAM frequency during the saturation phase. The initial adiabatic regime, defined as the regime where the time derivative of the EGAM frequency is much smaller than the squared EP bounce frequency, has been shown to have a transition to a non-adiabatic regime at the saturation. The exact analytical model of this phenomenon is still a work in progress and will be the first expansion of this work. The simplified picture is that the mode will chirp to maintain maximal exchange of power with the energetic particles , while balancing the Landau damping. The mode frequency is at the same time bound by the wave dispersion relation and plasma equilibrium, making their simultaneous examination with the power exchange equation necessary for determining the rate and direction of chirping. The non-adiabatic behavior is enabled by the non-perturbative response of the energetic particles, which consists of significant distortion of the distribution function in the resonant region, as well as of the particle trapping and detrapping . Differently, the adiabatic chirping is connected with either external driven fluctuations, or slow energetic particle redistribution, provided that the drive is sufficiently weak .
Interesting discussions with N. Carlevaro, G. Montani, F. Zonca on the nonlinear dynamics of the energetic particles are acknowledged. Interesting discussions with A. Mishchenko and B. Finney McMillan on the treatment of EGAMs with ORB5 are also acknowledged. Part of this work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053, within the framework of the Nonlinear energetic particle dynamics (NLED) European Enabling Research Project, WP 15-ER-01/ENEA-03. The views and opinions expressed herein do not necessarily reflect those of the European Commission. Part of this work has been funded by the Centre de Coopération Universitaire Franco-Bavarois - Bayerisch-Französisches Hochschulzentrum, for the grant FK 30_15. Simulations were performed on the Marconi supercomputer within the framework of the OrbZONE and ORBFAST projects. Part of this work was done while two of the authors (A. Biancalani and I. Novikau) were visiting LPP-Palaiseau, whose team is acknowledged for the hospitality. One of the authors (I.C.) wants to thank the Max Planck Princeton Center for the support during this work.
-  A. Hasegawa, C. G. Maclennan, and Y. Kodama, Phys. Fluids 22, 2122 (1979)
-  M.N. Rosenbluth and F.L. Hinton, Phys. Rev. Lett. 80,4 724 (1998)
-  P.H. Diamond, S.-I. Itoh, K. Itoh, and T. S. Hahm, Plasma Phys. Controlled Fusion 47, R35 (2005)
-  N. Winsor, J. L. Johnson, and J. M. Dawson, Phys. Fluids 11, 2448, (1968)
-  F. Zonca and L. Chen, Europhys. Lett. 83, 35001 (2008)
-  G. Fu et al. Phys. Rev. Lett. 101, 185002 (2008)
-  Z. Qiu, F. Zonca and L. Chen, Plasma Phys. Controlled Fusion 52, 095003 (2010)
-  Z. Qiu, F. Zonca, L. Chen Plasma Science and Tech. 13, 257 (2011)
-  Z. Qiu, F. Zonca, L. Chen Phys. Plasmas 19 082507 (2012)
-  D. Zarzoso et al. Phys. Plasmas 19, 022102-1 (2012)
-  H. Wang, Y. Todo and C. C. Kim Phys. Rev. Lett. 110, 155006 (2013)
-  G. B. Girardo, et al. Phys. Plasmas 21 092507 (2014)
-  D. Zarzoso et al. Phys. Rev. Lett. 110, 125002 (2013)
-  D. Zarzoso et al. Nucl. Fusion 57 072011 (2017)
-  P. Lauber, ITPA technical meeting on energetic particles (2014)
-  L. Horváth et al. Nucl. Fusion 56 112003 (2016)
-  Z. Qiu, I. Chavdarovski, A. Biancalani, and J. Cao, Phys. Plasmas 24, 072509 (2017)
-  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).
-  A. Bottino, T. Vernay, B. D. Scott, S. Brunner, R. Hatzky, S. Jolliet, B. F. McMillan, T. M. Tran, and L. Villard Plasma Phys. Controlled Fusion 53, 124027 (2011)
-  A. Bottino and E. Sonnendrücker, J. Plasma Phys. 81, 435810501 (2015).
-  A. Biancalani, A. Bottino, C. Ehrlacher, V. Grandgirard, G. Merlo, I. Novikau, Z. Qiu, E. Sonnendrücker, X. Garbet, T. Görler, S. Leerink, F. Palermo, and D. Zarzoso, Phys. Plasmas 24, 062512 (2017)
-  A. Biancalani, et al. Nucl. Fusion 54, 104004 (2014)
-  D. Zarzoso, A. Biancalani, A. Bottino, Ph. Lauber, E. Poli, J.-B. Girardo, X. Garbet and R.J. Dumont, Nucl. Fusion 54, 103006 (2014)
-  T. M. O’Neil Phys. Fluids 8, 2255 (1965)
-  T. M. O’Neil and J. H. Malmberg Phys. Fluids 11, 1754 (1968)
-  Z. Qiu, L. Chen and F. Zonca 41st EPS Conference on Plasma Physics, P4.004 (2014)
-  S. Briguglio, X. Wang, F. Zonca, G. Vlad, G. Fogaccia, C. Di Troia and V. Fusco, Phys. Plasmas 21, 112301 (2014).
-  G. J. Morales and T. M. O’Neil, Phys. Rev. Letters 28, 417 (1972).
-  T. Armon, L. Friedland, J. Plasma Phys. 82, 705820501 (2016)
-  M. B. Levin, M. G. Lyubarskii, I. N. Onishchenko, V. D. Shapiro, and V. I. Shevchenko, Sov. Phys. JETP 35, 898 (1972)
-  G. Fu and J. W. Van Dam, Phys. Fluids B 1, 1949 (1989)
-  H. Sugama, Phys. Plasmas 7, 466 (2000)
-  A. J. Brizard and T. S. Hahm, Rev. Mod. Phys 79, 421 (2007)
-  M. Lesur, Y. Idomura, and X. Garbet, Phys. Plasmas 16, 092305 (2009).
-  N. Carlevaro, et al, to be submitted (2017)
-  H.L. Berk, C.J. Boswell, D. Borba, A.C.A. Figueiredo, T. Johnson, M.F.F. Nave, S.D. Pinches, S.E. Sharapov and JET EFDA contributors, Nucl. Fusion 46 S888-S897 (2006)
-  X. Wang, S. Briguglio, L. Chen, C. Di Troia, G. Fogaccia, G. Vlad, and F. Zonca Phys. Rev. E 86, 045401(R) (2012)
-  H. S. Zhang, Z. Lin, and I. Holod Phys. Rev. Lett. 109, 025001 (2012)
-  L. Chen and F. Zonca, Rev. Mod. Phys. 88 1, (2016)
-  B. N. Breizman, H. L. Berk, M. S. Pekker, F. Porcelli, G. V. Stupakov, and K. L. Wong, Phys. Plasmas 4(5) 1559-1568, (1997)