Alpha particle driven Alfvénic instabilities in ITER post-disruption plasmas

The subject of runaway electron (RE) mitigation is of crucial importance to the success of reactor-relevant tokamaks such as ITER [1–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–10]. A multi-megaampere runaway beam has the potential to inflict significant damage to plasma-facing components [11, 12]. This paper discusses a phenomenon which could potentially act to passively mitigate RE generation at ITER.

The interaction of runaways with plasma waves has been investigated both in theory [13–20] and through observations in multiple tokamaks [21–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 REs—due 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–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 REs as a possible drive for CAEs (or possibly GAEs) on the DIII-D tokamak in the context of a massive gas injection 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–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 REs 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 REs—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–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.

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, 42, 55–58]. High-current scenarios are also expected to produce the largest and most energetic populations of REs [59]. The plasma background consists of a 1:1 mixture of deuterium and tritium. Main parameters of the pre-disruption scenario are shown in table 1 and the temperature profiles in figure 2(a).

Table 1. The plasma parameters of the pre-disruption ITER scenario.

Parameter name	Notation	Value
Major radius	R0	6.195 m
Minor radius	a	2 m
Magnetic field on axis	B	5.26 T
Effective charge	Zeff	1.0
Normalized flux	ψ	Ψ(r)/Ψ(a)
Normalized radius	r/a	$r/a\simeq \sqrt{\psi }\equiv s$
Plasma current	Ip	15 MA
Electron density	ne(ψ)	1020 m−3
Ion density	ni(ψ)	1020 m−3
Electron temperature on axis	${T}_{0,\mathrm{e}}^{\text{pre}}$	24.7 keV
Ion temperature on axis	${T}_{0,\mathrm{i}}^{\text{pre}}$	21.2 keV
q on axis	q0	1.07
q on edge	qa	3.67

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 [60], the Dreicer [61] 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 [9, 10]) 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:

$$\begin{equation}T\left(r,t\right)={T}_{\mathrm{f}}+\left[T\left(r,0\right)-{T}_{\mathrm{f}}\right]\mathrm{exp}\left\{-{t}_{\mathrm{N}}\right\},\end{equation} \tag{ 1 }$$

where T_f is the final temperature, t_N = t/t₀ where t = 0 marks the instant of the thermal quench starting, t₀ 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.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 equation (1) and in a scanned range of t₀ = 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.

Zoom In Zoom Out Reset image size

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

$$\begin{equation}\frac{\partial {f}_{\alpha }}{\partial t}+\frac{e{Z}_{\alpha }}{{m}_{\alpha }}E\left(\xi \frac{\partial }{\partial v}+\frac{1-{\xi }^{2}}{v}\frac{\partial }{\partial \xi }\right){f}_{\alpha }=\sum\limits _{s}{\mathcal{C}}_{\alpha ,s}\left\{{f}_{\alpha }\right\}+{S}_{\alpha },\end{equation} \tag{ 2 }$$

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 ${\mathcal{C}}_{\alpha ,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 [62])

$$\begin{equation*}{\langle \sigma v\rangle }_{\text{DT}}=3.68{\times}1{0}^{-12}{T}_{\mathrm{i}}^{-2/3}\enspace \mathrm{exp}\left\{-19.94{T}_{\mathrm{i}}^{-1/3}\right\}\enspace {\mathrm{c}\mathrm{m}}^{-3}\enspace {\mathrm{s}}^{-1},\end{equation*}$$

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 center-of-mass frame to the lab frame results in the following energy dependence of the alpha particle source [63]:

$$\begin{equation}{S}_{\alpha }={S}_{0}\enspace \mathrm{exp}\left\{-\frac{5}{16}\frac{{\left(\langle {E}_{\alpha }\rangle -{E}_{\alpha }\right)}^{2}}{{T}_{\mathrm{i}}\langle {E}_{\alpha }\rangle }\right\}.\end{equation} \tag{ 3 }$$

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 [64] 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 steady-state alpha distribution is established self-consistently by providing the initial profiles of temperature T_s(r) (see figure 2(a)) and densities n_s(r). We take the electron density as radially flat and constant at a value of n_e = 10²⁰ m⁻³ 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

$$\begin{equation}{p}_{\alpha }\left(r,t\right)=\frac{4\pi {m}_{\alpha }}{3}\int {v}^{4}{f}_{\alpha }\left(r,v,t\right)\mathrm{d}v,\end{equation} \tag{ 4 }$$

