The resonance broadened quasilinear (RBQ) model for the problem of relaxing the hot ion distribution function in constant-of-motion 3D space [Gorelenkov *et al.*, Nucl. Fusion **58**, 082016 (2018)] is presented with the self-consistent evolution of multiple Alfvén eigenmode amplitudes. The RBQ model represents the generalization of the earlier published model [Berk *et al.*, Nucl. Fusion **35**, 1661 (1995)] by carefully examining the wave particle interaction in the presence of realistic Alfvén eigenmode (AE) structures and pitch angle scattering with the help of the guiding center code ORBIT. One aspect of the generalization is that the RBQ model goes beyond the local perturbative-pendulumlike approximation for the wave particle dynamics near the resonance. An iterative procedure is introduced to account for eigenstructures varying within the resonances. It is found that a radially localized mode structure implies a saturation level 2–3 times smaller than that predicted by an earlier bump-on-tail quasilinear model that employed uniform mode structures. We apply the RBQ code to a DIII-D plasma with an elevated *q*-profile where the beam ion profiles exhibit stiff transport properties [Collins *et al.*, Phys. Rev. Lett. **116**, 095001 (2016)]. The properties of AE driven fast ion distribution relaxation are studied for validations of the applied RBQ model in DIII-D discharges. Initial results show that the model is robust, is numerically efficient, and can predict fast ion relaxation in present and future burning plasma experiments.

## I. INTRODUCTION

At the moment, the ITER (International Thermonuclear Experimental Reactor), a prototype of future DT fusion reactors, faces several physics challenges which are being investigated in the fusion research community. One challenge is the confinement of energetic charged products (or hot ions) of the DT fusion reaction, which are prone to be lost in the presence of Alfvénic instabilities.^{1} In present devices, strong Alfvénic activity has been shown to significantly flatten the fast-ion profile, reducing the fusion performance by up to 50%. Alfvénic eigenmodes (AEs) are driven by gradients in fast-ion phase space through wave-particle interactions. A major goal is to develop accurate and efficient models in order to predict and design fusion plasmas that avoid undesirable fast-ion transport.

Here, we continue to develop the Reduced Broadened Quasilinear (RBQ) model^{2} to allow the relaxation of the hot ion distribution function (DF) in the presence of multiple Alfvénic modes. The RBQ methodology targets the dimensionality reduction of the full kinetic phase space dynamics and fits the whole problem to the quasilinear (QL) methodology. We note here that many elements of the RBQ model with its realistic approach are adopted from the earlier proposed Line Broadened Quasilinear (LBQ) model.^{3} Several improvements are included in the RBQ model, which we describe in this paper. In particular, we resolve the fast ion dynamics near the wave particle resonance to account for the radial variation of the mode amplitude. This procedure allows us to correct the pendulum motion of the resonant hot ion and change its saturated mode amplitude (or, equivalently, the square of the nonlinear bounce frequency), considerably reducing it by a factor of 2–3 in some cases. Another improvement is that RBQ neglects the contribution of the growth rate to the resonance broadening of the Alfvénic modes, which is based on the results of Refs. 4 and 5. Such improvement allows us to avoid spuriously growing modes in simulations.

The RBQ model is a significant improvement of the critical gradient model (CGM) published earlier^{6,7} which is a limiting case of the further reduced quasilinear model. In the CGM, as well as in the stiff transport model,^{8} there is no velocity space resolution of Energetic Particles (EP) distribution during EP relaxation. This significantly limits the accuracy of the predictions of the saturated state of the system to the level of 50% of beam ion density.^{1} Some of the predictions, such as the beam driven current drive, require much better knowledge about the distribution function in the velocity space.

This paper, together with an earlier publication,^{2} addresses the arguments existing in the literature, implying that LBQ should not be the choice for “viable modeling” of EP losses in the presence of multiple AEs^{9} (also see Ref. 10). We demonstrate that the RBQ formulation inherits all the elements which were argued to be impossible to treat within the LBQ methodology, such as the details of plasma nonuniformities, equilibrium geometry via multidimensional resonance conditions, finite mode structures, and many others which make the chosen approach an appropriate computationally efficient one for the problem of interest in realistic conditions. Here, we also demonstrate that RBQ validation is very encouraging as well as numerically effective for the problem of the whole device modeling.^{2} It is true that more validation exercises are required, but RBQ should be first improved to include the energy variation in the fast resonance ion dynamics within the resonance island.

