Electron acceleration in a JET disruption simulation

In tokamak plasmas, fast electrons are said to run away when collision effects are not capable to compensate the electron acceleration induced by a driving electric field [4]. Highly energetic runaway electron (RE) beams are sometimes experimentally observed during plasma disruptions with possible harmful consequences for the reactor first wall [2]. Indeed, a plasma facing component (PFC) struck by RE might suffer damages due to the deposition of high heat loads [2]. The energy carried by these suprathermal electrons can considerably increase with the plasma current (${I_{\rm p}}$). Therefore, their presence during high ${I_{\rm p}}$ discharges in ITER (${I_{\rm p} \approx 15}$ MA) has to be considered as a serious threat to PFCs [5–7]. For this reason, RE prevention and mitigation systems for the ITER tokamak are under advanced state of design [5, 6]. Ideally, their development should be based on the complete understanding of the physics underlying the formation and dissipation of REs. Unfortunately, complete answers to questions concerning the mechanisms and plasma configurations allowing the generation of a disruptive initial (primary) RE seed have not been achieved yet.

Generally, two consequent phases characterise a tokamak disruption: the first one is called thermal quench (TQ) and consists of a fast (few milliseconds) and almost complete loss of the plasma thermal energy. The second one, known as current quench (CQ), is identified by a decrease of ${I_{\rm p}}$ which leads to discharge termination. This ${I_{\rm p}}$ reduction is imputable to the very large post-TQ plasma resistance [8] and induces a toroidal driving electric field. At the end of the CQ, disruptive runaway beams may be observed via a slowly decaying ${I_{\rm p}}$ which correlates with different radiation measurements such as: synchrotron, soft and hard x-rays, gamma and neutron emissions [2, 9]. The disruption runaway phase, which is not systematically seen in experiments, is known as runaway plateau [2, 9] and can last up to few hundreds of milliseconds.

The generation mechanism of disruptive REs can be decomposed into two different families: the first one is the primary generation which consists of all the processes capable to produce REs without requiring an already existing relativistic electron population [10]. On the other hand, the second one, called the secondary generation or electron avalanche [11, 12], exponentially amplifies an already existing primary RE seed via knock-on collisions between thermal and relativistic charge carriers [11, 12]. In ITER, the electron avalanche is thought to be one of the dominant generation processes but, unfortunately, estimations of the maximum runaway beam current are difficult to obtain. One of the reasons of such a high task complexity is the significant secondary generation sensitivity to the primary RE current and by large incertitudes on the prediction of these last. This is justified by the ITER large avalanche amplification factor which exponentially magnify the errors coming from RE seed current estimations to RE beam scales. As an example of the expected exponential growth, a  ∼$2\cdot 10^{-8}$ MA seed is foreseen to produce RE beams having currents up to  ∼2 MA by electron avalanche in a ITER 15 MA disruption being characterised by a CQ time of  ∼50 ms and mitigated by argon mixed with 7 kPa · m³ of deuterium [10]. These considerations highlight the importance to better understand the processes underlying the primary RE generation in order to achieve a more effective disruption mitigation system design.

At the moment, four different mechanisms are identified as belonging to the primary generation: Dreicer, Hot Tail, tritium β-decay and the Compton scattering [10]. The Dreicer mechanism [1, 4], which is strictly related to the subject of this paper, consists in the acceleration of thermal electrons up to relativistic energies due to the presence of a driving electric force stronger than the average Coulomb collision drag. In disruptions, this electric field is generally associated to the one caused by the CQ ${I_{\rm p}}$ decay [10] but it may also be induced by the virulent MHD activity taking place during the TQ. The CQ Dreicer generation is foreseen not to significantly contribute to the formation of REs in ITER [10] while the TQ one was not taken into account in previous work. Differently, the Hot Tail mechanism involves the acceleration of a supra-thermal electron distributions emerging from the TQ phase which is thought to be caused by the incomplete thermalisation of the pre-TQ ones [13, 14]. The Hot Tail is foreseen to be one of the main primary generation mechanisms in ITER mitigated disruptions [10]. The last two processes, which are respectively the hot electron emission due to tritium β-decay and the thermal electron acceleration caused by collisions with energetic photons emitted by activated wall components, have not been experimentally observed yet but they are thought to contribute to the generation of primary REs in ITER [10].

As discussed above, in ITER the ${I_{\rm p}}$ decay induced electric field is foreseen to be small enough to prevent the CQ Dreicer generation [10] but little is known about the possibility of obtaining RE due to the strong MHD activity taking place during the TQ. This possibility requires, at least, the satisfaction of two conditions: the presence of an electric field strong enough to drive electrons up to high energies and a sufficiently long electron confinement time. The second condition was already addressed theoretically in [15, 16] and numerically in [3]. Indeed, in [15, 16] it is shown that if the plasma magnetic topology during the TQ presents respectively residual closed flux surfaces at the edge or flux tubes non intercepting the first wall, electrons can be confined long enough to become RE. On the other hand, in [3] a numerical study of the electron confinement properties in a JET-like disruption simulation is presented. The latter concerns JET pulse 86887 which is an $I_{\rm p}=2~{\rm MA} - {B_0=2}$ T Ohmic discharge where a disruption was obtained via ${\rm D_2}$ massive gas injection (MGI). In this simulation, the MGI destabilises an MHD 'modal cascade' [17] from a large 2/1 mode up to the complete magnetic field stochastisation obtained setting q₀  >  1. In addition to the full magnetic field chaoticity, the 'artificial' q₀  >  1 setting has also the effect to suppress the internal kink mode observed in [17]. It was found in [3] that electrons are able to 'survive' the TQ for a wide range of initial energies and radial positions. All these elements sustain the possibility that (fast) particles are not totally deconfined by the strong MHD activity of the TQ. A pioneering study of the REs dynamics in tokamak disruption is presented in [18] where the confinement properties of electrons having initial energy above the RE threshold are analysed via combined NIMROD MHD-particle tracking simulations. The results presented in [18] also confirm the incomplete fast electron deconfinement during the TQ indeed, the RE confinement time is found to increase with  ∼R³ (where R is the tokamak major radius). Despite the presence of both accelerating and decelerating terms in the NIMROD particle model, investigations on the electron dynamics and, more precisely, on possible RE acceleration mechanisms during a disruption TQ have not been reported yet. From the experimental side, losses of RE within energies from 1 MeV to 3 MeV were observed before the CQ phase in DIII-D killer pellet-induced disruptions [19]. James et al [19] relates these 'prompt' losses to the RE generation due to the existence of a high loop voltage in the TQ probably induced by a significant increase of hyper-resistivity. Moreover, [19] shows that the spatial deposition of RE losses seen during these experiments is consistent with NIMROD disruption simulations.

For these reasons, the present work explores the electron dynamics during the TQ phase of the JOREK simulated JET pulse 86887 disruption via the fast particle tracker already introduced and exploited in [3]. In contrast with [3], we focus our investigation on assessing the possibility for thermal electrons to become RE due to the presence of acceleration mechanisms all along the disruption TQ phase.

