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

E. Nardon, K. Särkimäki, F. J. Artola, S. Sadouni

► To cite this version:

E. Nardon, K. Särkimäki, F. J. Artola, S. Sadouni. On the origin of the plasma current spike during a tokamak disruption and its relation with magnetic stochasticity. Nuclear Fusion, 2023, 63 (5), pp.056011. 10.1088/1741-4326/acc417. cea-04061073

HAL Id: cea-04061073

https://cea.hal.science/cea-04061073

Submitted on 6 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# On the origin of the plasma current spike during a tokamak disruption and its relation with magnetic stochasticity

E. Nardon¹, K. Särkimäki², F.J. Artola³, S. Sadouni¹, the JOREK team‡ and JET contributors§

¹CEA, IRFM, F-13108 Saint-Paul-lez-Durance, France

²Max Planck Institute for Plasma Physics, Boltzmannstr. 2, 85748 Garching b. M., Germany

³ITER Organization, Route de Vinon sur Verdon, 13115 St Paul Lez Durance, France

**Abstract.** 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 ($I_p$) spike [1], is analyzed in order to investigate the origin of the $I_p$ spike and its relation with magnetic stochasticity. The $I_p$ spike is associated to a current density ($j_\phi$) profile relaxation which appears to result from Shear Alfvén Wave (SAW) propagation along stochastic field lines, as proposed by Boozer [2][3], 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 [2][3], which connects the $I_p$ 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_\phi$ 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_\phi$ flattening. During the first phase, strong E×B flows develop that play a key role in mixing impurities into the core.

## 1. 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_\phi$) profile driven by intense magnetohydrodynamic (MHD) activity. It can indeed be shown [5] that magnetic helicity, $H \equiv \int \vec{B} \cdot \vec{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.

‡ See the author list of M. Hoelzl et al. 2021 Nucl. Fusion 61 065001

§ See the author list of ‘Overview of JET results for optimising ITER operation’ by J. Mailloux et al 2022 Nucl. Fusion 62 042026
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_0N_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.|| 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, 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

<footnote>|| It is interesting to note that a comparable $I_p$ 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.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Evolution of the plasma current: global view showing the experiment and 3D JOREK simulation [1] (left), and close-up view also showing the 2D mean field JOREK simulation (see Section 3) (right). Vertical dashed lines in the right figure indicate the different phases of the mean field simulation.</figcaption></figure>

define a reference case in which we set $\eta_h$ so as to match the 3D simulation in terms of $I_p$ and $j_\phi$ profile evolution. We study the effect of various modifications in the $\eta_h$ settings to gain insight. In order to test Boozer's theory, we then compare the $\eta_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 $\tau_{SAW}$. We summarize and discuss our results in Section 4.

## 2. 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_\phi$ profile dynamics and following with a discussion of the mechanisms at play.

### 2.1. $I_p$ and $j_\phi$ 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_\phi$ 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 line-integrated 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, $\langle j_{\phi} \rangle$, 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 $\langle j_{\phi} \rangle$ takes place between 5.84 and 6.00 ms in the region $\psi_N^{1/2} \lesssim 0.7$ ($\psi_N$ is the normalized poloidal magnetic flux, $\psi_N = 0$ corresponding to the magnetic axis and $\psi_N = 1$ to the last closed flux surface in the initial equilibrium). This corresponds to $q \lesssim 2$ initially. The evolution of the $j_{\phi}$ cross-section in Figure 3 (right) is rather spectacular (compare 5.72 and 5.99 ms). The large $\langle j_{\phi} \rangle$ peak visible in the very center at 6.00 ms in Figure 2 corresponds to the central $j_{\phi}$ ‘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 well-defined anymore at this time due to magnetic stochasticity, as can be seen in the Poincaré cross-sections overlayed on the $j_{\phi}$ 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 $\langle j_{\phi} \rangle$ 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 ($\psi_N^{1/2} \lesssim 0.7$). This violent relaxation is associated to a burst in magnetic and kinetic energies in the various toroidal Fourier harmonics $n \ge 1$, as can be seen in Figure 4. Figure 3 (left) shows that the amplitude of magnetic fluctuations normalized to the toroidal field, $\delta 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 $\langle j_{\phi} \rangle$ 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 $\langle j_{\phi} \rangle$ flattening in this region. The fact that fluctuations subsequently grow again is consistent with the reformation of a $\langle j_{\phi} \rangle$ gradient in the core as the (slower) more global $\langle j_{\phi} \rangle$ profile relaxation takes place (compare the magenta and green dashed profiles in Figure 2).

