Minority relativistic electron populations can occur in a range of complex plasmas. Of specific interest is when runaway electrons form among the presence of high-atomic-number ion species in a tokamak plasma discharge. It has been recently demonstrated that ion charge state distributions and radiation losses at low bulk electron temperatures can be dominated by relativistic electrons, even though their density is orders of magnitude lower. This was attributed to the relativistic enhancement of electron impact inelastic cross sections. In this work, we provide a closer inspection of the atomic physics underpinning this effect. We also demonstrate the consequences of runaway enhanced scattering on post-disruption tokamak fusion discharges with neon and argon impurities present. Effects on charge state distributions, radiation and spectral characteristics, and reduced-order modeling considerations are discussed.

## I. INTRODUCTION

For safe and viable fusion power production in any future tokamak device, the fusion community must address the outstanding issue of major disruption events.^{1,2} During a major disruption in a tokamak reactor, the rapid release of significant amounts of energy has the potential to cause massive damage to plasma-facing components of the containment vessel and to also exert unbearably high electromagnetic forces on the structure itself leading to a failure.^{3,4} Additionally, the formation of runaway electrons (REs) poses another threat in the event that runaway beams impact plasma-facing components inside the tokamak.^{1,5} While impact of runaway beams with low-energy densities can be tolerated to a point, the event of a high-energy runaway beam, carrying even a small fraction of the pre-disruption plasma current, impacting the plasma-facing components locally could be catastrophic. This type of runaway beam impact could result in surface damage and penetration through the vessel to cooling, tritium breeding, and instrumentation systems.

Ideally, disruption events would be avoided in tokamak operation, which is indeed an active avenue of current research,^{2} but in reality, there must be disruption mitigation strategies to minimize the impact of losing plasma stability and to salvage the vessel in the worst-case scenario.^{1–3} One of the primary disruption strategies being explored is the deposition of impurity material into the tokamak discharge following detection of a disruption event. In this mitigation strategy, a variety of materials have been proposed, such as argon, neon, boron, and hydrogen/deuterium, in a number of different mixtures and formats, such as gas, frozen pellets, dust, or encapsulated dust inside metallic shells.^{1,2,6,7}

Injection of these impurities can provide thermal quench mitigation by dissipating the plasma power load over the reactor first wall through radiative cooling of the thermal plasma.^{8–10} In the current quench phase, where plasma current is a mixture of thermal and RE currents, secondary impurity injection can provide another means for increasing the current quench rate.^{5,9,11}

In many of the proposed impurity injection strategies, a common outcome is a resulting plasma mixture of hydrogenic fuel, helium ash, and high Z impurities, with a total atomic number density on the order $1020\u221221$ m^{−3}. This ultimately leads to a radiatively cooled bulk electron population, accompanied by a minority RE component with a number density of approximately $1016\u221217$ $m\u22123,$ which is quite sufficient to carry a plasma current of a few to over 10 MA in potential future ITER discharges. While the current study of this complex type of plasma is indeed for tokamak fusion contexts, similar plasmas have been observed in a range of non-fusion examples such as planetary radiation belts^{12} and atmospheric lightning.^{13}

Some recent studies have begun exploring modeling approaches to account for high Z impurities and a minority RE component being present in a post-disruption plasma.^{14,15} In a mitigated argon-doped tokamak post-disruption discharge where $nimpurity\u2248nD=1020$ $m\u22123,$ it was shown via collisional-radiative (CR) modeling that RE-driven atomic processes can have a considerable impact in modulating ion charge state balance and radiation through line emission at lower plasma temperatures, which are to be expected following thermal quench and collapse of the bulk electron temperature to on the order of 1–10 eV.

The underlying cause for the observed RE influence on macroscopic plasma properties was discussed to be a quantum electrodynamics (QED) effect dating back to the seminal works of Møller, Breit, and Bethe.^{16–18} In conventional plasma modeling approaches, it is normally assumed that at electron energies significantly above the threshold energy, the excitation and ionization integrated cross section (ICS) used as fundamental input data decrease with increasing electron impact energy. However, this trend in fact reverses when relativistic energies are reached. The observed enhanced scattering cross section is a lowest-order QED correction associated with the Møller or generalized Breit interaction (GBI),^{19–23} which will be expanded on later in this work. As a very simple illustration of this, we refer to Fig. 1, which demonstrates the overlapping energy window in which REs are expected to populate, but also where enhanced inelastic scattering is observed. It is the interaction of these two relativistic features which we seek to probe and understand in this study.

While there are many established CR models used to model magnetic confinement fusion (MCF) plasmas,^{24–26} at this time there are a few that account for the enhanced inelastic scattering of relativistic electrons. Inclusion of such scattering is crucial to understanding and modeling the impact of minority runaway electron populations that may form in disruption scenarios. Further, the requirement for a CR model to span a wide range of thermal temperatures, such as $100\u2264Te\u2264104$ eV in tokamak disruption scenarios, is not a straightforward one, as CR codes are generally tailored to model either low temperature or atmospheric discharges^{27,28} over a limited range of low plasma temperatures, often with detailed excited-state level descriptions, or at the other end of the spectrum, highly ionized, and dense plasmas such as in high-energy density physics,^{29,30} where approximate level descriptions are often employed. As such, we have sought to produce a CR model, which can account for the necessary physics and plasma parameters of potential post-disruption tokamak discharges in a computationally efficient manner.

Given the progress made in forming an improved CR model, we will be able to demonstrate the qualitative and quantitative differences that arise when treating runaway and thermal bulk electrons through a coupled CR model. This will be contrasted against a common method of treating runaway and thermal electrons in a decoupled manner,^{31,32} where generally only the background Maxwellian population at given temperature is used to determine ionization balance and radiation losses, such as line emission. Additionally, we will address differences between a CR treatment and a reduced-physics model, often termed a coronal-like model,^{24,33,34} and what considerations practitioners should take into account when developing models of similar dilute MCF discharges.

In this paper, we will first outline the motivation for our recently formulated CR model in Sec. II, developed as a fork of the widely used FLYCHK CR code.^{29} We then demonstrate electron scattering benchmarks in Sec. III to help justify the application of this model to near-neutral targets of importance to post-thermal-quench tokamak plasmas, where significant densities of added impurity material may be present. We then outline the conditions and setup for CR simulations of post-thermal quench plasmas with added impurities in Sec. IV. This is followed by presentation and analysis of the results from these simulations in Sec. V A where we demonstrate the effect of minority RE populations on the charge state distributions (CSDs), radiation losses, and spectra characteristics for neon and argon. These runaway-driven modifications to charge state distributions are analyzed to highlight impacts on electron collision models that find use in modern kinetic and particle-based plasma simulations to highlight the importance of suitably accounting for relativistic effects in tokamak discharge CR modeling. Finally, we summarize our observations in Sec. VI and present an outlook on where to extend and implement our findings toward improving simulation models of the vital problem of tokamak disruption mitigation.

