2.1. Ions and neutrals

The tool STREAM allows for an arbitrary number of ion species to be considered in the simulation and evolves the density of each charge state separately. The fact that the ionization mean-free path of neutrals decreases as the plasma develops, so that neutrals can be screened out of the plasma core, is accounted for by using the two-volume model introduced by Lloyd, Carolan & Warrick (Reference Lloyd, Carolan and Warrick1996). In the two-volume model, illustrated in figure 1, the plasma is assumed to occupy a volume $V_{p}$, with neutrals of species $i$ only being able to penetrate a subvolume $V_{n,i}$ of the plasma. In addition, neutrals are assumed to be homogeneously spread out in the volume consisting of $V_{n,i}$ and region outside of the plasma. This means that the total volume occupied by the neutrals of species $i$ is $\gamma _{n,i} V$, where the neutral volume coefficient is given by

(2.1)\begin{equation} \gamma_{n,i} = 1 -\frac{V_{p} - V_{n,i}}{V}. \end{equation}

It is only the neutrals within the plasma which contribute to atomic processes.

Figure 1. The two-volume model assumes that neutral particles move freely inside the vacuum vessel of volume $V$, except in the centre of the plasma volume $V_{p}$ where plasma formation has come sufficiently far for neutrals to be fully screened out. The region of the plasma where neutrals and ions of species $i$ coexist has volume $V_{n,i}$.

The geometry of the plasma is specified using an analytic magnetic field model described in appendix A of Hoppe et al. (Reference Hoppe, Embreus and Fülöp2021). The model includes the effects of elongation and triangularity on the plasma. From the shaping profiles prescribed, the plasma volume is computed numerically by

(2.2)\begin{equation} V_{p} = \int_0^a \,\mathrm{d} r \int_0^{2 {\rm \pi}} \,\mathrm{d}\phi \int_{-{\rm \pi}}^{\rm \pi} \mathcal{J}\,\mathrm{d}\theta, \end{equation}

where $a$ is the (possibly time evolving) plasma minor radius and

(2.3)\begin{equation} \mathcal{J} = \frac{1}{|\boldsymbol{\nabla} \phi \boldsymbol{\cdot} (\boldsymbol{\nabla} \theta \times \boldsymbol{\nabla} r )|} \end{equation}

is the spatial Jacobian, $r$ is the minor radius coordinate and $\phi$ and $\theta$ are the toroidal and poloidal angles, respectively.

The volume occupied by the neutrals of species $i$ is determined from a formula accounting for the plasma elongation $\kappa$ and triangularity $\delta$ by subtracting the ionization mean-free path for species $i$ from the plasma radius, so that

(2.4)\begin{equation} V_{n,i} = 2{\rm \pi}^2 R_0\kappa\left[ a^2 - \left(a-\lambda_i\right)^2 \right] + 2\kappa\delta\left(8-3{\rm \pi}^2\right)\left[ a^3 - \frac{\left(a-\lambda_i\right)^3}{3} \right]. \end{equation}

The ionization mean-free path for neutrals of species $i$ is determined by (Lloyd et al. Reference Lloyd, Carolan and Warrick1996)

(2.5)\begin{equation} \lambda_i = \frac{v_{{\rm th},i}}{n_{e} I_i^{(0)}}, \end{equation}

with $v_{{\rm th},i} = \sqrt {2T_i/m_i}$ the thermal ion speed; $T_i$ and $m_i$ the temperature and mass, respectively, for species $i$; $n_{e}$ the electron density; and $I_i^{(0)}$ the rate at which neutrals of species $i$ are ionized. The ionization rate $I_i^{(0)}$ is taken from the Atomic Data and Analysis Structure (ADAS) (Summers Reference Summers2004).

2.1.1. Particle balance

The density of ions of species $i$ in charge state $j$ evolves through ionization, recombination and charge exchange processes according to

(2.6)\begin{gather} \frac{\mathrm{d} n_i^{(j)}}{\mathrm{d} t} = \frac{1}{V_i^{(j)}}\left[ \hat{V}_i^{(j-1)} I_i^{(j-1)} n_{e} n_i^{(j-1)} - \hat{V}_i^{(j)}I_i^{(j)} n_{e} n_i^{(j)}\right.\nonumber\\ \left.+\hat{V}_i^{(j+1)} R_i^{(j+1)} n_{e} n_i^{(j+1)} - \hat{V}_i^{(j)} R_i^{(j)} n_{e} n_i^{(j)} + \hat{V}_\star^{(0)} n_\star^{(0)} A_{i,{\rm cx}}^{(j)} \right] + S_i^{(j)}. \end{gather}

Here, $I_i^{(j)}$ is the rate of ionization of species $i$ from charge state $j$ to $j+1$ and $R_i^{(j)}$ the rate of recombination of species $i$ from charge state $j$ to $j-1$,Footnote ² both of which are taken from ADAS. The total volume $V_i^{(j)}$ occupied by ions of species $i$ in charge state $j$ is given by

(2.7)\begin{equation} V_i^{(j)} = \begin{cases} \gamma_{n,i} V, & j=0,\\ V_{p}, & j \geq 1, \end{cases} \end{equation}

while the volume inside the plasma $\hat {V}_i^{(j)}$ occupied by ions of species $i$ and charge state $j$ is

(2.8)\begin{equation} \hat{V}_i^{(j)} = \begin{cases} V_{n,i}, & j=0,\\ V_{p}, & j \geq 1. \end{cases} \end{equation}

The interaction term $A_{i,{\rm cx}}^{(j)}$, for charge exchange with neutrals of the main ion species (denoted with a star) and an impurity ion, is

(2.9)\begin{equation} A_{i,{\rm cx}}^{(j)} = \begin{cases} \left({-}1\right)^{(j+1)} \sum_{k,l\geq 1} R_{ik,{\rm cx}}^{(l)} n_k^{(l)}, & \quad i={\star}, k\neq{\star},\\ R_{ik,{\rm cx}}^{(j+1)} n_i^{(j+1)} - R_{ik,{\rm cx}}^{(j)}n_i^{(j)}, & \quad i\neq{\star}, k={\star}, \end{cases} \end{equation}

