Impact of supra-thermal particles on plasma performance at ASDEX Upgrade with GENE-Tango simulations

To support the quest for sustainable and efficient fusion energy production, it is crucial to develop sophisticated tools that can accurately capture the main physical effects occurring in plasma experiments. These tools will enable reliable predictions of plasma profiles, facilitating the investigation of various plasma regimes and assessing optimal conditions for maximizing fusion performance while avoiding unexpected confinement degradation.

In recent years, there has been significant progress in developing reduced turbulence models that efficiently capture the essential features of plasma micro-instabilities and turbulence. Notable examples of such models include TGLF (trapped-gyro-Landau-fluid) [1, 2] and QuaLiKiz [3–6]. These models offer the advantage of fast computation times, enabling turbulent flux calculations—based on saturation rules—within seconds on a single CPU. They have been successfully integrated with transport codes like TGLF-ASTRA [7, 8], TGLF-TGYRO [9], and QuaLiKiz-JINTRAC [10]. However, despite the ability of these reduced turbulence models to accurately reproduce experimental profiles across a wide range of plasma conditions, they still lack the capability to capture nonlinear energetic particle effects on core plasma turbulence [11–14]. They should be able to capture some linear fast ion effects [15]. As a result, they tend to underestimate the on-axis ion temperature in regimes where energetic particles play a significant role. This limitation might be particularly severe in future burning plasma experiments, impacting both the initial ramp-up phase, which relies on auxiliary systems to reach the optimal conditions for fusion reactions, as well as the subsequent operational phase where heating is solely reliant on MeV alpha particles.

Nonlinear gyrokinetic codes have proven effective in capturing energetic particle effects on turbulence in various plasma regimes [16–20], offering valuable insights into the underlying physical mechanisms that lead to observed ion-scale turbulence suppression in experiments. However, these codes are computationally significantly more demanding than quasi-linear models, despite extensive work to improve code performances, e.g. via GPU porting [21]. As a result, gyrokinetic simulations are often limited to single flux-tubes, and even in the few global simulations conducted so far, plasma profiles are not allowed to freely evolve due to the combined effect of turbulence and external sources, but are kept fixed by using artificial Krook-like operators. Recent advancements have successfully addressed these limitations [22–25]. A notable example is the coupling of the global gyrokinetic code GENE with the transport solver Tango [23, 25]. This coupling takes advantage of the distinct time-scale separation between microscopic and macroscopic physics, resulting in a substantial reduction in computational costs allowing, —for the first time—high-fidelity simulations that extend up to the transport time scale. In this coupling, the turbulence code determines the turbulence fluxes, while the transport code extracts transport coefficients based on the plasma profiles and evolves the plasma profiles themselves considering sources and turbulent fluxes. This approach ensures that the heavily computational gyrokinetic code is no longer constrained by the micro-turbulence time scale, which is several orders of magnitude shorter than the macroscopic scale at which background profiles evolve. Currently, there is an ongoing effort [26] to investigate potential differences between stand-alone gyrokinetic simulations up to the transport time-scales and gyrokinetic codes coupled to transport codes to speed-up the calculations, such as GENE-Tango. However, due to the immense computational cost of stand-alone gyrokinetic codes up to the transport time scale, this work has been confined to scenarios involving adiabatic electrons and simplified geometries, neglecting any energetic particle species. Future endeavors aim to broaden this investigation to include kinetic electrons.

In our study, we present the first GENE-Tango simulations of an ASDEX Upgrade H-mode discharge with large fast ion content, demonstrating a substantial improvement in plasma performance attributed to energetic particles. Remarkably, we find that only when the energetic particle species is consistently included in the radially global nonlinear electromagnetic GENE-Tango simulations, we observe an excellent agreement with the experimental plasma profiles. Removing the fast ions from the model leads to a significant flattening of the ion temperature profile, aligning with profiles computed using TGLF-ASTRA, which significantly under-predicts the on-axis ion temperature. Detailed analysis reveals that the thermal ion temperature gradient (ITG) increases markedly in correspondence with Kinetic Ballooning Modes (KBMs)/Alfvén eigenmodes (AEs). These modes enhance zonal flow activity and significantly enhance shearing rate levels compared to GENE-Tango simulations without energetic particles, resulting in turbulence suppression consistent with previous flux-tube simulations. When these modes become more unstable during the GENE-Tango simulation, turbulent fluxes increase for all species and Tango reacts relaxing thermal profiles.

Additionally, we observe that while the global GENE simulations can tolerate marginal instability of these modes and still achieve power balance agreement, flux-tube simulations impose stricter constraints on the linear stability of these modes. As a consequence, flux-tube simulations are likely to undergo a more pronounced flattening of plasma pressure profiles compared to global simulations to reduce the drive of these modes until they reach a state of marginal stability. These findings suggest that a global gyrokinetic treatment might be necessary to accurately model the effects of energetic particles on turbulence in regimes where energetic particles and electromagnetic effects are expected to play a non-negligible role on plasma turbulence.

