Effect of partially ionized impurities and radiation on the effective critical electric field for runaway generation

We derive a formula for the effective critical electric field for runaway generation and decay that accounts for the presence of partially ionized impurities in combination with synchrotron and bremsstrahlung radiation losses. We show that the effective critical field is drastically larger than the classical Connor-Hastie field, and even exceeds the value obtained by replacing the free electron density by the total electron density (including both free and bound electrons). Using a kinetic equation solver with an inductive electric field, we show that the runaway current decay after an impurity injection is expected to be linear in time and proportional to the effective critical electric field in highly inductive tokamak devices. This is relevant for the efficacy of mitigation strategies for runaway electrons since it reduces the required amount of injected impurities to achieve a certain current decay rate.


Introduction
When a plasma carrying a large electric current is suddenly cooled, as happens in tokamak disruptions, a large toroidal electric field is induced due to the dramatic increase of the plasma resistivity. If this electric field is larger than a certain critical electric field, a relativistic runaway electron beam can be generated [1,2]. Such runaway beams can damage the plasma facing components on impact due to localized energy deposition. Therefore, runaway electrons constitute a significant threat to large tokamak experiments (e.g. ITER) [3][4][5].
To minimize the risk of damage, it is crucial to understand the runaway electron dynamics. Disruption mitigation by material injection is motivated by the strong influence of partially ionized atoms, as observed in experiments [3,4,6]. It is therefore important to have accurate models of the interaction between fast electrons and the partially screened nuclei of heavy ions. Fast electrons are not simply deflected by the Coulomb interaction with the net charge of the ion, but probe its internal electron structure, so that the nuclear charge is not completely screened. Energetic electrons can therefore be expected to experience higher collision rates against partially ionized impurities compared to a fully ionized plasma with the same effective charge, leading to a more efficient damping. There has been a considerable effort to produce a detailed theoretical description of this process [7][8][9][10].
A recent paper presented a generalized collision operator which describes the interaction between fast electrons and partially screened impurities via analytic modifications to the collision frequencies [9]. The elastic electron-ion collisions were modeled quantum-mechanically in the Born approximation as in [7,8], however, to obtain the required electron density distribution of the impurity ions [7,8] used the Thomas-Fermi model. In [9] we used fitted results from density functional theory (DFT) thereby providing a more accurate description. To describe inelastic collisions with bound electrons, we employed Betheʼs theory for the collisional stopping power [11], with mean ionization energies for ions calculated in [12]. Our results show that, already at subrelativistic electron energies, the deflection and slowing-down frequencies are increased significantly compared to standard collisional theory [9].
The quantity that is arguably the most important for runaway generation and decay is the threshold, or critical, electric field, which in a fully ionized plasma without radiation losses is given by the Connor-Hastie field E c = n e mc ln 4 e 3 0 2 e 2  p L ( ) [2], where n e and m e are the electron density and mass, ln L is the Coulomb logarithm, ò 0 is the vacuum permittivity and c is the speed of light. Below the threshold field no new runaway electrons are produced and all preexisting runaways eventually thermalize. There is a wealth of experimental evidence that the critical electric field is much higher than E c given above [13][14][15][16][17][18]. Well-diagnosed and reproducible experiments in quiescent plasmas on a wide range of tokamaks show that measured threshold electric fields can be approximately an order of magnitude higher than predicted by the Connor-Hastie threshold [13,18]. Furthermore, it has been shown that the runaway electron current decays much faster after high-Z particle injection than expected from conventional theory [2], in contrast to low-Z particle injection which results in a current decay rate only slightly below that expected [14]. From a theoretical point of view, the threshold electric field is expected to be higher than E c , as can be influenced by synchrotron [19,20] and bremsstrahlung radiation losses, and also, as we will show here, by the presence of partially ionized atoms. The value of the critical electric field is not only interesting theoretically-it is of immense practical importance as it determines the amount of material that has to be injected in disruption mitigation schemes [21].
In this paper we derive an analytical expression for the effective critical field for runaway generation and decay that takes into account the presence of partially screened impurities, using the generalized collision operator derived in [9]. We present a formula that accounts for arbitrary ion species in combination with synchrotron and bremsstrahlung losses. We show that the effect of partially screened impurities is captured by replacing the plasma density in the critical electric field with an effective density n n n free bound k = + , where κ is typically in the range 1-2 which implies that the effect of bound electrons is significantly larger than suggested by previous studies [22]. Furthermore, using a kinetic equation solver with a 0 D inductive electric field, we verify the prediction from [21], that the runaway current in highly inductive tokamak devices after impurity injection will decay linearly with time at a rate proportional to the effective electric field. We expect these findings will facilitate future comparisons with experimental observations of runaway-current decay, however such analysis is beyond the scope of the present paper.
The structure of the paper is as follows. In section 2 we describe the kinetic model accounting for the effect of partial screening in both the generalized collision operator and the bremsstrahlung operator. Then we proceed in section 3 to derive analytical expressions for the effective critical electric field in the presence of partially ionized impurities. This calculation generalizes the results in [20], in which the critical electric field was calculated by assuming rapid pitch-angle dynamics in the Fokker-Planck equation. In contrast to [20], our study includes the effect of partially ionized impurities and bremsstrahlung losses. We demonstrate how the presence of partially screened impurities affects both synchrotron losses (through pitch-angle scattering) and bremsstrahlung (as partial screening affects the bremsstrahlung cross-section). In section 4 we discuss the decay of a runaway current when heavy impurities are injected. Through kinetic simulations, we demonstrate the accuracy of the analytical expressions for the effective critical electric field and the current decay. Finally in section 5 we summarize our conclusions.