with $R_{ik,{\rm cx}}^{(j)}$ denoting the charge exchange rate, which is also taken from ADAS.

The external particle source/sink $S_i^{(j)}$ models the influx of neutrals from the wall, as well as the outflux of ions due to transport,

(2.10)\begin{equation} S_i^{(j)} = \begin{cases} \varGamma_{i,{\rm in}}^{(0)}/V_i^{(0)}, & j=0,\\ -n_i^{(j)}/\tau_i, & j\geq 1, \end{cases} \end{equation}

where $\tau _i$ is the confinement time for species $i$. The neutral influx $\varGamma _{i,{\rm in}}^{(0)}$ is generally given by

(2.11)\begin{equation} \varGamma_{i,{\rm in}}^{(0)} = V_{p}\sum_k\sum_{l\geq 1} \frac{Y_k^i n_k^{(l)}}{\tau_k}. \end{equation}

The sputter yield $Y_k^i$ for species $i$ due to the bombardment of incident species $k$ is prescribed by the user. The value of $Y_k^i$ depends on the underlying physical mechanism. The sputtering can be due to chemical sputtering (e.g. in the case of carbon wall (Kim et al. Reference Kim, Fundamenski and Sips2012)), in which case $Y_k^i$ is approximately constant with a value based on plasma surface experiment data, or physical sputtering, when $Y_k^i$ depends on the incident ion energy (i.e. ITER wall) (Kim, Sips & Contributors Reference Kim, Sips and Contributors2013).

The number of free electrons in the plasma is constrained via quasi-neutrality to

(2.12)\begin{equation} n_{e} = \sum_{ij} Z_{0,i}^{(j)} n_i^{(j)}, \end{equation}

where $Z_{0,i}^{(j)}$ denotes the net charge number of ion species $i$ in charge state $j$, i.e. $Z_{0,i}^{(j)}=j$.

2.1.2. Particle confinement

The confinement time varies greatly during start-up and is determined by different mechanisms during different phases of the start-up. In the early stages, before closed flux surfaces have formed, the confinement is set by transport along magnetic field lines to the wall. With an effective connection length $L_{f}$, the thermal particle confinement time can be estimated as (Kim et al. Reference Kim, Fundamenski and Sips2012)

(2.13)\begin{equation} \tau_{i,\parallel} = \frac{L_{f}}{C_{s}}, \end{equation}

where $C_{s} = \sqrt {(T_{e}+T_{\rm \star })/m_\star }$ is the main ion sound speed. It was shown by Kim et al. (Reference Kim, Fundamenski and Sips2012) that the same confinement time can also be used for impurities when they are trace.

At low plasma currents, the effective connection length will depend on the magnitude of stray magnetic fields. As the plasma current increases and gradually exceeds the eddy currents in surrounding conducting structures, the connection length increases and approaches infinity as closed flux surfaces form. We therefore model the effective connection length as in DYON (Mineev et al. Reference Mineev, Belyakov, Gribov, Kavin, Lobanov, Belov, Lamzin and Sytchevsky2014; Kim et al. Reference Kim, Mineev, Ricci, Lee and Na2020),

(2.14)\begin{equation} L_{f} = \frac{3a}{4}\frac{B_\phi}{B_z}\exp\left(\frac{I_{p}}{I_{\rm ref}}\right), \end{equation}

where $B_\phi$ is the toroidal magnetic field strength, $I_{p}$ the total plasma current and $I_{\rm ref}$ denotes the plasma current where closed flux surfaces form and the effect of stray magnetic fields greatly reduces. The stray field $B_z$ is composed of the vertical magnetic field $B_{v}$ and a component $B_{\rm eddy}$ resulting from eddy currents in conducting structures surrounding the plasma,

(2.15)\begin{equation} B_{\rm eddy} = \frac{\mu_0}{2{\rm \pi} l_{\rm wall}}I_{\rm wall}, \end{equation}

where $\mu _0$ is the permeability of free space and $l_{\rm wall}$ denotes the distance between the centre of the plasma and the conducting structure.

During later stages of the discharge, the thermal particle confinement is expected to be dominated by turbulent transport. We can therefore estimate the confinement time during the later stages using the Bohm scaling

(2.16)\begin{equation} \tau_{i,\perp} = \frac{a^2}{2D_{\rm Bohm}}, \end{equation}

where the Bohm diffusion coefficient $D_{\rm Bohm} = T_{e}/(16eB_\phi )$. To allow the confinement time to smoothly transition from (2.13) to (2.16) during the start-up, we let the total confinement time $\tau _i$ satisfy

(2.17)\begin{equation} \frac{1}{\tau_i} = \frac{1}{\tau_{i,\parallel}} + \frac{1}{\tau_{i,\perp}}. \end{equation}

2.2. Heat

The temperature is evolved separately for each plasma species. The bulk electron temperature $T_{e}$ is related to the electron thermal energy density $W_{e} = 3n_{e} T_{e}/2$, and the thermal energy density is in turn evolved according to

(2.18)\begin{equation} \frac{\mathrm{d} W_{e}}{\mathrm{d} t} = j_\varOmega E_\parallel{-} n_{e} \sum_i\sum_{j=0}^{Z_i} \frac{\hat{V}_i^{(j)}}{V_i^{(j)}} n_i^{(j)} L_i^{(j)} - \sum_i Q_{{\rm e}i} - \frac{W_{e}}{\tau_{e}}. \end{equation}

The first term represents the ohmic heating provided by the electric field, and the second term represents energy losses via inelastic atomic processes with

(2.19)\begin{equation} L_i^{(j)} = L_{\rm line} + L_{\rm free} + \Delta W_i^{(j)}\left( I_i^{(j)} - R_i^{(j)} \right), \end{equation}

