On the origin of the plasma current spike during a tokamak disruption and its relation with magnetic stochasticity

A JOREK 3D non-linear MHD simulation of a disruption triggered by an argon massive gas injection in JET, which quantitatively reproduces the plasma current ( Ip ) spike (Nardon et al 2021 Plasma Phys. Control. Fusion 63 115006), is analyzed in order to investigate the origin of the Ip spike and its relation with magnetic stochasticity. The Ip spike is associated to a current density (j φ ) profile relaxation which appears to result from Shear Alfvén Wave (SAW) propagation along stochastic field lines, as proposed by Boozer (2019 Plasma Phys. Control. Fusion 61 024002; 2020 Phys. Plasmas 27 102305), possibly complemented by a macroscopic E×B flow structure. Using axisymmetric JOREK simulations involving a mean field Ohm’s law, we verify that the level of hyper-resistivity associated to SAWs is consistent with the prediction made in (Boozer 2019 Plasma Phys. Control. Fusion 61 024002; Boozer 2020 Phys. Plasmas 27 102305), which connects the Ip spike with the level of stochasticity. The relaxation comprises two main phases, the first one corresponding to a fast (0.1 ms) and almost complete j φ flattening in the q < 2 region, while the second one is longer (0.5 ms) and corresponds to a more gradual, global and incomplete j φ flattening. During the first phase, strong E×B flows develop that play a key role in mixing impurities into the core.


Introduction
A common observation during tokamak disruptions is that the plasma current I p rises by about 10% (order of magnitude for conventional aspect ratio tokamaks) over a ∼1 ms timescale before starting its decay to zero over a longer timescale [4]. The present work aims at improving our understanding of this so-called 'I p spike'. It is thought to be due to a relaxation of the current density (j ϕ ) profile driven by intense magnetohydrodynamic (MHD) activity. It can indeed be shown [5] that magnetic helicity, H ≡´⃗ B · ⃗ AdV, is approximately conserved during a fast MHD relaxation and that the release of magnetic energy at fixed H leads to an increase of the current in the relaxed region.
In a series of recent papers [2,3], Boozer hypothesizes that the j ϕ profile relaxation is mediated by Shear Alfvén Waves (SAWs) which propagate along stochastic field lines, transporting helicity with them. Based on this picture, it is argued that the timescale of the relaxation should be comparable to the time required for SAWs to travel across the plasma along stochastic field lines, τ SAW = 2π R 0 N t /V A , where V A is the Alfvén speed, R 0 is the major radius and N t is the number of toroidal turns needed for a stochastic field line to travel across the plasma. Such an insight appears potentially very useful since it may allow estimating N t , a quantity which cannot be measured directly, from the experimental I p (t) signal. Credible estimates of N t would in particular help understand and predict runaway electron (RE) generation as well as the fraction of the thermal energy conducted onto plasma facing components during the thermal quench, which both strongly depend on electron transport and losses along stochastic field lines [6,7]. One of the main aims of the present work is therefore to assess whether the insight from [2,3] is supported by 3D non-linear MHD simulations.
Of course, if one disposes of 3D non-linear MHD simulations, one may calculate electron losses directly by tracing test electrons, like done in recent years [8][9][10]. However, the credibility of such an approach depends on the credibility of the MHD simulations themselves, which is often not fully established. One particular aspect casting doubt is precisely the I p spike. Past publications, e.g. [11,12], have shown simulations producing an I p spike, but the latter had a tendency to be much smaller than what is observed. However, in a recent publication [1], we presented a JOREK simulation of a disruption triggered by an argon massive gas injection (MGI) in JET which produced a quantitatively realistic I p spike 1 . This constitutes an ideal case for our investigations.
In [13], a mean field Ohm's law which may be used to simulate fast MHD relaxations with 1D or 2D codes was proposed. In this equation, MHD fluctuations are accounted for by a hyper-resistivity term which conserves magnetic helicity. Several works have discussed the effect of this hyper-resistivity, 1 It is interesting to note that a comparable Ip spike was also obtained in the 'JET shot 1' case from [18] (see figure 4 in this paper), which is a JOREK simulation using the same JET target plasma as in the present paper but studying an 'imaginary' argon shattered pellet injection instead of the MGI used in the experiment.
which is sometimes called 'anomalous electron viscosity' (this terminology is justified by the fact that 'actual' electron viscosity gives rise to a hyper-resistivity term, see p. 204 in [14]; the latter is however usually neglected due to its smallness), see e.g. [15,16]. This model has been used for example in recent work investigating the effect of an MHD relaxation on RE generation with the DREAM code [17]. In [2,3], a formula relating the hyper-resistivity η h to N t is proposed, based on the above-mentioned argument that the relaxation timescale should correspond to τ SAW .
The paper is constructed as follows. In section 2, we analyze the JOREK 3D simulation and discuss the mechanisms responsible for the I p spike in this simulation. In section 3, we present JOREK 2D (axisymmetric) mean field model simulations. We define a reference case in which we set η h so as to match the 3D simulation in terms of I p and j ϕ profile evolution. We study the effect of various modifications in the η h settings to gain insight. In order to test Boozer's theory, we then compare the η h from the reference case with the one predicted by the above-mentioned formula when calculating N t from the 3D simulation. We also directly compare the relaxation timescale observed in the 3D simulation with τ SAW . We summarize and discuss our results in section 4.

