On the relativistic large-angle electron collision operator for runaway avalanches in plasmas

Abstract

Large-angle Coulomb collisions lead to an avalanching generation of runaway electrons in a plasma. We present the first fully conservative large-angle collision operator, derived from the relativistic Boltzmann operator. The relation to previous models for large-angle collisions is investigated, and their validity assessed. We present a form of the generalized collision operator which is suitable for implementation in a numerical kinetic equation solver, and demonstrate the effect on the runaway-electron growth rate. Finally we consider the reverse avalanche effect, where runaways are slowed down by large-angle collisions, and show that the choice of operator is important if the electric field is close to the avalanche threshold.

1 Introduction

Large-angle collisions are associated with large momentum transfers, but their influence can often be ignored in plasma physics, as the cumulative effect of many small-angle deflections are larger by a factor of the Coulomb logarithm, $\ln \unicode[STIX]{x1D6EC}$ (Rosenbluth, MacDonald & Judd 1957; Trubnikov 1965). In many plasmas, e.g. magnetic fusion plasmas and astrophysical plasmas, $\ln \unicode[STIX]{x1D6EC}$ is typically of the order of 10–30. This allows collisions to be accurately accounted for using a Fokker–Planck equation, originally derived for Coulomb interactions by Landau as the small-momentum-transfer limit of the Boltzmann equation (Landau 1965).

A unique situation occurs in runaway acceleration of electrons, where large-angle collisions can play a dominant role even for large $\ln \unicode[STIX]{x1D6EC}$ , as they cause an exponential growth of the runaway density – a runaway avalanche (Sokolov 1979). Runaway is the acceleration of particles in the presence of an electric field which exceeds the critical field $E_{\text{c}}=n_{\text{e}}\ln \unicode[STIX]{x1D6EC}e^{3}(4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}^{2}m_{\text{e}}c^{2})^{-1}$ (Connor & Hastie 1975), where $n_{\text{e}}$ is the electron density, $e$ is the elementary charge, $\unicode[STIX]{x1D700}_{0}$ is the vacuum permittivity, $m_{\text{e}}$ is the electron rest mass and $c$ is the speed of light. Since the collisional drag for superthermal electrons is given by $F_{\text{c}}=eE_{\text{c}}(v/c)^{-2}$ , any electrons with speed greater than the critical speed $v_{\text{c}}=c\sqrt{E_{\text{c}}/E}$ will be accelerated indefinitely, and are hence referred to as runaway electrons (Wilson 1925). Electron runaway occurs in a wide range of plasmas, e.g. in atmospheric discharges (Gurevich, Milikh & Roussel-Dupre 1992), in solar flares (Holman 1985) and in tokamak disruptions when the plasma current changes quickly and a strong electric field is induced (Gill 1993; Jaspers 1993). Due to the large plasma current they would carry, reactor-scale tokamaks such as ITER will be particularly susceptible to the conversion of plasma current to relativistic runaway-electron current by large-angle collisions during disruptions (Rosenbluth & Putvinski 1997). The subsequent uncontrolled loss of a runaway-electron beam could damage plasma-facing components, and the runaways therefore pose a critical threat to the viability of nuclear fusion for energy production (Hollmann et al. 2015).

In a plasma, runaways are mainly generated by two separate mechanisms. When the electrons in the runaway region $v>v_{c}$ are being accelerated, collisional velocity-space diffusion will feed thermal electrons into the runaway region at a steady rate. This primary runaway generation, or Dreicer mechanism (Dreicer 1960), generates new runaways at a rate which is exponentially sensitive to the electric field. The runaway population growth rate was derived in Connor & Hastie (1975), Cohen (1976) and is

(1.1) $$\begin{eqnarray}\left(\frac{\text{d}n_{\text{RE}}}{\text{d}t}\right)_{\text{prim}}\approx \unicode[STIX]{x1D705}\frac{n_{\text{e}}}{\unicode[STIX]{x1D70F}_{\text{c}}}\left(\frac{E}{E_{\text{D}}}\right)^{-(3/16)(1+Z_{\text{eff}})h}\exp \left[-\unicode[STIX]{x1D706}\frac{E_{\text{D}}}{4E}-\sqrt{\unicode[STIX]{x1D702}\frac{(1+Z_{\text{eff}})E_{\text{D}}}{E}}\right],\end{eqnarray}$$

where the undetermined constant factor $\unicode[STIX]{x1D705}$ is of order unity. Here $n_{\text{RE}}$ is the number density of runaways, $\unicode[STIX]{x1D70F}_{\text{c}}=4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}^{2}m_{\text{e}}^{2}c^{3}/(n_{\text{e}}e^{4}\ln \unicode[STIX]{x1D6EC})$ is the relativistic-electron collision time, $E_{\text{D}}=m_{\text{e}}^{2}c^{3}/(e\unicode[STIX]{x1D70F}_{\text{c}}T_{\text{e}})$ is the so-called Dreicer field and $Z_{\text{eff}}=\sum _{i}n_{i}Z_{i}^{2}/\sum _{i}n_{i}Z_{i}$ is the effective ion charge (with the sum taken over all ion species $i$ ). The parameters $h$ , $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D702}$ (not given here) depend on $E/E_{\text{c}}$ and approach unity as $E/E_{\text{c}}$ becomes large (in the non-relativistic limit), and ensure that the growth rate vanishes as $E\rightarrow E_{\text{c}}$ .

A secondary runaway generation mechanism is provided by large-angle collisions, whereby an electron with kinetic energy $\unicode[STIX]{x1D716}=(\unicode[STIX]{x1D6FE}-1)m_{\text{e}}c^{2}>2\unicode[STIX]{x1D716}_{\text{c}}$ can send a stationary target electron into the runaway region in a single collision event while remaining a runaway itself, where $\unicode[STIX]{x1D716}_{\text{c}}$ is the kinetic energy corresponding to the critical speed. Secondary generation, also referred to as avalanche generation due to the resulting exponential growth of the runaway population, generates new runaways at a rate calculated by Rosenbluth & Putvinski (1997) to be approximately

(1.2) $$\begin{eqnarray}\left(\frac{\text{d}n_{\text{RE}}}{\text{d}t}\right)_{\text{ava}}\approx C\frac{n_{\text{RE}}}{2\ln \unicode[STIX]{x1D6EC}\unicode[STIX]{x1D70F}_{\text{c}}}\left(\frac{E}{E_{\text{c}}}-1\right).\end{eqnarray}$$

The function $C=C(E,\,Z_{\text{eff}})$ was shown to be $C=1$ when collisional diffusion is neglected (formally by setting $Z_{\text{eff}}=-1$ ).

While the avalanche growth rate is formally of order $1/\text{ln}\unicode[STIX]{x1D6EC}$ smaller than the primary generation rate, the more favourable scaling with electric field makes it the dominant source of new runaways for sufficiently large runaway populations $n_{\text{RE}}$ or sufficiently small $E/E_{\text{D}}$ , i.e. at low temperature. In the presence of a constant electric field $E$ , with no initial runaway population (apart from a small primary seed), the secondary generation rate will exceed the primary one after approximately one avalanche e-folding time $t_{\text{ava}}\approx 2\ln \unicode[STIX]{x1D6EC}m_{\text{e}}c/[Ce(E-E_{\text{c}})]$ . This corresponds to the time for an initially slow electron to be accelerated to a kinetic energy $E_{\text{k}}\approx \ln \unicode[STIX]{x1D6EC}/C$ MeV (Jayakumar, Fleischmann & Zweben 1993) (neglecting the weak electric-field dependence of $C$ ). Numerically, $t_{\text{ava}}\approx 3.4\ln \unicode[STIX]{x1D6EC}/[C(E-E_{\text{c}})]~\text{ms}$ , with $E$ and $E_{\text{c}}$ in $\text{V}~\text{m}^{-1}$ . If the electric field decreases in magnitude with time, avalanche will become important even earlier. In many practical runaway scenarios, the runaway process will last for multiple $t_{\text{ava}}$ (Gurevich, Milikh & Roussel-Dupre 1994; Gurevich & Zybin 2001; Helander, Eriksson & Andersson 2002), and secondary generation will therefore be the dominant runaway mechanism.

