Low-frequency turbulence suppression by amplification of synchronized zonal flow with energetic particle driven modes

We investigate the suppression of the turbulent transport associated to the emergence of spontaneous internal transport barriers, due to the combined collective effects played by trapped electrons and energetic passing ions. Numerical experiments performed with a ‘particle mode’ model based on a double gyro-average over the fast cyclotron phase and over the bounce (or transit) phases are used to show the role played by energetic particles in the suppression of the ion-temperature-gradient driven turbulence. We show this occur via phase locking in a Kuramoto-type synchronization process.


Introduction
The important role of zonal flows (ZFs) in regulating turbulence and transport in tokamaks is now broadly accepted. However the impact of fast particles in the ZF amplification and in the associated turbulence suppression is an open question. Instabilities driven by energetic particles have been observed in many tokamaks and helical devices. In the Large Helical Device (LHD), the resistive interchange mode can be * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. destabilized by the neutral beam injection. A new mode is this way generated and referred to as the 'energetic ion driven resistive interchange mode' [1][2][3][4]. This mode is characterized by a global structure on the poloidal and toroidal numbers, m and n, respectively, which here satisfy m/n = 1/1. One of the distinguished features of this mode is the down-chirping of its frequency during the intermittent burst which characterizes its temporal dynamics. Furthermore, the mode is determined by the precession frequency of the energetic particles. The energetic ion driven interchange mode can be therefore included in the family of the energetic particles modes [5], similarly to the oscillating fishbones observed in tokamaks. In the latter, different physical mechanisms leading to the suppression of turbulence involving the dynamics of ZFs have been reported. A first scenario, recently proposed in [6][7][8], is the possibility to generate a new type of internal transport barrier (ITB), induced by fast ions driven by ion-cyclotron resonance heating when the ion magnetic drift frequency is close to the frequency of underlying ion temperature gradient (ITG) instabilities. This seems to be a necessary condition for the suppression of turbulence. Modes driven by energetic particles can be thus excited nonlinearly, depleting this way the turbulent energy by enhancing the energy transfer to ZFs [9] (which remain the dominant saturation channel of turbulence). The high frequency component of ZFs, generally attributed to the presence of geodesic acoustic modes (GAMs), can be this way affected by the energetic particles. This leads to the emergence of energetic GAM modes, named EGAM [10], at a frequency typically comprised between 0.40ω GAM and 2ω GAM (see [11][12][13]). In this frequency interval the ITG frequency may resonate with the ion precession frequency [14]. However, the role of the low-frequency component of the zonal flow (LFZF) remains little studied. From global gyrokinetic simulations of toroidal ITG-driven turbulence, another physical mechanism has been predicted in [15], which takes place when the unstable region width is larger than the ITG correlation length, and phase synchronization leads to a modification of the turbulent transport (displaying a gyro-Bohm scaling) on a large spatial scale. Furthermore, the synchronization of GAMs with magnetic fluctuations has been identified in the edge plasmas of the HL2A tokamak [16]. The experiments presented in [16] were conducted in a deuterium plasma heated via electron-cyclotron resonance. This suggests that energetic electrons have an impact on the low-frequency zonal flow (LFZF) component. While a mechanism of synchronization of GAMs with magnetic fluctuations has been clearly demonstrated, the amplification of the LFZF mode by synchronization effects remains an open question, since it involves fluctuations of electrostatic nature at much lower frequencies.
Recently, in [17,18] we have extended the trapped particle model, initially developed to study the impact of ITG-type instabilities on trapped particle modes, so to include circulating and counter-circulating ions. All these classes of particles are described by a separate distribution function in an actionangle formulation of the corresponding Vlasov equation, and are coupled one each other via the Poisson equation. In this context we name here 'particle mode' each of these classes of particles, depending on the general form of their precession frequency, and regardless of the stable or unstable nature of the linear modes which can be associated to them (see section 2, below). The resulting interaction between all these particle modes can lead to an amplification of the LFZF component by synchronization with the interchange type modes, which can exhibit a hybrid fluid-kinetic character in the presence of low-frequency ITG instabilities, that is, trapped ion modes (TIM) and trapped electron modes (TEM). While in the high-frequency regime the synchronization process of GAMs can be driven by magnetic fluctuations, the situation can be somewhat different in the low-frequency regime, where it is the low frequency of the ZF that can be affected by trapped particles modes having frequencies close to the magnetic drift frequency.
In this letter, we show that a population of fast particles and a synchronization mechanism mediated by phase locking can affect the ITG-driven turbulence in the very low-frequency regime, in which the LFZF can be affected. Its amplitude can be thus enhanced up to a critical value, leading to the suppression of turbulence: synchronization of LFZFs with lowfrequency trapped particles modes thus provides a new fundamental process in the suppression of ITG turbulence.

