Magnetohydrodynamic simulations of runaway electron beam termination in JET

Disruptions in tokamaks occur as a result of magnetic field stochastization caused by large scale magnetohydrodynamic (MHD) instabilities [1, 2]. This leads to a loss of most of the plasma thermal energy on a millisecond timescale, cooling it to a temperature ${\sim}{10}\ \mathrm{eV}$ [3]. While initial thermal losses are dominated by transport along stochastic field lines, the post thermal quench temperature evolution is determined by the relative magnitudes of Ohmic heating and impurity radiation [4]. The temperature drop increases the electrical resistivity of the plasma by several orders of magnitude, causing the plasma current to decay on a resistive timescale, referred to as the current quench (CQ). Increased resistivity also causes a large toroidal electric field, which can accelerate suprathermal electrons to relativistic velocities and energies of a few tens of MeV. The post thermal quench temperature via the impurity density and species play a crucial role here. Such runaway electrons (REs) could eventually form a beam carrying a large fraction of the pre-disruption current in fusion grade devices such as ITER [5–7]. Uncontrolled loss of REs can cause localized wall damage [8–10] and significant machine downtime. Hence the necessity to devise robust RE avoidance/mitigation strategies, motivating an improved understanding of plasma dynamics during the lifetime of the RE beam. The interplay between REs and plasma instabilities becomes important in this respect.

Several strategies have been proposed over the years for RE avoidance or mitigation [5, 11–13]. A potentially promising one has been demonstrated in a series of recent experiments at JET [12]. It has been observed that second injection of deuterium into an established post-disruption RE beam (created via high-Z impurity injection earlier in the discharge) can lead to benign terminations, i.e. without any wall damage. This can be attributed to a significant flush-out of high-Z impurities in the stationary beam phase followed by a large MHD instability causing RE losses. The lost REs are supposedly not replenished via reacceleration, mainly due to the deficiency of high-Z impurities [12]. In order to be able to potentially replicate similar benign terminations in ITER, it would be important to understand the MHD activity causing the RE losses. This is the focus of the present work.

In this paper, we present 3D non-linear MHD simulations of one such experimental shot at JET leading to a benign RE beam termination. Simulations are performed using the JOREK code [14–16]. JOREK is an extended 3D MHD code based on fully-implicit time stepping, and a spatial discretization with two-dimensional (2D) isoparametric Bezier Galerkin finite-elements in the poloidal plane and a Fourier decomposition in the toroidal direction. REs are treated via a fluid model self-consistently coupled to MHD within JOREK [17] (a kinetic model [18, 19] too is available in JOREK, but without back coupling to the plasma so far). The aim here is to shed light on the conditions causing the large instability and the ensuing non-linear dynamics until the final RE loss.

The paper is organized as follows. Section 2 summarizes the essential details of the experiment being simulated (refer to [12, 20] for further details), followed by a description of the MHD model and equilibrium in section 3. Simulation results and analysis is provided in section 4, followed by concluding remarks in section 5.

The experimental shot of interest in this work is #95135 at JET, in which a disruption is triggered by the injection of massive amounts of argon impurities in the form of shattered pellets. This causes a thermal quench and a prompt conversion to a runaway current beam over the CQ (see figure 1). During the plateau phase of the RE beam with a plasma current $I_p\approx {0.5}\ \mathrm{MA}$, deuterium shattered pellets are injected into the beam (second injection). After the second shattered pellet injection (SPI), it has been observed that most of the argon impurity species are flushed-out from the plasma due to recombination [12, 21]. The flush out of argon impurities causes a significant drop in the hard x-ray (HXR) signal as shown in figure 1(c), as Bremsstrahlung from REs is a strong function of the high-Z impurity content. The drop in HXR causes an associated drop in the rate of photoneutrons released by the interaction of HXR with Beryllium wall material as shown in figure 1(b). Expulsion of significant amount of high-Z impurities also reduces the effective resistivity and leads to a gradual increase in the plasma current [12, 21]. After about ${240}\ \mathrm{ms}$ from the second SPI trigger, when the plasma current reaches ${\approx}{0.75}\ \mathrm{MA}$, a fast and near-complete loss of REs from the plasma is observed. This is followed by a benign termination of the discharge in a few tens of milliseconds, without any localized damage on the first wall surface. The sudden loss of REs causes a large spike in the neutron count and HXR signals, as can be seen in figure 1. After the fast loss of REs, further on, there is no indication of the presence of any REs as seen from the infrared synchrotron images. This indicates that as impurity radiation is reduced significantly after the flush-out, the resulting lower toroidal electric field renders avalanche regrowth of REs from remnants ineffective.