In this work we derive a conservative large-angle (also known as ‘close’ or ‘knock-on’) collision model from the high-energy limit of the linearized relativistic Boltzmann collision integral. We will show how the operators used to model large-angle collisions in previous studies are obtained through various approximations of the Boltzmann collision operator, and how our more general operator resolves issues with previous models and allows the study of new physical effects. In particular, we resolve the issue of double counting of large-angle and small-angle collisions, and show that this development is essential to accurately capture the dynamics. We find that the change to the runaway growth rate due to the new operator is largest during the early stages of the runaway acceleration process, and the likelihood of a given runaway seed transforming into a serious runaway beam can thus potentially be affected. Furthermore, we consider the effect of the inverse knock-on process, where a runaway is slowed down in a single large-angle collision. This effect was recently shown by Aleynikov & Breizman (2015) to be significant for runaway in a near-threshold electric field.

The rest of the paper is organized as follows. In § 2 we introduce the theoretical models describing the large-angle collisions. After giving an overview of the existing models, we present a derivation of the new conservative operator. In § 3 we investigate the effect of the new operator on the runaway growth rate numerically, using the kinetic equation solver CODE (Landreman, Stahl & Fülöp 2014; Stahl et al. 2016). Finally, we summarize our conclusions in § 4.

2 Theoretical models for runaway generation due to large-angle collisions

One of the earliest models for avalanche runaway generation was introduced by Rosenbluth & Putvinski (1997). Due to its simple form, suitable for analytical development, it has been widely used to study the dynamics of an avalanching runaway population (Eriksson & Helander 2003; Smith et al. 2005; Fülöp et al. 2006; Nilsson et al. 2015). Rosenbluth and Putvinski proposed a kinetic equation for the electron distribution of the form

(2.1) $$\begin{eqnarray}\frac{\text{d}f_{\text{e}}}{\text{d}t}=C_{\text{FP}}(f_{\text{e}})+S(f_{\text{e}}),\end{eqnarray}$$

where $\text{d}f_{\text{e}}/\text{d}t$ represents the advective part of the motion, $C_{\text{FP}}$ is the Fokker–Planck collision operator and $S$ a source term representing ‘secondary high-energy electrons knocked out of their orbits by close collisions of a primary relativistic electron with low-energy electrons from the background plasma’ (Rosenbluth & Putvinski 1997). Assuming all existing runaways to be infinitely energetic and having zero pitch angle, they obtained (here adapting their more general result to a homogeneous plasma)

(2.2) $$\begin{eqnarray}S_{\text{RP}}(p,\unicode[STIX]{x1D709},\unicode[STIX]{x1D711})=\frac{n_{\text{RE}}}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}_{\text{c}}\ln \unicode[STIX]{x1D6EC}}\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D709}-\unicode[STIX]{x1D709}_{0})\frac{m_{\text{e}}^{3}c^{3}}{p^{2}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}p}\left(\frac{1}{1-\unicode[STIX]{x1D6FE}}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D709}=\cos \unicode[STIX]{x1D703}=p_{\Vert }/p$ is the pitch-angle cosine, $\unicode[STIX]{x1D6FE}=\sqrt{1+(p/m_{\text{e}}c)^{2}}$ is the Lorentz factor, $\unicode[STIX]{x1D709}_{0}=\sqrt{(\unicode[STIX]{x1D6FE}-1)/(\unicode[STIX]{x1D6FE}+1)}$ and the momentum-space volume element is $p^{2}\,\text{d}p\,\text{d}\unicode[STIX]{x1D709}\,\text{d}\unicode[STIX]{x1D711}$ , with $\unicode[STIX]{x1D711}$ the azimuthal angle of the momentum (the gyroangle). The delta function ensures that secondary electrons are only born on the parabola $p_{\bot }^{2}=2p_{\Vert }m_{\text{e}}c$ in momentum space. In the non-relativistic limit, $p\ll m_{\text{e}}c$ , secondaries are born at perpendicular angles, $p_{\Vert }\approx 0$ , and are prone to trapping in an inhomogeneous magnetic field (Rosenbluth & Putvinski 1997). Away from the magnetic axis of a tokamak, this can lead to a strong reduction in the avalanche growth rate, as recently also demonstrated in detailed numerical simulations by Nilsson et al. (2015).

A more general model was later described by Chiu et al. (1998) (from now on referred to as the Chiu–Harvey operator), which has also been used in runaway studies (Chiu et al. 1998; Harvey et al. 2000; Stahl et al. 2016). Allowing runaway-electron energies to be finite but assuming the runaway pitch angle to be zero, they obtained a knock-on source term

(2.3) $$\begin{eqnarray}S_{\text{CH}}(p,\unicode[STIX]{x1D709},\unicode[STIX]{x1D711})=\frac{1}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}_{\text{c}}\ln \unicode[STIX]{x1D6EC}}\frac{p_{1}^{2}}{m_{\text{e}}cp\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D709}}F(p_{1},t)\unicode[STIX]{x1D6F4}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FE}_{1}^{\prime }),\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F4}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FE}_{1}) & = & \displaystyle \frac{\unicode[STIX]{x1D6FE}_{1}^{2}}{(\unicode[STIX]{x1D6FE}_{1}^{2}-1)(\unicode[STIX]{x1D6FE}-1)^{2}(\unicode[STIX]{x1D6FE}_{1}-\unicode[STIX]{x1D6FE})^{2}}\left((\unicode[STIX]{x1D6FE}_{1}-1)^{2}-\frac{(\unicode[STIX]{x1D6FE}-1)(\unicode[STIX]{x1D6FE}_{1}-\unicode[STIX]{x1D6FE})}{\unicode[STIX]{x1D6FE}_{1}^{2}}\right.\nonumber\\ \displaystyle & & \displaystyle \left.\times \,[2\unicode[STIX]{x1D6FE}_{1}^{2}+2\unicode[STIX]{x1D6FE}_{1}-1-(\unicode[STIX]{x1D6FE}-1)(\unicode[STIX]{x1D6FE}_{1}-\unicode[STIX]{x1D6FE})]\vphantom{\frac{(\unicode[STIX]{x1D6FE}-1)(\unicode[STIX]{x1D6FE}_{1}-\unicode[STIX]{x1D6FE})}{\unicode[STIX]{x1D6FE}_{1}^{2}}}\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D6F4}=(2\unicode[STIX]{x03C0}r_{0}^{2})^{-1}\,\text{d}\unicode[STIX]{x1D70E}/\text{d}\unicode[STIX]{x1D6FE}$ is the normalized Møller differential cross-section for free–free electron–electron scattering (Møller 1932), $r_{0}=e^{2}/(4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}m_{\text{e}}c^{2})$ is the classical electron radius and $\unicode[STIX]{x1D6FE}_{1}$ is connected to $p$ and $\unicode[STIX]{x1D709}$ by the relation

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D709}\equiv \unicode[STIX]{x1D709}^{\ast }(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FE}_{1})=\sqrt{\frac{\unicode[STIX]{x1D6FE}_{1}+1}{\unicode[STIX]{x1D6FE}_{1}-1}\frac{\unicode[STIX]{x1D6FE}-1}{\unicode[STIX]{x1D6FE}+1}}\quad \Leftrightarrow \quad p_{1}=\frac{2p\unicode[STIX]{x1D709}}{1+\unicode[STIX]{x1D709}^{2}-\unicode[STIX]{x1D6FE}(1-\unicode[STIX]{x1D709}^{2})},\end{eqnarray}$$

where a misprint in the original paper incorrectly replaced the $\unicode[STIX]{x1D6FE}-1$ factor with $\unicode[STIX]{x1D6FE}_{1}$ . Since the authors work under the assumption that the runaway pitch angles are negligible, the distribution only appears in the angle-averaged form

