Evaluation of the Dreicer runaway generation rate in the presence of high-$Z$ impurities using a neural network

Abstract

Integrated modelling of electron runaway requires computationally expensive kinetic models that are self-consistently coupled to the evolution of the background plasma parameters. The computational expense can be reduced by using parameterized runaway generation rates rather than solving the full kinetic problem. However, currently available generation rates neglect several important effects; in particular, they are not valid in the presence of partially ionized impurities. In this work, we construct a multilayer neural network for the Dreicer runaway generation rate which is trained on data obtained from kinetic simulations performed for a wide range of plasma parameters and impurities. The neural network accurately reproduces the Dreicer runaway generation rate obtained by the kinetic solver. By implementing it in a fluid runaway-electron modelling tool, we show that the improved generation rates lead to significant differences in the self-consistent runaway dynamics as compared to the results using the previously available formulas for the runaway generation rate.

1 Introduction

If the electric field exceeds a critical field in plasmas, the accelerating force on fast electrons overcomes the collisional friction. The electrons are then accelerated to relativistic energies – they run away (Dreicer 1959) – until the electric force is balanced by another mechanism such as synchrotron radiation. Runaway acceleration occurs in solar flares (Holman 1985), lightning discharges (Dwyer 2007) and in magnetic fusion devices (Helander, Eriksson & Andersson 2002; Breizman et al. 2019). In tokamaks, runaway electrons form when large electric fields are present: during startup of the discharge, radio-frequency current drive or disruptions – a sudden termination of a tokamak discharge. In high-current, reactor-scale machines, a disruption has the potential to convert a major part of the plasma current to a relativistic runaway-electron beam. Such a beam would severely damage plasma-facing components if it is not well controlled (Lehnen et al. 2015).

An accurate modelling capacity of runaway electrons is essential to evaluate the different methods aimed at mitigating their effects. To describe the full evolution of the temperature, plasma composition, electric field, magnetic equilibrium and runaway current, it would be necessary to simultaneously solve for the relativistic momentum-space dynamics, the magnetic field evolution (including breakup of magnetic surfaces) and the transport (including both collisional and turbulent processes). This is, however, unfeasible with currently available computational resources and will likely remain so in the foreseeable future.

Until simulations of the full plasma evolution during a disruption can be realized, as an intermediate step, transport codes and equilibrium solvers could be coupled with analytical formulas for runaway generation rates. This means that instead of evolving the full runaway-electron distribution, only key quantities such as the runaway number density would be considered and computed from the instantaneous electric field and background plasma parameters. In such fluid models, the runaway-electron density evolves by analytical generation rates describing Dreicer, hot-tail and avalanche generation, as well as tritium decay and Compton scattering of $\unicode[STIX]{x1D6FE}$ -rays (which can be emitted by the activated wall in the nuclear phase of tokamak operation). This approach has been used in the past to gain insight into the runaway-electron dynamics and electric-field diffusion; some examples are the go code (Smith et al. 2006; Gál et al. 2008; Fehér et al. 2011; Papp et al. 2013) and the work by Martín-Solís, Loarte & Lehnen (2017).

Although runaway fluid models are useful to understand runaway dynamics, present tools use runaway generation rates that lack several important effects, including the effect of synchrotron radiation (Stahl et al. 2015), bremsstrahlung (Embréus, Stahl & Fülöp 2016) and screening effects in partially ionized plasmas (Hesslow et al. 2017). The need for amendments is particularly pronounced in scenarios involving disruption mitigation by massive material injection. The injected impurity ions will radiate the thermal energy content of the plasma during the thermal quench, and will become weakly ionized in the resulting cold plasma, implying that the nuclei will be partially screened by the bound electrons in interactions with fast electrons. Recent studies indicate substantial differences in the runaway dynamics compared to the fully ionized case: the effect of partial screening might give order of magnitude differences in both the Dreicer generation rate (Hesslow et al. 2018) and the avalanche multiplication factor (Hesslow et al. 2019).

