Alpha particle driven Alfv\'enic instabilities in ITER post-disruption plasmas

Fusion-born alpha particles in ITER disruption simulations are investigated as a possible drive of Alfv\'enic instabilities. The ability of these waves to expel runaway electron (RE) seed particles is explored in the pursuit of a passive, inherent RE mitigation scenario. The spatiotemporal evolution of the alpha particle distribution during the disruption is calculated using the linearized Fokker-Planck solver CODION coupled to a fluid disruption simulation. These simulations are done in the limit of no alpha particle transport during the thermal quench, which can be seen as a most pessimistic situation where there is also no RE seed transport. Under these assumptions, the radial anisotropy of the resulting alpha population provides free energy to drive Alfv\'enic modes during the quench phase of the disruption. We use the linear gyrokinetic magnetohydrodynamic code LIGKA to calculate the Alfv\'en spectrum and find that the equilibrium is capable of sustaining a wide range of modes. The self-consistent evolution of the mode amplitudes and the alpha distribution is calculated utilizing the wave-particle interaction tool HAGIS. Intermediate mode number ($n=7-15,~22-26$) Toroidal Alfv\'en Eigenmodes (TAEs) are shown to saturate at an amplitude of up to $\delta B /B \approx 0.1$\% in the spatial regimes crucial for RE seed formation. We find that the mode amplitudes are predicted to be sufficiently large to permit the possibility of significant radial transport of runaway electrons.


Introduction
The subject of runaway electron (RE) mitigation is of crucial importance to the success of reactor-relevant tokamaks such as ITER [1][2][3][4][5]. Generation of REs is most concerning during disruptions, as the avalanche mechanism [6,7] is expected to convert a significant portion of the plasma current into runaway current on a large tokamak [8][9][10]. A multimegaampere runaway beam has the potential to inflict significant damage to plasmafacing components [11,12]. This paper discusses a phenomenon which could potentially act to passively mitigate runaway electron generation at ITER.
The interaction of runaways with plasma waves has been investigated both in theory [13][14][15][16][17][18][19][20] and through observations in multiple tokamaks [21][22][23][24][25][26][27]. Several tokamaks (such as DIII-D or TEXTOR ) have reported runaway suppression in correlation with increased wave activity in the current quench of the disruption or the plateau phase of the runaway beam. Wave activities in a tokamak plasma can lead to a variety of instabilities with different effects on particle confinement to which runaway electronsdue to their high velocity -can be extremely susceptible to.
In fusion, the umbrella term "shear Alfvén wave" collects an important type of transverse, electromagnetic plasma waves characterized by their low (Alfvénic) frequency range and propagation along the magnetic field. Frequency gaps in the continuum damping allow the existence of these Alfvén Eigenmodes (AEs). Those gaps can occur through an extremum in the safety-profile (Reversed Shear AEs [28], Global AEs [29] and Beta-induced AEs [30]) or through geometric coupling of two poloidal harmonics (Toridicity-induced AEs [31][32][33], Ellipticity-induced AEs [34]). Compressional AEs (CAEs, [35,36]) are high frequency kinetic instabilities with both a perpendicular and a parallel component.
Instabilities in the Alfvénic frequency range are routinely driven in fusion plasmas by energetic ions. However, in the case of such modes observed during the quench or post-disruption, the source of the drive is less obvious. Former analytical work on runaway ions [37] was inconclusive and it was later deduced [38] that ion runaway formation even in reactor-sized tokamaks is unlikely at disruption time scales. Spontaneous ion acceleration caused by internal magnetic reconnection events in the MAST tokamak were observed [39]. Nevertheless, the experimental observation of these modes necessitates a drive to exist. A recent publication [40] identified runaway electrons as a possible drive for CAEs (or possibly GAEs) on the DIII-D tokamak in the context of a Massive Gas Injection (MGI) triggered disruption and the consecutive suppression of a runaway plateau formation.
For the ITER 15 MA scenario it is predicted [41,42] that marginally unstable modes can already be present in the quiescent burning phase. Decisive for the presence of an instability however is not only the mode drive but also the strength of the competing damping. Energetic beam experiments at the AUG tokamak revealed the importance of the heavily temperature dependent Landau damping and its role in allowing strongly unstable modes to exist in a cold plasma [43][44][45].
In this paper we consider the active phase of ITER, where suprathermal alpha particles born through the fusion process in the burning plasma exist at significantly higher energies than present day experiments. This work aims at a scenario, where the post-thermal-quench healing of magnetic surfaces is fast enough to keep both runaway electrons and alpha particles confined. Good alpha particle confinement in the thermal quench stage is a necessary condition for the scenario described in the paper. This may not be universally true in all disruptions. However, if the breakup of magnetic surfaces is sufficiently strong for a sufficiently long time to cause significant alpha particle losses, then the losses of seed runaway electrons -which possess larger thermal speeds -will be even larger. Investigating the possible suppression of runaways via alpha-driven modes is interesting in scenarios where the runaways are reasonably confined [46][47][48], and if runaways are confined then alphas are likely confined as well. If there are some runaways being generated, the runaway current can provide the subsequent magnetic equilibrium to confine the alphas post-quench.
We show that the alpha particle distribution remains sufficiently energetic during the disruptions considered to drive TAEs in the current quench, where the plasma temperature (hence Landau damping) has dropped significantly. We also show the presence of these modes to cause significant transport to a runaway seed population. The main reason for energetic alphas to exist in this stage of a discharge is that the alpha suprathermal collision time is long compared to the thermal quench time scale.
The paper is structured as follows. In section 2 we model the evolution of the main plasma parameter profiles during the disruption using GO [10,49,50], which are also used to construct the magnetic equilibrium. In section 3 we describe the calculation of the spatiotemporal evolution of the alpha distribution function using the CODION tool [38], and we show that a substantial suprathermal alpha population exists in the current quench. In section 4 we describe the LIGKA model [51] used to identify the Alfvén spectrum and mode structures. As evidenced by simulations carried out by the relativistic version of HAGIS [52] introduced in section 5, the presence of the alpha population drives these modes to amplitudes of up to δB/B ≈ 0.1%. Finally, in section 6 we use HAGIS to determine the transport of seed runaways in the presence of these modes, and show seed transport.