where $L_{\rm line}$ is the radiated power by line radiation, $L_{\rm free}$ by recombination radiation and bremsstrahlung, and the other terms represent the change in potential energy due to excitation and recombination, with the same rate coefficients as in (2.6). Note that $R_i^{(j)}$ and $L_{\rm free}$ vanish for $j=0$, while $I_i^{(j)}$ vanishes for $j=Z_i$. The ionization threshold $\Delta W_i^{(j)}$ is retrieved from the NIST database (Kramida, Ralchenko & Reader Reference Kramida, Ralchenko and Reader2020) while the other rate coefficients are taken from ADAS. The collisional heat transfer $Q_{kl}$, which appears in the term for heat exchange with ions in (2.18), is generally given for two arbitrary plasma species $k$ and $l$ by

(2.20)\begin{equation} Q_{kl} = \frac{\left\langle nZ^2 \right\rangle_k \left\langle nZ^2\right\rangle_l e^4\ln\varLambda_{kl}} {\left(2{\rm \pi}\right)^{3/2}\epsilon_0^2 m_km_l} \frac{T_k - T_l}{\left(\dfrac{T_k}{m_k}+\dfrac{T_l}{m_l}\right)^{3/2}}, \end{equation}

with the weighted charge $\langle nZ^2\rangle _i=\sum _{j=0}^{Z_i}n_i^{(j)}(Z_{0,i}^{(j)})^2$. Finally, the electron confinement time $\tau _{e}$ is taken to be the same as that of the main ions, given in (2.17), due to the requirement of ambipolarity. The Coulomb logarithms are modelled using (Wesson Reference Wesson1997)

(2.21a)\begin{gather} \ln\varLambda_{\rm ee} = \ln\varLambda_{{\rm e}i} = 14.9 + \ln\left(\frac{T_{e}}{1\,\text{keV}}\right) - 0.5\ln\left(\frac{n_{e}}{10^{20}\,\text{m}^{{-}3}}\right), \end{gather}

(2.21b)\begin{gather}\ln\varLambda_{ii} = 17.3 + \frac{3}{2}\ln\left(\frac{T_{e}}{1\,\text{keV}}\right) - 0.5\ln\left(\frac{n_{e}}{10^{20}\,\text{m}^{{-}3}}\right), \end{gather}

for electron–ion and ion–ion collisions, respectively. The thermal energy density for each ion species evolves according to

(2.22)\begin{equation} \frac{\mathrm{d} W_i}{\mathrm{d} t} = \sum_k Q_{ik} - \frac{3}{2}\frac{\hat{V}_{{\star}}^{(0)}}{V_{p}} n_\star^{(0)}\left( T_i - T_0 \right) R_{i,{\rm cx}}^{(1)}n_i^{(1)} - \frac{W_i}{\tau_i}, \end{equation}

where $T_0 = 300\,{\rm K}$ is the temperature of the main species neutrals, $\tau _i$ is given by (2.17), and the sum in the first term runs over all particle species, including electrons.

2.3. Electric field

The toroidal electric field inside the plasma depends on the externally applied loop voltage as well as the time rate of change of the plasma current. If the plasma is surrounded by conducting structures, such as a metallic wall, the plasma will be inductively coupled to these structures and the electric field and current dynamics will be correspondingly affected. In STREAM, we use the same model for the electric field and current dynamics as in DYON (Kim et al. Reference Kim, Fundamenski and Sips2012), described by one circuit equation for the plasma and one for the conducting structure surrounding it,

(2.23a)\begin{gather} 2{\rm \pi} R_0 E_\parallel + L_{p}\frac{\mathrm{d}I_{p}}{\mathrm{d} t} + M\frac{\mathrm{d}I_{\rm wall}}{\mathrm{d} t} = V_{{\rm loop},{\rm ext}}, \end{gather}

(2.23b)\begin{gather}V_{{\rm loop},{\rm wall}} + L_{\rm wall}\frac{\mathrm{d}I_{\rm wall}}{\mathrm{d} t} + M\frac{\mathrm{d}I_{p}}{\mathrm{d} t} = V_{{\rm loop},{\rm ext}}, \end{gather}

where $V_{{\rm loop},{\rm ext}}$ is the externally applied loop voltage, $R_0$ is the tokamak major radius, $E_\parallel$ the parallel electric field inside the plasma, $I_{p}$ the total plasma current, $I_{\rm wall}$ the current in the conducting structure and $V_{{\rm loop},{\rm wall}}=R_{\rm wall}I _{\rm wall}$ is the loop voltage in the conducting structure. The inductances $L_{p}$ (plasma inductance), $L_{\rm wall}$ (wall inductance) and $M$ (plasma–wall mutual inductance), and the wall resistance $R_{\rm wall}$, are free parameters in the model.

Due to STREAM's heritage from DREAM, the total plasma current $I_{p}$ is obtained from the total current density $j_{\rm tot}$ in the system by multiplying with the plasma cross-sectional area. An important difference between STREAM and DYON is that in STREAM the current density in turn consists of an ohmic component $j_{\varOmega } = \sigma E_\parallel$, where $\sigma$ is the conductivity of the plasma as determined by Redl et al. (Reference Redl, Angioni, Belli and Sauter2021), and a runaway electron component $j_{\rm re} = ecn_{\rm re}$, with $e$ the elementary charge, $c$ the speed of light and $n_{\rm re}$ the runaway electron density. The calculation of the runaway electron density $n_{\rm re}$ is described in detail in § 2.4.

2.4. Runaway electrons

Electrons are said to ‘run away’ when the collisional friction acting on them is weaker than other accelerating forces. In most cases, the accelerating force is an electric field, and in order for it to provide net acceleration to electrons it must exceed the Connor–Hastie threshold (Connor & Hastie Reference Connor and Hastie1975)

(2.24)\begin{equation} E > E_{c} = \frac{e^3n_{e}\ln\varLambda_{\rm ee}}{4{\rm \pi}\epsilon_0^2m_{e} c^2}, \end{equation}