While the avalanche growth rate calculation (Rosenbluth & Putvinski 1997) was recently generalized to include the effect of partial screening and radiation reaction (Hesslow et al. 2019), a similar generalization of the Dreicer generation rate calculation by Connor & Hastie (1975) appears intractable. This is because the Dreicer rate is exponentially sensitive to the plasma properties at near-thermal energies, where the collision frequencies have a complicated energy dependence in the presence of cold impurities. Instead, Martín-Solís, Loarte & Lehnen (2015) suggested to replace the Dreicer field $E_{\text{D}}$ and the effective charge $Z_{\text{eff}}$ in the Connor & Hastie (1975) formula, but otherwise keep the same form of the expression. This approach has however not been validated by solutions of the kinetic equation.

A different approach to improving the Dreicer generation rate is to use a large database of numerical solutions of the kinetic equation to train a neural network. Neural networks have proved to be useful to fit complicated, high-dimensional results to multi-parameter models, which is the case here since the density of each impurity species, in combination with a wide range of plasma parameters, gives a large number of free parameters. Neural networks are widely used in many areas of physics, including fusion plasma physics; see e.g. Svensson, von Hellermann & König (1999), Clayton et al. (2013), Citrin et al. (2015) and Boyer, Kaye & Erickson (2019).

In this paper, we construct a neural network that determines the Dreicer generation rate including the effect of partial screening as well as an energy-dependent Coulomb logarithm. After detailing the kinetic model and its validity (§ 2), we use a neural network to create a model of the Dreicer generation rate based on kinetic simulations (§ 3). The resulting network successfully reproduces the runaway generation rates predicted by the original kinetic solver, and can be directly implemented into integrated models. As a proof of concept, we implement the network into the runaway fluid simulation tool go, and demonstrate that partial screening has a significant effect on runaway generation (§ 4). Finally, we discuss the applications of the model as well as possible improvements (§ 5).

2 Kinetic model

The Dreicer (1959, 1960) mechanism for runaway generation originates from the interplay between collisional energy diffusion and electric-field acceleration. In a time-independent background plasma, the Dreicer mechanism causes a constant particle flow through any momentum-space boundary beyond the runaway separatrix ¹ ; this flow rate defines the steady-state Dreicer generation rate,

(2.1) $$\begin{eqnarray}\unicode[STIX]{x1D6FE}\equiv \frac{\text{d}n_{\text{RE}}}{\text{d}t},\end{eqnarray}$$

where $n_{\text{RE}}$ is the runaway-electron number density.

To be useful in runaway fluid modelling, the system must be well described by a quasi-steady-state approximation, so that $\unicode[STIX]{x1D6FE}(E,T,n_{\text{e}},\ldots ,t)\approx \unicode[STIX]{x1D6FE}[E(t),T(t),n_{\text{e}}(t),\,\ldots ]$ , where $E,T,n_{\text{e}},\ldots \,$ are the plasma parameters which influence the Dreicer generation rate in steady state, such as the electric field $E$ , the electron temperature $T$ and the electron density $n_{\text{e}}$ . The quasi-steady-state generation rate thus requires sufficiently slowly varying parameters.

The most accurate analytical treatment of the Dreicer problem was carried out by Connor & Hastie (1975), who extended the results of Kruskal & Bernstein (1962) to include relativistic effects, resulting in

(2.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D6FE}=C{\displaystyle \frac{n_{\text{e}}}{\unicode[STIX]{x1D70F}_{\text{ee}}}}\left({\displaystyle \frac{E}{E_{\text{D}}}}\right)^{-(3/16)(1+Z_{\text{eff}})h}\exp \left[-\unicode[STIX]{x1D706}{\displaystyle \frac{E_{\text{D}}}{4E}}-\sqrt{\unicode[STIX]{x1D702}{\displaystyle \frac{(1+Z_{\text{eff}})E_{\text{D}}}{E}}}\right],\\ \unicode[STIX]{x1D706}=8{\displaystyle \frac{E^{2}}{E_{\text{c}}^{2}}}\left[1-{\displaystyle \frac{1}{2}}{\displaystyle \frac{E_{\text{c}}}{E}}-\sqrt{1-{\displaystyle \frac{E_{\text{c}}}{E}}}\right],\\ \unicode[STIX]{x1D702}={\displaystyle \frac{1}{4}}{\displaystyle \frac{E^{2}}{E_{\text{c}}(E-E_{\text{c}})}}\left[{\displaystyle \frac{\unicode[STIX]{x03C0}}{2}}-\arcsin \left(1-{\displaystyle \frac{2E_{\text{c}}}{E}}\right)\right]^{2},\\ h={\displaystyle \frac{1}{3}}{\displaystyle \frac{1}{(E/E_{\text{c}})-1}}\left[{\displaystyle \frac{E}{E_{\text{c}}}}+2\left({\displaystyle \frac{E}{E_{c}}}-2\right)\sqrt{{\displaystyle \frac{E}{E-E_{\text{c}}}}}-{\displaystyle \frac{Z_{\text{eff}}-7}{Z_{\text{eff}}+1}}\right].\end{array}\right\}\end{eqnarray}$$

