Development and application of a hybrid MHD-kinetic model in JOREK

In nuclear fusion reactors, energetic particles (EPs), with characteristic energies much larger than the thermal energy of the bulk plasma, can arise due to fusion reactions or external heating systems. Confining these energetic particles is crucial for sustaining a burning fusion reaction. These EPs can interact resonantly with magnetohydrodynamic (MHD) waves or instabilities of the bulk plasma, leading to outward transport and possibly deconfinement of EPs. In present-day devices, EPs normally only provide a small fraction of the total pressure, while in future burning reactors, the EP pressure is high compared to the bulk plasma pressure. Thus, for predicting the performance and optimizing the design of future fusion devices, understanding the interaction of EPs with the bulk plasma in a regime with high EP pressure (but not necessarily high EP density) is a key area of research.

For simulating the interaction of EPs with a bulk plasma, a common technique is the hybrid MHD-kinetic model,¹ employed by, for example, MEGA,² (X)HMGC,^3,4 XTOR-K,⁵ M3D-K,⁶ and M3D-C1-K.⁷ Here, the bulk plasma is treated using an MHD model, while the EPs are treated kinetically. This approach saves computational time compared to a fully kinetic treatment such as implemented in ORB5,⁸ but can still reproduce the relevant physics.⁹ 

In this work, a hybrid MHD-kinetic extension of the non-linear extended MHD code JOREK¹⁰ is introduced, capable of simulating EP-driven instabilities in a high EP pressure discharge using realistic, experimental plasma parameters and EP distribution functions. A full-f formulation is employed for the EPs, and an anisotropic pressure coupling to the MHD fluid is used. The JOREK code contains a broad range of different MHD models, can simulate up to the first wall, and has proven capabilities for challenging, highly non-linear MHD scenarios. It also has a versatile kinetic particle module, capable of simulating many different types of particles, varying from slow impurities to relativistic electrons. This module is currently being ported to graphical processing units (GPUs).¹¹ These features make it an attractive option for the simulation of a wide variety of EP instabilities. Earlier work¹² used the JOREK kinetic extension for toroidal Alfvén eigenmodes (TAEs)¹³ in the presence of EPs, but the model used therein is generalized to anisotropic EP distribution functions in the present work.

This hybrid extension of JOREK is benchmarked to other codes for TAEs and energetic particle modes (EPMs),¹⁴ with good agreement regarding mode structures, frequencies, growth rates, and EP phase-space resonances in the linear regime. The code is then applied to a high EP pressure discharge in the ASDEX-Upgrade (AUG) tokamak to validate the model and show its potential for simulating scenarios relevant to ITER and DEMO.

The paper is structured as follows. In Sec. II, the JOREK code and the EP simulation model are introduced. In Sec. III, the phase-space diagnostic in JOREK is explained. Linear benchmarks with other EP-simulation capable codes are shown in Sec. IV. Results from the application of the developed code to a high EP pressure discharge in the ASDEX-Upgrade (AUG) tokamak are shown in Sec. V, and a summary and outlook is finally given in Sec. VI.

The JOREK code¹⁰ is a 3D non-linear extended MHD code, capable of simulating tokamak plasmas in realistic X-point geometry using a broad range of models. It uses implicit time-stepping, a finite-element discretization in the poloidal plane and a Fourier series expansion in the toroidal angle. Reduced and full MHD models with various extensions are available. A comprehensive overview of the models, numerics, and applications is available in Ref. 10.

A kinetic extension of JOREK was developed initially for the simulation of edge impurities during edge localized modes¹⁵ but has been extended for other applications such as runaway electrons,¹⁶ edge physics, and ion temperature gradient turbulence studies.¹¹ It uses a particle-in-cell scheme to solve the kinetic Boltzmann equation for a particular species, which can be coupled to the MHD fluid if desired. The kinetic particle module is coupled to the MHD fluid by projecting moments of the distribution function on the finite-element representation of this MHD fluid and using these projections in the timestepping of the MHD fluid. The specific terms coupled to the MHD fluid can depend on the application (e.g., collisionless EPs have different coupling terms than heavy impurities).