ITER natural disruption scenario
In this paper we consider unmitigated ITER disruptions, as our aim is to investigate the possibility of inherent, "natural" runaway suppression mechanisms. Furthermore, the lack of mitigation significantly reduces the dimensionality of the parameter space to be considered when selecting a disruption scenario. Mitigated ITER disruptions [10] in a similar context are left for a future study.
The pre-disruption scenario is that of the 15 MA inductive burning plasma "scenario #2" described by Polevoi et al. [53,54], which has been extensively studied in the literature [41,[55][56][57][58][59]. High-current scenarios are also expected to produce the largest and most energetic populations of runaway electrons [60]. The plasma background consists of a 1:1 mixture of deuterium and tritium. Main parameters of the predisruption scenario are shown in table 1 and the temperature profiles in figure 2a.
In order to study the post-disruption evolution of runaways and Alfvén waves, we need a disruption scenario and a magnetic equilibrium. As first step in constructing these, we use the GO-code [10,49,50] to solve the 1D induction equation and to obtain the time evolution of the induced electric field E, the current density j for both the Ohmic current and the different runaway generation mechanisms. In this paper the hot-tail [61], the Dreicer [62] and the avalanche [6] mechanisms are considered. Since we are investigating a mechanism that would act on runaway seed particles, we are focusing on a scenario that is dominated by primary generation through the hot-tail mechanism. For this reason runaway seed sources coming from the active phase of ITER (i.e. tritium beta decay and inverse Compton scattering [10,63]) are not considered, as they would contribute a comparatively small seed in the selected scenario. While GO is capable of self-consistent simulation of the thermal collapse in disruptions, in order to reduce the dimensionality of the parameter space and the necessary run times for a parameter scan we are externally forcing a thermal quench: where T f is the final temperature, t N = t/t 0 where t = 0 marks the instant of the thermal quench starting, t 0 is the exponential decay time and T (r, t) is the time evolving temperature profile, equal for both electrons and ions. Since t N is also a temperature axis for the background species, it becomes convenient to normalize the time this way and will be used throughout the sections to come. Within the GO framework we neglect impurities and alphas for simplicity (which would only have a secondary effect on the evolution through conductivity and by slightly modifying the effective charge in the Dreicer and avalanche sources). Since the thermal collapse is externally forced, we also set all plasma species at the same temperature, i.e. assuming fast equilibration between the electron and ion thermal distributions. Figure 1 depicts the GO-simulation output of a disruption identified by T f = 3 eV and t 0 = 0.7 ms. A significant fraction of the current is converted to runaway through the hot-tail mechanism and due to this, little flux is available for the Dreicer mechanism. With the induced electric field setting in, the avalanche effect multiplies this initial seed. Note that with eq. (1) and in a scanned range of t 0 = 0.1 − 1 ms the hot-tail current always generates around t N ≈ 5 and the electric field is induced around t N ≈ 8. The induced electric field will be discussed in more detail in section 3 in the context of its influence on the alpha particle distribution.

Spatiotemporal evolution of the alpha particle distribution
In this section we are going to discuss what happens with the fusion-born alpha particle population in a disruption where a significant runaway current is still providing confinement. This question may also be important for the accurate calculation of wall heat loads during the disruption, but in this paper we are focusing on the drive to Alfvénic instabilities through fast ion resonances. Since mode drive can manifest both through momentum-space and real-space anisotropies, it is necessary to calculate the 1D+2V alpha particle distribution function during the disruption.
The gyro-averaged kinetic equation for the alpha particles in a homogeneous plasma with a Fokker-Planck collision operator with the source term S α for the fusion of alpha particles is where f α is the alpha particle distribution function, e the elementary charge, Z α = 2 is the alpha particle charge number, m α the alpha particle mass and C α,s the linearized collision operator describing collisions with species s. The particle pitch ξ = v /v is defined with respect to the equilibrium magnetic field lines. Equation 2 is a reduction of the ion kinetic equation, computing the time evolution of a distribution function in velocity space (v, ξ) under the influence of electric field acceleration, the Lorentz force and an accumulation of small-angle Coulomb collisions.
Trading off spatial features such as neoclassical particle transport is justified by the charged particle motion in a tokamak being dominated by parallel dynamics because of the preferred direction of the conductivity (σ σ ⊥ ). Perpendicular dynamics will be discussed later on. The magnetic field strength does not enter as a quantity into the equation, as at the energies investigated Bremsstrahlung and synchrotron emission losses can be neglected for the alpha particles.
As discussed earlier, we consider a 1:1 mixture of deuterium and tritium. The fusion process births alpha particles isotropically, at an energy of E α = 3.5 MeV. The empirically derived reaction rate for Deuterium and Tritium is (NRL [64]) where T i is the temperature of both deuterium and tritium given in keV and must be below 25 keV for the formula to be valid. The reaction rate for deuterium-deuterium fusion is not included, as it is generally two orders of magnitude lower. With equal D-T ion densities n D,T and temperatures T i , transforming from the fusion reaction centerof-mass frame to the lab frame results in the following energy dependence of the alpha particle source [65]: We solve the kinetic equation (2) numerically with the tool CODION (COllisional Distribution of IONs) [38]. CODION was originally developed to study highly energetic ion runaway mechanisms [66] in fusion and astrophysical plasmas. It can be considered an extension of previously existing analytical models [39] for the cases of low electric field and trace impurities. CODION uses a linearized collision operator, which is valid as long as the density of alpha particles is sufficiently small compared to the background.
In order to compute the radial distribution, independent CODION calculations are executed on a radial grid spanning the plasma radius with 101 points. The initial steadystate alpha distribution is established self-consistently by providing the initial profiles of temperature T s (r) (see figure 2a) and densities n s (r). We take the electron density as radially flat and constant at a value of n e = 10 20 m −3 and the 50:50 fuel ion densities n D,T to fulfill quasi-neutrality. The alpha density n α and pressure p α , which is necessary for later calculations, are evaluated (assuming isotropy and applying gyro-averaging) with the use of moments Isotropy is assumed as the birth process is isotropic, and the pre-disruptive electric field is small, therefore for the initialization simulation set to zero. The integrals are computed with the use of Simpson quadrature weights. For later use, we further specify . The alpha particles are born at v α ≈ 13·10 6 m/s with a spread ∆v corresponding to the pre-disruptive electron background temperature T pre 0,e = 25 keV. The temperature dependent velocities v EP above which particles are considered energetic are marked as a triangle of the corresponding color. Also included are the Alfvén velocities v A for the plasma composition chosen. Note that the total density integral is conserved throughout the disruption as the fusion process comes to a hold and particle losses are ignored. the energetic part of the alpha population by limiting the lower bound of the integrals to v EP = 10 2T (r, t)/m α and thus define energetic particle density n α,EP and energetic particle pressure p α,EP . This definition is not suited for pre-disruption analysis (as the cut-off velocity would be too high), however, it is purposeful for our post-disruptive analysis. Starting with initially negligible alpha particle densities n α (for numerical reasons) we simulate the fusion process S α (eq. (3)) self-consistently until a desired n α (r = 0) = 0.01n e is reached. The alpha profile established by our simulation is shown in figure 2a. We do not remove a fused D/T atom from the ion density profile nor do we add the born alpha particle to it, as it remains a minority species. Targeted for ITER is an optimal 50:50 D:T mixture and the efficiency of the pump-out of He-ash is not yet clear. Radial particle transport is not captured by the code. However, since we consider situations of good confinement in which the alpha particle loss time is approximately three orders of magnitude higher than its slowing down time, this seems to be a good approximation.
In velocity space, after sufficient simulation time, the alpha particles form a slowingdown distribution [67] where C is a constant proportional to S α , Erfc(x) is the complementary error function, v α = 2E α /m α is the birth velocity, v c is the crossover velocity and ∆v is the velocity spread of the fusion reactants at birth. In the absence of a particle sink, the alphas born at v α eventually thermalize into a Maxwellian f M of background temperature (helium ash). The total alpha distribution therefore consists of a slowing-down part and a thermal Maxwellian f α = f M + f SD . Using this, we fit the CODION simulated timedependent distributions and determine v c , ∆v and C for further analytical processing of the energetic tail (v > v EP ), which we assume to be the described fully by the slowingdown part f SD . In steady state these coefficients can be determined through the given plasma parameters.
In simulating the effect of the disruption, the initial steady-state alpha particle distribution is subjected within CODION to the time varying background profile of temperature as described in section 2. Effects of the induced electric field will be neglected since we only calculate f α up to t N = 6. Including the electric field in this duration is more expensive numerically, and only leads to a sub-1% change in the fast particle pressure. The high-energy alphas exist during the current quench due to their relatively slow deceleration, rather than due to electric field acceleration. For the disruption the fusion source is disabled, conserving the total density and density profiles. Background populations are assumed to maintain thermal equilibrium, while the collisional cooling of the alpha particles is consistently calculated by CODION. The disrupting alpha distribution function for the chosen example case T f = 3 eV and t 0 = 0.7 ms is shown in figure 2b. The energetic tail of alphas is largely conserved during simulation time, due to the short background cooling time compared to the collision time of the fast alphas. The analytically derived [68] slowing-down time for an alpha particle colliding with the steady state background on ITER is on the order of a second [69] and drops to the order of milliseconds for T ≈ 1 keV. Additionally marked in figure 2b are the velocities v EP above which the particles are considered energetic and also the Alfvén velocity v A = B/ √ µ 0 ρ, with the mass density ρ, the magnetic field strength B and µ 0 the permeability. Note how v EP (t N = 3) separates the Maxwellian from the energetic tail, which we assume to be fully described by f SD . We conduct a scan over various disruption scenarios defined by T f = [1 − 15] eV and t 0 = [0.1 − 1] ms. The presence of an energetic population we quantify as the ratio of energetic to total pressure Π ≡ p α,EP /p α and indicate the ability to drive modes via its gradient in real space. Figure 3a shows the time-slices of (normalized) core pressures and uses colorcode for the energetic particle pressure gradient p α,EP,grad = p α,EP (r/a = 0) − p α,EP (r/a = 0.6). The general behaviour of Π core (t) and its gradient in this plot can be imagined as a wave propagating towards higher t 0 with the effect of lower T f forcing earlier drops in pressure, as indicated by the time-slice for t = 0.73 ms. The peak of this "wave" gradually drops after reaching up to 80% for the quickest thermal quench simulated (note that the denominator in Π does not contain electron and background ion pressure). Figure 3b depicts the t 0 parameter space for T f = 3 eV and the energetic pressure in the core in absolute numbers. The general rise in p α,EP is due to the definition of v EP ∝ √ T and what is considered by us to be "energetic" in reference to the background (compare to v EP in figure 2b). Its magnitude being a product of particle densities and energies serves as an indication to the remnant of the tail throughout the disruption. The interesting result is found in difference in rise (3 < t N )  1), the x-axis can be interpreted as a background temperature. and the lack of difference in rise (0 < t N < 3). First, pressure evolution shares nearly identical behavior for the scanned range and under the assumptions made (exponential temperature decay, pure D-T composition, instantly thermalizing background species). In the context of eq. (1) t N becomes a temperature axis of the background species, thus a collision time-scale. Exemplary shown in figure 2b (for t 0 = 0.7 ms) was the barely changing alpha distribution until t N = 3. The qualitatively different evolution afterwards suggests a deviation caused by a growing discrepancy of elapsed time t = t N t 0 to collisionality. As expected, the quicker and more violent thermal quenches sustain a more significant energetic alpha particle tail in reference to the background temperature. This similarity for t < t N = 3 for the core pressure hold qualitatively true up to a radial point of r/a ≈ 0.5 in our simulations. In the outer half of the plasma the background temperature to begin with is low enough to collisionally drag the energetic tail early on.
In the absence of a bump-on-tail drive (due to the electric field not yet setting in) we evaluate the possibility of drive through the pressure gradient of the alphas. For a preferably general result we opt to evaluate the distributions and the pressure they exert at t N = 3. Our definition of v EP , which separates f M and f SD , makes sure that the energetic tail can optimally be described by the slowing down distribution equation for the time-point chosen. As is shown in figure figure 2b important TAE resonance regions for energies above 100 keV [70] are populated at this point in time. Figure 4a shows the individual species' pressures and temperature profile at t N = 3 (for the selected case of t 0 = 0.7 ms and T f = 3 eV) and is going to be used to reconstruct the equilibrium. The total pressure is determined as p tot = p e + p D,T + p α , where electron pressure p e and deuterium/tritium pressure p D,T are calculated from the ideal gas equation. Though the alpha particles are a minority species in the plasma, they contribute significantly to the pressure and its gradient due to their large kinetic energy. In initial steady state, the alpha pressure contributed to ≈ 10% of the total pressure on axis. With the thermal background cooling rapidly its relative relevance grows. It even briefly dominates in the core before vanishing together with the energetic pressure at t > 6t N (compare to figure 3b). Note that the Alfvén mode drive we are investigating in the following sections is determined by the actual distribution functions f SD (r) and not simply by their integral moments.

