Fast ion relaxation in ITER mediated by Alfvén instabilities

We address the critical issue for future burning plasmas of whether high-energy fusion products or auxiliary heating-beam ions will be confined for a sufficiently long time to compensate for thermal plasma energy losses. This issue can be mitigated by one of the most deleterious collective phenomena—the instability of low, sub-cyclotron frequency Alfvén eigenmodes (AEs), such as toroidicity-induced AEs and reversed-shear AEs in the ITER steady-state scenario. Using a revised quasi-linear (QL) theory applied to energetic particle (EP) relaxation in the presence of AEs, we find that the AE instabilities can affect both neutral beam ions and alpha particles, although the resulting fast ion transport is expected to be modest if classical particle slowing down is assumed. On the other hand, the QL theory predicts that the AE amplitudes will be enhanced by the background microturbulence, although this topic remains outside our scope due to the significant numerical effort required to evaluate these effects. We report our results for EP relaxation dynamics obtained utilizing several tools: (i) a comprehensive linear stability study of the sub-cyclotron Alfvénic spectrum as computed by ideal magnetohydrodynamic NOVA simulations for the AE eigenproblem, (ii) drift kinetic NOVA-C calculations for wave–particle interaction and AE growth/damping rates, and (iii) predictive QL modeling coupled with the global transport code TRANSP to assess the EP relaxation on the equilibrium timescale.


Introduction
The problem of energetic particle (EP, also referred to here as energetic or fast ion) confinement in tokamaks is critical to achieving a successful self-sustained controlled thermonuclear reactor.For example, collective effects due to Alfvén eigenmodes (AEs) in ITER are expected to play a significant 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.
role in the relaxation of 3.52 MeV fusion products, alpha particles, and 1 MeV injected deuterium beam ions.
It has been recently noted that regimes of enhanced-AEinduced fast ion transport can occur, where the background microturbulence mediates this EP relaxation [1].This highlights a novel way to explain fast ion relaxation losses that goes beyond the scenarios described in [2].It has been numerically demonstrated that microturbulence-driven EP pitch angle scattering can increase the AE amplitude to a level significantly greater than the level likely to damage the ITER first wall [3].
Different reduced models can be used to simulate multimode EP relaxation under burning plasma conditions.For example, the critical gradient model [4,5], which does not include velocity space resolution, can nevertheless be used for the fast evaluation of EP relaxation.The main problem with this model is that it requires validations which are not available in future devices, such as ITER, nor can it be generalized for varying plasma regimes in present-day devices.Another model is a more detailed resonance broadened quasi-linear (RBQ) model that evaluates EP relaxation in both energy space and canonical toroidal momentum space in the presence of Fokker-Planck collisions and can include the pitch angle effective scattering and zonal flows caused by microturbulence [6].RBQ is general enough and ready for applications in planned future fusion devices.Both models rely on AE linear instability growth rate calculations, e.g.those provided by NOVA/NOVA-C [7], whereas RBQ can implement effective EP pitch angle scattering, if available, which in the case of classical coulomb collisions comes from NOVA-C simulations but needs to be enhanced by the contribution due to the background microturbulence [8].
It is important to highlight that RBQ simulations consider the wave-particle nonlinearity arising from fast ion interaction with AEs and the resulting alteration of the EP distribution function in the constants of motion space.
In ITER, super Alfvénic ions include the auxiliary deuterium heating-beam ions and fusion-born alpha particles needed in the plasma to replenish thermal ion losses.The alpha particle distribution function is far from that of a typical beam or ion cyclotron resonance heating ions which normally have narrow widths in pitch angles, χ = v ∥ /v, where v ∥ is the component of the ion velocity parallel to the magnetic field and v is its absolute value.For example, tangentially injected beam ions in ITER are predicted to have δχ ≲ 0.1 [9].Comprehensive linear and nonlinear studies of AE stability were recently performed for the ITER baseline scenario [10], where several important kinetic damping mechanisms were accounted for and their effects on EP relaxation were discussed in detail.These sources of damping should be considered standard for reliable predictions of EP confinement in fusion-grade reactor devices.One of the damping mechanisms mentioned in this reference is the trapped electron collisional damping, which sets boundary conditions at the plasma edge but is often ignored in benchmarks [11].
We note that an essential feature of EPs in fusion reactors that include fusion charged products is their bootstrap current, which can provide the additional current drive needed in tokamak reactors, especially in spherical tokamaks (STs) [12].The fusion product bootstrap current has a nonzero finite value at the plasma center which is important for STs.The fast ion current drive is extremely sensitive to the details of the energetic ion distribution function in the constant of motion (COM) space, requiring accurate modeling of AE-driven diffusive and convective transport [13].These issues are vital for making credible projections for fusion reactors.

