Effect of partially ionized impurities and radiation on the effective critical electric field for runaway generation

When a plasma carrying a large electric current is suddenly cooled, as happens in tokamak disruptions, a large toroidal electric field is induced due to the dramatic increase of the plasma resistivity. If this electric field is larger than a certain critical electric field, a relativistic runaway electron beam can be generated [1, 2]. Such runaway beams can damage the plasma facing components on impact due to localized energy deposition. Therefore, runaway electrons constitute a significant threat to large tokamak experiments (e.g. ITER) [3–5].

To minimize the risk of damage, it is crucial to understand the runaway electron dynamics. Disruption mitigation by material injection is motivated by the strong influence of partially ionized atoms, as observed in experiments [3, 4, 6]. It is therefore important to have accurate models of the interaction between fast electrons and the partially screened nuclei of heavy ions. Fast electrons are not simply deflected by the Coulomb interaction with the net charge of the ion, but probe its internal electron structure, so that the nuclear charge is not completely screened. Energetic electrons can therefore be expected to experience higher collision rates against partially ionized impurities compared to a fully ionized plasma with the same effective charge, leading to a more efficient damping. There has been a considerable effort to produce a detailed theoretical description of this process [7–10].

A recent paper presented a generalized collision operator which describes the interaction between fast electrons and partially screened impurities via analytic modifications to the collision frequencies [9]. The elastic electron–ion collisions were modeled quantum-mechanically in the Born approximation as in [7, 8], however, to obtain the required electron density distribution of the impurity ions [7, 8] used the Thomas–Fermi model. In [9] we used fitted results from density functional theory (DFT) thereby providing a more accurate description. To describe inelastic collisions with bound electrons, we employed Bethe's theory for the collisional stopping power [11], with mean ionization energies for ions calculated in [12]. Our results show that, already at sub-relativistic electron energies, the deflection and slowing-down frequencies are increased significantly compared to standard collisional theory [9].

The quantity that is arguably the most important for runaway generation and decay is the threshold, or critical, electric field, which in a fully ionized plasma without radiation losses is given by the Connor–Hastie field ${E}_{{\rm{c}}}\,={n}_{{\rm{e}}}{e}^{3}\mathrm{ln}\,{\rm{\Lambda }}/(4\pi {\epsilon }_{0}^{2}{m}_{{\rm{e}}}{c}^{2})$ [2], where n_e and m_e are the electron density and mass, $\mathrm{ln}\,{\rm{\Lambda }}$ is the Coulomb logarithm, ₀ is the vacuum permittivity and c is the speed of light. Below the threshold field no new runaway electrons are produced and all preexisting runaways eventually thermalize. There is a wealth of experimental evidence that the critical electric field is much higher than E_c given above [13–18]. Well-diagnosed and reproducible experiments in quiescent plasmas on a wide range of tokamaks show that measured threshold electric fields can be approximately an order of magnitude higher than predicted by the Connor–Hastie threshold [13, 18]. Furthermore, it has been shown that the runaway electron current decays much faster after high-Z particle injection than expected from conventional theory [2], in contrast to low-Z particle injection which results in a current decay rate only slightly below that expected [14]. From a theoretical point of view, the threshold electric field is expected to be higher than E_c, as can be influenced by synchrotron [19, 20] and bremsstrahlung radiation losses, and also, as we will show here, by the presence of partially ionized atoms. The value of the critical electric field is not only interesting theoretically—it is of immense practical importance as it determines the amount of material that has to be injected in disruption mitigation schemes [21].

In this paper we derive an analytical expression for the effective critical field for runaway generation and decay that takes into account the presence of partially screened impurities, using the generalized collision operator derived in [9]. We present a formula that accounts for arbitrary ion species in combination with synchrotron and bremsstrahlung losses. We show that the effect of partially screened impurities is captured by replacing the plasma density in the critical electric field with an effective density $n={n}_{\mathrm{free}}+\kappa {n}_{\mathrm{bound}}$, where κ is typically in the range 1–2 which implies that the effect of bound electrons is significantly larger than suggested by previous studies [22]. Furthermore, using a kinetic equation solver with a 0 D inductive electric field, we verify the prediction from [21], that the runaway current in highly inductive tokamak devices after impurity injection will decay linearly with time at a rate proportional to the effective electric field. We expect these findings will facilitate future comparisons with experimental observations of runaway-current decay, however such analysis is beyond the scope of the present paper.

The structure of the paper is as follows. In section 2 we describe the kinetic model accounting for the effect of partial screening in both the generalized collision operator and the bremsstrahlung operator. Then we proceed in section 3 to derive analytical expressions for the effective critical electric field in the presence of partially ionized impurities. This calculation generalizes the results in [20], in which the critical electric field was calculated by assuming rapid pitch-angle dynamics in the Fokker–Planck equation. In contrast to [20], our study includes the effect of partially ionized impurities and bremsstrahlung losses. We demonstrate how the presence of partially screened impurities affects both synchrotron losses (through pitch-angle scattering) and bremsstrahlung (as partial screening affects the bremsstrahlung cross-section). In section 4 we discuss the decay of a runaway current when heavy impurities are injected. Through kinetic simulations, we demonstrate the accuracy of the analytical expressions for the effective critical electric field and the current decay. Finally in section 5 we summarize our conclusions.

In a uniform, magnetized plasma, the kinetic equation for relativistic electrons can be written as follows:

$$\begin{eqnarray}&&\displaystyle \frac{\partial f}{\partial \tau }+\mathop{\underbrace{\displaystyle \frac{E}{{E}_{{\rm{c}}}}\left(\xi \displaystyle \frac{\partial f}{\partial p}+\displaystyle \frac{1-{\xi }^{2}}{p}\displaystyle \frac{\partial f}{\partial \xi }\right)}}\limits_{\mathrm{electric}\,\ \mathrm{field}}\\ &&\quad =\mathop{\underbrace{{C}_{\mathrm{FP}}\{f\}+{S}_{\mathrm{ava}}}}\limits_{\mathrm{collisions}}+\mathop{\underbrace{{C}_{\mathrm{br}}\{f\}-\displaystyle \frac{\partial }{\partial {\boldsymbol{p}}}\cdot ({{\boldsymbol{F}}}_{\mathrm{syn}}f)}}\limits_{\mathrm{radiation}\ \,\ \mathrm{reaction}},\end{eqnarray} \tag{ 1 }$$

where f is the electron distribution function, ${C}_{\mathrm{FP}}\{f\}$ is the partially screened Fokker–Planck collision operator as described in section 2.1, which accounts for ionizing as well as elastic collisions. The avalanche source is denoted S_ava and E is the component of the electric field which is antiparallel to the magnetic field ${\boldsymbol{B}}$. Radiation losses are modeled by C_br (the bremsstrahlung collision operator) and ${{\boldsymbol{F}}}_{\mathrm{syn}}$ (the synchrotron radiation reaction force), which are described in section 2.2. The normalized momentum is defined as p = γv/c is (with γ the Lorentz factor), $\xi ={\boldsymbol{p}}\cdot {\boldsymbol{B}}/({pB})$ is the cosine of the pitch-angle, and the time variable τ is normalized to the relativistic collision time

$$\begin{eqnarray*}&&{\tau }_{{\rm{c}}}=4\pi {\epsilon }_{0}^{2}{m}_{{\rm{e}}}^{2}{c}^{3}/({n}_{{\rm{e}}}{e}^{4}\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}),\end{eqnarray*}$$

where we introduced a relativistic Coulomb logarithm

$$\begin{eqnarray}&&\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}=\mathrm{ln}\,{{\rm{\Lambda }}}_{0}+\displaystyle \frac{1}{2}\mathrm{ln}\displaystyle \frac{{m}_{{\rm{e}}}{c}^{2}}{T}\approx 14.6+0.5\mathrm{ln}\,({T}_{\mathrm{eV}}{/n}_{{\rm{e}}20}).\end{eqnarray} \tag{ 2 }$$

Here, T_eV is the temperature in eV, ${n}_{{\rm{e}}20}$ is normalized to ${10}^{20}\,{{\rm{m}}}^{-3}$ and $\mathrm{ln}\,{{\rm{\Lambda }}}_{0}=14.9-0.5\mathrm{ln}{n}_{{\rm{e}}20}+\mathrm{ln}{T}_{\mathrm{keV}}$ is the thermal electron–electron Coulomb logarithm [23]. The temperature dependence of $\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}$ is reduced compared to $\mathrm{ln}\,{{\rm{\Lambda }}}_{0}$ as it describes collisions between thermal particles and relativistic electrons; (2) corresponds to evaluating the energy-dependent electron–ion Coulomb logarithm $\mathrm{ln}\,{{\rm{\Lambda }}}^{\mathrm{ee}}$ at γ = 2. For future reference, the superthermal Coulomb logarithms are given by [24]

$$\begin{eqnarray}&&\mathrm{ln}\,{{\rm{\Lambda }}}^{\mathrm{ee}}=\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}+\mathrm{ln}\sqrt{(\gamma -1)}\end{eqnarray} \tag{ 3 }$$

and