## II. COLLISIONAL RADIATIVE MODEL

CR modeling has been used to computationally model plasma composition, radiation, and spectroscopy for many years.^{35–38} Many detailed and thorough reviews and texts on modern CR modeling can be found in the literature,^{20,29,38–40} and as such, here we simply present a brief summary of the mechanics of how to describe and implement the CR model used in this work. The standard CR model can be written as

where ** n** is the state vector of excited-state populations of each ion stage being considered for the target atom(s), and $R$ is the CR rate matrix where each element

*R*is the sum of transition rates between state

_{ij}*i*to state

*j*due to all considered collisional and radiative processes in the plasma, such as electron impact (de-)excitation and ionization (recombination), or spontaneous radiative decay for example.

Depending on the atomic number of the target being studied, *Z*, and the level of refinement chosen for excited-state energy levels, the size of the state vector, ** n**, can vary massively between roughly $O(101)\u2212O(109)$. As such, the trade-off between physics accuracy in the energy level refinement and the computational time needed to solve the system must be considered for the application at hand. First, we may recall common energy level labeling via the following quantum numbers: principal quantum number

*n*, orbital angular momentum quantum number

*l*, and total angular momentum quantum number

*j*. Typical energy level refinement schemes used in CR models range from fine structure (

*n*,

*l*,

*j*) resolved, to configuration averaged (

*n*,

*l*) resolved, and finally averaging over

*l*to obtain superconfiguration (

*n*) resolved.

For this work, we continue the development of a fork of the popular FLYCHK CR model,^{14,29} albeit with refinements targeted for use in the MCF plasmas of interest. In this CR model formulation, a superconfiguration assumption is used to describe excited states with a principal quantum number (*n*). The FLYCHK model has been available online with hosting and support from the US National Institute of Standards and Technology^{41} and has found wide use among experimental and computational physicists. This code has been widely used in high-energy dense plasma applications,^{29} due to its demonstrated applicability of a screened hydrogenic superconfiguration model to replicate essential physics in a computationally efficient manner.^{30,42} In development of our compact CR model, we also make use of the Los Alamos suite of atomic physics codes^{20} to calibrate and verify our CR model. This is of particular benefit as a source of guidance and benchmarking of our model in sufficiently ionized plasmas, where the LANL atomic physics codes represent a state-of-the-art CR modeling suite that has been extensively benchmarked over the last decade.^{20,30,42,43}

A particular benchmarking example performed during development^{14} of the present CR model was to verify the approximate treatment for the contribution from $\Delta n=0$ transitions in line emissions. Given that a superconfiguration description does not naturally account for possible $\Delta n=0$ transitions, which can have significant contributions to the radiative power loss (RPL) via transitions within Layzer complexes of mid- to high-*Z* targets,^{44,45} a correction should be made to account for these potentially large emission lines. In particular, for conditions relevant to magnetically confined fusion plasmas (i.e., relatively low *n _{e}*), it has been shown that contributions from the $\Delta n=0$ lines of intermediate Z ions, such as neon or argon, between approximately $5\u2264Te\u226450$ eV can be equal to, or greater than, RPL contributions by $\Delta n>0$ lines.

^{34,38,46}To approximate this contribution, a partial local thermodynamic equilibrium (pLTE) assumption is made on the excited-state populations that would correspond to the $\Delta n=0$ fine structure transitions within each ion stage.

^{34}With these population estimates, and pre-computed transition energies and emission rates,

^{47}an approximation of the line power due to $\Delta n=0$ bound–bound transitions within a given ion stage can be computed.

^{29}The assumption of pLTE for transitions within the same principal quantum is grounded in the observation that electron impact excitation cross sections for $\Delta n=0$ transitions are significantly higher than cross sections for $\Delta n>0$ transitions. The collisional dominance of the $\Delta n=0$ transition pushes these excited states to an effective LTE at a lower temperature than the bulk

*T*, as the spontaneous radiative de-excitation is less dominant than collisional de-excitation.

_{e}^{34,38}The effect of this approximation for $\Delta n=0$ line emission was compared with thermal plasma results of Fournier

*et al.*

^{46}and the Los Alamos suite of codes' fine structure CR model

^{20}to ensure the approximations were indeed reasonable.

While FLYCHK has found wide use in applications involving dense, highly ionized plasmas, its applicability to near-neutral plasmas does require some care. One reason for this is the use of a screened hydrogenic approximation (SHA) to describe superconfiguration energy levels of excited states in near-neutral targets, where excited states are distinguished by (*n*) only. While accurate binding energies are used for ground-state levels of ion stages, we note that the SHA is indeed but an estimate when compared to results of detailed atomic physics calculations, such as HULLAC,^{47} FAC,^{48} or ATOMIC.^{20} While there are shortfalls of averaging over energy levels, it does vastly reduce the size of the rate-matrix system and thus the computational expense to study physics effects. Further, it has been found suitable for modeling plasmas when one is primarily concerned with macroscopic quantities, such as charge state populations and total radiative cooling power,^{49} especially at the plasma densities of interest to MCF.^{50} While noting this shortcoming of the model, we believe the general qualitative trends observed and reported in this study are indeed valid. The trade-off between fine spectroscopic accuracy vs computational tractability and ease of studying new physics effects is a balance that we seek to maintain in this work and subsequent studies.

While out of the scope of the current study, in which we are focusing on exploring the effects of relativistic scattering on CR modeling, refinement of the (*n*) resolved superconfiguration energy levels of selected near-neutral targets to (*n*, *l*) resolved configuration averaged levels may be explored in the future. This would have precedent in the way FLYCHK originally employed more accurate energy levels for Li-like, He-like, and H-like ions to allow more accurate K-shell spectral analysis.^{29}

## III. ELECTRON SCATTERING INPUT DATA

While the overarching machinery of the superconfiguration CR model employed in this work is indeed similar to that contained in FLYCHK^{29} and SCFLY,^{51} we now note some specific differences between those and the model discussed in this present study. Essentially, the key differences in the current model are based on prescribing fit-for-purpose atomic inputs for the fusion-relevant problems under consideration. These inputs appear as inelastic electron impact cross sections employed in our model. These changes were made to address physics requirements for modeling the impact of relativistic electron populations in relatively cold post-thermal-quench tokamak discharges. We first address some scaling for neutral and near-neutral targets, followed by relativistic scattering considerations.

### A. Neutral and near-neutral inelastic electron scattering