3D simulation
The JOREK simulation considered here is referred to as 'Case D' in [1], where it is described in substantial detail already. Here we provide additional information, beginning with a description of the I p and j ϕ profile dynamics and following with a discussion of the mechanisms at play.

Ip and j ϕ profile dynamics
Figure 1 (left) shows I p (t) in the experiment and simulation. One can see that the height and rise time of the I p spike are quantitatively matched. Unfortunately, a numerical issue prevents the simulation from running beyond 6.76 ms. It can be noticed that the moderate drop in I p before the beginning of the spike (i.e. in the first 6 ms), which comes from j ϕ decaying at the edge of the plasma due to cooling by the injected argon, is under-estimated by the simulation. Other simulations (Cases A, B and C in [1]) with a larger argon source (the latter being set ad hoc in the absence of a self-consistent model for gas dynamics in JOREK) reproduce the early I p (t) evolution much better and are also more consistent with the measured lineintegrated electron density [1]. Unfortunately, none of these simulations produced a quantitatively realistic I p spike, partly due to numerical issues which prevented from pushing them far enough in time. Figure 2 shows the profile of the flux-surface-averaged current density, ⟨j ϕ ⟩, at different times in the simulation, starting at 5.84 ms, i.e. just before I p starts rising, and ending at 6.54 ms, near the time when I p reaches its maximum. Note that for the moment, we are only interested in the plain lines, which correspond to the 3D simulation. It can be seen that a fast flattening of ⟨j ϕ ⟩ takes place between 5.84 and 6.00 ms in  the region ψ 1/2 N ≲ 0.7 (ψ N is the normalized poloidal magnetic flux, ψ N = 0 corresponding to the magnetic axis and ψ N = 1 to the last closed flux surface in the initial equilibrium). This corresponds to q ≲ 2 initially. The evolution of the j ϕ crosssection in figure 3 (right) is rather spectacular (compare 5.72 and 5.99 ms). The large ⟨j ϕ ⟩ peak visible in the very center at 6.00 ms in figure 2 corresponds to the central j ϕ 'ribbon' at 5.99 and 6.03 ms in figure 3. As explained in [1], we interpret this ribbon as a consequence of the (m = 2, n = 1) (where m and n are the poloidal and toroidal mode numbers) 'ghost island' (we employ this term because the island is not welldefined anymore at this time due to magnetic stochasticity, as can be seen in the Poincaré cross-sections overlayed on the j ϕ cross-sections in figure 3) 'running into itself' in the center, which is in turn caused by the fact that the q = 2 'ghost surface' moves towards the center as ⟨j ϕ ⟩ flattens. The evolution of the q profile is shown in figure 18 of [1], where it can be seen that the relaxation leads to q being very flat and slightly above 2 in the relaxed region (ψ 1/2 N ≲ 0.7). This violent relaxation is associated to a burst in magnetic and kinetic energies in the various toroidal Fourier harmonics n ⩾ 1, as can be seen in figure 4. Figure 3 (left) shows that the amplitude of magnetic fluctuations normalized to the toroidal field, δB/B t , is of the order of several% and reaches more than 4% in certain regions during the burst. After this phase, the dynamics become less violent but MHD activity continues, leading to a global 'smoothing' of the ⟨j ϕ ⟩ profile between 6.00 and 6.54 ms, as can be seen in figure 2. An interesting feature visible in figure 4 is that magnetic and kinetic fluctuations have a dip between 6.0 and 6.1 ms. This dip corresponds to a transient reduction of fluctuations in the core, which is consistent with their drive being removed after the strong ⟨j ϕ ⟩ flattening in this region. The fact that fluctuations subsequently grow again is consistent with the reformation of a ⟨j ϕ ⟩ gradient in the core as the (slower) more global ⟨j ϕ ⟩ profile relaxation takes place (compare the magenta and green dashed profiles in figure 2).