This paper is organized as follows. The ASDEX Upgrade discharge is described in section 2. While the GENE-Tango coupling is detailed in section 3, we outline the numerical setup employed for the GENE-Tango simulations in section 4. In section 5, we present the steady-state profiles obtained from TGLF-ASTRA, highlighting a significant under-prediction, particularly of the on-axis ion temperature. Building upon this, section 6 shows the plasma profiles resulting from various GENE-Tango simulations that retain or neglect energetic particles and electromagnetic effects in GENE. Notably, we demonstrate that the only case that successfully matches the experimental plasma profiles is when both fast particles and electromagnetic effects are retained in the GENE-Tango calculations. To shed light on the physical mechanisms driving the observed on-axis peaking of the temperature profile in the GENE-Tango simulations with fast ions and electromagnetic effects, we provide a detailed analysis in section 7. Additionally, stability analyses on the final steady-state GENE-Tango profiles that successfully reproduce the experimental results are conducted in section 8. The impact of strongly unstable high frequency electromagnetic modes on the thermal profiles is discussed in section 9. Section 10 presents a comparison between flux-tube and global GENE simulations for the steady-state GENE-Tango profiles. Finally, conclusions are drawn in section 11.

The plasma discharge under investigation in this study corresponds to the ASDEX Upgrade pulse $\#39230$. This is a deuterium plasma with an on-axis magnetic field strength of approximately $B_0 = 2.5$ T. The plasma current is $I_p \approx 0.8$ MA, and the safety factor at the flux surface containing 95$\%$ of the poloidal flux is $q_{95} \approx 5$. The heating power is comprised of $P_\mathrm{NBI} \approx 7.5$ MW of Neutral Beam Injection (NBI) with perpendicular beams and $P_\mathrm{ECRH} \approx 1.6$ MW of Electron Cyclotron Resonance Heating (ECRH). Figure 1 illustrates the time evolution of various parameters, including the heating powers, on-axis plasma temperatures and density. Notably, the application of NBI is closely associated with a significant rise in ion temperature. This effect is most prominent in the time range of $t = [2.5\text{-}3.2]$ s, i.e. at the maximum NBI power in the discharge, where the ion temperature reaches its peak values at the center, with relatively smaller impacts on electron temperature and density. Consequently, our analysis focuses on a specific time slice at t = 2.7 s, during which the maximum external NBI power is employed, resulting in a pronounced effect on the on-axis ion temperature. The electron temperature and density measurements are obtained using IDA [27], while the ion temperature and plasma toroidal rotation are determined through IDI [28].

Zoom In Zoom Out Reset image size

Considering the substantial NBI external heating power applied in this plasma discharge, there exists a non-negligible fraction of supra-thermal ions. In figure 2, we show the temperature and density profiles of these energetic particle species, computed via RABBIT [29] using the experimental profiles.

Zoom In Zoom Out Reset image size

The numerical simulations presented in this paper were performed using the newly developed tool GENE-Tango, which combines the gyrokinetic code GENE [30, 31] with the transport solver Tango [32]. This coupling enables reliable profile predictions with significant speed-ups compared to standalone flux-driven gyrokinetic simulations on the transport timescale [25]. The underlying scheme follows the approach employed in reduced turbulence models and transport codes, where turbulence and transport phenomena are simulated only on their natural timescales [22].

The coupling starts with a global GENE simulation using an initial guess for the temperature and density profiles, running until saturation. Afterwards, the turbulent fluxes (averaged over the saturated phase) are extracted and transferred to Tango. Tango, in turn, evaluates new plasma (temperature and density for each species) profiles consistent with the given turbulence levels and experimental sources, evolving temperatures and densities for each species. These updated profiles are then fed back into GENE, restarting from the previous turbulent phase and running until saturation. Given the restart from the previous iteration, the system quickly converges to a new saturated phase within a short time domain (approximately $t = 150 c_s/a$ for the simulations performed in this paper, where c_s represents the sound speed defined in section 4 and a the minor radius). This iterative process continues until the turbulent fluxes match the volume-averaged energy and particle injected by external sources. Achieving a steady-state solution is followed by further verification by performing a stand-alone GENE simulation with the final plasma profiles computed by Tango. The coupled GENE-Tango simulation is considered fully converged only when this stand-alone GENE simulation yields the same fluxes as the last GENE-Tango iteration. The coupling scheme is illustrated in figure 3.

Zoom In Zoom Out Reset image size

The use of GENE-Tango significantly reduces the computational expense by several orders of magnitude, making first-principle core plasma simulations over the confinement time feasible, even for devices such as ITER, using currently available computing resources [25]. The validity of GENE-Tango has been demonstrated through validation studies on multiple ASDEX Upgrade discharges with moderate external heating powers, where the effects of energetic particles on plasma profiles were negligible [25]. Recently, the application of GENE-Tango has been extended to optimized stellarator devices [33]. However, the results presented in this paper are based only on radially global GENE-Tango analyses. A radially global treatment is indeed required to describe energetic particle dynamics and energetic-particle-driven modes correctly [34].