One aspect of the original FLYCHK model, which may be less applicable to near-neutral plasmas, is the formula used to approximate electron scattering cross sections. These approximations are necessarily employed to provide a complete set of scattering data between the possible state populations in the CR system. The assumptions invoked in order to arrive at simpler approximate formulas are often sufficiently accurate for ionized targets, where electron scattering is dominated by a stronger interaction between the positively charged nucleus and projectile electron, compared to the interaction between target electrons and the projectile electron. Conversely, for neutral or near-neutral targets, the projectile electron generally experiences a more complex scattering interaction with the target nucleus as well as the bound target electrons. To address this, we now detail modifications of the electron impact excitation and ionization cross sections employed in the original version of FLYCHK, in an effort to improve the applicability of our model to near-neutral plasmas.

Electron impact excitation in FLYCHK is approximated by the form of van Regemorter,^{52} where for an excitation from a state *i* to state *j,* the cross section is

where $U=E/\Delta Eij$ is the scaled incoming electron energy, *I _{H}* is the Rydberg energy,

*f*is the oscillator strength, and

_{ij}*a*

_{0}is the Bohr radius. The Gaunt factor,

*g*(

*U*), is given by

and coefficients for each transition are determined through empirical fits of (3) to excitation cross sections computed by Chung *et al.*^{51} via plane wave Born (PWB) calculations.

While the use of PWB excitation cross sections is well motivated for providing data for the many transitions required for CR modeling of charged targets, such as Ar^{+5} or Ne^{+6} and so on, it is known to be less accurate for describing excitation of neutral and near-neutral targets.^{22} In an effort to remedy this, we employ scaling methods for cross sections outlined by Kim^{53,54} to modify the cross sections for neutral and singly ionized targets. In the approach of Kim, PWB and Coulomb Born (CB) cross sections are used as a basis for rapidly obtaining values of the required data, and then, scaling for neutral and singly ionized targets is applied to scale and shift the peak of the cross sections, which was shown to be successful in producing results comparable to accurate calculations.^{53,54} Given that the PWB method is used to obtain the excitation cross sections employed in (2) and (3), it is sensible to then apply the scaling rules of Kim^{53,54} in an effort to improve the applicability of our CR model for near-neutral targets.

The Binary Encounter Bethe (BEB) scaling rule^{53} used for scaling a neutral target electron impact excitation cross section, *σ _{ij}*, from (2), for a transition with threshold energy $\Delta Eij$, is given by

where $\Delta EI$ is the ionization energy of the target.

The Binary Encounter (BE) scaling rule^{54} for a singly ionized target is given by

where we note the absence of the ionization energy when compared to (4). We also note that Kim originally employed scaling of CB cross sections for singly ionized targets, not PWB cross sections. However, given the approximate nature of PWB cross sections used in the formulation of the present CR model, we find that using the prescribed scaling formula is beneficial nonetheless as a step toward improving the quality of input data.

For electron impact ionization (EII), the non-relativistic cross section is modeled by the semi-empirical expression of Burgess and Chidichimo^{55}

where $\Delta EI$ is the ionization threshold energy for the target of charge state *Z*, $U=E/\Delta EI$, *ξ* is the number of electrons in the shell being ionized, *C* is an arbitrary constant nominally around 2,^{55} and for this work is fixed at a constant value of 2.3 as is done in FLYCHK,^{29} and

where *Q _{n}* is the screened charge felt by a target electron, in the

*n*th shell, to be ejected following ionization.

^{29}

To scale the electron impact ionization cross sections of neutral and singly ionized targets, we make an empirical observation that applying the same scaling rule of Kim used for scaling excitation cross sections produces a beneficial correction to the ionization cross section of the valence shell when compared to accurate calculated results, see Sec. III C. Here, we simply replace the excitation energy found in Eq. (4) with the ionization energy such that for neutral and singly ionized targets the non-relativistic electron impact ionization cross section of valence shell ionization is scaled as

where the empirical cross section of Burgess and Chidichimo $\sigma IBC$ is given by Eq. (6).

We should note that while the original intention of the scaling formulas was to correct excitation cross sections, and not for ionization as we have done here, this style of scaling method has been previously demonstrated to be useful.^{56} This can be further understood considering that when Burgess and Chidichimo arrived at the empirical formula of Eq. (6), it was designed for application to doubly charged ions and above, and then benchmarked against ionization cross sections calculated under the CB approximation for a binary encounter.^{57} Similar to the PWB approach, this approximation is also well suited to charged targets, but is less applicable to near-neutral targets.^{58} Since dipole allowed excitations to continuum states are the dominant contribution to the total ionization cross section, we argue that the same scaling correction of PWB cross sections, which was designed to correct dipole allowed excitations,^{53} is also reasonably justified for scaling ionization cross sections grounded in the binary encounter limit, as is the case for Eq. (6).

We reiterate that the empirical scaling of ionization cross sections in Eq. (9) was only applied to ionization of valence shell electrons, and all inner shell electron impact ionization transitions remain governed by Eq. (6).

For both excitation and ionization cross sections discussed in this section, results of benchmarking the proposed scalings for near-neutral targets are shown in Sec. III C for key targets relevant to tokamak plasma studies.

### B. Relativistic inelastic scattering

A further requirement of our model is that we must accommodate for scattering of a minority RE population present at relativistic energies, as summarized in Garland *et al.*^{14} The physical origins of enhanced inelastic scattering at relativistic energies can be understood when considering a QED treatment of the scattering problem,^{59} accounting for relativistic effects. In doing so, two corrections are made to the familiar Coulomb interaction potential, which is otherwise well suited to describing electron–electron scattering in a non-relativistic case. These corrections account for a magnetic potential and a retarded potential that become non-negligible in relativistic scattering. The magnetic correction can be understood as the interaction of the target electron's magnetic dipole moment with the dynamic magnetic field produced by the moving projectile electron. The retardation correction is applied to account for the fact that the Coulomb interaction between two electrons does not occur instantaneously, due to the finite speed of light. This effect produces a time-dependent delay in the electric force between the electrons and can become significant when describing the interaction under relativistic conditions.

It is instructive to note that there are two primary, but equivalent,^{59} approaches of how one may include the aforementioned QED corrections to more accurately compute inelastic scattering cross sections. If one expresses the lowest-order Feynman diagram in the Coulomb gauge, then both magnetic and retarded terms are included in the generalized Breit interaction^{18,60} (GBI), which must be then added to the standard Coulomb interaction. If the Lorenz gauge is used instead, all three aforementioned interactions are described in one expression known as the Møller interaction.^{17}

In our approach, we approximate the enhanced inelastic scattering by employing a Møller–Bethe-like functional form^{21,23,61} at high electron kinetic energies. For non-relativistic energies, the existing electron impact excitation and ionization cross-sectional formulas employed in FLYCHK are used, with some modification as will be shortly discussed.

Using (2) as a basis, a Møller–Bethe-like relativistic correction^{16,17,61} is applied, so the total excitation ICS is

where

and the sigmoid smoothing function is defined as