where $\epsilon _0$ is the vacuum permittivity and $m_{e}$ the electron rest mass. Above this threshold, all electrons with relativistic momentum $p > p_{c}\approx 1/\sqrt {E/E_{c}-1}$ will be freely accelerated and run away. Electrons can find themselves above $p=p_{c}$, in the so-called runaway region, in a number of ways. For example, they can enter this region in momentum space through the collisionally diffusive leak from the thermal population at a steady rate – the so-called Dreicer generation (Dreicer Reference Dreicer1959), that is exponentially sensitive to the electric field normalized to the Dreicer field,

(2.25)\begin{equation} E_{D}= E_{c}\frac{m_{e} c^2}{T_{e}} = \frac{e^3n_{e}\ln\varLambda_{\rm ee}}{4 {\rm \pi}\epsilon_0^2 T_{e}}. \end{equation}

When $E\approx 0.215E_{D}$, the electric field acceleration exceeds the collisional friction for all electrons, leading to slide-away and a distribution that is far from thermal. Slide-away can also occur at lower values of $E$ since the Dreicer generation mechanism will gradually drain the thermal bulk of particles and pull the electrons into the runaway region.

Existing runaway electrons can also create new ones through close collisions with thermal electrons (Jayakumar, Fleischmann & Zweben Reference Jayakumar, Fleischmann and Zweben1993). This leads to an exponential growth of the number of runaway electrons – an avalanche. The avalanche mechanism has caused much concern for the potentially highly aggressive growth of the number of runaway electrons during tokamak disruptions (Boozer Reference Boozer2015; Breizman et al. Reference Breizman, Aleynikov, Hollmann and Lehnen2019), which is seen as almost inevitable in future reactors due to the presence of additional primary sources of runaway electrons. The first of these, labelled the ‘hot-tail’ mechanism, occurs during the thermal quench of a disruption as the fastest electrons of the original hot thermal distribution take longer to slow down. They may therefore be accelerated by the electric field which is induced as the plasma current drops due to the increased resistivity, before they have time to thermalize.

During tritium operation, energetic electrons resulting from tritium decay and Compton scattering of $\gamma$ photons originating from the activated wall are predicted to also provide a significant seed electron population to drive the avalanche multiplication (Martın-Solıs, Loarte & Lehnen Reference Martın-Solıs, Loarte and Lehnen2017; Fülöp et al. Reference Fülöp, Helander, Vallhagen, Embreus, Hesslow, Svensson, Creely, Howard and Rodriguez-Fernandez2020; Vallhagen et al. Reference Vallhagen, Embreus, Pusztai, Hesslow and Fülöp2020).

During tokamak start-up, the Dreicer and avalanche mechanisms are expected to dominate the generation of runaway electrons, while the other generation mechanisms do not contribute significantly. Hot-tail generation is only expected to occur in cooling plasmas, and is not relevant to the conditions of the initiating plasma. Tritium decay and Compton scattering sources of runaway electrons are typically much lower than Dreicer generation due to the large normalized electric field $E/E_{D}$ and low electron densities that are typical during plasma initiation.

2.4.1. Runaway electron model

Runaway electrons primarily influence the plasma evolution during start-up by contributing a relativistic current component to the plasma current. Since the distribution of runaways is usually strongly beamed along the magnetic field lines, with speeds close to the speed of light, the runaways contribute a current density $j_{\rm re}=ecn_{\rm re}$. This current density will in turn affect the evolution of the loop voltage in (2.23) and, by extension, the ohmic current and heating.

The runaway electron density is evolved according to

(2.26)\begin{equation} \frac{\partial n_{\rm re}}{\partial t} = \gamma_{\rm Dreicer} + \varGamma_{\rm ava}n_{\rm re} - \frac{n_{\rm re}}{\tau_{\rm re}}, \end{equation}

where $\gamma _{\rm Dreicer}$ is the rate at which electrons are generated by the Dreicer mechanism (Dreicer Reference Dreicer1959), $\varGamma _{\rm ava}$ is the avalanche growth rate by which runaway electrons exponentially multiply (Rosenbluth & Putvinski Reference Rosenbluth and Putvinski1997) and $\tau _{\rm re}$ is the runaway electron confinement time. The Dreicer generation rate is evaluated using the neural network developed by Hesslow et al. (Reference Hesslow, Unnerfelt, Vallhagen, Embreus, Hoppe, Papp and Fülöp2019b), while the avalanche growth rate is given by the semianalytical formula (Hesslow et al. Reference Hesslow, Embréus, Vallhagen and Fülöp2019a)

(2.27)\begin{equation} \left. \begin{gathered} \varGamma_{\rm ava} = \frac{e}{m_{e} c\ln\varLambda_{c}} \frac{n_{e}^{\rm tot}}{n_{e}} \frac{E{-}E_{c}^{\rm eff}} {\sqrt{4+\bar{\nu}_{s}\left(p_\star\right)\bar{\nu}_{D}\left(p_\star\right)}},\\ p_\star{=} \frac{\sqrt[4]{\bar{\nu}_{s}\left(p_\star\right)\bar{\nu}_{D}\left(p_\star\right)}} {\sqrt{E{/}E_{c}}},\\ \ln\varLambda_{c} = 14.6 + 0.5\ln\left( \frac{T_{e}\,[{\rm keV}]} {n_{e}\,[10^{20}\,{\rm m}^{{-}3}]} \right), \end{gathered} \right\} \end{equation}

where $n_{e}^\textrm {tot}$ is the total density of electrons (bound and free), $\ln \varLambda _{c}$ is a generalized Coulomb logarithm evaluated at relativistic energies, and $\bar {\nu }_{s}$ and $\bar {\nu }_{D}$ are the normalized slowing-down and deflection frequencies as defined by Hesslow et al. (Reference Hesslow, Embréus, Hoppe, DuBois, Papp and Fülöp2018b). The effective critical electric field $E_{c}^\textrm {eff}$ is calculated as described in appendix C.2 of Hoppe et al. (Reference Hoppe, Embreus and Fülöp2021), and takes into account the effect of bremsstrahlung and synchrotron radiation based on the approach of Hesslow et al. (Reference Hesslow, Embréus, Wilkie, Papp and Fülöp2018a). In partially ionized plasmas, both the Dreicer generation rate and the avalanche growth rate depend on the extent to which fast electrons can penetrate the bound electron cloud around the impurity ion. This effect of partial screening is taken into account in both the neural network for Dreicer generation and in (2.27) for the avalanche growth rate.

