Verification and validation of integrated simulation of energetic particles in fusion plasmas

This paper reports verification and validation of linear simulations of Alfvén eigenmodes in the current ramp phase of DIII-D L-mode discharge #159243 using gyrokinetic, gyrokinetic-MHD hybrid, and eigenvalue codes. Using a classical fast ion profile, all simulation codes find that reversed shear Alfvén eigenmodes (RSAE) are the dominant instability. The real frequencies from all codes have a coefficient of variation of less than 5% for the most unstable modes with toroidal mode number n  =  4 and 5. The simulated RSAE frequencies agree with experimental measurements if the minimum safety factor is adjusted, within experimental errors. The simulated growth rates exhibit greater variation, and simulations find that pressure gradients of thermal plasmas make a significant contribution to the growth rates. Mode structures of the dominant modes agree well among all codes. Moreover, using a calculated fast ion profile that takes into account the diffusion by multiple unstable modes, a toroidal Alfvén eigenmode (TAE) with n  =  6 is found to be unstable in the outer edge, consistent with the experimental observations. Variations of the real frequencies and growth rates of the TAE are slightly larger than those of the RSAE. Finally, electron temperature fluctuations and radial phase shifts from simulations show no significant differences with the experimental data for the strong n  =  4 RSAE, but significant differences for the weak n  =  6 TAE. The verification and validation for the linear Alfvén eigenmodes is the first step to develop an integrated simulation of energetic particles confinement in burning plasmas incorporating multiple physical processes.