(2.6) $$\begin{eqnarray}F(p_{1},t)=\int \,\text{d}\unicode[STIX]{x1D709}_{1}\,\text{d}\unicode[STIX]{x1D711}_{1}\,p_{1}^{2}f_{\text{e}}(p_{1},\unicode[STIX]{x1D709}_{1},t).\end{eqnarray}$$

Both models for large-angle collisions presented above suffer from several defects. In particular, they do not conserve particle number, energy or momentum. In addition, the Rosenbluth–Putvinski model assumes that the incoming particle momentum is infinite, which has the consequence that particles can be created with an energy higher than any of the existing runaways. This assumption is not made in the model derived by Chiu et al. (1998), where the electron energy distribution is properly taken into account, but all incident runaways are still assumed to have zero pitch angle. The magnitudes of both sources ( $S_{\text{RP}}$ and $S_{\text{CH}}$ ) increase rapidly with decreasing momenta and the sources are thus sensitive to the choice of cutoff momentum (introduced to avoid double counting of small-angle collisions).

In the following we will derive a knock-on collision model from the Boltzmann collision integral. As we will show, the model takes into account the full momentum dependence of the primary distribution, and conserves particle number, momentum and energy, while also consistently distinguishing between small- and large-angle collisions, therefore avoiding double counting.

2.1 The Boltzmann collision integral and the Fokker–Planck limit

The Boltzmann collision operator gives the time rate of change of the distribution function due to binary collisions, described by an arbitrary differential cross-section. It can be derived with the following heuristic argument (Montgomery & Tidman 1964; Cercignani & Kremer 2002). The collision operator can be defined as $C_{ab}(f_{a})=(\text{d}n_{a})_{\text{c},ab}/\text{d}t\,\text{d}\boldsymbol{p}$ , where $(\text{d}n_{a})_{\text{c},ab}$ is the differential change in the density of a species $a$ due to collisions with species $b$ , and is defined in terms of the differential cross-section $\text{d}\unicode[STIX]{x1D70E}$ by Lifshitz & Pitaevski (1981) and Cercignani & Kremer (2002)

(2.7)

The first term on the right-hand side, the gain term, describes the rate at which particles $a$ of momentum $\boldsymbol{p}_{1}$ will scatter to momentum $\boldsymbol{p}$ . The second term, the loss term, is the rate at which particles $a$ scatter away from momentum $\boldsymbol{p}$ . Here, we introduced the Møller relative speed and the differential cross-section $\text{d}\unicode[STIX]{x1D70E}_{ab}$ for scattering events $\boldsymbol{p}$ , $\boldsymbol{p}^{\prime }\rightarrow \boldsymbol{p}_{1}$ , $\boldsymbol{p}_{2}$ . The barred quantities are defined likewise, but with $\boldsymbol{p}$ exchanged for $\boldsymbol{p}_{1}$ and $\boldsymbol{p}^{\prime }$ for $\boldsymbol{p}_{2}$ . Since the interactions are viewed as instantaneous, the time labels of the distribution functions have been suppressed for clarity of notation.

The elastic differential cross-section satisfies the symmetry property (known as the principle of detailed balance (Weinberg 2005)), allowing the collision operator to be cast in the commonly adopted symmetric form

(2.8)

where $\boldsymbol{p}_{1}$ and $\boldsymbol{p}_{2}$ (six degrees of freedom) are uniquely determined in terms of $\boldsymbol{p}$ and $\boldsymbol{p}^{\prime }$ by two scattering angles and four constraints by the conservation of momentum and energy,

(2.9) $$\begin{eqnarray}\displaystyle & \boldsymbol{p}_{1}+\boldsymbol{p}_{2}=\boldsymbol{p}+\boldsymbol{p}^{\prime }, & \displaystyle\end{eqnarray}$$

(2.10) $$\begin{eqnarray}\displaystyle & m_{a}\unicode[STIX]{x1D6FE}_{1}+m_{b}\unicode[STIX]{x1D6FE}_{2}=m_{a}\unicode[STIX]{x1D6FE}+m_{b}\unicode[STIX]{x1D6FE}^{\prime }. & \displaystyle\end{eqnarray}$$

From this collision operator, the Fokker–Planck operator, which is often used in plasma physics, can be obtained by a Taylor expansion to second order in the momentum transfer $\unicode[STIX]{x0394}\boldsymbol{p}=\boldsymbol{p}_{1}-\boldsymbol{p}$ (Landau 1936; Akama 1970), motivated by the fact that the cross-section for Coulomb collisions is singular for small deflections. It is then seen that the contribution of small-angle collisions is larger than those of large-angle collisions by a factor of the Coulomb logarithm,

(2.11) $$\begin{eqnarray}\ln \unicode[STIX]{x1D6EC}=\int _{\cot (\unicode[STIX]{x1D703}_{\max }/2)}^{\cot (\unicode[STIX]{x1D703}_{\min }/2)}\frac{\text{d}\unicode[STIX]{x1D706}}{\unicode[STIX]{x1D706}}=\ln \left(\cot \frac{\unicode[STIX]{x1D703}_{\min }}{2}\right),\end{eqnarray}$$

where the maximum centre-of-mass deflection angle for self-collisions is $\unicode[STIX]{x1D703}_{\max }=\unicode[STIX]{x03C0}/2$ (not $\unicode[STIX]{x03C0}$ as it is for unlike-species collisions, or collisions would be double counted) and $\unicode[STIX]{x1D703}_{\min }$ is a cutoff required to regularize the expression, typically chosen as the scattering angle corresponding to impact parameters of order the Debye length, ¹ beyond which particles will not interact because of Debye screening.

Note that by using a total collision operator $C_{\text{FP}}+C_{\text{Boltz}}$ as prescribed by (2.1), the Boltzmann operator has effectively been added twice, although different approximations are used to evaluate the two terms. A subset of collisions will therefore be double counted. One way to resolve this issue is to apply the Fokker–Planck operator only to collision angles smaller than some $\unicode[STIX]{x1D703}_{\text{m}}$ , and the knock-on (Boltzmann) operator for $\unicode[STIX]{x1D703}>\unicode[STIX]{x1D703}_{\text{m}}$ . The Coulomb logarithm used in the Fokker–Planck operator then ought to be changed from (2.11) to

(2.12) $$\begin{eqnarray}\overline{\ln \unicode[STIX]{x1D6EC}}=\ln \unicode[STIX]{x1D6EC}-\ln \left(\cot \frac{\unicode[STIX]{x1D703}_{\text{m}}}{2}\right).\end{eqnarray}$$

When an incident electron of momentum $p$ knocks a stationary electron to momentum $p_{\text{m}}$ , the corresponding centre-of-mass scattering angle $\unicode[STIX]{x1D703}_{\text{m}}$ is given by

(2.13) $$\begin{eqnarray}\cot \frac{\unicode[STIX]{x1D703}_{\text{m}}}{2}=\sqrt{\frac{\unicode[STIX]{x1D6FE}-\unicode[STIX]{x1D6FE}_{\text{m}}}{\unicode[STIX]{x1D6FE}_{\text{m}}-1}}.\end{eqnarray}$$

By using this energy-dependent modification to the Coulomb logarithm, no collisions will be double counted. Indeed, by taking the energy moment of the test-particle collision operator (the sum of Fokker–Planck and Boltzmann), it can be verified that with this choice, the average energy loss rate experienced by a test particle becomes independent of the cutoff $p_{\text{m}}$ , when $p_{\text{m}}\ll m_{\text{e}}c$ .

The number of collisions that are double counted can often be significant when this effect is unaccounted for. Assuming $v_{\text{m}}/c\sim \sqrt{E_{\text{c}}/E}$ to be located at a non-relativistic energy (that is, we assume $E\gg E_{\text{c}}$ ), the modification to the Coulomb logarithm is approximately given by $\ln \sqrt{2(E/E_{c})(\unicode[STIX]{x1D6FE}-1)}$ . For highly energetic electrons with $\unicode[STIX]{x1D6FE}\sim 50$ and $E/E_{\text{c}}\sim 100$ , this corresponds to a change of approximately 5, which – depending on plasma parameters – typically constitutes a relative change to the Coulomb logarithm of 25–50 %.

