Nonlinear electromagnetic interplay between fast ions and ion-temperature-gradient plasma turbulence

In strong electromagnetic regimes, gyrokinetic simulations have linked a substantial ion-scale turbulence stabilization to the presence of supra-thermal particles, capturing qualitatively well the experimental observations in different devices worldwide. An explanation for the underlying physical mechanism responsible for the fast-ion-induced turbulent transport reduction observed in the numerical simulations has been proposed only recently by Di Siena et al. (Nucl. Fusion, vol. 59, 2019, p. 124001; Nucl. Fusion, vol. 60, 2020, p. 089501). It involves a nonlinear cross-scale coupling (nonlinear interaction involving different modes at different wavenumbers) between ion-temperature-gradient and marginally stable Alfvén eigenmodes, which in turn increases zonal flow activity. In view of an optimization of this turbulence-stabilizing effect, the key parameters controlling the nonlinear cross-scale coupling are here identified. At the same time, these findings provide useful insights for reduced-turbulence models and integrative approaches, which might be trained on the results presented in this paper to grasp the underlying physics and the parameter scaling of the beneficial effects of fast particles on plasma turbulence.


Introduction
A long-standing challenge in magnetic fusion research is the identification of mechanisms able to suppress turbulent transport -the drive in free energy is provided by steep gradients -with the ultimate goal of improving reactor performances. An important breakthrough along these lines is given by local gyrokinetic simulations showing a beneficial impact of supra-thermal particles -generated through auxiliary heating schemes -on ion-temperature-gradient (ITG) turbulence transport (Romanelli, Zocco & Crisanti 2010;Citrin et al. 2013;Garcia et al. 2015;Doerk et al. 2017;Di Siena et al. 2018, 2019a. Signatures of this turbulence regulation via fast particles -generated via auxiliary heating systems -have also been observed in experiments at AUG (see e.g. Tardini et al. (2007), Bock et al. (2017), JET (see e.g. Mantica et al. (2009Mantica et al. ( , 2011, It is worth mentioning that the range of supra-thermal particle parameters employed throughout this paper reflects the one typically achieved in present-day experiments with auxiliary heating systems (e.g., neutral-beam injection and ion cyclotron resonance heating). For such cases, the energetic-particle temperature is sub-critical (with energy lower than the critical energy, where the slowing-down rate on electrons is equal to the scattering rate on thermal ions). This is an important remark to be made, since the ion-scale turbulence regulation observed with such fast-ion parameters might not extrapolate to fusion-born alpha particles. To describe alpha particles correctly, radially global simulations might be required and the underlying interaction between them and ITG turbulence might change. However, sub-critical energetic-particle scenarios remain of high importance for ITER and DEMO ramp-up/-down scenarios.
This paper is organized as follows. In § 2, we describe the plasma parameters and numerical set-up employed for the gyrokinetic simulations. In § 3 we compare the impact of nonlinear electromagnetic effects in simulations with and without energetic particles, providing additional insights into the destabilization of high-frequency modes due to wave-wave nonlinear coupling. The nature of these modes is further investigated in § 4 with linear gyrokinetic simulations, while in § 5 the range of validity of the flux-tube approximation to study them is discussed. In § 6 the main dependencies of the nonlinear energetic-particle stabilizing effects on ITG turbulence are presented. Conclusions are drawn in § 7.

Numerical simulation set-up
The impact of energetic particles on turbulent transport is investigated via gyrokinetic simulations performed with the turbulence code GENE (Jenko et al. 2000). To reduce the otherwise prohibitive computational cost of these analyses the flux-tube approach is applied. It assumes periodic boundary conditions over the radial and bi-normal directions, thus enabling a Fourier decomposition of the perturbed quantities (Dannert & Jenko 2005). The range of validity of the flux-tube approximation in capturing the energetic-particle dynamics is substantiated in § 5. The numerical simulations are performed with a radial box size of 175ρ s and with a minimum k y ρ s of 0.05. Here, ρ s = c s /Ω i , where c s = (T e /m i ) 1/2 represents the sound speed, Ω i the ion gyro-frequency. The grid resolution employed in the radial (x), bi-normal ( y) and field-aligned (z) directions respectively is 192 × 48 × 32, while 32 grid points are used in the velocity parallel to the background magnetic field and 24 for the magnetic moment. Furthermore, for simplicity and in order to keep the computational expenses at a reasonable level, the fast particle species have been modelled with an equivalent Maxwellian distribution function. Future studies will take advantage of the non-Maxwellian GENE extensions discussed in detail in Di Siena et al. (2016Siena et al. ( , 2018a. However, the impact of the velocity space asymmetries and anisotropies in the background distributions is not expected to play a significant role in the range of parameters employed throughout this paper (Di Siena et al. 2018b;Di Siena 2019).
The reference plasma parameters are the same as those previously employed in Di Siena et al. (2019a) and are in the experimental range of the L-mode JET discharge studied in Citrin et al. (2013) and in Bravenec et al. (2016) at the radial position ρ tor = 0.33. Thermal deuterium, electrons and supra-thermal deuterium, with typical values of externally injected neutral beam fast ions, are considered. Collisions are retained in the simulations and modelled by a linearized Landau-Boltzmann collision operator (Florian 2008). The toroidal rotation is neglected in this paper, in order to isolate the impact of the energetic particles on ITG turbulence due only to the wave-wave nonlinear cross-scale interaction. For the same reason, the magnetic geometry, described by an analytical Miller equilibrium Here, T denotes the temperature normalized to the electron one, ω T,n = a/L T,n the normalized logarithmic temperature and density gradients, R 0 the major radius, a the minor radius, s = (ρ tor /q)(dq/dρ tor ) the magnetic shear, ρ * i = ρ s /a and β e = 8πn e T e /B 2 0 the ratio between the thermal electron and magnetic pressure. (Miller et al. 1998), has been kept fixed to the nominal one regardless of the change in the plasma parameters.
The reference values are summarized in table 1. We note that for the specific energetic-particle parameters selected -describing fast ions from neutral-beam injection heating -the wave-particle resonance interaction of Di Siena et al. ( , 2019b is here expected to be negligible. The condition requiring a larger energetic-particle temperature gradient compared to the density one is not satisfied.

Electromagnetic fast-ion turbulence suppression
Let us start by considering in figure 1 the effect that finite electromagnetic fluctuations have -for this particular set-up -on turbulent transport in the absence of energetic particles. The plasma parameters are the same as those in table 1, except for the ion density and its gradient, here set equal to the electron ones to fulfil quasi-neutrality. The nonlinear heat fluxes of each species -in gyro-Bohm units Q gB = T 2.5 i n i m 0.5 i /e 2 B 2 0 R 2 0 , with e the effective ion charge -are illustrated in figure 1 for β e = 0 (electrostatic) and β e = 0.012. They are obtained as the sum of the electrostatic and (when present) the electromagnetic (flutter) flux contributions. Interestingly, while only a mild impact of finite β e -effects is observed in figure 1(a) on the ion heat flux (reduced by 30 %), we note a destabilization of the electron heat flux (by roughly 60 %), which leads to overall larger turbulence levels (figure 1b). This increase is driven by the electromagnetic electron flux contribution, with its electrostatic counterpart only slightly affected by the presence of finite β e . These findings are consistent with the theoretical analyses reported in Kim, Horton & Dong (1993), Zonca, Chen & Santoro (1996) and Zonca et al. (1998Zonca et al. ( , 1999, where the stabilization of ITG micro-instabilities has been linked to an electromagnetic coupling between ITGs and Alfvénic ITGs/kinetic ballooning modes. Recent numerical studies have investigated this topic further, revealing an interesting interaction between stable and unstable modes and zonal flows that increases the triplet correlation time and enhances ITG transport reduction (Whelan, Pueschel & Terry 2018).
When supra-thermal particles (modelled with the parameters of table 1) are added to the nonlinear simulations, figure 2 reveals small deviations from the previous results (without energetic particles) only for the case β e = 0. A substantial reduction (up to 95 %) of the heat fluxes in each turbulent channel is instead observed when electromagnetic effects are included. In these simulations, similarly as discussed in Di Siena et al. (2019a), two distinct nonlinear phases are identified, marked by the vertical black line in figure 2. During the so-called phase I, a high-frequency modulation of the heat fluxes and slowly decaying transport levels are observed, which decrease until phase II, where they enter a stationary state at reduced levels. This transition between these two separate nonlinear phases is  systematically observed in all the turbulence simulations characterized by strong nonlinear and electromagnetic energetic-particle effects (but see also the discussion in § 6.1). The behaviour of the gyrokinetic simulations at low flux levels and/or close to marginality presented in this paper is discussed in the Appendix. The nature of the high-frequency oscillations in the turbulent transport has been already partially investigated in Di Siena et al. (2019a) and Di Siena (2019), which reveal a progressive destabilization of linearly-stable energetic-particle-driven modes with β e . In the following we extend upon the previous analyses by studying in figure 3 the frequency spectra of the electrostatic potential for each bi-normal mode number k y ρ s . The numerical diagnostic employed to extract the frequency applies the so-called windowed Fourier transform of the time trace (interpolated to equidistant time steps) of the electrostatic potential φ 1 in the time domain of the nonlinear simulations, which goes from 50 to 350 (phase I) in units of c s /a. The windowed Fourier transform enables the study of non-periodic time signals by applying smoothing functions (in this case Hanning apodization functions) and removing spurious high-frequency components.
The (k y ρ s , ω) structure of the perturbed electrostatic potential is shown in figure 3 after an average over the radial wavevectors and along the field-aligned coordinate z. The ratio between the electron and magnetic pressure is fixed at β e = 0 (electrostatic) and β e = 0.012. In all the plots of figure 3, a similar pattern of frequencies, increasing nearly linearly with k y , can be observed. This corresponds to the dispersion of the ITG mode. Only for finite β e and in the presence of fast ions (panel (c)), a second branch at higher frequency -with amplitude comparable to the ITG one -can be identified. This secondary peak is found at ω ∼ 2c s /a. It is localized in a narrow region in k y ρ s , which goes from 0.075 < k y ρ s < 0.15. In the case where energetic particles are neglected, its amplitude is strongly reduced by several orders of magnitude (note the logarithmic scale). It is important to mention that as the secondary peak arises the low-frequency ITG branch is weakened, supporting the interpretation of a nonlinear cross-scale coupling between low-and highfrequency modes. A critical requirement to observe a turbulence reduction via electromagnetic fast-particle effects in flux-tube simulations is that the high-frequency mode be linearly stable. If this condition is not fulfilled and the mode is linearly unstable, an overall increase in turbulent transport is typically observed. This is illustrated in figure 4, where the ion, electron and fast ion heat fluxes are shown in gyro-Bohm units for β e = 0.013, i.e. when the high-frequency mode is linearly unstable (the linear threshold is found at β e 0.012). A striking observation is that the electrostatic thermal (ion and electron) turbulent channels start oscillating at large amplitude around an average value close to the typical range of heat conductivities measured in present-day experiments, i.e. χ i [m 2 s −1 ] = 3.2 ± 5.1 and χ e [m 2 s −1 ] = 9.2 ± 7.5. The electromagnetic (flutter) counterpart, on the other hand, strongly increases for the electrons, while remaining at moderate values for the thermal ions. The large oscillations in the electrostatic fluxes lead to substantial uncertainties (large error bars), making any physical interpretation of these numerical simulations particularly questionable. Moreover, the energetic-particle fluxes are significantly enhanced, leading to an overall turbulence destabilization of several orders of magnitude. If the electron plasma beta is increased further, the time-averaged value of the electrostatic and electromagnetic thermal fluxes is strongly increased. These (flux-tube) findings are consistent with the results of Citrin et al. (2015) and recent radially global analyses ).

Linear stability analyses
The high-frequency mode responsible for the modulation of the turbulent fluxes observed during the nonlinear phase I is identified -for the specific range of parameters employed throughout this paper -as a toroidal Alfvén eigenmode (TAE) (Chen & Zonca 2012;Zonca & Chen 2014). Here, we provide additional evidence supporting this interpretation and demonstrate that such TAEs are linearly stable for the case under consideration. This is achieved by performing linear flux-tube GENE simulations for the bi-normal wavevector k y ρ s = 0.1 (corresponding to a toroidal mode number of 17), which lies at the centre of the high-frequency mode observed in figure 5.
Given the electromagnetic nature of TAEs we perform a scan over β e retaining the energetic particle species. The results are illustrated in figure 5, where the growth rates and frequencies of the dominant and subdominant (for β e < 0.013) modes are shown. The latter have been extracted by applying a low-frequency filter to the perturbed electrostatic potential, thus removing the slowly-varying ITG component, and measured by means of linear regression of the filtered (high-frequency) logarithmic time trace. Figure 5(a) reveals that the TAE is linearly stable for β e < 0.013, where the ITG mode represents the dominant micro-instability. However, the linear TAE damping decreases with β e until its drive term overcomes the damping and the TAE becomes the most unstable instability. This appears evident in figure 5(b) at β e ≈ 0.013, where a discontinuity in the linear frequency ω/[c s /a] is observed, indicating a mode transition. FIGURE 6. Linear growth rates (a) and frequencies (b) of the TAE for different β e and energetic-particle temperature at k y ρ s = 0.1.
Interestingly, the linearly stable TAE mode at β e = 0.012 has the same frequency as the high-frequency components observed in the turbulent spectra of the electrostatic potential in figure 3, thus enforcing the physical interpretation of Di Siena et al. (2019a) of marginally stable TAEs destabilized via nonlinear mode-to-mode coupling.
Given the key role played by the supra-thermal particles in the linear dynamics of TAEs, we extend the previous analyses by performing scans over β e for different values of the fast-particle temperature, going from T f /T e = 1 (thermal species) up to T f /T e = 13. The results are summarized in figure 6, where the linear growth rates and frequencies are shown for the range of β e where the TAEs are linearly stable. Figure 6(b) reveals that the frequency of the TAE depends not on the energetic-particle temperature but instead on the value of β e (or equivalently β i = β e n i T i /(n e T e )). Moreover, we observe very good agreement between the numerical results and the predicted values for the TAE frequency, which -in GENE normalized units -reads Here, R 0 is expressed in units of the minor radius, and β i = β e n i T i /(n e T e ) represents the thermal ion kinetic over magnetic pressure.
On the other hand, the TAE linear damping is found to decrease as the energetic-particle temperature T f /T e increases, thus moving the threshold for the linear destabilization of this mode to smaller β e (see figure 6a). Therefore, at a fixed value of β e , a more effective nonlinear coupling between TAEs and ITGs is expected (and later on confirmed in § 6) for larger values of the energetic-particle temperature T f /T e , i.e. when the TAE is still linearly stable but with a negligible damping. These results are consistent with the TAE scaling with T f /T e known in the literature (Hu & Chen 2004;Bierwage, Chen & Zonca 2009;Mishchenko, Könies & Hatzky 2009;Zonca & Chen 2014;Könies et al. 2018). An important remark here is that these findings might not extrapolate well to the range of fast-particle temperatures expected for fusion-born alpha particles, where the flux-tube approximation might break down.

Flux-tube approximation of energetic-particle modes
The degree of approximation in treating energetic-particle modes, such as TAEs, in the local flux-tube limit (1/ρ * → ∞) and its range of validity is discussed throughout this section. Results on energetic-particle modes obtained with the flux-tube gyrokinetic code GYRO (Candy & Waltz 2003a,b) have already been reported in previous papers, e.g. Bass & Waltz (2010) and Sheng, Waltz & Staebler (2017). These studies demonstrated the feasibility of local linear and nonlinear simulations of TAEs/energetic-particle modes for the DIII-D standard case parameters (Waltz, Kerbel & Milovich 1994). Here, we present a comparison between results obtained with the local version of the code GENE and well established linear global results on TAEs (related to the International Tokamak Physics Activity (ITPA) benchmark case) obtained with gyrokinetic global codes at the toroidal mode number n = 6 (Mishchenko et al. 2009;Könies et al. 2018). A detailed description of the plasma parameters can be found in Mishchenko et al. (2009) and Könies et al. (2018). They are summarized in table 2. The linear TAE growth rates and frequencies as functions of the fast-particle temperature and density are shown in figure 7. The flux-tube results are obtained at ρ tor = 0.5, which corresponds to the radial position of the peak of the global electrostatic and magnetic potential structures. Good qualitative agreement is observed for each of these simulations in both the frequencies and the growth rates at T f /T e < 400. The quantitative differences between the local and global codes (which are rather small in the frequencies) might be attributed to the absence of 'profile-averaging' as major finite size effects in the flux-tube simulation, which therefore experiences a stronger drive. Clearly, the flux-tube TAE description becomes less and less accurate as the ratio β f /β thermal = T f n f /(T e n e + T i n i ) increases. This is consistent with the interpretation of a progressive mode conversion of the TAE into an energetic-particle-mode (EPM) with β f as discussed in Mishchenko et al. (2009). In particular, a TAE has a dominant Landau damping Vannini et al. 2020), present in our flux-tube model, whereas an EPM has a dominant continuum damping (Chen & Hasegawa 1974;), which is not present here. The EPM is destabilized as the energetic-particle drive contribution overcomes the stabilizing effect provided by the shear Alfvén wave (SAW) continuum. The TAE/EPM mode conversion occurs when β f /(β e + β i ) > 1. In correspondence to this threshold, the TAE/EPM frequency exceeds the (upper) toroidicity continuum-gap frequency and starts to strongly interact with the SAW continuum. This interaction can only be captured by retaining a radially-global description, i.e. keeping phase-mixing effects. Therefore, the local flux-tube approximation breaks down as the EPM is destabilized.
However, we note that the plasma scenario considered for this paper has β f /β thermal < 0.29; i.e. the TAE mode is linearly stable and only marginally destabilized through nonlinear coupling with the ITG turbulence. Moreover, the high-frequency modes  (Mishchenko et al. 2009;Könies et al. 2018) at ρ tor = 0.5. Here, T denotes the temperature normalized to the electron one, ω T,n = a/L T,n the normalized logarithmic temperature and density gradients, R 0 the major radius, a the minor radius, s = (ρ tor /q)(dq/dρ tor ), ρ * i = ρ s /a the magnetic shear and β e = 8πn e T e /B 2 0 the ratio between the thermal electron and magnetic pressure. observed in figure 7(b) and in figure 7(d) are consistent with the TAE linear frequency of (4.1), which is significantly below the upper SAW continuum frequency and located well into the TAE gap. This is shown in detail in Mishchenko et al. (2009). Therefore, the local flux-tube gyrokinetic approximation appears to be a reasonably accurate description of energetic-particle modes which lie in the SAW gaps for the ITPA benchmark case. However, further comparisons are required to assess the range of validity of flux-tube simulations in modelling the linear and nonlinear dynamics of energetic-particle-driven modes. This will be addressed in future studies.

Nonlinear scans
In the previous sections, we have shown linear and nonlinear results supporting the physical interpretation of a nonlinear wave-wave interaction between marginally stable TAEs and ITG turbulence. Such interplay is able to explain the substantial turbulence stabilization observed in flux-tube numerical simulations (for fast ion species at sub-critical energies) and can qualitatively reproduce the experimental observations. According to the results contained in this paper and in Di Siena et al. (2019a), energetic particles can marginally destabilize modes at frequencies larger than the ITG one, which deplete the energy content of the ion-scale turbulence and act as a catalyst for an increase in the zonal-flow activity. As a result, energetic particles strongly suppress outward transport, potentially leading to improved plasma confinement. It is important to mention that this beneficial impact of energetic particles on turbulent transport requires that the linear damping of the TAEs be larger than the corresponding drive term. If this condition is not fulfilled, i.e. the EPMs are linearly unstable, a significant turbulence destabilization is observed ) (see figure 4). Typically, these modes become the dominant instabilities when the total plasma pressure and its radial gradient exceeds a particular threshold. The parameter controlling this transition is defined as α = −q 2 R dβ/dρ tor (Romanelli et al. 2010). Here, q represents the safety factor, R the major radius, β = s 8πn s T s /B 2 0 the total plasma beta (summed over all the species s) and B 0 the magnetic field on-axis. When α passes a critical value α c these modes strongly increase particle/heat fluxes, thus causing degradation of the overall plasma confinement (see figure 4). It is, therefore, essential for further exploitation of this stabilizing effect to explore the impact of the various plasma parameters on the turbulence stabilization observed previously and in the literature under the condition α < α c .
The main goal of this section is to identify the most effective parameters controlling the nonlinear turbulence stabilization. Moreover, these findings provide guidance for maximizing the turbulence suppression by nonlinear and electromagnetic fast-ion effects and well-converged data for training of reduced models to incorporate this nonlinear fast-particle effect. It is important to mention here that we consider in our analyses only supra-thermal particles generated via auxiliary heating systems, which are essentially different from alpha particles. Therefore, it may not be possible to extrapolate these results to fusion-born fast ions.
Given the nonlinear nature of this mechanism, we perform 30 nonlinear GENE simulations, each requiring approximately 100 000 CPU · h on the Marconi Skylake partition.
6.1. Energetic-particle temperature The first parameter dependence studied in this section is the energetic-particle temperature. It distinguishes the supra-thermal particles from the bulk species and, as previously discussed in § 4, directly affects the linear damping/drive of the TAE modes. The role of T f /T e on the nonlinear electromagnetic turbulence suppression observed in figure 2 is analysed by performing nonlinear GENE flux-tube simulations at β e = 0.006, keeping the other main plasma parameters and geometry fixed at the reference values summarized in table 1. In figure 8(a) we observe a substantial reduction of the main-ion turbulent fluxes as the fast-ion temperature increases, with a relative stabilization by about 40 % at T f = 13T e . The same reduction is observed for the electron heat flux (not shown here), consistently with the results of figure 2(b). We note a linear dependence (inside the error bars) of Q i /Q gB with T f /T e , which scales with a coefficient m = −5.1. Interestingly, as the main-ion heat flux decreases with T f /T e , a progressive destabilization of high-frequency modes (TAEs) is observed in the frequency spectra of the electrostatic potential at k y ρ s = 0.1 (see figure 8b). This is obtained by applying a Fourier transform to the electrostatic potentialφ 1 during the saturated steady-state time domain t[a/c s ] = [50-350]. A meaningful observation is that, consistently with the linear results shown previously in § § 4 and 5, the frequency of the TAE mode does not change significantly with T f /T e and is well captured by (4.1). We note that deviations in the TAE frequency with T f /T e are only expected for larger values of T f /T e (as shown in figure 7).
It is worth mentioning that we do not observe the transition to the second nonlinear phase -characterized by an increased in the zonal flow activity -for these specific simulations. This is attributed to the relatively small value of β e = 0.006 employed throughout these analyses, which does not allow the linear TAE drive to get close to the damping of the mode. Therefore, the nonlinearly excited TAE -sustained only via cross-scale coupling with the ITG turbulence -does not reach a meaningful amplitude to interact with the zonal flow.
6.2. Energetic-particle temperature gradient The second energetic-particle parameter analysed is the logarithmic temperature gradient ω T,f . It enters the definition of the fast-ion pressure gradient and thus affects the normalized pressure α and the TAE linear drive. It is varied from ω T,f = 0, i.e. a flat temperature profile, to ω T,f = 1.5, i.e. 50 % steeper than the reference value. The other plasma parameters are summarized in table 1, with the sole exception of the electron plasma beta, which is β e = 0.009 (increased with respect to the previous nonlinear scans to move the TAE closer to the marginal stability threshold). The magnetic geometry is kept fixed in all of these analyses, so that the results within this section arenot related to geometric stabilizing effects enhanced by the increase of the fast-ion pressure gradient (Bourdelle et al. 2005).
The behaviour of the bulk ion turbulent heat flux (Q i /Q gB ) in relation to the energetic-particle logarithmic temperature gradient is illustrated in figure 9(a). A progressive enhancement of the nonlinear fast-ion stabilization is observed when ω T,f is increased, leading to 60 % lower turbulent fluxes at ω T,f = 1.5 (with respect to ω T,f = 0). We note a linear dependence of Q i /Q gB on ω T,f , which scales with a coefficient m = −12.6. The growth of the nonlinear fast-ion stabilization with ω T,f is correlated with the destabilization of linearly stable TAEs. This is shown in figure 9(b), where the frequency spectra of the electrostatic potential are plotted for the bi-normal mode number k y ρ s = 0.1 -at the middle of the k y ρ s scales where the TAE is nonlinearly destabilized. Interestingly, ω T,f does not affect the TAE frequency, but rather the linear TAE drive/damping, thus moving the threshold for linearly destabilizing this mode -at fixed β e -to smaller values of ω T,f , or equivalently α. This result is similar to what has been observed previously for the energetic-particle temperature. Therefore, T f /T e and ω T,f play similarly critical roles in the dynamics of the electromagnetic cross-scale interaction studied here. However, the logarithmic temperature gradient is more effective than the temperature in suppressing ITG turbulence by more than a factor of two (m ω T,f /m T = 2.5), as estimated by looking at the coefficient of the linear fit of figures 8(a) and 9(a).
The characteristic transition to the second nonlinear phase is observed in all the simulations presented in this section. This feature was not found during the scan over the energetic-particle temperature because of the smaller value of β e employed, making the nonlinear electromagnetic fast-particle effects on turbulence weaker. This can also be noted by looking at the relatively larger ion heat fluxes observed during the scan over T f /T e .

Energetic-particle density gradient
Similarly to ω T,f , the energetic-particle logarithmic density gradient affects the normalized pressure gradient α and thus the TAE linear drive. In contrast to ω T,f , a change in ω n,f imposes a variation in the logarithmic density gradient of the bulk ions to ensure quasi-neutrality. The latter is always fulfilled in our analyses on both density and density gradients. The plasma parameters employed here are summarized in table 1, with the sole exception of the electron plasma beta, which is β e = 0.009. The dependence of the main-ion turbulent fluxes on the fast-particle density gradient is shown in figure 10(a). Again, an increase in ω n,f reflects into an enhancement of the nonlinear electromagnetic fast-particle stabilization, which is always accompanied by a destabilization of a linearly stable TAE mode. However, we note that the energetic-particle stabilization scales with the fast-ion logarithmic density gradient with a coefficient m ω n,f = −7.4, twice as small as the temperature gradient one. This result can be explained by an interplay between the wave-particle resonance interaction and the nonlinear wave-wave coupling. The increase in the supra-thermal ion temperature gradient leads to an enhancement of the resonant interaction between fast particles and the underlying ITG micro-instability. However, since ω n,f > ω T,f this resonant interaction destabilizes the bulk ion-driven ITG, leading to a larger ITG drive and thus larger turbulent fluxes. As a result, the critical threshold to destabilize the TAEs (α c ) is expected to move to larger values when ω n,f > ω T,f . The results shown in figures 9(a) and 10(a) seem to confirm this theoretical prediction, although the nonlinear nature of this fast-particle turbulence stabilization makes any extrapolation on the α c questionable.
These results suggest that ω n,f is less effective than ω T,f in enhancing the overall fast-particle turbulence suppression, or, equivalently, that the wave-particle interaction weakens the TAE drive. A more detailed analysis of the possible interaction between these two newly discovered fast-particle effects on turbulence is left for future study.

Safety factor
Another important parameter controlling the stability of the sub-marginal TAEs is the safety factor q. In a radially global set-up, it strongly affects the continuous spectrum of SAW and, hence, the radial position of the continuum gaps. Therefore, it impacts both the TAE frequency (see (4.1)) and the normalized pressure α (as α ∝ q 2 ), simultaneously affecting frequency and drive/damping. Assuming that the critical threshold to destabilize these modes α c is approximately constant with respect to the safety factor, an increase of q will lead to a strong enhancement of the nonlinear fast particle turbulence suppression at smaller values of β e . This behaviour is confirmed by the nonlinear GENE simulations in figure 11(a), where a scan over β e at q = 1.2 is compared to the one obtained at q = 1.7. These results show that the nonlinear electromagnetic turbulence suppression by supra-thermal ions is particularly sensitive to the change in the safety factor q. More precisely, an increase in the safety factor from q = 1.2 to q = 1.7 leads to a downshift of 54 % in the critical value of β e,c : it moves from β e,c = 0.013 at q = 1.2 to β e,c = 0.006 at q = 1.7. These findings suggest a quadratic dependence of the critical electron beta (β e,c ) on the safety factor, namely β e,c (q = 1.7)/β e,c (q = 1.2) ∼ (q = 1.2) 2 /(q = 1.7) 2 . A more accurate quantification of the scaling power of β e,c with q will be carried out in the near future by increasing the number of nonlinear simulations performed. The enhanced turbulence suppression obtained by increasing the safety factor q at fixed β e is again attributed to the nonlinear destabilization of linearly stable TAEs, as shown by the frequency spectra of the electrostatic potential ( figure 11b). An increase in q leads to a larger TAE drive and to a more effective cross-scale nonlinear coupling, as demonstrated in figure 11 by the behaviour of the amplitude of the high-frequency mode with q. The change in the safety factor also affects the marginally stable TAE frequency, which scales as predicted by (4.1), as 1/q, going from ω = 2.4[c s /a] at q = 1.2 to ω = 1.7[c s /a] at q = 1.7.
These findings highlight the importance of the safety factor in controlling the nonlinear fast-particle turbulence stabilization and make q a key parameter for the design of experiments with improved confinement.

Magnetic shear
The magnetic shear (defined as s = (ρ tor /q)(dq/dρ tor )) -similarly to the safety factor q -is another important geometrical parameter affecting the energetic-particle nonlinear electromagnetic turbulence suppression. While the safety factor q affects the normalized pressure gradient α, the magnetic shear s affects the stability threshold for the destabilization of AEs α c ; see e.g. Fu (1995). Therefore, at a fixed α, variations in the magnetic shear can enhance or weaken the nonlinear electromagnetic fast-ion turbulence suppression, depending on the difference between α and α c . The smaller this difference, the closer the sub-marginally stable fast-particle-driven TAE is to the linear stability threshold, leading to a more effective cross-scale coupling and to a stronger turbulence suppression. However, if α surpasses α c , a strong turbulence destabilization is observed, which is attributed to linearly unstable energetic-particle-driven modes ) (see figure 4). These predictions are consistent with the results presented in figure 12(a), where the main-ion turbulent fluxes are shown for different values of β e and magnetic shear. More precisely, figure 12(a) reveals that if s is increased from s = 0.1 to s = 0.52 (nominal) and s = 0.75, the nonlinear fast-particle turbulence suppression is progressively weakened, although α is unaffected. We note that the increase in the magnetic shear affects the critical β c (and hence α c ), shifting it to larger values. An increase in the magnetic shear of 44 % moves β c (or equivalently α c ) by 12 %, namely from β c (s = 0.52) = 0.006 to β c (s = 0.75) = 0.0067. These findings reveal that the magnetic shear has a minor impact on the critical threshold α c . However, it still has a large effect in enhancing the nonlinear wave-particle stabilization, which increases by 85 % at β e = 0.006 for the value of the magnetic shear s = 0.52 (as shown in figure 12a), relative to s = 0.75. Furthermore, consistently with the previous results, figure 12(b) shows that the amplitude of the TAE mode and the ITG are strictly connected; specifically, as the amplitude of the TAE increases, the ITG peak is consequently reduced. This again suggests a nonlinear cross-scale coupling between these different plasma instabilities. It is worth mentioning that for the case s = 0.1, a third peak (ω[c s /a] ∼ 2.31) arises at frequency larger than the TAE one (ω[c s /a] ∼ 1.38). This finding suggests that the nonlinear cross-scale between ITGs and marginally stable modes is not restricted to TAEs only. Similar effects are expected for other marginally stable modes driven by thermal and/or energetic ions.
The findings presented throughout this section are consistent with previous modelling results (Citrin et al. 2013 and experimental observations (Mantica et al. 2011), showing enhanced confinement in hybrid scenarios at JET due to ITG electromagnetic stabilization at low magnetic shear (Mantica et al. 2011;Citrin et al. 2013).
Hybrid scenarios (high beta, low magnetic shear in the inner half-radius and significant fast-ion fraction) are high-confinement regimes achieved in present-day devices (e.g. JET, AUG, DIII-D, KSTAR) as well as in ITER extrapolations. In such regimes, electromagnetic stabilization effects are considered to be an important contribution to the observed improved core confinement regime. Therefore, our findings, showing an improved turbulence stabilization via fast ions at low magnetic shear (e.g. s < 0.2), might be particularly valuable for such conditions.

Conclusions
Traditionally, fusion theory has largely been subdivided into separate research topics treated independently. More recently, however, we witness the emergence of a growing number of cases in which two or more fundamental processes with disparate time and/or space scales are inseparably connected, adding new facets to the qualitative and quantitative understanding of the respective phenomena. The present paper addresses one such case. Recent studies of magnetically confined plasmas have highlighted that turbulent transport can be controlled by the presence of energetic ions, often neglected in turbulence studies. This calls for a more integrative modelling approach, with the goal of grasping the underlying physics and identifying the range of parameters that maximize the beneficial effect of fast particles on magnetic confinement.
In the present contribution, significant steps along these lines are presented through flux-tube gyrokinetic simulations of realistic plasma conditions. More precisely, the impact of energetic particles on ITG turbulence is investigated in detail. It has previously been found, in a number of studies, that supra-thermal particles can strongly reduce turbulence levels, leading to a significant suppression in the outward turbulence fluxes for each plasma species. The mechanism responsible for such a substantial effect has been attributed to the cross-scale nonlinear interaction between linearly stable Alfvén eigenmodes (driven by the large fast-ion pressure and pressure gradients) and the ITG turbulence (Di Siena et al. 2019a).
To further corroborate these results and understand the range of validity of the flux-tube approximation in modelling EPMs, we have presented a linear benchmark of TAEs (the ITPA benchmark case) between the local version of the code GENE and well established linearly global results. We have found that, for low energetic-particle temperature and far from the SAW continuum gap, the flux-tube code successfully captures the TAE linear growth rates and frequencies. On the other hand, it breaks when the supra-thermal particle temperature increases and the interaction with the continuum is not negligible, which, however, is well beyond the parameter range of interest for the experimental parameters at hand. Furthermore, we have performed a series of nonlinear simulations to address the scaling of the electromagnetic fast-particle turbulence suppression with different plasma parameters. Among the fast-particle ones (T f /T e , ω T,f and ω n,f ), we have found that the most effective in enhancing the ITG turbulence suppression is the logarithmic temperature gradient, which leads to a sharp drop in the main-ion turbulent fluxes. The magnetic geometry is also found to lead to strong modifications in the turbulence suppression. In particular, we have studied the role of the safety factor and the magnetic shear in controlling the levels of the turbulent fluxes. Both these parameters have a strong impact, moving the marginal stability threshold of the energetic-particle-driven TAEs and thus enhancing their beneficial effect on turbulence (acting as an effective sink of free energy from the ITG micro-turbulence).
Beside providing physical understanding of the nonlinear interaction between energetic particles and ITG micro-turbulence, these findings provide useful insights for reduced turbulence models and integrative approaches, which still miss this particular effect .

Acknowledgements
The authors would like to thank P. Lauber for useful discussions.
Editor Francesco Califano thanks the referees for their advice in evaluating this article.

Funding
This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement no. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The simulations presented in this work were performed at the Marconi CINECA. Therefore, we acknowledge the CINECA award under the ISCRA initiative for the provision of high-performance computing resources and support.

Declaration of interests
The authors report no conflict of interest.

Appendix. Long-term zonal flow evolution in near-marginal turbulence
Certain turbulent systems that approach the nonlinear ITG threshold have recently been shown to exhibit a long-time evolution of the zonal flow amplitude Rath et al. 2016). In particular, a minimum main-ion heat flux value of Q i = 10 in gyro-Bohm units has been identified for the collisionless electrostatic CYCLONE base case parameters with adiabatic electrons. Below this value, Rath et al. (2016) showed that no turbulent state can exist, and periodic bursts of the main-ion heat flux have been observed in association with a secular evolution of the zonal flow amplitude. Therefore, the nonlinear simulations presented throughout this paper and in Di Siena et al. evolution of the zonal flow to a significant increase in the non-Gaussian features of the heat flux p.d.f.s. In particular, values of 3.4 and 18, respectively, for the skewness and excess kurtosis were observed in Rath et al. (2016) in the low-heat-flux simulations, characterized by non-physical zonal flow evolution and bursty transport. Therefore, the results presented in this section reveal that only minor deviations from the Gaussian p.d.f. distribution are observed in this paper for each value of β e . These observations suggest that the zonal flow increase between the so-called phase I and phase II of Di Siena et al. (2019a) cannot be explained in terms of a non-physical behaviour of the long-time zonal flow evolution.
Moreover, a stable turbulent 'stationary phase' is observed for a rather long time domain, namely t ∼ 400-1000[a/c s ], which is in contrast with the interpretation of a secular/intermittent evolution of the heat flux and zonal flow. Moreover, as discussed in Di Siena et al. (2019a), after a Fourier decomposition in time of the zonal component of the electrostatic potential, fast oscillations are observed at the specific TAE frequency, corroborating the physical interpretation given throughout this paper and in Di Siena et al. (2019a). It is remarked that the magnitude of the heat flux during the second nonlinear phase depends on the simulation parameters and does not necessarily imply proximity to marginality, while the nonlinear coupling to a marginally stable high-frequency mode remains a ubiquitous observation in all these scenarios.
It is worth mentioning that the results of this paper were obtained in a realistic plasma scenario, which includes the effects of inter-species collisions, kinetic electrons, electromagnetic effects and complex magnetic plasma geometry. Therefore, each of these terms can affect the minimum heat-flux threshold for the non-physical results discussed in Peeters et al. (2016) and Rath et al. (2016). In particular, collisions can act on the zonal flow evolution as a natural damping mechanism, suppressing or delaying its long-time secular growth.