Physical model
The description of particle modes is based on a Hamiltonian-Jacobi formalism, where each class of particles is represented by a distribution function f s = f s,Es,κ (ψ, α, t) in the ψ − α phase space plane, where ψ is the poloidal flux and α = φ − q 0 θ its corresponding phase.
In the Hamilton-Jacobi formalism, and for each particle 'species', namely ions, electrons and classes of trapped particles, one introduces a set of action-angle variables (J, α), where J = (J 1 , J 2 , J 3 ) and α = (α 1 , α 2 , α 3 ), the particle motion being integrable and quasi-periodic. The pair (J 1 , α 1 ) is related to the cyclotron motion of frequency ω cs where α 1 ∼ ω cs (J) t stands for the gyrophase and J 1 is proportional to the magnetic moment µ s . The pair (J 2 , α 2 ), where α 2 ∼ ω bs (J) t is the bounce phase and J 2 depends on the energy E s , is related to the poloidal motion. Finally, the pair (J 3 , α 3 ) corresponds to the toroidal motion, where J 3 ≡ P φ s = m s Rv φ s + e s ψ is the toroidal kinetic momentum, v φ s is the toroidal velocity and α 3 ∼ ω ds (J) t is the precession phase. Assuming an equilibrium magnetic field of the form B = I (ψ) ∇φ + ∇φ × ∇θ = ∇ (φ − q (ψ) θ) × ∇ψ, the trajectory of a particle can be written in terms of (see [19,20]) eψ G = J 3 + eψ (J, α 2 ), θ G = ϵ c α 2 +θ (J, α 2 ) and φ G = α 3 + q (J 3 )θ (J, α 2 ) +φ (J, α 2 ). The parameter ε c is taken to be one (zero, respectively) for passing (trapped, respectively) particles. The functionsθ (J, α 2 ) andφ (J, α 2 ) are periodic with respect to α 2 and describe the bounce (or transit) motion and the deviation from the regular precession motion, respectively. The low-frequency response for a given particle class s is obtained by gyroaveraging the Vlasov equation, since ω 2 , ω 3 ≪ ω 1 ; J 1 thus becomes an adiabatic invariant. Furthermore, for trapped (or circulating) particles, we have ω 3 ≪ ω 2 and also the bounce (or transit) motion can be gyroaveraged. So, also J 2 becomes an adiabatic invariant (also known ad the bounce-transit action) with respect to the bounce period, even if the guiding center deviates from the magnetic field line −see [21]. It is thanks to this reduction that, for the sake of notation, we could denote with the symbol α the only surviving variable α 3 in the model. Such hamiltonian formalism was used in [22,23] to recover the expression for the Rosenbluth-Hinton residual ZF [24] in axisymmetric tokamaks. Formulas for trapped-particle and passing particle guiding-center orbits were obtained in terms of the Jacobi elliptic functions in [25] or were generalized for finite inverse aspect ration in [26]. In general, such a scale separation is not always valid for the transit motion (e.g. in the case of GAMs) but it remains valid if we look at the very low frequency phenomena close to the ion precession frequency, which, can affect the LFZF. This model therefore excludes the study of GAMs. It can be used to describe the dynamics of low-frequency TIM modes induced by ITG-type instabilities, as well as TEM modes or interchange-type instabilities which may have a hybrid fluid-kinetic character. However, it does not allow us to study the dynamics of 'high-frequency' ITGrelated modes.
The model consists of four classes of populations of particles and allows the description of the associated modes. Each population s ∈ {i, e, ±} verifies N κ N E Vlasov equations, each cluster being obtained after a 'double gyro-average', i.e. by performing a time average over fast scales corresponding to the cyclotron plus bounce motion for trapped electron (s = e) and ion (s = i) populations, and to the cyclotron plus transit motion for circulating ions (s = +) and countercirculating ions (s = −). The model of the particle modes is based on the solution of N gyrokinetic-type Vlasov equations, nonlinearly coupled by the quasi-neutrality condition: Here e s is the species charge, where the apex c = ±1 stands for circulating and counter-circulating ions and is the fraction of trapped particles (taken identical for both electrons and ions). The polarization terms ('(Pol. terms)' in equation (1)) take into account the Laplacian of each species where C e = e e /T e = e/T e and C i = e i /T i . In equation (1) (1) can be quantified as usual using Pade approximates, which makes it depend on the population energy E s , on the 'banana' orbit width and on the Larmor radius of the population, ρ s , which, once normalized to the reference radius of the magnetic surface r 0 , we write as ρ * ≡ ρ s /r 0 (cf e.g. [17,27,28]).
Each distribution function f s (ψ, α, t; κ, E s ) also depends on the two adiabatic invariants κ and E s which play the role of two parameters. Since a precession frequency of the type ω ds = ω ds (k) E s is associated to each of these classes of particles, it is possible to characterize each of them as an eigenmode of the Hamiltonian. It is in this sense that here (and in [17,18]) we call each class of particles so defined a 'particle mode'.
As an illustrative example, let us consider the simplest situation, say of a class of TIMs, with a population of adiabatic electrons. Without loss of generality, it is possible to find a marginal solution of the Vlasov-Poisson system (1) and (2), for which the Hamiltonian writes H = ω d (κ) E i ψ + J 0i ϕ. Looking for stationary solutions of the form f = 2 nω d (κ) defines the class of (collisionless, interchange-type [20,29]) TIM eigenmodes, at the threshold of the ion-temperaturegradient (ITG) instability. Like for ITG modes, these TIMs are subject to an interchange instability driven by the magnetic curvature. The 'marginal solution' is defined by F κ, where the quantity E i is normalized to the ion temperature T i . Performing a linear Landau-type analysis of the Vlasov-Poisson system by considering a perturbation on the distribution function and the electric potential in the form δf n (ψ) exp i (nα − ωt) and δϕ n (ψ) exp i (nα − ωt), around the marginal equilibrium F κ,Es , one obtains , and with the usual Landau prescription on the imaginary part. Equation (3) allows us to define the threshold of the ITG-instability, together with the dispersion relation of collisionless (kinetic) TIM mode in the form ω ≡ ω TIM = 3 2 nω d (κ) as previously obtained from the research of the marginal solution. Such an linear analysis can be extended to other classes of eigenmodes, namely TEM, co-(or counter) circulating modes. However, these trapped-ion modes differ from those discovered by Kadomtsev and Pogutse [30], the dissipative (or collisional) trapped ion modes (DTIM). It is however found that for sufficiently large ion temperature gradients (i.e. above the critical threshold), the usual picture of dissipative modes is modified, and a transition from DTIMs to collisionless TIMs takes place, when the curvature drift effects dominate.