$$\begin{equation}{n}_{\alpha }\left(r,t\right)=4\pi \int {v}^{2}{f}_{\alpha }\left(r,v,t\right)\mathrm{d}v.\end{equation} \tag{ 5 }$$

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 energetic part of the alpha population by limiting the lower bound of the integrals to ${v}_{\text{EP}}=10\sqrt{2T\left(r,t\right)/{m}_{\alpha }}$ 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_α (equation (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 2(a). 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.

Zoom In Zoom Out Reset image size

In velocity space, after sufficient simulation time, the alpha particles form a slowing-down distribution [65]

$$\begin{equation}{f}_{\mathrm{SD}}\left(v\right)=\frac{C}{{v}_{\mathrm{c}}^{3}+{v}^{3}}\enspace \mathrm{e}\mathrm{r}\mathrm{f}\mathrm{c}\left(\frac{v-{v}_{\alpha }}{{\Delta}v}\right),\end{equation} \tag{ 6 }$$

where C is a constant proportional to S_α , erfc(x) is the complementary error function, ${v}_{\alpha }=\sqrt{2{E}_{\alpha }/{m}_{\alpha }}$ 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 time-dependent 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 slowing-down 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.7 ms is shown in figure 2(b). 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 [66] slowing-down time for an alpha particle colliding with the steady state background on ITER is on the order of a second [67] and drops to the order of milliseconds for T ≈ 1 keV. Additionally marked in figure 2(b) are the velocities v_EP above which the particles are considered energetic and also the Alfvén velocity ${v}_{\mathrm{A}}=B/\sqrt{{\mu }_{0}\rho }$, with the mass density ρ, the magnetic field strength B and μ₀ 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.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 3(a) 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 behavior of Π_core(t) and its gradient in this plot can be imagined as a wave propagating towards higher t₀ 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 3(b) depicts the t₀ 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}_{\text{EP}}\propto \sqrt{T}$ and what is considered by us to be 'energetic' in reference to the background (compare to v_EP in figure 2(b)). 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) 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 equation (1) t_N becomes a temperature axis of the background species, thus a collision time-scale. Exemplary shown in figure 2(b) (for t₀ = 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₀ 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.

Zoom In Zoom Out Reset image size

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 2(b) important TAE resonance regions for energies above 100 keV [68] are populated at this point in time.

Figure 4(a) shows the individual species' pressures and temperature profile at t_N = 3 (for the selected case of t₀ = 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 3(b)). 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.

Zoom In Zoom Out Reset image size

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 4(a)) and the current density profiles from GO (figure 4(b)) using the VMEC code [69]. 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 1(a)).

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 5(a) 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 [70]. We chose the nominal value of the core safety-factor to be q₀ = 1.071 with the toroidal harmonic ranges n = [7–15, 22–26] and respective poloidal harmonics $m=\left[\left\{\left(n-2\right)-\left(n+4\right)\right\},\enspace \left\{\left(n-2\right)-\left(n+6\right)\right\}\right]$ as justified by the core-localization. The center of the mode frequency gap is located at ${\omega }_{\text{TAE}}={v}_{\mathrm{A}}/\left(4\pi qR\right)\propto 1/\left(q\sqrt{{n}_{\mathrm{e}}}\right)$ and varies mainly due to the q-profile, since the electron profile is flat. The corresponding (normalized) radial structures are shown in figures 5(b) and (c). Our modeling shows that most of the modes possibly sustained by the post-disruption equilibrium are core localized, which is beneficial for interaction with the mainly core localized RE seed. The Alfén velocity v_A has a value of v_A ≃ 7.3 × 10⁶ 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 2(b)).

Zoom In Zoom Out Reset image size

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 [71, 72], while more recent research [73] 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, 74, 75]. 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.7 ms, T_f = 3 eV, t_N = 3 are fitted with the formula shown in equation (6) to facilitate implementing the CODION distributions into HAGIS. The alpha particles are represented by N = 10⁶ 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 equation (2) from Gorelenkov and Sharapov [76] to estimate the collisional damping rate for our disrupting plasma. Damping will be given in $\left[\gamma /\omega \right]=\left[{s}^{-1}/{s}^{-1}\right]$. With the electron–electron collision rate ν_e(T) and the plasma beta β_e(T) its value

$$\begin{equation*}\frac{{\gamma }_{\mathrm{e}}}{{\omega }_{\text{TAE}}}=-2.1{\beta }_{\mathrm{e}}{\left(\frac{{\nu }_{\mathrm{e}}}{{\omega }_{\text{TAE}}}\right)}^{1/2}\left[\mathrm{ln}\enspace 8{\left(\frac{2r{\omega }_{\text{TAE}}}{R{\nu }_{\mathrm{e}}}\right)}^{1/2}\right]\end{equation*}$$

reaches 1.2% at T = 60 eV for the n = 8 mode calculated and only varies slightly for the other modes calculated. Hence we neglect it for t < 6t_N ≃ 64 eV. Using the resistive MHD code CASTOR (Complex Alfvén Spectrum of TORoidal plasmas) [77] we also calculate fluid damping. Up to a resistivity of ≃0.56 × 10⁻⁴ Ω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 figures 2(b) and 3(b), 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⁻¹⁰, 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⁻⁵ to 10⁻⁴. On the other hand, the HAGIS code is known to overestimate the saturation amplitude due to lack of zonal-flow physics [68], these however are typically only relevant at perturbation magnitudes well above the ones discussed in this work.

Figure 6(a) shows single mode and figure 6(b) 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.

Zoom In Zoom Out Reset image size

We want to look at the situation at t_N = 6 into the disruption, since this is right before the avalanching begins (figure 1(a)) and the energetic particle pressure decays due to collisional cooling (figure 3(b)). Depending on the disruption scenario chosen, the amplitudes reached vary, as depicted in figure 6(b). 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.5, 0.7, 0.9] with a seed of REs.

Table 2. HAGIS simulation results: respectively maximum and average (⟨⋅⟩) mode amplitudes δB/B at t_N = 6 (figure 6(b)) 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 REs. Additional simulations are run with mode amplitudes of the t₀ = 0.9 ms case multiplied by 1.3 and 2 (figure A2).