$$\begin{eqnarray}&&\mathrm{ln}\,{{\rm{\Lambda }}}^{\mathrm{ei}}=\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}+\mathrm{ln}(\sqrt{2}p).\end{eqnarray} \tag{ 4 }$$

The parallel electric field E is thus most naturally compared to the critical electric field E_c defined with the relativistic Coulomb logarithm $\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}$ (rather than the thermal $\mathrm{ln}\,{{\rm{\Lambda }}}_{0}$):

$$\begin{eqnarray*}&&{E}_{{\rm{c}}}=\displaystyle \frac{{n}_{{\rm{e}}}{e}^{3}\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}}{4\pi {\epsilon }_{0}^{2}{m}_{{\rm{e}}}{c}^{2}}=\displaystyle \frac{{m}_{{\rm{e}}}c}{e{\tau }_{{\rm{c}}}}.\end{eqnarray*}$$

2.1. Collision frequencies with partially ionized impurities

When acting on relativistic electrons and $T\ll {m}_{{\rm{e}}}{c}^{2}$, the linearized Fokker–Planck collision operator ${C}_{\mathrm{FP}}\{f\}$ can be simplified to

$$\begin{eqnarray*}&&{C}_{\mathrm{FP}}\{f\}={\nu }_{{\rm{D}}}{\mathscr{L}}\{f\}+\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{\partial }{\partial p}({p}^{3}{\nu }_{{\rm{s}}}f),\end{eqnarray*}$$

where ${\mathscr{L}}=\tfrac{1}{2}\tfrac{\partial }{\partial \xi }(1-{\xi }^{2})\tfrac{\partial }{\partial \xi }$ is the Lorentz scattering operator. The slowing-down frequency ${\nu }_{{\rm{s}}}={\nu }_{{\rm{s}}}^{\mathrm{ee}}$ and the deflection frequency ${\nu }_{{\rm{D}}}={\nu }_{{\rm{D}}}^{\mathrm{ee}}+{\nu }_{{\rm{D}}}^{\mathrm{ei}}$ are well known in the limits of complete screening (i.e. the electron interacts only with the net ion charge) and no screening (the electron experiences the full nuclear charge). The generalized expressions for ${\nu }_{{\rm{D}}}^{\mathrm{ei}}$ and ${\nu }_{{\rm{s}}}^{\mathrm{ee}}$ taking into account partial screening are given in [9].

Focusing on the effective critical electric field E_c^eff in this paper, the following equations are specialized to the superthermal momentum region, in which the critical momentum p_c corresponding to E_c^eff is found. Thus all of the following expressions are given for superthermal electrons.

The generalized deflection frequency is, in units of ${\tau }_{{\rm{c}}}^{-1}$, given by

$$\begin{eqnarray}{\nu }_{{\rm{D}}} & = & \displaystyle \frac{\gamma }{{p}^{3}}{\bar{\nu }}_{{\rm{D}}},\\ {\bar{\nu }}_{{\rm{D}}} & = & \displaystyle \frac{1}{\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}}\phantom{\rule{}{3.5ex}}[\mathrm{ln}\,{{\rm{\Lambda }}}^{\mathrm{ee}}+\mathrm{ln}\,{{\rm{\Lambda }}}^{\mathrm{ei}}{Z}_{\mathrm{eff}}\\ & & \left.+\displaystyle \sum _{j}\displaystyle \frac{{n}_{j}}{{n}_{{\rm{e}}}}\left(({Z}_{j}^{2}-{Z}_{0,j}^{2})\mathrm{ln}({\bar{a}}_{j}p)-\displaystyle \frac{2}{3}{N}_{{\rm{e}},j}^{2}\right)\right].\end{eqnarray} \tag{ 5 }$$

Here, Z_0,j is the ionization state, Z_j is the charge number and ${N}_{{\rm{e}},j}={Z}_{j}-{Z}_{0,j}$ is the number of bound electrons of the nucleus for species j, ${Z}_{\mathrm{eff}}={\sum }_{j}{n}_{j}{Z}_{0,j}^{2}/{n}_{{\rm{e}}}$, where n_j is the density of species j, and n_e represents the density of free electrons. The parameter ${\bar{a}}_{j}$ was determined from DFT calculations, and is an effective ion size which depends on the ion species j. These constants are given for argon and neon in table A1 in appendix A. In (5), we have assumed $p\gg 1/{\bar{a}}_{j}\simeq {10}^{-2}$. Figure 1(a) shows the enhancement of the deflection frequency for singly ionized argon and neon. At typical runaway energies in the MeV range, the enhancement is more than an order of magnitude compared to taking the limit of complete screening and neglecting the variation of the Coulomb logarithm, which would give ${\bar{\nu }}_{{\rm{D}}}=1+{Z}_{\mathrm{eff}}$.

Zoom In Zoom Out Reset image size

In the limit of p ≫ 1, the deflection frequency (5) can be approximated by

$$\begin{eqnarray}&&{\bar{\nu }}_{{\rm{D}}}\approx ({\bar{\nu }}_{{\rm{D}}0}+{\bar{\nu }}_{{\rm{D}}1}\mathrm{ln}p),\end{eqnarray} \tag{ 6 }$$

where the constants are given by

$$\begin{eqnarray}&&{\bar{\nu }}_{{\rm{D}}0}=1+{Z}_{\mathrm{eff}}+\displaystyle \frac{1}{\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}}\displaystyle \sum _{j}\displaystyle \frac{{n}_{j}}{{n}_{{\rm{e}}}}\left(({Z}_{j}^{2}-{Z}_{0,j}^{2})\mathrm{ln}{\bar{a}}_{j}-\displaystyle \frac{2}{3}{N}_{{\rm{e}},j}^{2}\right),\quad \end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\bar{\nu }}_{{\rm{D}}1}=\displaystyle \frac{1}{\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}}\displaystyle \sum _{j}\displaystyle \frac{{n}_{j}}{{n}_{{\rm{e}}}}{Z}_{j}^{2}.\,\end{eqnarray} \tag{ 8 }$$

For the superthermal slowing-down frequency, we obtain, in units of τ_c⁻¹,

$$\begin{eqnarray}{\nu }_{{\rm{s}}} & = & \displaystyle \frac{{\gamma }^{2}}{{p}^{3}}{\bar{\nu }}_{{\rm{s}}},\\ {\bar{\nu }}_{{\rm{s}}} & = & \displaystyle \frac{1}{\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}}\left(\mathrm{ln}\,{{\rm{\Lambda }}}^{\mathrm{ee}}+\displaystyle \sum _{j}\displaystyle \frac{{n}_{j}}{{n}_{{\rm{e}}}}{N}_{{\rm{e}},j}(\mathrm{ln}\,{h}_{j}-{\beta }^{2})\right).\end{eqnarray} \tag{ 9 }$$

Here, ${h}_{j}=p\sqrt{\gamma -1}/{I}_{j}$ and I_j is the mean excitation energy of the ion, normalized to the electron rest energy [12]; see table A1 in appendix A. As ${\nu }_{{\rm{s}}}$ given in (9) is based on the Bethe stopping power formula matched to the low-energy asymptote corresponding to complete screening, we refer to it as the Bethe-like model. As shown in figure 1(b), the slowing-down frequency is enhanced significantly compared to the completely screened limit with constant Coulomb logarithm, where ${\bar{\nu }}_{{\rm{s}}}=1$. The enhancement is also significantly different from a widely used rule of thumb that is mentioned in passing by Rosenbluth and Putvinski [22], which suggests that inelastic collisions with bound electrons can be taken into account by adding half the number of bound electrons to the number of free electrons. As shown in figure 1, the Rosenbluth–Putvinski (RP) model overestimates the slowing-down frequency at low energies and is a significant underestimation at high runaway energies. The weak energy-dependence of the RP model is due to the energy-dependence in the electron–electron Coulomb logarithm in (3).

In the ultra-relativistic limit p ≫ 1, the slowing-down frequency  (9) is approximately

$$\begin{eqnarray}&&{\bar{\nu }}_{{\rm{s}}}\approx ({\bar{\nu }}_{{\rm{s}}0}+{\bar{\nu }}_{{\rm{s}}1}\mathrm{ln}p),\end{eqnarray} \tag{ 10 }$$

where

$$\begin{eqnarray}&&{\bar{\nu }}_{{\rm{s}}0}=1+\displaystyle \frac{1}{\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}}\displaystyle \sum _{j}\displaystyle \frac{{n}_{j}}{{n}_{{\rm{e}}}}{N}_{{\rm{e}},j}(\mathrm{ln}\,{I}_{j}^{-1}-1),\end{eqnarray} \tag{ 11 }$$

$$\begin{eqnarray}&&{\bar{\nu }}_{{\rm{s}}1}=\displaystyle \frac{1}{2}\displaystyle \frac{1}{\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}}\left(1+\displaystyle \sum _{j}3\,\displaystyle \frac{{n}_{j}}{{n}_{{\rm{e}}}}{N}_{{\rm{e}},j}\right).\,\end{eqnarray} \tag{ 12 }$$

2.2. Radiation losses

