Wall heating by subcritical energetic electrons generated by the runaway electron avalanche source

Subcritical energetic electrons (SEEs) produced by the runaway electron (RE) avalanche source at energies below the runaway threshold are found to be the primary contributor to surface heating of plasma-facing components (PFCs) during final loss events. This finding is supported by theoretical analysis, computational modeling with the Kinetic Orbit Runaway electrons Code (KORC), and qualitative agreement with DIII-D experimental observations. The avalanche source generates significantly more secondary electrons below the runaway threshold, which thermalize rapidly when well-confined. However, during a final loss event, the RE beam impacts the first wall, and SEEs are deconfined before they can thermalize. Additionally, because the energy deposition length decreases faster than energy, the deposited energy density, and thus the maximum PFC surface temperature change, is larger for SEEs than REs. KORC simulations employ an analytic first wall to model particle deconfinement onto a non-axisymmetric wall composed of individual tiles. PFC surface heating is calculated using a 1D model extended to include an energy-dependent deposition length scale. Simulations of DIII-D qualitatively agree with infrared (IR) imaging only when SEEs from the avalanche source are included. These results demonstrate that SEEs are the dominant contributor to PFC surface heating and indicate that the avalanche source plays a critical role in the PFC damage caused during final loss events. The prominence of SEEs also has important implications for interpreting IR imaging, one of the primary diagnostics for RE-wall interaction diagnosis, despite REs dominating the energy and current density. This result improves predictions of wall damage due to post-disruption REs to estimate material lifetime and design RE mitigation systems for ITER and future reactors.

States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, worldwide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes.The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (www.energy.gov/downloads/doe-public-access-plan).
* * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence.Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
An initial runaway electron (RE) can undergo a large-angle collision with a thermal electron where enough energy is transferred to the latter to move it over the runaway threshold [1,2], where electric field acceleration and collisional friction are balanced [3,4].Large-angle collisions of REs comprise the avalanche source of secondary electrons and lead to population exponentiation with a rate proportional to the plasma current [5].This exponential growth makes the RE avalanche mechanism perhaps the most dangerous consideration for ITER and reactor-level tokamaks [6,7].Additionally, compared to an initial RE generating a secondary RE during a large-angle collision, there is a higher probability for the initial RE to generate a subcritical energetic electron (SEE) that is born below the runaway threshold but significantly more energetic than the thermal background electrons.The avalanche source is the most significant RE generation mechanism during tokamak disruptions, where large Ohmic electric fields result from the rapid increase of resistivity during the thermal quench and large inductive electric fields originate from the change in current during the current quench, e.g.[6].Furthermore, the baseline disruption mitigation system designed for ITER plans to use a massive injection of impurities, which provides an additional inventory of free and bound electrons to undergo large-angle collisions [8].While previous studies have expanded our understanding of the avalanche threshold and amplification of seed RE populations during the current quench phase of a disruption, the avalanche mechanism also plays a role during RE mitigation.
References [1,2] were the first to consider the avalanche effect of REs, and [5] was the first to numerically verify the model, which assumed a primary RE with infinite energy and zero pitch angle.Restrictions on the primary REs were relaxed in [9] that allowed for finite energy primaries, and [10] and [7] that removed the assumption of zero pitch angle.Embréus et al [11] derived the operator from the relativistic Boltzmann operator, thereby rectifying the issue of double counting by including distribution function sinks in addition to the source, which [10] shows is important for continuous kinetic calculations, as opposed to previous studies using Monte Carlo kinetic calculations.These operators have been primarily used in the modeling studies of RE generation, e.g.[12] and references therein.
The avalanche source will also play a significant role during the final loss event of RE mitigation, where the RE beam impacts the wall and any remaining current rapidly quenches [13,14].During the final loss event, magnetic surfaces open and terminate on plasma-facing components (PFCs), which results in particles deconfining on rapid time scales.With confining magnetic fields, generated SEEs quickly lose energy due to collisional friction.Conversely, they can be rapidly lost when magnetic fields are open and significantly contribute to PFC heating.Additionally, the avalanche mechanism will evolve the energy distribution of the RE population.The volumetric heating is affected because the deposition length scale into plasma-facing components (PFCs) is energy-dependent.
Lastly, the amplification of the RE beam via the avalanche source can enhance the magnetic to kinetic energy conversion [15] due to the electric field induced by the current quench during the final loss event.Ultimately, the deconfinement of SEEs due to open magnetic fields, the evolution of the RE distribution, and increased magnetic to kinetic energy conversion due to the avalanche source have consequences for PFC heating to consider.
The present manuscript aims to elucidate the contribution of secondary electrons generated via large-angle collisions comprising the avalanche source to wall heating during RE final loss events, particularly the role of SEEs.We organize the remainder of this paper as follows.In section 2, we will review the avalanche source and the distribution of secondary electrons it produces and investigate their role in PFC surface heating by using an energy-dependent heat transport model.In section 3, we introduce the experimental case of DIII-D #177031 and modeling elements used in the Kinetic Orbit Runaway electrons Code (KORC) code [16] to support our findings in section 2. In section 4, we present the results of simulations with the KORC, detailing the evolution of the secondary electron distribution during the final loss event and connecting this to their role in wall heating.Lastly, in section 5, we offer concluding remarks and discuss the results in the context of RE mitigation efforts.

Theoretical background
In this section, we begin by reviewing the avalanche source, including partially ionized impurities.We then include RE energy dependence into a previously used PFC surface heating model, enabling us to explore the significance of SEEs generated by the avalanche source on the PFC surface heating.

Large-angle collision operator
The avalanche source of secondary electrons is due to largeangle collisions of existing REs with the background thermal electron population.References [7,10] provide the continuous operator that evolves the electron distribution due to the addition of secondary electrons as where subscripts 0 and 1 refer to quantities corresponding to the initial and secondary distributions from a collision, n tot is the total electron density from free and partiallyionized impurity bound electrons, p is the magnitude of the momentum, η is the pitch angle, and φ is the gyro-angle.dσ/dγ is the Møller differential cross-section for free-free electron-electron scattering [17] given by with γ = √ 1 + (p/m e c) 2 the Lorentz factor, r e = e 2 /(4πϵ 0 m e c 2 ) the classical electron radius, n α with α = 1, 0 are the unit vectors along the initial and secondary momenta making n 1 • n 0 the deflection angle, and cos η γ encapsulates energy and momentum conservation laws in a collision, yielding Note that cos η γ is > 1 for γ 1 > γ 0 and is not an accessible region.Because the gyro-angle appears in the momentum unit vectors, we integrate the Dirac-delta function over the gyroangle, which yields While equation (1) only accounts for the change in the RE distribution function due to the addition of secondary REs, the RE distribution will also change because energy and momentum must be conserved in the large-angle collision creating a secondary RE [11].Including these conservation properties is straightforward in the Monte Carlo formalism and is discussed in appendix A. Lastly, to exclude double counting in the Fokker-Planck collision operator, we include alterations to the electron-electron Coulomb logarithm [11,18] by employing the form from [18] lnΛ LA ee = lnΛ ee + ln where γ min = √ 1 + p min /(m e c) is the minimum energy a secondary electron can have in a given simulation.
Because Monte Carlo modeling involves discrete particles, we can calculate the source of secondary REs due to each initial RE by using the single-particle distributionJ function for each initial RE i in equation ( 1) to yield Note that from the 3D Dirac-delta function of the singleparticle distribution function in equation ( 6), a secondary RE has the same location as the primary RE.Ar /2 = n +0 Ar /2 = n e /3, T e = 1.5 eV, and Z eff = 1.The solid black line borders the allowable secondary region.It is given by where the argument of the radical in the denominator in equation ( 7) equals zero and is equivalent to the region where energy and momentum are conserved.This can be compared to similar operators plotted in [10] and [11] in (K, sin 2 η) and (p ∥ , p ⊥ ) phase-space.The solid red line indicates the approximate RE threshold from [19] using the effective critical momentum p crit in the presence of partially-ionized impurities from [8], which yields where ξ = cos η, E = 9.45 V m −1 is a representative value of the inductive electric field from DIII-D #177031, and E CH = n e e 2 lnΛ/(4πϵ 2 0 m e c 2 ) is the Connor-Hastie field.The effective critical momentum in equation ( 8) is given by where νS and νD are the slowing down and pitch angle scattering collision frequencies shown in [23] normalized to their values in the absence of partially-ionized impurities.The black, dot-dashed, vertical line in figure 1(a) shows the critical energy corresponding to the momentum in equation (9).Note that p crit /(m e c) ∼ 1 for the present example, making it approximately the same as that used in [24] where E CH,tot = (n tot /n e )E CH = 3.78 V m −1 is the Connor-Hastie field including all electrons, free and bound to partiallyionized impurities.Because this form is more straightforward to evaluate than equation ( 9), we will use this form in the calculations presented in section 4. Lastly, we set the maximum p 1 from energy conservation, namely that the energy of the secondary RE can at most be equal to the energy of the initial RE after the collision, which yields The blue trace multiplies the PFC surface temperature change by the differential of the density from the secondary RE source to quantify the energy dependence of the heating effect of the secondary RE source.The green trace (right axis) shows an analytic expression of the CSDA range [20] found by fitting the high energy ESTAR [21] and low energy IMFP [22] databases.
The black, dashed, vertical line in figure 1(a) shows the maximum energy of a secondary electron.

