Nonlinear dynamics and transport driven by energetic particle instabilities using a gyro-Landau closure model

Energetic particle (EP) destabilized Alfvén eigenmode (AE) instabilities are simulated for a DIII-D experimental case with a pulsed neutral beam using a gyro-Landau moments model which introduces EP phase-mixing effects through closure relations. This provides a computationally efficient reduced model which is applied here in the nonlinear regime over timescales that would be difficult to address with more complete models. The long timescale nonlinear evolution and related collective transport losses are examined including the effects of zonal flow/current generation, nonlinear energy cascades, and EP profile flattening. The model predicts frequencies and mode structures that are consistent with experimental observations. These calculations address issues that have not been considered in previous modelling: the EP critical gradient profile evolution in the presence of zonal flows/currents, and the dynamical nature of the saturated state. A strong level of intermittency is present in the predicted instability-driven transport; this is connected to the zonal flow growth and decay cycles and nonlinear energy transfers. Simulation of intermittent AE-enhanced EP transport will be an important issue for the protection of plasma facing components in the next generation of fusion devices.


Introduction
Energetic particle (EP) populations are essential ingredients for fusion plasmas. Plasma heating and current drive require supra-thermal tail populations (neutral beams, ion/electron cyclotron heating) and the DT fusion reaction creates energetic (3.5 MeV) alpha particles that must be sufficiently well confined so that they transfer sufficient energy to the thermal plasma for sustainment of the burning plasma state. These energetic populations can interact with and destabilize a complex collection of plasma waves. It is critical to model and understand such instabilities since they can enhance EP losses above classical levels, leading to significantly degraded heating efficiencies and damage to plasma-facing components. The study of the nonlinear state of EP-destabilized Alfvén modes and associated transport has been of significant recent interest. A variety of approaches are under development, including critical gradient models [1,2], quasilinear broadening methods [3,4], averaged particle displacement (kick) models [5], and nonlinear hybrid kinetic-MHD methods [6,7]. Additionally, several recent nonlinear simulations [8][9][10] have simultaneously accounted for both the effects of the thermal plasma microturbulence and Alfvén instabilities. Here a model is applied that can address several new issues related to EP instability-driven transport. First, while there is general agreement that EP profiles will evolve to some type of critical gradient profile, this has not been examined for DIII-D in a longer timescale, dynamic environment where the effects of zonal flows/currents, and nonlinear energy transfers are active. Zonal flow structures are thought not to have a strong influence on saturation near marginal stability; for example, in [10] a case with γ/ω < 3% was analysed, indicating only weak effects of zonal flows on saturation. The simulation described here considers a multi-toroidal-mode case further above marginal stability with growth rates for the different toroidal modes falling in the range 6% < γ/ω < 18%. The closest to a simulation of this type was that done with the MEGA multi-phase model for the tokamak fusion test reactor experiment [7]. Next, the time interval for evolving from one set of EP profiles to another as the instability drive is changed has not been examined; this is critical for assessing the degree to which profile stiffness effects will govern EP transport. Finally, a self-contained evaluation of intermittency effects driven by EP nonlinear dynamics is important and is demonstrated here. This calculation will focus on the EP driven Alfvén instabilities and not the thermal plasma turbulence; however, we note that global nonlinear gyrofluid models for thermal plasma ITG turbulence [11] based upon similar computational methods as those used here have been developed. In the future this type of model could in principle be coupled with the model described here to study combined Alfvén/thermal plasma turbulence. A reduced dimensionality model (FAR3d) is described for simulating these instabilities and applied to observations from the DIII-D tokamak; its application to the characteristics of nonlinear EP instability saturation and dynamics is demonstrated. This model is based on an extension of gyro-Landau closure techniques [12] to a set of global electromagnetic moments equations that contain the phase-mixing effects [13,14] required to destabilize shear Alfvén eigenmodes (AEs). This model and its predecessor (TAEFL) have been verified/validated with other models and experiments [15,16] and applied to a number of recent EP instability observations in both tokamaks and stellarators [17,18]. Reduced models offer computational efficiency, allowing a more rapid prototyping of new mechanisms and understanding the effects of parameter and profile variations on EP-driven instabilities. A simulation model that is feasible to apply to many different plasma equilibria can be valuable both for discharge optimization and profile uncertainty analysis. It should be noted that such models are intended more for physical insight and analysis of trends than for accurate detailed modelling of specific cases. The application of nonlinear gyro-Landau closure models for tokamak Alfvén instabilities initially [19] focused on the roles of nonlinearly induced E × B velocity shearing (i.e. zonal flows), zonal currents, and self-organization mechanisms on the saturated state of EP-driven Alfvén instabilities. Using a second-order predictor-corrector time stepping algorithm, these effects were then examined over longer simulation times [20]. The version of FAR3d used here employs a 3rd order accurate predictor-corrector algorithm; Adams-Bashforth is used for the predictor step and Adams-Moulton for the corrector step [21]. Although the time-stepping is explicit, this has allowed the simulation to be run with reasonable time steps and extended to longer timescales of increasing experimental relevance.