At the high densities typical of post-disruption scenarios, bremsstrahlung may be an important energy loss mechanism compared to synchrotron radiation reaction [25, 26]. In a fully ionized plasma, the required density for bremsstrahlung dominance is [27]

$$\begin{eqnarray}&&{n}_{{\rm{e}},20}\gtrsim {B}_{{\rm{T}}}^{2},\end{eqnarray} \tag{ 13 }$$

with B_T in units of Tesla and n_e,20 normalized to 10²⁰ m⁻³. In a partially ionized plasma, both bremsstrahlung and synchrotron losses will be enhanced, the latter through the increased pitch-angle scattering. Both radiative energy loss channels can therefore be significant at densities characteristic of disruptions and are included in this paper.

The synchrotron radiation reaction force is given by [28, 29]

$$\begin{eqnarray}\displaystyle \frac{\partial }{\partial {\boldsymbol{p}}}\cdot ({{\boldsymbol{F}}}_{\mathrm{syn}}f) & = & -\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{\partial }{\partial p}\left(\displaystyle \frac{{p}^{3}\gamma }{{\tau }_{\mathrm{syn}}}(1-{\xi }^{2})f\right)\\ & & +\displaystyle \frac{\partial }{\partial \xi }\left(\displaystyle \frac{\xi (1-{\xi }^{2})}{{\tau }_{\mathrm{syn}}\gamma }f\right),\end{eqnarray} \tag{ 14 }$$

where ${\tau }_{\mathrm{syn}}$ is the synchrotron radiation-damping time scale normalized to τ_c:

$$\begin{eqnarray}&&{\tau }_{\mathrm{syn}}^{-1}=\displaystyle \frac{{\tau }_{{\rm{c}}}{e}^{4}{B}^{2}}{6\pi {\epsilon }_{0}{m}_{{\rm{e}}}^{3}{c}^{3}}\approx \displaystyle \frac{1}{15.44\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}}\displaystyle \frac{{B}_{{\rm{T}}}^{2}}{{n}_{{\rm{e}},20}}.\end{eqnarray} \tag{ 15 }$$

We model partially screened bremsstrahlung with a Boltzmann operator as presented in [26], using the model that neglects the angular deflection due to the bremsstrahlung process:

$$\begin{eqnarray*}{C}_{\mathrm{br}}(p,\xi ) & = & \displaystyle \int {v}_{1}f({p}_{1},\xi )\displaystyle \frac{\partial {\sigma }^{\mathrm{br}}(p,{p}_{1})}{\partial p}{\rm{d}}{p}_{1}\\ & & -{vf}(p,\xi ){\sigma }^{\mathrm{br}}(p),\end{eqnarray*}$$

where $\partial {\sigma }^{\mathrm{br}}(p,{p}_{1})/\partial p$ is the normalized cross-section for an incident electron with momentum p₁ to end up with momentum p after emitting a bremsstrahlung photon carrying the energy difference, and σ^br is the total bremsstrahlung cross-section for an incident electron of momentum p. The integration is taken over $\sqrt{{(\gamma +{k}_{{\rm{c}}})}^{2}-1}\leqslant {p}_{1}\lt \infty$, where, following [26], photon energies are cut off at 0.1% of the kinetic energy of the outgoing electrons in order to resolve the infrared divergence, i.e. k_c = (γ − 1)/1000. The partially screened bremsstrahlung cross-section is given in [30, 31]:

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\sigma }^{\mathrm{br}}}{\partial p}(p,{p}_{1})=\displaystyle \frac{\alpha }{\pi \mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}}\displaystyle \frac{1}{k}\displaystyle \sum _{j}\displaystyle \frac{{n}_{j}}{{n}_{{\rm{e}}}}\left[\left(1+\displaystyle \frac{{\gamma }^{2}}{{\gamma }_{1}^{2}}\right)\right.\\ &&\quad \times \,\left({Z}_{j}^{2}+{\displaystyle \int }_{{q}_{0}}^{1}{[{Z}_{j}-{F}_{j}(q)]}^{2}\displaystyle \frac{{(q-{q}_{0})}^{2}}{{q}^{3}}{\rm{d}}q\right)\\ &&\quad -\,\displaystyle \frac{2}{3}\displaystyle \frac{\gamma }{{\gamma }_{1}}\left(\displaystyle \frac{5}{6}{Z}_{j}^{2}+{\displaystyle \int }_{{q}_{0}}^{1}{[{Z}_{j}-{F}_{j}(q)]}^{2}\right.\\ &&\quad \left.\left.\times \,\displaystyle \frac{{q}^{3}\,+\,3{{qq}}_{0}^{2}[1-2{{qq}}_{0}^{2}\mathrm{ln}\,(q/{q}_{0})]-4{q}_{0}^{3}}{{q}^{4}}{\rm{d}}q\right)\right],\end{eqnarray} \tag{ 16 }$$

where k is the photon momentum and q₀ = p₁ − p − k. We use the form factor F(q) for partially ionized atoms presented in [9],

$$\begin{eqnarray*}&&{F}_{j}(q)=\displaystyle \frac{{N}_{{\rm{e}},j}}{1+{(q{\bar{a}}_{j})}^{3/2}}.\end{eqnarray*}$$

In order to get an analytically tractable problem when deriving the effective critical electric field, a simplified bremsstrahlung mean-force stopping power will be used in section 3. Although a mean-force model has been shown to significantly alter the steady-state electron distribution compared to the full Boltzmann model, it captures the mean energy accurately [26], and is therefore sufficient for the purpose of deriving the effective critical electric field. This assumption is verified with numerical calculations using the full Boltzmann operator in section 4.

For the mean-force model, we have

$$\begin{eqnarray}&&{C}_{\mathrm{br}}\{f\}\approx -\displaystyle \frac{\partial }{\partial {\boldsymbol{p}}}\cdot ({{\boldsymbol{F}}}_{\mathrm{br}}f)=\displaystyle \frac{1}{{p}^{2}}\displaystyle \frac{\partial }{\partial p}({p}^{2}{F}_{\mathrm{br}}f),\end{eqnarray} \tag{ 17 }$$

where the bremsstrahlung mean-force is given by ${F}_{\mathrm{br}}(p)\,=\int \,k(\partial {\sigma }^{\mathrm{br}}({p}_{1},p)/\partial {p}_{1})\,{\rm{d}}{p}_{1}$, the integral taken over all allowed outgoing momenta p₁. For argon and neon, a numerical investigation of (16) shows that ${F}_{\mathrm{br}}$ is well approximated by

$$\begin{eqnarray}{F}_{\mathrm{br}} & \approx & p({\phi }_{\mathrm{br}0}+{\phi }_{\mathrm{br}1}\mathrm{ln}p)\\ & \equiv & \displaystyle \frac{p\alpha }{\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}}\displaystyle \sum _{j}\displaystyle \frac{{n}_{j}}{{n}_{{\rm{e}}}}{Z}_{j}^{2}(0.35+0.20\mathrm{ln}\,p).\end{eqnarray} \tag{ 18 }$$

The critical electric field is a central parameter for both generation of a runaway current and for its decay rate in a highly inductive tokamak; in the latter case, it is predicted that once the Ohmic current has dissipated, the induced electric field will be close to the critical electric field so that the current decays according to ${\rm{d}}I/{\rm{d}}t=2\pi {{RE}}_{{\rm{c}}}^{\mathrm{eff}}/L$ [21], where L ∼ μ₀R is the self-inductance and R is the major radius of the tokamak. The physical argument is that the runaway avalanche time scale is much faster than the inductive time scale, and therefore the electric field must be close to the critical electric field to prevent rapid current variations.

We calculate the effective electrical field due to collisions with partially screened ions by finding the minimum electric field ${E}_{{\rm{c}}}^{\mathrm{eff}}$ that satisfies the pitch-angle averaged force balance equation

$$\begin{eqnarray*}&&\langle {eE}\xi -F\rangle =0,\end{eqnarray*}$$

where F denotes the collisional and radiation forces on a runaway electron.

In order to find ${E}_{{\rm{c}}}^{\mathrm{eff}}$, we assume rapid pitch-angle dynamics compared to the time scale of the energy dynamics [20, 32]. In the kinetic equation (1), this amounts to requiring that the pitch-angle flux vanishes. Since ${\tau }_{\mathrm{syn}}^{-1}\ll 1$ from (15), we can neglect the effect of radiation on the pitch-angle distribution (term marked as 'neglect' below) as well as the effect of the avalanche source, which is slower than both pitch-angle scattering and collisional friction. We demonstrate the validity of these assumptions in appendix B by comparing the resulting critical electric field and angular distribution to kinetic simulations. Inserting the collision frequencies (6) and (10) as well as the radiation terms (14) and (17), the kinetic equation (1) can be rewritten