Toroidicity-induced AE (TAE)/reversed-shear AE (RSAE) linear stability in the ITER steady-state scenario
We start with a comprehensive linear stability analysis performed with the help of the ideal magnetohydrodynamic (MHD) code NOVA [14] and its hybrid drift kinetic extension NOVA-C4 [7].NOVA ideal MHD analysis uncovers the Alfvénic eigenmodes in the frequency range spanning from the geodesic acoustic-mode frequency up to the ellipticityinduced Alfvénic eigenmode gap.In figure 1 we present the ITER steady-state plasma profiles prepared by the ASTRA code [15] and made available for simulations through the IMAS framework [16].These were summarized in an IAEA presentation [17].In the selected ITER steady-state scenario, the only planned auxiliary heating and current drive schemes are neutral beam injection (NBI) and electron cyclotron heating (ECH) [15].
ITER steady-state operation is characterized by the deuterium NBI power P NBI = 33 MW and the ECH and current drive power P EC = 20 MW.Unlike the baseline case scenario considered earlier [10], this AE stability analysis addresses an ITER steady-state scenario characterized by a fusion αparticle pressure profile twice as large and a beam ion beta profile ten times larger.Developing plasma control in the steadystate scenario of interest involves the use of various techniques to optimize the NBI current drive in particular [15].Unsurprisingly, the NBI contribution to AE instabilities was ignored in the follow-up studies [18,19].We stress here that the NBI current drive can be quite an important factor for ITER operation even though we show here that the effect of AEs is benign for EP ion confinement.
Both fast ion species, i.e. fusion α particles and beam ions, are included in our simulations.Their distribution contour maps are shown in figure 2 and are taken near the injection (birth) energy in the plane of canonical angular momentum, P φ , and normalized to the EP energy adiabatic moment, λ = µB 0 /E.Given the ratio of their super-Alfvénic velocities to the Alfvén speed, it is expected that the AE instability drive is maximized relative to the v α,b0 /v A ratio.Furthermore, one can see jumps at the separatrix between the passing and trapped alpha particles in the right figure.These jumps are physical, and are due to jumps in the drift orbit precession times associated with the transitions from trapped to passing ion orbits and vice versa.If multiplied by the precession times of the fast ions, their 'distribution function times the precession time' quantity becomes a smooth function near the transition points and does not have those jumps.
We note that figures 2(a) and (b) are plotted on the same plane as that of figure 19 'Confinement loss domains in µ/E, P φ ' of [20], where various particle orbits are labeled.For instance, it is evident that the fast ions injected by NBI are exclusively composed of co-passing particles.
NOVA has found around 600 Alfvén modes of interest, out of which NOVA-C identified 42 unstable or marginally stable eigenmodes for the subsequent RBQ runs.The AE stability is addressed by the kinetic NOVA-C code, which incorporates rich physics including background damping and advanced fast particle representation that achieves favorable benchmarks Plasma profiles of ITER steady-state scenario used in simulations, showing the radial profiles of the beam pressure, P b ; the fusion product alpha particle pressure, Pα; the thermal ion pressure P i plotted in Pascals; the electron density, ne in units of cm −3 ; the safety multiplied by 1/3, i.e. q/3; and the magnetic shear s = q ′ r/q.The horizontal axis corresponds to the square root of the toroidal flux normalized to its value at the last closed surface.compared to the main available stability codes [11].We show the AE growth rates in figure 3.
Based on the linear AE stability results, several important observations can be made.First, our linear AE stability analysis is consistent with earlier results [18,19], especially the toroidal mode number range of the unstable AEs and their characteristic growth rates.The growth rates have a maximum or a rollover point at around n = 20 ÷ 30 in those calculations.Second, two regions in the radius have the strongest contributions to the growth rates: near the plasma center and near the q min location.Third, the multiplicity of AEs and their resonances likely results in overlaps of the resonances in the COM space due to the nonlinear amplitude broadening.Fourth, the dominant damping rate that controls the EP profiles near the edge, i.e. the rate that sets the boundary conditions, is the trapped electron collisional damping.For reliable predictions, such damping must be included in simulations.Ignoring it would change the radial domain of EP confinement and, as a result, would overestimate the loss fraction of fast ions as discussed in [21], where the critical gradient model (CGM) [5] was revisited.We would like to note that the abovementioned CGM approach has limitations.One of them is that it does not predict the AE amplitudes near saturation.It also does not allow for the slowing down of fast ions nor for EP pitch angle scattering over time.The CGM can be considered to a limiting case of the QL approach if the number of unstable AEs goes to infinity.This was shown in [4].We include most of the linear dampings relevant to AE stability, including thermal ion and electron Landau dampings, radiative damping, as well as the nonlinear (or QL) damping that emerges during the AE amplitude saturation.One of the dampings neglected in our calculations is the continuum damping.This approximation is not expected to change our conclusions concerning the fast ion relaxation, since the most unstable AEs have a high-n toroidal mode number whose continuum damping is negligible [22].We discuss this approximation in detail in appendix.