t0	$\mathrm{max}\left(\delta B/B\right)$ (%)	⟨δB/B⟩ (%)	$\mathrm{max}\left({\Delta}{O}_{\varphi }/{P}_{\varphi }\right)$ (%)	⟨ΔPφ /Pφ ⟩ (%)
0.5	0.006	0.002	3	1
0.7	0.056	0.027	11	8
0.9	0.083	0.042	18	10
0.9	0.083 × 1.3	0.042 × 1.3	39	13
0.9	0.083 × 2.0	0.042 × 2.0	48	25

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 RE 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 6(b), table 2).

In order to extract the interaction Green functions, 10 000 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 sub-relativistic 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⁷ m s⁻¹ ≃ 8v_A (for purely parallel motion).

In figure 7(a) we show the results of the passive simulation for the amplitudes extracted for the chosen t₀ = 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 6(b) 2.7 ms, corresponding to ultimately 6t_N for t₀ = 0.9 ms), this picture changes significantly (figure 7(b)). 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.

Zoom In Zoom Out Reset image size

Average/maximum values for the change in P_φ are collected in table 2 for those simulations, including a reference run for the t₀ = 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.9 ms case (figure A2, table 2). This fulfills 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 RE dynamics, but serve as an indication to the possibility to transport runaways via alpha particle driven modes. Recent studies [73] 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 [58, 72] could further aid the removal of runaways.

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 RE 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–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₀ = 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 analyzed 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.7 ms. RE test particles were used in the HAGIS simulations to analyze 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 REs can be subjected to significant displacement due to the presence of the TAEs, which may contribute to a reduction of the RE 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 REs—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.

The authors are grateful to E. Strumberger, Th. W. Hayward-Schneider, M. Schneller and S. Newton for fruitful discussions. Ph. W. Lauber acknowledges support from the Enabling research projects NAT (CfP-AWP17-ENR-MPG-01) and MET (CfP-AWP19-ENR-01-ENEA-05). This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under Grant Agreement No. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

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%.

Zoom In Zoom Out Reset image size

In figure A2 a scan is conducted where we take the respective t₀ = 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.

Zoom In Zoom Out Reset image size

In this appendix, we derive the relativistic equations of motion for REs in HAGIS. The numerical implementation has been verified [22] through extensive comparisons with the ANTS code [58, 72, 78, 79] by tracking the same test particle ensembles in the same background equilibria. The following is a relativistic extension to the derivations in the PhD thesis of Pinches [80]. The constants used are e for the electric charge, c for the speed of light, m for the particle mass. The relativistic guiding-centre Lagrangian reads (in CGS units) [81]

$$\begin{equation*}L\left(\mathbf{x},\dot {\mathbf{x}},t\right)=\left[\frac{e}{c}\mathbf{A}\left(\mathbf{x},t\right)+\frac{{p}_{{\Vert}}}{B}\mathbf{B}\left(\mathbf{x},t\right)\right]\cdot \dot {\mathbf{x}}+\frac{mc}{e}\mu \dot {\vartheta }-{H}_{\mathrm{g}\mathrm{c}}\end{equation*}$$

with the relativistic guiding-centre Hamiltonian

$$\begin{equation*}{H}_{\mathrm{g}\mathrm{c}}=\gamma m{c}^{2}+e\phi \left(\mathbf{x},t\right)\end{equation*}$$

and the relativistic factor

$$\begin{equation*}\gamma =\sqrt{1+\frac{2}{m{c}^{2}}\mu B+\frac{{p}_{{\Vert}}^{2}}{{m}^{2}{c}^{2}}}.\end{equation*}$$

A and ϕ are the vector and electric potential, respectively, x is the guiding-centre position, ϑ is the gyroangle, B the magnetic field strength, and the parallel momentum p_|| and the magnetic moment μ are defined as follows:

$$\begin{align*}{p}_{{\Vert}}& =m\gamma \mathbf{v}\cdot \mathbf{b},\\ \mu & =m{\gamma }^{2}\frac{{\mathbf{v}}_{\perp }^{2}}{2B}\end{align*}$$

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

$$\begin{equation*}\mathbf{B}=I\left({\psi }_{\mathrm{p}}\right)\nabla \theta +g\left({\psi }_{\mathrm{p}}\right)\nabla \varphi +\delta \left({\psi }_{\mathrm{p}},\theta \right)\nabla {\psi }_{\mathrm{p}},\end{equation*}$$

and the vector potential as

$$\begin{equation*}\mathbf{A}={\psi }_{\mathrm{t}}\nabla \theta -{\psi }_{\mathrm{p}}\nabla \varphi .\end{equation*}$$

ψ_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

$$\begin{align*}L& =\left(\frac{e}{c}\mathbf{A}+\frac{{p}_{{\Vert}}}{B}\mathbf{B}\right)\cdot \dot {\mathbf{x}}+\frac{mc}{e}\mu \dot {\vartheta }-{H}_{\mathrm{g}\mathrm{c}}\\ & =\left(\frac{e}{c}\left({\psi }_{\mathrm{t}}\nabla \theta -{\psi }_{\mathrm{p}}\nabla \varphi \right)+\frac{{p}_{{\Vert}}}{B}\left(\delta \nabla {\psi }_{\mathrm{p}}+I\nabla \theta +g\nabla \varphi \right)\right)\cdot \dot {\mathbf{x}}\\ & \quad +\frac{mc}{e}\mu \dot {\vartheta }-{H}_{\mathrm{g}\mathrm{c}}\\ & =\left(\frac{e}{c}{\psi }_{\mathrm{t}}+\frac{{p}_{{\Vert}}I}{B}\right)\dot {\theta }+\left(\frac{{p}_{{\Vert}}}{B}g-\frac{e}{c}{\psi }_{\mathrm{p}}\right)\dot {\varphi }+\frac{{p}_{{\Vert}}}{B}\delta {\dot {\psi }}_{\mathrm{p}}\\ & \quad +\frac{mc}{e}\mu \dot {\vartheta }-{H}_{\mathrm{g}\mathrm{c}}.\end{align*}$$