Using state-of-the-art models for the effect of partial screening is the main difference between the runaway generation models used in the SCENPLINT code and STREAM. As shown by Hesslow et al. (Reference Hesslow, Unnerfelt, Vallhagen, Embreus, Hoppe, Papp and Fülöp2019b) these effects can lead to orders of magnitude differences in the Dreicer generation rate, in plasmas containing partially ionized atoms. Also, the generalized expression for the avalanche growth rate has been shown to lead to large differences in the final runaway current, compared with previously used expressions (e.g. in Martın-Solıs et al. (Reference Martın-Solıs, Loarte and Lehnen2017)), as demonstrated in Hesslow et al. (Reference Hesslow, Embréus, Vallhagen and Fülöp2019a) and Vallhagen et al. (Reference Vallhagen, Embreus, Pusztai, Hesslow and Fülöp2020).

As for ions, the runaway electron confinement time varies significantly during start-up (Kavin, Mineev & Lobanov Reference Kavin, Mineev and Lobanov2017). In the early stages, before flux surfaces have formed, the transport will be dominated by parallel transport, and the distance traversed by a runaway before leaving the plasma is the same as for the thermal particles, i.e. the connection length given in (2.14). Starting from the equation of motion for the electron,

(2.28)\begin{equation} \frac{\mathrm{d} p}{\mathrm{d} t} = eE, \end{equation}

where $p$ is the electron momentum, and letting $E$ be positive in the cocurrent direction, we can estimate its confinement time. By assuming $E$ to be constant in time, integrating equation (2.28) and solving for the electron speed $v$, we obtain

(2.29)\begin{equation} v \equiv \frac{\mathrm{d} s}{\mathrm{d} t} = \frac{eE t/m_{e}}{\sqrt{1+\left(\dfrac{eE t}{m_{e} c}\right)^2}}. \end{equation}

This can in turn be integrated over the time period $\tau _\textrm {re}^{\parallel }$ it takes for the electron to travel a distance $L_{f}$, to obtain the relation

(2.30)\begin{equation} L_{f} = \frac{m_{e} c^2}{eE}\left[ \sqrt{1+\left(\frac{eE\tau_{\rm re}^{{\parallel}}}{m_{e} c}\right)^2} - 1 \right], \end{equation}

which can be solved for the runaway electron confinement time

(2.31)\begin{equation} \tau_{\rm re}^{{\parallel}} = \frac{m_{e} c}{eE}\sqrt{\left(\frac{eE L_{f}}{m_{e} c^2}+1\right)^2-1}. \end{equation}

After flux surfaces have formed, the confinement time is rather set by the rate at which the runaways gain energy and drift out of the plasma. Assuming that a runaway electron drifts out of the plasma when reaching a relativistic energy $m_{e} c^2\gamma _\textrm {max}$, and assuming the electric field to be roughly constant during acceleration, the confinement time in this second phase can be related to the runaway electron energy via

(2.32)\begin{equation} p_{\rm max} = \int_0^{\tau_{\rm re}^{\rm drifts}} eE\,\mathrm{d} t = eE\tau_{\rm re}^{\rm drifts}. \end{equation}

The maximum energy $\gamma _\textrm {max}=\sqrt {1+p_\textrm {max}^2/m_{e}^2c^2}$ for an electron was estimated by Knoepfel & Spong (Reference Knoepfel and Spong1979) as

(2.33)\begin{equation} \gamma_{\rm max}\approx\frac{56R_0}{a}I_{p}\,[\mathrm{MA}], \end{equation}

where $a$ and $R_0$ are the plasma minor and major radii, respectively. The runaway electron confinement time due to drifts can therefore be taken as

(2.34)\begin{equation} \tau_{\rm re}^{\rm drifts}\,[\mathrm{s}]\approx\frac{R_0}{10a} \frac{I_{\rm p}\,[\mathrm{MA}]}{E\,[\mathrm{V}\,\mathrm{m}^{{-}1}]}. \end{equation}

Since the transition from the first to the second stage is determined by when closed flux surfaces form, which depends on the development of the plasma current, we interpolate between the confinement times (2.31) and (2.34) using

(2.35)\begin{equation} \frac{1}{\tau_{\rm re}} = \frac{\exp\left({-}I_{p}/I_{\rm ref}\right)}{\tau_{\rm re}^{{\parallel}}} + \frac{1-\exp\left({-}I_{p}/I_{\rm ref}\right)}{\tau_{\rm re}^{\rm drifts}}, \end{equation}

with $I_\textrm {ref}$ the current at which flux surfaces form. Here $I_\textrm {ref}$ is assumed to be approximately $100\,\textrm {kA}$ in ITER (Kim et al. Reference Kim, Mineev, Ricci, Lee and Na2020).

3.1. ITER scenario

We start by considering an idealized ITER ohmic burn-through scenario with constant input parameters, similar to the one considered in § 2 of Kim et al. (Reference Kim, Mineev, Ricci, Lee and Na2020). We assume a circular plasma and a constant loop voltage corresponding to $0.3\,\textrm {V}\,\textrm {m}^{-1}$. The prefilled fuel gas pressure is $0.8\,\textrm {mPa}$, corresponding to an initial hydrogen density of $3.84\times 17\,\textrm {m}^{-3}$. The gas consists of pure hydrogen (no impurities) at $T_i=0.026\,\textrm {eV}$ (room temperature) in the beginning of the simulation. The plasma inductance is calculated from $L=\mu _0 R_{0} [\ln {(8R_{0}/a)}+l_{i}/2-2]$, where $l_{i}=0.5$ is assumed for the internal inductance and we neglect all currents in surrounding passive structures. The wall recycling coefficient for hydrogen $Y_{H}^{H}$ is set equal to one (no external gas fuelling). The input parameters are summarized in table 2. Note that while hydrogen is used as the main ion species, STREAM has been temporarily modified in these simulations to use charge-exchange rates for deuterium, which should correspond more closely to the charge-exchange rates used for this scenario in Kim et al. (Reference Kim, Mineev, Ricci, Lee and Na2020).

