{'fraction_non_alphanumeric': 0.03776210799777324, 'fraction_numerical': 0.014311560586379663, 'mean_word_length': 4.380381879445901, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 14, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On the relativistic large-angle electron collision operator for runaway avalanches in plasmas

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.


Introduction
Large-angle collisions are associated with large momentum transfers, however 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 Λ (Rosenbluth et al. 1957;Trubnikov 1965). In many plasmas, e.g. magnetic fusion plasmas and astrophysical plasmas, ln Λ is typically of order 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 Λ, as they cause an exponential growth of the runaway density -a runaway avalanche (Jayakumar et al. 1993). Runaway is the acceleration of particles in the presence of an electric field which exceeds the critical field E c = n e ln Λe 3 (4πε 2 0 m e c 2 ) −1 (Connor and Hastie 1975), where n e is the electron density, e is the elementary charge, ε 0 is the vacuum permittivity, m e is the electron rest mass, and c is the speed of light. Since the collisional drag for superthermal electrons is given by F c = eE c (v/c) −2 , any electrons with speed greater than the critical speed v c = c √ E c /E will be accelerated indefinitely, and are hence referred to as runaway electrons (Wilson 1925). Runaway electrons occur in a wide range of plasmas, e.g. in atmospheric discharges (Dwyer 2007), 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 and Putvinski 1997). The subsequent uncontrolled loss of a runaway-electron beam could damage plasmafacing 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 and Hastie 1975;Cohen 1976) and is where the undetermined constant factor κ is of order unity. Here n RE is the number density of runaways, τ c = 4πε 2 0 m 2 e c 3 /n e e 4 ln Λ is the relativistic-electron collision time, E D = m 2 e c 3 /(eτ c T e ) is the so-called Dreicer field, and Z eff = i n i Z 2 i / i n i Z i is the effective ion charge (with the sum taken over all ion species i). The parameters h, λ and η (not given here) depend on E/E c and approach unity as E/E c becomes large (in the non-relativistic limit), but ensure that the growth rate vanishes as E → E c .
A secondary runaway generation mechanism is provided by large-angle collisions, whereby an electron with kinetic energy = (γ−1)m e c 2 > 2 c can send a stationary target electron into the runaway region in a single collision event while remaining a runaway itself, where 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 and Putvinski (1997) to be approximately dn RE dt ava ≈ C n RE 2 ln Λτ c E E c − 1 . (1.1) The function C = C(E, Z eff ) was shown to be C = 1 when collisional diffusion is neglected (formally by setting Z eff = −1). While the avalanche growth rate is formally of order 1/ ln Λ smaller than the primary generation rate, the more favorable scaling with electric field makes it the dominant source of new runaways for sufficiently large runaway populations n RE or sufficiently small E/E D , i.e. at sufficiently 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 ava ≈ 2 ln Λm e c/[Ce(E − E c )]. This corresponds to the time when the fastest runaway has been accelerated to a kinetic energy E k ≈ ln Λ/C MeV (Jayakumar et al. 1993) (neglecting the weak electric-field dependence of C). Numerically, t ava ≈ 3.4 ln Λ/[C(E − E c )] ms, with E and E c in V/m. 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 ava (Gurevich et al. 1994;Gurevich and Zybin 2001;Helander et al. 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 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 and Breizman (2015) to be significant for runaway in a near-threshold electric field.
The rest of the paper is organized as follows. In Section 2 we introduce the theoretical mod-els describing the large-angle collisions. After giving an overview of the existing models, we present a derivation of the new conservative operator. In Section 3 we investigate the effect of the new operator on the runaway growth rate numerically, using the kinetic-equation solver CODE (Landreman et al. 2014;Stahl et al. 2016). Finally, we summarize our conclusions in Section 4.

Theoretical models for runaway generation due to large-angle collisions
One of the earliest models for avalanche runaway generation was introduced by Rosenbluth and 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 and 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 where d f e /dt represents the advective part of the motion, C 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 and 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) where ξ = cos θ = p /p is the pitch-angle cosine, γ = 1 + (p/m e c) 2 is the Lorentz factor, ξ 0 = (γ − 1)/(γ + 1) and the momentum-space volume element is p 2 dpdξdϕ, with ϕ the azimuthal angle of the momentum (the gyroangle). The delta function ensures that secondary electrons are only born on the parabola p 2 ⊥ = 2p m e c in momentum space. In the non-relativistic limit, p m e c, secondaries are born at perpendicular angles, p ∼ 0, and are prone to trapping in an inhomogeneous magnetic field. Away from the magnetic axis of a tokamak, this can lead to a strong reduction in the avalanche growth rate, as recently shown 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 where Σ = (2πr 2 0 ) −1 dσ/dγ is the normalized Møller differential cross-section for free-free electronelectron scattering (Møller 1932), r 0 = e 2 /(4πε 0 m e c 2 ) is the classical electron radius, and γ 1 is connected to p and ξ by the relation where a misprint in the original paper incorrectly replaced the γ − 1 factor with γ 1 . Since the authors work under the assumption that the runaway pitch-angles are negligible, the distribution only appears in the angle-averaged form F(p 1 , t) = dξ 1 dϕ 1 p 2 1 f e (p 1 , ξ 1 , t). (2.6) 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 RP and S CH ) increase rapidly with decreasing momenta and the sources are thus sensitive to the choice of cut-off momentum (introduced to avoid double-counting 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.

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 (Cercignani and Kremer 2002;Montgomery and Tidman 1964). The collision operator can be defined as C ab ( f a ) = (dn a ) c,ab /dtdp, where (dn a ) 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 dσ by (Cercignani and Kremer 2002;Lifshitz and Pitaevski 1981) (dn a ) c,ab = f a (p 1 ) f b ( p 2 )ḡ ø dσ ab d p 1 d p 2 dt − f a ( p) f b ( p )g ø dσ ab dp d p dt. (2.7) The first term on the right-hand side, the gain term, describes the rate at which particles a of momentum p 1 will scatter to momentum p. The second term, the loss term, is the rate at which particles a scatter away from momentum p. Here, we introduced the Møller relative speed g ø = (v − v ) 2 − (v × v ) 2 /c 2 and the differential cross-section dσ ab for scattering events p, p → p 1 , p 2 . The barred quantities are defined likewise, but with p exchanged for p 1 and p for 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ḡ ø dσ ab dp 1 d p 2 = g ø dσ ab d pd p (known as the principle of detailed balance (Weinberg 2005)), allowing the collision operator to be cast in the commonly adopted symmetric form where p 1 and p 2 (six degrees of freedom) are uniquely determined in terms of p and p by two scattering angles and four constraints by the conservation of momentum and energy, (2.10) 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 ∆p = p 1 − 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, where the maximum center-of-mass deflection angle for self collisions is θ max = π/2 (not π as it is for unlike-species collisions, or collisions would be double counted) and θ min is a cut-off 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 FP + C Boltz as prescribed by Eq. (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 θ m , and the knock-on (Boltzmann) operator for θ > θ m . The Coulomb logarithm used in the Fokker-Planck operator then ought to be changed from Eq. (2.11) to ln Λ = ln Λ − ln cot θ m 2 . (2.12) When an incident electron of momentum p knocks a stationary electron to momentum p m , the corresponding center-of-mass scattering angle θ m is given by (2.13) 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 cut-off p m , when p m m e c. The number of collisions that are double-counted can often be significant when this effect is unaccounted for. Assuming v m /c ∼ √ E c /E to be located at a non-relativistic energy (that is, we assume E E c ), the modification to the Coulomb logarithm is approximately given by ln 2(E/E c )(γ − 1). For highly energetic electrons with γ ∼ 50 and E/E c ∼ 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 θ m approaches the cut-off 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 ln Λ = 0. However, this corresponds to a cut-off momentum smaller than thermal, p m p 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 m , and it is therefore desirable to choose p 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 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 Boltz- † In the quantum-mechanical treatment, it is rather the de-Broglie wavelength of the center-of-mass momentum transfer λ = /| p * 1 − p * | that cuts off at the Debye length. mann operator given by Eq. (2.8), but instead use the alternative given directly from Eq. (2.7), where σ ab ( p, p ) = dp 1 ∂σ ab /∂ p 1 is the total cross-section.

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 e ( p) = f Me ( p) + δ f e (p), where the runaway distribution δ f e is much smaller than the bulk distribution f Me , ||δ f e || || f Me || (for example in terms of number densities n RE n e ). We may then linearize the bilinear Boltzmann operator by ignoring terms quadratic in δ f e , obtaining 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 particles. The second term, the field-particle term, describes the reaction of the bulk 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 Te = √ 2T e /m e , so that we may approximate the bulk population with a Dirac delta function: f Me ( p) ≈ n e δ( p). The collision operator then takes the form The total cross-section σ ee ( p) is given in Eq. (A 9) in Appendix A. The limiting momenta q * 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 γ 1 = γ + γ − 1, where γ and γ are the Lorentz factors of the two electrons after the collision. The conditions γ > γ or γ < γ determines whether γ 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 done in order to avoid double counting). We therefore obtain q * from setting γ = γ in the conservation law, giving γ * = 2γ − 1, which corresponds to q * = m e c γ * 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 m p Te . Setting γ = γ m then yields the lower limit γ 0 = γ + γ m − 1, corresponding to q 0 = m e c γ 2 0 − 1. Note that for the total operator C Boltz , the two gain terms in Eqs. (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 * 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 ∂σ ee /∂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 θ s is the deflection angle, Σ is defined in Eq. (2.4), and ξ * is defined in Eq. 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 and Putvinski (1997) and Chiu et al. (1998) -Eqs. (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 and 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 and 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 and Pankratov (1986).
There are multiple ways of carrying out the integration over the delta function; if we assume a distribution function independent of gyro-angle, f e ( p 1 ) = f e (p 1 , cos θ 1 ), a few convenient expressions are given by where we have introduced the quantities dξ f e (p, ξ)P L (ξ).
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 Boltz ( p) = hibits 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 ϕ s were also found by Gurevich and Zybin (2001) for the so-called ionization integral in neutral gases. The form of the integral taken over ξ 1 was obtained by Helander et al. (1993) in the analogous problem of elastic nucleon-nucleon scattering, and an equivalent formulation was also recently given by Boozer (2015).
The Legendre modes of the collision operator are explicitly given by (2.23) Note further that since we only consider those collisions where both the incident and outgoing particles have momenta p > p m , the gain terms must only be applied for γ > γ m , while the test-particle loss term is applied for γ > 2γ m − 1. In Appendix A it is explicitly demonstrated that this collision operator conserves density, momentum and energy. A qualitative illustration of the large-angle collision operators discussed here is shown in Figure 1. A test runaway distribution (Fig. 1a) was generated by applying a constant electric field E = 15E c for a short time t ≈ 0.5τ c with Z 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 Fig. 1c and Fig. 1d 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 1b 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 Fig. 1b and 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 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 Te v c limits the electric-field values to √ E √ E D /2 (effectively forming a lower limit in density and an upper limit in temperature for a given E). This limitation would be resolved by accounting for the velocity distribution of the target population in Eq. (2.7), resulting in a significantly more complicated operator. However, in most scenarios this is not an issue; the critical velocity tends to be significantly larger than thermal. If this was not the case, the entire electron population would run away within a few collision times by the primary generation mechanism regardless.

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 where F L is the Lorentz force, and F S is the radiation reaction force associated with synchrotron radiation and the brackets denote averaging over the azimuthal (gyro) angle. C ei and C ee are the gyroaveraged Fokker-Planck collision operators for electron-ion and electron-electron collisions, respectively. First, we will study the sensitivity of the avalanche dynamics to the arbitrary cut-off parameter p 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.

Sensitivity to the cut-off parameter p 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 cut-off momentum p m , as long as it is chosen small enough. To determine the sensitivity of the solutions to p m , we will consider the instantaneous runaway growth rate when the primary runaway population is described by a shifted Maxwellian runaway distribution f RE ∝ exp − ( p − p 0 ) 2 /q 2 . For this test we have chosen the momentum p 0 ≈ 6 m e c in the parallel direction, with width q ≈ 0.6m e c. Two electric-field strengths are investigated, a low-field case where E = 3E c and a high-field case E = 100E c . The resulting growth rates are shown in Fig. 2, as a function of the cut-off p m after a short time 0.03τ c . The growth rate obtained using the field-particle operator alone is nearly independent of p m as long as it is smaller than p c , indicating that secondary particles created  Figure 2. Runaway growth rate as function of momentum cut-off parameter p 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, Eq. (2.23), and (solid) the full operator including the test-particle operator when (red) ln Λ is held fixed or (black) modified according to Eq. (2.12), which is the physically most correct model. Plasma parameters: thermal electron density n e = 10 20 m −3 ; temperature T e = 100 eV.
with momentum p < p c are unlikely to run away. For the Rosenbluth-Putvinski operator, this behavior was also observed by Nilsson et al. (2015). When the test-particle operator is added, but the Coulomb logarithm ln Λ 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 m → 0. When ln Λ is modified (black line in Fig. 2, representing the most physically accurate model), the mean energy-loss rate of a runaway becomes independent of p m . The growth rate, however, is found to increase with decreasing p m , settling to a constant value in the limit p m → 0. The underlying mechanism for this behavior 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 percent, the effect being weaker for smaller electric fields. The effect is, however, directly proportional to 1/ ln Λ, as it depends on the relative importance of small-angle 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 runawayavalanche 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 p c .

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.1) in 1997. After an initial transient, the distribution function tends to approach the asymptotic quasi-steady state behavior f (t, p, ξ) ∼ n RE (t)f (p, ξ), where f dp = 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 Γ 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 c to E = 30E c and Z 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 c . At larger electric fields the error is insignificant. Interestingly, the full Boltzmann operator (solid black line) yields a correction of only a few percent 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 cut-off p m , diverging logarithmically as it approaches zero, and the result is included primarily for illustrative purposes.

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, Γ(E a ) = 0 (Breizman 2014). At these low electric fields, radiation losses have a large impact, and can not be ignored in the calculation of the runaway growth rate. In this section, we will include the effect of synchrotron-  (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" Γ ava . Plasma parameters: thermal electron density n e = 10 20 m −3 ; temperature T e = 1 keV; effective charge Z eff = 5, B = 1.81 T. radiation losses and investigate runaway generation when E ∼ E a . A model for this was recently presented by Aleynikov and Breizman (2015) (referred to as A&B), using a simplified kinetic equation following a method used by Lehtinen et al. (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 E a , existing runaways will be slowed down to v < v c in single largeangle collision events. This leads to a negative avalanche-growth rate, which previous largeangle 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 ∼ E a , accounting for radiation losses.
In Fig. 4 we show how the quasi-steady state growth rate Γ depends on the electric-field strength, similar to Fig. 4 of A&B. We use the same plasma parameters Z eff = 5 andτ rad = 3m e n e ln Λ/(2ε 0 B 2 ) = 70 (corresponding to B ≈ 1.81 T at n e = 10 20 m −3 ), although a slight discrepancy occurs due to our ln Λ = 14.9 -consistent with the background parameters chosen -compared to their ln Λ = 18. In this scenario, the A&B threshold electric field is E a ≈ 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.
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 and Breizman 2015). Unlike the A&B model, which predicts a significant negative growth rate due to this effect when E E a , we find that the largeangle collision operator always adds a positive contribution to the total growth rate compared to the no-avalanche case (see Fig. 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 behavior of the distribution shape functionsf = f /n, defined before Eq. (3.2), which are illustrated in Fig. 5. 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 popu- lation. 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.
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 ∼ 200 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 Γ(E) will determine its evolution, which previous theoretical models fail to describeparticularly for electric fields E E a .
Finally, we show the effective critical field E a calculated numerically by CODE for a wide range of Z eff and magnetic-field strength parameters τ r = 6π 2 0 m 2 e c 3 /e 4 B 2 τ c . This is shown in Fig. 6, along with the values given by the A&B model by determining the roots of their Eq. (11). It is seen that the predictions of Aleynikov and 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 pitchangle equilibration time (large Z eff ) and weak magnetic field (large τ r ) in order to reduce the kinetic equation to a tractable form.

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 behavior significantly different from that of the widely used Rosenbluth and 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 and 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-angle 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.
The full electron-electron Boltzmann operator C is known to satisfy conservation of density, momentum and energy, expressed by the relations d pC( f e ) = 0, dp pC( f e ) = 0, (A 1) dp m e c 2 (γ − 1)C( f e ) = 0, or in our case of a cylindrically symmetric plasma, in terms of the Legendre modes of the collision operator, dp p 2 C 0 (p) = 0, dp p 3 C 1 (p) = 0, (A 2) dp p 2 (γ − 1)C 0 (p) = 0.
We will show that our explicit form of the knock-on operator, accounting only for collisions involving electrons with momenta p > p m , satisfy the same conservation laws. Taking the full operator Eqs. (2.22)-(2.23), noting that the gain term only applies for γ > γ m and the loss term for γ > 2γ m − 1, we find upon integration (changing momentum integrals to energy integrals by vdp = m e c 2 dγ) In the first term, the integration order can be interchanged by using The integrals over γ 1 will vanish for all cross-sections that respect the indistinguishability of the electrons, i.e. for which Σ(γ 1 , γ) = Σ(γ 2 , γ) where γ 2 = γ + 1 − γ 1 . This follows directly from the observation that confirming that our operator indeed satisfies the conservation laws.

A.1. Total cross-section
For our case of the Møller cross-section, the differential cross-section of Eq. 2.4 can be integrated analytically to produce the total cross-section. One obtains

(B 2)
We also utilized cos θ s = cos θ when cos θ 1 = 1, and kinematics constrain the incident momentum p 1 according to Eq. (2.5). This result agrees exactly with the Chiu-Harvey source S CH of Eq. (2.3). In terms of an expansion in Legendre polynomials, the Chiu-Harvey avalanche source is obtained from the general field-particle operator in Eq. (2.23) simply by replacing f L (p) by (2L + 1) f 0 (p), corresponding to the delta-function approximation. In this representation, however, the approximation holds limited appeal as it does not provide a significant simplification of the collision operator; indeed, compared to the full operator it requires a larger number of Legendre polynomials to be retained since the true f L decreases rapidly with L for sufficiently large L. By the addition of the sink terms in (2.22) and (2.23), and extending the integration limit down from q * to q 0 , the Chiu-Harvey operator can be made conservative. However, the delta-function assumption in pitch-angle causes incorrect momentum dynamics, and the total momentum of the distribution will not be conserved in this treatment. This can be corrected by treating the L = 1 mode exactly, corresponding to a total conservative knock-on operator in the Chiu-Harvey approximation of the form C (cons) CH =S CH − 1 4τ c ln Λ v f e ( p)σ(p) − δ( p) 4τ c ln Λ p >q 0 (p m ) dp v f e ( p )σ(p ) − 3(m e c) −3 8πτ c ln Λ ξ γ(γ + 1) p 1 >q 0 d p 1 γ 1 + 1 γ 1 Σ(γ, γ 1 )(1 − ξ 1 ) f e ( p 1 ) + 3 8τ c ln Λ ξvσ(p)