Runaway dynamics in tokamak disruptions with current relaxation

Abstract

The safe operation of tokamak reactors requires a reliable modelling capability of disruptions, and in particular the spatio-temporal dynamics of associated runaway electron currents. In a disruption, instabilities can break up magnetic surfaces into chaotic field line regions, causing current profile relaxation, as well as a rapid radial transport of heat and particles. Using a mean-field helicity transport model implemented in the disruption runaway modelling framework Dream, we calculate the dynamics of runaway electrons in the presence of current relaxation events. In scenarios where flux surfaces remain intact in parts of the plasma, a skin current is induced at the boundary of the intact magnetic field region. This skin current region becomes an important centre concerning the subsequent dynamics: it may turn into a hot ohmic current channel, or a sizeable radially localized runaway beam, depending on the heat transport. If the intact region is in the plasma edge, runaway generation in the countercurrent direction can occur, which may develop into a sizeable reverse runaway beam. Even when the current relaxation extends to the entire plasma, the final runaway current density profile can be significantly affected, as the induced electric field is reduced in the core and increased in the edge, thereby shifting the centre of runaway generation towards the edge.

1. Introduction

Recent progress in fusion science has increased confidence in that the energy confinement requirements of upcoming experiments aiming to demonstrate a scientific break-even of the power balance, such as ITER and SPARC, can be met during normal operation (Rodriguez-Fernandez et al. 2020). However, off-normal events leading to the rapid loss of the energy content of a plasma, called disruptions (Hender et al. 2007; Boozer 2012; Lehnen et al. 2015), and associated generation of highly energetic runaway electron beams, represent an outstanding challenge of the otherwise most successful tokamak concept for magnetic confinement.

In a disruption, the thermal energy of the plasma is lost during the rapid thermal quench, followed by a resistive decay of the ohmic plasma current, called current quench. Parallel electric fields, well in excess of the Connor–Hastie critical electric field (Connor & Hastie 1975) for runaway generation, may then be induced, generating a runaway electron seed (Dreicer 1959; Harvey et al. 2000; Martın-Solıs, Loarte & Lehnen 2017) which can, in a collisional avalanche process, exponentiate into a macroscopic runaway current (Sokolov 1979; Rosenbluth & Putvinski 1997; Hesslow et al. 2019a). The amplification factor of the runaway current in the avalanche process is exponentially sensitive to the initial plasma current (Hender et al. 2007; Boozer 2018a); thus high-current, reactor-relevant devices, such as ITER, are predicted to be prone to convert a large fraction of their ohmic current into runaway electron current (Boozer 2012; Breizman et al. 2019; Vallhagen et al. 2020). The resulting multi-${\rm MA}$ runaway beam poses a major threat to the structural integrity of the plasma-facing components. Methods to mitigate the effects of a disruption (Lehnen et al. 2015) focus on limiting the exposure of the wall to localized heat losses and to the impact of high-current runaway electron beams, as well as on avoiding excessive forces on structural elements. The currently promoted disruption mitigation methods employ magnetic perturbations (Tinguely et al. 2021), and/or a rapid massive material injection, such as through shattered pellet injection (Commaux et al. 2016; Jachmich et al. 2021; Vallhagen et al. 2022), to tailor the spatio-temporal properties of dissipating the thermal and magnetic energy content of the plasma.

Disruption is a strongly nonlinear multi-physics process, with wide ranges of temporal and spatial scales and kinetic physics at play. Thus progress in developing robust disruption-mitigation strategies must inevitably rely on employing a broad range of numerical tools in terms of the trade-off between physical accuracy and numerical cost; these include three-dimensional (3-D) magnetohydrodynamic (MHD) simulation codes (Sovinec et al. 2004; Ferraro et al. 2016; Hoelzl et al. 2021), flux surface-averaged transport solvers (Linder et al. 2020), local kinetic models (Stahl et al. 2016), and sometimes a combination thereof (Hoppe et al. 2021b). The thermal quench often starts with an instability that breaks up the confining magnetic surfaces, leading to a highly elevated transport along multiple channels (Hender et al. 2007), including electron heat losses, transport of energetic particles and a redistribution of current density. This is a very consequential process concerning the rest of the disruption, while it is perhaps the most difficult phase of the disruption to model. Models that compromise on spatial dimensionality cannot self-consistently describe this process, while even in 3-D MHD simulations it is a challenge (Nardon et al. 2021; Boozer 2018b) to resolve the process at physically low resistivity.

A key mechanism of the destruction of flux surfaces during the thermal quench is fast magnetic reconnection (Boozer 2019), a rapid topological change of the magnetic field occurring on an ideal MHD time scale, independently of resistivity. During this process, field lines connect regions initially belonging to different magnetic surfaces with different current density and pressure, and the variation of $j_\|/B$ (parallel current density over magnetic field strength) along the field line is relaxed, mediated by the propagation of shear Alfvén waves. The effects of such current relaxation events – which do not necessarily extend over the entire plasma volume – are experimentally observed through a sizeable temporary increase of the total plasma current $I_p$ (Wesson, Ward & Rosenbluth 1990; de Vries et al. 2015), as well as the reduction of the internal inductance of the plasma $\ell _{i}=W_{\rm p}/(\mu _0 R_0 I_p^2/4)$, such that the magnetic energy in the poloidal magnetic field $W_p$ is only slightly reduced. Here $\mu _0$ is the vacuum permeability and $R_0$ is the major radius at the magnetic axis of the unperturbed magnetic field. The low resistivity and the ideal MHD time scale imply that the process approximately conserves magnetic helicity (Taylor 1974; Berger 1999), a topological measure of magnetic field linkage that can only change on global resistive time scales. Magnetic helicity $H^{M}$ is the volume integral of $\boldsymbol {A}\boldsymbol {\cdot }\boldsymbol {B}$ with the magnetic vector potential and field, $\boldsymbol {A}$ and $\boldsymbol {B}$, which can be expressed in tokamak geometry as $H^{M}=-2\int \psi \, {\rm d} \psi _t$, with the poloidal and toroidal magnetic fluxes, $\psi$ and $\psi _t$, defined as in figure 2 of Boozer (2017), and the entire plasma being the integration domain.