Zoom In Zoom Out Reset image size

We use an ansatz based reduced MHD model in JOREK for the simulations, with a single fluid representation of the background plasma consisting of thermal ions and electrons. For simplicity, since the thermal pressure is too low to affect the dynamics, the background plasma density ρ and temperature T are assumed to be time invariant. The magnetic field $\boldsymbol{B}$ and electric field $\boldsymbol{E}$ are, respectively, represented as

$$\begin{align} \boldsymbol{B} & = R^{-1} \left[ \nabla \psi \times \boldsymbol{e_{\phi}} + F_0 \boldsymbol{e_{\phi}} \right], \end{align} \tag{ 1 }$$

$$\begin{align} \boldsymbol{E} & = -F_0 \nabla u - R^{-1} \partial_t \psi \end{align} \tag{ 2 }$$

where, R is the major radial coordinate, ψ is the poloidal magnetic flux, F₀ u the electric potential, $\boldsymbol{e_\phi}$ a unit vector in the toroidal direction and F₀ a constant. The ion fluid velocity $\boldsymbol{v}$ consists of the leading order $\boldsymbol{E}\times\boldsymbol{B}$ drift and is given by

$$\begin{align} \boldsymbol{v} = \frac{\boldsymbol{E}\times\boldsymbol{B}}{B^2} \approx -R \nabla u \times \boldsymbol{e_\phi} \end{align} \tag{ 3 }$$

where B is the magnitude of the magnetic field. REs are represented as a fluid species with number density n_r , that interacts with the background plasma through a current coupling. The total current density $\boldsymbol{j}$ can be seen as decomposed into a thermal current density $\boldsymbol{j_\mathrm{th}}$ and an RE current density $\boldsymbol{j_r}$ as

$$\begin{align} \boldsymbol{j} = \boldsymbol{j_\mathrm{th}} + \boldsymbol{j_r}, \quad \boldsymbol{j_r} = -e \ c \ n_r \boldsymbol{\hat{b}}, \end{align} \tag{ 4 }$$

where e and c are electron charge and speed of light, respectively, and $\boldsymbol{\hat{b}}$ is the unit vector along the magnetic field. Transport of REs is described here by a large parallel diffusion $D_{r,||}$ as an ad-hoc to the computationally more demanding parallel advection, which is also an option in the JOREK RE fluid model. As described in [17], use of parallel diffusion that ensures transport timescale equivalent to advection at speed of light, can adequately capture the stochastic loss of REs along the field lines with some limitations (these aspects are discussed later in this paper). Furthermore, REs are assumed to have a constant parallel momentum, and RE sources are neglected.

For our simulations, we choose a starting point that represents the plasma state a few milliseconds before the final crash. As indicated earlier, most of the argon impurities are experimentally observed to be lost by this time. Therefore we do not model impurities in the present simulations. Additionally, presence of neutrals in the plasma are not considered. The equations governing the coevolution of MHD and REs are given in normalized form by

$$\begin{align} \frac{1}{R^2} \partial_t \psi &= \frac{\eta}{R^2} \left(j -cn_{r}\frac{F_0}{BR}\right) - \frac{F_0}{R^2} \partial_\phi u \nonumber \\ &\quad- \frac{1}{R} \left[u,\psi\right] - \frac{\eta_h}{R}\nabla^2\left(\frac{j-cn_{r}F_0/\left(BR\right)}{R}\right) \end{align} \tag{ 5a }$$

