{'fraction_non_alphanumeric': 0.032594185701572236, 'fraction_numerical': 0.010605161673094037, 'mean_word_length': 4.570956413964057, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical calculation of the runaway electron distribution function and associated synchrotron emission

Synchrotron emission from runaway electrons may be used to diagnose plasma conditions during a tokamak disruption, but solving this inverse problem requires rapid simulation of the electron distribution function and associated synchrotron emission as a function of plasma parameters. Here we detail a framework for this forward calculation, beginning with an efficient numerical method for solving the Fokker-Planck equation in the presence of an electric field of arbitrary strength. The approach is continuum (Eulerian), and we employ a relativistic collision operator, valid for arbitrary energies. Both primary and secondary runaway electron generation are included. For cases in which primary generation dominates, a time-independent formulation of the problem is described, requiring only the solution of a single sparse linear system. In the limit of dominant secondary generation, we present the first numerical verification of an analytic model for the distribution function. The numerical electron distribution function in the presence of both primary and secondary generation is then used for calculating the synchrotron emission spectrum of the runaways. It is found that the average synchrotron spectra emitted from realistic distribution functions are not well approximated by the emission of a single electron at the maximum energy.

ions [3] and positrons [4,5] can also be produced. Relatively large electric fields are required for runaway production, and in tokamaks these can arise during disruptions or in sawtooth events. Understanding of runaway electrons and their generation and mitigation is essential to planning future large experiments such as ITER. For this purpose it is important to know the electron distribution function accurately.
In this work, we demonstrate a new numerical tool named CODE (COllisional Distribution of Electrons) that can be used for computing the electron distribution function in general and the production rate of runaway electrons in particular. These quantities are determined by a balance between acceleration in the electric field and collisions with both electrons and ions. The calculation is fully relativistic, using a collision operator valid for both low and high velocities [6] and it includes both primary and secondary runaway electron generation. If only the primary runaway electron generation is of interest, CODE can be used in both time-dependent and time-independent modes. The latter mode of operation, in which a long-time quasi-equilibrium distribution function is calculated, is extremely fast in that it is necessary only to solve a single sparse linear system. Thus, the time-independent CODE is highly suitable for coupling within larger more expensive simulations. If secondary (avalanche) generation is also considered, CODE must be run in the time-dependent mode.
Several simplifying assumptions are made in CODE. The most restrictive assumption is that spatial inhomogeneity and the associated magnetic trapping are neglected. This approximation may be justified near the magnetic axis, where most runaways are produced, since magnetic trapping becomes negligible there. The neglect of spatial inhomogeneity also means that radial transport of runaways is not included. To include these effects of spatial nonuniformity, a code such as LUKE [7,8] or CQL3D [9,10] must be used. Furthermore, we neglect the effect of radiation (synchrotron and bremsstrahlung) on the distribution function.
We anticipate the distribution functions computed in CODE can be used as inputs to other calculations. In this paper, we present a first example of this type, a computation of the synchrotron radiation spectrum of a runaway beam. Another possible application could be a study of instabilities driven by the anisotropy of the electron distribution function. We leave this and other possible applications for future work.
The remainder of the paper is organized as follows. In Sec. II we present the kinetic equation and the collision operator used. Section III details the discretization scheme and calculation of the primary runaway production rate, with typical results shown in section IV. The avalanche source term and its implementation is described in Sec. V and examples of applications are given in Sec. VI. We conclude in Sec. VII.

II. KINETIC EQUATION AND NORMALIZATIONS
We begin with the kinetic equation Here, −e is the electron charge, E = |E| is the magnitude of the electric field E, e = −E/E is a unit vector antiparallel to the electric field, ∇ p is the gradient in the space of relativistic c is the speed of light, C is the electron collision operator, and S represents any sources. All quantities refer to electrons unless noted otherwise. Spatial inhomogeneity is ignored. We may write e · ∇ p f in (1) in terms of scalar variables using where p = |p|, and ξ = p · e/p is the cosine of the pitch angle relative to the electric field.
The distribution function is defined such that the density n is given by n = d 3 p f , so f has dimensions of (length × momentum) −3 , and we assume the distribution function for small momentum to be approximately the Maxwellian f M = nπ −3/2 (mv e ) −3 exp(−y 2 ) where v e = 2T /m is the thermal speed, and y = p/(mv e ) = γv/v e is the normalized momentum.
We use the collision operator from Appendix B of Ref. [6]. This operator is constructed to match the usual nonrelativistic test-particle operator in the limit of v ≪ c, and in the relativistic limit it reduces to the operator from Appendix A of Ref. [11]. The collision operator is where is identical to the Γ defined in Refs. [6,11,12], ν ee = 4 √ 2πe 4 n ln Λ/(3 √ mT 3/2 ) is the usual Braginskii electron collision frequency, φ(x) = 2π −1/2 x 0 exp(−s 2 ) ds is the error function, and is the Chandrasekhar function. In the nonrelativistic limit δ → 0, then y → x, and (3) reduces to the usual Fokker-Planck test-particle electron collision operator.
The collision operator (3) is approximate in several ways. First, it originates from the Fokker-Planck approximation in which small-angle collisions dominate, which is related to an expansion in ln Λ ≫ 1. Consequently, the infrequent collisions with large momentum exchange are ignored, so the secondary avalanche process is not included at this stage, but will be addressed later in Sec. V. Also, the modifications to the Rosenbluth potentials associated with the high-energy electrons are neglected, i.e. collisions with high-energy field particles are ignored.
The kinetic equation is normalized by multiplying through with m 3 v 3 e π 3/2 /(ν ee n), and defining the normalized distribution function so that F → 1 at p → 0. We also introduce a normalized electric field which, up to a factor of order unity, is E normalized by the Dreicer field. The normalized time ist = ν ee t and the normalized source isŜ = Sm 3 v 3 e π 3/2 /(ν ee n). We thereby obtain the dimensionless equation Notice that this equation has the form of a linear inhomogeneous 3D partial differential equation: for a linear time-independent differential operator M. If a time-independent equilibrium solution exists, it will be given by F = M −1Ŝ .
Since both the electric field acceleration term and the collision operator in the kinetic equation (1) have the form of a divergence of a flux in velocity space, the total number of particles is constant in time in the absence of a source: (d/dt) d 3 p f = d 3 p S. However, runaway electrons are constantly gaining energy, so without a source at small p and a sink at large p, no time-independent distribution function will exist. From another perspective, a nonzero source is necessary to find a nonzero equilibrium solution of (11), because when To find a solution, we must either consider a time-dependent problem or include a nonzero S. In reality, spatial transport can give rise to both sources and sinks, and a sink exists at high energy due to radiation. When included, secondary runaway generation (considered in Sec. V) also introduces a source. To avoid the added complexity of these sinks and sources and simultaneously avoid the intricacies of time dependence, when restricting ourselves to primary generation we may formulate a time-independent problem as follows. We takê S = αe −y 2 for some constant α, representing a thermal source of particles. Equation (11) for ∂/∂t = 0 may be divided through by α and solved for the unknown F/α. Then α may be determined by the requirement F (p = 0) = 1, and F is then obtained by multiplying the solution F/α by this α.
The constant α represents the rate at which particles must be replenished at low energy to balance their flux in velocity space to high energy. Therefore, α is the rate of runaway production. As we do not introduce a sink at high energies, F will have a divergent integral over velocity space.
CODE can also be run in time-dependent mode. Once the velocity space coordinates and the operator M are discretized, any implicit or explicit scheme for advancing a system

III. DISCRETIZATION
We first expand F in Legendre polynomials P L (ξ): Then the operation 2L + 1 2 is applied to the kinetic equation. Using the identities in the appendix, we obtain . Note that the collision operator is diagonal in the L index, and the electric field acceleration term is tridiagonal in L.
It is useful to examine the L = 0 case of (15), which corresponds to (half) the integral of the kinetic equation over ξ: Applying 4π −1/2 ∞ y b dy y 2 ( · ) for some boundary value y b , and assuming the source is negligible in this region, we obtain where n r is the number of runaways, meaning the number of electrons with y > y b , so that −1 dξ f . The runaway rate calculated from (17) should be independent of y b in steady state (as long as y b is in a region ofŜ 0 = 0), which can be seen by applying y b2 y b1 dy y 2 ( · ) to (16). We find in practice it is far better to compute the runaway production rate using (17) than from the source magnitude α, since the latter is more sensitive to the various numerical resolution parameters.
To discretize the equation in y, we can apply fourth-order finite differences on a uniform grid. Alternatively, a coordinate transformation can be applied so grid points are spaced further apart at high energies. The y coordinate is cut off at some finite maximum value y max . The appropriate boundary conditions at y = 0 are dF 0 /dy = 0 and F L = 0 for L > 0.
For the boundary at large y, we impose F L = 0 for all L. This boundary condition creates some unphysical grid-scale oscillation at large y, which may be eliminated by adding an  figure 1.a, excellent convergence is achieved for the parameters used here. As expected, the distribution function is increased in the direction opposite to the electric field (p || > 0). While the distribution function is reduced in the direction parallel to the electric field (p || < 0) for y < 5, F is actually slightly increased for y > 5 due to pitch-angle scattering of the high-energy tail electrons [13]. This interpretation can be confirmed by suppressing the pitch-angle scattering term in CODE, in which case F is reduced in the direction parallel to the electric field for all y. For comparison with previously published results, we show in Figure 3 results by Kulsrud et al [14], who considered only the nonrelativistic case δ → 0. The agreement with CODE is exceptional. The runaway production rate in CODE is computed using (17) for y b = 10. (Any value of y b > 5 gives indistinguishable results.) Ref. [14] uses a different normalized electric field E K which is related toÊ by E K = 2(3 √ π) −1Ê , and in Ref. [14] the runaway rate is also normalized by a different collision frequency ν K = 3 π/2 ν ee . It should also be noted that the Kulsrud Table 1 of Ref. [14].
the density of thermal electrons which may be accelerated in a close collision. The production of runaways through close collisions is known as secondary production, or as avalanche production since it may occur with exponential growth. To simulate secondary generation of energetic electrons, we use a source term derived in [15], starting from the Møller scattering cross-section in the w ≫ 1 limit, with w = p/(mc) = δy a normalized momentum. In this limit, the trajectories of the primary electrons are not much deflected by the collisions. The source then takes the form where 1/τ = 4πne 4 ln Λ/(m 2 c 3 ) is the collision frequency for relativistic electrons, n r is the density of the fast electrons and ξ 2 = w/(1 + √ 1 + w 2 ) is the cosine of the pitch angle at which the runaway is born. (Our (18) differs by a factor m 3 c 3 compared to the source in Ref.
[15] since we normalize our distribution function as n = d 3 p f instead of n = d 3 w f .
There is also a factor of 2π difference due to the different normalization of the distribution function.) Due to the approximations used to derive S, care must be taken in several regards.
First, to define n r in (18), it is not clear where to draw the dividing line in velocity space between runaways and non-runaways. One possible strategy for defining n r is to compute the separatrix in velocity space between trajectories that will have bounded and unbounded energy in the absence of diffusion, and to define the runaway density as the integral of f over the latter region [16]. This approach may somewhat overestimate the true avalanche rate, since it neglects the fact that some time must elapse between an electron entering the runaway region and the electron gaining sufficient energy to cause secondary generation.
As most runaways have ξ ≈ 1, we may approximate the separatrix by setting dw/dt = 0 where dw/dt = eE/(mc) − (1 + 1/w 2 )/τ defines the trajectory of a particle with ξ = 1, neglecting diffusion in momentum and pitch angle. The runaway region is therefore w > w c where w c = [(E/E c ) − 1] −1/2 and E c = eτ /(mc) is the critical field, and so we take n r = 2πm 3 c 3 1 −1 dξ ∞ wc dw w 2 f . (We cannot define n r by the time integral of (17), since (17) is no longer valid when S is nonzero away from p ≈ 0.) A second deficiency of (18) is that S is singular at w → 0, so the source must be cut off below some threshhold momentum. Following Ref. [10], we choose the cutoff to be w c . Neither of the cutoffs discussed here would be necessary if a less approximate source term than (18) were used, but derivation of such an operator is beyond the scope of this paper.
Normalizing and applying (14) as we did previously for the other terms in the kinetic equation, the source included in CODE becomeŝ When secondary generation is included, CODE must be run in time-dependent mode.
To benchmark the numerical solution of the kinetic equation including the above source term by CODE, we use the approximate analytical expression for the avalanche distribution function derived in Section II of Ref. [17]: where k is a constant. The quantityγ is the growth rateγ = (1/f )∂f /∂t, which must be independent of both time and velocity for (20) to be valid. Equation (20) is also valid only where p || ≫ p ⊥ and in regions of momentum space where S is negligible. (This restriction is not a major one since S = 0 everywhere except on the ξ = ξ 2 curve.) If most of the runaway distribution function is accurately described by (20), then we may approximate where n r e −γt is constant. (Equation (21) may be inaccurate in some situations even if (20) is accurate in part of velocity-space, because (21) requires (20) to apply in all of velocity-

space.) Figures 4 and 5 show comparisons between distributions from CODE and (20)-(21)
for two different sets of parameters. More precisely, the quantity plotted in figures 4-5 is log 10 (m 3 c 3 f /n r ). To generate the figures, CODE is run for a sufficiently long time that (1/f )∂f /∂t becomes approximately constant. The resulting numerical value of (1/n r )dn r /dt is then used asγ when evaluating (20)- (21). For a cleaner comparison between CODE and analytic theory in these figures, we minimize primary generation in CODE in these runs by initializing f to 0 instead of to a Maxwellian. For both sets of physical parameters, the agreement between CODE and (20) is excellent in the region where agreement is expected: where p || ≫ p ⊥ and away from the curve ξ = ξ 2 .

VI. APPLICATIONS
As an application of the distribution functions calculated by CODE, we consider the synchrotron radiation emitted by runaway electrons in a magnetized plasma as a result of their gyromotion. Synchrotron emission is a useful diagnostic for runaways, since due to the energy dependence of the emitted synchrotron power, the emission from runaways completely dominates that of the thermal particles. The emission also depends strongly on the pitch-angle of the particle. In a cylindrical plasma geometry, the emitted synchrotron power per wavelength at wavelength λ from a single highly energetic particle is given by [18] P(γ, γ , λ) = 4π √ 3 where the two-dimensional momentum of the particle is determined by γ and is a modified Bessel function of the second kind, and where B is the magnetic field strength.
Using CODE we will demonstrate that the synchrotron radiation spectrum from the entire runaway distribution is substantially different from the spectrum obtained from a single particle approximation. By transforming to the more suitable coordinates w and ξ, related to γ and γ through γ 2 = 1+w 2 and γ 2 = (1−w 2 ξ 2 /(1+w 2 )) −1 , and integrating (22) over the runaway region R in momentum space, we obtain the total synchrotron emission from the runaway distribution. Normalizing to n r , we find that the average emitted power per runaway particle at a wavelength λ is given by The per-particle synchrotron spectra generated by the CODE distributions in Figures   4 and 5 were calculated using this formula, and are shown in Figure 6, together with the spectra radiated by electron distributions for other electric field strengths. For the physical parameters used, we note that the peak emission occurs between 10 and 20 µm. The synchrotron spectra show a decrease in per-particle emission with increasing electric field strength. Even though a stronger electric field leads to more particles with high energy (and thus high average emission), it also leads to a more narrow distribution in pitch-angle. This reduction in the number of particles with large pitch-angle leads to a decrease in average emission. Both figures confirm that the average emission is reduced for higher electric fields, implying that the latter mechanism has the largest impact on the spectrum.
In calculating the spectra, the runaway region of momentum space, R, was defined such that the maximum particle momentum was w max = 50 (which translates to y max = 500 and y max = 1000 respectively for the cases shown in Figures 4 and 5), corresponding to a maximum particle energy of ≃ 25 MeV. Physically the cutoff at large energy can be motivated by the finite life-time of the accelerating electric field and the influence of loss mechanisms such as radiation. Since the radiated synchrotron power increases with both particle energy and pitch, this truncation of the distribution is necessary to avoid infinite emission, although the precise value for the cutoff depends on the tokamak and on dischargespecific limitations to the maximum runaway energy. For the lower bound of R, w min = w c = [(E/E c ) − 1] −1/2 was used, and all particles with ξ ∈ [0, 1] where included. Although no explicit cutoff was imposed in ξ, the distribution decreases rapidly as this parameter decreases from 1 (as can be seen in Figures 4 and 5) and there are essentially no particles below some effective cutoff value. Figure 6 also shows the synchrotron spectrum from a single particle with momentum corresponding to the maximum momentum of the distributions (w = 50), and pitch-angle y ⊥ /y = 0.2, which is qualitatively consistent with measurements of the maximum runaway pitch on Textor [19] and DIII-D [20]. This choice is equivalent to using a 2D δ-function approximation of the distribution, as was done in Refs. [19,21] (and with some modification in [20]). The figure shows that using this approximation overestimates the synchrotron emission per particle by several orders of magnitude. Note that in the figure, the values for the emitted power per particle were divided by a large number to fit in the same scale.
The overestimation is not surprising, since the δ-function approximation assumes that all particles emit as much synchrotron radiation as the most strongly emitting particle in the actual distribution. The figure also shows that the wavelength of peak emission is shifted towards shorter wavelengths when using the δ-function approximation. In order to obtain an accurate runaway synchrotron spectrum, it is thus crucial to use the full runaway distribution in the calculation.
In principle, we may also use the synchrotron spectra from distributions calculated through CODE to estimate the maximum energy of the runaways in existing tokamaks.
However, due to the region of sensitivity of the available detectors, there is only a limited wavelength range in which calculated spectra can be fitted to experimental data in order to determine the maximum runaway energy. The available range often corresponds to the short wavelength slope of the spectrum, where the emitted power shows an approximately linear dependence on wavelength. Indeed, the short-wavelength spectrum slope has been used to estimate the maximum runaway energy in experiments [19]. If the runaway distribution function is approximated by a δ-function at the maximum available energy and pitch angle, there is a monotonic relationship between the short-wavelength spectrum slope and the maximum particle energy (at fixed pitch angle). Such a relationship holds because increasing the particle energy leads to more emission at shorter wavelengths, resulting in a shift of the wavelength of peak emission towards shorter wavelengths, and a corresponding change in the spectrum slope.
Using an integrated synchrotron spectrum from a CODE distribution is much more accurate than the single particle approximation, but also introduces additional parameters (Ê, δ, Z). If the physical parameters are well known, a unique relation still holds between the spectrum slope and the maximum particle energy. During disruptions however, many parameters (like the temperature and the effective charge) are hard to measure with accuracy.
As the shape of the underlying distribution depends on the values of the parameters, the synchrotron spectrum will do so as well. This complexity is apparent in Figure 6, where the single particle approximation produces identical results in the two cases, whereas the spectra from the complete distributions are widely different. The dependence on distribution shape makes it possible in principle for two sets of parameters to produce the same spectrum slope for different maximum energies. Given this insight, using the complete runaway distribution when modelling experimentally obtained spectra is necessary for an accurate analysis and reliable fit of the maximum particle energy. In this context, CODE is a very useful tool with the possibility to contribute to the understanding of runaways and their properties.
In providing the electron distribution functions (and thus knowledge of a variety of quantities through its moments), the applicability of CODE is wide, and the potential in coupling CODE to other software, e.g. for modelling of runaway dynamics in disruptions, is promising. For a proper description of the runaways generated in disruptions it is important to take into account the evolution of the radial profiles of the electric field and fast electron current self-consistently. This can be done by codes such as GO, initially described in Ref. [22] and developed further in Refs. [23,24]. GO solves the equation describing the resistive diffusion of the electric field in a cylindrical approximation coupled to the runaway generation rates. In the present version of GO, the runaway rate is computed by approximate analytical formulas for the primary and secondary generation. Using CODE, the analytical formulas can be replaced by a numerical solution for the runaway rate which would have several advantages. One advantage would be that Dreicer, hot-tail and secondary runaways could all be calculated with the same tool, avoiding the possibilities for double-counting and difficulties with interpretations of the results. Also, in the present version of GO, it is assumed that all the runaway electrons travel at the speed of light, an approximation that can be easily relaxed using CODE, which calculates the electron distribution in both energy and pitch-angle. Most importantly, the validity region of the results would be expanded, as the analytical formulas are derived using various assumptions which are often violated in realistic situtions. The output would be a self-consistent time and space evolution of electric field and runaway current, together with the electron distribution function. This information can then be used for calculating quantities that depend on the distribution function, such as the synchrotron emission or the kinetic instabilities driven by the velocity anisotropy of the runaways.

VII. CONCLUSIONS
The CODE code provides a fast computation of the runaway electron distribution function. Both primary (Dreicer) and secondary (avalanche) generation are included. A Legendre spectral disctretization is applied to the pitch-angle coordinate, with high-order finite differences applied to the speed coordinate. A nonuniform speed grid allows high resolution of thermal particles at the same time as a high maximum energy without a prohibitively large number of grid points. If secondary generation is negligible, the long-time distribution function may be calculated by solving a single sparse linear system. The speed of the code makes it feasible to couple to other codes for integrated modelling of complex processes such as tokamak disruptions. CODE has been benchmarked against previous analytic and numerical results in appropriate limits, showing excellent agreement.
As an example application, we investigated the synchrotron radiation spectra from runaway electron distributions. We found that the radiation spectrum from a single electron can differ substantially from the overall spectrum generated by a distribution of electrons.
Therefore, experimental estimates of maximum runaway energy based on the single-particle synchrotron spectrum are likely to be inaccurate. A detailed study of the distributionintegrated synchrotron spectrum and its dependencies on physical parameters will be presented in a future publication.
Similarly, to evaluate the ξ integral of the ∂F/∂ξ term in (11)