2 Kinetic equation for positrons

In this paper we consider the dynamics of positrons during a relativistic electron runaway avalanche (Jayakumar, Fleischmann & Zweben Reference Jayakumar, Fleischmann and Zweben1993). Due to the non-monotonic dynamical friction acting on a charged test particle in a plasma, in an electric field larger than a critical value $E_{c}$ fast electrons may experience a net force that can rapidly accelerate them to energies in the range of tens of MeV. In a fully ionized plasma, the critical field is $E_{c}=\ln \unicode[STIX]{x1D6EC}n_{e}e^{3}/(4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}^{2}m_{e}c^{2})$ (Connor & Hastie Reference Connor and Hastie1975), where $\ln \unicode[STIX]{x1D6EC}\approx 14.6+0.5\ln (T[\text{eV}]/n_{e}[10^{20}~\text{m}^{-3}])$ is the Coulomb logarithm (Solodov & Betti Reference Solodov and Betti2008). We neglect a logarithmic energy dependence in $\ln \unicode[STIX]{x1D6EC}$ , and use the value for relativistic electrons at 1 MeV for simplicity. Here, $n_{e}$ is the electron density, $e$ the elementary charge, $\unicode[STIX]{x1D700}_{0}$ the vacuum permittivity, $m_{e}$ the electron rest mass and $c$ the speed of light. The background plasma is assumed to be nearly Maxwellian for all species with the same temperature $T$ . In a neutral gas, $\ln \unicode[STIX]{x1D6EC}$ depends on the mean excitation energy of the medium instead of the temperature, and corresponds to $\ln \unicode[STIX]{x1D6EC}\approx 11$ in air (Gurevich & Zybin Reference Gurevich and Zybin2001). In this case the electron density refers to the density of bound electrons.

A sufficiently energetic electron can produce new runaway electrons through elastic large-angle collisions. The result is an exponentially growing number of runaway electrons, a so-called runaway avalanche. Each $e$ -folding of the number density takes a time $t_{\text{ava}}=c_{Z}/[4\unicode[STIX]{x03C0}n_{e}r_{0}^{2}c(E/E_{c}-1)]$ where $c_{Z}$ is only weakly dependent on electric field, and can be approximated by $c_{Z}\approx \sqrt{5+Z_{\text{eff}}}$ in a fully ionized plasma (Rosenbluth & Putvinski Reference Rosenbluth and Putvinski1997), where the effective charge is $Z_{\text{eff}}=\sum n_{i}Z_{i}^{2}/\sum _{i}n_{i}Z_{i}$ with the sum taken over all ion species $i$ . We shall find that several results in the paper are insensitive to the details of $c_{Z}$ , assuming only that it is independent of $E$ . As such, more accurate models of the avalanche process can in principle be implemented by inserting for $c_{Z}$ the value characterizing any particular scenario of interest.

Since the electrons are ultra-relativistic, they will create positrons which are predominantly co-moving; these are created either directly in collisions or indirectly through the hard X-rays emitted in collisions (Heitler Reference Heitler1954), which can produce a pair in a subsequent interaction. Since the positrons experience an acceleration by the electric field in the direction opposite to the runaway-electron motion, they will immediately start decelerating. A fraction of these positrons will slow down to thermal speeds where they eventually annihilate, whereas the remainder obtain sufficiently large momenta perpendicular to the acceleration direction that they become runaway accelerated along the electric field, moving anti-parallel to the runaway electrons. Annihilation – which occurs at a rate that decreases with positron energy – does not have a significant effect on the dynamics of the energetic positrons since the avalanche rate is typically much faster, which is demonstrated in § 2.2.

Throughout this paper, we shall assume that the plasma is fully ionized. In a partially ionized plasma or a neutral gas, screening effects due to the bound electrons would enter into all binary interactions. In the 10 MeV energy range, these are however largely negligible for the pair-production mechanisms as well as for the emission of bremsstrahlung, meaning that they are to be calculated using the full nuclear charge of the target. The screening effects become significant when $p/m_{e}c\gtrsim 137/Z^{1/3}$ (Heitler Reference Heitler1954). Elastic Coulomb collisions are to a greater extent affected by screening effects, where the pitch-angle scattering rates may be reduced by approximately up to two thirds and energy loss rates by one third (Hesslow et al. Reference Hesslow, Embréus, Stahl, Dubois, Papp, Newton and Fülöp2017) in the energy range of interest, compared to the results obtained treating the medium as fully ionized. This would modify primarily two important quantities that affect our results: the avalanche growth rate factor $c_{Z}$ , as well as the critical field $E_{c}$ (Hesslow et al. Reference Hesslow, Embréus, Wilkie, Papp and Fülöp2018), which can here be assumed to be accurate only up to an order-of-unity factor in partially ionized plasmas. While the results we present are strictly valid for a fully ionized plasma, we expect to capture the correct order of magnitude also in a partially ionized plasma or neutral gas, if the effective charge and electron densities appearing in the formulas are always evaluated using the fully ionized values. We denote these by

(2.1) $$\begin{eqnarray}\displaystyle & \displaystyle n_{\text{tot}}=\mathop{\sum }_{i}Z_{i}n_{i}, & \displaystyle\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & \displaystyle Z_{\text{tot}}=\frac{1}{n_{\text{tot}}}\mathop{\sum }_{i}n_{i}Z_{i}^{2}, & \displaystyle\end{eqnarray}$$

where $Z_{i}$ is the atomic number of species $i$ . Thus, the density is always to be taken as the total density of free plus bound electrons, and in a single-component gas or plasma $Z_{\text{tot}}$ is the atomic number of the ion species regardless of ionization degree.

The dynamics described above can be most lucidly captured in a two-dimensional model. The distribution function of positrons with momentum $\boldsymbol{p}=m_{e}\boldsymbol{v}/\sqrt{1-v^{2}/c^{2}}$ , where the positron velocity is denoted $\boldsymbol{v}$ , at a time $t$ is denoted $f_{\text{pos}}(t,\boldsymbol{p})$ . In a homogeneous cylindrically symmetric plasma in the presence of an electric field $\boldsymbol{E}$ it satisfies the kinetic equation

(2.3) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}f_{\text{pos}}}{\unicode[STIX]{x2202}t}+eE\left[\unicode[STIX]{x1D709}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}p}+\frac{1-\unicode[STIX]{x1D709}^{2}}{p}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D709}}\right]f_{\text{pos}}=C_{\text{pos}}+S_{\text{pos}}+S_{\text{an}},\end{eqnarray}$$

where $E=|\boldsymbol{E}|$ , $p=|\boldsymbol{p}|$ , $\unicode[STIX]{x1D709}\equiv \cos \unicode[STIX]{x1D703}=\boldsymbol{p}\boldsymbol{\cdot }\boldsymbol{E}/pE$ is the pitch-angle cosine, $C_{\text{pos}}$ is the positron collision operator, $S_{\text{pos}}$ denotes the source term of positrons generated in collisions between relativistic runaway electrons and field particles of the plasma, as well as positron production by highly energetic photons, and $S_{\text{an}}$ denotes the annihilation term. In a magnetized plasma, the equation is valid for an axisymmetric positron distribution if $E$ is replaced by the component of the electric field parallel to the magnetic field, and the pitch angle is instead defined relative to the magnetic field.

In the limit of small-energy transfers, the elastic positron–electron and positron–ion differential scattering cross-sections coincide with the electron–electron and electron–ion cross-sections, respectively (Landau & Lifshitz Reference Landau and Lifshitz1983). Consequently, the positron collision operator $C_{\text{pos}}$ equals the electron collision operator $C_{e}$ up to terms small in the Coulomb logarithm $\ln \unicode[STIX]{x1D6EC}$ . Large-angle collisions, which are primarily important for avalanche generation when $\ln \unicode[STIX]{x1D6EC}$ is large, can be neglected since the thermal positron population will always be small in number. The positron distribution therefore satisfies the same kinetic equation as the electron distribution, except for the electric field accelerating them in the opposite direction (with these definitions positrons are accelerated towards $\unicode[STIX]{x1D709}=1$ , and electrons towards $\unicode[STIX]{x1D709}=-1$ ), and the presence of the terms $S_{\text{pos}}$ and $S_{\text{an}}$ describing their creation and annihilation, respectively.