$$\begin{align} \nabla\cdot\left[\rho R^2\nabla_\perp \frac{\partial u}{\partial t}\right] &= \frac{1}{2{R}}\left[R^2 |\nabla_\perp u|^2,R^2 p\right] \nonumber \\ &\quad+ \frac{1}{R}\left[R^4 \rho\omega,u\right] + \frac{1}{R} \left[\psi,j\right]\nonumber \\ &\quad- \frac{F_0}{R^2}\partial_\phi j -\frac{1}{R}\left[R^2,\rho T\right] \nonumber \\ &\quad+ R\mu_\perp\left(T\right)\nabla^2\omega + R\mu_{\perp,h}\left(T\right)\nabla^4\omega \end{align} \tag{ 5b }$$

$$\begin{align} j &= \Delta^{*}\psi \end{align} \tag{ 5c }$$

$$\begin{align} \omega &= \nabla\cdot\nabla_{\perp}u \end{align} \tag{ 5d }$$

$$\begin{align} \left[\partial_t \rho \hspace{4mm} \partial_t T\right] &= \left[0 \hspace{4mm} 0 \right] \end{align} \tag{ 5e }$$

$$\begin{align} \partial_t n_r &= \left(i_D -1\right) \left[-\frac{c}{B R} \left( \left[n_r,\psi\right] + \frac{F_0}{R} \partial_\phi n_r \right)\right] \nonumber \\ &\quad+ \nabla\cdot\left(i_D D_{r,\parallel}\nabla_\parallel n_{r} + D_{r, \perp}\nabla_\perp n_{r}\right). \end{align} \tag{ 5f }$$

We retain the same names for the normalized variables for simplicity. Details of normalization are provided in appendix A. In the above, η is the electrical resistivity, η_h the hyperresistivity, j is the toroidal current density, ω the toroidal vorticity and p the plasma pressure. Perpendicular viscosity and hyperviscosity are given by $\mu_\perp$ and $\mu_{\perp,h}$, respectively, while i_D is a boolean value to switch between parallel advection and parallel diffusion for RE density. A small perpendicular diffusion is used for RE density for numerical stability, via the diffusivity $D_{r,\perp}$. For η_h and $\mu_{\perp,h}$, small enough values are used such that the dynamics remain unaffected.

As mentioned earlier, we start with an equilibrium corresponding to several milliseconds before the crash, with $Ip = {0.747}\ \mathrm{MA}$ and q₀ = 5.3. It is a challenge to obtain direct information regarding the central safety factor q₀ and the shape of the q-profile from the experiment. This however, can be partly resolved through indirect evidence related to MHD mode structure and synchrotron emission before the crash. This is described in the following.

Excellent data on synchrotron radiation acquired during the SPI RE experiments at JET, especially using the IR camera (KLDT-E5WC, λ = 3–3.5 $\mu\mathrm{m}$) allows to discriminate some of the RE beam properties. Despite the fact that synchrotron radiation pattern and intensity is a function of multiple variables—local magnetic field direction and intensity, RE energy, pitch angle and local RE density on one side and the camera properties on the other—it is possible to infer various qualitative observations from the image data and forward synchrotron pattern modeling using, e.g. SOFT code [22]. Simulations were performed using SOFT Green's function tool to reconstruct the synchrotron images that would result from different RE density profiles (including JOREK equilibrium profile) and different regions of phase space allowed. Figure 2 shows the original camera image with subtracted background and the reconstructed image from JOREK equilibrium. Clearly, a qualitatively good match is observed with the JOREK profile. The results here corresponds to a moderate pitch angle range θ = 0.1–0.3 radians and lower energies (1–15 MeV). However, it is observed that the conclusions hold true also in the case when contributions from higher energy or higher/lower pitch angle populations are considered. SOFT simulations were repeated (for the set of RE density profiles described in appendix B) with exclusive contribution from each of the following four regions of the phase space, namely

• (a)  
  Low pitch angle (0.1–0.3 rad), low energy 0.5–15 MeV
