1. Introduction
Disruptions in tokamak plasmas may lead to the formation of a beam of relativistic, so-called runaway electrons (RE), which has the potential to severely damage plasma-facing components (Hender et al. Reference Hender, Wesley, Bialek, Bondeson, Boozer, Buttery, Garofalo, Goodman, Granetz and Gribov2007). For this reason, much effort is directed towards the development of schemes to avoid, limit or mitigate the formation of such a beam. One proposed measure is using massive material injection (MMI) in the form of gas (massive gas injection – MGI) or frozen pellets (shattered pellet injection) to avoid or dissipate the runaway electrons (Hollmann et al. Reference Hollmann, Aleynikov, Fülöp, Humphreys, Izzo, Lehnen, Lukash, Papp, Pautasso, Saint-Laurent and Snipes2015a; Lehnen et al. Reference Lehnen, Aleynikova, Aleynikov, Campbell, Drewelow, Eidietis, Gasparyan, Granetz, Gribov and Hartmann2015), and its efficiency has been demonstrated in medium-sized tokamaks (Hollmann et al. Reference Hollmann, Parks, Commaux, Eidietis, Moyer, Shiraki, Austin, Lasnier, Paz-Soldan and Rudakov2015b; Reux et al. Reference Reux, Plyusnin, Alper, Alves, Bazylev, Belonohy, Boboc, Brezinsek, Coffey and Decker2015; Pautasso et al. Reference Pautasso, Bernert, Dibon, Duval, Dux, Fable, Fuchs, Conway, Giannone and Gude2016; Esposito et al. Reference Esposito, Boncagni, Buratti, Carnevale, Causa, Gospodarczyk, Martin-Solis, Popovic, Agostini and Apruzzese2017; Paz-Soldan et al. Reference Paz-Soldan, Cooper, Aleynikov, Pace, Eidietis, Brennan, Granetz, Hollmann, Liu and Lvovskiy2017; Carnevale et al. Reference Carnevale, Ariola, Artaserse, Bagnato, Bin, Boncagni, Bolzonella, Bombarda, Buratti and Calacci2018; Coda et al. Reference Coda, Agostini, Albanese, Alberti, Alessi, Allan, Allcock, Ambrosino, Anand and Andrèbe2019; Mlynar et al. Reference Mlynar, Ficker, Macusova, Markovic, Naydenkova, Papp, Urban, Vlainic, Vondracek and Weinzettl2019; Pautasso et al. Reference Pautasso, Dibon, Dunne, Dux, Fable, Lang, Linder, Mlynek, Papp and Bernert2020). However, the plasma currents, temperatures and densities of future devices such as ITER will be significantly larger than what can be achieved in current experiments, and simulations are necessary to foresee the effectiveness of massive material injections for disruption mitigation under such conditions.
As MGI is widely used in current devices, cases where runaways are formed in MGI-induced disruptions provide a valuable dataset to properly understand the physics of runaway electron formation and dissipation in such scenarios. To gain this understanding, theoretical models have been formulated and implemented in computational tools which can be used to model disruptions (Breizman & Aleynikov Reference Breizman and Aleynikov2017). To be applicable for predictions for ITER and beyond, theoretical tools must first be validated against existing experimental data to ensure that they capture the relevant physics.
One such kinetic modelling tool is code (COllisional Distribution of Electrons), briefly described in § 2 and in detail in the paper by Stahl et al. (Reference Stahl, Embréus, Papp, Landreman and Fülöp2016). This tool includes modelling of many phenomena important for the studied scenarios, in particular partial screening of partially ionized impurities, and was therefore chosen for the present investigations. The present paper discusses the comparison of code simulations with experimental data obtained from the ASDEX Upgrade tokamak (AUG, § 3.1). Runaway-generating disruptions are deliberately triggered in AUG to obtain relevant data for studies of the connected phenomena, and parameters important for modelling are measured by plasma diagnostics, as described in §§ 3.2 and 3.3.
In a disruption, the induced electric field accelerates electrons above a certain critical velocity to relativistic energies. The Dreicer runaway generation is the result of velocity space diffusion of the electrons into the runaway region due to small-angle collisions (Dreicer Reference Dreicer1959). Existing runaways can create new fast electrons through close collisions with bulk electrons (secondary generation). The latter leads to an exponential growth of the number of REs – an avalanche (Sokolov Reference Sokolov1979; Rosenbluth & Putvinski Reference Rosenbluth and Putvinski1997).
In the case of sudden cooling (a thermal quench, TQ), fast electrons do not have time to thermalize, and a so-called hot-tail forms in the electron distribution. Hot-tail generation is the dominant primary generation when the TQ duration is shorter than the collision time at the runaway threshold velocity (Helander et al. Reference Helander, Smith, Fülöp and Eriksson2004). In future fusion devices, due to the higher initial plasma temperature, the collision time is not much shorter than the expected duration of the TQ, and therefore a sizable hot-tail RE seed is likely to be produced (Smith et al. Reference Smith, Helander, Eriksson and Fülöp2005). A simplified analytical model for hot-tail generation has been formulated by Smith & Verwichte (Reference Smith and Verwichte2008) but comparisons with kinetic simulations show that this model underestimates the runaway density by an order of magnitude (Stahl et al. Reference Stahl, Embréus, Papp, Landreman and Fülöp2016), so to gain a quantitative understanding of hot-tail generation, kinetic simulations are needed.
The future scenario we desire to understand is a spontaneous disruption, mitigated by MMI. However, since spontaneous disruptions do not occur reproducibly, we instead model disruptions which were deliberately triggered by an MGI, resulting in a scenario which is similar to the desired scenario in the important aspect that a runaway current is formed and dissipated in the presence of partly ionized high-Z materials.
Such scenarios have been considered also in the recent paper by Linder et al. (Reference Linder, Fable, Jenko, Papp and Pautasso2020), where the astra-strahl 1.5-dimensional (1.5-D) transport code (Dux et al. Reference Dux, Peeters, Gude, Kallenbach and Neu1999; Fable et al. Reference Fable, Angioni, Ivanov, Lackner, Maj, Medvedev, Pautasso, Pereverzev and Treutterer2013) was used, including reduced kinetic models for Dreicer and avalanche runaway generation (Hesslow et al. Reference Hesslow, Embréus, Vallhagen and Fülöp2019a,Reference Hesslow, Unnerfelt, Vallhagen, Embréus, Hoppe, Papp and Fülöpb). The modelling results presented here intend to assess the kinetic part of the process, namely to which extent code captures the central aspects of runaway current formation and dissipation and what role hot-tail generation plays. For a full understanding of the disruption scenario, kinetic simulation must be combined with other tools, most importantly modelling spatial dynamics, 3-D MHD (magnetohydrodynamic) evolution and the atomic physics needed to determine ionization states.
2. Kinetic modelling
The relativistic finite-difference Fokker–Planck solver code (Stahl et al. Reference Stahl, Embréus, Papp, Landreman and Fülöp2016) simulates electron dynamics in plasmas in 2-D momentum space. The plasma is assumed to be homogenous, i.e. it is a 0-D simulation in real space and radial transport or instabilities are not modelled. We focus on modelling the evolution of the electron momentum distribution function on axis, and will use experimental data representative of the conditions on the plasma axis to the extent possible.
code calculates the time-evolving electron distribution function under the influence of collisions, synchrotron radiation reactions and electric field acceleration. In the simulations presented here, collisions are modelled by a relativistic test-particle operator (Braams & Karney Reference Braams and Karney1989, detailed in Hesslow et al. Reference Hesslow, Unnerfelt, Vallhagen, Embréus, Hoppe, Papp and Fülöp2019b) and a simplified large-angle collision operator (Rosenbluth & Putvinski Reference Rosenbluth and Putvinski1997). Screening of partially ionized impurities is taken into account according to the model described by Hesslow et al. (Reference Hesslow, Embréus, Hoppe, DuBois, Papp, Rahm and Fülöp2018a). Bremsstrahlung radiation losses were found to be negligible in the scenarios considered here, as well as the differences between the fully conservative and the simplified avalanche operator (Embréus, Stahl & Fülöp Reference Embréus, Stahl and Fülöp2018).
The electric field 
$E$ and the plasma current density 
$j$ are calculated self-consistently throughout the simulation, according to (Hesslow et al. Reference Hesslow, Embréus, Wilkie, Papp and Fülöp2018b)
\begin{equation} E =-\frac{L a^2}{2 R}\frac{\textrm{d} j}{\textrm{d} t}, \end{equation}
where the inductance 
$L$ is given by
\begin{equation} L \approx \mu_0 R \left[ \ln \left(\frac{8 R}{a} \right)-2 +\frac{l_i}{2}\right]. \end{equation}
The major and minor radii, 
$R$ and 
$a$ respectively, are given in table 1 and 
$l_i$ is the internal inductance. The value of 
$l_i$ differs between discharges, but has a negligible effect on the simulation results. The difference in calculated final runaway current between the cases using 
$l_i$ = 0 and 
$l_i$ = 1.5 (which is a common estimate for 
$l_i$ in AUG) was less than 1 % for AUG discharge #35408. This is expected, partly since the inductance differs by only 3 % between the two cases, and partly since, as will be shown later, the primary RE generation is dominated by the hot-tail mechanism which is not sensitive to the induced electric field. A more accurate expression for the inductance is 
$L = \lvert \Psi _p/I_p \rvert$ (Boozer Reference Boozer2018), where 
$\Psi _p$ is the poloidal magnetic flux and 
$I_p$ is the plasma current. To calculate the on-axis magnetic flux 
$\Psi _p$ requires the plasma current density profile, which we do not know with an accuracy that motivates the use of the more accurate expression.
Common parameters of all modelled discharges. The initial current is slightly lower, 0.71 MA, in discharge #31318.
Using the test-particle collision operator is computationally efficient but leads to the underestimation of the ohmic current by about a factor of two (Helander & Sigmar Reference Helander and Sigmar2005) and in a self-consistent calculation, this needs to be compensated for. As the conductivity obtained by the test-particle operator 
$\sigma _{\textrm {CODE},\textrm {tp}}$ is proportional to the fully relativistic electric conductivity 
$\sigma _{\textrm {BK}}$ (Braams & Karney Reference Braams and Karney1989) for a wide range of effective charges and temperatures (Hoppe et al. Reference Hoppe, Hesslow, Embreus, Unnerfelt, Papp, Pusztai, Fülöp, Lexell, Lunt and Macusova2020), we can model the plasma current density as
\begin{equation} j = j_{\textrm{CODE},\textrm{tp}} + \Delta \sigma E, \end{equation}
where 
$j_{\textrm {CODE},\textrm {tp}}$ is the current density calculated by code using the test particle operator, and 
$\Delta \sigma (Z_{\textrm {eff}},T_e)=\sigma _{\textrm {BK}}-\sigma _{\textrm {CODE},\textrm {tp}}$. The calculation of the effective ionic charge 
$Z_{\textrm {eff}}$ is described in § 3.2.3 and the free electron temperature 
$T_e$ in § 3.3.
3. Experimental data
3.1. ASDEX Upgrade
In this paper we model ASDEX Upgrade discharges specifically designed for the study of runaway electron dynamics (Pautasso et al. Reference Pautasso, Bernert, Dibon, Duval, Dux, Fable, Fuchs, Conway, Giannone and Gude2016). ASDEX Upgrade is a medium-sized tokamak at the Max Planck Institute for Plasma Physics in Garching, Germany. The typical runaway electron scenario uses a low density (initial free electron density 
$n_{e0}\approx 3 \times 10^{19}\ \mathrm {m}^{-3}$), inner wall limited, circular (elongation 
$\kappa \approx 1.1$), ECRH (electron cyclotron resonance heating) heated plasma terminated with an argon MGI. The argon is held in a 0.1 l chamber, at room temperature, before the injection valve is triggered at 1.000 s into the discharge. In the considered discharges, the pressure in the valve ranged from 0.11 to 0.9 bar. During the simulated time span, the plasma position remained stable radially as well as vertically.
Eleven discharges with similar plasma parameters but with different amounts of argon injected were selected for modelling. As the discharges – which represent a scan of injected argon quantity – were selected from a database spanning multiple years, the electron temperature 
$T_{e0}$ before the disruption varies in the dataset (
$T_{e0} = 8.0 \pm 2.9\ \textrm {keV}$), because of varying experimental conditions and occasional ECRH gyrotron trips over the years. In ten of the modelled discharges, a runaway current was formed, and discharge #35400, in which no runaway current was formed, was added for comparison. Parameters that are common to all modelled discharges are shown in table 1 and an overview of the basic parameters for the eleven selected discharges is presented in table 2.
Basic on-axis parameters of the eleven simulated discharges and the notation used in this paper. The injection pressure is expressed in bars in the 0.1 l injection volume. The injected number of Ar atoms is estimated assuming a gas temperature of 300 K. The initial free electron density is the value given by 
$\textrm {CO}_2$ interferometry (average over the first 1.5 ms after the argon valve trigger, i.e. before Ar penetration into the plasma) and the initial electron temperature by the electron cyclotron emission measurements.
The initial on-axis current density is used as input to code when creating an initial electron momentum distribution function. The on-axis current density 
$j$ can be estimated using current density profiles obtained from equilibrium reconstructions by CLISTE (CompLete Interpretive Suite for Tokamak Equilibria) (McCarthy Reference McCarthy1999). For the modelled discharges, the initial on-axis current density 
$j_0$ was estimated to approximately 
$1.2\ \textrm {MA}\,\textrm {m}^{-2}$. The initial plasma current 
$I_{p0}$ is very similar between all discharges, and equals 0.76 MA as listed in table 1. The conversion factor between the estimated 
$j_0$ and the measured 
$I_{p0}$ is thus 
$0.76/1.2 = 0.63\ \textrm {m}^{2}$. Since this conversion factor has the unit of 
$\textrm {m}^{2}$, it will be referred to as an ‘effective area’, 
$A_{\textrm {eff}}$. Application of this conversion factor results in an initial on-axis current density of 
$1.21\ \textrm {MA}\,\textrm {m}^{-2} \pm <1\,\%$ for all discharges except the very early discharge #31318 in which the initial on-axis current density is 
$1.13\ \textrm {MA}\,\textrm {m}^{-2}$. This estimated conversion factor was applied in the simulations of all discharges and throughout the simulated time period. Since the pre-disruption plasma is inner wall limited, and the post-disruption runaway beam is surrounded by a low temperature companion plasma, the radial extent of the two stages is comparable.
3.2. Densities
The disruptions were triggered by injection of argon into the plasma (Fable et al. Reference Fable, Pautasso, Lehnen, Dux, Bernert and Mlynek2016; Pautasso et al. Reference Pautasso, Bernert, Dibon, Duval, Dux, Fable, Fuchs, Conway, Giannone and Gude2016). When penetrating into the plasma, the injected argon is partly ionized. The density of free electrons, as well as the density and charge states of the argon atoms, directly affect the collision operator, and hence the evolution of the electron momentum distribution function. Thus, these parameters must be estimated and used as input to the simulations.
3.2.1. Free electron density
The free electron density is measured by 
$\textrm {CO}_2$ interferometry. This method yields the line integrated free electron density along the line of sight of the interferometer. The measured value can be divided by an estimate of the chord length, i.e. the portion of the line of sight that passes through the plasma, and the resulting density value is an average over the estimated chord length. Thus, the density value obtained by this method is not fully representative of the on-axis conditions we intend to simulate.
Previous work by Fable et al. (Reference Fable, Pautasso, Lehnen, Dux, Bernert and Mlynek2016) indicates that the free electron density increases rapidly on the plasma edge, but remains constant on-axis until the MHD mixing event that also causes the plasma current 
$I_p$ to increase for about a millisecond just before starting to decay due to the increased plasma resistivity. The time of the onset and end of this current spike are referred to as 
$t_{\textrm {onset}}$ and 
$t_{\textrm {end}}$ respectively, where 
$t_{\textrm {end}}$ is defined by the peak of the current spike, and 
$t_{\textrm {onset}}$ is the time when the measured 
$I_p$ starts to increase, just before the current spike. We assume that the on-axis free electron density remains constant at the pre-injection value until 
$t_{\textrm {onset}}$, and that after 
$t_{\textrm {end}}$ the plasma is homogenous enough for the measured line-averaged free electron density to be representative of the on-axis value. The data are smoothed using the rloess algorithm in matlab (matlab2017) to avoid numerical difficulties caused by the signal noise. Between 
$t_{\textrm {onset}}$ and 
$t_{\textrm {end}}$, the free electron density is assumed to increase linearly. Simulations were run with different assumptions on the density increase rate, but the calculated current density was insensitive to these variations. The value for the initial on-axis free electron density is taken as the average measured free electron density during the first 1.5 ms after the argon valve trigger, since this is the time period before the measured chord-length-averaged density starts to increase due to the argon injection. The resulting time evolution of the free electron density is shown in figure 1(a), along with the measured free electron density.
$(a)$ The free electron density measured by interferometry (
${\cdots \cdots \cdots }$) and the assumed on-axis values for the densities of argon and free electrons (———) for AUG discharge #35408. The end of the measured 
$I_p$ spike is also marked with a vertical line for reference. 
$(b)$ The dissipation rates for the experimentally measured total plasma current 
$I_p$ and the calculated current densities for three different values of the assimilation fraction 
$f$. Note that the 
$I_p$ dissipation rate has been scaled with 
$A_{{\textrm {eff}}} = 0.63\ \textrm {m}^{2}$.
3.2.2. Argon density
The injected amount of argon is quantified by the pressure 
$p_{\textrm {Ar}}$ in the MGI chamber holding the argon gas before injection. This quantity is listed in table 2, as well as the corresponding number of injected argon atoms 
$N_{\textrm {Ar}}$, assuming a valve volume of 0.1 l and a temperature of 300 K. The average argon density inside the tokamak vacuum chamber after the injection can be calculated, but does not necessarily equal the on-axis argon density. Thus, some assumptions have to be made regarding which fraction of the injected argon assimilates in the plasma (referred to as the assimilation factor, 
$f$), and also the time dependence of the assimilation. 
$f$ is defined as the fraction of the total injected argon which, after the MHD mixing, resides within the plasma region defined by the major and minor radii listed in table 1.
As will be shown later, in § 4.1, the plasma current drops suddenly during the disruption (referred to as the current quench, CQ) and then dissipates more slowly during the so-called plateau phase. As shown in figure 1(b), experimental data show that there is a linear correlation between the injected argon amount and the current dissipation rate in this phase, as long as there is some RE generation but not full conversion of the initial current into RE current. The assimilation factor was estimated by comparison of measured and simulated current density dissipation rates, 
$\textrm {d}j/\textrm {d}t$, during the plateau phase. For comparison with the calculated current density dissipation rate, the measured current dissipation rate has been scaled with the effective area 
$A_{{\textrm {eff}}} = 0.63\ \textrm {m}^{2}$ explained in § 3.1. The current dissipation rate has been calculated, for each discharge except the no-RE discharge #35400, from the measured current, as an average over the time period from 20 to 30 ms after the argon valve trigger. This time period is well into the plateau phase for all simulated discharges.
The current density dissipation rate is calculated from the current density given by the kinetic simulations as the average over the same time span as the current dissipation rate (from 20 to 30 ms after the argon valve trigger). The calculated dissipation rates are plotted in figure 1(b) for 
$f = 10\,\%$, 
$f = 20\,\%$ and 
$f = 40\,\%$, along with a linear fit for each value of 
$f$. As the figure shows, the slope of the linear fit for 
$f = 20\,\%$ approximately reproduces the experimentally deduced slope. Thus, the assimilation factor 
$f = 20\,\%$ was used for all simulated discharges, as a best estimate. The argon densities given by this estimate are in agreement with results available in the literature (Papp et al. Reference Papp, Pautasso, Decker, Hesslow, Bernert, Blanchard, Bock, Bolzonella, Carnevale and Cavedon2019; Pautasso et al. Reference Pautasso, Dibon, Dunne, Dux, Fable, Lang, Linder, Mlynek, Papp and Bernert2020). The positive offset of the experimental current dissipation rates comes from the fact that the runaway plateau current is in controlled ramp down in the experiment, which in one case gives rise to a slight current increase during the early plateau phase. The electric field in the simulations for the runaway electron generation is developing self-consistently, without the comparatively small external electric field, as the inclusion of such field would require a full self-consistent simulation of the AUG automatic control system.
Similarly to the free electron density, the on-axis argon density has been shown to remain approximately constant (zero) until MHD mixing occurs at 
$t_{\textrm {onset}}$ (Fable et al. Reference Fable, Pautasso, Lehnen, Dux, Bernert and Mlynek2016). After 
$t_{\textrm {end}}$, the argon density is assumed to be constant, at a level given by
\begin{equation} n_{\textrm{Ar}} = \frac{N_{\textrm{Ar}}f}{V_{\textrm{plasma}}} = \frac{N_{\textrm{Ar}}f}{2{\rm \pi}^2Ra^2}. \end{equation}
Here, 
$R$ and 
$a$ are the major and minor radii of the plasma given in table 1. The assimilation fraction 
$f = 20\,\%$ was used in the simulations for all modelled discharges. Between 
$t_{\textrm {onset}}$ and 
$t_{\textrm {end}}$, the argon density is assumed to increase linearly. The calculated current was shown neither to be sensitive to the detailed time evolution of the density within the time interval between 
$t_{\textrm {onset}}$ and 
$t_{\textrm {end}}$, nor to changes of time interval length between 
$(t_{\textrm {end}}-t_{\textrm {onset}})/2$ and 
$(t_{\textrm {end}}-t_{\textrm {onset}})\times 2$.
3.2.3. Charge states
For the purpose of the kinetic simulations, the average charge state 
$Z_{\textrm {eff}}$ of the argon is inferred from measurements, i.e. not calculated from atomic physics; 
$Z_{\textrm {eff}}$ is estimated by dividing the density of free electrons attributable to argon by the density of argon atoms. All free electron densities above the initial value 
$n_{e0}$ are attributed to argon. The corresponding distribution over the discrete ionization states is calculated by linear interpolation between the two integers closest to the calculated average charge state, e.g. if the average charge state is calculated to be 4.5, then half of the states are assumed to be 4 and the other half 5. The thus inferred 
$Z_{\textrm {eff}}$ generally shows a distinct peak of around 6 immediately after the end of the current spike and then fluctuates at lower values during the remaining simulation time.
Detailed atomic physics modelling yields a broader distribution over multiple ionization states (Linder et al. Reference Linder, Fable, Jenko, Papp and Pautasso2020). Simulations were also performed using a distribution over multiple charge states given by equilibrium between excitation and recombination rates at the given temperature. It is, however, unlikely that this equilibrium is reached during the rapid thermal quench, and the difference in final calculated current was less than 
$1\,\%$ between the cases run with equilibrium assumption and linear interpolation respectively.
3.3. Thermal quench
The presence of impurities in the plasma abruptly increases radiative energy losses, which leads to the rapid TQ. The simulation results are very sensitive to the duration of the TQ.
3.3.1. Measured free electron temperature
Previous work by Fable et al. (Reference Fable, Pautasso, Lehnen, Dux, Bernert and Mlynek2016) indicates that the on-axis free electron temperature 
$T_e$ remains almost constant until MHD mixing occurs in MGI induced disruptions in ASDEX Upgrade. The free electron temperature is indirectly measured through electron cyclotron emission (ECE), which also yields a temperature profile in the plasma. The thus measured free electron temperature, averaged over a circular 
$15\ \textrm {cm}^{2}$ area surrounding the plasma axis, is shown in figure 3 (
${\cdots \cdots \cdots }$). However, the ECE signal is cut off from approximately 1.5 ms after the beginning of the argon injection until 
$t_{\textrm {end}}$, due to the high free electron density at the plasma edge (Fable et al. Reference Fable, Pautasso, Lehnen, Dux, Bernert and Mlynek2016). Thus, we calculate the initial free electron temperature 
$T_{e0}$ as the measured value averaged over the central 
$15\ \textrm {cm}^2$ area and the time interval between 1 and 1.5 ms after the argon injection valve trigger. This time interval was chosen to exclude both any initial temperature variations due to the heating system being shut off shortly before the disruption, and the beginning of the TQ; 
$T_e$ is then assumed to decrease exponentially with time 
$t$, as argued by Smith & Verwichte (Reference Smith and Verwichte2008). We thus assume that 
$T_e = T_{\exp }$, where 
$T_{\exp }$ is given by (3.2), from the beginning of the simulation until 
$T_{\exp } < T_{\textrm {EQ}}$, where 
$T_{\textrm {EQ}}$ is given by (3.4).
\begin{equation} T_{\exp} = T_{e0} - \frac{T_{e0}-T_{e,\textrm{final}}}{1+\exp\left(\dfrac{t_{\textrm{end}}-t_{\textrm{TQ}}-t}{t_{\textrm{TQ}}}\right)}. \end{equation}
Here, 
$t_{\textrm {TQ}}$ is determined as described in § 3.3.2. The final temperature 
$T_{e,\textrm {final}}$ can be any low value, because after the TQ, the values given by (3.2) are no longer used, but instead the equilibrium temperature is determined using the equilibrium assumption described in § 3.3.3. The final equilibrium temperature is approximately 1 eV for all the simulated discharges. The condition for using the equilibrium temperature is that it is higher than the temperature given by (3.2), so 
$T_{e,\textrm {final}}$ is fixed to 0.5 eV to make sure that the temperature given by (3.2) falls below the equilibrium temperature in all simulations.
3.3.2. Thermal quench time
The thermal quench time is quantified by the parameter 
$t_{\textrm {TQ}}$, whose significance is shown in (3.2). At the fastest phase of the TQ, the temperature drops from 
$62\,\%$ to 
$38\,\%$ of its initial value during the time span 
$t_{\textrm {TQ}}$, as a consequence of the formulation of (3.2). The choice of 
$t_{\textrm {TQ}}$ proves to be very important for the RE generation. For too large 
$t_{\textrm {TQ}}$, there is no RE generation in the simulation, whereas for too small 
$t_{\textrm {TQ}}$, the entire current density is converted to RE current density due to the exponential sensitivity of the hot-tail generation to this parameter (Smith & Verwichte Reference Smith and Verwichte2008). In the experiment, a RE current was formed in all modelled discharges except #35400, but full conversion was not observed in any of the discharges, and thus 
$t_{\textrm {TQ}}$ was chosen to reproduce this result. Using the same value of 
$t_{\textrm {TQ}}$ in each discharge did not yield the desired result – a low 
$t_{\textrm {TQ}}$ resulted in no RE formation in the high-injection cases, whereas a high 
$t_{\textrm {TQ}}$ resulted in full conversion in the low-injection cases. Also, a constant value of 
$t_{\textrm {TQ}}$ would not be expected, since the time scale of the disruption is affected by, among other parameters, the number of injected argon atoms. To quantify this dependence, the time between the argon valve trigger and the end of the 
$I_p$ spike was studied as a function of the injected argon pressure. As shown in figure 2(a), an inverse relationship is found between these parameters indicating that the onset of the CQ is faster for higher injected argon pressures 
$p_{\textrm {Ar}}$. Therefore in the modelling we use the assumption that the thermal quench time is shorter for higher 
$p_{\textrm {Ar}}$. An interpolation (fitting) formula describing the assumed relationship between 
$t_{\textrm {TQ}}$ and 
$p_{\textrm {Ar}}$ was used:
\begin{equation} t_{\textrm{TQ}} = \frac{A}{p_{\textrm{Ar}}-B} + C, \end{equation}
with 
$A=7 \times 10^{-4}\ \textrm {bar}\,\textrm {s}$, 
$B=0.1\ \textrm {bar}$ and 
$C=0.65 \times 10^{-4}\ \textrm {s}$. The parameters 
$A$, 
$B$ and 
$C$ were varied and the simulation results for the modelled discharges (represented by 
$p_{\textrm {Ar}}$ in the respective cases) are shown in figure 2(b), where simulations resulting in full conversion are represented by 
${\circ}$ and simulations resulting in no RE generation are represented by 
$\triangle$. As shown in the figure, choosing 
$t_{\textrm {TQ}}$ according to (3.3) resulted in some RE generation in all cases except the no-RE case #35400. In this case, 
$t_{\textrm {TQ}}$ according to (3.3) may be unphysically long, but it has the desired effect of preventing any RE formation in the simulation. In general, our simulations indicate that 
$t_{\textrm {TQ}} < 0.03\ \textrm {ms}$ always results in full conversion, and 
$t_{\textrm {TQ}} > 0.35\ \textrm {ms}$ does not result in any RE generation in any of the simulated discharges. For all the presented simulations, the respective values of 
$t_{\textrm {TQ}}$ were estimated using (3.3), inserting 
$p_{\textrm {Ar}}$ for the respective discharge.
$(a)$ The timing of the current spike end, relative to the argon valve trigger time, is used as an indication of the time scale of the disruption dynamics, here plotted against 
$p_{\textrm {Ar}}$, showing an approximately inverse relationship. 
$(b)$ Assuming that the duration of the TQ also follows an inverse relationship, the corresponding parameter 
$t_{\textrm {TQ}}$ can be chosen so that none of the simulations of the discharges results in neither full conversion (
$\circ$) nor no RE generation (
$\triangle$).
3.3.3. Post-TQ equilibrium temperature
After the thermal quench and the related MHD mixing the ECE signal is still, in many discharges, blocked due to a high free electron density, and in addition the signal noise at the low post-TQ temperatures (some tens of eV) exceeds the signal by an order of magnitude in the few discharges where the density is low enough not to cut off the ECE signal. Lacking a reliable measurement, we thus need to estimate the post-TQ temperature. A reasonable estimate can be made by assuming equilibrium between ohmic heating and line radiation losses (Martin-Solis, Loarte & Lehnen Reference Martin-Solis, Loarte and Lehnen2017), so that the equilibrium temperature 
$T_{\textrm {EQ}}$ must satisfy
\begin{equation} E^2\sigma(T_{\textrm{EQ}},Z_{\textrm{eff}}(T_{\textrm{EQ}}))=\sum_in_{e}(T_{\textrm{EQ}})n_iL_i(T,n_e(T_{\textrm{EQ}})). \end{equation}
The electric field is denoted 
$E$, 
$\sigma$ is the plasma conductivity and 
$Z_{\textrm {eff}}$ is the effective argon ion charge. The sum goes over all possible charge states 
$i$ of argon, 
$n_i$ is the density of charge state 
$i$ and 
$L_i$ is the corresponding line radiation coefficient, obtained by interpolating data from the open-ADAS database (Summers Reference Summers2004). As before, 
$n_e$ and 
$T_e$ are the free electron density and temperature, respectively. The equation is solved iteratively. The equilibrium temperature during the plateau phase was slightly above 1 eV for all the modelled discharges.
The calculated equilibrium temperature for discharge #35408 is indicated with crosses in figure 3(b), which also shows the ECE-measured temperature (
${\cdots \cdots \cdots }$). As shown, the temperature used as input for the code simulations (———) is taken as that given by (3.2) until this falls below the calculated equilibrium temperature, after which the calculated equilibrium temperature from (3.4) is used.
The ECE free electron temperature 
$T_e$ signal used in the simulations (———) and the measured temperature (
${\cdots \cdots \cdots }$) during 
$(a)$ the entire simulated time span and 
$(b)$ around the end of the TQ. The end of the measured 
$I_p$ spike is also marked by a vertical line for reference. In 
$(b)$, the equilibrium temperature according to (3.4) (
$\times$) is also shown. The plot of the ECE signal in 
$(b)$ demonstrates the issue with the signal noise. All data for AUG discharge #35408.
4. Results
Using the kinetic model described in § 2 and experimental data as described in § 3, the electron distribution function and the associated current density are calculated. After the disruption, the ohmic current density drops rapidly in agreement with the measured 
$I_p$ while an RE current is formed that completely comes to dominate the total current density. Then follows a plateau phase during which the current density dissipates at a lower rate due to the interaction between the runaway electrons and the bulk plasma. For all discharges, we model the first 30 ms after the injection valve trigger, which ensures that the modelled time interval covers the entire CQ and part of the plateau phase for all the modelled cases.
4.1. Calculated current densities
The total and RE current densities are calculated by code by integration of the distribution function in the direction parallel to the magnetic field. These current densities are shown in figure 4(a), and the measured total plasma current 
$I_p$ is plotted for comparison, divided by the conversion factor 
$A_{\textrm {eff}} = 0.63\ \textrm {m}^{2}$ that was introduced in § 3.1. The motivation for this value of 
$A_{\textrm {eff}}$, however, may not be valid for the post-disruption phase, so the comparison between the scaled total current and the calculated current density is only indicative.
$(a)$ The calculated total and runaway current densities (———) are shown along with the total current divided by 
$A_{\textrm {eff}} = 0.63\ \textrm {m}^{2}$ (
${\cdots \cdots \cdots }$) for AUG discharge #35408. The end of the CQ (defined as 
$j_{\textrm {RE}}/j_{\textrm {tot}} = 0.99$) as marked with vertical lines. 
$(b)$ The calculated total current density for #35408 for cases with (———) and without (
$-\;-\;-\;-$) screening, as well as for different values of 
$t_{\textrm {TQ}}$.
As described in § 3.3.2, the TQ time scale strongly affects the resulting RE generation. This fact is demonstrated in figure 4(b), where the calculated total current density is displayed again, and compared with the same quantity for 
$t_{\textrm {TQ}} = 0.03\ \textrm {ms}$ and 
$t_{\textrm {TQ}} = 0.30\ \textrm {ms}$. For the reference case #35408, 
$t_{\textrm {TQ}} = 0.08 \ \textrm {ms}$. The profound effect of the screening of partially ionized impurities, mentioned in § 2, is also shown in figure 4(b), where we also, for comparison, have plotted the total current density calculated for the reference case #35408 with 
$t_{\textrm {TQ}} = 0.08\ \textrm {ms}$, but turning off the partial screening effects in the simulation, i.e. the simulation models full screening of all impurities irrespective of ionization state.
4.2. Momentum distributions
The time evolution of the calculated current densities can be more thoroughly understood by observing the time evolution of the electron momentum distribution function in the parallel (toroidal) direction. This is shown in figure 5(a) for AUG discharge #35408, where each line represents a separate time step. The nonlinear scale of the colour axis reflects the fact that the momentum distribution changes rapidly during the CQ but very little during the plateau phase.
$(a)$ The simulated parallel electron momentum distribution for AUG discharge #35408. The parallel momentum distribution functions at 5.2 ms and 6.5 ms (end of CQ) are marked with 
$-\;\cdot \;-\;\cdot \;-\;\cdot \;-\;\cdot $ and 
$-\;-\;-\;-$ respectively. An arrow marks the hot-tail ‘bump’. 
$(b)$ The corresponding simulated two-dimensional electron momentum distribution at 5.2 ms (upper panel, corresponding to 
$-\;\cdot \;-\;\cdot \;-\;\cdot \;-\;\cdot $ in figure 5a) and at 6.5 ms (lower panel, corresponding to 
$-\;-\;-\;-$ in figure 5a). The momentum 
$p$ is given, in both figures, in units of 
$mc$, with 
$m$ being the electron mass and 
$c$ the speed of light.
During the first 5 ms, the low-momentum part of the initially Maxwellian distribution narrows due to the rapid cooling while the high-momentum tail remains almost constant, i.e. these electrons do not lose momentum. After approximately 5 ms, the high-momentum tail (or ‘hot-tail’) starts to gain momentum and forms a ‘bump’ (marked with an arrow in figure 5a), separate from the thermal Maxwellian which is now indistinguishable from the vertical axis. The ‘bump’ is then gradually accelerated and simultaneously the avalanche mechanism gives rise to an approximately exponentially decreasing distribution of fast electrons at lower momenta.
A contour plot of the complete 2-D momentum distribution is shown in figure 5(b), for the two time instances marked by different line styles in figure 5(a). The upper panel displays the momentum distribution when the hot-tail generated seed is most pronounced, at 5.2 ms. The lower panel displays the momentum distribution at 6.5 ms, i.e. the end of the CQ when 
$j_{\textrm {RE}}/j_{\textrm {tot}} = 0.99$, marked in figure 4(a).
4.3. RE generation mechanisms
Figure 7 shows the runaway generation rate 
$(1/n_e)(\textrm {d}n_{\textrm {RE}}/\textrm {d}t)$ for two different discharges with different initial temperatures but similar injected argon pressures 
$p_{\textrm {Ar}}$. Runaway electrons are, in this case, defined as having 
$p > 0.75\ \textrm {m}_{e}c$ (with 
$m_{e}$ the electron rest mass and 
$c$ the speed of light) and the total RE generation rate is directly derived from the simulation output as the increase in the fraction of RE electrons during each time step, divided by the length of the time step.
code calculates the evolution of the electron distribution, and as such does not categorize REs by generation mechanisms. For a detailed analysis of the RE current formation mechanisms, the Dreicer RE generation rate was also calculated using a neural network described in Hesslow et al. (Reference Hesslow, Unnerfelt, Vallhagen, Embréus, Hoppe, Papp and Fülöp2019b), and the avalanche growth rate was also calculated using the semi-analytical formula developed in Hesslow et al. (Reference Hesslow, Embréus, Vallhagen and Fülöp2019a). As shown in figure 7, the analytically calculated avalanche growth rate describes the total simulated growth rate well, except for in the beginning, where the hot-tail seed can be seen as a peak approximately half a ms after the end of the 
$I_p$ spike.
The Dreicer RE seed generation is seen as a peak immediately after the hot-tail peak. This timing is expected since the hot-tail generation is directly connected with the TQ whereas the Dreicer generation is a consequence of the electric field generated, with some delay, after the TQ. The Dreicer generation rate was in general very small as compared with the hot-tail and avalanche generation, and smaller for larger 
$p_{\textrm {Ar}}$, which is clearly illustrated by figure 6. The maximal Dreicer RE seed generation rate (
$G_{\textrm {Dreicer}}$) is never larger than 
$10^{-5}$ times the maximal hot-tail RE seed generation rate, and the relative importance ranges all the way down to 
$10^{-18}$ for the largest Ar injection pressure (0.9 bar, in discharge #31318). The hot-tail RE seed generation was evaluated by subtracting the calculated avalanche generation from the total RE generation rate given by the simulations.
The ratio between the maximal Dreicer RE seed generation rate 
$G_{\textrm {Dreicer}}$ and the maximal hot-tail RE seed generation rate 
$G_{\text {hot-tail}}$, as a function of injected Ar pressure 
$p_{\textrm {Ar}}$.
RE generation rates, total and specified per generation mechanism, for discharges #34084 and #34140. The Dreicer RE generation rate is scaled to be visible in the respective panels – note the different scaling factors given in the plots. The hot-tail RE seed is seen as a peak in the total RE generation rate approximately half a ms after the end of the 
$I_p$ spike.
Generation rates for #34084 and #34140 are shown to display the effect of the initial temperature. These two discharges have very similar 
$p_{\textrm {Ar}}$, but different initial temperatures 
$T_{e0}$ (5.2 and 6.9 keV, respectively), and the hot-tail RE generation is found to be significantly higher for 6.9 keV, i.e. it is found to increase with increasing temperature. The Dreicer generation rates were scaled to be distinguishable in the respective plots.
4.4. Plateau phase current and dissipation rate
In our simulations, the fraction of the ohmic current which is converted to RE current during the CQ is sensitive to 
$t_{\textrm {TQ}}$ and the argon assimilation factor 
$f$. However, on AUG the measured post-CQ plasma current shows weak correlation with individual plasma parameters such as temperature, density or injected argon quantity; as long as these parameters are within the range of the RE generation window. The post-CQ calculated current density is plotted against the measured post-CQ total plasma current in figure 8(a), which shows that the measured post-CQ total plasma current is fairly constant for all RE discharges, whereas the simulation results show a much larger variability and no correlation with the measurement data. This would suggest that there are negative feedback effects (such as runaway seed transport) which are not captured by the 0-D kinetic modelling.
The calculated post-CQ current density plotted against 
$(a)$ the measured post-CQ total plasma current and 
$(b)$ the injected argon pressure 
$p_{\textrm {Ar}}$. In 
$(b)$, a linear fit has been added to visualize the trend above the injection threshold for significant RE generation. The linear correlation excludes the outliers #34084 (
$\triangledown$) and #31318 (
$\circ$) as well as the discharges #35400 and #35401 in which 
$p_{\textrm {Ar}}$ falls below the threshold.
However, figure 8(b) shows a correlation between the post-CQ calculated current density and 
$p_{\textrm {Ar}}$. The experimental observations show that there is a certain threshold 
$p_{\textrm {Ar}}$ under which no significant RE current is formed, since it does not lead to a TQ and a CQ quick enough for RE generation. Above this threshold, the post-CQ calculated current density is generally smaller for larger 
$p_{\textrm {Ar}}$. Two outliers are noted, however. In discharge #34084 (
$\triangledown$) the comparatively low initial free electron temperature (5.2 keV) results in a low hot-tail RE generation, as also shown in figure 7(b), and hence a low post-CQ current density. In discharge #31318 (
$\circ$), the initial temperature was instead comparatively high (9.3 keV), resulting in high hot-tail generation and an accordingly high post-CQ current density.
Given that we calculate a current density but only have access to an integrated total plasma current measurement, we cannot assess directly how well the simulation reproduces the on-axis current density development during the initial rapid current decrease. The current density profile changes drastically already during the MHD phase, and thus, the previously discussed conversion factor 
$A_{\textrm {eff}}$ is no longer valid. However, if we assume that the current density profile remains roughly constant throughout the plateau phase, i.e. that the conversion factor between current and current density remains constant during the plateau phase, we can make a meaningful comparison between the calculated current density dissipation rate and the measured current dissipation rate during the plateau phase. As described in § 3.2.2, it has been observed that the plateau phase current dissipation rate scales linearly with injected Ar quantity, under the condition that there is some RE generation but not full conversion of the initial current into RE current. As discussed in § 3.3.2, the thermal quench time parameter 
$t_{\textrm {TQ}}$ was tuned to fulfil this condition, but the dissipation rate is not sensitive to 
$t_{\textrm {TQ}}$ when the condition is fulfilled. As shown in figure 1(b), the dependence of the dissipation rate on the injected amount of Ar is reproduced by the simulations.
5. Discussion and conclusions
We have presented kinetic modelling of runaway generation and dissipation in argon-induced disruptions in ASDEX Upgrade, where the initial on-axis free electron temperature, the free electron density, the amount of injected argon and the initial current density have been based on experimental data. The timing and fraction of the argon assimilation and the time scale of the thermal quench were estimated. The fraction of argon assimilated in the plasma volume was fixed by fitting the calculated current density dissipation rates with the experimentally observed current dissipation rates during the plateau phase. The argon was assumed to assimilate in the on-axis plasma during the MHD phase, which coincides with the spike in the measured plasma current preceding the CQ. The time scale of the TQ was assumed to be inversely proportional to the injected argon amount, similarly to the observed timing of the current spike relative to the argon injection valve trigger time. These parameters could be estimated using the same assumptions for a set of eleven discharges with varying initial temperatures and injected argon amounts, yielding reasonable simulation results for all cases, i.e. neither full conversion nor complete CQ except for the no-RE discharge #35400. A 0D kinetic simulation cannot capture certain effects (such as e.g. runaway seed transport) and therefore the simulation results are not expected to give quantitative agreement with the experimental results.
Simulations show that, above a threshold injected argon quantity, a larger injection leads to a lower post-CQ current density, which is expected since the presence of argon increases energy loss from the plasma. The simulations also show that hot-tail RE generation is the most important RE generation mechanism in all the modelled discharges, having a significant impact on the post-CQ current density. For quantitatively accurate predictions of the plasma current, more elaborate models including transport phenomena could be used, such as the one by Linder et al. (Reference Linder, Fable, Jenko, Papp and Pautasso2020). code could also be coupled with a transport code such as go, as done by Hoppe et al. (Reference Hoppe, Hesslow, Embreus, Unnerfelt, Papp, Pusztai, Fülöp, Lexell, Lunt and Macusova2020). Such simulations are underway, following up on the simulations presented herein.
Acknowledgements
The authors are grateful for fruitful discussions with M. Hoppe. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under grant agreement No 633053 and from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No 647121. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The work was also supported by the Swedish Research Council (Dnr. 2018-03911) and the EUROfusion – Theory and Advanced Simulation Coordination (E-TASC).
Editor Paolo Ricci thanks the referees for their advice in evaluating this article.
Declaration of interests
The authors report no conflict of interest.
Open access
