An IMAS-integrated workflow for energetic particle stability

The confinement of energetic particles (EPs) generated by fusion reactions and external heating methods is crucial for the performance of future fusion devices. However, EP transport can occur due to their interaction with electromagnetic perturbations, affecting heating efficiency and overall performance. Robust reduced models are needed to analyze stability and transport. This paper presents an automated IMAS-based workflow for analyzing the time-dependent stability of EP-driven modes, focusing on the linear properties of Toroidal Alfvén Eigenmodes (TAEs) in general tokamak geometry. The workflow utilizes efficient computational methods and reduced models to deliver fast and reproducible results. A demonstration of the workflow’s effectiveness was performed, identifying key linear properties of TAEs in various simulated ITER scenarios. This approach represents a critical step toward developing tools for analyzing EP transport and optimizing the performance of future fusion reactors.


Introduction
Energetic particles (EPs) generated by fusion reactions and by external heating methods such as neutral beam injection need to be well confined in future fusion devices such as ITER or DEMO. However, the EPs' resonant and non-resonant interaction with electromagnetic perturbations can lead to EP transport which affects the heating efficiency and thus the overall performance of the burning plasma which is of concern for future reactors. Previous analyses show that in order to be able to determine the stability limits of weakly damped Alfvénic perturbations and to model the related transport, a * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. non-linear, global, kinetic treatment is necessary [1,2]. But in order to study a whole scenario including ramp-up, flattop and ramp-down that e.g. at ITER is foreseen to last >500 s, one needs to find a compromise between computational fidelity and speed. Comprehensive codes with high fidelity (e.g. ORB5,-global electromagnetic gyrokinetic particle-incell code [3,4], GTC [5] or GENE [6] can be used only for relatively short timescales (several ms), due to their extensive computational cost. Non-linear (gyro-) fluid models [7] can cover longer timescales, but besides the lack of kinetic effects, the deep non-linear evolution poses still some unresolved challenges. Also, direct coupling schemes of high-fidelity simulations and transport codes are presently developed [8], which however are also extremely costly. Therefore, robust reduced models for the analysis of stability and transport are needed as an additional tool. Presently, the following implementations based on different physics models are available: the critical gradient model [9], the kick model [10], the 'resonance broadening quasi-linear model' (RBQ) [11], and the phase space zonal structure model (PSZS) [12]. Although they are intended to be faster than comprehensive codes, significant effort is still required to identify the fundamental components (such as linear mode features including frequency, damping/growth rate, mode structure, and associated saturation rules) for the reduced models. In addition, to perform multiple linear computations at different time slices of a given scenario (ramp-up, flat-top and ramp-down) a high degree of automatization is required.
In this work, the first automated, time-dependent IMAS (Integrated Modeling & Analysis Suite [13]-partially based on earlier developments [14]) workflow (EP-Stability-WF) is presented and used to provide insight into the stability of Toroidal Alfvén Eigenmodes (TAEs) in different phases of various scenarios. It is based on the equilibrium code HELENA [15] and different hierarchical models of the linear gyrokinetic eigenvalue code LIGKA [16] that are connected via a centralized IMAS database. A complex workflow (EP-Stability-WF [17]) was created to connect the different stages of the analysis with the database to ensure the reproducibility and consistent interpretation of the different analysis steps. With this tool, both time-dependent predictive scenario simulations and experimental data can be analyzed.
To demonstrate the capabilities of this workflow, several ITER scenarios for deuterium-tritium plasmas are analyzed: a time-independent scenario generated by the ASTRA [18] transport code, a time-dependent predictive scenario generated by the METIS [19,20] [22] transport workflow which combines a free boundary equilibrium evolution code with the 1.5D core/2D SOL code. Here time-dependent means that an automated linear stability analysis is performed for different time points of a simulated ITER scenario. These time points are sufficiently separated in time (typically 0.5-5 s), i.e. it is assumed that the equilibrium evolution time is much slower than the Alfvén time scales connected to the linear and non-linear TAE physics.
The paper is organized as follows: section 2 contains a technical description of the models and numerical tools, section 3 gives an overview of the scenarios and parameters chosen for the computation. The presentation of the results (section 4) is split into three subsections: the evolution of damped TAEs without EPs, what happens to them when EPs are introduced, and a convergence test for a whole scenario. In the end, section 5 presents the conclusions and outlook of the paper.