A “full-f” particle-in-cell scheme is used for the solution of the distribution function of the EPs, such that the distribution function f for N_p marker particles is

where w_i indicates the number of physical particles the ith marker particle represents. The marker particles are pushed using one of the available pushers, ranging from full-orbit to relativistic gyro-center. Collisions with the bulk fluid and the kinetic particles can be used, and ionization, recombination, and radiation models are included. The particle pushing is parallelized on central processing units (CPUs) using hybrid OpenMP-MPI (Message passing interface) while for GPU acceleration OpenACC is used. Marker particles can be initialized using analytical profiles for spatial and velocity space distribution functions, or from arbitrary, numerical distribution functions, such as a realistic distribution function produced by the NUBEAM code.¹⁷ 

For EPs, collisions are neglected and a full-orbit Boris pusher¹⁸ is used to retain all finite orbit width (FOW) and finite Larmor radius (FLR) effects. Guiding and gyro-center pushers are available and may be used in future applications. The coupling between EPs and the bulk MHD fluid is provided by the pressure coupling scheme¹ in the MHD momentum equation as

where

is the pressure tensor of the EPs calculated from the EP distribution function f_h. ρ is the bulk fluid mass density, V is the bulk fluid velocity, and p is the bulk fluid pressure. $J t$ is the total current carried by the bulk plasma and EPs, and B is the magnetic field. In the latter equation, m_h is the EP mass and v the velocity coordinate of the distribution function f_h. The symbol $⊥$ denotes the component perpendicular to the magnetic field. The current coupling scheme has been implemented as well, see Ref. 12, but was too noisy to obtain useful results. In a previous simplified implementation,¹² only scalar pressure was considered. Although the full tensor was implemented, it proved to be unstable and noisy, such that for practical purposes the full tensor is in the rest of this work simplified to

with

This tensor is calculated using Eq. (1) as

where $b i$ is the magnetic field at the location of marker particle i. The pressure tensor is then projected onto the fluid finite elements using smoothing terms to mitigate noise that are always chosen small enough not to affect physical results. For details of this projection procedure, see Refs. 10 and 15, while the effect of the smoothing terms and their impact on the physical results is investigated in Ref. 12. This projection is averaged over several orbits (commonly the MHD time step) to remove high-frequency noise and then entered as an explicit momentum source term into the bulk fluid timestepping. The bulk MHD fluid timestepping remains implicit, and it is only the contribution of the EPs that is considered explicitly. This is a common strategy across different hybrid codes.^3,7,19 This coupling has been implemented in both the full²⁰ and reduced MHD models in JOREK.

To analyze details of the EP dynamics, a versatile diagnostic has been developed which can provide insight regarding real space as well as velocity space dynamics. Similar kind of diagnostics exist in other codes.^21–24 Consider some single-particle quantity g_i, for example, the density by choosing g_i = 1, or the energy transfer between bulk and EPs by using $g i = Δ E i$⁠, where $Δ E i$ denotes the change of particle kinetic plus potential energy. The distribution of g_i throughout the real and velocity space is given by the following expression:

Although this yields all information about the distribution of g_i, it is an inconvenient representation. Some coordinates might not be of interest (e.g., gyrophase or toroidal angle), and some coordinates are not natural coordinates for particles moving in electromagnetic fields (e.g., v is inconvenient compared to coordinates like the magnetic moment μ and parallel velocity $v ∥ )$⁠. Therefore, the idea is to apply a coordinate transformation and integrate over coordinates that are ignorable for the considered analysis. Thus, for the distribution G of g_i in the coordinates $ξ$⁠, integrate over the ignorable coordinates to obtain

An example usage is the investigation of resonances by visualizing the distribution of energy loss or gain (integrated over some time) in terms of $( v ∥ , μ )$⁠, yielding

Delta-functions are not useful for a diagnostic, as they consume a large amount of memory (for large number of particles), cannot be represented graphically, and are noisy. To solve these issues while keeping the useful integration properties of the delta-function, the diagnostic replaces the delta-function by finite-support kernels. These finite-support kernels K_l are defined by