The number of positrons created with momentum $\boldsymbol{p}$ in time $\text{d}t$ has two main contributions: (i) the collisions between stationary ions of species $i$ with density $n_{i}$ and the number $\text{d}n_{\text{RE}}$ of runaways at momentum $\boldsymbol{p}_{1}$ and speed $v_{1}$ , and (ii) the pair production of $\text{d}n_{\unicode[STIX]{x1D6FE}}$ photons in the field of ions $i$ :

(2.4) $$\begin{eqnarray}\text{d}n_{\text{pos}}=\mathop{\sum }_{i}[n_{i}v_{1}\,\text{d}\unicode[STIX]{x1D70E}_{ci}^{+}\,\text{d}n_{\text{RE}}\,\text{d}t+n_{i}c\,\text{d}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FE}i}^{+}\,\text{d}n_{\unicode[STIX]{x1D6FE}}\,\text{d}t].\end{eqnarray}$$

Here, $\text{d}\unicode[STIX]{x1D70E}_{ci}^{+}$ is the differential cross-section for producing a positron in a collision between an electron and a stationary ion, and similarly $\text{d}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FE}i}^{+}$ for a photon interacting with stationary ions, and are given in appendix A. We use the Madgraph 5 tool (Alwall et al. Reference Alwall, Frederix, Frixione, Hirschi, Maltoni, Mattelaer, Shao, Stelzer, Torrielli and Zaro2014) for obtaining the pair-production cross-sections throughout this paper. In appendix A, figure 6, we compare the total cross-section as a function of electron energy with a previously published formula by Gryaznykh (Reference Gryaznykh1998). The latter is shown to be approximately a constant factor of four times as large as the cross-section that we have obtained.

Using $\text{d}n_{\text{RE}}(\,\boldsymbol{p}_{1})=f_{\text{RE}}(\,\boldsymbol{p}_{1})\text{d}\boldsymbol{p}_{1}$ , where $f_{\text{RE}}$ is the distribution function of runaway electrons, and similarly for the positron distribution $f_{\text{pos}}(\boldsymbol{p})=\text{d}n_{\text{pos}}/\text{d}\boldsymbol{p}$ and photons $\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D6FE}}(\boldsymbol{k})=\text{d}n_{\unicode[STIX]{x1D6FE}}/\text{d}\boldsymbol{k}$ where $\boldsymbol{k}/c$ is the photon momentum, we find the following form for the positron source $S_{\text{pos}}$ :

(2.5) $$\begin{eqnarray}S_{\text{pos}}\equiv \left(\frac{\unicode[STIX]{x2202}f_{\text{pos}}}{\unicode[STIX]{x2202}t}\right)_{\text{pp}}=\mathop{\sum }_{i}n_{i}c\left[\int \text{d}\boldsymbol{p}_{1}\,\frac{v_{1}}{c}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{ci}^{+}}{\unicode[STIX]{x2202}\boldsymbol{p}}f_{\text{RE}}(\,\boldsymbol{p}_{1})+\int \text{d}\boldsymbol{k}\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FE}i}^{+}}{\unicode[STIX]{x2202}\boldsymbol{p}}\unicode[STIX]{x1D719}_{\unicode[STIX]{x1D6FE}}(\boldsymbol{k})\right].\end{eqnarray}$$

In an avalanching runaway scenario, the photon distribution can be eliminated in favour of an expression involving only the runaway distribution because of the relatively slow evolution of the photon energy spectrum. The runaway-electron population grows exponentially in time on the time scale (Rosenbluth & Putvinski Reference Rosenbluth and Putvinski1997)

(2.6) $$\begin{eqnarray}t_{\text{ava}}=\frac{c_{Z}}{4\unicode[STIX]{x03C0}n_{\text{tot}}r_{0}^{2}c(E/E_{c}-1)}.\end{eqnarray}$$

The photons on the other hand evolve on the Compton-scattering time scale (Heitler Reference Heitler1954)

(2.7) $$\begin{eqnarray}t_{\text{Co}}=\frac{k}{\unicode[STIX]{x03C0}n_{\text{tot}}r_{0}^{2}cm_{e}c^{2}\ln [2k/(m_{e}c^{2})]},\end{eqnarray}$$

where the photon energies $k=|\boldsymbol{k}|$ are larger than the pair-production threshold $2m_{e}c^{2}$ , and $r_{0}=e^{2}/(4\unicode[STIX]{x03C0}\unicode[STIX]{x1D700}_{0}m_{e}c^{2})\approx 2.82\times 10^{-15}$  m is the classical electron radius. Comparing the two time scales shows that the photons do not have time to change significantly from the distribution in which they are created whenever

(2.8) $$\begin{eqnarray}\frac{k/m_{e}c^{2}}{\ln (2k/m_{e}c^{2})}\gg \frac{c_{Z}}{4(E/E_{c}-1)}.\end{eqnarray}$$

Since the right-hand side is typically smaller than unity, this is generally well satisfied in an avalanching runaway scenario. The photon distribution is then given by

(2.9) $$\begin{eqnarray}\unicode[STIX]{x1D719}(\boldsymbol{k})=t_{\text{ava}}\mathop{\sum }_{i}n_{i}\int \text{d}\boldsymbol{p}_{1}\,v_{1}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\text{br,i}}}{\unicode[STIX]{x2202}\boldsymbol{k}}(\boldsymbol{k},\boldsymbol{p}_{1})f_{\text{RE}}(\boldsymbol{p}_{1}),\end{eqnarray}$$

where $\text{d}\unicode[STIX]{x1D70E}_{\text{br,i}}$ is the differential bremsstrahlung cross-section for interactions between electrons and particle species $i$ .

Since the cross-sections appearing in these formulas depend on target species only through $Z_{i}^{2}$ , the target charge squared (Heitler Reference Heitler1954), these may be factored out when screening effects are neglected, yielding a factor of the effective plasma charge $Z_{\text{tot}}$ when summed over $i$ . We shall therefore suppress the indices $i$ of the cross-sections by writing $\sum _{i}n_{i}\unicode[STIX]{x1D70E}_{ci}^{+}=n_{\text{tot}}Z_{\text{tot}}\unicode[STIX]{x1D70E}_{c}^{+}$ for collisional pair production, and $\sum _{i}n_{i}\unicode[STIX]{x1D70E}_{\text{br,i}}=n_{\text{tot}}(Z_{\text{tot}}+1)\unicode[STIX]{x1D70E}_{\text{br}}$ (and likewise for $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FE}}^{+}$ ) for the photon pair-production cross-sections. Here we have added the contribution from electron–electron bremsstrahlung in the approximation that the electron–electron and electron–proton bremsstrahlung cross-sections are the same, which has satisfactory accuracy since the majority of interactions occur with negligible momentum transfer to the target particle (Haug Reference Haug1975). Conversely, for collisional pair production the electron–electron interactions are negligible, which was verified with MadGraph 5 simulations (Alwall et al. Reference Alwall, Frederix, Frixione, Hirschi, Maltoni, Mattelaer, Shao, Stelzer, Torrielli and Zaro2014) which indicated that the $e$ – $e$ cross-section is 10 %–20 % of the $e$ – $i$ cross-section when the incident electron laboratory-frame energy ranges over 10–20 MeV and $Z_{i}=1$ .

The positron source can then be written