In principle, as $\unicode[STIX]{x1D703}_{\text{m}}$ approaches the cutoff imposed by Debye screening (or the binding energy of atoms in the case of electrons in neutral media), the Boltzmann operator will account for all collisions and $\overline{\ln \unicode[STIX]{x1D6EC}}=0$ . However, this corresponds to a cutoff momentum smaller than thermal, $p_{\text{m}}\ll p_{\text{Te}}$ , and the assumption of stationary targets is violated when evaluating the operator in the bulk region. In addition, to numerically resolve the Boltzmann operator in a finite-difference scheme, the grid spacing in momentum must be much smaller than $p_{\text{m}}$ , and it is therefore desirable to choose $p_{\text{m}}$ as large as allowed while having a well converged description of the secondary generation rate. The sensitivity of the result to the choice of $p_{\text{m}}$ is investigated in the next section.

In the following we will find it more useful not to work with the symmetric form of the Boltzmann operator given by (2.8), but instead use the alternative given directly from (2.7),

(2.14)

where $\unicode[STIX]{x1D70E}_{ab}(\,\boldsymbol{p},\boldsymbol{p}^{\prime })=\int \,\text{d}\,\boldsymbol{p}_{1}\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ab}/\unicode[STIX]{x2202}\boldsymbol{p}_{1}$ is the total cross-section.

2.2 Derivation of a conservative knock-on operator

For the avalanche problem, one is concerned with the electron–electron Boltzmann operator. We consider the scenario where a small runaway population has been accelerated by an electric field (or other mechanism), leaving a largely intact thermal bulk population. We may then write our electron distribution as $f_{\text{e}}(\,\boldsymbol{p})=f_{\text{Me}}(\,\boldsymbol{p})+\unicode[STIX]{x1D6FF}f_{\text{e}}(\,\boldsymbol{p})$ , where the runaway distribution $\unicode[STIX]{x1D6FF}f_{\text{e}}$ is much smaller than the bulk distribution $f_{\text{Me}}$ , $\Vert \unicode[STIX]{x1D6FF}f_{\text{e}}\Vert \ll \Vert f_{\text{Me}}\Vert$ (for example in terms of number densities $n_{\text{RE}}\ll n_{\text{e}}$ ). We may then linearize the bilinear Boltzmann operator by ignoring terms quadratic in $\unicode[STIX]{x1D6FF}f_{\text{e}}$ , obtaining

(2.15) $$\begin{eqnarray}C_{\text{ee}}^{B}\{f_{\text{e}},f_{\text{e}}\}\approx C_{\text{ee}}^{B}\{f_{\text{e}},f_{\text{Me}}\}+C_{\text{ee}}^{B}\{f_{\text{Me}},f_{\text{e}}\}\equiv C_{\text{boltz}}(\,\boldsymbol{p}),\end{eqnarray}$$

where terms $C_{\text{ee}}\{f_{\text{Me}},f_{\text{Me}}\}$ vanish since $f_{\text{Me}}$ is chosen as an equilibrium distribution. The first term, the test-particle term, describes the effect of large-angle collisions on the runaway electrons as they collide with the thermal population. The second term, the field-particle term, describes the reaction of the bulk electrons as they are being struck by the runaways. Intuitively, one could expect this field-particle term to constitute the avalanche knock-on source. We shall show below that this is indeed the case.

Before giving the explicit forms of the collision operator, we will make one final approximation. We assume that both the incident and outgoing electrons in the large-angle collisions are significantly faster than the thermal speed $v_{\text{Te}}=\sqrt{2T_{\text{e}}/m_{\text{e}}}$ , so that we may approximate the bulk population with a Dirac delta function: $f_{\text{Me}}(\,\boldsymbol{p})\approx n_{\text{e}}\unicode[STIX]{x1D6FF}(\,\boldsymbol{p})$ . The collision operator then takes the form

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle C_{\text{ee}}^{B}\{f_{\text{e}},f_{\text{Me}}\}=n_{\text{e}}\int _{q^{\ast }>p_{1}>q_{0}}\,\text{d}\,\boldsymbol{p}_{1}v_{1}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70E}}_{\text{ee}}}{\unicode[STIX]{x2202}\boldsymbol{p}}f_{\text{e}}(\,\boldsymbol{p}_{1})-n_{\text{e}}v\unicode[STIX]{x1D70E}_{\text{ee}}(\,\boldsymbol{p})f_{\text{e}}(\,\boldsymbol{p}), & \displaystyle\end{eqnarray}$$

(2.17) $$\begin{eqnarray}\displaystyle & \displaystyle C_{\text{ee}}^{B}\{f_{\text{Me}},f_{\text{e}}\}=n_{\text{e}}\int _{p_{1}>q^{\ast }}\text{d}\boldsymbol{p}_{1}v_{1}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70E}}_{\text{ee}}}{\unicode[STIX]{x2202}\boldsymbol{p}}f_{\text{e}}(\,\boldsymbol{p}_{1})-n_{\text{e}}\unicode[STIX]{x1D6FF}(\,\boldsymbol{p})\int _{p_{1}>q_{0}}\text{d}\,\boldsymbol{p}_{1}v_{1}\unicode[STIX]{x1D70E}_{\text{ee}}(\,\boldsymbol{p}_{1})f_{\text{e}}(\,\boldsymbol{p}). & \displaystyle\end{eqnarray}$$

The total cross-section $\unicode[STIX]{x1D70E}_{\text{ee}}(\,\boldsymbol{p})$ is given in (A 9) in appendix A. The limiting momenta $q^{\ast }$ and $q_{0}$ are determined from constraints imposed by conservation laws. For the gain term, i.e. the first term in each equation, energy conservation in each collision reads $\unicode[STIX]{x1D6FE}_{1}=\unicode[STIX]{x1D6FE}+\unicode[STIX]{x1D6FE}^{\prime }-1$ , where $\unicode[STIX]{x1D6FE}$ and $\unicode[STIX]{x1D6FE}^{\prime }$ are the Lorentz factors of the two electrons after the collision. The conditions $\unicode[STIX]{x1D6FE}^{\prime }>\unicode[STIX]{x1D6FE}$ or $\unicode[STIX]{x1D6FE}^{\prime }<\unicode[STIX]{x1D6FE}$ determines whether $\unicode[STIX]{x1D6FE}$ refers to the bulk particle or runaway particle after the collision, respectively (note that the electrons are in fact indistinguishable, but an artificial distinction like this must be performed in order to avoid double counting). We therefore obtain $q^{\ast }$ from setting $\unicode[STIX]{x1D6FE}^{\prime }=\unicode[STIX]{x1D6FE}$ in the conservation law, giving $\unicode[STIX]{x1D6FE}^{\ast }=2\unicode[STIX]{x1D6FE}-1$ , which corresponds to $q^{\ast }=m_{\text{e}}c\sqrt{\unicode[STIX]{x1D6FE}^{\ast 2}-1}$ .

Similarly, we cannot account for all collisions, since we have assumed the bulk particles to be much slower than the outgoing particles. We therefore choose to account only for those collisions where incident and outgoing particles have momenta larger than some $p=p_{\text{m}}\gg p_{\text{Te}}$ . Setting $\unicode[STIX]{x1D6FE}^{\prime }=\unicode[STIX]{x1D6FE}_{\text{m}}$ then yields the lower limit $\unicode[STIX]{x1D6FE}_{0}=\unicode[STIX]{x1D6FE}+\unicode[STIX]{x1D6FE}_{\text{m}}-1$ , corresponding to $q_{0}=m_{\text{e}}c\sqrt{\unicode[STIX]{x1D6FE}_{0}^{2}-1}$ . Note that for the total operator $C_{\text{Boltz}}$ , the two gain terms in (2.16) and (2.17) combine into one integral, taken over all momenta $p_{1}>q_{0}$ . The full expression is thus independent of the parameter $q^{\ast }$ which distinguishes the two outgoing particles (as is expected, since the distinction is not physically relevant for scattering of identical particles).