Numerical tools and physical model
The aim of the WF is to perform an automated linear stability analysis on different time slices of a simulated scenario or reconstructed experimental equilibrium. This is achieved by connecting various numerical tools with the data infrastructure (IMAS), facilitating the retrieval/saving of data from/to the database (DB) through a series of IDS files and fast configuration of numerical tools via XML files. As described in [13] IMAS relies on IDSs (integrated data structures) that combine a set of well-defined physical quantities into the complete description of a complex physical object such as an equilibrium or the linear stability information related to this equilibrium. In figure 1 one can see how the IDSs are being passed throughout the workflow between the codes and the DB. The inputs and outputs of each actor are well defined in terms of IDSs and the static configuration parameters that define the respective simulation are given in XML. In order to facilitate the interaction with the code and the underlying infrastructure, a graphical user interface (GUI) is also available, this is made possible by the standardization of code-specific parameters in XML files. Figure 1 explains how the individual actors are combined on a conceptual and technical level: (i) Take the equilibrium_ids from the transport code (given by ASTRA, METIS or DINA-JINTRAC) and pass it to HELENA or CHEASE [23] to perform high-resolution MHD equilibrium calculations. HELENA is used in the workflow as a coordinate transformation tool as LIGKA requires a certain representation of the equilibrium: straight field line coordinates, metric tensor representation and, mapping to the R-Z grid. As output, they will produce another equilibrium_ids/0. From here on, * /0 means occurrence 0, one can have several instances of the same IDS in one workflow/database. (ii) Together with equilibrium_ids/0, core_profiles_ids is being taken from the transport code and given as inputs for model 5 of LIGKA (see next subsection). Model 5 performs analytical estimates of TAE (and other AEs such as e.g. BAEs, EAEs, RSAEs, as explained later in the LIGKA model) and creates as an output mhd_linear_ids/0 which in turn it is saved in the local database.
(iii) Based on this result from model 5, the next model of LIGKA hierarchy-model 4-can be run: as input, it uses equilibrium_ids/0, core_profiles_ids/0 and mhd_linear_ids/0 written previously by model 5 and creates a new occurrence of the mhd_linear_ids/1. Model 4 returns a local estimate of the linear AE properties, i.e. it solves the local kinetic dispersion relation, using the approximations described in detail in the following papers [1,24,25]. (iv) Finally, models 1 and 2 of LIGKA, the standard global solver can be started. As necessary input IDSs, we have equilibrium_ids/0, core_profiles_ids/0 and mhd_linear_ids/1 (output from model 4). Again, it creates a new occurrence of mhd_linear_ids/2 or mhd_linear_ids/6.
This WF is not only driven by the need to analyze the stability problem with different levels of complexity and speed that is reflected by the different LIGKA models but also, each LIGKA model profits from the results of the previously run model: the local dispersion relation solver (model 4) needs a good starting guess for the complex contour integration in the frequency plane that is provided by the analytical estimate (model 5). The global solver (model 1/2) in turn uses model 4 results to set up an efficient 'antenna' scan based on the knowledge of the kinetic continuum structure. Consistent mapping of the results of each model into the DB allows us to construct a consistent hierarchy of results and eventually perform uncertainty quantification for each model. The modular setting and the usage of standard IMAS data structures facilitates also the usage of various other codes than LIGKA under the same umbrella in the future. LIGKA models could be replaced by other codes that perform different computations but arrive at the same physics result or, more specifically, a sufficiently populated data structure (mhd_linear) that can be used as input for the next model hierarchy (i.e. the output of replacing code includes everything necessary to run the next LIGKA model).
In practice, the analytical (model 5) and local solvers (model 4) are used to obtain a fast overview of the scenario before attempting global, more expensive runs. The global solver (model 1) can then be used to validate-and in most cases-to improve the results obtained by the local runs.
By solving the linearized gyrokinetic equations LIGKA outputs eigenvalues and eigenfunctions such as frequency, damping and mode structures (electrostatic and electromagnetic components, and derived from that also the parallel electric field). Models of LIGKA used in this work are: (i) (1 s/mode) Model 5: local analytical estimates of various AEs properties: frequency, estimated mode structure, rational surface, next and previous gap information based on the analytical formulae of [26,27] (in case we are dealing with AEs related to gaps in the shear Alfvén continuum). (ii) (10 s/mode) Model 4: based on model 5 results, the local analytical dispersion relation for one AE (n, m-pair) is calculated. The output is the complex frequency of the maximum or minimum of the continuum closest to the AE's  figure 2 (last)). (iv) (2-5 min/mode) Model 2: uses the same solver as model 1, but model 2 tries to accurately find the phase jump during the frequency scan, such that a more accurate growth/damping rate can be obtained, or the runtime can be reduced by starting the scan close to the AE's frequency as determined by model 1.
In this paper, we focus exclusively on TAEs for simplicity, although also other types of AEs such as RSAEs (reversed shear AEs), BAEs (beta-induced AEs) and EAEs (ellipticity induces AEs) are implemented. TAE are the prototype of all Alfvénic Eigenmodes (AE) [28]. They have their frequency inside the shear Alfvén continuum gap and their mode structure is given by coupling of two counterpropagating regular shear Alfvén waves. TAEs can be excited through resonance by fast particles and fusion products because their parallel group velocity in toroidal devices is typical of the same order as the parallel particle velocity. It is well known that these modes can lead to significant radial particle transport and EP losses [29][30][31]. This workflow is also equipped to work with other numerical tools such as HAGIS [32] 1 and 2 for the non-linear stability of EPs, but this surpasses the scope of this work.