(2.10) $$\begin{eqnarray}S_{\text{pos}}=Z_{\text{tot}}n_{\text{tot}}\int \text{d}\boldsymbol{p}_{1}\,v_{1}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}^{+}}{\unicode[STIX]{x2202}\boldsymbol{p}}f_{\text{RE}}(\,\boldsymbol{p}_{1}),\end{eqnarray}$$

where the effective pair-production cross-section $\text{d}\unicode[STIX]{x1D70E}^{+}$ , accounting for both direct pair production in collisions as well as by X-rays, is given by

(2.11) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}^{+}}{\unicode[STIX]{x2202}\boldsymbol{p}}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{c}^{+}}{\unicode[STIX]{x2202}\boldsymbol{p}}+\frac{(Z_{\text{tot}}+1)^{2}}{Z_{\text{tot}}}t_{\text{ava}}n_{\text{tot}}c\int \text{d}\boldsymbol{k}\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FE}}^{+}}{\unicode[STIX]{x2202}\boldsymbol{p}}(\boldsymbol{p},\boldsymbol{k})\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\text{br}}}{\unicode[STIX]{x2202}\boldsymbol{k}}(\boldsymbol{k},\boldsymbol{p}_{1}).\end{eqnarray}$$

Positrons are created with a significant fraction of the energy of the incident electron that created them, but with a momentum perpendicular to the direction of the incident electron of order (Heitler Reference Heitler1954; Landau & Lifshitz Reference Landau and Lifshitz1983) $p_{\bot }\approx m_{e}c$ . This means that the differential cross-section for their production is strongly peaked in the direction of the incident electron; throughout this work we assume that it is delta distributed in the scattering angles, and write

(2.12) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{c}^{+}}{\unicode[STIX]{x2202}\boldsymbol{p}}=\frac{\unicode[STIX]{x1D6FF}(\cos \unicode[STIX]{x1D703}-\cos \unicode[STIX]{x1D703}_{1})}{2\unicode[STIX]{x03C0}(m_{e}c)^{2}p\unicode[STIX]{x1D6FE}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{c}^{+}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}}(p,p_{1}),\\ \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FE}}^{+}}{\unicode[STIX]{x2202}\boldsymbol{p}}=\frac{\unicode[STIX]{x1D6FF}(\cos \unicode[STIX]{x1D703}-\cos \unicode[STIX]{x1D703}_{k})}{2\unicode[STIX]{x03C0}(m_{e}c)^{2}p\unicode[STIX]{x1D6FE}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FE}}^{+}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}}(p,k),\\ \displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\text{br}}}{\unicode[STIX]{x2202}\boldsymbol{k}}=\frac{\unicode[STIX]{x1D6FF}(\cos \unicode[STIX]{x1D703}_{k}-\cos \unicode[STIX]{x1D703}_{1})}{2\unicode[STIX]{x03C0}k^{2}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\text{br}}}{\unicode[STIX]{x2202}k}(k,p_{1}),\end{array}\right\} & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6FE}=\sqrt{1+(p/m_{e}c)^{2}}$ is the Lorentz factor. The angles $\unicode[STIX]{x1D703}$ , $\unicode[STIX]{x1D703}_{1}$ and $\unicode[STIX]{x1D703}_{k}$ are the angles between the accelerating electric field $\boldsymbol{E}$ (or in a magnetized plasma the magnetic field $\boldsymbol{B}$ ) and $\boldsymbol{p}$ , $\boldsymbol{p}_{1}$ and $\boldsymbol{k}$ , respectively. With an axisymmetric runaway distribution $f_{\text{RE}}(\,\boldsymbol{p}_{1})=f_{\text{RE}}(p_{1},\cos \unicode[STIX]{x1D703}_{1})$ , we then obtain the approximated positron source term

(2.13) $$\begin{eqnarray}S_{\text{pos}}(\unicode[STIX]{x1D6FE},\cos \unicode[STIX]{x1D703})=\frac{n_{\text{tot}}Z_{\text{tot}}m_{e}c^{2}}{p\unicode[STIX]{x1D6FE}}\int _{\unicode[STIX]{x1D6FE}+2}^{\infty }\text{d}\unicode[STIX]{x1D6FE}_{1}\,(\unicode[STIX]{x1D6FE}_{1}^{2}-1)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}}^{+}f_{\text{RE}}(\unicode[STIX]{x1D6FE}_{1},\cos \unicode[STIX]{x1D703}),\end{eqnarray}$$

where the effective cross-section now takes the form

(2.14) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}^{+}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FE}_{1})=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{c}^{+}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FE}_{1})+\frac{(Z_{\text{tot}}+1)^{2}}{Z_{\text{tot}}}t_{\text{ava}}n_{\text{tot}}c\int _{\unicode[STIX]{x1D6FE}+1}^{\unicode[STIX]{x1D6FE}_{1}-1}\text{d}k\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FE}}^{+}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D6FE},k)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\text{br}}}{\unicode[STIX]{x2202}k}(k,\unicode[STIX]{x1D6FE}_{1}),\end{eqnarray}$$

depending only on the simpler integrated cross-sections that are only differential in the energy of the outgoing particle of interest. In (2.13) it is explicit that positrons are only generated in the direction of the incident energetic electrons – the electron distribution is sampled at the same pitch angle as the source. Further details of the positron source term in the presence of an avalanching electron distribution are given in appendix A, where we present the differential cross-sections used as well as illustrate typical shapes of the source term in (2.13).

Figure 1. Critical field $E_{\text{pp}}$ above which collisional positron production is the dominant pair-production mechanism in a uniform plasma, normalized to the avalanche threshold field $E_{c}$ , calculated from the two expressions given in (2.23) with $c_{Z}=\sqrt{5+Z_{\text{tot}}}$ .

The annihilation source takes the simpler form

(2.15) $$\begin{eqnarray}S_{\text{an}}(\boldsymbol{p})=-n_{\text{tot}}v\unicode[STIX]{x1D70E}_{\text{an}}(p)f_{\text{pos}}(\boldsymbol{p}),\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{\text{an}}$ is the cross-section for free positron–electron two-quanta annihilation against stationary target electrons (Heitler Reference Heitler1954)

(2.16) $$\begin{eqnarray}\unicode[STIX]{x1D70E}_{\text{an}}=\frac{\unicode[STIX]{x03C0}r_{0}^{2}}{(p/m_{e}c)(\unicode[STIX]{x1D6FE}+1)}\left[\frac{\unicode[STIX]{x1D6FE}^{2}+4\unicode[STIX]{x1D6FE}+1}{p/m_{e}c}\ln \left(\unicode[STIX]{x1D6FE}+\frac{p}{m_{e}c}\right)-\unicode[STIX]{x1D6FE}-3\right].\end{eqnarray}$$

2.1 Relative importance of pair production by collisions and by photons

A peculiar phenomenon occurs when considering pair production in the presence of a strong electric field, where the number of energetic electrons grows exponentially in time. Because there is a delay between the emission of photons and their subsequent pair production, if the electron population has time to grow by a significant amount during one such photon pair-production time, the direct positron generation in collisions may contribute with a relatively larger production of pairs. We will now proceed to derive the threshold electric field above which pair production in collisions is dominant due to this effect.

In order to evaluate the pair-production source terms we need an expression for the runaway-electron distribution. In a spatially uniform fully ionized plasma with constant electric field, when the runaway generation is dominated by the avalanche mechanism – i.e. by multiplication through large-angle collisions – it is given by (Fülöp et al. Reference Fülöp, Pokol, Helander and Lisak2006)