During the thermal quench the diagnostics used to measure the current density are unreliable due to the excessive photon background; in practice the total plasma current and – in the case of shaped plasmas – the internal inductance are the quantities that can be experimentally obtained during a disruption. Thus, for practical reasons, several disruption runaway electron studies use the pre-disruption current profile (Linder et al. 2020) or assume a flat electric field profile (Insulander Björk et al. 2021; Hoppe et al. 2021b) as their initial condition, which is then evolved in an inductive–resistive fashion. However, the rearrangement of the current density in a fast magnetic reconnection event, and in particular the development of skin currents (Wesson et al. 1990) at boundaries between chaotic field line regions and intact flux surfaces (Boozer 2017), is expected to have a significant effect on the development of the runaway current. It has already been shown that the transport of runaway electrons – also caused by the chaotic magnetic field – has a major impact on the evolution of the runaway current density (Svenningsson et al. 2021), but studies concerning the effect of the current relaxation are lacking. To capture current relaxation through fast magnetic reconnection in a model that only retains the radial (i.e. across the unperturbed magnetic flux surfaces) variations, one may employ the mean-field model of Boozer (2018a), which describes the current redistribution as a hyperdiffusion of the poloidal flux, and is constructed to conserve magnetic helicity.

In this paper we study the effect of current relaxation in the thermal quench phase of tokamak disruptions on the evolution of the ohmic and runaway current components. The plasma parameters considered are ITER-like, and deuterium and neon are assumed to be injected as part of the disruption mitigation. We consider both full-radius and partial reconnection events, with a particular focus on the fate of skin current regions in the latter case. We employ the recently developed Dream code (Hoppe, Embreus & Fülöp 2021a), equipped with an implementation of the helicity transport model of Boozer (2018a), which allows the self-consistent simulations of plasma cooling and associated runaway electron (RE) dynamics during disruptions. We find that full-radius reconnection tends to shift the centre of the runaway generation radially outward. Possible outcomes include a reduced core-localized runaway current, or significantly increased runaway currents with hollow radial profile. We also show that skin current regions, when developed, play a central role in the subsequent current evolution, and could turn into long-lived hot ohmic channels or runaway hotspots. As a possible, somewhat exotic, outcome, ‘reverse’ (countercurrent) runaway beams may develop in intact edge regions.

The rest of the paper is organized as follows. First, in § 2.1, we provide the equations relevant for the current evolution in the numerical model, then in § 2.2 we describe the baseline simulation set-up. We first consider full-radius current relaxations; the possible outcomes – core-localized and edge-localized runaway generation – are discussed in §§ 3.1 and 3.2, respectively. The role of heat transport concerning the fate of a core-localized skin current region is discussed in § 3.3. Finally, we consider the possibility of an intact edge and associated reverse runaway generation in § 3.4. Limitations and future directions are touched upon in § 3.5, and the conclusions are summarized in § 4.

2. Methods

2.1. Numerical approach

For our simulations we use the Dream tool, which is a finite-volume fluid-kinetic framework developed to model runaway electron dynamics in disruptions. For a detailed description of the code we refer the reader to Hoppe et al. (2021a), while here we overview some aspects particularly relevant for our analysis. The code resolves one spatial dimension – a radial coordinate $r$ – and it can resolve the entire gyro-averaged momentum space of electrons, or parts of it, parametrized by the magnitude of the momentum $p$ and $\xi =p_\|/p$, with the component of the momentum along the magnetic field $p_\|$ (taken at the lowest magnetic field position along a collisionless orbit). It has various reduced models for the momentum-space dynamics as well; here we predominantly use a fluid model that retains only a thermal bulk, characterized by a density $n_e$, a temperature $T_e$ and an ohmic current density $j_{\varOmega }$, and a runaway electron population, characterized by the current density they carry $j_{\rm re}$. In the fluid model it is assumed that the runaways, which have density $n_{\rm re}$, move with the speed of light $c$ parallel to the magnetic field, hence $j_{\rm re} =e n_{\rm re} c$, where $e$ is the elementary charge.

The total current density is computed through evolving the poloidal flux:

(2.1)\begin{equation} \frac{\partial \psi}{\partial t}={-}\mathcal{R}+\mu_0\frac{\partial}{\partial\psi_{t}}\left( \psi_{t} \varLambda_{m} \frac{\partial}{\partial \psi_{t}}\frac{j_{\rm tot}}{B}\right); \end{equation}

in this equation the toroidal flux $\psi _{t}(r)=(2{\rm \pi} )^{-1}\int _0^r V' \langle \boldsymbol {B} \boldsymbol {\cdot } \boldsymbol {\nabla } \varphi \rangle \, \mathrm {d} r'$ is used as a radial coordinate, with $\langle \cdot \rangle$ denoting flux surface average, $V'$ is the spatial Jacobian (incremental volume enclosed by two infinitesimally close flux surfaces) and $\mathcal {R}=2 {\rm \pi}\langle \boldsymbol {E}\boldsymbol {\cdot }\boldsymbol {B}\rangle /\langle \boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\varphi \rangle$, with the toroidal angle $\varphi$. This flux function term is equal to the loop voltage when the second term on the right-hand side of (2.1) vanishes, and it describes ohmic dissipation, as we may express the term $\langle \boldsymbol {E}\boldsymbol {\cdot }\boldsymbol {B}\rangle =j_{\varOmega } \langle B^2\rangle /(\sigma B)$, with $\sigma$ the electric conductivity and $j_{\varOmega }$ the ohmic current density.¹ The total current density $j_{\rm tot}$ includes the ohmic and the runaway current densities, $j_{\varOmega }$ and $j_{\rm re}$, and in some kinetic simulations $j_{\rm hot}$ as well, the current density carried by electrons at an intermediate energy range between the thermal bulk and the highly energetic runaways. The total current density is related to the poloidal flux through