## 2.2. Underlying mechanisms

As explained in the introduction, Boozer proposed that the $\langle j_{\phi} \rangle$ 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 $\simeq 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_{\phi}$ (right) at 6.03 ms. It can also be noticed in Figure 5 that there is
<figure><img src="image_3.png" /><figcaption>Figure 2. Flux-surface-averaged current density profile in the 3D (plain) and 2D (dashed) simulations. The blue profile at 5.84 ms is common to both simulations.</figcaption></figure>

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

### 3. 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 $\psi$ reads (in JOREK normalized units [19]):

$$ \frac{\partial \psi}{\partial t} = R[\psi, u] + \eta j - R^2 \nabla \cdot \left( \frac{\eta_h}{R} \nabla j \right), \quad (1) $$

where $j \equiv \Delta^* \psi \equiv R^2 \nabla \cdot (\frac{1}{R^2} \nabla \psi)$. The first term on the right-hand side is an advection term ([, .] is a Poisson bracket and $u$ is the flow potential), the second one is
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Poloidal cross-sections showing: left column: the magnitude of the non-axisymmetric part of the magnetic field $\delta B$ normalized to the toroidal magnetic field $B_t$; middle column: the $E \times B$ flow potential (colors) and streamlines (black lines); right column: the toroidal current density $j_\phi$ (colors) as well as Poincaré cross-sections (white dots) (right) at different times in the 3D simulation. Note that cross-sections showing the evolution of the axisymmetric part of $j_\phi$, which make the flattening more</figcaption></figure>

clear.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Normalized magnetic and kinetic energies in the various toroidal Fourier harmonics in the 3D simulation.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Normalized magnetic and kinetic energies versus toroidal Fourier mode number at 5.98 ms in the 3D simulation. The normalization factor is the same for magnetic and kinetic energies.</figcaption></figure>

the resistive term and the third is the hyper-resistive term.

We start these 2D simulations at 5.84 ms, i.e. just before $\langle j_\phi \rangle$ 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 $\eta$ and hyper-resistivity $\eta_h$, which allow reproducing the evolution of the $\langle j_\phi \rangle$ 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 $\eta$ and $\eta_h$, which helps gain intuition. Finally, in Section 3.3, we calculate the value of $\eta_h$ predicted by Boozer's formula and compare it to our reference case.

### 3.1. Reference settings and comparison to 3D simulation

We found that in addition to prescribing $\eta_h$, it is also necessary to prescribe $\eta$ in these 2D simulations since calculating $\eta$ 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 pre-disruption $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 $\eta = 10^{-7}$ in the bulk during this phase, i.e. between 5.84 and 6.59 ms (we give $\eta$ in JOREK units here and below). This value corresponds to $4.2 \times 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 $\eta = 2 \times 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:

$$ \eta(\psi_N) = 10^{-7} + (2 \times 10^{-5} - 10^{-7}) \times 0.5 \times (1 + \tanh((\psi_N - 0.95)/0.05)) $$

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 $\eta$ profile with a value of $2 \times 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 $\eta(\psi_N)$ profile in the different phases of the simulation is shown in Figure 6.

Regarding $\eta_h$, based on observations from Section 2, we define 3 phases for which we use 3 different profiles. During Phase 1, between 5.84 and 6.09 ms, we use a large $\eta_h$ in the core ($\psi_N \le 0.5$): $\eta_h = 2.4 \times 10^{-6}$ in JOREK units, in order to reproduce the fast $\langle j_\phi \rangle$ flattening there, while we use a 6 times lower value, $\eta_h = 4 \times 10^{-7}$, in the outer part of the bulk plasma, and we cut-off $\eta_h$ in the SOL. The complete hyper-resistivity profile in Phase 1 is given by

$$ \eta_{h,1}(\psi_N) = 2 \times 10^{-6} \times 0.5 \times (1 - \tanh((\psi_N - 0.5)/0.05)) + 4 \times 10^{-7} \times 0.5 \times (1 - \tanh((\psi_N - 0.95)/0.05)) $$

and is shown in Figure 6 (left). In Phase 2, between 6.09 and 6.59 ms, $\eta_h$ is reduced to $8 \times 10^{-7}$ in the core and maintained at $4 \times 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