• (b)  
  Low pitch angle (0.1–0.3 rad), high energy only 20–25 MeV
• (c)  
  High pitch angle (0.3–0.6 rad), low energy 0.5–15 MeV
• (d)  
  High pitch angle (0.3–0.6 rad), high energy only 20–25 MeV.

Zoom In Zoom Out Reset image size

In all of these cases, hollow profiles provide synchrotron images that were very distinct from those obtained from any peaked profile in the IR camera spectral sensitivity band—crescent-like shape (as the one observed in experiment) extended to the top and bottom parts of the circular cross-section could be achieved only by a hollow profile. Images from the peaked profile were far from experimental observation throughout the studied phase space. Furthermore, comparison with other RE density profiles considered excludes the probability of non-hollow profiles (see appendix B). Therefore, as summarized here and also discussed in [12], it can be concluded that only a hollow current density profile allows to recover the experimental radiation distribution.

Furthermore, m = 4 magnetic islands visible in the synchrotron images before the crash (not shown here) provide an additional constraint on the location of the q = 4 surface, that can only be satisfied by a hollow profile. In fact, the equilibria with hollow profiles feature two q = 4 surfaces and the location of the inner one is able to explain the experimental structures. JOREK MHD test simulations (not shown here for brevity) with a monotonous current profile for the equilibrium obtained from EFIT did not reproduce the correct dynamics as in the experiment. For example, these simulations showed a dominant (m, n) = (2, 1) mode as compared to the dominant (m, n) = (4, 1) mode observed in the experiment, and also did not result in any stochastization of the plasma. On the other hand, simulations with different hollow profiles recover these dynamics in a very robust way as shown later i.e. the results are not strongly sensitive to the value of q₀. To summarize, the SOFT simulations, observation of m = 4 islands and JOREK test simulations provide sufficient evidence for the existence of a hollow current density profile before the crash. Hence only hollow profiles are considered henceforth.

Further constraints are necessary in order to completely define the current/q-profile. This is provided by information on the MHD modes as observed from the infrared synchrotron data over a time span of a few tens of milliseconds starting before the crash. This is shown in figure 3 which shows the approximate minor-radial location of the m = 6, m = 5 and m = 4 islands during various time points until the crash, where m is the poloidal mode number. The islands shown here correspond to the n = 1 toroidal mode, so that m = q, the safety factor. It has been observed that the m = 6 islands move radially outward and eventually leave the plasma in ${\sim}{50}\ \mathrm{ms}$, as shown in figure 3. Following this, the m = 5 islands are visible which again move radially outward. They are invisible after around ${48.58}\ \mathrm{s}$, but are still within the last closed flux surface. Finally, before the crash (indicated by the red vertical dotted line), m = 4 islands are visible at much lower minor-radial locations. This provides an approximate location for the inner q = 4 surface for the equilibrium being sought.

Zoom In Zoom Out Reset image size

A scan of q profiles with the q = 4 surface in the vicinity of $r = {0.2}\ \mathrm{m}$ have been performed to determine the non-linearly dominant toroidal mode in the JOREK simulations (shown in figure 4). It can be seen that with an increase in q₀, beyond a threshold, there is a clear transition in the most unstable MHD mode, from an n = 2 dominant behavior to an n = 1 dominant behavior. In line with the aforementioned experimental observation of n = 1 dominant mode before the crash, one of the q-profiles (green colored lines) would satisfy now all the known constraints. Furthermore, it is seen from the sensitivity test that the final results of interest are not very sensitive to the precise profile chosen among the plausible ones. This enables us to choose the dotted profile (in figure 4) as representative of the equilibrium before the crash.

Zoom In Zoom Out Reset image size