The Alfvén spectrum in the current quench
In order to determine the Alfvén spectrum, it is necessary to reconstruct the equilibrium. The characteristic time point chosen during the disruption is set at t N = 3. On the one hand there is still a substantial population of alpha particles, on the other hand the core background temperature has fallen to ≈ 1232 eV, which means that various mode damping mechanisms are also lowered. The plasma equilibrium is constructed using the pressure profiles (figure 4a) and the current density profiles from GO (figure 4b) using the VMEC code [71]. As aforementioned, the pressure profiles are very similar for the various disruption scenarios up to this time-point and since no runaway current is generated yet (see section 2) the same holds true for the current density profiles. The total time-window for the mode-evolution is going to be 3t N < t < 6t N , during which we will assume the plasma equilibrium to be constant. The low post-disruption pressure has diminishing influence on shaping the equilibrium and the current density profile is essentially constant in this time-frame, since the current quench is just about to begin (see figure 1a).
In order to evaluate the Alfvén spectrum, as well as the frequency and mode structure of modes possibly supported by this equilibrium, we employ the LIGKA code (LInear GyroKinetic Alfvén physics) [51]. We carry out a scan for the absolute scaling of the safety factor to account for deviations in the scenario and the time evolution of the plasma current, while maintaining the shape (determined by the plasma background). Figure 5a shows the results of this scan and reveals frequency gaps for TAEs of toroidal mode numbers 1 < n < 30. Due to the alpha particle orbit width in ITER the toroidal mode numbers are high compared to present day tokamaks. The solutions shown are from the even TAE branch that have been found to be the most unstable AEs in previous ITER analyses [41]. Also, Beta-induced Alfvén Eigenmodes (BAEs) are found to be present in the steep pressure region, however not further considered in this work. TAEs and BAEs are often deemed dominant considering the particle transport they cause due to their generally low frequency hence higher (potential) particle displacement per wave-particle energy transfer [72]. We chose the nominal value of the core safetyfactor to be q 0 = 1.
n e ) and varies mainly due to the q-profile, since the electron profile is flat. The corresponding (normalized) radial structures are shown in figure 5b-c. Our modeling shows that most of the modes possibily sustained by the post-disruption equilibrium are core localized, which is beneficial for interaction with the mainly core localized runaway electron seed. The Alfén velocity v A has a value of v A 7.3 · 10 6 m/s for the plasma composition chosen. For TAEs the most fundamental resonances occur at v = v A and v = v A /3 [33], which is still populated by the energetic alpha tail in velocity space at t N < 6 (figure 2b).