Kinetic equation including partially screened impurities
In a uniform, magnetized plasma, the kinetic equation for relativistic electrons can be written as follows: where f is the electron distribution function, C f FP { } is the partially screened Fokker-Planck collision operator as described in section 2.1, which accounts for ionizing as well as elastic collisions. The avalanche source is denoted S ava and E is the component of the electric field which is antiparallel to the magnetic field B. Radiation losses are modeled by C br (the bremsstrahlung collision operator) and F syn (the synchrotron radiation reaction force), which are described in section 2.2. The normalized momentum is defined as p=γv/c is (with γ the Lorentz factor), Here, T eV is the temperature in eV, n e20 is normalized to 10 m    n taking into account partial screening are given in [9].
Focusing on the effective critical electric field E c eff in this paper, the following equations are specialized to the superthermal momentum region, in which the critical momentum p c corresponding to E c eff is found. Thus all of the following expressions are given for superthermal electrons.
The generalized deflection frequency is, in units of c = å , where n j is the density of species j, and n e represents the density of free electrons. The parameter a j was determined from DFT calculations, and is an effective ion size which depends on the ion species j. These constants are given for argon and neon in table A1 in appendix A. In(5), we have assumed p a 1 10 Figure 1(a) shows the enhancement of the deflection frequency for singly ionized argon and neon. At typical runaway energies in the MeV range, the enhancement is more than an order of magnitude compared to taking the limit of complete screening and neglecting the variation of the Coulomb logarithm, which would give and I j is the mean excitation energy of the ion, normalized to the electron rest energy [12]; see table A1 in appendix A. As s n given in(9) is based on the Bethe stopping power formula matched to the low-energy asymptote corresponding to complete screening, we refer to it as the Bethe-like model. As shown in figure 1(b), the slowingdown frequency is enhanced significantly compared to the completely screened limit with constant Coulomb logarithm, where 1 s n = . The enhancement is also significantly different from a widely used rule of thumb that is mentioned in passing by Rosenbluth and Putvinski [22], which suggests that inelastic collisions with bound electrons can be taken into account by adding half the number of bound electrons to the number of free electrons. As shown in figure 1, the Rosenbluth-Putvinski (RP) model overestimates the slowing-down frequency at low energies and is a significant underestimation at high runaway energies. The weak energy-dependence of the RP model is due to the energy-dependence in the electron-electron Coulomb logarithm in(3).
In the ultra-relativistic limit p?1, the slowing-down frequency (9) is approximately p ln , 10

