Gyrokinetic study of transport suppression in JET plasmas with MeV-ions and toroidal Alfvén eigenmodes

The impact of fast ions, generated in the MeV-range through the efficient application of the three-ion scheme in JET plasmas, on the turbulence properties is presented through complex numerical simulations. The suppression of the ion-scale turbulent transport is studied by means of in-depth gyrokinetic numerical analyses. Such a suppression is demonstrated to be achieved in the presence of toroidal Alfvén eigenmodes (TAEs) destabilized by the highly energetic ions. Details on the TAE excitation are also provided with a multi-code analysis. The inherently nonlinear and multi-scale mechanism triggered by the fast ions, also involving the high-frequency modes and the large-scale zonal flows, is deeply analyzed. Such mechanism is thus demonstrated, with experimental validating studies, to be the main cause of turbulence suppression and improvement of ion thermal confinement. Additional simulations address the implications of reversed shear magnetic equilibrium on the turbulent transport.


Introduction
The success of magnetic confinement fusion as an energy source crucially relies on reaching simultaneously high temperatures and high densities for the fuel D and T ions. Fusionborn alpha particles are the main source of central plasma heating in ITER and future fusion power plants. Yet, these highly energetic alpha particles heat primarily electrons rather than bulk ions through Coulomb collisions [1]. The physics of plasma heating by alpha particles is complex and its extrapolation to future devices is not straightforward, partly because of the mutual interplay between turbulence and energetic ions. This can lead to additional nonlinearities, further increasing the complexity of plasma dynamics. As an example, a significant reduction of the ITG-induced transport was recently predicted for ITER D-T plasmas in the presence of fusion-born alpha particles [2].
Simultaneously with plasma heating, alpha particles can also excite instabilities, in particular, the non-damped Alfvén eigenmodes (AEs) [3]. The excitation of these modes by energetic ions can be related to the particular geometry of toroidal fusion plasmas, such as the toroidicity (toroidicity-induced AEs or TAEs) [4,5]. Expected to be destabilized also in ITER [6], the presence of TAEs is commonly considered detrimental for plasma confinement in fusion plasmas, as they can cause increased transport of energetic ions [7][8][9][10].
Earlier studies have already demonstrated that the thermal ion confinement is improved in the presence of fast ions introduced through neutral beam injection (NBI) or ion cyclotron heating resonance (ICRH) systems [11][12][13][14][15][16] (for comprehensive reviews on the topic, the reader can consult [17,18]). Insights on such a turbulence reduction have been provided [19], pointing to the important role of marginally stable fastion-driven modes acting as intermediate actors between the ion-scale turbulent energy and the zonal perturbation of the electrostatic potential. Nevertheless, a possible complex threewave interaction between MeV-range ions, fast-ion driven AEs and microturbulence has not been systematically studied at ITER-relevant conditions, neither theoretically nor experimentally.
Controlling turbulence in a fusion power plant will finally lead to a more economical operation and ultimately to a reduction in the cost of fusion electricity. This paper reports on the three-wave mechanism of turbulence suppression in the presence of MeV range fast ions and fast-ion driven TAEs, identified during detailed transport analysis of JET experiments in D- 3 He plasmas. The improved confinement in JET plasmas with MeV-range fast ions was confirmed with state-ofthe-art turbulence gyrokinetic modeling. Electrostatic thermal ion energy fluxes were suppressed due to intense poloidally directed shear flows, known as zonal flows (ZFs) [20]. These ZFs act as a barrier by shearing and de-correlating [21] for the outward radial transport of particles and energy, thus increasing the thermal insulation of the central plasma. Complex gyrokinetic simulations show that ZFs are generated as a result of the non-linear interaction between fast ions, TAEs and microturbulence.
The results reported in this paper highlight that AEs excited by highly energetic ions can be, in fact, beneficial to improve thermal energy confinement, a concept that has not been extensively studied before. A detailed understanding of this complex interplay paves the way toward controlling turbulence and enhanced performance of future fusion reactors with strong alpha particle heating, such as the forthcoming D-T campaign in ITER.
Section 2 is dedicated to present the experimental outcomes useful to validate the following gyrokinetic numerical analyses, which are reported in sections 3 and 4. Whereas in the former a linear study of the selected plasma stability with the identification of the destabilized modes is carried out, in the latter section the results of the nonlinear analyses are described. In section 5, additional studies on the impact of the negative magnetic shear on the transport are performed. The conclusions and future perspectives are drawn in section 6.