Energy-dependent PFC heating model
To understand how the avalanche distribution will change the surface temperature of the PFCs, we employ a model derived from the one-dimensional solution of the heat diffusion equation in a semi-infinite solid [25].Lehnen et al [26] first used this model, assuming an exponential decay of the RE energy deposition into the PFC with a constant heat flux density inferred from experimental observations.An evolving heat flux density from a zero-dimensional model was incorporated by [13] by using Duhamel's theorem [25].We can extend these previous efforts by calculating the energy-dependent deposition length scale for single REs at different energies.This results in an equation for the change in PFC surface temperature where we use parameters corresponding to ATJ graphite at a constant temperature of 800 K: κ = K/(ρc p ), K = 55.1W(m • K) −1 is the heat conductivity, c p = 1640 J (kg • K) −1 the specific heat capacity, and ρ = 1730 kg m −3 the density.q is the kinetic energy deposited across the PFC interface per unit time and unit area by a single RE, and has the same units as heat flux density 4 .The energy-dependent deposition length scale δ dep corresponds to an analytic expression of the continuousslowing-down approximation (CSDA) range from [20] that fits the high energy ESTAR [21] and low energy inelastic mean free path (IMFP) [22] databases to semi-empirical models.We plot δ dep as a function of particle kinetic energy in the green trace in figure 1(b) with values on the right axis.The CSDA range assumes the rate of energy loss along the path is equal to the total stopping power from [27]; this is generally not consistent with the exponential decay of the RE energy deposition into the PFC.Because KORC evolves individual particles' locations and momenta, equation (12) will be directly applied to simulation results in section 4.3.2.The red trace in figure 1(b) shows the maximum PFC surface temperature change due to a single particle as a function of kinetic energy with values on the left axis.Along with the RE kinetic energy, q i also requires a deposition area and time.Because the time interval and surface area of a single RE impact are infinitesimal, we normalize the deposition time and area out of equation (12) when plotting in figure 1, which can be observed in the units on the left axis.There is a nonmonotonic change in PFC surface temperature coinciding with the turning point of the IMFP.Physically, for kinetic energy in the range K ∈ [ 10 2 , 10 6 ] eV the energy deposition length scale increases faster than the kinetic energy, leading to a decreased maximum PFC temperature as energy increases in this range.Note that there is also an increase as the kinetic energy increases past the plotting range.However, this energy range is not expected to be significant for DIII-D [28,29].
We then multiply the maximum PFC temperature change by the normalized integral of the avalanche source over the secondary pitch, found by integrating over η 1 and normalizing by the value at 100 eV.The blue plotted line shows the effective PFC surface heating due to the avalanche source for an initial RE with kinetic energy K 0 = 10 MeV and pitch angle η 0 = 10 • , where all the secondary REs are rapidly deconfined.The change in PFC surface temperature decreases at higher energies, reaching a maximum value substantially below the approximate runaway threshold in figure 1(a).
The equation for the change in PFC surface temperature equation (12) captures the evolution in time, and figure 1 shows the maximum of ∆T surf .Figure 2 shows the dynamics of PFC surface temperature spanning the energy range plotted in figure 1.We have used a time range that is uniformly spaced in the decimal logarithm of time to resolve the smalltime scale dynamics.The scaled complementary error function in equation (12) captures the heat diffusion of the exponential energy deposition, and it becomes more transient as the initial energy decreases.Thus, while the lower energy SEEs provide increased maximum PFC surface temperature changes, this heating is temporary.
The PFC surface heating model equation (12) used in this work advances over previous studies [13,26].However, future studies need to address several shortfalls.First, the present model uses an exponential volumetric energy deposition to be analytically tractable, but this precludes the more physical Bragg curve [30].Next, there are several shortfalls in modeling the heating of graphite tiles.While the one-dimensional solution applies to the semi-infinite plane, this assumption is least valid at the leading edge, where most of the heating occurs.The model of volumetric energy deposition does not account for the generation of secondary electrons within the graphite.As tiles heat up, the thermal conductivity will decrease while the specific heat capacity increases, but this is not considered in the present modeling.While infrared (IR) imaging observes graphite tile sublimation and an explosion event from one of the leading edges (not shown), we limit the focus of the present work to surface heating.While graphite sublimates, metal PFCs will undergo a phase transition and melt, requiring a more sophisticated model to consider the flow of molten metal.Future work will aim to bridge these gaps by using higher fidelity models for volumetric energy deposition and surface melting dynamics, building on the work of [31].

