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

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.

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 . † Email address for correspondence: hesslow@chalmers.se 2 L. Hesslow

and others
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 momentumspace 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 γ -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 , bremsstrahlung ) 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 .
While the avalanche growth rate calculation (Rosenbluth & Putvinski 1997) was recently generalized to include the effect of partial screening and radiation reaction , 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 D and the effective charge Z 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 Runaway-electron generation rate using a neural network 3 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).

Kinetic model
The Dreicer (1959Dreicer ( , 1960 mechanism for runaway generation originates from the interplay between collisional energy diffusion and electric-field acceleration. In a timeindependent background plasma, the Dreicer mechanism causes a constant particle flow through any momentum-space boundary beyond the runaway separatrix 1 ; this flow rate defines the steady-state Dreicer generation rate, where n 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 γ (E, T, n e , . . . , t) ≈ γ [E(t), T(t), n e (t), . . .], where E, T, n e , . . . 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 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) 1 Above the second separatrix, where fast particles are decelerated due to for example radiation reaction losses, the particles will slow down again. This eventually creates a bump-on-tail Decker et al. 2016;Guo, McDevitt & Tang 2017), which prevents further runaway growth through the Dreicer mechanism, but the involved energies and time scales makes this effect irrelevant except for near the threshold E c . An exception is with strong synchrotron radiation, temperatures of several keV and an electric field approaching the effective critical electric field; in this case a steady-state Dreicer generation rate cannot clearly be determined and a full kinetic simulation may be necessary. , 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.
Here, E D = n e e 3 ln Λ 0 /(4π 2 0 T) is the Dreicer field, E c = E D T/(m e c 2 ) is the critical electric field, Z eff = i n i Z 2 i /n e is the effective plasma charge and τ ee = 4πε 2 0 m 2 e v 3 Te /(n e e 4 ln Λ 0 ) is the thermal electron collision time, with the thermal speed v Te = √ 2T/m e , and where ln Λ 0 is the Coulomb logarithm evaluated at the thermal speed (in the original work, ln Λ 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 ≈ 1.0 with our time normalization (Kruskal & Bernstein 1962;Jayakumar, Fleischmann & Zweben 1993), and the parameters λ, η and h constitute the relativistic generalization of the generation rate; they approach unity as E/E c → ∞.
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 Λ, 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 1b), 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 c 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 Runaway-electron generation rate using a neural network 5 even at relatively synchrotron-favouring parameters with T = 10 keV, B = 5 T and n D = 10 20 m −3 . In these simulations, radiation reaction mainly affects the generation rate close to the effective critical field, which is consistent with previous work . We also note that Dreicer generation predominantly takes place at low temperatures during disruption scenarios, since the value of E/E D ∼ j/(σ E D ) ∼ 1/ √ 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 E 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 E D . Otherwise, a major part of the thermal electron distribution runs away and a nonlinear Fokker-Planck equation, solved by for example NORSE , 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 Zα, where α ≈ 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 → 0. As the Dreicer generation rate is most sensitive to the dynamics near the critical momentum given by p c (E/E c − 1) −1/2 , the model is only strictly valid for E/E D (%) T/(20 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 eff and E D by expressions involving the increased collision rates evaluated at an approximate value of the critical momentum p 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.

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 ratē γ ≡ τ ee n e γ .  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 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. Ar +n for n = 0 · · · 18, and Ne +m for m = 0 · · · 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 Ar + = 219.4 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.
Rather than using the full set of impurity ion densities as input to the neural network, we use six derived parameters, ln n e , Z eff , n e n tot , Here, Z i and Z 0,i are the atomic number and charge number of species i, respectively, and n tot = 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 e and Z eff ) parameters that naturally appear in the partially screened collision frequencies (n e /n tot and i (n i /n tot )(Z 2 i − Z 2 0,i )), 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γ , was randomly split into a training and validation set with a 4 : 1 ratio.
Runaway-electron generation rate using a neural network 7 FIGURE 2. Comparison between normalized runaway generation rates obtained from CODE and those from the neural network regression.
The neural network 2 is designed as a multilayer perceptron described by the following series of matrix multiplications and function applications: Here, x is the input vector described above, the matrices W (k) describe the weights between the layers and the vectors 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. W 1 ∈ R 20×8 , {W 2 , W 3 , W 4 } ∈ R 20×20 and W 5 ∈ R 1×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 γ , 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 eV and deuterium density n D = 10 20 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 8 L. Hesslow and others (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.

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 RE . It solves the coupled equations for the runaway generation and resistive diffusion of the electric field. In elongated plasmas, the latter reads where E is the electric field, σ is the Spitzer conductivity with a neoclassical correction and κ 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 × 10 20 argon atoms and a runaway-electron plateau of 590 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 m; minor radius a = 0.88 m; radius of the conducting wall b = 1.3 m; initial plasma current I p = 1.93 MA; magnetic field on axis B = 2 T; elongation κ = 1.3; density n e (r) = n 0 (1 − 1.27 · r 2 ) 0.43 , with n 0 = 2.59 × 10 19 m −3 ; and temperature T(r) = T 0 (1 − 1.03 · r 2 ) 2 with T 0 = 2.17 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 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 k Z i , k = 0 · · · Z i ), where I k denotes the electron impact ionization rate for the kth 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 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 Ar /n e,ini = 0.6, n 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 5b). As shown figure 5(c), this time period coincides with the largest values of E/E 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.

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 v ≈ 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 v ≈ c may significantly overestimate the runaway current. As the mean velocity evolves in time during runaway generation in such scenarios, v 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 where f e is the electron distribution function, E is the component of the electric field that is antiparallel to the magnetic field B and ξ = p · B/(pB) is the cosine of the pitch angle. Collisions are modelled by the Fokker-Planck collision operator (as this work focuses on Dreicer generation, a large-angle collision operator describing avalanche generation is not included). Radiation losses are modelled by C br (the bremsstrahlung collision operator) and F syn (the synchrotron radiation reaction force). The Fokker-Planck solver CODE includes several options for the collision operator. In this work we use the test-particle collision operator given by Braams & Karney (1989) and Pike & Rose (2014), with corrections for partial screening according to Hesslow et al. (2018), where the slowing-down, parallel and deflection frequencies are given by, respectively, Here,p = p/(m e c) is the normalized momentum, Θ = T/(m e c 2 ) and τ = 4πε 2 0 m 2 e c 3 / (n e e 4 ln Λ 0 ) is a relativistic collision time. The energy-dependent Coulomb logarithm is modelled by matching the thermal Coulomb logarithm ln Λ 0 = 14.9-0.5 ln(n e [10 20 m −3 ]) + ln(T[keV]) from Wesson (2011)  wherep Te = 2T/(m e c 2 ) is the normalized thermal momentum and the parameter k = 5 is chosen to give a smooth transition. The contributions from the fully ionized plasma are captured by the functions Ψ s (p, Θ) = γ 2 Ψ 1 − ΘΨ 0 + (Θγ − 1)pe −γ /Θ p 2 K 2 (1/Θ) , Ψ D (p, Θ) = (p 2 γ 2 + Θ 2 )Ψ 0 + Θ(2p 4 − 1)Ψ 1 + γ Θ[1 + Θ(2p 2 − 1)]pe −γ /Θ 2γp 3 K 2 (1/Θ) , N 2 e,i (pā i ) 3/2 (pā i ) 3/2 + 1 , Runaway-electron generation rate using a neural network 13 where Z 0,i is the charge number, Z i is the atomic number and N e,i = Z i − Z 0,i is the number of bound electrons of the nucleus for species i. The length parameterā i is given in Hesslow et al. (2018) and the mean excitation energy I i in Sauer et al. (2015). The collision operator given above is a linearized relativistic test-particle operator, i.e. the field-particle part of the operator is neglected. To assess the discrepancy due to the combined effect of the field-particle operator and nonlinear effects, we compared the runaway generation rates computed with CODE, using the test-particle operator given above, and simulations with the fully nonlinear kinetic equation solver NORSE , which uses the collision operator derived by Braams & Karney (1989). The simulations were performed at temperatures in the range 10 eV-10 keV, and for electric fields E/E D ≈ 1 %-3.5 %, corresponding to E/E c ≈ 2-2700. In all simulations, we find agreement between CODE and NORSE simulations within a factor of two, which is enough to capture the order of magnitude of the Dreicer seed. Below 1 keV, where E/E c 1, the two tools agreed within 15 %. In contrast, when the temperature increases beyond 1 keV, the differences were up to 50 %. This is because E/E c decreases with temperature at constant E/E D , and since the Dreicer generation rate is dramatically sensitive to the electric field at near-critical values, even a negligible difference in the electric field can lead to a significant difference in the generation rate. In other words, the observed differences may be smaller than any reasonable uncertainty in the electric field. Given the orders-of-magnitude variation of the Dreicer generation rate with electric field, the errors here are deemed acceptable and sufficient to capture the magnitude of the Dreicer seed.