Hamiltonian phase synchronization aspects
By following the same kind of analysis recently developed in references [27,28,31] for the Vlasov-Poisson system, we can identify Λ s (ψ, α, t) ≡ − [J 0s ϕ, f s ], as a complex order parameter. Equation (1), which thus reads can be written using the Fourier representation of f s (ψ, α, t) = n f s,n (ψ, t) e inα . Thus, by introducing f s,n = | f s,n | e iφs,n and Λ s,n = |Λ s,n | e iΘs,n , and by separating the real and imaginary parts, equation (4) becomes Equations (5) and (6) describe the system in the form of a Kuramoto model, in which each 'particle mode' (i.e. the quantities labeled with s) plays the role of an oscillator of phase φ n = φ s,n , of frequency ω n = −nω ds (E s , κ), whose dynamics is governed by a Kuramoto-type equation: In the standard Kuramoto model [32][33][34][35], which describes the weak interaction of N < ∞ oscillators, the coupling term K is a constant and ω n is a random frequency chosen from a distribution with a probability density function g (ω). In the mean field Kuramoto approach, which corresponds to the N → +∞ coupling limit, the collective behavior of an ensemble of oscillators is quantified by the order parameter r = (1/N) N n=1 e iφn (following the notation of [21]) and undergoes a phase transition from incoherence (described by the r = 0 state) to synchronization (r ∼ 1) as the coupling coefficient K trespasses a critical value K c . Most of the research on the Kuramoto model focused on the case where the distribution g (ω) is unimodal (see [5]). Specifically, g is usually assumed to be symmetric about a maximum value at a given frequency ω 0 and to decrease monotonically to zero outside of ω ∼ ω 0 . In that case, it was found in [34,35] that the transition in the continuum limit (where N → +∞) takes the form of a bifurcation towards a discrete mode emerging out from a continuum spectrum [36]. However, a natural question arises at this stage: what is the impact of a more general frequency distribution on the bifurcation mechanism, and, more precisely, what does it happen when we take into account several classes of particle modes, which are characterized by very different frequency distributions ω d,s (E s , κ). To explore the impact of several particles modes on the incoherence-synchronization transition, we first solve numerically the discrete Kuramoto equation (7) for N < ∞.

