Stability analysis of runaway-driven waves in a tokamak

Runaway electrons (REs) have attracted increased attention in recent years as a potential risk factor for ITER where they can develop easily during plasma disruptions [1]. A very high RE current of several MA can be produced due to a large magnetic energy stored in the plasma. To minimize the detrimental impact of disruptions, such runaways have to be mitigated to avoid massive losses to the wall. Micro-instabilities are potentially able to cause enhanced scattering of the runaways and thereby facilitate the mitigation process. This points to the need to examine plausible micro-instabilities in ITER systematically.

In the early tokamaks, such as TM-3, T-6, TFR and others, the 'fan' instability [2, 3] was observed frequently in the presence of REs. Parail and Pogutse [4, 5] have attributed this phenomenon to excitation of magnetized electron plasma waves. Their local linear and quasilinear analysis has been instrumental in understanding experimental data including enhanced scattering of the REs. Yet, the initial local theory does not appear to be fully sufficient for runaway stability analysis. The point is that it does not cover the effect of plasma non-uniformity on the excited waves. This includes any system size constraints on the wave growth. Also, the analysis of magnetized plasma waves needs to be extended to other potentially unstable modes. References [6–8] reflect the efforts motivated by these unsettled issues. In particular, they emphasize the role of whistler waves and suggest an estimate for their convective damping, in addition to collisional dissipation. However, the expression for the instability threshold obtained in [6–8] is problematic for two reasons. First, the collisional damping rate for whistlers is overestimated. Second, the explanation of convective damping misses the possibility of wave ducting due to internal reflection that enables amplification of the wave over multiple radial passes.

It is appropriate to mention that the characteristic wavelengths of the runaway-driven modes are typically smaller than the plasma radius, which suggests the use of the Wentzel–Kramers–Brillouin (WKB) approximation rather than a full-wave description of the waves of interest. We implement this approach in a ray-tracing code COIN (convective instability) that is designed to examine kinetic instabilities of a runaway beam in a tokamak for any given equilibrium configuration of the plasma and any distribution function of the RE. The code evaluates an amplification factor of the wave by integrating the kinetic drive and collisional damping rate along the wave-packet trajectory in the plasma. It incorporates refined expressions for the drive and damping and covers appropriately convective aspects of the instability as well as the effect of linear mode coupling in non-uniform plasmas. The COIN code can thereby serve as a convenient computational module that is supposed to be used together with independent calculations of the plasma and runaway beam equilibria.

This paper describes the physics content of the code, verification tests, and illustrative calculations of kinetic instability thresholds for ITER-relevant parameters. The paper is organized in the following way.

In section 2 we consider the dispersion relation for high frequency electron plasma waves and identify the candidate instabilities. In sections 3 and 4 we discuss collisional damping and kinetic drive of the wave, respectively. Section 5 presents an evaluation of the kinetic drive required for code benchmarking and verification. In section 6, we compare our local stability calculation with experimental data and discuss the local instability threshold for ITER-relevant parameters. Sections 7 and 8 describe the ray-tracing procedure and discuss internal reflection of the waves from the plasma boundary and the transformation of whistlers into magnetized plasma waves. In section 9 we discuss instability thresholds in DIII-D experiments where REs form during plasma disruptions initiated by Ar pellet injection [9]. Section 10 introduces a message passing interface (MPI) version of the COIN code that allows fast statistical analysis of large number of waves. It also contains an illustrative analysis for ITER-relevant parameters. Section 11 is a brief summary of the results.

In this section we consider waves that can be excited by REs in cold magnetized plasmas. The dispersion relation $\omega(\bi k)$ for any such wave must allow wave–particle resonance. The energy and momentum conservation laws in the process of wave emission [10] show that only two resonances need to be considered for runaways with a small angular spread: the anomalous Doppler resonance,

$$\begin{equation} \omega -{{k}_{\parallel}}V_\|+\frac{\left| {{\omega}_{\rm c}} \right|}{\gamma}=0, \end{equation} \tag{ 1 }$$

and the Cherenkov resonance,

$$\begin{equation} \omega -{{k}_{\parallel}}V_\|=0, \end{equation} \tag{ 2 }$$

where k_|| is the parallel (to the unperturbed magnetic field B) component of the wave vector, V is the particle velocity, $\gamma \equiv 1/\sqrt{1-V^2/c^2}$ is the relativistic factor, and |ω_c| ≡ |eB/mc| is the electron cyclotron frequency.

Of these two, the anomalous Doppler resonance is the one that can produce significant beam scattering, because it involves partial conversion of the electron parallel motion into gyro-motion.

The wave equation and the corresponding dispersion equation in cold magnetized plasma have the form [11, 12]

$$\begin{equation} \left[ {{k}_{\alpha}}{{k}_{\beta}}{{c}^{2}}-{{\delta}_{\alpha \beta}}{{k}^{2}}{{c}^{2}}+{{\omega}^{2}}{{\varepsilon}_{\alpha \beta}} \right]{{E}_{\beta}}=0, \end{equation} \tag{ 3 }$$

$$\begin{equation} \det \left[ {{N}_{\alpha}}{{N}_{\beta}}-{{\delta}_{\alpha \beta}}{{N}^{2}}+{{\varepsilon}_{\alpha \beta}} \right]=0, \end{equation} \tag{ 4 }$$

where

$$\begin{equation} {{\varepsilon}_{\alpha \beta}}(\omega )\equiv \varepsilon {{\delta}_{\alpha \beta}}+ig{{e}_{\alpha \beta \gamma}}{{b}_{\gamma}}+\left( \eta -\varepsilon \right){{b}_{\alpha}}{{b}_{\beta}} \end{equation} \tag{ 5 }$$

is a dielectric tensor, k_α is the wave vector, E_β is the polarization vector, N_α ≡ k_αc/ω is the refractive index, and b_α ≡ B_α/B is the unit vector along the magnetic field. In the absence of collisions, the dielectric tensor is Hermitian, and its components are [12]:

$$\begin{eqnarray} \varepsilon = \varepsilon^{\rm H} &\equiv& 1- \sum_{e;i} \frac{\omega _{\rm p}^{2}}{{{\omega}^{2}}-\omega _{\rm c}^{2}}\nonumber\\ g = g^{\rm H} &\equiv& - \sum_{e;i} \frac{{{\omega}_{\rm c}}}{\omega}\frac{\omega _{\rm p}^{2}}{{{\omega}^{2}}-\omega _{\rm c}^{2}}\\ \eta = \eta^{\rm H} &\equiv& 1- \sum_{e;i} \frac{\omega _{\rm p}^{2}}{{{\omega}^{2}}},\nonumber \end{eqnarray} \tag{ 6 }$$