$$\begin{eqnarray}\displaystyle \frac{\partial \bar{f}}{\partial \tau } & = & \displaystyle \frac{\partial }{\partial p}\left[\left(-\displaystyle \frac{\xi E}{{E}_{{\rm{c}}}}+p{\nu }_{{\rm{s}}}+{F}_{\mathrm{br}}+\displaystyle \frac{p\gamma }{{\tau }_{\mathrm{syn}}}(1-{\xi }^{2})\right)\bar{f}\right]\\ & & +\displaystyle \frac{\partial }{\partial \xi }[(1-{\xi }^{2})\mathop{\underbrace{\left(-\displaystyle \frac{1}{p}\displaystyle \frac{E}{{E}_{{\rm{c}}}}\bar{f}+\displaystyle \frac{1}{2}{\nu }_{{\rm{D}}}\displaystyle \frac{\partial \bar{f}}{\partial \xi }\right)}}\limits_{=0}\\ & & -\mathop{\underbrace{\displaystyle \frac{\xi (1-{\xi }^{2})}{{\tau }_{\mathrm{syn}}\gamma }}}\limits_{\mathrm{neglect}}\bar{f}]\,,\end{eqnarray} \tag{ 19 }$$

where $\bar{f}={p}^{2}f$.

Following the method and notation of [20], the condition that the pitch-angle flux vanishes yields the following form for the angular distribution:

$$\begin{eqnarray}&&\bar{f}=G(t,p)A\exp (A\xi )/2\sinh A,\end{eqnarray} \tag{ 20 }$$

where the parameter A is defined as

$$\begin{eqnarray*}&&A(p)\equiv \displaystyle \frac{2E}{p{\nu }_{{\rm{D}}}{E}_{{\rm{c}}}}.\end{eqnarray*}$$

Then, (19) integrated over pitch-angle yields a continuity equation

$$\begin{eqnarray*}&&\displaystyle \frac{\partial G}{\partial \tau }+\displaystyle \frac{\partial }{\partial p}[U(p)G]=0,\end{eqnarray*}$$

where

$$\begin{eqnarray}U(p) & = & \displaystyle \frac{E/{E}_{{\rm{c}}}}{\tanh A}-\left[p{\nu }_{{\rm{s}}}+{F}_{\mathrm{br}}+\displaystyle \frac{p{\nu }_{{\rm{D}}}}{2}\right.\\ & & +\left.\displaystyle \frac{{p}^{2}\gamma {\nu }_{{\rm{D}}}}{{\tau }_{\mathrm{syn}}E/{E}_{{\rm{c}}}}\left(\displaystyle \frac{1}{\tanh A}-\displaystyle \frac{1}{A}\right)\right].\end{eqnarray} \tag{ 21 }$$

As the sign of U(p) determines if the distribution at p is accelerated or decelerated, the effective critical electric field is the minimum electric field for which force balance is possible:

$$\begin{eqnarray}&&{E}_{{\rm{c}}}^{\mathrm{eff}}\equiv \mathop{\min }\limits_{p}[E| U(p,E)=0].\end{eqnarray} \tag{ 22 }$$

The minimum can be found analytically if A ≫ 1 (so that $\tanh A\approx 1$) and the critical momentum fulfills ${p}_{{\rm{c}}}({E}_{{\rm{c}}}^{\mathrm{eff}})\gg 1$, which are consistent with our final solution if partially ionized impurities dominate. Hence (6), (10) and (18) may be used, and (22) is approximately solved by (see appendix C for more details):³

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{{\rm{c}}}^{\mathrm{eff}}}{{E}_{{\rm{c}}}}\approx {\bar{\nu }}_{{\rm{s}}0}+{\bar{\nu }}_{{\rm{s}}1}\left[\left(1+\displaystyle \frac{{\bar{\nu }}_{{\rm{D}}1}}{{\bar{\nu }}_{{\rm{D}}0}}\right)\mathrm{ln}\displaystyle \frac{{\bar{\nu }}_{{\rm{D}}0}}{2{\bar{\nu }}_{{\rm{s}}1}}+\sqrt{2\delta +1}\right],\end{eqnarray} \tag{ 23 }$$

where the constants are given in (2), (7), (8), (11), (12), (15), and (18), and δ, which is a measure of the effect of radiation losses, is given by

$$\begin{eqnarray}&&\delta ({E}_{{\rm{c}}}^{\mathrm{eff}})=\displaystyle \frac{{\bar{\nu }}_{{\rm{D}}0}}{{\bar{\nu }}_{{\rm{s}}1}^{2}}\left(\displaystyle \frac{{\bar{\nu }}_{{\rm{D}}0}{\tau }_{\mathrm{syn}}^{-1}}{{E}_{{\rm{c}}}^{\mathrm{eff}}/{E}_{{\rm{c}}}}+{\phi }_{\mathrm{br}0}+{\phi }_{\mathrm{br}1}\mathrm{ln}({\bar{\nu }}_{{\rm{D}}0}/2{\bar{\nu }}_{{\rm{s}}1})\right).\end{eqnarray} \tag{ 24 }$$

Since δ depends on ${E}_{{\rm{c}}}^{\mathrm{eff}}$, (23) is not in a closed form, and therefore (23) and (24) are evaluated iteratively starting at ${E}_{{\rm{c}}}^{\mathrm{eff}}={E}_{{\rm{c}}}^{\mathrm{tot}}$, where ${E}_{{\rm{c}}}^{\mathrm{tot}}$ is the critical electric field including the density of both bound and free electrons:

$$\begin{eqnarray}&&{E}_{{\rm{c}}}^{\mathrm{tot}}\equiv \displaystyle \frac{{n}_{{\rm{e}}}^{\mathrm{tot}}}{{n}_{{\rm{e}}}}{E}_{{\rm{c}}}=\displaystyle \frac{{n}_{{\rm{e}}}^{\mathrm{tot}}{e}^{3}\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}}{4\pi {\epsilon }_{0}^{2}{m}_{{\rm{e}}}{c}^{2}},\end{eqnarray} \tag{ 25 }$$

with ${n}_{{\rm{e}}}^{\mathrm{tot}}={n}_{{\rm{e}}}+{\sum }_{j}{n}_{j}{N}_{{\rm{e}},j}$. Here, we iterate once so that ${\delta }_{0}\,=\delta ({E}_{{\rm{c}}}^{\mathrm{eff}}={E}_{{\rm{c}}}^{\mathrm{tot}})$ and $\delta \approx {\delta }_{1}=\delta [{E}_{{\rm{c}}}^{\mathrm{eff}}({\delta }_{0})]$. Equation (23) was found to be accurate to within 10% for magnetic fields in the range ${B}_{{\rm{T}}}^{2}\lesssim 100\,{n}_{20}^{\mathrm{tot}}$ for all considered impurity species and plasma compositions.

Figure 2 shows the effective critical electric field normalized to ${E}_{{\rm{c}}}^{\mathrm{tot}}$. Our model, corresponding to (23), is shown in black and compared to the full numerical solution to (22) (using the algorithm in [33], implemented as fmincon in matlab) for three different values of the magnetic field: B = 0 T in solid line, B = 2 T dashed and B = 5 T in dotted line. These are shown for singly ionized argon in figure 2(a) and singly ionized neon in 2(b). The behavior is only weakly dependent on ionization states; this is illustrated with neutral argon and Ar⁴⁺ in figure 3. In addition, we find that the background deuterium density has a negligible effect on ${E}_{{\rm{c}}}^{\mathrm{eff}}$ when Zn_Z ≫ n_D.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 2–3 also show that with weakly ionized impurities,

$$\begin{eqnarray*}&&{E}_{{\rm{c}}}^{\mathrm{eff}}\gtrsim {E}_{{\rm{c}}}^{\mathrm{tot}}\gg {E}_{{\rm{c}}}.\end{eqnarray*}$$

Hence, it is more accurate to include all electrons in the critical electric field, than to count for instance half of the bound electrons as done in the RP model (${E}_{{\rm{c}}}^{\mathrm{RP}}\,={E}_{{\rm{c}}}^{\mathrm{tot}}({n}_{{\rm{e}}}\,+0.5{n}_{\mathrm{bound}})/{n}_{{\rm{e}}}^{\mathrm{tot}}$). This underestimation of the effective critical field by the RP model is a result of using a simplistic form of the inelastic collision rate as well as neglecting the effect of pitch-angle scattering and radiation losses. To further explore the scaling of ${E}_{{\rm{c}}}^{\mathrm{eff}}$ with magnetic field strength and impurity content, we approximate (23) in the case where one weakly ionized state j dominates:

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{{\rm{c}}}^{\mathrm{eff}}}{{E}_{{\rm{c}}}^{\mathrm{tot}}}\approx \displaystyle \frac{{n}_{{\rm{e}}}}{{n}_{{\rm{e}}}^{\mathrm{tot}}}+\displaystyle \frac{{N}_{{\rm{e}},j}{n}_{Z}}{{n}_{{\rm{e}}}^{\mathrm{tot}}}\displaystyle \frac{1}{\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}}\left({S}_{j}+{R}_{j}\sqrt{\displaystyle \frac{{B}_{{\rm{T}}}^{2}}{{n}_{20}^{\mathrm{tot}\ }}+\displaystyle \frac{0.9}{\mathrm{ln}\,{\bar{a}}_{j}}}\right),\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{S}_{j}=\left[\mathrm{ln}\,{I}_{j}^{-1}-1+\displaystyle \frac{3}{2}\left(1+\displaystyle \frac{1}{\mathrm{ln}\,{\bar{a}}_{j}}\right)\mathrm{ln}\,\left(\displaystyle \frac{{Z}_{j}}{3}\mathrm{ln}\,{\bar{a}}_{j}\right)\right],\end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray}&&{R}_{j}=0.09(Z+{Z}_{0})\mathrm{ln}\,{\bar{a}}_{j}.\,\end{eqnarray} \tag{ 28 }$$