Interaction of alpha particles and Alfvén waves
The time evolution of the modes and their saturation amplitude is a critical question to determine their potency for runaway transport. Earlier studies showed that a magnetic perturbation with an amplitude of about δB/B ≈ 0.1% is sufficient to suppress runaway avalanche [73,74], while more recent research [75] decreases this threshold by about a factor of 2.
The interaction of the modes with the alpha particles (providing the drive) and the runaways is calculated using HAGIS (HAmiltonian GuIding center System) [42,76,77]. HAGIS is a perturbative, non-linear wave-particle interaction model which allows the modes to evolve in the presence of EPs, and the EPs to redistribute in phase space due to the interaction with the modes self-consistently.
The equations of motions in HAGIS are written in Boozer coordinates, thus we assume for the radial coordinates to be related to the normalized poloidal flux ψ 1/2 ≡ s ≈ r/a when transferring the distributions between codes. Calculations are supplied with the equilibrium, fast alpha population and the mode structures introduced in the previous sections. HAGIS uses a δf formalism, which allows us to omit the Maxwellian bulk and use the energetic part of the distribution f SD for analysis. The numerically calculated distributions for t 0 = 0.7 ms, T f = 3 eV, t N = 3 are fitted with the formula shown in eq. (6) to facilitate implementing the CODION distributions into HAGIS. The alpha particles are represented by N = 10 6 markers in phase space, which are initialized isotropically in pitch with velocities and positions in the tokamak according to f SD (r, t N = 3). Movement along the plasma equilibrium is dictated by their resonant interaction with the calculated modes. The interaction redistributes the particles and transfers energy to the modes, evolving their amplitude. The individual mode growth rates are competing against various damping mechanisms. LIGKA readily provides the damping rates caused by the ion/electron Landau mechanism and radiative damping and are of the order of γ/ω TAE ≈ 0.1%. This damping is typically 10 times smaller than reported for TAEs in the pre-disruption phase [41]. The reason is the Landau damping vanishing exponentially as the ion temperature drops to ≈ 1 keV [44].
As the temperature drops further, collisional damping by trapped electrons becomes increasingly important. We use eq. (2) from Gorelenkov et al. [78] to estimate the collisional damping rate for our disrupting plasma. Damping will be given in  [79] we also calculate fluid damping. Up to a resistivity of 0.56 · 10 −4 Ωm which corresponds to 6 eV background temperature, the damping is below 1%, hence will also not be included in our mode evolution simulations. Because the HAGIS model does not include collisional cooling of the driving alpha particles, its driving force only changes due to the redistribution of particles as dictated by their interaction with the modes. The loss of drive due to the thermal quenching however is captured by the CODION simulations. As indicated by figure 2b) and figure 3b), the loss of resonance with the Alfén velocity (first harmonic) occurs around t N = 6. As this is also the time-point up to which we calculated the collisional damping to be negligible, we restrict the window of mode evolution to 3t N < t < 6t N . Considering the fact that the dominating (ion Landau) damping mechanism drops exponentially with the temperature, it can be assumed that the mode excitation actually begins earlier than in our computations. Every initial mode amplitude in our simulation is set to δB/B = 10 −10 , though previous studies [41,42] have shown TAEs to be only marginally stable in ITER steady state with amplitudes of the order of δB/B = 10 −5 to 10 −4 . On the other hand, the HAGIS code is known to overestimate the saturation Figure 6. Evolution of mode amplitudes δB/B as caused by resonant interaction with the energetic alpha particle population f SD , 3t N into an ITER disruption in a) single mode and b) multi mode simulation. All available poloidal modes are included. Horizontal lines in b) mark t N = 6 for t 0 = [0.5, 0.7, 0.9] from which we extract the individual mode amplitudes for further use.
amplitude due to lack of zonal-flow physics [70], these however are typically only relevant at perturbation magnitudes well above the ones discussed in this work. Figure 6a shows single mode and figure 6b multi mode results. Both are qualitatively different because in multi-mode the energy transfer to waves and the subsequent redistribution of particles may push the modes through multiple resonance regions, as seen by the non-monotonic behavior. In single mode this redistribution may lead to a loss in drive and the damping taking over. The multi-mode simulation is the more realistic one, its linear growth rate of the most pronounced modes n = 23, 24, 25, 26 (in this time window) saturates at δB/B ≈ 0.05% after roughly 2 ms, a timescale which is sufficiently short even during a current quench. The effective growth rate in the linear phase for those modes is ≈ 14%. A saturation of the n = 8, 22 modes is observed approximately 2.5 ms into the simulation. The slowing-down distributions are isotropic in pitch and monotonically decreasing in velocity space, hence the mode drive comes mainly from gradients in real-space and is saturated by its flattening. In Appendix A we conduct the multi mode evolution with alpha particle densities increased to 130% and 200%, resulting in a growth rate in the linear phase of respectively 33% and 100%, indicating a high sensitivity. The saturated amplitude however is not strongly affected.
We want to look at the situation at t N = 6 into the disruption, since this is right before the avalanching begins (figure 1a) and the energetic particle pressure decays due to collisional cooling (figure 3b). Depending on the disruption scenario chosen, the amplitudes reached vary, as depicted in figure 6b. The maximum amplitude as well as the root mean square of all mode amplitudes is made note of in table 2. In the next section we simulate the interaction of said modes at their respective amplitudes at 6t N for t 0 = [0.5, 0.7, 0.9] with a seed of runaway electrons.  Table 2. HAGIS simulation results: Respectively maximum and average ( · ) mode amplitudes δB/B at t N = 6 (figure 6b) and maximum individual (max) and maximum ensemble-averaged ( · ) RE toroidal momentum changes ∆P ϕ /P ϕ (figure 7) caused by said modes after additional 2t N to an initial seed of runaway electrons. Additional simulations are run with mode amplitudes of the t 0 = 0.9 ms case multiplied by 1.3 and 2 (figure A2).