Underlying mechanisms
As explained in the introduction, Boozer proposed that the ⟨j ϕ ⟩ profile relaxation and I p spike come from SAW activity [2,3]. It is thus important to look for signs of such activity in the 3D simulation. Identifying individual SAWs is however not trivial. A number of SAWs propagating in either toroidal direction probably coexist, and stochasticity makes the situation complicated. This said, it is clear from figure 3 that during the period ≃5.8-6.5 ms, small scale fluctuations are excited, which we interpret as a sign of SAW dynamics. These are visible for example in the cross-sections of both the E×B flow pattern (left) and j ϕ (right) at 6.03 ms. It can also be noticed in figure 5 that there is a trend towards equipartition between magnetic and kinetic energies for large toroidal mode numbers, which we interpret as a sign of dominant SAW dynamics at small scales.
It appears possible, however, that SAW propagation is not the only mechanism at play. Indeed, during the fast central relaxation phase (≃5.8-6.0 ms), a macroscopic flow structure can also be observed, characterized by 'streamers', i.e. long, densely packed, almost parallel streamlines (see figure 3 left at 5.99 ms). These streamers, whose velocity is of several tens of km s −1 , are aligned with the above-mentioned current ribbon, which is remindful of the classical Sweet-Parker reconnection picture [14], although things are more complicated here because the magnetic field is strongly stochastic. An important observation from JOREK simulations ( [1,18] as well as unpublished work) is that such a flow structure appears only in simulations producing a large (i.e. realistic) I p spike. One may thus speculate that the I p spike is to some extent related to this flow.

2D mean field model simulations
We now move to JOREK 2D (axisymmetric) mean field model simulations. In this model, the equation for the evolution of the poloidal flux ψ reads (in JOREK normalized units [19]): The first term on the righthand side is an advection term ([., .] is a Poisson bracket and u is the flow potential), the second one is the resistive term and the third is the hyper-resistive term.
We start these 2D simulations at 5.84 ms, i.e. just before ⟨j ϕ ⟩ in the core begins to flatten and I p begins to rise in the 3D simulation. We use a restart file from the 3D simulation, discarding the non-axisymmetric part of the fields. In section 3.1, we describe our reference settings for the resistivity η and hyper-resistivity η h , which allow reproducing the evolution of the ⟨j ϕ ⟩ profile and I p from the 3D simulation. In section 3.2, we discuss the sensitivity of 2D simulations to various changes in the settings of η and η h , which helps gain intuition. Finally, in section 3.3, we calculate the value of η h predicted by Boozer's formula and compare it to our reference case.