$$ \eta_{h,2}(\psi_N) = 4 \times 10^{-7} \times 0.5 \times (1 - \tanh((\psi_N - 0.5)/0.05)) + 4 \times 10^{-7} \times 0.5 \times (1 - \tanh((\psi_N - 0.95)/0.05)) $$
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Profiles of the normalized resistivity $\eta$ and hyper-resistivity $\eta_h$ during the 3 successive phases of the ‘best match’ 2D mean field model simulation. Note that the $\eta_h$ curve is not visible in the rightmost plot because $\eta_h = 0$ during Phase 3.</figcaption></figure>

and is shown in Figure 6 (middle). Finally, in Phase 3, after 6.59 ms, we simply set $\eta_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.

<table><caption>**Table 1.** Settings of the ‘best match’ 2D mean field model simulation (see also Figure 6 and equations given in the main text).</caption><thead><tr><th>Phase</th><th>1</th><th>2</th><th>3</th></tr></thead><tbody><tr><td>Time (ms)</td><td>5.84-6.09</td><td>6.09-6.59</td><td>6.59-end</td></tr><tr><td>η(ψ<sub>N</sub> ≲ 0.9)</td><td>10<sup>-7</sup></td><td>10<sup>-7</sup></td><td>2 × 10<sup>-5</sup></td></tr><tr><td>η(ψ<sub>N</sub> ≳ 1.0)</td><td>2 × 10<sup>-5</sup></td><td>2 × 10<sup>-5</sup></td><td>2 × 10<sup>-5</sup></td></tr><tr><td>η<sub>h</sub>(ψ<sub>N</sub> ≲ 0.45)</td><td>2 × 10<sup>-6</sup></td><td>8 × 10<sup>-7</sup></td><td>0</td></tr><tr><td>η<sub>h</sub>(0.55 ≲ ψ<sub>N</sub> ≲ 0.9)</td><td>4 × 10<sup>-7</sup></td><td>4 × 10<sup>-7</sup></td><td>0</td></tr><tr><td>η<sub>h</sub>(ψ<sub>N</sub> ≳ 1.0)</td><td>0</td><td>0</td><td>0</td></tr></tbody></table>

Figure 1 (right) shows that the 2D simulation reproduces well the evolution of $I_p$ from the 3D simulation. The evolution of the $\langle j_\phi \rangle$ 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 $\eta_h$ (e.g. with an inward-expanding large $\eta_h$ region), but we prefer keeping the setting simple.

## 3.2. Sensitivity studies