The screening constant S_j is given for all argon and neon species in table A1 in appendix A. For typical magnetic fields, the terms inside the brackets tend to be roughly 1–2 times $\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}$. As ${n}_{{\rm{e}}}+{N}_{{\rm{e}},j}{n}_{j}={n}_{{\rm{e}}}^{\mathrm{tot}}$ with only one impurity species j, one obtains ${E}_{{\rm{c}}}^{\mathrm{eff}}\gtrsim {E}_{{\rm{c}}}^{\mathrm{tot}}$. From (26), we thus conclude that the effect of partially stripped impurities scale approximately linearly with impurity density; more specifically, ${E}_{{\rm{c}}}^{\mathrm{eff}}\,={E}_{{\rm{c}}}^{\mathrm{tot}}({n}_{{\rm{e}}}+\kappa {n}_{\mathrm{bound}})/{n}_{{\rm{e}}}^{\mathrm{tot}}\approx \kappa {E}_{{\rm{c}}}^{\mathrm{tot}}$, where κ is between 1 and 2. Consequently, our calculated of ${E}_{{\rm{c}}}^{\mathrm{eff}}$ is up to 4E_c^RP in typical tokamak scenarios.

The radiation term R_j quantifies the effect of bremsstrahlung and synchrotron losses; these are dominated by synchrotron radiation reaction if

$$\begin{eqnarray*}&&{B}_{{\rm{T}}}^{2}\gtrsim 0.2{n}_{20}^{\mathrm{tot}},\end{eqnarray*}$$

which is lower than the fully ionized estimation (13). In this case, ${E}_{{\rm{c}}}^{\mathrm{eff}}$ depends linearly on ${B}_{{\rm{T}}}/\sqrt{{n}_{20}^{\mathrm{tot}\ }}$. This agrees with the scaling found in [20] for the fully ionized case. In contrast, for low magnetic fields, bremsstrahlung can increase the effective critical field by up to 20% for argon. This number is insensitive to the plasma density and depends only on its ionic composition.

The critical electric field, especially as modified by the effects of partially screened nuclei and radiation losses, plays an important role during the relaxation of runaway electrons. In this section, we demonstrate with kinetic simulations that (23) well characterizes the threshold between runaway growth and decay under these modifications. Then, when the electric field evolves self-consistently, we show that it remains tied to ${E}_{{\rm{c}}}^{\mathrm{eff}}$ under certain assumptions during the current decay phase of a tokamak disruption.

If the current is carried by runaway electrons and the shape of the runaway distribution is constant in time, the time derivative of the current is related to the steady-state runaway growth rate

$$\begin{eqnarray}&&{\rm{\Gamma }}(E)\equiv \displaystyle \frac{1}{{n}_{\mathrm{RE}}}\displaystyle \frac{{\rm{d}}{n}_{\mathrm{RE}}}{{\rm{d}}t}\approx \displaystyle \frac{1}{I}\displaystyle \frac{{\rm{d}}I}{{\rm{d}}t}.\end{eqnarray} \tag{ 29 }$$

The scaling of the growth rate with impurity content may be estimated from the RP formula [22] by replacing E_c with ${E}_{{\rm{c}}}^{\mathrm{eff}}$ and the density by the total electron density due to the fact that bound and free electrons have equal probability of becoming runaway electrons through knock-on collisions:

$$\begin{eqnarray}&&{\rm{\Gamma }}(E)\sim \displaystyle \frac{1}{\mathrm{ln}\,{{\rm{\Lambda }}}_{{\rm{c}}}}\displaystyle \frac{1}{{\tau }_{{\rm{c}}}^{\mathrm{tot}}}\left(\displaystyle \frac{E}{{E}_{{\rm{c}}}^{\mathrm{eff}}}-1\right),\end{eqnarray} \tag{ 30 }$$

with ${\tau }_{{\rm{c}}}^{\mathrm{tot}}=({n}_{{\rm{e}}}/{n}_{{\rm{e}}}^{\mathrm{tot}}){\tau }_{{\rm{c}}}$. The qualitative scaling of the analytic growth rate is confirmed in figure 4, where the growth rate is numerically calculated using code [34, 35], which directly solves the kinetic equation (1). These simulations employed the general field-particle knock-on operator of [36–38] and a Boltzmann operator for partially screened bremsstrahlung losses as described in section 2.2. The vertical lines denote the analytic prediction in (23) for when one would expect the transition between growth and decay of an existing runaway population. Radiation losses affect where this threshold lies and the analytic model ${E}_{{\rm{c}}}^{\mathrm{eff}}$ accurately and robustly captures this effect. In particular, we note that the mean-force bremsstrahlung model employed in the analytical derivation of ${E}_{{\rm{c}}}^{\mathrm{eff}}$ agrees with the Boltzmann-type bremsstrahlung operator used in the simulations within a few percent.

Zoom In Zoom Out Reset image size

The electric field is hypothesized to remain close to ${E}_{{\rm{c}}}^{\mathrm{eff}}$ during the current decay phase of a tokamak disruption [21]. The mechanism by which this occurs is the fast time scale of the avalanche generation in relation to the inductive time scale of the system. A toroidal electric field is induced when there is a time-changing magnetic flux through a current loop such as a runaway beam. This magnetic flux is proportional to the total current through the loop. The induced electric field is therefore related to the rate of change of the current:

$$\begin{eqnarray}&&E=-\displaystyle \frac{L}{2\pi R}\displaystyle \frac{{\rm{d}}I}{{\rm{d}}t},\end{eqnarray} \tag{ 31 }$$

where R is the major radius of the tokamak. This inductance model has recently been implemented in code to calculate the electric field self-consistently with the evolution of the electron velocity distribution. In general, the exact value of the inductance L will depend on the spatial distribution of current, which will change in time. For a large-aspect ratio current loop (such as a runaway beam), L can be approximated by [39]

$$\begin{eqnarray}&&L\approx {\mu }_{0}R\left[\mathrm{ln}\left(\displaystyle \frac{8R}{a}\right)-2+\displaystyle \frac{{l}_{i}}{2}\right].\end{eqnarray} \tag{ 32 }$$

Here, R is the major radius of the tokamak, a is the radius of the runaway beam, and l_i parametrizes the distribution of current within the beam. We have chosen l_i = 1.5 as a representative mid-plateau value, based on experimental results from European medium sized tokamaks.

When $E\approx {E}_{{\rm{c}}}^{\mathrm{eff}}$, the growth rate can be expanded according to

$$\begin{eqnarray*}&&{\rm{\Gamma }}={\rm{\Gamma }}^{\prime} ({E}_{{\rm{c}}}^{\mathrm{eff}})[E-{E}_{{\rm{c}}}^{\mathrm{eff}}]+...,\end{eqnarray*}$$

which allows (31) to be solved analytically:

$$\begin{eqnarray}&&E\approx {E}_{{\rm{c}}}^{\mathrm{eff}}\left(1-\displaystyle \frac{2\pi R}{{{LI}}_{\mathrm{RE}}{\rm{\Gamma }}^{\prime} ({E}_{{\rm{c}}}^{\mathrm{eff}})}\right).\end{eqnarray} \tag{ 33 }$$

This yields a condition under which the electric field remains close to ${E}_{{\rm{c}}}^{\mathrm{eff}}$:

$$\begin{eqnarray*}&&{{LI}}_{\mathrm{RE}}{\rm{\Gamma }}^{\prime} ({E}_{{\rm{c}}}^{\mathrm{eff}})\gg 2\pi R.\end{eqnarray*}$$

With the estimate of ${\rm{\Gamma }}^{\prime} (E)$ from the numerical results of figure 4 (at B = 0 T) and estimating R/a ≈ 5 we find that the minimum required current for $E\approx {E}_{{\rm{c}}}^{\mathrm{eff}}$ is approximately

$$\begin{eqnarray}&&{I}_{\mathrm{RE}}\,\gg 60\,\mathrm{kA}.\end{eqnarray} \tag{ 34 }$$

This value is substantially lower than the estimation of 250 kA in [21], which did not include the effect of partial screening or radiation losses. Since this threshold current is inversely proportional to the inductance, the estimate (34) is only weakly dependent on the details of the spatial current distribution. Therefore, the exact value of the instantaneous inductance does not affect the primary result of this section: for large enough inductance, the electric field remains approximately tied to ${E}_{{\rm{c}}}^{\mathrm{eff}}$ during the current decay phase, leading to a predictable decay time scale.

