{'fraction_non_alphanumeric': 0.030476082777338008, 'fraction_numerical': 0.01624637189490746, 'mean_word_length': 4.6343846235531485, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Electron acceleration in a JET disruption simulation

Runaways are suprathermal electrons having sufficiently high energy to be continuously accelerated up to tens of MeV by a driving electric field [1]. Highly energetic runaway electron (RE) beams capable of damaging the tokamak first wall can be observed after a plasma disruption [2]. Therefore, it is of primary importance to fully understand their generation mechanisms in order to design mitigation systems able to guarantee safe tokamak operations. In a previous work, [3], a test particle tracker was introduced in the JOREK 3D non-linear MHD code and used for studying the electron confinement during a simulated JET-like disruption. It was found in [3] that relativistic electrons are not completely deconfined by the stochastic magnetic field taking place during the disruption thermal quench (TQ). This is due to the reformation of closed magnetic surfaces at the beginning of the current quench (CQ). This result was obtained neglecting the inductive electric field in order to avoid the unrealistic particle acceleration which otherwise would have happened due to the absence of collision effects. The present paper extends [3] analysing test electron dynamics in the same simulated JET-like disruption using the complete electric field. For doing so, a simplified collision model is introduced in the particle tracker guiding center equations. We show that electrons at thermal energies can become RE during or promptly after the TQ due to a combination of three phenomena: a first REs acceleration during the TQ due to the presence of a complex MHD-induced electric field, particle reconfinement caused by the fast reformation of closed magnetic surfaces after the TQ and a secondary acceleration induced by the CQ electric field.

Abstract. Runaways are suprathermal electrons having sufficiently high energy to be continuously accelerated up to tens of MeV by a driving electric field [1]. Highly energetic runaway electron (RE) beams capable of damaging the tokamak first wall can be observed after a plasma disruption [2]. Therefore, it is of primary importance to fully understand their generation mechanisms in order to design mitigation systems able to guarantee safe tokamak operations. In a previous work, [3], a test particle tracker was introduced in the JOREK 3D non-linear MHD code and used for studying the electron confinement during a simulated JET-like disruption. It was found in [3] that relativistic electrons are not completely deconfined by the stochastic magnetic field taking place during the disruption thermal quench (TQ). This is due to the reformation of closed magnetic surfaces at the beginning of the current quench (CQ). This result was obtained neglecting the inductive electric field in order to avoid the unrealistic particle acceleration which otherwise would have happened due to the absence of collision effects. The present paper extends [3] analysing test electron dynamics in the same simulated JET-like disruption using the complete electric field. For doing so, a simplified collision model is introduced in the particle tracker guiding center equations. We show that electrons at thermal energies can become RE during or promptly after the TQ due to a combination of three phenomena: a first REs acceleration during the TQ due to the presence of a complex MHD-induced electric field, particle reconfinement caused by the fast reformation of closed magnetic surfaces after the TQ and a secondary acceleration induced by the CQ electric field.

Introduction
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 p ). Therefore, their presence during high I p discharges in ITER (I p ≈ 15MA) has to be considered as a serious threat to PFCs [5][6][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 p which leads to discharge termination. This I 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 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 · 10 −8 MA seed is foreseen to produce RE beams having currents up to ∼ 2MA by electron avalanche in a ITER 15MA disruption being characterised by a CQ time of ∼ 50ms and mitigated by argon mixed with 7kPa · m 3 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 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 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, [3] numerically studied the electron confinement properties of a JET-like disruption simulation. The latter concerns JET pulse 86887 which is an I p = 2MA -B 0 = 2T Ohmic discharge where a disruption was obtained via 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 0 > 1. In addition to the full magnetic field chaoticity, the 'artificial' q 0 > 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 3 (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 1MeV to 3MeV were observed before the CQ phase in DIII-D killer pellet-induced disruptions [19]. [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).

A collision drag model for the JOREK Guiding Center tracker
The representation of the electric field used in the JOREK code is the following [20]: where Φ = R 0 B 0 u is the electric potential as a function of the stream function u, the magnetic axis major radius R 0 and the reference toroidal magnetic field B 0 , ψ is the poloidal magnetic flux, e φ is the unit vector in the geometrical toroidal direction and R is the major radius. In a previous work [3], the inductive term ( ∂ψ ∂t ) of Eq.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 ∂ψ ∂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: where X is the GC position vector, p is the GC momentum parallel to the magnetic is the magnetic moment [21], B is the magnetic field intensity, b = B B is the magnetic field direction, B * = p ∇ × b + qB 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 (D 2 ) [23]: where p is the particle momentum in 3D momentum space, 0 is the vacuum permittivity, E 0 = mc 2 is the electron rest energy, n = n ef = n if and n D 2 are respectively the background plasma and molecular deuterium impurity number densities, Z nucl,D 2 is the deuterium nuclear charge which is set equal to 2 under the assumption of simultaneous collision with the two D 2 nuclei, λ D is the Debye length (λ D = 0 k B TeT i q 2 n(Te+T i ) ), r e = q 2 /(4π 0 E 0 ) is the classical electron radius [24] and I z = 15.5eV is the D 2 ionisation energy taken from [25]. It has to be remarked that in Eq.3c and 3d the quasi-neutrality assumption n = n i = n e of the JOREK MHD model is used (where n i and n e are respectively the plasma ion and electron densities).
The introduction of the drag force (Eq.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 Eq.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: where α e and α i are defined by Eq.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 Eq.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.

Effective electric field
In order to assess whether particles are accelerated or decelerated, the critical quantity is what we call here the effective electric field, which represents the net parallel force due to the action of the parallel electric field (E ) and collision drag (for which the particle kinetic energy and pitch angle are uniform and constant free parameters), divided by the particle charge |q|: where B * is the effective magnetic field defined in Section 2, E denotes the electric field and F ,coll is the collision drag defined in Equation 4b. It should be remarked that the parallel momentum equation 4a can be rewritten in terms of effective electric field as follows: We point out that the B * B * ·b · p ∂b ∂t and B * B * ·b · µ∇B γ GC terms have small effects on the electron distribution spreading in momentum space during the TQ. Thus, they are not included in the E eff definition. In Tables 1 the E eff for energies and pitch angle respectively of [1, 10, 100]keV and 170 • are reported at different times in the disruption simulation.  A column-wise reading of Tables 1 shows the E eff evolution for the times: [3.55, 3.83, 4.03, 6.94]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 eff time slices. On the other hand, a row-wise scan allows comparisons between different energies. A Poincaré plot is also provided for each time. For a good understanding of the E eff plots, it should be mentioned that in JET the plasma current (I 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 eff cause respectively the acceleration or deceleration of runaways. Let us describe the temporal evolution of E eff . The first row of Table 1 shows E eff 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 eff evolves from a fully decelerating to an almost fully accelerating condition when the kinetic energy is augmented from 1keV to 100keV. 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 100keV 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 1keV collision drag (left plot) and the parallel electric field (right plot) for the same simulation time (3.55ms). 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 eff . The parallel electric field variation with the kinetic energy has important effects on the electron population dynamics. If electrons are initialised in thermal conditions (E kin 1keV) the collision drag will prevent their acceleration. On the other hand, electrons from the far tail of the thermal distribution (E kin ≥ 10keV) 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 effective electric field is completely dominated by the E 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.55ms (i.e. during the pre-TQ phase).
The third row of figures composing Table 1 is dedicated to the E eff 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] and in [30] but they do not describe the E evolution during the TQ. As before, solutions for [1, 10, 100] keV do not differ significantly so, also in this case, the E term is dominant. The strongest electric activity is around mid-radius which is where E eff fluctuations up to ∼2kV 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.83ms) but their intensities are similar. This particular E eff 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.94ms) 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: 1keV electrons are always decelerated due to high collision braking while particles having a kinetic energy ≥ 10keV and confined in the plasma core will be accelerated and become runaway. As remarked for the 1keV 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 is dominated by the ∂ψ ∂t (inductive) term, which is related to the plasma current decay caused by the increase of plasma resistivity. At this point of our discussion, we present a summary of Table 1 results. In the pre-TQ (first row of Table 1) and during the CQ (second row of Table 3) phases, E eff is strongly dependent on the kinetic energy. The dynamics of a 1keV (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 10keV, regions of accelerating E appear at the plasma core allowing the generation of RE. Further increase of E 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 eff is given by E . 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.03ms) and then decay. The presence of large E eff fluctuations at the TQ and the possibility to generate RE if E kin > 10keV 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.

Electron acceleration during the TQ phase
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.01ms, i.e., just before the time of the Table 1 third row, up to the beginning of the CQ (t=6.94ms), for a total simulated time of ∼ 3ms. The late initialisation in the disruption simulation (t=4.01ms) is necessary to avoid a significant decrease of p due to the intense drag force typical of the E kin = 1keV 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 3 GC is randomly initialised on each n=0 surface identified by a specific normalised poloidal flux label (ψ) denotedψ init (n is the toroidal mode number). It has to be remarked that theψ 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 kin = [1, 5, 10, 25, 50]keV while the pitch angle is set to 170 • (counter-current passing particle) for each run. Then, the electron distribution is evolved in the disruption simulation using a time step of 14·T gyro (where T 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 3 particles initialised in the core region with an initial kinetic energy of 1keV are displayed for the first 0.3ms 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 1MeV. 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, 1keV particles which are not accelerated during the TQ cannot reach runaway conditions during the CQ (last row-left plot of Table 1). This means that the simulated runaway electrons interact with regions of counter-I p accelerating E eff 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.
The intense p 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 < 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.
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.015ms (top left), t=4.06ms (top middle), t=4.11ms (top right) and t=4.21ms (bottom middle). Red, blue and green dots denote respectively electrons having final kinetic energy above 1MeV, below 1MeV 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 φ = 180 • plane as reference, a toroidal angle interval of ±30 • and a time window of ±0.005ms. 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. Figures 7 and 8 report the RE final radial distribution for a set of initial kinetic energies and initial radial positions respectively. 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 (ψ ≤ 0.2). This beam-like focusing happens independently from the particle initial position as shown in Figure 8 for the 1keV case. The RE near-to-magnetic-axis positioning found in our simulation seems to be in agreement with observations obtained from DIII-D [31] 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 kin are discussed. Figure 9 upper and lower plots present respectively the fraction of electrons having final kinetic energy above 1MeV 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 Figure 7: Radial positions at t=6.94ms for electrons having E kin > 1MeV (averaged over allψ init ). Black, blue, red, magenta and green lines are respectively associated to initial kinetic energies of [1,5,10,25,50]keV. Figure 8: Radial positions at t=6.94ms for electrons having E kin > 1MeV and being initialised at E kin = 1keV. Each line colour is associated to a specific initial radial position. electron population reaches runaway conditions. In particular, about ∼ 1% of the initial thermal (1keV) population (we recall that the pre-TQ central electron temperature is 1keV) 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.
A second remark on Figure 9 is on the augmentation of the runaway and lost electron fractions with E kin . While the first one is related to the decrease of collision drag Figure 9: Fraction of electrons having final kinetic energy above 1MeV (top plot) and lost to the PFCs (bottom plot). Black, blue, red, magenta and green lines are respectively associated to initial kinetic energies of [1,5,10,25,50]keV. 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 ψ init = 0.05, E kin = 10keV 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 eff which negatively afflict the particle transport. Finally, it has to be recognised that for E kin < 10keV 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 10keV due to the electron smaller sensitivity to the collisions.

Possible factors causing a runaway electron over-generation in the JOREK simulated disruptions
Three different hypotheses to explain the RE over-generation observed in Section 4 are discussed in this Section. The first one (Subsection 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] report 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 Subsection 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 Subsection 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.

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 (Eq.3a of Section 2). In order to take into account electron-W interactions, the following terms have to be added respectively to Eq.3c and Eq.3d of Section 2 [22]: where n eb,W and n W are respectively the W bounded electron and ion number densities, < Z W > is the average W charge state, Z nucl,W = 74 is the W nuclear charge [25], I W = 7.86403eV is the W ground ionisation energy [25], λ D , r e and E 0 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 kin ) varying from 1keV to 10keV and a pitch angle of θ = 170 • on a constant and uniform W distribution (n W = 10 16 (m −3 )) is calculated using Eq.4b of Section 3 (where α e and α i are substituted with α eb,W and α i,W ) and then used to evaluate the W-induced E eff (E eff,W ). The plasma parameters used for calculating λ D are the pre-TQ ones (n = 3 · 10 19 (m −3 ) and T = 1.26keV) for the TQ phase and n = 6 · 10 19 (m −3 ) and T = 270eV for the CQ [23]. The bounded electron number density is estimated as n eb = n w (Z nucl,W − < Z >). A scan in < Z > from neutral (< Z >= 0) to fully ionised (< Z >= 74) impurities shows that the largest W contribution to the E eff is of: E eff,W = 9.6(V/m) for the E kin = 1keV case and of E eff,W = 1.11(V/m) for 10keV electrons. A comparison of these estimates with the E eff 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 eff (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 eff is orders of magnitude higher than E eff,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 60keV (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 eff , 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.

Estimations of the electron drag force reduction due to slow MGI-induced plasma density increase
In this Subsection, we analyse the slower-than-in-experiment plasma density increase as a possible explanation for the RE over-production. This is justified by the evidences reported in [17] (Figure 2) which show that the measured augmentation of the line integrated density (related to the assimilation of the MGI neutral gas) is higher in experiments than in JOREK simulations. For doing so, a scan of the background plasma contribution to the E eff as a function of its density is performed using Eqs.3a, 3c and 3d of Section 2 (the TQ and CQ reference plasma temperatures are respectively of 1.26keV and 270eV while the plasma densities are of 3 · 10 −19 (m −3 ) and of 6 · 10 −19 (m −3 ), the deuterium density is set to be zero). This scan reveals that a complete suppression of the TQ driving E eff (Table 1 third row) for an electron kinetic energy of 1keV and a pitch angle of 170 • can be attained adding a drag force corresponding to approximately a hundred times higher plasma density than the pre-TQ one (taken at the magnetic axis). On the other hand, the addition of a roughly 7.5 times higher density during the CQ (and with respect to the CQ one) is sufficient to completely suppress the RE generation (this estimate is obtained considering an electron kinetic energy of 60keV which is higher than the maximum E kin reported in the Figure 4 upper plot). While a hundred times increase of the core plasma density is unlikely to happen, it is reasonable to think that an increase of the CQ plasma density would significantly reduce the fraction of the initial electron population running away.

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 eff and RE generation. This is motivated by the impossibility of obtaining JET disruption simulations with realistic input parameters. Indeed, the plasma resistivity (η) 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 η 0 (noting that in JOREK the Spitzer-like plasma resistivity model η(T ) = η 0 (T /T 0 ) −1.5 is used, where T 0 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 η 0 = 3.85 · 10 −7 (Ωm) to η 0 = 3.85 · 10 −6 (Ωm) and η 0 = 3.85 · 10 −5 (Ω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 1ms for η 0 = 3.85 · 10 −6 (Ωm) and 0.28ms for η 0 = 3.85 · 10 −5 (Ωm)).
In Table 2 the Poincaré and E eff plots for the JOREK simulations obtained using plasma resistivities of 3.85 · 10 −5 (Ωm), 3.85 · 10 −6 (Ωm) and 3.85 · 10 −7 (Ω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 eff (second column of Table 2). Indeed, in all three cases E eff has a cellular topology composed of poloidally alternating accelerating and decelerating cells. However, while the E eff magnitude is similar for the η 0 = 3.85 · 10 −6 (Ωm) and η 0 = 3.85 · 10 −7 (Ωm) cases, it is much smaller for η 0 = 3.85 · 10 −5 (Ωm). One can also note the presence of an accelerating electric field at the plasma core for the η 0 = 3.85 · 10 −5 (Ω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 η 0 = 3.85 · 10 −6 (Ωm) and of η 0 = 3.85 · 10 −7 (Ωm). The magnetic and E eff topologies during the CQ, respectively third and fourth columns of Table 2, vary significantly with η 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 eff . 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 p decay induced by the lower resistivity, which is visible in Figure 10 displaying the experimental and JOREK simulated I p traces. Summarising, Table 2 shows that, in this case of study, the CQ E eff is strongly sensitive to the choice of η 0 . Conversely, initial resistivity variations weakly influence the CQ magnetic configuration and the TQ fields, Poincaré TQ E eff = 1keV TQ Poincaré CQ E eff = 1keV CQ Table 2: Poincaré and E eff ( V m ) (E kin = 1keV, θ = 170 • ) plots of the TQ and CQ phases for simulations with different plasma resistivities: top, middle and bottom plots refer respectively to η 0 of 3.85 · 10 −5 (Ωm), 3.85 · 10 −6 (Ωm) and 3.85 · 10 −7 (Ωm). especially if we compare the η 0 = 3.85 · 10 −6 (Ωm) and η 0 = 3.85 · 10 −7 (Ωm) cases.
Analysing Figure 10, one can see that the experimental I 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 η 0 = 7.7 · 10 −8 (Ωm). For this purpose, we estimated the current decay rate via linear regression of the CQ I p profiles exhibited in Figure 10 and, then, we extrapolated them to η 0 = 7.7 · 10 −8 (Ωm) via logarithmic fitting. The result is reported in Figure 11 which suggests that the experimental I p decay rate would be recovered if the realistic resistivity could be used. As a side remark on Figure 11, the JOREK dIp dt (black line) is not directly proportional to η 0 . We presume that this is due to the CQ plasma temperature increasing with η 0 due to a larger Ohmic heating.
The fraction of electrons becoming RE (averaged over the initial positions) as a function of η 0 for different initial E kin is given in Figure 12. This figure shows that Figure 12: Fraction of electrons having final kinetic energy above 1MeV as a function of the initial plasma resistivity for different initial kinetic energies (averaged over all initial radial positions). Black, blue, red, magenta and green lines correspond respectively to initial E kin of 1keV, 5keV, 10keV, 25keV and 50keV a decrease of η 0 reduces the number of produced RE. In particular, the fraction of RE generated from the 1keV populations drops from ∼ 37% for η 0 = 3.85 · 10 −5 (Ωm) to ∼ 2.6% and ∼ 0.6% respectively for η 0 = 3.85 · 10 −6 (Ωm) and η 0 = 3.85 · 10 −7 (Ωm). However, it seems that for η 0 = 7.7 · 10 −8 (Ω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 η 0 ∈ [3.85 · 10 −7 , 3.85 · 10 −6 ](Ωm) weakens when the E kin is increased.

Summary and conclusions
In order to study the generation of fast electrons, a drag force, modeling collisions between relativistic electrons and an 'impure' (presence of D 2 neutrals) background plasma, is introduced in the GC pusher of the JOREK fast particle tracker.
After having introduced the effective electric field (electric force plus drag force), we analysed its evolution during the treated disruption simulation for kinetic energies of E kin = [1, 10, 100]keV. Before (t=3.55ms) and after (t=6.94ms) the TQ, the 1keV E eff is dominated by the drag force whereas at higher energies a transition towards an E dominated E eff is observed. In contrast, throughout the TQ, E dominates the drag force independently from the initial E kin . During this phase, the E eff fluctuations reach intensities up to ∼2kV/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 activity taking place during the TQ causes an important spreading of the momentum space particle distribution in counter and coplasma current directions. Considering the counter-I p accelerated particles, a few % of them reach kinetic energies at which E eff 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 kin are generally lost to the wall while at low E 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 η 0 conducted in order to asses its importance on both MHD and particle dynamics. Indeed, simulations show that, while the TQ E eff is weakly affected by η 0 variations (in agreement with the discussions reported in [30] and references therein), the CQ E eff (at E kin = 1keV) varies significantly with η 0 , i.e., for the very high resistivity of η 0 = 3.85 · 10 −5 (Ωm) a strong accelerating E is the dominant contribution to E eff whereas for η 0 = 3.85 · 10 −7 (Ωm) the collision drag dominates. An extrapolation of the JOREK results towards a realistic resistivity suggests that the experimental CQ dIp dt , 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 present a weaker MHD activity than in experiment as can be deduced by the small I 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 simulation and experiment.
Despite the quantitative mismatch between simulations and experiment, 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 vs limiter) and on the type of MGI gas. In addition, test electron studies will be also repeated for the JOREK shattered pellet injection simulations [32] 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]. 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 [33] but it will also allow the evaluation of the electron transport in physical space due to collision scattering.

Acknowledgments
The authors wish to thank Daan van Vugt, Gergely Papp and Allen H. Boozer for the fruitful discussions.
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.