In order to achieve the natural Lagrangian form and read the canonical momenta straight away, we have to reabsorb ${\dot {\psi }}_{\mathrm{p}}$. This is done using the same argumentation as on p 47/48 of Pinches PhD thesis [80] for the derivation for the non-relativistic equations: we transform the guiding center velocity $\dot {x}\to \dot {x}+w$ with

$$\begin{equation*}\mathbf{w}=-\frac{\left(\frac{{p}_{{\Vert}}}{B}\delta +\frac{e}{c}{\tilde {A}}_{{\psi }_{\mathrm{p}}}\right){\dot {\psi }}_{\mathrm{p}}}{\frac{e}{c}\mathbf{A}\cdot \mathbf{B}+\frac{{p}_{{\Vert}}}{B}{B}^{2}}\mathbf{B}.\end{equation*}$$

We now add perturbations to both potentials of the form

$$\begin{align*}& \mathbf{A}\left(\mathbf{x},t\right)={\tilde {A}}_{{\psi }_{\mathrm{p}}}\nabla {\psi }_{\mathrm{p}}+{\tilde {A}}_{\theta }\nabla \theta +{\tilde {A}}_{\varphi }\nabla \varphi \quad \text{and}\\ & \tilde {\phi }=\tilde {\phi }\left(\mathbf{x},t\right),\end{align*}$$

and are now able to express the canonical momenta for the set of canonical variables (θ, φ, ϑ) as

$$\begin{equation}{P}_{\theta }=\frac{{p}_{{\Vert}}}{B}I+\frac{e}{c}{\psi }_{\mathrm{t}}+\frac{e}{c}{\tilde {A}}_{\theta }\end{equation} \tag{ B.1 }$$

$$\begin{equation}{P}_{\varphi }=\frac{{p}_{{\Vert}}}{B}g-\frac{e}{c}{\psi }_{\mathrm{p}}+\frac{e}{c}{\tilde {A}}_{\varphi }\end{equation} \tag{ B.2 }$$

$$\begin{equation*}{P}_{\vartheta }=\frac{mc}{e}\mu .\end{equation*}$$

B.1. Equations of motion

Combining canonical momenta (B.1) and (B.2) yields

$$\begin{equation*}g\left({P}_{\theta }-\frac{e}{c}{\psi }_{\mathrm{t}}-\frac{e}{c}{\tilde {A}}_{\theta }\right)=I\left({P}_{\varphi }+\frac{e}{c}{\psi }_{\mathrm{p}}-\frac{e}{c}{\tilde {A}}_{\varphi }\right).\end{equation*}$$

We differentiate with respect to θ, φ, P_θ and P_φ and using q = ∂ψ_t/∂ψ_p to obtain

$$\begin{align*}\frac{\mathrm{d}}{\mathrm{d}\theta }& :{g}^{\prime }\frac{\partial {\psi }_{\mathrm{p}}}{\partial \theta }\frac{{p}_{{\Vert}}I}{B}+g\left(-q\frac{e}{c}-\frac{e}{c}\left(\frac{\partial {\tilde {A}}_{\theta }}{\partial \theta }+{\tilde {A}}_{\theta }^{\prime }\frac{\partial {\psi }_{\mathrm{p}}}{\partial \theta }\right)\right)\\ & \quad -{I}^{\prime }\frac{\partial {\psi }_{\mathrm{p}}}{\partial \theta }\frac{{p}_{{\Vert}}g}{B}-I\left(\frac{e}{c}-\frac{e}{c}\left(\frac{\partial {\tilde {A}}_{\varphi }}{\partial \theta }+{\tilde {A}}_{\varphi }^{\prime }\frac{\partial {\psi }_{\mathrm{p}}}{\partial \theta }\right)\right)=0\\ & {\Leftrightarrow}\frac{\partial {\psi }_{\mathrm{p}}}{\partial \theta }=\frac{1}{{D}_{r}}\left(I\frac{\partial {\tilde {A}}_{\varphi }}{\partial \theta }-g\frac{\partial {\tilde {A}}_{\theta }}{\partial \theta }\right)\end{align*},$$$

where

$$\begin{equation*}{D}_{r}=\frac{c}{e}\frac{{p}_{{\Vert}}}{B}\left({I}^{\prime }g-I{g}^{\prime }\right)+I+qg+g{\tilde {A}}_{\theta }^{\prime }-I{\tilde {A}}_{\varphi }^{\prime }.\end{equation*}$$

Analogously,

$$\begin{equation*}\frac{\partial {\psi }_{\mathrm{p}}}{\partial \varphi }=\frac{1}{{D}_{r}}\left(I\frac{\partial {\tilde {A}}_{\varphi }}{\partial \varphi }-g\frac{\partial {\tilde {A}}_{\theta }}{\partial \varphi }\right).\end{equation*}$$

Thus the next differentiation becomes