Physics model
This section provides an overview of the model elements contained in KORC that this study uses.The full details of implementing the relativistic guiding-center equations of motion in cylindrical (R, ϕ, Z) coordinates with synchrotron and bremsstrahlung radiation, incorporating spatiotemporal experimental reconstructions of magnetic and toroidal electric fields, and Fokker-Planck linearized Coulomb collision operator appear in [23].We have also implemented the large-angle collision operator comprising the avalanche source of secondary electrons from section 2.1.We describe the implementation of this operator in appendix A and provide benchmarking with previous studies in appendix B.
To explore the role of PFC surface heating by SEEs detailed in section 2, we have chosen to examine the final loss event in DIII-D discharge #177031.In this discharge, an injected cryogenic Ar pellet triggers the disruption, and the resulting RE beam is formed and controlled.Secondary deuterium shattered pellet injection (SPI) does not cool the plasma to the point of fully recombining and is followed by a ramp down of the central solenoidal field, resulting in a rapid loss event to the high field side (HFS) inner wall limiter.This discharge is the 'ionized' case in [32].Figure 3(a) shows the current evolution (blue trace), the maximum toroidal electric field from experimental reconstructions as calculated in [23] (green trace), and the average electron density measured by an interferometer in a horizontal chord across the midplane (red trace) during the final loss event.We observe a current spike near t = 1.6015 s that suggests a magnetohydrodynamic (MHD) event, but we ignore its effects on the present work.During the RE final loss event, the average electron density increases transiently.This is likely due to the wall heating by the REs, causing the graphite to sublimate and carbon to ionize.Note that [32] also saw this for the 'recombined' case, where an unstable MHD mode causes a final loss event with a large wetted area.Following the final RE loss event, a small remnant current slowly decays by 1.63 s.Furthermore, there is a transient increase in the average electron density near the end of the final loss event when the RE beam scrapes off on the wall.However, to simplify the present modeling, we use a constant and uniform electron density with magnitude n e = 2.25 × 10 20 m −3 .
The three plots to the right in figures 3(b)-(d) indicate the induced toroidal electric field in the colormap overlaid with contours of poloidal magnetic flux for times before (b), during (c), and after (d) the final loss event.The three times show that the magnetic surfaces open up during the final loss event, leaving no closed surfaces afterward.The RE beam, which self-consistently provides the poloidal magnetic field, advects toward and scrapes off on the inner wall, inducing a large toroidal electric field.The induced toroidal electric field is larger than the numerically calculated critical electric field consistent with equation ( 9) We note that the ratio of E crit /E CH,tot ≈ 1.15 is consistent with that found by [33].From any REs not deconfined when the RE beam scrapes off and opens magnetic surfaces to the wall during the final loss event, the avalanche source will generate SEEs and additional REs.
Lehnen et al [26] suggests that REs could deposit their energy inhomogeneously on the upper divertor of JET due to tile misalignment, which would be significant at shallow angles of incidence to the wall.For the present case in DIII-D, REs and SEEs are being deconfined on the inner wall, comprised of 15 cm wide graphite tiles [34].To investigate the deposition of REs and SEEs onto a more realistic inner wall compared to the axisymmetric case in [23], a faceted wall has been implemented in KORC to model the inner wall limiter on DIII-D.We approximate the wall as a regular polygon, where the major radius of the wall as a function or toroidal angle is given by where R circ = 1.016 m is the circumradius at the leading edges where individual PFC tiles meet, and n tiles = 42 is the approximate number of 15 cm wide graphite tiles.Note that tile gaps are not included in this model because they would add substantial complexity to this initial study.However, adding tile gaps would introduce more leading edges, and the heating at tile corners could be different than the 1D model used in the present work.We show a schematic of the analytic inner wall used in KORC calculations in figure 4. Using the chosen circumradius and number of tiles results in a leading-edge located at a radius ∼3 mm greater than that at the center of the tile.The inset of figure 4 shows the offset in wall location at a smaller major radius compared to the circumradius except at leading edges.In section 4.3.2,we will see that nearly all deconfined electrons deposit along the leading edges.

Results
This section will present the results of KORC simulations of initial RE and secondary electron distribution evolution in DIII-D discharge #177031.The presented results include reduced-model calculations used to choose neutral content and set the lower bound on energy for computed electron orbits and the system evolution using the full physics model.They will culminate in the analysis of PFC heating due to all deconfined electrons.The simulations in this study, unless otherwise noted, employ an initial RE distribution with N RE = 1.088 × 10 4 particles that are monoenergetic with K = 10 MeV, monopitch with η = 10 • , and spatially uniform within the last closed flux surface from JFIT reconstruction at t = 1.5947 s as shown in figure 3. We choose this initial distribution following the modeling from [23] that showed only quantitative differences between using such a distribution compared to the experimentally-inferred distribution found in [28].
To characterize the evolution of the ensemble of particles at a macroscopic level for comparing with DIII-D current measurements, we define the simulated current I as derived in [23] where b ϕ,i is the unit vector of the magnetic field in the toroidal direction and H RE effectively splits the ensemble of particles into two populations, confined REs and thermalized or deconfined electrons.Note that this definition is consistent with the discussion of the lower energy bound in the following section.

Estimating neutral content and lower energy bound
Using experimental data from the horizontal midplane interferometer shown in figure 3, we set the electron density to be constant and uniform with a magnitude of 2.25 × 10 20 m −3 .However, because there was a cryogenic Ar pellet injected to cause the disruption and additional D 2 SPI attempting to mitigate the RE beam, the ratio of Ar +1 and D +1 in the post-disruption companion plasma is unknown.Additionally, because the companion plasma is ∼1 eV [35,36], there is a relatively low ionization factor and a significant level of Ar +0 neutrals.Beidler et al [37] found that for DIII-D, setting a neutral impurity density equal to the ionized impurity density led to an agreement between KORC modeled and observed current reduction during post-disruption secondary high-Z impurity injection with Ne SPI.We adopt the same assumption here.We performed KORC calculations without the avalanche source turned on (not shown), where different ratios of Ar +1 (Ar +0 ) and D +1 are used.The different ratios affect the Fokker-Planck collision operators, and the KORC simulation using a ratio of n +1 Ar /n D +1 = 2 provides the best match to the experimentally observed evolution of current.While this choice ignored higher charge states for simplicity, efforts are ongoing to use higher fidelity, experimentally-inferred plasma and impurity profiles from the model developed in [38].
During the final loss event, a SEE can only deposit its energy to PFCs if deconfined before collisional friction thermalizes it.This balance effectively sets the lower energy bound that KORC calculations must resolve to model SEE dynamics.To determine where the deconfinement and thermalization times are equal, KORC calculations that do not evolve spatial degrees of freedom are used to evaluate thermalization time scales.In contrast, we use collisionless calculations to assess deconfinement time scales.Each simulation starts with 2176 electrons and runs until the distributions achieve a constant number of electrons after particles are confined or thermalized.Spatially independent simulations are run with the plasma parameters discussed in the previous paragraph and the maximum electric field inferred from the experiment.Collisionless simulations use the experimental magnetic and electric field configuration at t = 1.6022 s during the RE final loss with open flux surfaces and maximum electric field.Particles in the collisionless simulation are initialized uniformly within the flux surface that encompasses the HFS limiter.The time scales are taken to be when the number of electrons dropped below halfway between the initial and final total.Figure 5 shows two sets of calculations varying the initial kinetic energy of a monoenergy/monopitch distribution with different pitch angles.While the thermalization time is nearly independent of initial pitch angles, the collisionless calculations indicate that the change in parallel velocity along the open magnetic field lines strongly affects deconfinement times.Figure 5 indicates that for electrons with η = 10 • , the deconfinement time scale is shorter than the thermalization time scale for ≲3keV, and only these electrons will contribute to the wall heating.While the collisionless simulations include many trapped electrons at higher energies for η = 80 • that do not deconfine, the effects of electric field acceleration and collisional pitch angle scattering cause detrapping on a faster time scale than collisional friction energy loss time scale.Both time scales are also well below the time resolution of the experimental reconstructions at 0.5 ms.A more thorough understanding of this lower energy bound on the SEE energy range is of interest, notably how it scales toward ITER and future reactors.A follow-up study will explore this topic further.
To avoid having to self-consistently calculate p crit according to equation ( 9), KORC calculates the approximate form according to equation (10) using the global maximum toroidal electric field and evaluating the Coulomb logarithm of equation ( 5) at 100 keV, which is consistent with a more precise calculation shown in figure 1(a).Then, we choose γ min = [γ crit (p crit ) + 1]/2, making p min (γ min ) the maximum possible momentum of a secondary electron generated by an initial electron at p crit .Lastly, we include the option to scale the minimum momentum lower as p ′ min = a pscale p min .For the following calculations, we use a pscale = 0.25 that corresponds to a lower energy bound of K min = 4.94 keV for the maximum toroidal electric field at t = 1.6017 s.Note that we have also used a lower scaling of a pscale = 0.125 corresponding to K min = 1.24 keV, and do not get qualitatively different results.The lack of variation using this minimum energy is consistent with the results shown in figure 5 that the lower bound on energy is where the deconfinement and thermalization time scales are equal.