Reference settings and comparison to 3D simulation
We found that in addition to prescribing η h , it is also necessary to prescribe η in these 2D simulations since calculating η from the electron temperature T e like done in the 3D simulation would not be appropriate. This is because 2D simulations behave very differently from the 3D simulation in terms of T e evolution: T e tends to collapse (due to radiation) down to the eV level in certain locations in 2D while it remains quite homogeneous at a few hundred eV in 3D due to thermal conduction along stochastic field lines (which is not accounted for in the 2D simulations). Based on the observation that in the 3D simulation, T e is homogeneous in the bulk of the plasma with a value of a few hundred eV during the I p rise phase (see figure 19 in [1]; note that this is much lower than the predisruption T e , which is in the keV range, because most of the thermal collapse takes place before the I p spike, see [1]), we use η = 10 −7 in the bulk during this phase, i.e. between 5.84 and 6.59 ms (we give η in JOREK units here and below). This  value corresponds to 4.2 × 10 −7 Ω.m, i.e. to the Spitzer resistivity at 160 eV. In the scrape-off layer (SOL), i.e. the edge region in which field lines are connected to the wall within a short distance and not connected to the bulk (in spite of magnetic stochasticity, a clear 'frontier' exists; this perhaps surprising statement can be understood by studying properties related to the stable and unstable manifolds of the X-point [20,21]), T e is much lower, in the range of a few eV, and we thus use η = 2 × 10 −5 , which corresponds to the Spitzer resistivity at 5 eV. A hyperbolic tangent is used to smoothly connect the bulk and SOL, the complete expression for the resistivity profile being: From 6.59 ms, we assume that a global radiative collapse has taken place (which is related to the efficient mixing of argon impurities taking place during the MHD relaxation) such that T e is now in the range of a few eV throughout the plasma, and we thus use a flat η profile with a value of 2 × 10 −5 . This value is supported by the fact that it allows reproducing the experimental I p decay rate during the early current quench (CQ) phase. The η(ψ N ) profile in the different phases of the simulation is shown in figure 6.
Regarding η h , based on observations from section 2, we define three phases for which we use 3 different profiles. During Phase 1, between 5.84 and 6.09 ms, we use a large η h in the core (ψ N ⩽ 0.5): η h = 2.4 × 10 −6 in JOREK units, in order to reproduce the fast ⟨j ϕ ⟩ flattening there, while we use a 6 times lower value, η h = 4 × 10 −7 , in the outer part of the bulk plasma, and we cut-off η h in the SOL. The complete hyperresistivity profile in Phase 1 is given by and is shown in figure 6 (left). In Phase 2, between 6.09 and 6.59 ms, η h is reduced to 8 × 10 −7 in the core and maintained at 4 × 10 −7 in the outer part of the bulk plasma, with still a cut-off in the SOL. The complete profile in Phase 2 is given by and is shown in figure 6 (middle). Finally, in Phase 3, after 6.59 ms, we simply set η h = 0 everywhere, based on the assumption that MHD activity decays, as suggested by the drop in magnetic fluctuations measured by Mirnov coils, see figure 2 in [1].
Settings of this reference 2D simulation are summarized in table 1. Figure 1 (right) shows that the 2D simulation reproduces well the evolution of I p from the 3D simulation. The evolution of the ⟨j ϕ ⟩ profile is also approximately matched, as can be seen by comparing the dashed (2D) and plain (3D) lines in figure 2, although it can be observed that the flattening in the central part during Phase 1 is not perfectly reproduced. A better match would probably be possible with a more sophisticated setting for η h (e.g. with an inward-expanding large η h region), but we prefer keeping the setting simple.  Time (ms) 5.84-6.09 6.09-6.59 6.59-end

Sensitivity studies
We now test different variations of the η h settings to see their influence on 2D simulation results. We begin with investigating the effect of η h in the outer part of the bulk plasma (0.5 ⩽ ψ N ⩽ 0.95) during Phase 1. First of all, if we set η h to 0 in this region, we can see in figure 7 (black curve) that I p does not rise. This is because, even though ⟨j ϕ ⟩ still flattens fast in the core, a negative skin current appears at the edge of the flattened region, compensating the rise in plasma current in the latter. This is probably what happens during sawteeth, which typically do not produce I p spikes. Conversely, if we set η h to a large value in the outer part of the bulk plasma, e.g. the same value as in the core (2.4 × 10 −6 ), I p rises much faster than in the reference simulation, as can again be seen in figure 7 (green curve). This is because now the fast ⟨j ϕ ⟩ profile flattening is too global.  Secondly, we test the effect of keeping η h during Phase 2 (6.09-6.59 ms) the same as in Phase 1 (5.84-6.09 ms), with a large η h (2.4 × 10 −6 ) in the core. Figure 8 shows that this results in a slightly larger rise in I p during Phase 2, degrading the match with the 3D simulation. This also degrades the match in terms of ⟨j ϕ ⟩ profile evolution, leading to too much flattening in the core during Phase 2.
A third test is to maintain η h during Phase 3 (after 6.59 ms) the same as it is in Phase 2 (6.09-6.59 ms). Figure 9 shows that this results in a more rounded I p spike, which is due to the fact that with such settings, the ⟨j ϕ ⟩ profile continues to relax in Phase 3, partly counter-balancing the effect of the CQ. This more rounded I p spike shape is clearly less similar to the experimental data.
It is also instructive to see in figure 10 that starting Phase 3 earlier or later results in a different I p spike height and duration, since after the switch to a large η, resistive diffusion becomes dominant, resulting in a 'truncation' of the I p spike. This suggests that the time at which I p reaches its peak value in the experiment is the time at which T e drops to a few eV throughout the plasma.
Finally, we add that 2D simulation results are not very sensitive to the precise location of the cut-off radius of the η and η h profiles at the edge during the first two phases. This radius was chosen to be ψ N = 0.95 for both η and η h in the reference case, but changing it (whether for η or η h ) by ±0.05 does not modify the results significantly. In particular, the precise degree of overlap between the large η region (the SOL) and the region where η h is applied does not have much influence on the results. This is probably due to the fact that no sharp ⟨j ϕ ⟩ gradients are produced during the process.