with *E* normalized with respect to 1 eV. The form of (12) is chosen to adequately describe known relativistic ICS for inelastic processes.

In a similar manner as was done for excitation, we use previously discussed non-relativistic ionization cross sections (6) as a basis to form a total ionization cross section

augmented by the relativistic correction^{21,62}

Figures 2 and 4 compare the proposed fits with calculated ICS for excitation and ionization in the literature. We note agreement within a factor of two of recent R-matrix theories at the peak of the cross section and that, unlike those reported by Wang *et al.*,^{63} normally such calculations are not extended into the relativistic regime of relevance here.

### C. Electron impact cross-sectional benchmarking

In Secs. III A and III B, modifications to the base electron impact excitation and ionization cross sections contained in the FLYCHK model were outlined: first to correct near-threshold behavior for neutral and singly ionized targets, and second to add a cross section component to accommodate enhanced inelastic scattering at relativistic electron energies. Here, we now demonstrate some results of our benchmarking exercises which we routinely perform to ensure validity of the approximations and assumptions invoked in our cross sections required for input to the CR model.

We first show an example of electron impact excitation of neutral neon and argon from the ground state to the first and second excited states of both targets. As a reference point for benchmarking, we compare against available data produced by B-spline R-matrix (BSR) calculations of Bartschat and Zatsarinny^{64–67} obtained via the LXCat database.^{68} We also compare against relativistic distorted wave (RDW) calculations from the Los Alamos suite of atomic physics codes^{20} to provide an estimated upper bound on neon cross sections. We note that a mostly complete set of benchmark BSR cross sections was available for neon to make a comparison. In contrast, the available data for argon are more sparse, and missing some transitions from initial states. Effectively, this comparison then acts as a lower bound on the cross section, and one may reasonably expect the true cross section to be marginally higher than what is plotted for the argon BSR result. Similarly, the available incomplete set of transitions highlighted in Bretagne *et al.*^{61} were used to provide a similar lower bound on the cross sections for argon.

Following the presentation of excitation cross sections, we now show benchmarking for electron impact ionization of fusion relevant targets in Figs. 3–5. In Fig. 3, we first demonstrate the ability of our approximate cross-sectional formulas to describe ionization from M, L, and K shells of Al-like argon ($Ar+5)$ by compared to scaled hydrogenic results of the GIPPER ionization code contained in the Los Alamos suite of atomic physics codes.^{20} Here, we see that for a significantly ionized target there is good agreement between our present approach and the cross sections produced from a first-principles solution of the scattering problem.

For electron impact ionization of neutral targets, we demonstrate applicability of our method in Figs. 4 and 5. In Fig. 4, we show ionization of the valence shell for helium, nitrogen, neon, and argon—which are all relevant atoms in various scenarios of the intended operation of ITER. It is evident to see the benefit of scaling the empirical formula of Burgess and Chidichimo^{55} so that our present formula can better replicate theory and experimental results.

For the conditions relevant to tokamak plasmas, the majority of electron impact ionization is generally considered to come via valence shell ionization. So far our studies have found that this is indeed a fair generalization, with inner shell electron impact ionization of neon and argon contributing approximately a few percent of total ionization rates for neon and up to 5%–10% for argon. For the purposes of completeness, and to enable application of the model presented in this study to other plasma conditions, we also benchmark inner shell electron impact ionization for targets of interested. In contrast to ionization of electrons in the valence shell of a neutral or near-neutral target, it was found that inner shell ionization cross sections using the empirical formula of Burgess and Chidichimo^{55} did not need to be scaled in the same manner to better agree with values in the literature, shown in Fig. 5. This was expected as inner-shell excitation and ionization processes have significantly less dependence on coupling to the continuum.

### D. Relativistic rate formulas

Finally, in order to thoroughly detail the process required to modify a typical CR model to account for relativistic electron populations, we now briefly highlight relativistic considerations in evaluating collision rates. The relativistically invariant definition of a typical rate, *R*, can be written^{70}

where *N _{e}* is the electron density,

*v*is the impact electron velocity,

*p*is the impact electron linear momentum, $fe(p)$ is the electron momentum distribution function, and $\sigma (E)$ is the electron impact cross section for the process being considered, generally given as a function of the electron kinetic energy,

*E*.

For convenience in implementing the CR model, we seek to rewrite (15) as a function solely of *E*, with the caveat of ensuring that a relativistically invariant expression is retained. By using the differential relation^{70,71}

and normalization condition on the distribution function

we can rewrite (15) as the very familiar expression

which is indeed relativistically invariant.

We note that while (18) is a very standard expression many modelers will have implemented, one should take care to ensure that a relativistically correct expression for the impact electron velocity,

is used when evaluating (18), in lieu of the familiar,

expression that is often used in the non-relativistic limit.

## IV. SIMULATION SETUP

Simulations presented in this study are characterized by a number of input parameters to the CR model outlined in Sec. II. First, the ion species must be specified, and in this study, we show results given some fixed density of deuterium *n _{D}* and another impurity, either neon or argon with a fixed density,

*n*.

_{I}Given a fixed atomic density of a species, *n _{I}*, with an atomic number

*Z*, and knowledge of the charge state populations, the electron density due to a given species present alongside a population of singly ionized deuterium, $nD+$, is thus

where *Z _{j}* is the charge of ion stage

*j*, and

*n*is the population density of the ion stage.

_{j}We note that for the typical MCF discharges of focus for this model, species densities of fuel hydrogen isotopes and associated injected impurities are generally on the order of $1013\u20131015$ cm^{−3}, depending on specific device operations. Additionally, the electron and ion densities considered in these contexts are generally said to be sufficiently low that radiation transport effects can be reasonably neglected.^{38} Therefore, in this study, the effects of plasma opacity are being neglected for now. It can be noted that in some MCF plasma conditions, UV lines can be trapped, while visible lines are not. Our future studies will explore this effect, but for now, we make clear the assumption of an optically thin plasma. In addition to species densities, the plasma electron energy distribution properties should be specified to allow evaluation of CR transition rates. This is now addressed in Sec. IV A.

### A. Parameterizing electron energy distribution

As noted earlier, to evaluate the transition rates between excited and/or ionic states, such as the standard rate integral (18), an electron energy distribution function (EEDF), $fe(E)$, is required. In practice, any form of $fe(E)$ can be used; however, for the purposes of this investigation, it is convenient to define a simple analytic expression where possible, to simplify the numerical integration procedures required while still retaining a reasonable distribution informed by physics of RE formation.

