2. Zero-dimensional fluid model for runaway electron generation

The density of the runaway electron population $n_{\mathrm{r}}$ is determined by the sum of the primary (seed) and secondary (avalanche) generation mechanisms as

(2.1) \begin{equation} \frac {\mathrm{d}n_{\mathrm{r}}}{\mathrm{d}t} = \sum _{\mathrm{seeds}} \gamma _{\mathrm{seed}} + \varGamma _{\mathrm{ava}} n_{\mathrm{r}}, \end{equation}

with the total seed generation rate (in principle consisting of Dreicer, hot-tail, tritium decay and Compton scattering sources) $\sum _{\mathrm{seeds}}\gamma _{\mathrm{seed}}$ and the avalanche growth rate $\varGamma _{\mathrm{ava}}$ . Here, it is assumed that no runaway electrons are present in the plasma at the time of the disruption onset (at $t = 0$ ). All runaway electron generation mechanisms depend on the parallel electric field $E_\parallel$ , which is inductive and is given by

(2.2) \begin{equation} E_{\parallel } = E_{\parallel 0} - \frac {L}{2 \pi R_0} \frac {\mathrm{d}I_{\mathrm{p}}}{\mathrm{d}t}, \end{equation}

where

(2.3) \begin{equation} L \approx \mu _0 R_0\left [ \ln {\left ( \frac {8R_0}{a} \right )} - 2\right ] \end{equation}

is the self-inductance of the plasma in the large aspect ratio limit, with the vacuum permeability $\mu _0$ , the plasma major radius $R_0$ and the plasma minor radius $a$ ; $E_{\parallel 0}$ is the contribution to the electric field due to the external transformer. The plasma current $I_{\mathrm{p}}$ consists of two contributions, namely the Ohmic current and the runaway electron current, and can be written as

(2.4) \begin{equation} I_{\mathrm{p}} = \sigma E_\parallel A + n_{\mathrm{r}} e c A, \end{equation}

where $\sigma$ is the Spitzer conductivity, and $A$ is the (effective) cross-sectional area of the plasma. The constants $e$ and $c$ denote the electron charge and the speed of light, respectively. Following the procedure proposed by Helander et al. (Reference Helander, Eriksson and Andersson2002), (2.1), (2.2) and (2.4) can be renormalised by introducing the dimensionless quantities

(2.5) \begin{align} n & = \frac {n_{\mathrm{r}} e c}{j_0}, \end{align}

(2.6) \begin{align} E & = \frac {E_\parallel }{E_{\mathrm{c}}}, \end{align}

(2.7) \begin{align} t' & = \frac {t}{\sqrt {5 + Z_{\mathrm{eff}}} \tau _{\mathrm{c}} \ln {\Lambda _{\mathrm{c}}}}, \end{align}

(2.8) \begin{align} s & = \frac {\sigma E_{\mathrm{c}}}{j_0}, \end{align}

where $\tau _{\mathrm{c}} = 4 \pi \varepsilon _0^2 m_{\mathrm{e}}^2 c^3 / (n_{\mathrm{e}} e^4 \ln {\Lambda _{\mathrm{c}}} )$ is the relativistic electron collision time, $m_{\mathrm{e}}$ is the electron mass, $\varepsilon _0$ is the vacuum permittivity, $Z_{\mathrm{eff}}$ is the effective charge, $E_{\mathrm{c}} = m_{\mathrm{e}} c / (e \tau _{\mathrm{c}})$ is the critical electric field introduced by Connor & Hastie (Reference Connor and Hastie1975), $\ln {\Lambda _{\mathrm{c}}} = \ln {\Lambda _{\mathrm{th}}} + 0.5\ln {(m_{\mathrm{e}} c^2 /T)} \approx 14.6 + 0.5 \ln {(T_{\mathrm{eV}}/n_{{\mathrm{e}} 20})}$ is the relativistic Coulomb logarithm, $\ln {\Lambda _{\mathrm{th}}} = 14.9 - 0.5\ln {n_{{\mathrm{e}} 20}} + \ln {(T_{\mathrm{eV}}/1000)}$ is the thermal electron–electron Coulomb logarithm, $T_{\mathrm{eV}}$ is the electron temperature in units of electronvolts, $n_{{\mathrm{e}} 20}$ is the electron density in units of $10^{20}\: \mathrm{m}^{-3}$ , and $j_0 = I_{\mathrm{p}0}/A$ and $I_{\mathrm{p0}}$ are the pre-disruption current density and current, respectively. The resulting dimensionless differential equations become