Testing Boozer's formula for η h
In [2,22], Boozer discusses the relaxation of the ⟨j ϕ ⟩ profile during a fast MHD relaxation, which he models with an equation equivalent to equation (1) above. As mentioned in section 1, [2,22] also propose an estimate for the level of hyper-resistivity based on the assumption that the timescale of the relaxation should be of the order of τ SAW = 2π R 0 N t /V A . Under certain approximations, this estimate can be translated into the following expression for the hyper-resistivity η h to be used in JOREK 2D mean field model simulations, as detailed in appendix: with ρ the mass density, ρ 0 the reference mass density used in the JOREK normalization (i.e. the initial mass density at the plasma center) [19], a the minor radius and B 0 = 3 T the toroidal magnetic field. In this simulation, ρ 0 = 7.1 × 10 −8 kg m −3 . As can be seen in figure 11, ρ evolves as a consequence of the MGI, increasing first in the outer part of the plasma and then homogenizing during the relaxation, reaching a level of the order of ρ ≃ 4 − 5ρ 0 . For simplicity, we will consider that √ ρ 0 /ρ = 1/2 in equation (2). We want to test whether equation (2) corresponds to the level of η h that allows matching the 3D simulation with 2D mean field model simulations as described in section 3.2. For this, we need to estimate N t . Various methods are possible. One of them consists in calculating the field line transport coefficients in terms of the advection-diffusion model introduced in [23,24]. Note that by 'transport' here we mean how field lines move and spread radially as one progresses along them at a fixed time. The advection-diffusion model considers that the radial density of field lines is transported according to a Fokker-Planck equation, the advection and diffusion coefficients of which can be obtained by methods based on field line tracing which are described in [23] and have been implemented as a possible post-treatment of JOREK simulations. The result is shown in figure 12, which displays the field line diffusion coefficient D FL (in m 2 m −1 ) (left) and advection coefficient (in m m −1 ) A FL (right). It can be seen that field line transport is rather limited until about 5.60 ms, except at the edge where A FL is large. This is consistent with the Poincaré crosssection at 5.26 ms in figure 3 (top right) which shows stochastization only at the edge. From 5.60 until 5.80 ms, D FL grows in a region around mid-radius which expands over time. This corresponds to the growth of the (m = 2, n = 1) (ghost) island. Between 5.80 and 6.00 ms, this region expands all the way to the center, and a burst in both A FL and D FL in the core region is visible. This corresponds to the fast relaxation phase described in section 2 and modeled as Phase 1 in the above-described 2D simulations. After a relatively quiet phase between 6.05 and 6.20 ms, which corresponds to the dip in the magnetic energies visible in figure 5, A FL and D FL grow again, with now a rather homogeneous and steady level across the whole plasma until the end of the simulation, i.e. until the end of Phase 2 in the 2D simulations.
Over the whole 'active' period, D FL ≃ 0.5 − 1 × 10 −3 m 2 m −1 and A FL ≃ 1 − 2 × 10 −3 m m −1 . Since a ≃ 1 m, this means that the Péclet number [23] aA FL /D FL ≃ 2, i.e. that advection dominates over diffusion. It is not fully clear why this is so, but this may be related to the fact that the magnetic  field structure is strongly dominated by the large (m = 2, n = 1) mode. This is different from a stochastic field that would result from the overlap of many small island chains, for which one would expect a diffusive behavior.
The typical distance along a field line required to travel across the minor radius is a 2 /D FL ≃ 1000-2000 m (resp. a/A FL ≃ 500-1000 m) when only diffusion (resp. advection) is taken into account. Combining diffusion and advection, we estimate the required distance to be ≃500 m, which corresponds to N t ≃ 25 since R 0 ≃ 3 m.
In order to cross-check this estimate, we also used a more direct method based on the following procedure. We initialized 200 field lines at a given (R, Z) position and at equally spaced toroidal positions, and tracked them over 40 turns. We then calculated their radial density after a given number of turns. For this, we defined 20 radial bins (i.e. ψ 1/2 N intervals) and calculated the approximate field line density in each bin by dividing the number of field lines in the bin by ψ 1/2 N , the latter being roughly proportional to the bin's volume. An example result is shown in figure 13 (left), where it can be seen that the field line distribution is initially peaked but quickly spreads over the whole domain. We performed such simulations for various initial positions and times. Figure 13 (center and right) shows how the average and variance of the radial position (ψ 1/2 N ) of remaining field lines depends on the number of toroidal turns for these different simulations. It can be seen that all curves essentially merge, i.e. initial conditions are forgotten, after 10-20 turns, at which point the variance saturates at ≃0.25, which corresponds to a spreading over the whole volume. This analysis thus suggests that N t ≃ 10-20, which is of the same order, although somewhat smaller, as our above estimate N t ≃ 25. It is not clear which estimate is more accurate, but we are mainly concerned with orders of magnitude here.
Plugging N t = 20 into equation (2), one obtains η h ≃ 6 × 10 −7 . Looking back at the settings of the reference 2D case described in section 3.2 (see table 1), we see that this is of the same order as the 4 × 10 −7 used in the outer part of the bulk plasma in Phases 1 and 2 and the 8 × 10 −7 used in the core in Phase 2, while it is 4 times smaller than the 2.4 × 10 −6 used in the core in Phase 1. However, as discussed above, N t may be somewhat smaller than 20, perhaps by a factor 2, especially during Phase 1 (see the blue curves in figure 13), which would reduce the mismatch. In addition, during the early part of Phase 1, the density has not yet increased much in the plasma core (see figure 11), such that √ ρ 0 /ρ is larger than the 1/2 we have assumed, improving the match to the theoretical prediction, although the latter still gives a too small value in this phase. Nonetheless, altogether, Boozer's estimate is consistent with our simulation results, at least in terms of orders of magnitude.