RE current and distribution function evolution
After performing reduced-model KORC calculations without the RE avalanche source, and separately collisionless or spatially-independent, to estimate the impurity content and lower bound on the energy needed to capture SEE deposition onto PFCs, we continue with KORC simulations that include the avalanche source.This section will focus on the evolution of experimentally-relevant macroscopic quantities, the initial and secondary electron distribution functions, and the spatial distribution of the generated secondary electrons.Because it would be computationally expensive to simulate all the secondary electrons that are generated and thermalized in regions of closed magnetic field lines, we only include the avalanche source when the initial RE beam begins to advect into the inner wall.We choose this time to be t = 1.5997 s, which coincides with when initial REs first deconfine (not shown).Additionally, we only use the avalanche source for particles where equation ( 10) is finite, effectively decreasing the number of secondaries to ensure computational tractability.This assumption likely leads to an underestimate of the surface heating.When the RE beam, which self-consistently provides its own confining poloidal magnetic field, advects into the HFS inner wall, its outer REs and flux surfaces are scraped off on the wall.The advection of the beam toward the HFS and the scraping-off causes a change in the poloidal magnetic field, and an inductive electric field is generated according to Faraday's law.The scraping-off generates an inductive electric field with a poloidally symmetric structure that accelerates REs.At the same time, the advection will generate a RE-accelerating electric field located on the trailing side of the beam toward the LFS.In the present case, because the largest induced fields are always on the trailing side of the RE beam when the magnetic field region connects to the wall, this methodology includes the most critical population of avalanche-generated secondary electrons.In future studies, we plan to relax these restrictions using a GPU-accelerated version of KORC that can efficiently simulate a significantly increased number of electrons.
Figure 6 shows high-level results of a full KORC simulation.In figure 6(a), we further subdivide the total current I (black trace) into the components from the initial REs with subscript 'init' (red trace) and secondary electrons with subscript 'sec' (blue trace).The violet trace shows the experimentally measured current from DIII-D #177031.The evolution of the simulated current is qualitatively similar to the experimental measurement.However, it is plausible that the MHD event around 1.615 s affects the confinement and pitch angle distribution to cause a difference in the current evolution.Figure 6(b) focuses on the secondary electrons generated during the final loss event.The black trace shows the total number of secondaries generated, and the solid blue and red traces correspond to the total thermalized and deconfined secondaries.At the beginning of the final loss event, we observe that most of the secondaries are thermalized before they are deconfined, and it is not until the end of the final loss event when the flux surfaces are fully opened to the wall, that secondaries are deconfined before they are thermalized.We also plot subsets of the rapidly thermalized or deconfined populations in the dotted traces.We define that an electron is rapidly thermalized or deconfined if these events happen in the same time interval between experimental reconstructions (0.5 ms) as the avalanche source that generates them.We note that this definition of 'rapid' is a misnomer given the time scales in figure 5.However, this is the available time resolution of experimental reconstructions.Using this definition, we observe that nearly all the secondary electrons that are thermalized or deconfined do it rapidly.
We next turn from the evolution of macroscopic quantities to the evolution of the electron distribution functions.Figure 7 shows the evolution of the initial RE (a), (b) and secondary electron (c), (d) distributions of kinetic energy (a), (c) and pitch angle (b), (d).Note that these distributions include only non-thermalized electrons while including confined and deconfined particles.The initial REs begin as a monoenergetic/monopitch distribution and are evolved primarily by collisional effects, electric field acceleration, and mirror forces along their orbits.Because of the mirror force, many deconfined particles are at pitch angles close to 90 • when the RE beam scrapes off on the HFS inner wall limiter, which appears as a local maximum in the observed pitch angle distribution of the initial REs.For the secondary electrons, we only plot distributions starting at t = 1.6007 s due to the lack of good statistics during the time interval just after including the avalanche source.While it can be shown that the avalanche source from equation ( 7) has an energy distribution ∼ K −1 to leading order, collisional friction and electric field acceleration effects will lead to energy distribution ∼ exp[−E/T] in asymptotic time state [5].However, particles deconfine before reaching a time-asymptotic state, resulting in an energy distribution more closely resembling a power law than an exponential one.Because the avalanche distribution includes a larger proportion of low-energy electrons, collisions play a more significant role for secondary electrons than for the initial REs, resulting in a flatter distribution in pitch for the secondary electrons.Lastly, we observe a hollow pitch angle distribution of the secondary electrons at 90 • .This is due to the plotted distribution only showing the non-thermalized secondary electrons.When we look at the thermalized population (not plotted), the pitch angle distribution is ∝ sin η.The nonthermalized electrons are deconfined to the first wall before they can equilibrate their pitch angle.Additionally, the secondaries near 90 degrees would not be accelerated by the electric field to the same extent due to trapping and would thermalize more readily.Note that trapping is not considered in the spatially-independent collisional calculations shown in figure 5.
The final result we present in this section is the evolution of the spatial generation of secondary electrons via the avalanche source shown in figure 8.As described at the beginning of this section, the present study only applies the avalanche source for a given electron if equation ( 10) is finite, which limits the localization to the trailing edge of the RE beam as it advects toward the HFS and scrapes off on the inner wall limiter.Figure 8 shows snapshots in time of the spatial distribution of secondary electrons generated by the avalanche source.The generation is limited early by the relatively low electric field and later due to the confinement loss of the initial REs.We observe the generation of secondary electrons in regions of open magnetic fields at the end of the RE final loss event.For a circular cross-section tokamak assuming constant safety factor, [16] states that electrons will contain an 'energy drift' according to where R 0 is the major radius of the magnetic axis, v ϕ is the component of the velocity in the toroidal direction, and Ω e = eB/γm e is the relativistic electron gyrofrequency.As the energy increases, v ϕ will approach c and γ will increase without bound, increasing the poloidal radius of the orbit R c and shifting from the magnetic axis a distance ∆.Because the initial and secondary REs have significant energy drifts, they can remain confined even when the magnetic axis has shifted within the first wall to the HFS and continue to generate additional secondary electrons.The lower energy, secondary SEEs have an insignificant energy drift and will subsequently travel along magnetic fields after being generated, deconfine rapidly, and deposit their energy into the first wall.

Spatially-dependent kinetic PFC wall heating
In the previous section, we detailed how REs in the initial RE beam generate secondary electrons via the avalanche source, which are then rapidly deconfined to the first wall in open magnetic regions during the RE final loss event.Section 2.2 discussed the change in PFC surface temperature for single REs as a function of energy and used the avalanche distribution to estimate the heating of the RE avalanche source.We begin this section by recasting equation ( 12) for the ensemble of REs and SEEs and use it to evaluate the PFC surface heating during a RE final loss event.Note that in this work, the effects of the plasma sheath on the RE distribution deposited on the PFC surface are neglected for simplicity.
Because the heat diffusion equation with constant thermal conductivity and specific heat is linear, we can extend equation (12) to calculate the energy-dependent deposition length scale for a spatially-dependent ensemble of electrons by adding each contribution to calculate the total change in surface temperature.We express this in integral form by where d 3 X is the 3-dimensional volume element and r w is the location of the wall, is the distribution function of particles impacting the wall with X ∈ R 3 denoting the spatial location of an electron and t w the time an electron strikes the wall, and with ∆t and the deposition time and area will be determined later.Note that in the present study, due to the computational unfeasibility of tracking low-energy electrons, that any simulated electron whose momenta drops below p min is considered thermalized and does not contribute to the wall heating calculations.
In the following subsections, we determine the proper angle of incidence of an electron modeled using a guiding center orbit for equation (17) and then look at the relative importance of the initial REs and the secondary electrons generated by the avalanche source in PFC surface heating and make comparisons with experimental observations.