(2.17) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle f_{\text{RE}}(p,\unicode[STIX]{x1D709},t)=\frac{n_{\text{RE}}(t)A(p)}{2\unicode[STIX]{x03C0}m_{e}c\unicode[STIX]{x1D6FE}_{0}p^{2}}\frac{\exp \left[-{\displaystyle \frac{\unicode[STIX]{x1D6FE}}{\unicode[STIX]{x1D6FE}_{0}}}-A(p)(1+\unicode[STIX]{x1D709})\right]}{1-\text{e}^{-2A}},\\ \displaystyle A(p)=\frac{E/E_{c}+1}{Z_{\text{tot}}+1}\unicode[STIX]{x1D6FE},\\ \displaystyle n_{\text{RE}}(t)=n_{\text{RE}}(0)\text{e}^{t/t_{\text{ava}}},\\ \displaystyle \unicode[STIX]{x1D6FE}_{0}=c_{Z}\ln \unicode[STIX]{x1D6EC}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Our choice for $A$ differs slightly from that in Fülöp et al. (Reference Fülöp, Pokol, Helander and Lisak2006), however agrees in the limit $E\gg E_{c}$ , $p\gg m_{e}c$ and $1+\unicode[STIX]{x1D709}\ll 1$ where the solution is expected to be valid, but is here generalized to also capture the near-threshold limit $E\rightarrow E_{c}$ (Lehtinen, Bell & Inan Reference Lehtinen, Bell and Inan1999; Hesslow et al. Reference Hesslow, Embréus, Wilkie, Papp and Fülöp2018). When pitch-angle averaged, the electron distribution is given by

(2.18) $$\begin{eqnarray}{\mathcal{F}}_{\text{RE}}(p,t)=2\unicode[STIX]{x03C0}\int _{-1}^{1}\text{d}\unicode[STIX]{x1D709}\,f_{\text{RE}}(p,\unicode[STIX]{x1D709},t)=\frac{n_{\text{RE}}(t)}{m_{e}c\unicode[STIX]{x1D6FE}_{0}}\text{e}^{-\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D6FE}_{0}},\end{eqnarray}$$

where the average runaway energy is given by $\unicode[STIX]{x1D6FE}_{0}m_{e}c^{2}\approx (c_{Z}\ln \unicode[STIX]{x1D6EC}/2)$  MeV, which is typically of the order of 10–30 MeV in most scenarios of interest. Runaway-electron populations with energies in this range have been observed in avalanching scenarios in tokamaks (Paz-Soldan et al. Reference Paz-Soldan, Cooper, Aleynikov, Pace, Eidietis, Brennan, Granetz, Hollmann, Liu and Lvovskiy2017). In our results we will indicate the dependence on average runaway energy through the parameters $\unicode[STIX]{x1D6FE}_{0}$ or $c_{Z}$ .

The total number of pairs created per unit time and volume is obtained by integrating the positron source function (2.13) over all momenta, yielding

(2.19) $$\begin{eqnarray}\frac{\text{d}n_{\text{pair}}}{\text{d}t}=n_{e}Z_{\text{tot}}m_{e}c^{2}\int _{3}^{\infty }\text{d}\unicode[STIX]{x1D6FE}_{1}\,\unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D6FE}_{1}){\mathcal{F}}_{\text{RE}}(\unicode[STIX]{x1D6FE}_{1})\equiv n_{e}n_{\text{RE}}Z_{\text{tot}}\langle v\unicode[STIX]{x1D70E}^{+}\rangle _{\text{RE}},\end{eqnarray}$$

(2.20) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D6FE}_{1})=\int _{1}^{\unicode[STIX]{x1D6FE}_{1}-2}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}^{+}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}}(\unicode[STIX]{x1D6FE},\unicode[STIX]{x1D6FE}_{1})\,\text{d}\unicode[STIX]{x1D6FE},\\ \displaystyle \langle v\unicode[STIX]{x1D70E}^{+}\rangle _{\text{RE}}=\frac{1}{n_{\text{RE}}}\int _{\sqrt{8}}^{\infty }\text{d}p_{1}\,v_{1}\unicode[STIX]{x1D70E}^{+}(\unicode[STIX]{x1D6FE}_{1}){\mathcal{F}}_{\text{RE}}(\unicode[STIX]{x1D6FE}_{1}).\end{array}\right\} & & \displaystyle\end{eqnarray}$$

With the analytic form of (2.18) for the electron distribution, the pair-production rate defined by the above equations is characterized by the two integrals

(2.21) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle \!\!\!\langle v\unicode[STIX]{x1D70E}_{c}^{+}\rangle _{\text{RE}}=\frac{m_{e}c^{2}}{n_{\text{RE}}}\int _{3}^{\infty }\text{d}\unicode[STIX]{x1D6FE}_{1}\,\unicode[STIX]{x1D70E}_{c}^{+}(\unicode[STIX]{x1D6FE}_{1}){\mathcal{F}}_{\text{RE}}(\unicode[STIX]{x1D6FE}_{1})\approx \unicode[STIX]{x1D6FC}^{2}r_{0}^{2}c\frac{\unicode[STIX]{x1D6FE}_{0}-6.7}{15},\\ \displaystyle \!\!\!\langle v\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FE}}\rangle _{\text{RE}}=\frac{m_{e}c^{2}}{n_{\text{RE}}}\int _{3}^{\infty }\text{d}\unicode[STIX]{x1D6FE}_{1}\,{\mathcal{F}}_{\text{RE}}(\unicode[STIX]{x1D6FE}_{1})\int _{2}^{\unicode[STIX]{x1D6FE}_{1}-1}\text{d}k\,\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FE}}^{+}(k)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\text{br}}}{\unicode[STIX]{x2202}k}(k,\unicode[STIX]{x1D6FE}_{1})\approx \unicode[STIX]{x1D6FC}^{2}r_{0}^{4}c(2.6\unicode[STIX]{x1D6FE}_{0}-14.8),\end{array}\!\right\} & & \displaystyle \nonumber\\ \displaystyle \qquad \qquad & & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FE}}^{+}=\int _{1}^{k-1}(\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FE}}^{+}/\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE})\,\text{d}\unicode[STIX]{x1D6FE}$ , and the approximate formulas are least-square fits on the interval of $\unicode[STIX]{x1D6FE}_{0}$ between 20 and 80, giving a maximum error of 2.5 %. Within an error of less than $3\,\%$ , the second expression differs from the first by a constant factor $40.75r_{0}^{2}$ , allowing the total pair-production rate to be written

(2.22) $$\begin{eqnarray}\frac{\text{d}n_{\text{pair}}}{\text{d}t}\approx Z_{\text{tot}}\unicode[STIX]{x1D6FC}^{2}n_{e}r_{0}^{2}c\frac{\unicode[STIX]{x1D6FE}_{0}-6.7}{15}\left(1+40.75\frac{(Z_{\text{tot}}+1)^{2}}{Z_{\text{tot}}}t_{\text{ava}}n_{e}cr_{0}^{2}\right).\end{eqnarray}$$

With $t_{\text{ava}}n_{e}cr_{0}^{2}=c_{Z}/[4\unicode[STIX]{x03C0}(E/E_{c}-1)]$ , it is clear that there is a threshold field $E=E_{\text{pp}}(Z_{\text{tot}})$ above which the collisional pair production (described by the first term within the parentheses) will be dominant. When $c_{Z}$ is independent of $E$ , this threshold field is given by

(2.23) $$\begin{eqnarray}\displaystyle \frac{E_{\text{pp}}}{E_{c}}-1 & = & \displaystyle \frac{(1+Z_{\text{tot}})^{2}}{Z_{\text{tot}}}\frac{c_{Z}}{4\unicode[STIX]{x03C0}r_{0}^{2}}\frac{\displaystyle \int _{3}^{\infty }\text{d}\unicode[STIX]{x1D6FE}_{1}\,\text{e}^{-\unicode[STIX]{x1D6FE}_{1}/\unicode[STIX]{x1D6FE}_{0}}\int _{2}^{\unicode[STIX]{x1D6FE}_{1}-1}\text{d}k\,\unicode[STIX]{x1D70E}_{\unicode[STIX]{x1D6FE}}^{+}(k)\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\text{br}}}{\unicode[STIX]{x2202}k}(k,\unicode[STIX]{x1D6FE}_{1})}{\displaystyle \int _{3}^{\infty }\text{d}\unicode[STIX]{x1D6FE}_{1}\,\text{e}^{-\unicode[STIX]{x1D6FE}_{1}/\unicode[STIX]{x1D6FE}_{0}}\unicode[STIX]{x1D70E}_{c}^{+}(\unicode[STIX]{x1D6FE}_{1})}\nonumber\\ \displaystyle & {\approx} & \displaystyle 3.25c_{Z}\frac{(1+Z_{\text{tot}})^{2}}{Z_{\text{tot}}}.\end{eqnarray}$$

