As mentioned in Section 1, several operational scenarios are foreseen for JT-60SA, ranging from full-Ip (5.5 MA) inductive scenarios at low (0.77 × 10²⁰ m⁻³) or high (1.23 × 10²⁰ m⁻³) central density, to advanced inductive (hybrid) scenarios and even high β_N full current drive scenarios with a strongly reversed magnetic shear [6]. In all scenarios, at least 30 MW of NBI will be used (the N-NBI of 500 keV at 10 MW being always used), and the analysis presented here regards scenario 3 (full-Ip inductive at high density) and advanced scenario 4 (hybrid). Scenario 3 operates at a plasma current of 5.5 MA, toroidal magnetic field of 2.25 T, and elongation and triangularity of 1.86 and 0.5, respectively, almost at full NBI power (30 MW out of 34 MW). In this scenario, the counter-current beams are not used, so only 20 MW is foreseen from 10 of the 12 P-NBI units. In terms of normalized performance, the scenario operates at high β _N = 2.6, n_e/n_GW = 0.8 (n_GW being the Greenwald density), and non-inductive current drive fraction of 0.36. Hybrid scenario 4 operates at a lower plasma current of 3.5 MA, toroidal magnetic field of 2.28T (very close to scenario 3), and elongation and triangularity of 1.8 and 0.4, respectively, with the same 30 MW of NBI power. However, it also obtains additional 7 MW of ECH, operates at higher β _N = 3, and has similar n_e/n_GW = 0.8 but a higher non-inductive current drive fraction of 0.58. The main ion species of the plasma in both scenarios is deuterium, and the only impurity present is carbon (JT-60SA first operates with a carbon wall) with a corresponding Z_eff on-axis of 1.6 (∼3 at the edge) for scenario 3 and 1.55 (∼2.6 at the edge) for hybrid scenario 4.

In this work, the equilibrium and kinetic profile (densities/temperatures) references for the two scenarios are taken from variants of previously predictive modeling results from the CRONOS suite [12, 33], which used full heating power (NEMO [34]/SPOT [35] for simulating NBI) for the inductive scenario at high density instead of the reduced 30 MW. However, the transport models used were the same, i.e., GLF23 [36] for particle/energy transport in scenario 3 and CDBM [37] for energy transport in hybrid scenario 4 ([12] provides a more detailed rationale for the appropriateness of such transport models). Figure 1 shows the plasma equilibrium shape and poloidal magnetic flux surfaces, as well as the plasma species densities (including fast ion density), temperatures, and safety factor profiles, for both scenarios. The radial coordinate ( ρ t o r _ n o r m ) used in CRONOS for densities, temperatures, and safety factor profiles is related to the normalized toroidal magnetic flux Φ_N as ρ tor _ norm = Φ N . The foreseen fast ion density for scenario 3 is clearly off-axis, which is easily justified by the beam deposition and slowing down, as confirmed later when analyzing the results from the Monte Carlo orbit-following code ASCOT [38]. The summary overview for hybrid scenario 4, also shown in Figure 1, shows that the safety factor is q∼0.83 on the axis and, contrary to scenario 3, features a finite magnetic shear in the plasma core. A clear ITB region ( Φ N ∼ 0.3 − 0.5 or ψ N ∼ 0.4 − 0.6 ) is identified and stems from the stabilizing effect on turbulence from fast ions, as accounted for in the CDBM transport model [12, 39].

In order to obtain the energetic particle distribution for each plasma scenario, the ASCOT code suite [38] was used. ASCOT is an orbit-following code, which uses the Monte Carlo approach to simulate the motion of charged particles in 3D magnetic fields and, together with appropriate collision operators, simulate the beam ionization (through the BBNBI code [40]), slowing down, and thermalization of the energetic ions toward a steady state. The simulations took into account the carbon impurity species in addition to the main deuterium ions from the plasma, and all energy components for positive-ion-based beams were included, i.e., E, E/2, and E/3. A thermalization “standard” end condition of 1.5 times the thermal energy was considered in the simulations. Separate runs were conducted for the positive-ion-based NBI (P-NBI) and negative-ion-based NBI (N-NBI) systems, powered at the foreseen nominal values (10 × 2 MW and 1 × 10 MW, respectively) to allow for an independent evaluation of the contribution of each NBI system to the Alfven eigenmode stability. P-NBI and N-NBI yield quite different deposition profiles, considering both the injection energies (85 keV for P-NBI and 500 keV for N-NBI) and geometries (predominantly perpendicular for P-NBI and only tangential for N-NBI). The different plasma densities in both scenarios under consideration also lead to evident differences, with the P-NBI deposition and slowing down clearly dominant closer to the plasma edge for scenario 3, as shown in Figure 2A (particle distribution integrated over velocity space). Since only two out of the 10 P-NBI units used have a tangential injection, the velocity pitch angle distribution is mostly symmetric in the parallel velocity, as shown in Figure 2C (particle distribution integrated over the spatial domain). The three different energy components are quite evident, showing a significant fraction of ions at their birth energies.