Scenario description
For the purpose of demonstrating the capabilities of the newly developed workflow, three scenarios were chosen from the many available in the ITER-IMAS database. The scenarios analyzed in this paper are the baseline deuterium-tritium predictive ITER scenarios (generated by different transport codes). The transport data was obtained from running the ASTRA, METIS or DINA-JINTRAC workflows, transport packages that are available at ITER through the Integrated Modeling & Analysis Suite (IMAS).
In order to give a clear overview, this chapter is separated into three sub-chapters, each describing one scenario.

ASTRA-131025/34
This specific ASTRA-generated scenario is time-independent, i.e. only one time slice at 208 s in the flat top phase has been exported into IMAS. ASTRA, similar to other transport codes, is a tokamak simulation tool based on 2D equilibria and a set of 1D diffusions equations for particle and heat transport [18]. For this particular scenario, a semi-empirical scaling law for the transport and a sawtooth mixing model were used [33]. In this scenario the plasma composition is a 50:50 deuteriumtritium mix.
The electron, deuterium and tritium profiles can be seen in figure 3. Electrons have a temperature on the axis of T e,0 ≈ 23 keV, and T e,p ≈ 4 keV at the top of the pedestal. For D and T, the on-axis temperature T D/T,0 ≈ 22 keV with the temperatures at the pedestal the same as for electrons. We have an almost flat density profile for electrons with n e,0 ≈ 1.2 × 10 20 m −3 . For the ion species considered in this scenario (D and T) the on-axis density is n D/T,0 ≈ 0.4 × 10 20 m −3 . The profiles peak slightly in the pedestal region n D/T,p ≈ 0.5 × 10 20 m −3 due to the presence of fast ions and helium ash in the core. However, since no information about the EP fractions of density and pressure is present in the output of the transport code, no fast particles were included in the analysis of this scenario. Instead, we will focus on the damping mechanisms.
In figure 4 the safety factor profile (with a zoom showing the values in the inner part of the tokamak and the minimum value of the profile) can be seen. q 0 = 1.00 in this case, meaning that in principle due to q TAE = (m + 0.5)/n all TAEs with n ⩾ 1 could be present. Also, the q-profile is not strictly increasing. This can lead to multiple TAE surfaces, as shown by the blue line that indicates q TAE = (14 + 0.5)/14 = 1.036 for the (n, m) pair (14,14).