In order to do so, a guiding center (GC) collision drag force is introduced into the GC model presented in the section 3 of [3]. This JOREK fast particle tracker development is described in section 2. The evolution of the parallel effective electric field (sum of the parallel electric force and collision drag) during the simulated disruption is studied in section 3. In section 4 an analysis of the particle behaviour during the TQ is given. In this phase, RE formation is observed due to electron acceleration by large local MHD-induced parallel electric fields. A study on scenarios and parameters of the MHD disruption simulations altering the particle acceleration processes is furnished in section 5. Conclusions are presented at the end of this paper (section 6).

The representation of the electric field used in the JOREK code is the following [20]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathbf{E} = -\nabla \Phi - \frac{\partial \psi}{\partial t}\frac{\mathbf{e}_{\phi}}{R} \label{chp4_eq_Efield} \nonumber \end{align} \tag{ 1 }$$

where ${\Phi=R_0B_0u}$ is the electric potential as a function of the stream function u, the magnetic axis major radius R₀ and the reference toroidal magnetic field B₀, ψ is the poloidal magnetic flux, $\mathbf{e}_{\phi}$ is the unit vector in the geometrical toroidal direction and R is the major radius. In a previous work [3], the inductive term ($\frac{\partial \psi}{\partial t}$) of equation (1) was neglected in order to avoid an unrealistic electron acceleration during the pre-TQ phase of a disruption. Indeed, the GC model used in [3] did not take into account the energy dissipation of test electrons due to their collisions with the background plasma, dissipation which tends to counteract the $\frac{\partial \psi}{\partial t}$-induced acceleration. In the present work, a simplified collision drag model is introduced in the GC equations presented in [3] which are reported below for sake of completeness:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \dot{\mathbf{X}} = {{\frac{1}{{\mathbf{b}} \cdot {\mathbf{B^{\ast}}}}} \left({q{\mathbf{E}}{\times}{\mathbf{b}}} - {{p_{\parallel}}{\frac{{\partial}{\mathbf{b}}}{\partial t}} \times {\mathbf{b}}} + {\frac{m{\mu}{\mathbf{b}}{\times}{\mathbf{{\nabla}B}} + {p_{\parallel}}{\mathbf{B^{\ast}}}}{m \gamma_{\rm GC}}}\right)} \label{gc_position} \nonumber \end{align} \tag{ 2a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{\dot{p_{\parallel}}} ={{\frac{\mathbf{B^{\ast}}}{\mathbf{b} \cdot \mathbf{B^{\ast}}}} \cdot \left(q{\mathbf{E}}-{{p_{\parallel}}{\frac{{\partial}{\mathbf{b}}}{\partial t}}}-{\frac{{\mu}{\nabla}{B}}{\gamma_{\rm GC}}}\right)}} \label{gc_momentum} \nonumber \end{align} \tag{ 2b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\gamma_{\rm GC}} = {\sqrt{1+{\left(\frac{p_{\parallel}}{mc}\right)^2} + {\frac{2{\mu}B}{mc^2}}}} \label{gamma_GC} \nonumber \end{align} \tag{ 2c }$$

where ${\mathbf{X}}$ is the GC position vector, $p_{\parallel}$ is the GC momentum parallel to the magnetic field, ${{\mu}=\frac{{\| \mathbf{p}-{p_{\parallel}}\mathbf{b} \|}^2}{2mB}}$ is the magnetic moment [21], ${B}$ is the magnetic field intensity, ${\mathbf{b}}={\frac{\mathbf{B}}{{B}}}$ is the magnetic field direction, ${{\mathbf{B^{\ast}}}={p_{\parallel}}{{\mathbf{\nabla}}\times {\mathbf{b}}}+q{\mathbf{B}}}$ is the so-called 'effective magnetic field', ${q}$ and ${m}$ are respectively the particle charge and mass while ${c}$ is the speed of light.

The drag force used in this work is the one given in [22] which is adapted to the JOREK MHD model including molecular deuterium (${\rm D_2}$) [23]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \mathbf{F} = -{\frac{q^4}{4\pi \epsilon^2_0 E_0} \frac{\gamma((\gamma + 1)\alpha_{\rm e} + \alpha_{\rm i})}{(\gamma^2 - 1)^{\frac{3}{2}}}} \frac{\mathbf{p}}{{mc}} \label{eq_chp4_drag_force_fussmann} \nonumber \end{align} \tag{ 3a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \gamma = \sqrt{1+\left(\frac{\mathbf{p}}{{mc}}\right)^2} \nonumber \end{align} \tag{ 3b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\alpha_{\rm e} = n{\;\ln{\left({\Lambda_{\rm ef}}\right)}} + n_{\rm D_2}{\;\ln{\left({\Lambda_{{\rm eb}}}\right)}} \label{chp4_eq_electron_collisional_density}} \nonumber \end{align} \tag{ 3c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\alpha_{\rm i} = n{\;\ln{\left({\Lambda_{\rm if}}\right)}} + n_{\rm D_2}{{\left({Z_{{\rm nucl},{\rm D_2}}}\right)}^2}{\;\ln{\left({\Lambda_{{\rm nucl},{\rm D_2}}}\right)}} \label{chp4_eq_ion_collisional_density}} \nonumber \end{align} \tag{ 3d }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\Lambda_{\rm ef} = \frac{(\gamma-1)\sqrt{\gamma+1}\lambda_{\rm D}}{2\gamma r_{\rm e}},} \quad {\Lambda_{{\rm eb}} = (\gamma-1)\sqrt{\gamma+1} \frac{E_0}{I_z}} \nonumber \end{align} \tag{ 3e }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\Lambda_{\rm if} = \frac{(\gamma^2 - 1)\lambda_{\rm D}}{\gamma r_{\rm e}},} \quad {\Lambda_{{\rm nucl},k} = \frac{(\gamma^2 - 1)E_0}{\gamma I_z}} \nonumber \end{align} \tag{ 3f }$$

where $\mathbf{p}$ is the particle momentum in 3D momentum space, $\newcommand{\e}{{\rm e}} {\epsilon_0}$ is the vacuum permittivity, ${E_0=mc^2}$ is the electron rest energy, ${n=n_{\rm ef}=n_{\rm if}}$ and ${n_{\rm D_2}}$ are respectively the background plasma and molecular deuterium impurity number densities, ${Z_{{\rm nucl}, {\rm D_2}}}$ is the deuterium nuclear charge which is set equal to 2 under the assumption of simultaneous collision with the two ${\rm D_2}$ nuclei, $\lambda_{\rm D}$ is the Debye length ($\newcommand{\e}{{\rm e}} {\lambda_{\rm D} = \sqrt{\frac{\epsilon_0 k_{\rm B} T_{\rm e} T_{\rm i}}{q^2n(T_{\rm e} + T_{\rm i})}}}$), $\newcommand{\e}{{\rm e}} {r_{\rm e}= q^2/(4\pi\epsilon_0E_0)}$ is the classical electron radius [24] and ${I_z=15.5~{\rm eV}}$ is the ${\rm D_2}$ ionisation energy taken from [25]. It has to be remarked that in equations (3c) and (3d) the quasi-neutrality assumption ${n = n_{\rm i} = n_{\rm e}}$ of the JOREK MHD model is used (where ${n_{\rm i}}$ and ${n_{\rm e}}$ are respectively the plasma ion and electron densities).