Transport of runaway electrons
In order to model the interaction of relativistic electrons with Alfvén waves, the HAGIS model had to be extended [22]. The derivation of the relativistic equations of motion in Boozer coordinates is provided in Appendix B. As the runaways are not expected to have a back-reaction to the wave evolution (due to the lack of suitable resonances) the calculations are run in a "passive" mode, where only the effect of the presence of the Alfvén modes is evaluated on the runaway electron test particles. The mode evolution is disabled and the amplitudes are set to their respective values at ultimately t N = 6 into the disruption ( figure 6b, table 2). In order to extract the interaction Green functions, 10000 test particles are launched on a phase space grid in energy, pitch and radial position. The radial region of interest is restricted to r/a = [0.05 − 0.45] and the energies are selected on a logarithmic grid ranging from 10 keV to 30 MeV. For each value of energy, pitch and radial position a total of 25 electrons are distributed evenly on the flux surface with uniform random distribution. With 5 radial positions, each phase space point is therefore represented by 125 electrons for statistical averaging.
We use the canonical toroidal momentum P ϕ (equation (B.2)) to quantify the displacement of the test particle orbits. The canonical toroidal momentum is a function of both parallel kinetic momentum and poloidal flux surface, i.e. a change in P ϕ in a sense represents the contribution to the change in the current profile of the given subrelativistic test particle. Changes to P ϕ are dominated by radial transport caused by magnetic perturbation of the modes when no resonant processes are happening. Those are unlikely given that the electron gyro-frequency on axis is 150 GHz and the parallel velocity of the lowest energy electron approximates to v || 6 · 10 7 m/s 8v A (for purely parallel motion).
In figure 7a we show the results of the passive simulation for the amplitudes extracted for the chosen t 0 = 0.7 ms case after 2t N integration time corresponding to 125 poloidal turns of each runaway particle. Each circle represents a radial starting position and its color shows the overall change in P ϕ , ensemble-averaged over the 25 REs that started on this flux surface. We observe minor changes of few % to P ϕ in few selected phase-space positions. However, when we take the mode amplitudes reached further into the thermal quench (figure 6b, 2.7 ms, corresponding to ultimately 6t N for t 0 = 0.9 ms), this picture changes significantly (figure 7b). As expected from the mode structures, the inner-and outermost particles initiated are not affected by the TAEs, however, the change is up to 10% for REs localized in between. The effect is stronger for higher particle pitches but applies for a wide range in the phase space.
Average/maximum values for the change in P ϕ are collected in table 2 for those simulations, including a reference run for the t 0 = 0.5 ms case. In addition, as a sensitivity scan we carried out simulations with 130% and 200% individual mode amplitude of the t 0 = 0.9 ms case ( figure A2, table 2). This fulfils the purpose of both a numerical sensitivity scan as well as it considers the experimental possibility of external amplitude amplification. At twice the amplitudes (max(δB/B) ≈ 0.17%) the average displacement on a flux surface is found to be as high as 25% for various points in initial phase space, including the most dangerous MeV REs.
These simulations do not evaluate the runaway electron dynamics, but serve as an indication to the possibility to transport runaways via alpha particle driven modes. Recent studies [75] suggest that the perturbation amplitudes and particle displacement caused by the modes discussed in this paper can lead to runaway avalanche mitigation (or even suppression). If the core transport is enhanced by for example alpha-driven modes, runaways are more easily transported to the edge, where other methods, such as Resonant Magnetic Perturbations [59,74] could further aid the removal of runaways.