In the numerical simulations conducted for this study, the NBI heating power deposited on electrons and ions, as well as the fast deuterium pressure, were calculated using RABBIT [29]. The ECRH power was computed using TORBEAM [35, 36], while the ohmic power density was determined from ASTRA [7] and remained fixed throughout the GENE-Tango loop. The power density of the collisional energy exchange between ions and electrons was self-consistently computed by Tango and varies at each iteration. Additionally, the radiated power density, which takes into account bremsstrahlung radiation, was calculated using the experimental effective charge profile and the tungsten density. While the effective charge profile and the tungsten density were kept constant during the GENE-Tango loop, the radiated power density is consistently adjusted at each iterations accounting for variations in the electron temperature profile.

The GENE simulations are performed employing the realistic ion-to-electron mass ratio. External flow shear and collisions are retained and the latter modeled using a linearized Landau-Boltzmann collision operator with energy and momentum conserving terms [37]. Neoclassical fluxes are expected to be negligible for the analyzed discharge and are therefore not included in the calculations [38, 39]. For simplicity, each species has been modeled with a Maxwellian distribution function. Specifically, for the energetic particle species, this involves a Maxwellian distribution with the zeroth and second-order moments equal to the fast ion density and energy as computed by RABBIT. In the near future, the simplification of an equivalent Maxwellian distribution will be relaxed. While the GENE-Tango loop allows us to evolve the energetic particle profiles through Rabbit, we have kept the energetic particle profiles fixed at the ones computed for the experimental plasma profiles for each GENE-Tango simulations. This approach facilitates the comparison between the various simulations performed that cover electrostatic and electromagnetic regimes considering energetic particles.

The use of a more realistic distribution functions for modeling the energetic particle species is not expected to yield significant differences in the numerical results [40]. While GENE-Tango has recently been coupled with CHEASE [41] to enable the self-consistent evolution of magnetic geometry as the plasma pressure profiles evolve within the Tango iterations, for simplicity, we have kept the magnetic equilibrium fixed to the one reconstructed via IDE [42] for the ASDEX Upgrade discharge $\#39230$ at the time t = 2.7 s.

In this paper, we perform GENE-Tango simulations using only the radially global version of GENE. These simulations employ a grid resolution of $(n_x \times n_{k_y} \times n_z) = (225 \times 36 \times 32)$ in the radial (x), bi-normal (y), and field-aligned (z) directions. The velocity space grids consist of $(n_{v_\shortparallel} \times n_\mu) = (48 \times 32)$ points, with $v_\shortparallel/v_{\mathrm{th},s} = [-3.5, 3.5]$ and $\mu B_0 / T_s = [0,12]$ for each species. Here, $v_\shortparallel$ denotes the velocity component parallel to the background magnetic field, and µ is the magnetic moment. These extended velocity ranges can capture for changes in pressure profiles within the radial box. The radial domain covered in the radially global GENE-Tango simulations is $\rho_\mathrm{tor} = [0.025, 0.7]$. Here, $\rho_\mathrm{tor}$ represents the radial coordinate based on the toroidal flux Φ, calculated as $\rho_\mathrm{tor} = \sqrt{\Phi / \Phi_\mathrm{LCFS}}$, where $\Phi_\mathrm{LCFS}$ represents the value of the toroidal flux at the last closed flux surface.

In GENE, $v_\shortparallel$ is normalized to the thermal velocity of each species s, $v_{\mathrm{th},s} = \left(2T_s/m_s\right)^{1/2}$, where T_s is the temperature and m_s the mass; and the magnetic moment is normalized to B₀, the on-axis magnetic field and T_s . In the electrostatic simulations, we set the minimum toroidal mode number $n_{0,\mathrm{min}} = 2$, while in the electromagnetic simulations, we used $n_{0,\mathrm{min}} = 1$. The toroidal mode number n is discretized as $n = n_{0,\mathrm{min}} \cdot j$, where j ranges from 0 to $n_{k_y}-1$. In the electromagnetic global simulations $B_{\shortparallel}$ fluctuations are neglected.

The global GENE simulations are performed in gradient-driven mode, using Krook-type particle and heat operators to maintain the plasma profiles close to the reference profiles provided by Tango during each iteration [31]. To achieve this, the Krook heat coefficient (γ_k ) and particle coefficient (γ_p ) are set to $\gamma_k = 0.01 c_s / a$ and $\gamma_p = 0.01 c_s /a$, respectively. Here, c_s represents the sound speed, calculated as $(T_e / m_i)^{1/2}$, where T_e is the electron temperature at the reference radial position, and m_i is the bulk ion mass in proton units. In order to mitigate the possible impact of unresolved scales and instabilities such as electron temperature gradient driven modes, a numerical fourth-order hyperdiffusion method is employed to dampen fluctuations [43].

Additionally, buffer regions covering 10% of the GENE simulation's radial domain are implemented. These buffers serve the purpose of suppressing fluctuations to zero near the boundaries of the domain and ensuring compatibility with the Dirichlet boundary conditions. Within these regions, a Krook operator with an amplitude of $\gamma_b = 1.0 c_s/a$ is applied.