Comparison of the relaxation timescale with τ SAW
A more direct test of Boozer's theory is to compare the observed relaxation timescale, τ relax , with τ SAW . Judging from the timescale of the ⟨j ϕ ⟩ profile evolution, we may estimate that τ relax ≃ 100 µs. On the other hand, if we consider that the I p spike has a shape of the type 1 − exp(−(t − t 0 )/τ relax ), we find a larger value, τ relax ≃ 200 µs. The discrepancy between these two values is related to the dynamics of the negative current induced near the edge during the relaxation, which is visible in figure 3 at 5.99 and 6.03 ms and makes I p evolve slower than the 'actual' relaxation timescale. We thus consider that τ relax ≃ 100 µs is the better estimate. Taking ρ = 4ρ 0 ≃ 3 × 10 −7 kg m −3 , we have V A ≃ 5 × 10 6 m s −1 , leading to τ SAW ≃ 75 µs for N t = 20, which is indeed close to τ relax .

Summary and discussion
In summary, it appears plausible that the ⟨j ϕ ⟩ profile relaxation (and hence the I p spike) observed in the above-described JOREK 3D non-linear MHD simulation of a disruption triggered by an argon MGI in JET is mediated to a large extent by SAWs propagating along stochastic field lines, as proposed by Boozer [2,3]. We find in particular that the relaxation timescale is comparable to τ SAW and that the level of hyperresistivity needed in 2D (axisymmetric) mean field JOREK simulations to match the 3D simulation is consistent with Boozer's prediction.
The relaxation observed in the 3D simulation comprises two main phases. During Phase 1 (5.84-6.09 ms), a fast and almost complete flattening of ⟨j ϕ ⟩ takes place in the q < 2 region, while during Phase 2 (6.09-6.59 ms), a slower, more global flattening takes place, which remains incomplete by the end of the simulation.
An important remark is that the existence of Phase 1, i.e. of a fast almost complete relaxation in the q < 2 region, is a characteristic feature of JOREK simulations which produce a 'large' (i.e. comparable to the experiment) I p spike, as compared to those which produce only a small I p spike.
We interpret the small scale fluctuations of the E×B velocity observed in the 3D simulation as a signature of SAW dynamics in the stochastic field. As noted in section 2, in addition to these small scale fluctuations, a macroscopic E×B flow pattern is visible in the core during Phase 1, which seems to also contribute to the ⟨j ϕ ⟩ profile relaxation. Transport coefficients calculated with test particles subject only to the E×B flow (with 0 parallel velocity) are shown in figure 14 (left). It can be seen that E×B transport is most active during Phase 1, with a typical diffusion coefficient of 1000 m 2 s −1 and advection velocity of 5000 m s −1 , corresponding to a transport timescale of the order of 0.1-0.2 ms. For comparison, figure 14 (center) shows field line transport coefficients, similar to those of figure 12, but here given in m 2 s −1 and m s −1 , having been calculated assuming a parallel propagation at the Alfvén speed. It can be seen that 'SAW transport' is generally much stronger than E×B transport, except during Phase 1 (5.84-6.09 ms), when the two are of the same order, although E×B transport is still somewhat weaker. This supports a contribution of the E×B flow to the relaxation during Phase 1, which could explain why equation (2) under-predicts the η h used in our reference 2D simulation in this phase.
An important remark is that the E×B flow, in addition to possibly contributing to the ⟨j ϕ ⟩ profile relaxation, also plays a major role in mixing particles into the core. In figure 11, it can be seen that the density profile, which is hollow just before the relaxation onset, quickly flattens during Phase 1 of the relaxation. This is essentially due to the E×B flow, while parallel advection along stochastic field lines plays a minor role, as can be seen by comparing the left and right plots in figure 14, which show the transport coefficients for test particles moving respectively at the E×B and plasma parallel velocity.
It would obviously be good to understand the precise mechanisms behind the relaxation. For this, numerical experiments (e.g. parameter scans) would be useful. Unfortunately, these are presently difficult because simulations, in addition to being costly (taking weeks to run on supercomputers), have a strong tendency to get stalled due to numerical issues as the MHD activity becomes violent in conjunction with large density and temperature gradients created by massive material injections. The precise nature of these numerical issues remains to be elucidated, but they appear to be related to negative densities and/or temperatures. Solving this 'numerical bottleneck' will be the object of future efforts.
In section 3.5.1 of [2] and appendix C of [22], it was observed that the experimentally measured I p spike is smaller than the one expected from a complete current profile relaxation at fixed helicity. It was hypothesized that this is due to resistive helicity dissipation taking place during the relaxation. According to our JOREK simulations, the explanation is rather that the relaxation is incomplete, while helicity is well conserved (consistently with the fact that T e remains on the order of 100 eV, and thus that resistive helicity dissipation is small, before I p has reached its peak value).
We mentioned in the introduction that equation (2) appears possibly very valuable because it may allow estimating N t from the measured I p (t) signal, without needing to run costly 3D non-linear MHD simulations. The fact that our results support the validity of this formula is thus encouraging. Our study is however based on a single case and further analyzes on other cases are needed to gain confidence.
In the 3D simulation discussed in this paper, we find N t ≃ 10-25, which is a useful order of magnitude to keep in mind. We however call attention on the fact that N t may be a misleading indicator regarding electron losses. Indeed, it can be seen in figure 15, which was produced from the same field line tracing simulations as figure 13 and shows the number of field lines remaining inside ψ N ⩽ 1 as a function of the number of turns, that it takes substantially more turns for field lines to escape the domain than to distribute evenly across it. In particular, at 5.99 ms (blue curve), it takes 33 turns to lose half the field lines while it takes only about 10 turns for field lines  to distribute evenly (see the blue curves in figure 13). This is related to the fact that field line transport is comparatively slow at the edge, as already observed in [8].
Finally, an open question with critical implications for RE generation is how fast flux surfaces reform after the relaxation. The magnetic field is still globally stochastic at the end of the 3D simulation studied in this paper, which cannot be pushed further because of a numerical issue. This highlights again the importance of solving the 'numerical bottleneck' problem.