The computational domain used is shown in figure 5(a) colored by the normalized RE density at equilibrium, along with the flux surfaces. The boundary of the domain closely approximates the first wall surface. The runaway beam is nearly circular and approximately ${0.7}\ \mathrm{m}$ in diameter. The plasma current $I_p = {0.747}\ \mathrm{MA}$ and the on-axis toroidal magnetic field $B_0 = {3.49}\ \mathrm{T}$. Being a cold plasma, the equilibrium is practically pressure-less. The corresponding equilibrium safety factor and current density profiles are summarized in figure 5(b). The background plasma density $n_e = {1 \times 10^{19}}\ {\mathrm{m}^{-3}}$. The resistivity and viscosity are spatially uniform, with values $\eta = {3.3 \times 10^{-5}}\ {\Omega\mathrm{m}}$ (which corresponds to a Spitzer value at $T = {10}\ \mathrm{eV}$) and $\mu = {1.6 \times 10^{-6}}\ \mathrm{kg\ m^{-1} s^{-1}}$. All the current is assumed to be carried by REs at the initial state. Runaway parallel diffusion is chosen to be $D_{r,||} = {10^{10}}\ \mathrm{m^{2}\ s^{-1}}$ which is ${\sim}L_{||} c$, where $L_{||}\sim {100}\ \mathrm{m}$ is an appropriate length scale for stochastic parallel transport or connection length. Simulations were performed including toroidal modes n = 0, ..., 8 on a poloidal grid $n_R\times n_Z = 161\times64$ with radially local clustering to resolve resonant surface regions. A dedicated grid sensitivity study showed no noticeable change in linear growth rates with further increase in grid resolution.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 6(a) shows the evolution of the energy in the poloidal magnetic field (in the toroidal modes n = 1, ..., 8) with time. It can be seen that the most unstable mode is the (m, n) = (4, 1) double tearing mode, dominant primarily at the outer q = 4 surface as shown in figure 6(b). Being a relatively low density plasma, the Alfvén timescale is much smaller ($t_a \sim a\sqrt{\rho \mu_0}/B \sim {0.2 \times 10^{-7}}\ \mathrm{s}$) and hence the tearing mode growth rate is relatively faster than in a typical density scenario. Via a test in the linear phase, we confirmed that the n = 1 growth rate using parallel diffusion is similar to that with advection at the speed of light (less than 2% difference). Running the full simulation with advection is possible, but computationally too expensive with the existing numerical methods. Additionally, the growth rate is found to increase with parallel transport in line with the prediction of the linear theory by Zhao et al [23]. Growth of the n = 1 mode triggers the linear growth of the subsequent toroidal modes in order until saturation. Subsequent non-linear interaction triggers stochastization of the magnetic field. Magnetic island structure at the early non-linear phase ($t = {0.97}\ \mathrm{ms}$) with the (4, 1) and (5, 1) islands is shown in the poincare plot in figure 7 (panel (a)). It is observed that stochastization starts initially around the outer q = 4 surface (panel (b)), and extends inwards toward the core (panel (c)). REs are lost as the field stochastizes. Since REs are lost much faster from the stochastic regions than from within the island regions, differential loss of REs can be observed from the localized high density of REs in the islands (panel (c)). Subsequent to this, the core (region still retaining the closed flux surfaces) shrinks until it reaches a size of a few centimeters of diameter (panels (d) and (e)). This observation matches well with the core-shrinking also seen during the final crash in the experiment. By this stage, about 94% of the REs are lost from the plasma. Stochastic loss of REs happens over a total time span of approximately ${100}\ \mu\mathrm{s}$. From the simulations, the core-shrinking is seen to occur as a series of instabilities of sequentially higher poloidal mode numbers. For example, between the span of a few tens of microseconds between panels (c) and (e) of figure 7, islands corresponding to n = 5, n = 6 and n = 7 are broken successively leading to the core-shrinkage.

Zoom In Zoom Out Reset image size

The corresponding evolution of total plasma current and RE current is shown in figure 8. The quick loss of RE current is accompanied by a reformation of thermal current during this duration via the induction effect. A small I_p spike ≈3% is observed in the simulation as a result of magnetic helicity conservation, which is similar to the I_p spike that appears after the thermal quench of a disruption. The time resolution of the current from the experiment is insufficient to conclude if such a short timescale I_p spike actually occurred in the shot. However, I_p spikes (of magnitudes between 1% and 8% of I_p ) have been observed in other discharges from the series of experiments described in [12, 20]. After the crash, the thermal plasma current decays at the resistive timescale. Accurate modeling of the downstream process requires the inclusion of Ohmic heating, radiative losses etc in the simulations and is not of interest for the present work.