We can now derive explicit expressions for the collision operator. Since there are only two degrees of freedom in the scattering process (for example two independent scattering angles), the differential cross-section $\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70E}}_{\text{ee}}/\unicode[STIX]{x2202}\boldsymbol{p}$ will invariably contain a delta function. In Møller scattering (relativistic electron–electron scattering), the cross-section is azimuthally symmetric (assuming the electrons to be spin unpolarized) and takes the form

(2.18) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70E}}_{\text{ee}}}{\unicode[STIX]{x2202}\boldsymbol{p}}=\frac{r_{0}^{2}}{(m_{\text{e}}c)^{2}p\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D6FF}(\cos \unicode[STIX]{x1D703}_{\text{s}}-\unicode[STIX]{x1D709}^{\ast })\unicode[STIX]{x1D6F4}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FE}_{1}), & \displaystyle\end{eqnarray}$$

(2.19) $$\begin{eqnarray}\displaystyle & \displaystyle \cos \unicode[STIX]{x1D703}_{\text{s}}=\frac{\boldsymbol{p}_{1}\boldsymbol{\cdot }\boldsymbol{p}}{p_{1}p}\equiv \unicode[STIX]{x1D709}_{s}, & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D703}_{\text{s}}$ is the deflection angle, $\unicode[STIX]{x1D6F4}$ is defined in (2.4) and $\unicode[STIX]{x1D709}^{\ast }$ is defined in (2.5). The delta function enforces the relation between scattering angle and energy transfer that follows from the conservation of 4-momentum. The gain term then takes the form

(2.20) $$\begin{eqnarray}\displaystyle n_{\text{e}}\int \,\text{d}\,\boldsymbol{p}_{1}\,v_{1}\frac{\unicode[STIX]{x2202}\bar{\unicode[STIX]{x1D70E}}_{\text{ee}}}{\unicode[STIX]{x2202}\boldsymbol{p}}f_{\text{e}}(\,\boldsymbol{p}_{1}) & = & \displaystyle \frac{1}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}_{\text{c}}\ln \unicode[STIX]{x1D6EC}}\frac{1}{p\unicode[STIX]{x1D6FE}}\int \,\text{d}p_{1}\frac{p_{1}^{3}}{\unicode[STIX]{x1D6FE}_{1}}\unicode[STIX]{x1D6F4}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FE}_{1})\nonumber\\ \displaystyle & & \displaystyle \times \,\int \,\text{d}\unicode[STIX]{x1D709}_{1}\,\text{d}\unicode[STIX]{x1D711}_{1}\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D709}_{s}-\unicode[STIX]{x1D709}^{\ast })f_{\text{e}}(\,\boldsymbol{p}_{1}).\end{eqnarray}$$

This expression (when choosing integration limits appropriate for the field-particle term) is the generalized ‘knock-on source term’ $S$ , which reduces to the expressions given by Rosenbluth & Putvinski (1997) and Chiu et al. (1998) – (2.2) and (2.3), respectively – using appropriate approximations, as shown in appendix B. This connection has not been acknowledged in previous studies, to the degree that Chiu et al. (1998) incorrectly ascribe the discrepancy between their result and that of Besedin & Pankratov (1986) by the fact that ‘The present expressions are simply a statement of the total rate at which electrons in different velocity-space elements of primary electrons knock a collection of cold bulk electrons into velocity-space elements of the secondary electrons. The expression in (Besedin & Pankratov 1986) uses a Boltzmann-like integral operator’. In fact, as we show in appendix B, the approaches are completely equivalent, and the discrepancy is the result of an error in the calculation of Besedin & Pankratov (1986).

There are multiple ways of carrying out the integration over the delta function; if we assume a distribution function independent of gyroangle, $f_{\text{e}}(\,\boldsymbol{p}_{1})=f_{\text{e}}(p_{1},\cos \unicode[STIX]{x1D703}_{1})$ , a few convenient expressions are given by

(2.21) $$\begin{eqnarray}\displaystyle \int _{-1}^{1}\,\text{d}\unicode[STIX]{x1D709}_{1}\int _{0}^{2\unicode[STIX]{x03C0}}\,\text{d}\unicode[STIX]{x1D711}_{1}\,\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D709}_{s}-\unicode[STIX]{x1D709}^{\ast })f_{\text{e}}(p_{1},\unicode[STIX]{x1D709}_{1}) & = & \displaystyle \int _{0}^{2\unicode[STIX]{x03C0}}\,\text{d}\unicode[STIX]{x1D711}_{\text{s}}\,\,f_{\text{e}}(p_{1},\unicode[STIX]{x1D709}_{1})\nonumber\\ \displaystyle & = & \displaystyle 2\int _{\cos (\unicode[STIX]{x1D703}+\unicode[STIX]{x1D703}^{\ast })}^{\cos (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}^{\ast })}\,\text{d}\unicode[STIX]{x1D709}_{1}\frac{f_{\text{e}}(p_{1},\unicode[STIX]{x1D709}_{1})}{\sqrt{1-\unicode[STIX]{x1D709}^{\ast 2}-\unicode[STIX]{x1D709}_{1}^{2}-\unicode[STIX]{x1D709}^{2}+2\unicode[STIX]{x1D709}^{\ast }\unicode[STIX]{x1D709}_{1}\unicode[STIX]{x1D709}}}\nonumber\\ \displaystyle & = & \displaystyle 2\unicode[STIX]{x03C0}\mathop{\sum }_{L}f_{L}(p_{1})P_{L}(\unicode[STIX]{x1D709})P_{L}(\unicode[STIX]{x1D709}^{\ast }),\end{eqnarray}$$

where we have introduced the quantities

(2.22) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\cos \unicode[STIX]{x1D711}_{\text{s}}={\displaystyle \frac{\unicode[STIX]{x1D709}_{1}-\unicode[STIX]{x1D709}^{\ast }\unicode[STIX]{x1D709}}{\sqrt{1-\unicode[STIX]{x1D709}^{\ast 2}}\sqrt{1-\unicode[STIX]{x1D709}^{2}}}},\\ f_{L}(p)={\displaystyle \frac{2L+1}{2}}\displaystyle \int _{-1}^{1}\,\text{d}\unicode[STIX]{x1D709}\,f_{\text{e}}(p,\unicode[STIX]{x1D709})P_{L}(\unicode[STIX]{x1D709}).\end{array}\right\}\end{eqnarray}$$

In particular the form involving Legendre polynomials $P_{L}$ is a powerful result, as it demonstrates that the linearized Boltzmann operator is diagonal in $L$ , in the sense that if $C_{\text{Boltz}}(\,\boldsymbol{p})=\sum _{L}C_{L}(p)P_{L}(\cos \unicode[STIX]{x1D703})$ , then $C_{L}$ depends only on $f_{L}$ (and not other $f_{l}$ with $l\neq L$ ). This behaviour exhibits the spherical symmetry inherent in scattering on stationary targets. Utilizing this property leads to significant practical gains in terms of numerical computation times. Analogous expressions in terms of Legendre polynomials and the integration over $\unicode[STIX]{x1D711}_{\text{s}}$ were also found by Gurevich & Zybin (2001) for the so-called ionization integral in neutral gases. The form of the integral taken over $\unicode[STIX]{x1D709}_{1}$ was obtained by Helander, Lisak & Ryutov (1993) in the analogous problem of elastic nucleon–nucleon scattering, and equivalent formulations were also recently given by Aleynikov et al. (2014) and Boozer (2015).

The Legendre modes of the collision operator are explicitly given by