$$\begin{align*}\frac{\mathrm{d}}{\mathrm{d}{P}_{\theta }}& :{g}^{\prime }\frac{\partial {\psi }_{\mathrm{p}}}{\partial {P}_{\theta }}\frac{{p}_{{\Vert}}I}{B}+g\left(1-\frac{e}{c}q\frac{\partial {\psi }_{\mathrm{p}}}{\partial {P}_{\theta }}-\frac{e}{c}{\tilde {A}}_{\theta }^{\prime }\frac{\partial {\psi }_{\mathrm{p}}}{\partial {P}_{\theta }}\right.\\ & \quad \left.-{I}^{\prime }\frac{\partial {\psi }_{\mathrm{p}}}{\partial {P}_{\theta }}\frac{{p}_{{\Vert}}g}{B}-I\left(\frac{e}{c}-\frac{e}{c}{\tilde {A}}_{\varphi }^{\prime }\frac{\partial {\psi }_{\mathrm{p}}}{\partial {P}_{\theta }}\right)\right)=0\\ & {\Leftrightarrow}\frac{\partial {\psi }_{\mathrm{p}}}{\partial {P}_{\theta }}=\frac{c}{e}\frac{g}{{D}_{r}}\end{align*}$$

and

$$\begin{equation*}\frac{\partial {\psi }_{\mathrm{p}}}{\partial {P}_{\varphi }}=-\frac{c}{e}\frac{I}{{D}_{r}}\end{equation*}$$

when performing the same operation with d/dP_φ . Differentiation of p_||/B with respect to θ, φ, P_θ and P_φ yields

$$\begin{align*}\frac{\mathrm{d}}{\mathrm{d}\theta }& :0=\frac{\partial }{\partial \theta }\left(\frac{{p}_{{\Vert}}}{B}\right)I+\frac{{p}_{{\Vert}}}{B}{I}^{\prime }\frac{\partial {\psi }_{\mathrm{p}}}{\partial \theta }+\frac{e}{c}q\frac{\partial {\psi }_{\mathrm{p}}}{\partial \theta }\\ & \quad +\frac{e}{c}\left(\frac{\partial {\tilde {A}}_{\theta }}{\partial \theta }+{\tilde {A}}_{\theta }^{\prime }\frac{\partial {\psi }_{\mathrm{p}}}{\partial \theta }\right)\\ & {\Leftrightarrow}\frac{\partial }{\partial \theta }\left(\frac{{p}_{{\Vert}}}{B}\right)=-\frac{1}{I{D}_{r}}\left(\left(\frac{{p}_{{\Vert}}}{B}{I}^{\prime }+\frac{e}{c}q+\frac{e}{c}{\tilde {A}}_{\theta }^{\prime }\right)\right.\\ & \qquad {\times}\left.\left(I\frac{\partial {\tilde {A}}_{\varphi }}{\partial \theta }-g\frac{\partial {\tilde {A}}_{\theta }}{\partial \theta }\right)+\frac{e}{c}\frac{\partial {\tilde {A}}_{\theta }}{\partial \theta }{D}_{r}\right)\\ & \qquad =-\frac{1}{{D}_{r}}\left(\frac{\partial {\tilde {A}}_{\varphi }}{\partial \theta }\left(\frac{{p}_{{\Vert}}}{B}{I}^{\prime }+\frac{e}{c}q+\frac{e}{c}{\tilde {A}}_{\theta }^{\prime }\right)\right.\\ & \qquad \left.+\frac{\partial {\tilde {A}}_{\theta }}{\partial \theta }\left(-\frac{{p}_{{\Vert}}}{B}{g}^{\prime }+\frac{e}{c}-\frac{e}{c}{\tilde {A}}_{\varphi }^{\prime }\right)\right)\\ & {\Leftrightarrow}\frac{\partial }{\partial \theta }\left(\frac{{p}_{{\Vert}}}{B}\right)=\frac{1}{{D}_{r}}{\times}\left(\frac{\partial {\tilde {A}}_{\theta }}{\partial \theta }\left(\frac{{p}_{{\Vert}}}{B}{g}^{\prime }-\frac{e}{c}+\frac{e}{c}{\tilde {A}}_{\varphi }^{\prime }\right)\right.\\ & \qquad \left.-\frac{\partial {\tilde {A}}_{\varphi }}{\partial \theta }\left(\frac{{p}_{{\Vert}}}{B}{I}^{\prime }+\frac{e}{c}q+\frac{e}{c}{\tilde {A}}_{\theta }^{\prime }\right)\right).\end{align*}$$

Analogously

$$\begin{align*}\frac{\partial }{\partial \varphi }\left(\frac{{p}_{{\Vert}}}{B}\right)& =\frac{1}{{D}_{r}}\left(\frac{\partial {\tilde {A}}_{\theta }}{\partial \varphi }\left(\frac{{p}_{{\Vert}}}{B}{g}^{\prime }-\frac{e}{c}+\frac{e}{c}{\tilde {A}}_{\varphi }^{\prime }\right)\right.\\ & \quad \left.-\frac{\partial {\tilde {A}}_{\varphi }}{\partial \varphi }\left(\frac{{p}_{{\Vert}}}{B}{I}^{\prime }+\frac{e}{c}q+\frac{e}{c}{\tilde {A}}_{\theta }^{\prime }\right)\right).\end{align*}$$