Summary and outlook
In this paper we simulate unmitigated disruptions in burning ITER plasmas and we show that alpha particles remain sufficiently energetic during the current quench to drive Alfvén modes sustained by the current quench equilibrium. In turn these modes can enhance the transport of the runaway electron seed prior to the onset of significant induced electric field, preceding the runaway avalanche.
The disruptions considered are characterized by an exponential temperature decay time t 0 = [0 − 1] ms, and we have conducted a parameter scan for the final temperature and speed of the thermal quench using the GO code. The output of these simulations was used in the Fokker-Planck solver CODION to track the collisional cooling of the alpha particle distribution. We find that the energetic tail is sustained into the current quench. Throughout the parameter range, the evolution of the simulated core fast particle pressure shows a general similarity up to t N ≡ t/t 0 = 3. The simulated current and pressure profiles are used as input to construct equilibria using the VMEC code. The Alfvén continuum of these equilibria is analysed using LIGKA, showing the possible existence of core-localized TAEs.
Using HAIGS we have calculated the growth and saturation amplitudes of these TAEs, driven by the energetic alpha particles. The modes reach a saturation level of up to δB/B ≈ 0.1% for cases of t 0 > 0.7 ms. Runaway electron test particles were used in the HAGIS simulations to analyse the impact of these TAEs on runaway transport. The onset and saturation of the modes occurs before the rise of the electric field induced in the current quench, before the start of significant runaway avalanche. We show that the runaway electrons can be subjected to significant displacement due to the presence of the TAEs, which may contribute to a reduction of the runaway electron avalanche. The parameter and sensitivity scans indicate that slower thermal quenches are beneficial from the perspective of the studied phenomenon, and that there is a high sensitivity of runaway transport to mode amplitude.
All of our simulations have considered good confinement of alpha particles after the thermal quench. We argue that this is a conservative estimate, i.e. when the breakup of magnetic surfaces is bad enough to lead to alpha particle losses, the runaway electrons -which have higher thermal speeds -are likely to get lost as well. A rising runaway current in the current quench may in fact contribute to alpha particle confinement and in turn the drive of TAEs. This paper analyses the "worst case scenario" of unmitigated disruptions. Future studies will have to be conducted to evaluate the possibility of alpha particle drive of Alfvénic instabilities in ITER disruptions mitigated by e.g. shattered pellet injection. or other means.
In conclusion, natural disruptions in D-T ITER plasmas may be able to provide a natural mechanism that contributes to runaway avalanche suppression. Following this proof-of-principle paper, further studies are necessary to identify the optimal disruption scenario which maximizes runaway transport, and the self-consistent evaluation of runaway dynamics will also be necessary. Figure A1 presents a sensitivity scan of the multi mode amplitude growth caused by the energetic alpha particle population as computed by HAGIS. While the original alpha particle density causes a growth rate of (at max) 14% in the linear phase, this growth rate increases to 33% (a) and 100% (b) for respectively 130% and 200% original density values. Though saturated earlier, the amplitudes reached after the linear phase are approximately the same. HAGIS is valid up to δB/B ≈ 1%. Figure A1. Multi mode amplitude evolution of available TAEs (including poloidal harmonics) during an ITER disruption as caused by an energetic alpha particle population. Its originally calculated density 3t N into the disruption is increased to a) 130% and b) 200% to establish a sensitivity estimation of figure 6b.
In figure A2 a scan is conducted where we take the respective t 0 = 0.9 ms mode amplitudes at a) 130% and b) 200% (see table 2). This serves both as a sensitivity scan and a means to explore the influence of additional external amplitude amplification. We observe a significant increase of runaway displacement as the TAE amplitude is raised, which suggests that a disruption scenario can be optimized to maximize runaway seed transport. p and the magnetic moment µ are defined as follows: where b is the unit vector along the magnetic field and v is the particle velocity vector. In Boozer coordinates, the magnetic field can be represented as B = I(ψ p )∇θ + g(ψ p )∇ϕ + δ(ψ p , θ)∇ψ p , and the vector potential as A = ψ t ∇θ − ψ p ∇ϕ. ψ t and ψ p are the toroidal and poloidal fluxes, I and g are the toroidal and poloidal currents, δ the radial covariant component of B, and θ and ϕ the poloidal and toroidal angles, respectively. Using Boozer coordinates for the guiding-centre Lagrangian, we find In order to achieve the natural Lagrarian form and read the canonical momenta straight away, we have to reabsorbψ p . This is done using the same argumentation as on p.47/48 of S. Pinches PhD thesis [84] for the derivation for the non-relativistic equations: we transform the guiding center velocityẋ →ẋ + w with We differentiate with respect to θ, ϕ, P θ and P ϕ and using q = ∂ψ t /∂ψ p to obtain