The introduction of the drag force (equation (3a)) in the JOREK fast particle tracker GC model is obtained neglecting its component acting on the particle perpendicular velocity (the magnetic moment remains an adiabatic invariant of motion) and the drag-induced drifts appearing in equation (2a) resulting from the GC expansion. Moreover, the plasma fields are approximated substituting the GC position to the particle one. Thus, the modified GC parallel momentum equation is:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\frac{{\rm d}p_{\parallel}}{{\rm d}t} = \frac{\mathbf{B^{\ast}}}{{\mathbf{B^{\ast}}}\cdot{\mathbf{b}}} \cdot \left({q{\mathbf{E}} - {p_{\parallel}}\frac{\partial {\mathbf{b}}}{\partial t} -\frac{\mu \nabla B}{\gamma_{\rm GC}}}\right) + F_{{\parallel},~{\rm coll}} \label{chp4_eq_parallel_momentum_drag}} \nonumber \end{align} \tag{ 4a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {F_{{\parallel},~{\rm coll}} = - {\frac{q^4}{4{\pi}{\epsilon^2_0}E_{0}}} {\frac{\gamma_{\rm GC}{\left({{\left({\gamma_{\rm GC} + 1}\right)}\alpha_{\rm e} + \alpha_{\rm i}}\right)}}{{\left({\gamma_{\rm GC} - 1}\right)}^{\frac{3}{2}}}} \frac{p_{\parallel}}{{m}c} \label{chp4_eq_collision_drag}} \nonumber \end{align} \tag{ 4b }$$

where ${\alpha_{\rm e}}$ and ${\alpha_{\rm i}}$ are defined by equations (3c) and (3d). This simplistic model is justified by the scope of the present work. Indeed, our main aim is to perform a first assessment of the running away possibility of a thermal electron due to the TQ electric field acceleration and not to precisely describe its phase space dynamics. For this reason, we preferred the faster and more intuitive drag model composed by equations (4a) and (4b) leaving as future work the implementation of a more accurate Monte Carlo solver such as the one presented in [26].

One of the drawbacks of this simple collision drag model is the absence of a 'thermal bath'. We minimise its consequences initialising particle populations just before the TQ. It has also to be remarked that neither the effects of the background plasma velocity nor the ones due to plasma fluxes are taken into account by this collision drag. Indeed, the plasma rotation was found to be significantly smaller than the electron thermal velocity in the considered discharge thus, its effects on the collision drag are neglected in the present work.

In order to assess whether electrons are accelerated or decelerated, the critical quantity to be considered is the net parallel force (${F_{\parallel}}$) acting on each particle. In the GC model used in this work, ${F_{\parallel}}$ is given by the right hand side of equation (4a):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {F_{\parallel} = \frac{\mathbf{B^{\ast}}}{{\mathbf{B^{\ast}}}\cdot{\mathbf{b}}} \cdot \left({q{\mathbf{E}} - {p_{\parallel}}\frac{\partial {\mathbf{b}}}{\partial t} -\frac{\mu \nabla B}{\gamma_{\rm GC}}}\right) + F_{{\parallel},~{\rm coll}} . \label{sec3_eq_parallel_force}} \nonumber \end{align} \tag{ 5 }$$

For the kinetic energy levels investigated hereafter, the $\frac{\mathbf{B^{\ast}}}{{\mathbf{B^{\ast}}}\cdot{\mathbf{b}}} \cdot {p_{\parallel}}\frac{\partial {\mathbf{b}}}{\partial t}$ and $\frac{\mathbf{B^{\ast}}}{{\mathbf{B^{\ast}}}\cdot{\mathbf{b}}} \cdot \frac{\mu \nabla B}{\gamma_{\rm GC}}$ terms have small effects on the electron distribution spreading in energy space during the TQ. Consequently, we focus our attention on the net force resulting from the parallel electric field and the collision drag:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {F_{\parallel,~{\rm acc}} = \left(\frac{\mathbf{B^{\ast}}}{\mathbf{B^{\ast}}\cdot \mathbf{b}} \cdot q \mathbf{E} + F_{\parallel,~{\rm coll}}\right)}. \label{sec3_eq_acc_parallel_force} \nonumber \end{align} \tag{ 6 }$$

It is worth remarking that ${F_{\parallel, {\rm acc}}}$ can be interpreted in terms of a parallel effective electric field (${E_{\parallel, {\rm ef\,\!f}}}$) defined as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {E_{\parallel,~{\rm ef\,\!f}} =\frac{1}{|q|}\left(\frac{\mathbf{B^{\ast}}}{\mathbf{B^{\ast}}\cdot \mathbf{b}} \cdot q \mathbf{E} + F_{\parallel,~{\rm coll}}\right) = \frac{F_{\parallel,~{\rm acc}}}{|q|}}. \label{chp4_effective_electric_field} \nonumber \end{align} \tag{ 7 }$$

This quantity is particularly convenient because it allows to rewrite the collisional GC parallel momentum equation (equation (4a)) into the same form of the non-collisional one (equation (2b)):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\frac{{\rm d}p_{\parallel}}{{\rm d}t} = |q|E_{\parallel,~{\rm ef\,\!f}} - \frac{\mathbf{B^{\ast}}}{{\mathbf{B^{\ast}}}\cdot{\mathbf{b}}} \cdot \left({{p_{\parallel}}\frac{\partial {\mathbf{b}}}{\partial t} +\frac{\mu \nabla B}{\gamma_{\rm GC}}}\right)}. \nonumber \end{align} \tag{ 8 }$$

Therefore, ${E_{\parallel, {\rm ef\,\!f}}}$ can be interpreted as the electric field accelerating an 'equivalent' non-collisional particle. This interpretation makes ${E_{\parallel, {\rm ef\,\!f}}}$ the most suitable quantity for evaluating the capability of a MHD field to accelerate or decelerate electrons thus, it will be extensively used in the remaining of this work.

In table 1 the ${E_{\parallel, {\rm ef\,\!f}}}$ for energies and pitch angle respectively of $[1, 10, 100]$ keV and $170^{\circ}$ are reported at different times in the disruption simulation. This pitch angle value is chosen within the typical experimental interval of $\theta \approx \left(5^{\circ}, 12^{\circ}\right)$ seen in various machines [31–33] and respects the RE counter-current motion [18, 27] (the JET plasma current and magnetic field are both in clockwise direction seen from above thus, electron counter current motions are characterised by a negative ${p_{\parallel}}$). It has also to be remarked that this setting allows comparisons between the analysis reported in the following of this paper and the one of [3].

Table 1. Parallel effective electric field (V m⁻¹) at a toroidal angle, ϕ, of $45^{\circ}$ for a pitch angle of $170^{\circ}$. From top to bottom different disruption instants are reported: pre-TQ (t  =  3.55 ms), TQ beginning (t  =  3.83 ms), fully developed TQ (t  =  4.03 ms) and CQ beginning (t  =  6.94 ms). Kinetic energies of $\left[1, 10, 100\right]$ keV are shown from left to right. Blue and red shades represent respectively regions of accelerating and decelerating ${E_{\parallel, {\rm ef\,\!f}}}$.