Tango employs different strategies to extend turbulent fluxes and sources in the GENE buffers, accounting for the undesired damping effects on plasma fluctuations. Within the inner buffer region, the physical sources are intentionally set to zero, with a minimal impact on the volume integral of the injected sources due to the reduced plasma volume near the magnetic axis. At the outer buffer, we employ extrapolation regions in Tango, with distinct radial domains chosen for temperatures and density. These regions are carefully selected in areas that remain unaffected by the Krook operator. Within these specified domains, Tango performs interpolation of the GENE turbulent fluxes. As we move beyond this region towards the outer boundary, the turbulent fluxes are replaced by a linear extrapolation until reaching the end of the GENE-Tango grid. We refer the reader to references [23, 25, 32] for further details.

To enhance numerical stability, Tango incorporates relaxation coefficients for plasma pressure $(\alpha_p)$ and turbulent fluxes $(\alpha_q)$, facilitating effective averaging over previous iterations and mitigating statistical variance arising from turbulence simulation [23, 25, 32]. In this study we used $\alpha_p = 0.1$ and $\alpha_q = 0.3$. Furthermore, Neumann boundary conditions are applied for the inner boundary, while Dirichlet boundary conditions are used for the outer boundary in Tango.

We begin our analysis of the ASDEX Upgrade discharge $\#39230$ at t = 2.7 s using TGLF-ASTRA covering the same GENE-Tango radial domain $(\rho_\mathrm{tor} = [0.025 - 0.7])$. In this analysis, we employed the SAT2 saturation rule for TGLF [44]. The steady-state profiles obtained from this analysis are shown in figure 4 and compared with the experimental profiles. While we observe a good agreement for the electron density, we note a significant underestimation of the temperature profiles, especially concerning the ion temperature predicted by TGLF-ASTRA, in comparison to the experimental measurements. This finding is consistent with previous observations indicating that TGLF-ASTRA generally fails to capture the positive effects of supra-thermal particles on turbulent transport often linked to nonlinear electromagnetic effects, resulting in larger turbulent fluxes and, thus flatter profiles [11–14]. It is worth mentioning that recently, wave-particle resonant effects between energetic particles and ion-scale micro-instabilities have been successfully reproduced with TGLF simulations [15]. As fast ions primarily impact the thermal ion turbulent fluxes, this discrepancy leads to a noticeable decrease in the on-axis T_i , with a minor effect on the electron temperature while still maintaining qualitative agreement with the density profile.

Zoom In Zoom Out Reset image size

We apply the global version of GENE interfaced with Tango to simulate this ASDEX Upgrade discharge, employing the numerical setup described in section 4. The plasma profiles were initialized to the ones reconstructed via IDA for electron temperature and density, and IDI for toroidal rotation and ion temperature. While the toroidal rotation remained fixed during the GENE-Tango iterations, the other profiles were self-consistently adjusted due to the combined impact of heating and particle sources. When included in the modeling, the temperature and density profiles of energetic particles were computed using RABBIT. Additionally, the thermal ion density is kept equal to the electron density profile to ensure quasi-neutrality in the simulations without energetic particles.

Four simulations were performed: (i) electrostatic without fast particles, (ii) electromagnetic without fast particles, (iii) electrostatic with fast particles, and (iv) electromagnetic with fast particles. In each case, GENE-Tango was run until the global turbulent fluxes computed by GENE matched the volume integral of the injected sources at each location, e.g. as shown in figure 5 for the electromagnetic simulation with fast particles. The final steady-state profiles are illustrated in figure 6 for each of the different cases and compared with the experimental profiles.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6 shows a particularly good agreement between GENE-Tango and experimental profiles for the electron temperature and plasma density. However, when focusing on ion temperature, it is found that achieving an agreement with experimental data requires the inclusion of both fast particles and electromagnetic effects in the modeling. This suggests that the joint influence of these factors is crucial in accurately reproducing the observed T_i profiles. It further highlights the importance of considering the complex interplay between fast particles and electromagnetic effects for a realistic description of plasma behavior in plasma regimes with high external heating power and energetic particle concentration.

The analysis also reveals that the magnetic equilibrium alone, characterized by features such as reversed shear and rational surfaces, is insufficient in explaining the T_i peaking, since it was kept fixed in all of the different simulations (figure 7).

Zoom In Zoom Out Reset image size

Electrostatic simulations, regardless of the inclusion of fast ions, predict excessive turbulent transport. As a consequence, the central T_i remains at a relatively constant value of $T_{i,0} = 4$ keV, with minimal effects on T_e and plasma density. On the other hand, the inclusion of electromagnetic effects without fast ions leads to a mild peaking of T_i . Although there is some improvement in reproducing the experimental profiles, the influence on the on-axis T_i value remains limited.

The central role played by electromagnetic effects and fast particles on the ion temperature profile peaking becomes evident when looking at the radial profile of the logarithmic gradient of the thermal ion temperature and the energetic particle heat flux. This is shown in figure 7, together with the results obtained from the electromagnetic GENE-Tango simulations without fast particles (as a reference).

Figure 7 illustrates a clear correlation between the regions experiencing a significant rise in the logarithmic gradient of the thermal ion temperature and the locations where the energetic particle heat flux increases. This occurs at the rational surface q = 1 and at s = 0. These locations correspond to the excitation of high frequency KBMs/AEs (as discussed in section 8) [45]. This indicates that the presence of these modes is likely to be linked to the turbulence suppression of the thermal ion turbulent fluxes, resulting in the peaking of the thermal ion temperature profiles.