Zoom In Zoom Out Reset image size

The RE loss during the crash is accompanied by a current profile flattening due to fast magnetic reconnection [24]. This is shown in figure 9(a) along with the corresponding q-profiles in figure 9(b) at various instances of time during the losses. It is interesting to note the minor-radially inward shift of the peak of the current profiles with time (or the location of high shear in q profile). This correlates to the previously mentioned sequence of instabilities with increasing poloidal mode number that causes the core-shrinking. We now turn to the impact of the stochastic RE losses on the wall.

Zoom In Zoom Out Reset image size

Figure 10(a) shows the distribution of toroidally-averaged RE flux ($\mathrm{m^{2}\ s^{-1}}$) on the wall during the time span of the RE loss, as a function of the poloidal location and time. It can be seen that most of the RE flux occurs in the near limiter region until ψ_N ≈ 1.2 (or a poloidal distance of ${\approx}{50}\ \mathrm{cm}$), which can be considered as a poloidal broadening of RE flux. This partially explains the reason for the low impact of REs on the first wall in this termination event, in addition to the absence of RE re-acceleration due to lack of high-Z impurities. It must be noted that the use of parallel diffusion as an ad-hoc leads to RE losses deposited on both sides of the limiter point (up and down), as opposed to the way deposition would occur with parallel advection of REs. In our analysis, this effect is accounted for by mapping the RE flux on both sides of the limiter to a single side, as per the ψ_N value of the location. Furthermore, it is observed that large RE losses occur within a duration of ${\sim}{20}\ \mu\mathrm{s}$, which matches very well with observation from the experiment. The distribution of magnetic field line connection length to the boundary, which quantifies the degree of stochasticity, also supports this observation (see appendix C). Also, at regions away from the limiter point, there is a slight time delay in the peak RE flux, which can be attributed to reduced stochasticity and longer connection time to the wall in that region.

Zoom In Zoom Out Reset image size

Figure 10(b) shows the distribution of RE fluence (particles per square meter) on the wall during the time period shown in figure 9(a), as function of poloidal and toroidal locations. In addition to the poloidally broad deposition, there is variation of RE deposition along the toroidal direction. The poloidal broadening is due to the stochasticity setting in first at the plasma edge, unlike a situation with predominantly core-stochastization wherein the deposition would be concentrated close to the limiter point. The toroidal distribution is smooth without any localized hot-spots, and this is clearly attributed to the n = 1 mode that was dominant even through the non-linear phase.

It was verified that the overall dynamics of this beam termination event with RE parallel diffusion is very similar to that in a simulation with RE parallel advection at the speed of light (fluence could not be obtained from the advection simulation though).

The first self-consistent non-linear MHD simulation of an RE beam termination event in JET has been presented. Mutual interaction between MHD and REs is treated by the use of an RE fluid model. Useful insights have been obtained on the dynamics of the crash and RE loss. A low density plasma with a hollow RE current profile leads to fast MHD growth dominated by n = 1, stochastization and corresponding RE losses. This is accompanied by current profile flattening by fast magnetic reconnection. Results show a good match with experiment. RE deposition shows significant toroidal variation (but without local hotspots), and poloidal broadening that potentially reduces the RE load on the wall. This explains in part why no first-wall damage is observed in the experiment. Results presented in this work improve our understanding of specific aspects of MHD behavior and therefore contribute to the development of a safe RE beam termination scenario in high current tokamaks. This opens up new avenues toward devising a robust RE mitigation scheme for ITER.

