2 Kinetic model

The Dreicer (Reference Dreicer1959, Reference Dreicer1960) 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 separatrixFootnote ¹ ; 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 (Reference Connor and Hastie1975), who extended the results of Kruskal & Bernstein (Reference Kruskal and Bernstein1962) 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 (Reference Connor and Hastie1975) but has been quantified in subsequent work to around $C\approx 1.0$ with our time normalization (Kruskal & Bernstein Reference Kruskal and Bernstein1962; Jayakumar, Fleischmann & Zweben Reference Jayakumar, Fleischmann and Zweben1993), 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. (Reference Kulsrud, Sun, Winsor and Fallon1973) 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. Reference Stahl, Hirvijoki, Decker, Embréus and Fülöp2015). 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. Reference Nilsson, Decker, Peysson, Granetz, Saint-Laurent and Vlainic2015). 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 Reference McDevitt and Tang2019).

The results in this work, including figure 1, were obtained with the linearized Fokker–Planck solver code (Landreman, Stahl & Fülöp Reference Landreman, Stahl and Fülöp2014; Stahl et al. Reference Stahl, Embréus, Papp, Landreman and Fülöp2016), 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. Reference Stahl, Landreman, Embréus and Fülöp2017), 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. Reference Hesslow, Embréus, Hoppe, DuBois, Papp, Rahm and Fülöp2018). 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 Reference Mott and Massey1965). 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. (Reference Martín-Solís, Loarte and Lehnen2017) – 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 (Reference Connor and Hastie1975). 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. (Reference Landreman, Stahl and Fülöp2014) 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 Reference Sauer, Oddershede and Sabin2015), 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. Reference Hesslow, Embréus, Hoppe, DuBois, Papp, Rahm and Fülöp2018), 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 networkFootnote ² 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 Reference Kingma and Ba2014), which is a stochastic gradient descent method. The training of the neural network was implemented in Python using the library PyTorch (Paszke et al. Reference Paszke, Gross, Chintala, Chanan, Yang, Devito, Lin, Desmaison, Antiga and Lerer2017). 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. Reference Smith, Helander, Eriksson, Anderson, Lisak and Andersson2006). 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. Reference Fülöp, Helander, Vallhagen, Embreus, Hesslow, Svensson, Creely, Howard and Rodriguez-Fernandez2019)

(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. (Reference Papp, Fülöp, Fehér, de Vries, Riccardo, Reux, Lehnen, Kiptily, Plyusnin and Alper2013). 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. (Reference Hesslow, Embréus, Vallhagen and Fülöp2019), 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 Reference Summers2004). 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 (Reference Connor and Hastie1975) 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 (Reference Connor and Hastie1975) 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 (Reference Connor and Hastie1975) 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. Reference Smith, Helander, Eriksson, Anderson, Lisak and Andersson2006), 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. Reference Fable, Pautasso, Lehnen, Dux, Bernert and Mlynek2016), jorek (Bandaru et al. Reference Bandaru, Hoelzl, Artola, Papp and Huijsmans2019) and ets (Pokol et al. Reference Pokol, Olasz, Erdos, Papp, Aradi, Hoppe, Johnson, Ferreira, Coster and Peysson2019).