where summation involves electrons and ions (note that ω_c is negative for the electrons). Equations (3) and (4) permit the following representation for the polarization vector:

$$\begin{eqnarray} &&{\bit E}=\frac{{{E}_{x}}\left( {\bit N}-{\bit b}\left( {\bit N}\cdot {\bit b} \right) \right)+{{E}_{y}}\left[ {\bit b}\times {\bit N} \right]}{\sqrt{{{N}^{2}}-{{\left( {\bit N}\cdot {\bit b} \right)}^{2}}}}+{{E}_{z}}{\bit b}\nonumber\\ &&{{E}_{x}}=1 \nonumber\\ &&{{E}_{y}}=i\frac{g}{\varepsilon -{{N}^{2}}}\\ &&{{E}_{z}}=-\frac{\left( {\bit N}\cdot {\bit b} \right)}{\eta -{{N}^{2}}+{{\left( {\bit N}\cdot {\bit b} \right)}^{2}}}\sqrt{{{N}^{2}}-{{\left( {\bit N}\cdot {\bit b} \right)}^{2}}}. \nonumber \end{eqnarray} \tag{ 7 }$$

Here and further the z index denotes parallel to the magnetic field component, whereas the x and y indices refer to the transverse components.

The dispersion relation (4), together with the resonance condition (1), determine the parallel (N_||) and transverse (N_⊥) components of the refractive index as a function of ω for given values of ω_p/ω_c and relativistic factor γ. The condition that N_⊥(ω) is a real quantity determines the range of allowable wave frequencies. Figure 1 shows a sample plot of N_⊥(ω) in the limit of immobile (infinitely heavy) ions. This figure illustrates the case of relatively strong magnetic field $(\omega_{\rm p}^2/\omega_{\rm c}^2\ll 1)$, in which the lower frequency branches represent magnetized plasma oscillations with

$$\begin{equation} \omega ={{\omega}_{\rm p}}\frac{\left| {{k}_{\parallel}} \right|c}{\sqrt{{{k}^{2}}{{c}^{2}}+\omega _{\rm p}^{2}}} \end{equation} \tag{ 8 }$$

and a whistler mode with

$$\begin{equation} \omega =\left| {{\omega}_{\rm c}} \right|\frac{\left| {{k}_{\parallel}} \right|k{{c}^{2}}}{\omega _{\rm p}^{2}\sqrt{1+{{k}^{2}}{{c}^{2}}\omega _{\rm c}^{2}/\omega _{\rm p}^{4}}}. \end{equation} \tag{ 9 }$$

Note that the usual dispersion relation for whistlers has been generalized here to the case of strong magnetic field $(\omega_{\rm p}^2/\omega_{\rm c}^2 \ll 1)$, in which the expression under the square root can differ from unity significantly.

Zoom In Zoom Out Reset image size

The higher frequency branch asymptotes to the upper-hybrid mode with

$$\begin{equation} \omega =\sqrt{\omega _{\rm c}^{2}+\omega _{\rm p}^{2}}. \end{equation} \tag{ 10 }$$

We note that both resonances (anomalous Doppler and Cherenkov) are present for the higher frequency branch, whereas the lower frequency modes can interact with a relativistic beam only via the anomalous Doppler resonance.

Accurate evaluation of collisional damping is essential for runaway driven waves, because this damping mechanism can be strong enough to compete with the runaway drive. It was shown in [13] that electron collisions can be appropriately included in the conductivity tensor by adding an imaginary contribution to the wave frequency ω. Reference [13] suggests that this contribution is

$$\begin{equation} \frac{i8\sqrt{2\pi}\Lambda {{e}^{4}}}{3{{m}^{1/2}}T_{\rm e}^{3/2}}\sum\limits_{Z}{{{Z}^{2}}{{n}_{Z}}}, \end{equation} \tag{ 11 }$$

where Λ is a Coulomb logarithm. However, this expression is not quite accurate, because of a misprint in the collisional operator given by equation (7) in [13], where the integral has a numerical coefficient 4π (after the substitution of equation (9) into (7) in [13]), whereas this coefficient should be equal to 2π, as can be found in the original work by Landau [14] as well as in the later review by Trubnikov [15]. This misprint translates into the overestimation of the imaginary part by a factor of 2. The correct substitution is

$$\begin{equation} \omega => \omega + i\nu_{\rm e}, \end{equation} \tag{ 12 }$$

where ν_e is an electron–ion collisional frequency

$$\begin{equation} {\nu_{\rm e}}=\frac{4\sqrt{2\pi}\Lambda {{e}^{4}}}{3{{m}^{1/2}}T_{\rm e}^{3/2}}\sum\limits_{Z}{{{Z}^{2}}{{n}_{Z}}}. \end{equation} \tag{ 13 }$$

The simple recipe, offered by [13], can be misinterpreted as an evidence that the collisional damping rate is always equal to ν_e. In particular, this issue has affected calculations of the runaway instability threshold in [6] and [8]. The missing point in those calculations is that the substitution ω => ω + iν_e is permissible in the conductivity tensor only but not in the entire wave dispersion relation³.

In the limit of ω ≫ ν_e, the collision frequency adds a small anti-Hermitian part $\varepsilon_{\alpha \beta}^{\rm A}$ to the dielectric tensor (5), so that the expressions for ε, g and η change from (6) to ε = ε^H + ε^A, g = g^H + g^A and η = η^H + η^A, where the anti-Hermitian corrections are given by

$$\begin{eqnarray} {{\varepsilon}^{\rm A}}&=&i\frac{{{\nu}_{\rm e}}}{\omega}\frac{\omega _{\rm pe}^{2}{{\omega}^{2}}+\omega _{\rm ce}^{2}\omega _{\rm pe}^{2}}{{{\left( {{\omega}^{2}}-\omega _{\rm ce}^{2} \right)}^{2}}}\nonumber\\ {{g}^{\rm A}}&=&-2i\frac{{{\nu}_{\rm e}}}{\omega}\frac{\omega _{\rm pe}^{2}\omega \left| {{\omega}_{\rm ce}} \right|}{{{\left( {{\omega}^{2}}-\omega _{\rm ce}^{2} \right)}^{2}}}\\ {{\eta}^{A}}&=&i\frac{{{\nu}_{\rm e}}}{\omega}\frac{\omega _{\rm pe}^{2}}{{{\omega}^{2}}}.\nonumber \end{eqnarray} \tag{ 14 }$$