Here, $E_{\text{D}}=n_{\text{e}}e^{3}\,\text{ln}\,\unicode[STIX]{x1D6EC}_{0}/(4\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}_{0}^{2}T)$ is the Dreicer field, $E_{\text{c}}=E_{\text{D}}T/(m_{\text{e}}c^{2})$ is the critical electric field, $Z_{\text{eff}}=\sum _{i}n_{i}Z_{i}^{2}/n_{\text{e}}$ is the effective plasma charge and $\unicode[STIX]{x1D70F}_{\text{ee}}=4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}^{2}m_{\text{e}}^{2}v_{\text{Te}}^{3}/(n_{\text{e}}e^{4}\ln \unicode[STIX]{x1D6EC}_{0})$ is the thermal electron collision time, with the thermal speed $v_{\text{Te}}=\sqrt{2T/m_{\text{e}}}$ , and where $\ln \unicode[STIX]{x1D6EC}_{0}$ is the Coulomb logarithm evaluated at the thermal speed (in the original work, $\ln \unicode[STIX]{x1D6EC}$ was assumed to be energy independent). The order-unity parameter $C$ is undetermined by Connor & Hastie (1975) but has been quantified in subsequent work to around $C\approx 1.0$ with our time normalization (Kruskal & Bernstein 1962; Jayakumar, Fleischmann & Zweben 1993), and the parameters $\unicode[STIX]{x1D706}$ , $\unicode[STIX]{x1D702}$ and $h$ constitute the relativistic generalization of the generation rate; they approach unity as $E/E_{\text{c}}\rightarrow \infty$ .

Figure 1. Effect of partial screening, radiation reaction and an energy-dependent Coulomb logarithm on the Dreicer generation rate. Panel (a) has singly ionized argon impurities with density $n_{\text{Ar}}=n_{\text{D}}=10^{20}~\text{m}^{-3}$ at a temperature of $T=10~\text{eV}$ , whereas panel (b) displays the effect of synchrotron and bremsstrahlung radiation reaction in a 5 T magnetic field plasma with $n_{\text{D}}=10^{20}~\text{m}^{-3}$ and $T=10~\text{keV}$ . In both (a) and (b), the solid black line shows the analytical formula from (2.2) with $C=1$ , which was derived in the completely screened limit (denoted ‘CS’ in the legend), using a constant Coulomb logarithm and neglecting the effect of radiation (denoted ‘no rad’). Black dots and red crosses represent the simplified ideal plasma model where screening and radiation effects are ignored; black dots show code simulations using a constant Coulomb logarithm, whereas red crosses account for its energy dependence. The blue solid line with plus markers shows results that account for all kinetic effects. For comparison, the black squares show the numerical results by Kulsrud et al. (1973) in (a), which were obtained in the non-relativistic limit.

Despite its apparent complexity, the generation rate in (2.2) neglects certain effects that are necessary for an accurate treatment of the problem. As illustrated in figure 1, the generation rate deviates significantly from (2.2) already when accounting for the energy dependence of the Coulomb logarithm $\ln \unicode[STIX]{x1D6EC}$ , an effect which is most pronounced at the lower temperatures of figure 1(a). By furthermore considering the effect of partial screening, the difference can reach several orders of magnitude.