$E_{\rm kin}=1$ keV	$E_{\rm kin}=10$ keV	$E_{\rm kin}=100$ keV	Poincaré plot	Time and phase
				t  =  3.55 ms: pre-TQ with large 2/1 magnetic island
				t  =  3.83 ms: TQ beginning
				t  =  4.03 ms: developed TQ
				t  =  6.94 ms: CQ beginning

A column-wise reading of table 1 shows the ${E_{\parallel, {\rm ef\,\!f}}}$ evolution for the times: ${[3.55, 3.83, 4.03, 6.94]~{\rm ms}}$. These time slices correspond respectively to the pre-TQ, TQ beginning, TQ and beginning of the CQ phases as visualised in figure 1 where solid lines represent the mode magnetic energies of the simulated disruption and dash-dot lines are associated to the ${E_{\parallel, {\rm ef\,\!f}}}$ time slices. On the other hand, a row-wise scan allows comparisons between different energies. A Poincaré plot is also provided for each time. As discussed above, in JET the plasma current (${I_{\rm p}}$) and the toroidal magnetic field are in the same direction while RE move always in the opposite direction [27]. Thus, regions of negative (in blue) and positive (in red) ${E_{\parallel, {\rm ef\,\!f}}}$ cause respectively the acceleration or deceleration of runaways. Let us describe the temporal evolution of ${E_{\parallel, {\rm ef\,\!f}}}$. The first row of table 1 shows ${E_{\parallel, {\rm ef\,\!f}}}$ when the presence of neutral gas has significantly destabilised an ${m=2, n=1}$ tearing mode (magnetic island) but before the TQ onset. A comparison among figures at different energies reveals that ${E_{\parallel, {\rm ef\,\!f}}}$ evolves from a fully decelerating to an almost fully accelerating condition when the kinetic energy is augmented from 1 keV to 100 keV. This evolution is related to the reduction of collisionality with the kinetic energy increase. As remarked in [17], the highest electron density rise is localised within the ${m=2, n=1}$ magnetic island, which causes an inevitable increase of drag force and plasma resistivity in this region. At higher kinetic energies, electron collisions become less probable implying an inevitable decrease of the drag force. At the same time, the 100 keV plot reveals the presence of an accelerating electric field mainly due to the increase of plasma resistivity. This double effect due to the MGI deposition in the ${m=2, n=1}$ tearing mode is more clearly visualised in figure 2 which presents the 1 keV collision drag (left plot) and the parallel electric field (right plot) for the same simulation time (3.55 ms). When figure 2 is juxtaposed to the first row of table 1 it becomes evident that the augmentation of particle energy causes a transition from a collision to an electric field dominated ${E_{\parallel, {\rm ef\,\!f}}}$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The parallel effective electric field variation with the kinetic energy has important effects on the electron population dynamics. If electrons are initialised in thermal conditions (${E_{\rm kin} \simeq 1~{\rm keV}}$) the collision drag will prevent their acceleration. On the other hand, electrons from the far tail of the thermal distribution (${E_{\rm kin} \geqslant 10~{\rm keV}}$) may run away before the TQ, especially inside the ${m=2, n=1}$ magnetic island.

We now turn our attention to the second row of table 1, which corresponds to the beginning of the magnetic field stochastisation (in the following, we will refer to this time instant as the TQ beginning). In contrast with the previous phase, at the TQ beginning the parallel effective electric field is completely dominated by the ${E_\parallel}$ term: the figures composing the second row of table 1 are indeed virtually impossible to distinguish. In addition, it has to be remarked that the parallel electric field activity is mainly focused at the plasma core with an intensity two orders of magnitude higher than at t  =  3.55 ms (i.e. during the pre-TQ phase).

The third row of figures composing table 1 is dedicated to the ${E_{\parallel, {\rm ef\,\!f}}}$ acting during the MHD-activity peak of the TQ, when closed magnetic surfaces are completely destroyed. The characterisation of the source and type of the observed MHD fluctuations is beyond the scope of this paper but, at this stage, a plausible hypothesis involves the generation of MHD turbulence by the magnetic field ergodisation as discussed in [28] chapter 8.2.5. Further theoretical works on the subject are reported in [29, 30] but they do not describe the ${E_{\parallel}}$ evolution during the TQ. As before, solutions for ${\left[1, 10, 100\right]~{\rm keV}}$ do not differ significantly so, also in this case, the ${E_{\parallel}}$ term is dominant. The strongest electric activity is around mid-radius which is where ${E_{\parallel, {\rm ef\,\!f}}}$ fluctuations up to  ∼2 kV are observed. During this phase, the electric field presents a cellular-like topology with an alternation of accelerating and decelerating regions in the poloidal direction which extends up to the plasma edge. These cells are smaller than the ones observed at the TQ beginning (t  =  3.83 ms) but their intensities are similar. This particular ${E_{\parallel, {\rm ef\,\!f}}}$ topology suggests that it is mainly due to the large MHD fluctuations taking place during this phase but a dedicated analysis would be required to understand the precise mechanisms at play.

The last row of table 1 corresponds to the beginning of the CQ (t  =  6.94 ms) which is characterised by the presence of large areas having good confinement properties and by the beginning of the plasma current decay. The CQ stage is distinguished by the return of the competition between collision drag and accelerating electric field: 1 keV electrons are always decelerated due to high collision braking while particles having a kinetic energy  ⩾10 keV and confined in the plasma core will be accelerated and become runaway. As remarked for the 1 keV plot of the table 1 first row, the drag force is stronger at the plasma edge due to the higher MGI-induced increase in electron density. In figure 3 the parallel electric field is decomposed into scalar and vector potential components. Figure 3 shows that the ${E_{\parallel}}$ is dominated by the ${\frac{\partial \psi}{\partial t}}$ (inductive) term, which is related to the plasma current decay caused by the increase of plasma resistivity.

Zoom In Zoom Out Reset image size

At this point of our discussion, we present a summary of table 1 results. In the pre-TQ and CQ phases (respectively first and last rows of table 1), ${E_{\parallel, {\rm ef\,\!f}}}$ is strongly dependent on the kinetic energy. The dynamics of a 1 keV (thermal) electron is everywhere dominated by the drag force; thus, a thermal population cannot reach runaway energies in these time periods. Conversely for kinetic energies of 10 keV, regions of accelerating ${E_{\parallel}}$ appear at the plasma core allowing the generation of RE. Further increase of ${E_{\rm kin}}$ implies a greater drag force reduction thus stronger accelerating electric fields which extend towards the plasma edge. In contrast, all along the TQ (second and third rows of table 1), the most prominent contribution to ${E_{\parallel, {\rm ef\,\!f}}}$ is given by ${E_{\parallel}}$. During this phase, the MHD activity generates cells of accelerating and decelerating electric fields which strengthen and reduce in size until the complete magnetic field stochastisation (t  =  4.03 ms) and then decay. The presence of large ${E_{\parallel, {\rm ef\,\!f}}}$ fluctuations at the TQ and the possibility to generate RE if $E_{\rm kin} >10$ keV during the CQ raise the question of whether the TQ electric fields are able to accelerate a fraction of an initially thermal electron population up to this critical energy level leading to the formation of RE. This question motivates the study presented in section 4.