This paper is organized as follows. In Sec. II, we discuss the system of quasilinear equations including the applied boundary conditions (BCs). In Sec. III, we summarize the implementation of the resonant frequency approach which is fundamental for RBQ methodology being inherited from earlier works by Refs. 11 and 12. Then, we show how the realistic mode structure affects the saturated amplitude of a single mode to be used within the multiple mode RBQ simulations in Sec. IV. We proceed with RBQ1D applications for DIII-D plasmas for the validation exercise in Sec. V. In Sec. VI, we summarize the results and present future plans. In the Appendixes, we offer verification exercises for the AE amplitudes of the expected saturation level of a single mode which is dependent on either the effective collisionality, Appendix B, or the proximity of the problem to the excitation threshold, Appendix C.

## II. RBQ SYSTEM OF EQUATIONS

For the consistency of our presentation, let us write the system of equations published earlier,^{2} assuming the ansatz for the perturbed quantities $e\u2212i\omega t\u2212in\phi +im\theta $. For a single mode, the QL diffusion of an ion is along the paths of constant values of the expression,

where *ω* and *n* are the angular frequency and the toroidal mode number of a mode, $P\phi $ is the canonical toroidal momentum, and $E$ is the unperturbed energy.

Following Ref. 2, we denote each considered mode with symbol *k* and the corresponding resonances with symbol ** l**. If the ensemble of modes is known, the distribution function evolves according to QL equations implemented in RBQ. The system of QL equations in its general form is given in Ref. 29. It was adopted to NOVA-K notations in Refs. 2 and 13 which is trivially generalized to the following 2D form:

where the diffusion coefficients are

and $\nu \u22a5$ is the 90° pitch angle scattering rate, *f* and *f*_{0} are the distribution functions at times *t* and *t *=* *0, respectively, $Fl$ is the resonance window function, and *G _{kmp}* are the wave particle interaction (WPI

**)**matrices. At present, 1D version RBQ uses the action variable $Ik=E/\omega k\u2212P\phi /nk$ just as a toroidal canonical momentum divided by

*n*ignoring the energy dependence. Also, here,

_{k}*I*is the action at a given resonance. Subindex

_{kr}*k*denotes that the action is related to the

*k*th mode. In general, three action variables can be used within RBQ, i.e., $P\phi ,E$, and the magnetic moment

*μ*. In the current version of the code, we consider the subcyclotron frequency range implying that

*μ*= const.

The equation for the amplitude of each mode can be written formally without explicit contributions from other modes although we should note here that 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 $\gamma L,k$ and $\gamma d,k$ are the linear growth and damping rates of the *k*th mode at time *t* and $Ck=\delta B\theta ,k/B$ is the amplitude of that mode.

The RBQ model is designed for applications having both isolated and overlapping modes, whereas the conventional quasilinear theory applies to the multiple overlapping modes. We use the same structure of the quasilinear equations for the ion distribution function evolution but with the resonance delta function broadened over some region in the direction of the action variable *I _{k}*. Resonance broadening is a key element of the proposed diffusion model whose parametric dependencies are verifiable against known analytical asymptotic behaviors for isolated modes, see Appendixes A–C.

If the diffusion is approximately along $P\phi $, i.e., at low mode frequency, the whole problem reduces to a set of one-dimensional equations, which is also discussed in Fitzpatrick's thesis.^{14} An appropriate rotation in phase space reduces the diffusion equation in the RBQ model to a set of one-dimensional equations that are straightforward to solve. The effect of several low-*n* modes on energetic particles cannot be accurately determined without the code that captures the diffusion in full $E,P\phi $, 2Dspace. This is due to a complex multidirectional diffusion ion dynamics in 2D space and a possible resonance overlap in the presence of multiple AEs. Thus, the generalization to the 2D case is a major extension target for the RBQ model in the near future. In this paper, we present the 2D RBQ model formulation although we do not implement it in simulations. The RBQ model is a valuable predictive tool needed at the design stages of burning plasma experiments.