At high temperatures (figure 1 b), the main deviation from the ideal behaviour is the energy dependence of the Coulomb logarithm. Since radiation reaction only becomes important at highly relativistic energies, these losses have a small effect on Dreicer generation except for where the critical momentum fulfils $p_{\text{c}}\gg 1$ , which is obtained at near-critical electric fields if the temperature is high. As shown in figure 1(b), the effect of synchrotron and bremsstrahlung radiation reaction is however modest even at relatively synchrotron-favouring parameters with $T=10~\text{keV}$ , $B=5~\text{T}$ and $n_{\text{D}}=10^{20}~\text{m}^{-3}$ . In these simulations, radiation reaction mainly affects the generation rate close to the effective critical field, which is consistent with previous work (Stahl et al. 2015). We also note that Dreicer generation predominantly takes place at low temperatures during disruption scenarios, since the value of $E/E_{\text{D}}\sim j/(\unicode[STIX]{x1D70E}E_{\text{D}})\sim 1/\sqrt{T}$ decreases with temperature at constant current density. For the remainder of this work, we therefore neglect the effect of radiation reaction on the Dreicer generation rate.

We also disregard the effect of toroidicity, which may have an appreciable effect on Dreicer generation off the magnetic axis: if the bounce time is much shorter than the detrapping time, the generation rate is reduced due to magnetic trapping (Nilsson et al. 2015). Conversely, at high densities and electric fields $E\gg E_{\text{c}}$ , which can be present during tokamak disruptions, the Dreicer generation is approximately local (as in this work), but will be spatially non-uniform since the induced electric field decreases in magnitude with major radius (McDevitt & Tang 2019).

The results in this work, including figure 1, were obtained with the linearized Fokker–Planck solver code (Landreman, Stahl & Fülöp 2014; Stahl et al. 2016), which models a spatially homogeneous, magnetized plasma. The test-particle collision operator in code is given in appendix A. Since the collision operator is linearized, it is only valid for weak electric fields $E\ll E_{\text{D}}$ . Otherwise, a major part of the thermal electron distribution runs away and a nonlinear Fokker–Planck equation, solved by for example norse (Stahl et al. 2017), should be used. Conversely, the linear and nonlinear operators should agree at weak electric fields, and we verify in appendix A that the test-particle collision operator in code gives a similar generation rate to norse.

code has a model of partial screening based on the Born approximation (Hesslow et al. 2018). This is the most significant limitation of this work, since the Born approximation is strictly valid only for electron speeds obeying $v/c\gg Z\unicode[STIX]{x1D6FC}$ , where $\unicode[STIX]{x1D6FC}\approx 1/137$ is the fine-structure constant and $Z$ is the atomic number. For argon and neon, the validity of the Born approximation has been experimentally verified to extend beyond this requirement, down to kinetic energies of approximately 1 keV (Mott & Massey 1965). Below this threshold, the model for partial screening is approximate, although it has been asymptotically matched to the correct behaviour as $p\rightarrow 0$ . As the Dreicer generation rate is most sensitive to the dynamics near the critical momentum given by $p_{\text{c}}\gtrsim (E/E_{\text{c}}-1)^{-1/2}$ , the model is only strictly valid for $E/E_{\text{D}}(\%)\lesssim T/(20~\text{eV})$ , implying that the accuracy of our screening model is compromised for certain parameters.

Notably, we found that the generalization of the generation rate suggested by Martín-Solís et al. (2017) – to replace $Z_{\text{eff}}$ and $E_{\text{D}}$ by expressions involving the increased collision rates evaluated at an approximate value of the critical momentum $p_{\text{c}}$ – gave poor agreement with kinetic simulations. A more involved amendment of the Dreicer generation rate (2.2) is therefore required, and would likely result in a significantly more involved analysis than was done by Connor & Hastie (1975). For this reason, we resorted to the use of a neural network, which will be described in the following section.

3 Neural network model for the Dreicer generation rate

We used code to determine the steady-state momentum-space distribution of the fast electrons, from which we determine the normalized generation rate

(3.1) $$\begin{eqnarray}\bar{\unicode[STIX]{x1D6FE}}\equiv \frac{\unicode[STIX]{x1D70F}_{\text{ee}}}{n_{\text{e}}}\unicode[STIX]{x1D6FE}.\end{eqnarray}$$