The N-NBI beams, on the other hand, penetrate deeper inside and yield roughly an annular distribution over the RZ plane (see Figure 2B) due to a significant fraction of a co-passing population (see Figure 2D). Figure 2C clearly shows that much like the P-NBI distribution, the distribution exhibits characteristic bump-on-tail features close to the birth beam energy of 500 keV, in the v_//>0 region, due to tangential injection. From these observations, one can already anticipate the potential effect on the stability of the MHD shear Alfven waves associated with this scenario. Although the P-NBI distribution is primarily relevant for edge-localized modes, where the spatial gradient of the distribution function (proxy for the derivative with respect to P_phi) is largest in magnitude, the N-NBI distribution can potentially destabilize a larger pool of modes yet dominantly located from midradius outward where ∂ F EP ∂ P p h i < 0 is anticipated. The maximum of the flux surface-averaged N-NBI distribution is, indeed, located at s = ψ N ∼ 0.42 , where ψ N is the normalized poloidal magnetic flux. Averaged over velocity space, this yields a fast ion density of the order 6 × 10 17 m − 3 , which is roughly 0.6% of the local bulk plasma density. The bump-on-tail features of the N-NBI distribution can potentially assist the drive of some modes unstable by co-passing particles through ∂ F EP ∂ E > 0 , something not as easily accessible to the P-NBI distribution except for deeply trapped particles.

For hybrid scenario 4, owing to the lower plasma density and similar electron/ion temperature from midradius outward, both P-NBI and N-NBI beams are expected to penetrate deeper inside the plasma volume, so a larger pool of modes, e.g., core localized modes, might be expected to interact with the energetic particle population. This is particularly relevant if we consider the ITB region ( Φ N ∼ 0.3 − 0.5 or ψ N ∼ 0.4 − 0.6 ) since fast ions are believed to play an important role in the stabilization of turbulence in this scenario [12], and fast ion redistribution resulting from destabilized Alfvén eigenmodes could potentially degrade ITB. Figure 3 shows the ASCOT distribution function over RZ and velocity spaces for the P-NBI (A, C) and N-NBI (B, D) settings, respectively. Although less evident, the singled-out bump on the tail of the N-NBI distribution close to the birth energy is again observed. As previously anticipated, the fast ions penetrate deeper into the plasma, yielding a non-negligible particle density on the axis.