Quasi-linear (QL) modeling of fast particle relaxation
The recently built numerically efficient, self-consistent QL code RBQ is applied to the ITER steady-state scenario [15].It is capable of addressing EP confinement in the presence of several or many Alfvénic modes by following the time evolution of their amplitudes and by advancing the EP distribution function in the space of the invariants of the unperturbed motion.This can be done either within the RBQ code itself or by providing the diffusion and convective coefficients for subsequent calculation by a whole-device modeling (WDM) package such as NUBEAM [23], which has rich physics.Furthermore, the RBQ code has been interfaced with the WDM TRANSP code and has been preliminarily shown to agree well with experimental measurements when applied to DIII-D critical gradient experiments [6].In the NUBEAM guiding-center orbit-following simulations, RBQ transport coefficients are applied while maintaining the same wave-particle nonlinearity utilized in the RBQ model.
The RBQ model solves a system of equations (published earlier for 2D, i.e. the R 2 case [6]) by assuming the ansatz for the harmonics of the perturbed quantities e −iωt−inφ+imθ .For a single mode of a tokamak, the QL diffusion of an ion occurs in the R 1 slanted direction near each resonance along the paths of constant values of the expression [24]: where ω and n are the angular frequency and the toroidal mode number of a mode, respectively.An immediate consequence of this expression is that at positive values of n, increasing P φ (the EP diffusion to the plasma edge) leads to a decrease in EP energy.In other words, this is how the unstable system transfers energy to the eigenmode.Following [25], each analyzed mode is denoted by the symbol k and corresponding resonances by the symbol l.The latter identifies the resonance diffusion region using the index of the poloidal harmonic and the index of the poloidal sidebands, hence it is a vector.
If the ensemble of modes is known, the distribution function evolves according to the following QL equations implemented in RBQ.The system of QL equations was given in its general form in [26].It was adapted for the NOVA-C notation in [25,27], which was generalized for the 2D case as:  where the diffusion coefficients are: and where ν ⊥ is the 90 0 pitch angle scattering rate frequency, f and f 0 are the distribution functions at times t and t = 0 respectively, F l is the resonance window function, G kmp are the wave-particle interaction (WPI) matrices, and Ω l is the resonance frequency of the WPI between the mode k and the resonant ion.In an earlier 1D version of RBQ, we used the action variable e. a toroidal canonical momentum divided by −n k , ignoring energy dependence [6,25].Also, here I kr is the action at a given resonance.The subindex k denotes action related to the kth mode.In general, three action variables can be used within the RBQ framework, i.e.L operates in the 2D plane, L : (P φ , E) ∈ R 2 at each value of µ.In the current version of RBQ2D, we consider the sub-cyclotron frequency range, implying that µ = const for the operator L. Also, unlike the recent revision of the QL theory [28], equation (2) includes only the diffusion terms and ignores the convective velocities which require the knowledge of zonal flows (ZFs).ZF modifications of the RB model are straightforward if the convective velocities are known but are beyond the scope of this paper.The equation for the amplitude of each mode can be formally written without the explicit contributions from other modes, though such contributions are mediated by the fast ion distribution function.The amplitudes of the modes of interest evolve according to the set of equations: where γ L,k and γ d,k are the linear growth and damping rates of the kth mode at time t, and C k = δB θ,k /B is the amplitude of that mode.The AE damping rate is considered to be constant in time but it can, in principle, change on a slow timescale during RBQ evolution.The RBQ model is designed for applications that use both isolated and overlapping modes, whereas the conventional quasilinear theory applies to multiple overlapping modes.We use the same structure of quasilinear equations for the ion distribution function evolution [26,29,30] but with the resonance delta function broadened over some region in the direction of the action I k variation: The resonance broadening function F l is the key novel element of the proposed diffusion model, whose parametric dependencies are verifiable against known analytic asymptotic behaviors for isolated modes (see the appendices of [6]).It satisfies the property in [31].
If the diffusion occurs approximately along P φ , i.e. in the low mode frequency approximation, the whole problem reduces to the set of 1D equations discussed in Fitzpatrick's thesis [29] (in addition, see [32] where the window function was introduced in its preliminary form; it was revisited in [31]).The effect of several low-n modes on EPs cannot be accurately determined without numerical simulations that capture the diffusion in 2D space, (P φ , E) ∈ R 2 .This is due to the complex multidirectional diffusion dynamics of the resonant ions in R 2 space and possible resonance overlap in the presence of multiple AEs.Thus, this 2D generalization is a major extension target for the RBQ model.
We applied the RBQ to an ITER steady-state plasma characterized by 42 unstable or marginally stable AEs and prepared the COM diffusion coefficients for WDM processing.Both beam ion and alpha particle plasma components contribute to the AE drive, although the linear stability analysis indicates that, on average, the growth rate of beam ions is twice as large.This is because beam ions are injected into the most unstable location in the COM space; almost all the beam ions are passing ions and resonate easily with AEs.
The RBQ precomputed diffusion coefficients turn out to be quite large and localized where most of the 42 tested AE modes are unstable or marginally unstable, i.e. near the q min point, reaching values of up to 50 m 2 s −1 in the COM space and going down to ∼ 1 m 2 s −1 at the periphery.This implies a local flattening of EP profiles near q min , but the overall EP redistribution and losses in particular can be modest.
Figure 4 shows the evolution of EP relaxation over time corresponding to AEs driven purely by injected beam ions; in other words, this is the case of a pure beam ion drive unassisted by fusion alpha particles.We can see that the amplitudes of the AE modes are quite small, reaching values of δB θ /B ∼ 10 −5 ÷ 10 −3 .As a result, no significant beam ion losses to the wall are present in NUBEAM simulations.At the same time, the analogous, unassisted AE driven by fusion alpha particles shows very weak amplitude modes at δB θ /B < 10 −4 .
In contrast to unassisted AE drives, let us consider a case more relevant to experimental expectations, namely AE instabilities in which both fusion alpha particles and beam ions contribute to the drive.The evolution of the AE amplitude over time is shown for this case in figure 5. We follow the prescription given by QL theory [26], which controls DF evolution according to the equations written for that species.It follows from these figures that the slowing down beam ions together with fusion alpha particles drive AEs to higher amplitudes in comparison to the overshoot case, up to δB θ /B ≲ 3 × 10 −3 .Using AE amplitudes prepared by the RBQ code, the diffusion coefficients were transferred to the NUBEAM Monte Carlo package, whose calculations are discussed in section 5.It is important to note that the initial saturation level we see in RBQ simulations at t = 1.5 ÷ 3 ms (see figure 4) corresponds to the case in which the EP profiles are still far from the relaxed state and thus the AE amplitudes are still converging to the saturated state.It is also interesting that both cases represented in figures 5(a) and (b) have approximately the AE amplitudes at saturation; that is, they are determined primarily by the values of the AE growth rates, which are the same in both cases.