Description of the nonlinear model
The model used here is based on reduced MHD with acoustic couplings for the thermal plasma and a two-pole (two moments) energetic ion closure model. This leads to the six evolution equations (1) to (6) given below. Here ψ is the perturbed poloidal magnetic flux, U ζ is the toroidal component of vorticity, p th is the perturbed thermal plasma pressure, v ,th is the thermal plasma parallel momentum moment, n fast is the perturbed EP density and v ,fast is the EP parallel momentum moment. These equations (expressed below in coordinate-free vector operator form) are then transformed to Boozer coordinates [22] and normalized; the details of the Boozer coordinates normalized version have been described in [18] and will not be repeated here. The nonlinearities include convective nonlinearities in the vorticity, plasma and EP moment equations, a J × B nonlinearity in the vorticity equation and a B·∇φ nonlinearity in the magnetic flux evolution equation. These nonlinearities generate effects such as zonal flows, zonal currents, and modification/transport of the EP density that are important in establishing a nonlinear saturated state. The destabilization of the Alfvén modes is introduced into these equations through the Landau closure terms in equation (6), which have the L 0 and L 1 coefficients and involve the absolute value of the parallel gradient of the parallel velocity moment. This introduces the phase-mixing associated with the parallel linear Landau resonance (ω = k v ) and allows EP destabilization of the normally stable AEs. The absolute value of the parallel gradient operator is easily obtained in this model due to the use of a Fourier series representation in the toroidal/poloidal angles for all dynamical variables; the parallel gradient then reduces to a multiplication by mode numbers and the absolute value can be readily applied. It should be noted that this is intended as a minimalist model designed simply to trigger the basic de-stabilizing effects of passing fast ions. It should be considered as a surrogate for more complex kinetic resonances in 5D phase space. The field and moments nonlinearities included in these equations for the saturation of the EP instabilities (convective, J × B, B·∇φ) should be universal, and are expected to be active independent of the details of how the instabilities are triggered. The main limitation is that there may be additional saturation effects of a more kinetic origin (velocity space gradient flattening, particle-wave trapping, resonance broadening, etc) that can play a role that are not fully included in this model. where T fast eBn f 0 ·b eq × ∇ L 0 = 2.72, L 1 = −1.31 Here equations (1) and (2) are the Ohm's law and toroidal component of vorticity. Equation (2) is derived from multiplying the momentum balance equation by the Jacobian √ g of the Boozer coordinates transformation and taking the toroidal component of the curl of this equation. ρ m is the ion mass density; T fast is the average energy of the EP component (matched to a slowing-down distribution); B eq is the equilibrium magnetic field; p th , v ,th , n fast , and v ,fast are the fluctuating thermal plasma pressure, thermal plasma parallel velocity moment, fast ion density, and fast ion parallel velocity moment, respectively. p th,eq and n f 0 are the equilibrium thermal plasma pressure and equilibrium fast ion density. v d is the fast ion guiding-center drift velocity, evaluated at the average energy. Each equation includes a diffusive term of the form D∇ 2 ⊥ (. . . ). These provide diffusion of the perturbed fields as needed for a sink of the nonlinear energy cascades to short wavelengths, but also aid in preventing numerical grid separation effects which can adversely affect numerical stability. With the exclusion of the resistivity (which is set based on plasma parameters), the diffusion coefficients are chosen by prior testing; they are chosen large enough to avoid grid separation, but small enough so as not to significantly lower the linear growth rates and quench the instability. The diffusion terms in the fast ion moment equations also act to transport the perturbed moment fields through resonance regions (which for this model are in real space rather than phase space) in a way that is somewhat analogous to the way collisional effects in more kinetic models act with velocity space resonances. Collisional effects have recently been identified [23,24] as an important mechanism in setting saturation levels for Alfvén instabilities. Gyro-Landau closure models have been developed that include collisional effects [25]; this can be an important area for future improvements in this type of model. Each equation also includes a convective term (not shown here) of the form −v ς,eq ∂(. . . )/∂ζ to incorporate the effects of toroidal plasma rotation. The above equations do not include ion/electron Landau damping, or FLR terms for the EP or thermal ion components. These effects will make minor changes in the linear growth rates for this case.
They have been included and tested in the linear version of FAR3d and will be implemented in the nonlinear version in the future. The model used here does include continuum damping effects.