where l is the so-called bandwidth. This controls the smallest visible structures and thus the amount of smoothing. Denote the amount of dimensions of $ξ$ by α. The α-dimensional kernel $K α ( ξ )$ can then be written as

where $l = { l j }$⁠. This multi-dimensional kernel still has the property

Then, by substitution with $δ ( ξ − ξ ) → K l α ( ξ )$ the diagnostic implemented in JOREK is obtained

Numerically, an α-dimensional grid is constructed. For every marker particle I, its contribution is added at every position $ξ$ where the kernel is non-zero. For some applications (such as resonance visualization), it is useful to integrate the projection over time to reduce noise, which in this case means adding contributions of many timesteps corresponding to many gyro-motions.

Applications include visualizing the full distribution function in various coordinates, for example, in terms of $( R , Z , E , λ = v ∥ / v )$ or $( E , μ , P ϕ )$⁠, visualizing resonances as a function of various coordinates, for example, in terms of $( v ∥ , μ )$ or $( E , P ϕ )$⁠, visualizing resonant particle density as a function of the minor radius r or poloidal magnetic flux ψ to probe EP transport or saturation mechanisms, etc. For several figures in this work, this diagnostic is already used (e.g., Figures 10, 4, and 6).

For resonance visualization, a complicating factor is that during a gyro-motion the particle may lose and gain energy, while there is no net energy transfer over the full gyro-motion. Resonant particles mainly lose net energy due to the interaction of the drift velocities with EP instabilities.²⁵ This net energy loss can be much less than the magnitude of the oscillation during the gyro-motion. For a full-orbit particle, $( v ∥ , μ )$ are oscillating during the gyro-motion, and the energy loss or gain is deposited at slightly different coordinates. This leads to the actual resonances potentially being obscured by an irrelevant oscillation (which might still be physical and not noise; i.e. particles all gain and lose energy at the same phase-space locations). The gyro-motion of the particle is then fitted to obtain only the trend of $v ∥$ and μ, without the oscillation, such that the energy loss and gain during the gyro-motion cancels. This procedure is shown in detail in the  Appendix. However, for resonance visualization in coordinates that are conserved for full-orbit particles (such as $E , P ϕ$⁠), this is not necessary.

In this section, two separate benchmarks are considered. The first is the well-known International Tokamak Physics Activity (ITPA) case²⁶ concerning a TAE in a high aspect ratio tokamak. The second is the ASDEX-Upgrade-NonLinear Energetic particle Dynamics (AUG-NLED)^9,27 case, which is based on a high EP pressure discharge in the realistic geometry of an ASDEX Upgrade²⁸ tokamak discharge. The results are compared to other codes in terms of mode structure, frequency, growth rates, and phase-space resonances.

The ITPA case and the results from other codes are described in Ref. 26. It concerns a high aspect ratio tokamak (major radius R = 10 m, minor radius a = 1 m), with a flat hydrogen bulk fluid density (⁠ $n b = 2 × 10 19$ m⁻³). The q-profile is $q = 1.71 + 0.16 ( r / a ) 2$⁠, such that an n = 6, m = 10, 11 TAE gap is expected at $r = 0.5 a$⁠. The EPs follow a Maxwellian temperature distribution with varying fast particle temperature T_f. The simulation is restricted to a single toroidal mode number n = 6. This benchmark was previously performed in JOREK for the reduced MHD model using isotropic pressure coupling in Ref. 12, but here the full MHD model is used with the anisotropic pressure coupling scheme.

The MHD time step was chosen as $3 τ A$ where $τ A ≡ a μ 0 ρ 0 / B$ is the Alfvén time (ρ₀ is central mass density, B is central magnetic field, and a is the minor radius), while the particle time step was chosen as $2 × 10 − 10 s ≈ τ g / 200$⁠, where τ_g is the time it takes for the EPs to complete one gyro-motion. The grid was chosen to feature 102 radial elements with 128 poloidal elements. 40 × 10⁶ numerical particles were used.