(2.23) $$\begin{eqnarray}\displaystyle C_{L}\{f_{\text{e}},f_{\text{Me}}\} & = & \displaystyle \frac{(m_{\text{e}}c)^{-3}}{2\unicode[STIX]{x1D70F}_{\text{c}}\ln \unicode[STIX]{x1D6EC}}\frac{1}{\unicode[STIX]{x1D6FE}p}\int _{q_{0}}^{q^{\ast }}\,\text{d}p_{1}\frac{p_{1}^{3}}{\unicode[STIX]{x1D6FE}_{1}}f_{L}(p_{1})P_{L}(\unicode[STIX]{x1D709}^{\ast })\unicode[STIX]{x1D6F4}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FE}_{1})\nonumber\\ \displaystyle & & \displaystyle -\,\frac{1}{4\unicode[STIX]{x1D70F}_{\text{c}}\ln \unicode[STIX]{x1D6EC}}\frac{v}{c}f_{L}(p)\int _{\unicode[STIX]{x1D6FE}_{\text{m}}}^{\unicode[STIX]{x1D6FE}+1-\unicode[STIX]{x1D6FE}_{\text{m}}}\,\text{d}\unicode[STIX]{x1D6FE}_{1}\unicode[STIX]{x1D6F4}(\unicode[STIX]{x1D6FE}_{1},\unicode[STIX]{x1D6FE}),\end{eqnarray}$$

(2.24) $$\begin{eqnarray}\displaystyle C_{L}\{f_{\text{Me}},f_{\text{e}}\} & = & \displaystyle \frac{(m_{\text{e}}c)^{-3}}{2\unicode[STIX]{x1D70F}_{\text{c}}\ln \unicode[STIX]{x1D6EC}}\frac{1}{\unicode[STIX]{x1D6FE}p}\int _{q^{\ast }}^{\infty }\,\text{d}p_{1}\frac{p_{1}^{3}}{\unicode[STIX]{x1D6FE}_{1}}f_{L}(p_{1})P_{L}(\unicode[STIX]{x1D709}^{\ast })\unicode[STIX]{x1D6F4}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FE}_{1})\nonumber\\ \displaystyle & & \displaystyle -\,\frac{(m_{\text{e}}c)^{-1}}{4\unicode[STIX]{x1D70F}_{\text{c}}\ln \unicode[STIX]{x1D6EC}}\unicode[STIX]{x1D6FF}_{L,0}\frac{\unicode[STIX]{x1D6FF}(p)}{p^{2}}\int _{q_{0}(p_{\text{m}})}^{\infty }\,\text{d}p^{\prime }\frac{{p^{\prime }}^{3}}{\unicode[STIX]{x1D6FE}^{\prime }}f_{0}(p^{\prime })\int _{\unicode[STIX]{x1D6FE}_{\text{m}}}^{\unicode[STIX]{x1D6FE}^{\prime }+1-\unicode[STIX]{x1D6FE}_{\text{m}}}\,\text{d}\unicode[STIX]{x1D6FE}_{1}\unicode[STIX]{x1D6F4}(\unicode[STIX]{x1D6FE}_{1},\unicode[STIX]{x1D6FE}^{\prime }).\qquad\end{eqnarray}$$

Note further that since we only consider those collisions where both the incident and outgoing particles have momenta $p>p_{\text{m}}$ , the gain terms must only be applied for $\unicode[STIX]{x1D6FE}>\unicode[STIX]{x1D6FE}_{\text{m}}$ , while the test-particle loss term is applied for $\unicode[STIX]{x1D6FE}>2\unicode[STIX]{x1D6FE}_{\text{m}}-1$ . In appendix A it is explicitly demonstrated that this collision operator conserves density, momentum and energy.

Figure 1. Illustration of the large-angle collision operators investigated in this study. Darker colours represent larger amplitudes (in arbitrary units), where white and black are separated by 3 orders of magnitude. (a) The distribution function $\log f_{\text{e}}$ with which we evaluate the large-angle collision operator; (c) the Chiu–Harvey operator $\log C_{\text{CH}}$ ; (d) the full field-particle operator $\log C_{\text{boltz}}^{\text{(fp)}}$ (dashed: the line $\unicode[STIX]{x1D709}=\sqrt{(\unicode[STIX]{x1D6FE}-1)/(\unicode[STIX]{x1D6FE}+1)}$ , where the Rosenbluth–Putvinski operator creates knock-ons); (b) the magnitude of the full test-particle operator $\log |C_{\text{boltz}}^{\text{(tp)}}|$ , where the dotted line separates the region of negative contributions (to the right) from the positive contributions (to the left).

A qualitative illustration of the large-angle collision operators discussed here is shown in figure 1. A test runaway distribution (figure 1 a) was generated by applying a constant electric field $E=15E_{\text{c}}$ for a short time $t\approx 0.5\unicode[STIX]{x1D70F}_{\text{c}}$ with $Z_{\text{eff}}=5$ , and the large-angle collision operators were evaluated in the final time step. The figures show a snapshot of where large-angle collisions between runaways and bulk particles create or remove electrons in phase space; comparing 1(c) and 1(d) shows that the Chiu–Harvey operator creates secondary runaways in a significantly smaller region in momentum space than the full field-particle operator, however the total number of secondary runaways created is equal between the models. Figure 1(b) shows the Boltzmann test-particle operator, illustrating the reaction of the already present runaways: they are removed at small pitch angles where the runaway distribution is largest, and placed at larger pitch angles and lower energy. The sum of figure 1(b,d) is the full Boltzmann operator which conserves particle number, momentum and energy.

Note finally that all of the knock-on models described in this paper share the assumption of a stationary bulk, which means that the operators can only be evaluated at speeds much larger than the thermal speed $v_{\text{Te}}$ . Since the sources must be applied for speeds smaller than the critical speed $v_{c}$ in order to accurately capture the runaway rate, the condition $v_{\text{Te}}\ll v_{c}$ limits the electric-field values to $\sqrt{E}\ll \sqrt{E_{\text{D}}/2}$ (effectively forming a lower limit in density and an upper limit in temperature for a given $E$ ). As a consequence, avalanche generation in electric fields large enough for Dreicer generation to be significant is not accessible by the models used here. This limitation would be resolved by accounting for the velocity distribution of the target population in (2.7), resulting in a significantly more complicated operator. However, in many scenarios of interest this is not an issue; the critical velocity tends to be significantly larger than thermal, or is comparable only for a relatively short period of time during which a runaway seed is generated, which then proceeds to grow primarily by avalanche generation.

3 Numerical study of the effect of large-angle collisions

We use the kinetic equation solver CODE (Landreman et al. 2014; Stahl et al. 2016) to compare the various models for the knock-on collision operator. We use CODE to solve the relativistic 0D $+$ 2P kinetic equation for the electron distribution

(3.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}f_{e}}{\unicode[STIX]{x2202}t}+\left\langle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\boldsymbol{p}}\boldsymbol{\cdot }[(\boldsymbol{F}_{L}+\boldsymbol{F}_{S})f_{e}]\right\rangle =C_{\text{ei}}+C_{\text{ee}}+C_{\text{boltz}},\end{eqnarray}$$

where $\boldsymbol{F}_{L}$ is the Lorentz force, and $\boldsymbol{F}_{S}$ is the radiation reaction force associated with synchrotron radiation and the brackets denote averaging over the azimuthal (gyro) angle. $C_{\text{ei}}$ and $C_{\text{ee}}$ are the gyroaveraged Fokker–Planck collision operators for electron–ion and electron–electron collisions, respectively. The thermal bulk population is resolved in the simulations, as well as the relativistic runaway tail. The collision operator $C_{\text{ei}}+C_{\text{ee}}$ and the numerical scheme used are described in Landreman et al. (2014). In the numerical solutions of (3.1) a Dirichlet boundary condition $f_{e}=0$ is imposed at the upper boundary $p_{\max }$ in momentum. The boundary is chosen large enough for the results to be insensitive to variations in $p_{\max }$ or details of the boundary condition; the value of the distribution is naturally negligibly small near the boundary, since it asymptotically tends to decrease exponentially with energy.

First, we will study the sensitivity of the avalanche dynamics to the arbitrary cutoff parameter $p_{\text{m}}$ and investigate the effects of adding the test-particle Boltzmann operator, which restores conservation laws in the knock-on collisions. We then focus on two scenarios: (i) we revisit the classical calculation of the steady-state avalanche growth rate in a constant electric field, (ii) we calculate the runaway growth rate in the near-critical field, accounting for synchrotron energy loss.

3.1 Sensitivity to the cutoff parameter $p_{\text{m}}$