In this section we study the capability of the TQ electric field to accelerate a thermal electron up to runaway conditions. For this purpose, we track multiple test GC populations from t  =  4.01 ms, i.e. just before the time of the table 1 third row, up to the beginning of the CQ (t  =  6.94 ms), for a total simulated time of  ∼3 ms. The late initialisation in the disruption simulation (t  =  4.01 ms) is necessary to avoid a significant decrease of ${p_{\parallel}}$ due to the intense drag force typical of the ${E_{\rm kin}=1~{\rm keV}}$ case. This momentum reduction is caused by the drag force model dissipation which does not preserve the thermodynamic equilibrium. The procedure is similar to the one exposed in [3]: the plasma minor radius, expressed in normalised magnetic flux coordinates, is divided into 10 nodes from the plasma core to the edge. A population consisting of 10³ GC is randomly initialised on each n  =  0 surface identified by a specific normalised poloidal flux label (${\bar{\psi}}$) denoted ${\bar{\psi}_{\rm init}}$ (n is the toroidal mode number). It has to be remarked that the ${\bar{\psi}}$ coordinate system is chosen to be the one used for the particle initialisation procedure and is kept constant all along the simulation. As done in [3], a mono-energetic mono-pitch angle electron beam is used. The chosen energies are: ${E_{\rm kin}=[1, 5, 10, 25, 50]~{\rm keV}}$ while the pitch angle is set to ${170^{\circ}}$ (counter-current passing particle) for each run. Then, the electron distribution is evolved in the disruption simulation using a time step of ${14{\cdot}T_{\rm gyro}}$ (where ${T_{\rm gyro}}$ is the non-relativistic gyration period) which was shown to be a good compromise between result accuracy and simulation computational cost [3].

In figure 4, the kinetic energy (upper plot) and the parallel momentum (lower plot) time profiles for 10³ particles initialised in the core region with an initial kinetic energy of 1 keV are displayed for the first 0.3 ms of the simulation. Lost particle profiles are shown using green lines while red and blue lines are associated to electrons having a final energy respectively above and below 1 MeV. It has to be noted that this set of initial conditions is representative of a core background electron population. Indeed, the particle initial kinetic energy is consistent with the background thermal one at this time of the simulated disruption. Moreover, the core magnetic field is ergodised just before the TQ MHD activity peak therefore, core electrons remains confined during both the pre-TQ and the TQ beginning phases justifying the late particle initialisation discussed in section 2. The first plot of figure 4 clearly shows that a fraction of the initial population (in blue) loses its kinetic energy until reaching the minimum energy level allowed by the drag operator. Contrarily, a few electrons (in red) see an increase in their kinetic energies up to relativistic conditions, positively answering to the question asked at the end of section 3. Recalling the interpretation also given in section 3, 1 keV particles which are not accelerated during the TQ cannot reach runaway conditions during the CQ because, for this energy level, the collision drag overwhelms the driving electric field at each radial position (last row-left plot of table 1). This means that the simulated REs interact with regions of counter-${I_{\rm p}}$ accelerating ${E_{\parallel, {\rm ef\,\!f}}}$ during the TQ. Despite the beam-type initialisation, which obliges all electrons to have equal parallel momentum, when the population enters the TQ phase a significant spread of the distribution function in velocity space is recorded. This spreading is mainly caused by the presence of accelerating and decelerating electric cells shown in the second and third rows of table 1.

Zoom In Zoom Out Reset image size

The intense ${p_{\parallel}}$ fluctuations in concert with the permanent presence of particles at the plasma core (discussed in [3]) allow the confinement of electrons (less than ${2\%}$ of the initial population) having energies high enough and the correct direction (${p_{\parallel}<0}$) for becoming runaway due to the CQ inductive electric field, as depicted by the red lines of figure 4.

In order to assess the importance of MHD fluctuations in the electron acceleration process described above, the same simulation presented in figure 4 was conducted using only the n  =  0 component of the background plasma fields. The electron kinetic energy time traces for this simulation are reported in figure 5.

Zoom In Zoom Out Reset image size

From figure 5 it is clear that the electron acceleration during the TQ is strictly related to the presence of MHD fluctuations indeed, when only the n  =  0 mode is used the maximum kinetic energy variation is a hundred times smaller than the one obtained using the full background plasma fields, preventing the generation of REs.

Figure 6 reports the electron Pseudo-Poincaré and the field line Poincaré plots for the times: t  =  4.015 ms (top left), t  =  4.06 ms (top middle), t  =  4.11 ms (top right) and t  =  4.21 ms (bottom middle). Red, blue and green dots denote respectively electrons having final kinetic energy above 1 MeV, below 1 MeV and lost particles while field line positions are identified with black dots. As introduced in [3], Pseudo-Poincaré plots represent the nearest particle positions to a specific poloidal plane within a given time interval. Figure 6 is obtained using the ${\phi=180^{\circ}}$ plane as reference, a toroidal angle interval of  ±$30^{\circ}$ and a time window of  ±0.005 ms. As observed in [3], the magnetic stochasticity destroys the initial particle torus (figure 6 top left plot) spreading electrons in the whole plasma volume (figure 6 top middle plot). In this simulation, electrons dispersed outside the plasma centre do not reach runaway energies. Indeed, they are either lost to the wall (green dots) or decelerated to low kinetic energies (blue dots) probably due to the MGI-induced increase of collisionality in the plasma outer region. Conversely, those remaining close to the magnetic axis are reconfined by the reformation of magnetic surfaces at the end of the TQ chaotic phase. In this case, electrons being accelerated before reconfinement become RE (red dots) otherwise they are slowed down by collisions (blue dots). Thereby, the probability of an electron to become runaway seems to be related to the combined effects of the momentum and physical space transports.

Zoom In Zoom Out Reset image size

Figures 7 and 8 report the RE final radial distribution for a set of initial kinetic energies and initial radial positions respectively.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 7 furnishes the final RE radial distribution averaged over all initial radii in terms of normalised poloidal flux, for different initial kinetic energies. Clearly, RE are focused in the plasma core region ($\bar{\psi} \leqslant 0.2$). This beam-like focusing happens independently from the particle initial position as shown in figure 8 for the 1 keV case.

The RE near-to-magnetic-axis positioning found in our simulation seems to be in agreement with observations obtained from DIII-D [34] and, less clearly, from TEXTOR [9] RE experiments, where the measured RE synchrotron radiation suggest that the beam is mainly localised at the core of the post-disruptive plasma.

As last part of this section, the dependencies between the total number of generated RE and the electron population initial radial position and ${E_{\rm kin}}$ are discussed.

Figure 9 upper and lower plots present respectively the fraction of electrons having final kinetic energy above 1 MeV and lost to the wall as a function of the initial radial position, for a range of initial kinetic energies. A first observation of the figure 9 upper plot reveals that, independently from the initial kinetic energy, a few % of each initial electron population reaches runaway conditions. In particular, about  ∼1% of the initial thermal (1 keV) population (we recall that the pre-TQ central electron temperature is 1 keV) runs away. This is one to two orders of magnitude larger than the RE density needed for carrying the whole plasma current. Since no signs of RE were observed during the particular experiment modelled here, it seems that our model strongly over-estimates the runaway seed production. Different possible reasons for this over-prediction will be studied in the next section.