In order to minimize the computational cost of the simulations, code was used in the steady-state mode, as described by Landreman et al. (2014) and detailed in appendix B.

code simulations were performed for a large number of points randomly sampled in the region described in table 1, where $n_{Z}$ is the impurity ion density, and $n_{\text{D}}$ is the density of deuterium (or other hydrogen species; the isotope does not affect the generation rate). For each ionization state of argon and neon, i.e. $\text{Ar}^{+n}$ for $n=0\cdots 18$ , and $\text{Ne}^{+m}$ for $m=0\cdots 10$ (one at a time), 8000 points were sampled. Additionally, 10 000 points were sampled without any impurities ( $n_{Z}=0$ ), and 10 000 points with a mix of different impurities with total density $n_{Z}$ . The maximum temperature was set to either 20 keV or twice the mean excitation energy, $2I_{Z}$ ; the latter for cases with a single impurity species. This is because a given charge state does not typically occur at higher temperatures. For example, $I_{\text{Ar}^{+}}=219.4~\text{eV}$ (Sauer, Oddershede & Sabin 2015), at which temperature argon would already be multiply ionized in equilibrium. Moreover, in our model for partial screening, the enhancement of the slowing-down collision frequency starts to extend into the thermal population at such temperatures (Hesslow et al. 2018), and the validity of the screening model starts to become questionable.

Table 1. Parameters used in the code simulations, their range and how they were sampled.

Rather than using the full set of impurity ion densities as input to the neural network, we use six derived parameters,

(3.2a-f ) $$\begin{eqnarray}\displaystyle \ln n_{\text{e}},\quad Z_{\text{eff}},\quad {\displaystyle \frac{n_{\text{e}}}{n_{\text{tot}}}},\quad \mathop{\sum }_{i}\frac{n_{i}}{n_{\text{tot}}}(Z_{i}^{2}-Z_{0,i}^{2}),\quad \mathop{\sum }_{i}\frac{n_{i}}{n_{\text{tot}}}Z_{0,i}Z_{i},\quad \text{and}\quad \mathop{\sum }_{i}\frac{n_{i}}{n_{\text{tot}}}{\displaystyle \frac{Z_{0,i}}{Z_{i}}}. & & \displaystyle \nonumber\\ \displaystyle & & \displaystyle\end{eqnarray}$$

Here, $Z_{i}$ and $Z_{0,i}$ are the atomic number and charge number of species $i$ , respectively, and $n_{\text{tot}}=\sum _{i}Z_{i}n_{i}$ is the total density of free and bound electrons. These derived parameters were chosen to both include the relevant parameters in the completely screened limit ( $n_{\text{e}}$ and $Z_{\text{eff}}$ ) parameters that naturally appear in the partially screened collision frequencies ( $n_{\text{e}}/n_{\text{tot}}$ and $\sum _{i}(n_{i}/n_{\text{tot}})(Z_{i}^{2}-Z_{0,i}^{2})$ ), as well as the last two parameters in (3.2), which vary significantly with plasma composition. This reduced the required size of the network, and allowed it to generalize better to other impurity species (for example, the neural network could accurately predict the Dreicer generation rate with carbon impurities although it was not trained for these). Accordingly, the input data were composed of 8 parameters: six density-related inputs, the normalized electric field and the natural logarithm of the temperature. This input, combined with the output $\bar{\unicode[STIX]{x1D6FE}}$ , was randomly split into a training and validation set with a 4 : 1 ratio.

Figure 2. Comparison between normalized runaway generation rates obtained from code and those from the neural network regression.

The neural network ² is designed as a multilayer perceptron described by the following series of matrix multiplications and function applications:

(3.3) $$\begin{eqnarray}\ln \bar{\unicode[STIX]{x1D6FE}}=\unicode[STIX]{x1D652}_{5}\tanh \{\unicode[STIX]{x1D652}_{4}\cdots \tanh (\unicode[STIX]{x1D652}_{1}\boldsymbol{x}+\boldsymbol{b}_{1})+\cdots +\boldsymbol{b}_{4}\}+\boldsymbol{b}_{5}.\end{eqnarray}$$