The presence of these modes is particularly important in the nonlinear GENE simulations to strongly suppress ion-scale turbulent transport and increase the logarithmic gradient of the thermal ion temperature profiles. This effect is illustrated in figure 8 where the time evolution of the radial profile of the flux-surface averaged radial electric field is shown for the electromagnetic simulation with supra-thermal particles. Interestingly, we observe the generation of zonal layers in the proximity of the radial location where the energetic particle heat fluxes are larger in the nonlinear GENE simulations. However, no such behavior is observed in the zonal currents (not shown here), as they do not exhibit a noticeable increase at the rational surface q = 1 or at the location s = 0. Additionally, figure 8 compares the time-averaged (over the saturated nonlinear phase) radial electric field computed in electromagnetic simulations with and without energetic particles during the saturated nonlinear phase. It is evident that the simulation without fast particles lacks the layer at s = 0 of the radial electric field, which is driven by the KBM/AE that is absent in the simulation without fast particles.

Zoom In Zoom Out Reset image size

Additionally, it is worth noting that while the electromagnetic simulation neglecting fast particles does indicate an increase in the radial electric field around q = 1, this effect is confined to a localized region. In contrast, in the simulation incorporating fast particles, the radial electric field exhibits a similar magnitude at q = 1, but with the notable difference of a more extensive zonal structure that spans a broader radial range from $\rho_\mathrm{tor}$ = [0.2–0.3]. This extended zonal structure may account for the larger amplification of the logarithmic temperature gradient in the thermal ions near the rational surface q = 1 when fast ions are present, as compared to the case without fast ions.

These results are consistent with previous findings observed in flux-tube simulations, showing that fast particles play a significant role in destabilizing marginally stable modes through nonlinear interactions [18, 46]. In particular, as proposed by [47], the mechanism leading to the nonlinear excitation is likely the inverse Compton scattering by energetic ions. These modes are shown to lead to an energy redistribution from the thermal modes to the fast ion-driven modes, initially depleting the energy content of ion-scale turbulence. When the drive from fast ions is sufficiently large, it was observed that these fast particle modes interact with zonal flows, leading to a substantial increase in shearing rate levels, ultimately resulting in a strong suppression of heat and particle fluxes. The impact of energetic particle modes in enhancing zonal flow and zonal currents has also been observed in global gyrokinetic [48, 49] simulations and predicted by analytic theory [47, 50, 51].

Finally, figure 7 also shows that the other rational surfaces do not significantly impact the fast ion heat flux. This can be attributed to the relatively low fast ion density and temperature at $\rho_\mathrm{tor} \gt 0.3$, which strongly reduce the drive of these modes. We observe that the simulations without fast ions—despite having same geometry—do not show clear improvements of the logarithmic gradient of the thermal ion temperature profile at s = 0 and q = 1.

Given the high relevance of the electromagnetic modes destabilized by fast particles on the T_i temperature peaking observed in GENE-Tango, it is crucial to characterize their linear and nonlinear stability. To this aim we analyze the electrostatic and electromagnetic ($A_\shortparallel$) potentials obtained from the nonlinear stand-alone GENE simulations retaining fast particles and electromagnetic fluctuations. The results, shown in figure 9, are obtained using the GENE-Tango steady-state plasma profiles. Figure 9 reveals that the mode with toroidal mode numbers n = 7 at $\rho_\mathrm{tor} = 0.13$—which corresponds to the location of s = 0—exhibits the largest amplitude. Additionally, significant contributions are also observed from modes with $n = [1,2]$ at $\rho_\mathrm{tor} = 0.23$, located at the rational surface q = 1. These modes have a considerable impact on the turbulent heat fluxes at these specific locations for each plasma species. This is shown in figure 10, which displays the ion heat flux spectra (electron heat flux is negligible at s = 0) at the location of s = 0 (left) and the electron spectra at q = 1 (right). In particular, the turbulent fluxes at s = 0 are dominated by the mode at n = 7 and its harmonics. Notably, figure 10 shows a particularly large electromagnetic contribution (flatter) to the turbulent fluxes at these specific locations of interest for the modes n = 1 and n = 7, thereby highlighting the fundamental electromagnetic nature of these modes. In order to provide a more comprehensive characterization of these modes, we present in figure 11 the frequency spectra of the electrostatic potential for different toroidal mode numbers obtained from the nonlinear electromagnetic radially global GENE simulation using the steady-state GENE-Tango profiles taken at the location of s = 0 (left) and q = 1 (right). Figure 11 reveals that at the location s = 0, the electrostatic potential spectrum is predominantly influenced by high-frequency modes at $n = [7, 14]$, with only a minor contribution from the low-frequency ITG branch. As we move to the rational surface at q = 1, the significance of the low-frequency branch increases. However, it is important to note that a non-negligible contribution from a high-frequency mode with ω ≈ 150 kHz, occurring within the range of toroidal mode numbers $n = [1\text{-}7]$, is still observed. Figure 11 also showcases slices of the electrostatic potential spectra for the toroidal mode numbers n = 7 and n = 1, respectively, at the locations of s = 0 and q = 1, revealing more clearly that the electromagnetic mode observed is characterized by a high-frequency component of approximately $\omega \sim [150\text{-}200]$ kHz.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To provide further insights into the nature of this high-frequency mode, we perform radially global linear GENE simulations using the final steady-state GENE-Tango profiles corresponding to the case retaining electromagnetic effects and fast particles. The linear results are shown in figure 12, respectively for the growth rate and frequency varying the plasma beta. This has been done by artificially re-scaling the plasma beta of each species by constant value, keeping the same shape as the nominal beta profile (that measures $\beta_e = 0.013$ at mid radius for the steady-state profiles). While the mode at n = 7 is found to be linearly unstable for the steady-state GENE-Tango profiles and gets marginally stable at $\beta_e = 0.008$, the modes at $n = [1,2]$ observed in the nonlinear spectra of the electrostatic potential (see figure 9) are found to be linearly stable for the steady-state profile since at these toroidal mode numbers we only see an electrostatic mode at low frequency. This finding suggest that the electromagnetic high-frequency mode with n = 1 is linearly stable for the final GENE-Tango profile and possibly destabilized nonlinearly.