## III. RESONANT FREQUENCY APPROACH JUSTIFIES THE RESONANCE BROADENING METHODOLOGY

Here, we expand the developed earlier model approach of the resonant ion dynamics^{3} to a realistic RBQ approach.^{2} A key element of the model is the resonant frequency in the presence of the *k*th mode which we express as

The corresponding resonance between the hot ion and the mode with index *k* is $\Omega =\xi \u0307=\omega k$. The relevant angular variable to describe the nonlinear hot ion dynamics is $\xi =n\phi \u2212p\theta $ for the case when the oscillation frequency is subcyclotronic.

When the displacement of the particles is much less than the spatial scale size of the eigenfunction, the resonant particles sufficiently close to the center of the resonance island can be approximately described by the nonlinear pendulum equation, i.e.,

where *ξ*_{0} is the initial value of the phase and *ω _{b}* is the nonlinear bounce (trapping) frequency.

The full effective energy $Er$ of the ion due to its motion around the resonant island is a combination of the effective potential energy

and the effective kinetic energy

In the Hamiltonian formalism, potential and kinetic energies are independent of each other as the former depends on angular variables and the latter depends on actions (*φ* and *P _{φ}* in 1D case). However, the nonlinear resonant island connects the potential and kinetic energy components to conserve the full effective energy along the resonant ion motion near the island.

If the ion is trapped inside the island, the trapping region can be expressed through the phase, $0\u2264\xi \u2212\omega t\u2212\xi 0<2\pi $. It implies the range of effective energies required for trapping, $\u2212\omega b2\u2264Er<\omega b2$, with the separatrix being determined by $Er=\omega b2$. For such a full effective energy, $Er$, we have the range of an effective kinetic energy $0\u2264\Omega 2/2<2\omega b2$, with the largest value occurring for $\xi \u2212\omega t\u2212\xi 0=0$. Hence, $\Omega min=\u22122\omega b$ and $\Omega max=+2\omega b$. Thus, the maximum frequency range within the separatrix is

where the broadening coefficient is *a *=* *4. To the best of our knowledge, this result was first verified numerically for a realistic mode structure in Ref. 15 using the ORBIT code. Recently, the broadening in the cases of collisional and collisionless relaxation was studied in Ref. 16.

Equation (9) with *a *=* *4 gives the maximum local broadening which needs to be phase averaged. For example, for the collisionless case, the broadening coefficient for the RBQ model can be calculated by enforcing the saturation level to match the simulation result $\omega b,sat\u22483.2\gamma L,0$^{17} (also see Ref. 3 and Appendix C). In the derivation of the quasilinear theory using action and angle variables, the fast angular dynamics is averaged out and only the slower, diffusive relaxation is kept, which represents the relaxation in the direction of an action. This implies that the resonance islands cannot be resolved in the conventional quasilinear approach, i.e., the resonance platform in the quasilinear theory is a rectangle of width in the action space *P _{φ}* by a length of 2

*π*in angle. However, the resonance island shape accounting for the mode radial structure is accurately resolved in RBQ in order to compute the corresponding broadening along the direction of the action variable. This is considered in detail in Sec. IV.

In contrast, in the full Vlasov approach, the resonance platform has an elliptic shape with the maximum width of $\Delta \Omega =4\omega b$ at the elliptic point (center of the island) and a zero width at the hyperbolic (or stagnating) points. The broadening coefficient *a* can be calculated by enforcing the wave-particle momentum exchange to be the same in both quasilinear and Vlasov approaches. The integration of the momentum density over the elliptic island leads to the relation $a=(3\pi \omega b,sat/\gamma L)1/3$, which for the collisionless case $\omega b,sat\u22483.2\gamma L$ implies that *a *≈* *3.1. We note, therefore, that in RBQ, the maximum extension of the resonance is effectively reduced by a factor of 0.775. For the collisional case, the broadening coefficients were derived in Ref. 2.

## IV. THE EFFECT OF THE REALISTIC MODE STRUCTURE ON THE SATURATED AMPLITUDE

