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 Reference Dreicer1959) – until the electric force is balanced by another mechanism such as synchrotron radiation. Runaway acceleration occurs in solar flares (Holman Reference Holman1985), lightning discharges (Dwyer Reference Dwyer2007) and in magnetic fusion devices (Helander, Eriksson & Andersson Reference Helander, Eriksson and Andersson2002; Breizman et al. Reference Breizman, Aleynikov, Hollmann and Lehnen2019). 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. Reference Lehnen, Aleynikova, Aleynikov, Campbell, Drewelow, Eidietis, Gasparyan, Granetz, Gribov and Hartmann2015).
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.
Reference Smith, Helander, Eriksson, Anderson, Lisak and Andersson2006; Gál et al.
Reference Gál, Fehér, Smith, Fülöp and Helander2008; Fehér et al.
Reference Fehér, Smith, Fülöp and Gál2011; Papp et al.
Reference Papp, Fülöp, Fehér, de Vries, Riccardo, Reux, Lehnen, Kiptily, Plyusnin and Alper2013) and the work by Martín-Solís, Loarte & Lehnen (Reference Martín-Solís, Loarte and Lehnen2017).
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. Reference Stahl, Hirvijoki, Decker, Embréus and Fülöp2015), bremsstrahlung (Embréus, Stahl & Fülöp Reference Embréus, Stahl and Fülöp2016) and screening effects in partially ionized plasmas (Hesslow et al. Reference Hesslow, Embréus, Stahl, DuBois, Papp, Newton and Fülöp2017). 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. Reference Hesslow, Embréus, Hoppe, DuBois, Papp, Rahm and Fülöp2018) and the avalanche multiplication factor (Hesslow et al. Reference Hesslow, Embréus, Vallhagen and Fülöp2019).
While the avalanche growth rate calculation (Rosenbluth & Putvinski Reference Rosenbluth and Putvinski1997) was recently generalized to include the effect of partial screening and radiation reaction (Hesslow et al.
Reference Hesslow, Embréus, Vallhagen and Fülöp2019), a similar generalization of the Dreicer generation rate calculation by Connor & Hastie (Reference Connor and Hastie1975) 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 (Reference Martín-Solís, Loarte and Lehnen2015) suggested to replace the Dreicer field
$E_{\text{D}}$
and the effective charge
$Z_{\text{eff}}$
in the Connor & Hastie (Reference Connor and Hastie1975) 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 (Reference Svensson, von Hellermann and König1999), Clayton et al. (Reference Clayton, Tritz, Stutman, Bell, Diallo, Leblanc and Podestà2013), Citrin et al. (Reference Citrin, Breton, Felici, Imbeaux, Aniel, Artaud, Baiocchi, Bourdelle, Camenen and Garcia2015) and Boyer, Kaye & Erickson (Reference Boyer, Kaye and Erickson2019).
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 (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 1 ; this flow rate defines the steady-state Dreicer generation rate,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn1.gif?pub-status=live)
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
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn2.gif?pub-status=live)
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$
.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_fig1g.gif?pub-status=live)
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
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn3.gif?pub-status=live)
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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_tab1.gif?pub-status=live)
Rather than using the full set of impurity ion densities as input to the neural network, we use six derived parameters,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn4.gif?pub-status=live)
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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_fig2g.gif?pub-status=live)
Figure 2. Comparison between normalized runaway generation rates obtained from code and those from the neural network regression.
The neural networkFootnote 2 is designed as a multilayer perceptron described by the following series of matrix multiplications and function applications:
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn5.gif?pub-status=live)
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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_fig3g.gif?pub-status=live)
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)
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn6.gif?pub-status=live)
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}$
),
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn7.gif?pub-status=live)
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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_fig4g.gif?pub-status=live)
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.
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_fig5g.gif?pub-status=live)
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).
Acknowledgements
The authors are grateful to S. Newton for fruitful discussions. This work was supported by the Swedish Research Council (Dnr. 2018-03911) and the Knut and Alice Wallenberg Foundation. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 647121. This work was supported by the EUROfusion – Theory and Advanced Simulation Coordination (E-TASC). This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under grant agreement no. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.
Appendix A. Collision operator used in the Fokker–Planck simulations
The kinetic equation solved by code (Landreman et al. Reference Landreman, Stahl and Fülöp2014; Stahl et al. Reference Stahl, Embréus, Papp, Landreman and Fülöp2016) describes a magnetized, homogeneous plasma,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn8.gif?pub-status=live)
where
$f_{\text{e}}$
is the electron distribution function,
$E$
is the component of the electric field that is antiparallel to the magnetic field
$\boldsymbol{B}$
and
$\unicode[STIX]{x1D709}=\boldsymbol{p}\boldsymbol{\cdot }\boldsymbol{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_{\text{br}}$
(the bremsstrahlung collision operator) and
$\boldsymbol{F}_{\text{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 (Reference Braams and Karney1989) and Pike & Rose (Reference Pike and Rose2014), with corrections for partial screening according to Hesslow et al. (Reference Hesslow, Embréus, Hoppe, DuBois, Papp, Rahm and Fülöp2018),
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn9.gif?pub-status=live)
where the slowing-down, parallel and deflection frequencies are given by, respectively,
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn10.gif?pub-status=live)
Here,
$\bar{p}=p/(m_{\text{e}}c)$
is the normalized momentum,
$\unicode[STIX]{x1D6E9}=T/(m_{\text{e}}c^{2})$
and
$\unicode[STIX]{x1D70F}=4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}^{2}m_{\text{e}}^{2}c^{3}/(n_{\text{e}}e^{4}\ln \unicode[STIX]{x1D6EC}_{0})$
is a relativistic collision time. The energy-dependent Coulomb logarithm is modelled by matching the thermal Coulomb logarithm
$\ln \unicode[STIX]{x1D6EC}_{0}=14.9{-}0.5\ln (n_{\text{e}}[10^{20}\,\text{m}^{-3}])+\ln (T[\text{keV}])$
from Wesson (Reference Wesson2011) and the high-energy formula from Solodov & Betti (Reference Solodov and Betti2008), according to
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn11.gif?pub-status=live)
where
$\bar{p}_{\text{Te}}=\sqrt{2T/(m_{\text{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
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn12.gif?pub-status=live)
with
$\unicode[STIX]{x1D6F9}_{0}=\int _{0}^{\bar{p}}\exp [-\sqrt{1+s^{2}}/\unicode[STIX]{x1D6E9}]/(\sqrt{1+s^{2}})\,\text{d}s$
,
$\unicode[STIX]{x1D6F9}_{1}=\int _{0}^{\bar{p}}\exp [-\sqrt{1+s^{2}}/\unicode[STIX]{x1D6E9}]\,\text{d}s$
and
$K_{2}$
is the second-order modified Bessel function of the second kind. Finally, the partial screening corrections are
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn13.gif?pub-status=live)
where
$Z_{0,i}$
is the charge number,
$Z_{i}$
is the atomic number and
$N_{\text{e},i}=Z_{i}-Z_{0,i}$
is the number of bound electrons of the nucleus for species
$i$
. The length parameter
$\bar{a}_{i}$
is given in Hesslow et al. (Reference Hesslow, Embréus, Hoppe, DuBois, Papp, Rahm and Fülöp2018) and the mean excitation energy
$I_{i}$
in Sauer et al. (Reference Sauer, Oddershede and Sabin2015).
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 (Stahl et al.
Reference Stahl, Landreman, Embréus and Fülöp2017), which uses the collision operator derived by Braams & Karney (Reference Braams and Karney1989). The simulations were performed at temperatures in the range 10 eV–10 keV, and for electric fields
$E/E_{\text{D}}\approx 1\,\%{-}3.5\,\%$
, corresponding to
$E/E_{\text{c}}\approx 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_{\text{c}}\gg 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_{\text{c}}$
decreases with temperature at constant
$E/E_{\text{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.
Appendix B. Calculation of the steady-state runaway generation rate
To rapidly calculate the Dreicer generation rate, code can be used in steady-state mode as described by Landreman et al. (Reference Landreman, Stahl and Fülöp2014). Equation (A 1) has the form
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn14.gif?pub-status=live)
which implies that the steady-state distribution can be obtained by
$f_{\text{e}}=M^{-1}S$
, where
$S$
is a source term that cancels the constant outflow of particles due to the Dreicer mechanism. The steady-state generation rate
$\unicode[STIX]{x1D6FE}$
can then be obtained by integrating equation (A 1) over
$2\unicode[STIX]{x03C0}\int _{p_{b}}^{\infty }p^{2}\,\text{d}p\int _{-1}^{1}(\cdot )\,\text{d}\unicode[STIX]{x1D709}$
, where
$p_{\text{b}}$
is taken so that the source
$S$
vanishes for
$p\geqslant p_{\text{b}}$
. If radiation reaction is neglected, this yields
![](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20191219131049522-0883:S0022377819000874:S0022377819000874_eqn15.gif?pub-status=live)
where
$f_{L}=(2L+1)2^{-1}\int _{-1}^{1}f_{\text{e}}P_{L}(\unicode[STIX]{x1D709})\,\text{d}\unicode[STIX]{x1D709}$
is the
$L$
th Legendre mode of the steady-state distribution. This generation rate is insensitive to the exact shape of
$S$
as long as it is limited to the bulk and does not affect the runaway population (Landreman et al.
Reference Landreman, Stahl and Fülöp2014); accordingly, the source is taken to be
$S=\unicode[STIX]{x1D6FC}\text{e}^{-p/p_{\text{Te}}^{2}}$
, with the magnitude
$\unicode[STIX]{x1D6FC}$
determined from the normalization of
$f_{\text{e}}$
. If (B 2) is numerically well resolved, the generation rate is also independent of
$p_{\text{b}}$
as long as the source
$S$
vanishes for
$p\geqslant p_{\text{b}}$
. Here, we took the average value of
$\unicode[STIX]{x1D6FE}$
over
$p_{\text{b}}/p_{\text{Te}}\in [10,50]$
. The lower limit was raised if there were large variations within the interval, which may occur due to floating point errors when evaluating (B 2) for weak electric fields.