Radiation losses
At the high densities typical of post-disruption scenarios, bremsstrahlung may be an important energy loss mechanism compared to synchrotron radiation reaction [25,26]. In a fully ionized plasma, the required density for bremsstrahlung dominance is [27]  In a partially ionized plasma, both bremsstrahlung and synchrotron losses will be enhanced, the latter through the increased pitch-angle scattering. Both radiative energy loss channels can therefore be significant at densities characteristic of disruptions and are included in this paper.
The synchrotron radiation reaction force is given by [28,29] We model partially screened bremsstrahlung with a Boltzmann operator as presented in [26], using the model that neglects the angular deflection due to the bremsstrahlung process: ) is the normalized cross-section for an incident electron with momentum p 1 to end up with momentum p after emitting a bremsstrahlung photon carrying the energy difference, and σ br is the total bremsstrahlung cross-section for an incident electron of momentum p. The integration is taken over k p 1 , where, following [26], photon energies are cut off at 0.1% of the kinetic energy of the outgoing electrons in order to resolve the infrared divergence, i.e.k c =(γ−1)/1000. The partially screened bremsstrahlung cross-section is given in [30,31]: where k is the photon momentum and q 0 =p 1 −p−k. We use the form factor F(q) for partially ionized atoms presented in [9], In order to get an analytically tractable problem when deriving the effective critical electric field, a simplified bremsstrahlung mean-force stopping power will be used in section 3. Although a mean-force model has been shown to significantly alter the steady-state electron distribution compared to the full Boltzmann model, it captures the mean energy accurately [26], and is therefore sufficient for the purpose of deriving the effective critical electric field. This assumption is verified with numerical calculations using the full Boltzmann operator in section 4.
For the mean-force model, we have where the bremsstrahlung mean-force is given by F p , the integral taken over all allowed outgoing momenta p 1 . For argon and neon, a numerical investigation of (16) shows that F br is well approximated by

Effective critical electric field
The critical electric field is a central parameter for both generation of a runaway current and for its decay rate in a highly inductive tokamak; in the latter case, it is predicted that once the Ohmic current has dissipated, the induced electric field will be close to the critical electric field so that the current decays according to I t , where L∼μ 0 R is the self-inductance and R is the major radius of the tokamak. The physical argument is that the runaway avalanche time scale is much faster than the inductive time scale, and therefore the electric field must be close to the critical electric field to prevent rapid current variations.
We calculate the effective electrical field due to collisions with partially screened ions by finding the minimum electric field E c eff that satisfies the pitch-angle averaged force balance equation where F denotes the collisional and radiation forces on a runaway electron.
In order to find E c eff , we assume rapid pitch-angle dynamics compared to the time scale of the energy dynamics [20,32]. In the kinetic equation (1), this amounts to requiring that the pitch-angle flux vanishes. Since 1 syn 1 t - from (15), we can neglect the effect of radiation on the pitch-angle distribution (term marked as 'neglect' below) as well as the effect of the avalanche source, which is slower than both pitch-angle scattering and collisional friction. We demonstrate the validity of these assumptions in appendix B by comparing the resulting critical electric field and angular distribution to kinetic simulations. Inserting the collision frequencies (6) and (10) as well as the radiation terms (14) and (17), the kinetic equation (1) can be rewritten Following the method and notation of [20], the condition that the pitch-angle flux vanishes yields the following form for the angular distribution: where the parameter A is defined as Then,(19) integrated over pitch-angle yields a continuity equation As the sign of U(p) determines if the distribution at p is accelerated or decelerated, the effective critical electric field is the minimum electric field for which force balance is possible: The minimum can be found analytically if A?1 (so that   ). This underestimation of the effective critical field by the RP model is a result of using a simplistic form of the inelastic collision rate as well as neglecting the effect of pitch-angle scattering and radiation losses. To further explore the scaling of E c eff with magnetic field strength and impurity content, we approximate (23) in the case where one weakly ionized state j dominates: The screening constant S j is given for all argon and neon species in table A1 in appendix A. For typical magnetic fields, the terms inside the brackets tend to be roughly 1-2 times ln c L . As n N n n tot . This agrees with the scaling found in [20] for the fully ionized case. In contrast, for low magnetic fields, bremsstrahlung can increase the effective critical field by up to 20% for argon. This number is insensitive to the plasma density and depends only on its ionic composition.