DINA-JINTRAC-134173/106
The second scenario taken into consideration is a timedependent D-T scenario produced with the DINA-JINTRAC  transport workflow. It assesses the entire scenario evolution from the early ramp-up phase (from X-point formation) until the late ramp-down phase (to X-point-limiter transition) by means of integrated simulations including core (1.5D), edge (2D), and SOL transport with time-dependent freeboundary plasma geometry and pedestal pressure determined by continuous self-consistent edge MHD stability analysis. Neoclassical transport is modeled with NCLASS [34] and turbulent transport is described by the L-mode version of the Bohm-gyroBohm model in L-mode [35] or by the GLF23 [36] model in H-mode [37].
In figure 5 the fusion power peaks at 150 s in the simulation and then slowly decays over the rest of the scenario. Unfortunately, at the time of writing this paper no fast particle information was available in the transport code output, and thus, no fast particles (EPs) were included in the analysis of this scenario. In

METIS-130012/02
The third scenario is a time-dependent scenario generated by METIS. To achieve faster speeds (real-time and faster) the physics model in METIS is simplified with respect to the DINA-JINTRAC workflow. It combines 0D scaling-law normalized heat and particle transport modeled with 1D current diffusion modeling and 2D equilibria. The details are described in reference [19]. It is again a D-T plasma, the baseline ITER scenario for Q = 10, and thus, one of the most important scenarios for fusion at ITER.
In figure 7 the evolution of heating and fusion powers for the METIS scenario are shown. Of special interest for our study is the P FUS (purple line). It peaks at approximately 90  s with P FUS ≈ 130 MW. This evolution allows us to test and validate the EP part of the EP-WF. Figure 8 shows, like before, the density and temperature profiles this time at t ≈ 90 s i.e. at the peak of P FUS . For electrons, the temperatures on the axis and at the pedestal peak at T e,0 ≈ 27 keV, and T e,p ≈ 4 keV. For background ions (D and T in this case) the temperature profile is the same with T D/T,0 ≈ 25 keV at the axis, and the same as the electrons at the pedestal. The electron and background ions densities are mostly flat, peaking at n e,0 ≈ 1 × 10 20 m −3 for electrons and n D/T,0 ≈ 0.45 × 10 20 m −3 for background ions. In this scenario, we have alpha particles from the fusion process. Their density is also shown in the same figure (100 ×), peaking on axis at n He4,0 ≈ 0.38 × 10 18 m −3 .
In figure 9 we have the safety profile at t = 201.96 s for this scenario and the minimum value of this profile is q min = 1.06 depicted by a green dot in the figure.