Zoom In Zoom Out Reset image size

The poloidal cross section of the electrostatic potential for the toroidal mode number n = 7 at $\beta_e = 0.013$ is shown in figure 13. We can observe that the two dominant poloidal mode numbers are m = 6 and m = 7. Additionally, it is worth mentioning that the mode is well localized at the location where s = 0 also in the linear simulations, which is consistent with the nonlinear results presented previously. Figure 13 displays the radial frequency spectrogram of the linear simulations performed for n = 7 at the nominal plasma beta $\beta_e = 0.013$, accompanied by the Alfvén continuum computed using the expression from [52]. Notably, figure 13 highlights that the high-frequency mode observed in the GENE linear simulations is destabilized between two distinct Toroidal AE gaps, with a frequency exceeding that of a reversed shear AE [45, 53, 54].

Zoom In Zoom Out Reset image size

Interestingly, this high frequency electromagnetic mode can be destabilized even in the absence of energetic particles. This is shown in figure 14, where linear growth rates and frequencies are compared for various toroidal mode numbers in simulations both with and without energetic particles. The thermal profiles were initialized to those obtained from the GENE-Tango simulation retaining simultaneously energetic particles and electromagnetic effects. Specifically, figure 14 indicates that the high-frequency mode at n = [6–7] persists even in the absence of fast ions, provided the plasma beta reaches sufficiently large values, akin to what was observed in the GENE-Tango simulations with fast ions and electromagnetic effects.

Zoom In Zoom Out Reset image size

Nonetheless, there are notable distinctions in the high frequency mode behavior with and without the presence of energetic particles. This is shown in figure 15 when looking at growth rate and frequency dependence on the total plasma beta (obtained as the sum of thermal ion, electron and fast ion contributions). Specifically, we observe that the linear growth rates are reduced, and the mode frequency is higher in the absence of fast ions. Moreover, figure 15 shows that the threshold for destabilizing this high frequency mode shifts to higher beta values when fast ions are incorporated into the model. The substantial increase in growth rates for both the high frequency mode and ITG, observed at the same total plasma beta in the absence of energetic particles, bears significant consequences for turbulent fluxes. Specifically, when conducting two stand-alone GENE electromagnetic global simulations with and without energetic particles initializing the same profiles as the GENE-Tango electromagnetic simulation with energetic particles and electromagnetic effects, we notice large differences in the turbulence levels. The simulation without energetic particles displays notably enhanced turbulent fluxes in each channel, as shown in figure 16. These findings imply that in the absence of energetic particles, both ITGs and high frequency modes are more unstable compared to scenarios involving energetic particles, resulting in overall larger turbulent fluxes. This aligns with the analyses presented in the previous sections.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

It is essential to note that the magnetic coil spectrogram analyses for this ASDEX Upgrade discharge do not show signatures of these high-frequency modes. This absence of the mode may be attributed to the weak linear drive observed in the GENE simulations for the steady-state plasma profiles. Specifically, the mode with n = 7 at the location s = 0 exhibits a very small linear growth rate, only slightly exceeding the ITG growth rate, as seen in figure 12. In addition, the other mode observed in figure 9 with n = 1 at the rational surface q = 1 remains linearly stable and is possibly destabilized nonlinearly. As a result, also this mode is characterized by a notably weak drive. Therefore, it is possible that the experimental measurements of these weak signals is particularly challenging. Future investigations will focus on designing similar plasma discharges at ASDEX Upgrade to evaluate the existence of these high-frequency modes through density fluctuation measurements. Unfortunately, such measurements were not available in the plasma core for the discharge analyzed in this paper.

