Runaway dynamics in disruptions with current relaxation

The safe operation of tokamak reactors requires a reliable modeling 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 modeling 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 center concerning the subsequent dynamics: It may turn into a hot ohmic current channel, or a sizable radially localized runaway beam, depending on the heat transport. If the intact region is in the plasma edge, runaway generation in the counter-current direction can occur, which may develop into a sizable 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 center of runaway generation towards the edge.


Introduction
Recent progress in fusion science has increased confidence in that the energy confinement requirements of up-coming experiments aiming to demonstrate a scientific breakeven 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 the 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 et al. 2017) which can, in a collisional avalanche process, exponentiate into a macroscopic runaway current (Sokolov 1979;Rosenbluth & Putvinski 1997;Hesslow et al. 2019a). Owing to the exponential sensitivity of the avalanche to the initial plasma current, high-current reactor relevant devices are predicted to be prone to convert a large fraction of their ohmic current into † Email address for correspondence: pusztai@chalmers.se arXiv:2206.00904v1 [physics.plasm-ph] 2 Jun 2022 runaway electron current (Boozer 2012;Breizman et al. 2019;Vallhagen et al. 2020a), thereby posing a major threat to the structural integrity of the plasma facing components. Methods to mitigate the effects of a disruption  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 a rapid massive material injection, such as through shattered pellet injection (Commaux et al. 2016;Jachmich et al. 2021;Vallhagen et al. 2022), or magnetic perturbations (Tinguely et al. 2021), 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 (3D) magnetohydrodynamic (MHD) simulation codes (Hoelzl et al. 2021;Sovinec et al. 2004;Ferraro et al. 2016), flux surface averaged transport solvers (Linder et al. 2020), local kinetic models (Stahl et al. 2016), 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 3D 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 timescale, 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 sizable temporary increase of the total plasma current I p (Wesson et al. 1990;de Vries et al. 2015), as well as the reduction of the internal inductance of the plasma i = W p /(µ 0 R 0 I 2 p /4), such that the magnetic energy in the poloidal magnetic field W p is only slightly reduced. Here µ 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 timescale 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 timescales. Magnetic helicity H M is the volume integral of A · B with the magnetic vector potential and field, A and B, which can be expressed in tokamak geometry as H M = −2 ψ dψ t , with the poloidal and toroidal magnetic fluxes, ψ and ψ t , defined as in Fig. 2 of (Boozer 2017), and the entire plasma being the integration domain.
The current density evolution during a thermal quench is too rapid to be experimentally directly observable; in practice the total plasma current and-in case of shaped plasmasthe 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 (Hoppe et al. 2021b;Insulander Björk et al. 2021) as their initial condition, which is then evolved in an inductive-resistive fashion. However, the rearrangement of the current density in the fast magnetic reconnection event, and in particular the development of skin currents 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), that describes the current redistribution as a hyper-diffusion 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 et al. 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 center 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 hot-spots. As a possible, somewhat exotic, outcome, "reverse" (counter-current) runaway beams may develop in intact edge regions.
The rest of the paper is organized as follows. First, in Sec. 2.1, we provide the equations relevant for the current evolution in the numerical model, then in Sec. 2.2 we describe the baseline simulation setup. We first consider full-radius current relaxations; the possible outcomes-core localized and edge localized runaway generation-are discussed in Secs. 3.1 and 3.2, respectively. The role of heat transport concerning the fate of a core localized skin current region is discussed in Sec. 3.3. Finally, we consider the possibility of an intact edge, and associated reverse runaway generation in Sec. 3.4. Limitations and future directions are touched upon in Sec. 3.5, and the conclusions are summarized in Sec. 4.

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 ξ = 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 will 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 Ω ; and a runaway electron population, characterized by the current density they carry j re . In the fluid model it is assumed that the runaways, which have the density n re , move with the speed of light c parallel to the magnetic field, hence j re = en re c, where e is the elementary charge.
The total current density is computed through evolving the poloidal flux (2.1) in this equation the toroidal flux ψ t (r) = (2π) −1 r 0 V B · ∇ϕ dr is used as a radial coordinate, with · denoting flux surface average, V is the spatial Jacobian (incremental volume enclosed by two infinitesimally close flux surfaces), and R = 2π E · B / B · ∇ϕ , with the toroidal angle ϕ. This flux function term is equal to the loop voltage when the second term in the right hand side of (2.1) vanishes, and it describes ohmic dissipation, as we may express the term E · B = j Ω B 2 /(σB), with σ the electric conductivity, and j Ω the ohmic current density †. The total current density j tot includes the ohmic and the runaway current densities, j Ω and j re , and in some kinetic simulations j 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 where R is the major radius of a given point. From the form of (2.2) we can see that the term involving Λ m in Eq. (2.1) includes a 4 th order radial derivative of ψ, and as such it describes a hyper-diffusion of ψ wherever Λ 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 Eq. (2.1) is constructed to respect this conservation property is most easily seen by considering d t H M = −2 dψ t ∂ t ψ. Clearly, the contribution of the second term on the r.h.s. of Eq. (2.1) to ∂ t ψ-which is of the form ∂ ψt (. . . )-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 1. In the simulations, the helicity transport coefficient Λ m can have an arbitrary prescribed spatio-temporal structure. Currents in passive structures are taken into account through the boundary condition on Eq. (2.1). The relevant equations are the following where a and b are the minor radii corresponding to the plasma edge and the wall. I(r) is the total plasma current within a given radius r. The edge-wall mutual inductance is M ew = (2π) 2 µ 0 b a dr(V |∇r| 2 /R 2 ) −1 , which in the cylindrical limit reduces to µ 0 R 0 ln(b/a). The external inductance is L ext = µ 0 R 0 ln(R 0 /b), with R = R 0 at the † In kinetic simulations jΩ is replaced by jΩ − δjcorr, with a correction that plays a role in the presence of rapid time variation, as explained at Eq. (29) of (Hoppe et al. 2021a). magnetic axis. The inputs to the wall model are the resistive timescale of the wall τ w = L ext /R wall , and b.