From Lvovskiy *et al.*,^{72} a recent indirect measurement of an experimentally produced runaway EEDF is reported for a RE beam, following the purge of higher Z impurities. In that study, the EEDF was obtained via inverting a measured hard x-ray spectra. We add a caveat that this inversion process may struggle to reconstruct a hollow shell, or bump on tail in energy space, accurately and so the downward dip at the low-energy side of the tail maximum therefore has a larger uncertainty associated with it, as shown in Fig. 6. Further, the larger uncertainty shown in measurements above 10 MeV may be attributed to detector inefficiency for those higher energies. While the measurement of runaway electron distributions found in experimental discharges is indeed still a preliminary science, the reported result indeed indicates the “bump-on-tail” form of distribution, at characteristic 1–10 MeV energies, often shown in theory and modeling studies.^{5,73–75}

Aiming to replicate the general bump-on-tail feature, we trialed and compared various analytic distributions, namely,

- Maxwellian EEDF(22)$fM(E)=2E\pi Te1.5exp\u2009(\u2212E/Te),$
for a given temperature

*T*in eV._{e} - Gaussian EEDF(23)$fG(E)=12\pi \sigma Eexp\u2009(\u2212(E\u2212\mu E)22\sigma E2),$
where

*μ*is the mean energy and_{E}*σ*is the standard deviation, which may be defined in terms of the full-width-half-maximum length via $\sigma E=EFWHM22\u2009log\u20092$._{E} - Maxwell–Jüttner EEDF(24)$fMJ(\gamma )=\gamma 2\beta \theta K2(1/\theta )exp\u2009(\u2212\gamma /\theta ),$
where for a given temperature

*T*, $\theta =Te/mec2$. Also, $\gamma =1+p2$ is the relativistic Lorentz factor, where_{e}*p*is electron momentum normalized to*m*, $\beta =1\u22121/\gamma 2$, and_{ec}*K*_{2}is the modified Bessel function of the second kind. - Boosted Maxwell–Jüttner EEDF: where a boosted (translated) Maxwell–Jüttner distribution is defined in $(p\u2225,p\u22a5)$ space as(25)$fBMJ(p,p\u2225)=14\pi \theta \gamma bK2(1/\theta )exp\u2009(\u2212(\gamma b\gamma (p)\u2212pbp\u2225)/\theta ),$
where $\gamma b=1+pb2$ and

*p*is the normalized momentum of the translation in reference frame._{b}

This distribution can be integrated over pitch given $p\u2225=p\xi $ to obtain

With these analytic expressions, we performed fits of trial EEDF forms against the reconstructed experimental EEDF to best align the RE tail peak. Results of this parameterization are shown in Fig. 6. It was found that both Maxwellian and Maxwell–Jüttner EEDFs were not able to well describe the features, while a Gaussian representation captured some of the bump-on-tail features at the maximum; however, the high-energy tail is not well described. Interestingly, a boosted Maxwell–Jüttner with a lower temperature of $Te=105$ eV and a translation in normalized momentum space of *p _{b}* = 8 appears to well describe both the bump-on-tail features and the decay of the high-energy tail.

As a result, for most parts of this study we chose to employ this boosted Maxwell–Jüttner parameterization for the runaway component of the EEDF. For the thermal component, a standard Maxwellian distribution was assumed. Given distinct analytic EEDFs, we may define the total EEDF as the addition of a majority thermal component Maxwellian distribution, with some analytic function defining the runaway tail population containing a much smaller electron density fraction

where $x=nRA/ne$ is a density fraction computed as the ratio of the desired runaway tail electron density to the total electron density, which is effectively the bulk electron density. Based on current-carrying capacity of a RE beam at low *T _{e}*, this fraction is often on the order of $x\u224810\u22124\u201310\u22123$ in order to carry a typical plasma current of present tokamaks as well as proposed future plasma currents in devices such as ITER.

## V. RESULTS AND DISCUSSION

### A. Effect of relativistic CR corrections on ionization balance and radiative power loss rate

To examine the effect of including the aforementioned relativistic QED corrections to a CR modeling approach, we present a very simple comparison. We present three sets of simulations, for both neon and argon atomic targets, using (i) a standard thermal plasma assumption with non-relativistic scattering assumed, (ii) a thermal plus runaway EEDF with non-relativistic scattering assumed, and (iii) a thermal plus runaway EEDF with relativistic scattering corrections. In this way, we can see the differences between using a somewhat typical model with a common Maxwellian EEDF assumption, using a typical model with a runaway tail added, and finally using a relativistically QED corrected CR model with a runaway tail added.

In formulating this comparison study, we aimed to think of the simplest approach one may first take to include a runaway population in a model. We concluded that was to simply add a Maxwellian distribution at a higher *T _{e}*, in addition to the thermal component as per (27). While this is indeed different from the boosted Maxwell–Jüttner EEDF discussed in Sec. IV A, our aim with this comparison is to emphasize the impact of scattering physics more so than being an exact replication of a tokamak RE scenario. For Secs. V B, D, E, and C of this study, the boosted Maxwell–Jüttner form is assumed, with relevance to more realistic scenarios.

We first present the average ion charge state, $Z\xaf$, and total radiative power losses (RPL) of neon as a function of thermal plasma temperature, *T _{e}*, in Fig. 7. These results are computed given a runaway population described by a Maxwellian with $TeRA=10$ MeV, and gas densities $nNe=nD=1014$ cm

^{−3}and $nRA=5\xd71010$ cm

^{−3}.

As seen in Fig. 7, there is an enhancement in average charge state and the associated RPL at very low *T _{e}* when only considering the addition of a runaway tail population. When the runaway tail population is accounted for, with the relativistically QED corrected scattering input to the CR model, there is a further major enhancement in $Z\xaf$ and RPL. This effect is a direct result of the relativistic electron population sampling inelastic scattering cross sections that have an increased magnitude above the electron rest mass energy.

Given the enhancement in the average charge state due to RE-driven ionization, it is worth considering the effect on not just the average charge, but also the spread and variety of ion charge states that may be found in a given plasma configuration. To examine this, we present three sample charge state distributions (CSDs) in Fig. 8, along with the standard deviation of the ion charge distribution as a function of *T _{e}*. The sample CSDs were taken at 2, 5, and 25 eV to demonstrate the variation in the ion spread at various stages.

From these figures, it is clear that the presence of runaway driven processes not only enhances a key quantity such as the average ion charge state but also increases the spread and variety of ion charge states present in a discharge at low bulk *T _{e}*. To further demonstrate this effect, we refer the reader to a supplementary animation that shows the increased spread of ion charge states with progressively lower values of

*T*(see the supplementary material). Acknowledging this effect is an important step in understanding the impact of relativistic scattering on CR modeling results, especially when considering the effects of electron scattering from different ion species for the input data to kinetic models as discussed in Sec. V D.

_{e}Following the results for neon, we now present the same quantities for an argon–deuterium discharge with the same density, $nAr=nD=1014$ cm^{−3}, and RE population properties as before. The average ion charge state, $Z\xaf$, and RPL are shown in Fig. 9, complemented by sampled CSDs and standard deviation of ion charge states in Fig. 10.