Table 2. Simulation parameters for the ITER and JET benchmark cases. Some parameters appear only in the models used in the JET case.

Figure 2 shows the time evolution of the plasma current, effective connection length, particle confinement time, electron density, electron temperature and ion temperature. After modifying STREAM to use the ADAS charge-exchange rates for deuterium also for hydrogen, as was done in (Kim et al. Reference Kim, Mineev, Ricci, Lee and Na2020), good agreement was found between the simulation results of DYON and STREAM. Only a slightly delayed burn-through is seen in STREAM which leads to a delayed rise of the plasma current, and a slightly faster heating rate for both ions and electrons in STREAM. These small differences observed could potentially be explained by the different conductivities used (the STREAM simulations use the conductivity of Redl et al. (Reference Redl, Angioni, Belli and Sauter2021) while the DYON simulations use the Spitzer conductivity), the use of different line radiation and charge exchange rates (the STREAM simulations use ADAS rates while the DYON simulations use analytical fits), and differences in the electron–ion equilibration term (2.20) (the DYON simulations use a constant value $\ln \varLambda =10$ for the Coulomb logarithm).

Figure 2. Time evolution of plasma parameters for the idealized ITER scenario: (a) plasma current; (b) effective connection length;(c) confinement time; (d) electron density; (e) electron temperature; and (f) ion temperature. The input parameters are given in table 2. The solid blue line is obtained by STREAM, dashed black line by DYON.

3.2. JET scenario

Next, we consider an advanced burn-through scenario of a deuterium plasma with multiple impurity species and a conducting structure near the plasma. As in Kim et al. (Reference Kim, Mineev, Ricci, Lee and Na2020), we use parameters from a JET discharge with a carbon wall, which was used to validate the DYON code. In this case, both the breakdown region and the loop voltage are time dependent, and the data is given in figure 10 of Kim et al. (Reference Kim, Mineev, Ricci, Lee and Na2020). The plasma minor radius is inferred from the plasma volume given in figure 10(a). The eddy current is calculated with the two-rings circuit model in (2.23). To obtain agreement it was necessary to modify the effective connection length according to $L_{f}\to L_{f}/3$ from (2.14) in STREAM to use the exact same model as was used in Kim et al. (Reference Kim, Mineev, Ricci, Lee and Na2020).

In the simulation, a time evolving deuterium recycling coefficient is used. The coefficient is given generally by

(3.1)\begin{equation} Y_{D}^{D}(t) = c_1 - c_2\left(1 - \mathrm{e}^{{-}t/c_3}\right), \end{equation}

and in this scenario we set $c_1=1.1$, $c_2=0.05$ and $c_3=0.1\,\textrm {s}$.

Figure 3 shows the time evolution of the plasma current, effective connection length, particle confinement time, electron density, electron temperature and ion temperature. Also in this case STREAM and DYON results agree well, with a slight deviation in the effective connection length $L_{f}$, and at early times in the deuterium confinement time $\tau _{D}$, possibly explaining differences between the two codes. The differences in $L_{f}$ and $\tau _{D}$ could in turn possibly be explained by differences in the evolution of the plasma volume and the plasma major radius. The plasma volume was determined from figure 10(a) of Kim et al. (Reference Kim, Mineev, Ricci, Lee and Na2020), while the plasma major radius was assumed constant. The plasma volume plays an important role for the density and temperature evolution, while the plasma major radius affects the value of the toroidal magnetic field used in the simulation, all of which in turn affect the effective connection length and the particle confinement time.

Figure 3. Time evolution of plasma parameters for the JET carbon wall discharge $\#77\,210$: (a) plasma current; (b) electron density;(c) electron temperature; (d) effective connection length; (e) ion temperature; and (f) confinement time. The input parameters are given in table 2. The solid blue line is obtained by STREAM, dashed black line by DYON and dotted line is the measured plasma current.

4.1. ITER first-plasma scenario

We begin by considering runaway generation in an ITER ohmic first plasma scenario, with the parameters given in the first column of table 2 and a deuterium–deuterium recycling coefficient $Y_{D}^{D}=1$. Note, these parameters were chosen for benchmarking purposes by Kim et al. (Reference Kim, Mineev, Ricci, Lee and Na2020), and they are quantitatively different from an actual plasma initiation scenario in ITER. In this baseline scenario, the plasma volume, position and loop voltage are held constant while in reality they should be evolving in time, and no impurities are considered. Furthermore, the vessel volume is assumed to be as in ITER with a fully completed first wall, i.e. much smaller than the vessel volume in 2025 ITER first plasma. Nevertheless, since we are primarily interested in qualitatively studying the prevalence of runaway electrons in start-up scenarios, we use these parameters to reduce the complexity of the simulations.

We focus on the tokamak start-up from just after breakdown when $0.2\,\%$ of deuterium atoms are ionized, the plasma current has reached $I_{p}=2.4\,\textrm {kA}$ and the electron temperature is $T_{e}=1\,\textrm {eV}$. In this situation, closed flux surfaces have yet to be formed and deuterium burn-through has not yet occurred. figure 4(a,d,g,j) illustrates the evolution of plasma current, temperature and electric field. Burn-through is achieved within the first $100\,\textrm {ms}$, after which the temperature continues to increase and the plasma current is ramped up linearly. In this case, no appreciable number of runaway electrons are generated, despite the fact that $E/E_{D}$ takes values as high as 85 % early in the start-up, and remains well above the threshold value $E_{c}$ throughout the start-up phase. The reason for this is the strong transport which expels all runaways that are generated before closed flux surfaces have formed. By the time the flux surfaces are closed, the free electron density and temperature have risen sufficiently for the ratio $E/E_{D}$ to be negligibly small. The subsequent increase in $E/E_{D}$ is due to the rising temperature, but never reaches values large enough for significant Dreicer generation to occur.

