Runaway electron generation during tokamak start-up

Abstract

Tokamak start-up is characterized by low electron densities and strong electric fields, in order to quickly raise the plasma current and temperature, allowing the plasma to fully ionize and magnetic flux surfaces to form. Such conditions are ideal for the formation of superthermal electrons, which may reduce the efficiency of ohmic heating and prevent the formation of a healthy thermal fusion plasma. This is of particular concern in ITER where engineering limitations put restrictions on the allowable electric fields and limit the prefill densities during start-up. In this study, we present a new 0D burn-through simulation tool called STREAM (STart-up Runaway Electron Analysis Model), which self-consistently evolves the plasma density, temperature and electric field, while accounting for the generation and loss of relativistic runaway electrons. After verifying the burn-through model, we investigate conditions under which runaway electrons can form during tokamak start-up as well as their effects on the plasma initiation. We find that Dreicer generation plays a crucial role in determining whether a discharge becomes runaway-dominated or not, and that a large number of runaway electrons could limit the ohmic heating of the plasma, thus preventing successful burn-through or further ramp-up of the plasma current. The runaway generation can be suppressed by raising the density via gas fuelling, but only if done sufficiently early. Otherwise a large runaway seed may have already been built up, which can avalanche even at relatively low electric fields and high densities.

1. Introduction

Operation in future fusion devices will have to ensure a start-up scenario that is reliable, reproducible and low-risk. To achieve that, a solid understanding of all aspects of plasma initiation is needed. If the operating parameters are not chosen with care, plasma initiation can sometimes lead to the formation of a beam of superthermal electrons, sometimes referred to as runaway electrons (Knoepfel & Spong 1979).

Runaway beam formation is usually a significant concern in tokamak disruptions, when the temperature falls and induces a large electric field that accelerates electrons (Breizman et al. 2019). Runaway acceleration was frequently observed during plasma initiation in the early days of tokamak operation, but has received less attention lately, as most current tokamaks can tune the operational parameters so that runaway electron discharges are avoided. However, it is not clear how the parameters should be chosen for future reactor-scale tokamaks, like ITER. To avoid the risk of failed start-up and corresponding delays or possible damage to in-vessel components caused by runaway electrons, careful investigation of the plasma initiation is needed.

In future tokamaks, the electric field applied for ionization and to ramp up the plasma current is limited due to engineering constraints related to the superconducting magnetic coils. In ITER, the maximum electric field available for plasma breakdown is $0.3\,{\rm V}\,{\rm m}^{-1}$ (Gribov et al. 2007). This is considerably lower than the typical value for the electric field available in current devices (${\sim }1\,{\rm V}\,{\rm m}^{-1}$). Tokamak start-up using a weak electric field requires a low prefill gas pressure (de Vries & Gribov 2019). Due to the low collisionality at low prefill pressure electrons are more prone to accelerate to high energies. Therefore, plasma initiation using low prefill pressure is often associated with discharges dominated by runaway electrons.

Tokamak plasma initiation consists of the breakdown phase, the burn-through phase and the subsequent ramp-up of current. Plasma breakdown by electron avalanche is induced by an applied electric field (Lloyd et al. 1991). In the breakdown phase, losses are dominated by transport along magnetic field lines. After plasma breakdown, the next phase is the burn-through, when ionization continues, as long as the heating (that can consist of both ohmic and auxiliary heating) overcomes losses due to ionization and radiation from the fusion fuel and impurities (Kim, Fundamenski & Sips 2012).

Burn-through allows the plasma to reach full ionization, high temperature and low resistance. After this phase, the plasma current can be ramped-up efficiently by an applied electric field. As the plasma current increases closed flux-surfaces are formed. Burn-through is successful if the plasma reaches high enough ionization for the heating power provided to be greater than the power lost due to ionization, radiation and transport. The runaway electron content strongly affects the efficiency of the burn-through and the electron density development (de Vries et al. 2020). Runaway electrons may prevent further increase in temperature by impairing the ohmic coupling between the electric field and plasma, and thereby hinder the ramp-up of the plasma current.

