scieee AI-readable full text Open interactive document viewer

Hybrid kinetic-MHD modeling of alpha-driven TAEs in the SPARC tokamak

Tinguely, R. Alexander; González Martín, Javier; Todo, Yasushi

Abstract

As the magnetic confinement fusion community prepares for the next generation of fusion devices and burning plasmas, there is still a question of whether fast ions (FIs) will drive MHD instabilities, causing significant redistribution or even loss of FIs, thereby leading to reduced plasma performance and possibly threatening the integrity of the first wall. In this paper, we explore the existence and stability of toroidicity-induced Alfvén eigenmodes (TAEs) in the > 100 MW , Q ∼ 9 -11 DT-fusion power ‘Primary Reference Discharge’ (PRD) of the SPARC tokamak; the PRD has a relatively low on-axis alpha pressure, β α 0 ≈ 0.6 % , due to the high magnetic field strength, B 0 = 12.2 T . A scan in toroidal mode number is performed in the vicinity of the estimated ‘most unstable’ modes, n ≈ 5-20, with the linear eigenvalue code NOVA-K and nonlinear initial-value code MEGA. Both codes identify the same (even) n = 10 TAE located near q = 1 with frequency f ≈ 360 kHz and alpha drive γ / ω ≈ + 0.6 % . While MEGA evaluates this mode to be marginally unstable for the nominal alpha pressure, NOVA-K instead identifies a higher frequency (odd) n = 10 TAE as marginally destabilized; different evaluations of radiative damping are likely the cause of this discrepancy. These results indicate that AEs may be only marginally unstable for the highest performing SPARC PRD, at least for the q profile explored here. They also serve as a starting point for further scans, inclusion of FIs from auxiliary heating systems, and exploration of AE-induced FI transport, as well as a guide for diagnostic measurements of these n ≈ 10 AEs.

Full text