Figure 4. Time evolution of plasma parameters for the ITER first plasma scenario. Panels (a,d,g,j) show the baseline case, when no significant runaway current is generated. Panels (b,e,h,k) show a case when a significant runaway current is generated. Panels (c,f,i,l) show the same case as panels (b,e,h,k), but with runaway generation turned off. Here, $I_{p}$ denotes the total plasma current, $I_\textrm {re}$ the runaway electron current and $I_{\varOmega }$ the ohmic current.

By lowering the prefill pressure $p_\textrm {prefill}$ by a factor of 10, to $0.08\,\textrm {mPa}$, so that the initial electron density is an order of magnitude lower than in the baseline case, we obtain a scenario in which burn-through succeeds and significant runaway electron generation occurs. In figure 4(b,e,h,k) the evolution in this lower-prefill case is illustrated. As in the baseline case (figure 4a,d,g,j), the plasma current rises at the expected rate, but the composition of the current is very different. The ohmic current saturates at a value of approximately $1\,\textrm {MA}$ while the current carried by the runaway electrons continues to rise unhindered. As with the ohmic current, the electron temperature also levels off since most of the energy transferred via the electric field is deposited to the runaway electrons which, due to their low rate of collisions, retain most of the energy without transmitting it to the thermal electrons. This also causes $E_{D}$ to reach a steady value, and since the externally applied loop voltage is held constant, the parameter $E/E_{D}$ saturates. It should be pointed out than in ITER, the loop voltage will eventually be forced to decrease, thus potentially greatly reducing the runaway generation. As long as $E>E_{c}$, runaways can continue to reproduce through the avalanche mechanism, requiring that the model be coupled to a simulation of the current in surrounding coils for a quantitative study.

Figure 5(a) shows the role played by the Dreicer and avalanche generation mechanisms, respectively, in the low-prefill pressure case. Throughout the simulation, Dreicer generation remains the dominant mechanism while avalanche only starts producing a noticeable number of runaways after a few seconds, when a sufficient seed population has built up. The reason for the strong Dreicer production is the large value of $E/E_{D}$ which is sustained for a long time. This is in contrast to the typical situation in disruptions where the electric field usually attains significant fractions of $E_{D}$ only for a brief period during the current quench, after which it quickly drops and approaches $E_{c}^\textrm {eff}$, driving further runaway production mainly via avalanche multiplication.

Figure 5. Time evolution of runaway generation and runaway electron confinement time for case B.

Transport plays an important role for the runaway suppression only at very early times, before closed flux surfaces have formed (which occurs when the current reaches $I_{p}\approx 100\,\textrm {kA}$ at approximately $t\approx 40\,\textrm {ms}$). Beyond this point, confinement rapidly improves and enables net runaway production. As shown in figure 5, the losses of runaways due to drifts remain negligible throughout the simulation, even as drift losses become the dominant loss mechanism.

In figure 4(c,f,i,l) we show the evolution of the same case as in figure 4(b,e,h,k), but with all runaway generation disabled (i.e. setting $\partial n _\textrm {re}/\partial t=0$ in (2.26)). In this case, burn-through succeeds and the plasma current and temperature are effectively ramped up. The large ratio $E/E_{D}$ is, however, a major warning sign, as it quickly rises beyond $E/E_{D}=5\,\%$ where significant runaway generation would generally be expected. This illustrates the need to not just ensure that burn-through is reached for successful start-up, but also to verify that $E/E_{D}$ remains sufficiently low for negligible runaways to be produced. Note that $E/E_{D}$ should be a better parameter to track in start-up simulations, rather than the threshold parameter $E/E_{c}$, due to the relative importance of the Dreicer mechanism in these scenarios, which is exponentially sensitive to $E/E_{D}$.

One of the most crucial uncertain parameters of the simulation is $\tau _\textrm {re}$, the runaway electron confinement time. In § 2.4.1 we derived a heuristic model for the runaway confinement time, and as observed in the simulations above, it is short enough to completely eliminate the seed runaways produced in the burn-through phase. The model is, however, not validated against experiment, and it could therefore be informative to investigate how sensitive our results are to variations in this parameter. To this end, we introduce a scale factor $\alpha _\tau$ in the transport term of the runaway density equation (2.26), such that $\tau _\textrm {re}\to \alpha _\tau \tau _\textrm {re}$. The result of varying $\alpha _\tau$ between $10^{-3}$ and 1 is shown in figure 6, illustrated by the time evolutions of the runaway electron density $n_\textrm {re}$ for different $\alpha _\tau$ and runaway loss fraction

(4.1)\begin{equation} f_{{\rm re},{\rm loss}}= \frac{\int_0^{t_{\rm max}}\,\mathrm{d} t\left(\gamma_{\rm Dreicer}+\varGamma_{\rm ava}n_{\rm re}\right)} {\int_0^{t_{\rm max}}\,\mathrm{d} t\,n_{\rm re}/\tau_{\rm re}}, \end{equation}

where $t_\textrm {max}=0.15\,\textrm {s}$ was used for these simulations. As previously observed, with $\alpha _\tau =1$ (transport exactly according to the heuristic model) all runaway electrons are lost from the plasma. This remains true for scale factors down to approximately $\alpha _\tau \approx 0.2$, where some runaway electrons are able to survive until the closed flux surfaces form. With $\alpha _\tau \approx 10^{-3}$, effectively all runaways generated in the early phase survive the simulation.

Figure 6. (a) Time evolution of the runaway electron density $n_\textrm {re}$ for different values of the runaway transport scale factor $\alpha _\tau$. (b) Runaway electron loss fraction $f_{\textrm {re},\textrm {loss}}$, defined according to (4.1), as a function of the scale factor $\alpha _\tau$.

4.2. Runaway electron generation after gas fuelling