NUBEAM gyro-center simulations of fast beam ions and the slowing down of fusion alpha particles
Ideally, the EP relaxation used in WDM simulations should include the evolution of the EP distribution function in COM space over a timescale that includes three characteristic times: the inverse linear AE growth rate, the damping rate, and the effective pitch angle scattering rate [33].Such dynamics have been demonstrated and are well captured by the QL theory implemented in RBQ simulations [1].Bearing this approach in mind in the longer term, we aim at steady-state plasma operation regimes but ignore such intermittent behavior.
This approach can unify different simulations of several initial value codes, such as those used in a recent comprehensive benchmark [11], by making use of the WDM prototype NUBEAM package [23].It requires the following diffusion and convective transport coefficients to be recorded in simulations [34] either by a reduced model, such as the RBQ model, or by the initial value codes discussed in [11]: where D xȳ is the diffusion coefficient and C x is the convective motion coefficient of fast ions due to either the AE activity, collisions, or the associated ZF.Here, the overbar symbol ( ¯) means that the value is normalized to either ψ θ w (in the case of canonical momentum) or to either the EP birth or the injection energy E 0 (in the case of EP kinetic energy).
We note that the coefficients D, C need a special rule for implementation in NUBEAM to describe the resonant ion motion in the COM space, since COM coordinates are specific for each species under consideration.They are given by a well-known expression if one AE is the only mediator, which follows from equation (1): keeping µ = const.The minus sign here implies that the resonant particle loses its energy when Pφ increases, i.e. the resonant ion moves toward the edge and loses its energy to drive the instability.The WPI is the only case relevant to fusion plasmas when other scattering mechanisms are negligible, such as when the turbulence-induced scattering is small and the AE/EP system goes into the frequency chirping regime [35].For Alfvénic oscillations, the eigenmode frequency is typically small, much smaller than the cyclotron frequency of plasma ions, so that it is quite safe to assume ∆ Ē/ Ē ≪ |∆ Pφ / Pφ | for analytic estimates, if needed.The above condition, equation ( 8), needs to be enforced for EP relaxation.One way to do this is relevant to cases in which the bins required to describe EP distribution function in the COM space are large, such as those extensively used in the kick model [36].In this case, a special probability density function (PDF) has to be introduced to differentiate resonant ions from non-resonant ones.This scheme adds dimensionality and complexity to the problem.
Another route exists if the dimensions of the problem are sufficient to resolve EP dynamics.Let us consider the COM grid bins required for NUBEAM simulations; these contain N points with Pφ and Ē grids that are equidistant, as in the NOVA-C/RBQ model.The dimensions fixed in our problem are N Pφ = 200 and N E = 40.We then construct the local bin PDF characterized by the COM location, as follows: where 0 < ρ < 1, so that at ρ = 0 the resonant ion experiences a 'kick' from cell i, j with equal probabilities of moving in the ( Pφ − Pφ0 ) /σP φ or Ē − Ē0 /σ Ē directions over a time step ∆t.The PDF function needs to be normalized: An example of such a PDF is graphically given in figure 6.If ρ = 0.9, the natural for the WPI slanted direction is enforced as shown in equation ( 8) which is close to (and it is exact at ρ = 1).In the near-threshold case, i.e. when the diffusion is primarily in the Pφ direction, it needs to be corrected by an additional Ē diffusion by the amount given in equation ( 8).
Thus, given the time step ∆t, the resonant ion is expected to have a diffusive change in the P φ direction given by 11) and a change in Ē: where the signs σ 1 = ±1 and σ 2 = ±1 are independent of each other and determined randomly.Another change it experiences is due to the convective motion in the P φ direction, given by CP φ ∆t, and in the E direction, given by C Ē ∆t.The random steps in the COM variables of NUBEAM are natural for the Monte Carlo technique.
In this paper, we follow the third method to implement the transport coefficients D, C in NUBEAM calculations, as described below.Since the EP diffusion and convection are found for each cell and fast ions diffuse (and convectively move) according to the precomputed values of D, C, each ion at a grid point (grid cell or narrow bin) Pφ , Ē diffuses with the coefficient D and at the same time experiences convection C.
We now prescribe the following changes in the COM variables for the particle kicks: This complicated recipe guarantees that the slanted path is used in COM EP diffusion, equation (8).It may not always be possible to distinguish EP diffusion from convection in initial value codes such as GTC or GYRO.In these cases, it seems reasonable to prescribe EP motion using convection, but more often, an exchange of EP convection will be required.The representation of the EP motion during the transport processes in terms of its diffusive and convective kicks can be considered to be a simplification if the characteristic timescale and the spacial grids are large.However, with sufficient grid resolution, the above representation, equations ( 6) and ( 7), seems to be adequate for the problems of interest.Another simplification of this problem is that the coefficients D, C are the sums of the contributions from several resonances, since AEs have global mode structures.
The optimization of this model with regard to the time step ∆t can also be done, for example, in the following way.The main limitation of Monte Carlo simulations arises from the requirement that the diffusion/convection steps should not be larger than the grid size ∆ Ē < N E , ∆ Pφ < N Pφ .Moreover, since the ion motion is primarily expected to take place in the P φ direction for the Alfvénic oscillations and associated WPI, one would expect that the choice of timestep can be made according to the requirement for the radial kick.Thus, we can write it as