Zoom In Zoom Out Reset image size

A second remark on figure 9 is on the augmentation of the runaway and lost electron fractions with ${E_{\rm kin}}$. While the first one is related to the decrease of collision drag at higher kinetic energies, the second is probably linked to the faster particle transport shown in [3]. However, in these simulations particle losses are lower than the ones given in [3] for each initial energy level, e.g. less than $10\%$ of the ${\left\{\bar{\psi}_{\rm init}=0.05, E_{\rm kin}=10~{\rm keV} \right\}}$ distribution is lost in figure 9 against the $50\%$ reported in [3] (figure 7 upper plot). These discrepancies are probably related to the presence of regions having decelerating ${E_{\parallel, {\rm ef\,\!f}}}$ which negatively afflict the particle transport. Finally, it has to be recognised that for $E_{\rm kin}<10$ keV the number of REs reduces considerably when the initial radius increases, supporting the idea that the plasma core has the most favourable conditions for the generation of runaways, essentially because the collision drag is smaller in the core due to the smaller density. This spatial dependency considerably reduces for initial energies above 10 keV due to the electron smaller sensitivity to the collisions.

Three different hypotheses to explain the RE over-generation observed in section 4 are discussed in this section. The first one (section 5.1) concerns the absence of high-Z impurities, such as tungsten (W), in the treated disruption simulation. Indeed, high-Z pollutants coming from the metallic walls have the effect to increase the collision drag possibly reducing the electron acceleration. The second hypothesis is related to the slower density rise observed in JOREK simulations than in experiments which would lead to an underestimated drag. Indeed, [17] reports that the MGI-induced density increase obtained using the JOREK code is approximately two times smaller than the measured one. This possibility is the subject of section 5.2. Finally, an evaluation of the influence of the plasma resistivity setting used in the MHD disruption simulation on the RE generation is presented in section 5.3. This study is justified by the higher plasma resistivity used in JOREK simulations than the JET estimated one which may increase the driving electric field thus, the production of REs.

5.1. Collision drag due to high-Z pollutants into the plasma discharge

As introduced above, the presence of high-Z pollutant in the discharge may reduce or avoid the run-away process. Indeed, if W impurities pollute a background plasma, fast electrons have a non-zero probability to collide with them, resulting in an increase of the test particle drag force (equation (3a) of section 2). In order to take into account electron-W interactions, the following terms have to be added respectively to equations (3c) and (3d) of section 2 [22]:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\alpha_{{\rm eb},W} = n_{{\rm eb},W}{\;}\ln\left(\Lambda_{{\rm eb},W}\right)} \nonumber \end{align} \tag{ 9a }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\alpha_{{\rm i},W} = n_W\left[{\langle Z_W\rangle }^2{\;}\ln\left(\Lambda_{\langle Z\rangle ,W}\right) +Z_{{\rm nucl},W}^2{\;}\ln\left(\Lambda_{Z_{{\rm nucl},W}}\right) \right]} \nonumber \end{align} \tag{ 9b }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\Lambda_{{\rm eb},W} = \left(\gamma-1\right)\sqrt{\gamma+1}\frac{E_0}{I_W}} \nonumber \end{align} \tag{ 9c }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\Lambda_{\langle Z\rangle } = \frac{\lambda_{\rm D}}{{\langle Z\rangle }r_{\rm e}}\frac{I_W}{E_0}}, \quad {\Lambda_{Z_{{\rm nucl},W}} = \frac{\gamma^2-1}{\gamma}\frac{E_0}{I_W}} \nonumber \end{align} \tag{ 9d }$$

where ${n_{{\rm eb}, W}}$ and n_W are respectively the W bounded electron and ion number densities, ${\langle Z_W\rangle }$ is the average W charge state, ${Z_{{\rm nucl}, W}=74}$ is the W nuclear charge [25], ${I_W = 7.864\, 03}$ eV is the W ground ionisation energy [25], ${\lambda_{\rm D}}$, ${r_{\rm e}}$ and E₀ are respectively the plasma Debye length, the electron classical radius and rest energy as defined in section 2. Using this model, the collision drag experienced by an electron having kinetic energy (${E_{\rm kin}}$) varying from 1 keV to 10 keV and a pitch angle of ${\theta=170^{\circ}}$ on a constant and uniform W distribution (${n_W = 10^{16} ~({\rm m}^{-3})}$) is calculated using equation (4b) of section 3 (where ${\alpha_{\rm e}}$ and ${\alpha_{\rm i}}$ are substituted with ${\alpha_{{\rm eb}, W}}$ and ${\alpha_{{\rm i}, W}}$) and then used to evaluate the W-induced ${E_{\parallel, {\rm ef\,\!f}}}$ (${E_{\parallel, {\rm ef\,\!f}, W}}$). The plasma parameters used for calculating ${\lambda_{\rm D}}$ are the pre-TQ ones (${n=3 \cdot 10^{19}~ ({\rm m}^{-3})}$ and ${T=1.26~{\rm keV}}$) for the TQ phase and ${n=6 \cdot 10^{19}~ ({\rm m}^{-3})}$ and ${T=270~{\rm eV}}$ for the CQ [23]. The bounded electron number density is estimated as ${n_{{\rm eb}}=n_w(Z_{{\rm nucl}, W}-\langle Z\rangle)}$. A scan in ${\langle Z\rangle }$ from neutral (${\langle Z\rangle =0}$) to fully ionised (${\langle Z\rangle =74}$) impurities shows that the largest W contribution to the ${E_{\parallel, {\rm ef\,\!f}}}$ is of: ${E_{\parallel, {\rm ef\,\!f}, W}=9.6~ ({\rm V}~{\rm m}^{-1})}$ for the ${E_{\rm kin}=1~{\rm keV}}$ case and of ${E_{\parallel, {\rm ef\,\!f}, W}=1.11 ~({\rm V}~{\rm m}^{-1})}$ for 10 keV electrons. A comparison of these estimates with the ${E_{\parallel, {\rm ef\,\!f}}}$ given in the table 1 of section 3 shows that a possible W pollution of the background plasma has the double effect of decreasing the TQ accelerating ${E_{\parallel, {\rm ef\,\!f}}}$ (table 1 third row) and of increasing the RE energy threshold at the CQ beginning (table 1 last row). While the latter can be reduced by a factor of two, the TQ driving ${E_{\parallel, {\rm ef\,\!f}}}$ is orders of magnitude higher than ${E_{\parallel, {\rm ef\,\!f}, W}}$ thus, the absence of a W impurity background in our simulation can reduce but not completely suppress the RE over-generation observed in figure 9 upper plot. This conclusion is also supported by the estimation of the W density required to completely suppress the CQ RE production which, for an electron kinetic energy of 60 keV (from figure 4 upper plot) is approximately 30 times higher than the plausible n_W value used above. Anyway, it has to be noted that these considerations are valid only from the electron kinetic point of view. Indeed, modifications of the TQ MHD dynamics, thus to the ${E_{\parallel, {\rm ef\,\!f}}}$, in disruption simulations caused by the introduction of a high-Z impurity background are difficult to predict and may completely alter the number of REs obtained in simulations.