ASTRA-131025/34
In this time-independent scenario, we perform a scan over a wide array of TAEs (n = 1, . . . 35, and m = n . . . n + 4) for the given time slice in order to check if all WF elements (including the kinetic multi-species part) are behaving as expected. For this case, we include only thermal background particles, i.e. electrons and a 50:50 D-T mix satisfying quasineutrality at each radial slice.
In figure 10 the radial location of the modes over the whole radial domain is presented. Also, the q ≈ 1 rational surface is marked as a red line. LIGKA model 5 uses as an analytical estimate the TAE condition q = m+ 1 2 n . Also, more sophisticated formulae based on [26] are implemented, taking into account elongation and β corrections (not shown here). However, it turns out that the simplest expression gives a good estimate for the center of the TAE gap (that is needed for the subsequent analysis) in all geometries and mode number ranges investigated so far. Part of the modes can be typically discarded in the following analysis (i.e. models 4 and 1) based on their position relative to the EP gradient: the alpha particle gradient is significant only in the region 0.3 . . . 0.65, thus the modes that do not fall in this range were excluded since they are expected to be linearly stable. Unfortunately, as explained above, details about EP pressures and densities are not available for this simulation, although they have been included in the calculation. The radial range of interest can also be conveniently defined via the input files (XML data structure, editable directly, or via a GUI), allowing the user to have flexible control over the considered mode spectrum.
For all selected toroidal mode numbers the shear Alfvén wave continuum (SAW) can be calculated (model 6), as shown in figure 11 for n = 14, 15, 16. All branches of TAEs in the radial range discussed above are added. It can be seen that the simple analytical guesses (model 5) lie typically in the TAE gap, as expected. The two solutions for n = 14 are shown (see figure 4), separated in frequency by the changes in the local Alfvén velocity. Now we calculate the global solutions with model 1. In figure 12 (round dots) the normalized mode structures (electrostatic potential) for the two modes n = 14, m = (14,15) found at different radial locations are displayed (as explained before). By taking these structures and adding them to the continuum (see figure 12) one can see that they intersect the reduced SAW continuum. This is not surprising since the  global solutions include various kinetic effects not accounted for in the reduced MHD model, in particular, the diamagnetic upshift of the continua in regions of high plasma pressure [24]. In figure 12 last subfigure, the local kinetic continuum-model 3-is added. For a detailed description of this local solution of the kinetic dispersion function see [25]. From the figure one can see that using the kinetic model consistently in both local and global calculations, the TAEs are, as expected, situated in the kinetic TAE gap, avoiding local continuum damping. It becomes apparent from this discussion that the kinetic corrections for these ITER scenarios and mode number ranges are significant and cannot be neglected as often done in presentday experiments where both the plasma beta and the excited TAE mode numbers are typically much lower.

DINA-JINTRAC-134173/106
In this subsection, the previously shown capabilities of the workflow are used to analyze the ramp-up phase of a DINA-JINTRAC scenario simulation (134173/106).
Similarly to what has been shown in the last section, a run was performed at one time point (t = 97.93 s), using model 5. For this run, TAEs with n = 1 to n = 35 with m = n − 2 to m = n + 2 as main poloidal harmonics were considered. Figure 13 shows the radial position of these modes together with the q ≈ 1 rational surface. As expected, modes with low toroidal mode numbers are shifted away from q = 1 towards the core, when the poloidal mode numbers are lower (m = n − 2) or towards the edge when m = n + 2 according to q TAE = m+ 1 2 n . Next, a specific TAE with (n, m) = (25, 24) was tracked during the ramp-up phase of 134173/106. This mode was chosen based on the fact that q min is slightly below 1 in most of the time points, meaning that the TAE with (n, m = n − 1) is present throughout the whole ramp-up and flat top phase. Figure 14 shows the frequency evolution. Plotted are the different models used by LIGKA. Starting with model 5 (blue line), continuing with the local solver, model 4 (orange line) and finally the global solver, model 1 (green line). Performing convergence tests as described above we found that the range of poloidal harmonics [−4,+32] is necessary to study this particular part of the scenario. Features from the computations regarding the frequency of the mode can be seen here as well in addition to the general Alfvén scaling trend: due to  the rising density, the frequency decreases due to f TAE ∝ 1 √ n . Model 5 results in the highest frequency (middle of the gap), model 4 the closest continuum (downwards, lower frequency) and finally model 1 determines the exact location of the TAE within the gap via a global analysis. Figure 15 shows the damping rate of the TAE, during the same period of time (ramp-up). Again, a considerable difference between models can be noticed, with model 1 being the most reliable, showing clearly the local damping rate underestimates the global one considerably, as explained above. It can be seen that the chosen TAE is strongly damped in the flat top, whereas during the ramp-up phase, the damping rate is much smaller. Clearly, this motivates future detailed studies including the drive from the neutral beams and the alpha particle drive that were not included in this calculation.
For completeness, the mode structure is presented in figure 16. One can see that despite the high toroidal mode number, the mode is quite global due to the flat shear around s = 0.5. However, as figure 17 shows, the continuum damping is negligible since the mode touches the continuum at s = 0.9 where all poloidal harmonics with significant relative amplitude have decayed. Thus, just other non-local damping  effects (ion/electron/radiative) are contributing via the wide mode structure to the overall mode damping.

