Inelastic collisions facilitating runaway electron generation in weakly-ionized plasmas

Introduction.– Runaway Electrons (REs) are electrons that are accelerated freely by the electric field due to a decreasing drag force with increasing velocity [1]. They are important not only in astrophysics [2, 3] and geophysics [4, 5] but also in tokamaks, where RE generation is undesirable due to the danger of triggering plasma instabilities [6, 7] and ability to damage the device [8, 9, 10]. In tokamaks REs can be generated during both plasma initiation (startup) [11, 12, 10] and disruptions [13, 14, 15]. Startup REs are more critical in ITER than in the existing fusion devices because the low plasma density required for plasma burn-through facilitates their generation [16, 17, 18, 19, 20, 21] and the large device size ensures their good confinement [10]. Development of a successful startup scenario for ITER and other future reactors thus requires a robust understanding of RE physics and reliable predictions of their generation rates in cold weakly ionised plasmas.

Several RE generation mechanisms are distinguished based on how electrons surpass the critical momentum $p_{c}$ (force-free) in phase space [22]. The most important during startup is the Dreicer mechanism, pioneered in Ref. [23] and elaborated in subsequent works [24, 25, 26]. According to this mechanism, a diffusive flux from the Maxwellian core across the critical momentum boundary is formed under the effect of electric field if the electric field exceeds the critical value, $E_{c}\equiv\frac{e^{3}n_{e}\ln\Lambda}{4\pi\varepsilon_{0}^{2}mc^{2}}$, i.e. $E/E_{c}>1$. Owing to the dominance of small angle collisions in a fully ionised plasma this process is diffusive and is adequately described by the Fokker-Planck model [27, 28].

In a weakly ionized plasma, such as a tokamak plasma during startup, electron-atom collisions need to be considered [29, 30, 31, 32, 33]. The effect of the weakly ionized plasma on highly energetic electrons may still be appropriately described with the Fokker-Planck operator, since the energy transfer in a typical small angle collision exceeds the ionization potential for highly energetic electrons [22]. Then, their stopping power is typically described with the Bethe stopping power [34], and elastic scattering can be accounted for with the help of the Thomas–Fermi model [35, 22]. However, in a typical startup the critical momentum may lie in the energy range which is comparable to (or exceeds only by one order of magnitude) the ionization potential, in which case the key role in determining the RE generation rate is played by the region of the phase space for which the small angle approximation is not appropriate.

Indeed, for instance in ITER startup the inductive electric field is $E\approx 0.2\ V/m$ [16, 36], the plasma temperature is increasing, but is expected to be below 10 eV during the plasma breakdown [32], where the free electron density is slowly growing and is initially extremely low. For instance, when $n_{e}=10^{17}m^{-3}$, the ratio of the inductive electric field to the Dreicer field, $E_{d}\equiv\frac{e^{3}n_{e}\ln\Lambda}{4\pi\varepsilon_{0}^{2}T_{e}}$, is $E/E_{d}\approx 0.06$, which is above the known threshold for effective Dreicer generation (0.02-0.025) [37], but the corresponding critical momentum, $p_{c}\approx(E/E_{c})^{-\frac{1}{2}}$, is only around 80 eV - a few times higher than the H ionization potential. Evidently, the non-differential nature of collisions within this energy range needs to be considered in order to correctly describe electron kinetics.

In this Letter we present numerical calculations of the Dreicer flow accounting for the non-differential nature of inelastic collisions using the Fokker-Planck-Boltzmann operator (FPB) that treats soft and hard collisions with the FP and Boltzmann operator, respectively. The double counting of collisions is avoided by an appropriate selection of the integration boundaries in the Coulomb logarithm definition for the FP operator [38]. We demonstrate a significant effect of the non-differential inelastic collisions, which can not be captured using the FP operator.

Kinetic model.– The kinetic equation for the electrons is

where $F_{e}(p,\mu)$ is the distribution function normalized such that the density $n_{e}=\int F_{e}(p,\mu)dpd\mu$, $p$ is the particle momentum normalized to $mc$ and $\mu=\cos\theta$ with the pitch angle $\theta$. The FPB operator for collisions between free electrons and hydrogen atoms consists of two parts

where $\mathcal{C}^{eH}_{FP}\{F_{e}\}$ is the FP operator and $\mathcal{B}^{eH}\{F_{e}\}=\mathcal{B}^{eH}_{iz}\{F_{e}\}+\mathcal{B}^{eH}_{ex}\{F_{e}\}$ is the Boltzman operator. The subscripts $iz$ and $ex$ denote ionization and excitation, respectively.