The poloidal mode structure and frequency for a simulation with T_f = 400 keV are shown in Fig. 1, exhibiting a clear TAE structure. The energy in the n = 6 harmonic is shown in Fig. 7, where also the linear phase is indicated. A simulation setup using T_f = 700 keV was repeated for varying particle time step and varying particle number to ensure sufficient convergence. These results are shown in Fig. 2, where all results agree quantitatively within 5%.

Quantitative growth rate and frequency comparisons are given in Fig. 3, showing reasonable agreement. In the other codes, the relaxation of the EP distribution function has been switched off, while this relaxation is intrinsically present in JOREK due to the full-f scheme. The initialized Maxwellian distribution function relaxes as the distribution function cannot be expressed solely in terms of conserved particle quantities, changing the real and velocity space distribution functions and thus the EP drive (Fig. 4). This relaxation is visualized for $T f = 400$ keV in 4. Therefore, it is difficult to compare results one-to-one in this benchmark, motivating future comparisons with fully stationary distribution functions. In Fig. 5, the phase space resonances are visualized, showing that the $v A / 3$ resonance is dominant, in agreement with results in Ref. 31. These theoretical resonances are derived using assumptions on the aspect ratio (for details, see Ref. 25) and therefore do not hold exactly.

For simple non-linear behavior, the 400 keV simulation has been continued into the non-linear phase. Using the results from Fig. 5, resonant particles are selected, and their density before and after saturation is shown in Fig. 6. Density flattening of the resonant particle density is observed in the vicinity of the mode location, as expected. Although the pressure coupling scheme does not conserve energy exactly, Fig. 7 shows that energy is conserved rather well until the saturation phase. This simulation was performed with fixed equilibrium (n = 0), and non-zero resistivity and viscosity. Therefore, the energy associated with the dissipation of the mode is not conserved, leading to a difference between EP energy and bulk energy at the later, saturation stage, as shown in Fig. 7.