To test the hypothesis that $E\approx {E}_{{\rm{c}}}^{\mathrm{eff}}$ when I_RE ≫ 60 kA, we generate a forward-beamed initial distribution obtained from a simulation with a large electric field; the initial average runaway energy in our simulation is 17.2 MeV. We then inject singly ionized argon with a density that is four times the deuterium density ${n}_{{\rm{D}}}={10}^{20}\,{{\rm{m}}}^{-3}$. Starting at an initial current density ${j}_{0}=12.9\,\mathrm{MA}\,{{\rm{m}}}^{-2}$, we let the electron distribution evolve with a self-consistent electric field in a strongly, intermediate or weakly inductive system. At a constant current density, varying ${I}_{0}^{\mathrm{RE}}L/({\mu }_{0}R)$ corresponds to varying L/(μ₀R) through the beam aspect ratio R/a or the initial current ${I}_{0}\,={j}_{0}\pi {a}^{2}$. The following values were chosen in the simulations: $\pi {a}^{2}L/({\mu }_{0}R)\,=({\rm{i}})\,4.30,({\rm{ii}})\,1.57\,\mathrm{and}\,({\rm{iii}})\,0.14$. If R/a = 5 and l_i = 1.5, these three values correspond to an initial current of $({\rm{i}})\,{I}_{0}^{\mathrm{RE}}\,=\,23\,{\rm{MA}};$ $({\rm{ii}})\,{I}_{0}^{\mathrm{RE}}\,=\,8.3\,{\rm{MA}};$ and $({\rm{iii}})\,{I}_{0}^{\mathrm{RE}}\,=\,0.75{\rm{MA}}$. As in the growth rate simulations, we include both synchrotron losses, the full bremsstrahlung model and a Chiu-Harvey type avalanche operator.

Figure 5(a) shows the current decay, which is linear (as expected) and faster in the low-inductance case. Figure 5(b) shows the electric field evolution. Clearly, in the high-inductance case, the electric field is close to the critical field after an initial transient. This means that, in highly inductive devices such as ITER, the current decay is to a very good approximation given by ${\rm{d}}{I}_{\mathrm{RE}}/{\rm{d}}t=-2\pi {{RE}}_{{\rm{c}}}^{\mathrm{eff}}/L$. Enhanced E_c^eff will lead to faster current decay, and (23) quantifies how fast the decay is.

Zoom In Zoom Out Reset image size

On the other hand, the induced electric field deviates by approximately 10% from ${E}_{{\rm{c}}}^{\mathrm{eff}}$ in the low-inductance case. Since the initial current I₀^RE = 750 kA is high in relation to many medium sized tokamak experiments, $E\approx {E}_{{\rm{c}}}^{\mathrm{eff}}$ gives an overestimation of the current decay rate in many of today's devices. The relative deviation from ${E}_{{\rm{c}}}^{\mathrm{eff}}$ observed in figure 5(b) is consistent with the estimation $1-E/{E}_{{\rm{c}}}^{\mathrm{eff}}\approx 60\,\mathrm{kA}/{I}_{\mathrm{RE}}$ from (33) and (34).

Although the predicted induced electric field obeys $E\leqslant {E}_{{\rm{c}}}^{\mathrm{eff}}$ with our assumptions, several effects could lead to a higher induced electric field in an actual experimental discharge. For example, a stronger electric field would be necessary to balance a runaway population with sub-relativistic energy, in which case the steady-state growth rate used here is inaccurate. Other effects such as transport [40–42], trapping [22, 43] and wave–particle interaction [10, 44–46] may also increase the runaway current decay rate and accordingly the induced electric field. Such complete modeling remains the subject of future work. Nevertheless, partial screening has a major effect on the critical electric field as demonstrated here, and therefore the results derived herein should be an important piece toward improved experimental comparison of the runaway current decay rate as well as the avalanche growth rate.

Finally, we note that the simulations with an inductive electric field validate the initial assumption of rapid pitch-angle dynamics in (19); we find that the resulting pitch-angle distribution in (20) is accurate for $E\approx {E}_{{\rm{c}}}^{\mathrm{eff}};$ see appendix B. The distribution function in (20) is consequently appropriate for determining the effective critical electric field, but not for describing runaway generation.

Recent experimental studies on several tokamaks show that the onset and decay of runaway electrons occurs for critical electric fields that are considerably higher than the Connor–Hastie field E_c. One reason is that there are other runaway loss mechanisms in addition to damping due to collisions in a fully ionized plasma that seem to dominate both in disruptive and quiescent cases. In this paper, we show that if there are heavy partially ionized impurities present in the plasma, the dominant effect on the critical electric field is the effect of partial screening. The effective critical field is further increased due to the enhanced radiation loss rates when partially ionized impurities are present.

We give analytical formulas for the effective critical electric field ${E}_{{\rm{c}}}^{\mathrm{eff}}$ including partial screening and radiation effects, derived under the condition of rapid pitch-angle dynamics. The validity of this assumption and the value of the effective critical electric field is demonstrated by numerical simulations with the kinetic equation solver code. The most complete expression for the critical electric field is given in (23). It has been shown to be valid for a wide range of magnetic fields, impurity species and plasma composition. To make the parametric dependencies more transparent, we also give an approximate expression in (26) that is valid when one weakly ionized state dominates, which is often the case in a cold post-disruption tokamak plasma.

As expected, we find that in the presence of large amounts of heavy impurities, the effective critical field can be drastically higher than ${E}_{{\rm{c}}}$ which is proportional to the density of free electrons: ${E}_{{\rm{c}}}^{\mathrm{eff}}$ even exceeds the value obtained by including the total density of both free and bound electrons. In contrast to RP [22], where the effective density includes half of the bound electrons, $n={n}_{{\rm{e}}}+0.5{n}_{\mathrm{bound}}$, our calculations show that the bound electrons are weighted by a factor of typically 1–2. This enhancement is attributed to the energy-dependent collisional friction, pitch-angle scattering as well as radiation losses. Bremsstahlung and synchrotron losses both increase the effective critical field, typically by tens of percent.

Using a 0 D inductive electric field we calculate the runaway current decay after impurity injection. Through kinetic simulations we confirm the accuracy of the formula for the effective critical field (23), and demonstrate that the electric field stays close to the effective critical field when the runaway current satisfies I_RE ≫ 60 kA, in which case ${\rm{d}}{I}_{\mathrm{RE}}/{\rm{d}}t\propto {E}_{{\rm{c}}}^{\mathrm{eff}}$. These findings are relevant for the efficacy of mitigation strategies for runaway electrons in tokamak devices: since the runaway current decay rate is typically 2–4 times higher than what is predicted by the RP formula, a lower quantity of assimilated material is required for successful mitigation. As screening significantly increases the critical electric field, we anticipate that this effect is of importance to include in experimental comparisons; however, accurate predictions may require the modeling of spatial effects which are not considered here.

The authors are grateful to E Hollmann, S Newton and A Stahl for stimulating discussions and to T DuBois and M Rahm for the simulations needed to determine the effective ion size. This work was supported by the Swedish Research Council (Dnr. 2014-5510), the Knut and Alice Wallenberg Foundation and the European Research Council (ERC-2014-CoG grant 647121). The work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Table A1 summarizes the constants needed to compute the value of the effective electric field in the presence of argon and neon. The effective ion size ${\bar{a}}_{j}$ is determined by DFT simulations and is related to a_j in [9] through ${\bar{a}}_{j}=2{a}_{j}/\alpha$ where α ≈ 1/137 is the fine-structure constant. The mean excitation energy I_j is taken from [12]. These give S_j from (27) according to

$$\begin{eqnarray*}&&{S}_{j}=\left[\mathrm{ln}{I}_{j}^{-1}-1+\displaystyle \frac{3}{2}\left(1+\displaystyle \frac{1}{\mathrm{ln}{\bar{a}}_{j}}\right)\mathrm{ln}\left(\displaystyle \frac{{Z}_{j}}{3}\mathrm{ln}{\bar{a}}_{j}\right)\right].\end{eqnarray*}$$

Table A1.  Constants to determine ${E}_{{\rm{c}}}^{\mathrm{eff}}$.

	$\mathrm{ln}{\bar{a}}_{j}$	$\mathrm{ln}{I}_{j}^{-1}$	Sj			$\mathrm{ln}{\bar{a}}_{j}$	$\mathrm{ln}{I}_{j}^{-1}$	Sj
Ar0	4.6	7.9	13.0		Ne0	4.7	8.2	12.2
Ar1+	4.5	7.8	12.8		Ne1+	4.6	8.0	12.0
Ar2+	4.4	7.6	12.6		Ne2+	4.5	7.9	11.8
Ar3+	4.4	7.5	12.5		Ne3+	4.4	7.7	11.6
Ar4+	4.3	7.3	12.3		Ne4+	4.3	7.5	11.4
Ar5+	4.2	7.2	12.2		Ne5+	4.1	7.3	11.2
Ar6+	4.1	7.0	12.0		Ne6+	4.0	7.0	10.8
Ar7+	4.0	6.8	11.8		Ne7+	3.7	6.6	10.4
Ar8+	3.9	6.6	11.5		Ne8+	3.2	5.9	9.5
Ar9+	3.8	6.5	11.4		Ne9+	3.1	5.8	9.5
Ar10+	3.7	6.4	11.3
Ar11+	3.6	6.2	11.1
Ar12+	3.6	6.1	11.0
Ar13+	3.5	5.9	10.8
Ar14+	3.3	5.7	10.5
Ar15+	3.1	5.3	10.1
Ar16+	2.6	4.7	9.4
Ar17+	2.5	4.7	9.4