Beam ions and fusion alpha particles are depleted in the plasma core by AEs
We apply the NUBEAM Monte Carlo guiding center-orbitfollowing package to compute the confined alpha particle and beam ion density profiles, which are typically monotonic in radius, as shown by the solid curves in figure 7. When the AE diffusion prepared by the RBQ model is turned on, NUBEAM finds that the EP profiles of both alpha particles and beam ions become depleted near the plasma center, as shown by the dashed lines in the graph.Thus far, this is perhaps the most important implication of the effect of the Alfvénic modes on fast ion confinement in ITER tokamak plasmas.According to our calculations, the AE-driven losses are expected to be negligible for beam ions, whereas for alpha particles, they are found to be around 1.7%.This is likely due to the significant fraction of trapped alpha particles in comparison with beam ions which are primarily passing; see a figure 2. Trapped alpha particles deviate by 2qρ Lh / r/R from the magnetic surface, whereas passing ions stay closer, deviating by qρ Lh in the radial direction [24] where ρ Lh denotes the fast ion Larmor radius.
Another important factor missing in our calculations in comparison with our own results applied to DIII-D [6] is the absence of robust evaluations of the additional, effective pitch angle scattering which needs to be included using dedicated effort for ITER.This seems to be a concern for future applications of many reduced simulations which need to eventually rely on the values of thermal electron and ion conductivities.Our study is based on TRANSP calculations which may be far from the expectations for ITER experiments.This concern is a further motivation to experimentally and theoretically study the effect of microturbulence on MeV ions in ITER.
Nevertheless, we will use what exists in the literature, and in particular in [37], where the required effective scattering frequency of fast ions is based on simulations performed using the GS2, GYRO, and GKW codes.If we follow this route, we make use of the thermal ion and electron conductivities found by TRANSP, which are χ i = 4.2 * 10 3 cm 2 s −1 and χ e = 4.9 * 10 3 cm 2 s −1 .Also, TRANSP has found that the thermal electron temperature T e = 27.6 keV at the q min location, which is r/a = 0.345 in the run of interest.We find that the effective pitch angle scattering rate is 0.026 ms −1 for beam ions and 0.013 ms −1 for alpha particles.Given this uncertainty, we fix it uniformly in radius to 0.03 ms −1 for both species for long-term RBQ simulations corresponding to ∆t = 100 ms, which seems to be sufficient to address the steady-state diffusion due to AE/EP interactions.Calculations using the RBQ model are still expensive and have to be run on a Princeton Plasma Physics Laboratory (PPPL) computer cluster for several days, whereas NUBEAM with the RBQ diffusion rates requires an order of magnitude less time, since no updates of the AE growth rates are required.
We have found that the fast ion slowing down relaxation using only classical, coulomb scattering leads to insignificant density flattening near the q min location.However, the most important effect is in the EP density depletion near the plasma center.This is true for both beam ion and alpha particle species that have comparable injection or birth energies per nuclei that are about twice as large as the energy corresponding to the Alfvén speed.
As previously mentioned, RBQ simulations find that saturated AEs with modest amplitudes, δB θ /B ∼ 10 −5 ÷ 10 −3 , account for strong local radial diffusion, which reaches up to ∼ 50 m 2 s −1 for resonant fast ions near q min .However, near the plasma center the EP density is depleted, likely due to the finite orbit width effect, with the result that it is comparable to the minor radius value in the central region.A similar depletion was reported in DIII-D beam ion profiles measured and simulated by the kick model [38].A preliminary interpretation of this feature ponts to finite orbit width effects, which are much larger in DIII-D than in ITER.This effect, coupled with strongly driven AE diffusion of several radially overlapping modes near the center and relatively weak coulomb pitch angle scattering, results in the EP profile inversion obtained by the NUBEAM package.
We see a somewhat similar effect in ITER simulations.To illustrate this conjecture, we write the WPI resonant condition for fast MeV energy ions that we simulate with the NOVA-C code [7].The WPI resonance condition reads: where v dr is the EP drift velocity and l is an integer.In conventional tokamak experiments on present-day devices, the dominant resonance contributions are due to two resonances, v ∥ = v A /3 and v ∥ = v A , corresponding to l = ±1 [7,39,40].This case corresponds to k ∥ v ∥ ≫ |k ⊥ v dr |, which is typical for DIII-D [38].When this inequality is reversed, the contributions of the sidebands to the AE growth rate become dominant.This is also the case for the modes we find in ITER simulations; see figures 5, n = 15 ÷ 20 in the saturated state.In other words, v α,b0 l/qR ≃ qnρ Lh v α,b0 /rR near the plasma center.This implies that r/a ≃ q 2 /l, which is true for higher-order sidebands expected at the plasma center, l > 1.
The consequences of the central region density depletion are not severe for the plasma power balance, since no significant losses are found.What may be more important for the plasma discharge from the point of view of AE excitation is their effect on the current drive.Indeed, the co-injection of beam ions plays an important role in creating and maintaining the reversed safety factor profile.As a result, beam ion depletion at the plasma center can pose severe limitations in the steady-state scenario.

