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

A. Ghizzo, D. Del Sarto

► To cite this version:

A. Ghizzo, D. Del Sarto. Low-frequency turbulence suppression by amplification of synchronized zonal flow with energetic particle driven modes. Nuclear Fusion, 2023, 63 (10), pp.104002. 10.1088/1741-4326/acee95. hal-04545536

HAL Id: hal-04545536

https://hal.science/hal-04545536

Submitted on 14 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Low-frequency turbulence suppression by amplification of synchronized zonal flow with energetic particle driven modes

A. Ghizzo, D. Del Sarto

► To cite this version:

A. Ghizzo, D. Del Sarto. Low-frequency turbulence suppression by amplification of synchronized zonal flow with energetic particle driven modes. Nuclear Fusion, 2023, 63 (10), pp.104002. 10.1088/1741-4326/acee95. hal-04545536

HAL Id: hal-04545536

https://hal.science/hal-04545536v1

Submitted on 14 Apr 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
LETTER • OPEN ACCESS

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

To cite this article: A. Ghizzo and D. Del Sarto 2023 *Nucl. Fusion* **63** 104002

View the [article online](#) for updates and enhancements.

You may also like

- A new formulation of gradient boosting
  Alex Wozniakowski, Jayne Thompson,
  Mile Gu et al.

- Systematic design of a multi-input multi-output controller by model-based decoupling: a demonstration on TCV using multi-species gas injection
  J.T.W. Koenders, A. Perek, C. Galperti et al.

- The planted k-factor problem
  Gabriele Sicuro and Lenka Zdeborová
LETTER • OPEN ACCESS

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

To cite this article: A. Ghizzo and D. Del Sarto 2023 *Nucl. Fusion* **63** 104002

View the [article online](#) for updates and enhancements.

You may also like

- A new formulation of gradient boosting
  Alex Wozniakowski, Jayne Thompson,
  Mile Gu et al.

- Systematic design of a multi-input multi-output controller by model-based decoupling: a demonstration on TCV using multi-species gas injection
  J.T.W. Koenders, A. Perek, C. Galperti et al.

- The planted k-factor problem
  Gabriele Sicuro and Lenka Zdeborová
Letter

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

A. Ghizzo* and D. Del Sarto

Institut Jean Lamour, UMR 7198, Université de Lorraine, BP 70239, Vandoeuvre les Nancy cedex,
54506, France

E-mail: alain.ghizzo@univ-lorraine.fr

Received 28 March 2023, revised 17 June 2023

Accepted for publication 9 August 2023

Published 23 August 2023

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

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.

Keywords: plasma physics, vlasov equation, gyrokinetic modelling, transport barrier, phase synchronization

(Some figures may appear in colour only in the online journal)

## 1. 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

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–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–8],

<figure><img src="image_3.png" /></figure>

Original Content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

<footnote>* Author to whom any correspondence should be addressed.</footnote>
Letter

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

A. Ghizzo* and D. Del Sarto

Institut Jean Lamour, UMR 7198, Université de Lorraine, BP 70239, Vandoeuvre les Nancy cedex,
54506, France

E-mail: alain.ghizzo@univ-lorraine.fr

Received 28 March 2023, revised 17 June 2023

Accepted for publication 9 August 2023

Published 23 August 2023

<figure><img src="image_2.png" /></figure>

CrossMark

## Abstract

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.

Keywords: plasma physics, vlasov equation, gyrokinetic modelling, transport barrier, phase synchronization

(Some figures may appear in colour only in the online journal)

## 1. 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

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–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–8],

<figure><img src="image_3.png" /></figure>

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.