(2.9) \begin{align} \frac {\mathrm{d}n}{\mathrm{d}t'} &= \sum _{\mathrm{seeds}} \bar {\gamma }_{\mathrm{seed}} + \bar {\varGamma }_{\mathrm{ava}} n, \end{align}

(2.10) \begin{align} \frac {\mathrm{d} }{\mathrm{d}t'} \left ( sE + n \right ) &= \frac {E_0 - E}{\alpha } , \end{align}

where

(2.11) \begin{align} E_0 & = \frac {E_{\parallel 0}}{E_{\mathrm{c}}}, \end{align}

(2.12) \begin{align} \alpha & = \frac {2}{\sqrt {5 + Z_{\mathrm{eff}}}}\frac {L}{\mu _0 R_0} \frac {I_{\mathrm{p}0}}{I_{\mathrm{A}} \ln {\Lambda _{\mathrm{c}}}}, \end{align}

(2.13) \begin{align} \bar {\gamma }_{\mathrm{seed}} & = \sqrt {5 + Z_{\mathrm{eff}}} \tau _{\mathrm{c}} \ln {\Lambda _{\mathrm{c}}} \frac {e c}{j_0} \gamma _{\mathrm{seed}}, \end{align}

(2.14) \begin{align} \bar {\varGamma }_{\mathrm{ava}} & = \sqrt {5 + Z_{\mathrm{eff}}} \tau _{\mathrm{c}} \ln {\Lambda _{\mathrm{c}}} \varGamma _{\mathrm{ava}}. \end{align}

Here, $I_{\mathrm{A}} = 4 \pi m_{\mathrm{e}} c /(\mu _0 e) \approx {0.017}\,\mathrm{MA}$ is the Alfvén current. Note that, using this normalisation, $n$ is equivalent to the conversion rate of the pre-disruption Ohmic current into a runaway electron current. Dividing (2.10) by (2.9), the system is reduced to the single nonlinear ordinary differential equation

(2.15) \begin{equation} s\frac {\mathrm{d}E}{\mathrm{d}n} = -1 - \frac {E - E_0}{\alpha \left [\sum _{\mathrm{seeds}} \bar {\gamma }_{\mathrm{seed}} + \bar {\varGamma }_{\mathrm{ava}}n\right ]}. \end{equation}

Assuming that the runaway electron generation is initially dominated by primary generation mechanisms, the avalanche source term can be neglected, leading to

(2.16) \begin{equation} s\frac {\mathrm{d}E}{\mathrm{d}n} = -1 - \frac {E - E_0}{\alpha \sum _{\mathrm{seeds}} \bar {\gamma }_{\mathrm{seed}} }. \end{equation}

At early times in the disruption, the electric field $E$ exceeds $E_0$ by a large margin. For example, in a high performance ITER plasma without material injection, $E_0 \approx 0.05$ , whereas $E$ at the very least is larger than unity. In addition, $\varSigma _{\mathrm{seeds}}\bar {\gamma }_{\mathrm{seed}}$ is assumed to be much smaller than $E/\alpha$ , i.e. the second term on the right hand side of equation (2.16) dominates. Using these assumptions, equation (2.16) reduces to

(2.17) \begin{equation} s\frac {\mathrm{d}E}{\mathrm{d}n} = - \frac {E}{ \alpha \sum _{\mathrm{seeds}} \bar {\gamma }_{\mathrm{seeds}} }, \end{equation}

which can be integrated to yield the seed density

(2.18) \begin{equation} n_{\mathrm{seed}} = s \alpha \int _{\bar {E}_{\mathrm{c}}^{\mathrm{eff}}}^{E_1} \frac {\sum _{\mathrm{seeds}}\bar {\gamma }_{\mathrm{seed}}}{E}\,{\mathrm{d}} E. \end{equation}

The upper integration limit is given by $E_1$ , which denotes the normalised electric field immediately after the thermal quench. The lower integration limit, $\bar {E}_{\mathrm{c}}^{\mathrm{eff}}$ , is the effective critical electric field (Hesslow et al. Reference Hesslow, Embréus, Wilkie, Papp and Fülöp2018b ) in units of $E_{\mathrm{c}}$ . Here, $\bar {E}_{\mathrm{c}}^{\mathrm{eff}}$ includes radiation forces and the enhanced collision frequency due to runaway electrons probing the internal structure of partially ionised ions.

The avalanche gain factor can, in principle, be estimated in a similar manner by assuming that after some short time (compared with the entire avalanche phase), the runaway electron density is sufficiently large for secondary generation to dominate, resulting in the avalanche gain factor

(2.19) \begin{equation} n = n_{\mathrm{seed}} \exp ({N_{\mathrm{ava}}}), \quad \quad N_{\mathrm{ava}} = s\alpha \int _{\bar {E}_{\mathrm{c}}^{\mathrm{eff}}}^{E_1} \frac {\bar {\varGamma }_{\mathrm{ava}}}{E} \, {\mathrm{d}} E. \end{equation}

However, as the final runaway electron density is exponentially sensitive to $N_{\mathrm{ava}}$ , it is important to take into account the radial diffusion of the parallel electric field. This means that $N_{\mathrm{ava}}$ has to be modified by a factor $2 a_{\mathrm{wall}}^2/(x_1^2 a^2)$ (see Hesslow et al. Reference Hesslow, Embréus, Vallhagen and Fülöp2019), where $a_{\mathrm{wall}}$ is the minor radius of the wall and $x_1 \approx 2.4$ is the first zero of the zeroth-order Bessel function of the first kind, $J_0(x)$ . Here, $N_{\mathrm{ava}}$ takes the form (Hesslow et al. Reference Hesslow, Embréus, Vallhagen and Fülöp2019)

(2.20) \begin{equation} N_{\mathrm{ava}} = \tau _{\mathrm{CQ}} \int _{\bar {E}_{\mathrm{c}}^{\mathrm{eff}}}^{E_1} \frac {\varGamma _{\mathrm{ava}}}{E} \, {\mathrm{d}} E, \end{equation}

where $\tau _{\mathrm{CQ}}$ is the characteristic time scale of the current decay in the cylindrical limit and is given by

(2.21) \begin{equation} \tau _{\mathrm{CQ}} = \frac {L}{\mu _0 R_0}\frac {\sigma \mu _0 a_{\mathrm{wall}}^2}{x_1^2}. \end{equation}

The exponent $N_{\mathrm{ava}}$ is thus sensitive to the parameter $a_{\mathrm{wall}}$ , and the most accurate choice for $a_{\mathrm{wall}}$ is not necessarily the minor radius of the first wall. Therefore, in this context, the wall refers to the toroidally continuous conducting structure closest to the plasma, which is normally at a larger radius than the first wall of the vacuum chamber. The parameter $a_{\mathrm{wall}}$ is made sufficiently large to contain a representative poloidal magnetic energy that is released during the disruption. Therefore, the values $\ a_{\mathrm{wall}} = {2.833}\,\mathrm{}$ and $ {0.621}\,\mathrm{m}$ are used in the ITER and SPARC simulations, respectively (Pusztai et al. Reference Pusztai, Ekmark, Bergström, Halldestam, Jansson, Hoppe, Vallhagen and Fülöp2023; Ekmark et al. Reference Ekmark, Hoppe, Tinguely, Sweeney, Fülöp and Pusztai2025).

For an exponentially increasing runaway population, the assumption that the first term on the right-hand side of (2.15) is negligible eventually becomes invalid, and the runaway electron current reaches unphysical values. However, the purpose of this model is not to determine the exact value of the runaway current, but rather to determine whether the exponential growth is sufficiently strong for the runaway current to become comparable to the pre-disruption Ohmic current.

3. Avalanche generation and partial screening

Runaway electron mitigation schemes usually involve massive material injection, where a mixture of deuterium and noble gases (primarily neon or argon) is injected into the plasma, either through gas puffing (Reux et al. Reference Reux2015; Pautasso et al. Reference Pautasso2020) or in the form of cryogenic pellets (Reux et al. Reference Reux2022; Halldestam et al. Reference Halldestam, Heinrich, Papp, Hoppe, Hoelzl, Pusztai, Vallhagen, Fischer and Jenko2025; Patel et al. Reference Patel2025). The main objective of the material injection is to enhance line radiation to reduce the localised thermal loads and electromagnetic forces on the device (Bandyopadhyay et al. Reference Bandyopadhyay2025), as well as to increase the collisional drag, which leads to a larger effective critical electric field and critical momentum. Not only do free electrons contribute to the enhanced fraction force, but bound electrons as well, since energetic electrons may probe the internal structure of an ion, so that the nucleus is only partially screened by its bound electrons (Hesslow et al. Reference Hesslow, Embréus, Wilkie, Papp and Fülöp2018b ). Accounting for this effect, the critical momentum becomes (Hoppe, Embréus & Fülöp Reference Hoppe, Embréus and Fülöp2021)

(3.1) \begin{equation} p_{\mathrm{c}}^2 = \frac {\sqrt {4 \bar {\nu }_{\mathrm{s}}(p_\star )^2 + \bar {\nu }_{\mathrm{s}}(p_\star ) \bar {\nu }_{\mathrm{D}}(p_\star ) }}{E - \bar {E}_{\mathrm{c}}^{\mathrm{eff}}}, \end{equation}

where $\bar {E}_{\mathrm{c}}^{\mathrm{eff}}$ is evaluated using expressions (23) and (24) in Hesslow et al. (Reference Hesslow, Embréus, Wilkie, Papp and Fülöp2018b ) and $\bar {\nu }_{\mathrm{D}}$ and $\bar {\nu }_{\mathrm{s}}$ denote the normalised deflection and slowing-down frequencies, respectively, and are defined as (Hesslow et al. Reference Hesslow, Embréus, Hoppe, DuBois, Papp, Rahm and Fülöp2018a )

(3.2) \begin{align} \bar {\nu }_{\mathrm{D}}(p) & = 1 + Z_{\mathrm{eff}} \frac {\ln {\Lambda ^{{\mathrm{e}} \mathrm{i}}}}{\ln {\Lambda ^{{\mathrm{e}}{\mathrm{e}}}}} \sum _j \frac {n_j}{n_{\mathrm{e}}} g_j(p), \end{align}

(3.3) \begin{align} \bar {\nu }_{\mathrm{s}}(p) & = 1 + \frac {1}{\ln {\Lambda ^{{\mathrm{e}}{\mathrm{e}}}}}\sum _j \frac {n_j}{n_{\mathrm{e}}} N_{{\mathrm{e}},j} \left [ \frac {1}{k} \ln {\left(1 + h_j^k\right)} - \beta ^2 \right ]\!, \end{align}

where the functions $g_j$ and $h_j$ are defined as

(3.4) \begin{align} g_j(p) & = \frac {2}{3} \left ( Z_j^2 - Z_{0,j}^2 \right ) \ln {\left [(p \bar {a}_j)^{3/2} + 1\right ]} - \frac {2}{3} \frac {N_{{\mathrm{e}},j}^2 (p \bar {a}_j)^{3/2}}{(p \bar {a}_j)^{3/2} + 1}, \end{align}

(3.5) \begin{align} h_j(p) & = p\sqrt {\gamma - 1} \frac {m_{\mathrm{e}} c^2}{I_j}. \end{align}

In the expressions above, $Z_j$ denotes the atomic number of ion species $j$ , $Z_{0,j}$ denotes the charge number and $N_{{\mathrm{e}},j} = Z_j - Z_{0,j}$ is the number of bound electrons. Furthermore, $p$ is the electron momentum in units of $m_{\mathrm{e}} c$ , $\beta$ is the electron velocity in units of $c$ and $\gamma$ is the Lorentz factor. The atomic data $\bar {a}_j$ and $I_j$ denote the normalised effective length scale and the mean excitation energy, respectively, and are specific to each ion and charge state. They have been tabulated by Hesslow et al. (Reference Hesslow, Embréus, Hoppe, DuBois, Papp, Rahm and Fülöp2018a ) and Sauer, Sabin & Oddershede (Reference Sauer, Sabin and Oddershede2018), respectively. The ad hoc parameter $k = 5$ is an interpolation parameter between the low and high energy regimes. The electron–electron and electron–ion Coulomb logarithms, $\ln {\Lambda ^{{\mathrm{e}} {\mathrm{e}}}}$ and $\ln {\Lambda ^{{\mathrm{e}} \mathrm{i}}}$ , are evaluated as

(3.6) \begin{align} \ln {\Lambda ^{{\mathrm{e}} {\mathrm{e}}}} & = \ln {\Lambda _{\mathrm{th}}} + \frac {1}{k}\ln {\left [ 1 + \left (\frac {2(\gamma - 1)}{p_{\mathrm{th}}^2}\right )^{k/2} \right ]}, \end{align}

(3.7) \begin{align} \ln {\Lambda ^{{\mathrm{e}} \mathrm{i}}} & = \ln {\Lambda _{\mathrm{th}}} + \frac {1}{k}\ln {\left [ 1 + \left (\frac {2p}{p_{\mathrm{th}}} \right )^k\right ]}, \end{align}

where $p_{\mathrm{th}} = \sqrt {2 T_{\mathrm{e}}/(m_{\mathrm{e}} c^2)}$ is the normalised thermal momentum. The normalised slowing-down and deflection frequencies are evaluated at $p = p_\star$ , which is defined implicitly as

(3.8) \begin{equation} p_\star ^2 = \frac {\sqrt {\bar {\nu }_{\mathrm{s}}(p_\star ) \bar {\nu }_{\mathrm{D}}(p_\star )}}{E}, \end{equation}

and has to be evaluated numerically. In a fully ionised plasma (or in the low energy limit), $\bar {\nu }_{\mathrm{s}} \rightarrow 1$ and $\bar {\nu }_{\mathrm{D}} \rightarrow 1 + Z_{\mathrm{eff}}$ . Implementing the model for partially screened ion nuclei derived by Hesslow et al. (Reference Hesslow, Embréus, Hoppe, DuBois, Papp, Rahm and Fülöp2018a ) would make the otherwise analytical model a semi-analytical model, in the sense that some simple numerical evaluations would be required. For many applications, this is not an obstacle. For this reason, two versions of the criterion will be presented in § 5, namely an analytical criterion that consists solely of analytical expressions, and a semi-analytical criterion of a similar nature but with a more detailed form of the seed densities and the avalanche gain factor. In the analytical model, the critical momentum is evaluated in its completely screened limit and reduces to $p_{\mathrm{c}}^2 = \sqrt {5 + Z_{\mathrm{eff}}}/(E - 1)$ .

Massive material injection not only increases collisional drag but also provides more target electrons during the avalanche, which could enhance the avalanche gain factor by many orders of magnitude if the composition of the injected material is chosen poorly. A model that captures this behaviour was developed by Hesslow et al. (Reference Hesslow, Embréus, Vallhagen and Fülöp2019) and later modified by Hoppe et al. (Reference Hoppe, Embréus and Fülöp2021) to improve the accuracy of the avalanche generation rate in nearly neutral low- $Z$ plasmas. The avalanche generation rate is given by

(3.9) \begin{equation} \varGamma _{\mathrm{ava}} = \frac {1}{\tau _{\mathrm{c}} \ln {\Lambda _{\mathrm{c}}}}\frac {n_{\mathrm{e}}^{\mathrm{tot}}}{n_{\mathrm{e}}} \frac {E - \bar {E}_{\mathrm{c}}^{\mathrm{eff}}}{\sqrt {4 \bar {\nu }_{\mathrm{s}}(p_\star )^2 + \bar {\nu }_{\mathrm{s}}(p_\star ) \bar {\nu }_{\mathrm{D}}(p_\star ) }}, \end{equation}

where $n_{\mathrm{e}}^{\mathrm{tot}}$ is the total number of electrons (free and bound) in the plasma. In the limit of a fully ionised plasma, the avalanche generation rate derived by Rosenbluth & Putvinski (Reference Rosenbluth and Putvinski1997) is recovered. To allow for analytical treatment, the avalanche generation rate of Rosenbluth and Putvinski can be adjusted slightly to approximate the effect of partially ionised ions by including the bound target electrons and by modifying the critical electric field (Putvinski et al. Reference Putvinski, Fujisawa, Post, Putvinskaya, Rosenbluth and Wesley1997) and the effective charge (Parks, Rosenbluth & Putvinski Reference Parks, Rosenbluth and Putvinski1999). To be consistent with the notation in Hesslow et al. (Reference Hesslow, Embréus, Hoppe, DuBois, Papp, Rahm and Fülöp2018a ), these will be denoted $E_{\mathrm{c}}^{\mathrm{RP}}$ and $Z_{\mathrm{eff}}^{\mathrm{RP}}$ , respectively. The approximate avalanche generation rate becomes

(3.10) \begin{equation} \varGamma _{\mathrm{ava}}^{\mathrm{RP}} = \frac {1}{\tau _{\mathrm{c}} \ln {\Lambda _{\mathrm{c}}}}\frac {n_{\mathrm{e}}^{\mathrm{tot}}}{n_{\mathrm{e}}} \frac {E/\bar {E}_{\mathrm{c}}^{\mathrm{RP}} - 1}{\sqrt {5 + Z_{\mathrm{eff}}^{\mathrm{RP}} }}, \end{equation}

where

(3.11) \begin{align} E_{\mathrm{c}}^{\mathrm{RP}} & = \frac {1}{2}\left ( 1 + \frac {n_{\mathrm{e}}^{\mathrm{tot}}}{n_{\mathrm{e}}}\right )E_{\mathrm{c}} \equiv \bar {E}_{\mathrm{c}}^{\mathrm{RP}} E_{\mathrm{c}}, \end{align}

(3.12) \begin{align} Z_{\mathrm{eff}}^{\mathrm{RP}} & = \sum _{\substack {j \: \text{part.} \\ \text{ionized}}} \frac {n_j}{n_{\mathrm{e}}} \frac {Z_j^2}{2} + \sum _{\substack {j \: \text{fully} \\ \text{ionized}}} \frac {n_j}{n_{\mathrm{e}}} Z_j^2. \end{align}

The exponent in the avalanche gain factor becomes

(3.13) \begin{equation} N_{\mathrm{ava}}^{\mathrm{RP}} = \frac {1}{\ln {\Lambda _{\mathrm{c}}} \sqrt {5 + Z_{\mathrm{eff}}^{\mathrm{RP}} }} \frac {\tau _{\mathrm{CQ}}}{\tau _{\mathrm{c}} }\frac {n_{\mathrm{e}}^{\mathrm{tot}}}{n_{\mathrm{e}}} \left [ \frac {E_1}{\bar {E}_{\mathrm{c}}^{\mathrm{RP}}} - \ln {\left (\frac {E_1}{\bar {E}_{\mathrm{c}}^{\mathrm{RP}}}\right )} - 1\right ] {\textrm {H}}{\left (E_1 - \bar {E}_{\mathrm{c}}^{\mathrm{RP}}\right )}, \end{equation}

where the Heaviside function $\textrm {H}$ reflects that avalanche generation occurs only when $E_1 \gt \bar {E}_{\mathrm{c}}^{\mathrm{RP}}$ .

In this work, the ion charge state distribution and electron temperature are determined assuming that the plasma is in a collisional–radiative equilibrium where the Ohmic heating is balanced by the line radiation and bremsstrahlung (see Martín-Solís et al. Reference Martín-Solís, Loarte and Lehnen2017; Vallhagen et al. Reference Vallhagen, Embréus, Pusztai, Hesslow and Fülöp2020) according to

(3.14) \begin{align} \frac {j_0^2}{\sigma (T_{\mathrm{e}},Z_{\mathrm{eff}})} & = \sum _{i,l} n_{\mathrm{e}} n_i^l L_i^l (T_{\mathrm{e}},n_{\mathrm{e}}), \end{align}

(3.15) \begin{align} n_i^{l} & = n_i\left ( 1+ \sum _{j=0}^{l-1}\prod _{k=j+1}^{l}\frac {R_i^{k}}{I_i^{k-1}} +\sum _{j=l+1}^Z\prod _{k=l}^{j-1}\frac {I_i^{k}}{R_i^{k+1}} \right )^{-1}, \end{align}

where $n_i^l$ is the density of ion species $i$ at charge state $l$ , $L_i^l(T_{\mathrm{e}}, n_{\mathrm{e}})$ is the combined rate coefficient for bremsstrahlung and line radiation from excitation and recombination, $I_i^{l} = I_i^{l}(T_{\mathrm{e}},n_{\mathrm{e}})$ denotes the ionisation rate of ion species $i$ from charge state $l$ to charge state $l + 1$ , $R_i^{l} = R_i^{l}(T_{\mathrm{e}},n_{\mathrm{e}})$ denotes the recombination rate of ion species $i$ from charge state $l$ to charge state $l - 1$ and $n_i$ is the total density of ion species $i$ . The ion charge state distribution (and thereby $n_{\mathrm{e}}$ and $Z_{\mathrm{eff}}$ ) and the electron temperature are initially evaluated self-consistently for the given plasma composition and pre-disruption current density. The temperature is then kept constant. Ionisation, recombination and radiation rates are obtained from the Atomic Data and Analysis Structure (ADAS) database.Footnote ¹ Determining the temperature from a collisional–radiative equilibrium – along with neglecting time-dependent impurity sources – implies that the detailed time evolution of the atomic physics is not accounted for, which can affect the accuracy of both the temperature evolution and the runaway electron generation process (Vallhagen et al. Reference Vallhagen, Embréus, Pusztai, Hesslow and Fülöp2020). This approximation is nevertheless required to maintain the semi-analytical tractability of the problem, as it avoids explicit time integration of the rate equations and provides a level of accuracy sufficient for the purposes of the model.

Figure 1.

Equilibrium electron temperature evaluated using (3.14)–(3.15) for (a) ITER and (b) SPARC. The effective charge $Z_{\mathrm{eff}}$ evaluated using a charge state distribution consistent with the temperature for (c) ITER and (d) SPARC. The quantities on the axes denote injected deuterium and neon densities. As the colour scale of $T_{\mathrm{e}}$ is saturated to highlight interesting ranges, additional contours are added outside these ranges. Note also, that in particular the neon concentration ranges are different in ITER and SPARC.

Equilibrium temperatures and the associated effective charge $Z_{\mathrm{eff}}$ for a range of deuterium and neon concentrations are shown in figure 1. The current density is set to $j_0 = I_{\mathrm{p0}}/(\pi a^2)$ , which becomes ${1.19}$ and ${10.0}\,\mathrm{MA m^{-2}}$ in ITER and SPARC, respectively. The value of $Z_{\mathrm{eff}}$ can be used as a measure of the degree of ionisation of the injected neon and is shown in figures 1(c) and 1(d). Although neon can have a large degree of ionisation in the hot lower left quadrant, it is too diluted to contribute substantially to $Z_{\mathrm{eff}}$ . This is primarily true for ITER, where neon densities down to ${1\times 10^{16}}\,\mathrm{m^{-3}}$ are considered. As the neon concentration increases, more power is radiated away, and the temperature drops, leading to a lower degree of ionisation. In SPARC (figure 1 d), $Z_{\mathrm{eff}}$ peaks around $n_{\mathrm{Ne}} = {1\times 10^{19}}\,\mathrm{m^{-3}}$ at $Z_{\mathrm{eff}} \approx 2.9$ , and this trade-off between plasma temperature and relative neon concentration is clearly visible. In ITER (figure 1 c), however, due to the lower current density, the temperature is lower, and $Z_{\mathrm{eff}}$ peaks in the top left corner at a value of $Z_{\mathrm{eff}} \approx 1.5$ .

4. Nuclear sources

In activated devices, seed runaway electrons can be generated from nuclear sources such as tritium beta decay and Compton scattering of gamma photons originating from the irradiated plasma-facing components. These sources can be significant, in the sense that both sources can (under certain circumstances, see § 5) produce enough seed electrons for them to avalanche into mega-ampere runaway electron currents. In this section, the nuclear seed generation rates are presented, and it is shown how they can be integrated to be compatible with an analytical criterion for significant runaway electron generation.

4.1. Tritium decay

Tritium has a half-life of $\tau _{\mathrm{T}} = {4500\pm 8}\,\mathrm{d}$ and decays as

(4.1) \begin{equation} \mathrm{T} \rightarrow {}^3_2\mathrm{He} + {\mathrm{e}}^- + \bar {\nu }_{\mathrm{e}}+ {18.6}\,\mathrm{k eV}. \end{equation}

On average, the electron receives an energy of 5.7 keV, but all energies from 0 to $W_{\mathrm{max}} = {18.6}\,\mathrm{k eV}$ are possible. The production rate of beta electrons is given by

(4.2) \begin{equation} \frac {\mathrm{d}n_\beta }{\mathrm{d}t} = \ln {2} \frac {n_{\mathrm{T}}}{\tau _{\mathrm{T}}}, \end{equation}

where $n_{\mathrm{T}}$ is the tritium density. In a post-disruption plasma with temperatures up to 100 eV, beta electrons can be considered relatively energetic, and some fraction $F_\beta (W_{\mathrm{c}})$ of these will be born above the critical energy $W_{\mathrm{c}} = m_{\mathrm{e}} c^2 \left(\sqrt {p_{\mathrm{c}}^2 + 1} - 1 \right)$ . The runaway electron generation due to tritium decay is then given by (Martín-Solís et al. Reference Martín-Solís, Loarte and Lehnen2017; Fülöp et al. Reference Fülöp, Helander, Vallhagen, Embréus, Hesslow, Svensson, Creely, Howard and Rodriguez-Fernandez2020)

(4.3) \begin{equation} \gamma _{\mathrm{T}}(W_{\mathrm{c}}) = \ln {2} \frac {n_{\mathrm{T}}}{\tau _{\mathrm{T}}} F_\beta (W_{\mathrm{c}}), \end{equation}

where $F_\beta (W_{\mathrm{c}})$ is given by

(4.4) \begin{equation} F_\beta (W_{\mathrm{c}}) = \frac {\int _{W_{\mathrm{c}}}^{W_{\mathrm{max}}} F(p,2)p\mathcal{W}(W_{\mathrm{max}} - W)^2{\mathrm{d}} W}{\int _{0}^{W_{\mathrm{max}}} F(p,2)p\mathcal{W}(W_{\mathrm{max}} - W)^2{\mathrm{d}} W}, \end{equation}

where, $F(p,Z)$ is the Fermi function, and $\mathcal{W} = m_{\mathrm{e}} c^2 \gamma$ is the total energy of the electron. Neglecting the interaction between the positive nucleus and the beta electron, the fraction of beta electrons born above the critical energy is given by (Martin & Shaw Reference Martin and Shaw2019; Fülöp et al. Reference Fülöp, Helander, Vallhagen, Embréus, Hesslow, Svensson, Creely, Howard and Rodriguez-Fernandez2020)

(4.5) \begin{align} F_\beta (W_{\mathrm{c}}) & \approx \frac {\int _{W_{\mathrm{c}}}^{W_{\mathrm{max}}} p(W_{\mathrm{max}} - W)^2{\mathrm{d}} W}{\int _{0}^{W_{\mathrm{max}}} p(W_{\mathrm{max}} - W)^2{\mathrm{d}} W} {\textrm {H}}{(W_{\mathrm{c}}-W_{\mathrm{max}})} \end{align}

(4.6) \begin{align} & = \left (1 - \frac {35}{8} w^{3/2} + \frac {21}{4} w^{5/2} - \frac {15}{8} w^{7/2}\right ) {\textrm {H}}{(w-1)}, \end{align}

where $w = W_{\mathrm{c}}/W_{\mathrm{max}}$ . $F_\beta (W_{\mathrm{c}})$ is plotted against the critical energy in figure 2(a).

Figure 2.

Fraction of $\beta$ -electrons born in the runaway regime, $F_\beta$ , (blue) and normalised Compton cross-section averaged over the gamma $\gamma$ -spectrum, $\bar {\sigma }_{\mathrm{eff}}$ , plotted against critical energy $W_{\mathrm{c}}$ (a). The effective cross-section is plotted for DT plasmas in both ITER (orange) and SPARC (green), and is obtained by averaging the cross-section (black) over the corresponding Compton spectrum $\varGamma _\gamma$ (b). The cross-section is plotted for critical energies 0, 1, 10 and 100 keV (solid, dashed, dot-dashed and dotted curves, respectively).

If the critical energy is larger than $W_{\mathrm{max}}$ , no runaway electrons are generated according to this fluid model. Note that in a kinetic model, energetic electrons born below the critical energy can end up in the runaway regime through collisional processes (Ekmark et al. Reference Ekmark, Hoppe, Fülöp, Jansson, Antonsson, Vallhagen and Pusztai2024). Using the same normalisation as in § 2, the dimensionless runaway electron generation rate due to tritium beta decay takes the form

(4.7) \begin{equation} \bar {\gamma }_{\mathrm{T}} = \ln {2}\sqrt {5 + Z_{\mathrm{eff}}} \ln {\Lambda _{\mathrm{c}}} \frac {ec n_{\mathrm{T}}}{j_0} \frac {\tau _{\mathrm{c}}}{\tau _{\mathrm{T}}}F_\beta (W_{\mathrm{c}}), \end{equation}

and the tritium seed density becomes

(4.8) \begin{equation} n_{\mathrm{seed,T}} = \frac {x_1^2 \ln {2}}{2} \left ( \frac {a}{a_{\mathrm{wall}}} \right )^2 \frac {ec n_{\mathrm{T}}}{j_0} \frac {\tau _{\mathrm{CQ}}}{\tau _{\mathrm{T}}} \int _{E_{\mathrm{min}}}^{E_1} \frac {F_\beta (W_{\mathrm{c}})}{E} \, {\mathrm{d}} E, \end{equation}

where $E_{\mathrm{min}}$ is the value of the normalised electric field $E$ that corresponds to $W_{\mathrm{c}} = W_{\mathrm{max}}$ . Assuming complete screening of the ion nuclei, $E_{\mathrm{min}}$ becomes

(4.9) \begin{equation} E_{\mathrm{min}} = \frac {\sqrt {5 + Z_{\mathrm{eff}}}}{w_{\mathrm{max}}^2 + 2 w_{\mathrm{max}}} + 1 \approx 34, \end{equation}

where $w_{\mathrm{max}} = W_{\mathrm{max}}/(m_{\mathrm{e}} c^2)$ is the maximum beta decay energy normalised to the electron rest energy, and the numerical evaluation assumes that $Z_{\mathrm{eff}} = 1$ .

In the limit of complete screening of the ion nuclei, $F_\beta (W_{\mathrm{c}})/E$ can be integrated analytically to become

(4.10) \begin{align} \int \frac {F_\beta (W_{\mathrm{c}})}{E} \, {\mathrm{d}} E &= \ln {E} - \sqrt {5+Z_{\mathrm{eff}}}w_{\mathrm{max}}^{-7/2} w_{\mathrm{c}}^{1/2} \left (21 w_{\mathrm{max}} + \frac {45}{2} - \frac {5}{2} w_{\mathrm{c}}\right ) \nonumber \\ &\quad + \sum _{i = 1}^3 c_n \left (\frac {1}{w_{\mathrm{max}}}\right )^{(2n + 1)/2} G_{2n + 1}\left (Z_{\mathrm{eff}};w_{\mathrm{c}}\right )\!, \end{align}

where $w_{\mathrm{c}} = W_{\mathrm{c}}/(m_{\mathrm{e}} c^2)$ is the critical energy normalised to the electron rest energy, and the coefficients $c_1$ through $c_3$ are given by 35/4, 21/2 and 15/4, respectively. The function $G_{n}(Z_{\mathrm{eff}};w_{\mathrm{c}})$ is defined as

(4.11) \begin{align} G_{n}\left (Z_{\mathrm{eff}};w_{\mathrm{c}}\right ) & = 2^{n/2}\arctan {\left ( \sqrt {\frac {w_{\mathrm{c}}}{2}}\right )} - (5 + Z_{\mathrm{eff}})^{n/8} T_n(u) \arctan {\left (\frac {2 (5 + Z_{\mathrm{eff}})^{1/8} u \sqrt {w_{\mathrm{c}}} }{(5 + Z_{\mathrm{eff}})^{1/4} - w_{\mathrm{c}}}\right )} \nonumber \\ & \quad - (5 + Z_{\mathrm{eff}})^{n/8} v U_{n - 1} (u) \textrm {arctanh}{\left (\frac {2 (5 + Z_{\mathrm{eff}})^{1/8} v \sqrt {w_{\mathrm{c}}}}{(5 + Z_{\mathrm{eff}})^{1/4} + w_{\mathrm{c}}}\right )}, \end{align}

where $T_n$ and $U_n$ are Chebyshev polynomials of the first and second kinds, respectively, and

(4.12) \begin{align} u & = \frac {1}{\sqrt {2}} \sqrt {1 + (5 + Z_{\mathrm{eff}})^{-1/4}}, \end{align}

(4.13) \begin{align} v & = \frac {1}{\sqrt {2}} \sqrt {1 - (5 + Z_{\mathrm{eff}})^{-1/4}}. \end{align}

4.2. Compton scattering

The activated wall radiates gamma photons up to the MeV range that can, through Compton scattering, launch both free and bound electrons into the runaway region of momentum space. In this work, the energy spectrum $\varGamma _\gamma (W_\gamma )$ is given as a closed-form expression by fitting the function

(4.14) \begin{align} \varGamma _\gamma (W_\gamma ) & = \varGamma _0 \exp {[-\exp {(-z)} - z + 1]}, \end{align}

(4.15) \begin{align} z & = \frac {\ln {W_\gamma \,[\mathrm{MeV}]} + C_1}{C_2} + C_3(W_\gamma \, [\mathrm{MeV}])^2, \end{align}

to radiation transport calculations. The gamma spectrum is machine- and shot-dependent, and fitting parameters $C_i$ for a few relevant machines and plasma compositions (Martín-Solís et al. Reference Martín-Solís, Loarte and Lehnen2017; Ekmark et al. Reference Ekmark, Hoppe, Tinguely, Sweeney, Fülöp and Pusztai2025) are reproduced in table 1. The normalisation constant $\varGamma _0$ ensures that the total photon flux becomes $\varGamma _{\mathrm{flux}}$ .

Table 1.

Fitting coefficients for gamma photon spectra for three reference scenarios (Martín-Solís et al. Reference Martín-Solís, Loarte and Lehnen2017; Ekmark et al. Reference Ekmark, Hoppe, Tinguely, Sweeney, Fülöp and Pusztai2025). The final column explicitly lists the value of $\bar {\sigma }_{\mathrm{eff}}(0)$ used in (4.22).

The rate at which electrons are scattered into the runaway region by gamma photons emitted with a certain energy $W_\gamma$ is proportional to the cross-section $\sigma (W_\gamma ,W_{\mathrm{c}})$ , which is obtained by integrating the Klein–Nishina differential cross-section (Klein & Nishina Reference Klein and Nishina1929)

(4.16) \begin{equation} \frac {\mathrm{d}\sigma }{\mathrm{d}\varOmega } = \frac {r_{\mathrm{e}}^2}{2} \frac {W_\gamma '^2}{W_\gamma ^2} \left ( \frac {W_\gamma }{W_\gamma '} + \frac {W_\gamma '}{W_\gamma } - \sin ^2{\theta }\right )\!, \end{equation}

where $r_{\mathrm{e}} = e^2/(4\pi \varepsilon _0 m_{\mathrm{e}} c^2)$ is the classical electron radius, and $\theta$ and $W_\gamma '$ are the deflection angle and the energy of the scattered photon, respectively. Here, $W_\gamma '$ can be related to the incident photon energy and the deflection angle through the kinematic relation

(4.17) \begin{equation} W_\gamma ' = \frac {W_\gamma }{1 + ({W_\gamma }/{m_{\mathrm{e}} c^2}) (1 - \cos {\theta })}. \end{equation}

By letting $W_\gamma -W_\gamma ' = W_{\mathrm{c}}$ , the kinematic relation above can be used to find the critical angle $\theta _{\mathrm{c}}$ , i.e. the minimum photon deflection angle required to scatter an electron into the runaway region for a given photon energy and critical energy. The critical angle becomes

(4.18) \begin{equation} \cos {\theta _{\mathrm{c}}} = 1 - \frac {w_{\mathrm{c}}}{w_\gamma } \frac {1}{w_\gamma - w_{\mathrm{c}}}, \end{equation}

yielding the integrated Compton cross-section

(4.19) \begin{align} \sigma (W_\gamma ,W_{\mathrm{c}}) &= 2\pi \int _{\theta _{\mathrm{c}}}^\pi \frac {\mathrm{d}\sigma }{\mathrm{d}\varOmega } \sin {\theta } \, {\mathrm{d}} \theta \nonumber\\& = \frac {3\sigma _{\mathrm{T}}}{8} \Bigg \{ \frac {w_\gamma ^2 -2w_\gamma -2}{w_\gamma ^3} \ln \left[{\frac {1 + 2w_\gamma }{1 + w_\gamma (1 - \cos {\theta _{\mathrm{c}}})}}\right] \nonumber \\ &\quad + \frac {1}{2w_\gamma } \left [ \frac {1}{\big [ 1 + w_\gamma (1 - \cos {\theta _{\mathrm{c}}}) \big ]^2} - \frac {1}{(1 + 2w_\gamma )^2} \right ] \nonumber \\ &\quad - \frac {1}{w_\gamma ^3}\left [ 1 - w_\gamma -\frac {1 + 2w_\gamma }{1 + w_\gamma (1-\cos {\theta _{\mathrm{c}}})} - w_\gamma \cos {\theta _{\mathrm{c}}}\right ]\Bigg \} , \end{align}

where $\sigma _{\mathrm{T}} = 8 \pi r_{\mathrm{e}}^2/3$ is the Thomson scattering cross-section and $w_\gamma = W_\gamma /(m_{\mathrm{e}} c^2)$ is the photon energy normalised to the electron rest energy (Martín-Solís et al. Reference Martín-Solís, Loarte and Lehnen2017).

Using the gamma spectrum and cross-section above, the runaway electron generation rate due to Compton scattering can be approximately evaluated as

(4.20) \begin{equation} \gamma _{\mathrm{C}}(W_{\mathrm{c}}) \approx n_{\mathrm{e}}^{\mathrm{tot}} \int _0^\infty \varGamma _\gamma (W_\gamma ) \sigma (W_\gamma ,W_{\mathrm{c}}) {\mathrm{d}} W_\gamma \equiv n_{\mathrm{e}}^{\mathrm{tot}} \varGamma _{\mathrm{flux}} \sigma _{\mathrm{eff}}(W_{\mathrm{c}}), \end{equation}

where $\sigma _{\mathrm{eff}}(W_{\mathrm{c}})$ is the Compton cross-section averaged over the energy spectrum. The normalised seed density becomes

(4.21) \begin{equation} n_{\mathrm{seed,C}} = \frac {x_1^2}{2} \left ( \frac {a}{a_{\mathrm{wall}}} \right )^2 \frac {e c n_{\mathrm{e}}^{\mathrm{tot}}}{j_0} \tau _{\mathrm{CQ}} \varGamma _{\mathrm{flux}} \sigma _{\mathrm{T}} \int _{\bar {E}_{\mathrm{c}}^{\mathrm{eff}}}^{E_1} \frac {\bar {\sigma }_{\mathrm{eff}}(W_{\mathrm{c}})}{E} \, {\mathrm{d}} E, \end{equation}

where $\bar {\sigma }_{\mathrm{eff}} \equiv \sigma _{\mathrm{eff}}/\sigma _{\mathrm{T}}$ is plotted against the critical energy in figure 2(a) for DT plasmas in both ITER and SPARC. An upper bound on the Compton seed can be established by assuming that all electrons struck by a gamma photon run away. This overestimates the Compton seed by approximately half an order of magnitude compared with expression (4.21), so a reasonable estimate of the integrated Compton seed becomes

(4.22) \begin{equation} n_{\mathrm{seed,C}} \approx \frac {x_1^2}{4} \left ( \frac {a}{a_{\mathrm{wall}}} \right )^2 \frac {e c n_{\mathrm{e}}^{\mathrm{tot}}}{j_0} \tau _{\mathrm{CQ}} \varGamma _{\mathrm{flux}} \sigma _{\mathrm{T}} \bar {\sigma }_{\mathrm{eff}}(0) \ln {\left (\frac {E_1}{\bar {E}_{\mathrm{c}}^{\mathrm{eff}}}\right )}, \end{equation}

where a factor of 1/2 has been included to compensate for (part of) the overestimation. For convenience, the value of $\bar {\sigma }_{\mathrm{eff}}(0)$ is given in table 1 for each Compton spectrum.

As the plasma disrupts and the fusion reactions cease, the energy spectrum is assumed to retain its shape but have a photon flux that is three to four orders of magnitude lower than the fluxes listed in table 1 (Martín-Solís et al. Reference Martín-Solís, Loarte and Lehnen2017; Vallhagen et al. Reference Vallhagen, Hanebring, Artola, Lehnen, Nardon, Fülöp, Hoppe, Newton and Pusztai2024; Ekmark et al. Reference Ekmark, Hoppe, Tinguely, Sweeney, Fülöp and Pusztai2025). The model presented in this work has no explicit temporal resolution and is therefore unable to include a time-dependent photon flux; thus, the photon flux is set to be a constant factor of 1000 smaller than the value listed in table 1. The criterion is only logarithmically sensitive to the seed densities, so using the exact value is not essential to the validity of the model. In addition, early seed electrons produced before the fusion reactions cease are prone to radial losses, increasing the relevance of the reduced photon flux.

5. Criterion for significant runaway electron generation

The runaway electron generation is assumed to be significant if $n \sim 1$ , or equivalently

(5.1) \begin{equation} \mathcal{Z} \equiv N_{\mathrm{ava}} + \ln {n_{\mathrm{seed}}} \gt 0, \end{equation}

where, again, $n$ is the normalised runaway electron density introduced in § 2. Note that $n$ is equivalent to the ratio between the runaway electron current and the pre-disruption current, and that $n$ is allowed to be larger than unity, as the saturation of the runaway electron current was neglected in § 2. This criterion can be used to assess whether a post-disruption plasma with a given set of parameters can be expected to produce a large number of runaway electrons. Here, $\mathcal{Z}$ can be evaluated using either the analytical expressions (3.13), (4.10) and (4.22), or the semi-analytical expressions (2.20), (4.8) and (4.21).Footnote ² This makes the criterion a useful tool for quickly determining favourable regions in parameter space. To verify the criterion, it is compared with zero-dimensional fluid simulations using the Dream code (Hoppe et al. Reference Hoppe, Embréus and Fülöp2021) for both ITER- and SPARC-relevant parameters. Note that only tritium and Compton seeds are considered, as similar criteria based on Dreicer and hot-tail generation have previously been studied elsewhere (see Helander et al. Reference Helander, Eriksson and Andersson2002; Fülöp et al. Reference Fülöp, Smith and Pokol2009).

Figure 3.

Dream simulations in zero dimensions (filled contours) compared with inequality (5.1) (green and grey contours for analytical and semi-analytical expressions, respectively) using (a) ITER-relevant and (b) SPARC-relevant parameters. The criterion $\mathcal{Z} = 0$ approximately delineates regions in parameter space where a significant fraction of the Ohmic current is converted into a runaway electron current. The grey dotted contour represents $\mathcal{Z} = 0$ , evaluated using only the tritium seed. In the SPARC case, the dotted contour overlaps with the solid contour. Note the nonlinearity in the lower part of the colour map.

The result of this comparison is illustrated in figure 3, where the filled contours represent the fraction of the Ohmic current converted into runaway electron current over a large range of injected deuterium and neon densities. The initial plasma current was set to $I_{\mathrm{p}0} = {15}\,\mathrm{MA}$ in ITER and $I_{\mathrm{p}0} = {8.7}\,\mathrm{MA}$ in SPARC, and the tritium density (as well as the initial deuterium density) was set to $n_{\mathrm{T}} = {5\times 10^{19}}\,\mathrm{m^{-3}}$ in ITER and $n_{\mathrm{T}} = {2\times 10^{20}}\,\mathrm{m^{-3}}$ in SPARC. The temperature was taken to be constant throughout the simulation and was determined by solving (3.14)–(3.15) for $T_{\mathrm{e}}$ for each combination of $n_{\mathrm{D}}$ and $n_{\mathrm{Ne}}$ , using the pre-disruption Ohmic current.

As shown in figure 3, the contour $\mathcal{Z} = 0$ approximately delineates the region in parameter space where the zero-dimensional fluid model in Dream predicts 1 % runaway electron current conversion in ITER and 10 %–30 % runaway electron current conversion in SPARC. A runaway electron current conversion of 1 % in ITER is consistent with the requirement to keep the runaway electron current below ${150}\,\mathrm{k A}$ (Lehnen Reference Lehnen2021), but the criterion is less conservative for SPARC-like parameters due to the exponential dependency on the initial plasma current. However, it is possible to compensate for this by explicitly making the criterion more conservative, i.e. defining significant runaway electron generation as e.g. $n\sim 0.1$ or $n\sim 0.01$ .

To better understand the contribution from each generation mechanism, inequality (5.1) can be decomposed into its constituents, namely the tritium decay seed, the Compton scattering seed and the avalanche gain factor $\exp {(N_{\mathrm{ava}})}$ , all of which are illustrated in figure 4 for both ITER and SPARC. For comparison, the ratio $N_{\mathrm{ava}}^{\mathrm{RP}}/N_{\mathrm{ava}}$ is plotted as well. The quantitative values differ between the two devices but are qualitatively similar for a given source. Very few runaway electrons are generated in the lower left quadrant, above temperatures of $\gtrsim\!{50}\,\mathrm{eV}$ . Note that additional physics mechanisms not included in the model presented here, such as those related to a vertical displacement event (including wall impurities entering the plasma), could lead to runaway production even in this parameter region. Furthermore, in addition to minimising the potential damage caused by runaway electrons, a disruption mitigation system must also manage thermal loads and electromagnetic forces – constraints that are not considered in the present work.

Moving up in either neon or deuterium concentration, the equilibrium temperature drops, and a large electric field is induced. At higher deuterium concentrations, the free electron density (and thus the critical electric field) increases, yielding a lower value of the normalised electric field $E$ . In the neon rich top left quadrant, the critical electric field remains relatively low, and $E$ thereby becomes large. This results in a low critical momentum that enables tritium beta electrons to run away. This also coincides with strong avalanching due to the combination of a strong electric field (relative to $E_{\mathrm{c}}$ ) and a large number of bound target electrons, which contribute only weakly to the collisional drag. The result is a substantial generation of runaway electrons due to the combination of a large tritium seed (see figures 4 a and 4 b) and strong avalanching (see figures 4 e and 4 f). Note that the low critical momentum and short thermal quench time can also give rise to a high hot-tail generation rate, which may be greater than the tritium seed generation. But again, hot-tail generation is not considered in this work.

In this region of parameter space, the Compton scattering seed (figures 4 c and 4 d) is a couple of orders of magnitude weaker than the tritium seed and does not contribute significantly. However, the Compton seed is relatively insensitive to the critical energy and is therefore present everywhere in the considered density space where the temperature is sufficiently low. In fact, going to high deuterium densities enhances the Compton seed due to the larger number of available target electrons.

Figure 4.

Tritium (a and b) and Compton (c and d) seeds evaluated using (4.8) and (4.21), respectively. The avalanche gain factor (e and f) is evaluated using the semi-analytical expression (2.20). Furthermore, the semi-analytical formulation is compared with the analytical expression (3.13) (g and h). The expressions are evaluated for ITER-like parameters in the left column, and SPARC-like parameters in the right column.

In ITER, the avalanche gain is sufficiently strong that the Compton seed alone can be amplified into a large runaway current. Therefore, the region in density space where significant runaway occurs is limited by the balance between the Compton seed and the avalanche gain. In SPARC, the avalanche gain is approximately five to ten orders of magnitude weaker than in ITER, and the Compton seed is not sufficiently large to give rise to a significant runaway electron current according to this model. However, this may change if a larger value is chosen for $\varGamma _{\mathrm{flux}}$ or if a different value for the current density is chosen. In this work, the current density is taken to be the average current density $I_{\mathrm{p0}}/(\pi a^2)$ , which is equivalent to a flat current density profile. The current density can also be chosen as the on-axis value of some peaked profile. Using the current density profile below in § 5.1 as an example, the on-axis current density is approximately 40 % higher than the average current density, yielding an avalanche gain that is approximately five orders of magnitude higher in ITER. This would be in line with the maximum avalanche gain reported by, for example, Wang et al. (Reference Wang, Nardon, Artola, Bandaru and Hoelzl2024).

The contribution from the tritium seed is illustrated in figure 3, where, in addition to the full nuclear seed, $\mathcal{Z}$ is also evaluated using only the tritium seed, indicating where large runaway currents can be expected if the Compton seed is neglected. In the ITER case in figure 3(a), this significantly reduces the extent of the region of significant runaway current, indicating that the dominant seed near the solid boundary is indeed the Compton seed. In the SPARC case in figure 3(b), however, the dotted contour overlaps with the solid contour, indicating that tritium beta decay is the dominant nuclear seed wherever large runaway currents can be expected.

In addition, it is noteworthy that for ITER-relevant parameters and deuterium densities $\gtrsim\!{7\times 10^{21}}\,\mathrm{m^{-3}}$ , the plasma temperature drops to such an extent that deuterium recombines, leading to a strong inductive electric field and a large runaway electron current. This is consistent with previous findings in e.g. Vallhagen et al. (Reference Vallhagen, Embréus, Pusztai, Hesslow and Fülöp2020) and McDevitt et al. (Reference McDevitt, Tang, Fontes, Sharma and Chung2023).

The analytical version and the slightly more detailed semi-analytical version of $\mathcal{Z}$ make similar predictions regarding where in density space large runaway currents can be expected, although the analytical formulation tends to be a little more conservative. This stems from $N_{\mathrm{ava}}^{\mathrm{RP}}$ being slightly larger than $N_{\mathrm{ava}}$ in a plasma consisting of mostly ionised hydrogen and small amounts of weakly ionised neon or argon, as shown in figures 4(g) and 4(h). The reason for the good agreement between the analytical and semi-analytical versions of the criterion is that where the criterion is valid, i.e. where $n\lesssim 1$ , the factor $\sqrt {4 \bar {\nu }_{\mathrm{s}}(p_\star )^2 + \bar {\nu }_{\mathrm{s}}(p_\star ) \bar {\nu }_{\mathrm{D}}(p_\star ) }$ is close to its completely screened limit, and the correction for partially ionised impurities is small. In a plasma dominated by weakly ionised high- $Z$ impurities on the other hand, effects of partial screening become significant and should be treated more carefully. However, injecting high concentrations of impurities is not attractive from a runaway electron mitigation point of view, as the large number of target electrons can give rise to very large runaway electron currents (Hesslow et al. Reference Hesslow, Embréus, Vallhagen and Fülöp2019). Hence, this approximation is not a significant limitation of the analytical model.

5.1. Comparison with one-dimensional fluid simulations

As shown above, the criterion (5.1) is a good surrogate model for a zero-dimensional fluid model. However, there are cases where the radial dynamics is important. To demonstrate this, one-dimensional Dream simulations with radially varying current density and temperature profiles that were allowed to evolve self-consistently were performed. The shapes of the initial current density and temperature profiles are given by

(5.2) \begin{align} \hat {j}(r) & = \left [1 - \left ( \frac {r}{a} \right )^2 \right ]^{0.41}, \end{align}

(5.3) \begin{align} T(r) & = T_0 \left [ 1 - 0.99\left (\frac {r}{a}\right )^2 \right ], \end{align}

where $T_0 = {20}\,\mathrm{k eV}$ is the on-axis temperature before the thermal quench. The current density profile is normalised such that the current density integrates to the total plasma current $I_{\mathrm{p0}}$ . At $t = 0$ , a uniform density of neutral neon and deuterium with a temperature of 1 eV is injected. To simulate heat transport losses across the stochastic flux surfaces during the thermal quench, a Rechester–Rosenbluth diffusion coefficient (Rechester & Rosenbluth Reference Rechester and Rosenbluth1978) is used. The diffusion coefficient is given by $D\sim v_{\mathrm{th}}R_0 (\delta B/B)^2$ , where $v_{\mathrm{th}}$ is the thermal speed, and $\delta B/B$ is a measure of the magnetic field fluctuations, which in this work was set to $3\times 10^{-3}$ (Vallhagen et al. Reference Vallhagen, Hanebring, Artola, Lehnen, Nardon, Fülöp, Hoppe, Newton and Pusztai2024; Ekmark et al. Reference Ekmark, Hoppe, Tinguely, Sweeney, Fülöp and Pusztai2025) to obtain a thermal quench time of the order of 1 ms (Sweeney et al. Reference Sweeney2020; Hu et al. Reference Hu, Nardon, Hoelzl, Wieschollek, Lehnen, Huijsmans, van Vugt and Kim2021).

Introducing radial profiles and dynamics complicates the comparison between Dream and the criterion (5.1), as the local runaway electron production rate is constrained by atomic physics, which, in turn, depends on the local current density. The current density is lower towards the plasma edge, leading to a lower equilibrium temperature and potentially to deuterium recombination, which can result in enhanced avalanche generation. This is exemplified by Vallhagen et al. (Reference Vallhagen, Embréus, Pusztai, Hesslow and Fülöp2020), where their ‘case 3’ corresponds to a scenario in which the temperature in the outer part of the plasma is initially substantially lower than that in the core, leading to an off-axis runaway current. This type of behaviour is not captured by the reduced zero-dimensional model presented in this work.

Figure 5.

Dream simulations including radial profiles (filled contours) compared with inequality (5.1) (green and grey contours for analytical and semi-analytical expressions, respectively) using (a) ITER-relevant and (b) SPARC-relevant parameters. In ITER, the criterion $\mathcal{Z} = 0$ approximately delineates regions in parameter space where a significant fraction of the Ohmic current is converted into a runaway electron current. In SPARC, the criterion does not capture the off-axis runaway electron generation in the high $n_{\mathrm{D}}$ regime. The orange dashed contour indicates the 50 eV isotherm to illustrate how the equilibrium temperature limits the runaway electron generation. Note the nonlinearity in the lower part of the colour map.

In light of this, the criterion can be adapted to the one-dimensional fluid simulations by using a lower value of $j_0$ in the power balance (3.14) to better represent the cold outer plasma. For example, the atomic physics can be evaluated using the average current density $I_{\mathrm{p0}}/(\pi a^2)$ , whereas the runaway electron dynamics is determined using the on-axis current density. Overall, this makes the criterion more conservative compared with evaluating both using the on-axis current density, since a lower current density leads to a lower equilibrium temperature, which in turn results in a higher electric field and thus a higher runaway electron generation.

In fact, in some cases, the runaway electron generation rate can be viewed as being directly limited by the temperature. This can be seen in figure 5, where the $T_{\mathrm{e}} = {50}\,\mathrm{eV}$ contour is included to illustrate this point. This contour is obtained by solving (3.14)–(3.15) for $n_{\mathrm{D}}$ and $n_{\mathrm{Ne}}$ with $T_{\mathrm{e}} = {50}\,\mathrm{eV}$ and $j_0 = I_{\mathrm{p0}}/(\pi a^2)$ . For deuterium densities below ${1.5\times 10^{21}}\,\mathrm{m^{-3}}$ in ITER and ${1\times 10^{21}}\,\mathrm{m^{-3}}$ in SPARC, the contour $\mathcal{Z} = 0$ coincides with the contour $T_{\mathrm{e}} = {50}\,\mathrm{eV}$ , indicating that for higher temperatures, the critical momentum is sufficiently large to prohibit any significant runaway electron generation. Increasing or decreasing the current density used for the evaluation of the equilibrium temperature therefore shifts the contour to higher and lower densities. Choosing an excessively low Ohmic current could therefore lead the criterion to predict significant runaway electron generation in a region where the plasma would otherwise be too hot to produce a runaway electron beam. The remainder of the contour is determined by the balance of the seed and the avalanche gain, and therefore other parameters such as the total electron density, photon flux, tritium density, and perhaps most importantly, the current density used to evaluate the avalanche gain.

For ITER-relevant parameters, a reasonable agreement between the criterion and one-dimensional Dream fluid simulations is found. The exception is for low neon and high deuterium content, where the off-axis runaway current is not captured (see figure 5(a)). For SPARC-relevant parameters, the contour $\mathcal{Z} = 0$ approximately delineates the region of significant on-axis runaway electron generation in the top left quadrant of figure 5(b). Again, the criterion does not capture the off-axis runaway current that arises at large deuterium concentrations ( $\gtrsim\!{1\times 10^{22}}\,\mathrm{m^{-3}}$ ). This is expected, considering that the criterion has no radial resolution.

Recent modelling efforts have highlighted the importance of runaway electron losses due to the scraping off of the plasma column during a vertical displacement event (Martín-Solís et al. Reference Martín-Solís, Mier, Lehnen and Loarte2022; Wang et al. Reference Wang, Nardon, Artola, Bandaru and Hoelzl2024; Vallhagen et al. Reference Vallhagen, Hanebring, Fülöp, Hoppe, Votta and Pusztai2025; Bandaru et al. Reference Bandaru, Hoelzl, Artola and Lehnen2025). This process has the potential to significantly reduce the final runaway current, especially if the current profile is strongly flattened during the thermal quench, making the plasma prone to producing runaways near the edge. The corresponding non-trivial radial dynamics would be difficult to account for within the current simplified model. However, knowing when large runaway currents are expected in the absence of such losses is nevertheless useful, not only for a quick conservative exploration of the parameter space, but also because scrape-off losses are more effective if the runaway current is not too large even without them; otherwise the current decay stops early, strongly braking the scrape-off process (Vallhagen et al. Reference Vallhagen, Hanebring, Fülöp, Hoppe, Votta and Pusztai2025).