Summary
We performed a comprehensive stability analysis of the ITER steady-state plasma using the ideal MHD code NOVA, its drift kinetic extension NOVA-C, and the recently developed 2D QL code RBQ, which employs a novel and revised methodology.This analysis helps to evaluate the AE saturation amplitudes and the relaxation dynamics of EPs present in planned ITER operations when both EPs, i.e. super-thermal beam ions and fusion α particles are included.
Calculations using NOVA and NOVA-C revealed 42 unstable and marginally stable AEs that were used for the subsequent RBQ processing.The modes thus identified included RSAEs, TAEs, and ellipticity-induce Alfven eigenmodes (EAEs) residing in the corresponding gaps of the Alfvén continuum.The results of the linear stability analysis showed that the most unstable AEs' toroidal numbers span from n = 1 to n = 40 with relatively modest growth rates, γ L /ω < 5%, which justifies the application of the perturbative analysis.These results also justify the subsequent QL analysis, since most AE growth rates are not strong enough to change the EP relaxation and the effect of ZF is not expected to change it [41].
The application of RBQ in its 2D version showed that the AE amplitudes remain relatively low, δB θ /B < 3 × 10 −3 , for all 42 analyzed modes.In its current version, RBQ provided the AE diffusion and convective coefficients for the subsequent analysis.In the final stage of our analysis, the NUBEAM package was applied to evolve the fast ion distribution functions more accurately in the COM space, which means that NUBEAM evolved EPs on a slowing-down timescale.This required NUBEAM to account for EP slowing down under realistic ITER conditions, that is to follow EP slowing down for several seconds.Both RBQ distribution evolution and NUBEAM calculations found that the EP confinement with coulomb scattering does not result in significant fast ion losses, which is consistent with earlier studies of the baseline calculations [18,19].However, we have previously found that background microturbulence can significantly boost the AE saturation amplitudes by broadening the phase-space locations near the WPI resonances [1,6].For example, the AE amplitudes can go up significantly with an increase in the anomalous scattering, δB θ /B ∼ ν 2 eff , and as a result, EP losses can increase.This dependence of the fast ion relaxation on the microturbulence intensity does not, however, allow us to conclusively predict the fast ion confinement in advanced steadystate ITER scenarios without making quantitative predictions for the background microturbulence levels.
We have identified a potentially important effect of AEs on EP confinement which is due to EP depletion near the plasma center.This effect is connected to the beam ion and fusion alpha particle current drives which are also depleted near the center, so that the generation of current drive is required for WDM simulations.A self-consistent analysis of the plasma discharge that includes this effect is needed to evaluate its consequences for the plasma scenario.
In summary, we have found that the beam ions injected at 1 MeV lead to stronger AE growth rates in comparison with the effect of fusion alpha particles, which are created as a source that is isotropic in pitch angle.This was not the case in earlier studies of the ITER baseline scenario [10], where NBI injected fast ions had a much smaller (around ten times smaller) beta.On the other hand, the background microturbulence can enhance EP losses in ITER plasmas, which deserves careful consideration.Current applications of RBQ and NUBEAM to the ITER steady-state case have shown a weak loss of fast ions to the wall (at the level of a few percent).
Discussions with Prof. G.-Y. Fu concerning the long-term evolution of AEs are appreciated.This work was supported by the U.S. Department of Energy under contract number DE-AC02-09CH11466.The United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.