The collisional damping rate Γ_ν can now be determined via the standard perturbation technique applied to the wave equation (3). This procedure gives the following expression for Γ_ν:

$$\begin{equation} {{\Gamma}_{\nu}}=-i\frac{E_{\alpha}^{*}{{E}_{\beta}}{{\omega}^{2}}\varepsilon _{\alpha \beta}^{A}}{E_{\alpha}^{*}{{E}_{\beta}}\frac{\partial}{\partial \omega}{{\omega}^{2}}\varepsilon _{\alpha \beta}^{H}}, \end{equation} \tag{ 15 }$$

where E_α is the wave polarization vector determined by (7)

Figure 2 shows the resulting plots of Γ_ν as a function of ω for the waves that meet the anomalous Doppler resonance condition. An immediate conclusion that can be drawn from this figure is that the collisional damping rate is much smaller than ν_e for the low-frequency branch of the resonant waves, which makes it easier for the REs to excite such waves. The waves with the weakest collisional damping (solid part of the low-frequency curve in figure 2) are in fact whistlers and the stronger damped low-frequency branch (dashed curve on the plot) represents magnetized electron plasma waves.

Zoom In Zoom Out Reset image size

The kinetic instability drive from runaways can be derived from a known general expression for the plasma dielectric tensor with a relativistic electron component. We will use the dielectric tensor in the form given in [16] which employs spherical coordinates for the unperturbed electron distribution function F(p; θ) in momentum space. The distribution function F is a sum of the cold plasma contributions from electrons and ions

$$\begin{equation} F_{\rm e}=n_\rme\delta({\bit p}),F_{\rm i}=n_{\rm i}\delta({\bit p}), \end{equation} \tag{ 16 }$$

and a kinetic contribution from runaways, F_b, normalized by the condition

$$\begin{equation} 2\pi \int{{{F}_{\rm b}}{{p}^{2}}}{\rm d}p\sin \theta {\rm d}\theta ={{n}_{\rm b}}, \end{equation} \tag{ 17 }$$

where n_e, n_i and n_b are the cold electron density, cold ion density and the RE density, respectively. Recalling expressions (5) and (6) for the components of the cold electron dielectric tensor, we can then present the complete dielectric tensor in the form