To investigate the linear interaction of the energetic beam particle populations with the possible pool of MHD Alfven eigenmodes in each scenario, a workflow [41] relying on the fixed-boundary high-resolution equilibrium code HELENA [42], the ideal incompressible MHD linear stability code MISHKA [43], and the hybrid MHD-drift kinetic code CASTOR-K [44, 45] is used. CASTOR-K uses a perturbative approach for the perturbed distribution function of the energetic ions associated with a given normal eigenmode of the MHD kernel (ideal/resistive) and, therefore, cannot be used to study energetic particle modes or any non-linear behavior of the mode amplitude and/or frequency. The perturbed Lagrangian used in the simulations assumes the MHD eigenmode to be incompressible; thus, this precludes the analysis of any modes where coupling with the sonic branch of compressible MHD occurs, e.g., BAEs. The model does, however, include finite particle orbit width effects and can, therefore, address all different types of orbits and orbit topological regions, e.g., trapped-passing boundary. Arbitrary distribution functions in the (P_phi,E,μ) constant of motion phase space (parametrized or from tabulated data) can be used in the model, which grants CASTOR-K the capability to address arbitrary fast ion sources such as fusion alphas or NBI/ICRH additional heating in any possible configuration, as well as thermal background ions. Since CASTOR-K code operates as an eigenvalue solver, it is much more adequate to analyze the damping rates of stable modes than initial value codes and sub-dominant modes at the same toroidal mode number and neighboring frequencies. It is also far less computationally intensive than any equivalent initial value code. The model only addresses kinetic effects from the thermal or energetic particle populations, which means that fluid damping due to conversion of the mode energy into strongly damped kinetic Alfvén waves (radiative damping) and collisional electron damping are not directly included in CASTOR-K but can be suitably calculated by a resistive MHD code (see Ref. 44). The identification of all AE solutions of the linearized MHD equations with MISHKA is carried out over a range of toroidal mode numbers (n = 1–25) and with normalized frequency ω/ω_A within 0.01 and 2.7, where ω_A is the Alfvén frequency on axis. The upper limit in frequency is a reasonable compromise, considering the maximum beam energy of 500 keV, the particle energy at the Alfven velocity (on axis) of E_A∼136 and 150 keV for scenario 3 and hybrid scenario 4, respectively, and the ability to identify not only TAEs but also EAEs and NAEs. The upper limit in the toroidal mode number is somewhat optimistic, considering the inherent limitations to the drift-kinetic model over which the perturbative analysis of CASTOR-K takes place. To be strictly valid, the perpendicular wavevector of the perturbation k_⊥ should satisfy k ⊥ r L < 1 , where r L is the typical Larmor radius of the energetic particle. This translates into a maximum toroidal mode number N max ∼   s / q m o d e / r L / a with s as the normalized radial coordinate and the lower script “mode” as the radial location, where a particular mode has the energy density of the highest mode (hereafter referred to as “mode location”). Taking an indicative average particle energy E∼250keV, midradius mode location, and q∼1 yields N_max∼15. As the mode locations and dominant particle resonances are largely unknown a priori, we consider this value of N_max as a more realistic upper limit than the assumed n = 25, although, in the following sections, the results also cover modes up to n = 25. The pool of modes within the prescribed toroidal mode number and frequency range is then reduced by discarding all modes that cross the Alfvén continuum as the interaction of the modes with the continuum is considered to strongly damp the modes [46]. In all CASTOR-K simulations, a 201 radial-by-256 poloidal plasma equilibrium grid was considered, and a 32-poloidal harmonic spectral resolution was considered when calculating the MHD eigenmodes. The minimum poloidal mode number depended on the toroidal mode number and on the safety factor q(s) on the axis to ensure that the most relevant resonance/gap modes, e.g., q=(m-1/2)/n up to s∼0.8, were covered for the modes of interest, i.e., below n∼15. CASTOR-K is executed independently over all this set using three different distribution functions: the thermal ion Maxwellian population, and the P-NBI and the N-NBI beam ion populations. When using the beam ion populations in CASTOR-K, these are first mapped appropriately to the (P_phi,E,μ) constant of motion space from the ASCOT output.

In this scenario, a total of 213 AEs were identified, following the criteria detailed in Section 2.3, covering most of the frequency range and yielding modes localized from deep inside the core up to the plasma boundary (see Supplementary Figure S1 in Supplementary Material). When considering both the effect of the neutral beam energetic ions (PNBI + NNBI sources) and the thermal ions on the mode stability, it is observed, as shown in Figure 4, that all modes are found to be stable independent of their frequency (A) or location (B). Thermal ion Landau damping dominates and, whenever beam ions drive the modes unstable, it is roughly one order of magnitude lower, i.e., γ ω A ≲ 0.2 % (see discussion in the following paragraph) for the N-NBI beams and γ ω A ≲ 0.02 % for the P-NBI beams. Although the figure may suggest that there are fewer than 213 modes, there is a multiplicity of modes that have frequencies so similar that they appear to be overlapping. This is to be expected in low-shear regions of the plasma, as is the case of the plasma core, where multiple modes can be found within a single Alfven gap and with increasing radial wavenumber on the mode harmonics as the mode frequency approaches either lower or upper ends of a given Alfven gap [47].