Appendix. AE continuum damping
Continuum damping calculations are important for relatively low toroidal mode number AEs [22].It has been shown that perturbative calculations of this damping may have convergence issues [42].A detailed study of finite-Larmour-radius and nonperturbative AE effects on the continuum damping has been published elsewhere [43].
As shown in [22], a rough formula for the continuum damping rate can be derived for the magnetic shear value |s| > 0.3 and when 1 < mϵ < 20.The following expression is valid within a factor of two: It follows from our calculations that the most unstable AEs which can contribute to EP relaxation exist near q min and toward the plasma center.The maximum magnetic shear value in that region is |s| ⩽ 1 (see figure 1) and thus −γ c /ω ≲ 2 × 10 −2 over the range of expected toroidal mode numbers.Recent analysis of Joint European Torus (JET) deuteriumtritium (DT) experiments has confirmed this assessment using kinetic CASTOR-K simulations [44].

Figure 1 .
Figure 1.Plasma profiles of ITER steady-state scenario used in simulations, showing the radial profiles of the beam pressure, P b ; the fusion product alpha particle pressure, Pα; the thermal ion pressure P i plotted in Pascals; the electron density, ne in units of cm −3 ; the safety multiplied by 1/3, i.e. q/3; and the magnetic shear s = q ′ r/q.The horizontal axis corresponds to the square root of the toroidal flux normalized to its value at the last closed surface.