We will now demonstrate that our complete knock-on model satisfies the essential property that the solutions to the kinetic equation are independent of the arbitrary cutoff momentum $p_{m}$ , as long as it is chosen small enough. To determine the sensitivity of the solutions to $p_{\text{m}}$ , we will consider the instantaneous runaway growth rate when the primary runaway population is described by a shifted Maxwellian runaway distribution $f_{\text{RE}}\propto \exp [-(\,\boldsymbol{p}-\boldsymbol{p}_{0})^{2}/q^{2}]$ . For this test we have chosen the momentum $p_{0}\approx 6m_{\text{e}}c$ in the parallel direction, with width $q\approx 0.6m_{\text{e}}c$ . Two electric-field strengths are investigated, a low-field case where $E=3E_{\text{c}}$ and a high-field case $E=100E_{\text{c}}$ . The resulting growth rates are shown in figure 2, as a function of the cutoff $p_{\text{m}}$ after a short time $0.03\unicode[STIX]{x1D70F}_{\text{c}}$ . The growth rate obtained using the field-particle operator alone is nearly independent of $p_{\text{m}}$ as long as it is smaller than $p_{c}$ , indicating that secondary particles created with momentum $p<p_{c}$ are unlikely to run away. For the Rosenbluth–Putvinski operator, this behaviour was also observed by Nilsson et al. (2015).

Figure 2. Runaway growth rate as function of momentum cutoff parameter $p_{\text{m}}$ for two different electric fields, normalized to the field-particle $p_{m}=0$ value. Lines correspond to (dotted blue) the field-particle Boltzmann operator, equation (2.24) and (solid) the full operator including the test-particle operator when (red) $\ln \unicode[STIX]{x1D6EC}$ is held fixed or (black) modified according to (2.12), which is the physically most correct model. Plasma parameters: thermal electron density $n_{\text{e}}=10^{20}~\text{m}^{-3}$ ; temperature $T_{\text{e}}=100~\text{eV}$ .

When the test-particle operator is added, but the Coulomb logarithm $\ln \unicode[STIX]{x1D6EC}$ is left unmodified, the growth rate is decreased. This can be understood from the fact that the test-particle operator represents a source of energy loss for the runaways, which diverges logarithmically as $p_{\text{m}}\rightarrow 0$ . When $\ln \unicode[STIX]{x1D6EC}$ is modified (black line in figure 2, representing the most physically accurate model), the mean energy loss rate of a runaway becomes independent of $p_{\text{m}}$ . The growth rate, however, is found to increase with decreasing $p_{\text{m}}$ , settling to a constant value in the limit $p_{\text{m}}\rightarrow 0$ . The underlying mechanism for this behaviour is that a fraction of all collisions are now accounted for with a Boltzmann operator rather than with a Fokker–Planck operator. This leads to an increase in the runaway probability for particles with $p<p_{c}$ , since the Boltzmann operator fully captures the stochastic nature of the collisions; instead of continuously experiencing the average energy loss, an electron is accelerated freely until it undergoes a collision, by which point it may have gained enough energy to enter the runaway region ( $p>p_{c}$ ). Note that this effect only appears to modify the growth rate with a few per cent, the effect being weaker for smaller electric fields. The effect is, however, directly proportional to $1/\ln \unicode[STIX]{x1D6EC}$ , as it depends on the relative importance of small- and large-angle collisions. This implies that for higher-density or lower-temperature plasmas, the effect can be expected to be more pronounced.

It should be remarked that the field-particle knock-on operator uses a constant Coulomb logarithm in the Fokker–Planck operator, yet is still well behaved when $p_{m}$ becomes small. We have pointed out that the field-particle knock-on operators, like those used in previous runaway avalanche studies, double count collisions with the Fokker–Planck operator. However, they do so only with the field-particle Fokker–Planck operator, and not the test-particle operator which describes the friction on runaways. Therefore, only the Coulomb logarithm in the field-particle operator should be modified when using such models. The field-particle Fokker–Planck operator is essential when considering the dynamics of the bulk population, however it does not significantly affect the avalanche growth rate, thereby explaining the insensitivity to $p_{m}$ for $p_{m}\lesssim p_{c}$ .

3.2 Steady-state avalanche growth rate at moderate electric fields

The steady-state avalanche growth rate in a constant electric field is a classical result; Rosenbluth and Putvinski derived the growth rate formula (1.2) in 1997. After an initial transient, the distribution function tends to approach the asymptotic quasi-steady-state behaviour $f(t,p,\unicode[STIX]{x1D709})\sim n_{\text{RE}}(t)\bar{f}(p,\unicode[STIX]{x1D709})$ , where $\int \bar{f}\,\text{d}\,\boldsymbol{p}=1$ . The kinetic equation, being linear in the runaway distribution, then prescribes that the runaway population will grow with a constant growth rate

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D6E4}=\frac{1}{n_{\text{RE}}}\frac{\text{d}n_{\text{RE}}}{\text{d}(t/\unicode[STIX]{x1D70F}_{c})}.\end{eqnarray}$$

In figure 3 we show the growth rate $\unicode[STIX]{x1D6E4}$ obtained from numerical solutions of the kinetic equation using various models for the knock-on operator, for moderate electric fields ranging from $E=1.5E_{\text{c}}$ to $E=30E_{\text{c}}$ and $Z_{\text{eff}}=1$ . We see that using the Rosenbluth–Putvinski knock-on operator leads to a significant error compared to the more accurate models when the electric field is near the critical – of order 30 % at $1.5E_{\text{c}}$ . At larger electric fields the error is insignificant. Interestingly, the full Boltzmann operator (solid black line) yields a correction of only a few per cent compared to the field-particle operator alone. This means that the test-particle part of the operator does not influence the growth rate significantly. This result is robust; it is not affected by changes in thermal electron density and temperature, and only slightly modified by changes in the effective charge. Note that a significant error is obtained if one fails to account for the double counting of small-angle collisions – the size of this error is sensitive to the cutoff $p_{\text{m}}$ , diverging logarithmically as it approaches zero, and the result is included primarily for illustrative purposes.

Figure 3. Steady-state runaway growth rate normalized to the diffusion-free result $\unicode[STIX]{x1D6E4}_{0}=(E/E_{\text{c}}-1)/2\ln \unicode[STIX]{x1D6EC}$ , equation (1.2), in the presence of a constant electric field, neglecting radiation losses. The red dot-dashed line represents the theoretical prediction, Rosenbluth & Putvinski (1997, equation (18)). Plasma parameters: thermal electron density $n_{\text{e}}=10^{20}~\text{m}^{-3}$ ; temperature $T_{\text{e}}=100~\text{eV}$ ; effective charge $Z_{\text{eff}}=1$ .

3.3 Avalanche generation in a near-threshold electric field with synchrotron radiation losses

In tokamaks, the runaway dynamics in electric fields near the runaway generation threshold is of particular interest. Due to the large self-inductance of tokamaks, after a transient phase during which the ohmic current of the background is dissipated, the electric field will tend towards that value $E_{a}$ – the threshold field – for which the runaway growth rate vanishes, $\unicode[STIX]{x1D6E4}(E_{a})=0$ (Breizman 2014).

At these low electric fields, radiation losses have a large impact, and cannot be ignored in the calculation of the runaway growth rate. In this section, we will include the effect of synchrotron radiation losses and investigate runaway generation when $E\sim E_{a}$ . A model for this was recently presented by Aleynikov & Breizman (2015) (referred to as A&B), using a simplified kinetic equation following a method used by Lehtinen, Bell & Inan (1999). An interesting prediction by the A&B model was that reverse knock-on can have a significant effect on the growth rate, where for electric fields $E\lesssim E_{a}$ , existing runaways will be slowed down to $v<v_{c}$ in single large-angle collision events. This leads to a negative avalanche growth rate, which previous large-angle collision models are incapable of describing, as this process is inherently a large-angle test-particle effect. Using the knock-on operator presented in this work, we will now assess the magnitude of the reverse knock-on effect, as well as determine the threshold field $E_{a}$ and the growth rate when $E\sim E_{a}$ , accounting for radiation losses.