Approximating angle of incidence to PFCs.
In equation ( 17), the energy-dependent deposition length scale is given by where δ dep,i (K i ) is the energy-dependent deposition length scale for a given kinetic energy K i as determined from the composite analytic expression fitting ESTAR and IMFP databases in [20], and sin Θ is the sine of the angle of incidence on the inner wall limiter.The KORC simulations in the previous section followed electron guiding center orbits, evolving two velocity dimensions rather than the full three by ignoring the evolution of the gyroangle.However, the gyroangle is important for calculating the angle of incidence, as seen in the schematic of figure 9.The full orbit given by the blue helix travels approximately parallel to the magnetic field as it intersects the PFC surface, but depending on the gyroangle, the angle of incidence between the PFC and the purple vector v, indicating the electron motion at the time of impact, can be significantly different than the guiding center angle of incidence θ GC .In the primed coordinate system rotated along the x-axis according to, we describe the full orbit as where v ⊥ = v sin η is the magnitude of the velocity perpendicular to the magnetic field, v ∥ = v cos η is the magnitude of the velocity parallel to the magnetic field, Ω is the gyrofrequency, and α the initial gyroangle.We define the angle of incidence of a particle hitting the wall as where χ = Ωt + α is the gyroangle, and sin Θ > 0 requires the particle to intersect the surface from the plasma region.This restriction constrains the gyroangle to a range of values depending on the pitch and guiding center angle of incidence.From our KORC simulations, we directly evolve the pitch angle η.KORC also retains trajectory locations when crossing the analytic wall surface, allowing us to do a linear interpolation afterward to calculate the deposition location and θ GC .Because KORC does not calculate the gyroangle for guiding center orbits, we use the definition sin Θ > 0 to constrain the possible gyroangles.In the calculation of PFC surface heating in the following section, for the allowable range of gyroangles, we effectively split every simulated electron into 10 electrons, each having 1/10 of the energy with a randomly chosen gyrophase.
Figure 10 shows distribution functions of the various angles calculated for the simulated ensemble of electrons that strike the wall for all (blue trace), initial RE (red trace), and secondary electron (green trace) populations.Note that in figure 10(a) of the distribution of sin Θ, we scale all the distributions to the distribution of the ensemble of all electrons.Because the limiting surface is on the HFS wall, particle trajectories that penetrate the wall are traveling mostly in the toroidal direction with R∆ϕ ∼ 20∆Z, and θ GC is similar for both initial RE and secondary electron populations as seen in figure 10(b).Thus, the qualitative difference between sin Θ for the initial RE and secondary electron populations is caused by the different pitch angle distributions as seen in figure 10(c).Note that the pitch angle distributions shown here are the timeasymptotic distributions of those shown in figure 7. Figure 10(d) shows the allowable range of gyroangles for particle deposition.From the convention used in equations (22a), (22c) for the gyroangle, χ = 0 is when an electron has a minimum value of z ′ and v z ′ > 0 for χ ∈ [0, π].It makes sense that sin Θ < 0 in this interval of χ because the electron is moving away from the z = 0 surface.

Comparing wall heating of primary and secondary
With a reasonable method for calculating the angle of incidence to include the gyrophase, we can continue by calculating equation ( 17) for the change in PFC surface temperature.First, because we are simulating many orders of magnitude (∼10 13 ) fewer particles than exist in the experiment, we compute a proportionality factor by taking the ratio of the simulated to experimental current at the initial time I DIII-D /I KORC = 4.5 × 10 12 .Next, we determine when a particle is lost to the wall and bin each particle within the 0.5 ms time intervals separating the experimental magnetic reconstructions.We also bin the surface of the inner wall to smooth out the individual particle deposition locations and to provide a measure of the unit area for calculating the heat flux density.
Because the simulated inner wall is a regular polygon and we use axisymmetric magnetic reconstructions, we would expect the heating pattern to be qualitatively similar across each tile.However, this is not the case (not shown) due to the finite number of particles used in our simulation.To look more closely at the evolution of the local heating at leading edges, we proceed by averaging over each tile in the toroidal direction.The spatial resolution of the tile average uses 50 bins across a simulated graphite PFC tile in ϕ and 30 bins in Z to yield an area of deposition (R∆ϕ × ∆Z) ≈ (0.003 m × 0.053 m) = 1.6 × 10 −4 m 2 .Note that there will also be several tile interfaces in the Z direction, but because the DIII-D inner wall has no curvature in the vertical direction other than near the lower hybrid current drive antenna, we do not treat these interfaces as leading edges.We then calculate each particle's contribution to the surface temperature change and sum all the contributions together.Figure 11 shows the evolution of the PFC surface temperature change for times spanning the final loss event, where the colormap is the change in surface temperature in a log scale.Additionally, this figure splits the PFC surface heating into separate contributions from the initial RE beam in (a)-(d) and the secondary electrons in (e)-(h).
Figure 11 indicates that the PFC surface temperature change due to secondary electrons is more significant than the contribution due to the initial RE beam.However, this larger heating is also more transient, consistent with figure 2.
The tile-averaging also allows us to focus on the spatial localization of the PFC surface heating.Because B ϕ < 0 and J ϕ > 0, electrons with pitch angle η < 90 • will be traveling in the 'forward' − φ-direction and will intersect the inner wall on the left leading edge traveling from the Z < 0 region, consistent with the simulation results.Interestingly, the secondary electron population also has significant energy deposition on the right leading edge traveling from the Z > 0 region.Due to the broad pitch distribution of the avalanche distribution around η = 90 • and the greater effect of pitch angle scattering at their lower energies, many secondary electrons will also travel in the 'backward' φ-direction.The last observation in figure 11 is that the deposition of the initial RE beam is relatively localized near the midplane while the secondary electron deposition is vertically elongated.Because the initial REs have significant drift orbit effects, they can generate secondary electrons on open flux surfaces that terminate at different vertical positions along the inner wall.While the large drift orbit effects of the initial REs allow them to travel from confined to open flux regions, the secondary electrons generated at lower energies do not have significant drift orbit effects and follow the magnetic fields closely.
In the KORC simulations, we do not retain the initial energy of secondary electrons.Thus we cannot get statistics on the SEE population compared to secondary REs.However, we can calculate the relative proportion of SEEs to secondary REs, which allows us to infer that the SEE population dominates the wall heating due to secondary electrons.For the parameters yielding the secondary electron source plotted in figure 1, integrating the source over secondary electron phase space for the entire kinetic energy range from K therm to the maximum secondary kinetic energy yields 3.3 × 10 6 secondary electrons generated by a single RE per second.Integrating in the range starting from the approximate balance between thermalization and deconfinement time scales from figure 5 of K ∼ 5 keV yields 1.0 × 10 3 secondary electrons per second.Lastly, integrating in the range starting from K crit yields only 36 REs per second.Note that the balance between thermalization and deconfinement times scales and K crit used here will be the minimum values because we are considering the largest inductive electric field in figure 1.Thus, in the following calculations, only a maximum of 3.5% of secondary electrons generated by the avalanche source are REs, while the remaining 96.5% are SEEs.