Since the results of the previous section show that the most crucial parameter for the generation of runaways is $E/E_{D}$, we can consider this parameter to try to understand how runaway electrons could potentially be prevented. The quantity $E/E_{D}$ depends on three physical parameters, namely the electric field strength $E$, the electron temperature $T_{e}$ and the electron density $n_{e}$. Since the goal of the start-up is to reach a target plasma current and temperature, limiting $T_{e}$ can only be a temporary measure and is as such unfeasible for preventing runaways. Limiting $E$ can be effective in preventing the formation of runaways, but it will also limit the rate at which the current and temperature can be increased. This leaves the electron density as the main control parameter. Since $E/E_{D}\propto n _{e}^{-1}$, increasing the electron density (and thereby the electron collision frequency) could potentially be used to limit the growth of runaway electrons.

As was already illustrated by case A in § 4.1, a high prefill pressure can be effective in limiting the runaway electron growth. However, too high a prefill pressure will also prevent successful burn-through. Since persistent runaways tend to be generated during later stages of start-up, after burn-through has occurred, one could imagine a start-up scenario in which the prefill pressure is kept low to guarantee burn-through, with the density being subsequently increased to lower, or at least maintain, $E/E_{D}$ at safe levels.

In figure 7, we consider the same scenario as case B of § 4.1, but with the neutral influx (2.11) modified to include a source function $S_{i,\textrm {fuel}}^{(0)}$:

(4.2)\begin{equation} \varGamma_{i,{\rm in}}^{(0)} = V_{p}\sum_k\sum_{l\geq 1} \frac{Y_k^i n_k^{(l)}}{\tau_k} + \frac{\hat{V}_i^{(0)}}{V_i^{(0)}} S_{i,{\rm fuel}}^{(0)}. \end{equation}

Here we take the source function for deuterium to be a box function

(4.3)\begin{equation} S_{D,{\rm fuel}}^{(0)}(t) = n_{D,0}\begin{cases} 1, & t_0\leq t\leq t_0+\Delta t,\\ 0, & \text{otherwise}, \end{cases} \end{equation}

with $n_{D,0}=5\times 17\,\textrm {m}^{-3}\,\textrm {s}^{-1}$, $t_0$ is the activation time and $\Delta t$ is the duration of the source. As shown in figure 7(a), in this study we only vary the onset $t_0$ of the source (i.e. its magnitude and duration are kept constant), and it results in the same number of injected electrons in all cases. In this section we consider injections starting at $0.5\,\textrm {s}$, $1\,\textrm {s}$ and $3\,\textrm {s}$ after breakdown, all with a duration of $2\,\textrm {s}$.

Figure 7. Evolution of (a) electron density, and (b) electron temperature when neutral deuterium is continuously injected for a period of two seconds. The initiation of the fuelling is varied between the three different cases, and it is started at $0.5\,\textrm {s}$ (black), $1.0\,\textrm {s}$ (red) and $3.0\,\textrm {s}$ (blue), respectively. The shaded regions correspond to the intervals during which the fuelling is active in the different cases.

Figure 8(a) shows the evolution of the total plasma current $I_{p}$ and runaway current $I_\textrm {re}$. While $I_{p}$ remains almost exactly the same in all three cases, the fraction of the current carried by runaways differs immensely. This also implies that the ohmic component of the current differs greatly, which will affect the ability of the electric field to heat the plasma. In the case with late fuelling ($3-5\,\textrm {s}$, blue), significant Dreicer generation occurs which causes runaway electrons to carry ${\sim} 75\,\%$ of the total current by the time the fuelling is initiated. The parameter $E/E_{D}$ plateaus near $E/E_{D}\approx 6\,\%$ due to the poor ohmic coupling, allowing the plasma current to continue to increase while the temperature and ohmic current remain constant. Once fuelling commences, the Dreicer generation is suppressed as illustrated in figure 8(c), but due to the large seed of runaways already created, and the difficulty of raising the threshold electric field $E_{c}^\textrm {eff}$ sufficiently high, the runaways can still multiply via the avalanche mechanism and raise the runaway current further.

Figure 8. Evolution of (a) plasma and runaway current, (b) electric field and (c) runaway generation rate $\mathrm {d}n_\textrm {re}/\mathrm {d} t$ in the gas fuelling scenarios. By activating the fuelling earlier, the quantity $E/E_{D}$ can be kept down to suppress the Dreicer generation, effectively delaying the growth of the runaway electron beam.

When the fuelling source is initiated earlier, at $1\,\textrm {s}$ after breakdown, the effect on the runaway current is significant. Again, the fuelling is successful in suppressing the Dreicer generation, but in contrast to the delayed fuelling case, in this case the runaway seed is much smaller. As a result, avalanche multiplication is much slower and gives the ohmic current time to rise without runaway electrons overtaking the current. Eventually, after the gas fuelling is finished and the temperature has risen further, $E/E_{D}$ recovers somewhat and again allows Dreicer generation to produce more runaways.

By initiating the fuelling source even earlier, just $0.5\,\textrm {s}$ after breakdown, the effect on the Dreicer generation becomes even more pronounced. In this case, almost all Dreicer runaways are suppressed, preventing them from multiplying. Also here, $E/E_{D}$ eventually recovers and allows the Dreicer generation to commence again, but only after the ohmic plasma current and temperature have been raised significantly. Adding to this the fact that, in reality, the loop voltage will have to decrease significantly after burn-through (de Vries & Gribov Reference de Vries and Gribov2019), $E/E_{D}$ (and by extension the runaway generation) will be much lower at later times than in our simulations.

From the above analysis we can conclude that gas fuelling may be effective in preventing significant runaway generation. Our simulations suggest that the timing of the fuelling is crucial, so that a higher density is reached as soon as possible after burn-through in order to reduce the parameter $E/E_{D}$ to which the Dreicer generation mechanism is exponentially sensitive. If the density can be raised sufficiently high early on during the start-up, the potential avalanching of runaways can be greatly reduced and bring the plasma into an essentially runaway-free operating state.