(2.2)\begin{equation} 2{\rm \pi} \mu_0 \langle \boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla} \varphi \rangle \frac{j_{\rm tot}}{B}=\frac{1}{V'} \frac{\partial}{\partial r}\left[ V' \left\langle \frac{|\boldsymbol{\nabla} r|^2}{R^2}\right\rangle \frac{\partial \psi}{\partial r} \right], \end{equation}

where $R$ is the major radius of a given point. From the form of (2.2) we can see that the term involving $\varLambda _{m}$ in (2.1) includes a fourth-order radial derivative of $\psi$, and as such it describes a hyperdiffusion of $\psi$ wherever $\varLambda _{m}$ is non-zero. This term describes a local transport of magnetic helicity, such that it conserves the total magnetic helicity. That the form of (2.1) is constructed to respect this conservation property is most easily seen by considering $d_t H^{M}=-2 \int \,\mathrm {d} \psi _t\,\partial _t \psi$. Clearly, the contribution of the second term on the right-hand side of (2.1) to $\partial _t \psi$ – which is of the form $\partial _{\psi _t}(\cdots )$ – integrates to zero if it vanishes at the boundaries; that is, helicity can only be transported across the plasma boundary, but not created or destroyed inside. The numerical conservation of $H^{M}$ is demonstrated in Appendix B. In the simulations, the helicity transport coefficient $\varLambda _{m}$ can have an arbitrary prescribed spatio-temporal structure.

Currents in passive structures are taken into account through the boundary condition on (2.1). The relevant equations are the following:

(2.3)\begin{gather} \psi_{\rm edge} \equiv \psi(r=a)= \psi_{\rm wall} - M_{\rm ew} I(a), \end{gather}

(2.4)\begin{gather} \psi_{\rm wall} \equiv \psi(r=b)={-}L_{\rm ext}\left[I(a) + I_{\rm wall}\right], \end{gather}

(2.5)\begin{gather} \mathcal{R}^{\rm wall} = R_{\rm wall}I_{\rm wall}, \end{gather}