Three-ion scheme scenario at JET: improved thermal ion confinement
The selected pulse #94701 of the recently developed threeion D − D NBI − 3 He scenario at JET [22] has been used as the background plasma for the following gyrokinetic simulations. The three-ion scheme [23][24][25] has been adopted in a novel combination to efficiently accelerate the NBI-generated deuterons up to MeV energy range with ICRH waves. Such deuterons are indeed injected with the NBI in the vicinity of the ion-ion hybrid layer, where the mode conversion and the energy absorption are enhanced. The radial position of such a layer can be determined by means of the density ratio of the three ion species in the plasma [23]. With 3 He concentration X[ 3 He] ≈ 25%, the mode conversion layer is located in the very deep core of the plasma [22,26], and an efficient energy absorption by the NBI deuterons takes place. Further analyses with the TOFOR [27] and the neutron camera [28] diagnostics demonstrate the effectiveness of the three-ion scheme and the core-localization of the MeV-range fast D [29] in this particular JET scenario. These particles at larger energy with respect to the critical energy (∼500 keV) lead to a predominant electron heating [1] in the plasma core, with electron temperature reaching ∼7 keV.
An observed direct effect of the efficient application of the three-ion scheme, and the subsequent generation of MeVrange ions in the plasma core, is the strong Alfvén activity systematically measured. Indeed, in the frequency range Mirnov coil spectrogram measurements in the frequency range 75-260 kHz for JET pulse #94701. Beyond t = 7.4 s, when the three-ion scheme heating is applied, a series of AEs are destabilized. Toroidal mode numbers are reported close to the corresponding signals, and the TAE antenna [33,34] measurement, useful to probe stable AEs, is also displayed. f = 185-225 kHz TAEs are destabilized, as shown in the magnetic spectrogram of figure 1 for pulse #94701. X-mode correlation reflectometer [30,31] measurements of JET pulse #95669, performed at very similar experimental conditions of #94701, determine that the radial localization of the unstable TAEs roughly corresponds to 0.21 ≲ ρ tor ≲ 0.37 (3.21m ≲ R ≲ 3.36 m) [32], where ρ tor is the square root of the normalized toroidal magnetic flux.
It must be however noted that the energy distribution of the suprathermal population computed by means of the integrated modeling TRANSP [35] departs from the isotropic Maxwellian distribution, as expected for the fusion-born alpha particles in burning plasmas. Indeed, co-passing fast ions (with v ∥ > 0) are mainly generated in the three-ion scenario at JET. Although ICRH systems mainly produce ions with large perpendicular velocities, the substantial deposition in the deep core where particles primarily experience passing trajectories creates these particular conditions [22]. This consideration may be also related to the destabilization of the reversed shear Alfvén eigenmodes (RSAEs), systematically observed in the D − D NBI − 3 He scenario at JET [36]. In fact, the plasma current generated by the MeV-ions may be the cause of the nonmonotonic q-profile, necessary to destabilize the RSAEs.
Another striking observation, together with the fast-iondriven AEs and the strong electron heating due to the MeVions, is the improved global confinement when the three-ion heating scheme is applied [22,26,32]. Several experimental indicators, such as the ion temperature and the plasma energy content, are strongly enhanced if compared to very Frequency spectrum of the density fluctuation measurements with the X-mode reflectometer for JET pulse #95669 at R = 3.25 m, corresponding to ρtor ≈ 0.24. In red, the fluctuations are measured in the phase where only NBI is applied as external heating system, whereas the blue curve illustrates the second phase of the discharge with ICRH combined in the three-ion scheme. similar experimental conditions with only NBI heating. This improved confinement from direct experimental measurements is exemplified in figure 2. The amplitude of the density fluctuations measured by means of the X-mode reflectometer is illustrated. A strong difference is evident comparing two phases of the same pulse, with a significant decrease of the fluctuation amplitude especially in the low-frequency range when MeV-ions destabilize the TAEs. The spikes at f ≈ 200-240 kHz are a clear demonstration of the TAE destabilization. Moreover, the outcomes of the TRANSP integrated modeling are consistent with such a picture. Indeed, it is clearly observed that concomitantly with the destabilization of the TAEs, the improvement of the ion thermal confinement, measured as a decrease of the ion heat diffusivity in the plasma core, is obtained [32]. For more details on this experimental comparison the reader is addressed to [22,26,32]. Such an unexpected observation of improved thermal confinement in the presence of unstable TAEs and predominant electron heating in the core motivated the in-depth transport analyses, whose results are reported in the following sections of this contribution.
The achieved experimental conditions in the three-ion scenario at JET, i.e. strong electron heating from MeV-range particles, T e ∼ T i and low plasma toroidal rotation are of paramount importance in view of ITER burning plasmas [26]. Indeed, the MeV-ions generated through the efficient threeion scheme in these ITER-relevant features allow to mimic the impact of the alpha particles on the turbulence regimes and properties. Moreover, the TAEs are expected to be destabilized also in ITER D-T plasmas [6], addressing thus also the possible effect of the TAEs on the global confinement. Although the direct extrapolation to ITER plasmas is extremely challenging, the experimental observations closely together with the detailed gyrokinetic numerical simulations presented in the following suggest an unforeseen beneficial impact of the alpha particles on the turbulent transport.