In a recent paper,^{18} the deviation of the resonant frequency from Eq. (9) is explored where it was found that at moderate mode amplitudes, $\delta B\theta /B\u223c3\xd710\u22123$, the contribution from the radial variation of the mode structure is strong as we will see. This variation was shown to cause the deviation of the nonlinear bounce frequency from the analytical nonlinear pendulum bounce frequency.

In a similar fashion to the approach introduced in Ref. 18, here called accurate iterative, the RBQ methodology is extended here to account for ion precession within the resonant island. The amplitude is treated as a function of the action angle *Q* via the poloidal flux, i.e., *ψ*(*Q*), which parametrizes the orbit around the island (as explained in Appendix A). Even with the pendulumlike dependence of *ω _{b}* on the local amplitude, the asymmetry can be found because of the

*ψ*(

*Q*) dependence. Another approach is proposed in Appendix A, called approximate first order Taylor expansion, which accounts for the next-order correction due to the next-order Taylor expansion of the mode amplitude with respect to the action, but all quantities are evaluated at the center of the island. It is a faster but a less accurate approach to obtain the correction to the resonant island nonlinear bounce frequency.

Let us compare the model pendulum,^{15} accurate iterative, and approximate Taylor expansion approaches for a simple Gaussian mode structure. For the purposes of illustration, in this section, we do this comparison in the limit of massless resonant particles so that $P\phi \u223c\psi $ over the particle orbit. Following the example of Ref. 18, we consider two Gaussian mode structures with distinctive widths (a broad and a narrow one), as shown in Fig. 1(a). Figure 1 shows the kinetic Poincaré plots of resonant orbits for two mode structures. The minor variation of the mode amplitude across the island in panel (b) makes the pendulum approximation valid to good accuracy. However, when the amplitude changes appreciably within the island (panel c), the island structure is considerably distorted, making the model pendulum framework a poor description of the system dynamics. The same is true for the first order Taylor expansion ion dynamics.

The equation that parameterizes the surfaces of the Poincaré plots in Fig. 1 is $(P\phi \u2212P\phi ,0)2\u221d\xi (\psi )(cos\u2009Q\u2212cos\u2009Q0)$, where the proportionality sign indicates that the square of *P _{φ}* excursion around its original central value $P\phi ,0$ is proportional to a mode amplitude that multiplies the mode structure $\xi (\psi )$ shown in Fig. 1(a). In Figs. 1(b) and 1(c), the $\xi (\psi )$ dependence is taken using three approaches: $\xi (\psi )|\psi =P\phi ,0$ (i.e., local pendulum approximation, in black), $\xi (\psi )|\psi =P\phi ,0+\u2009(\psi \u2212P\phi ,0)\xi \u2032(\psi )|\psi =P\phi ,0$ (first-order Taylor expansion, in blue), and $\xi (\psi )|\psi =P\phi $ (exact value of the mode structure, in red).

The RBQ1D code includes the option to model realistic resonant islands for each resonance using iterations [Fig. 1(c), red curve] and to account for a proper reduction of the predicted RSAE amplitude. We depict the nonlinear amplitude evolution of the *n *=* *5 Reversed Shear Alfven Eigenmode (RSAE) mode (localized at $q(\psi )|\psi =0.45\psi w=qmin$) in Fig. 2 for different numbers of iterations: zero, two, and six iterations as indicated in the figure. Iterations are included following the approach introduced in Ref. 18 and accounting for the fact that, for a narrow mode structure, resonant particles experience a significant amplitude variation of the perturbed fields as they rotate around the island. The iterations convert the radial location of the mode, represented by *ξ*(*ψ*), into a shift of resonant particle orbits in *P _{φ}* space, while maintaining $E=const$.

Within RBQ simulations, each resonance island is resolved and integrated efficiently to find the nonlinear bounce frequency and the resonance broadening. This is done by using the grid which is adapted to the integral Jacobian which is singular at the island separatrix. To avoid the singularity, a small deviation from the separatrix is allowed. The described procedure provides extremely efficient simulation convergence with just 5 integration points over 1/4 of the island phase.

## V. MULTIPLE MODE RELAXATION

### A. Boundary conditions