In figure 4 we show how the quasi-steady-state growth rate $\unicode[STIX]{x1D6E4}$ depends on the electric-field strength, similar to figure 4 of A&B. We use the same plasma parameters $Z_{\text{eff}}=5$ and $\unicode[STIX]{x1D70F}_{r}=3m_{\text{e}}n_{\text{e}}\ln \unicode[STIX]{x1D6EC}/(2\unicode[STIX]{x1D700}_{0}B^{2})=70$ (corresponding to $B\approx 1.81~\text{T}$ at $n_{\text{e}}=10^{20}~\text{m}^{-3}$ ), although a slight discrepancy occurs due to our $\ln \unicode[STIX]{x1D6EC}=14.9$ – consistent with the background parameters chosen – compared to their $\ln \unicode[STIX]{x1D6EC}=18$ . In this scenario, the A&B threshold electric field is $E_{a}\approx 1.71E_{c}$ . Several models for the knock-on operator are included in the comparison, in addition to the no-avalanche case since we are now interested in the sub-threshold dynamics. The simulations are run for approximately 300 relativistic collision times $\unicode[STIX]{x1D70F}_{c}$ , upon which the growth rates have settled to the asymptotic steady-state value, corresponding to approximately 6 s with the plasma parameters given above.

Figure 4. Steady-state runaway growth rate in the presence of a constant electric field, accounting for synchrotron radiation losses and using various models for the large-angle collision operator: the Rosenbluth–Putvinski operator, equation (2.2) (green, dashed); the Boltzmann operator equations (2.23)–(2.24) (black, solid) and without any large-angle collision operator (black, dash-dotted). For comparison we have included equation (11) of Aleynikov & Breizman (2015) (blue, dotted). In (b), the avalanche-free growth rate (black, dotted line in (a)) has been subtracted to yield a pure ‘avalanche growth rate’ $\unicode[STIX]{x1D6E4}_{\text{ava}}$ . Plasma parameters: thermal electron density $n_{\text{e}}=10^{20}~\text{m}^{-3}$ ; temperature $T_{\text{e}}=1~\text{keV}$ ; effective charge $Z_{\text{eff}}=5$ , $B=1.81~\text{T}$ .

It is interesting to observe that the test-particle operator, which allows runaways to be thermalized in a single large-angle collision, does not significantly modify the dynamics, in contrast to the theoretical prediction by Aleynikov & Breizman (2015). Unlike the A&B model, which predicts a significant negative growth rate due to this effect when $E\lesssim E_{a}$ , we find that the large-angle collision operator always adds a positive contribution to the total growth rate compared to the no-avalanche case (see figure 4(b) where the no-avalanche growth rate has been subtracted). It can be concluded that the negative growth rates in the sub-threshold regime is a result primarily of the Fokker–Planck dynamics, rather than of large-angle collisions. The reason for this discrepancy to the A&B model can be understood by considering the behaviour of the distribution shape functions $\bar{f}=f/n$ , defined before (3.2), which are illustrated in figures 5 and 6. A&B predicted the distribution to be a delta function in momentum, located at the point of force balance, $p_{\max }$ . When this occurs near the critical speed $v_{c}$ , runaways cannot produce knock-ons with sufficient energy to become runaway. Large-angle collisions then act only to slow down the existing population. In numerical solutions of the full kinetic equation, conversely, it is found that the runaway population takes on a wide energy spectrum, and there will always be sufficiently many runaways with the energy required to produce new runaways to counter the reverse knock-on effect.

Figure 5. Steady-state runaway momentum distributions $\int p^{2}\bar{f}\,\text{d}\unicode[STIX]{x1D6FA}$ (defined to have unit area under the shown curves). Included in the panels is the maximum runaway momentum $p_{\max }$ predicted by Aleynikov & Breizman (2015). Plasma parameters: background electron density $n_{\text{e}}=10^{20}~\text{m}^{-3}$ ; temperature $T_{\text{e}}=1~\text{keV}$ ; effective charge $Z_{\text{eff}}=5$ , $B=1.81~\text{T}$ .

Figure 6. Steady-state normalized runaway momentum distributions $\log _{10}\bar{f}$ from the $E=E_{a}$ case of figure 5.

A notable difference between the A&B model and full solutions of the kinetic equation considered here, which can play an important role when considering the decay of the runaway current in tokamaks, is that the growth rate is not as sensitive to variations in electric field close to (but below) the effective critical field $E_{a}$ as predicted by A&B. It is known that the decay rate is determined primarily by the value of the effective critical field $E_{a}$ when the self-inductance can be considered large (roughly when the total runaway current is much larger than ${\sim}200~\text{kA}$ ) (Breizman 2014). Since the threshold field $E_{a}$ given by the A&B model is reasonably accurate in many cases, it is likely that it may be used to describe runaway current decay in large-current scenarios. However, for moderate runaway currents in the range of hundreds of kA, the overall shape of $\unicode[STIX]{x1D6E4}(E)$ will determine its evolution, which previous theoretical models fail to describe – particularly for electric fields $E\lesssim E_{a}$ .

Finally, we show the effective critical field $E_{a}$ calculated numerically by CODE for a wide range of $Z_{\text{eff}}$ and magnetic-field strength parameters $\unicode[STIX]{x1D70F}_{\text{r}}=6\unicode[STIX]{x03C0}\unicode[STIX]{x1D716}_{0}^{2}m_{e}^{2}c^{3}/e^{4}B^{2}\unicode[STIX]{x1D70F}_{c}$ . This is shown in figure 7, along with the values given by the A&B model by determining the roots of their equation (11). It is seen that the predictions of Aleynikov & Breizman (2015) are typically accurate unless the effective charge is very large, and are most accurate for sufficiently small or large $B$ . The observed trend in the accuracy of their model is unexpected, since they have utilized fast pitch-angle equilibration time (large $Z_{\text{eff}}$ ) and weak magnetic field (large $\unicode[STIX]{x1D70F}_{\text{r}}$ ) in order to reduce the kinetic equation to a tractable form.

Figure 7. Threshold electric field determined numerically from solutions of the kinetic equation, as a function of normalized magnetic-field strength $\unicode[STIX]{x1D70F}_{\text{r}}$ for various values of the effective charge. Predictions by the theoretical model of Aleynikov & Breizman (2015) are included for comparison.

4 Conclusions

Predictions indicate that a major part of the initial plasma current in large tokamaks can be converted to runaway-electron current. This is partly due to the large plasma size limiting the loss of runaway-electron seeds, but more importantly, it is due to the avalanche mechanism which leads to an exponential growth of runaways. The runaway-electron growth rate due to avalanching is exponentially sensitive to the plasma current, and avalanche runaway generation is therefore expected to be a serious issue in ITER and other high-current reactor-scale tokamaks. As the plasma current in present devices cannot be increased above a few megaamperes, full experimental simulation of high-current tokamak disruptions is not possible. Therefore it is very important to develop accurate theoretical models from first principles, to test the validity of approximative models.

In this paper we have developed a fully conservative knock-on collision operator derived from the relativistic Boltzmann operator, and compared it to existing models. Close to the critical electric field, the new model leads to behaviour significantly different from that of the widely used Rosenbluth & Putvinski (1997) avalanche model. This influences the predictions for the transformation of a runaway seed to an avalanching population; fortunately the new operator predicts a lower growth rate than Rosenbluth & Putvinski (1997) and therefore the implications for ITER should be positive, although the difference between models is marginal for high electric fields. We have also described how to resolve the issue of double counting the small- and large-angle collisions, and have illustrated the importance of this issue. The new operator includes both the test-particle and field-particle parts of the collision operator, however we have shown that the test-particle part does not influence the growth rate significantly.

Using kinetic simulations we have performed a careful study of the runaway growth rate in the presence of synchrotron radiation losses and several different avalanche operators. Again, we find a significant difference in runaway rates close to the critical field, however, the effective critical field appears to be well reproduced by simplified models unless the effective charge is very large.