Simulation setup
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 m −3 , the temperature is parabolic with 20 keV on-axis, and the current density corresponds to a total plasma current of 15 MA. The simulations consider an ITER-like magnetic geometry with R 0 = 6 m, a = 2 m, b = 2.15 m, B(r = 0) = 5.3 T, and a resistive wall time of τ w = 0.5 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 breakup of flux surfaces. At t = 6 ms, when the maximum plasma temperature has dropped to ≈ 100 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 by 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;Vallhagen et al. 2020a), the hot-tail seed is calculated using the model in Sec. 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 Eq. (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 re E c , with E c = e 3 n e ln Λ c /(4π 0 m e c 2 ) the critical electric field, 0 the vacuum permittivity, and m e the electron mass. In the definition of E c we take ln Λ c , the momentum dependent Coulomb logarithm, defined in Eq. (18) of (Hoppe et al. 2021a), at the momentum p = 20 m e c; in addition ln Λ 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 where v te = 2T e /m e is the local electron thermal speed, and δB/B is the relative magnetic perturbation amplitude †. This 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 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 δB/B is applied after the thermal quench, throughout the simulation, while physically electron heat transport might then be caused by electrostatic turbulence.
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 RE = πR 0 c(δB/B) 2 . We employ here the same δ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 ∆t = 6.6 µs.

Full-radius current relaxation; core localized runaways
First we briefly discuss the temperature evolution through our baseline case; in most cases a similar qualitative behavior is observed, with the exception of cases with strong ohmic skin currents, to be discussed in Sec. 3.3. After the instantaneous and homogeneous deposition of neutral deuterium and neon of densities n D,inj = 7 · 10 20 m −3 and n Ne,inj = 1.5·10 19 m −3 , the plasma temperature drops on a timescale of 0.01-0.1 µs as these species get ionized and the plasma diluted. The peak temperature drops from 20 to ∼ 2 keV. Then the electron heat transport corresponding to a magnetic perturbation amplitude of δB/B = 3.5·10 −3 cools the plasma further on a ms time scale †, until T e reaches ∼ 100 eV, at which point radiative energy losses become dominant and cool the plasma further to the 1-10 eV range very rapidly. The injected densities are chosen to produce a disruption in the parameter region deemed favorable in (Vallhagen et al. 2020b), 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 δB/B = 4 · 10 −4 (at 6 ms, indicated by the dashed vertical line in Fig. 1a). Then the temperature evolves on the ∼ 10 ms timescale 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.
Without current relaxation (thin light lines in Fig. 1), the total ohmic current (solid line in Fig. 1a) decreases monotonically. After a short period of hot-tail runaway generation that is concluded at 5 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 ∼ MA level through avalanche, the ohmic current drops on the time scale of the avalanche growth, and the runaway current saturates around 5.3 MA. Then the RE current keeps increasing only very slowly as poloidal magnetic filed is diffusing back into the vacuum chamber through the resistive wall.
In case a spatially homogeneous hyperdiffusivity of Λ = 3 · 10 −2 Wb 2 m/s is applied in the first 6 ms (solid dark lines in Fig. 1), when the flux surfaces are assumed to be † It is in fact this characteristic cooling timescale that was targeted by choosing this specific value of δB/B. broken up, we find that the ohmic current first exhibits a peak of 18.5 MA at 3 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 Fig. 1c (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 Fig. 1b, showing the electric field normalized to the Dreicer field E D = (e 3 n e ln Λ c )/(2π 2 0 m e v 2 te ). 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 Fig. 1c) are everywhere lower than the initial ohmic current density. In the case with current relaxation the same observation can be done, 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 MA to 3.9 MA).
In the first 6 ms of the case with current relaxation runaway electron transport is also active, corresponding to δB/B = 3.5 · 10 −3 (consistently with the electron heat transport). This reduces the runaway electron seed, as seen in Fig. 1a. 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.

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 setup is such that it favors 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 behavior 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 Fig 2a to Fig 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 Fig 2b, 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 Λ m = 1.5 · 10 −2 Wb 2 m/s is applied in the first 10 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 Fig 1. We may conclude then, that a current relaxation that extends over the entire plasmareducing 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 favorable in the core, and increase it if the plasma is more prone to develop edge-localized runaway currents. Panels c and 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 ms.

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 3D 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 favor island overlap.
To model such a scenario, we assume that the magnetic surfaces are broken up in the first 6 ms, when we employ a runaway and an electron heat transport corresponding to δB/B = 3.5 · 10 −3 , as well as a hyperdiffusivity of Λ = 3 · 10 −2 Wb 2 m/s, 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 RE and Λ, now denoted by X, jump between their core and edge values X core and X edge according to the following radial variation with the error function Erf, and here X core = 0, X edge corresponding to the values given above, a transition radius r m = 0.3 m and a characteristic transition width of w = 0.01 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-TQ) heat transport has a major role in the dynamics. First we consider a lower remnant electron heat transport at δB/B = 4 · 10 −4 . Similarly to Fig. 1c and 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 Fig. 3a. 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 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 ms, reaching a maximum runaway rate of 2.7 · 10 17 m −3 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.
The ohmic current channel in the intact region broadens inward in time, and it locally heats the plasma above 100 eV, as seen in Fig. 3b. Once the runaway current sheet increased 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 Fig. 3c, taken at t = 60 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 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 the extremely large temperature gradients at its edges are likely to be unstable to microinstabilities, 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 being 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 Fig. 3a-c, but corresponding to a remnant electron heat diffusivity at δB/B = 2 · 10 −3 (compared to δB/B = 4 · 10 −4 ) are shown in Fig. 3d-f. This time the ohmic current quench completes within 60 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, Fig. 3d and 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 MA before 20 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 to the lower heat diffusivity case; compare Fig. 3c and f. As the ohmic current have decayed to small values, the contribution of ohmic heating (black dotted) is negligible compared to 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 modeled by an inward expansion of the region with finite Λ m . We note that in simulations with prescribed inward propagating Λ 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 behavior is indeed observed in 3D MHD simulations, e.g. in Fig. 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 †.

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 saw-teeth 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 , but may also happen in shattered pellet injection cases, depending on the detailed dynamics of impurity deposition and the radiative collapse.
The simulation setup is similar to the intact core case considered in Sec. 3.3, but now the spatial variation of transport coefficients, according to Eq. (3.1), uses a transition radius of r m = 1.9 m, X edge = 0, and X core corresponding to δB/B = 3.5 · 10 −3 for runaway particle transport, and a hyperdiffusivity of Λ = 3 · 10 −2 Wb 2 m/s. These transport coefficients are active for t < 6 ms, along with a radially constant electron heat diffusion at the same magnetic perturbation amplitude; the latter is reduced to δB/B = 4 · 10 −4 afterwards.
As the current density flattens in the chaotic field region, as seen in Fig. 4b, a strong reverse skin current is generated in the intact edge region. Since the temperature is moderate (∼ 1 keV) in this edge region already in the beginning and it drops below ∼ 50 eV within 2 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 ms. Its evolution during the current relaxation is shown in Fig. 4c. 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 mA/m 2 . Nevertheless, the reason for the runaway current curve (dash-dotted) in Fig. 4a to go outside the plotted logarithmic   Figure 4. Intact edge case yielding a negligible reverse runaway current. Panel 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 ms, indicated by dotted line. Panel b: Current density profiles at various times, with ohmic (solid) and runaway (dashed) components. Panel c: Current density profiles in the skin layer during current relaxation.   Figure 5. Intact edge case in an alternative scenario yielding a sizable reverse runaway beam (kinetic simulation). Panel a: Radial profiles of the electric field normalized to the local Dreicer field. Panel b: Radial profiles of the current density, with ohmic (solid) and runaway (dashed) components. Panel c: Time evolution of the parallel electric field (orange dotted line) and current density components (solid lines) taken at r = 1.944 m. Blue jΩ -ohmic, black jRE -total runaway, red j − RE -reverse runaway, green j + RE -forward runaway.
range in the first few milliseconds is the presence of a negative runaway current in the skin layer (reaching only −2.4 mA).
Due to the short magnetic diffusion timescale in the edge, changes in the poloidal magnetic flux during the reconnection is 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 Fig. 4a (solid line). However, while the I p spike starts with a clear finite positive dI p /dt in case there is no intact region in the edge, see e.g. Fig. 1a, now there is no appreciable change of I p in the first ms, then the dI p /dt gradually increases before the I p spike. Such delay of the I p spike onset, 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 values the reverse electric field reaches in the skin layer, and the timescale 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 (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 Fig. 2, with an intact edge region outside r m = 1.9 m. In this case a fully kinetic simulation is performed, which is practically necessary in case 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 MA/m 2 , and the highest negative value of the total runaway current is 3.2 MA. The runaway current density is localized to a ∼ 1 cm thin layer, as seen in Fig. 5b. At the location of the highest negative runaway current density, r = 1.944 m, the electric field stays negative in the first 4.7 ms, and then it changes to the forward direction †, see the orange dotted curve in Fig. 5c.
The electric field first draws a reverse ohmic skin current (blue line in Fig. 5c), which starts to diffusively decay already after 0.5 ms. A macroscopic reverse runaway current density is generated a bit after 1 ms, which starts to dominate the local current density already before 2 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.

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 magnetohydrodynamic 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 ∆ -is exceeded, or 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 analyzed by higher fidelity models. † The sign of the electric field is defined such that a positive electric field drives a positive ohmic current that is directed as the initial plasma current.

Conclusions
We have analyzed 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 modeling 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 center 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 behavior 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 hot-spot for runaway current generation. Revisiting such situations by the micro and macro stability 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 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. Not only does the skin current affect the time evolution of the I p spike, that is reliably measured experimentally, but it can act as a hot-spot for runaway electron evolution, but now in the counter-current 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 to 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.  Figure 6. Panels a-c: Initial plasma profiles in the baseline (solid) and the alternative (dashed) case; a: Total density of hydrogen isotopes present, b: electron temperature, c: current density (the total current in both cases is 15 MA, which corresponds to a lower j in the elongated baseline case). Panel d: Shaping parameters in the baseline case.
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.

Appendix A. Initial plasma and magnetic geometry profiles
The initial plasma parameter profiles are shown in figure 6a-c in the baseline case (solid curves) and the alternative case (dashed). The current density j is taken at the outboard mid-plane, which is the definition used throughout the paper (and in Dream). In the baseline case the plasma is shaped with magnetic geometry parameters shown in figure 6d. The geometric quantity G determines the magnitude of the toroidal magnetic field B t = |G∇ϕ|, with the toroidal angle ϕ; its on-axis value is B 0 = 5.3 T. Furthermore the elongation κ and the triangularity δ are defined as in the Miller model equilibrium (Miller et al. 1998), and the Shafranov shift parameter, ∆ = R(r) − R 0 , is defined here to be zero on-axis. The magnetic equilibrium is not evolved self-consistently in the simulation, thus these shaping parameters are held fixed.
The alternative case also represents an ITER-sized plasma, but the physics content of the simulation is strongly reduced compared to that of the baseline. In this case we use a cylindrical model of the plasma (e.g., all toroidicity effects neglected and there is no plasma shaping) with minor radius a = 2 m and major radius R 0 = 6 m. A perfectly conducting wall is located at the radius b = 1.35a, and the on-axis magnetic field is B 0 = 5.3 T. The plasma composition is prescribed (charge states are not evolved), and it consists of fully ionized deuterium with density, as shown in 6a, and Ar 5+ of density n Ar = 0.2n D . The temperature evolution is prescribed as where the final temperature T f = 3 eV is radially constant, and the characteristic time of the temperature decay is t TQ = 1 ms. This choice of T f is somewhat arbitrary; it is smaller than typical early post-thermal quench temperatures tend to be (5-10 eV), and it favors a large seed generation. This circumstance is not crucial for the behavior we intend to illustrate with this example (for instance, favoring runaway generation towards the edge). Here, the plasma conductivity is calculated using the collisionless limit of the model of (Redl et al. 2021). Only Dreicer and avalanche runaway generation mechanisms are active, using the same models as in the baseline case. Runaway electrons are not transported radially in these simulations.
In the kinetic simulation presented in Fig. 5 we set Λ m = 1.5 · 10 −2 Wb 2 m/s inside of the radius r m = 1.9 m transitioning to zero outside, over a characteristic distance of w = 5 mm. This simulation uses 40 radial grid cells, 25 of which are packed around the skin current region. The momentum space of the thermal and superthermal regions are resolved by 25 cells in ξ = p /p, where p is the magnitude of momentum and p is its component along the magnetic field (taken at the lowest magnetic field point along the orbit). The thermal region extends up to p = 0.07m e c and is resolved by 60 cells in p, and the superthermal region extends to p = 2m e c and is resolved by 80 momentum cells. The runaway region extends up to p = 40m e c and has 50 cells in p and 100 in ξ.
a similar simulation without current relaxation (Λ m = 0). Thus the size of the helicity change we observe is consistent with being caused by the finite resistivity of the plasma. The radial resolution used is 400 grid cells for 0.1 < r/a < 0.95, and 50 cells distributed in the remaining radial domain.