Important for our simulations are the boundary conditions (BCs) that we impose on the distribution function (DF) during its phase space relaxation. They are derived from the physical constrains that DF remains fixed at zero value at the plasma wall, the Dirichlet BC, and that any particle reaching the plasma center returns back to the plasma which implies that the Neumann BC, zero DF derivative, has to be enforced in the center. The Neumann boundary conditions in the present RBQ version were applied for the first two points on the *P _{φ}* grid. If the resonances of the modes extend to the center (which we found in simulations) without sources, the number of particles in the central region goes down. This also means that the growth rates cannot recover their initial values. This is why in some verification studies, the Dirichlet BC has to be used, as discussed in Appendix B.

### B. RBQ1D simulations

The relaxation of the fast ion distribution in velocity space is simulated for a DIII-D discharge that was part of a dedicated fast-ion transport experiment. In the experiment, the AE activity during the current ramp of L-mode plasmas was varied shot-to-shot using neutral beam injected (NBI) power. The total AE amplitude increased approximately linearly with the applied NBI power. One particular aspect to point out here is that above a threshold in beam power, the fast ion profiles exhibited resiliency, or remained mainly unchanged, despite increased beam power.

The chosen DIII-D discharge, #159243, used 6.4 MW NBI and had a reversed shear magnetic safety factor profile with a minimum of *q _{min}* = 2.9 at the time of analysis,

*t*=

*805 ms. The same discharge was studied in our earlier paper.*

^{2}The spectrogram of temperature fluctuation data in Fig. 3 shows that multiple RSAEs and Toroidicity-induced Alfven Eigenmodes (TAEs) are present.

An unexpected hollow profile of beam ions was measured using the Fast Ion D_{α} (FIDA) diagnostics.^{19} When the kick model^{20} was applied to the same discharge, a similar feature was simulated. The hollow fast ion profiles were understood as a result of strong localized transport of the primarily passing ions near the plasma center.^{2} An application of a limiting QL model, critical gradient model,^{7} to the same plasma conditions did not find a hollow EP profile behavior but rather monotonic profiles.

Recent improvements to the RBQ model allowed us to study this effect further by modeling the amplitudes of AEs self-consistently. That is, the effect of each mode is included in the distribution function and is processed to other modes via contributions to the growth rate. Figure 4 illustrates such simulations, where all 11 modes evolved together with the resonances overlapped in radial directions which mediate the nonlinear interaction between the modes. As one can see at small amplitudes, all the modes evolve with the growth rates which are constant for about 0.5 ms. Then, the wave particle nonlinearity plays a more important role in modifying the distribution function which leads to the saturation of the modes. At the end of the RBQ1D run, the maximum values of the diffusion coefficients are collected in the form of the probability density function (PDF).^{21}

We note that the accuracy of computing the growth rates and the damping rates is an important factor in RBQ simulations as we observed. This motivates thorough investigation of both damping and growth rates. Within RBQ, the growth rate is favorably compared with the values found by the NOVA-K code. The most dominant damping mechanisms that we have found for the case analyzed here are electron collisional damping and the thermal electron and ion Landau damping. For some modes, the radiative damping could be important.^{22}

The amplitudes of 4 (mostly localized near *q _{min}* surface) modes were found in the predictive RBQ1D simulations to be higher than those reported earlier.

^{19}Partially, this is due to 1D relaxation which evolves EP distribution in the direction of

*P*only, leaving the stabilizing energy gradient unchanged. Another reason for this is that the RSAE damping rates in the present version of NOVA-K do not include the radiative damping and as a result underestimate the total damping rate. According to Eq. (4), the amplitudes have an exponential dependence on the rates. The effect of those localized modes is limited to the low shear region and is the likely contributor to the deviation of the predicted neutron flux deficit from the measured value.

_{φ}The set of PDFs is then used within TRANSP NUBEAM package^{23} to prescribe the fast ion diffusion. The distribution function is computed using the same Monte Carlo package as TRANSP does for the kick model computed diffusion. As a result, TRANSP computes the distribution function for further incorporation into the plasma analysis. For the multimode case analyzed here, the same procedure was implemented and the results are shown in Fig. 5. The distribution function is computed integrating its values from *θ* = −*π* to *θ* = *π*. The results from TRANSP/NUBEAM simulations are further averaged radially over several radial zones in the simulation to reduce the Monte Carlo statistical noise. The average is performed using a Gaussian weight centered at *ψ* = 0.5*ψ _{w}* and with a width of $\delta \psi =0.3\psi w$ so that the results shown in Fig. 5 are representative of a broad region in terms of the normalized minor radius.