Here, $\boldsymbol{x}$ is the input vector described above, the matrices $\unicode[STIX]{x1D652}_{(k)}$ describe the weights between the layers and the vectors $\boldsymbol{b}_{(k)}$ are the biases. The activation function, $\tanh$ , is applied element-wise to the four hidden layers which were of size 20 (i.e. $\unicode[STIX]{x1D652}_{1}\in \mathbb{R}^{20\times 8}$ , $\{\unicode[STIX]{x1D652}_{2},\unicode[STIX]{x1D652}_{3},\unicode[STIX]{x1D652}_{4}\}\in \mathbb{R}^{20\times 20}$ and $\unicode[STIX]{x1D652}_{5}\in \mathbb{R}^{1\times 20}$ ). To determine the optimized weight and bias values, we used the Adam algorithm (Kingma & Ba 2014), which is a stochastic gradient descent method. The training of the neural network was implemented in Python using the library PyTorch (Paszke et al. 2017). The validation set was used to determine when the optimization was sufficiently converged and to avoid overfitting.

Due to feedback between the current and the electric field, an order-unity error in the Dreicer generation rate will have a marginal impact on the maximum runaway current at the end of a current quench. This implies that some errors can be accepted, although the dominating factors must be correctly modelled to capture the final runaway-electron current. A comparison between the regression neural network and code outputs for a set of 6500 points (different from both the validation and training set) is shown in figure 2. We found that 99.6 % of the neural network predictions were within a factor of two of the correct value of $\unicode[STIX]{x1D6FE}$ , and the mean absolute logarithmic (base 10) error was 0.0283. This indicates that the training data provided sufficient coverage of the parameter space. The evaluation time of the neural network is approximately five orders of magnitude faster than a steady-state code simulation. The difference is even larger if the neural network is compared to a full time-dependent simulation with varying parameters and simultaneously accounting for avalanche multiplication, which is the type of simulation the generation rates can replace when used in integrated modelling.

The typical quality of the fits is shown in figure 3, where we display the normalized generation rate for a temperature of $T=10~\text{eV}$ and deuterium density $n_{\text{D}}=10^{20}~\text{m}^{-3}$ . Figure 3(a) shows the runaway generation rate as a function of normalized electric field calculated by the neural network, together with code simulations, showing excellent agreement. For reference, the analytical formula (2.2) is also included (dashed line). Figure 3(b) shows the generation rate as a function of the density of triply ionized neon normalized to the density of hydrogen. Again, the agreement between the simulations and neural network is excellent, whereas the disagreement with the analytical expression is substantial.

Figure 3. Comparison between the normalized generation rate obtained by the neural network (solid lines), simulations with the code Fokker–Planck solver (blue dots) and the analytical formula from (2.2) with $C=1$ (dashed lines). The temperature was set to 10 eV, and $n_{\text{D}}=10^{20}~\text{m}^{-3}$ . In (a) $n_{\text{Ne}^{3+}}/n_{\text{D}}=1$ , and in (b) $E/E_{\text{D}}=0.04$ .

4 Application in runaway current modelling

To demonstrate the impact of the modified Dreicer generation rates, we use the neural network in a self-consistent simulation of the electric field and current profile evolution performed by the go numerical tool (Smith et al. 2006). Instead of modelling the velocity space dynamics for the electrons in the runaway region, go considers only their total density $n_{\text{RE}}$ . It solves the coupled equations for the runaway generation and resistive diffusion of the electric field. In elongated plasmas, the latter reads (Fülöp et al. 2019)

(4.1) $$\begin{eqnarray}\frac{1+\unicode[STIX]{x1D705}^{-2}}{2r}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}r}\left(r\frac{\unicode[STIX]{x2202}E}{\unicode[STIX]{x2202}r}\right)=\unicode[STIX]{x1D707}_{0}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}t}(\unicode[STIX]{x1D70E}_{\Vert }E+n_{\text{RE}}ec),\end{eqnarray}$$

where $E$ is the electric field, $\unicode[STIX]{x1D70E}_{\Vert }$ is the Spitzer conductivity with a neoclassical correction and $\unicode[STIX]{x1D705}$ is the elongation.