When $Z_{\text{tot}}=1$ , the threshold field is $E_{\text{pp}}\approx 33E_{c}$ , but grows rapidly with $Z_{\text{tot}}$ . With an air-like value of $Z_{\text{tot}}=14.5$ , one obtains $E_{\text{pp}}\approx 240E_{c}$ .

The role of direct pair production in lightning has been investigated by Vodopiyanov et al. (Reference Vodopiyanov, Dwyer, Cramer, Lucia and Rassoul2015) using comprehensive Monte Carlo simulations of the runaway dynamics to investigate the degree to which positrons contribute to the relativistic feedback mechanism (Dwyer Reference Dwyer2012). Their numerical study revealed a similar qualitative dependence on the electric field when comparing the contributions from photon-produced and directly produced positrons; direct production becomes relatively more important for larger electric fields. Their study, however, uses the Gryaznykh (Reference Gryaznykh1998) cross-section for pair production as well as a simplified description of the energy spectrum of newly created positrons, and they find a significantly lower threshold field than predicted here. Spatial dynamics along the discharge may also play a role in determining the feedback gain, which prevents a simple direct comparison between our results.

In the above we assumed an infinitely large homogeneous system. When runaway acceleration occurs only over a finite distance of length $L$ of constant background parameters, the threshold field calculated above is valid when $L\gg L_{\text{ava}}=ct_{\text{ava}}$ , that is, when the system is significantly longer than one avalanche mean-free path. When this is not satisfied, i.e. when $L\lesssim L_{\text{ava}}$ , a threshold condition for the length of the system is obtained instead, taking the form

(2.24) $$\begin{eqnarray}L_{\text{pp}}\approx \frac{Z_{\text{tot}}}{(1+Z_{\text{tot}})^{2}}\frac{3\times 10^{7}\,\text{m}}{n_{e}[10^{20}\,\text{m}^{-3}]}.\end{eqnarray}$$

When $L_{\text{pp}}\lesssim L\lesssim L_{\text{ava}}$ , photon pair production will be the dominant positron-generation mechanism.

In the remainder of this work, we will focus on scenarios where either $E\gtrsim E_{\text{pp}}$ or $L\lesssim L_{\text{pp}}$ , so that pair production by photons is negligible. This is typical of runaway scenarios in tokamaks, where $L/L_{\text{pp}}\ll 10^{-5}$ due to the small size of the device.

2.2 Distribution function of fast positrons

Equipped with the kinetic equation for positrons in a runaway scenario, we can now characterize its solutions. When the electric field is sufficiently large for the average pitch angle to be small, typically well satisfied when $A=(E/E_{c}+1)\unicode[STIX]{x1D6FE}/(Z_{\text{tot}}+1)\gtrsim 1$ , the distribution function of fast positrons can be readily calculated analytically. The kinetics are then essentially one-dimensional, with the pitch-angle dynamics playing a peripheral role in the evolution of the energy spectrum.

We introduce a half-plane pitch-angle-averaged positron distribution function ${\mathcal{F}}$ as

(2.25) $$\begin{eqnarray}{\mathcal{F}}(p)=2\unicode[STIX]{x03C0}p^{2}\times \left\{\begin{array}{@{}ll@{}}\displaystyle \int _{0}^{1}\text{d}\unicode[STIX]{x1D709}\,f_{\text{pos}}(p,\unicode[STIX]{x1D709}),\quad & p\geqslant 0\\ \displaystyle \int _{-1}^{0}\text{d}\unicode[STIX]{x1D709}\,f_{\text{pos}}(|p|,\unicode[STIX]{x1D709}),\quad & p<0,\end{array}\right.\end{eqnarray}$$

where the coordinate $p$ now ranges from $-\infty$ to $\infty$ . This distribution is defined so that $\int _{p_{c}}^{\infty }{\mathcal{F}}\,\text{d}p=n_{\text{RP}}$ equals the total runaway-positron density, with $p_{c}$ a superthermal threshold in momentum distinguishing thermal positrons from runaways. In the same way, the thermal number density of positrons is $n_{\text{TP}}=\int _{-p_{c}}^{p_{c}}{\mathcal{F}}\,\text{d}p$ .

In appendix B we solve the positron kinetic equation (2.3) in the limit $(p/m_{e}c)^{2}\gg 1$ assuming small pitch angles $1-|\unicode[STIX]{x1D709}|\ll 1$ . The resulting positron distribution is given by

(2.26) $$\begin{eqnarray}{\mathcal{F}}(p)=\frac{Z_{\text{tot}}}{4\unicode[STIX]{x03C0}\ln \unicode[STIX]{x1D6EC}r_{0}^{2}}\frac{n_{\text{RE}}(t)}{\unicode[STIX]{x1D6FE}_{0}m_{e}c}\text{e}^{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D6FE}_{0}}\int _{\unicode[STIX]{x1D6FE}}^{\infty }\text{d}\unicode[STIX]{x1D6FE}^{\prime }\int _{\unicode[STIX]{x1D6FE}^{\prime }+2}^{\infty }\text{d}\unicode[STIX]{x1D6FE}_{1}\,\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}^{+}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FE}^{\prime }}(\unicode[STIX]{x1D6FE}^{\prime },\unicode[STIX]{x1D6FE}_{1})\exp \left(-\frac{\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FE}^{\prime }+\unicode[STIX]{x1D6FE}_{1}}{\unicode[STIX]{x1D6FE}_{0}}\right)\end{eqnarray}$$

for $p<0$ , which describes the slowing-down distribution of the newly created positrons, where $\unicode[STIX]{x1D70C}=(E/E_{c}-1)/(E/E_{c}+1)$ , and

(2.27) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle {\mathcal{F}}(p,t)=\frac{n_{\text{RP}}(t)}{m_{e}c\unicode[STIX]{x1D6FE}_{0}}\text{e}^{-\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D6FE}_{0}},\\ \displaystyle n_{\text{RP}}(t)=n_{\text{RP}}(0)\text{e}^{t/t_{\text{ava}}}\end{array}\right\} & & \displaystyle\end{eqnarray}$$

for $p>0$ , describing the runaway-positron population that is undergoing acceleration in the electric field. Note that the prefactor, including the runaway-positron density evaluated at $t=0$ , is not determined in this derivation, but must instead be calculated in a more comprehensive kinetic equation accounting for the dynamics near $p\lesssim m_{e}c$ .

We see that when $p>0$ , the runaway-positron distribution satisfies ${\mathcal{F}}(p)=(n_{\text{RP}}/n_{\text{RE}}){\mathcal{F}}_{\text{RE}}(-p)$ . Indeed, for the full positron distribution, since the kinetic equation (2.3) is identical to the runaway-electron equation for $\unicode[STIX]{x1D709}>0$ where the pair-production source vanishes, we would expect

(2.28) $$\begin{eqnarray}f_{\text{RP}}(p,\unicode[STIX]{x1D709},t)\approx \frac{n_{\text{RP}}(t)}{n_{\text{RE}}(t)}f_{\text{RE}}(p,-\unicode[STIX]{x1D709},t).\end{eqnarray}$$

The expressions given above are valid for collisional as well as for photon pair production during runaway scenarios.

We can now accurately evaluate the annihilation rate of runaway positrons, obtaining in the ultra-relativistic limit,