Comparing modeled and experimental wall heating.
One of the main reasons we use DIII-D #177031 in this study is the availability of IR thermography during the final loss event.Figure 12 compares the IR camera data (looking opposite the forward-beamed direction) in the (c) with the KORC modeling results for the initial and secondary populations of electrons in (a), (b).The modeled PFC surface temperature is plotted by taking a single toroidal tile-averaged region from figure 11 and repeating it to make a better connection with the multiple toroidal tiles in the IR observation.We show IR data at a time later than the KORC results because an explosive event of one of the tile leading edges causes a dust ejection (not shown) that obscures the image during the final loss event.Because the modeling shows the change from the initial PFC surface temperature, comparing directly with the IR is challenging because the IR starts at different values before the final loss event.The tile edge starts at around 243 • C, as shown in figure 13, suggesting it would be more correct to have the cutoff be at 45 K, while the tile face starts at around 178 • C, suggesting the cutoff should be 110 K. Thus, any cutoff you choose will not give you a direct comparison to the IR.We used 100 K in the present analysis, but the resulting plots are similar for different cutoffs.The IR observation indicates a localized impact near the midplane and a lower magnitude of heating in elongated bands vertically along the toroidal leading edges, qualitatively agreeing with KORC results only when including the avalanche source of secondary electrons.Moreover, it is also possible to discern from the IR camera image that there are both forward and backward leading-edge deposition, notably in the middle leading edge and the parity between Z > 0, < 0 is consistent with KORC modeling.There is a minor discrepancy in the spreading of the localized impact site in the KORC calculations.Because the integration time and frame rate of the IR camera was 0.2 ms and 1004.62Hz, while the experimental reconstructions are spaced by 0.5 ms, we posit that the time resolution of experimental reconstructions leads to this vertical spreading of the modeled RE/SEE deposition during horizontal advection of the RE beam during the final loss event.We envision that the initial REs would be more locally deconfined vertically along the inner wall, but snapshots would lead to rapid deconfinement that may be vertically elongated.Follow-up work will aim to use more highly resolved reconstructions or interpolations of the field configuration between snapshots.
In addition to the spatial structure of the PFC surface heating at a fixed time, we also compare the PFC surface temperature evolution from KORC and IR imaging at fixed locations.Figure 13 shows the time evolution of the calibrated IR data from DIII-D #177031 (blue and red traces) and modeled KORC PFC surface temperature (green trace).The blue trace is the temperature at a PFC edge with location (ϕ, Z) = (109.21• , −0.1221 m), the red trace is at a PFC face with location (ϕ, Z) = (105.67• , −0.1226 m), and the grayed region is where the IR signal is saturated.The locations on the wall are indicated by the symbols of matching color in figure 12(c).The KORC modeled value is the evolution of the toroidally tile-averaged PFC surface temperature change at the location of the maximum PFC temperature change at t = 1.6012 s, (ϕ, Z) = (0.492 • , −0.293 m) is close to the leading edge of the forward direction in the coordinate frame of figures 11 and 12.While the KORC modeling results predict a significant transient change in the PFC surface temperature that reaches the approximate DIII-D graphite tile sublimation threshold of 1000 • C-2000 • C [34], the subsequent modeled thermal transport occurs much more rapidly than observed  12(c).The temperature using KORC RE trajectories reaches the approximate heating range needed for graphite ablation, but only transiently before decreasing rapidly compared to DIII-D observations.Additionally, there is no simulated temperature change on the tile face because all modeled REs deposit into the leading edges and the 1D heat transport model assumes constant specific heat and conductivity.
experimentally.Additionally, because the 1D heating model does not consider the 2D thermal diffusion, the model cannot predict heating at the tile face.We also note that the confinement loss during the RE final loss event will also deconfine the remaining companion plasma, which is not considered in the current model.
While there are other sources of wall heating during the RE final loss, we argue that RE energy deposition is the dominant source.The RE beam kinetic energy can be roughly estimated as W kin ∼ KI p 2π R/ec [13,39,40] where I p = 250 kA and K ∼ 1 MeV, giving W kin ∼ 10 kJ.This energy will be deposited locally at the plasma-wetted area, which the IR camera observed.The RE beam magnetic energy is more significant, and is estimated as W mag ∼ 0.5LI 2 p ∼ 50 kJ [13,40], where L p is the plasma inductance.This will primarily go into line radiation following the RE loss to the wall [28].The line radiation will uniformly heat the center post, as compared to the localized RE deposition.Much of the current will transfer into the toroidal wall current and heat the vacuum vessel wall inside the tiles on the center post, which the IR camera does not observe.The background plasma thermal energy W th ∼ 3T e n e V p is small because T e is only 2 eV, W th ∼ 1 kJ.Additionally, this energy will likely stay in the plasma.After the REs are lost, the remaining thermal plasma is only ∼1 eV, having a low parallel thermal conductivity, and most of this energy will be lost to radiation.There is little previous work on the halo currents of the RE final loss event.However, the halo currents are a fraction (∼1/10) of the toroidal current, so we estimate the dissipated energy is 5 kJ or lower.During the current decay following the RE loss, the plasma is ∼1 eV, so the resistivity is much larger than the wall resistivity.This can be observed by looking at the I p decay, which is ∼1 ms compared with the wall time ∼8 ms.Thus, the halo current energy is primarily dissipated in the plasma, going into line radiation.The small fraction dissipated in the wall will volumetrically heat the conducting wall, not give localized leading edge heating.
The heating due to the final RE loss can also be compared to the wall heating during the thermal and current quenches.As discussed, the present experiment uses a small, cryogenic Ar pellet to trigger a disruption and form a RE beam [41,42].This is different than the planned DMS for ITER, which uses primarily low-Z material with a small fraction of a high-Z impurity.This method is utilized because, in present experiments, a RE beam does not form unless a large amount of high-Z material is included.This results in roughly >90% of the pre-disruption thermal energy radiating via line radiation [43,44].In DIII-D 177 031, the pre-disruption internal energy is ∼ 100 kJ.Thus, ∼90 kJ will volumetrically heat the entire wall surface due to radiation, and ∼10 kJ will heat the wall due to the advection of thermal energy.Predicting the heat loads due to the advection of thermal energy is challenging because it depends on impurity transport and MHD response.Roughly, the MHD response will lead to large wetted areas on the wall compared to localized RE strikes.The pre-disruption internal energy is significantly lower than experiments discussed in [44] that had internal energy of 0.78 MJ and resulting wall temperatures under 300 • C at the divertor inner strike point following the thermal quench.The pre-disruption current is ∼1.2 MA making the magnetic energy ∼1 MJ, assuming 10% is the halo current.Due to the ∼1 eV temperature of the plasma after the thermal quench and the impurities from Ar injection and carbon from DIII-D tiles, the dissipation of magnetic energy via line radiation will follow similar processes to those acting during the final RE loss event [36,44].However, the significantly increased magnetic energy during the current quench, as compared to the final RE loss event, will likely lead to observable heating.A detailed analysis of the wall heating similar to that in, e.g.[45] would yield insight into this process.