Current decay
The critical electric field, especially as modified by the effects of partially screened nuclei and radiation losses, plays an important role during the relaxation of runaway electrons. In this section, we demonstrate with kinetic simulations that (23) well characterizes the threshold between runaway growth and decay under these modifications. Then, when the electric field evolves self-consistently, we show that it remains tied to E c eff under certain assumptions during the current decay phase of a tokamak disruption. If the current is carried by runaway electrons and the shape of the runaway distribution is constant in time, the time derivative of the current is related to the steady-state runaway growth rate E n n t I I t The scaling of the growth rate with impurity content may be estimated from the RP formula [22] by replacing E c with E c eff and the density by the total electron density due to the fact that bound and free electrons have equal probability of becoming runaway electrons through knock-on collisions: ) . The qualitative scaling of the analytic growth rate is confirmed in figure 4, where the growth rate is numerically calculated using CODE [34,35], which directly solves the kinetic equation (1). These simulations employed the general field-particle knock-on operator of [36][37][38] and a Boltzmann operator for partially screened bremsstrahlung losses as described in section 2.2. The vertical lines denote the analytic prediction in(23) for when one would expect the transition between growth and decay of an existing runaway population. Radiation losses affect where this threshold lies and the analytic model E c eff accurately and robustly captures this effect. In particular, we note that the mean-force bremsstrahlung model employed in the analytical derivation of E c eff agrees with the Boltzmann-type bremsstrahlung operator used in the simulations within a few percent.
The electric field is hypothesized to remain close to E c eff during the current decay phase of a tokamak disruption [21].
The mechanism by which this occurs is the fast time scale of the avalanche generation in relation to the inductive time scale of the system. A toroidal electric field is induced when there is a time-changing magnetic flux through a current loop such as a runaway beam. This magnetic flux is proportional to the total current through the loop. The induced electric field is therefore related to the rate of change of the current: where R is the major radius of the tokamak. This inductance model has recently been implemented in CODE to calculate the electric field self-consistently with the evolution of the electron velocity distribution. In general, the exact value of the inductance L will depend on the spatial distribution of current, which will change in time. For a large-aspect ratio current loop (such as a runaway beam), L can be approximated by [39] L R R a l ln Here, R is the major radius of the tokamak, a is the radius of the runaway beam, and l i parametrizes the distribution of current within the beam. We have chosen l i =1.5 as a representative mid-plateau value, based on experimental results from European medium sized tokamaks.
When E E c eff » , the growth rate can be expanded according to which allows(31) to be solved analytically: This yields a condition under which the electric field remains close to E c eff : With the estimate of E G¢( ) from the numerical results of figure 4 (at B=0 T) and estimating R/a≈5 we find that the minimum required current for E E c eff » is approximately This value is substantially lower than the estimation of 250 kA in [21], which did not include the effect of partial screening or radiation losses. Since this threshold current is inversely proportional to the inductance, the estimate(34) is only weakly dependent on the details of the spatial current distribution. Therefore, the exact value of the instantaneous inductance does not affect the primary result of this section: for large enough inductance, the electric field remains approximately tied to E c eff during the current decay phase, leading to a predictable decay time scale. To test the hypothesis that E E c eff » when I RE ?60 kA, we generate a forward-beamed initial distribution obtained from a simulation with a large electric field; the initial average runaway energy in our simulation is 17.2 MeV. We then inject singly ionized argon with a density that is four times the deuterium density n 10 m . As in the growth rate simulations, we include both synchrotron losses, the full bremsstrahlung model and a Chiu-Harvey type avalanche operator. Figure 5(a) shows the current decay, which is linear (as expected) and faster in the low-inductance case. Figure 5(b) shows the electric field evolution. Clearly, in the highinductance case, the electric field is close to the critical field after an initial transient. This means that, in highly inductive devices such as ITER, the current decay is to a very good approximation given by . Enhanced E c eff will lead to faster current decay, and(23) quantifies how fast the decay is.
On the other hand, the induced electric field deviates by approximately 10% from E c eff in the low-inductance case. Since the initial current I 0 RE =750 kA is high in relation to many medium sized tokamak experiments, E E c eff » gives an overestimation of the current decay rate in many of today's devices. The relative deviation from E c eff observed in figure 5 (33) and (34). Although the predicted induced electric field obeys E E c eff  with our assumptions, several effects could lead to a higher induced electric field in an actual experimental discharge. For example, a stronger electric field would be necessary to balance a runaway population with sub-relativistic energy, in which case the steady-state growth rate used here is inaccurate. Other effects such as transport [40][41][42], trapping [22,43] and wave-particle interaction [10,[44][45][46] may also increase the runaway current decay rate and accordingly the induced electric field. Such complete modeling remains the subject of future work. Nevertheless, partial screening has a major effect on the critical electric field as demonstrated here, and therefore the results derived herein should be an important piece toward improved experimental comparison of the runaway current decay rate as well as the avalanche growth rate.
Finally, we note that the simulations with an inductive electric field validate the initial assumption of rapid pitchangle dynamics in (19); we find that the resulting pitch-angle distribution in (20) is accurate for E E ; c eff » see appendix B. The distribution function in (20) is consequently appropriate for determining the effective critical electric field, but not for describing runaway generation.