(2.29) $$\begin{eqnarray}\displaystyle \frac{1}{\unicode[STIX]{x1D70F}_{\text{aR}}}=\frac{n_{e}}{n_{\text{RP}}}\int _{p_{c}}^{\infty }\text{d}p\,v{\mathcal{F}}(p)\unicode[STIX]{x1D70E}_{\text{an}}(p) & {\approx} & \displaystyle \frac{1}{4\ln \unicode[STIX]{x1D6EC}(E/E_{c}-1)t_{\text{ava}}}\int _{1}^{\infty }\text{d}\unicode[STIX]{x1D6FE}\frac{\ln 2\unicode[STIX]{x1D6FE}-1}{\unicode[STIX]{x1D6FE}}\text{e}^{-\unicode[STIX]{x1D6FE}/\unicode[STIX]{x1D6FE}_{0}}\nonumber\\ \displaystyle & {\approx} & \displaystyle \frac{\ln ^{2}(\unicode[STIX]{x1D6FE}_{0}/2.42)+1.55}{8\ln \unicode[STIX]{x1D6EC}(E/E_{c}-1)t_{\text{ava}}},\end{eqnarray}$$

the final approximation having an error less than $2\,\%$ for $\unicode[STIX]{x1D6FE}_{0}>20$ , and where the annihilation cross-section $\unicode[STIX]{x1D70E}_{\text{an}}$ was given in (2.16). We find that typically $t_{\text{ava}}/\unicode[STIX]{x1D70F}_{\text{aR}}\lesssim 0.1/(E/E_{c}-1)$ , showing that annihilation has negligible impact on the avalanche time scale dynamics except for very close to the threshold field $E_{c}$ . At that point, however, most of the created positrons will become thermalized, and only a negligible fraction will have time to annihilate before reaching thermal energies.

For $v\ll c$ the annihilation cross-section takes the simple form $\unicode[STIX]{x1D70E}_{\text{an}}\sim \unicode[STIX]{x03C0}r_{0}^{2}c/v$ , so that the thermal annihilation time $\unicode[STIX]{x1D70F}_{\text{aT}}$ for a thermal positron population of temperature $T\ll 511$  keV is given simply by

(2.30) $$\begin{eqnarray}\frac{1}{\unicode[STIX]{x1D70F}_{\text{aT}}}=\unicode[STIX]{x03C0}n_{e}r_{0}^{2}c=\frac{1}{4t_{\text{ava}}}\frac{c_{Z}}{E/E_{c}-1}.\end{eqnarray}$$

In the presence of partially ionized or neutral gases, however, the cold positrons may annihilate also through the formation of positronium, which has significantly shorter, sub- $\unicode[STIX]{x03BC}\text{s}$ lifetime. The annihilation time of thermal positrons is then rather set by the positronium formation rate, which is of the order of $n_{i}va_{0}^{2}$ , with $a_{0}$ the Bohr radius (Charlton & Humberston Reference Charlton and Humberston2001).

2.3 Numerical distribution function

The positron Fokker–Planck equation, equation (2.3), can be solved as an initial value problem to give the evolution of the positron distribution function in the presence of an accelerating electric field. By adding the source and annihilation terms to the code (Landreman et al. Reference Landreman, Stahl and Fülöp2014; Stahl et al. Reference Stahl, Embréus, Papp, Landreman and Fülöp2016) numerical kinetic solver, we calculate the distribution function for various electric fields and effective ion charges. code uses a continuum-spectral discretization scheme and has been used extensively to calculate runaway-electron distributions including partial screening effects (Hesslow et al. Reference Hesslow, Embréus, Stahl, Dubois, Papp, Newton and Fülöp2017), synchrotron radiation (Hirvijoki et al. Reference Hirvijoki, Pusztai, Decker, Embréus, Stahl and Fülöp2015; Stahl et al. Reference Stahl, Hirvijoki, Decker, Embréus and Fülöp2015), bremsstrahlung (Embréus, Stahl & Fülöp Reference Embréus, Stahl and Fülöp2016) and close collisions (Embréus, Stahl & Fülöp Reference Embréus, Stahl and Fülöp2018).

Figure 2. Pitch-angle-averaged distribution functions ${\mathcal{F}}$ after 10 avalanche times $t_{\text{ava}}$ , with an initial runaway-electron density $n_{\text{RE,0}}=10^{10}~\text{m}^{-3}$ , $n_{e}=5\times 10^{19}~\text{m}^{-3}$ and $T=100$  eV. Runaway electrons in red and positrons in black and blue. Dashed lines denote the theoretical predictions of (2.26) and (2.27); in the ( $ii$ )  $E=10E_{c}$ , $Z_{\text{tot}}=1$ case it fully overlaps with the numerical solution.

Figure 2 illustrates the angle-averaged positron distribution for two cases: (i) with $E=2E_{c}$ and $Z_{\text{tot}}=10$ , and (ii) $E=10E_{c}$ and $Z_{\text{tot}}=1$ . The analytic solution, equations (2.26) and (2.27), is nearly indistinguishable from the numerical solution for $p<0$ for both cases, and for $p>0$ in case ( $ii$ ) with the higher electric field (shown in blue). The analytic solution fails to fully capture the low energy behaviour in case ( $i$ ) with low electric field and high plasma charge (black), where pitch-angle dynamics becomes important. The accuracy of the analytic solution at high electric field further motivates the neglect of annihilation in the dynamics of fast positrons. The sharp peaks at $p=0$ in the numerical positron energy spectra contain the thermalized positron populations, which we do not consider the detailed dynamics of here.

3 Rate equations for runaway positrons

From the kinetic description of § 2 we can find a reduced set of fluid equations which govern the evolution of the number densities of runaway positrons as well as thermal positrons. We introduce the runaway-positron density $n_{\text{RP}}$ and thermal positron density $n_{\text{TP}}$ in the same way as in the previous section. These then satisfy the equations

(3.1) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}n_{\text{RP}}}{\unicode[STIX]{x2202}t}=Z_{\text{tot}}n_{e}n_{\text{RE}}\unicode[STIX]{x1D705}(E,Z_{\text{tot}})\langle v\unicode[STIX]{x1D70E}_{c}^{+}\rangle _{\text{RE}}-n_{\text{RP}}/\unicode[STIX]{x1D70F}_{\text{aR}} & \displaystyle\end{eqnarray}$$

(3.2) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}n_{\text{TP}}}{\unicode[STIX]{x2202}t}=Z_{\text{tot}}n_{e}n_{\text{RE}}\unicode[STIX]{x1D702}(E,Z_{\text{tot}})\langle v\unicode[STIX]{x1D70E}_{c}^{+}\rangle _{\text{RE}}-n_{\text{TP}}/\unicode[STIX]{x1D70F}_{\text{aT}} & \displaystyle\end{eqnarray}$$

(3.3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\unicode[STIX]{x2202}n_{\text{RE}}}{\unicode[STIX]{x2202}t}=n_{\text{RE}}\unicode[STIX]{x1D6E4}_{\text{ava}}(E,Z_{\text{tot}}), & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D705}$ denotes the fraction of created positrons that are accelerated as runaways, $\unicode[STIX]{x1D702}$ the fraction that is thermalized and $\unicode[STIX]{x1D6E4}_{\text{ava}}=1/t_{\text{ava}}$ is the avalanche growth rate of runaway electrons.

Figure 3. (a) Positron runaway fraction $\unicode[STIX]{x1D705}$ defined by (3.1), for various electric-field strengths and plasma effective charge $Z_{\text{tot}}$ . (b) Positron thermalization fraction $\unicode[STIX]{x1D702}$ defined by (3.2), for various electric-field strengths and plasma effective charge $Z_{\text{tot}}$ .