Application to pulsed beam DIII-D case
This model is applied to a DIII-D discharge (#176523) shown in figure 1 where a pulsed tangential beam was used to periodically turn Alfvénic activity [26] on and off. Our modelling will be based on the central pulse outlined by the dashed line box in figure 1(a) where a single dominant instability was present with a frequency near 100 kHz. The decay phase of this pulse after the beam is turned off will be modelled here over about a 5 ms interval. As indicated, the frequency was relatively constant, although the fluctuations in the AE amplitude envelope were large. The virtue of modelling an isolated pulse of this type is that the beam input can be considered as an isolated impulse, precluding the need to develop detailed source and sink models; also, the amplitude decay rate of the fast ion driven instability is measurable both in the experiment and the simulation, allowing a basis for comparison. For the time interval simulated here, the neutral beam fuelling rate on the magnetic axis is ∼4 × 10 19 ions m −3 s −1 ; this is balanced off against the loss of beam ions to slowing-down into the thermal population of ∼1.3 × 10 19 ions m −3 s −1 . Over a 5 ms interval this would lead to about a 3% net increase in the central beam ion density (+1.33 × 10 17 ions m −3 ), which will not substantially modify the simulation results. So, although the model applied here (without a source or sink) is formally designed for the beam turn-off phase, the short time of the simulation interval allows approximate comparison during the beam heated phase as well. As longer-term simulations become possible, modelling of more steady-state regimes and source/sink models is required-this will be a topic for future research. The plasma/beam profiles are shown in figures 1(d) and (e) and the q-profile is shown in figure 2(a). For figures 1 and 2, and subsequent plots, ρ is the square root of the toroidal magnetic flux and it will be normalized to its edge value ρ edge . Figures 2(b) and (c) show the linear growth rates for this case and frequencies of the unstable Alfvén modes for toroidal mode numbers n = 1 to 6 and the location of these frequencies on the Alfvén continuum plot. Here for simplicity only the n = 3 continuum is shown; superimposing the continua for the higher and lower toroidal mode numbers will add further frequency lines but does not significantly change the width of the open gap structures.
For the nonlinear simulation a collection of 277 Fourier mode m/n pairs have been used. Both cos(mθ + nζ) and sin(mθ + nζ) components are included for each m/n pair due to the symmetry-breaking effects introduced by the fast ion gyrofluid equations. This includes n = 0 (m = −14 to 14), n = 1 (m = 1 to 10), n = 2 (m = 3 to 20), n = 3 (m = 7 to 30), n = 4 (m = 10 to 40), n = 5 (m = 13 to 27), n = 6 (m = 15 to 40); a radial resolution of 400 points is used. This collection of modes is chosen based on linear stability convergence studies for each toroidal mode group. The n = 0 modes are not linearly unstable but are required to allow the fast ion density flattening, and zonal flow/current generation that arise from the nonlinear products of the n > 0 modes and are essential in this model for saturating the exponential growth of the instabilities. Although, as indicated in figure 2(b), all the toroidal modes are unstable based on an initial value analysis, a more detailed examination of the eigenvalue spectrum shows that there are many damped modes present at similar scale lengths as the unstable modes; these provide nonlinear damping channels for the growing mode energies. A Lundquist number of S = τ resistive /τ Alfvén = 10 7 is used, which is close to the physical value for the parameters of this experiment, and diffusivities  of Dτ Alfvén /a 2 = 3 × 10 −6 are used in each of the dynamical equations; a time step of Δt/τ Alfvén = 0.1 is used. The fast ion central β EP (0) is 0.007. A reduced value for the adiabatic index Γ = 0.02 is used here; higher values in the range 0.7 to 1.35 have been tried, but these move the dominant nonlinear frequency up into the ellipticity Alfvén eigenmode (EAE) range (∼230 kHz). This effect is described in more detail in appendix A. The typical linear growth followed by nonlinear saturation is shown in figure 3 for the volume-integrated magnetic energy. The n = 0, 3 and 6 toroidal mode amplitudes reach a transition to saturation at around t = 0.5 ms. This causes a change in growth rate for the n = 4 and 5 mode families; n = 1, 2, 4, and 5 then reach saturation together at around 2 ms. These differences can be attributed to the different growth rates and mode amplitude initializations (only n = 1, 2, and 3 are initialized. After all the n > 0 modes reach saturation, the n = 0 eventually dominates due to its role There is also a coupling to an m = 3 structure at inner radii apparent in figures 4(c) and ( f ). Looking at other times in the full simulation indicates that the mode structures follow a complex evolution with cyclic radial expansion and contraction as the EP density profile locally flattens in one region while steepening nearby; the instability opportunistically follows the steeper gradient regions. Also, the zonal flows provide shearing effects, and the zonal currents alter the continuum damping and mode couplings. The fluctuating poloidal magnetic fields are stored at each time step for several different radial positions. An example is shown in figure 5(a) at ρ/ρ edge = 0.2 and indicates both the rapid Alfvén time scale oscillations as well as periodic bursts. This time series data can then be Fourier analysed using a moving time window, resulting in a spectrogram as shown in figure 5(b), which can be compared with the experimental spectrogram shown in figure 1(a). A similar dominant frequency line at around 100 kHz is seen as highlighted by the dashed box in similarity with the experiment. The simulation includes the toroidal rotation profile of the plasma which is introduced through a convective propagator term in each dynamical equation. At earlier times the simulation shows higher frequency components, but these are likely start-up transients related to the fact that the low amplitude initializations used in the simulation do not reflect the (unknown) state of the fields in the experimental plasma. In addition to the 100 kHz activity there is also a continuous low frequency ( about 10 kHz) component; this is interpreted as being a product of the nonlinear beatings of the higher frequency modes but could also contain lower frequency instabilities that will not be specifically analysed here.
Based on the simulation data given in figure 5(a), a similar time-windowed averaging can be performed as was done in figure 1(c). These results are presented in figure 6, where an extension of the above simulation out to 12 ms has been used and the amplitude has been adjusted to a similar level as the arbitrary scaling of the experimental data. The time axis has been shifted so that the simulation matches with the low amplitude data just before the beam pulse. It should be noted that the rise in the experimental fluctuations is caused by the beam turn on and fueling while the rise in the simulation is due to the natural growth from low initial amplitudes. There is at least a similarity in the level and period of the amplitude fluctuations. An absolute calibration of the experimental δB θ /B levels and comparison with the simulation has not been available. However, electron cyclotron emission (ECE) measurements have provided an electron temperature fluctuation level of δT e ∼ 2 eV near the location of the peak mode amplitude. From the simulation δT e can be inferred from the  radial averages around the location where the mode structure peaks at the times indicated by the arrows on figure 5(a) indicates δT e ∼ 3.4 eV at t = 4.65 ms and δT e ∼ −1.7 eV at t = 4.74 ms. Improved comparisons should be possible in the future by designing temporal and spatial averaging in the simulation that more closely mimic those used by the ECE diagnostic.
These instabilities drive EP transport. This can be analysed either in detail by following tracer particles in the time-varying fields or more directly at a macroscopic level by calculating the moments-based radial flux. The latter approach is used here and based on recording the n = 0 nonlinear products of the EP density with the fluctuating radial E × B and v ,EP δB ρ /B velocities leading to a radial particle flux, as indicated below.
From this radial flux, an effective diffusion coefficient can be inferred by dividing the m/n = 0/0 radial flux by the local density gradient χ = −Γ ρ /(dn fast /dρ). The time variation of this transport coefficient is shown in figures 7(a) and (b) at two radial positions over the full simulation time period. As can be seen, the transport coefficient transiently reaches very large values and shows a large degree of intermittency. Also, there is a rapid variation from ρ/ρ edge = 0.2 to 0.3, indicating that the transport likely has strong nonlocal characteristics. AE driven EP transport intermittency effects have been observed in DIII-D [27] for other discharges than considered here based upon fast ion loss detector measurements. These discharges have recently been modelled using a quasilinear critical gradient model that examines [28] intermittency which is driven by a different mechanism (microturbulent noise effects on the AE thermal plasma damping) than that considered here. Figures 7(a) and (b) show the flux surface averaged radial transport rate; however, the radial transport also varies within a flux surface. Figure 7(c) shows a reconstruction the twodimensional radial flux at the final time (t = 5.3 ms) of the simulation. As indicated, the radial flux has a dominant m = 1 structure near the magnetic axis; it is outward (red) in the upper half and outer regions of the plasma and inward (blue) in a small region below the magnetic axis. The flux surface averaged transport fluxes, as indicated in figure 7, exhibit very bursty/intermittent features. These characteristics are further analysed in figure 8 where a spectrogram is plotted in figure 8(a), showing that the transport fluxes have a broadband frequency spectrum up to about 200 kHz, due to the fact that they involve products of quantities (fast ion density, potential, and magnetic field) that have frequency peaks around 100 kHz, as shown in the spectrogram of figure 5. Associated with the bursts, there are also periodic excursions to higher frequency. The spectral kurtosis [29,30] is used as a measure of this intermittency and is shown in figure 8(b) based on the transport fluxes at four different radii. Spectral kurtosis is an appropriate measure for oscillatory signals such as characterize Alfvén instabilities. This is obtained by calculating the kurtosis using time averages of the short-time Fourier transform of the radial transport fluxes and is a measure of non-stationary or non-Gaussian behavior in the frequency domain. It increases with frequency indicating that the higher frequencies present during the transport bursts are characterized by higher levels of intermittency and non-Gaussianity. In addition to this statistical analysis of intermittency, its large peak to average ratio (5 to 10) and bursty nature (periods of ∼0.1 ms) can also have important implications for wall protection. A detailed evaluation of these effects will depend on the materials and geometry of the plasma facing components and is beyond the scope of this paper.
The nonlinearly generated m = 0, n = 0 EP density components manifest themselves through a depletion and flattening of the equilibrium EP density profile. This is illustrated in figure 9(a) where the initial and final EP density profiles are given. The gradient erosion effects from the Alfvén instabilities are shown in figure 9(b), where the derivatives of the profiles from figure 9(a) are displayed. While the current simulation is limited to about 5 ms, longer simulations are underway in order to better test the decay phase of the fluctuations shown in figure 1(a). As indicated, the peak of the fast ion density is depleted and the gradient reduced over the simulation interval; however, this is a very dynamic process, and not just a smooth relaxation of the profile. The nonlinearly generated m/n = 0/0 components of the potential (φ) and poloidal magnetic flux (ψ) are also active in the nonlinear regime. The effects of these fields are plotted in figure 10. In figure 10(a) the zonal flow derived from φ 00 is plotted vs radius for several times and shows an increase in magnitude and a radial broadening during the simulation time. These flows lead to a shearing of the mode structure. In figure 10(b) the effects of ψ 00 on the q-profile are shown, indicating a sort of low-level corrugation near the inner radii. The exact impact of this on the instability has not yet been determined but could lead to radial shifts in gap locations and altered continuum damping effects.
In addition to the above profiles of the m/n = 0/0 fields, the n = 0 modes also include a range of poloidal harmonics. These allow reconstruction of the 2D structures for zonal flows and currents. In figure 11 the contours of the poloidal flow velocity are plotted for two different times in figures 11(a) and (b) and the perturbed toroidal current density in figures 11(c) and (d). These are for the times t = 4.65 ms (high amplitude) and t = 4.74 ms (lower amplitude) indicated by arrows near the end of the time history in figure 5(a). There is a clear zonal structure evident; plots similar to these done over the full simulation interval show similar features, except for the very early growth phase. The main changes in these structures with time are an amplitude variation and a periodic oscillation in the radial location (as for example is present in the (a) to (b) sequence). The amplitude level of the zonal flows in figure 11(a) is ∼1.9 times larger than those in figure 11(b) and the zonal currents in figure 11(c) are ∼1.7 times larger than those in figure 11(d). These radially varying zonal structures support the hypothesis that the zonal flows regulate the nonlinear saturation level of this instability by providing shearing of the growing AE mode structure; as the mode structure is dissipated by the shearing, the amplitude of the instability and the zonal flow level decreases until the mode can again begin growing, thus repeating the cycle.
The global effects of these 0/0 components on the nonlinear evolution can be tested by selectively turning them off  and re-running the simulation. This type of study is shown in figures 12(a) and (b). The largest change occurs when the 0/0 fast ion density component is not present (red curves). In effect, this freezes the fast ion density profile at its initial equilibrium value. Another way to characterize this limit is that it effectively adds a fast ion density source which exactly fills in for any EP density flattening caused by the turbulence. As can be seen, without the feedback from the n fast (0, 0) component (which causes the flattening and density loss shown in figure 9) the amplitude of the fluctuations rises to a much higher level, and in the case of this particular run, resulted in a numerical instability that stopped the simulation at around 1.1 ms.
The intermittency in the transport fluxes shown in figure 7 can be analysed by performing cross-correlations with other quantities. The main candidates that have been considered for this are the nonlinearly generated m/n = 0/0 volume-averaged components of the field variables. In figure 13 these are plotted vs time for the potential, poloidal flux, and fast ion density. The only component that shows a similar rapid time variation as the radial fluxes is the φ(0, 0) component, which is connected with the zonal flow velocity generation. The zonal current levels (related to ψ(0, 0) ) and fast ion density profile changes follow a more gradual monotonic change with time. As indicated, the integrated m/n = 0/0 fast ion density is negative and decreasing in time in consistency with the dropping central fast ion density and profile flattening shown in figure 8(a). The edge boundary conditions for this simulation have been chosen so that the EP radial transport flux (equation (7)) is finite at the plasma edge, allowing fast ion density to be transported out of the system by the EP instability. In figure 13(d) the φ(0, 0) field is cross correlated with the radial transport fluxes for several radial positions over the full-time interval of the simulation; strong correlations are indicated with the fluxes at ρ/ρ edge = 0.4 and 0.5. This suggests that the rapid time variations in the transport fluxes are related to time variations in the zonal flows which play a role in regulating the nonlinear saturation levels. A possible interpretation for this connection relates to predator-prey cycles between the instability and the   zonal flows. i.e. as the instability amplitude grows, it drives sheared flows that break up the radial structures and suppress the coherency of the mode, the instability amplitude drops until the sheared flows are no longer large enough to have an effect, the instability amplitude begins to grow again, until sheared flows again quench the growth, and so on. However, it can be difficult to see these effects directly from the mode structure plots since it occurs at different radial locations for the different toroidal mode groups, and other phenomena (i.e. density flattening, nonlinear energy transfers) are going on at the same time.

Conclusion
The FAR3d gyro-Landau closure model has been applied here to a pulsed-beam DIII-D discharge using both its nonlinear and linear implementations. The nonlinear simulation of discharge #176523 demonstrates the ability of this model to follow long-time scale nonlinearly saturated EP-destabilized Alfvén instabilities without limitations from numerical instability or perturbed EP amplitudes becoming too large. Results are obtained that are consistent with experimental observations with respect to frequency spectra and mode type (TAE).
The reduced model used here is simplified in several aspects (isotropic two-moment fast ion dynamics, only continuum damping, reduced MHD), but provides a foundation for more advanced versions in the future. Linear versions of FAR3d now include electron/ion Landau damping, EP/thermal ion FLR effects, and two fast ion populations. Examination of the full fast ion density decay and profile relaxation will require longer simulations and a more detailed level of comparison with the data which is beyond the scope of the current paper. The nonlinear saturation is most directly regulated by the fast ion density profile flattening and gradient depletion. However, achieving this state is a very dynamic process with strong intermittent variations in transport fluxes. The fast ion density profile in this model does not simply smoothly relax to a critical gradient profile but experiences oscillations about this state. These rapid time variations have been correlated with zonal flow variations, implying that predator-prey phenomena are active between the instability drive and the nonlinearly driven sheared flows.
While the simulation times presented here remain short compared to experimentally observed modes, no limitations have been uncovered yet as to the length of time over which this model can be applied. This motivates further developments and applications, including introduction of models for fast ion sources and sinks, coupling to other types of fast ion transport modelling, extension to higher order moments, and new closure relations that capture kinetic effects and resonances in a more exact way. Applications to fusion reactor scale devices such as ITER which will have much denser collections of unstable Alfvén toroidal modes than current experiments should be feasible. This type of reduced model can ultimately provide scans with more comprehensive signatures of how the fast ion transport characteristics change as plasma conditions vary and can lead to enhanced understanding as compared to looking at simulations of only individual cases.

Disclaimer
This report was prepared as an account of work sponsored by an agency of the United States Government. Neither the United States Government nor any agency thereof, nor any of their employees, makes any warranty, express or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or otherwise does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States Government or any agency thereof.  calculations have been done to determine the appropriate value of Γ for a magnetically confined plasma. A value that is often used has been given in [31] as: For the case studied here this would result in a value for Γ around 1.3 in the region where the instability is dominant. However, the derivation of this value for Γ and its verification have been carried only in the context of linear instability analysis. As mentioned above, the nonlinear simulations done in this paper show best agreement with the observed mode frequencies if smaller values of Γ are used. This conclusion may be an artifact specific to this regime or model and is not necessarily presumed to be generally applicable. In this appendix results are presented which demonstrate that the response of this model to increasing Γ is quite different in the linear regime as compared to the nonlinear regime.
First in figure A.1 the influence of increasing Γ on the frequencies determined from linear stability calculations of n = 1 to 6 mode numbers like those presented in figure 2(b). In some cases, the curves do not extend all the way to right because modes stabilized at the higher values of Γ, and in others (n = 3a, n = 4a, and n = 5a) a jump occurred in the eigenmode. As can be seen there is a gradual upward trend in frequency with increasing Γ, as expected for linear stability effects of increasing Γ, but even at the highest values of Γ, the modes which remain unstable (n = 3, 4, 5, 6) remain in a frequency range (100 to 130 kHz) that is not significantly different from the observations. For the nonlinear simulations, a similar study has been made of the effects of increasing Γ, but with quite different results than for the linear regime. The same parameters and initial conditions are used as for figure 5; results for the fluctuating poloidal field and its spectrogram similar to figures 5(a) and (b) are calculated. These are given in figure A.2. As indicated, both the structure of the time-history waveform and the dominant frequency are quite sensitive to Γ. Increasing Γ even up to 0.35 raises the frequency into the 200 kHz range and it remains there for Γ = 0.7 and 1.35. This is in the frequency range of the elliptical Alfvén gap (EAE) and is an effect of the nonlinear couplings since the linear stability study of figure A.1 does not indicate this level of frequency upshift. Since the Γ = 0.02 results are closer to the observations, this was chosen for further analysis. However, this type of deviation in sensitivity to Γ between the frequency of the linearly most unstable mode and the dominant mode in the nonlinear simulation warrants further study.