During the electromagnetic GENE-Tango convergence loop performed with energetic particles, in a few iterations Tango increased the plasma pressure at the point to destabilize modes at the rational surface q = 1 above their linear stability threshold. When this happened in the GENE-Tango loop (i.e. when the high frequency mode became strongly unstable), we observed a significant increase in the turbulent fluxes of the energetic particle species. As a result, there is an increase in turbulent fluxes for each thermal species. Figure 17 illustrates this, comparing the turbulent fluxes of energetic particles and thermal ions in two iterations: (i) when the mode at the rational surface q = 1 is marginally stable, and (ii) when Tango increases the pressure beyond the linear stability threshold, leading to additional destabilization. These findings demonstrate the negative impact of these modes on plasma profiles and overall performance when left uncontrolled consistently with previous studies [46, 48, 49, 55].

Zoom In Zoom Out Reset image size

Figure 17 shows the impact of the increased fluxes on the thermal profiles, revealing how Tango responded by reducing the pressure gradients. This adjustment was intended to regulate and mitigate the destabilizing effects induced by energetic particles at q = 1 by reducing their drive. These observations indicate that the simulation aims to achieve a critical equilibrium where the modes at q = 1 are maintained marginally stable or weakly unstable, as these modes corresponds to an increase in the thermal turbulent fluxes. In other words, the profile stiffness in these regimes and regions is determined by the energetic particles if the total plasma beta and its radial gradient reach certain thresholds.

In this section, we compare the GENE-Tango turbulent fluxes with the ones obtained running the flux-tube version of GENE at different locations using the steady-state profiles for the case retaining both fast ions and electromagnetic effects. Specifically, we seek to determine if the flux-tube simulations can appropriately represent the energetic particle stabilizing effects, leading to a matching of the volume integral of the injected sources at different locations, similar to what is observed in the global GENE-Tango simulations. This comparison is particularly relevant when considering computational cost of the simulations, especially when looking at future fusion reactors, where the value of $1/\rho^*$ becomes very large. Under these conditions, local gyrokinetic simulations require significantly lower resolution to accurately represent the relevant bi-normal and radial scales together with the velocity space. However, it is important to note that fusion reactors will have fundamentally very different fast particles (alpha particles) compared to those with moderate energies generated through heating schemes. Therefore, the findings of this section should be applied specifically to fast particle effects of present devices.

The local results are obtained by running GENE flux-tube simulations at seven different radial locations, covering a range of $\rho_\mathrm{tor} = [0.1\text{-}0.7]$. We used a grid resolution of $(n_{k_x} \times n_{k_y} \times n_z) = (256 \times 48 \times 32)$ in the radial (x), bi-normal (y), and field-aligned (z) directions. The velocity space grids consist of $(n_{v_\shortparallel} \times n_w) = (32 \times 16)$ points, covering the velocity range $v_\shortparallel/v_{\mathrm{th},s} = [-3.0, 3.0]$ and $\mu B_0 / T_s = [0,9]$ for each species (s). Periodic boundary conditions are employed along the radial direction allowing us a Fourier decomposition along the radial direction as well. The plasma geometry and energetic particle profiles are kept the same as the ones used for the global GENE-Tango simulation with energetic particles and electromagnetic effects.

In figure 18, we present a comparison between the turbulent fluxes obtained from the flux-tube version of GENE and those from the global version (considering the last five GENE-Tango iterations), alongside the volume integral of the experimental sources. Strikingly, we find a significant discrepancy across all channels, with the flux-tube simulations greatly over-predicting turbulent fluxes for $\rho_\mathrm{tor}\gt 0.2$. Specifically, at $\rho_\mathrm{tor} = 0.3$, where the linearly unstable KBM/AE modes are present, there is a pronounced increase in the turbulent fluxes, particularly for the electrons (both heat and particles), and to a lesser extent, also for the thermal ions.

Zoom In Zoom Out Reset image size

Interestingly, the mismatch between the global and flux-tube fluxes is not limited to the locations with high frequency modes but also extends to areas without these modes. This observation aligns with other analyses demonstrating that energetic particles can influence turbulent fluxes even in radial locations where locally fast ions play a less prominent role [34]. Reference [34] shows differences on the turbulent fluxes qualitatively similar to the ones observed in this section in presence of unstable high frequency modes. In cases where a KBM/AE mode is destabilized in neighboring locations, global modifications of the underlying electric field can impact turbulent fluxes. Such global effects are not adequately captured in flux-tube models, where fluxes at each radial location are treated independently of neighboring positions.

It is also interesting to perform a linear stability analysis for each flux-tube selected for the nonlinear comparison. The resulting linear growth rates and frequencies are displayed in figure 19. As shown in this figure, we observe at $\rho_\mathrm{tor} = 0.3$ a linearly unstable high frequency with frequencies similar to those observed in the linear and nonlinear global GENE simulations. Surprisingly, another linearly unstable mode is found at the location $\rho_\mathrm{tor} = 0.4$ not present in the global simulations, but with a notably lower linear drive compared to the former case. This finding suggests that capturing the radial changes in plasma pressure and safety factor profiles is essential to reproduce the correct drive/damping mechanisms for these high frequency modes. Figure 19 reveals that the flux-tube approximation has a strong impact also on the linear properties of these high frequency modes. One notable limitation—among others—is the inability of local modes to correctly describe the interaction with the shear Alfvén wave continuum, which significantly impact the growth rate of the mode, thus strongly affecting in turn also the nonlinear dynamics that govern the interaction between high frequency modes, ITG turbulence, and zonal flows [45].