The simulations with an inductive electric field (figure 5) can be used to validate the initial assumption of rapid pitch-angle dynamics in (19) leading to the pitch-angle distribution in (20). Expanding $\bar{f}$ in Legendre polynomials

$$\begin{eqnarray*}&&\bar{f}=\displaystyle \sum _{L}{\bar{f}}_{L}(p){P}_{L}(\xi ),\end{eqnarray*}$$

we relate the predicted analytical distribution in (20) to the ratio between the zeroth and the first Legendre modes of the distribution:

$$\begin{eqnarray}&&\displaystyle \frac{1}{3}\displaystyle \frac{{\bar{f}}_{1}}{{\bar{f}}_{0}}=\left(\displaystyle \frac{1}{\tanh A}-\displaystyle \frac{1}{A}\right).\end{eqnarray} \tag{ B.1 }$$

The ratio given in (B.1) quantifies the narrowness of the electron distribution: ${\bar{f}}_{1}/3{\bar{f}}_{0}=0$ corresponds to an isotropic distribution while the ${\bar{f}}_{1}/3{\bar{f}}_{0}\to 1$ for a narrow, beam-like distribution. Figure B1 compares the numerical value of ${\bar{f}}_{1}/3{\bar{f}}_{0}$ as computed in code in solid black line, to the analytical prediction (B.1) in dashed green line. The analytical formula accurately predicts the distribution width on the entire interval from a fully isotropic distribution at p = 0 to a narrow beam for p ≫ 1. This validates our assumptions on the rapid pitch-angle dynamics in (19). In contrast, for larger electric fields ($E/{E}_{{\rm{c}}}^{\mathrm{eff}}\gtrsim 5$), we find that the distribution rather follows the formula in Fülöp et al [47], which is derived in the limit of $E\gg {E}_{{\rm{c}}}^{\mathrm{eff}}$.

Zoom In Zoom Out Reset image size

The effective critical field can be found analytically noting that the critical momentum fulfills ${p}_{{\rm{c}}}^{\star }\equiv {p}_{{\rm{c}}}({E}_{{\rm{c}}}^{\mathrm{eff}})\gg 1$. Moreover, we assume that A, which is defined in (20), fulfills A ≫ 1 (so that $\tanh A\approx 1$). These two assumptions are consistent with our final solution if partially ionized impurities dominate. Hence (6), (10) and (18) may be used in the expression for the effective critical field (22), and the requirement U(p) = 0 [with U given in (21)] results in a quadratic equation in E/E_c:

$$\begin{eqnarray}&&{\left(\displaystyle \frac{E}{{E}_{{\rm{c}}}}\right)}^{2}-\displaystyle \frac{E}{{E}_{{\rm{c}}}}h(p)-\epsilon (p)=0,\end{eqnarray} \tag{ C.1 }$$

where h(p) and (p) are both positive functions of p within the assumption p_c ≫ 1:

$$\begin{eqnarray*}h(p) & \equiv & {\bar{\nu }}_{{\rm{s}}0}+{\bar{\nu }}_{{\rm{s}}1}\mathrm{ln}p+\displaystyle \frac{1}{2p}({\bar{\nu }}_{{\rm{D}}0}+{\bar{\nu }}_{{\rm{D}}1}\mathrm{ln}p)\\ & & +p({\phi }_{\mathrm{br}0}+{\phi }_{\mathrm{br}1}\mathrm{ln}p),\\ \epsilon (p) & \equiv & p({\bar{\nu }}_{{\rm{D}}0}+{\bar{\nu }}_{{\rm{D}}1}\mathrm{ln}p){\tau }_{\mathrm{syn}}^{-1}.\end{eqnarray*}$$

Consequently, finding the effective critical field amounts to evaluating the positive solution to (C.1)

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{{\rm{c}}}^{\mathrm{eff}}}{{E}_{{\rm{c}}}}\approx \displaystyle \frac{1}{2}\phantom{\rule{}{3.00ex}}[h({p}_{{\rm{c}}}^{\star })+\sqrt{h{({p}_{{\rm{c}}}^{\star })}^{2}+4\epsilon ({p}_{{\rm{c}}}^{\star })}\phantom{\rule{}{3.00ex}}]\end{eqnarray} \tag{ C.2 }$$

at the minimum ${p}_{{\rm{c}}}^{\star }$, the critical momentum which minimizes ${E}_{{\rm{c}}}^{\mathrm{eff}}$ in (C.1), which is determined by

$$\begin{eqnarray}&&\displaystyle \frac{{E}_{{\rm{c}}}^{\mathrm{eff}}}{{E}_{{\rm{c}}}}h^{\prime} ({p}_{{\rm{c}}}^{\star })+\epsilon ^{\prime} ({p}_{{\rm{c}}}^{\star })=0.\end{eqnarray} \tag{ C.3 }$$

The derivatives of h(p) and (p) are given by

$$\begin{eqnarray*}h^{\prime} ({p}_{{\rm{c}}}^{\star }) & \approx & \displaystyle \frac{{\bar{\nu }}_{{\rm{s}}1}}{{p}_{{\rm{c}}}^{\star }}-\displaystyle \frac{1}{2{({p}_{{\rm{c}}}^{\star })}^{2}}[{\bar{\nu }}_{{\rm{D}}0}+{\bar{\nu }}_{{\rm{D}}1}(\mathrm{ln}{p}_{{\rm{c}}}^{\star }-1)]\\ & & +{\phi }_{\mathrm{br}0}+{\phi }_{\mathrm{br}1}(\mathrm{ln}{p}_{{\rm{c}}}^{\star }+1),\\ \epsilon ^{\prime} ({p}_{{\rm{c}}}^{\star }) & \approx & {\tau }_{\mathrm{syn}}^{-1}[{\bar{\nu }}_{{\rm{D}}0}+{\bar{\nu }}_{{\rm{D}}1}(\mathrm{ln}{p}_{{\rm{c}}}^{\star }+1)],\end{eqnarray*}$$

and thus (C.3) is solved by

$$\begin{eqnarray*}&&{p}_{{\rm{c}}}^{\star }\approx \displaystyle \frac{2{p}_{{\rm{c}}0}^{\star }}{1+\sqrt{1+2\delta }},\end{eqnarray*}$$

where

$$\begin{eqnarray}\delta & \equiv & \displaystyle \frac{{\bar{\nu }}_{{\rm{D}}0}+{\bar{\nu }}_{{\rm{D}}1}(\mathrm{ln}{p}_{{\rm{c}}}^{\star }-1)}{{\bar{\nu }}_{{\rm{s}}1}^{2}}[{\phi }_{\mathrm{br}0}\\ & & +{\phi }_{\mathrm{br}1}(\mathrm{ln}{p}_{{\rm{c}}}^{\star }+1)]({x}_{\mathrm{rad}}+1),\end{eqnarray} \tag{ C.4 }$$

$$\begin{eqnarray}&&{p}_{{\rm{c}}0}^{\star }\equiv \displaystyle \frac{{\bar{\nu }}_{{\rm{D}}0}+{\bar{\nu }}_{{\rm{D}}1}(\mathrm{ln}{p}_{{\rm{c}}}^{\star }-1)}{2{\bar{\nu }}_{{\rm{s}}1}},\,\end{eqnarray} \tag{ C.5 }$$

$$\begin{eqnarray}&&{x}_{\mathrm{rad}}\equiv \displaystyle \frac{[{\bar{\nu }}_{{\rm{D}}0}+{\bar{\nu }}_{{\rm{D}}1}(\mathrm{ln}{p}_{{\rm{c}}}^{\star }+1)]{\tau }_{\mathrm{syn}}^{-1}}{[{\phi }_{\mathrm{br}0}+{\phi }_{\mathrm{br}1}(\mathrm{ln}{p}_{{\rm{c}}}^{\star }+1)]{E}_{{\rm{c}}}^{\mathrm{eff}}/{E}_{{\rm{c}}}}.\,\end{eqnarray} \tag{ C.6 }$$

Here, ${x}_{\mathrm{rad}}$ describes the relative importance of synchrotron radiation compared to bremsstrahlung.

To evaluate (C.2), we first simplify $h({p}_{{\rm{c}}}^{\star })$ using ${(1+\sqrt{1+2\delta })}^{-1}=(\sqrt{1+2\delta }-1)/2\delta$:

$$\begin{eqnarray}h({p}_{{\rm{c}}}^{\star }) & = & {\bar{\nu }}_{{\rm{s}}0}+{\bar{\nu }}_{{\rm{s}}1}\mathrm{ln}{p}_{{\rm{c}}}^{\star }\\ & & +\displaystyle \frac{{\bar{\nu }}_{{\rm{s}}1}}{2}\displaystyle \frac{({\bar{\nu }}_{{\rm{D}}0}+{\bar{\nu }}_{{\rm{D}}1}\mathrm{ln}{p}_{{\rm{c}}}^{\star })}{{\bar{\nu }}_{{\rm{D}}0}+{\bar{\nu }}_{{\rm{D}}1}(\mathrm{ln}{p}_{{\rm{c}}}^{\star }-1)}(1+\sqrt{2\delta +1})\\ & & +\displaystyle \frac{{\bar{\nu }}_{{\rm{s}}1}}{2}\displaystyle \frac{{\phi }_{\mathrm{br}0}+{\phi }_{\mathrm{br}1}\mathrm{ln}{p}_{{\rm{c}}}^{\star }}{{\phi }_{\mathrm{br}0}+{\phi }_{\mathrm{br}1}(\mathrm{ln}{p}_{{\rm{c}}}^{\star }+1)}\displaystyle \frac{\sqrt{1+2\delta }-1}{{x}_{\mathrm{rad}}+1}\\ & \approx & \mathop{\underbrace{{\bar{\nu }}_{{\rm{s}}0}+{\bar{\nu }}_{{\rm{s}}1}\left(\mathrm{ln}{p}_{{\rm{c}}}^{\star }+1+\displaystyle \frac{{\bar{\nu }}_{{\rm{D}}1}}{{\bar{\nu }}_{{\rm{D}}0}}\mathrm{ln}{p}_{{\rm{c}}}^{\star }\right)}}\limits_{\equiv {h}_{0}({p}_{{\rm{c}}}^{\star })}\\ & & +\displaystyle \frac{{\bar{\nu }}_{{\rm{s}}1}}{2}(\sqrt{1+2\delta }-1)\left(\displaystyle \frac{{x}_{\mathrm{rad}}+2}{{x}_{\mathrm{rad}}+1}\right),\end{eqnarray} \tag{ C.7 }$$

where we assumed ${\bar{\nu }}_{{\rm{D}}0}\gg {\bar{\nu }}_{{\rm{D}}1}(\mathrm{ln}{p}_{{\rm{c}}}^{\star }-1)$ since ${\bar{\nu }}_{{\rm{D}}0}\gg {\bar{\nu }}_{{\rm{D}}1}$ typically; see (7) and (8). Furthermore, we assumed ${\phi }_{\mathrm{br}1}\ll {\phi }_{\mathrm{br}0}+{\phi }_{\mathrm{br}1}\mathrm{ln}{p}_{{\rm{c}}}^{\star }$. To simplify $\epsilon ({p}_{{\rm{c}}}^{\star })$, we approximate ${E}_{{\rm{c}}}^{\mathrm{eff}}/{E}_{{\rm{c}}}\approx h({p}_{{\rm{c}}}^{\star })$ and assume ${\bar{\nu }}_{{\rm{D}}1}\ll {\bar{\nu }}_{{\rm{D}}0}+{\bar{\nu }}_{{\rm{D}}1}\mathrm{ln}{p}_{{\rm{c}}}^{\star }$:

$$\begin{eqnarray}\epsilon ({p}_{{\rm{c}}}^{\star }) & = & \displaystyle \frac{{\bar{\nu }}_{{\rm{s}}1}}{2}{E}_{{\rm{c}}}^{\mathrm{eff}}(\sqrt{1+2\delta }-1)\displaystyle \frac{{\bar{\nu }}_{{\rm{D}}0}+{\bar{\nu }}_{{\rm{D}}1}\mathrm{ln}{p}_{{\rm{c}}}^{\star }}{{\bar{\nu }}_{{\rm{D}}0}+{\bar{\nu }}_{{\rm{D}}1}(\mathrm{ln}{p}_{{\rm{c}}}^{\star }+1)}\\ & & \times \left(\displaystyle \frac{{x}_{\mathrm{rad}}}{{x}_{\mathrm{rad}}+1}\right)\\ & \approx & \displaystyle \frac{{\bar{\nu }}_{{\rm{s}}1}}{2}h({p}_{{\rm{c}}}^{\star })(\sqrt{1+2\delta }-1)\left(\displaystyle \frac{{x}_{\mathrm{rad}}}{{x}_{\mathrm{rad}}+1}\right).\end{eqnarray} \tag{ C.8 }$$

Then,

$$\begin{eqnarray}&&\begin{array}{l}\sqrt{h{({p}_{{\rm{c}}}^{\star })}^{2}+4\epsilon ({p}_{{\rm{c}}}^{\star })}\\ \quad \approx \,\sqrt{{h}_{0}({p}_{{\rm{c}}}^{\star })+\displaystyle \frac{{\bar{\nu }}_{{\rm{s}}1}}{2}(\sqrt{1+2\delta }-1)\displaystyle \frac{{x}_{\mathrm{rad}}+2}{{x}_{\mathrm{rad}}+1}}\\ \qquad \times \,\sqrt{{h}_{0}({p}_{{\rm{c}}}^{\star })+\displaystyle \frac{{\bar{\nu }}_{{\rm{s}}1}}{2}(\sqrt{1+2\delta }-1)\displaystyle \frac{5{x}_{\mathrm{rad}}+2}{{x}_{\mathrm{rad}}+1}}\\ \quad \approx \,\left({h}_{0}({p}_{{\rm{c}}}^{\star })+\displaystyle \frac{{\bar{\nu }}_{{\rm{s}}1}}{2}(\sqrt{1+2\delta }-1)\displaystyle \frac{3{x}_{\mathrm{rad}}+2}{{x}_{\mathrm{rad}}+1}\right),\end{array}\quad \end{eqnarray} \tag{ C.9 }$$

where the last approximation is a matching between the behavior at ${x}_{\mathrm{rad}}\gg 1$ and ${x}_{\mathrm{rad}}\ll 1$ for $2{h}_{0}({p}_{{\rm{c}}}^{\star })\,\gg {\bar{\nu }}_{{\rm{s}}1}(\sqrt{1+2\delta }-1)$, i.e. screening effects dominate over radiation reaction effects. This assumption also motivates the approximation

$$\begin{eqnarray}&&\mathrm{ln}{p}_{{\rm{c}}}^{\star }\approx \mathrm{ln}{p}_{{\rm{c}}0}^{\star }\approx \mathrm{ln}({\bar{\nu }}_{{\rm{D}}0}/2{\bar{\nu }}_{{\rm{s}}1}).\end{eqnarray} \tag{ C.10 }$$

Finally, the effective critical field (C.2) is the mean of (C.7) and (C.9):

$$\begin{eqnarray}\displaystyle \frac{{E}_{{\rm{c}}}^{\mathrm{eff}}}{{E}_{{\rm{c}}}} & \approx & {h}_{0}({p}_{{\rm{c}}}^{\star })+{\bar{\nu }}_{{\rm{s}}1}(\sqrt{1+2\delta }-1)\\ & \approx & {\bar{\nu }}_{{\rm{s}}0}+{\bar{\nu }}_{{\rm{s}}1}\left[\left(1+\displaystyle \frac{{\bar{\nu }}_{{\rm{D}}1}}{{\bar{\nu }}_{{\rm{D}}0}}\right)\mathrm{ln}\displaystyle \frac{{\bar{\nu }}_{{\rm{D}}0}}{2{\bar{\nu }}_{{\rm{s}}1}}+\sqrt{2\delta +1}\right].\end{eqnarray} \tag{ C.11 }$$

For δ in equation (C.4), we again approximate $\mathrm{ln}{p}_{{\rm{c}}}^{\star }$ using (C.10) but also neglect the ${\bar{\nu }}_{{\rm{D}}1}$ terms compared to ${\bar{\nu }}_{{\rm{D}}0}$, which is motivated both by the smallness of ${\bar{\nu }}_{{\rm{D}}1}$ compared to ${\bar{\nu }}_{{\rm{D}}0}$ and the fact that (C.10) overestimates $\mathrm{ln}{p}_{{\rm{c}}}^{\star }$ if the effect of radiation reaction is significant. Accordingly, we obtain

$$\begin{eqnarray}&&\delta \approx \displaystyle \frac{{\bar{\nu }}_{{\rm{D}}0}}{{\bar{\nu }}_{{\rm{s}}1}^{2}}\left(\displaystyle \frac{{\bar{\nu }}_{{\rm{D}}0}{\tau }_{\mathrm{syn}}^{-1}}{{E}_{{\rm{c}}}^{\mathrm{eff}}/{E}_{{\rm{c}}}}+{\phi }_{\mathrm{br}0}+{\phi }_{\mathrm{br}1}\mathrm{ln}\displaystyle \frac{{\bar{\nu }}_{{\rm{D}}0}}{2{\bar{\nu }}_{{\rm{s}}1}}\right).\end{eqnarray} \tag{ C.12 }$$

Equation (C.11) is a not in a closed form since δ depends on ${E}_{{\rm{c}}}^{\mathrm{eff}}$, but an accurate approximation is obtained after one iteration of (C.11) and (C.12). This is shown in a comparison with the full numerical solution to (22) in figures 2 and 3.