The model employed in go has several limitations, e.g. it neglects radial losses due to magnetic perturbations and coupling to the coils (i.e. go has a perfectly conducting wall at radius $r=b$ , which gives an electric-field boundary condition at the plasma edge $r=a$ by matching to the vacuum solution). Therefore, the numerical results are not expected to match the experimental values exactly, and the following examples are shown only as an illustration of the magnitude of the effect expected in an experimentally relevant scenario.

As an example of a scenario where the effect of partially ionized atoms is expected to be important, we consider JET discharge no. 79423, in which a disruption was triggered with injection of $7.4\times 10^{20}$ argon atoms and a runaway-electron plateau of $\simeq 590~\text{kA}$ was observed. The experimental parameters and the details of the discharge are described by Papp et al. (2013). The pre-disruption parameters in the simulation were: major radius $R=3~\text{m}$ ; minor radius $a=0.88~\text{m}$ ; radius of the conducting wall $b=1.3~\text{m}$ ; initial plasma current $I_{\text{p}}=1.93~\text{MA}$ ; magnetic field on axis $B=2~\text{T}$ ; elongation $\unicode[STIX]{x1D705}=1.3$ ; density $n_{\text{e}}(r)=n_{0}(1-1.27\cdot r^{2})^{0.43}$ , with $n_{0}=2.59\times 10^{19}~\text{m}^{-3}$ ; and temperature $T(r)=T_{0}(1-1.03\cdot r^{2})^{2}$ with $T_{0}=2.17~\text{keV}$ and where $r$ is the normalized radial distance from the magnetic axis.

In the go simulations, we only included Dreicer and avalanche runaway sources (no hot-tail generation). For the avalanche growth rate, we use the formula derived in Hesslow et al. (2019), which has been shown to give accurate results compared to numerical solutions of the kinetic equation. The temperature evolution was taken from the experiment, but the post-disruption measurements exhibit a high degree of noise. To correct for this artefact, we set the central temperature to $T_{0}=20~\text{eV}$ after the thermal quench, which gives a current quench time that agrees with the experiment. During the thermal quench, the ionization of the impurities was determined by calculating the density of each charge state for every ion species ( $n_{Z_{i}}^{k},k=0\cdots Z_{i}$ ),

(4.2) $$\begin{eqnarray}\frac{\text{d}n_{Z_{i}}^{k}}{\text{d}t}=n_{\text{e}}(I_{k-1}n_{Z_{i}}^{k-1}-(I_{k}+R_{k})n_{Z_{i}}^{k}+R_{k+1}n_{Z_{i}}^{k+1}),\end{eqnarray}$$

where $I_{k}$ denotes the electron impact ionization rate for the $k$ th charge state and $R_{k}$ is the radiative recombination rate, which were both taken from the ADAS database (Summers 2004). After the disruption, the exact amount of argon was unknown as the assimilation is highly uncertain, and we therefore scanned over the argon density, which was assumed to be uniformly spread throughout the plasma. It is reasonable to assume that in connection with the disruption some additional carbon will penetrate into the plasma, and therefore we present results with both 0 % and 20 % additional carbon.

Figure 4. Plateau runaway current as a function of argon density normalized to the initial electron density. Blue dashed lines are simulations with the Connor & Hastie (1975) formula and solid black lines are simulations with the neural network. Simulations with squares include 20 % additional carbon, $n_{\text{C}}=0.2n_{\text{e,ini}}$ . The experimental value (590 kA) is shown with the dotted horizontal line, and the vertical dotted line shows the value for 100 % assimilation, corresponding to $n_{\text{Ar}}/n_{\text{e,ini}}=0.88$ .

Figure 4 shows the plateau runaway current given by go as a function of argon density, using the Connor & Hastie (1975) analytical formula or the neural network, with or without additional carbon. With the analytical formula, argon densities corresponding to more than 100 % assimilation are required to match the experimental value, whereas lower assimilation is sufficient with the neural network.