Although the combined effect of the thermal ions and energetic beam ions leads to stabilization of all identified Alfvén eigenmodes, it is important to investigate which particular modes are destabilized by the beam ions and identify the underlying features in the beam ion distribution function or wave–particle resonances that cause such destabilization. Since the P-NBI beam ions have a much less relevant impact on mode stability compared to N-NBI, here, we focus only on the driven/damped modes by the N-NBI beam ions. The overview results on the growth/damping rates stemming from N-NBI fast ions alone are shown in Figure 5, and one observes that, at most, one obtains a normalized growth rate γ ω A ∼ 0.2 % . In order to obtain the net drive when considering both thermal ion damping and potential NBI drive, we found that the threshold for fast ion density (n_fast-crit) would need to be five times than expected in the scenario reference (n_fast-ref), i.e., n_fast-crit> 5n_fast-ref, translating to 3% in the relative density of fast ions to thermal ions. From Figure 5, a first inspection confirms, as anticipated in Section 2.2, that the modes driven unstable are located at the s > 0.5 plasma domain where the fast ion distribution function has a negative radial gradient. However, this is a rather simplistic assessment since a) the magnetic flux of a magnetic surface is just a proxy for the canonical toroidal angular momentum P_phi and b) the drive for the mode excitation can equally arise from a positive gradient in energy of the distribution function if relevant wave–particle resonances are found at the bump-on-tail regions of phase space of the distribution (see Figure 2). Likewise, mode damping may well arise from local positive gradients of the distribution in P_phi, i.e., ∂ F EP ∂ P p h i > 0 , in addition to the usual ∂ F EP ∂ E < 0 contribution.

As shown in the following sections for some representative modes, in the majority of the cases where the beam ions drive the modes unstable, wave–particle resonances located around the 500-keV birth energy of the beam ions play a relevant role in mode destabilization. In addition, the drive coming from ∂ F EP ∂ E > 0 may be significant, as anticipated in bump-on-tail distributions (see Figure 2B when particle energy is close to maximum). Mode damping is also observed not to arise exclusively from ∂ F EP ∂ E < 0 , with wave–particle resonances crossing regions of phase space, where ∂ F EP ∂ P phi > 0 , particularly relevant in the plasma core (see Figure 2B).

A typical example of net mode drive by the energetic particle distribution of the AE spectra is evidenced by the n = 9 TAE mode with ω/ω_A = 0.302. As shown as follows, it features most of the potentially different stability trends associated with the beam distribution from the N-NBI system. The mode is not radially localized but extends from midradius up to the plasma edge, owing to broad poloidal harmonic spectra as is usual when the magnetic shear is not small (see Figure 6A). For dominantly co-passing particles, the wave–particle resonances are located in the ∂ F EP ∂ P p h i < 0 region of phase space, and thus, the EP distribution destabilizes the mode, with noticeable energy exchange around the birth energy of the injected neutral atoms (500 keV). This is shown in Figure 6B, where the distribution function and the dominant wave–particle resonance (magenta dashed curve) for the co-passing particle population (labeled as sigma = 1) for normalized magnetic moment λ = μ B E = 0.4725 are shown. As shown in Figure 6C, the highest energy transfer between particles and waves for the co-passing particle population (labeled as sigma = 1) occurs close to the bump-on-tail features at E∼500 keV, which clearly yields a destabilizing contribution stemming both from ∂ F EP ∂ E > 0 and ∂ F EP ∂ P p h i < 0 contributions. We note that in all the figures showing the distribution function over phase space (E,P_phi) for a given λ, the apparently coarse resolution of the domain is just a visual artifact arising from the symbol size used in the scatter plot and from the adaptive grid refinement used by the CASTOR-K code to calculate the total number of interacting particles, i.e., the grid need not have the same resolution in all phase space regions but rather be refined only where it matters the most.

As λ is increased (increased relative weight of perpendicular velocity to the total particle velocity or energy), the wave–particle resonance gradually shifts to the lower end (in P_phi) of the distribution where a net stabilizing effect ensues, assisted by a positive gradient in P_phi. Similar to the case at λ = 0.4725, the bump-on-tail features contribute to the energy exchange (see Figures 6D, E for the resonance map, distribution function, and energy exchange at λ = 0.7665). The guiding center orbits exchanging the most energy with the mode are located within the innermost regions of the mode eigenfunction, coinciding with the highest EP density from N-NBI, as evidenced in Figure 7. It should be noted that as λ increases, the shape and magnitude of the distribution function also change, and the net result, once integrated over λ, turns out to be destabilizing.