$$\begin{eqnarray} \fl {{\varepsilon}_{xx}}&=&\varepsilon +\frac{4\pi {{e}^{2}}}{\omega}\int{{{\rm d}^{3}}p\sum\limits_{n=-\infty}^{\infty}{{{A}_{n}}}}\frac{{{n}^{2}}\omega _{\rm cb}^{2}}{k_{\bot}^{2}} J_n^2\nonumber\\ \fl {{\varepsilon}_{xy}}&=&\! -{{\varepsilon}_{yx}}\!=\! ig\!+\!i\frac{4\pi {{e}^{2}}}{\omega}\int{{{\! \rm d}^{3}}p\sum\limits_{n=-\infty}^{\infty}{{\!\! {A}_{n}}}}\frac{n{{\omega}_{\rm cb}}V\sin \theta}{{{k}_{\bot}}} J_n {J'}_n \nonumber\\ \fl {{\varepsilon}_{yy}}&=&\varepsilon +\frac{4\pi {{e}^{2}}}{\omega}\int{{{\rm d}^{3}}p\sum\limits_{n=-\infty}^{\infty}{{{A}_{n}}}}{{V}^{2}}{{\sin}^{2}}\theta {J'}_n^2\\ \fl {{\varepsilon}_{xz}}&=&{{\varepsilon}_{zx}}=\frac{4\pi {{e}^{2}}}{\omega}\int{{{\rm d}^{3}}p\sum\limits_{n=-\infty}^{\infty}{{{B}_{n}}}}\frac{n{{\omega}_{\rm cb}}V\cos \theta}{{{k}_{\bot}}} J_n^2\nonumber\\ \fl {{\varepsilon}_{yz}}&=&-{{\varepsilon}_{zy}}=-i\frac{4\pi {{e}^{2}}}{\omega}\int{{{\rm d}^{3}}p\sum\limits_{n=-\infty}^{\infty}{{{B}_{n}}}}{{V}^{2}}\sin \theta \cos \theta J_n {J'}_n\nonumber\\ \fl {{\varepsilon}_{zz}}&=&\eta +\frac{4\pi {{e}^{2}}}{\omega}\int{{{\rm d}^{3}}p\sum\limits_{n=-\infty}^{\infty}{{{B}_{n}}}}{{V}^{2}}{{\cos}^{2}}\theta J_{n}^{2}.\nonumber \end{eqnarray} \tag{ 18 }$$

The argument of the Bessel function J_n and its derivative J'_n in (18) is k_⊥V sin θ/ω_cb, where V = p/mγ is the particle velocity and ω_cb = ω_c/γ is the gyro-frequency. The quantities A_n and B_n are defined as

$$\begin{eqnarray} \fl {{A}_{n}}=&\frac{1}{\left( \omega -{{k}_{\parallel}}V\cos \theta -n{{\omega}_{\rm cb}} \right)}\nonumber\\ \fl & \times \left[ \frac{1}{V}\frac{\partial {{F}_{\rm b}}}{\partial p}+\frac{\omega \cos \theta -{{k}_{\parallel}}V}{\omega pV\sin \theta}\frac{\partial {{F}_{\rm b}}}{\partial \theta} \right] \nonumber\\ \fl {{B}_{n}}= &\frac{1}{\left( \omega -{{k}_{\parallel}}V\cos \theta -n{{\omega}_{\rm cb}} \right)}\nonumber\\ \fl & \times \left[ \frac{1}{V}\frac{\partial {{F}_{\rm b}}}{\partial p}+\frac{1}{pV\cos \theta}\left( \frac{n{{\omega}_{\rm cb}}}{\omega \sin \theta}-\sin \theta \right)\frac{\partial {{F}_{\rm b}}}{\partial \theta} \right] \end{eqnarray} \tag{ 19 }$$

Note that there is a typo in the expression for B_n in [16]. That typo has been corrected in (19). Equation (18) agrees with the expression for the dielectric tensor given in [12], except for the incorrect sign of ε_xy, ε_yx, ε_yz, and ε_zy in [12], which is apparently a misprint there.

Given that the runaway density is relatively small, we treat their contribution to the dielectric tensor as a perturbation and express the instability drive Γ_b through the anti-Hermitian part of the dielectric tensor similarly to (15):

$$\begin{equation} {{\Gamma}_{\rm b}}=i\frac{E_{\alpha}^{*}{{E}_{\beta}}{{\omega}^{2}}\varepsilon _{\alpha \beta}^{\rm A}}{E_{\alpha}^{*}{{E}_{\beta}}\frac{\partial}{\partial \omega}{{\omega}^{2}}\varepsilon _{\alpha \beta}^{\rm H}}, \end{equation} \tag{ 20 }$$

where $\varepsilon_{\alpha \beta}^{\rm H}$ and the polarization vector are still given by (6) and (7), but $\varepsilon_{\alpha \beta}^{\rm A}$ is now determined by (18). The anti-Hermitian terms in (18) arise from the poles in A_n and B_n. Combining the contributions from the poles in accordance with the Landau rule we obtain

$$\begin{eqnarray} \fl {{\Gamma}_{\rm b}}=& \left[4{{\pi}^{2}}{{e}^{2}}\int{{\rm d}^{3}}p\sum\limits_{n=-\infty}^{\infty}{{Q}_{n}}\right. \nonumber\\ \fl & \times \left[ V\frac{\partial {{F}_{\rm b}}}{\partial p}+\frac{V}{p}\frac{n{{\omega}_{\rm cb}}-\omega {{\sin}^{2}}\theta}{\omega \cos \theta \sin \theta}\frac{\partial {{F}_{\rm b}}}{\partial \theta} \right]\nonumber\\ \fl & \times\delta \left( \omega -{{k}_{\parallel}}V\cos \theta -n{{\omega}_{\rm cb}} \right)\bigg] \left[\left( 1-E_{y}^{2} \right)\frac{1}{\omega}\frac{\partial}{\partial \omega}{{\omega}^{2}}\varepsilon \right. \nonumber\\ \fl & \left. +2i{{E}_{y}}\frac{1}{\omega}\frac{\partial}{\partial \omega}{{\omega}^{2}}g+E_{z}^{2}\frac{1}{\omega}\frac{\partial}{\partial \omega}{{\omega}^{2}}\eta \right]^{-1}\nonumber\\ \fl {{Q}_{n}}\equiv & {{\left\{\frac{n{{\omega}_{\rm cb}}}{{{k}_{\bot}}V}{{J}_{n}}+{{E}_{z}}\cos \theta {{J}_{n}}+i{{E}_{y}}\sin \theta {{{{J}'}}_{n}} \right\}}^{2}} \end{eqnarray} \tag{ 21 }$$

The anomalous Doppler resonance is represented by the n = 1 term in the sum in (21). We simplify the drive in the small-angle limit by keeping the lowest order terms in θ to obtain

$$\begin{eqnarray} &&\fl {{\Gamma}_{\rm b}}=\frac{{{\pi}^{2}}{{e}^{2}}}{m}\nonumber\\ &&\fl \times\frac{\int{{{\rm d}^{3}}p{{\left\{1+i{{E}_{y}}+{{E}_{z}}\gamma \frac{{{k}_{\bot}}V}{{{\omega}_{\rm c}}} \right\}}^{2}}\delta \left( \omega -{{k}_{\parallel}}V-\frac{{{\omega}_{\rm c}}}{\gamma} \right)\frac{1}{{{\gamma}^{2}}}\frac{{{\omega}_{\rm c}}}{\omega}\theta \frac{\partial F}{\partial \theta}}}{\left( 1-E_{y}^{2} \right)\frac{1}{\omega}\frac{\partial}{\partial \omega}{{\omega}^{2}}\varepsilon +2i{{E}_{y}}\frac{1}{\omega}\frac{\partial}{\partial \omega}{{\omega}^{2}}g+E_{z}^{2}\frac{1}{\omega}\frac{\partial}{\partial \omega}{{\omega}^{2}}\eta} \nonumber\\ \end{eqnarray} \tag{ 22 }$$

Integration by parts over θ transforms (22) to

$$\begin{eqnarray} &&\fl {{\Gamma}_{\rm b}}=\frac{2{{\pi}^{2}}{{e}^{2}}}{m}\nonumber\\ &&\fl \times \frac{\int_{0}^{\infty}{{\rm dp}{{\left\{1+i{{E}_{y}}+{{E}_{z}}\gamma \frac{{{k}_{\bot}}V}{{{\omega}_{\rm c}}} \right\}}^{2}}\delta \left( \omega -{{k}_{\parallel}}V-\frac{{{\omega}_{\rm c}}}{\gamma} \right)\frac{1}{{{\gamma}^{2}}}\frac{\left| {{\omega}_{\rm c}} \right|}{\omega}G}}{\left( 1-E_{y}^{2} \right)\frac{1}{\omega}\frac{\partial}{\partial \omega}{{\omega}^{2}}\varepsilon +2i{{E}_{y}}\frac{1}{\omega}\frac{\partial}{\partial \omega}{{\omega}^{2}}g+E_{z}^{2}\frac{1}{\omega}\frac{\partial}{\partial \omega}{{\omega}^{2}}\eta}\nonumber\\ \end{eqnarray} \tag{ 23 }$$

where

$$\begin{equation} G(p)=2\pi \int{{{p}^{2}}F\theta \rmd\theta} \end{equation} \tag{ 24 }$$

is the distribution function of runaways over absolute values of their momenta. The presence of a δ-function allows further simplification of (23) via integration over momentum p:

$$\begin{eqnarray} \fl {{\Gamma}_{\rm b}}=&\frac{2{{\pi}^{2}}{{e}^{2}}}{m\gamma} \frac{{{\left\{1+i{{E}_{y}}+{{E}_{z}}\frac{{{k}_{\bot}}c}{{{\omega}_{\rm c}}}\sqrt{{{\gamma}^{2}}-1} \right\}}^{2}}mc\gamma G}{\left( 1-E_{y}^{2} \right)\frac{\partial}{\partial \omega}{{\omega}^{2}}\varepsilon +2i{{E}_{y}}\frac{\partial}{\partial \omega}{{\omega}^{2}}g+E_{z}^{2}\frac{\partial}{\partial \omega}{{\omega}^{2}}\eta}\nonumber\\ &\times \frac{1}{\left| \sqrt{\frac{{{\gamma}^{2}}-1}{{{\gamma}^{2}}}}-\frac{{{k}_{\parallel}}c}{{{\omega}_{\rm c}}\gamma} \right|} \end{eqnarray} \tag{ 25 }$$

The value of γ in this expression is determined by the resonance condition

$$\begin{equation} \omega -{{k}_{\parallel}}c\sqrt{1-1/{{\gamma}^{2}}}-\frac{{{\omega}_{\rm c}}}{\gamma}=0, \end{equation} \tag{ 26 }$$

which gives

$$\begin{equation} \gamma =\frac{k_{\parallel}^{2}{{c}^{2}}+\omega _{\rm c}^{2}}{\omega {{\omega}_{\rm c}}+\sqrt{{{\omega}^{2}}\omega _{\rm c}^{2}+\left( k_{\parallel}^{2}{{c}^{2}}-{{\omega}^{2}} \right)\left( k_{\parallel}^{2}{{c}^{2}}+\omega _{\rm c}^{2} \right)}}. \end{equation} \tag{ 27 }$$

As suggested in [17], the energy spectrum of REs is exponential, so that the distribution function of runaways with an average momentum p₀ and negligible angular spread has the form:

$$\begin{equation} F(p,\theta) = \frac{n_{\rm b}}{2\pi} \frac{exp(-\frac{p}{p_0})}{p_0 p^2} 2\delta(\theta^2). . \end{equation} \tag{ 28 }$$

Equations (28) and (24) with p₀ ≈ mcγ give

$$\begin{equation} mc\gamma G \approx \frac{n_{\rm b}}{e} \end{equation} \tag{ 29 }$$

and

$$\begin{eqnarray} &&\fl \frac{{{\Gamma}_{\rm b}}}{\left| {{\omega}_{\rm c}} \right|} \frac{n}{n_{\rm b}}= \frac{\pi}{2e}\frac{\omega _{\rm p}^{2}}{\gamma \left| {{\omega}_{\rm c}} \right|}\nonumber\\ &&\fl \times \frac{{{\left\{1+i{{E}_{y}}+{{E}_{z}}\frac{{{k}_{\bot}}c}{{{\omega}_{\rm c}}}\sqrt{{{\gamma}^{2}}-1} \right\}}^{2}}}{\left( 1-E_{y}^{2} \right)\frac{\partial}{\partial \omega}{{\omega}^{2}}\varepsilon +2i{{E}_{y}}\frac{\partial}{\partial \omega}{{\omega}^{2}}g+E_{z}^{2}\frac{\partial}{\partial \omega}{{\omega}^{2}}\eta}\frac{1}{\left| \sqrt{\frac{{{\gamma}^{2}}-1}{{{\gamma}^{2}}}}-\frac{{{k}_{\parallel}}c}{{{\omega}_{\rm c}}\gamma} \right|} \nonumber\\ \end{eqnarray} \tag{ 30 }$$

Figure 3 shows the dependence of the normalized kinetic drive $\frac{{{\Gamma}_{\rm b}}}{| {{\omega}_{\rm c}} |} \frac{n}{n_{\rm b}}$ on the wave frequency for γ = 25 for the magnetized plasma wave branch (green solid curve) and the whistler branch (red solid curve). Although the magnetized plasma branch exhibits larger values of the drive than the whistler branch in figure 3, the whistler mode may in fact prevail because it has significantly smaller collisional damping and it is also less susceptible to finite Larmor radius of the beam due to its longer wavelength. Indeed, the dotted lines in figure 3 represent the normalized drive calculated for a beam with a small angular spread (we replace the δ-function in (28) by a Gaussian distribution $\exp(-\frac{\theta^2}{\theta_0^2})/2\theta_0^2$ with θ₀ = 0.1 and use (21)). The difference between the solid and dotted curves in figure 3 demonstrates that even a small angular spread of the beam is already sufficient to suppress the magnetized plasma mode dramatically, whereas the drive for whistlers exhibits only a minor change. These additional factors tend to facilitate excitation of whistlers. Note that [4, 5] assume the dominance of magnetized plasma waves whereas later works [6, 8] assume the dominance of whistlers.

Zoom In Zoom Out Reset image size

The linear growth rate Γ for each wave is the difference between the kinetic drive Γ_b (21) and the collisional damping rate Γ_ν (15).

The damping rate is a function of ω, N_∥, ω_p(n_e), ω_c(B), ν_e(T) that is straightforward to calculate. The kinetic drive involves the fast electron distribution function that has to be provided. In what follows, we use an exponential spectrum of runaway momenta (in the form (28)) with a Gaussian angular distribution:

$$\begin{equation} F_{\rm b}(p,\theta) = \frac{n_{\rm b}}{2\pi} \frac{exp(-\frac{p}{p_0})}{p_0 p^2}\frac{exp(-\frac{\theta^2}{\theta_0^2})}{2\theta_0^2}, \end{equation} \tag{ 31 }$$

where the density n_b is determined by the experimental value of the runaway current.

We compare our calculations with observations of the runaway driven instability reported in [2, 3]. The so-called 'fan' instability, described in these papers, develops in a low-density plasma (n_e≈10¹⁹ m⁻³), leading to deceleration of the RE beam. The corresponding theoretical model developed in [4] and [5] attributes this effect to the excitation of magnetized plasma waves. The x-ray measurements show that the REs are accelerated to 0.5–2 MeV before the instability onset. In [3], the transverse energy of the 0.5 MeV electron beam is estimated to be in the range of 7–60 keV, which gives θ₀ ≈ 0.1–0.3.

The Carbon line radiation measurements in [3] show that the temperature of the main body of electrons is in the range of 10–15 eV. The toroidal magnetic field was approximately 3 T in those experiments.

In figure 4, the data-points denote the RE current density and the bulk plasma electron density at which the instability was recorded in different tokamaks [2].

Zoom In Zoom Out Reset image size

To find the instability threshold numerically, we vary the runaway current density for every value of the plasma density, and find the values of ω and N_∥ for the most unstable wave. The threshold runaway current density, plotted with colour curves in figure 4, represents the case when the kinetic drive balances the damping rate.

Figure 4 shows the results for different RE energies for the bulk plasma temperature of 10 eV and the pitch angle θ₀ = 0.1. We observe that the calculated threshold agrees with the experimental data when the average energy is ≈1–2 MeV (i.e., the energy range relevant to the experiments). The same calculation for the runaway beam with broader pitch-angle (θ₀ = 0.3) puts a higher constraint on the runaway current density. In this case, the threshold matches the experimental results if the bulk plasma temperature is set to 15 eV (which is still consistent with experimental data).

Similar calculations of the maximum possible growth rate performed for ITER-relevant parameters reveal excitation of whistlers rather than magnetized plasma waves, which is a noteworthy distinction from what occurred in previous experiments [2] and [3]. The maximum growth rate for whistlers as a function of plasma density and temperature is shown in figure 5. The magnetic field is taken to be 5.2 T, the distribution function for REs is taken in the form (31) with the current density of 1 MA m⁻², average energy of 15 MeV (as estimated in [17]) and the pitch-angle θ₀ of 0.1.

Zoom In Zoom Out Reset image size

The instability of whistlers is apparently convective, which raises the question whether the wave can be amplified significantly along its path in the plasma. We will address this issue by using a WKB description of the wave packets.

In the WKB limit, a wave packet in an inhomogeneous plasma can be viewed as a particle with a Hamiltonian ω(k, r) determined by the local dispersion relation. The wave-vector k and the wave packet position r are the canonical momentum and spatial coordinate, respectively. This approach is commonly known as the ray-tracing and it was used before in many different codes, for example [18, 19] and more recently [20].

The axial symmetry of a tokamak suggests usage of the cylindrical coordinates. The canonical moments in this case are (P_r = k_r, P_φ = r k_φ, P_z = k_z) and the Hamiltonian ω(k_r; P_φ; k_z; r, z) is independent of the azimuthal angle φ. The resulting equations of motion of the packet on the poloidal (r; z) plane are

$$\begin{eqnarray} &&\frac{\rmd k_r}{\rmd t}=-\frac{\partial \omega (k_r; P_\phi; k_z; r, z)}{\partial r} \\ &&\frac{\rmd k_z}{\rmd t}=-\frac{\partial \omega (k_r; P_\phi; k_z; r, z)}{\partial z}\\ &&\frac{\rmd r}{\rmd t}=\frac{\partial \omega (k_r; P_\phi; k_z; r, z)}{\partial k_r}\\ &&\frac{\rmd z}{\rmd t}=\frac{\partial \omega (k_r; P_\phi; k_z; r, z)}{\partial k_z}. \end{eqnarray} \tag{ 32 }$$

The Hamiltonian here is given implicitly by the dispersion relation for cold electron plasma (4)–(6).

We have developed a ray-tracing code COIN to solve (32) and calculate the amplification factor for the wave.

The code uses cylindrical coordinates. The wave packet is characterized by its location in the poloidal plane (r, z) and the three components (N_r, N_φ, N_z) of its refractive index vector N ≡ k c/ω. The magnetic field data are specified by the EQDSK file. The field components are expressed through the poloidal flux, Ψ(r, z):

$$\begin{equation} {{B}_{r}}=-\frac{1}{r}\frac{\partial \Psi}{\partial z};{{B}_{z}}=\frac{1}{r}\frac{\partial \Psi}{\partial r};{{B}_{\phi}}=\frac{{{B}_{0}}{{R}_{0}}}{r}. \end{equation} \tag{ 33 }$$

The stored poloidal flux, Ψ(r, z), is splined with a 2D cubic spline at the initialization stage of the wave trajectory calculation. This procedure is applicable to an arbitrary 2D equilibrium. It automatically provides the values of the cyclotron frequency and its spatial derivatives.

The 1D profiles of electron density, electron temperature, and runaway current density fraction are specified by another data-file. This data is stored as a function of normalized poloidal flux or minor radius. At the initialization of the calculation, a 2D cubic spline for the plasma frequency is used to evaluate the required spatial derivatives.

Initial parameters for wave packet tracing are

• (r₀, z₀) position in meters (same coordinate system as in the EQDSK file [21]);
• ω/ω_c wave frequency normalized to the local value of the electron cyclotron frequency at (r₀, z₀);
• N_∥ parallel component of the refractive index;
• Θ angle (in degrees) that specifies direction of the transverse component of the wave vector. This angle needs to be introduced because the wave dispersion relation specifies only the absolute value of N_⊥ but not the direction of k_⊥. We assign Θ = 0 to the horizontal direction of k_⊥, so that k_⊥ has no z component at Θ = 0 (see figure 6).

Zoom In Zoom Out Reset image size

Axial symmetry of the equilibrium ensures conservation of the canonical angular momentum, rk_φ. The toroidal component of the refractive index is therefore evaluated as $N_\phi = N_\phi^0 r_0/r$, where $N_\phi^0$ and r₀ are the initial values of N_φ and r. In order to advance other quantities (r; z; N_r;N_z), the 5th order Cash–Karp Runge–Kutta method is implemented.

The amplification factor K of the propagating wave is a time integral of the local growth rate (Γ_b − Γ_ν) along the wave trajectory, i.e.

$$\begin{equation} K = \int (\Gamma_{\rm b}-\Gamma_\nu) \,\rmd t, \end{equation} \tag{ 34 }$$

where Γ_b and Γ_ν are defined by (21) and (15), where an arbitrary number of resonances can be taken into account in equation (21). All calculations presented in this paper include 4 resonances with n ranging from 0 to 3 in equation (21). An arbitrary distribution function of REs can be used as an input in these calculations. The time integration in (34) takes into account variations of the plasma and runaway parameters along the wave trajectory as well as the evolution of the wave vector that may break the wave–particle resonance condition.

The dispersion relations for whistlers and magnetized plasma waves preclude their escape from the plasma in the WKB limit. Consequently, the wave packets should bounce back and forth radially inside the plasma due to internal reflection from the low-density outer layers.

Considering a radially non-uniform plasma cylinder, we observe from the dispersion relation for whistlers (9) that $kc/\omega_{\rm p}^2 = {\rm const}$ along the wave-packet trajectory. On the other hand, radial propagation of the wave packet requires k_r to be real, so that $k=\sqrt{k_\parallel^2+k_\phi^2+k_r^2} \geq |k_\parallel|$. This means that a wave packet launched with k = k₀ and k_∥ = k_∥0 at an inner magnetic surface where the plasma density is n₀ will be reflected back before it reaches an outer surface where the density is n₀|k_∥0|/k₀. As a result, the low-density outer area can act as a reflection layer for the potentially unstable runaway-driven waves. An example of a wave packet trajectory is shown in figure 7. The wave experiences multiple reflection from the plasma boundary. It should, however, be mentioned that there is also partial absorption of the wave at the boundary, due to collisional damping. Sensitivity of collisional absorption to the profiles of electron density and temperature near the boundary calls for a careful treatment of the boundary layer. In particular, the existence of halo currents may alter collisional damping significantly. Also, density fluctuations at the boundary may need to be accounted for as they effectively change the mirror-like reflection of the wave packet to diffusive reflection with a resulting randomness of the wave-packet trajectory.

Zoom In Zoom Out Reset image size

The dispersion relation under consideration (4) is biquadratic. It gives two values of $N_\perp^2$ for a given parallel refractive index N_||. The lower root represents whistlers; the higher root - magnetized plasma waves. The situation, when the discriminant of the dispersion relation vanishes and the two roots coincide does not produce any singularity on the trajectory, but the two branches become indistinguishable at this point. As a result, a whistler wave coming to the zero discriminant location can continue its path as a magnetized plasma wave and vice versa. We note that these transitions take place without any violation of the WKB approximation, which distinguishes them from many other mode conversion processes described in [22]. The processes described in [22] typically involve an intersection of two dispersion curves at some location on the wave trajectory. This is not the case in our problem. The point of zero discriminant is not a brunching point on the trajectory in our case. And it is not a singular point in the WKB description. All derivatives (32) are uniquely defined there, and there is only one possibility for wave evolution (represented by the numerical solution). An example of such transition is shown in figure 7(b).

In this example, we follow the trajectory of a wave with initial parameters corresponding to the whistler branch of the dispersion relation. The solid line represents the results of the Runge–Kutta integration for N_⊥ along the trajectory. We also calculate N_|| on this trajectory (figure 7(a)) and use it to plot the two roots of (4). These two roots are shown by the dashed and the dotted curves in figure 7(b). These curves do not necessarily represent any real trajectory because they are not obtained from the evolution equations self-consistently. In this particular example, the wave evolves along the solid curve and the second value of N_⊥ corresponds to the second 'possible' wave for the given plasma parameters, ω and N_||. The presented example shows that both roots of the dispersion relation have the same value at, for example, t ≈ 0.85 µs and that the trajectory switches from one to another branch of the dispersion relation. Symmetry of the trajectory with respect to time reversal shows that the transformation can go in both directions. It is therefore possible for magnetized plasma waves to transform into whistlers (as seen in figure 7 at t ≈ 0.75 µs). In fact, some waves can experience multiple transformations along their path.

Note that magnetized electron plasma waves have higher collisional damping than whistlers (see figure 2), which indicates that conversion of whistlers into magnetized waves can limit the lifetime of whistlers.

In DIII-D, a beam of REs forms during a plasma disruption initiated by Ar pellet injection [9]. This beam can be sustained for an extended time (≈40 ms). Spectroscopic measurements provide the ionization state of impurities, bulk plasma electron density and temperature profiles. Measurements of the runway electron energy show a mean energy of several MeV and a maximum energy of 30–40 MeV. The anisotropy of the high-energy part of the beam is estimated to be very high, with an average pitch angle θ₀ ≈ 0.2. These data provide sufficient information for runaway stability analysis. Micro-instabilities have never been observed in the DIII-D device, and it is interesting to find out whether theory predicts an experimentally accessible parameter range where these instabilities should develop.

The required post-disruption profiles are shown in figure 8. These profiles (obtained via 1D transport modelling) are consistent with those measured in DIII-D [9].

Zoom In Zoom Out Reset image size

It is important that the post-disruption temperature and electron density profiles depend significantly on the presence of the runaway current. The runaways heat the bulk electrons and ionize neutrals in the core after the decay of the bulk plasma current. In contrast, the outer area with low runaway current has very low electron density (as seen from the n_e profile in figure 8). The resulting density gradient provides an efficient reflection for whistler waves.

Reference [9] provides experimental data on the runaway energy spectrum in DIII-D. For the stability analysis, we approximate this spectrum by an analytical function shown in figure 9.

Zoom In Zoom Out Reset image size

Although the pitch-angle of most energetic REs was estimated as ∼0.2 (10°) in DIII-D experiments, it is impossible for the majority of runaways to have such a small pitch angle, because they would then carry an excessively large current ≈37 MA m⁻². In other words, the angular distribution of the runaways appears to be energy-dependent with a roughly isotropic majority. It should be noted, however, that this distribution function was measured during the plateau phase of the discharge, which is characterized by a relatively low electric field. The electric field tends to be much stronger during the avalanche development, causing much stronger anisotropy of the runaways at that phase.

It is instructive to examine the stability of the runaway beam with the experimental energy spectrum (figure 9), implying that the average pitch-angle of the whole beam is θ₀ = 0.2. This distribution function does not represent the DIII-D reality, but it sets a constraint on the kinetic drive. If the drive for the narrow beam is smaller than collisional damping, one should not expect any instability from the realistic (nearly isotropic) distribution.

Due to the small RE beam radial size and small variation of plasma parameters inside the beam the wave vector remains nearly constant within the beam. Consequently, the local calculations and the ray-tracing calculations give very close results because the profiles shown in figure 8 lock the wave packet within the beam. We thus use only a local analysis in the DIII-D case.

In DIII-D, magnetized waves are driven primarily by moderately energetic electrons (≈1 MeV). One might therefore expect a larger kinetic drive for these waves than in the case of an exponential distribution function (31), due to a larger number of resonant electrons (see figure 9). Nevertheless, our calculation shows no positive growth rate for the magnetized waves at any background plasma parameters. The reason is that the positive kinetic drive from the anomalous Doppler resonance is overcompensated by strong absorption from the Cherenkov resonance.

This stabilizing effect should be even stronger for the actual distribution function where the low energy electrons are nearly isotropic, which allows us to conclude that the runaway distribution in DIII-D prohibits excitation of the magnetized waves. Note that it is the magnetized wave that is responsible for the instability described in [2] and [3]. The reason is that the distribution function in those experiments is apparently different from the distribution function measured in [9].

For the whistler waves, the role of the Cherenkov resonance is much less significant. As a result, whistler waves driven by resonant electrons with energies under 10–15 MeV are close to the instability threshold but the collisional damping rate at the experimental temperature is still too high to allow the development of the wave. The whistlers growth rate becomes positive only if we artificially increase the background plasma temperature in our calculations to ∼5 eV (see figure 10).

Zoom In Zoom Out Reset image size

These waves are primarily driven by electrons with 3–5 MeV. Conversion of whistlers into magnetized plasma waves may limit this growth, but the measured distribution function is still likely to generate whistler waves if the high-energy part (>3 MeV) of this distribution function is strongly anisotropic. For a more realistic beam, even higher background electron temperature is needed for the instability to appear.

The avalanche-type distribution function (see figure 9) exhibits a lower density of resonant electrons for whistler waves. Predictably, the calculated threshold for the avalanche-type distribution function shows that even higher background plasma temperature is needed for the whistler wave instability (>10 eV). On the other hand, the damping effect from the Cherenkov resonance is much lower for the avalanche-type distribution. As a result, the magnetized waves become unstable if we lower the density to 0.04 × 10²⁰ m⁻³.

We observe that amplification of a whistler wave packet is likely to require multiple bounces in the plasma. Also, the wave packet trajectory presented in figure 7 suggests that the lifetime of whistlers may be limited by their conversion into magnetized plasma waves. Because collisional damping of the magnetized wave is very high, it is likely that the wave will decay soon after conversion. It is apparently essential to predict the integrated effect of such conversion. The wave trajectory is very sensitive to plasma parameters. These parameters can fluctuate, especially near the edge, which virtually prohibits an exact calculation of the trajectory over an extended time.

The uncertainty introduced by plasma fluctuations calls for a statistical analysis of many trajectories to predict the instability onset reliably.

To perform such analysis, the COIN code has an MPI version. This enables fast analysis of several million wave trajectories to determine the most favourable conditions for instabilities. We take into account that kinetic drive and collisional damping are additive quantities for a particular wave trajectory. The collisional damping is proportional to the electron collision frequency ν that scales with temperature as T^−3/2. This allows us to calculate contributions to the amplification factor from the kinetic drive and collisional damping separately for a reference temperature and then find the minimal temperature required for instability using the scaling:

$$\begin{equation} {{T}_{\min}}={{T}_{\rm ref}}{{\left( \frac{\int{{{\Gamma}_{\rm b}}\rmd t}}{\int{\Gamma _{\nu}^{\rm ref}\rmd t}} \right)}^{-\frac{2}{3}}}. \end{equation} \tag{ 35 }$$

As an illustration of the statistical analysis, we analyze an ITER RE beam that may be generated during a disruption with the COIN code. We use plasma profiles during a vertical displacement event (VDE) with REs calculated by the disruption version of the DINA code [23]. An ITER 15 MA inductive scenario was chosen as an initial condition for the disruption. A flat Argon profile of 0.4 × 10²⁰ m⁻³ was introduced at the thermal quench phase to keep the plasma temperature low. An arbitrary moment during a VDE is picked to illustrate the analysis. The plasma shape changes slightly by the time most of the current is carried by REs. A better understanding of the RE distribution function is still required for conclusive predictions. For this demonstration, an exponential RE spectrum with average energy of 15 MeV was selected. The density profile was modified at the plasma edge to provide a reflection layer. The temperature profile was flattened for better understanding of the temperature dependence of the instability onset. The resulting plasma profiles are shown in figure 11.

Zoom In Zoom Out Reset image size

A set of random initial coordinates (r, z, ω, N_||, Θ) is generated for every whistler wave. r and z are generated uniformly inside the plasma boundary. Θ is a random number in a range [0 : 2π]. ω and N_|| are generated within ranges broad enough to cover the area of spectrum with a positive growth rate in every (r, z) of the plasma. For every given wave, the code determines a part of the trajectory with a maximum amplification factor:

$$\begin{equation} \forall {{t}_{1}},{{t}_{2}}\in ({{t}_{1}}>{{t}_{2}})\to {{K}_{\max}}=\max \left( K \bigg|_{{{t}_{2}}}^{{{t}_{1}}} \right). \end{equation} \tag{ 36 }$$

The wave is recorded for subsequent analysis to determine whether this factor is higher than a certain value. The results of such calculations are presented in figures 12–14. Figure 12 shows the time integrated kinetic drive and a reference collisional damping for the calculated trajectories versus the duration of the wave packet growth. The maximum value of the time averaged growth rate is comparable to the value of the maximum local linear growth rate for certain waves. Using (35), we find that the maximum amplification factor is around 10 at the plasma background temperature of T_min = 22 eV (figure 13).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 14 shows a histogram of frequencies and initial values of the parallel refractive indices for the amplified waves.

The local calculations predict that the first positive growth rate at the plasma core in these conditions should appear at ∼11 eV background plasma temperature (figure 5). However this temperature is too low for significant amplification of any wave packet along its trajectory. Using convective analysis, we find that the minimal background temperature for the amplification of the whistler waves is s2 times higher than predicted by local analysis.

Figure 15 shows an example of strongly unstable waves in ITER at T_e = 22 eV.

Zoom In Zoom Out Reset image size

A ray-tracing code has been developed and benchmarked for linear stability analysis of REs in realistic tokamak geometry. The code enables fast stability assessment for any given input distribution function of REs. Whistlers and magnetized plasma waves are two primary candidates for high frequency kinetic instabilities. Radial non-uniformity of the tokamak plasma creates a cavity for these waves. A wave packet trapped in this cavity exhibits multiple conversions of these two modes into each other. Calculations of the instability threshold are consistent with previous experimental observations of the runaway-driven instability in several tokamaks (T-3, T-6, TFR, and T-10). The most unstable mode in this case is the magnetized plasma wave.

Based on the experimentally measured parameters of the RE beam in DIIID, our code predicts robust stability of such a beam with respect to high-frequency kinetically driven modes. Preliminary calculations for ITER indicate that the minimal temperature for the whistler instability onset is ∼22 eV in the case of a perfectly reflective plasma boundary. The calculations were done for a flat temperature profile, runaway current of 12 MA (j_re ≈ 1 MA/m²), n_e ≈ 1.3 × 10²⁰ m⁻³ and average runaway energy of 15 MeV. Conclusive predictions for ITER require: (a) a refined description of the runaway distribution function together with the self-consistent decrease of the inductive electric field during saturation of the avalanche; (b) a refined understanding of the plasma temperature evolution during runaway production; and (c) an improved description of wave reflection from the plasma boundary.

This work was supported by the Principality of Monaco/ITER Postdoctoral Fellowship, by ITER under contract ITER-CT-12-4300000273, and by the US Department of Energy Contract No. DEFG02-04ER54742.

The views and opinions expressed herein do not necessarily reflect those of the ITER Organization.