Kinetic modelling of runaway electron generation in argon-induced disruptions in ASDEX Upgrade

Massive material injection has been proposed as a way to mitigate the formation of a beam of relativistic runaway electrons that may result from a disruption in tokamak plasmas. In this paper we analyse runaway generation observed in eleven ASDEX Upgrade discharges where disruption was triggered using massive gas injection. We present numerical simulations in scenarios characteristic of on-axis plasma conditions, constrained by experimental observations, using a two-dimensional momentum space description of the runaway dynamics with self-consistent electric field and temperature evolution. We describe the evolution of the electron distribution function during the disruption, and show that the runaway seed generation is dominated by hot-tail in all of the simulated discharges. We reproduce the observed dependence of the current dissipation rate on the amount of injected argon during the runaway plateau phase. Our simulations also indicate that above a threshold amount of injected argon, the current density after the current quench is generally lower for higher argon densities. This trend is not observed in the experiments due to the control system which we do not model.


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 plasmafacing components (Hender et al. 2007). 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 -SPI) to avoid or dissipate the runaway electrons (Hollmann et al. 2015a;Lehnen et al. 2015), and its efficiency has been demonstrated in medium-sized tokamaks (Hollmann et al. 2015b;Pautasso et al. 2016;Paz-Soldan et al. 2017;Carnevale et al. 2018;Reux et al. 2015;Esposito et al. 2017;Coda et al. 2019;Mlynar et al. 2019). 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 MGIinduced disruptions provide a valuable dataset to properly understand the physics of runaway electron formation and dissipation in such scenarios. To gain such understanding, theoretical models have been formulated and implemented in computational tools which can be used to model disruptions (Breizman & Aleynikov 2017). 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 section 2 and in detail in the paper by Stahl et al. (2016). 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, section 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 sections 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 1959). 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 1979;Rosenbluth & Putvinski 1997).
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. 2004). 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. 2005). A simplified analytical model for hot-tail generation has been formulated by Smith & Verwichte (2008) but comparisons with kinetic simulations show that this model underestimates the runaway density by an order of magnitude (Stahl et al. 2016), 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. (2020), where the astra-strahl 1.5D transport code (Fable et al. 2013;Dux et al. 1999) was used, including reduced kinetic models for Dreicer and avalanche runaway generation (Hesslow et al. 2019b,a). 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, 3D MHD evolution and the atomic physics needed to determine ionization states.

Kinetic Modelling
The relativistic finite-difference Fokker-Planck solver code (Stahl et al. 2016) simulates electron dynamics in plasmas in 2D momentum space. The plasma is assumed to be homogenous, i.e. it is a 0D 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 (1989), detailed in Hesslow et al. (2019b) and a simplified largeangle collision operator (Rosenbluth & Putvinski 1997). Screening of partially ionized impurities is taken into account according to the model described by Hesslow et al. (2018a). 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 et al. 2018).
The electric field E and the plasma current density j is calculated self-consistently throughout the simulation, according to Hesslow et al. (2018b): where the inductance L is given by (2. 2) 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 = |Ψ p /I p | (Boozer 2018), where Ψ p is the poloidal magnetic flux and I p is the plasma current. To calculate the on-axis magnetic flux Ψ p requires the plasma current density profile, which we do not know with an accuracy that motivates the use of the more accurate expression. 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 2005) and in a self-consistent calculation, this needs to be compensated for. As the conductivity obtained by the test-particle operator σ CODE,tp is proportional to the fully relativistic electric conductivity σ BK (Braams & Karney 1989) for a wide range of effective charges and temperatures, we can model the plasma current density as where j CODE,tp is the current density calculated by code using the test particle operator, and ∆σ(Z eff ) = σ BK − σ CODE,tp . The calculation of the effective ionic charge Z eff is described in section 3.2.3.
Magnetic field B = 2.5 T Major radius R = 1.65 m Minor radius a = 0.50 m Initial plasma current Ip0 = 0.76 MA Table 1: Common parameters of all modelled discharges. The initial current is slightly lower, 0.71 MA, in discharge #31318.