Summary and conclusions
In this manuscript, we discuss the contribution of secondary electrons generated via large-angle collisions comprising the avalanche source to wall heating during RE final loss events, particularly the role of SEEs.In section 2.1, we present the form of the avalanche source and, with the energy-dependent PFC surface heating model for a single particle equation ( 12) in section 2.2, describe how the SEE population at energies below the runaway threshold play a dominant role in the PFC surface heating.An example using parameters from DIII-D discharge #177031 was shown in figure 1 to illustrate this finding.In section 3, we discuss the final loss event in DIII-D discharge #177031 that we use to explore wall heating by SEEs. Figure 3 shows that the RE beam advects horizontally and scrapes off on the inner wall limiter, which induces a large toroidal electric field.Additionally, in this section, we approximate the inner wall comprised of individual graphite tiles as a regular polygon to model the inhomogeneous RE and SEE energy deposition onto tile leading edges.The following appendices provide details on implementing and benchmarking the large-angle collision operator in KORC.Results using the KORC code began by performing calculations without large-angle collisions in section 4.1 to estimate the neutral content that subsequent calculations use.Per separate collisionless and spatially-independent calculations, we found the kinetic energy where collisional and deconfinement effects were roughly equal during the final RE loss event in DIII-D #177031, implying the lower energy required to model wall heating effects due to SEEs adequately.In section 4.2, we discuss the evolution of KORC simulations of DIII-D discharge #177031, where figure 6 showed the generation of many secondary electrons during the final RE loss event.During the initial portion of the final RE loss event, most of the secondary electrons thermalize while field lines are still closed and confining.However, near the end of the final RE loss event, magnetic fields open to the wall and secondary electrons are rapidly lost, consistent with figure 8.We then extend the single particle heating model to incorporate a spatially-dependent ensemble of electrons equation ( 17) in section 4.3 followed by providing a method for approximating the angle of incidence to the PFC surface equation (23) given the GC angle of incidence and pitch angle in section 4.3.1 for use in PFC surface heating calculations.PFC surface heating calculations are performed in section 4.3.2 that find secondary electrons generated by the avalanche source were responsible for most of the heating on the leading edge of PFC tiles and that the surface temperature change was of the order needed for the sublimation of graphite tiles in DIII-D.Lastly, in section 4.3.3,we found qualitative agreement with IR imaging of PFC in DIII-D discharge #177031, showing the forward and backward electron deposition along leading edges.
Moving toward the reactor level plasma currents planned for ITER, SPARC, and other future tokamaks, optimizing RE final loss events to disperse RE deposition while limiting the generation of SEEs and better understanding how the RE beam interacts with the first wall will be required to ensure successful RE mitigation.We first comment that the present work on PFC surface heating by SEEs could contribute to the low PFC surface heating observed in [46,47], where low-Z material injection into a RE beam precipitates unstable MHD modes.In addition to the MHD mode spreading out the deposited electrons onto a larger wetted area, the final RE loss events occur over microsecond time scales rather than the millisecond scales in RE mitigation scenarios with high-Z material injection, which would significantly limit the generation of additional SEEs.Second, part of the ITPA DivSOL/MDC collaboration B1.3: characterization of power deposition to PFCs by runaway electrons is using the methodologies of the present study.These studies use KORC and analytic first wall capabilities to track REs in 3D magnetic fields from MARS-F [48] and simulate where they are deposited into a graphite dome inserted into DiMES [49,50] during a controlled RE impact.The resulting spatial, energy, and pitch distribution of deposited electrons will be the initial conditions of calculations using the successor to the MEMOS-U code [51] to compare against experimental measurements.
Lastly, we discuss the parametric trends of the secondary electron surface heating to assess the scalability of the simulated DIII-D observations to reactor-size devices.The RE beam current scales proportional to the plasma current [39].The RE current is proportional to the number of particles, while the velocity and pitch should be relatively unchanged.This would lead to a proportional increase in the number of secondaries generated.The RE beam current scales proportional to the magnetic field strength [39].This would lead to a proportional increase in the number of secondaries generated.Also, higher energy REs remain confined for increasing magnetic fields rather than being lost to drift orbit effects [37,39].However, 10 MeV to 20 MeV initial RE energy with the same lower threshold in the avalanche source leads to only a ∼0.8% increase in the number of secondaries generated.Scaling with size, deconfinement time increases while thermalization stays the same, leading to a higher lower bound in surface-heating secondaries and more trapped secondaries.This, coupled with a larger first wall due to the size increase, ultimately leads to less surface heating.Nuclear operation in a reactor provides a consistent source of primary REs from the activated wall and tritium [52] throughout the RE final loss event.While the PFC shaping includes bevels to avoid leading edges for steady-state heat fluxes [53], these bevels create leading edges for RE final losses that are directed forward and backward relative to the magnetic field.
To determine at which (p 1 , η 1 ) to insert the new RE, we employ a 2D rejection sampling algorithm.The sampling begins by choosing two uniformly distributed random numbers that are scaled to have trial values satisfying (p 1,t , η 1,t ) ∈ ([p min , p i ], [0, π]).Then, we check that these trial values yield a positive argument to the radical in the denominator of the RHS of equation ( 4), indicating that the trial values are physically possible; if not, we generate new trial values.Next, we use (p 1,t , η 1,t ) to calculate S Aval.from equation (7).The ratio of the trial calculation of S Aval. to its maximum value over all (p 1 , η 1 ) is then compared to another uniformly distributed random number.The algorithm cycles until the ratio is larger than the random number.
Recall that equation ( 7) is the change in the RE distribution function due to the addition of secondary REs due to largeangle collisions of a single initial RE.However, the RE distribution will also change because the energy and momentum of the initial RE must be conserved in the large-angle collision creating a secondary RE [11].We account for this in our Monte Carlo implementation by enforcing that For the present case, where the thermal background is cold with T e ∼ eV, we assume that the secondary RE initially has insignificant energy and parallel momentum.In practice, this is a valid approximation because p min is over two orders of magnitude larger than p th (v th ).
Because the avalanche source is exponentially growing, the Monte Carlo implementation often leads to a problem with the load balancing of simulated REs across a distributed memory parallelization.We initialize calculations with approximately 1/20 the maximum number of simulated REs on a given computational node to combat this issue.Once the avalanche source has increased the number of REs on any node to the maximum value, we stop the calculation.The calculation restarts with additional nodes, where the REs from the previous calculation are reordered and evenly spread across the new nodes.

Appendix B. Large-angle collision operator benchmarking
Next, we discuss verifying the implementation of the largeangle collision operator.We begin by running KORC without updating the spatial information of the REs but retaining the effects of electric field acceleration, synchrotron radiation, collisional friction, pitch angle scattering, and the large-angle collision operator.We start calculations with a monoenergetic/monopitch distribution with 10 MeV/10 • and let the distribution functions evolve until ∼30τ c . Figure 14 shows an example of such a calculation in the left plot of for E/E CH = 5, τ c /τ R = 0.5, Z eff = 2, and ln Λ = 20.Because these calculations ignore spatial effects, the electric and magnetic fields and electron density can be chosen arbitrarily as long as the ratios for E/E CH and τ c /τ R are maintained.The solid blue line shows the evolution of the total number of REs evolved in the calculation, where the solid red line shows the number of REs in the calculation that have not been thermalized (i.e.satisfy p > p min ).We find the RE avalanche growth rate by taking an exponential fit to the last 75% of the calculation, which the dashed lines show.The middle plot of figure 14 shows the excellent agreement of KORC results with the benchmark in figure B3 of [18], where the electric field was varied to find the critical electric field corresponding to zero RE avalanche growth rate.
We perform a second benchmark to verify the proper computation of toroidicity and particle trapping effects on the avalanche source [56,57] with KORC.The right plot of figure 14 shows KORC guiding center calculations similar to those done in [55], with a major radius R 0 = 6 m, an aspect ratio a/R 0 = 1/3, a toroidal field at the magnetic axis B 0 = 53 T, and a poloidal magnetic field with safety factor q = 2.1 + 2(r/a) 2 .The toroidal electric field scales as E ϕ ∼ 1/R for a constant loop voltage.The electron density is n e = 1 × 10 20 m −3 , electron temperature T e = 2 eV, Z eff = 1, and Coulomb logarithm lnΛ = 10 are all uniform and constant.We distribute the initial REs on a disk for given r/a with uniformly distributed energy ∈ [5,10] MeV and pitch ∈ [0, 36.7]• .We choose the toroidal electric field to have a value near and much larger than the critical electric field, and the radial location of the initial REs is varied.Results from KORC compare well with those from [55].Results of calculations ignoring spatial evolution, using electric field acceleration, collisional friction, pitch angle scattering, synchrotron radiation, and large-angle collision source.Calculations are run for >30τc and converge to a steady growth rate, (a) shows the evolution of all (blue trace) and a subpopulation of active (red trace) REs included in a calculation with E/E CH = 5, where the dashed lines show an exponential fit to the last 75% of the evolution.The results compare favorably with [18,54] as the electric field is varied (b), using τc/τ R = 0.5 and ln Λ = 20.The toroidal and trapped particle effects on the avalanche growth rate are also consistent with simulations in [55] (c), capturing the physical effects discussed in [56,57].Data from [18,54,55] is reconstructed using an online program, WebPlotDigitizer [58].