Conclusion
Recent experimental studies on several tokamaks show that the onset and decay of runaway electrons occurs for critical electric fields that are considerably higher than the Connor-Hastie field E c . One reason is that there are other runaway loss mechanisms in addition to damping due to collisions in a fully ionized plasma that seem to dominate both in disruptive and quiescent cases. In this paper, we show that if there are heavy partially ionized impurities present in the plasma, the dominant effect on the critical electric field is the effect of partial screening. The effective critical field is further increased due to the enhanced radiation loss rates when partially ionized impurities are present.
We give analytical formulas for the effective critical electric field E c eff including partial screening and radiation effects, derived under the condition of rapid pitch-angle dynamics. The validity of this assumption and the value of the effective critical electric field is demonstrated by numerical simulations with the kinetic equation solver CODE. The most complete expression for the critical electric field is given in (23). It has been shown to be valid for a wide range of magnetic fields, impurity species and plasma composition. To make the parametric dependencies more transparent, we also give an approximate expression in (26) that is valid when one weakly ionized state dominates, which is often the case in a cold post-disruption tokamak plasma. As expected, we find that in the presence of large amounts of heavy impurities, the effective critical field can be drastically higher than E c which is proportional to the density of free electrons: E c eff even exceeds the value obtained by including the total density of both free and bound electrons. In contrast to RP [22], where the effective density includes half of the bound electrons, n n n 0.5 e bound = + , our calculations show that the bound electrons are weighted by a factor of typically 1-2. This enhancement is attributed to the energydependent collisional friction, pitch-angle scattering as well as radiation losses. Bremsstahlung and synchrotron losses both increase the effective critical field, typically by tens of percent.
Using a 0 D inductive electric field we calculate the runaway current decay after impurity injection. Through kinetic simulations we confirm the accuracy of the formula for the effective critical field (23), and demonstrate that the electric field stays close to the effective critical field when the runaway current satisfies I RE ?60 kA, in which case . These findings are relevant for the efficacy of mitigation strategies for runaway electrons in tokamak devices: since the runaway current decay rate is typically 2-4 times higher than what is predicted by the RP formula, a lower quantity of assimilated material is required for successful mitigation. As screening significantly increases the critical electric field, we anticipate that this effect is of importance to include in experimental comparisons; however, accurate predictions may require the modeling of spatial effects which are not considered here. framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 under grant agreement No 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.
Appendix A. Constants for the effective electric field Table A1 summarizes the constants needed to compute the value of the effective electric field in the presence of argon and neon. The effective ion size a j is determined by DFT simulations and is related to a j in [9] through a a 2 j j a = where α≈1/137 is the fine-structure constant. The mean excitation energy I j is taken from [12]. These give S j from (27)   The simulations with an inductive electric field (figure 5) can be used to validate the initial assumption of rapid pitch-angle dynamics in (19) leading to the pitch-angle distribution in (20). Expanding f¯in Legendre polynomials we relate the predicted analytical distribution in (20) to the ratio between the zeroth and the first Legendre modes of the The ratio given in(B.1) quantifies the narrowness of the electron distribution: f f 3 0 1 0 =¯corresponds to an isotropic distribution while the f f 3 1 1 0 ¯for a narrow, beam-like distribution. Figure B1 compares the numerical value of f f 3 1 0¯a s computed in CODE in solid black line, to the analytical prediction(B.1) in dashed green line. The analytical formula accurately predicts the distribution width on the entire interval from a fully isotropic distribution at p=0 to a narrow beam for p?1. This validates our assumptions on the rapid pitch-angle dynamics in (19). In contrast, for larger electric fields (E E 5 c eff  ), we find that the distribution rather follows the formula in Fülöp et al [47], which is derived in the limit of E E c eff  .