We introduce a soft-hard separation factor $h$ to determine if an inelastic collision is categorized as soft or hard. When the energy loss fraction of the incident electron is higher than $h$, such collisions are included in the Boltzmann operator; otherwise, they are described by the FP operator. For high $h$ such as $h=0.99$, the test particle operator is approximated as the FP operator. Following the conventional practice, we assume that a scattered electron possesses more energy than an ejected one. The slowing-down frequency for free-bound collisions comprises the soft ionization and excitation collision parts, from $1\ s\rightarrow 2\ p$ to $1\ s\rightarrow 10\ p$ transitions.

Here, $\nu_{iz}^{sti,soft}=n_{H}(\int_{0}^{W_{bnd}}dW\frac{W+B}{mc^{2}}\frac{d\sigma_{iz}}{dW})v$ and $\nu_{ex}^{1n}=n_{H}\sigma_{ex}^{1s\rightarrow np}v$, where $\frac{d\sigma_{iz}}{dW}$ is the singly differential cross section described by the Relativistic Bineary Encounter Dipole (RBED) model [39, 40] and $\sigma_{ex}^{1s\rightarrow np}$ is the total excitation cross section with $1s\rightarrow np$ transition by the BE-scaled plane wave Born approximation model [41, 42]. $W_{bnd}=\min\{\max(0,hT-B),(T-B)/2\}$, $n_{H}$ is the hydrogen atom density, $T$ is the incident electron energy, $B$ is the binding energy, $W$ is the secondary electron energy, $E^{\prime}_{1n}=E_{1n}/mc^{2}$ is the $1s\rightarrow np$ transition energy $E_{1n}$ normalized to $mc^{2}$ and $\mathcal{H}$ is the heaviside function. The $W_{bnd}$ is chosen taking into account two conditions: the absence of soft inelastic collisions for $0\geq hT-B$ and the requirement that the ejected electron energy be less than $(T-B)/2$.

The electron-hydrogen atom soft collisions are collisions against stationary targets. Their FP operator takes the form [43]

There are the deflection frequencies due to collisions with bound electrons and hydrogen atoms, $\nu_{D}^{e,eb}=(1/\gamma)\nu_{S}^{e,eb}$ and $\nu_{D}^{eH}=(e^{4}\gamma/4\pi\varepsilon_{0}^{2}m^{2}c^{3}p^{3})(\chi n_{H}I_{2}^{H}(y))$, respectively, where $\gamma$ is the relativistic factor and $I_{2}^{H}(y)$ is the second screening coefficient based on the TF theory [22], which we fit as $I_{2}(y)\approx I_{2}(y_{*})+\ln[(y/y_{*})^{4}+\exp(-I_{2}(y_{*}))^{4}]/4$ with $y=274p$ and $y_{*}=26$. The TF model is only valid for the energetic particles, we extrapolate the operator to low energies using experimental data [44, 45] and the NIST Standard Database [46, 47]. The correction factor $\chi=1+a^{el}(511000(\gamma-1))^{b^{el}}/(1+c^{el}(511000(\gamma-1))^{b^{el}})$ is found by fitting, where $a^{el}=79.837201$, $b^{el}=-1.0992754$ and $c^{el}=1.6387662$.

Figure 1 (a) shows the refined scattering frequency $\nu_{D}^{eH}$ as a function of the incident electron energy. It is consistent with the TF model at high energy and approximates the experimental data in the low energy region. The soft collisions slowing down frequency $\nu_{S}^{e,eb}$ (Fig. 1 (b)) agrees with the Bethe model for high energy particles. In low energy range, the soft collisions represent only a fraction of all collisions and the frequency is therefore lower. In $h=0.1$ case, for instance, $\nu_{S}^{e,eb}=0$ and $\mathcal{B}^{eH}$ treats all inelastic collisions below $102.043\ eV$, this energy corresponds to $1/h$ of the minimal excitation energy for H. We assume that hydrogen is in the ground state. This treatment allows the electrons that didn’t experience inelastic collisions (accounted by $\mathcal{B}^{eH}_{iz}$) not to lose any energy during acceleration to $102.043\ eV$, the ‘‘continuous’’ soft collisions energy loss is significantly reduced even above $102.043\ eV$. Such accounting of non-differential nature of inelastic collisions changes dramatically the character of critical energy crossing, i.e. it allows some electrons to accelerate completely freely up to a certain energy.