5.2. Estimations of the electron drag force reduction due to slow MGI-induced plasma density increase

5.3. MHD fields and RE generation dependency on the plasma resistivity used in JOREK

As concluding subject, we present a first analysis of the influence of the plasma resistivity used in disruption simulation on the magnetic field, ${E_{\parallel, {\rm ef\,\!f}}}$ and RE generation. This is motivated by the impossibility of obtaining JET disruption simulations with realistic input parameters. Indeed, the plasma resistivity ($\newcommand{\e}{{\rm e}} {\eta}$) used for obtaining the results presented above is believed to be up to an order of magnitude higher than the estimated JET one. We focus our attention on the central plasma resistivity $\newcommand{\e}{{\rm e}} \eta_0$ (noting that in JOREK the Spitzer-like plasma resistivity model $\newcommand{\e}{{\rm e}} \eta(T) = \eta_0 (T/T_0){}^{-1.5}$ is used, where T₀ is the central plasma temperature). This is one of the key parameters ruling the electric field dynamics and the reformation of closed magnetic surfaces after the TQ. For doing this, we analyse two more JOREK disruption simulations in which the resistivity is increased from $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-7}$ (${\Omega}$m) to $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-6}$ (${\Omega}$m) and $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-5}$ (${\Omega}$m), all other parameters being left unchanged. As before, electrons are initialised just before the magnetic field complete stochastisation and followed until closed magnetic surfaces are reformed (total simulation time of ≃1 ms for $\newcommand{\e}{{\rm e}} {\eta_0=3.85\cdot10^{-6}}$ (${\Omega}$m) and ≃0.28 ms for $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-5}$ (${\Omega}$m)).

In table 2 the Poincaré and ${E_{\parallel, {\rm ef\,\!f}}}$ plots for the JOREK simulations obtained using plasma resistivities of $3.85\cdot10^{-5}$ (${\Omega}$m), $3.85\cdot10^{-6}$ (${\Omega}$m) and $3.85\cdot10^{-7}$ (${\Omega}$m) (respectively upper, middle and bottom rows) during the disruption TQ (first and second columns) and CQ (third and fourth columns) phases are reported. The TQ Poincaré plots (fist column) show that the magnetic field is globally chaotic independently from the resistivity. Similarities are found in ${E_{\parallel, {\rm ef\,\!f}}}$ (second column of table 2). Indeed, in all three cases ${E_{\parallel, {\rm ef\,\!f}}}$ has a cellular topology composed of poloidally alternating accelerating and decelerating cells. However, while the ${E_{\parallel, {\rm ef\,\!f}}}$ magnitude is similar for the $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-6}$ (${\Omega}$m) and $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-7}$ (${\Omega}$m) cases, it is much smaller for $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-5}$ (${\Omega}$m). One can also note the presence of an accelerating electric field at the plasma core for the $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-5}$ (${\Omega}$m) case which is not present in the other cases. These results show that disruption simulations obtained using the JOREK code recover the insensitivity of the TQ MHD activity with respect to the plasma resistivity already observed in [30] and references therein, at least for the reasonable values of $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-6}$ (${\Omega}$m) and of $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-7}$ (${\Omega}$m). The magnetic and ${E_{\parallel, {\rm ef\,\!f}}}$ topologies during the CQ, respectively third and fourth columns of table 2, vary significantly with $\newcommand{\e}{{\rm e}} \eta_0$. Indeed, while all simulations display the reformation of closed magnetic surfaces at the center, differences are visible in the Poincaré cross sections and, more importantly, on ${E_{\parallel, {\rm ef\,\!f}}}$. The latter transits from a strong accelerating to a decelerating configuration with the decrease of plasma resistivity. This changeover is explained by the slower ${I_{\rm p}}$ decay induced by the lower resistivity, which is visible in figure 10 displaying the experimental and JOREK simulated ${I_{\rm p}}$ traces. Summarising, table 2 shows that, in this case of study, the CQ ${E_{\parallel, {\rm ef\,\!f}}}$ is strongly sensitive to the choice of $\newcommand{\e}{{\rm e}} \eta_0$. Conversely, initial resistivity variations weakly influence the CQ magnetic configuration and the TQ fields, especially if we compare the $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-6}$ (${\Omega}$m) and $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-7}$ (${\Omega}$m) cases.

Zoom In Zoom Out Reset image size

Table 2. Poincaré and ${E_{\parallel, {\rm ef\,\!f}}}$ (V m⁻¹) (${E_{\rm kin}=1~{\rm keV}, \theta = 170^{\circ}}$) plots of the TQ and CQ phases for simulations with different plasma resistivities: top, middle and bottom plots refer respectively to $\newcommand{\e}{{\rm e}} \eta_0$ of $3.85\cdot10^{-5}$ (${\Omega}$m), $3.85\cdot10^{-6}$ (${\Omega}$m) and $3.85\cdot10^{-7}$ (${\Omega}$m).

Poincaré TQ	${E_{\parallel, {\rm ef\,\!f}}=1~{\rm keV}}$ TQ	Poincaré CQ	${E_{\parallel, {\rm ef\,\!f}}=1~{\rm keV}}$ CQ

Analysing figure 10, one can see that the experimental ${I_{\rm p}}$ decay is slower than the simulated one, as could be expected from the artificial increase of resistivity in the simulations. This probably partly explains the RE overproduction observed in section 4. Thus, it is interesting to try to extrapolate the JOREK results to the realistic level of $\newcommand{\e}{{\rm e}} \eta_0 = 7.7 \cdot 10^{-8}$ (${\Omega}$m). For this purpose, we estimated the current decay rate via linear regression of the CQ ${I_{\rm p}}$ profiles exhibited in figure 10 and, then, we extrapolated them to $\newcommand{\e}{{\rm e}} \eta_0=7.7 \cdot 10^{-8}$ (${\Omega}$m) via logarithmic fitting. The result is reported in figure 11 which suggests that the experimental ${I_{\rm p}}$ decay rate would be recovered if the realistic resistivity could be used.

Zoom In Zoom Out Reset image size

As a side remark on figure 11, the JOREK ${\frac{{\rm d} I_{\rm p}}{{\rm d} t}}$ (black line) is not directly proportional to $\newcommand{\e}{{\rm e}} \eta_0$. We presume that this is due to the CQ plasma temperature increasing with $\newcommand{\e}{{\rm e}} \eta_0$ due to a larger Ohmic heating.

The fraction of electrons becoming RE (averaged over the initial positions) as a function of $\newcommand{\e}{{\rm e}} \eta_0$ for different initial ${E_{\rm kin}}$ is given in figure 12. This figure shows that a decrease of $\newcommand{\e}{{\rm e}} \eta_0$ reduces the number of produced RE. In particular, the fraction of RE generated from the 1 keV populations drops from  ∼$37\%$ for $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-5}$ (${\Omega}$m) to  ∼2.6% and  ∼0.6% respectively for $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-6}$ (${\Omega}$m) and $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-7}$ (${\Omega}$m). However, it seems that for $\newcommand{\e}{{\rm e}} \eta_0=7.7 \cdot 10^{-8}$ (${\Omega}$m) the fraction of RE would still be significant compared to the $10^{-1} - 10^{-2}\%$ required for carrying the whole plasma current. It is worth remarking that the dependence between the RE number and the initial plasma resistivity for the range $\newcommand{\e}{{\rm e}} \eta_0 \in [3.85\cdot10^{-7}, 3.85\cdot10^{-6}]$ (${\Omega}$m) weakens when the ${E_{\rm kin}}$ is increased.