The AUG-NLED case is described in Refs. 9 and 27. It is based on an experimental AUG discharge (#31213 at t = 0.84 s) with high EP pressure. Here, only the case with off-axis peaking of the EP distribution is considered and the reduced MHD model is used as results agree with full MHD, while computational costs are lower. The EP distribution is assumed to be an isotropic Maxwellian at T = 93 keV, while the EP density is varied. The poloidal mode structure and frequency spectra for the nominal case (as described in Ref. 9) are shown in Fig. 8. From the mode structure and the mode location in the continuum, this mode can be identified as an EPM. The Maxwellian initial distribution is not stationary as the distribution function cannot be expressed in terms of conserved particle quantities. Thus, the distribution function undergoes a fast relaxation in the full-f treatment applied here. This case in realistic geometry relaxes much more than the ITPA case such that results cannot be directly compared anymore without compensating for the relaxation in some way. In order to still be able to compare growth rates and frequencies, the pressure gradient after relaxation was varied using two different methods: first, by simply increasing the weight of the particles and thus scaling initial EP total density, and second by modifying the initial profile such that the gradient at the mode location after relaxation is higher. The equilibrium (and thus the Alfvén continuum) is not changed in this procedure though (as increasing the density artificially would lead to a very low bulk mass density).

In all cases, 10 × 10⁶ particles were used. 51 radial elements were used with 64 poloidal elements. Timesteps of $0.26 μ$s $≈ 4 τ A$ were used for the MHD fluid and $2 × 10 − 10$ s for the EPs.

The comparison with other codes is shown in Fig. 9. As in the ITPA case, the comparison is not directly one-to-one anymore due to the real and velocity space relaxation, but the agreement is good in terms of the dependency of EP drive on the pressure gradient, in terms of the mode structure excited, and in terms of the mode location in the continuum. The frequency does not match as closely. In the original benchmark results, the bulk ion density was different for different EP densities, while the bulk ion density was constant in the JOREK runs (as noted above). Therefore, the Alfvén continua in the JOREK runs will be slightly different compared to the original benchmark results. The frequencies of the waves that can be excited will thus differ as well, providing an explanation for the difference in frequencies between JOREK and the original results for higher EP densities.

Phase-space resonances are also compared with MEGA for the nominal case in Fig. 10, showing similar features but quantitative differences. These differences can arise due to the very different models (full-f and full-orbit in JOREK compared to δf and drift-kinetic that was used in MEGA), and the relaxation of the distribution function in real and velocity space. Furthermore, note from Fig. 9 that MEGA finds a frequency near 165 kHz in the nominal case, while JOREK (without correcting for gradient) obtains about 150 kHz, which also impacts resonances. As EPMs, even in a linear description, require a nonperturbative analysis of the EP response,¹⁴ these differences in the distribution function and the frequency of the mode could explain the quantitative differences.

To validate the code with realistic parameters, it is applied to a high EP pressure discharge using experimental plasma profiles and a realistic EP distribution function. The AUG discharge considered is, as in the AUG NLED benchmark, discharge #31213. The spectrogram is shown in Fig. 11. The goal of this section is to reproduce some characteristics of the spectrum, such as emergence of modes and frequency sweeping of other modes. To this end, several timepoints are simulated. These time slices are taken at ${ 0.6 , 0.65 , 0.7 , 0.75 , 0.8 , 0.84 , 0.9 , 0.93 }$ and are indicated in Fig. 11. The simulations are restricted to n = 1 for simplicity.

AUG is equipped with a large number of diagnostics, which are used via an integrated data analysis (IDA) framework to obtain accurate equilibria profiles.^33–37 The equilibria and plasma profiles are directly imported into JOREK for these simulations. As the spectrum is expected to consist of core-located modes, the simulation domain does not extend beyond the separatrix. The bulk ion charge density can be obtained by subtracting the fast ion (charge) density from the electron (charge) density. The bulk mass density can be obtained by dividing the ion charge density by the effective charge of the ion population and multiplying by the effective mass of the ion population. The main impurity is assumed to be Boron (atomic mass of 10–11) as in the IDA equilibrium reconstruction³⁴ and the rest of the plasma is deuterium, such that the effective atomic mass per unit charge is approximately two. The choice was made to treat the plasma as a fully deuterium plasma for the purposes of the bulk mass density (as in AUG NLED benchmark setup).

The MHD temperature $T = T i + T e$ is set as to reproduce the pressure of the IDA equilibrium. The $F F ′$ profile is imported directly. Then, the built-in Grad–Shafranov solver in JOREK is used to obtain a discretely accurate force balance in the initial conditions. These JOREK-calculated equilibria can also be used to calculate the Alfvén continua with the HELENA and CASTOR codes. The MHD temperature is not set to an experimental profile but to reproduce the equilibrium pressure, and thus, temperature dependencies in viscosity and resistivity terms are not desired. These dependencies are switched off, leading to a fixed value of the viscosity and resistivity.

A highly anisotropic, realistic distribution function is used, obtained from Ref. 27. This distribution was calculated with the NUBEAM¹⁷ NBI code and is shown in Fig. 12. This distribution function was found to relax much less in both velocity and real space, as can be seen in Fig. 13 

As the NBI beam was injected co-current, the particles are almost solely counterpassing (i.e., not trapped and propagating in the direction opposite to the magnetic field). A slight complication is that the EPs (in a stationary state like used here) are not uniformly distributed along the poloidal angles, such that in principle the EP density is not solely a function of ψ. However, this is not taken into account at present, and the flux-surface averaged density is used to obtain the bulk mass density. The EP marker particle weight is then set such that the total amount of particles is equal to the integrated flux-surface averaged density from the IDA diagnostics.

20 × 10⁶ particles were used for all timepoints. 51 radial elements were used with 64 poloidal elements. Timesteps of $≈ 4 τ A$ were used for the MHD fluid and $2 × 10 − 10$ s for the EPs.

Before the results are shown, it is worthwhile to shortly discuss the evolution of the discharge. The q-profile evolution obtained from the IDA data is shown in Fig. 14. The discharge is characterized by a reverse shear q-profile, which decreases in time. In the IDA diagnostic, the q-profile is also quite flat in the core (s < 0.4) for $t ≥ 0.75$⁠. The EP pressure (and density) increases, but the total (EP + bulk) pressure stays consistent throughout the discharge. The bulk mass density slowly decreases in time.

The growth rates, frequency, total amount of EPs, and type of mode for the most unstable modes at the considered times are given in Table I. The frequency spectra and poloidal harmonic structure for the most unstable modes are shown in Fig. 15. It is clear that three regimes can be identified: while q > 2.5 (⁠ $t = 0.6 , 0.65$⁠), the m = 3 dominated EPMs are present. When q crosses the 2.5 mark, a core TAE emerges (t = 0.7). Finally, when q < 2.5 for s < 0.6 (⁠ $t ≥ 0.75$⁠), the q-profile is flat in the core and m = 2 dominated EPMs or reversed shear Alfvén eigenmodes (RSAEs) emerge. Superimposing the most unstable modes in the simulations on the experimental spectrogram yields Fig. 16. Good agreement is found for $t ≤ 0.7$ and $t ≥ 0.9$⁠, but the RSAEs at ${ 0.75 , 0.8 , 0.84 }$ are not present in the spectrogram. To assess the cause for this discrepancy, consider the uncertainty in the core q-profile. The error bars on the q-profile lead to uncertainties in the continua as shown in Fig. 17. Since the observed modes are all consistent with the Alfvén spectra, this rather large uncertainty in the Alfvén spectra could explain the differences. Furthermore, a modification of the core q-profile such that it is monotonic for s < 0.4 could lead to these low-shear EPMs or RSAEs at $t ≥ 0.75$ s no longer being the most unstable modes, such that possibly modes could emerge that are visible in the experimental spectrogram. Also, non-linear effects modifying the EP distribution function could play a role.

Confining energetic particles is crucial for sustaining a burning fusion reactor, but EP-driven instabilities can prove a threat to this confinement. Understanding these instabilities is important for developing next-generation devices and operational scenarios. In this work, the non-linear MHD code JOREK has been extended by anisotropic pressure coupling of EPs to the full and reduced MHD models such that it can simulate these EP instabilities using a wide variety of possible bulk physics models, powerful numerics, and grids up to the first wall.

The JOREK model for EP simulations was explained, including a description of a versatile projection diagnostic that has been developed to investigate EP dynamics in more detail in velocity space.

The JOREK code was then benchmarked against other codes using two cases: a TAE benchmark in simple geometry and an EPM benchmark in realistic geometry with high EP pressure. Although other codes used a different model for the EPs, which does not include the relaxation of the initial Maxwellian distribution, good results have been obtained in terms of mode structure, frequency, growth rates, and phase-space resonances.

The code is then applied to a high EP pressure discharge in Sec. V using experimental plasma geometry and parameters and a realistic distribution function for the EPs. Three regimes could be identified from these simulations and the calculated Alfvén continua, namely, an m = 3 EPM regime, a m = 2, 3 TAE, and finally a RSAE regime. Observed frequencies do match well with the Alfvén continua calculated from the experimental equilibria. Some of the obtained frequencies are present in the experimental spectrogram, but some are not. The impact of uncertainties in the q-profile was shown to be large for the core Alfvén continuum, providing a possible explanation for the RSAEs that are not present in the experimental spectrogram.

In the future, codes could be benchmarked in the linear and non-linear regime using realistic distribution functions and realistic or experimental plasma profiles in order to get closer to experimentally relevant scenarios. Further investigations into the AUG discharge #31213 (possibly using multiple codes) would be of interest as well, since this discharge was designed to mimic reactor-relevant EP conditions.

This work is only a first step regarding EP instability studies using JOREK, and not all features of JOREK have been used. For example, collisions can be used to self-consistently evolve a EP distribution function with sources and sinks, requiring, however, particular attention to energy conservation. Gyro-center particles can be used to improve the performance after detailed comparisons to full-orbit particles. Fully consistent equilibria can be generated by evolving the equilibrium and the EPs axisymmetrically. The projection diagnostic can be used to investigate non-linear behavior in depth. The full MHD capabilities could be used for investigating fishbone-type instabilities in realistic geometry.

We thank Philipp Lauber for providing the spectrograms of AUG shot No. 31213. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200—EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

The authors have no conflicts to disclose.

Timo J Bogaarts: Conceptualization (equal); Data curation (lead); Formal analysis (lead); Investigation (lead); Methodology (lead); Software (lead); Validation (lead); Visualization (lead); Writing – original draft (lead). Matthias Hoelzl: Conceptualization (equal); Data curation (supporting); Formal analysis (supporting); Funding acquisition (lead); Investigation (supporting); Methodology (supporting); Project administration (lead); Resources (lead); Software (supporting); Supervision (lead); Validation (supporting); Visualization (supporting); Writing – review & editing (equal). Guido T.A. Huijsmans: Conceptualization (equal); Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Methodology (supporting); Software (supporting); Supervision (equal); Validation (supporting); Visualization (supporting); Writing – review & editing (equal). Xin Wang: Data curation (supporting); Formal analysis (supporting); Investigation (supporting); Validation (supporting); Visualization (supporting); Writing – review & editing (supporting). The JOREK Team: Funding acquisition (equal); Software (equal).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

As mentioned in Sec. III, the oscillating $( v ∥ , μ )$ can be troubling for a resonance visualization diagnostic. This is not a problem for resonance visualization for quantities that are fully conserved during full-orbit motion such as the energy E and the toroidal canonical momentum $P ϕ$⁠. However, to compare with other codes and theory, it is useful to have the option to visualize resonances in terms of $( v ∥ , μ )$ as well. Therefore, in this appendix, a method is introduced to eliminate this oscillation. This method is then illustrated with results from the ITPA benchmark of Sec. IV A. During the gyro-motion, the energy of a full-orbit particle varies with a larger oscillation and a smaller trend, shown in Fig. 18. As the energy oscillation is of higher magnitude than the trend in the energy, the diagnostic will hide this trend in the oscillation, leading to the power exchange diagnostic producing dipole-like results in both (R, Z) space and $( μ , v ∥ )$ space, shown in Fig. 19. Although these plots can still show the resonances (and might be interesting results in their own right), it would be more useful if the diagnostic only produces the net power exchange, without the oscillating contributions.

To convert a full-orbit particle with oscillating $v ∥$ and μ to a gyro-kinetic particle with fixed μ and non-oscillating $v ∥$⁠, these quantities can be modeled as a linear function of time modulated with an oscillation at the gyro-frequency (within short enough timescales). The function to be fitted is then

Considering a set of values ${ μ , v ∥ } i = y i$ at time t_i, the difference (or residual) r_i between the values from the full orbit quantity y_i and the function $f ( t i , α , β , γ , δ , ω )$ at that time is

The sum of the residual squares is

such that the partial derivative of the sum of squares by the parameter $β j ∈ { α , β , γ , δ , ω }$ is

These can be calculated easily as

Without the last equation, it is a linear problem. As the gyrofrequency $ω = q B / m$ only depends on the value of B during its full orbit, the gyrofrequency can be estimated by taking the mean value of B during the orbit.

The expression of the function can then be rewritten as

with

Assuming N observations, the observation matrix $X i j : = ϕ j ( x i )$ is then

The quantity $X T X = X j i X i k$ can then be calculated as

The inverse of this can be calculated using Cramer's rule. The quantity $X T Y$⁠, where $Y = y i$⁠, is

The solution for the parameter vector $β$⁠, which minimizes the sum of squares of residuals can now be obtained by³⁸ 

Then, the linear part of this gyro-motion can be used to project the power exchange at the $( μ , v ∥ )$ location of the non-oscillating gyro-center particle corresponding to the full-orbit particle (either at a few points if the $v ∥$ trend is steep, or only at the average, which proved to be sufficient in the cases considered in this work). The effect of this procedure is shown in Fig. 20, where now only actual resonances are visible.

A final remark is that it is very useful to sum (or average) the power exchange over many fluid timesteps, as a single fluid time step only contains about five gyro-motions, which is insufficient for the power exchange diagnostic.