For $\mathcal{B}^{eH}_{iz}$, we symmetrize the collisional cross section to address the asymmetry in the RBED model originating from its dipole contribution, i.e. $(d\sigma_{iz}/dW)(W,T)=(d\sigma_{iz}/dW)(T-W-B,T)\ if\ W>(T-B)/2$. In addition, a scattering angle distribution is obtained by interpolating the isotropic and Moller scattering distribution [49] as low and high energy asymptotes, respectively, since pitch-angle equilibration is rapid in the low energy region and free-bound collisions can be approximated as free-free collisions in the high energy region. After gyro-averaging,

Here, $\mu^{\ast}=\sqrt{\frac{(\gamma_{1}+1)(\gamma-1)}{(\gamma_{1}-1)(\gamma+1)}}$, $p_{bnd}$ is the momentum of $\min\{\max(B+W,W/(1-h)),2W+B\}$, essential for the birth of secondary electron with $p$, $\nu_{iz}^{hard}=n_{H}(\int_{W_{bnd}}^{(T-B)/2}dW\frac{d\sigma_{iz}}{dW})v$ and $1/(1+0.1/p)$ is an interpolation factor.

For $\mathcal{B}^{eH}_{ex}$, we assume the scattering angle to be isotropic:

where $p^{+}_{1n}=p(T+E_{1n})$ and $F^{0}_{e}=\int d\mu F_{e}(p,\mu)$.

The full operator $\mathcal{C}^{eH}\{F_{e}\}$ retains invariance in particle and energy loss of electrons with respect to $h$, i.e. $\frac{\partial}{\partial h}[\int dpd\mu\ \mathcal{C}^{eH}\{F_{e}\}]=\frac{\partial}{\partial h}[\int dpd\mu\ (\gamma-1)\mathcal{C}^{eH}\{F_{e}\}]=0$:

where $\nu_{iz}=n_{H}(\int_{0}^{(T-B)/2}dW\frac{d\sigma_{iz}}{dW})$ and $B^{\prime}=B/mc^{2}$.

The Coulomb collisions with free electrons and ions, $C^{e,ei}$ in Eqn. 1, are accounted for using the linearized FP operator as formulated in Ref [50] with the Coulomb logarithms listed in [22]. This collisional integral assumes that the majority of free electrons are Maxwellian. This can be guaranteed if $n_{e}\ln\Lambda_{free}\gg n_{H}\ln\Lambda_{b}$, with $\ln\Lambda_{free}$ describing free-free electron collisions and $\ln\Lambda_{b}$ free-bound collisions. Because $\ln\Lambda_{b}$ is vanishingly small below $1s\to 2p$ transition energy, this condition is satisfied and we therefore expect that the lower energy part of the population is predominantly Maxwellian.

Numerical simulation.– We implement the FPB operator, given by Eqn. (2), in the self-consistent kinetic simulation [51] using FiPy Finite Volume Partial Differential Equation (PDE) Solver [52].

The non-differential nature of inelastic collisions facilitates electron acceleration in weakly-ionized plasma. This is demonstrated in Fig. 2, where the steady-state solutions are shown. For the demonstration, we only took the test particle part to eliminate RE avalanche effect [53] by setting upper boundary of integral in Eqn. (5) to $p$ corresponding to $2W+1$. The blue curve in Fig. 2 shows the solution where most collisions are described with the FP operator (h=0.99). The effective $E/E_{d}^{eff}<0.02$, with $E_{d}^{eff}\equiv\frac{e^{3}(n_{e}\ln\Lambda+n_{H}\ln\Lambda_{b})}{4\pi\varepsilon_{0}^{2}T_{e}}$ ($\ln\Lambda_{b}\approx 5$) and the resulting $F_{e}^{0}<10^{3}m^{-3}$ with electron energy above $100\ eV$. Indeed, no Dreicer flow is expected in such conditions. The red dashed curve shows the solution with the FPB operator ($h=0.1$). This solution is close to Maxwellian in lower energy where free-free collisions dominate. But it differs significantly at higher energy where binary Boltzmann collisions allow for virtually free acceleration for some particles. These freely accelerated particles become runaways.