Plasma initiation is an inherently dynamic situation in which many plasma parameters (temperature, density, current, electric field, etc.) evolve simultaneously and depend on each other, and is often described using a set of coupled, nonlinear differential equations. The situation is complicated by the fact that the presence of runaway electrons affects the plasma conductivity and the ionization rate of atoms, hence altering the evolution of the plasma parameters. It is therefore essential that start-up modelling tools include runaway electron physics. However, most of the previously available numerical solvers for the tokamak start-up problem do not include the effect of energetic electrons (Kim et al. 2020). The code currently used for ITER start-up development – the SCENPLINT code (Belyakov et al. 2003a, b) – does include effects of runaway electrons, but it uses simplified models for the runaway generation rates.

In this paper we investigate the effect of runaway electrons on plasma initiation. In § 2 we present the theoretical model underlying the simulation tool STREAM¹ (STart-up Runaway Electron Analysis Model) that includes the essential physical processes necessary for investigating start-up scenarios. The tool STREAM computes the runaway electron generation self-consistently with the plasma density, temperature and ion charge-state evolution, as well as the electric field evolution, models for coupling to the conducting structures in the wall, and transport along magnetic field lines. In § 3 we show that the burn-through model in STREAM is in excellent agreement with the predictions of DYON (Kim et al. 2012), both for ITER and JET parameters, when runaway electrons are not considered. In § 4 we demonstrate the effect of runaway electrons on ITER burn-through scenarios. Finally, in § 5 we summarize our conclusions.

2. Burn-through model

During the burn-through and ramp-up phases of tokamak start-up the applied electric field is responsible for driving an ohmic plasma current and thereby heating the plasma. As the temperature evolves, the resistivity changes, which impacts the plasma current. The atoms in the plasma chamber are initially primarily neutrals, but as the burn-through proceeds, they rapidly ionize and contribute to radiation loss processes, significantly affecting the temperature evolution. Superthermal electrons can be accelerated during many stages of the start-up, and their generation depends sensitively on background plasma parameters such as the electric field strength, electron temperature and density.

Building on the self-consistent disruption simulation tool DREAM (Hoppe, Embreus & Fülöp 2021), a new tool STREAM has been developed. The new tool introduces some physics essential to tokamak start-up which is not included in DREAM, and implements a 0D plasma model similar to the DYON code (Kim et al. 2012). Specifically, STREAM evolves the quantities listed in table 1, including the temperatures of electrons and an arbitrary number of ion species, the thermal, runaway electron and ion charge state number densities, as well as the parallel electric field and the plasma current through the ohmic and runaway electron current densities. In the following section we describe the main elements of the model implemented in STREAM.

Table 1. Background plasma quantities evolved by STREAM.

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 (1996). 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. (2021). 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. 1996)

(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 2004).

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$,² 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. 2012)), 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 2013).

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. 2012)

(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. (2012) 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. 2014; Kim et al. 2020),

(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 2020) 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 1997)

(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. 2012), 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. (2021), 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 1975)

(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 1959), 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 1993). 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 2015; Breizman et al. 2019), 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 2017; Fülöp et al. 2020; Vallhagen et al. 2020).

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 1959), $\varGamma _{\rm ava}$ is the avalanche growth rate by which runaway electrons exponentially multiply (Rosenbluth & Putvinski 1997) 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. (2019b), while the avalanche growth rate is given by the semianalytical formula (Hesslow et al. 2019a)

(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. (2018b). The effective critical electric field $E_{c}^\textrm {eff}$ is calculated as described in appendix C.2 of Hoppe et al. (2021), and takes into account the effect of bremsstrahlung and synchrotron radiation based on the approach of Hesslow et al. (2018a). 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. (2019b) 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. (2017)), as demonstrated in Hesslow et al. (2019a) and Vallhagen et al. (2020).

As for ions, the runaway electron confinement time varies significantly during start-up (Kavin, Mineev & Lobanov 2017). 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 (1979) 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. 2020).

3. Verification of the burn-through model