Acknowledgment
We thank Allen Boozer for inspiring this work and for many helpful discussions. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200 -EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.
Simulations were run on the Marconi-Fusion supercomputer hosted at CINECA and using HPC resources from GENCI-TGCC (Grant A8-gen11358).

Disclaimer
ITER is the Nuclear Facility INB No. 174. This paper explores physics processes during the plasma operation of the tokamak when disruptions take place; nevertheless the nuclear operator is not constrained by the results of this paper. The views and opinions expressed herein do not necessarily reflect those of the ITER Organization.

Appendix. Boozer's 'prescription' for η h to be used in JOREK
In [2,22], Boozer derives and analyzes a 1D mean field equation for the evolution of the poloidal flux ψ p : .
Here, the toroidal flux ψ t is used as a flux surface label and it is considered that ψ p = ψ p (ψ t , t). I = I(ψ t , t) is the toroidal current inside the flux surface labeled by ψ t at time t, and R ψ is a resistance. In appendix C of [22], assuming an initially parabolic current density profile, the following estimate for the hyper-resistivity coefficient Λ m is given: where V A is the Alfvén speed, Ψ t is the total toroidal flux inside the plasma, N t is the number of toroidal turns needed for a stochastic field line to travel across the plasma, and κ 0 is the vertical elongation (which is assumed not to depend on the flux surface). This estimate comes from the consideration that the current profile should relax on a timescale comparable to the time required for SAWs to travel across the plasma along stochastic field lines, which is τ SAW = 2π R 0 N t /V A , where R 0 is the major radius. The factor 1/144 in equation (A.2), which may seem surprising, finds its origin in the fact that a parabolic profile relaxes under equation (A.1) on a timescale τ parabolic which is 144 times shorter than the 'naive' characteristic timescale of this equation, τ naive = 2κ0 1+κ 2 0 µ 0 R 0 Ψ 2 t /Λ m (note that 2κ0 1+κ 2 0 µ 0 R 0 is the plasma inductance). In [2], this equation and its eigenmodes are analyzed and it is indeed found that already the slowest eigenmode (excluding the mode corresponding to a homogeneous j || /B) has a timescale of only τ naive /56.6. Since we request τ parabolic = τ SAW , not τ naive = τ SAW , a factor 1/144 is needed in equation (A. 2). Note that the current density profile before the relaxation in the JOREK simulation studied here is not exactly parabolic, so equation (A.2) should only be considered as an order of magnitude estimate. Equation (A.1) is obtained from the 3D mean field Ohm's law (equation A.5 in [22]), by performing flux surface averages.
In the present paper, we want to convert Boozer's estimate for Λ m (equation (A.2)) into an estimate for the hyperresistivity η h to be used in JOREK. The conversion is not direct since Λ m appears in a 1D equation (equation (A.1)) while η h appears in a 3D (or 2D under axisymmetry) equation (equation (1)). We thus use λ, which appears in Boozer's 3D equation (equation (A.3)), to make the link.
The equation solved by JOREK (equation (1)) is equivalent to equation (A.3) under a large aspect ratio reduced MHD approximation, with: where B 0 the toroidal field and ρ 0 the reference mass density used in the JOREK normalization [19]. Now, λ and Λ m are related by equation (A.14) in [22]: where the integral is performed over the flux surface labeled by ψ t , θ and ϕ are poloidal and toroidal angles, and J ≡ 1/((∇ψ t × ∇θ) · ∇ϕ) = 1/(2π B · ∇ϕ) is the Jacobian of (ψ t , θ, ϕ) coordinates (see appendix A1 in [22]). In the limit of high aspect ratio and assuming a circular cross-section for simplicity, we have ψ t = π r 2 B 0 and J = R 0 /(2π B 0 ), and equation (A.5) becomes: This leads to: Taking Λ m from equation (A.2), assuming again a circular cross-section and injecting V A = B 0 / √ µ 0 ρ where ρ is the mass density, we obtain: (A.8)