The n = 12 mode considered is localized at the core, has a normalized frequency ω/ω_A = 0.248, has a dominant m = 13 poloidal harmonic, and is the mode (n < 15) with the highest damping by N-NBI. The electrostatic potential is shown in Figure 8A. Owing to the quite central localization of the mode, wave–particle resonances above moderate energy ( E ≳ 50 keV are expected to be accessible primarily to counter-passing particles only. From the CASTOR-K calculations, stability is indeed dominated by counter-passing particles, with ∂ F EP ∂ P p h i > 0 contributing significantly to the overall damping of the mode. Figure 8B shows the distribution of counter-passing particles (labeled as sigma = −1) for the normalized magnetic moment λ = 0.6825, where the absolute value of the differential in energy exchange (dW_hot/dλ) has its maximum in the absolute value. The dominant resonance region (damping) is also highlighted in magenta and occurs at P_phi∼0.28 eVs and E∼90keV, corresponding to the v_A/3 resonance. When looking for the highest energy exchange over phase space, at the dominant wave–particle resonance, one easily confirms that the positive gradient in P_phi term of Eq. 2 (Figure 8D) contributes much higher to mode damping than the negative gradient in energy term (Figure 8C). A cross section of the electrostatic potential of this n = 12 TAE over the (R,Z) space is shown in Figure 8E, together with the two possible orbits for the constant of motion (P_phi = 0.28 eVs, E = 90 keV, λ = 0.6825), where damping from the energy gradient is largest. As shown in Figure 8B, these represent a small fraction of phase space for the energetic particle distribution, and thus, the overall damping is unsurprisingly smaller compared to the N-NBI drive of the n = 9 mode presented in Section 3.2. At lower values of λ, a similar pattern of wave–particle resonances is observed, although at slightly lower energy and toroidal angular momentum and with v_A/3 still acting as single/dominant resonance.

To model a possible uncertainty in the q-profile, a variation of ±1% was considered on the total plasma current around the reference value. In the ensuing analysis, only the plasma equilibrium is recalculated, i.e., the same thermal density/temperature profiles and energetic particle distribution functions as the reference scenario were considered. Although these slightly different scenarios are, thus, not fully self-consistent, the emphasis of the analysis is on the qualitative impact on the mode spectra and possible impact of the drive/damping nature of the interaction between the energetic distribution functions and the modes. The calculated equilibria for the three scenarios have, as expected, negligible differences in the flux surfaces mapping in (R,Z), considering the small variation in plasma current and/or in the equilibrium pressure, i.e., no change in the Shafranov shift is observed and the relative change in the q-profile is almost hardly noticeable (see Figure 9).

Given the shape of the q-profile, which exhibits very low shear in the core, and the fact that Alfvén eigenmode spectra are highly dependent on the local q-profile features, it is anticipated that the effect of the small changes in plasma current on the eigenmode spectrum is most prominent for modes localized at the plasma core. This is confirmed when the total spectra of modes that do not cross the Alfvén continuum are plotted, as shown in Figure 10. Although modes located outside the low-magnetic shear region have very small variations in frequency, core modes show significant variations in frequency and, in some cases, are eliminated from the analysis since, albeit small, the change in the safety factor can cause the frequency of some modes to touch neighboring branches of the modified Alfvén continuum.

Regarding MHD stability, there is little change overall since thermal ion landau damping still dominates the wave–particle interaction, and both 1% increase and decrease in plasma current still yield stabilized modes (see Supplementary Figure S2). Similar to the reference equilibrium, the variants of 1% in plasma current also evidence mode destabilization by the N-NBI beam ions. It was verified in the simulation results that the unstable modes typically have mode numbers n = 3,7–10, their poloidal harmonic content is localized mostly at s > 0.6, and the mode drive comes primarily from the co-passing population with normalized magnetic moment λ = [0,0.65]. In the (P_phi,E) space, the drive arises once more from wave–particle resonances at the downslope region of the distribution function, i.e., where ∂ F EP ∂ P p h i < 0 , extending up to the birth energy of the beam ions. The bump-on-tail features of the distribution assist in the driving/damping of the mode whether particle energy is slightly smaller/higher than the birth energy, similar to the results already obtained for the reference scenario.

In this scenario, a total of 358 AEs were identified, following the criteria detailed in Section 2.3. Similar to the full Ip inductive scenario 3, very few modes (here, only EAEs with a frequency of the order of the Alfvén frequency) were found in the region where the ITB is located, i.e., s = ψ N ∼ 0.4 − 0.6 , as shown in Figure 11B. One notices that there is a reasonable number of modes localized very close to the axis (see Figure 11B). Such modes are, in some cases, ill-resolved radially and not considered in the MHD stability analysis. The scarcity of modes within the ITB region is encouraging as the absence of Alfvén eigenmodes in that region eliminates the possibility of fast ion transport induced by a resonant interaction with these modes. However, it also raises the question of why no modes are present within ITB. Is the bulk plasma pressure gradient so large that no shear Alfvén wave solutions are permitted in the linearized MHD equations [48], or do some of the existing solutions cross the Alfvén continuum and are, thus, excluded from the analysis?

In order to explore the former hypothesis, the ideal MHD spectra were re-calculated using two scaled variants of the reference equilibrium where the plasma beta (pressure profile) was scaled down to 70% and 50%, respectively. The analysis evidenced that there is no increase in the number of modes in the radial range of ITB (which do not eventually cross the continuum) when decreasing plasma beta (see Figure 12 showing only modes with a normalized frequency below 1). Moreover, all modes located at s > 0.6 with a frequency close to ω/ω_A = 1, corresponding to the EAEs shown in Figure 11, are still observed. One should note that the imposed beta scaling of the equilibria does impact the safety factor q(s) profile even though the total plasma current only shows a variation below 0.1%. This causes slight variations in the mode locations, similar to the plasma current scan performed for scenario 3 (see Figure 10).

In order to confirm that the local absence of TAE modes is a continuum effect, it is instructive to plot the Alfvén continuum for a representative subset of the toroidal mode spectra, e.g., n < 10, to avoid cluttering the plot. As shown in Figure 13, since magnetic shear is not small and considering the equilibrium mass density and q-profile, most TAE gaps in the range s ≳ 0.25 are relatively close to each other and do not have gap extrema aligned in frequency. This implies that gap modes either touch the continuum at inner radii, e.g., 0.25 < s < 0.6 gap modes, or touch the continuum at s∼0.9. Two representative examples (n = 5 and n = 7) are shown in Figure 14. When n = 5 with ω/ω_A = 0.525 (Figure 14A), although the frequency matching the Alfvén continuum occurs only near the edge, it is unavoidable that the relevant (“resonant”) poloidal harmonics of the eigenfunction exhibit sharp local variations over radius and are, thus, discarded. When n = 7 with ω/ω_A = 0.401 (see Figure 14B), the continuum crossing occurs at a radial location where the poloidal components of the eigenfunction still hold significant energy, the sharp radial variations in the eigenfunction are much more noticeable, and thus, the mode is trivially discarded. One notes that there are other eigenfunction solutions for n = 7 (extensible to other toroidal mode numbers) which feature dominant m=(n-1),n poloidal harmonics, which resonate strongly with a continuum close to the TAE gap centered at s∼0.2–0.3. In such cases, these harmonics overshadow the remaining poloidal harmonics (barely noticeable in the eigenfunction radial profile). Figure 13 shows that as the toroidal mode number increases, at 0.25 < s < 0.6, it only allows for sufficiently localized (high-n) of odd parity modes at the top of the gap ω/ω_A ∼[1,1.3], in agreement with the results shown in Figure 11B.

When the contribution from N-NBI, P-NBI, and thermal ions to the mode stability is calculated, similar to scenario 3, all modes are observed to be stable (Figure 15A). However, in scenario 4, there is a more sizeable number of modes which are driven unstable by N-NBI, as shown in Figure 15B, showing the stability overview results including only the contribution from N-NBI fast ions. Such modes are primarily in the TAE and EAE frequency ranges (not shown) and are all located in the outer regions (s > 0.6) of the plasma. Although thermal ion landau damping is smaller than at the core, it is enough to overcome the beam drive.