$$\begin{align*}\frac{\mathrm{d}}{\mathrm{d}{P}_{\theta }}:1& =\frac{\partial }{\partial {P}_{\theta }}\left(\frac{{p}_{{\Vert}}}{B}\right)I+\frac{{p}_{{\Vert}}}{B}{I}^{\prime }\frac{\partial {\psi }_{\mathrm{p}}}{\partial {P}_{\theta }}+\frac{e}{c}q\frac{\partial {\psi }_{\mathrm{p}}}{\partial {P}_{\theta }}+\frac{e}{c}{\tilde {A}}_{\theta }^{\prime }\frac{\partial {\psi }_{\mathrm{p}}}{\partial {P}_{\theta }}\\ & {\Leftrightarrow}\frac{\partial }{\partial {P}_{\theta }}\left(\frac{{p}_{{\Vert}}}{B}\right)=\frac{1}{I{D}_{r}}\left({D}_{r}-\left(\frac{{p}_{{\Vert}}}{B}{I}^{\prime }+\frac{e}{c}q+\frac{e}{c}{\tilde {A}}_{\theta }^{\prime }\right)\frac{c}{e}g\right)\\ & =\frac{1}{{D}_{r}}\left(1-{\tilde {A}}_{\varphi }^{\prime }-\frac{c}{e}\frac{{p}_{{\Vert}}}{B}{g}^{\prime }\right)\end{align*}$$

and

$$\begin{equation*}\frac{\partial }{\partial {P}_{\varphi }}\left(\frac{{p}_{{\Vert}}}{B}\right)=\frac{1}{{D}_{r}}\left(q+{\tilde {A}}_{\theta }^{\prime }+\frac{c}{e}\frac{{p}_{{\Vert}}}{B}{I}^{\prime }\right).\end{equation*}$$

Combining above expressions for the equations of motion

$$\begin{align*}\frac{\partial H}{\partial {P}_{\theta }}& =\dot {\theta }=\frac{\partial H}{\partial {\psi }_{\mathrm{p}}}\frac{\partial {\psi }_{\mathrm{p}}}{\partial {P}_{\theta }}+\frac{\partial H}{\partial \left(\frac{{p}_{{\Vert}}}{B}\right)}\frac{\partial \left(\frac{{p}_{{\Vert}}}{B}\right)}{\partial {P}_{\theta }}\\ \frac{\partial H}{\partial {\psi }_{\mathrm{p}}}& =\frac{{B}^{\prime }}{\gamma }\left(\mu +\frac{B}{m}{\left(\frac{{p}_{{\Vert}}}{B}\right)}^{2}\right)+e{\phi }^{\prime }\\ \frac{\partial H}{\partial \left(\frac{{p}_{{\Vert}}}{B}\right)}& =\frac{{B}^{2}}{\gamma m}\frac{{p}_{{\Vert}}}{B}\end{align*}$$

yields

$$\begin{align*}{\Leftrightarrow}\dot {\theta }& =\frac{1}{{D}_{r}}\left[\left(\frac{{B}^{\prime }}{\gamma }\left(\mu +\frac{B}{m}{\left(\frac{{p}_{{\Vert}}}{B}\right)}^{2}\right)+e{\phi }^{\prime }\right)\frac{c}{e}g\right.\\ & \quad \left.+\frac{{B}^{2}}{\gamma m}\frac{{p}_{{\Vert}}}{B}\left(1-{\tilde {A}}_{\varphi }^{\prime }-\frac{c}{e}\frac{{p}_{{\Vert}}}{B}{g}^{\prime }\right)\right]\\ \dot {\varphi }& =\frac{1}{{D}_{r}}\left[-\left(\frac{{B}^{\prime }}{\gamma }\left(\mu +\frac{B}{m}{\left(\frac{{p}_{{\Vert}}}{B}\right)}^{2}\right)+e{\phi }^{\prime }\right)\frac{c}{e}I\right.\\ & \quad \left.+\frac{{B}^{2}}{\gamma m}\frac{{p}_{{\Vert}}}{B}\left(q+{\tilde {A}}_{\theta }^{\prime }+\frac{c}{e}\frac{{p}_{{\Vert}}}{B}{I}^{\prime }\right)\right]\\ {\dot {P}}_{\theta }& =-\frac{1}{\gamma }\left(\mu +\frac{B}{m}{\left(\frac{{p}_{{\Vert}}}{B}\right)}^{2}\right)\frac{\partial B}{\partial \theta }-e\frac{\partial \phi }{\partial \theta }\\ & \quad -\frac{1}{{D}_{r}}\left(\frac{{B}^{\prime }}{\gamma }\left(\mu +\frac{B}{m}{\left(\frac{{p}_{{\Vert}}}{B}\right)}^{2}+e{\phi }^{\prime }\right)\left(I\frac{\partial {\tilde {A}}_{\varphi }^{\prime }}{\partial \theta }-g\frac{\partial {\tilde {A}}_{\theta }}{\partial \theta }\right)\right.\\ & \quad \left.-\frac{1}{{D}_{r}}\frac{{B}^{2}}{\gamma m}\frac{{p}_{{\Vert}}}{B}\left(\left(\frac{{p}_{{\Vert}}}{B}{g}^{\prime }-\frac{e}{c}+\frac{e}{c}{\tilde {A}}_{\varphi }^{\prime }\right)\frac{\partial {\tilde {A}}_{\theta }}{\partial \theta }\right.\right.\\ & \quad \left.\left.-\left(\frac{{p}_{{\Vert}}}{B}{I}^{\prime }+\frac{e}{c}q+\frac{e}{c}{\tilde {A}}_{\theta }^{\prime }\right)\frac{\partial {\tilde {A}}_{\varphi }^{\prime }}{\partial \theta }\right)\right)\\ {\dot {P}}_{\varphi }& =-\frac{1}{\gamma }\left(\mu +\frac{B}{m}{\left(\frac{{p}_{{\Vert}}}{B}\right)}^{2}\right)\frac{\partial B}{\partial \varphi }-e\frac{\partial \phi }{\partial \varphi }\\ & \quad -\frac{1}{{D}_{r}}\left(\frac{{B}^{\prime }}{\gamma }\left(\mu +\frac{B}{m}{\left(\frac{{p}_{{\Vert}}}{B}\right)}^{2}+e{\phi }^{\prime }\right)\left(I\frac{\partial {\tilde {A}}_{\varphi }}{\partial \varphi }-g\frac{{\tilde {A}}_{\theta }}{\partial \varphi }\right)\right.\\ & \quad \left.-\frac{1}{{D}_{r}}\frac{{B}^{2}}{\gamma m}\frac{{p}_{{\Vert}}}{B}\left(\left(\frac{{p}_{{\Vert}}}{B}{g}^{\prime }-\frac{e}{c}+\frac{e}{c}{\tilde {A}}_{\varphi }^{\prime }\right)\frac{\partial {\tilde {A}}_{\theta }}{\partial \varphi }\right.\right.\\ & \quad \left.\left.-\left(\frac{{p}_{{\Vert}}}{B}{I}^{\prime }+\frac{e}{c}q+\frac{e}{c}{\tilde {A}}_{\theta }^{\prime }\right)\frac{\partial {\tilde {A}}_{\varphi }}{\partial \varphi }\right)\right).\end{align*}$$