Kinetic simulations must be used to determine the runaway fraction $\unicode[STIX]{x1D705}$ and thermalization fraction $\unicode[STIX]{x1D702}$ (note that they will not sum to unity, since the population of fast newly born positrons with $\unicode[STIX]{x1D709}<0$ also grows in time). Results from numerical simulations for a variety of electric fields and plasma charges are shown in figure 3, obtained for constant electric fields and plasma charge. These are applicable to scenarios where $E$ and $Z_{\text{tot}}$ vary slowly in time compared to the avalanche time. When $E\gg E_{c}$ , the runaway fraction is near unity, but decreases exponentially in magnitude when the electric field approaches the threshold value $E_{c}$ . As only a small fraction of positrons are annihilated before slowing down (or entering the runaway region) (Heitler Reference Heitler1954), the effect of annihilation on the positron runaway-generation dynamics can be ignored, and we can assume that positrons are only annihilated after being either thermalized or runaway accelerated.

In the presence of a constant electric field, background density and charge, the rate equations have a simple analytic solution given by (after a short transient phase on the scale of $t_{\text{ava}}$ )

(3.4) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}\displaystyle n_{\text{RE}}(t)=n_{0}\exp (\unicode[STIX]{x1D6E4}_{\text{ava}}t),\\ \displaystyle n_{\text{RP}}(t)=n_{\text{RE}}(t)\frac{Z_{\text{tot}}\unicode[STIX]{x1D705}n_{e}\langle v\unicode[STIX]{x1D70E}_{c}^{+}\rangle }{\unicode[STIX]{x1D6E4}_{\text{ava}}+\unicode[STIX]{x1D70F}_{\text{aR}}^{-1}},\\ \displaystyle n_{\text{TP}}(t)=n_{\text{RE}}(t)\frac{Z_{\text{tot}}\unicode[STIX]{x1D702}n_{e}\langle v\unicode[STIX]{x1D70E}_{c}^{+}\rangle }{\unicode[STIX]{x1D6E4}_{\text{ava}}+\unicode[STIX]{x1D70F}_{\text{aT}}^{-1}}.\end{array}\right\} & & \displaystyle\end{eqnarray}$$

Note that the positron populations grow in time despite annihilation; this occurs due to the ever-increasing amplitude of the positron source, since the runaway electrons are avalanching.

When the electric field is significantly above threshold one finds that $\unicode[STIX]{x1D6E4}_{\text{ava}}\gg \unicode[STIX]{x1D70F}_{\text{aR}}^{-1}$ meaning that annihilation is negligible, so that

(3.5) $$\begin{eqnarray}\frac{n_{\text{RP}}}{n_{\text{RE}}}\approx \frac{\unicode[STIX]{x1D705}(E,Z_{\text{tot}})}{E/E_{c}-1}\frac{Z_{\text{tot}}\langle v\unicode[STIX]{x1D70E}_{c}^{+}\rangle }{4\unicode[STIX]{x03C0}cr_{0}^{2}/c_{Z}}\approx \unicode[STIX]{x1D6FC}^{2}c_{Z}Z_{\text{tot}}\frac{\unicode[STIX]{x1D705}(E,Z_{\text{tot}})}{E/E_{c}-1}\frac{\unicode[STIX]{x1D6FE}_{0}-6.7}{60\unicode[STIX]{x03C0}},\end{eqnarray}$$

where again $\unicode[STIX]{x1D6FE}_{0}=c_{Z}\ln \unicode[STIX]{x1D6EC}$ . The electric-field dependence is fully captured in the factor $\unicode[STIX]{x1D705}/(E/E_{c}-1)$ , which takes its maximal value ${\approx}0.2$ near $E\approx 2E_{c}$ , only weakly dependent on the charge $Z_{\text{tot}}$ . With $\ln \unicode[STIX]{x1D6EC}=15$ , we then find that the maximal ratio of runaway positrons to electrons is $n_{\text{RP}}/n_{\text{RE}}\lesssim 8.5\cdot 10^{-7}Z_{\text{tot}}c_{Z}(c_{Z}-0.45)$ , which for a low- $Z$ plasma with $Z_{\text{tot}}=1$ is approximately $4\times 10^{-6}$ , and for a high- $Z$ plasma with $Z_{\text{tot}}=20$ of the order of $4\times 10^{-4}$ . This means that the runaway-positron synchrotron and hard X-ray (HXR) emission may be challenging to distinguish from the radiation emitted by the runaway electrons in a tokamak, since even a small fraction of reflected or scattered radiation from electrons or noise from other sources could drown out the positron signal.

4 Radiation from positrons in tokamak plasmas

In the previous section we found that runaway positrons are less numerous than the runaway electrons by a factor smaller than approximately $10^{-4}$ . This causes a direct measurement of runaway positrons in a laboratory plasma to be challenging, and an appealing option is instead to detect the annihilation radiation of the positrons that have slowed down, which is distinctly peaked around photon energies of 511 keV. The annihilation radiation from slow positrons is emitted approximately isotropically, whereas runaway electrons emit radiation primarily along their direction of motion, which when the electric field is large is along the electric field, or along the magnetic field in a magnetized plasma. This means that when measuring perpendicularly to the direction of runaway acceleration, even though the positrons are much fewer, their annihilation radiation may be detected through the X-ray background of runaway electrons for which only a small fraction is emitted at a $\unicode[STIX]{x03C0}/2$ angle and nearFootnote ¹   $511$  keV. Furthermore, coincidence measurement techniques can be employed to carry out measurements in poor signal-to-noise ratio cases (Guanying et al. Reference Guanying, Liu, Xie and Li2017).

We can make the heuristic discussion above stricter by the following arguments. The number density of bremsstrahlung photons emitted per unit solid angle, time and photon energy is given by

(4.1) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}n_{\text{HXR}}}{\unicode[STIX]{x2202}t\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}\unicode[STIX]{x2202}k}=n_{e}Z_{\text{tot}}\int _{\unicode[STIX]{x1D6FE}>k+1}v\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{\text{br}}}{\unicode[STIX]{x2202}k\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}f_{\text{RE}}(\,\boldsymbol{p})\,\text{d}\boldsymbol{p}.\end{eqnarray}$$

This can be compared to the number density of annihilation photons emitted per unit time and solid angle due to the thermal positrons $n_{\text{TP}}$ annihilating against the cold background,

(4.2) $$\begin{eqnarray}\frac{\unicode[STIX]{x2202}n_{\text{an}}}{\unicode[STIX]{x2202}t\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}=\frac{n_{\text{TP}}}{4\unicode[STIX]{x03C0}\unicode[STIX]{x1D70F}_{\text{aT}}}\approx Z_{\text{tot}}\frac{n_{\text{RE}}n_{e}\unicode[STIX]{x1D702}\langle v\unicode[STIX]{x1D70E}_{c}^{+}\rangle _{\text{RE}}}{4\unicode[STIX]{x03C0}},\end{eqnarray}$$

where we have assumed the thermal positron-annihilation rate to be much larger than the avalanche growth rate, $\unicode[STIX]{x1D6E4}_{\text{ava}}\unicode[STIX]{x1D70F}_{\text{aT}}\ll 1$ .

The annihilation radiation will have a line profile in photon energy with a characteristic width comparable to the background temperature. We consider the case where the profile is not resolved in the measurement, and the full line is captured in one channel. In this case, since the hard X-rays have a broad spectrum, we find it useful to characterize the visibility of the annihilation line with the parameter

(4.3) $$\begin{eqnarray}\unicode[STIX]{x0394}k=\frac{\unicode[STIX]{x2202}n_{\text{an}}/\unicode[STIX]{x2202}t\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}{\unicode[STIX]{x2202}n_{\text{HXR}}/\unicode[STIX]{x2202}t\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}\unicode[STIX]{x2202}k},\end{eqnarray}$$