ASDEX Upgrade
In this paper we model ASDEX Upgrade discharges specifically designed for the study of runaway electron dynamics (Pautasso et al. 2016). ASDEX Upgrade is a medium sized tokamak at the Max Planck Institute of Plasma Physics in Garching, Germany. The typical runaway electron scenario uses a low density (initial free electron density n e0 ≈ 3·10 19 m −3 ), inner wall limited, circular (elongation κ ≈ 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.
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 slightly in the dataset (T e0 = 7.1 ± 1.1 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.
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 (McCarthy 1999). For the modelled discharges, the initial on-axis current density j 0 was estimated to approximately 1.2 MA/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 m 2 . Since this conversion factor has the unit of m 2 , it will be referred to as an "effective area", A eff . Application of this conversion factor results in an initial on-axis current density of 1.21 MA/m 2 ± <1% for all discharges except the very early discharge #31318 in which the initial on-axis current density is 1.13 MA/m 2 .

Densities
The disruptions were triggered by injection of argon into the plasma (Pautasso et al. 2016;Fable et al. 2016). 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 distribution function. Thus, these parameters must be estimated and used as input to the simulations.  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 CO 2 interferometry and the initial electron temperature by the ECE measurements.

Free electron density
The free electron density is measured by 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. (2016) 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 o and t e respectively. We assume that the on-axis free electron density remains constant at the pre-injection value until t o , and that after t e the plasma is homogenous enough for the measured line-averaged free electron density to be representative of the on-axis value. The data is smoothed using the rloess algorithm in matlab (MATLAB 2017) to avoid numerical difficulties caused by the signal noise. Between t o and t e , 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 1a, along with the measured free electron density.

Argon density
The injected amount of argon is quantified by the pressure p 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 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 section 4.1, the plasma current drops suddenly during the disruption and then dissipates more slowly during the so-called plateau phase. As shown in figure 1b, experimental data shows that there is a linear correlation between the injected argon amount and the current dissipation rate in this phase, and the assimilation factor was estimated by comparison of measured and simulated current density dissipation rates, dj/dt, 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 eff = 0.63 m 2 explained in section 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 entire plateau phase.
The current density dissipation rate is calculated from the current density given by the kinetic simulations at 30 ms after the argon valve trigger. This time point is well into the plateau phase for all simulated discharges. The calculated dissipation rates are plotted in figure 1b for f = 20%, along with a linear fit which approximately reproduces the experimentally deduced slope. The positive offset of the experimental current dissipation rates comes from the fact that the current in the experiment is also driven by an external electric field, generated by the control system. The electric field in the simulations is developing self-consistently, without the external electric field, as it 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 zero until MHD mixing occurs at t o (Fable et al. 2016). After t e , the argon density is assumed to be constant, at a level given by Between t o and t e , the argon density is assumed to increase linearly. The calculated current was shown not to be sensitive to the detailed time evolution of the density within the time range between t o and t e . 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 modeled discharges.

Charge states
For the purpose of the kinetic simulations, the average charge state Z eff of the argon is inferred from measurements, i.e. not calculated from atomic physics. Z eff is estimated by dividing the density of free electrons attributable to argon by the density of argon atoms. All free electron density above the initial value n e0 is 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.
Detailed atomic physics modelling yields a broader distribution over multiple ionization states (Linder et al. 2020). 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 (a) (b) Figure 1: (a) The free electron density measured by interferometry ( ) 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 measured total plasma current I p and the calculated current densities. Note that the measured current dissipation rate has been scaled with A eff = 0.63 m 2 .
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.

Thermal quench
The presence of impurities in the plasma abruptly increases radiative energy losses, which leads to the rapid thermal quench (TQ). The simulation results are very sensitive to the duration of the TQ.

Measured free electron temperature
Previous work by Fable et al. (2016) 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 measured free electron temperature, averaged over a circular 15 cm 2 area surrounding the plasma axis, is shown in figure 3 ( ). However, the ECE signal is cut off from approximately 1.5 ms after the beginning of the argon injection until t e , due to the high free electron density at the plasma edge (Fable et al. 2016). Thus, we calculate the initial free electron temperature T e0 as the measured value averaged over the central 15 cm 2 area and the time period of 1 -1.5 ms after the argon injection valve trigger. T e is then assumed to decrease exponentially with time t according to equation 3.2, as argued by Smith & Verwichte (2008).
where t T Q is determined as described in section 3.3.2. The final temperature T e,final can be any low value, because after the TQ, the values given by equation 3.2 are no longer used, but instead the equilibrium temperature is determined using the equilibrium assumption described in section 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 equation 3.2, so T e,final is fixed to 0.5 eV to make sure that the temperature given by equation 3.2 falls below the equilibrium temperature in all simulations.

Thermal quench time
The thermal quench time is quantified by the parameter t T Q , whose significance is shown in equation 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 T Q . The choice of t T Q proves to be very important for the RE generation. For too large t T Q , there is no RE generation in the simulation, whereas for too small t T Q , 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 2008). 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 T Q was chosen to reproduce this result. Using the same value of t T Q in each shot did not yield the desired result -a low t T Q resulted in no RE formation in the high-injection cases, whereas a high t T Q resulted in full conversion in the low-injection cases. Also, a constant value of t T Q 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  2a, an inverse relationship is found between these parameters indicating that the onset of the current quench is faster for higher injected argon pressures p Ar . Therefore in the modelling we use the assumption that the thermal quench time is shorter for higher p Ar . The relationship between t TQ and p Ar is given by with A = 7 · 10 −4 bar·s, B = 0.1 bar and C = 0.65 · 10 −4 s. The parameters A, B and C were varied and the simulation results for the modelled discharges (represented by p Ar in the respective cases) are shown in figure 2b, where simulations resulting in full conversion are represented by and simulations resulting in no RE generation are represented by △. As shown in the figure, choosing t T Q according to equation 3.3 resulted in some RE generation in all cases except the no-RE case #35400. In general, our simulations indicate that t T Q < 0.03 ms always results in full conversion, and t T Q > 0.35 ms does not result in any RE generation in any of the simulated discharges.

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 shots 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 , so that the temperature must satisfy: j 2 Ω σ(T e , Z eff (T e )) = i n e (T e )n i L i (T, n e (T e )).
( 3.4) The ohmic current density is denoted j Ω , σ is the plasma conductivity and Z 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 (a) (b) Figure 3: The free electron temperature T e used in the simulations ( ) and the measured temperature ( ) 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 equation 3.4 (×) is also shown. All data for AUG discharge #35408. by interpolating data from the open-ADAS database (Summers 2004). 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 shot #35408 is indicated with crosses in figure 3b, which also shows the ECE-measured temperature ( ). As shown, the temperature used as input for the code simulations ( ) is taken as that given by equation 3.2 until this falls below the calculated equilibrium temperature, after which the calculated equilibrium temperature from equation 3.4 is used.  Figure 4: (a) The calculated total and runaway current densities ( ) are shown along with the total current divided by A eff = 0.63 m 2 ( ) for AUG discharge #35408. The end of the CQ (defined as j RE /j 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 T Q .

Results
Using the kinetic model described in section 2 and experimental data as described in section 3, the electron distribution function and associated current density is 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.

Calculated current densities
The calculated total and RE current densities are shown in figure 4a, and the measured total current I p is plotted for comparison, divided by the conversion factor A eff = 0.63 m 2 explained in section 3.1. The motivation for this value of A eff , however, is only valid for the pre-disruption phase, so the comparison between the scaled total current and the calculated current density is only indicative.
As described in section 3.3.2, the TQ time scale strongly affects the resulting RE generation. This fact is demonstrated in figure 4b, where the calculated total current density is displayed again, and compared with the same quantity for t T Q = 0.03 ms and t T Q = 0.30 ms. For the reference case #35408, t T Q = 0.08 ms. The profound effect of the screening of partially ionized impurities, mentioned in section 2, is also shown in figure  4b, where we also, for comparison, have plotted the total current density calculated for the reference case #35408 with t T Q = 0.08 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.

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 5a for AUG discharge #35408, where each line represents a separate time step. The non-linear scale of the color axis reflects the fact that the momentum distribution changes rapidly during the CQ but very little during the plateau phase.
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 about 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 2D momentum distribution is shown in figure 5b, for the two time instances marked by different line styles in figure 5a. 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 RE /j tot = 0.99, marked in figure 4a. Figure 6 shows the runaway generation rate (1/n e )(dn RE /dt) for two different discharges with different intial temperatures but similar injected argon pressures p Ar . Runaway electrons are, in this case, defined as having p > 0.75 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.

RE generation mechanisms
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. (2019b), and the avalanche growth rate was also calculated using the semi-analytical formula developed in Hesslow et al. (2019a). As  Figure 6: 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 plots -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.
shown in figure 6, 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 Ar . In the discharge with the largest p Ar , 0.9 bar for #31318, the Dreicer generation rate was 1.8·10 −14 times the maximum total generation rate, whereas for the discharge with the smallest p Ar leading to RE formation, 0.15 bar for #35401, the corresponding ratio was 0.083, i.e. a difference by almost 13 orders of magnitude. This shows that there are qualitative differences between the generation mechanisms, depending on the quantity of injected argon.
Generation rates for #34084 and #34140 are shown to display the effect of the initial temperature. These two discharges have very similar p 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 #34140. The Dreicer generation rates were scaled to be distinguishable in the respective plots. The scaling factors for discharges #34084 and #34140 fall in between the extremes mentioned.

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 T Q 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 7a, 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. However, figure 7b shows a correlation between the post-CQ calculated current density and p Ar . As expected and experimentally observed, there is a certain threshold p Ar under which no RE current is formed, since it does not lead to a thermal quench quick enough for RE generation. Above this threshold, however, the post-CQ calculated current density is smaller for larger p Ar . Two outliers are noted, however. In discharge #34084 (▽) the comparatively low initial free electron temperature (5.2 keV) results in a low hot-tail RE generation, as also shown in figure 6b, and hence a low post-CQ current density. In discharge #31318 ( ), 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 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 section 3.2.2, it has been observed that the plateau phase current dissipation rate scales linearly with injected Ar quantity. As shown in figure 1b, this dependence is reproduced by the simulations.

Discussion and conclusions
We have presented kinetic modelling of runaway generation and dissipation in argoninduced 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. However, the electric field generated by the control system in the experiment is not included in the simulations, 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. (2020).