Figure 5(a,b) shows the time evolution of the plasma current for a case when the current density nearly matches the experimental value ( $n_{\text{Ar}}/n_{\text{e,ini}}=0.6$ , $n_{\text{C}}=0$ ) comparing the analytical prediction with the neural network. The Dreicer generation mainly occurs in a short interval around 1 ms (although the Dreicer seed is barely visible in figure 5 b). As shown figure 5(c), this time period coincides with the largest values of $E/E_{\text{D}}$ , which is expected as the Dreicer generation rate is highly sensitive to this normalized electric field. Comparing figures 5(a) and 5(b), the neural network predicts a significantly reduced Dreicer seed, which is only partially compensated by the increased avalanche multiplication in the self-consistent electric field in figure 5(c). The result is an order-unity reduction in the plateau runaway current.

Figure 5. Time evolution of the plasma current for the case when the current approximately matches the experimental value ( $n_{\text{Ar}}/n_{\text{e,ini}}=0.6$ , $n_{\text{C}}=0$ ) with (a) the Connor & Hastie (1975) formula and (b) the neural network. The Dreicer current is shown by dotted lines, and dashed lines show the total runaway current (Dreicer $+$ avalanche). In panel (c), the left $y$ -axis shows the induced electric field (Connor–Hastie formula in solid line, neural network in dashed line), and the right $y$ -axis shows the temperature in dash-dotted blue line. Both quantities were evaluated at the magnetic axis. Note the shorter time scale in (c) compared to (a) and (b).

5 Discussion and conclusions

Runaway acceleration of particles occurring in plasmas with strong electric fields has been studied for more than a century, but only recently has it become possible to perform kinetic simulations in complex scenarios. However, coupling such kinetic calculations to a self-consistent simulation of the evolution of the background plasma parameters still poses a significant computational challenge. It is therefore useful to have runaway generation rates that can be used in so-called runaway fluid models, i.e. integrated runaway simulations where generation rates are used to evolve the runaway current instead of solving the full kinetic problem.

In this work, we presented a neural network that determines the steady-state Dreicer runaway generation rate accounting for collisions with partially ionized atoms. The network was trained on a large number of kinetic simulations and gives accurate results for the Dreicer runaway generation rate in plasmas consisting of hydrogen isotopes, neon and argon. This tool can therefore be used for rapid evaluation of such generation rates in any laboratory, space or astrophysical plasma where runaway electrons are produced. In particular, the tool should be valuable in simulations of tokamak disruptions involving massive material injection, in which case partially ionized impurities are expected to play an important role in the dynamics. Together with previously developed generation rates for avalanche, hot-tail, tritium decay and Compton scattering, the Dreicer generation rate neural network offers an improved runaway fluid model.

In their present form, runaway fluid models have some limitations. Most notably, these models typically relate the number density generation rate to the associated current density by approximating the mean parallel runaway velocity $\langle v_{\Vert }\rangle \approx c$ . This assumption is expected to be more accurate in avalanche-dominated scenarios, where the runaway seed population carries a negligible part of the current. In such scenarios, the runaway number density is amplified through the avalanche mechanism while the average speed approaches the speed of light. Conversely, if a significant part of the current is converted into a runaway current through the Dreicer or hot-tail mechanism, the assumption $\langle v_{\Vert }\rangle \approx c$ may significantly overestimate the runaway current. As the mean velocity evolves in time during runaway generation in such scenarios, $\langle v_{\Vert }\rangle$ cannot be determined by steady-state simulations like those performed here, but would require time-dependent modelling including the history of the electric field. Nevertheless, this generation rate tool offers an improvement of the fluid models already without accounting for such effects, especially because runaway generation is generally avalanche dominated in the present and future large tokamaks that carry multi-megaampere plasma currents.

As an illustration, we implemented this model in go (Smith et al. 2006), and presented numerical solutions of the coupled equations of runaway generation and electric-field diffusion in a JET-like disruptive scenario. In this scenario, we observed that the plateau runaway current was significantly reduced when using the neural network instead of the Connor–Hastie formula, which demonstrates the need to account for partially ionized atoms for realistic modelling of Dreicer generation. The neural network presented here can therefore be useful to improve runaway-electron modelling in tokamak simulation codes such as astra (Fable et al. 2016), jorek (Bandaru et al. 2019) and ets (Pokol et al. 2019).