which (when $\unicode[STIX]{x0394}k\ll k$ ) can be interpreted as the photon-energy interval $\unicode[STIX]{x0394}k$ around $k=m_{e}c^{2}$ within which the total HXR emission equals the annihilation photon flux. From a detection point of view, $\unicode[STIX]{x0394}k$ would approximate the energy resolution required for the annihilation peak to appear with twice the amplitude of the continuous X-ray background. The finite line width of the annihilation peak would need to be accounted for when the plasma temperature satisfies $T\gtrsim \unicode[STIX]{x0394}k$ .

In figure 4 we show $\unicode[STIX]{x0394}k$ for detection at a $\unicode[STIX]{x03C0}/2$ angle relative to the direction of runaway acceleration, using the analytic runaway distribution of (2.17). We observe a relatively weak dependence on electric field where the main trend is approximately captured, within approximately 25 %, by

(4.4) $$\begin{eqnarray}\unicode[STIX]{x0394}k\approx \frac{7\,\text{keV}}{\sqrt{Z_{\text{tot}}+1}}.\end{eqnarray}$$

This means that in order for the annihilation peak to be clearly distinguishable from the X-ray background due to runaway electrons, an energy resolution better than or comparable to $(7/\sqrt{Z_{\text{tot}}+1})$  keV is desirable. The contrast of X-rays to annihilation radiation, quantified through $\unicode[STIX]{x0394}k$ , is largely insensitive to other parameters of the scenario, since the cross-sections for the two processes scale in the same way with the background parameters. Commercial high purity germanium (HPGe) radiation detectors are available with energy resolutions near 1 keV. Therefore, our calculations indicate that sufficient signal-to-noise ratios can be achieved in the presence of most magnetic-fusion relevant plasma compositions.

Figure 4. The photon-energy resolution parameter $\unicode[STIX]{x0394}k=(\unicode[STIX]{x2202}n_{\text{an}}/\unicode[STIX]{x2202}t\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA})/(\unicode[STIX]{x2202}n_{\text{HXR}}/\unicode[STIX]{x2202}t\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}\unicode[STIX]{x2202}k)$ for perpendicular detection of annihilation radiation from thermalized positrons and hard X-rays from runaway electrons.

There are two main competing effects which are sensitive to $E$ and $Z_{\text{tot}}$ that determine the observed behaviour in $\unicode[STIX]{x0394}k$ . When $E$ increases, the thermalization fraction $\unicode[STIX]{x1D702}$ of positrons rapidly decreases, as illustrated in figure 3(b), which reduces the amount of 511 keV annihilation radiation. At the same time the runaway-electron population becomes more anisotropic, which sharply reduces the amount of bremsstrahlung emitted at a perpendicular angle. In the parameter range shown in the figure, these effects are found to approximately cancel, leaving only a weak $E$ dependence. On the other hand, an increase in charge $Z_{\text{tot}}$ causes the electron population to become more isotropic, increasing the amount of bremsstrahlung emitted at a perpendicular angle, however, it also increases the average runaway-electron energy which increases the number of positrons created per electron. The former effect is significantly stronger, which causes a net $1/\sqrt{Z_{\text{tot}}+1}$ dependence.

Finally we note that, in the post-disruption runaway plateau where the runaway current slowly dissipates on the inductive time scale of the device, the analytical avalanche runaway distribution that we have used here is not valid, as it tends to significantly overestimate the average energy of the distribution. Due to its experimental accessibility, we consider this scenario separately for a singly ionized argon-dominated plasma (Pautasso et al. Reference Pautasso, Bernert, Dibon, Duval, Dux, Fable, Fuchs, Conway, Giannone and Gude2016). For the runaway electron distribution we use the self-consistent slowing-down distribution of Hesslow et al. (Reference Hesslow, Embréus, Wilkie, Papp and Fülöp2018) obtained from a numerical solution of the kinetic equation with an inductive electric field and accurate modelling of screening effects on collisions, at a plasma temperature $T=10$  eV. Using such a numerical distribution function in evaluating the rate of pair production $\langle v\unicode[STIX]{x1D70E}_{c}^{+}\rangle _{\text{RE}}$ and the bremsstrahlung photon flux yields $\unicode[STIX]{x0394}k=0.31$  keV, which is approximately $20\,\%$ of the value predicted by the rule-of-thumb given in (4.4) evaluated with $Z_{\text{tot}}=18$ . One may therefore expect that during the runaway plateau following a disruption, the positron radiation is overwhelmed by X-rays more easily than during the transient avalanche phase.

As well as being distinguishable from the runaway X-rays, it is required that the total number of annihilation photons reaching a detector is sufficiently large. While this is highly sensitive to the details of the set-up, we can provide a rough estimate in the following way. The total number of annihilation photons per unit time reaching a detector with area $A_{\text{det}}$ placed a distance $R$ from the plasma detecting emission within an opening angle $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}$ , is given approximately by

(4.5) $$\begin{eqnarray}\displaystyle \frac{\unicode[STIX]{x2202}N_{\text{an}}}{\unicode[STIX]{x2202}t} & {\approx} & \displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D703}\frac{A}{R}A_{\text{det}}\frac{\unicode[STIX]{x2202}n_{\text{an}}}{\unicode[STIX]{x2202}t\unicode[STIX]{x2202}\unicode[STIX]{x1D6FA}}\approx \unicode[STIX]{x0394}\unicode[STIX]{x1D703}A_{\text{det}}AZ_{\text{tot}}\frac{n_{\text{RE}}n_{e}\unicode[STIX]{x1D702}\langle v\unicode[STIX]{x1D70E}_{c}^{+}\rangle _{\text{RE}}}{4\unicode[STIX]{x03C0}R}\nonumber\\ \displaystyle & {\approx} & \displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D703}A_{\text{det}}\frac{An_{\text{RE}}ec}{e}Z_{\text{tot}}\frac{n_{e}r_{0}^{2}}{4\unicode[STIX]{x03C0}137^{2}R}\unicode[STIX]{x1D702}\frac{\unicode[STIX]{x1D6FE}_{0}-6.7}{15}\nonumber\\ \displaystyle & {\approx} & \displaystyle (1.4\times 10^{6}\,\text{s}^{-1})n_{20}I_{\text{RE}}[100\,\text{kA}]\unicode[STIX]{x0394}\unicode[STIX]{x1D703}\frac{A_{\text{det}}[\text{dm}^{2}]}{R[\text{m}]}Z_{\text{tot}}(\unicode[STIX]{x1D6FE}_{0}-6.7)\unicode[STIX]{x1D702}.\end{eqnarray}$$

Here, the cross-sectional area $A$ of the plasma is assumed to be completely within the detector field-of-view. Then, discharges with higher plasma charge, background density and runaway current are seen to yield higher total annihilation photon fluxes. Note that a strong decrease in total photon flux is found when the electric field increases above the threshold value $E_{c}$ , due to the change in the thermalization fraction $\unicode[STIX]{x1D702}$ . As an example, inserting values typical of a disruption in a medium-sized tokamak with $R=1.5$  m, $Z_{\text{tot}}=10$ , $n_{\text{tot}}=10^{20}~\text{m}^{-3}$ , $I_{\text{RE}}=400$  kA, $E=2E_{c}$ , $A_{\text{det}}=1~\text{dm}^{2}$ , $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}=0.5$  rad and with $\ln \unicode[STIX]{x1D6EC}=15$ , one obtains a detected 511 keV annihilation photon count of $\unicode[STIX]{x2202}N_{\text{an}}/\unicode[STIX]{x2202}t\approx 7\times 10^{8}~\text{s}^{-1}$ .

In poor signal-to-noise ratio cases coincidence measurements can be employed, where only positrons annihilated between two detectors are counted. This can be approximately accounted for in the previous formula by using an opening angle $\unicode[STIX]{x0394}\unicode[STIX]{x1D703}=\sqrt{A_{\text{det}}}/R$ if two identical detectors are placed on either side of the plasma, which reduces the number of counts by another factor $0.1\,\sqrt{A_{\text{det}}}[\text{dm}]/R[\text{m}]$ .