We now test different variations of the $\eta_h$ settings to see their influence on 2D simulation results. We begin with investigating the effect of $\eta_h$ in the outer part of the bulk plasma ($0.5 \le \psi_N \le 0.95$) during Phase 1. First of all, if we set $\eta_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 $\langle j_\phi \rangle$ still flattens fast in the core, a negative skin current appears at the edge of the flattened
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Effect of changing $\eta_h$ in the outer part of the bulk plasma ($0.5 \le \psi_N \le 0.95$) during Phase 1 (5.84-6.09 ms).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Effect of keeping $\eta_h$ in Phase 2 (6.09 - 6.59 ms) the same as in Phase 1 (5.84 - 6.09 ms).</figcaption></figure>

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 $\eta_h$ to a large value in the outer part of the bulk plasma, e.g. the same value as in the core ($2.4 \times 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 $\langle j_\phi \rangle$ profile flattening is too global.

Secondly, we test the effect of keeping $\eta_h$ during Phase 2 (6.09-6.59 ms) the same as in Phase 1 (5.84-6.09 ms), with a large $\eta_h$ ($2.4 \times 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 $\langle j_\phi \rangle$ profile evolution, leading to too much flattening in the core during Phase 2.

A third test is to maintain $\eta_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 $\langle j_\phi \rangle$ 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.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Effect of keeping $\eta_h$ in Phase 3 (after 6.59 ms) the same as in Phase 2 (6.09-6.59 ms).</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Effect of changing the onset time of Phase 3.</figcaption></figure>

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 $\eta$, 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 $\eta$ and $\eta_h$ profiles at the edge during the first two phases. This radius was chosen to be $\psi_N = 0.95$ for both $\eta$ and $\eta_h$ in the reference case, but changing it (whether for $\eta$ or $\eta_h$) by $\pm 0.05$ does not modify the results significantly. In particular, the precise degree of overlap between the large $\eta$ region (the SOL) and the region where $\eta_h$ is applied does not have much influence on the results. This is probably due to the fact that no sharp $\langle j_\phi \rangle$ gradients are produced during the process.

### 3.3. Testing Boozer’s formula for $\eta_h$

In [2][22], Boozer discusses the relaxation of the $\langle j_\phi \rangle$ profile during a fast MHD relaxation, which he models with an equation equivalent to Equation 1 above. As mentioned in Section 1, references [2][22] also propose an estimate for the level of hyper-
<figure><img src="image_12.png" /><figcaption>Figure 11. Midplane cut of the axisymmetric part of the mass density ρ.</figcaption></figure>

resistivity based on the assumption that the timescale of the relaxation should be of the order of $\tau_{SAW} = 2\pi R_0 N_t / V_A$. Under certain approximations, this estimate can be translated into the following expression for the hyper-resistivity $\eta_h$ to be used in JOREK 2D mean field model simulations, as detailed in Appendix A:

$$ \eta_h = \sqrt{\frac{\rho_0}{\rho} \frac{B_0 a^4}{4608 \pi R_0^2 N_t}}, \qquad (2) $$

with $\rho$ the mass density, $\rho_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, $\rho_0 = 7.1 \times 10^{-8}$ kg/m$^3$. As can be seen in Figure 11, $\rho$ 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 $\rho \simeq 4 - 5\rho_0$. For simplicity, we will consider that $\sqrt{\rho_0/\rho} = 1/2$ in Equation 2.

We want to test whether Equation 2 corresponds to the level of $\eta_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) (left) and advection coefficient (in m/m) $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é cross-section 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
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Field line diffusion coefficient $D_{FL}$ and drift (advection) velocity $A_{FL}$ as a function of time and minor radius in the 3D simulation.</figcaption></figure>

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 modelled 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} \approx 0.5 - 1 \times 10^{-3} \text{ m}^2/\text{m}$ and $A_{FL} \approx 1 - 2 \times 10^{-3} \text{ m/m}$. Since $a \approx 1 \text{ m}$, this means that the Péclet number [23] $aA_{FL}/D_{FL} \approx 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 behaviour.

The typical distance along a field line required to travel across the minor radius is $a^2/D_{FL} \approx 1000 - 2000 \text{ m}$ (resp. $a/A_{FL} \approx 500 - 1000 \text{ m}$) when only diffusion (resp. advection) is taken into account. Combining diffusion and advection, we estimate the required distance to be $\approx 500 \text{ m}$, which corresponds to $N_t \approx 25$ since $R_0 \approx 3 \text{ 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. $\psi_N^{1/2}$ intervals) and calculated the approximate field line density in each bin by dividing the number of field lines in the bin by $\psi_N^{1/2}$, 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 ($\psi_N^{1/2}$) of remaining
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Left: Radial density of test field lines as a function of the number of toroidal turns for a set of 200 lines initialized at $R = 3.0$ m, $Z = 0.2$ m and equally spaced toroidal positions at $t = 5.99$ ms. Center and right: Statistics on such simulations for various initial $R$ (plain = 3.0 m, dashed = 3.3 m, dash-dotted = 3.6 m) and $t$ (blue = 5.99 ms, red = 6.34 ms), showing the average (center) and variance (right) of the radial position ($\psi_N^{1/2}$) of remaining field lines.</figcaption></figure>

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 $\simeq 0.25$, which corresponds to a spreading over the whole volume. This analysis thus suggests that $N_t \simeq 10 - 20$, which is of the same order, although somewhat smaller, as our above estimate $N_t \simeq 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 $\eta_h \simeq 6 \times 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 \times 10^{-7}$ used in the outer part of the bulk plasma in Phases 1 and 2 and the $8 \times 10^{-7}$ used in the core in Phase 2, while it is 4 times smaller than the $2.4 \times 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 $\sqrt{\rho_0/\rho}$ 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.

### 3.4. Comparison of the relaxation timescale with $\tau_{SAW}$

A more direct test of Boozer's theory is to compare the observed relaxation timescale, $\tau_{relax}$, with $\tau_{SAW}$. Judging from the timescale of the $\langle j_\phi \rangle$ profile evolution, we may estimate that $\tau_{relax} \simeq 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)/\tau_{relax})$, we find a larger value, $\tau_{relax} \simeq 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 $\tau_{relax} \simeq 100~\mu\text{s}$ is the better estimate. Taking $\rho = 4\rho_0 \simeq 3 \times 10^{-7} \text{ kg/m}^3$, we have $V_A \simeq 5 \times 10^6 \text{ m/s}$, leading to $\tau_{SAW} \simeq 75~\mu\text{s}$ for $N_t = 20$, which is indeed close to $\tau_{relax}$.