Numerical results for the Kuramoto model
Let us first discuss some numerical results obtained by integrating equations (5) and (6) for a finite number N of oscillators represented by particle modes. In a first simulation, whose results are shown in figure 1, we consider the case of two clusters of trapped particles oscillating as TIM and TEM. Each cluster is made of N κ N Es = 22 2 modes, which corresponds to an ensemble of 968 'oscillators'. The electron and ion temperatures are both equal to a reference temperature T 0 . For trapped particle populations s ∈ {e, i}, the precession fre- Here B(κ) = is the magnetic shear, assumed to be a constant in this model; K(κ) and E(κ) are the complete elliptic integral of the first and second kind, respectively. A characteristic drift frequency ω d0 = q 0 T 0 /(eB 0 r 0 R 0 ); T 0 can be defined, T 0 , B 0 and R 0 being  (7) is first solved numerically, for two clusters of TEMs and TIMs. In the upper frame, the order parameter r = ∑ n e iφn is shown versus the coupling parameter K. As expected, a transition from incoherence (r ∼ 0) to synchronization (r ∼ 1) takes places when we increase K. In the lower frame: the case of four clusters of particle modes: TIM/TEM (with κ < 1 and with temperature of T i = T 0 and Te = 2.5T 0 ) but now including also the two clusters of circulating (σ ∥ = 1, κ > 1 and T i = 0.5T 0 ) and counter-circulating (σ ∥ = −1, κ > 1, and T i = 0.5T 0 ). The introduction of two new clusters of passing ions of smaller ion temperature, together with a TEM of higher temperature can affect the bifurcation (lower frame). a reference temperature, magnetic field amplitude and radius values, respectively.
At the top of figure 1. The value of the order parameter r is shown as a function of K in the upper frame: the expected fast transition to a full synchronization, similar to a bifurcation, occurs for K c ∼ 2 π g(0) = 2×968 π×44 ≃ 14. In the lower frame, in figure 1, we display a similar diagnostics for a simulation initialized with four clusters of particle modes, TIM and TEM (having κ < 1), and two clusters of circulating (σ ∥ = 1; κ > 1) and counter-circulating (σ ∥ = −1 and κ > 1) ions. This time N κ N E = 15 2 and N + N − = 2 2 , for a total of N = 900 modes. The physical parameters are T i = T 0 , T e = 2.5T 0 , T ± = 0.5T 0 , q 0 = 2 and ρ * = 0.10. To show the role of energetic particles, we consider the following repartition in frequency for passing particles with σ ∥ = ±1: We thus find that the synchronization degree is highly dependent on the chosen distribution in frequency, but that for large K the transition is robust and the dynamics leads to a bifurcation. Moreover, the simultaneous inclusion of both passing ions and hot electrons breaks the symmetry of the distribution g (ω) in frequency (not shown here), in contrast with the previous case without energetic particles. The transition begins at a lower critical value K c , here found to be close to 2×900 π×g(ω0) ∼ 8.81 for a value of g (ω 0 ) ≃ 65. However the larger domain in the Kspace (see bottom panel of figure 1), facilitates the transition to a global phase synchronization.