Furthermore the change in poloidal flux can now be written as

$$\begin{align*}{\dot {\psi }}_{\mathrm{p}}& =\frac{\partial {\psi }_{\mathrm{p}}}{\partial \theta }\dot {\theta }+\frac{\partial {\psi }_{\mathrm{p}}}{\partial \varphi }\dot {\varphi }+\frac{\partial {\psi }_{\mathrm{p}}}{\partial {P}_{\theta }}{\dot {P}}_{\theta }+\frac{\partial {\psi }_{\mathrm{p}}}{\partial {P}_{\varphi }}{\dot {P}}_{\varphi }\\ & =\frac{1}{{D}_{r}}\left[\left(I\frac{\partial {\tilde {A}}_{\varphi }}{\partial \theta }-g\frac{\partial {\tilde {A}}_{\theta }}{\partial \theta }\right)\dot {\theta }+\left(I\frac{\partial {\tilde {A}}_{\varphi }}{\partial \varphi }-g\frac{\partial {\tilde {A}}_{\theta }}{\partial \varphi }\right)\dot {\varphi }\right.\\ & \quad \left.+\frac{c}{e}g{\dot {P}}_{\theta }-\frac{c}{e}I{\dot {P}}_{\varphi }\right]\end{align*}$$

and the change in parallel momentum as

$$\begin{align*}\dot {\left(\frac{{p}_{{\Vert}}}{B}\right)}& =\frac{1}{{D}_{r}}\left[\left(\left(\frac{{p}_{{\Vert}}}{B}{g}^{\prime }-\frac{e}{c}+\frac{e}{c}{\tilde {A}}_{\varphi }^{\prime }\right)\frac{\partial {\tilde {A}}_{\theta }}{\partial \theta }\right.\right.\\ & \quad \left.\left.-\left(\frac{{p}_{{\Vert}}}{B}{I}^{\prime }+\frac{e}{c}q+\frac{e}{c}{\tilde {A}}_{\theta }^{\prime }\right)\frac{\partial {\tilde {A}}_{\varphi }}{\partial \theta }\right)\dot {\theta }\right.\\ & \quad \left.+\left(1-{\tilde {A}}_{\varphi }^{\prime }-\frac{c}{e}\frac{{p}_{{\Vert}}}{B}{g}^{\prime }\right)\dot {{P}_{\theta }}\right.\\ & \quad \left.+\left(\left(\frac{{p}_{{\Vert}}}{B}{g}^{\prime }-\frac{e}{c}+\frac{e}{c}{\tilde {A}}_{\varphi }^{\prime }\right)\frac{\partial {\tilde {A}}_{\theta }}{\partial \varphi }\right.\right.\\ & \quad \left.\left.-\left(\frac{{p}_{{\Vert}}}{B}{I}^{\prime }+\frac{e}{c}q+\frac{e}{c}{\tilde {A}}_{\theta }^{\prime }\right)\frac{\partial {\tilde {A}}_{\varphi }}{\partial \theta }\right)\dot {\varphi }\right.\\ & \quad \left.+\left(q+{\tilde {A}}_{\theta }^{\prime }+\frac{c}{e}\frac{{p}_{{\Vert}}}{B}{I}^{\prime }\right){\dot {P}}_{\varphi }\right].\end{align*}$$

Setting a constraint to the perturbation of the vector potential

$$\begin{equation*}\tilde {\mathbf{A}}=\tilde {\alpha }\left(\mathbf{x},t\right)\mathbf{B},\end{equation*}$$

one finds that