## 4. Summary and discussion

In summary, it appears plausible that the $\langle j_\phi \rangle$ 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 $\tau_{SAW}$ and that the level of hyper-resistivity 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 $\langle j_\phi \rangle$ 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 $\langle j_\phi \rangle$ 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 \text{ m}^2\text{/s}$ and advection velocity of $5000 \text{ m/s}$, 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 $\text{m}^2\text{/s}$ and m/s, 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 $\eta_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 $\langle j_\phi \rangle$ 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
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Radial drift (advection) (top) and diffusion (bottom) coefficients associated to E×B flows (left), SAWs (center) and parallel advection (right). Note that E×B and SAWs share the same colour scale.</figcaption></figure>

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 \approx 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
<figure><img src="image_16.png" /><figcaption>**Figure 15.** Number of field lines remaining inside $\psi_N \le 1$ in the same test field line simulations as in Figure 13</figcaption></figure>

15, which was produced from the same field line tracing simulations as Figure 13 and shows the number of field lines remaining inside $\psi_N \le 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.

**5. Acknowledgements**

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 A. Boozer’s ‘prescription’ for ηh to be used in JOREK

In [22] and [2], Boozer derives and analyzes a 1D mean field equation for the evolution of the poloidal flux ψp :

$$
\frac{\partial \psi_p}{\partial t} = R_{\psi} \frac{\partial I}{\partial \psi_t} - \frac{\partial}{\partial \psi_t} (\psi_t \Lambda_m \frac{\partial^2 I}{\partial \psi_t^2}). \quad (A.1)
$$

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 labelled 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:

$$
\Lambda_m = \frac{1}{144} \frac{2\kappa_0 \mu_0 V_A \Psi_t^2}{1 + \kappa_0^2 4\pi N_t}, \quad (A.2)
$$

where $V_A$ is the Alfvén speed, $\Psi_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 $\kappa_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 $\tau_{SAW} = 2\pi 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 $\tau_{parabolic}$ which is 144 times shorter than the ‘naive’ characteristic timescale of this equation, $\tau_{naive} = \frac{2\kappa_0}{1+\kappa_0^2}\mu_0 R_0 \Psi_t^2 / \Lambda_m$ (note that $\frac{2\kappa_0}{1+\kappa_0^2}\mu_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_\parallel/B$) has a timescale of only $\tau_{naive}/56.6$. Since we request $\tau_{parabolic} = \tau_{SAW}$, not $\tau_{naive} = \tau_{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]),

$$
\mathbf{E} \cdot \mathbf{B} = \eta \mathbf{j} \cdot \mathbf{B} - \nabla \cdot \left( \lambda \nabla \frac{j_{\|}}{B} \right), \quad (A.3)
$$

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 hyper-resistivity η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:

$$
\eta_h = \sqrt{\frac{\rho_0}{\mu_0 R_0 B_0^2}}, \qquad (A.4)
$$
where $B_0$ the toroidal field and $\rho_0$ the reference mass density used in the JOREK normalization [19].

Now, $\lambda$ and $\Lambda_m$ are related by Equation A.14 in [22]:

$$ \Lambda_m(\psi_t, t) = \frac{\oint \lambda(\nabla\psi_t)^2 J d\theta d\phi}{\psi_t}, \qquad (A.5) $$

where the integral is performed over the flux surface labelled by $\psi_t$, $\theta$ and $\phi$ are poloidal and toroidal angles, and $J \equiv 1/((\nabla\psi_t \times \nabla\theta) \cdot \nabla\phi) = 1/(2\pi\mathbf{B} \cdot \nabla\phi)$ is the Jacobian of $(\psi_t, \theta, \phi)$ coordinates (see Appendix A1 in [22]). In the limit of high aspect ratio and assuming a circular cross-section for simplicity, we have $\psi_t = \pi r^2 B_0$ and $J = R_0/(2\pi B_0)$, and Equation A.5 becomes:

$$ \Lambda_m = 8\pi^2 R_0 \lambda. \qquad (A.6) $$

This leads to:

$$ \eta_h = \sqrt{\frac{\rho_0}{\mu_0}} \frac{\Lambda_m}{8\pi^2 R_0^2 B_0^2}. \qquad (A.7) $$

Taking $\Lambda_m$ from Equation A.2, assuming again a circular cross-section and injecting $V_A = B_0/\sqrt{\mu_0\rho}$ where $\rho$ is the mass density, we obtain:

$$ \eta_h = \sqrt{\frac{\rho_0}{\rho}} \frac{B_0 a^4}{4608\pi R_0^2 N_t}. \qquad (A.8) $$

[1] E. Nardon et al., ‘Thermal quench and current profile relaxation dynamics in massive material-injection-triggered tokamak disruptions’, Plasma Phys. Control. Fusion 63 115006 (2021)

[2] A. Boozer, ‘Magnetic surface loss and electron runaway’, Plasma Phys. Control. Fusion 61 024002 (2019)

[3] A. Boozer, ‘Flattening of the tokamak current profile by a fast magnetic reconnection with implications for the solar corona’, Phys. Plasmas 27 102305 (2020)

[4] T. Hender et al., ‘Progress in ITER Physics Basis, Chapter 3: MHD stability, operational limits and disruptions’, Nucl. Fusion 47 S128-202 (2007)

[5] D. Biskamp, ‘Nonlinear magnetohydrodynamics’, Cambridge University Press (1997)

[6] I. Sveningsson et al., ‘Hot-Tail Runaway Seed Landscape during the Thermal Quench in Tokamaks’, Phys. Rev. Lett. 127 035001 (2021)

[7] P. Svensson et al., ‘Effects of magnetic perturbations and radiation on the runaway avalanche’, J. Plasma Phys. 87 905870207 (2021)

[8] C. Sommariva et al., ‘Test particles dynamics in the JOREK 3D non-linear MHD code and application to electron transport in a disruption simulation’, Nucl. Fusion 58 016043 (2017)

[9] C. Sommariva et al., ‘Electron acceleration in a JET disruption simulation’, Nucl. Fusion 58 106022 (2018)

[10] K. Särkimäki et al., ‘Confinement of passing and trapped runaway electrons in the simulation of an ITER current quench’, to be submitted

[11] V. Izzo and P. Parks, ‘Comment on Plasma current spikes due to internal reconnection during tokamak disruptions’, Nucl. Fusion 50 058001 (2010)

[12] E. Nardon et al., ‘Progress in understanding disruptions triggered by massive gas injection via 3D non-linear MHD modelling with JOREK’, Plasma Phys. Control. Fusion 59 014006 (2017)

[13] A. Boozer, ‘Ohm’s law for mean magnetic fields’, J. Plasma Phys. 35 133 (1986)

[14] D. Biskamp, ‘Magnetic Reconnection in Plasmas’, Cambridge University Press (2000)
[15] A. Aydemir, ‘Magnetohydrodynamic modes driven by anomalous electron viscosity and their role in fast sawtooth crashes’, Phys. Fluids B: Plasma Phys. **2** 2135 (1990)

[16] P. Kaw et al., ‘Tearing Modes in a Plasma with Magnetic Braiding’, Phys. Rev. Lett. **43** 1398 (1979)

[17] I. Pusztai et al., J. Plasma Phys. **88** 905880409 (2022)

[18] D. Hu et al., ‘Radiation asymmetry and MHD destabilization during the thermal quench after impurity shattered pellet injection’, Nucl. Fusion **61** 026015 (2021)

[19] M. Hoelzl et al., ‘The JOREK non-linear extended MHD code and applications to large-scale instabilities and their control in magnetically confined fusion plasmas’, Nucl. Fusion **61** 065001 (2021)

[20] A. Lichtenberg and M. Lieberman, ‘Regular and stochastic motion’. Springer (2013)

[21] T. Evans et al., ‘Homoclinic tangles, bifurcations and edge stochasticity in diverted tokamaks’, Contrib. Plasma Phys. **44** 235 (2004)

[22] A. Boozer, ‘Pivotal issues on relativistic electrons in ITER’, Nucl. Fusion **58** 036006 (2018)

[23] K. Särkimäki et al., ‘An advection-diffusion model for cross-field runaway electron transport in perturbed magnetic fields’, Plasma Phys. Control. Fusion **58** 125017 (2016)

[24] A. Boozer and G. Kuo-Petravic, ‘Monte Carlo evaluation of transport coefficients’, Phys. Fluids **5** 851 (1981)