Gyrokinetic numerical experiments
We now discuss a set of numerical results obtained with a selfconsistent integration of equations (1) and (2) by means of an extended version of the TERESA code [17,18]. This uses a semi-Lagrangian scheme [, 28, 37, 38] to integrate equation (1) for each class of solutions defined by the set of adiabatic invariants (κ, E s ). We use normalized quantities: time is normalized to ω −1 d0 ; ψ is given in ∆ψ units. The electric potential ϕ, together with the constants C e and C i , are expressed in ω d0 ∆ψ units. The initial distribution function is Here the function h 1 = ∆τ ω d Es Ts − 3 2 ψ, with ω d = ω d (κ) ω d0 , allows us to start an ITG instability from the initial temperature gradient ∆τ , provided the latter is chosen larger than the critical threshold ∆τ th = C e /ξ where ξ = 1 − 3 4 δ 2 bs κ + 15 64 δ 4 bs κ . The function Es Ts − 5 2 ϕ max sin(2πψ), instead, allows us to excite a shear flow by introducing a potential perturbation of amplitude ϕ max . If the species s does not include an energetic population, the function F 0s = e −Es/Ts is usually a Maxwellian distribution. In presence of a fast ion beam of average energy E b and density n b , the standard Maxwellian is replaced by In [39,40], it was shown that it is possible to decompose the total particle energy of each species s into three different contributions involved in the energy transfer in ITG -type instabilities: the kinetic energy ϵ kin,s = ⟨ϖ kin ⟩ α,ψ , the energy ϵ turb = ⟨ϖ turb ⟩ α,ψ contributing to turbulence and the energy ϵ ZF = ⟨ϖ ZF ⟩ α,ψ associated to ZFs. The energy densities are ϖ kin = T s ψ ⟨ω d (κ)E s ⟩ κ,E , ϖ turb = Aδϕ 2 + s B s ∇ s ϕ 2 and ϖ ZF = s D s (∂ ⟨ϕ⟩ α /∂ψ), respectively.
Here ⟨.⟩ α,ψ = 1 4π´2 π 0 dα´1 0 dψ introduces an average over the phase space variables and the average ⟨.⟩ κ,E is performed over the pitch angle κ and over the energy, depending on the type of particle mode: for trapped particles where C κ is a constant that grants the integral to be equal to one and  (1) and (2)). On top, an initial population of energetic co-circulating ions is taken into account, the other species being described by a Maxwellian. At the bottom, the initial energetic ion beam is in the population of counter-circulating ions. The same physical parameters were used in both simulations.
K is the elliptic function of the first kind. The other constants are defined as A = C e (1 − R p ), B s = R p C s for trapped particle modes and B σ ∥ = 1−Rp 2 C i for circulating or countercirculating ion modes; D s = s B s δ 2 s and δ s ∈ δ be , δ bi , δ cσ ∥ where the subscript b denotes the banana width and c the transit particle width. To define ϖ turb we have introduced the normalized gradient in the ψ − α phase space in the form ∇ s = ρ * ∂ α e α + δs ∆ψ ∂ ψ e ψ . It should be noted that the presence of an important ZF energy component does not necessarily imply turbulence suppression, which can occur instead only in presence of a sheared ZF (see numerical results below).
Numerical simulations have been performed using a phase space sampling of N α N ψ = 1025 × 257 grid points, a time step of ∆tω d0 = 0.005 for more than N κ N E = 16 × 120 adiabatic invariants in the pitch angle (κ) and energy (E s ) directions. All numerical simulations were performed for the same set of physical parameters, namely for a magnetic shear of s m = 0.80, a maximum shear flow amplitude ϕ max / (ω d0 ∆ψ) of 2.50, an ion temperature of T i = T 0 , an electron temperature of T e = 10T 0 , C e = 0.40 and a polarization term C i = 0.60 in normalized units. Note that in the LHD experiments [1][2][3][4], the electron temperature varies from T i to ∼10T i . Figure 2 shows the temporal evolution of the energies E ZF and E turb , which correspond to the contributions of the zonal flow (including a sheared component) and of the ITG-turbulence respectively, in the presence of a population of energetic passing particles. In figure 2, on top, the energetic ion population introduced in the simulation corresponds to a passing ion distribution with E b = 10T 0 and with a density of n b /n 0 of 5%. At times smaller than 5ω −1 d0 we clearly observe a strong amplification of the turbulence, of more than a 3/0.0006 ≃ 5000 factor with respect to the initial phase of the dynamics. This corresponds to the achievement of a self-organized turbulence state. At later times, however, a complete suppression of turbulence, visible at times tω d0 ⩾ 5 and tω d0 ⩾ 10 in the upper and lower panel of figures 2, respectively, takes place under the combined action of energetic ions and of the stabilization by a ('rotational-like') shear flow. This results in a strong increase of the zonal flux, which is estimated to be of the order of 1/0.0005 ≃ 2000 (compare figures 2 and 4, below).
In the absence of a shear flow (shown in figure 3), we observe again the strong growth of the LFZF energy component (at a slightly higher level than in the case of figure 2) and we observe a significant reduction of turbulence, as well, although the latter this time does not disappear completely.
Finally, in figure 4 we show the results of a simulation performed using the same physical parameters of figure 2, but now in the absence of energetic ions: no turbulence suppression is this time observed.
Several points should be stressed: • The effect discussed above is strongly related to the presence of the hot electrons, although they, alone, are not sufficient to completely reduce the level of ionic turbulence shown in figure 4. Indeed, at a temperature of T e ∼ T 0 (not shown here), the decrease in turbulence is not as much efficient, and a complete suppression of the latter is not observed. The top panel in figure 2 corresponds to a synchronization process rather than to a resonant process, since the copassing particle mode and the TEM propagate in opposite directions. The energy redistribution occurs between fast particles and TEM turbulence when the magnetic -drift frequency of fast ions is close to the frequency of the underlying TEM instability, i.e. when both propagate in the same direction and ω d fast ions ∼ q 0 T ehot /(eB 0 r 0 R 0 ). • The turbulence suppression also occurs if energetic ions are introduced into the counter-circulating ion population (case shown at the bottom of figure 2). In this case, a slightly larger amplitude of the ZF energy is measured. Furthermore, the dynamics is somewhat faster than that shown in the top panel of figure 2 but, as a whole, it remains similar. This regime corresponds to that observed in [6] when a resonance occurs between (counter-passing) energetic ions and electron-driven TEM turbulence.  figure 2, the turbulence becomes predominant on long scales and is associated with interchange-ITG type turbulence (see figure 5). Lower frame: the beginning of the instability is zoomed, by showing that the ZF is initially excited when turbulence is growing. However, the turbulence level is not high enough to maintain a large amplitude ZF.
• The ITB associated to the ZF formation is very stable. Once the turbulence suppression is obtained, the plasma dynamics is quasi-stationary, at least over scales of the order of tω d0 ∼ 100. No mechanisms of destabilization of the internal transport barrier have been observed in this time interval. • The bifurcation to a synchronous state, which leads to the suppression of turbulence, is related to the emergence of the global mode n = 1 (see below). This process has been observed for a wide range of physical parameters, for example at lower electron temperatures (T e ∼ T i ) but at a higher beam density n b = 0.1n 0 in [18]. Figure 5 show the contour plot of the distribution function of trapped electrons in the (ψ, α) phase space at three different times. The case when fast co-passing ions are taken into account is shown in the left column (the corresponding plots of E turb and E ZF versus time are shown on top panel of Figure 2), whereas, for comparison, the corresponding case with no energetic ions is shown in the right panels (cf figure 4 for the corresponding E turb and E ZF ). On top left panel, a global toroidal mode n = 1 emerges as a result of a global transition from partial to global synchronization. Such an LFZF structure corresponds to the global E × B flow pattern, sometimes referred to as a ZF staircase, which is associated to the formation of the internal transport barrier. Here it seems to occur when the global structure with n = 1 is formed, similarly to what was previously observed in [29,41], in the case of a fishbone in presence of fast passing ions. An initial mean shear flow (here provided by the function h 2 in equation (10)) is needed to build-up the staircase solution, although, alone, it does not suffice for the emergence of the ITB. Nevertheless, it facilitates the phase locking of ZF and ITG-turbulence at a frequency close to ω d fast ions ∼ q 0 T ehot /(eB 0 r 0 R 0 ). Thus the LFZF Figure 5. Plot of the electron distribution function in the (ψ, α) phase space: including fast energetic ions (left panels) and without energetic ions (right panels). While in the right panels the electron distribution clearly shows the emergence of a streamer structure driven by an interchange-ITG instability, turbulence suppression is observed in presence of fast co-passing ions, in the left panels. mode characterized by the toroidal number n = 1 is determined by the precession frequency of the energetic particles. The fact that this is the same behavior observed for the energetic interchange mode in the LHD device suggests that this is a general mechanism. The appearance of fine phase-space filaments in the distribution f TEM during the saturation phase (cf the middle left panel of figure 5) is reminiscent of the (Landau) phase-mixing mechanism associated to the suppression of turbulence. By contrast, in the absence of a fast ion beam (right column of figure 5), the interchange-ITG driven turbulence, plays a dominant role due to the relative weakness of ZFs.

Conclusion
Thanks to the self-consistent integration of gyro-averaged Vlasov Equations written in action-angle variables, which we have re-formulated and interpreted in terms of a Kuramototype model, we have shown that a global transition is produced by a phase synchronization mechanism between TEM modes and circulating ions. This synchronization is induced by a population of fast ions, which makes the different particle modes to adopt similar frequencies close to the precession frequency of the fast ions, even if the modes are initially distributed over a wide spectrum of eigenfrequencies. An asymmetry is thus introduced by the fast ion population on the distribution of the natural frequencies of the mode populations present in the system. The LFZF is strongly affected by the strong growth of the initial interchange instability. This leads to the formation of a global structure, a ZF staircase of mesoscopic size, characterized by the formation of a global structure on the fundamental toroidal number (n = 1), resulting in the formation of an internal transport barrier. This result is quite general and, although it depends on the electron-to-ion temperature ratio and on the ratio between the fast ion beam and the trapped particle density, is confirmed for different values of these parameters: numerical experiments performed with the gyroaveraged Vlasov model including circulating ions have shown that the lower T e is, making T e /T i close to unity, the larger n b must be, making n b /n 0 close to a decimal fraction of unity.