Linear stability analysis
In this section, the linear stability of the selected JET pulse #94701 of the three-ion D − 3 He − D NBI scenario is analyzed. The state-of-the-art gyrokinetic code Gene [37] is employed in its local gradient-driven version. In figure 3, the experimentally measured electron density and both electron and ion temperature profiles of the JET pulse #94701 at the time t ≈ 9.6 s are shown in the radial domain. The input parameters for the gyrokinetic modeling are provided by the integrated modeling TRANSP, and are reported in table 1.
In the local version, the Gene code simulates only a fluxtube of the entire plasma volume. In this particular case, the simulation domain is centered at the radial position ρ tor = 0.23. This location has been chosen since (i) the combination of experimental measurements [29] and integrated modeling calculations [38] reveals that a substantial population of core-generated MeV ions is still present, (ii) reliable charge exchange measurements of the ion temperature profile are obtained [22], and (iii) TAEs are expected to be destabilized in this radial range, as shown in the previous section. Regarding the numerical setup, the distribution functions of up to four different kinetic species are evolved in time in a 5D numerical grid of 256 radial (k x ) and 48 binormal (k y ) wavenumbers, 32 points in the parallel direction (z), 48 points in the parallel velocity (v ∥ ) dimension and 64 points in the magnetic moment (µ) direction. The different species are electrons, thermal deuterium (majority ions), 3 He (minority ions) and fast deuterium from NBI as fast ions. The high resolution employed in µ is necessary to capture possible resonant interactions occurring at large energies. This is even more important in the case of the employed non-equidistant Gauss-Legendre discretization of the magnetic moment dimension in Gene. This consideration allows to relax the resolution in µ down to 16 when the fast ions are not retained in the simulations. Convergence tests, also in the nonlinear regime, have been carried out to ensure the correctness of the grid resolution. The minimum binormal wavenumber is k y,min ρ s = 0.025, with ρ s the thermal ion Larmor radius at the sound speed. The magnetic geometry is a EFIT equilibrium constrained by the total pressure computed by TRANSP, and a Miller parametrization [39] of such equilibrium is employed in the flux-tube gyrokinetic simulations. Fully electromagnetic simulations are performed, including hence both parallel and perpendicular fluctuations of the vector potential, and collisions are modeled through a Landau-Boltzmann operator. The externally induced plasma rotation is neglected due to low-NBI input power. An equivalent-Maxwellian fast-ion distribution function has been employed in Gene. Earlier studies showed that a non-negligible impact of realistic fast particle distribution on the ion-scale turbulent transport level may occur. Nevertheless, the qualitative results are comparable [40,41]. It is to be noted that the following analyses are firstly carried out for a time window before the onset of the RSAEs. A further study, reported in the following section 5, addresses the impact of the reversed shear on the transport.
As a final remark, it must be stressed that for this kind of studies the local (or flux-tube) approximation may not be entirely accurate due to the large value of the fast deuterium Larmor radius compared to the minor radius, i.e. ρ * FD ≈ 1/90. For this reason, this study is intended to be as a first step, whereas global gradient-driven simulations are planned to be performed in the near future to assess and verify the results here reported.
In figure 4, the linear spectra computed by Gene are illustrated. The growth rate is normalized to c s /a, where c s = T e /m p , with T e the electron temperature and m p the proton mass, and a is the minor radius. The blue curve represents the case without including the fast ions. The ITG instability peaking around k y ρ s = 0.45 is dominant, since the positive (resp. negative) frequencies mean that the mode is propagating in the ion (resp. electron) diamagnetic direction. In the low-k y region fast-ion-driven modes are destabilized beyond a threshold in the normalized fast ion pressure gradient R/L pFD . For the set of input parameters described in table 1, the threshold is identified to be around R/L pFD = 10.7. This latter case is represented by the green curve in figure 4, and it will also be named marginally stable case in the remainder. For the cases with R/L pFD = 13.6 and 16.2, the fast-ion-driven modes are fully destabilized (in figure 4, for the sake of simplicity, only the case with R/L pFD = 16.2 is illustrated). Although not dominant over the whole binormal spectrum, their growth rates are not negligible. Notably, the frequencies of these modes are within the experimentally measured range of frequencies of the TAEs ω ≈ 185-235 kHz [22], as also highlighted in  Table 1. Employed plasma parameters in Gene simulations modeling JET pulse #94701 at ρtor = 0.23 and t ≈ 9.6 s. Here, R 0 is the major radius, B 0 the on-axis magnetic field, Ip the plasma current, P NBI the NBI injected power, P ICRH the ICRH injected power, ne and Te respectively the local values of electron density and temperature, ϵ represents the inverse aspect ratio, R/L n,T the normalized logarithmic density and temperature gradient, βe the electron-beta, and ν * the normalized collision frequency. The ratios n i /ne and T i /Te, with i indexing the species, are the relative density and temperature with respect to the electron ones. The reported input parameters are common to all the numerical Gene simulation cases. The various cases, however, differ essentially in the fast ion pressure gradient R/Lp FD .  figure 4 with a gray shaded area. More details on the identification of these fast-ion-driven modes are provided in section 3.1. In the small scale region, ETG modes are found unstable up to k y ρ s ≈ 30 due to a large value of electron temperature gradient, as illustrated by the blue curve in figure 5. Although electron-scale instabilities may affect the ion-scale transport [42][43][44], recent gyrokinetic studies on JET L-mode plasmas showed a negligible impact of ETG modes on the ITG turbulence [45][46][47]. Moreover, the rule of thumb proposed in [48] is calculated and compared to the growth rate spectrum in figure 5. The criterion is based on the evaluation of the ratio between the maximum of the ETG and ITG linear growth rates, both divided by the binormal wavelength of the peaks. This parameter is here defined as represents the peak of the ETG (ITG) spectrum. With τ crit > 1, ETG instability is expected to affect the total electron transport, whereas for τ crit < 1 a negligible effect is generally measured. Although the maximum ETG growth rate is larger than the ITG one, the ETG region peak for the parameter of merit γ/k y is lower than the one for the ITG region, with τ crit = 0.42, as shown in figure 5. This justifies the choice to set the upper boundary of the binormal wavenumbers to k y ρ s = 1.2 and to exclude, thereby, the ETG modes in the following nonlinear simulations. Further studies are planned to be performed also retaining the electron-scale dynamics in the near future.
As can be seen in figure 4, the spectra without and with fast ions (also for different values of R/L pFD ) almost overlap. This suggests that no linear fast-ion resonant effect on the ion-scale instability is present. This is consistent with earlier analyses [15] showing a negligible effect for large values of T FD /T e , similarly to the value in this three-ion JET scenario.