While the runaway enhancement is not as prominent as in the neon discharge scenarios, the same effects on the macroscopic plasma properties can be observed in argon. An interesting comparison can be identified between the two atomic targets considered in this study, neon and argon. A direct result of the higher binding energy of electrons in neon is that the influence of runaway-driven ionization is demonstrated over a wider window of plasma temperature, *T _{e}*. This can be considered as an inability of low

*T*thermal electron populations to excite, and eventually ionize, the neon targets, allowing the runaway-driven inelastic collisions to stand out. In contrast, the lower binding energy of argon means a cooler thermal population can ionize argon targets much more effectively than it could neon for the same

_{e}*T*. As a result, the influence of the RE population gives way to the thermal electrons at a much lower

_{e}*T*than neon.

_{e}The differences between RE interactions with impurity species are a useful comparison to make, considering the current plans for using neon in ITER's disruption mitigation systems, and the fact that most current tokamak experiments use argon in RE experimental studies. As such, any theory or model that then seeks to accurately describe the injection of these impurities should therefore be constructed with an awareness of these differences, and that what works for one target may not necessarily apply to the other.

The trends identified here naturally will be important also when considering other species of impurities potentially expected in a tokamak plasma, such as helium ash or nitrogen as a buffer gas. Given the high ionization energy of neon, 21.57 eV, we can thus expect similar CR model results for runaway interaction with helium, given the higher ionization energy, 24.59 eV, of helium. Likewise, the comparatively lesser influence of runaway-driven CR processes in argon, with an ionization energy of 15.76 eV, would be expected to persist in CR modeling of nitrogen, where the ionization energy is 14.53 eV.

In addition to the impact on macroscopic parameters, such as average charge state of ion species and the total RPL experienced by the plasma, it is useful to examine the differences in line emission spectra when the RE population is accounted for in the CR model. An interesting comparison between the spectra found in a basic thermal plasma scenario and the spectra found when runaways are included is shown in Fig. 11 for a bulk plasma temperature of *T _{e}* = 2 eV. Here, the qualitative differences are immediately clear. A few dominant spectral lines in the thermal case are replaced by a dense series of smaller but abundant lines owing to the wider variation of ion species present in the plasma.

From our observations, the qualitative structure of these spectral results is representative of impurity doped discharges at lower thermal electron temperatures, approximately $Te<10$ eV, when a current carrying RE population is present. Further, we note the somewhat isolated spectral series at higher photon energies to the right of the lower panel. These emission lines are from Ne-like Ar^{+8}, which has considerably large excitation energies owing to the closed-shell structure of the Ne-like ion.^{76} As a result of the large energy gap, a noticeable spectral line contribution can be seen even when small fractional populations are present. As such, we believe this isolated, high-energy (short wavelength) qualitative spectral signature may be a potential option for diagnosing the presence of REs in cooler post-disruption plasmas.

### B. Impact on bulk electron power balance

An intended application of the improved CR model described in this paper is to help quantify the power balance of the bulk or background electron population, which along with the runaway evolution to be discussed in Sec. V D, holds the key in setting the current quench dynamics that involve both Ohmic and runaway currents. Although our eventual approach in addressing this problem is a coupled CR and plasma transport calculation, here we give a brief discussion on the various energy transfer channels associated with bulk electron population. We are primarily concerned with a post-thermal-quench plasma that has the runaway current fully formed. The issue is the range of thermal electron temperatures that could be sustained given sources of energy and accompanying mechanisms of energy loss. In Fig. 12, we show the power loss via excitation and ionization collisions of both thermal and runaway electron populations in the prototype neon-doped plasma simulation as outlined earlier in this section prior to Fig. 7. In separating these populations, we may estimate ranges of *T _{e}* where the energy loss through inelastic collisions may be balanced by heating mechanisms acting on the thermal electron population.

For RE collisional losses, we can clearly see the dominant ionization contribution at lower *T _{e,}* which reflects the RPL observations seen in Fig. 7. This contribution naturally decreases with increasing

*T*as the thermal electron population begins to dominate. For the thermal population, we note rapid decreases of collisional losses with decreasing

_{e}*T*due to the inability of lower energy populations to overcome threshold energies required to undergo the inelastic collisions. At approximately 5 eV, we see a hand-over between thermal and RE populations in terms of which population is the dominant driver of radiative power loss.

_{e}In terms of delivering energy to a thermal electron population, we can consider heating via runaway-thermal Coulomb collisions. The energy deposition rate can be approximated in the limit of ultra-relativistic electrons^{5,75}

where *m _{e}* is the electron mass,

*c*speed of light,

*r*classical electron radius, $log\u2009\Lambda $ is free electron Coulomb logarithm,

_{e}*n*is the background free electron density, and $nRE$ is runaway electron density.

_{e}By comparing the curve produced by (28) with the thermal electron cooling mechanisms, we can see a non-negligible heating contribution below approximately *T _{e}* = 5 eV. This indicates that a practical floor in the thermal temperature exists, where a balance between heating and cooling mechanisms may be found. While this study presents long-time steady-state results, we anticipate the timescales of possible heating and cooling mechanisms to be varied, and thus, a dynamic evolution of the EEDF is necessitated to obtain an accurate description of runaway and thermal electron dynamics in these post-thermal-quench scenarios.

We also consider an additional heating source by secondary electron formation following RE driven ionization. At relativistic impact energies, the electron energy sharing probability distribution between the high-energy primary, *E _{p}*, and low-energy progeny electron,

*E*, can be approximated by the expression of Opal

_{s}*et al.*

^{77}

where $E\xaf$ is a target-specific parameter proposed by Opal *et al.*,^{77} for example, 24.2 eV for neon. Employing this energy distribution to determine an average secondary electron energy, we may compute the energy transfer to the thermal population by RE impact ionization. Here, we see it is relatively insignificant, until lower *T _{e}* where we see this contribution approaching a similar order of magnitude as the RE-thermal electron Coulomb heating. This demonstrates that when any subsequent models seek to dynamically evolve an EEDF through disruption-like scenarios, the contribution of secondary electron formation should be considered alongside Coulomb heating.

An accurate quantification of the bulk electron temperature would need to account for additional energy transfer mechanisms such as the plasma particle and energy transport in space, and the energy exchange between bulk electrons and the ion/neutral population due to elastic collisions. This is the task of a coupled CR and plasma transport calculation, the fidelity of which is anticipated to be significantly improved by accounting for the relativistic CR corrections reported in this paper.

### C. Sensitivity to EEDF form