Zoom In Zoom Out Reset image size

Notably, our results suggest that global GENE simulations exhibit greater tolerance for marginal instability of these modes, achieving power balance agreement even in the presence of such instabilities. On the other hand, flux-tube simulations impose stricter constraints on the linear stability of these modes. The flux-tube approach, in particular, shows a notable increase in overall injected heat and particle fluxes when a linearly unstable high frequency electromagnetic mode occurs, far exceeding the levels from external sources. These observations strongly advocate for the necessity of a global gyrokinetic treatment to accurately model the effects of energetic particles on turbulence. The limitations of flux-tube simulations hinder their ability to fully capture the intricate interplay between energetic particles and turbulence. Sometimes they yield qualitatively correct results [56].

This paper presents the first gyrokinetic simulations of an ASDEX Upgrade discharge up to the transport time scale retaining energetic particles and electromagnetic effects. Specifically, the study focuses on the H-mode pulse $\#39230$, which exhibits a pronounced peaking of the on-axis ion temperature. The numerical simulations are performed using the newly developed tool GENE-Tango, which integrates the radially global version of the gyrokinetic code GENE with the transport solver Tango. This coupling ensures reliable profile predictions—based on global gyrokinetic calculations—while significantly reducing computation time in comparison to standalone flux-driven simulations.

To investigate the influence of energetic particles on the observed improvement in plasma performance, we conducted four distinct GENE-Tango simulations: (i) electrostatic simulation without fast particles, (ii) electromagnetic simulation without fast particles, (iii) electrostatic simulation with fast particles, and (iv) electromagnetic simulation with fast particles.

We found that the inclusion of energetic particle species in the radially global nonlinear electromagnetic GENE-Tango simulations led to an excellent agreement with experimental plasma profiles. On the contrary, when fast ions and electromagnetic effects were simultaneously neglected in the GENE modeling we observed a significant flattening of the ion temperature profile, consistent with profiles computed using TGLF-ASTRA, which notably under-predicts the on-axis ion temperature.

Thorough analyses conducted in our study have revealed the critical role of high-frequency KBMs/AEs ($\omega \sim [150\text{-}200]$ kHz) in suppressing ion-scale plasma turbulence. These modes play a critical role in the global GENE simulations enhancing zonal flow activity and increasing the shearing rate levels. These findings are consistent with recent experimental observations at DIII-D [57].

Interestingly, as the GENE-Tango simulation advanced, there were a few iterations when Tango raised the plasma pressure beyond a critical threshold leading to a destabilization of these high frequency modes. However, this led to a strong increase of the turbulent fluxes of each species within the GENE simulations. As a result, Tango responded by flattening the thermal profiles to achieve a match of the turbulence fluxes with the experimental power balance leading—potentially—to a large confinement degradation. These findings suggest that the simulation seeks a delicate balance by keeping these modes marginally stable or only weakly unstable.

It is worth emphasizing that the magnetic coil spectrogram analyses of this ASDEX Upgrade discharge did not reveal any indications of the high-frequency modes observed in the GENE simulations. However, it is possible that the particularly weak linear drive of these quite central modes may impose a considerable challenge to detect their signatures through experimental measurements. Future investigations will focus on designing similar plasma discharges at ASDEX Upgrade to evaluate the existence of these high-frequency modes through density fluctuation measurements, not available in the plasma core for the discharge analyzed in this paper.

The mechanism responsible for the observed turbulence suppression bears qualitative similarity to that previously identified in flux-tube simulations where marginally stable energetic particle-driven modes were nonlinearly destabilized, leading to a depletion of the free energy of the ITG turbulence and a subsequent increase in the levels of zonal flows [18, 46]. However, we noted fundamental differences when comparing the global GENE-Tango simulations with results obtained from stand-alone (without Tango) simulations using the flux-tube code at various radial locations. Specifically, global GENE simulations exhibited a higher tolerance for marginal instability of these modes while maintaining power balance agreement. In contrast, flux-tube simulations imposed stricter constraints on the linear stability of these modes, limiting their ability to capture the full complexity of the interplay between energetic particles and turbulence. As a consequence, the plasma pressure profiles in flux-tube simulations would lead to a more pronounced flattening if coupled to a transport code compared to the one observed when running global gyrokinetic simulations to reduce the drive of these modes until they reached a state of marginal stability. These findings underscore the necessity of a global gyrokinetic treatment to accurately model the effects of energetic particles on turbulence, especially in scenarios where energetic particle modes are linearly unstable [56].

The authors would like to acknowledge insightful discussions with P. Lauber, T. Hayward-Schneider, C. Angioni and E. Fable. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200—EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. Numerical simulations were performed at the Marconi and Leonardo Fusion supercomputers at CINECA, Italy. This work was carried out partially using supercomputer resources provided under the EU-JA Broader Approach collaboration in the Computational Simulation Centre of International Fusion Energy Research Centre (IFERC-CSC).