The distribution shown here includes the effects of all modes evolved self-consistently within RBQ1D simulations, whereas, previously,^{2} we showed them using only one RSAE mode localized near the *q _{min}* surface. Multiple modes inflict more fast ion transport in comparison with the single mode case as expected [cf. Fig. 4(a) of Ref. 2]. Finally, TRANSP provides the results which are shown in Fig. 6. Here, the diffusion rates used by TRANSP code correspond to only one time,

*t*=

*805 ms, RBQ1D runs, and were kept constant during its simulations, i.e., no amplitude modulation of AEs was prescribed for the computations of the distribution function of the fast ions shown in Fig. 5 and for the neutron rate and density profiles shown in Fig. 6. Visible oscillations in the neutron rate dependence vs time are due to the modulation of the injection power.*

Another feature reproduced by the RBQ1D code is the hollow beam ion density profile, Fig. 6(b). In both cases, kick and RBQ1D, the beam ion density profiles are hollow. This is the feature which was not seen in the critical gradient model simulations,^{2} but is found readily in the RBQ1D runs. Consistent with the previous study, this feature is due to strong diffusion of the RSAEs near the *q _{min}* region and the fact that NBI delivers mostly resonant passing ions with close to single pitch angle velocity distribution. That is, almost all the resonant ions are affected by AEs. A better agreement of the RBQ and kick model profiles with the neutron rate deficit shown in Fig. 5 of Ref. 2 in comparison with Fig. 6 can be understood by the fact that the RBQ and kick models were run in the interpretive regimes. Amplitudes of all AEs used in those simulations were inferred from measurements around

*t*=

*805 ms and then rescaled as a function of time at earlier/later times to match the measured neutron rate.*

## VI. SUMMARY

We have built a numerically efficient quasilinear resonance broadened code RBQ1D in its 1D version. By resolving the nonlinear wave particle interaction dynamics, the RBQ1D self-consistently computes the evolution of the Alfvénic modes. The saturation of AE amplitudes accounts for resonance overlapping in the presence of multiple modes. As a result, the RBQ1D code provides a matrix of the diffusion coefficients for subsequent whole-device-modeling simulations performed by the TRANSP code. TRANSP simulations show that RBQ1D captures many features observed in the experiments including the hollow beam ion distribution radial profile.

The RBQ1D code was rigorously verified against several known or newly developed methods. As we elucidate above, the effect of multiple Alfvénic modes on energetic ions cannot be determined accurately without numerical evaluations that analyze multiple resonances and capture the full 2D diffusion in the constant of motion $E,P\phi $ space. Proper treatment of fast ion transport in phase-space can significantly affect the accuracy of several calculated quantities, including neutrons which are produced by fusion reactions. The cross section of Deuterium-Deuterium reactions is sensitive to the fast ion energy.^{24} Ignoring the energy change will not allow for an accurate evaluation of the energy loss of a single ion and will lead to the miscalculation of the neutron source. This is due to the complex particle dynamics in 2D space and the likely resonance overlap which was already demonstrated in the 1D case. Another point to make is that the modifications of the velocity space gradients can be accurately computed only with the 2D diffusion properly addressed and are important (i) for the stability calculations, (ii) for correct saturation levels, and as a consequence (iii) for correct computation of the EP confinement. We conclude that the generalization to the 2D case is a major extension of RBQ code for the near term.

## ACKNOWLEDGMENTS

We acknowledge the enlightening and stimulating discussions we had with Professor H. L. Berk and Dr. R. Nazikian. This manuscript has been authored in part by Princeton University under Contract No. DE-AC02-09CH11466 and is based upon the work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, using the DIII-D National Fusion Facility, a DOE Office of Science user facility under Contract No. DE-FC02-04ER54698. The United States Government retains and the publisher, by accepting this article for publication, acknowledges that the United States Government retains a nonexclusive, 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 A: NONLINEAR BOUNCE FREQUENCY FOR THE MODE STRUCTURE THROUGH THE FIRST DERIVATIVE IN RADIUS