In Sec. IV A, a comparison was made between potential analytic EEDF forms and some experimental results of a reconstructed RE spectra. From the efforts of the community to understand RE distributions from experimental studies,^{10,72} as well as theoretical and computational studies,^{7,32,73,78–80} we know there are many different forms and shapes that a RE tail distribution may take. The shape of a RE EEDF will be heavily influenced by electric fields experienced, collisionality (with both electrons and impurity ions), and/or radiation such as synchrotron or Bremsstrahlung. Despite the many variations in potential RE distribution, one common concept in generalizing the description of a RE population is to describe a peak location or bump-on-tail of the electron population at relativistic energies. With this idea in mind, we present a brief discussion of trialing different analytic EEDFs that all have the same RE peak location, with the finding that the subsequent CR solutions yield very similar results.

Given $nAr=nD=1014$ cm^{−3} and the RE number density fraction set to $10\u22123$, we trialed Maxwellian (22), Gaussian (23), and boosted Maxwell–Jüttner (25) EEDFs to describe a RE electron population with a peak at 5 MeV. These EEDFs are shown in Fig. 13. Steady-state CR solutions were found over a range of thermal Maxwellian *T _{e}*, and the resulting average ion charge state curve is shown in Fig. 13.

As can be observed in the average ion charge state Fig. 13, we found that despite the variation of the trialed RE EEDFs around their peaks at 5 MeV, the outcome of the steady-state CR solutions was very similar for all three cases. While not shown here for brevity, the ion charge state distributions and subsequent RPL curve as a function of *T _{e}* were practically indistinguishable, as was the case for average ion charge state shown in Fig. 10. In our studies, multiple RE peak energies other than 5 MeV were also trialed, with these cases also reinforcing the observation presented in this work. Ultimately, from this trial, it appears that the key lever in the way RE populations may impact CR results is through the energy of the bump-on-tail peak at relativistic energies. This result appears to indicate that if one does not know the exact RE EEDF profile relevant to a particular problem at hand, in a theoretical or computational model of the problem, a suitable analytic form with the appropriate bump-on-tail location may serve as an adequate surrogate to obtain desired discharge properties in steady-state CR modeling.

### D. Impact on electron–ion collision rates in runaway electron modeling

So far in this study, our CR modeling approach has been solely based on inelastic collisions in plasmas. In addition to inelastic collisions, electron–ion *elastic* collisions can also have a critical role in runaway electron evolution through deflection, or pitch angle scattering, following collisions. For partially ionized impurities, partial screening of highly energetic projectile electrons can expose more of the nuclear charge than the fully screened situation, which is described by the ion charge. These can greatly enhance the pitch angle scattering rate of runaway electrons by weakly ionized high-atomic number impurities in a mitigated tokamak disruption.^{32,73,79} The relativistically corrected CR results, which address the inelastic collisions, can affect these electron–ion elastic collision rates through a modified charge state distribution. This effect will be quantified next.

Following Hesslow *et al.*,^{32} the collision frequency for elastic pitch angle scattering (deflection) of electrons is

where *p* is electron momentum normalized to *m _{ec}*, $\gamma =1+p2$ is the relativistic Lorentz factor, and $\tau c=(4\pi necr02\u2009ln\u2009\Lambda c)\u22121$ is the relativistic collision period. The relativistic Coulomb logarithm is given as

where $ln\u2009\Lambda 0$ is the thermal electron–electron Coulomb logarithm. Superthermal Coulomb logarithms for electron–electron and electron–ion collisions are

and

Further, $Zeff=\Sigma jZj2njne,$ with *Z _{j}* being the charge of an ion

*j*, and

*n*is the population density of said ion charge state. The atomic number of the impurity species considered is $Za,j$, that is, 18 for argon, 10 for neon. $Ne,j=Za,j\u2212Zj$ is the number of bound electrons of a given ion

_{j}*j*. Finally, an effective ion size $a\xafj$ is a fundamental input data quantity, taken from tabulated computation results.

^{81}

The inelastic collisions can also help increase the slowing down of impacting electrons. Hesslow *et al.*^{32} also devised a simple-to-use formula to account for the collisional energy loss due to inelastic collisions that are explicitly treated in our CR model. This Bethe-like slowing down frequency, which appears in the Fokker–Planck model as a collisional slowing down operator, has the form

where $I\xafj$ is the mean excitation energy of an ion, *j*, normalized to the electron rest energy.^{82} Again, this simplified model for inelastic collisions between energetic electrons and partially ionized impurities has an explicit dependence on the ion charge state distribution, which itself is strongly affected by the relativistic corrections in the CR model.

Rosenbluth and Putvinski^{78} proposed an often used rule of thumb for the energy loss rate, being that the effect of inelastic collisions can be represented by adding half of the bound electrons to the free electron density

As this was long thought of as a baseline approach for describing the effects of inelastic electron–ion collisions, we thus include a comparison with this expression in Fig. 14, where we present a comparison of scaled collision frequencies obtained by evaluating (30) and (31) using ion populations of a neon–deuterium discharge, given varying plasma temperatures. For both collision frequency types, we use the following scalings to simply aide comparison within a single figure, while still capturing the implicit electron density dependence of (30) and (31)

Comparison of the collision rates for pitch angle scattering and drag of electrons in Fig. 14 demonstrates clear differences due to different ion populations at the plasma temperatures shown. These results at intermediate and low temperatures, which may be representative of post-thermal quench plasma scenarios, indicate that modeling of RE formation and dynamics may be heavily influenced by the model chosen to provide charge state populations and thus collision rates as input data.

In these figures, agreement is found between the different CR model descriptions at high normalized momentum, *p* > 10, where ultra-relativistic electrons can be described. Through low to intermediate values of *p,* we see notable deviations.

At lower values of *p,* large differences are obtained, with the differences reducing with increasing *p*. While these large differences may not directly affect the scattering of a relativistic RE at say, 5–10 MeV, this will impact the modeling of how one describes the formation of RE populations, as initially low-energy electrons accelerate through phase space, and are ultimately captured in a runaway vortex phase-space flow.^{73} This window of sub-relativistic electron momentum/energy that connects the runaway population to the thermal bulk is crucial to accurately describing the conditions for runaway formation and subsequent avalanche, and therefore, provision of accurate input data to particle or kinetic models is a clear use case for the relativistically QED-corrected CR modeling outlined in this study.

### E. Effects of reduced order modeling

So far, we have discussed the role of including REs in a CR model, and how one may correct said CR model to more accurately model the enhanced inelastic electron scattering that occurs at relativistic energies. It was clear that differences arise compared to a simpler scenario in which only the thermal electrons are considered, or if a relativistic population is included but non-relativistic inelastic scattering cross sections are employed. Both of these approaches are standard assumptions that are often included when modeling tokamak plasmas with or without runaway electron component. As a bookend to this study, we briefly examine the effects of another standard assumption in MCF plasma modeling, use of a coronal, or coronal-like, model in lieu of a CR model.^{34,83} For this comparison, we compare the results of using the CR model discussed in this work to that from using a coronal-like model similar to what is popularly used in literature.^{24,33,34}