Dreicer generation rate as a function of the ionization fraction $\alpha=n_{e}/(n_{e}+n_{H})$ is shown in Fig. 3. The parameters for this scan correspond roughly to standard Ohmic discharge in KSTAR tokamak [32, 36] and are similar to that of ITER plasma initiation [17]. In this case, a critical role played by the non-differential nature of collisions is expected for $E/E_{d}>0.025$ but $E/E_{d}^{eff}<0.02$. The Dreicer generation without accounting for the non-differential nature of collisions is ineffective for these parameters. However, it reaches $10^{13}\ m^{-3}s^{-1}$ for $\alpha=0.01$ in the case of Boltzmann collisions. This is consistent with the analytical Dreicer rates since $E/E_{d}\approx 0.43\gg 0.02$ but $E/E_{d}^{eff}\approx 0.015$. The observed rates are much lower than the analytic prediction by Connor-Hastie [26] if collisions with bound electrons are ignored in the analytic calculation (green curve in Fig. 3.). We also note that that our calculations using FP operator (blue squares) agree with the Connor-Hastie formula if $E_{d}^{eff}$ with $\ln\Lambda_{b}\sim 5$ deduced from Eqn. (3) is used in the formula (dashed curve).

Dreicer generation rate is known to be exponentially sensitive to the value of the electric field [24]. Yet, strong effect of the non-differential nature of collisions should be considered separately: the sensitivity to the electric field is mainly the result of the exponential dependence of the distribution function on $p_{c}$, wheres in the case of non-differential collisions ($\mathcal{B}^{eH}$) friction-less acceleration ‘‘through’’ $p_{c}$ of some electrons facilitates the Dreicer flow. The effect of the $\pm 10\%$ deviation in $E$ is demonstrated with the error bars in Fig. 3.

The Closed Flux Surface Formation (CFSF), before which we assume a complete loss of fast particles, follows the main-speceis burn-through in small-size tokamaks [32]. Nevertheless, enhanced Dreicer generation in weakly ionized plasma can have significant implications for the startup scenarios in ITER [18] and other current [17, 10] and future tokamaks [54, 55, 56, 57] where the flux surfaces may be formed before the complete plasma ionization.

Figure 4 shows kinetic simulation results for the RE density and current during an ITER-relevant startup [17, 10], where plasma parameters $n_{e}$, $T_{e}$ and $E$ are prescribed by the ‘‘fluid’’ startup code DYON [58, 59, 60, 36](black curves). $E$ is deduced from the resistive voltage. In this calculation only the phase after the global CFSF is simulated, i.e. when the maximum ratio of the stray field to the toroidal field is $10^{-3}$ [17] and therefore the initial plasma current $I_{p}=21.1kA$. The kinetic model considers the full FPB operator in Eqn. (2).

Dreicer generation is most effective early in the simulation, since the rapidly growing electron density $n_{e}$, which exceeds $n_{H}$ at $4\ ms$, makes the generation ineffective quickly as shown in Fig. 4. Accounting for the non-differential nature of collisions enhances the Deicer generation increasing the ultimate RE density $n_{RE}$ (4) and current $I_{RE}$ (4) by a factor of 3. In this particular example the resulting $I_{RE}$ overtakes a significant fraction ($7.5\%$) of the total current, posing a danger of a complete startup failure, due to a reduction of the Ohmic heating (not modeled self-consistently here). Indeed, simulation show that we would already have the complete current takeover if E is only slightly higher as shown in Fig. 4. This simulation is not shown in Fig. 4, but the ultimate $n_{RE}=5.6\cdot 10^{13}m^{-3}$ in this case.

Conclusions.– Reliable startup RE prediction in large-scale future tokamaks such as ITER [17, 10] and DEMOs [54, 55, 56, 57] require appropriate description of Dreicer generation during early ionization stage. In this work we have demonstrated an essential role of the non-differential nature of ionizing collision, especially when $E/E_{d}\geq 0.025$ but $E/E_{d}^{eff}\leq 0.02$. We show that such enhanced Dreicer generation may lead to a complete failure of tokamak startup if the Ohmic current is taken over by the RE current.

This research was supported by National R&D Program through the National Research Foundation of Korea (NRF) funded by Ministry of Science and ICT (2021M3F7A1084419). This research was supported by National R&D Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Science & ICT (NRF-2019R1A2C1010757). This research was supported by R&D Program of "High Performance Fusion Simulation R&D (code No. EN2441-1)"through the Korea Institute of Fusion Energy(KFE) funded by the Government funds, Republic of Korea.