The plasma initiation model presented in § 2 involves coupled equations, most of which are non-trivial in their numerical implementation. To verify that the implementation of STREAM is correct, we must compare its predictions with previously established results. The runaway electron generation models included in STREAM are inherited from DREAM, for which extensive benchmarks to previous results have been made (some of which are documented in Hoppe et al. (2021)).

For STREAM, the crucial part to verify is the burn-through model, excluding the runaway electron physics. Such models have been implemented in several codes before, including in the DYON (Kim et al. 2012), SCENPLINT (Belyakov et al. 2003a, b) and BKD0 (Granucci et al. 2015) codes, which were carefully benchmarked against each other recently in Kim et al. (2020). In this section we will reproduce the results of the two burn-through scenarios considered in Kim et al. (2020) with STREAM. In the first scenario, a pure hydrogen plasma start-up in an ITER-like setting is considered, assuming a perfectly insulating vacuum vessel. In the second scenario, a more advanced JET case is studied, involving multiple impurity species and a conducting structure which affects the current and electric field evolution.

Although the burn-through models in STREAM and DYON are similar regarding all the important physics components: the particle and heat balance, the electrical circuit model, the neutral screening and the impurity charge-state evolution, there are slight differences. While DYON assumes all the ions to have the same temperature, in STREAM we allow for different ion species to have different temperatures and only assume that the temperature is the same in every charge state. Furthermore, there are also differences in how the conductivity is calculated as described in § 2.3. However, as we shall see, in spite of these subtle differences, the agreement between the results of the two codes is very good.

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. (2020). 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. (2020).

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. 2020), 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. (2021) 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. (2020), 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. (2020). 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. (2020).

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. (2020), 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. ITER burn-through simulations with runaway electrons

The role of runaway electrons during tokamak start-up has been studied both experimentally (Esposito et al. 1996; Yoshino & Seki 1997; Esposito et al. 2003; de Vries et al. 2020) and theoretically (Sharma & Jayakumar 1988; de Vries & Gribov 2019), and it is well known that the presence of superthermal or runaway electrons during start-up can lead to burn-through failing or the creation of a relativistic electron beam. Here, we use ‘superthermal’ to mean electrons which are moving at speeds much faster than the thermal speed, but not undergoing free acceleration. By ‘runaway electrons’, on the other hand, we mean electrons which are sufficiently fast for the electric field acceleration to dominate collisional slowing-down, and which therefore accelerate freely. To study the dynamics of the former, one generally needs to solve a kinetic equation which accounts for the momentum dynamics of the electrons, while the latter can often be studied using conceptually simpler fluid models (see e.g. Vallhagen et al. 2020). In this paper, we will only consider runaway electron discharges and leave considerations concerning superthermal discharges for a future publication.

In this section we consider runaway electron generation in an ITER ohmic first-plasma scenario. We investigate conditions for the appearance of runaway electrons in the plasma as well as their effect on the plasma start-up.

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. (2020), 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 2019), $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.

5. Discussion and conclusions

Runaway electrons produced during plasma initiation can have a significant effect on the evolution of plasma parameters, in particular in future tokamaks, where due to restrictions on the allowed toroidal electric field, the prefill pressure needs to be kept low in order for burn-through to be achievable. The low prefill pressure, and subsequent plasma density, leads to low values for the runaway electron threshold electric field $E_{c}$, which makes significant runaway generation much more likely. Since most plasma parameters vary greatly during tokamak start-up, it is crucial to evolve the plasma parameters self-consistently in simulations of the start-up. In this study, we have therefore developed the new simulation tool STREAM to self-consistently couple models for the background plasma to models for the runaway electrons. The background plasma models have been benchmarked to the results of DYON presented in Kim et al. (2020), which were in turn compared with simulations of the same scenarios using BKD0 and SCENPLINT, and show good agreement.