<footnote>* Author to whom any correspondence should be addressed.</footnote>
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.4\omega_{GAM}$ and $2\omega_{GAM}$ (see [11–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 action-angle 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 low-frequency trapped particles modes thus provides a new fundamental process in the suppression of ITG turbulence.

## 2. 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,E_s,\kappa}(\psi, \alpha, t)$ in the $\psi-\alpha$ phase space plane, where $\psi$ is the poloidal flux and $\alpha = \varphi - q_0\theta$ 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, \alpha)$, where $J = (J_1, J_2, J_3)$ and $\alpha = (\alpha_1, \alpha_2, \alpha_3)$, the particle motion being integrable and quasi-periodic. The pair $(J_1, \alpha_1)$ is related to the cyclotron motion of frequency $\omega_{cs}$ where $\alpha_1 \sim \omega_{cs}(J)t$ stands for the gyrophase and $J_1$ is proportional to the magnetic moment $\mu_s$. The pair $(J_2, \alpha_2)$, where $\alpha_2 \sim \omega_{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, \alpha_3)$ corresponds to the toroidal motion, where $J_3 \equiv P_{\varphi s} = m_s R v_{\varphi s} + e_s \psi$ is the toroidal kinetic momentum, $v_{\varphi s}$ is the toroidal velocity and $\alpha_3 \sim \omega_{ds}(J)t$ is the precession phase. Assuming an equilibrium magnetic field of the form $B = I(\psi) \nabla\varphi + \nabla\varphi \times \nabla\theta = \nabla(\varphi - q(\psi)\theta) \times \nabla\psi$, the trajectory of a particle can be written in terms of (see [19, 20]) $e\psi_G = J_3 + e\hat{\psi}(J, \alpha_2)$, $\theta_G = \epsilon_c\alpha_2 + \hat{\theta}(J, \alpha_2)$ and $\varphi_G = \alpha_3 + q(J_3)\hat{\theta}(J, \alpha_2) + \hat{\phi}(J, \alpha_2)$. The parameter $\epsilon_c$ is taken to be one (zero, respectively) for passing (trapped, respectively) particles. The functions $\hat{\theta}(J, \alpha_2)$ and $\hat{\phi}(J, \alpha_2)$ are periodic with respect to $\alpha_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 $\omega_2, \omega_3 \ll \omega_1$; $J_1$ thus becomes an adiabatic invariant. Furthermore, for trapped (or circulating) particles, we have $\omega_3 \ll \omega_2$ and also the bounce (or transit) motion can be gyroaveraged. So, also $J_2$ becomes an adiabatic invariant (also known as 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 $\alpha$ the only surviving variable $\alpha_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 ratio in [26]. In general, such a scale separation is not always valid for the transit motion (e.g. in the
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.4\omega_{GAM}$ and $2\omega_{GAM}$ (see [11–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 action-angle 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 low-frequency trapped particles modes thus provides a new fundamental process in the suppression of ITG turbulence.

## 2. 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,E_s,\kappa}(\psi, \alpha, t)$ in the $\psi-\alpha$ phase space plane, where $\psi$ is the poloidal flux and $\alpha = \varphi - q_0\theta$ 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, \alpha)$, where $J = (J_1, J_2, J_3)$ and $\alpha = (\alpha_1, \alpha_2, \alpha_3)$, the particle motion being integrable and quasi-periodic. The pair $(J_1, \alpha_1)$ is related to the cyclotron motion of frequency $\omega_{cs}$ where $\alpha_1 \sim \omega_{cs}(J)t$ stands for the gyrophase and $J_1$ is proportional to the magnetic moment $\mu_s$. The pair $(J_2, \alpha_2)$, where $\alpha_2 \sim \omega_{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, \alpha_3)$ corresponds to the toroidal motion, where $J_3 \equiv P_{\varphi s} = m_s R v_{\varphi s} + e_s \psi$ is the toroidal kinetic momentum, $v_{\varphi s}$ is the toroidal velocity and $\alpha_3 \sim \omega_{ds}(J)t$ is the precession phase. Assuming an equilibrium magnetic field of the form $B = I(\psi) \nabla\varphi + \nabla\varphi \times \nabla\theta = \nabla(\varphi - q(\psi)\theta) \times \nabla\psi$, the trajectory of a particle can be written in terms of (see [19, 20]) $e\psi_G = J_3 + e\hat{\psi}(J, \alpha_2)$, $\theta_G = \epsilon_c\alpha_2 + \hat{\theta}(J, \alpha_2)$ and $\varphi_G = \alpha_3 + q(J_3)\hat{\theta}(J, \alpha_2) + \hat{\phi}(J, \alpha_2)$. The parameter $\epsilon_c$ is taken to be one (zero, respectively) for passing (trapped, respectively) particles. The functions $\hat{\theta}(J, \alpha_2)$ and $\hat{\phi}(J, \alpha_2)$ are periodic with respect to $\alpha_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 $\omega_2, \omega_3 \ll \omega_1$; $J_1$ thus becomes an adiabatic invariant. Furthermore, for trapped (or circulating) particles, we have $\omega_3 \ll \omega_2$ and also the bounce (or transit) motion can be gyroaveraged. So, also $J_2$ becomes an adiabatic invariant (also known as 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 $\alpha$ the only surviving variable $\alpha_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 ratio 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’ ITG-related modes.

The model consists of four classes of populations of particles and allows the description of the associated modes. Each population $s \in \{i, e, \pm\}$ verifies $N_\kappa 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 counter-circulating ions ($s=-$). The model of the particle modes is based on the solution of N gyrokinetic-type Vlasov equations,

$$ \frac{\partial f_s}{\partial t} + \omega_{ds}(E_s, \kappa) \frac{\partial f_s}{\partial \alpha} + [J_{0s}\phi, f_s] = 0, \quad (1) $$

nonlinearly coupled by the quasi-neutrality condition:

$$ C_e(1 - \mathcal{R}_p)(\phi - \langle\phi\rangle_\alpha) + (\text{Pol. terms}) = \sum_s \text{sgn}(e_s)n_s. \quad (2) $$

Here $e_s$ is the species charge, $\sum_s \text{sgn}(e_s)n_s = \mathcal{R}_p(n_i^{(t)} - n_e^{(t)}) + \frac{1-\mathcal{R}_p}{2}(\sum_{c=\pm 1} n_i^{(c)}) - (1-\mathcal{R}_p)$, where the apex $c = \pm 1$ stands for circulating and counter-circulating ions and $\mathcal{R}_p = \frac{2\sqrt{2e}}{\pi}$ 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 in the form $-\mathcal{R}_p C_i \bar{\Delta}_i^{(t)} \phi - \mathcal{R}_p C_e \bar{\Delta}_e^{(t)} \phi - (1-\mathcal{R}_p) C_i \bar{\Delta}_i^{(c)} \phi$, where $C_e = e_e/T_e = e/T_e$ and $C_i = e_i/T_i$. In equation (1), $[\dots, \dots]$ denotes the standard Poisson brackets defined by $[A, B] = \frac{\partial A}{\partial \psi} \frac{\partial B}{\partial \alpha} - \frac{\partial A}{\partial \alpha} \frac{\partial B}{\partial \psi}$. The gyrokinetic operator $J_{0s}$ in the Poisson bracket of equation (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, $\rho_s$, which, once normalized to the reference radius of the magnetic surface $r_0$, we write as $\rho^* \equiv \rho_s/r_0$ (cf e.g. [17, 27, 28]).

Each distribution function $f_s(\psi, \alpha, t; \kappa, E_s)$ also depends on the two adiabatic invariants $\kappa$ and $E_s$ which play the role of two parameters. Since a precession frequency of the type $\omega_{ds} = \omega_{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 = \omega_d(\kappa)E_i\psi + J_{0i}\phi$. Looking for stationary solutions of the form $f=$

$F_{\kappa, E_i}(\psi, \alpha - \omega t)$ and $\phi = \phi(\psi, \alpha - \omega t)$, equation (1) reads $[(\omega_d(\kappa)E_i - \omega)\psi + J_{0i}\phi, F_{\kappa, E_i}] = 0$ and $\omega \equiv \omega_{TIM} = \frac{3}{2}n\omega_d(\kappa)$ defines the class of (collisionless, interchange-type [20, 29]) TIM eigenmodes, at the threshold of the ion-temperature-gradient (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_{\kappa, E_s} = T_i^{-3/2} \exp(-E_i) (1 + \Delta\tau\omega_d(\kappa)(E_i - \frac{3}{2}))$ 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 $\delta f_n(\psi)\exp i(n\alpha - \omega t)$ and $\delta\phi_n(\psi)\exp i(n\alpha - \omega t)$, around the marginal equilibrium $F_{\kappa, E_s}$, one obtains

$$ C_e \delta \phi_n = \left\langle n \Delta \tau \exp(-E_i) \frac{E_i - \frac{3}{2}}{E_i - \frac{\omega}{n\omega_d(\kappa)}} J_{0i}^2 \delta \phi_n \right\rangle_{\kappa, E_i} \quad (3) $$

where $\langle\cdot\rangle_{\kappa, E_i} = \frac{2}{\sqrt{\pi}} \int_0^{+\infty} dE \sqrt{E} \int_0^1 d\kappa \kappa K(\kappa)$, 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 $\omega \equiv \omega_{TIM} = \frac{3}{2}n\omega_d(\kappa)$ 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.

### 3. 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 $\Lambda_s(\psi, \alpha, t) \equiv -[J_{0s}\phi, f_s]$, as a complex order parameter. Equation (1), which thus reads

$$ \frac{\partial f_s}{\partial t} + \omega_{ds}(E_s, \kappa) \frac{\partial f_s}{\partial \alpha} = \Lambda_s(\psi, \alpha, t), \quad (4) $$

can be written using the Fourier representation of $f_s(\psi, \alpha, t) = \sum_n f_{s,n}(\psi, t) e^{in\alpha}$. Thus, by introducing $f_{s,n} = |f_{s,n}| e^{i\varphi_{s,n}}$ and $\Lambda_{s,n} = |\Lambda_{s,n}| e^{i\Theta_{s,n}}$, and by separating the real and imaginary parts, equation (4) becomes

$$ \frac{\partial |f_{s,n}|}{\partial t} = |\Lambda_{s,n}| \cos(\Theta_{s,n} - \varphi_{s,n}) \quad (5) $$

$$ \frac{\partial \varphi_{s,n}}{\partial t} = -n\omega_{ds}(E_s, \kappa) + \frac{|\Lambda_{s,n}|}{|f_{s,n}|} \sin(\Theta_{s,n} - \varphi_{s,n}). \quad (6) $$

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
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’ ITG-related modes.

The model consists of four classes of populations of particles and allows the description of the associated modes. Each population $s \in \{i, e, \pm\}$ verifies $N_\kappa 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 counter-circulating ions ($s=-$). The model of the particle modes is based on the solution of N gyrokinetic-type Vlasov equations,

$$ \frac{\partial f_s}{\partial t} + \omega_{ds}(E_s, \kappa) \frac{\partial f_s}{\partial \alpha} + [J_{0s}\phi, f_s] = 0, \quad (1) $$

nonlinearly coupled by the quasi-neutrality condition:

$$ C_e(1 - \mathcal{R}_p)(\phi - \langle\phi\rangle_\alpha) + (\text{Pol. terms}) = \sum_s \text{sgn}(e_s)n_s. \quad (2) $$

Here $e_s$ is the species charge, $\sum_s \text{sgn}(e_s)n_s = \mathcal{R}_p(n_i^{(t)} - n_e^{(t)}) + \frac{1-\mathcal{R}_p}{2}(\sum_{c=\pm 1} n_i^{(c)}) - (1-\mathcal{R}_p)$, where the apex $c = \pm 1$ stands for circulating and counter-circulating ions and $\mathcal{R}_p = \frac{2\sqrt{2e}}{\pi}$ 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 in the form $-\mathcal{R}_p C_i \bar{\Delta}_i^{(t)} \phi - \mathcal{R}_p C_e \bar{\Delta}_e^{(t)} \phi - (1-\mathcal{R}_p) C_i \bar{\Delta}_i^{(c)} \phi$, where $C_e = e_e/T_e = e/T_e$ and $C_i = e_i/T_i$. In equation (1), $[\dots, \dots]$ denotes the standard Poisson brackets defined by $[A, B] = \frac{\partial A}{\partial \psi} \frac{\partial B}{\partial \alpha} - \frac{\partial A}{\partial \alpha} \frac{\partial B}{\partial \psi}$. The gyrokinetic operator $J_{0s}$ in the Poisson bracket of equation (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, $\rho_s$, which, once normalized to the reference radius of the magnetic surface $r_0$, we write as $\rho^* \equiv \rho_s/r_0$ (cf e.g. [17, 27, 28]).

Each distribution function $f_s(\psi, \alpha, t; \kappa, E_s)$ also depends on the two adiabatic invariants $\kappa$ and $E_s$ which play the role of two parameters. Since a precession frequency of the type $\omega_{ds} = \omega_{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 = \omega_d(\kappa)E_i\psi + J_{0i}\phi$. Looking for stationary solutions of the form $f=$

$F_{\kappa, E_i}(\psi, \alpha - \omega t)$ and $\phi = \phi(\psi, \alpha - \omega t)$, equation (1) reads $[(\omega_d(\kappa)E_i - \omega)\psi + J_{0i}\phi, F_{\kappa, E_i}] = 0$ and $\omega \equiv \omega_{TIM} = \frac{3}{2}n\omega_d(\kappa)$ defines the class of (collisionless, interchange-type [20, 29]) TIM eigenmodes, at the threshold of the ion-temperature-gradient (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_{\kappa, E_s} = T_i^{-3/2} \exp(-E_i) (1 + \Delta\tau\omega_d(\kappa)(E_i - \frac{3}{2}))$ 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 $\delta f_n(\psi)\exp i(n\alpha - \omega t)$ and $\delta\phi_n(\psi)\exp i(n\alpha - \omega t)$, around the marginal equilibrium $F_{\kappa, E_s}$, one obtains

$$ C_e \delta \phi_n = \left\langle n \Delta \tau \exp(-E_i) \frac{E_i - \frac{3}{2}}{E_i - \frac{\omega}{n\omega_d(\kappa)}} J_{0i}^2 \delta \phi_n \right\rangle_{\kappa, E_i} \quad (3) $$

where $\langle\cdot\rangle_{\kappa, E_i} = \frac{2}{\sqrt{\pi}} \int_0^{+\infty} dE \sqrt{E} \int_0^1 d\kappa \kappa K(\kappa)$, 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 $\omega \equiv \omega_{TIM} = \frac{3}{2}n\omega_d(\kappa)$ 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.

### 3. 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 $\Lambda_s(\psi, \alpha, t) \equiv -[J_{0s}\phi, f_s]$, as a complex order parameter. Equation (1), which thus reads

$$ \frac{\partial f_s}{\partial t} + \omega_{ds}(E_s, \kappa) \frac{\partial f_s}{\partial \alpha} = \Lambda_s(\psi, \alpha, t), \quad (4) $$

can be written using the Fourier representation of $f_s(\psi, \alpha, t) = \sum_n f_{s,n}(\psi, t) e^{in\alpha}$. Thus, by introducing $f_{s,n} = |f_{s,n}| e^{i\varphi_{s,n}}$ and $\Lambda_{s,n} = |\Lambda_{s,n}| e^{i\Theta_{s,n}}$, and by separating the real and imaginary parts, equation (4) becomes

$$ \frac{\partial |f_{s,n}|}{\partial t} = |\Lambda_{s,n}| \cos(\Theta_{s,n} - \varphi_{s,n}) \quad (5) $$

$$ \frac{\partial \varphi_{s,n}}{\partial t} = -n\omega_{ds}(E_s, \kappa) + \frac{|\Lambda_{s,n}|}{|f_{s,n}|} \sin(\Theta_{s,n} - \varphi_{s,n}). \quad (6) $$

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
$\varphi_n = \varphi_{s,n}$, of frequency $\omega_n = -n\omega_{ds}(E_s, \kappa)$, whose dynamics is governed by a Kuramoto-type equation:

$$ \dot{\varphi}_n = \omega_n + \sum_{j=1}^{N} \frac{K}{N} \sin(\varphi_j - \varphi_n). \quad (7) $$

In the standard Kuramoto model [32–35], which describes the weak interaction of $N < \infty$ oscillators, the coupling term $K$ is a constant and $\omega_n$ is a random frequency chosen from a distribution with a probability density function $g(\omega)$. In the mean field Kuramoto approach, which corresponds to the $N \to +\infty$ coupling limit, the collective behavior of an ensemble of oscillators is quantified by the order parameter $r = (1/N) \sum_{n=1}^{N} e^{i\varphi_n}$ (following the notation of [21]) and undergoes a phase transition from incoherence (described by the $r=0$ state) to synchronization ($r \sim 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(\omega)$ is unimodal (see [5]). Specifically, $g$ is usually assumed to be symmetric about a maximum value at a given frequency $\omega_0$ and to decrease monotonically to zero outside of $\omega \sim \omega_0$. In that case, it was found in [34, 35] that the transition in the continuum limit (where $N \to +\infty$) 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 $\omega_{d,s}(E_s, \kappa)$. To explore the impact of several particles modes on the incoherence- synchronization transition, we first solve numerically the discrete Kuramoto equation (7) for $N < \infty$.

## 4. 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_\kappa N_E = 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 \in \{e, i\}$, the precession frequency is $\omega_{ds}(E_s, \kappa) = \pm T_s \bar{\omega}_d(\kappa) E_s$, where

$$ \bar{\omega}_d(\kappa) = B(\kappa) \left[ \frac{2\mathcal{E}(\kappa)}{\mathcal{K}(\kappa)} - 1 + 4s_m \left( \frac{\mathcal{E}(\kappa)}{\mathcal{K}(\kappa)} + \kappa^2 - 1 \right) \right]. \quad (8) $$

Here $B(\kappa) = \frac{1}{1-\epsilon+2\epsilon\kappa^2}$, $\kappa$ being the pitch angle satisfying $0 \le \kappa < 1$ for trapped particles; $s_m = \frac{r_0}{q_0} \left(\frac{dq}{dr}\right)_{r_0}$ is the magnetic shear, assumed to be a constant in this model; $\mathcal{K}(\kappa)$ and $\mathcal{E}(\kappa)$ are the complete elliptic integral of the first and second kind, respectively. A characteristic drift frequency $\omega_{d0} = q_0 T_0 / (eB_0 r_0 R_0)$; $T_0$ can be defined, $T_0$, $B_0$ and $R_0$ being

<figure><img src="image_4.png" /><figcaption>**Figure 1.** The discrete Kuramoto equation (7) is first solved numerically, for two clusters of TEMs and TIMs. In the upper frame, the order parameter $r = \sum_n e^{i\varphi_n}$ is shown versus the coupling parameter $K$. As expected, a transition from incoherence ($r \sim 0$) to synchronization ($r \sim 1$) takes places when we increase $K$. In the lower frame: the case of four clusters of particle modes: TIM/TEM (with $\kappa < 1$ and with temperature of $T_i = T_0$ and $T_e = 2.5T_0$) but now including also the two clusters of circulating ($\sigma_\parallel = 1$, $\kappa > 1$ and $T_i = 0.5T_0$) and counter-circulating ($\sigma_\parallel = -1$, $\kappa > 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).</figcaption></figure>

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 \sim \frac{2}{\pi g(0)} = \frac{2 \times 968}{\pi \times 44} \approx 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 $\kappa < 1$), and two clusters of circulating ($\sigma_\parallel = 1$; $\kappa > 1$) and counter-circulating ($\sigma_\parallel = -1$ and $\kappa > 1$) ions. This time $N_\kappa 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_\pm = 0.5T_0$, $q_0 = 2$ and $\rho^* = 0.10$. To show the role of energetic particles, we consider the following repartition in frequency for passing particles with $\sigma_\parallel = \pm 1$:

$$ \omega_{di}(E, \kappa) = \frac{T_{\text{circ}}}{T_0} \left( \bar{\omega}_d(\kappa) E_i + \sigma_\parallel \frac{\bar{\omega}_i(\kappa)}{q_0^2 \rho^*} \sqrt{E_i} \right). \quad (9) $$
$\varphi_n = \varphi_{s,n}$, of frequency $\omega_n = -n\omega_{ds}(E_s, \kappa)$, whose dynamics is governed by a Kuramoto-type equation:

$$ \dot{\varphi}_n = \omega_n + \sum_{j=1}^{N} \frac{K}{N} \sin(\varphi_j - \varphi_n). \quad (7) $$

In the standard Kuramoto model [32–35], which describes the weak interaction of $N < \infty$ oscillators, the coupling term $K$ is a constant and $\omega_n$ is a random frequency chosen from a distribution with a probability density function $g(\omega)$. In the mean field Kuramoto approach, which corresponds to the $N \to +\infty$ coupling limit, the collective behavior of an ensemble of oscillators is quantified by the order parameter $r = (1/N) \sum_{n=1}^{N} e^{i\varphi_n}$ (following the notation of [21]) and undergoes a phase transition from incoherence (described by the $r=0$ state) to synchronization ($r \sim 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(\omega)$ is unimodal (see [5]). Specifically, $g$ is usually assumed to be symmetric about a maximum value at a given frequency $\omega_0$ and to decrease monotonically to zero outside of $\omega \sim \omega_0$. In that case, it was found in [34, 35] that the transition in the continuum limit (where $N \to +\infty$) 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 $\omega_{d,s}(E_s, \kappa)$. To explore the impact of several particles modes on the incoherence- synchronization transition, we first solve numerically the discrete Kuramoto equation (7) for $N < \infty$.

## 4. 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_\kappa N_E = 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 \in \{e, i\}$, the precession frequency is $\omega_{ds}(E_s, \kappa) = \pm T_s \bar{\omega}_d(\kappa) E_s$, where

$$ \bar{\omega}_d(\kappa) = B(\kappa) \left[ \frac{2\mathcal{E}(\kappa)}{\mathcal{K}(\kappa)} - 1 + 4s_m \left( \frac{\mathcal{E}(\kappa)}{\mathcal{K}(\kappa)} + \kappa^2 - 1 \right) \right]. \quad (8) $$

Here $B(\kappa) = \frac{1}{1 - e^{2\epsilon\kappa^2}}$, $\kappa$ being the pitch angle satisfying $0 \le \kappa < 1$ for trapped particles; $s_m = \frac{r_0}{q_0} \left(\frac{dq}{dr}\right)_{r_0}$ is the magnetic shear, assumed to be a constant in this model; $\mathcal{K}(\kappa)$ and $\mathcal{E}(\kappa)$ are the complete elliptic integral of the first and second kind, respectively. A characteristic drift frequency $\omega_{d0} = q_0 T_0 / (eB_0 r_0 R_0)$; $T_0$ can be defined, $T_0$, $B_0$ and $R_0$ being

<figure><img src="image_4.png" /><figcaption>**Figure 1.** The discrete Kuramoto equation (7) is first solved numerically, for two clusters of TEMs and TIMs. In the upper frame, the order parameter $r = \sum_n e^{i\varphi_n}$ is shown versus the coupling parameter $K$. As expected, a transition from incoherence ($r \sim 0$) to synchronization ($r \sim 1$) takes places when we increase $K$. In the lower frame: the case of four clusters of particle modes: TIM/TEM (with $\kappa < 1$ and with temperature of $T_i = T_0$ and $T_e = 2.5T_0$) but now including also the two clusters of circulating ($\sigma_\parallel = 1, \kappa > 1$ and $T_i = 0.5T_0$) and counter-circulating ($\sigma_\parallel = -1, \kappa > 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).</figcaption></figure>

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 \sim \frac{2}{\pi g(0)} = \frac{2 \times 968}{\pi \times 44} \approx 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 $\kappa < 1$), and two clusters of circulating ($\sigma_\parallel = 1; \kappa > 1$) and counter-circulating ($\sigma_\parallel = -1$ and $\kappa > 1$) ions. This time $N_\kappa 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_\pm = 0.5T_0$, $q_0 = 2$ and $\rho^* = 0.10$. To show the role of energetic particles, we consider the following repartition in frequency for passing particles with $\sigma_\parallel = \pm 1$:

$$ \omega_{di}(E, \kappa) = \frac{T_{\text{circ}}}{T_0} \left( \bar{\omega}_d(\kappa) E_i + \sigma_{\parallel} \frac{\bar{\omega}_i(\kappa)}{q_0^2 \rho^*} \sqrt{E_i} \right). \quad (9) $$
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 $\frac{2 \times 900}{\pi \times g(\omega_0)} \sim 8.81$ for a value of g(ω₀) ≈ 65. However the larger domain in the K-space (see bottom panel of figure 1), facilitates the transition to a global phase synchronization.

## 5. Gyrokinetic numerical experiments

We now discuss a set of numerical results obtained with a self-consistent 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 ω_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

$$f_{\kappa, E_s}(\psi, \alpha, t=0) = (1+h_1(\psi)+h_2(\psi))F_{0s}(E_s). \quad (10)$$

Here the function $h_1 = \Delta\tau\omega_d(\frac{E_s}{T_s} - \frac{3}{2})\psi$, with $\omega_d = \bar{\omega}_d(\kappa)\omega_{d0}$, allows us to start an ITG instability from the initial temperature gradient $\Delta\tau$, provided the latter is chosen larger than the critical threshold $\Delta\tau_{th} = C_e/\xi$ where $\xi = 1 - \frac{3}{4}\langle\delta_{bs}^2\rangle_\kappa + \frac{15}{64}\langle\delta_{bs}^4\rangle_\kappa$. The function $h_2 = 2\pi C_i\delta_b^2|\frac{E_s}{T_s} - \frac{5}{2}|\phi_{\max}\sin(2\pi\psi)$, instead, allows us to excite a shear flow by introducing a potential perturbation of amplitude $\phi_{\max}$. If the species s does not include an energetic population, the function $F_{0s} = e^{-E_s/T_s}$ 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

$$F_{0s} = (1 - n_b)e^{-\frac{E_s}{T_s}} + n_b e^{-\left(\sqrt{\frac{E_s}{T_s}} - \sqrt{\frac{E_b}{T_s}}\right)^2}. \quad (11)$$

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 $\epsilon_{kin,s} = \langle\bar{\omega}_{kin}\rangle_{\alpha,\psi}$, the energy $\epsilon_{turb} = \langle\bar{\omega}_{turb}\rangle_{\alpha,\psi}$ contributing to turbulence and the energy $\epsilon_{ZF} = \langle\bar{\omega}_{ZF}\rangle_{\alpha,\psi}$ associated to ZFs. The energy densities are $\bar{\omega}_{kin} = T_s\psi\langle\bar{\omega}_d(\kappa)E_s\rangle_{\kappa,E}$, $\bar{\omega}_{turb} = A\delta\phi^2 + \sum_s B_s|\bar{\nabla}_s\phi|^2$ and $\bar{\omega}_{ZF} = \sum_s D_s(\partial\langle\phi\rangle_\alpha/\partial\psi)$, respectively. Here $\langle\cdot\rangle_{\alpha,\psi} = \frac{1}{4\pi}\int_0^{2\pi} d\alpha \int_0^1 d\psi$ introduces an average over the phase space variables and the average $\langle\cdot\rangle_{\kappa,E}$ is performed over the pitch angle κ and over the energy, depending on the type of particle mode: for trapped particles $\langle\cdot\rangle_{\kappa,E} = \frac{2}{\sqrt{\pi}}\int_0^{+\infty} dE\sqrt{E}\int_0^1 d\kappa\kappa K(\kappa)$; for circulating particles $\langle\cdot\rangle_{\kappa,E} = \frac{2}{\sqrt{\pi}}\int_0^{+\infty} dE\sqrt{E}\int_1^{+\infty} d\kappa C_\kappa K(\kappa^{-1})$, where $C_\kappa$ is a constant that grants the integral to be equal to one and

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Results of numerical simulations of the reduced Hamiltonian model (equations (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.</figcaption></figure>

$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_{\sigma_\parallel} = \frac{1 - R_p}{2}C_i$ for circulating or counter-circulating ion modes; $D_s = \sum_s B_s\delta_s^2$ and $\delta_s \in \{\delta_{be}, \delta_{bi}, \delta_{c\sigma_\parallel}\}$ where the subscript b denotes the banana width and c the transit particle width. To define $\bar{\omega}_{turb}$ we have introduced the normalized gradient in the $\psi-\alpha$ phase space in the form $\bar{\nabla}_s = \rho^*\partial_\alpha e_\alpha + \frac{\delta_s}{\Delta\psi}\partial_\psi e_\psi$. 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_\alpha N_\psi = 1025 \times 257$ grid points, a time step of $\Delta t\omega_{d0} = 0.005$ for more than $N_\kappa N_E = 16 \times 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 $\phi_{max}/(\omega_{d0}\Delta\psi)$ 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–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
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 $\frac{2 \times 900}{\pi \times g(\omega_0)} \sim 8.81$ for a value of g(ω₀) ≈ 65. However the larger domain in the K-space (see bottom panel of figure 1), facilitates the transition to a global phase synchronization.

## 5. Gyrokinetic numerical experiments

We now discuss a set of numerical results obtained with a self-consistent 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 ω_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

$$f_{\kappa, E_s}(\psi, \alpha, t=0) = (1+h_1(\psi)+h_2(\psi))F_{0s}(E_s). \quad (10)$$

Here the function $h_1 = \Delta\tau\omega_d(\frac{E_s}{T_s} - \frac{3}{2})\psi$, with $\omega_d = \bar{\omega}_d(\kappa)\omega_{d0}$, allows us to start an ITG instability from the initial temperature gradient $\Delta\tau$, provided the latter is chosen larger than the critical threshold $\Delta\tau_{th} = C_e/\xi$ where $\xi = 1 - \frac{3}{4}\langle\delta_{bs}^2\rangle_\kappa + \frac{15}{64}\langle\delta_{bs}^4\rangle_\kappa$. The function $h_2 = 2\pi C_i\delta_b^2|\frac{E_s}{T_s} - \frac{5}{2}|\phi_{max}\sin(2\pi\psi)$, instead, allows us to excite a shear flow by introducing a potential perturbation of amplitude $\phi_{max}$. If the species s does not include an energetic population, the function $F_{0s} = e^{-E_s/T_s}$ 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

$$F_{0s} = (1 - n_b)e^{-\frac{E_s}{T_s}} + n_b e^{-\left(\sqrt{\frac{E_s}{T_s}} - \sqrt{\frac{E_b}{T_s}}\right)^2}. \quad (11)$$

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 $\epsilon_{kin,s} = \langle\bar{\omega}_{kin}\rangle_{\alpha,\psi}$, the energy $\epsilon_{turb} = \langle\bar{\omega}_{turb}\rangle_{\alpha,\psi}$ contributing to turbulence and the energy $\epsilon_{ZF} = \langle\bar{\omega}_{ZF}\rangle_{\alpha,\psi}$ associated to ZFs. The energy densities are $\bar{\omega}_{kin} = T_s\psi\langle\bar{\omega}_d(\kappa)E_s\rangle_{\kappa,E}$, $\bar{\omega}_{turb} = A\delta\phi^2 + \sum_s B_s|\bar{\nabla}_s\phi|^2$ and $\bar{\omega}_{ZF} = \sum_s D_s(\partial\langle\phi\rangle_\alpha/\partial\psi)$, respectively. Here $\langle\cdot\rangle_{\alpha,\psi} = \frac{1}{4\pi}\int_0^{2\pi} d\alpha \int_0^1 d\psi$ introduces an average over the phase space variables and the average $\langle\cdot\rangle_{\kappa,E}$ is performed over the pitch angle κ and over the energy, depending on the type of particle mode: for trapped particles $\langle\cdot\rangle_{\kappa,E} = \frac{2}{\sqrt{\pi}}\int_0^{+\infty} dE\sqrt{E}\int_0^1 d\kappa\kappa K(\kappa)$; for circulating particles $\langle\cdot\rangle_{\kappa,E} = \frac{2}{\sqrt{\pi}}\int_0^{+\infty} dE\sqrt{E}\int_1^{+\infty} d\kappa C_\kappa K(\kappa^{-1})$, where $C_\kappa$ is a constant that grants the integral to be equal to one and

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Results of numerical simulations of the reduced Hamiltonian model (equations (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.</figcaption></figure>

$\mathcal{K}$ is the elliptic function of the first kind. The other constants are defined as $A = C_e(1 - \mathcal{R}_p)$, $B_s = \mathcal{R}_p C_s$ for trapped particle modes and $B_{\sigma_\parallel} = \frac{1 - \mathcal{R}_p}{2}C_i$ for circulating or counter-circulating ion modes; $D_s = \sum_s B_s\delta_s^2$ and $\delta_s \in \{\delta_{be}, \delta_{bi}, \delta_{c\sigma_\parallel}\}$ where the subscript b denotes the banana width and c the transit particle width. To define $\bar{\omega}_{turb}$ we have introduced the normalized gradient in the $\psi-\alpha$ phase space in the form $\bar{\nabla}_s = \rho^*\partial_\alpha e_\alpha + \frac{\delta_s}{\Delta\psi}\partial_\psi e_\psi$. 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_\alpha N_\psi = 1025 \times 257$ grid points, a time step of $\Delta t\omega_{d0} = 0.005$ for more than $N_\kappa N_E = 16 \times 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 $\phi_{max}/(\omega_{d0}\Delta\psi)$ 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–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
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Top frame: Results similar to those shown in figure 2, but in absence of a shear flow (i.e. with $\phi_{max} = 0$). Bottom frame: a zoom of the turbulence energy versus time.</figcaption></figure>

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\omega_{d0}^{-1}$ we clearly observe a strong amplification of the turbulence, of more than a 3/0.0006 $\simeq$ 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 $\omega_{d0} \gg 5$ and $\omega_{d0} \gg 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 $\simeq$ 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 \sim 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 co-passing 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 $\omega_{d\text{fastions}} \sim q_0 T_{\text{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><img src="image_6.png" /><figcaption>**Figure 3.** Top frame: Results similar to those shown in figure 2, but in absence of a shear flow (i.e. with $\phi_{max} = 0$). Bottom frame: a zoom of the turbulence energy versus time.</figcaption></figure>

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\omega_{d0}^{-1}$ we clearly observe a strong amplification of the turbulence, of more than a 3/0.0006 $\simeq$ 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 $\omega_{d0} \gg 5$ and $\omega_{d0} \gg 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 $\simeq$ 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 \sim 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 co-passing 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 $\omega_{d\text{fastions}} \sim q_0 T_{\text{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><img src="image_7.png" /><figcaption>**Figure 4.** Upper frame: results of a numerical simulation without energetic ions. The physical parameters are identical to those used in the results of figure 2. Although the energy levels are much lower, compared to those in 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.</figcaption></figure>

* 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\omega_{d0} \sim 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 \sim 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 $(\psi, \alpha)$ 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_{\text{turb}}$ and $E_{\text{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_{\text{turb}}$ and $E_{\text{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 \times 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 $\omega_{\text{dfastions}} \sim q_0 T_{\text{ehot}} / (eB_0 r_0 R_0)$. Thus the LFZF
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Upper frame: results of a numerical simulation without energetic ions. The physical parameters are identical to those used in the results of figure 2. Although the energy levels are much lower, compared to those in 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.</figcaption></figure>

* 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\omega_{d0} \sim 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 \sim 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 $(\psi, \alpha)$ 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_{\text{turb}}$ and $E_{\text{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_{\text{turb}}$ and $E_{\text{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 \times 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 $\omega_{\text{dfastions}} \sim q_0 T_{\text{ehot}} / (eB_0 r_0 R_0)$. Thus the LFZF
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Plot of the electron distribution function in the $(\psi, \alpha)$ 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.</figcaption></figure>

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.

## 6. 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 Kuramoto-type 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
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Plot of the electron distribution function in the $(\psi, \alpha)$ 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.</figcaption></figure>

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.

## 6. 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 Kuramoto-type 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<sub>e</sub> is, making T<sub>e</sub>/T<sub>i</sub> close to unity, the larger n<sub>b</sub> must be, making n<sub>b</sub>/n<sub>0</sub> close to a decimal fraction of unity.

## Acknowledgment

The authors are grateful to the IDRIS computational centre, Orsay, France, for computer time allocation on their computers. This work was granted access to the HPC resources (Grant A014054028) by GENCI (Grand Equipement National de Calcul Intensif).

## ORCID iDs

A. Ghizzo https://orcid.org/0000-0003-3870-9229

D. Del Sarto https://orcid.org/0000-0001-7489-8843

## References

[1] Du X.D. et al (LHD Experiment Group) 2015 *Phys. Rev. Lett.* **114** 155003

[2] Du X.D. et al (LHD Experiment Group) 2016 *Nucl. Fusion* **56** 016002

[3] Ogawa K. et al (LHD Experiment Group) 2020 *Nucl. Fusion* **60** 112011

[4] Idouakass M., Todo Y., Wang H., Wang J., Seki R. and Sato M. 2021 *Phys. Plasmas* **28** 080701

[5] Cheng L. and Zonca F. 2016 *Rev. Mod. Phys.* **88** 015008

[6] Di Siena A. et al (The ASDEX Upgrade Team) 2021 *Phys. Rev. Lett.* **127** 025002

[7] Di Siena A., Banon Navarro A. and Jenko F. 2020 *Phys. Rev. Lett.* **125** 105002

[8] Di Siena A., Gorler T., Poli E., Banon Navarro A., Biancalani A. and Jenko F. 2019 *Nucl. Fusion* **59** 124001

[9] Diamond P.H., Itoh S.I., Itoh K. and Hahm T.S. 2005 *Plasma Phys. Control. Fusion* **47** R35

[10] Fu G.Y. 2008 *Phys. Rev. Lett.* **101** 185002

[11] Nazikian R. et al 2008 *Phys. Rev. Lett.* **101** 185001

[12] Ido T. et al (LHD Experiment Group) 2016 *Phys. Rev. Lett.* **116** 015002

[13] Wang H., Todo Y., Ido T. and Suzuki Y. 2018 *Phys. Rev. Lett.* **120** 175001

[14] Zarzoso D. et al 2013 *Phys. Rev. Lett.* **110** 125002

[15] Yi D.S., Jhang H., Kim S.S. and Kwon J.M. 2022 *Nucl. Fusion* **62** 076037

[16] Zhao K.J. et al 2016 *Phys. Rev. Lett.* **117** 145002

[17] Ghizzo A. and Del Sarto D. 2022 *Phys. Plasmas* **29** 042506

[18] Ghizzo A. and Del Sarto D. 2022 *Phys. Plasmas* **29** 042507

[19] Garbet X., Laurent L., Mourgues F., Roubin J.P. and Samain A. 1990 *J. Comput. Phys.* **87** 249

[20] Depret G., Garbet X., Bertrand P. and Ghizzo A. 2000 *Plasma Phys. Control. Fusion* **42** 949

[21] Northrop T.G. 1963 *The Adiabatic Motion of Charged Particles* (Interscience Publishers)

[22] Duthoit F.X., Brizard A. and Hahm T.S. 2014 *Phys. Plasmas* **21** 122510

[23] Wang L. and Hahm T.S. 2009 *Phys. Plasmas* **16** 062309

[24] Rosenbluth M.N. and Hinton F.L. 1998 *Phys. Rev. Lett.* **80** 724

[25] Brizard A. 2011 *Phys. Plasmas* **18** 022508

[26] Stephens C.D., Garbet X. and Jenko F. 2020 *Phys. Plasmas* **27** 052504

[27] Guo Z.B., Wang Y.H. and Xu S.K. 2020 *Phys. Rev. E* **101** 03201(R)

[28] Ghizzo A., Del Sarto D., Garbet X., Sarazin Y. 2010 *Phys. Plasmas* **17** 092051

[29] Ghizzo A. and Del Sarto D. 2019 *Plasma* **2** 229

[30] Kadomtsev B.B. and Pogutse O.P. 1971 *Nucl. Fusion* **11** 67

[31] Xu S., Guo Z.B. and Gurcan A.D. 2021 *Phys. Rev. E* **103** 023208

[32] Kuramoto Y. 1984 *Chemical Oscillations, Waves and Turbulence* (Springer)

[33] Acebron J.A., Bonilla L.L., Perez Vicente C.J., Ritort F. and Spligler R. 2005 *Rev. Mod. Phys.* **77** 137

[34] Strogatz S.H., Mirollo R.E. and Matthews P.C. 1992 *Phys. Rev. Lett.* **68** 2730

[35] Martens E.A., Barreto E., Strogatz S.H., Ott E., So P. and Antonsen T.M. 2009 *Phys. Rev. E* **79** 026204

[36] Balmforth N.J. and Sassi R. 2000 *Physica D* **143** 21

[37] Sonnendrucker E., Roche J.R., Bertrand P. and Ghizzo A. 1999 *J. Comput. Phys.* **149** 201

[38] Ghizzo A., Bertrand P., Shoucri M., Fijalkow E. and Feix M.R. 1993 *J. Comput. Phys.* **108** 105

[39] Ghizzo A. and Palermo F. 2015 *Phys. Plasmas* **22** 082303

[40] Ghizzo A. and Palermo F. 2015 *Phys. Plasmas* **22** 082304

[41] Brochard G., Dumont R., Garbet X., Lutjens H., Nicolas T. and Orain F. 2018 *J. Phys. : Conf. Ser.* **1125** 012003
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<sub>e</sub> is, making T<sub>e</sub>/T<sub>i</sub> close to unity, the larger n<sub>b</sub> must be, making n<sub>b</sub>/n<sub>0</sub> close to a decimal fraction of unity.

## Acknowledgment

The authors are grateful to the IDRIS computational centre, Orsay, France, for computer time allocation on their computers. This work was granted access to the HPC ressources (Grant A014054028) by GENCI (Grand Equipement National de Calcul Intensif).

## ORCID iDs

A. Ghizzo https://orcid.org/0000-0003-3870-9229

D. Del Sarto https://orcid.org/0000-0001-7489-8843

## References

[1] Du X.D. et al (LHD Experiment Group) 2015 *Phys. Rev. Lett.* **114** 155003

[2] Du X.D. et al (LHD Experiment Group) 2016 *Nucl. Fusion* **56** 016002

[3] Ogawa K. et al (LHD Experiment Group) 2020 *Nucl. Fusion* **60** 112011

[4] Idouakass M., Todo Y., Wang H., Wang J., Seki R. and Sato M. 2021 *Phys. Plasmas* **28** 080701

[5] Cheng L. and Zonca F. 2016 *Rev. Mod. Phys.* **88** 015008

[6] Di Siena A. et al (The ASDEX Upgrade Team) 2021 *Phys. Rev. Lett.* **127** 025002

[7] Di Siena A., Banon Navarro A. and Jenko F. 2020 *Phys. Rev. Lett.* **125** 105002

[8] Di Siena A., Gorler T., Poli E., Banon Navarro A., Biancalani A. and Jenko F. 2019 *Nucl. Fusion* **59** 124001

[9] Diamond P.H., Itoh S.I., Itoh K. and Hahm T.S. 2005 *Plasma Phys. Control. Fusion* **47** R35

[10] Fu G.Y. 2008 *Phys. Rev. Lett.* **101** 185002

[11] Nazikian R. et al 2008 *Phys. Rev. Lett.* **101** 185001

[12] Ido T. et al (LHD Experiment Group) 2016 *Phys. Rev. Lett.* **116** 015002

[13] Wang H., Todo Y., Ido T. and Suzuki Y. 2018 *Phys. Rev. Lett.* **120** 175001

[14] Zarzoso D. et al 2013 *Phys. Rev. Lett.* **110** 125002

[15] Yi D.S., Jhang H., Kim S.S. and Kwon J.M. 2022 *Nucl. Fusion* **62** 076037

[16] Zhao K.J. et al 2016 *Phys. Rev. Lett.* **117** 145002

[17] Ghizzo A. and Del Sarto D. 2022 *Phys. Plasmas* **29** 042506

[18] Ghizzo A. and Del Sarto D. 2022 *Phys. Plasmas* **29** 042507

[19] Garbet X., Laurent L., Mourgues F., Roubin J.P. and Samain A. 1990 *J. Comput. Phys.* **87** 249

[20] Depret G., Garbet X., Bertrand P. and Ghizzo A. 2000 *Plasma Phys. Control. Fusion* **42** 949

[21] Northrop T.G. 1963 *The Adiabatic Motion of Charged Particles* (Interscience Publishers)

[22] Duthoit F.X., Brizard A. and Hahm T.S. 2014 *Phys. Plasmas* **21** 122510

[23] Wang L. and Hahm T.S. 2009 *Phys. Plasmas* **16** 062309

[24] Rosenbluth M.N. and Hinton F.L. 1998 *Phys. Rev. Lett.* **80** 724

[25] Brizard A. 2011 *Phys. Plasmas* **18** 022508

[26] Stephens C.D., Garbet X. and Jenko F. 2020 *Phys. Plasmas* **27** 052504

[27] Guo Z.B., Wang Y.H. and Xu S.K. 2020 *Phys. Rev. E* **101** 03201(R)

[28] Ghizzo A., Del Sarto D., Garbet X., Sarazin Y. 2010 *Phys. Plasmas* **17** 092051

[29] Ghizzo A. and Del Sarto D. 2019 *Plasma* **2** 229

[30] Kadomtsev B.B. and Pogutse O.P. 1971 *Nucl. Fusion* **11** 67

[31] Xu S., Guo Z.B. and Gurcan A.D. 2021 *Phys. Rev. E* **103** 023208

[32] Kuramoto Y. 1984 *Chemical Oscillations, Waves and Turbulence* (Springer)

[33] Acebron J.A., Bonilla L.L., Perez Vicente C.J., Ritort F. and Spligler R. 2005 *Rev. Mod. Phys.* **77** 137

[34] Strogatz S.H., Mirollo R.E. and Matthews P.C. 1992 *Phys. Rev. Lett.* **68** 2730

[35] Martens E.A., Barreto E., Strogatz S.H., Ott E., So P. and Antonsen T.M. 2009 *Phys. Rev. E* **79** 026204

[36] Balmforth N.J. and Sassi R. 2000 *Physica D* **143** 21

[37] Sonnendrucker E., Roche J.R., Bertrand P. and Ghizzo A. 1999 *J. Comput. Phys.* **149** 201

[38] Ghizzo A., Bertrand P., Shoucri M., Fijalkow E. and Feix M.R. 1993 *J. Comput. Phys.* **108** 105

[39] Ghizzo A. and Palermo F. 2015 *Phys. Plasmas* **22** 082303

[40] Ghizzo A. and Palermo F. 2015 *Phys. Plasmas* **22** 082304

[41] Brochard G., Dumont R., Garbet X., Lutjens H., Nicolas T. and Orain F. 2018 *J. Phys. : Conf. Ser.* **1125** 012003