Figure 2 .
Figure 2. Beam ion (left) and alpha particle (right) distribution functions shown as contour maps in the plane of the toroidal canonical momentum Pφ and normalized to the ion energy adiabatic moment λ.Both distribution functions (DFs) are the sums over different signs of the EP parallel component of particle velocities, v ∥ / v ∥ .Notably, in the right-hand panel, the co-passing and counter-passing ions show a discontinuity in DF values which is due to their precession time discontinuity.

Figure 3 .
Figure 3.Figure (a) shows the sum of AE growth and damping rates for each mode computed for the unstable and marginally stable AEs as blue open circles.Net growth rates due to fusion α particles only are shown as smaller red circles.AE net growth rates include the beam ion drive.Figure (b) represents the AE growth rate as a function of the location of its maximum amplitude value on the mode's minor radius.

Figure 4 .
Figure 4.The initial, overshooting time period of AE evolution of unstable and marginally unstable modes prepared by the NOVA/NOVA-C suite of codes and processed by the RBQ2D quasi-linear code over a 3 ms time window.The toroidal mode numbers of some of the most unstable modes are indicated on the figure.

Figure 5 .
Figure 5.The same as figure 4 but characterized by longer RBQ simulation times in order to see the AE saturation.The Alfvén modes are driven by both beam ions and fusion alpha particles.Figure (a) corresponds to the case of relaxed beam ions, whereas the AE growth rates include the alpha particle drive.Figure (b) corresponds to alpha particles' DF relaxation assisted by the beam ion drive.

Figure 7 .
Figure 7. Beam ion (a) and fusion alpha particle (b) density profiles before (solid curves) and after AE saturation.The profiles are computed by the NUBEAM package with constant diffusion due to AEs over 4 s.For the EP diffusion coefficients taken as corresponding to the overshoot region, i.e. at around 2 ms (see figures 4 and 5), the EP profiles are shown as red dashed curves.For the EP profiles corresponding to the long-term saturated phase, i.e. at around 100 ms on figure 5, the NUBEAM predicted profiles are shown as blue dashed curves.