Introduction
Energetic particle (EP) confinement is a key physics issue for future burning plasma experiments, since ignition relies on self-heating by energetic fusion products (α-particles). EP transport can affect plasma profiles, beam deposition, and current drive, and can erode reactor walls [1]. Due to the strong coupling of EPs with burning thermal plasmas, plasma confinement properties in the ignition regime are some of the most uncertain factors when extrapolating from existing tokamaks to the international thermonuclear experimental reactor (ITER). Fully self-consistent simulations of EP transport and EP coupling with thermal plasmas must incorporate microturbulence and magnetohydrodynamic (MHD) instabilities with kinetic effects of both EPs and thermal plasmas on an equal footing, which requires an integrated kinetic-MHD simulation model based on the gyrokinetic formalism [2]. Coordinated efforts in verification and validation (V&V) are needed to develop integrated simulation tools for EP transport due to mesoscale Alfvénic instabilities primarily excited by EPs and EP coupling with microturbulence and macroscopic MHD modes mostly driven by thermal plasmas.
The first-principles simulations and reduced transport models are built upon a hierarchical construction of EP transport prediction based on more fundamental constituents by the progression from linear dispersion relation to nonlinear dynamics and eventually to EP transport. Nonlinear V&V will take on an increased importance as gyrokinetic and kinetic-MHD hybrid simulation models progress from linear to nonlinear simulations for understanding EP confinement properties regarding instability saturation mechanisms, interactions between mesoscale EP turbulence with microturbulence and MHD modes, and EP transport statistics. While it is unlikely that different models will agree in all situations, the regimes of deviation will need to at least be characterized and understood. This is a continuous process since models and computational methods evolve in time. As updated results become available from the first-principles models, they will provide new calibration points for the reduced EP transport models and stimulate their further development.
The V&V studies should use a hierarchical approach, starting with test cases from existing experiments and quantities that are well-diagnosed. For this purpose, an NBIheated low-confinement (L-mode) plasma (DIII-D discharge #159243) with many small-amplitude RSAEs and toroidal Alfvén eigenmodes (TAEs), significant flattening of the EP profile, and strong microturbulence [3,4] has been selected as the reference case for V&V studies by the Integrated Simulation of Energetic Particle (ISEP) project, part of the Scientific Discovery through Advanced Computing (SciDAC) initiative. High quality data for the AE structure, frequency, and amplitude as well as the EP distribution, phase-space flows, and intermittent losses are all available from comprehensive DIII-D diagnostics. Taking advantage of this recent experimental progress, the early linear V&V studies [5] have been extended to nonlinear V&V studies of EP transport by using more newly available EP simulation codes and new EP reduced transport models. Linear and nonlinear simulations of AE and microturbulence in this reference case have been carried out by gyrokinetic, kinetic-MHD hybrid, and eigenvalue codes. Modeling of EP transport have also been carried out by reduced transport models. These V&V studies will proceed from linear simulation of instabilities, to nonlinear simulation of saturation mechanisms, to coupling of mesoscale turbulence with microturbulence and MHD modes, and finally to reduced EP transport models. The V&V for the linear simulations of Alfvén eigenmodes reported in this paper is the first step to develop an integrated simulation of energetic particles confinement in burning plasmas.
In this paper, we present linear simulations of RSAEs and TAEs observed in shot #159243 by using five initial value gyrokinetic codes (EUTERPE [6], GEM [7], GTC, GYRO, ORB5 [8]), two initial value gyrokinetic-MHD codes (FAR3D [9], MEGA [10]), and a perturbative eigenvalue code (NOVA-K [11]). Since fast ion profiles have the biggest uncertainty among all equilibrium profiles measured in the experiment, we use the fast ion profiles both from the kinetic EFIT reconstruction [12], which subtracts the thermal from the total plasma pressure, and from the more realistic kick model [13], which takes into account EP transport by the RSAEs and TAEs. The energetic particle distribution function, which is expected to be an anisotropic slowing-down in the experiment, is approximated by a Maxwellian in this V&V for all simulation codes. This approximation may cause some differences between simulations and experiments regarding the AE dispersion relation, especially the growth rate.
Using the EFIT fast ion profile, all simulation codes find that a RSAE is the dominant instability. The real frequencies from all eight codes have a coefficient of variation (CV ) less than 5% for the most unstable modes with toroidal mode number n = 4 and 5. The simulated growth rates of these two RSAE exhibit greater variations with a CV up to 17% for the five gyrokinetic codes, and a CV up to 26% for all eight codes. Mode structures of the dominant modes agree well among all seven non-perturbative codes regarding radial eigenmodes, 2D shape on poloidal plane, ballooning characteristics, radial extent, and radial symmetry breaking. The TAE observed in the outer edge of the DIII-D experiment is not found in these initial value simulations using the EFIT fast ion profile, indicating that TAE is either linearly stable or sub-dominant to the RSAE. Using the outward-shifted fast ion profile from the kick model, GTC simulations find the n = 6 TAE to be the dominant instability in the outer edge, consistent with the ECE data. Variations of the real frequencies and growth rates for this TAE from seven simulation codes are slightly larger than those of the RSAE, partially due to the co-existence of multiple radial eigenmodes with similar frequencies and growth rates. Finally, GTC simulation data, which has been processed by the Synthetic Diagnostic Platform (SDP) [14] to produce electron temperature fluctuations and radial phase shifts, is compared to the corresponding n = 4 and 6 ECE data for the experimental time of interest. The comparisons show no significant differences in radial mode structure for the strong n = 4 RSAE, but significant differences for the weak n = 6 TAE. These linear results provide a necessary foundation for the next step of nonlinear V&V studies.
The rest of the paper is organized as follows. In section 2, we describe the RSAE and TAE observations in the DIII-D experiment, and the equilibrium and profiles of this experiment as used in all simulation codes. In section 3, we compare different physics models and numerical parameters used in this V&V by all simulation codes. In section 4, we quantify the agreements and differences in RSAE and TAE linear dispersion from these eight independently developed simulation codes. In section 5, we process GTC data by a synthetic diagnostic to compare simulation results with the experimental ECE and ECEI data. Conclusions and discussions are presented in section 6.

DIII-D EP experiment for verification and validation
This work uses profiles and magnetic equilibrium obtained from DIII-D shot #159243 during L-mode current ramp phase at t = 805 ms, which has a safety factor, q, with reversed shear and q min = 2.9. Multiple unstable RSAEs and TAEs are excited using early deuterium beam power injected at 70-81 keV, with 4.0 MW of co-current, on-axis NBI, 0.7 MW of cocurrent, off-axis NBI, and 1.7 MW of counter-current, on-axis NBI. This discharge has excellent diagnostic coverage and was examined extensively in studies [3,4,15] of AE-induced fast-ion transport in critical gradient experiments. Figure 1(a) shows the spectrogram of electron cyclotron emission (ECE) data during the current ramp for shot #159243, along with calculated RSAE frequency evolution from an ad hoc model [16]. The model was used to aid in toroidal mode number identification and to constrain the value of q min for the kinetic EFIT equilibrium reconstruction. In the zero-pressure limit, the RSAE frequency is where V A is the Alfvén speed and R is the major radius. The sensitive dependence on q min causes the RSAEs to chirp up in frequency over 20-40 ms as the q-profile evolves. The observed modes are also Doppler shifted due to toroidal rotation, with f lab = f + nf rot , where f rot is the toroidal rotation frequency. In figure 1(a), RSAEs with multiple toroidal mode numbers appear simultaneously at the integer value q min = 3 near t = 770 ms. The relatively constant frequency modes are TAE modes. The approximate TAE frequency near q min is plotted as a dashed line, given by

Experiment fast-ion profile
While diagnostics provide measurements of portions of the fast-ion distribution, we are unfortunately not able to reconstruct the entire experimental fast-ion distribution function. Instead, an estimate of the experimental fast-ion pressure profile is obtained by subtracting the measured thermal pressure from the computed total pressure from the kinetic EFIT equilibrium reconstruction, which is constrained by Motional Stark Effect measurements, external magnetics data, and knowledge of q min from AE behavior. The kinetic EFIT fast ion density profile is shown in figure 2.
As described in [15], a more realistic fast-ion pressure profile was obtained using the time-dependent kick model of AE transport [17]. The kick model computes the probability for AE-induced change in energy and toroidal angular momentum throughout fast-ion phase space and evolves the fast-ion distribution function through the TRANSP-NUBEAM code. In this case, AE mode structures were first computed by the NOVA code and then scaled to match experiment measurements at a single timeslice. The kick model then evolved the mode amplitudes in time ( figure 1(b)) so that the modeled neutron rate matched the measured value. The resulting fastion profile agrees well with experimental measurements [4]. The kick model fast ion pressure profile is used in section 4.2.

Equilibrium and profile comparison
All benchmarking codes use magnetic equilibria calculated from EFIT [18], except for FAR3D and EUTERPE which use the same equilibrium calculated from VMEC [19]. Profiles are obtained using kinetic EFIT calculations. Figure 3 shows the equilibrium and profiles for all codes, as outputted from each code, after the experimental inputs have been internally processed. Note, while this may seem trivial, with a verification exercise of this magnitude it is an absolutely essential first step. Figure 3(a) depicts ten magnetic flux surfaces ranging from ρ = 0.1-1.0 in uniform steps, where ρ is the square root of the toroidal flux normalized by its separatrix value. Figure 3(b) shows the magnetic field magnitude on the low and high field sides of the mid-plane, respectively. Magnetic surfaces align within tolerance and magnetic field magnitudes are almost identical with the exception that it, on the low field side at ρ = 1.0, is 3% larger in NOVA. Figures 3(c) and (d) show that the q profile, density profiles, and temperature profiles are almost identical for all codes.

Simulation models
This paper presents linear AE simulations of the DIII-D EP experiment described in section 2 by using a perturbative eigenvalue code (NOVA-K), and seven non-perturbative initial value codes including five gyrokinetic codes (EUTERPE, GEM, GTC, GYRO, ORB5) and two gyrokinetic-MHD hybrid codes (FAR3D, MEGA). A tabulated comparison summary of the different codes is presented in table 1. Detailed comparisons of all codes can be found in appendix.

RSAE
Using the equilibrium and profiles from figure 3, linear electromagnetic simulations from all codes find an unstable RSAE, peaked at the q min = 2.94 surface, to be the dominant linear instability for DIII-D shot #159243 at 805 ms. The RSAE linear dispersion has been obtained and is presented in figure 4. In the figure, codes are grouped according to physics model via the plot marker used. Namely, diamond, star, and circle markers are used for the gyrokinetic, gyrokinetic-MHD hybrid, and perturbative eigenvalue codes, respectively.
All models show excellent agreement in real frequency. The coefficients of variation of real frequency values, CV ω = σ ω /µ ω , where σ and µ are the standard deviation and mean, respectively, for all data points per toroidal mode number in the dispersion are presented in table 2, which shows that CV ω < 5% throughout the dispersion, with the exception of the subdominant toroidal mode, n = 3, where CV ω = 8.4%. Figure 4 also contains ECE measured frequency values, which are shifted into the plasma frame, and are grouped with square plot markers. The ECE values are plotted for two experimental times in the discharge, 790 ms and 805 ms, to give a qualitative estimation of rate of change of the experimental q min value.  Simulation results agree better with the ECE frequency value at 790 ms, which is due to limitations in the accuracy of the q min calculation in the equilibrium reconstruction. For context, errors as small as 1% in q min can lead to variations in frequency as high as 18%. q min can be manually changed until simulation frequency values agree with the ECE data at 805 ms,but here we accept the EFIT calculations as they are.
Growth rates exhibit greater variations than those in the real frequency comparison. Nonetheless, there is agreement in the general trend of the dispersion, with n = 4 or 5 being the dominant mode. FAR3D and NOVA show exceptions to this trend, as the growth rates are found to increase monotonically. For NOVA, this is expected as some damping mechanisms are ignored, such as the radiative damping for RSAEs which is expected to increase strongly with the toroidal mode number. For TAE modes the radiative damping is added via the analytic expression developed earlier [20]. FAR3D uses Padé approximate fits to the Bessel functions entering into the fast ion moment equations for the energetic particle FLR effects and a first order expansion for the thermal ion FLR effects (in the vorticity equation). Taking into account the perpendicular wave numbers that can be inferred from the dominant components in the calculated mode structures, and the fast and thermal ion energies, these approximations have been checked and should be valid. However, as indicated in figure 4(b), there are deviations in the FAR3D growth rates from the gyrokinetic results particularly at n = 6, and, to a lesser degree at n = 5. This may be due to the subdominant modes having stronger spatial variations that create effective wave-numbers somewhat beyond the range of the FLR approximations. Moreover, the fast-ion gyro-fluid model used in FAR3D uses two moment equation, yet more may be needed. Also, the n = 6 data for GYRO is not presented here, as a coexistent ITG mode was observed to affect the numerical properties of the simulation. The coefficient of variation for n = 4, 5 is CV γ = 16%, 17%, for the gyrokinetic codes, and CV γ = 18%, 26% for all codes. Even with physics model differences, mode structures agree well between all codes. Figure 5 shows the radial structures of the poloidally-rms-averaged dominant poloidal harmonic, and two accompanying side bands, for the n = 4 RSAE. All codes show maximum mode intensity localized near the q min surface, ρ = 0.44, with the FWHM mean value and the coefficient of variation of the dominant poloidal harmonic being ∆ρ FWHM = 0.13, and 0.12, where ρ is the normalized square root of the toroidal flux, CV FWHM = 7.7%, and 15% for the gyrokinetic codes and all of the codes, respectively. The distance between the two q = 3 surfaces is comparable to RSAE radial mode width (∆ρ ∼ 0.16) as shown in figure 5. Figure 6 shows the two dimensional RSAE eigenfunction structures, which show agreement in 2D shape on the poloidal plane, ballooning characteristics, radial extent and radial symmetry breaking. In the MHD limit the Alfvén mode structure is updown symmetric; the presence of an EP component breaks this symmetry and leads to the teardrop/triangular shaped mode structures. These drift effects on TAE mode structures were first presented in figure 14 of [21]; and later discussed for RSAE modes in [5,22,23], and experimentally in figures 6 and 7 of [24]. The differences in linear mode structures are not significant and most likely not detectable in experiments since     other effects, such as nonlinear effects, could cause larger differences than those between codes. Furthermore, the significance of the effects of thermal ion and electron gradients on the n = 4 RSAE instability drive has been examined by several codes, which found that there is a large thermal plasma contribution to the destabilization of the AEs in this case, consistent with theoretical expectations [25]. In the absence of any thermal plasma density or temperature gradients, GEM, GTC, and GYRO growth rates for the n = 4 RSAE are found to be reduced by 100%, 85%, and 62%, respectively.

TAE
In addition to RSAEs, experimental observations also find unstable TAEs at 805 ms, as shown in figure 1. Spatial analysis of the ECE data shows that the TAEs are localized near ρ ≈ 0.75. To see if simulation can find an unstable TAE, the radial domain is restricted to the range ρ = [0.564-0.902] in a GTC simulation to avoid the dominant RSAE. In doing so, linear simulations do show a marginally unstable TAE with ω = 99 kHz and γ/ω = 0.0121. The approximation of energetic particle distribution as Maxwellian may cause some errors in the TAE growth rate. However, previous simulations in other work using an anisotropic slowing-down found only small differences in the TAE growth rate.
To see if a linearly unstable TAE appears without artificial domain restrictions, we use a more realistic fast-ion density profile taking into account transport caused by the AE, calculated from the kick model [17]. The kick model is used, as the actual fast-ion distribution is not able to be reconstructed from measurements. Figure 2 compares the fast-ion density profiles from the kinetic EFIT and the kick model. The figure shows that the kick model calculation predicts a higher density and larger gradients beyond ρ = 0.4 than the kinetic EFIT prediction. These two profiles are used in GTC n = 4 and n = 6  Figure 7 shows the obtained 2D modes strucutres of the perturbed electrostatic potential for these four simulations. The top row shows the n = 4 and 6 mode structures, both of which show an unstable RSAE, with no TAE, obtained using the kinetic EFIT profile. The bottom row shows a transition of the dominant mode from RSAE to TAE as the toroidal mode number increases from n = 4 to n = 6. From figure 7(c), it can be seen that the dominant n = 4 RSAE is accompanied by a lower amplitude TAE at larger radius. The real frequency for the n = 4 RSAE shows almost no difference if using the kinetic EFIT or kick model profiles, but the growth rate is 30% lower when using the kick model profile.
After observing the differences in the fast ion density profiles between the kinetic EFIT and kick model in GTC AE simulations, a verification test is carried out for the n = 6 TAE mode, using the kick model fast ion density profile scaled upwards by a factor of 1.5 times. This increase of the fast-ion density is done so as to ensure that all codes yield unstable results. Table 3 tabulates the calculated TAE frequency and growth rates for n = 6, from seven codes, along with the measured ECE frequency (shifted into the plasma frame). The mean of the calculated real frequencies shows a 6.0% difference from the experimental ECE value, and variation between codes is characterized by CV ω = 7.9%. Growth rates vary much more significantly, however, with CV γ = 33% for the gyrokinetic codes. This discrepancy correlates with different observed mode structures between codes. Figure 8 shows the 2D mode structures for the perturbed electrostatic potentials in poloidal cross sections. Here, it can be seen that there are three patterns of structures, each with one, two, or three local peaks of mode amplitude. The radial eigenfunctions can be Table 3. n = 6 real frequencies (in the plasma frame) and growth rates for DIII-D 159293.00805 simulations using a fast-ion density profile calculated using the kick model, scaled upwards by a factor of 1.5 times, and the corresponding coefficient of variation of the results. ECE frequency data at 805 ms is also included.  seen in figure 9, which shows the corresponding three radial eigenstates with zero, one, and two crossings of the zero value for the electrostatic potential.
The discrepancy between the mode structures between codes is consistent with the discrepancy in the real frequency. This is seen in a frequency scan using NOVA simulations, FAR3D's eigenvalue solver and the CKA-EUTERPE code.   solution space is frequency dependent, but the dependency is different between the different codes.

Comparison to experimental values
Experimental data from the DIII-D ECE radiometer [26] has been obtained for validation purposes. The data corresponds to the lowest radial harmonic of the n = 4 mode, which is averaged over 11 steps in the time range [791.5, 802.5] ms of the DIII-D shot #158243. Corresponding data for the n = 6 mode is also obtained. The diagnostic spans the frequency range [83.5, 129.5] GHz in 40 channels, providing diagnostic coverage in the radial range R = [148.1, 228.7] cm, where R is the major radius, 4 cm above the midplane. We compare the measured and simulated radial structure of the magnitude of the mode amplitude, |δT e |/T e0 , where δT e and T e0 are the perturbed and equilibrium electron temperatures, respectively, and phase profile relative to a specified radial location.
For this comparison, GTC has been interfaced with the open source Synthetic Diagnostic Platform (SDP) [14], where GTC simulation data is processed through SDP. Figure 11(a) shows |δT e |/T e0 , obtained from GTC via SDP, for both the kinetic EFIT (black) and kick model (magenta) fast-ion density profiles and the experimental data (red). The GTC data corresponds to the two n = 4 RSAE-dominated cases in figures 7(a) and (c). All three structures show peak amplitude near R ≈ 198 cm. The full width half max are nearly the same, with that from the kick-model being slightly larger. These results show there is no significant difference in |δT e |/T e0 of the RSAE between simulations and the experimental data, when using either the kinetic EFIT or kick model fast ion density profiles. The experimental data may indicate the presence of radially increasing fluctuations between R = [210, 220] cm, which may correspond to TAE activity; however, the uncertainty in the data is large in that region. Figure 11(b) shows the mode's phase difference for different radial locations, relative to R = 195.0 cm, for the experimental data and the GTC simulations with the kinetic EFIT and kick model fast ion density profiles. The disagreement between the phase values for the GTC simulations with the kinetic EFIT and kick model fast ion density profiles in the outer radial regions is due to the presence of a subdominant TAE near R ≈ 215 cm in the simulation using the kick model fast ion density.
The comparison of the n = 4 simulated RSAE has been extended to 2D ECEI data for the n = 4 mode of DIII-D shot #159243 at 807.3 ms. Figures 12(a) and (b) show GTC n = 4 RSAE data of δT e , processed through SDP, and the n = 4 mode filtered ECEI data, respectively. The comparison shows that the simulated and experimental data agree well in mode location, radial extent, and shape. In both cases, the mode is peaked at the q min = 2.94 location.
Comparison of GTC n = 6 simulation data, via SDP, with experimental ECE data does show a significant difference between experimental data and GTC simulations using the kinetic EFIT or kick model fast-ion density calculations. Figure 13(a) shows the |δT e |/T e0 profiles for GTC simulation results and the experimental ECE data. The ECE data shows peak magnitude near R ≈ 226 cm, and the peak amplitudes from GTC simulations are R = 201.1 cm and R = 210.0 cm for the kinetic EFIT and kick model results, respectively. Qualitatively, there is a better agreement of experimental data  with the kick model result than the kinetic EFIT simulation result. The large discrepancy in location of the peak magnitude between the kick model simulation result and experiment may be attributed to the kick model's prediction of the outward shift of the fast ion density profile gradients being too modest, but further testing is needed to confirm this. Another reason for the discrepancy may be that the experimentally observed TAE is simply nonlinearly generated, which cannot be reproduced in linear simulations. Figure 13 In both figures 11(b) and 13(b), there is a large difference between the simulation and experimental phase values, although the overall radial trends of the phases qualitatively agree. This may indicate that there are dominant nonlinear physics present in the experiment that cannot be reproduced in these linear simulations.

Conclusion and discussion
Using kinetic EFIT equilibrium data from DIII-D shot #159243, gyrokinetic, gyrokinetic-MHD hybrid, and perturbative eigenvalue simulations have obtained the RSAE linear dispersion of toroidal mode numbers n = [3,6], for verification and validation purposes. The simulations are conducted using five initial value gyrokinetic codes (EUTERPE, GEM, GTC, GYRO, ORB5), two initial value kinetic-MHD codes (FAR3D , MEGA), and a perturbative eigenvalue code (NOVA-K). All simulation results predict a linearly unstable RSAE and find excellent agreement in mode structure and real frequency. Simulated RSAE frequencies agree well with experimental ECE values for the experimental time of 790 ms, rather than 805 ms from which time profiles are taken. This discrepancy between simulation results and experimental data is due to small error in the reconstructed equilibrium q min value. Growth rates are found to show a larger variance, Figure 13. Comparison of GTC simulation data, after being processed through the synthetic-diagnostic-platform, to experimental ECE data for DIII-D #159243 at 805 ms. (a) Radial structure of |δT e |/T e0 . (b) The phase profile relative to R = 221.4 cm. The measured n = 6 TAE ECE frequency is 98.9 kHz (plasma frame) and the GTC calculated value is 95.2 kHz. with a coefficient of variation of 18% for the dominant mode number.
Moreover, experimental measurements observe the presence of TAE modes in the outer edge. Therefore linear simulations are repeated with a more realistic fast ion density profile obtained using the kick model, which takes AE induced transport into account. Using this fast ion profile, GTC simulations show that the observed instability transitions from RSAE to TAE as the toroidal mode number is increased from n = 4 to n = 6, whereas no TAE is observed when using the EFIT fast ion profile. TAE simulations from seven codes find variations of the real frequencies and growth rates are slightly larger than those of the RSAE, partially due to the co-existence of multiple radial eigenmodes with similar frequencies and growth rates.
Further validations are obtained by comparing GTC simulation data, processed through the synthetic-diagnosticplatform, for n = 4 and n = 6, using both the kinetic EFIT and kick model predicted fast ion density profiles, to experimental ECE measurements of |δT e |/T e0 and phase profiles. The comparisons shows excellent agreement in radial mode structure of n = 4, for both fast ion density profiles. The n = 4 comparison is also extended to ECEI data, where agreement is found for radial location and the qualitative structure of the mode. The n = 6 comparison shows better agreement with experimental data when using the kick model fast ion density profile.
Previously, V&V studies of the linear gyrokinetic simulations of reversed shear Alfvén eigenmodes (RSAEs) excited by fast ions from neutral beam injection (NBI) in the DIII-D tokamak have been carried out [5] by using a gyrokinetic particle code GTC [27], a gyrokinetic continuum code GYRO [28], and a gyro-Landau fluid code TAEFL [21]. Good agreement in RSAE frequency, growth rate, and mode structure have been obtained among these initial value simulations, and between simulation results and experimental measurements using electron cyclotron emission imaging (ECEI) [29]. The successful linear V&V lends some degree of confidence to nonlinear gyrokinetic simulations [30][31][32] that provide new kinetic insights on nonlinear Alfvén eigenmode dynamics and EP transport, and help the construction of reduced EP transport models [13,33,34] which are needed for fast parameter scans, shot-to-shot analysis, and optimization of ITER experiments.
The linear RSAE V&V study presented here expands on the earlier V&V study [5] with a larger sets of simulation codes and with a co-existing TAE when using a more realistic fast ion profile from the kick model. Nevertheless, robust comparisons of theory and experiment would require nonlinear integrated kinetic-MHD simulations, which can investigate the effects of mesoscale Alfvénic instabilities on EP transport as well as nonlinear couplings of Alfvén eigenmodes with microturbulence and MHD instability. To this end, nonlinear verification studies are needed to converge theoretical calculations, and build a reliable computational toolbox to understand EP transport and aid optimizations of ITER experiments. 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.

Appendix. Simulation models and setups
This section presents a detailed comparison of the simulation models used in each participating code in this study, as well as the simulation set ups used.

A.1. Gyrokinetic model
A.1.1. EUTERPE EUTERPE is three dimensional, full volume, and electromagnetic gyrokinetic particle in cell (PIC) code. It solves the gyrokinetic Vlasov-Maxwell system neglecting B || perturbations. To avoid numerical difficulties associated with the so-called 'cancellation problem', the gyrokinetic equations are formulated using mixed variables [35] and the 'pullback transformation scheme' [6]. It can be interpreted as an explicit reset of the time integrator bringing the system back to the v || scheme [36]. The spatial directions are discretized with B-splines (here B-splines of order two have been used). The code uses Fourier filtering of the perturbations in the angular directions. Furthermore, the sparse matrices resulting from the finite element decomposition have been Fourier transformed and filtered to construct a Fourier solver guaranteeing high accuracy. The spatial resolution is provided by using 150 radial and 128 poloidal splines. It is possible to provide a leading Fourier factor (∼e i(mpθ+npφ) ), which is called 'phase factor' in code terminology. It allows to single out the toroidal direction [37] and allows a lower resolution in poloidal direction.
The Vlasov equation is solved using the so-called δf -ansatz [38], i.e. the distribution function is split into an equilibrium part and the perturbation. The number of marker particles is 64 million for the ions, 256 million for the electrons and 64 million for the fast ions. The equilibrium provided by a mapping from the computational domain extends from r = 0 to r = a with Dirichlet boundary conditions at the outer boundary and natural boundary conditions for the radial finite elements at the inner. Lost particles are re-inserted such that their weight is zero and the constants of motion are preserved.
Although there are several models of different physical complexities installed in the EUTERPE code, such as fluid electrons and/or ions (FLU-EUTERPE) [39] or perturbative kinetic MHD model (CKA-EUTERPE) [40], here always the full gyrokinetic model with a realistic electron/ion mass ratio is used. The electrons are drift kinetic i.e. their gyro-radius is zero, while for the ion species the gyro-averages resulting from the theory are performed with n g -point averages where n g ranges between 4 and 32 and adapts to the size of the gyroradii [41].
A.1.2. GEM. GEM is a gyrokinetic delta-f PIC code that was originally developed for the study of tokamak core plasma microturbulence and anomalous transport. GEM uses the field-aligned coordinates in general magnetic equilibria. Electromagnetic perturbations are included using the parallel canonical momentum formalism. The split-weight scheme [7,42] is used to enhance the time step otherwise limited by the fast electron motion along the magnetic field. GEM also has a fluid electron model for studying the long wavelength energetic particle-driven modes [43]. The fluid electron model consists of the electron continuity equation, the isothermal condition for the electron temperature and the Ohm's law for determining the parallel electric field. Both the kinetic electron model and the fluid electron model are used for the RSAE simulation, and the two models agree. Only the kinetic electron model is used for the TAE simulations. The simulation domain is 0.1 < r/a < 0.9, with a grid resolution of (N x , N y , N ) = (256, 64,64) in the field-line-following coordinates (x, y, z), with 16 particles/cell per ion species. The ion FLR effect in particle motion and weight evolution is treated with four-point averaging. The gyrokinetic Poisson equation is discretized using a full spectral representation of the ion FLR effect in the ion polarization density.
A.1.3. GTC. The gyrokinetic code (GTC) [27,44] is a full torus particle code using both the δf and full-f methods. Thermal and fast ions are simulated using gyrokinetic equations [45]. The electron drift kinetic equation can be solved either exactly using a conservative scheme for both tearing and non-tearing parity [46], or approximately using a fluidkinetic hybrid electron model that removes the tearing parity [47]. The perturbed electromagnetic field is solved for from the gyrokinetic Poisson's equation [48] and Ampère's Law. GTC has been widely used to study microturbulence [49], Alfvén eigenmodes [50], and other low frequency MHD modes in tokamaks, stellarators [51], and field-reversed configuration [52].
For the present work, marker particles are loaded in velocity space according to a Maxwellian distribution, f 0 , and the plasma perturbation is described via the δf method. The electrons are modeled according to the fluid-kinetic hybrid model [47,53]. In the lowest order, or adiabatic limit, electrons are described via the fluid continuity equation. In the higher orders, kinetic effects are solved by the particle method for the non-adiabatic part of the electrons distribution using the drift kinetic equation. In this study, we neglect δB effects [54]. Finite larmor effects are implemented via the fourpoint average method [55]. The simulation domain used is ρ = [0.12, 0.9], the time step size is ∆t = 0.14R 0 /v A0 , where v A0 and R 0 are the on axis Alfvén speed and major radius, ∆r/ρ i ∼ r∆θ/ρ i ∼ 1.7, where ∆r and r∆θ are the radial and poloidal grid spacing, and ρ i is the local thermal ion gyroradius. Each particle species has 20 particles per cell.
A. 1.4. GYRO. GYRO [28,56] is a continuum gyrokinetic code with field-aligned coordinates in kinetic phase space: flux-surface label r = r a , where r is the flux-surface half width at the centroid height and a is the plasma minor radius; normalized parallel orbit time τ ; pitch angle variable where v ⊥ and v are the perpendicular and total velocities and B is the magnetic field normalized to a reference value; and kinetic energy K normalized to the local temperature T: Ê = K T . The toroidal degree of freedom is treated spectrally. The rapid cross-field variation is factored out through a k || = 0 eikonal, with the solver finding the slowly-varying envelope function. The slow θ variation of the envelope is neglected in the equations of motion and in the gyro-average. Trapped and passing particles lie on separate phase-space grid points, with τ normalized as is appropriate. The gyrokinetic equations are solved for three kinetic species: electrons, thermal deuterium ions, and hot beam deuterium ions modeled with an equivalent high temperature Maxwellian distribution. Coupling between species occurs through a Poisson-Ampere field solver performed at each time step. The solver considers electrostatic potential φ and perpendicular magnetic field fluctuations (through perturbed parallel vector potential A || ), but parallel magnetic field fluctuations are neglected in the present study. Both ion species are treated gyrokinetically and the electrons are treated with drift kinetics. The Eulerian time step uses a hybrid implicit-explicit scheme, with the electron dynamics treated implicitly to avoid tracking stiff electron parallel motion.
The present linear simulations are performed over a radial domain r = [0.23, 0.83] or ρ = [0.20, 0.81], with a grid resolution (Nr, N τ , N λ , NÊ) = (150, 20,8,8) at each value of toroidal mode number n. The λ grid includes four passing and four trapped values. Radial grid points are nonuniformly spaced to optimize resolution of flux surfaces. This gives a mean radial grid spacing of nearly two thermal ion Larmor radii, inadequate for resolving ion-scale drift-wave turbulence but well converged for the presented EP-driven Alfvén eigenmodes. The required time step of cs∆t a = 0.01 is smaller than typical values used for simulating microturbulence due to the faster EP dynamics.
Finite Larmor radius effects are accounted for through a pseudospectral gyro-average of relevant potentials, considering a finite radial stencil around the interest point (20 radial gridpoints here). The stencil is wide enough to adequately treat the largest orbits in the simulation. See [28] for details.
A.1.5. ORB5. ORB5 is a nonlinear gyrokinetic PIC code [8] with extension to electromagnetic physics and multiple species [57][58][59]. The p formulation is used and the adjustable control variate method is implemented in ORB5 [57,60], in order to avoid the 'cancellation problem'. Recent development in ORB5 allows larger time step size due to the implementation of the 'pullback transformation scheme' [6,35]. The linear, quadratic and cubic splines are implemented in the code for discretization in radial and poloidal directions and the cubic spline is used in this work. The Fourier representation is used in the toroidal direction. Fourier filters are applied in poloidal and toroidal directions in addition to a field-aligned filter which keeps the poloidal harmonics in the range nq − ∆m m nq + ∆m and in this work, ∆m = 5 for RSAE simulation and ∆m = 16 for TAE simulation.
The simulation is performed in the radial domain ρ = [0, 1.0], with a grid resolution of (N ρ , N θ , N φ ) = (256,192,48) for RSAE simulations and (N ρ , N θ , N φ ) = (256, 256, 64) for TAE simulations, where θ and φ are poloidal-like and toroidal angles. The number of marker particles is 16 million thermal ions, 64 million for electrons and 16 million fast ions for n = 4, 5 RSAEs. Doubled marker numbers are adopted for n = 3, 6 RSAEs. Gyro-averages are included for all ion species and the points number for gyro-averaging is determined by using the gyro-adaptive method [41]. The time step size is dt = 0.065R 0 /v A0 (the permissible dt in TAE case is larger), where v A0 and R 0 are the on axis Alfvén velocity and major radius. While the traditional δf method and the direct δf method [61] are both implemented in the code, the former one is adopted in the simulation.
The large fast ion pressure profile flattening and the electron temper ature fluctuations due to the TAEs observed in a DIII-D experiment were successfully reproduced by comprehensive MEGA simulations [64]. In the MEGA code, the bulk plasma is described by the nonlinear MHD equations, and the energetic particles are simulated with the gyrokinetic PIC method. The electromagnetic field is given by the MHD model. The effects of the energetic particles on the MHD fluid is taken into account through the perpendicular energetic particle current in the MHD momentum equation. Either the full f method or the δf method can be applied to the energetic particles. The thermal ion diamagnetic drift is considered in the MHD momentum equation, and the finite Larmor radius effect is retained for the energetic particle dynamics, using the four-point average method. The spatial derivatives in the MHD equations are calculated with a fourth order finite difference method, and the fourth order Runge-Kutta scheme is employed for the time integration of both the MHD equations and the particle dynamics. For the benchmark results presented in this paper, (128,16,256) grids are used for cylindrical coordinates (R, ϕ, z), respectively, with 0 ϕ < π/2 for the study of the n = 4 RSAE.The number of computational particles is one million.
A.2.2. FAR3D. Gyro-Landau closure techniques [65] allow excitation of Alfvén instabilities within a hybrid (fluidkinetic) global model; this technique was originally implemented and applied in the TAEFL model [21,66], and more recently extended to 3D configurations with the FAR3D model [9,67]. The motivations for such models are: computational efficiency; the fact that the equations can be cast into a matrix eigenmode form, allowing examination of both growing/ damped modes, and the capability to follow long-time scale nonlinear phenomena [68]. For the calculations reported here, the FAR3D model was used; this model is based on VMEC equilibria and can treat both 2D and 3D configurations as well as up-down asymmetric tokamaks. The initial equilibrium is obtained from EFIT; this is converted to a VMEC input file and then recalculated using VMEC. The VMEC data are transformed to Boozer coordinates [69], which are the native coordinates used in the code. The version of FAR3D used for this study included two moment equations (density and parallel momentum) for the fast ion component; options are available with three and four moments, which allow extension of the model [22] to non-Maxwellian distribution functions, such as slowing-down distributions. The fast ion moment equations include four scalar closure coefficients [69] that are selected via calibration against analytical results for Alfvén instability growth rates. The MHD component of the FAR3D model is based on the reduced MHD approximation; a poloidal flux evolution equation (Ohm's law), a toroidal component of the vorticity equation, and a pressure and parallel velocity evolution equation for the thermal plasma (to include sound wave couplings) are included. Toroidal rotation is included, but not used for this paper. Finite Larmor radius effects are introduced into the fast ion equations using Padé approximate fits to the Bessel functions and for the thermal ions using a perturbative approach. Ion and electron Landau damping is included through perturbative terms in the vorticity equation [70]. Since separate equations for thermal electrons and ions are not currently implemented, an ω * ion correction is added to the real frequencies of the modes analyzed in this paper. The equations are solved using Fourier series representations for the poloidal and tororidal angle dependencies; the radial variable is the square-root of the normalized toroidal magnetic flux and a finite difference grid is used in this coordinate. The equations can either be integrated in time, using a semi-implicit stepping procedure, or as a single eigenmode solution, based on a targeted Jacobi-Davidson algorithm. For the calculations reported in this paper 400 radial surfaces were used, with 22 to 30 Fourier modes for the perturbed fields and ten Fourier modes for the equilibrium fields. In most of the cases in this paper, the eigensolver option was used instead of the initial value option, since the instabilities studied here had growth rates that were subdominant to other AE and MHD modes.

A.3. Perturbative eigenvalue NOVA simulations
NOVA and NOVA-K codes are linear hybrid MHD/kinetic codes to study EP driven MHD eigenmode instabilities. NOVA solves ideal MHD equations and finds eigenmodes, such as TAEs [11], including such effects as plasma compressibility and realistic geometry. NOVA-K evaluates the wave particle interaction of the eigenmodes of interest such as TAEs or RSAEs by employing the quadratic form with the perturbed distribution function of energetic ions coming from the drift kinetic equations [71]. NOVA-K is able to predict various kinetic growth and damping rates perturbatively, such as the phase space gradient drive from energetic particles, continuum damping, radiative damping, ion/electron Landau damping and trapped electron collisional damping. NOVA is routinely used for AE structure computations and comparisons with the experimentally observed instabilities [29,72]. The main limitations of the NOVA code are caused by neglecting thermal ion FLR, toroidal rotation, and drift effects in the eigenmode computations. Therefore it can not describe accurately radiative damping for example. Finite element methods are used in radial direction and Fourier harmonics are used in poloidal and toroidal directions. In the results reported here we used the uniform in ψ radial grid with 201 and 256 points in radial and poloidal directions respectively, and poloidal harmonics ranging from 7 to 32.