In the coronal regime, the main collisional and radiative processes in a plasma are generally assumed ahead of time given a low-electron density, *n _{e}*, in the plasma. In this limit, some assumptions

^{34,83}are often made: (i) the populations of ion excited states remain small compared with ground-state populations, (ii) collisional recombination is negligible compared to radiative recombination, and (iii) collisional de-excitation is negligible compared with radiative decay. As a result, a smaller and simpler model that does not explicitly resolve excited-state populations is produced. With a reduced-state vector of

*Z*+ 1, ground states of the possible ion stages of a given atomic target of atomic number

*Z*computation time are reduced at the expense of the fidelity of potential CR pathways that could occur in the plasma.

To simulate a coronal-like model, we solve our CR system with restricted CR pathways that may occur in a typical CR model used in our work. First, we ignore all three-body recombination rates and electron impact de-excitation rates. We assume that electron impact excitation and ionization transitions can only occur from the ground state of each ion stage, not excited states. Radiative recombination rates and radiative decay, or spontaneous emission, rates are left unchanged. Finally, autoionization, and thus electron capture, rates are left unchanged.

In this comparison, we show results for both neon and argon discharges, with equal parts deuterium density, $nI=nD=1014$ cm^{−3}, over a range of thermal temperatures. For the runaway component, we prescribe the same boosted Maxwell–Jüttner runaway component as used in Sec. V D with a maximum in the runaway distribution at approximately 6 MeV and a density equivalent to a 5 MA runaway current. We present some results of using a CR vs a coronal-like model in Fig. 15, where we show variation in the averaged charge state and standard deviation of ion charge states as a function of the temperature of the thermal plasma population alongside the RE population.

In general, the trends of coronal-like model results of the average charge state and standard deviation are in agreement with the results produced from the full CR model. We note a clear increased average charge state at low *T _{e}*, which can be found to be a result of two main things, the first being the limitation of three-body recombination in the coronal-like model. This process nominally acts to counter ionization in the presence of the RE component. As a result of this pathway being removed, the net ionization balance shifts higher. Second, the elimination of collisional de-excitation results in higher excited-state populations found in the coronal-like model than the CR model. As a result of this, autoionizing excited states are more populated, and thus, a net increase in the ionization balance is produced, due to the removal of competing collisional de-excitation processes.

In terms of the standard deviation, there are temperature windows where the CR result is greater than the coronal-like result and vice versa. As the interplay of the scattering and radiative pathways evolves over changing bulk *T _{e}*, the net effect of a broader or narrow distribution in the coronal-like model, compared to the CR model, is not as straightforward to trace as it is for the average charge state result. While the general trend of the standard deviation is comparable between the two models, nonetheless, the differences are apparent in the presented figures.

As such, the demonstration of such differences between a CR and coronal-like model when a RE population is present should highlight the fact that implementation of such reduced-order models can indeed have a considerable impact on the macroscopic plasma properties used to serve plasma modeling or interpretation of experimental measurements. Indeed, while the use of a coronal model can be very accurate for fusion plasmas where electron densities are sufficiently low, and the CR pathways in the plasma can be reasonably predicted,^{34} we however issue caution when considering non-trivial discharges, such as any containing a considerable relativistic electron component, because application of tried-and-tested reduced order models may not be as accurate as expected.

## VI. CONCLUSION

In this study, we have presented a systematic approach to analyzing the behavior of relativistic electron inelastic scattering and how this physics can impact the macroscopic properties of tokamak relevant discharges containing neon and argon impurities. We have extended upon recent findings^{14} to demonstrate appropriate benchmark results of electron scattering approximate formulas that can be implemented in CR models in a straightforward manner while still retaining the enhanced relativistic scattering. We subsequently have shown the effect that a minority RE population can have in discharge conditions applicable to modern and future tokamaks. Specifically, we observe consistent enhancement of average charge state on the order of 50% and broadening of charge state distributions. As a consequence of this effect, we observe clear changes to the associated radiation and spectral characteristics that can be attributed to REs, with variation of up to an order of magnitude from RPL values produced with a standard thermal plasma assumption. Along with the numerical change in the observed line radiation produced by the plasma, we have demonstrated that the spectral lines observed at low *T _{e}*, during regimes of RE-dominated scattering, are qualitatively different from those normally observed. These spectral properties may be detected via spectroscopy and provide a potential avenue for characterization and detection of when REs are present in cool, post-disruption discharges.

Given the potential for variation in the ion charge state distributions when REs are properly accounted for, compared to simple assumptions of thermal-only EEDFs, we demonstrated how electron scattering models used in modern kinetic and particle-based plasma simulation^{32} can produce significantly different pitch-angle and energy loss collision rates, often varying by factors of up to 4 or 5 in certain conditions.

Finally, we have presented some differences produced from a CR modeling approach, and a reduced-order coronal-like model often employed in modeling tokamak relevant discharges. We note that such reduced-order models are indeed sufficiently accurate in well-understood, thermal plasma scenarios but advise caution and careful application of such models when dealing with complex plasmas, such as those discussed in this work featuring non-negligible populations of relativistic electrons. Ultimately, we hope to have demonstrated the benefits of accounting for often neglected relativistic scattering considerations when facing the task of modeling fusion plasma scenarios with minority populations of relativistic electrons. Further, we endeavor to generally promote expanded collaboration and research opportunities between plasma and atomic physics communities to achieve greater outcomes in both fields, so that we may together gain momentum in solving crucial problems, such as the runaway mitigation challenge outlined in this study.

## SUPPLEMENTARY MATERIAL

See the supplementary material for this article is a brief animation of neon charge state distributions for the conditions presented in Sec. V A. In this animation, we cycle through the ion charge state distributions found with decreasing values of *T _{e}* given the runaway electron tail described in Sec. V A. Key features of the animation include the following: (i) increased average ion charge at lower

*T*due to runaway driven processes and (ii) an accompanying broad trailing tail of ion states populated by the runaway electron collisions.

_{e}## ACKNOWLEDGMENTS

This work was jointly supported by the U.S. Department of Energy through the Magnetic Fusion Theory Program and the Tokamak Disruption Simulation (TDS) SciDAC project at Los Alamos National Laboratory (Contract No. 89233218CNA000001). Additional support is provided by the Laboratory Directed Research and Development Program of Los Alamos National Laboratory under Project No. 20200356ER. This work was supported by the R&D Program of the Korea Institute of Fusion Energy (KFE), which is funded by the Ministry of Science and ICT of the Republic of Korea (No. KFE-EN2141–7).

## AUTHOR DECLARATIONS

### Conflict of Interest

No conflicts of interest to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## References

*NIST Atomic Spectra Database (Ver. 5.6.1)*(