METIS-130012/02
Using the METIS scenario, the workflow demonstrates all the capabilities of a complete and automatic tool for linear stability inside a tokamak. Starting from an overview of the TAE locations and frequencies (via model 5) to a more indepth analysis via the local (model 4) and then a global solver (model 1/2) was performed. First alpha particles are excluded to determine the damping rates, and then the effects of alpha particles are added. Also, a convergence test was performed in order to investigate and optimize the computational resources for the various analysis steps.
Due to the ability of model 5 to perform basic analytical AE calculations in a short amount of time (<1 s-mode-time step), a scan over the whole METIS scenario (106-time points) was performed taking n = 1 . . . 35 and m = n + 1 . . . m = n + 5 for a total of 12014 modes, averaging about 137 modes per time point. In figure 18 a plot of all these modes' positions over time is given. Pointed in the graph is also the q ≈ 1 value for each time point. The slight inversion of the q-profile around

√
ψ N ≈ 0.5 in the flat top part of the scenario (>90 s) leads to two different TAE branches with the same mode numbers. Also, after (>90 s) no TAEs are found in the core (s < 0.4). This is a consequence of the extremely small magnetic shear in the core assumed by the METIS-given equilibrium.
In figure 19 the reduced SAW continuum and the kinetic continuum for n = 10 are plotted. The results for one TAE (m = 11, 12) using the model hierarchy are added. As in the previous case, the different models give slightly different TAE frequencies (colored dots in figure 19), in line with the discussion above.
Since the SAW continuum in figure 19 shows an almost open TAE gap, it is investigated how many poloidal harmonics are needed for the convergence of this mid-n (n = 10) TAE. This test is also of significant importance for solver validation inside LIGKA and the EP-WF. Several runs were done by looking at the same mode (n, m) = (10, 11) using a different number of poloidal harmonics.  In figure 20 the mode frequencies calculated using different numbers of poloidal harmonics are shown. The frequency differences between the different runs are not significant. During the ramp-up (50 <= t <= 100 s) and ramp-down (300 <= t <= 400 s), the variation in frequency is determined by the evolution of the density, temperature and safety profile. Note that the time step size for equilibria and profiles in this scenario is ≈ 3 s. For the flat-top phase of the scenario with stable profiles, almost no changes in the TAE frequency are observed.
In figure 21 the damping rate of the same mode is plotted, similarly to figure 20. Figure 22 shows the mode structure of this mode for the same number of harmonics. In this case, if the number of poloidal harmonics is varied, significant differences are observed. As it can be seen, the number of poloidal harmonics for convergence depends on the equilibrium, i.e. particularly at the end of the ramp-up phase nonlocal continuum damping towards the plasma edge requires many poloidal harmonics to be included. This difference is particularly pronounced in the most unstable (profile-wise) regions of the scenario: ramp-up and ramp-down. As depicted in the graph with blue we use 21 harmonics (from m = n − 4 to m = n + 16). One can see that as more harmonics are added, the results become more and more weakly damped. In the flat top region of the scenario, the difference is not that significant ±0.5%, whereas, in the ramp-up phase, the damping strongly depends on the number of poloidal harmonics (≈ − 25%). Adding more harmonics leads eventually to a converged damping rate, especially in the ramp-up phase (90-110 s). Obviously, the mode number n = 10 was chosen as the most challenging case for this scenario, since the global TAE touches the continuum close to the plasma edge, explaining the slow convergence behavior: as the number of poloidal harmonics is increased, the interaction with the continuum is radially further and further separated from the main TAE location, decreasing the damping rate (see figure 21). The continuum intersection moves to the edge and thus smaller and smaller poloidal harmonics interact with the (not yet resolved) continuum until it is fully resolved (−4/+28 and −4/+32). In order to demonstrate convergence three extra runs were performed: two with the higher number of harmonics (+36 and +40) and one for a higher radial resolution (512 radial grid points). This can be seen for 1-time point (t = 250 s) in the flat top region of the scenario (figure 23). Based on this, one can conclude that −4/+28 and −4/+32 harmonics are enough to achieve convergence for the particularly demanding TAE.
Next, alpha particles coming from the fusion reactions were added, as given by the scenario simulation. For simplicity, an effective hot Maxwellian was chosen, consistent with the EP pressure and density as given by the transport code. In figures 24 and 25, the frequencies and damping rates for the same n = 10 and (m = 11, 12) TAE for the whole time domain (t = 45-360 s) are given. Based on the previous analysis, the poloidal harmonics included in this simulation was −4/+32 which was found to be the necessary range needed for convergence. Again, the various models are run and compared (dashed lines: no alphas, solid: alphas). The blue vertical line on the graphs indicates the moment when a significant population of alpha particles starts to be present in the plasma due to fusion processes, As expected, the frequency did not shift by a significant margin (figure 24) between runs with alphas or without.    Figure 25 shows the damping rate for both cases, with and without alpha particles taken into consideration. Using model 4 (local solver-orange line) one can see that there is a small difference between the cases with and without EPs. Moreover, this local model suggests that in the ramp-down, the TAE becomes unstable in the presence of alphas. Using the global solver (model 1/2-green line) one can see that this is not the case, due to the non-local damping mechanisms. At t ≈ 90 s when alphas are turned on, the damping rates start to differ. However, while the alphas do reduce damping, the TAE does not become unstable throughout the whole scenario. The same is true for other n's (tested for n = 10 . . . 30) showing that this extremely (unrealistically) flat q-profile results in a decoupling of TAE locations and steep EP gradient regions, leading to a stable TAE spectrum. Small equilibrium deviations can change this, particularly assuming different models for density peaking. Using the established workflow sensitivity scans based on slightly differing assumptions in the scenario modeling process can be systematically investigated in the future.

Summary and outlook
Our study aimed to develop an automated workflow for timedependent EP transport analysis. The focus of this paper was on the linear properties of TAEs in a torus, specifically in ITER predictive scenarios. A demonstration of the effectiveness of this workflow in identifying the key linear properties of modes such as damping rates, frequency, location and structure was performed.
The workflow together with the associated physics codes, have demonstrated the ability to generate hierarchical physics results. By utilizing efficient computational methods, reduced models and ensuring that all code and data are versioncontrolled (via IMAS) and documented we were able to produce physically consistent and reproducible results. As an overview, table 1 contains detailed information about the runtime of the test cases discussed above is presented.
The successful application of this approach represents a fundamental and challenging first step toward developing more sophisticated tools for analyzing EP transport. These types of tools will be critical for optimizing the performance of fusion reactors, such as ITER. This work will contribute directly to various reduced EP transport models. Table 1. Run times for different cases presented in this paper. The computational resources are also shown in terms of processes being used for each run. The total run time for the equilibrium code at each timeslice is also included.

Appendix. Reproducibility of the results
The codes used to generate the results in this paper are controlled via a release process. The codes are also available as modules on the ITER Linux cluster ('SDCC').
In this paper, the following versions of the codes were used, with the Data Dictionary (DD) version 3.