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 Reference Montgomery and Tidman1964; Cercignani & Kremer Reference Cercignani and Kremer2002). 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 (Reference Lifshitz and Pitaevski1981) and Cercignani & Kremer (Reference Cercignani and Kremer2002)

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 Reference Weinberg2005)), allowing the collision operator to be cast in the commonly adopted symmetric form

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,

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 Reference Landau1936; Akama Reference Akama1970), 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,

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,Footnote ¹ 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

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

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),

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

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

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

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

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 (Reference Rosenbluth and Putvinski1997) and Chiu et al. (Reference Chiu, Rosenbluth, Harvey and Chan1998) – (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. (Reference Chiu, Rosenbluth, Harvey and Chan1998) incorrectly ascribe the discrepancy between their result and that of Besedin & Pankratov (Reference Besedin and Pankratov1986) 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 Reference Besedin and Pankratov1986) 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 (Reference Besedin and Pankratov1986).

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

where we have introduced the quantities

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 (Reference Gurevich and Zybin2001) 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 (Reference Helander, Lisak and Ryutov1993) in the analogous problem of elastic nucleon–nucleon scattering, and equivalent formulations were also recently given by Aleynikov et al. (Reference Aleynikov, Aleynikova, Breizman, Huijsmans, Konovalov, Putvinski and Zhogolev2014) and Boozer (Reference Boozer2015).

The Legendre modes of the collision operator are explicitly given by

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.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. (Reference Nilsson, Decker, Peysson, Granetz, Saint-Laurent and Vlainic2015).

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

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 (Reference Rosenbluth and Putvinski1997, 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 Reference Breizman2014).

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 (Reference Aleynikov and Breizman2015) (referred to as A&B), using a simplified kinetic equation following a method used by Lehtinen, Bell & Inan (Reference Lehtinen, Bell and Inan1999). 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 (Reference Aleynikov and Breizman2015) (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 (Reference Aleynikov and Breizman2015). 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 (Reference Aleynikov and Breizman2015). 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 Reference Breizman2014). 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 (Reference Aleynikov and Breizman2015) 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 (Reference Aleynikov and Breizman2015) are included for comparison.