Figure 1 (
a) shows a colormap of log 1 0(S Aval. ) for an initial RE with kinetic energy K 0 = 10 MeV and pitch angle η 0 = 10 • indicated by the red dot.We also choose plasma and impurity values representative of the post-disruption background plasma in DIII-D #177031, to be discussed in section 3, with n e = 2.25 × 10 20 m −3 , n +1 D = n +1

Figure 1 .
Figure 1.(a) Avalanche source due to large-angle collisions in equation (7) is shown by colormap bounded by the solid black trace for a chosen initial RE energy and pitch at the red dot.The vertical dash-dotted and dashed black lines indicate the critical and maximum secondary RE energies and the vertical dotted line indicates the thermal energy of the cold companion plasma.The red trace indicates the approximate runaway region from[19] using the critical RE energy in the presence of impurities from[8].(b) The red trace (left axis) is the maximum change in graphite PFC surface temperature multiplied by unit area and unit time as a function of single-particle electron kinetic energy according to equation(12).The blue trace multiplies the PFC surface temperature change by the differential of the density from the secondary RE source to quantify the energy dependence of the heating effect of the secondary RE source.The green trace (right axis) shows an analytic expression of the CSDA range[20] found by fitting the high energy ESTAR[21] and low energy IMFP[22] databases.

Figure 2 .
Figure 2. Evolution of surface ∆T for single electrons with different energies.

Figure 3 .
Figure 3.Time interval of the RE final loss event in DIII-D discharge #177031 simulated by KORC.In (a), the blue trace indicates the measured plasma current, the red trace indicates the average midplane density measured by the midplane interferometer, and the green trace indicates the maximum toroidal electric field as determined from the time evolution of JFIT reconstructions.Plots (b)-(d) show the toroidal electric field normalized to the critical electric field from equation (13) in the colormap with overlaid contours of poloidal flux before (b), during (c), and after (d) the final loss event.

Figure 4 .
Figure 4. Schematic of the periodic inner wall used in KORC simulations of DIII-D.The inset figure displays a single facet of two adjacent tiles, where the solid, red line is the wall and the dotted, blue trace is the circumradius.

Figure 5 .
Figure 5.Time scale of thermalization (blue traces) and deconfinement (red traces) as a function of initial kinetic energy using KORC simulation parameters for DIII-D discharge #177031.Solid traces/dots indicate an initial pitch angle of η = 10 • in dotted traces/open dots indicate an initial pitch angle of η = 80 • .Vertical violet lines indicate where thermalization and deconfinement time scales are equal.

Figure 6 .
Figure 6.(a) Evolution of normalized current.The red (blue) traces indicate contributions due to initial (secondary) electrons in a KORC simulation, and the black trace indicates the total.The violet trace corresponds to the plasma current in DIII-D discharge #177031.(b) The evolution of the number of secondary electrons generated each time step is indicated by the black trace, total thermalized (solid trace) and rapidly thermalized (dotted trace) by the blue traces, and lost (solid trace) and rapidly lost (dotted trace) by the red traces.The green trace indicates the maximum E ϕ from figure 3.

Figure 7 .
Figure 7. Evolution of kinetic energy (a), (c) and pitch angle (b), (d) distributions of initial RE (a), (b) and secondary electrons generated by the avalanche source (c), (d).

Figure 8 .
Figure 8. Histograms showing the spatial distribution of the secondary REs born during each time interval of the experimental reconstructions.

Figure 9 .
Figure 9. Schematic of the PFC geometry with an incident RE.The normal to the PFC is R and the guiding center of the RE intersects the PFC approximately along − φ ′ .The blue helix illustrates the electron motion before impact, while the purple vector v indicates the electron motion at the time of impact.The difference between the angle of incidence Θ and the guiding center angle of incidence θ GC depends on pitch angle and gyro-phase, which we choose randomly in the present analysis.

Figure 10 .
Figure10.Distribution functions of the sine of the angle of incidence, guiding center angle of incidence, pitch angle, and gyro-phase for all (blue traces), initial (red traces), and secondary (green traces) REs.Note that the (a) normalizes the distribution functions to the total number of REs.

Figure 11 .
Figure 11.Change in PFC surface temperature calculated by equation (17) using RE trajectories from KORC simulation.We split REs into initial (a)-(d) and secondary (e)-(h) populations, and average over the symmetry of the regular polygon in the toroidal direction.

Figure 12 .
Figure 12.Comparison of temperature change using RE trajectories from KORC simulation (a), (b) with IR on DIII-D discharge 177 031 (c).The color map for KORC temperature change has a ceiling of 100 K to compare with the saturated IR image from DIII-D.Vertically extended lobes directed forward and backward relative to the RE beam are qualitatively consistent with heating due to secondary REs.In (c), the blue and red symbols indicate the locations where the evolving temperature is shown in figure 13.

Figure 13 .
Figure 13.Comparison of the evolution of PFC surface temperature at facet edge (blue trace) and face (red trace)from IR on DIII-D discharge 177 031 with the evolution of maximum temperature change using RE trajectories from KORC simulation.The locations on the wall are indicated by the symbols of matching color in figure12(c).The temperature using KORC RE trajectories reaches the approximate heating range needed for graphite ablation, but only transiently before decreasing rapidly compared to DIII-D observations.Additionally, there is no simulated temperature change on the tile face because all modeled REs deposit into the leading edges and the 1D heat transport model assumes constant specific heat and conductivity.

Figure 14 .
Figure14.Results of calculations ignoring spatial evolution, using electric field acceleration, collisional friction, pitch angle scattering, synchrotron radiation, and large-angle collision source.Calculations are run for >30τc and converge to a steady growth rate, (a) shows the evolution of all (blue trace) and a subpopulation of active (red trace) REs included in a calculation with E/E CH = 5, where the dashed lines show an exponential fit to the last 75% of the evolution.The results compare favorably with[18,54] as the electric field is varied (b), using τc/τ R = 0.5 and ln Λ = 20.The toroidal and trapped particle effects on the avalanche growth rate are also consistent with simulations in[55] (c), capturing the physical effects discussed in[56,57].Data from[18,54,55] is reconstructed using an online program, WebPlotDigitizer[58].