In § 4.1 we considered an ITER first plasma-like scenario and varied the prefill pressure, as well as whether or not runaway electron generation was accounted for in the model. We found that a low prefill pressure would lead to significant runaway generation, which would also impact the ohmic coupling between the electric field and the plasma, thus limiting the amount of heat which could be provided to the plasma in the presence of runaways. The clearest sign of whether significant runaway generation would occur in our simulations was when $E/E_{D}$, the electric field normalized to the Dreicer field, reached values of more than ${\sim }3\,\%$ during a significant amount of time after successful burn-through. Prior to burn-through, $E/E_{D}$ usually takes extremely high values, but any runaway electrons which are generated during this phase are almost immediately lost to the wall due to the open field line configuration. Considering that significant runaway electron generation can impact the plasma properties during start-up, monitoring the evolution of $E/E_{D}$ and ensuring that it is kept sufficiently low could be a straightforward way to verify the applicability of other existing burn-through models.

The confinement time for runaway electrons is a highly uncertain parameter in our simulations, and they suggest that the result that all runaway electrons are lost is robust to within a factor of ${\sim }10$ in the confinement time. This assumes that the generation is dominated by relativistic electrons, something which is not necessarily the case in an experimental setting. It is possible that electrons only attain moderate superthermal energies, likely leading to improved confinement of these electrons which are prone to continue accelerating and eventually turning into relativistic runaway electrons.

We also studied the effect of gas fuelling on the generation of runaway electrons and found that injection of deuterium after burn-through can be effective in suppressing much of the runaway generation. The large values of $E/E_{D}$ obtained early on in low prefill scenarios gives rise to significant Dreicer generation. This provides a large seed of runaway electrons which will avalanche even after $E/E_{D}$ has been reduced by the increased density, thus allowing a significant runaway current to build up. By fuelling the plasma early after burn-through, the Dreicer generation can be suppressed, leaving little or no runaway electrons available to avalanche when $E/E_{D}$ eventually recovers. By increasing the density even further than we have done in this paper it might also be possible to raise the runaway threshold $E_{c}$ above the applied electric field $E$, thus also suppressing the avalanche multiplication. However, note that if fuelling is increased too much or too early, plasma burn-through will fail for the same reason as high prefill gas pressure cases. This calls for careful optimization of the prefill gas and fuelling to avoid the risk of both runaways and failed burn-through.

It is important to note that our studies are qualitative rather than quantitative, and that we have not considered the response of the tokamak control system in these simulations. Since the plasma resistivity decreases with temperature, the electric field should eventually be limited to maintain a constant plasma current when the target temperature is reached. When this happens, it should limit the avalanche generation mechanism, stop further runaway growth and instead reduce the number of runaways, if any, via collisions.

The focus in § 4.2 has been on preventing the generation of runaway electrons. A natural follow-up question to ask would be if it is possible to suppress already existing runaways. In our model, the only reasonable means for reducing the number of runaway electrons is to reduce the threshold parameter $E/E_{c}^\textrm {eff}$ below unity to allow the runaways to avalanche ‘in reverse’ and collisionally dissipate their energy to the bulk electrons. This might be feasible, and a natural result of successful start-up, if the fraction of runaway current is low. If, however, the runaway current is significant, this process would be difficult to reconcile with maintaining the plasma at a specified current and temperature.

In this study we have only considered the effect of relativistic electrons. During start-up it is also possible for the electron distribution function to be significantly distorted from thermal equilibrium due to the large values of $E/E_{D}$, and for a large fraction of superthermal (albeit not relativistic) electrons to form. Such electrons will also affect the heating properties and general evolution of the plasma, but are not captured by the model presented here. An analysis of these electrons would require the solution of the Fokker–Planck equation with nonlinear collision coefficients, since linearizations of the collision operator typically assume that most electrons are in thermal equilibrium. While DREAM, and thus by extension STREAM, supports the solution of a Fokker–Planck equation simultaneously with the evolution of the (fluid) background plasma, all collision operators available in the code are linearized. A kinetic treatment of the runaway problem during start-up could also allow us to study the effect of runaways on the ionization of atoms, as well as the effect of electron cyclotron heating – the latter which is anticipated to be used in ITER to assist start-up – for the generation of fast electrons. This could be done by, for example, coupling STREAM to the bounce-averaged Fokker–Planck code LUKE (Decker & Peysson 2004) which also solves for wave–particle interactions. We therefore leave for future studies a more detailed analysis of the momentum–space dynamics of fast electrons during start-up.