The starting equilibrium used in our simulations is clearly representative rather than fully empirical, but sufficiently adequate for the objective of this work. Sensitivity studies confirm a weak dependence on the exact current profile. Furthermore, as previously mentioned, parallel diffusion was used as an ad-hoc to relativistic advection of REs in order to avoid excessive computational costs. Nevertheless, from test simulations with advection, it was confirmed that the dynamics remain very similar in both cases. For example, linear growth rates differed by less than 2% and the RE loss rates due to stochasticity were very similar as well. Additionally, while using parallel diffusion, sufficient care was taken to map the RE flux onto one side of the limiter, to mimic the way REs would be deposited on the first-wall in case of advection. Hence all the important physical effects (if not exact details) are still captured adequately with the use of ad-hoc parallel diffusion.

Our simulations further substantiate the broad theme proposed by Reux et al [12] on the potential utility of deuterium second injection. From this study, it is seen that an MHD instability causing strong stochastization first in the edge region (achieved by a hollow current density profile in this case) can be a very conducive ingredient to proactively enable scattered RE losses to avoid wall damage. However, it is vital to understand the physics behind the formation of a hollow profile to be able to reliably obtain it and assess the applicability of the D₂ second SPI scheme in ITER. We plan to study this aspect in the future.

Furthermore, in elongated plasma scenarios such as in ITER, possible vertical instability after the first or second injection makes the applicability of the scheme challenging. Hence it is important for future studies to perform RE termination studies in elongated shapes with free-boundary computations in order to incorporate vertical dynamics realistically. Major MHD codes have demonstrated the ability to do so, including 3D instabilities in a benchmark exercise [25].

The authors acknowledge help from G Szepesi for providing equilibrium data. Fruitful discussions with F J Artola and E Nardon are also acknowledged. We also acknowledge useful suggestions on the manuscript from S Gerasimov. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under Grant Agreement No. 633053. This work was also supported by the EUROfusion—Theory and Advanced Simulation Coordination (E-TASC). The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Normalization used for the variables shown in the governing equations (5a )–(5f ). The factors µ₀, n₀, ρ₀ refer to the magnetic permeability of free space, central number density and the central mass density of the background plasma, respectively, and k_B is the Boltzmann constant.

Quantity	Normalization
Time, t	$t^\mathrm{SI} = t \sqrt{\mu_0 \rho_0}$
RE density, nr	$n_{r}^\mathrm{SI} = n_{r} \sqrt{\rho_0/\mu_0} / \left(eR\right)$
Speed of light, c	$c^\mathrm{SI} = c/\sqrt{\mu_0 \rho_0}$
Electric potential, u	$u^\mathrm{SI} = u/ \sqrt{\mu_0 \rho_0}$
Toroidal vorticity, ω	$\omega_\phi^\mathrm{SI} = u/ \sqrt{\mu_0 \rho_0}$
Current density, j	$j_\phi^\mathrm{SI} = -j/\left(R \mu_0 \right)$
Resistivity, η	$\eta^\mathrm{SI} = \eta \sqrt{\mu_0/\rho_0}$
Viscosity, $\mu_{\perp}$	$\mu_{\perp}^\mathrm{SI} = \mu_{\perp} \sqrt{\rho_0/\mu_0}$
Diffusivity, $D_{r,\parallel or \perp}$	$D_{r}^\mathrm{SI} = D_{r} /\sqrt{\mu_0\rho_0}$
Density, ρ	$\rho^\mathrm{SI} = \rho \rho_0$
Temperature, T	$T\,^\mathrm{SI} = T/\left(k_B \mu_0 n_0\right)$

Some of the RE number density profiles (normalized) tested with the Green's functions are displayed in figure 11(a). Corresponding to the displayed profiles a series of synthesized images is shown in figure 11(b). Although only qualitative comparison has been conducted, it is obvious that any kind of peaked profile cannot create the observed radiation pattern given the diagnostics parameters.

Zoom In Zoom Out Reset image size

The evolution of the connection length distribution during the stochastic RE loss phase is shown in figure 12. Smaller connection length corresponds to a higher degree of stochasticity. It can be observed that stochasticity in the radially inner regions only last for ${\sim}{20}\ \mu\mathrm{s}$, as also indicated by the flux deposition.

Zoom In Zoom Out Reset image size