$$\begin{align*}\tilde {\alpha }\left(\delta \nabla {\psi }_{\mathrm{p}}+I\nabla \theta +g\nabla \varphi \right)& !{\enspace =\enspace }{\tilde {A}}_{{\psi }_{\mathrm{p}}}\nabla {\psi }_{\mathrm{p}}+{\tilde {A}}_{\theta }\nabla \theta +{\tilde {A}}_{\varphi }\nabla \varphi \\ {\Rightarrow}{\tilde {A}}_{{\psi }_{\mathrm{p}}}& =\delta \tilde {\alpha },\\ {\tilde {A}}_{\theta }& =I\tilde {\alpha },\\ {\tilde {A}}_{\varphi }& =g\tilde {\alpha }.\end{align*}$$

Using these expressions in the equations of motion and transferring to SI units

$$\begin{align*}c& \to \frac{1}{\sqrt{{\mu }_{0}{{\epsilon}}_{0}}},\hspace{25.0pt}e\to \frac{e}{\sqrt{{\mu }_{0}{{\epsilon}}_{0}}},\\ E& \to \sqrt{4\pi {{\epsilon}}_{0}}E,\hspace{25.0pt}B\to \sqrt{\frac{4\pi }{{\mu }_{0}}}B,\end{align*}$$

we get the following final expressions for the equations of motion:

$$\begin{align*}\dot {\theta }& =\frac{1}{{D}_{r}}\left[\left(\frac{{B}^{\prime }}{\gamma }\left(\mu +\frac{B}{m}{\left(\frac{{p}_{{\Vert}}}{B}\right)}^{2}\right)+e{\phi }^{\prime }\right)\frac{g}{e}\right.\\ & \quad \left.+\frac{{B}^{2}}{\gamma m}\frac{{p}_{{\Vert}}}{B}\left(1-{g}^{\prime }\left(\frac{{p}_{{\Vert}}}{eB}+\tilde {\alpha }\right)-{\tilde {\alpha }}^{\prime }g\right)\right]\\ \dot {\varphi }& =\frac{1}{{D}_{r}}\left[-\left(\frac{{B}^{\prime }}{\gamma }\left(\mu +\frac{B}{m}{\left(\frac{{p}_{{\Vert}}}{B}\right)}^{2}\right)+e{\phi }^{\prime }\right)\frac{I}{e}\right.\\ & \quad \left.+\frac{{B}^{2}}{\gamma m}\frac{{p}_{{\Vert}}}{B}\left(q+{I}^{\prime }\left(\frac{{p}_{{\Vert}}}{eB}+\tilde {\alpha }\right)+{\tilde {\alpha }}^{\prime }I\right)\right]\\ {\dot {\psi }}_{\mathrm{p}}& =\frac{1}{{D}_{r}}\left[-\frac{1}{\gamma e}\left(\mu +\frac{B}{m}{\left(\frac{{p}_{{\Vert}}}{B}\right)}^{2}\right)\left(g\frac{\partial B}{\partial \theta }-I\frac{\partial B}{\partial \varphi }\right)\right.\\ & \quad \left.-\left(g\frac{\partial \phi }{\partial \theta }-I\frac{\partial \phi }{\partial \varphi }\right)+\frac{{B}^{2}}{\gamma m}\frac{{p}_{{\Vert}}}{B}\left(g\frac{\partial \tilde {\alpha }}{\partial \theta }-I\frac{\partial \tilde {\alpha }}{\partial \varphi }\right)\right]\\ \dot {\left(\frac{{p}_{{\Vert}}}{B}\right)}& =\frac{1}{{D}_{r}}\left[\left(I\frac{\partial \tilde {\alpha }}{\partial \varphi }-g\frac{\partial \tilde {\alpha }}{\partial \theta }\right)\left(\frac{{B}^{\prime }}{\gamma }\left(\mu +\frac{B}{m}{\left(\frac{{p}_{{\Vert}}}{B}\right)}^{2}\right)+e{\phi }^{\prime }\right)\right.\\ & \quad \left.-\frac{1}{\gamma }\left(\mu +\frac{B}{m}{\left(\frac{{p}_{{\Vert}}}{B}\right)}^{2}\right)\left(\left(1-{g}^{\prime }\tilde {\alpha }-g{\tilde {\alpha }}^{\prime }-\frac{{p}_{{\Vert}}}{eB}{g}^{\prime }\right)\frac{\partial B}{\partial \theta }\right.\right.\\ & \quad \left.\left.+\left(q+{I}^{\prime }\tilde {\alpha }+I{\tilde {\alpha }}^{\prime }+\frac{{p}_{{\Vert}}}{eB}{I}^{\prime }\right)\frac{\partial B}{\partial \varphi }\right)\right.\\ & \quad \left.-e\left(\left(1-{g}^{\prime }\tilde {\alpha }-g{\tilde {\alpha }}^{\prime }-\frac{{p}_{{\Vert}}}{eB}{g}^{\prime }\right)\frac{\partial \phi }{\partial \theta }\right.\right.\\ & \quad \left.\left.+\left(q+{I}^{\prime }\tilde {\alpha }+I{\tilde {\alpha }}^{\prime }+\frac{{p}_{{\Vert}}}{eB}{I}^{\prime }\right)\frac{\partial B}{\partial \varphi }\right)\right]\end{align*},$$$

where

$$\begin{equation*}{D}_{r}=\frac{{p}_{{\Vert}}}{eB}\left({I}^{\prime }g-I{g}^{\prime }\right)+I+qg+g{I}^{\prime }\tilde {\alpha }-I{g}^{\prime }\tilde {\alpha }.\end{equation*}$$