Identification of the fast-ion-driven mode nature and excitation mechanism
A detailed study on the identification of the fast-ion-driven modes destabilized at large fast ion pressure gradient in the low-k y region is described in this subsection. The analyses focus on the binormal wavenumbers k y ρ s = 0.025 and 0.05, corresponding to the toroidal mode numbers n = 4 and 7. It must be noted, for the sake of clarity, that in order for the lowest-mode rational surface to correspond to an integer toroidal mode number n, a slight radial shift of the binormal wavenumbers is considered in the linear simulations. Hence, the actual binormal wavenumbers corresponding to n = 4 and 7 are, respectively, k y ρ s = 0.0269 and 0.0474. However, in the following, for consistency with the numerical binormal grid employed in the nonlinear simulations, we will refer to these wavenumbers as k y ρ s = 0.025 and k y ρ s = 0.05. Including the lower toroidal mode numbers would have required an amount of CPU resources out of the scope of the work. Moreover, as can be seen in figure 1, the n = 1 and 2 TAE modes are expected to be stable.
Motivated by the good matching of the Gene linear mode frequencies of the fast-ion-driven modes with the experimental TAE range of frequencies, the Alfvén continuum is calculated using equation (8) of [49] and thus taking into account only the impact of the toroidicity. This continuum is calculated fixing the plasma parameters to those for ρ tor = 0.23, not varying along the radial coordinate. However, consistently with the flux-tube approximation employed in this gyrokinetic study, the safety factor profile is linearized around the magnetic surface of interest in the following fashion [50]: with x the radial coordinate, q 0 ,ŝ 0 and x 0 respectively the values of the q-profile, the magnetic shear and the radial coordinate at ρ tor = 0.23. In this approximation, the mode frequencies computed by Gene lie very close to the the toroidicity-induced gaps in the continua. The poloidal mode numbers of these gaps are respectively m = [−5, −4] and m = [−8, −7] for n = 4 and n = 7. Indeed, the gaps induced by the toroidicity are due to a coupling between consecutive poloidal mode numbers [49], and hence two m are given for each gap. This latter result is consistent with the poloidal structure of the modes in Gene by projecting the flux-tube domain on the poloidal section for the simulation with R/L pFD = 16.2. Moreover, as already reported in an earlier publication [32], the radial localization of the unstable TAEs has been measured in the three-ion scenario at JET to be in the range 3.22 m ≲ R ≲ 3.36 m by means of the X-mode reflectometer. Consistently, the outer mid-plane of the Gene flux-tube domain is within this radial range, exactly at R = 3.27 m.
For all the evidence provided above, the fast-ion-driven modes can be firmly identified as TAEs, with a qualitative agreement with the experimental measurements. Although a more accurate description of these large-scale instabilities would require a global approach, the flux-tube approximation was demonstrated to be adequate for a good qualitative evaluation of their impact on the microturbulence properties [19,[51][52][53].
The focus is now moved to the excitation mechanism of the TAEs, which can be addressed by analyzing the evolution of the free energy in the system. The free energy E f can be written as the sum of the kinetic E k and the potential energy E w , and their time derivatives represent the energy balance equation, determining the direction of the energy flow in the system [54]: with B 0 the magnetic field, n 0,s and T 0,s the relative density and temperature of the species s, F 0,s and f s the equilibrium and the perturbed distribution functions of the species s and ϕ the electrostatic potential. Hence, the linear growth rate can also be expressed as [55]: where the contribution of each species s is thus separated. For more details on such a derivation from the Vlasov-Poisson system of equations and the application to the Gene code, the reader can refer to [56][57][58]. For the sake of clarity, the potential energy can be associated with the energy stored in the perturbed part of the electrostatic potential fluctuations. Therefore, E w is associated with the instability energy. In Gene simulations, the free energy balance and subsequently the contribution of each species to the linear growth rate can be evaluated in the velocity space [15,57], allowing to determine resonant structures stabilizing or destabilizing the mode depending on their sign (negative values mean stabilization and positive values destabilization). In the following, the fast-ion contribution to the instability is referred to as γ FD . Thus, γ FD is evaluated at the outer mid-plane (z = 0 in Gene convention) and then plotted in the velocity space. This is done in figures 6(a) and (b), respectively for the toroidal mode numbers n = 4 and 7, for the case with R/L pFD = 16.2. Similar structures are observed also for the case with R/L pFD = 13.6. The energy flowing from the particles to the wave is highlighted with positive values (destabilizing), whereas negative values are denoting the opposite energy transfer (stabilizing). For this particular analysis, the parallel velocity resolution has been increased to n v ∥ = 120 in order to resolve any possible resonance more accurately. Two dominant positive structures at large perpendicular velocity and roughly extended in the low-v ∥ domain are clearly visible in the trapping region, where the trapped/passing boundary is illustrated by the black dashed line. Together with these resonances, some other minor structures related to the barely passing and fully passing fast ions are also induced.
To check the validity of such dominant resonances within the trapping cone, the TAPaS code [59,60] has been employed. The TAPaS code integrates the motion's equations of particles, tracing their trajectory and calculating thus the characteristic frequencies of their motion. For this particular study, TAPaS is employed in its guiding-center approximation, integrating the unperturbed equations within the gyro-kinetic framework for a substantial number of tracers of fast deuterium. The tracers are initialized on the equatorial mid-plane of the low-field-side at the magnetic surface ρ tor = 0.23 of JET pulse #94701. The fixed magnetic equilibrium is the same used in Gene, and consistently the safety factor profile has been linearized around the selected magnetic surface as done in Relation 1. In figure 6, the resonant condition f TAE − nΩ prec = 0 between the TAE frequency computed by Gene (f TAE ) and the toroidal precession frequency (Ω prec ) calculated by TAPaS is illustrated with solid black curves for both n = 4 and n = 7 wavenumbers. It is to be noted that the tracers initialized with v ∥ > 0 are all passing, when considering their unperturbed trajectories. This asymmetry in the parallel velocity space is due to the odd parity of the toroidal canonical momentum P φ with respect to the toroidal velocity v φ , which can be approximated with the velocity parallel to the magnetic field (v φ ∼ v ∥ ). As within the TAPaS the poloidal magnetic flux varies along the radial direction, such an asymmetry is captured by the code; thereby, the particles with v ∥ > 0, which are passing, do not fulfill the resonant condition. Instead, the local approach of the Gene simulations does not allow to include such a radial variation of the poloidal magnetic flux, and therefore the structures of the destabilizing fastion contribution γ FD are nearly symmetric in the parallel velocity space.
As can be seen, a quite good agreement is achieved with the resonant structures within the trapping region. This demonstrates that the TAE instabilities identified in Gene simulations are excited through a wave-particle resonant mechanism induced by the trapped fast ions. It is to be noted that the three-ion scheme predominantly generated co-passing fast ions in this JET scenario [22,26]. Therefore, this consideration suggests that the TAEs are excited by the fast ions remained trapped in magnetic wells in their slowing-down process from the original co-passing motion. Yet, deeper analyses on the experimental data are needed to further validate this latter statement.

Ion-scale turbulent transport suppression
In this section, the effects of the suprathermal population on the turbulent transport are investigated by means of nonlinear simulations with the Gene code. A detailed study on the underlying complex mechanisms is reported, showing the beneficial impact of the high-frequency Alfvénic mode on the turbulence stabilization.
In order to assess the impact of the fast ions on the ionscale turbulent transport, the electrostatic heat diffusivity χ ES of the thermal species is plotted as a function of the normalized fast ion pressure gradient R/L pFD in figure 7. The fluxsurface-averaged heat diffusivities are computed with nonlinear Gene simulations. A time average is then performed over a sufficiently large time window to ensure the flux convergence in each simulation (around t ≈ 1000 a/c s ). Comparing to the simulation without fast ions (horizontal dashed lines), χ ES are only slightly or partially reduced when the fast ion pressure gradient is not large enough to excite TAEs. This is represented by the cases with R/L pFD = 5.9 and 8.1. On the other hand, when the TAEs are destabilized for R/L pFD ⩾ 10.7, the electrostatic heat diffusivities are strongly decreased, down to values around χ ES ≈ 0.5-1 m 2 s −1 . A reduction of around 95% with respect to the case without fast ions is achieved in these latter cases. Moreover, such a suppression is necessary also to match the power balances computed by the integrated modeling TRANSP, and illustrated with an horizontal black dashed line also in figure 7 for the thermal ions.
It is important to note that such a suppression is achieved in conditions of fully destabilized TAEs, as shown by the frequency spectrum of ϕ fluctuations for the cases with R/L pFD = 13.6 and R/L pFD = 16.2 in figure 8, in good agreement with the experimental observations. These results, therefore, are in line with what was previously reported [12,19], where the partial reduction was obtained in conditions of marginal stability of the fast-ion-driven modes.
It is to be stressed that in the current work mitigated electron transport is obtained also in conditions of fully destabilized AEs. Whereas for the case with R/L pFD = 16.2, the electron electromagnetic heat diffusivity is χ e EM = 3.19 m 2 s −1 , in the case with R/L pFD = 13.6 χ e EM = 0.72 m 2 s −1 , approaching the experimental power balance value of χ e = 0.5 m 2 s −1 . It is eventually reported that with a sensitivity scan on the ion and electron temperature gradient R/L Ti,e within the experimental uncertainty (∼15% of reduction) the consistency of both χ i  and χ e is obtained. Indeed, for the case with R/L pFD = 13.6 and reduced R/L Ti,e , it is obtained χ e = 0.69 m 2 s −1 (note that here the total, electrostatic plus electromagnetic, electron diffusivity χ e is reported). In these latter configurations, the ion heat diffusivities are not strongly modified.
Therefore, it is implied that to obtain a validated gyrokinetic result for pulse #94701, a precise balance of TAE amplitude and background turbulence must be fulfilled. This latter condition reveals necessary due to the strong increment of the perturbed field fluctuations in case of larger R/L pFD , thereby strongly affecting the transport levels. As demonstrated by the substantial stabilization of the turbulence-driven density fluctuations (figure 2), which is in qualitative good agreement with the gyrokinetic modeling results, these conditions of proper balance are obtained in the three-ion scenario at JET. Furthermore, it is to be noted that the Gene nonlinear simulations are performed retaining the n = 4 and n = 7 wavelengths, nonetheless excluding the n ⩽ 3 modes, as already reported in section 3. This constraint is due to the computational affordability of the numerical simulations, since including the n ⩽ 3 wavenumbers would have required to more than double the grid resolution in the binormal direction, affecting also the interplay with large-scale radial wavenumbers; therefore, extreme computational resources would have been necessary, which were not available for this study. Although this approximation is quite remarkable, the agreement between the numerical and the experimental outcomes. This may suggest that the n = 1 and n = 2 wavenumbers, clearly stable in the experiment (see figure 1), are not playing a crucial role in the mechanism; whereas the n = 3, expected to be unstable, may not impact the qualitative picture here described.

Strong ZF activity triggered by TAEs
In the following, the main causes of the ion-scale turbulent transport suppression are analyzed. Reminiscent of the mechanism proposed in [19], the large reduction of the electrostatic heat fluxes is related to an increase of the ZF activity, which is well-known to have a beneficial effect on the ITG-driven turbulent transport [20,21]. This is demonstrated by inspecting the spectrum of the ZF shearing rate γ E×B,zonal ≡ |k 2 x ϕ(k x , k y = 0)| in the radial direction, illustrated in figure 9 for different Gene simulations. It can be observed that a strong increase of γ E×B,zonal is measured in the case of fully destabilized TAEs, i.e. for R/L pFD = 13.6 and 16.2. Such an increase with respect to the case without fast ions and with stable TAEs is especially noteworthy in the low-k x region, where the ITG-induced transport is mainly concentrated. Therefore, the increase of the ZF activity can be related to the suppression of the ion-scale turbulent transport observed in conditions of fully destabilized TAEs.
Hence, the present work indicates that a multi-scale synergy between the ZFs and the high-frequency TAEs is underlying. This is consistent with previous numerical results [19,61,62] and theoretical studies [63], which demonstrated a ZF generation due to an energy transfer from the TAE spatio-temporal scale. In order to enforce such a consideration, a multimode analysis has been performed. The employed bispectral method based on the wavelet decomposition [64][65][66][67] allows the detection of the nonlinear interaction between the considered wavevectors. This highlights the multi-scale coupling between the TAEs and the ZFs. The wavelet transform of the perturbed electrostatic potential ϕ reads W ϕ kx,ky (ω, t) =ˆR ϕ kx,ky (τ )ψ ω (t − τ )dτ (5) where the selected Morlet wavelet function ψ ω is with ω the frequency of the characteristic wavelet scale. Thus, the bispectral analysis identifies the intensity of the nonlinear coupling by measuring the modulus of the bispectrum, which is defined as with the summation to be performed on both binormal and radial wavenumbers, † indicating the complex conjugate and the brackets ⟨·⟩ z representing the average over the parallel direction. Given the conditions ω = ω ′ + ω ′ ′ , k x = k ′ x + k ′ ′ x and k y = k ′ y + k ′ ′ y to be fulfilled for a nonlinear interaction, the considered triplets of wavevectors are limited in this analysis to the most destabilized TAE mode (k ′ x , k ′ y ) = (0, 0.05), the ZF components (k ′ ′ x , k ′ ′ y ) = (a, 0) and the missing mode to complete the triplets (k x , k y ) = (−a, 0.05), with a ∈ [−k x,max , +k x,max ]. For the sake of clarity, in the following the expression of the specific bispectrum calculated in this study is reported: This calculation is done for the case with R/L pFD = 16.2, in which TAEs are fully destabilized, in the TAE growing phase of the simulation. In figure 10(a), the modulus of the bispectrum displays a dominant structure at the intersection between the TAE frequency range, represented by the white dashed lines, and the ZF frequency, i.e. f ′ ′ = 0. It is to be noted that the modulus of the bispectrum is here reported as a function of the frequency in standard units instead of the normalized ones (f ≡ ω/(2π)). Thus, the analysis clearly shows that a threewave nonlinear interaction occurs between the TAE and the ZF spatio-temporal scales. Considering the strong increase of the ZF activity in correspondence of the TAE destabilization, this analysis suggests that a net energy transfer from the TAE to the zonal wavelengths is underlying, as already demonstrated by analytical studies [63,68] and shown to be occurring in numerical simulations [19].

Additional nonlinear interactions destabilizing low-frequency modes
Furthermore, as shown in theφ spectrum of figure 8, nondominant spikes are present in the range f ≈ 20-40 kHz, representing modes destabilized at k y ρ s = 0.025 and 0.05. Since an accurate identification is still missing, for the sake of simplicity, the modes are called low-frequency modes (LFMs) in the following. Thus, the same multi-mode analysis has been employed to study the dynamics of such nonlinearly excited instabilities in the low-frequency range. The bispectral analysis is performed for the wavevectors (k ′ x , k ′ y ) = (0, 0.05), the LFMs (k ′ ′ x , k ′ ′ y ) = (a, 0.025) and (k x , k y ) = (−a, 0.025), and hence the bispectral calculation in this case is performed as: As done in the previous section and depicted in panel (a) of figure 10, also here the bispectrum is calculated for the case with R/L pFD = 16.2. In panel (b), the modulus of the bispectrum for the case with R/L pFD = 16.2 illustrates an effective nonlinear interaction between the TAE and the LFM ranges of frequency. The same analysis carried out for R/L pFD = 13.6 shows very similar results, not reported here for the sake of simplicity. These bispectral analyses have been performed in a time window different from the one displayed in panel (a). Indeed, while panel (a) refers to the TAE linear growing phase, in panel (b) the analysis is related to the TAE saturated phase. For the sake of clarity, if the same study on the LFM dynamics is performed in the TAE growing phase, or for lower fast ion pressure gradient case, the nonlinear coupling between LFMs and TAEs results absent or very weak [69]. These considerations suggest that the LFMs, are destabilized through an energy redistribution from the TAE to the LFM components when the TAE amplitude is larger than a critical threshold.
Interestingly, earlier experimental studies reported a very similar three-wave nonlinear interaction between TAEs and fishbone-like instabilities in the NSTX device [70]. This may indicate that the LFMs are destabilized when largeamplitude TAE-induced fluctuations are present. These modes may also be a fingerprints of the beneficial complex mechanism triggered by the highly energetic fast ions, since similar outcomes have been obtained in preliminary gyrokinetic studies on a specific TCV pulse which shows a turbulence reduction in the presence of highly energetic ions [69].

Role of cross phase in the transport suppression
The impact of the TAE destabilization on the cross-phase of the heat-flux-relevant quantities is now analyzed. Indeed, the transitioning of different turbulence regimes can be detected by inspecting the cross phase between fluctuating physical parameters, such as the electrostatic potential and the electron/ion density or pressure [71]. This technique is widely used in experiments (see, e.g. [72]), and it could be also employed as an effective validating feature in gyrokinetic studies [73,74]. In the present work instead, it is exploited to gain more insights on the turbulent transport suppression mechanism found to be active in the gyrokinetic simulations. As no cross-phase measurements are available in the three-ion scenario at JET, a comparison with experiments of such results is not possible to date.
Since the cross phase and the amplitude of the fluctuating quantities determines together the transport level, the phase angle α between the quantities A and B and defined as can indicate a possible underlying mechanism in the turbulent transport suppression. The phase angle distribution is illustrated in figure 11 for three different cross phase, i.e. ϕ × n, ϕ × T ∥ and ϕ × T ⊥ (for the 3 He), and for two different cases, namely the simulation without fast ions (top row) and with R/L pFD = 16.2 (bottom row). The cross phase is measured for each radial x and parallel z grid-point, and then collected into 62 bins spacing from −π to +π. Eventually this is shown in histograms as a function of the binormal wavenumbers, after being averaged over a significant time window of the nonlinear Gene simulation and weighted by the amplitude of the product of the respective quantities absolute values at each k y ρ s . The three combinations displayed in figure 11 are the kernel of the computation of the electrostatic thermal ion heat flux. Note that for the sake of simplicity only the quantities for the 3 He are reported, but a very similar behavior is observed also for the thermal deuterium and for the electron quantities.
In the case without retaining the fast ion component, the turbulence is dominated by the ITG-induced fluctuations. This is confirmed by monitoring the cross phase, as illustrated by the top row of figure 11. The dominant peaks of the distributions are localized at the binormal wavenumber k y ρ s = 0.2, in the ITG region, and the phase α is close to π/2 [75]. In these conditions of phase angle distribution, the heat transport is favorable [76]. Nevertheless, when the TAE regime is dominant, likewise in the case with R/L pFD = 16.2, a noteworthy shift in the cross phase is observed. The dominant peaks in the contour plots of the bottom row of figure 11 are indeed localized in the TAE spatial scales (k y ρ s = 0.025) and at α = ±π. This clearly demonstrates a mode transition from ITG-dominated to TAE-dominated regime. Moreover, although the TAEs increase the amplitude of the field and physical parameter fluctuations, the heat flux is strongly reduced since the angle distribution is mainly in-phase. Indeed, when the relative phase between heat-flux-relevant quantities is a multiple of π, the transport cannot be developed. It should be stressed that, not only is the phase angle of the dominant peaks modified, but also the spatial scales at which the outward electrostatic heat fluxes are mainly generated are shifted. A similar phase angle shift has already been observed in earlier local gyrokinetic studies for JET plasmas [77], moving from an ITG-dominated regime to a KBM/BAE one. Eventually, this result of cross-phase shift from out-of-phase to in-phase conditions is consistent with the turbulent transport suppression measured with unstable TAEs. The preferentially cross-phase shift toward the values α = ±π, instead of α = 0 which would be equivalent in term of contribution to the resultant heat flux, remains yet elusive. Future efforts will be dedicated to study this cross-phase shift in more details.
Therefore, the described picture may indicate a crucial role of the cross-phase angle distribution in the complex mechanism leading to the turbulent transport suppression. Indeed, a remarkable correlation between the heat flux reduction and the modification of the dominant structures in the cross phase distribution has been shown. As the ZF activity is hugely increased in conditions of fully destabilized TAEs, the cross-phase shift may be a subsequent effect of the strong sheared flows generated by such high-frequency TAEs [78,79]. Nevertheless, a definitive answer to this issue is still elusive and further works are required.

Effects of the reversed shear on the fast-ion transport
In this section, the fast-ion transport induced by the largeamplitude TAE fluctuations of the fields is assessed, and possible actions to be considered for its control are suggested. Indeed, the total fast ion heat diffusivity in the case with fully destabilized TAEs (R/L pFD = 16.2) is + χ FD = 76.1 m 2 s −1 , with the superscript + indicating the simulation with positive magnetic shear. This large value is caused by the perturbations of the electrostatic and magnetic potentials induced by the TAE destabilization. Even in the case of consistency with experimentally-derived power balances, i.e. with R/L pFD = 13.6 and decreased R/L Ti,e by 15%, the fast-ion transport is far from being negligible ( + χ FD = 18.9 m 2 s −1 ). Nevertheless, the experimental outcomes imply that a good confinement for the fast ions generated through the three-ion scheme is achieved [22,26,80]. This consideration suggests that an additional mechanism improving the fast-ion confinement is present. A suitable candidate is represented by the modification of the magnetic equilibrium. Indeed, a systematic observation of chirping-up RSAEs destabilized in the threeion scenario [22,26,36] demonstrates that the safety factor presents a non-monotonic radial profile. The physical mechanism leading to the reversal of the q-profile is, however, still not identified, although the substantial population of MeV-ions generated as co-passing particle in the plasma core may drive non-inductive current inducing, thus, the reversed magnetic shear [22,26,36]. An important concurrent observation with the RSAE destabilization is about the neutron rate and plasma stored energy signals. In fact, an increase of both these latter physical parameters is measured in the same time window where RSAEs are found unstable. This consideration suggests a correlation between the two phenomena, and also motivates the following gyrokinetic studies performed in negative magnetic shear geometry.
The simulations presented in the previous sections have been performed with a positive magnetic shear, as also displayed in table 1. This is consistent with the EFIT equilibrium computed for pulse #94701 at t ≈ 9.6 s. Moreover, as shown in the Mirnov coil spectrogram of figure 1, only right after t ≈ 9.6 s the chirping-up RSAE harmonics are clearly visible.
As already pointed out, the detection of such instabilities distinctly indicates that a non-monotonic safety factor profile is obtained even in JET pulse #94701. The studies on the radial localization of the RSAEs place such Alfvén eigenmodes in the deep core [36]. Nevertheless, large uncertainties affect the calculation of the safety factor profile in the inner core. Thus, to possibly investigate and determine the shape of the q-profile in the core, detailed sensitivity scans with numerical simulations must be performed.
Therefore, in order to study the impact of a reversed q-profile, both linear and nonlinear gyrokinetic simulations have been carried out in different configurations. The only difference in the numerical setup reported in table 1 is on the local value of the magnetic shear. This is justified by the will of disentangling other possible physical effects, and to focus only on the impact of the reversed shear on the turbulence regime and transport characteristics. In figure 12, the linear growth rate and the mode frequency are plotted for three different binormal wavenumbers as a function of the magnetic shear valueŝ for the case with fast ion pressure gradient R/L pFD = 16.2. The selected k y represent the two TAE modes retained in the nonlinear simulations (k y ρ s = 0.025 and k y ρ s = 0.05) and the ITG peaking mode (k y ρ s = 0.4). The ITG growth rate is modified by the magnetic shear, presenting lower values in the negative magnetic shear region [81]. It is important to note that the ITG is the dominant instability for the selected space parameter −0.5 <ŝ < 0.85, as shown by the absence of sharp transitions in the mode frequency spectrum. For what concerns the wavenumbers k y ρ s = 0.025 and k y ρ s = 0.05, instead, the modification of the magnetic shear affects more strongly the TAE linear stability. As can be seen, the TAE is not the dominant instability in the range 0.1 ≲ŝ ≲ 0.4. Importantly, the negative magnetic shear already delineates a reduction of the linear growth rate, more evident for the wavelength k y ρ s = 0.05. Therefore, a lower turbulent transport is expected to develop in the nonlinear phase, however maintaining the TAE definitely unstable.
The negative magnetic shear has a stabilizing effect on the ballooning-type instabilities, by twisting the magnetic field lines along the particle motion and thus by reducing the curvature drive of the modes [82][83][84]. In the negative magnetic shear configuration, the ballooned perturbations are stabilized, and this explains the decrease of both ITG and TAE linear growth rate shown in figure 12. An additional mechanism consistent with the decrease of the TAE drive in this particular case is likely to be represented by the stabilizing effect of the negative magnetic shear on the trappedparticle-driven instabilities [82]. Indeed, as shown in figure 6, the main contribution to the TAE excitation is given by the trapped fast ions. The stabilization of trapped-particledriven instabilities has also been validated in experimental devices [85].
Therefore, an additional analysis in the linear regime has been carried out to enforce these latter considerations. The free energy exchanged in the velocity space has been analyzed for a numerical setup withŝ = −0.3, following the method already described in section 3.1. In figure 13, the fast-ion contribution to the instability for the binormal   k y ρ s = 0.05, the resonance is more localized towards v ∥ = 0, whereas for k y ρ s = 0.025 the structure has even disappeared, leaving the passing fast ions to mainly contribute to the TAE destabilization.
The nonlinear simulations then illustrate a coherent picture, where the negative magnetic shear configuration drastically reduces the fast-ion transport. This is shown in figure 14, where the time evolution of the fast ion heat flux for both the configuration withŝ = 0.63 (solid curves) andŝ = −0.3 (dotted curves) are compared. After the modification of the magnetic equilibrium, occurring around t = 200 a/c s , the fluxes settle at considerably lower levels. The horizontal black lines illustrate the time-averaged values, starting from which is possible to determine the heat diffusivities. Hence, the heat diffusivities move from + χ FD ES = 71.9 m 2 s −1 and + χ FD EM = 4.19 m 2 s −1 in the case withŝ = 0.63 to − χ FD ES = 1.3 m 2 s −1 and − χ FD EM = 0.4 m 2 s −1 , respectively for the electrostatic and electromagnetic contributions. Such a decrease implies that a better fast-ion confinement is achieved in conditions of negative magnetic shear, in which the TAEs are however found unstable. Indeed, since the TAEs are still dominating the frequency spectrum, the turbulent transport of the thermal ions is suppressed. Very slight variations of the thermal ion diffusivities are indeed measured moving from positive to negative shear configuration.
Hence, the decrease of the fast-ion transport in conditions of reversed magnetic shear measured in these gyrokinetic local simulations indicate an improved confinement of the suprathermal component of the plasma. This may unveil a relation with the increased neutron rate and plasma energy content experimentally measured in conditions of destabilized RSAEs [22,26]. However, although promising, further studies with more specific numerical tools are required to better assess the consistency of these results with the experimental outcomes.

Conclusions
The improvement of bulk-ion thermal confinement was observed in recent JET fast-ion experiments in D-3 He plasmas with strong core electron heating from MeV-range D ions. The detailed transport analysis of these studies reveals a novel mechanism of turbulence suppression in the presence of highly energetic ions and unstable toroidal AEs.
In first place, linear gyrokinetic simulations, performed in the flux-tube approximation, successfully identify the unstable TAEs in good agreement with the experimental evidence. In addition, in-depth multi-code analyses describe the TAE excitation mechanism by trapped fast ions in their slowing-down process. Subsequently, the nonlinear simulations demonstrate a clear suppression of the electrostatic turbulent transport developed at the ion-scale only in conditions of unstable TAEs. A bispectral analysis establishes that a strong nonlinear coupling between the TAE and the ZF dynamics exists, suggesting that the strong stabilizing zonal activity is due to a net energy transfer from the TAE scales consistently with earlier studies [19,61,63]. An additional underlying mechanism involving the shift of the cross phase of the heat-flux-relevant quantities from out-of-phase to in-phase is also demonstrated to have a crucial role in the turbulent transport suppression.
Numerical analyses about the impact of the non-monotonic q-profile on the fast-ion confinement are also performed. This reveals necessary due to the difficulty of measuring the qprofile in the inner core. A beneficial effect on the fast-ion heat flux is measured in the configuration with negative magnetic shear. This finding is in line with the experimental observations of increased neutron rate and total plasma energy content in conditions of reversed shear in the deep plasma core [22,26]. It is to be noted that additional scans on the local values of q andŝ may be required to assess the shape of the safety factor profile in the plasma core.
It is also to be noted that the validity of the flux-tube approximation used in the gyrokinetic numerical analyses has been carefully checked, also with experimental validating studies. Nevertheless, global simulations and detailed analyses on this specific topic are planned to be reported in future contributions, as non-local effects may also be important for a quantitative comparison with experimental diagnostic measurements.
These promising results indicate that similar conditions might be realized in ITER and future fusion power plants, where fusion-born alpha particles provide a strong source of core electron heating [86] and can potentially destabilize TAEs [6]. The latter effect is usually considered as potentially detrimental for plasma confinement. Nonetheless, the results of the present studies show that this is not necessarily so. TAEs can contribute to efficient ion heating through the onset of intense ZFs, providing an improved thermal insulation of the central plasma. Earlier predictive gyrokinetic studies already pointed out a beneficial impact of the fusion-born alpha particles on the ITG-driven turbulent transport in ITER hybrid plasmas [2], but the AE activity was purposely avoided. Indeed, the current work shows that AE activity is not necessarily a phenomenon to be avoided in D-T plasmas, but actually one that could be purposely tailored using external actuators [87]. The identification of this mechanism for turbulence suppression has important implications for fusion research, as it holds promise to enhance the performance of future fusion devices with strong alpha particle heating.
It is finally noted that the concept of anomalous ion heating by alpha particles was introduced to explain the observed increase in T i during the past full-scale D-T experiments on JET [88]. Consistent with previous linear analyses [89,90], the novel mechanism identified in this study points out that this anomalous ion heating could well be identified with the suppression of microturbulence when AEs are marginally stable or even unstable, resulting in a very effective heating of the thermal ions in plasmas with a significant population of MeVrange fast ions.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.