Zoom In Zoom Out Reset image size

In order to study the generation of fast electrons, a drag force, modeling collisions between relativistic electrons and a background plasma containing ${\rm D_2}$ neutrals, is introduced in the GC pusher of the JOREK fast particle tracker.

After having introduced the parallel effective electric field (electric force plus drag force), we analysed its evolution during the treated disruption simulation for kinetic energies of $E_{\rm kin}=[1, 10, 100]$ keV. Before (t  =  3.55 ms) and after (t  =  6.94 ms) the TQ, the 1 keV ${E_{\parallel, {\rm ef\,\!f}}}$ is dominated by the drag force whereas at higher energies a transition towards an ${E_{\parallel}}$ dominated ${E_{\parallel, {\rm ef\,\!f}}}$ is observed. In contrast, throughout the TQ, ${E_{\parallel}}$ dominates the drag force independently from the initial ${E_{\rm kin}}$. During this phase, the ${E_{\parallel, {\rm ef\,\!f}}}$ fluctuations reach intensities up to  ∼2 kV m⁻¹ and have a topology characterised by poloidally alternated accelerating and decelerating cells. Moreover, these cells shrink in size and extend from the plasma core to the edge with time. The origin of this field is related to the strong MHD activity taking place during the TQ but, at the moment, the precise mechanisms remain to be investigated.

Then, we used test particles simulations in order to analyse the generation of fast electrons. Results show that the ${E_{\parallel}}$ activity taking place during the TQ causes an important spreading of the momentum space particle distribution in counter and co-plasma current directions. Considering the counter-${I_{\rm p}}$ accelerated particles, a few % of them reach kinetic energies at which ${E_{\parallel, {\rm ef\,\!f}}}$ remains dominated by the electric field after the TQ while remaining within the plasma core region. After the TQ, these electrons are confined by the reformation of closed magnetic surfaces and driven to RE energies during the CQ by the inductive electric field. The fate of the non-RE electrons strongly depends on the population initial energy, i.e. particles having high initial ${E_{\rm kin}}$ are generally lost to the wall while at low ${E_{\rm kin}}$ electron thermalisation is the dominant process.

In the JET 86887 disruption experiment, no RE were observed. In contrast, the JOREK particle simulation indicates a strong generation of RE even for initially thermal electron populations. Three different possible reasons explaining this discrepancy have been addressed: the absence of high-Z (tungsten) impurities in the simulation, a plasma density increase slower in simulations than in experiments and the fact that the JOREK disruption simulations are run with a plasma resistivity significantly higher than the JET estimated one. While the tungsten concentration required to completely suppress the CQ RE production is found to be unrealistically elevated, a 7.5 times higher plasma density should be sufficient to avoid the run-away process. This last result has to be considered in concert with the scan in $\newcommand{\e}{{\rm e}} {\eta_0}$ conducted in order to asses its importance on both MHD and particle dynamics. Indeed, simulations show that, while the TQ ${E_{\parallel, {\rm ef\,\!f}}}$ is weakly affected by $\newcommand{\e}{{\rm e}} {\eta_0}$ variations (in agreement with the discussions reported in [30] and references therein), the CQ ${E_{\parallel, {\rm ef\,\!f}}}$ (at $E_{\rm kin}=1$ keV) varies significantly with $\newcommand{\e}{{\rm e}} {\eta_0}$, i.e. for the very high resistivity of $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-5}$ (${\Omega}$m) a strong accelerating ${E_{\parallel}}$ is the dominant contribution to ${E_{\parallel, {\rm ef\,\!f}}}$ whereas for $\newcommand{\e}{{\rm e}} \eta_0=3.85\cdot10^{-7}$ (${\Omega}$m) the collision drag dominates. An extrapolation of the JOREK results towards a realistic resistivity suggests that the experimental CQ ${\frac{{\rm d}I_{\rm p}}{{\rm d}t}}$, thus the ${E_{\parallel}}$, would be recovered. Therefore, it is reasonable to hypothesise that the combined effect of both higher plasma density and lower ${E_{\parallel}}$ would prevent the generation of REs and, in particular, their acceleration during the CQ. Anyway, it has to be mentioned that the analysed disruption simulation probably presents a weaker MHD activity than in experiment as can be deduced by the small ${I_{\rm p}}$ spike associated to the simulated TQ (figure 10). The consequences of this are not clear and, as a consequence, further efforts are needed to reconcile simulations and experiments.

Despite the quantitative mismatch between simulations and experiments, the present work suggests that a kind of Dreicer generation might take place during the TQ and at the CQ beginning. Indeed, electrons can be accelerated by large parallel electric field associated to the TQ MHD activity and, after a prompt reconfinement due to the reformation of magnetic surfaces, become RE thanks to the subsequent acceleration induced by the CQ inductive electric field. This mechanism, which was not reported before at the best knowledge of the authors, may strongly influence the primary RE seed estimates performed for ITER. Thus, further studies are advisable in order to understand and characterise the nature of the actual TQ MHD activity of a mitigated disruption and its capabilities to generate supra-thermal electron populations.

Considering future developments, the present work suggests a multiplicity of research axes. More advanced theoretical investigations should be performed (possibly with the help of simplified numerical models) in order to better understand the physics underlying the parallel electric field dynamics and electron deconfinement during the TQ. Further numerical experiments using the JOREK code will be also performed for improving the quantitative match with the experiments, in particular, with this JET MGI case. For example, it has been found recently that simulations including an impurity background display larger plasma current spike than the one presented in this paper. Another axis of research consists in testing the ability of codes such as JOREK to qualitatively reproduce robust experimental trends like the RE existence domain in JET as a function of the toroidal magnetic field or the quantity of injected impurities [2], the dependencies on the magnetic configuration (divertor versus limiter) and on the type of MGI gas. In addition, test electron studies will be also repeated for the JOREK shattered pellet injection simulations [35] in order to assess differences and similarities with respect to the MGI cases.

Finally, the JOREK fast particle tracker could be improved by implementing a Monte Carlo operator for simulating the electron Coulomb collisions such as the one reported in [26] or in [36]. It has to be remarked that a Monte Carlo approach to collisions will not only improve the phase space dynamics description for particles having nearly thermal energies [37] but it will also allow the evaluation of the electron transport in physical space due to collision scattering.

The authors wish to thank Daan van Vugt, Gergely Papp and Allen H. Boozer for the fruitful discussions. In addition, we thank the two anonymous reviewers whose comments and suggestions helped improve and clarify this manuscript.

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 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Simulations were performed on Curie and Marconi supercomputers respectively at CEA-TGCC and CINECA (MARCONI-Fusion project). We thank these institutions for the resources granted.