where $a$ and $b$ are the minor radii corresponding to the plasma edge and the wall. Parameter $I(r)$ is the total plasma current within a given radius $r$. The edge–wall mutual inductance is $M_{\rm ew}=(2{\rm \pi} )^2\mu _0\int _a^b\,{\rm d}r(V'\langle |\boldsymbol {\nabla } r|^2/R^2 \rangle )^{-1}$, which in the cylindrical limit reduces to $\mu _0R_{0} \ln (b/a)$. The external inductance is $L_{\rm ext}=\mu _0R_{0} \ln (R_{0}/b)$, with $R=R_{0}$ at the magnetic axis. The inputs to the wall model are the resistive time scale of the wall $\tau _w=L_{\rm ext}/R_{\rm wall}$ and $b$.

2.2. Simulation set-up

The disruption simulations assume an initially ($t<0$) pure deuterium–tritium plasma (with even isotope concentrations² ). Specifically, the initial electron density is spatially constant $10^{20}\,{\rm m}^{-3}$, the temperature is parabolic with $20\,{\rm keV}$ on-axis and the current density corresponds to a total plasma current of $15\,{\rm MA}$. The simulations consider an ITER-like magnetic geometry with $R_0=6\,{\rm m}$, $a=2\,{\rm m}$, $b=2.15\, {\rm m}$, $B(r=0)=5.3\,{\rm T}$ and a resistive wall time of $\tau _w=0.5\,{\rm s}$, as well as a Miller model equilibrium (Miller et al. 1998), with the radially varying realistic shaping parameters given in Appendix A.

At $t=0$ an instantaneous and homogeneous deposition of additional neutral deuterium and neon is done. At the same time an elevated transport of electron heat and energetic electrons is activated, along with a current profile relaxation, as described by the appropriate transport coefficients, to emulate fast magnetic reconnection and associated break-up of flux surfaces. At $t=6\,{\rm ms}$, when the maximum plasma temperature has dropped to $\approx 100\,{\rm eV}$ and the current density has already relaxed, the transport of energetic electrons and magnetic helicity is switched off, and the electron heat transport is strongly reduced. This is to account for the reformation of flux surfaces as the drive of the instability is removed. A weaker electron heat diffusivity remains active until the end of the simulation. Toroidicity effects are accounted for in the conductivity through the model of Redl et al. (2021) and in the runaway generation mechanisms.

The Dream simulations are performed in fully fluid mode, unless stated otherwise. The Dreicer runaway generation rate is calculated using a neural network (Hesslow et al. 2019b), Compton scattering and tritium decay seed sources are accounted for as in Martın-Solıs et al. (2017) and Vallhagen et al. (2020), the hot-tail seed is calculated using the model in § 4.2 of Svenningsson (2020) and the avalanche growth rate accounts for partial screening (Hesslow et al. 2019a).

The bulk electron temperature evolution is calculated from the time-dependent energy balance throughout the simulation, according to equation (43) in Hoppe et al. (2021a), accounting for ohmic heating, line and recombination radiation and bremsstrahlung, as well as a radial heat transport. Recombined deuterium is assumed to be opaque to Lyman line radiation, which can have a non-negligible effect on the post-thermal quench plasma temperature and indirectly on the avalanche gain (Vallhagen et al. 2022). Note that here we do not evolve ion temperatures separately and do not have a kinetic runaway population; thus the electron–ion heat exchange term and the kinetic term describing heating by runaway electrons are zero. However, the latter process is approximately accounted for by a term $j_{\rm re} E_{c}$, with $E_{c}=e^3n_e \ln \varLambda _c/(4{\rm \pi} \epsilon _0 m_e c^2)$ the critical electric field, $\epsilon _0$ the vacuum permittivity and $m_e$ the electron mass. In the definition of $E_{c}$ we take $\ln \varLambda _c$, the momentum-dependent Coulomb logarithm, defined in equation (18) of Hoppe et al. (2021a), at the momentum $p=20\, m_e c$; in addition, $\ln \varLambda _c$ depends on the bulk electron density and temperature, $n_e$ and $T_e$.

In the early phase of the thermal quench electron heat losses are dominated by a transport along the chaotic magnetic field lines. To capture this, we use a Rechester–Rosenbluth-type model in the collisionless limit (Rechester & Rosenbluth 1978), with a heat diffusivity given by $D_W\approx 2\sqrt {{\rm \pi} } R_0 v_{te} (\delta B/B)^2$, where $v_{te}=\sqrt {2 T_e/m_e}$ is the local electron thermal speed and $\delta B/B$ is the relative magnetic perturbation amplitude³; specific numerical values of $\delta B/B$ are provided in § 3.1. This model is based on the assumption of chaotic magnetic field line topology, where $\delta B$ is a representative value for the radial perturbation amplitude. Note, however, that this heat transport coefficient is applied spatially homogeneously even in the cases where part of the plasma is assumed to have intact magnetic surfaces (where no current relaxation and runaway transport is active). This is so as to be able to reach sufficiently low temperatures to initiate the radiative temperature collapse without the need to modify the impurity content, thereby making it easier to compare various cases. Also, technically the same form of transport coefficient with reduced $\delta B/B$ is applied after the thermal quench, throughout the simulation, while physically electron heat transport might then be caused by electrostatic turbulence. The value of the post-thermal quench heat diffusivity is only important in the presence of intact regions in the plasma, as discussed in § 3.3.

During the thermal quench we also account for a diffusive transport of runaway electrons using a diffusion coefficient of similar form but with parallel streaming along the perturbed field lines at the speed of light⁴ $D_{\rm RE}={\rm \pi} R_0 c (\delta B/B)^2$. We employ here the same $\delta B/B$ as for electron heat transport, but only in the chaotic field line regions.

The simulations use $200$ radial grid cells, and in cases with skin current regions, $150$ of these are packed around the radial region with sharp current density variation. The first microsecond of the simulation, when the injected neutral material ionizes, is resolved by $1000$ time steps. The rest of the simulation typically uses a time step⁵ of $\Delta t= 6.6\,\mathrm {\mu }{\rm s}$.

3. Results

3.1. Full-radius current relaxation; core-localized runaways

First we briefly discuss the temperature evolution through our baseline case; in most cases a similar qualitative behaviour is observed, with the exception of cases with strong ohmic skin currents, discussed in § 3.3. After the instantaneous and homogeneous deposition of neutral deuterium and neon of densities⁶ $n_{{\rm D,inj}}=7\times 10^{20}\,{\rm m}^{-3}$ and $n_{\rm Ne,inj}=1.5\times 10^{19}\,{\rm m}^{-3}$, the plasma temperature drops on a time scale of $0.01\unicode{x2013}0.1\,\mathrm {\mu }{\rm s}$ as these species get ionized and the plasma diluted. The peak temperature drops from $20$ to ${\sim }2\,{\rm keV}$. Then the electron heat transport corresponding to a magnetic perturbation amplitude of $\delta B/B = 3.5\times 10^{-3}$ cools the plasma further on a millisecond time scale, until $T_e$ reaches $\sim 100\,{\rm eV}$, at which point radiative energy losses become dominant and cool the plasma further to the $1$–$10\,{\rm eV}$ range very rapidly. This specific value of $\delta B/B$ is chosen such to target a millisecond characteristic cooling time scale while radial transport dominates the heat losses. The injected densities are chosen to produce a disruption in the parameter region deemed favourable in Vallhagen et al. (2020), producing relatively low maximum runaway currents and current quench time scales within tolerable limits.

Once the thermal quench is complete the electron heat transport is reduced to a level corresponding to $\delta B/B = 4\times 10^{-4}$ (at $6\,{\rm ms}$, indicated by the dashed vertical line in figure 1a). Then the temperature evolves on the ${\sim }10\,{\rm ms}$ time scale of the current quench, and it is determined by an approximate balance of ohmic heating (later the friction of energetic electrons on the bulk) and radiative losses.

Figure 1. Comparison of simulations without (thin light curves) and with (thick dark curves) a full-radius current relaxation. The latter case employs spatially homogeneous helicity and runaway electron transport coefficients in the first $6\,{\rm ms}$ (marked by dotted vertical line in a). (a) Time evolution of the total plasma current (dashed) and its ohmic (solid) and runaway (dash-dotted) components. (b) Radial profiles of the electric field normalized to the local Dreicer field at $t=6\,{\rm ms}$. (c) Radial profiles of the current density, with ohmic (solid) and runaway (dashed) components. Profiles at $t=\{0,\,6,\, 80\}\,{\rm ms}$ are shown with orange, purple and black lines. Note that the runaway component is significant only in the third time point, and the ohmic component in the first two time points.

Without current relaxation (thin light lines in figure 1), the total ohmic current (solid line in figure 1a) decreases monotonically. After a short period of hot-tail runaway generation that is concluded at $5\,{\rm ms}$, the dominant seed generation mechanism is tritium decay, while the Dreicer and Compton seeds remain negligible in comparison. Once the runaway current (dash-dotted) reaches $\sim {\rm MA}$ level through avalanche, the ohmic current drops on the time scale of the avalanche growth, and the runaway current saturates around $5.3\,{\rm MA}$. Then the runaway electron current keeps increasing only very slowly as poloidal magnetic field is diffusing back into the vacuum chamber through the resistive wall.

In the case where a spatially homogeneous hyperdiffusivity of $\varLambda _{m} = 3\times 10^{-2} \, \,{\rm Wb}^2\,{\rm m}\,{\rm s}^{-1}$ is applied in the first $6\,{\rm ms}$ (solid dark lines in figure 1), when the flux surfaces are assumed to be broken up, we find that the ohmic current first exhibits a peak of $18.5\,{\rm MA}$ at $3\,{\rm ms}$. This is the ‘$I_p$ spike’ regularly observed in plasma disruptions, as the current density is radially redistributed at approximately constant magnetic helicity. The ohmic current density indeed flattens,⁷ as seen in figure 1(c) (thick purple line).

The ohmic current is sustained by the induced electric field, which thus needs to adjust itself to the current redistribution. This means that the electric field must increase at the edge and reduce in the core compared to the case without current relaxation, which is clearly seen in figure 1(b), showing the electric field normalized to the Dreicer field $E_D=(e^3 n_e \ln \varLambda _c)/(2{\rm \pi} \epsilon _0^2 m_e v_{te}^2)$. We note that the changes in $E/E_D$ are more affected by differences in $E$ rather than $E_D$, as changes in the temperature profiles are modest.

That current relaxation can lead to strongly disparate outcomes concerning the runaway current evolution can be better understood from the observation that the runaway electron density profiles often grow as to replace the ohmic current; thus the initial ohmic current density acts as an envelope to the final runaway current density (even if there are exceptions from this rough rule of thumb, due to trapping, finite wall time and electric field diffusion). Indeed, we find that the final runaway current densities (dashed lines in figure 1c) are everywhere lower than the initial ohmic current density. In the case with current relaxation the same observation can be made, but now with the relaxed current profile playing the role of the envelope. As in both cases the runaway density peaks in the core, this translates to a reduced final runaway electron current in the case with current relaxation (from $5.3$ to $3.9\,{\rm MA}$).

In the first $6\,{\rm ms}$ of the case with current relaxation runaway electron transport is also active, corresponding to $\delta B/B = 3.5\times 10^{-3}$ (consistently with the electron heat transport). This reduces the runaway electron seed, as seen in figure 1(a). However, after the flux surfaces are re-formed (following the dotted vertical line) the runaway current growth quickly recovers, and it is not as much due to the reduced seed, rather due to a slightly reduced avalanche growth rate, that the runaway current reaches macroscopic values later in the case with current relaxation.

3.2. Full-radius current relaxation; edge-localized runaways

To illustrate that the current flattening could potentially lead to a very different outcome, we consider another case with an ITER-sized plasma (with initial profiles different from those of the baseline). In this case the initial profiles and the simulation set-up is such that they favour Dreicer seed generation and the development of a higher electric field in the edge. Such differences include a higher current density, a relatively low post-thermal quench temperature and a lower edge electron density. This case has a reduced physics fidelity, and as such, it is likely less representative of the behaviour in ITER. For instance, it uses a prescribed temperature variation, the evolution of ion charge states and runaway electron transport are disabled, and no shaping and toroidicity effects are included. More details of corresponding settings are provided in Appendix A.

Unlike in our baseline case, where $E/E_D$ went from centrally peaked to flat when current relaxation was applied, in this case it goes from mostly flat to peaked towards the edge (compare figure 2a with figure 1b). Note also the significantly higher values of $E/E_D$, resulting in a rapid, and almost full, ohmic-to-runaway current conversion. In the case without current relaxation the final runaway current is peaked in the core, as shown by the thin lines in figure 2(b), similarly to the baseline. Now the initial current density does not quite envelope the final runaway current, due to an electric field diffusion into regions where the runaway current is already high. More importantly, when a radially constant $\varLambda _{m} =1.5\times 10^{-2} \,{\rm Wb}^2\,{\rm m}\,{\rm s}^{-1}$ is applied in the first $10\,{\rm ms}$ of the simulation (solid lines), the runaway current profile becomes edge-localized and hollow. In fact, the runaway current first grows at the edge, and once it replaces the ohmic current density it continues to grow radially inward. Since the weight of the current density towards the total current increases radially, the final total runaway current corresponding to the case with current relaxation is higher than that without current relaxation. This is unlike the situation in our baseline case shown in figure 1.

Figure 2. Alternative scenario that yields a hollow runaway profile with current relaxation. (a) Radial profiles of the electric field normalized to the Dreicer field at $t=4\,{\rm ms}$. (b) current density profiles with ohmic (solid) and runaway (dashed) components. Thinner lines with lighter colour correspond to a simulation without current relaxation; thicker darker lines use a spatially homogeneous $\varLambda _{m}$, applied in the first $10\,{\rm ms}$. Profiles at $t=\{0,\,4,\, 15\}\,{\rm ms}$ are shown with orange, purple and black lines. Note that the runaway component is significant only in the third time point, and the ohmic component in the first two time points.

We may conclude, then, that a current relaxation that extends over the entire plasma – reducing the electric field in the core and increasing it in the edge – tends to decrease the final runaway current as long as the conditions for runaway generation are more favourable in the core, and increase it if the plasma is more prone to develop edge-localized runaway currents.

3.3. Intact core and the role of heat transport

We now consider the possibility of a partial current relaxation, where a range of flux surfaces in the deep core remain intact, while the magnetic field becomes chaotic in most of the plasma up to the edge. As often observed in 3-D MHD disruption simulations (Sommariva et al. 2018; Bandaru et al. 2021; Artola et al. 2022), the deep plasma core tends to be a region where flux surfaces are not completely destroyed, or are broken up only for a very short period within the thermal quench. This may be related to the typically low magnetic shear in the deep core, which does not favour island overlap.

To model such a scenario, we assume that the magnetic surfaces are broken up in the first $6\,{\rm ms}$, when we employ a runaway and an electron heat transport corresponding to $\delta B/B = 3.5\times 10^{-3}$, as well as a hyperdiffusivity of $\varLambda _{m} = 3\times 10^{-2} \, \,{\rm Wb}^2\,{\rm m}\,{\rm s}^{-1}$, as in the case considered previously. However, now the current relaxation and the runaway electron transport are only active in the outer part of the plasma. The transport coefficients, $D_{\rm RE}$ and $\varLambda _{m}$, now denoted by $X$, jump between their core and edge values $X_{\rm core}$ and $X_{\rm edge}$ according to the following radial variation:

(3.1)\begin{equation} X= \dfrac{X_{\rm core}\left\{{\rm Erf}\left[\dfrac{a-r_m}{w}\right] + {\rm Erf}\left [\dfrac{r_m-r}{w}\right] \right\} +X_{\rm edge}\left\{{\rm Erf}\left[\dfrac{r_m}{w}\right] + {\rm Erf}\left[\dfrac{r-r_m}{w}\right] \right\} }{{\rm Erf}\left[\dfrac{a-r_m}{w}\right] + {\rm Erf}\left[\dfrac{r_m}{w}\right] }, \end{equation}

with the error function ${\rm Erf}$, and here $X_{\rm core}=0$, $X_{\rm edge}$ corresponding to the values given above, a transition radius $r_m=0.3\,{\rm m}$ and a characteristic transition width of $w=0.01\,{\rm m}$. Again, for easier comparison with other cases, we keep $D_W$ radially constant.

With intact regions, and associated skin currents, remnant (i.e. pre-thermal quench) heat transport has a major role in the dynamics. First we consider a lower remnant electron heat transport at $\delta B/B=4\times 10^{-4}$. Similarly to figure 1(c,a), we observe a rapid flattening of the current density throughout most of the plasma, and an associated $I_p$ spike. At the same time, a strongly peaked skin current appears at the boundary of the intact region, shown in the zoomed-in logarithmic plot of figure 3(a). This skin current broadens out with time due to electric field diffusion, and turns into a long-lived ohmic current channel, such that the ohmic decay of $I_p$ halts. In the meantime avalanche, slowly but steadily, leads to an increasing runaway current. This runaway current grows in a very narrow layer around $r=0.25\,{\rm m}$, where $E/E_D$ spikes first during the formation of the skin current. In this region the hot-tail seed generation peaks around $0.5\,{\rm ms}$, reaching a maximum runaway rate of $2.7 \times 10^{17}\,{\rm m}^{-3}\,{\rm s}^{-1}$, which is seven orders of magnitude higher than the maximum of the tritium seed runaway rate that has the second highest value among the seed generation mechanisms. As the seed generation is exponentially sensitive to the electric field, the seed in the skin layer is many orders of magnitude higher than elsewhere. Since there is no transport of runaways in the intact region, the generated seed does not broaden radially. Consequently, essentially all runaway generation in the intact region happens in this very narrow peak throughout the simulation, even though $E/E_D$ is not the highest in this region after the current relaxation.

Figure 3. Dynamics of the intact core case with low (a–c) and high (d–f) remnant heat diffusivity (zoom-in around skin current region; note the logarithmic scales). (a,d) Radial profiles of the current density at various time points, with ohmic (solid) and runaway (dashed) components. (b,e) Radial profiles of temperature at various time points. (c,f) Power balance showing ohmic heating (dotted), heating by runaways (orange solid), radiative losses (red solid) and net heat transport out from (blue dashed) and into (green dashed) a given radial location; taken at $t=60\,{\rm ms}$.

The ohmic current channel in the intact region broadens inward in time, and it locally heats the plasma above $100\,{\rm eV}$, as seen in figure 3(b). Once the runaway current sheet increases to macroscopic values at the outer edge of the ohmic channel, the ohmic current density peak moves inward, along with the temperature peak. It is interesting to consider the inner structure of the hot ohmic channel in terms of power balance, shown in figure 3(c), taken at $t=60\,{\rm ms}$. The main heat source inside the temperature peak is ohmic heating (black dotted line). As the temperature in the middle of the peak is too high for radiative losses to be dominant, ohmic heating is balanced by an outward diffusion of heat there (blue dashed). Towards the sides of the temperature peak this transported heat is being deposited (dashed green), and as the temperature drops rapidly outward of the peak – along with the ohmic heating – it is the divergence of the outward heat flux that becomes the dominant heat source, balanced by radiative heat losses (solid red) that are more effective below $100\,{\rm eV}$. We can also see the contribution of the friction of runaway electrons on the bulk to the heating (orange solid), but at this time point this contribution is subdominant to the ohmic heating.

The hot current channel exhibits strong similarities to the hot current filaments discussed in detail by Putvinski et al. (1997). Such a formation with extremely large temperature gradients at its edges are likely to be unstable to microinstabilities, and thus its existence would likely be of transient nature, if it could be formed at all. Nevertheless, they are acceptable solutions within the model employed here, where the heat diffusivity is simply prescribed. Heat transport driven by microinstabilities tends to be stiff, with fluxes rapidly increasing above a critical gradient; however, in this work any excess of instability thresholds is not monitored.

To account for the activity of such instabilities, we may simply increase the remnant heat diffusivity, which can essentially change the dynamics of the skin current region. Plots similar to figure 3(a–c), but corresponding to a remnant electron heat diffusivity at $\delta B /B = 2 \times 10^{-3}$ (compared with $\delta B/B=4\times 10^{-4}$), are shown in figure 3(d–f). This time the ohmic current quench completes within $60\, {\rm ms}$, since the formation of a long-lived hot ohmic current channel is inhibited in the presence of the higher heat diffusivity. Clearly the hot current channel does not form, as seen in the $j$ and $T_e$ plots (figure 3d,e). Once the current relaxation is over the electric field also drops to small values, while it still remains significant enough to drive the avalanche, with the runaway current exceeding $1\,{\rm MA}$ before $20\,{\rm ms}$. Most of the runaway current is carried by a narrow current channel where the seed formation was the strongest, similarly to the case with lower heat diffusivity. However, the power balance has an essentially different character from the lower-heat-diffusivity case; compare figure 3(c,f). As the ohmic current has decayed to small values, the contribution of ohmic heating (black dotted) is negligible compared with the heating from runaway friction (orange solid). This heating is balanced by a transport outward from the very narrow layer where the runaways are localized. Further out the heat deposited (green dashed) is being removed by radiation (red solid). We note that on a longer time scale the lower-heat-diffusivity case would also develop into a similar state.

Since no MHD stability is being monitored in these simulations, such a runaway – or ohmic – current layer can develop enormous current densities, so that the total current it carries becomes comparable with the current in the rest of the plasma. Once that happens, or perhaps even earlier, the current sheet would become MHD unstable, which would again increase all transport channels, as well as would lead to a flattening of the current profile. It is plausible that this process can be modelled by an inward expansion of the region with finite $\varLambda _{m}$. We note that in simulations with prescribed inward-propagating $\varLambda _{m}$ (not shown here), we observe that the skin current layer does not disappear; instead it moves inward along with the boundary between the intact and chaotic field line regions. Such behaviour is indeed observed in 3-D MHD simulations, e.g. in figure 9 of Bandaru et al. (2021). This inward expansion of the chaotic region may proceed until the field becomes chaotic down to the magnetic axis, or the reformation of flux surfaces may stop this process earlier, in which case the skin current layer survives.⁸

3.4. Intact edge; reverse skin current and runaways

After studying the dynamics in the presence of an intact core we now consider the possibility of the edge region remaining intact while the core undergoes a current profile relaxation. This situation may be representative of an internal plasma instability, which can also arise even in non-disruptive plasmas, such as during sawteeth activity. In addition, it can also be relevant for scenarios exhibiting ‘inside-out thermal quench’, which can occur in case of a shell pellet injection (Hollmann et al. 2019), but may also happen in shattered pellet injection cases, depending on the detailed dynamics of impurity deposition and radiative collapse.

The simulation set-up is similar to that of the intact core case considered in § 3.3, but now the spatial variation of transport coefficients, according to (3.1), uses a transition radius of $r_m=1.9\,{\rm m}$, $X_{\rm edge}=0$ and $X_{\rm core}$ corresponding to $\delta B/B = 3.5\times 10^{-3}$ for runaway particle transport, and a hyperdiffusivity of $\varLambda _{m} = 3\times 10^{-2} \,{\rm Wb}^2\,{\rm m}\,{\rm s}^{-1}$. These transport coefficients are active for $t< 6\,{\rm ms}$, along with a radially constant electron heat diffusion at the same magnetic perturbation amplitude; the latter is reduced to $\delta B/B = 4\times 10^{-4}$ afterwards.

As the current density flattens in the chaotic field region, as seen in figure 4(b), a strong reverse skin current is generated in the intact edge region. Since the temperature is moderate (${\sim }1\,{\rm keV}$) in this edge region already in the beginning and it drops below ${\sim }50\,{\rm eV}$ within $2 \,{\rm ms}$, the resistive diffusion in this region is significant. This leads to a rapid diffusive decay of the reverse skin current; the current density in the skin layer changes direction already at $t=3.5\,{\rm ms}$. Its evolution during the current relaxation is shown in figure 4(c). The peak electric field magnitude in the skin layer ($0.26\,\% E_D$) exceeds that in the rest of the plasma during the course of the entire simulation ($0.13\,\% E_D$); however, it decays too rapidly to lead to a significant reverse runaway current generation. The maximum reverse runaway current density is only $2.5\,{\rm mA}\,{\rm m}^{-2}$. Nevertheless, the reason for the runaway current curve (dash-dotted) in figure 4(a) going outside the plotted logarithmic range in the first few milliseconds is the presence of a negative runaway current in the skin layer (reaching only $-2.4\,{\rm mA}$).

Figure 4. Intact edge case yielding a negligible reverse runaway current. (a) Time evolution of the total plasma current (dashed) and its ohmic (solid) and runaway (dash-dotted) components; enhanced transport is applied at $t<6\,{\rm ms}$, indicated by dotted line. (b) Current density profiles at various times, with ohmic (solid) and runaway (dashed) components. (c) Current density profiles in the skin layer during current relaxation.

Due to the short magnetic diffusion time scale in the edge, changes in the poloidal magnetic flux during the reconnection are not trapped inside the plasma too long, and an $I_p$ spike is observed already during the time period of the current relaxation, as seen in figure 4(a) (solid line). However, while the $I_p$ spike starts with a clear finite positive ${\rm d}I_p/{\rm d}t$ in the case where there is no intact region in the edge (see e.g. figure 1a), now there is no appreciable change of $I_p$ in the first millisecond, then ${\rm d}I_p/{\rm d}t$ gradually increases before the $I_p$ spike. Such delay of the $I_p$ spike onset, which is often observed experimentally, along with the negative skin current, has already been discussed by Wesson et al. (1990), although without the possibility to account for runaway electrons.

How high the values that the reverse electric field reaches in the skin layer, and the time scale it lasts, determine whether a significant reverse runaway beam could form. These factors, in turn, depend strongly on the resistive diffusion time, and thus ultimately on the temperature of the skin layer. In the presence of internal instabilities which are localized deeper in the core, thus allowing for a higher temperature, experiments have shown the generation of superthermal electrons in measurable quantities (Savrukhin 2001; Klimanov et al. 2007; Kamleitner et al. 2015; Mai et al. 2021). To illustrate that it is conceivable to get significant reverse runaway beams in disruptions with an intact edge region, we return again to the alternative scenario detailed in Appendix A, but now, unlike in the case studied in figure 2, with an intact edge region outside $r_m=1.9\,{\rm m}$. In this case a fully kinetic simulation is performed, which is practically necessary in the case where a major reverse runaway population develops during the simulation.

In this scenario $-E/E_D$ reaches $8.1\,\%$ compared to $0.2\,\%$ in the baseline case. This difference translates to a dramatic disparity for the runaway current densities reached; the maximum negative runaway current density reached in the simulation is $36\,{\rm MA}\,{\rm m}^{-2}$, and the highest negative value of the total runaway current is $3.2\,{\rm MA}$. The runaway current density is localized to an $\sim 1\,{\rm cm}$ thin layer, as seen in figure 5(b). At the location of the highest negative runaway current density, $r=1.944\,{\rm m}$, the electric field stays negative in the first $4.7\,{\rm ms}$, and then it changes to the forward direction⁹ (see the orange dotted curve in figure 5c).

Figure 5. Intact edge case in an alternative scenario yielding a sizeable reverse runaway beam (kinetic simulation). (a) Radial profiles of the electric field normalized to the local Dreicer field. (b) Radial profiles of the current density, with ohmic (solid) and runaway (dashed) components. (c) Time evolution of the parallel electric field (orange dotted line) and current density components (solid lines) taken at $r=1.944\,{\rm m}$. Blue $j_{\varOmega }$, ohmic; black $j_{\rm RE}$, total runaway; red $j_{\rm RE}^{-}$, reverse runaway; green $j_{\rm RE}^{+}$, forward runaway.

The electric field first draws a reverse ohmic skin current (blue line in figure 5c), which starts to diffusively decay already after $0.5\,{\rm ms}$. A macroscopic reverse runaway current density is generated slightly after $1\,{\rm ms}$, which starts to dominate the local current density already before $2\,{\rm ms}$, and keeps growing as long as the electric field is negative (see red solid curve, overlaid here by the black curve, the total runaway current density). Soon after the electric field becomes positive, a positively directed runaway current component arises (green). Note that the maximum positive electric field at this location reaches only $0.2\,\%$, which would not be sufficient to generate a macroscopic runaway current under such a short time. However, energetic superthermal electrons which are present due to the reverse runaway population can drift over to the positive direction, and turn into an avalanching forward runaway population, without first slowing back to the bulk. We note that the avalanche source term we employ takes into account the momentum direction of the runaways that generate new ones by close collisions. For a while the backward and forward runaway populations coexist in the remnant of the skin current region.

As mentioned previously, there are observations of energetic electrons in the presence of internal instabilities (Savrukhin 2001; Klimanov et al. 2007; Kamleitner et al. 2015; Mai et al. 2021), such as large sawteeth. However, we are not aware of any clear observation of a significant reverse runaway electron beam. In this literature the propagation direction of the most energetic electrons is either ambiguous, or no directional preference was observed (Kamleitner et al. 2015). This may be explained by the relativistic beaming of the synchrotron or bremsstrahlung emitted by the electrons – that could provide direct evidence for the propagation direction – not yet being strong at these subrelativistic energies. Nevertheless, the energetic electrons are radially localized in the vicinity of the inversion radius, as seen for instance from elevated electron cyclotron emission from this region in Klimanov et al. (2007); this is consistent with the picture of the energetic electrons being generated in the skin layer. In order to guide experimental observations, further high-fidelity numerical analysis will be necessary to explore the parameter regions where significant reverse runaway electron beams are to be expected to develop.

3.5. Outlook

The analysis presented here is not fully self-consistent as the transport coefficients are prescribed. As we have shown, there is a fair degree of robustness to changes in transport coefficients, in particular to the transport of runaways during the thermal quench, or the remnant electron heat transport after the thermal quench. However, particularly in the presence of skin current regions, it would be desirable to capture the dynamic nature of current- and pressure-driven instabilities. A fully self-consistent treatment would require nonlinear MHD simulations with runaway dynamics included. While such numerical models with fluid runaways have become available recently (Bandaru et al. 2019; Liu et al. 2020) they have their own limitations and represent an incomparably larger computational expense than the approach pursued here. There are various ways to move towards self-consistency, ranging from implementing an automatic enhancement of transport once some specified threshold – in for instance pressure gradient or the tearing parameter $\Delta '$ – is exceeded, to utilizing more complex instability criteria and prescriptions based on wisdom gained from nonlinear MHD simulations.

There is a large degree of freedom in prescribing the spatial and temporal structure of helicity transport and the cases shown here only represent a very limited number of possible scenarios. The approach can be informed by close comparison and even fitting to nonlinear MHD simulations. The machinery to turn simulated magnetic perturbations into transport coefficients and using these in Dream simulations is in place (Tinguely et al. 2021). A natural next step is to extend this workflow to helicity transport. The information gained can be used for exploring larger parameter regions with Dream, which can provide interesting scenarios to be analysed by higher-fidelity models.

4. Conclusions

We have analysed the runaway electron dynamics in disruptions of ITER-like plasmas in the presence of current relaxation that extends over the entire plasma, or only to limited radial domains. We use a one-dimensional helicity transport model to capture the current relaxation due to a fast magnetic reconnection, and employ the disruption runaway electron modelling framework Dream. During the prescribed reconnection event, besides the magnetic helicity transport, heat and energetic electrons also undergo spatial diffusion, and heat losses are also enhanced by an injected deuterium–neon mixture.

We find that the current relaxation event reduces the efficiency of runaway generation in the core and increases it towards the edge, thus shifting the centre of runaway generation towards the edge. This may lead to disparate outcomes depending on the scenario, including the possibility of reduced core-localized runaway generation, or edge-localized hollow runaway profiles with larger total runaway current.

In scenarios where the magnetic field does not become chaotic all the way to the magnetic axis, the strong skin current developing at the boundary of the intact region becomes a dynamically important location. If the remnant electron heat diffusivity after the thermal quench is small, such regions may develop into a hot ohmic current channel, with an interesting internal structure. Not accounting for the development of instabilities, the long-time behaviour of these structures is not well represented within the current analysis. Allowing for larger heat diffusivity – that, in scenarios without skin currents, would not have a significant effect – inhibits the development of hot ohmic channels, and instead the skin current region turns into a hotspot for runaway current generation. Revisiting such situations by the micro- and macrostability of the plasma monitored and acted upon appears to be a fruitful path forward.

Finally, in intact edge scenarios, where magnetic flux surfaces do not break up all the way to the separatrix, a reverse skin current layer develops. Intact edge scenarios can be relevant in the case of internal MHD instabilities, in an ‘inside-out thermal quench’, or if helicity transport is inhibited or strongly reduced towards the edge for some other reason. The skin current flowing in the countercurrent direction affects the time evolution of the $I_p$ spike, which is reliably measured experimentally; indeed the $I_p$ spike is often subsequent to the temperature collapse in experiments, consistently with the existence of a layer at the plasma edge with reduced helicity transport. In addition the reverse skin current can act as a hotspot for runaway electron evolution, but now in the countercurrent direction. Since field diffusion tends to be faster at the less hot plasma edge, the conditions for runaway generation tend to decay more quickly compared with a skin current region in the core. However, kinetic simulation results indicate that it may be possible to develop a non-negligible reverse runaway electron current. Even if a macroscopic reverse electron beam would not have time to develop, once the electric field becomes positive, the reverse runaway seed in the skin region can be accelerated into the forward direction and provide a seed for runaway generation even at electric field strengths that would otherwise be too weak for a significant seed generation. The dynamics of reverse skin current regions is worth further investigation, and can be connected to experiments of superthermal electron observations in the presence of large sawtooth crashes.