PAPER • OPEN ACCESS Hybrid kinetic-MHD modeling of alpha-driven TAEs in the SPARC tokamak To cite this article: R.A. Tinguely et al 2025 Nucl. Fusion 65 036021 View the article online for updates and enhancements. You may also like Simultaneous measurements of unstable and stable Alfvén eigenmodes in JET R.A. Tinguely, J. Gonzalez-Martin, P.G. Puglia et al. - Alfvén eigenmode stability analysis and energetic particle transport prediction for CFETR hybrid scenario Yunpeng ZOU, , Minyou YE et al. - Measurement and calculation of Alfvén eigenmode damping and excitation over a full toroidal spectrum J. Sears, R.R. Parker, J.A. Snipes et al. - This content was downloaded from IP address 150.214.182.230 on 06/03/2025 at 14:19 International Atomic Energy Agency Nuclear Fusion Nucl. Fusion 65 (2025) 036021 (15pp) https://doi.org/10.1088/1741-4326/adaf40 Hybrid kinetic-MHD modeling of alpha-driven TAEs in the SPARC tokamak R.A. Tinguely1,a,∗, J. Gonzalez-Martin2,aand Y. Todo3 1Plasma Science and Fusion Center, Massachusetts Institute of Technology, Cambridge, MA, United States of America 2Universidad de Sevilla, Seville, Spain 3National Institute for Fusion Science, Toki, Gifu 509-5292, Japan E-mail: [email protected] Received 16 August 2024, revised 20 November 2024 Accepted for publication 28 January 2025 Published 17 February 2025 Abstract As the magnetic confinement fusion community prepares for the next generation of fusion devices and burning plasmas, there is still a question of whether fast ions (FIs) will drive MHD instabilities, causing significant redistribution or even loss of FIs, thereby leading to reduced plasma performance and possibly threatening the integrity of the first wall. In this paper, we explore the existence and stability of toroidicity-induced Alfvén eigenmodes (TAEs) in the >100MW, Q∼9–11 DT-fusion power ‘Primary Reference Discharge’ (PRD) of the SPARC tokamak; the PRD has a relatively low on-axis alpha pressure, βα0≈0.6%, due to the high magnetic field strength, B0=12.2T. A scan in toroidal mode number is performed in the vicinity of the estimated ‘most unstable’ modes, n≈5–20, with the linear eigenvalue code NOVA-K and nonlinear initial-value code MEGA. Both codes identify the same (even) n=10 TAE located near q=1 with frequency f≈360kHz and alpha drive γ/ω ≈+0.6%. While MEGA evaluates this mode to be marginally unstable for the nominal alpha pressure, NOVA-K instead identifies a higher frequency (odd) n=10 TAE as marginally destabilized; different evaluations of radiative damping are likely the cause of this discrepancy. These results indicate that AEs may be only marginally unstable for the highest performing SPARC PRD, at least for the qprofile explored here. They also serve as a starting point for further scans, inclusion of FIs from auxiliary heating systems, and exploration of AE-induced FI transport, as well as a guide for diagnostic measurements of these n≈10 AEs. aShared first authorship. ∗Author to whom any correspondence should be addressed. Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. 1741-4326/25/036021+15$33.00 Printed in the UK 1 © 2025 The Author(s). Published by IOP Publishing Ltd on behalf of the IAEA Nucl. Fusion 65 (2025) 036021 R.A. Tinguely et al Keywords: Alfvén eigenmodes, stability, SPARC, NOVA-K, MEGA (Some figures may appear in colour only in the online journal) 1. Introduction It has been observed in many tokamaks that fast ions (FIs), with velocities of order the Alfvén speed, can destabilize Alfvén eigenmodes (AEs) and that, in turn, these AEs can be correlated with FI transport and deconfinement; for examples, see [1–3] and others. Such FI redistribution in phase space, caused by wave-particle resonances and energy exchange, could thus affect—and possibly degrade—plasma heating, which is usually the goal of the FIs in the first place. Of particular interest is the 3.5MeV alpha particle born in the deuterium– tritium (DT) fusion reaction, which slows down from its birth velocity v0≈1.3×107ms−1via collisions. Both TFTR and JET, the only two tokamaks to have operated with DT fuel, have some evidence of alpha-driven AEs [4–6]. However, in those experiments, the alpha populations—and power—were relatively small, with plasma fusion gains Q<1. An area of active study, therefore, is the prediction of AE stability, alpha drive and transport in the burning plasmas (Q>5) of future DT devices. Many analyses have been carried out for ITER (see [7,8] and others); this paper focuses on the SPARC tokamak [9], which is currently under construction by Commonwealth Fusion Systems. SPARC is a compact device with major and minor radii R0=1.85m and a=0.57m, respectively; its high magnetic field strength B0=12.2T and plasma current IP=8.7MA are enabled by high temperature superconducting technology [10]. These strong fields, together with relatively large plasma densities O(1020 m−3), yield a high Alfvén velocity vA≈9×106ms−1, which yet is still less than the alpha birth velocity. Previous works have considered AE stability and related transport in SPARC: in [11], the authors investigated trends of parameters relevant to AE stability with the magnetic field strength B: alpha pressure; AE mode structures, growth rates, resonance conditions; and more. A scan in Bwas performed using a suite of codes: HELENA [12] to compute the magnetic equilibrium, MISHKA [13] to identify AEs, and CASTOR-K [14] to evaluate AE growth (and damping) rates. Importantly, a notional set of profiles (scaled from an Alcator C-Mod discharge) was used, simply because the paper was published before the SPARC Physics Basis [9]. In addition, only alpha drive and ion Landau damping were included in the stability calculation, without any finite Larmor radius (FLR) effects. A drift kinetic theory of alpha transport was derived in [15] and applied to cases of toroidal magnetic field ripple and a generic Toroidicity-induced AE (TAE) in SPARC. Scott et al [16] complemented this work with a numerical investigation of ripple-induced alpha transport, comparing the ASCOT5 [17] and SPIRAL [18] orbitfollowing codes. Further numerical analysis of TAE-induced alpha transport is underway, but will be left for a future publication. This work presents the most thorough analysis to date of alpha-driven TAEs in SPARC. We improve upon [11] by analyzing a realistic magnetic geometry and set of plasma profiles for SPARC [19], given in section 2. Furthermore, section 3 complements the stability analysis of [11] using the linear eigenvalue solver and stability code NOVA-K [20–22] which includes additional drive and damping mechanisms as well as FLR effects. We then compare linear results to those from the nonlinear code MEGA [23] in section 4. Finally, sections 5 and 6provide a discussion and summary. 2. Modeling inputs The magnetic geometry, plasma profiles, and FI distribution functions are common inputs to both codes. Specifically, the scenario of interest is SPARC’s planned highest performing ‘Primary Reference Discharge’ (PRD), a double-null H-mode DT plasma with maximum fusion power Pfus =140MW [9, 19,24]. Of course, a plasma with maximal DT fusion power also has the highest alpha particle generation rate, which is why the PRD is the most relevant scenario for initial exploration of alpha-driven AEs in SPARC. Other plasma scenarios could lead to even greater alpha drive or a wider variety of FIrelated instabilities, as discussed later; however, we decide to focus on the PRD, the SPARC scenario most explored to date. The first high fidelity simulations of core plasma performance for the PRD were carried out in [19]; figure 7 within shows the time-evolution of the standard baseline discharge, including the sawtooth instability. During the quasi-quiescent period between sawtooth crashes, there is minimal change in the on-axis and volume-averaged electron densities, while the on-axis electron and ion temperatures vary by roughly ∼10%. Of course, the sawtooth crash affects the temperatures greatly, reducing them by almost a factor of 2, along with the requisite change in the core safety factor (q) profile, i.e. the existence and movement of the q=1 surface. It is primarily computational resources that limit the following simulations of sections 3and 4to a single time slice, and we simply choose a time shortly before a sawtooth crash, which is explored further in [19] and in many follow-on publications. It is difficult to quantify uncertainties introduced by selecting this single time; however, we posit a rough ∼10% uncertainty in relevant plasma parameters (neglecting the effect of sawteeth) and discuss the resulting uncertainty propagation in later results. Electron and ion density and temperature profiles from TRANSP [25–27] simulations are shown in figure 1for the single time selected, along with the qprofile before the sawtooth crash. In the figures, ψNis the normalized poloidal flux, and √ψNis approximately the normalized minor radius (i.e. √ψN=ρpol ∼r/a). A poloidal cross-section of √ψNis depicted in figure 2, along with the reduced-resolution, field-aligned 2 Nucl. Fusion 65 (2025) 036021 R.A. Tinguely et al Figure 1. (a) Density and (b) temperature profiles, from TRANSP, for electrons (dotted) and ions (dashed) in the SPARC primary reference discharge (PRD), plotted vs normalized minor radius √ψN, with ψNthe normalized poloidal flux. Note that ni=nD+nTis the combined D+T ion density. The safety factor profile (solid) is also shown in (a). Figure 2. A poloidal cross-section of the SPARC PRD magnetic equilibrium, showing the square root of normalized poloidal flux √ψNtogether with the field-aligned grid used for post-processing in MEGA, with a reduced resolution for visualization. grid used in MEGA. Uniform Zeff =1.5 and (nD+nT)/ne= 0.85 were assumed in the TRANSP simulations, with an impurity mix that consists of the Helium-3 (He3) minority, tungsten (W), and a lumped low-Z impurity (F). An effective density, incorporating dilution as well as the two main ion species of different masses, is used in the following simulations of sections 3and 4. We note, however, that the correction to the Alfvén speed and frequencies due to this dilution is <5%. Here, it is important to discuss the qprofile used for modeling. From the analysis of [19,24], the sawtooth period is ∼1s, similar to the energy confinement time of ∼0.8s and ∼4×the alpha slowing down time, meaning that the plasma is approximately in steady state before the sawtooth crash. At the selected time slice, the rational surface q=1 is located at mid-radius, around √ψN≈0.5, as seen in figure 1(a). Because this time is just before the sawtooth crash, the q=1 surface is likely near the largest radius at which it will exist, at least during the flat-top plasma current of the SPARC PRD scenario. For this reason, the chosen time slice has an ‘extreme’ q=1 location, which affects the existence of AEs and their stability. We discuss the impact of this choice in the following sections, and other qprofiles should be investigated in future work. The SPARC PRD will be externally heated by ∼11MW of ion cyclotron radio frequency (RF) power, utilizing He3 as the minority ion species [28]. In the TRANSP simulation used here, the total fusion power is ∼111MW and Ohmic heating is ∼1MW, meaning that alpha power will contribute ∼22MW of self-heating to this Q≈9 plasma. The normalized pressure (β) profiles for the main and FI species are shown in figure 3(a). All values are relatively low due to the high magnetic pressure. The central β0is actually larger for the He3 RF minority than alphas, i.e. ∼1%compared to ∼0.6%, respectively. However, the alpha profile is much broader, with a fairly constant radial gradient to mid-radius (√ψN≈0.1–0.5), whereas the He3 FI gradient is steepest within the core (√ψN<0.2), as shown in figure 3(b). Note that the small positive gradient for the alphas at the plasma center (√ψN<0.1) is not real, but due to poor marker statistics in the NUBEAM simulation [29]. The alpha population, computed by NUBEAM, is isotropic in pitch angle and exhibits the usual slowing down distribution in velocity. The He3 population from TORIC-FPP [30, 31], on the other hand, is quite anisotropic. The total, parallel, and perpendicular ‘effective’ temperature profiles are plotted in figure 3(c). Note the similarity to figure 13 in [16]. These 3 Nucl. Fusion 65 (2025) 036021 R.A. Tinguely et al Figure 3. (a) Normalized pressure (beta) profiles, from TRANSP, for electrons (dotted), ions (dashed), Helium-3 RF minority ions (dot-dashed), and alphas (solid) in the SPARC PRD. (b) Gradients of He3 and alpha betas from (a). (c) Total (solid), parallel (dotted), and perpendicular (dashed) effective temperatures for the He3 RF minority population vs normalized minor radius √ψN. data are used in section 3to evaluate the RF contribution to linear AE stability with NOVA-K where, at each radial location, the RF tail is modeled with an anisotropic equivalent temperature [32]. As will be seen, these RF-accelerated FIs are not expected to destabilize mid-radius TAEs, so the follow-on nonlinear MEGA calculations in section 4only simulate the alphas. The study of RF drive, perhaps of energetic particle modes, will be pursued in future work and is further discussed in section 5. In both NOVA-K and the single-nversion of MEGA—as well as most, if not all, other AE stability codes, the toroidal mode number nis also an input value. Therefore, we are interested in modeling the ‘most unstable’ n. As in other works [7, 11,15,16,33], this is computed from the resonance condition, equating the mode and particle orbit widths. The mode width is approximated as ∆m≈rm m≈rm nq(rm),(1) where mis the poloidal mode number and q(rm)is the safety factor at the mode’s minor radial location rm. The orbit width is ∆o≈qv∥ ωFI =qmFIv∥ ZFIeB .(2) Here, the FI is characterized by its parallel velocity v∥, gyrofrequency ωFI, electric charge ZFIe, and mass mFI. From equations (1) and (2), the most resonant toroidal mode number is given by n∗≈ZFIe mFI Brm q2v∥ .(3) There are a few free choices here: First is the mode location, which we can estimate to be at mid radius, rm≈a/2. For TAEs, the primary resonance is at the Alfvén speed, v∥=vA. Plugging in values for alpha-driven TAEs in SPARC gives n∗≈8–18 for q=1−3/2. While this is a relatively wide range, it at least provides a starting point for simulation scans. 3. Linear modeling with NOVA-K 3.1. TAE existence and stability The eigenvalue solver NOVA-K [20–22] is first used to compute all possible AE mode structures and eigenfrequencies among the Alfvén continua. Profiles from section 2are provided as inputs. A fit to the magnetic equilibrium is performed internally in NOVA-K, with the magnetic axis, last closed flux surface, and qprofile serving as constraints (see figures 1(a) and 2). No toroidal rotation is considered here, but its expected effect would only be to add a Doppler shift to the frequency. A coarse scan in toroidal mode number is performed: n=5,10,15 and 20. This is chosen to cover the range of most unstable mode number n∗predicted in section 2; however, the upper bound is ultimately limited by the need for finer radial resolution with increasing nand resulting computational expense. Of all AE solutions, a subset is chosen based on (i) eigenfrequencies within the TAE gap and (ii) mode structures minimally intersecting with the Alfvén continuum. Representative low and high frequency AEs—i.e. even and odd modes— are selected for stability analysis, near the bottom and top of the TAE gap, respectively; their poloidal mode structures are shown in figure 4. Note that the modes are localized near the flux surface q=1, i.e. with dominant poloidal harmonics m,m+1≈n. Only the n=5 TAEs, with broadest mode widths, seem to have significant interactions with the Alfvén continuum. Other TAEs at higher q>1, i.e. √ψN>0.5 (see figure 1(a)), would intersect the continuum and likely have too high damping; thus, they are not assessed here. NOVA-K is then used to calculate the linear stability of each AE from a variety of drive and damping mechanisms. Table 1gives a breakdown of all contributions, with growth rates γ > 0 normalized to each eigenfrequency ω=2πf. (Note that γ/ω < 0 indicates damping.) As expected, the n=5 lowfrequency TAE exhibits large continuum damping compared to the n⩾10 modes; however, the high-frequency n=5 mode interestingly does not. Radiative damping [34,35] dominates for most modes and has a greater impact on the even vs 4 Nucl. Fusion 65 (2025) 036021 R.A. Tinguely et al Figure 4. NOVA: Poloidal mode structures (m, solid and dot-dashed) vs normalized radius for AEs with toroidal mode numbers n=5,10,15,and 20. For each n, low and high frequency modes (horizontal lines) are shown within TAE gaps of the Alfvén continua (thin lines). odd modes (i.e. at the bottom vs top of the TAE gaps), as discussed in [36]. Electron and ion Landau damping play a lesser role, with electron collisional damping the smallest. The net damping rate, without the contribution of FIs, is highest for the low-frequency TAEs, −γ/ω ∼1%–3%, while the highfrequency TAEs are only marginally damped, −γ/ω < 0.5%. 5 Nucl. Fusion 65 (2025) 036021 R.A. Tinguely et al Table 1. NOVA-K: Eigenfrequencies and normalized growth rates γ/ω [%] (damping <0) for AEs with toroidal mode numbers n=5,10,15,and 20. A breakdown of damping and drive mechanisms is provided, with and without finite Larmor radius (FLR) effects. Net growth rates are given with (bolded) and without (italicized) the contributions from fast ions (FI) including FLR effects: RF-accelerated He3 FIs and DT alphas. n5 10 15 20 f(kHz) 292 350 362 482 361 481 371 487 Continuum −0.06 −0.01 −0.02 −0.00 −0.02 −0.01 −0.00 −0.00 Radiative −2.69 −0.29 −1.35 −0.00 −0.95 −0.01 −2.63 −0.03 Electron collisional −0.01 −0.01 −0.02 −0.01 −0.01 −0.01 −0.01 −0.01 Electron Landau −0.08 −0.13 −0.01 −0.01 −0.01 −0.02 −0.00 −0.01 Ion Landau −0.02 −0.08 −0.08 −0.03 −0.05 −0.02 −0.08 −0.01 He3 FIs w/o FLR −0.18 −0.07 −0.20 −0.06 −0.19 −0.08 −0.74 −0.14 He3 FIs w/ FLR −0.17 −0.07 −0.19 −0.06 −0.17 −0.07 −0.35 −0.11 Alphas w/o FLR +0.06 +0.01 +0.75 +0.34 +1.01 +0.32 +1.85 +0.22 Alphas w/ FLR +0.06 +0.01 +0.61 +0.29 +0.79 +0.25 +0.71 +0.13 Total w/o FI −2.86 −0.44 −1.48 −0.05 −1.04 −0.06 −2.72 −0.06 Total w/ FI −2.98 −0.49 −1.05 +0.18 −0.42 +0.12 −2.35 −0.04 NOVA-K models the anisotropic distribution function of RF minority ions using the pressure profile in figure 3(a) and effective tail temperature in figure 3(c) [32]. (FLR effects can also be turned on and off in NOVA-K and will be discussed in the next section.) Interestingly, the He3 minority population is predicted to damp all modes. There are several possible explanations for this: First, the He3 FI pressure gradient is steepest in the plasma core, around √ψN∼ 0.1, whereas the TAEs are located at mid radius. Another reason could be that, given the effective parallel temperature O(75keV), very few He3 FIs will have parallel velocities comparable to the Alfvén speed. This situation could be significantly different for a H-minority species at planned lower-B0(8T) SPARC discharges; that is, H FIs could be ∼70% faster while the Alfvén speed would be ∼30% lower. Alphas are modeled in NOVA-K again using their pressure profile in figure 3(a), but now implementing a slowing down distribution in velocity space. Since the alpha birth velocity exceeds vA, the alpha population is calculated to drive all modes, as expected. In addition, the alpha pressure gradient is low, but fairly constant across the plasma (see figure 3(b)), including the mode locations. An important observation is that the alpha drive is stronger for the n⩾10 TAEs compared to n=5. This agrees with our earlier prediction of the most unstable mode number, from section 2. The total growth rate, including the contributions of both He3 FIs and alphas with FLR effects, is listed in table 1. Only the high-frequency n=10,15 TAEs are assessed to be marginally destabilized, γ/ω < +0.2%, even though the alpha drive is larger for the low-frequency TAEs. Odd TAEs, located at the top of the TAE gap, have been observed before in experiment [37]. However, if the contribution from radiative damping is overestimated by NOVA-K, we might expect low-frequency n=10,15 TAEs to also be driven unstable. 3.2. Uncertainties and FLR effects Regarding stability, it should be noted that absolute uncertainties ±0.1%are expected in NOVA-K’s calculation of continuum damping, larger than the values reported in table 1. In addition, as described in [38,39], the relative uncertainty in radiative damping is of order ∼10% when assuming 10% uncertainties in qand Tprofiles. Electron and ion Landau damping are exponentially sensitive [40] to the ratio of the Alfvén and thermal electron and ion velocities, respectively, making it more difficult to propagate uncertainties here; however, we note that their values in table 1are small compared to radiative damping and/or alpha drive for each mode. Finally, electron collisional damping depends on the electron beta and collision frequency (see equation (2) in [41]) such that O(10%) relative uncertainty may also be expected from the already low predicted value. Thus, even including the estimated uncertainties in damping, the overall stability would not change for each mode, with the exception of the highfrequency n=20 TAE, which is already near the marginal stability threshold. Furthermore, as discussed in section 2, the sawtooth instability will modify the location of the q=1 surface and AE parameters as a result. At the time of the sawtooth crash, any destabilized TAE near q=1 would likely disappear as q0∼qmin ∼1. Yet as the temperature increases in the core, energetic ion populations re-equilibrate, and the q=1 surface moves outward to mid radius, the TAEs identified by NOVAK could re-form (which has been observed in experiments before, such as in [42] and others). From figure 3, we might actually expect the TAEs to be more strongly driven in the core where the alpha and He3 FI betas and gradients are largest. The time evolution of the qprofile and its effect on AE stability will be explored in future work; interestingly, the present work may have identified the time when TAEs near q=1 are least driven by FIs. 6 Nucl. Fusion 65 (2025) 036021 R.A. Tinguely et al Lastly, we consider FLR effects. The maximum Larmor radii, at the approximate mode location √ψN∼r/a∼0.5, for 1MeV He3 FIs and 3.5MeV alphas are 1.2cm and 2.6cm, respectively, which correspond roughly to ∆r/a∼2.1%and 4.5%. The differences in drive and damping with and without FLR effects are shown in table 1. Very little effect is seen for He3 FIs’ damping of n=5–15 TAEs, while the damping decreases by ∼20%–50% for the n=20 TAE when including FLR effects. This is likely due to the mode width decreasing with n(see figure 4) and becoming more comparable to the Larmor radius, thereby enhancing FLR effects. In a similar way, FLR effects grow with nfor the alphas, although more significantly than for He3 FIs due to the alpha’s larger Larmor radius. For the n=20 TAE, a ∼40%–60% decrease in alpha drive is calculated when adding FLR effects. Thus, we would have overestimated the contributions from both He3 FIs and alphas by neglecting FLR effects in our stability calculations; this may have impacted the results of [11] as well. Even so, this would not have changed the conclusions about total growth rates for the n=5–20 TAEs, except the high frequency n=20 TAE, which would be just marginally unstable, instead of marginally stable. 4. Nonlinear modeling with MEGA 4.1. Numerical methods The MEGA code [23] is used to self-consistently solve the evolution of kinetic particles and the bulk plasma. The kinetic population modeled by Monte Carlo markers can be electrons [43], thermal ions [44], or FIs, as in the simulations discussed in this manuscript. In these runs, the kinetic population is simulated using the δfmethod, including FLR effects. Realistic injections of particles are available when using the full-fversion of the code, which is left for future work and will simultaneously include both alphas and RF-accelerated FIs. The bulk plasma is described by the non-linear, full MHD equations, including diamagnetic drift and toroidal flows. MEGA assumes quasineutrality on the MHD grid, providing a single density, velocity, and pressure for both thermal ions and electrons. The coupling between the MHD grid and the kinetic species is included through the energetic particle current density term in the MHD momentum equation. More details on the MHD module of the MEGA code can be found in [45–47]. MEGA has been extensively validated against experimental data. Some examples include reproducing the observation of ‘abrupt large-amplitude events’ in JT60-U [48], the visualization of the AE-induced FI flow using imaging neutral particle analyzers [49–51], the impact of the bump-on-tail of the NBI FI distribution on TAE growth rate [52], the spectrum of externally applied perturbations on AE stability [53,54], as well as measurements of stable and unstable TAEs using the Alfvén Eigenmode Active Diagnostic in JET [38]. The simulations described here are single-n, in which only a portion of the toroidal geometry of the tokamak (ϕ∈ [0,2π/n]) is simulated, including periodic boundary conditions. While these simulations do not capture the interaction between modes of different toroidicities, they are computationally efficient, as only Nϕ=32 toroidal grid points are required to resolve the instability for n⩽16. The poloidal resolution, NR=NZ=256, is set to resolve n=10 perturbations, using the radial profile of the TAEs calculated by NOVAK as a mock-up and comparing against the grid distributions on the poloidal plane (see figures 2and 4). In these simulations, the frequency of the destabilized modes will depend on the background mass density. As these DT plasmas are expected to have a non-negligible impurity (see figure 1(a)) and nHe3/ne=5%concentration, a correction of 1.15 is applied to the simulated electron density to account for the ratio of the effective density to the bulk deuterium simulated by MEGA. While the typical on-axis electron temperature for AUG and DIII-D plasmas is about 5keV [53], the on-axis temperature for this SPARC case is about 20keV (see figure 1(b)), resulting in a much lower Spitzer resistivity ηS. Therefore, the simulated value of normalized resistivity is reduced to ηMEGA/(vAR0µ0) = 5×10−8, which is 10 ×smaller when compared to similar simulations of AUG and DIII-D plasmas [47,53]. Yet this ensures that the same ratio of the simulated resistivity with respect to the Spitzer resistivity is maintained, about ηMEGA/ηS∼2000. This reduction is admissible as the poloidal resolution is increased compared to the AUG and DIII-D simulations, where numerical convergence has been already studied [54]. The values of viscosity and diffusivity are maintained at ν=χ=5×10−7vAR0, which is the same value reported in [46]. The kinetic population is an isotropic slowing down distribution resulting from alphas generated at 3.5MeV. The different mass and charge with respect to the bulk plasma is carefully considered when normalizing the inputs. The simulated distribution is isotropic in pitch angle and has the sloweddown energy profile depicted in figure 5(b), which is determined by the space-dependent critical velocity calculated by MEGA based on the input temperature and density profiles (see figure 1). The radial profile of the alpha pressure (βα) is determined by the quasi-analytic built-in distribution [55], βα(ψN) = βα0exp[−(ψN/∆ψ)a],(4) where the pressure on axis (βα0=0.65%), radial gradient scale length (∆ψ=0.16), and exponential factor (a=1.00) are adjusted to resemble the TRANSP/NUBEAM outputs, as depicted in figure 5(a). The on-axis value of normalized alpha pressure (βα0) is nearly half of the typical values used for ITER simulations; for instance, βα0≈1.2%was used in [45]. This small value of the normalized pressure, despite the >100MW of fusion power, is explained by the high magnetic field strength B0=12.2T. 7 Nucl. Fusion 65 (2025) 036021 R.A. Tinguely et al Figure 5. MEGA: (a) Radial profile of the alpha normalized pressure βαsimulated in MEGA, which is adjusted to match the NUBEAM distribution (see figure 3(a)), as a function of normalized radius ρpol =√ψN. (b) Energy distribution of the alpha particle population. Figure 6. MEGA: (a) Temporal evolution of the n=10 radial velocity perturbation, δvr, at the midplane during the linear phase. The mode amplitude is normalized by exp(−γtk)so that the mode location can be observed throughout the entire linear phase. (b) Fast Fourier Transform of the fields with the Shear Alfvén Wave continuum overplotted, showing the mode near the bottom of the TAE gap. 4.2. Marginally unstable n =10 TAE The single-nsimulations described in the previous section are run for longer than t=0.22ms =163τA0, with τA0= 2πR0/vA0the Alfvénic period. The increased poloidal resolution requires increased computational resources, so each simulation takes an average of 72h using 512 CPUs. The evolution of the n=10 radial velocity perturbation δvr, once the mode is destabilized, is plotted in figure 6(a). In this figure, the amplitude at each time step tkis multiplied by a factor of exp(−γtk), with γan ad-hoc parameter adjusted to visualize the evolution of the perturbation; the mode is be observed to oscillate around ρpol =0.38 for the entire linear phase. Figure 6(b) depicts the Fast Fourier Transform (FFT) of the δvrfield, together with the Shear Alfvén Wave (SAW) continuum calculated using the ALCON code [56]. As the linear growing phase is relatively long, the FFT has a good frequency resolution, and we observe the mode located near f=370kHz, exactly at the location where the m=9 and m=10 even coupling is located [37]. This even coupling of m=9,10 is even more evident when looking at the poloidal harmonics during the linear phase (figure 7), which allows us to classify this 8 Nucl. Fusion 65 (2025) 036021 R.A. Tinguely et al [51] Du X., Heidbrink W., Van Zeeland M., Gonzalez-Martin J., Austin M., Yan Z. and McKee G. 2023 Visualization of fast ion phase-space flow in plasmas well-below, near and well-above Alfvén eigenmode stability threshold in tokamak Nucl. Fusion 63 046020 [52] Van Zeeland M.A. et al 2021 Beam modulation and bump-on-tail effects on Alfvén eigenmode stability in DIII-D Nucl. Fusion 61 066028 [53] Gonzalez-Martin J. 2023 Active control of Alfvén eigenmodes by externally applied 3D magnetic perturbations Phys. Rev. Lett. 130 035101 [54] Gonzalez-Martin J. et al 2024 Active control of Alfvén eigenmodes by external magnetic perturbations with different spatial spectra Nucl. Fusion 64 076022 [55] Wang H. and Todo Y. 2013 Linear properties of energetic particle driven geodesic acoustic mode Phys. Plasmas 20 012506 [56] Deng W., Lin Z., Holod I., Wang Z., Xiao Y. and Zhang H. 2012 Linear properties of reversed shear Alfvén eigenmodes in the DIII-D tokamak Nucl. Fusion 52 043006 [57] Pinches S.D., Chapman I.T., Lauber P.W., Oliver H.J.C., Sharapov S.E., Shinohara K. and Tani K. 2015 Energetic ions in ITER plasmas Phys. Plasmas 22 021807 [58] Wallace G., Migliore C., Wright J., Brookman M. and Garrett M. 2023 Retiring risk for ion cyclotron range of frequency heating in SPARC through modeling APS Division of Plasma Physics Meeting Abstracts (Denver, Colorado,30 October 3–November 2023) vol 2023 p NO05–008 (available at: https://meetings.aps.org/Meeting/ DPP23/Session/NO05.8) [59] Tinguely R.A., Puglia P.G., Dowson S., Porkolab M., Douai D., Fasoli A., Frassinetti L., King D. and Schneider P. (JET Contributors) 2024 Isotope effects and Alfvén eigenmode stability in JET H, D, T, DT and He plasmas Nucl. Fusion 64 096002 [60] Sorbom B.N. et al 2015 ARC: a compact, high-field, fusion nuclear science facility and demonstration power plant with demountable magnets Fusion Eng. Des. 100 378–405 [61] Kuang A.Q. et al 2018 Conceptual design study for heat exhaust management in the ARC fusion pilot plant Fusion Eng. Des. 137 221–42 [62] Creely A. et al 2022 Demonstration of fusion pilot plant physics in SPARC APS Division of Plasma Physics Meeting Abstracts (Spokane, Washington,17–21 October 2022) vol 2022 p NO03–003 (available at: https://meetings.aps.org/ Meeting/DPP22/Session/NO03.3) [63] Hillesheim J. et al 2023 ARC physics basis status APS Division of Plasma Physics Meeting Abstracts (Denver, Colorado,30 October 3–November 2023) vol 2023 pp J11–115 (available at: https://meetings.aps.org/Meeting/ DPP23/Session/NO05.8) [64] Fu G.Y. and Van Dam J.W. 1989 Excitation of the toroidicity-induced shear Alfvén eigenmode by fusion alpha particles in an ignited tokamak Phys. Fluids B11949 [65] Cheng C.Z. 1991 Alpha particle destabilization of the toroidicity-induced Alfvén eigenmodes Phys. Fluids B 32463–71 [66] Breizman B.N. and Sharapov S.E. 1995 Energetic particle drive for toroidicity-induced Alfven eigenmodes and kinetic toroidicity-induced Alfven eigenmodes in a low-shear tokamak Plasma Phys. Control. Fusion 37 1057 15