Consider the following resonant particle Hamiltonian, where for simplicity, we ignored the subscript *k* for the action variable and mode frequency,

where $|H1/H0|\u226a1$ and *I _{r}* is the action at the center of the resonance. For the purpose of evaluating the bounce frequency and the resonance broadening more realistically, we are interested in analyzing the effect of particle action deviating from its value at the central resonance point. It is therefore convenient to choose a canonical transformation that leads to a new action being $\delta I=I\u2212Ir$. The old set of coordinates is

the new set is chosen as

This transformation can be accomplished by defining a generating function *F*_{3} (Goldstein's notation was adopted, in which it is a function of the new coordinate and the old momentum^{25}) defined as

A new Hamiltonian can be defined as

Let us expand the old Hamiltonian around the resonance center up to second order in the perturbation,

Noting that, in the absence of the perturbation, we have $\xi \u0307=\u2202H0\u2202I=\Omega (I)$, at the center of the resonance, we have $\Omega (I=Ir)=\omega $, and therefore, $\u2202H0\u2202I|Ir=\omega $. This allows the new Hamiltonian to be written as

where the prime denotes the partial derivative with respect to *I*. From Hamilton's equations,

In order to calculate the bounce frequency *ω _{b}*, we can use

The integration is along the path of constant *I* so that $\omega b=\omega b(I)$. When *δI* = 0, the bounce frequency becomes the inverse of a complete elliptic integral of the first kind, which is the solution of the problem under the pendulum approximation. If $\delta I\u22600$, we can suppress *δI* in favor of *K* and write

Therefore, the lowest order correction to the pendulum approximation involves first-order derivatives of the matrix elements via the term $H1\u20322(Ir)$.

### APPENDIX B: COLLISIONAL EFFECTS ON THE SATURATION

In the quasilinear framework, the effective Fokker-Planck collisional operator is exhibited primarily through the effective pitch angle scattering in the regimes of interest since it is usually much larger than the effective drag-induced transport in the resonance, in the quasilinear regime. Expressions for the effective drag and scattering are shown in Eqs. (6) and (7) of Ref. 26. Within RBQ1D simulations, we neglect the Coulomb drag and assume that the hot ion relaxation in the velocity space occurs due to WPI and the collisional scattering. Here, we focus on the verification through the saturation amplitude dependence on the scattering frequency.

The Petviashvili formula for the saturation amplitude predicts that, regardless of the closeness to the instability threshold, one obtains that, at saturation, the value of *ω _{b}* is linearly proportional to $\nu eff=[\nu \u22a5R2\u27e8v2\u2212v||2\u27e9]1/3(\u2202\Omega l\u2202I|Ir)2/3$.

^{27,28}This means that $A\u223c\nu eff2$ scaling is expected. We are demonstrating in Fig. 7 that RBQ exhibits close scaling.

An important element of this RBQ verification exercise is the way that the boundary conditions are set. They are set different from the nominal physical BCs used in RBQ runs as discussed in Subsection V A. This is done by specifying the Dirichlet BC with $f|P\phi =P\phi w=0$ and $f|P\phi =P\phi min=const,$ where $P\phi w$ is the maximum value of the canonical momentum ions present when they are at the plasma wall and $P\phi min$ is its minimum value corresponding to the canonical momentum near the axis. We note that the distribution function is computed by the TRANSP code and used within NOVA-K as described in Ref. 6.

### APPENDIX C: VERIFICATION OF THE EXPECTED SATURATION LEVEL

The following three cases are relevant to the verification of quasilinear simulations depending on the value of the ratio of the growth to the damping rates.

##### 1. Collisionless, undamped case

When damping is absent, the relaxation of the fast ion distribution was modeled earlier.^{17} For this case, the expected saturation level for the nonlinear bounce frequency *ω _{b}* is

The corresponding dependences are plotted in Fig. 8.

##### 2. Collisional case far from marginal stability

For this case, the analytical prediction for the saturation of a single mode is

The mode saturation is achieved with one overshoot oscillation as shown in Fig. 9.

##### 3. Collisional case close to marginal stability

For this case, the analytical prediction for the saturation of a single mode is

This is shown in Fig. 10.