The Impact of Collisionality on the Runaway Electron Avalanche during a Tokamak Disruption

The exponential growth (avalanching) of runaway electrons (REs) during a tokamak disruption continues to be a large uncertainty in RE modeling. The present work investigates the impact of tokamak geometry on the efficiency of the avalanche mechanism across a broad range of disruption scenarios. It is found that the parameter $\nu_{*, crit}$ describing the collisionality at the critical energy to run away delineates how toroidal geometry impacts RE formation. In particular, utilizing a reduced but self-consistent description of plasma power balance, it is shown that for a high-density deuterium-dominated plasma, $\nu_{*, crit}$ is robustly less than one, resulting in a substantial decrease in the efficiency of the RE avalanche compared to predictions from slab geometry. In contrast, for plasmas containing a substantial quantity of neon or argon, $\nu_{*, crit} \gtrsim 1$, no reduction of the avalanche is observed due to toroidal geometry. This sharp contrast in the impact of low-versus high-Z material results primarily from the relatively strong radiative cooling from high-Z impurities, enabling the plasma to be radiatively pinned at low temperatures and thus large electric fields, even for modest quantities of high-Z material.

The exponential growth (avalanching) of runaway electrons (REs) during a tokamak disruption continues to be a large uncertainty in RE modeling.The present work investigates the impact of tokamak geometry on the efficiency of the avalanche mechanism across a broad range of disruption scenarios.It is found that the parameter ν * ,crit describing the collisionality at the critical energy to run away delineates how toroidal geometry impacts RE formation.In particular, utilizing a reduced but selfconsistent description of plasma power balance, it is shown that for a high-density deuterium-dominated plasma, ν * ,crit is robustly less than one, resulting in a substantial decrease in the efficiency of the RE avalanche compared to predictions from slab geometry.In contrast, for plasmas containing a substantial quantity of neon or argon, ν * ,crit ≳ 1, no reduction of the avalanche is observed due to toroidal geometry.This sharp contrast in the impact of low-versus high-Z material results primarily from the relatively strong radiative cooling from high-Z impurities, enabling the plasma to be radiatively pinned at low temperatures and thus large electric fields, even for modest quantities of high-Z material.

I. INTRODUCTION
Disruptions pose a significant threat to the success of the tokamak concept in magnetic fusion.During a tokamak disruption, nearly all of the plasma's thermal energy is lost, leading to a drastic increase in the plasma resistivity and a subsequent spiking of the electric field.The force exerted on an electron by the resulting electric field can greatly exceed the electron's collisional drag, thus enabling the acceleration of electrons to relativistic energies 1 .
The number of these so-called runaway electrons (REs) can in turn be amplified by a process referred to as 'avalanching' 2 .Here, a pre-existing RE undergoes a large-angle collision with a low-energy electron.If the energy transferred during this collision is sufficient to scatter the initially cold electron to high enough energy to run away (often referred to as a secondary electron), this will lead to the growth of the initial RE population.This process has been observed in existing tokamaks to enable the transfer of plasma current to the RE population, with the potential for mega-Ampere RE currents being inadvertently generated in reactorscale tokamak plasmas 3 .Moreover, the energy of RE beams is typically in the MeV range, which implies the possibility of subsurface damage to plasma-facing components, further raising the importance of RE mitigation.
Understanding the efficiency of the avalanche mechanism is important for predicting the amount of RE current that can be expected to form during a tokamak disruption.While the impact of partially ionized impurities on avalanche generation has been extensively studied [4][5][6] , substantial uncertainties remain regarding the impact of tokamak geometry on the avalanche mechanism.In particular, the magnetic field strength in a tokamak device scales as B ∝ 1/R, where R is the major radius.Electrons with modest values of pitch ξ ≡ p ∥ /p, are susceptible to magnetic trapping, leading to the formation of banana orbits.In the presence of a parallel electric field, electrons undergoing banana orbits will be accelerated on one leg of the orbit and decelerated on the returning leg, leading to a negligible acceleration of the electron after a complete banana orbit.Since a substantial fraction of secondary electrons are born with small values of the pitch ξ ≡ p ∥ /p, the presence of a trapped region can substantially reduce the efficiency of the avalanche mechanism 3,[7][8][9] .
A caveat to the above picture, however, is that it requires the electron to have sufficiently low collisionality such that the electron is able to complete a full bounce orbit before being collisionally de-trapped.While MeV electrons are expected to be in the collisionless regime for parameters characteristic of a tokamak disruption, the collisionality at the critical energy for an electron to run away is expected to be far larger.In particular, defining the critical energy for an electron to run away by the energy at which the drag and electric field terms balance, this energy can be approximated by m e v 2 crit /2 ∼ (m e c 2 /2) E c /E ∥ , where we have assumed E ∥ /E c ≫ 1, E c ≡ m e c/(eτ c ) is the Connor-Hastie electric field 1 , and / (e 4 n e ln Λ) is the relativistic collision time.During a tokamak disruption, this energy can obtain values ranging from several to tens of keV, depending on the specific plasma composition present during the disruption.Noting that the collisionality scales inversely with energy, electrons near v crit often will not be in the banana collisionality regime, potentially negating electron trapping effects on the avalanche growth rate 10 .
The present paper aims to employ an idealized but self-consistent model of a tokamak disruption to identify the collisionality regime characterizing runaway electron formation across a broad range of disruption scenarios, with particular emphasis on the avalanche mechanism of RE formation.Specifically, defining the dimensionless collisionality parameter by ν * ≡ τ bounce /τ dt , where τ bounce = qR 0 /(ϵ 1/2 v) and τ dt = ε/ν D (v) are the bounce and collisional de-trapping times, respectively 11 .The collisionality ν * can be further expressed as where a is the minor radius, R 0 is the major radius, ε = r/R 0 is the inverse aspect ratio, q is the safety factor, and ν D (v) is the pitch-angle scattering rate.Noting that cτ c /a ∝ 1/(n e ln Λ), where n e is the free electron density and ln Λ is the Coulomb logarithm, the collisionality at the critical speed v crit for an electron to run away scales as where we have omitted the geometric parameter qR 0 / aε 3/2 for simplicity.Here, ν * ,crit depends sensitively on the plasma density, pitch-angle scattering frequency, and v crit .During a tokamak disruption, these three quantities can span a broad range of values, depending on the disruption mitigation scheme pursued and the device under consideration.For example, the large inductive electric field present during a disruption is the result of a sharp decrease in the plasma temperature T e driving an abrupt rise in the plasma resistivity.The decrease in T e is dependent on the plasma composition (n impurity , n D ), which directly relates E ∥ to the quantity and species of atoms injected into the plasma as part of a disruption mitigation scheme or that is released from the wall due to plasma-material interaction.The primary motivation of this paper will thus be to quantify how distinct plasma compositions impact the value of ν * ,crit and thus the efficiency of the RE avalanche.
The remainder of this paper is organized as follows: Section II describes the self-consistent formulation for evaluating ν * ,crit and evaluates it across a broad range of plasma composi- During the thermal quench (TQ) phase of a tokamak disruption, where large regions of magnetic field stochasticity are expected, the plasma primarily cools by heat conduction along the open magnetic field lines to a temperature of roughly a hundred eV 12 , after which radiation further cools the plasma.Losses from this latter mechanism will eventually be balanced by the Ohmic heating of the resistive plasma, resulting in a saturated post-TQ plasma temperature T e .In order to self consistently evaluate T e for a given plasma composition, a power balance model must be employed.Here, data from the collisional-radiative models FLYCHK 13 and ATOMIC 14 are used, where we consider a steady state 0-D power balance between the Ohmic heating of the bulk plasma (j 2 Ω ) and the radiative cooling (S rad ).In addition to radiative losses, these atomic data sets are used to evaluate the average charge state of the impurity Z impurity (neon or argon) and deuterium Z D ions, and n e is evaluated through quasineutrality.Two representative cases are shown in Fig. 1 for singly charged neon [panels  (d)], the ∼ 1 eV root enables the deuterium species to begin to recombine, thus defining a critical deuterium density above which the plasma is expected to be in a low ionization state.For these high deuterium density cases, where multiple roots between ∼ 1 eV and ∼ 10 eV exist, the larger of the two roots is used since we are considering scenarios where the plasma is cooling from higher temperatures.If the deuterium density is further increased or if the current density is reduced, the higher temperature root will be removed, leaving only the root at one to two eV.
Identifying the precise criteria for accessing this low-temperature root requires the incorporation of physics beyond the idealized power balance analysis discussed above.In particular, the power balance curves shown in Fig. 1 were evaluated under the assumption of an optically thin plasma.While this is expected to be a good approximation for the majority of plasma compositions under consideration, recent work has noted that for sufficiently high deuterium densities, the plasma will become opaque to Lyman radiation, leading to a substantial reduction of radiative losses S rad for deuterium 15 .While this reduction of radiative losses postpones the density at which deuterium recombines, this effect only impacts plasmas with exceptionally high deuterium densities, which make up a small portion of the cases treated in this analysis.We will thus employ the assumption of an optically thin plasma throughout this paper, where finite opacity effects will be left for future work.

B. Electric Field Evaluation
Once T e , n e , and Z impurity are evaluated across a range of plasma compositions, E ∥ is obtained from Ohm's law: E ∥ = η s j ∥ , where η s is the Spitzer resistivity 16 and j ∥ is the local current density.Here we are interested in scenarios during the early phase of a disruption where the RE current j RE is negligible, hence we neglect its impact on Ohm's law. Figure 2 shows the electric fields evaluated from the power balance for argon and neon impurities and two different values of j ∥ .The yellow, orange, and red contours are for T e = 10, 17, and 20 eV, respectively.While the current decay time in ITER will depend on a range of factors 17 , plasmas with temperatures of approximately ten to twenty eV will have current decay times near the ITER target range of 50-150 ms 18 , and will thus be the primary focus of this analysis.We can understand the behavior of second large E ∥ /E c region at large deuterium densities is no longer present since the critical deuterium density at which the plasma recombines is increased.The impact of the choice of impurities on E ∥ /E c is also shown in Figure 2. Comparing Figures 2 (a) and (c), it can be seen that for argon, the T e = 10-20 eV channel is generally at lower deuterium and argon densities than neon, and the width is narrower for argon.As a result, argon impurities produce larger E ∥ /E c in the region where 10 eV < T e < 20 eV, since E c ∝ n e .Furthermore, argon impurities globally increase E ∥ /E c across all deuterium and argon densities, resulting from more efficient radiative cooling than neon.Specifically, comparing Figures 1 (a) and (c), the radiative cooling for argon is stronger than that for neon, where at a given deuterium and argon density, the temperature root is lower for argon, thus increasing E ∥ in comparison to neon.

C. Collisionality at the RE Critical Energy
The specific form used to evaluate the collisionality at the critical energy v crit is where the collisionality was evaluated at mid-radius (r/a ≈ 0.5), and the q profile used for this analysis is q = 2.1 + 2(r/a) 2 .Due to the strong dependence of the pitch-angle scattering rate ν D on v crit we will need to employ an accurate description of how v crit varies with the system parameters.While the energy at which collisional drag is balanced by electric field acceleration provides an order of magnitude estimate of the critical speed for electron acceleration, a more refined estimate can be made by identifying the X-point energy 19 .Specifically, defining the energy and pitch fluxes by where C F is the collisional drag, γ is the Lorentz factor, α ≡ τ c /τ s , and τ s = 6πϵ 0 m 3 e c 3 /(e 4 B 2 ) defines the synchrotron radiation timescale.The form of the collisional drag and pitch-angle scattering coefficients employed in this work, including partial screening effects, is given in Ref. 20 .For cases when the electric field is above the avalanche threshold, two nulls of the momentum space flux (i.e.Γ p = Γ ξ = 0) can be identified.The first is present at high energy and is distinct from, but correlated with, the bump in the steady state primary distribution 19,21 .This null in the momentum space flux, often referred to as an Opoint, represents the energy at which the combination of drag, synchrotron radiation, and pitch-angle scattering balance electric field acceleration and thus provides the characteristic energy of the saturated primary distribution.The second null, referred to as the X-point, occurs at much lower energy for E ∥ ≫ E av , and represents the separatrix between runaway and non-runaway orbits.Since the evaluation of the location of the X-point incorporates pitch-angle scattering and synchrotron radiation in addition to collisional drag, it provides a more accurate estimate of the threshold energy for electron acceleration compared to simply balancing collisional drag against electric field acceleration.We will thus use the X-point energy to estimate v crit in Eq. ( 3), where the analytic model derived in Ref. 22 will be employed to efficiently estimate the location of the X-point energy.Section IV A will utilize direct drift-kinetic simulations as a means of verifying this approximate analytic model.
The collisionality at the critical energy to run away ν * ,crit for plasmas with different compositions and j ∥ is shown in Figure 3.We see that ν is present, and we can anticipate that impurities robustly produce a collisional regime near the critical energy, thus reducing toroidal trapping effects.As the plasma composition approaches n D ≫ n impurity , ν * ,crit approaches ν * ,crit ≪ 1, suggesting that initially trapped electrons will complete many bounce orbits before collisionally de-trapping.We also see in the critical deuterium density for recombination is increased, and the recombination region is also shifted to higher n D .Larger j ∥ increases ν * ,crit across the entire range of plasma compositions, which arises from E ∥ ∝ j ∥ at a given plasma composition.

III. DRIFT-KINETIC SOLVER FOR RE SIMULATION
This section provides a brief description of the drift-kinetic code RAMc that will be used to directly evaluate the avalanche growth rate at different radial locations and for different plasma compositions [further details of the code can be found in Ref. 5 ].The governing electron description used is where time is normalized to τ c and the relativistic guiding center from Ref. 23 are used.Synchrotron radiation's effect on runaway physics 19,21 is accounted for in C rad (f ).The collisions Other parameters used were B = 5.3 T, q = 2.1 + 2(r/a) 2 , R 0 = 6 m, and a = 2 m. are modeled in C(f ), where C(f ) = C SA + C LA .Here C SA is the small-angle collision operator and C LA is the large-angle collision operator.The small-angle collisions are described by the monte-carlo equivalents of the Lorentz and energy scattering operators 5 .Large-angle collisions are captured through a Möller source term 24 , where both free and bound electrons are considered in the bulk electron density when generating secondary runaway electrons.

A. Runaway Probability Function
The runaway probability function (RPF) provides a convenient means of assessing how tokamak geometry modifies the efficiency of the avalanche mechanism across a range of plasma compositions.Specifically, the RPF [25][26][27][28] indicates the probability that an electron with a given initial phase space location will run away at a later time.By considering how this function varies across different plasma compositions this will provide direct insight into which regions of phase space contribute to the avalanche of runaway electrons.To evaluate the RPF, we will utilize the drift-kinetic solver described in Sec.III, but with the largeangle collision operator turned off, and a large number of marker particles used to sample the region of interest.Over time marker particles will either accelerate to relativistic energies, or slow down to the bulk plasma.By keeping track of the marker particle's initial location, and its final energy, the probability of electrons running away can then be computed.
The RPFs for two different plasma compositions and two different radii are shown in One immediate result of the RPF calculation is the verification of the critical energy expected reduction in the RPF is present for trapped electrons as r/a increases (compare the r/a ≈ 0 and r/a ≈ 0.8 cases).However, in the collisional regime, those electrons that would typically be trapped at ν * ,crit ≪ 1 are now not trapped, as shown in Figure 4 (d).In fact, the RPF is minimally impacted at ν * ,crit ≳ 1, highlighting that electrons in the trapped region indeed collisionally de-trap and become REs.These collisionally de-trapped electrons can become runaways, which is anticipated to enhance the avalanche efficiency.Moreover, in the ν * ,crit ≳ 1 regime, the energy where RPF = 0.5 and ξ = −1 decreases by over an order of magnitude, allowing more electrons to enter the runaway regime and generate secondary electrons via avalanche.
The impact of ν * ,crit can also be seen in spatial coordinates, where R and Z will be the radial and axial coordinates.In the collisionless limit, electrons can be described by their invariants f e ≈ f e (K, µ, p ϕ ), which are the electrons' energy K, magnetic moment µ, and toroidal canonical momentum p ϕ .For small ρ * e ≡ ρ e /a, where ρ e is the electron gyroradius, p ϕ becomes dominated by the poloidal flux function ψ.Therefore, the electron distribution function can be approximated as f e ≈ f e (γ, µ, ψ).The RPF at different ν * ,crit is shown in Figure 5 for argon impurities, where electrons are initialized near the critical energy and with µ = 0 (ξ ≈ -0.9 to -1.0 on the weak field side).In Figure 5 where R 0 is the major radius and E 1 is the electric field at the magnetic axis (R = R 0 ).The modified RPF structure thus follows E ∥ in the limit of ν * ,crit > 1, where electrons across the R-Z plane are largely passing.Since E ∥ is strongest on the inboard side, electrons thus have the largest probability of running away there, and the RPF is the largest.

B. Avalanche Growth Rate Evaluation
Quantifying the impact of ν * ,crit on the avalanche mechanism is done by directly evaluating the avalanche growth rate γ av as a function of r/a using the drift kinetic solver described in Sec.III.Here, a small initial seed population of REs is used to initiate the exponential  I.
growth.The avalanche growth rate is evaluated as a function of r/a at varying ν * ,crit for both argon and neon and is shown in Fig. 6, with the tabulated parameters in Table I.
In the collisionless (ν * ,crit ≈ 10 −2 ) region, the avalanche growth rate is reduced by almost 50% at r/a = 0.8.However, as ν * ,crit increases to the collisional regime, the reduction of the avalanche growth rate at large r/a becomes smaller, up to the point where it can be considered negligible, contrary to widely used formulations of the avalanche growth rate 3,7 .
The largest impact from impurities on the avalanche growth rate is between ν * ,crit ≈ 0.1 − 1, where the radial reduction varies anywhere from ∼ 40% to ∼ 10%.
One feature present is in the case of argon, where ν * ,crit ≈ 7. It is seen that the avalanche growth rate increases slightly with radius.The increased avalanching on the outboard side has been noted in the analysis of Dreicer RE generation in Ref. 10 , however here it is seen in the context of avalanche generation, albeit with only a very modest increase evident.This slight increase is due to the average electric field seen by electrons are large radii slightly exceeding the on-axis value.In particular, if the electric field E ∥ is expanded at the inboard for the outboard side.While the first order terms cancel, implying that the electric field an electron sees averages out as it does a poloidal transit, the second order terms do not.As a result electrons at large r/a can see an average electric field enhanced by a factor of (r/R 0 ) 2 .This feature can be seen by evaluating the average electric field an electron sees as it does many poloidal transits.Figure 7 illustrates the instantaneous and average electric fields an electron located at r/a = 0.8 sees as it completes several poloidal transits.Here, the oscillations in the instantaneous electric field seen by the electron are due to the 1/R dependence of the electric field; however, as indicated by the orange line, the average value of the electric field is slightly greater than one (⟨E⟩ /E 1 ≈ 1.024).While this enhancement of the electric field is modest, the RE avalanching deviates from a linear scaling with electric field strength for resulting in enhanced avalanching at large r/a, where the enhancement at r/a ≈ 0.8 was ∼ 4%.

C. Efficiency of the Avalanche Mechanism
A powerful means of estimating the maximum number of RE exponentiations possible for a given disruption scenario is by predicting the amount of poloidal flux that must be consumed to effect an order of magnitude amplification of a given seed RE population.Such an estimate directly measures the efficiency through which the avalanche process converts poloidal flux into increases in the RE population, often a more useful metric than specific values of the avalanche growth rate at given parameters.In particular, taking the avalanche growth rate to scale linearly with the electric field strength 3 , and assuming the electric field is much greater than the avalanche threshold, allows the avalanche growth rate to be expressed as: where γ exp is a constant coefficient that determines the efficiency of the avalanche.Noting that the inductive electric field is related to the change in the poloidal flux via the relation , where ψ is the poloidal flux function, Eq. ( 7) can be rewritten as Integrating this expression over the time interval t i to t f , yields the total number of exponentiations of the RE population, i.e.
where ∆ψ ≡ ψ (t f ) − ψ (t i ) and we have assumed the plasma composition to be unchanged such that γ exp / (E c R 0 ) is constant.The number of exponentiations of the RE population is thus directly linked to the change of the poloidal flux, the coefficient γ exp , and the parameters E c and R 0 .In addition, it is evident that the number of expected avalanche exponentiations does not depend on the specifics of the disruption trajectory, only the change of the poloidal flux so long as the coefficient γ exp /E c R 0 remains fixed.If we take the familiar expression for γ exp derived in Ref. 3 , , where φ ≈ (1 + 1.46 √ ε + 1.72ε), this allows Eq. ( 9) to be written as Here, we have introduced the quantity ψ exp defined by which is a constant that determines the number of exponentiations a RE population undergoes for a given change of poloidal flux.Equation ( 11) implies the efficiency of the RE avalanche is only sensitive to a modest number of quantities, primarily the Coulomb logarithm ln Λ, the effective charge Z ef f and the radial location through ε = r/R 0 .The major radius dependence indicated in Eq. ( 11) will cancel, due to the poloidal flux scaling with the major radius.For example, in circular geometry the poloidal flux function can be expressed in terms of the plasma current by where I p (r) indicates the plasma current inside a radius r.A convenient normalization of ψ exp will thus be ψexp ≡ 2πψ exp / (µ 0 R 0 ), which removes the major radius scaling, and has units of current.The quantity ψexp thus provides a rough measure of the drop of the plasma current required to effect one exponentiation of the RE population.A closely related, but more convenient quantity, is ψ10 = ln 10 ψexp , which is directly related to the drop in plasma current required to induce an order of magnitude increase of the RE population via avalanching 29 .Defining the Alfven current by I A ≡ 4πm e c/ (µ 0 e) ≈ 0.017 MeV, ψ10 can be expressed as: For a pure deuterium plasma, taking the RE beam to be tighly localized about the magnetic axis such that φ ≈ 1, and noting that the Coulomb logarithm seen by relativistic electrons will exceed the Coulomb logarithm for thermal electrons 30 , such that we may take ln Λ = 20, Eq. ( 13) implies ψ10 ≈ 0.94 MA, indicating that a tokamak with 15 MA of current, would have an avalanche gain of roughly ∼ 10 15/0.94 ≈ 10 16 , which is the often cited expected avalanche gain of ITER 17 .
While Eq. ( 13) provides a useful baseline for estimating the avalanche potential for a range of tokamak concepts, several uncertainties in the appropriate value of ψ10 remain.In particular, the toroidal factor φ implies that the avalanche growth rate for a given E ∥ can be reduced by roughly a factor of two at large minor radii.While it may often be the case that RE beams are largely localized to the inner half of the plasma, even a thirty percent decrease in γ exp due to toroidal corrections, would reduce the avalanche gain in ITER to 10 15/1.22 ≈ 2 × 10 12 , substantially reducing the importance of nuclear seeding mechanisms 31 .
As indicated in Fig. 6, this reduction due to toroidicity is only expected for plasmas with ν * ,crit < 1, which tend to be deuterium dominated plasmas.The expected ψ10 for ν * ,crit < 1  I.
plasmas is shown in Fig. 8 (blue and black curves) as a function of minor radius.Here it is evident that a substantial increase of ψ10 is present as r is increased, though the on-axis value is substantially less than the often invoked value of ψ10 ≈ 0.94 MA.This is due to Ref. 17 choosing a Coulomb logarithm of ln Λ = 20, which is larger than the Coulomb logarithm estimated from the more recent expressions derived in Refs. 20,30, particularly for the high deuterium densities that are required to maintain a plasma temperature of roughly ten to twenty eV in a plasma with modest impurity content.
In addition to toroidal corrections, Ref. 32 showed that the presence of partially ionized impurities led to a strongly nonlinear scaling of the avalanche growth rate, leading to a substantial increase in the efficiency of the RE avalanche at large electric fields.The presence of this nonlinear scaling of the avalanche growth rate prevents a unique values of ψ10 from being evaluated since the time integration indicated in Eq. ( 9) can no longer be analytically computed.Instead, the effective value of ψ10 will depend on the strength of the electric field, where larger electric fields will lead to a smaller value of ψ10 E ∥ .In Ref. 5 values near ψ10 E ∥ ≳ 0.5 MA were found on-axis for electric fields of tens of times the Connor-Hastie threshold.For the cases indicated in Table I, γ exp can be evaluated via the expression: where γ av is the avalanche growth rate previously evaluated.The resulting values of ψ10 E ∥ are shown in Fig. 8 for different minor radii.Perhaps the most striking feature in Fig. 8 lies in the region where ν * ,crit ≳ 1 for both argon and neon, resulting in ψ10 values as low as ψ10 ≈ 0.25.The highly efficient avalanche at these parameters arises from the non-linear dependence of the avalanche growth rate on the electric field.From Table I, the nonnormalized electric field ranged from 10 V/m to 60 V/m (hundreds to over a thousand times the Connor-Hastie threshold), which is well within the non-linear region of the avalanche growth rate 32 , resulting in the efficiency of the RE avalanche being substantially enhanced compared to Eq. ( 13).An additional striking feature is that depending on the value of ν * ,crit , ψ10 can vary from ψ 10 ≈ 0.25 − 1.1 at large minor radii, a far larger scatter in the efficiency of the RE avalanche than has been reported previously.

V. CONCLUSIONS
This work addresses the net impact of impurities and toroidal geometry on the RE avalanche mechanism for a range of parameters typical of tokamak disruptions.As RE avalanching is expected to be a large source of RE generation in future tokamaks such as ITER 17 , resolving uncertainties in this mechanism is imperative to identifying appropriate disruption mitigation strategies.The net impact from competing effects on the critical energy for RE generation produced by impurities is evaluated using a self-consistent power balance relation and an electric field estimated via Ohm's law.It was found that the strong radiative cooling resulting from substantial quantities of neon and argon impurities, and the subsequent large electric fields (see Fig. 2), leads to a decrease in the critical energy for an electron to run away.Consequently, electrons were found to be in the collisional regime (ν * ,crit ≳ 1) for cases where a substantial inventory of impurities was present (n impurity /n D ≳ 1).For these plasmas, the efficiency of the RE avalanche mechanism was found to be largely insensitive to magnetic trapping by the tokamak magnetic field.In contrast, as the plasma composition approaches a deuterium-dominant plasma, it is found that ν * ,crit ≪ 1.For this latter case the efficiency of the RE avalanche mechanism was found to drop as the radius was increased, consistent with previous results.The analysis in this work also shows that ν * ,crit serves as a guiding parameter to dictate the efficiency of the avalanche mechanism, in addition to indicating when toroidal trapping effects are prominent.
(a) and (b)] and argon [panels (c) and (d)] impurities, where Figs.1(a) and (c) represent a scenario where a significant amount of impurities are present (n impurity /n D ≈ 0.1) and Figs.1(b) and (d) represent a deuterium dominant plasma (n D /n impurity ≫ 1).Depending on the current density and plasma composition, stable roots may exist at ∼ 1 eV, ∼ 5-30 eV, and∼ 100-300 eV; however, we are only interested in roots with temperatures below the regime where heat conduction along stochastic magnetic field lines is subdominant.Hence, we will
n e and considering the case of neon at a given j ∥ and T e [Fig. 2 (b)].As n D decreases along a given temperature contour (say, the yellow 10 eV contour), E ∥ /E c increases significantly, which is due to the drop in the free electron density n e and the increase in Z ef f .Another feature present is a second large E ∥ /E c region at deuterium densities of n D ≈ 10 16 cm −3 in Figs.2(a) and (c).The corresponding temperature at these deuterium densities is roughly an eV, thus driving recombination and lowering n e which decreases E c .Once the current density is increased from 1 MA/m 2 to 2 MA/m 2 , this

FIG. 2 :
FIG. 2: Post-TQ electric fields, normalized to E c for neon [panels (a) and (b)] and argon [panels (c) and (d)], where (a) and (c) are for j ∥ = 1 MA/m 2 and (b) and (d) are for j ∥ = 2 MA/m 2 .Other parameters used were an inverse aspect ratio of ϵ ≡ a/R 0 = 1/3, a safety factor of q = 2.1 + 2(r/a) 2 , a minor radius of a = 200 cm, and a magnetic field of B = 5.3 T. The yellow, orange, and red contour lines represent T e = 10, 17, and 20 eV, respectively.The white contour represents electric fields below the Connor Hastie threshold electric field.
* ,crit peaks similarly to E ∥ /E c .As the plasma composition transitions from n impurity /n D ≈ 1 to n D ≫ n impurity in the region where 10 eV < T e < 20 eV (the red, yellow, and orange contours), ν * ,crit spans the entire collisionality regime.In the regime where n impurity /n D ≈ 1 the collisional regime (ν * ,crit ≳ 1)

Fig. 4 .
Fig. 4. One million electrons are initialized with energies near the critical energy obtained from the O-X merger model and ξ ∈ (−1, 1).Here, the black contours are the 50% contour and the white contours indicate the trapped-passing boundary.We have selected a plasma composition with a substantial quantity of impurities and thus a high ν * ,crit ≳ 1 [Figs.4(c) and (d)] and a deuterium-dominated plasma with ν * ,crit ≪ 1 [Figs.4(a) and (b)].
we see in panel (a) (ν * ,crit ≈ 10 −2 ) that the RPF contours indeed follow ψ, however, moving from panel (a) to panel (d), where ν * ,crit increases from 10 −2 to 7, the RPF structure changes significantly.Between panels (a) and (b), electrons continue to largely follow the flux surfaces, but as ν * ,crit increases from 10 −1 to 1 [panels (b) and (c)], the RPF structure departs from following flux surfaces.The modified RPF structure results from electrons on the outer flux surfaces collisionally de-trapping and becoming passing electrons.For larger ν * ,crit [panel (d)], electrons have the largest probability of running away on the inboard side, and the RPF structure now follows roughly a 1/R function, which results from the form of the electric field

FIG. 5 :
FIG. 5: RPF projected on to the R-Z plane for argon at varying ν * ,crit .10 6 maker particles were used and initialized with energy K = K crit and µ = 0 (ξ ≈ -0.9 to -1.0 on the weak field side.Simulation parameters can be found in Table I, and other parameters used are B = 5.3 T, j ∥ =1.5 MA/m 2 , q = 2.1 + 2(r/a) 2 , R 0 = 6 m, and a = 2 m.

FIG. 6 :
FIG. 6: Avalanche growth rates, normalized to the on axis value for neon (a) and argon (b).The dashed lines represent cases where synchrotron radiation was turned off, and B = 53 T. The solid lines represent cases where synchrotron radiation was kept on, and B = 5.3T.The full list of input parameters can be found in TableI.

FIG. 7 :
FIG. 7: Simulation of the toroidal electric field seen by an electron as a function of time, where time is normalized to τ c .Simulation parameters used here were B = 5.3 T, q = 2.1 + 2(r/a) 2 , R 0 = 6 m, and a = 2 m.The blue curve represents the electric field the electron sees at a given poloidal location normalized to E 1 , and the orange curve represents the average of this quantity over time.

FIG. 8 :
FIG. 8: Efficiency of the avalanche mechanism as a function of radius at different collisionalities for neon (a) and argon (b).The parameters are in TableI.

TABLE I :
Summary of parameters used for the drift-kinetic simulations in Figs.4,6, and 8 for a range of neon (a) and argon (b) densities, with j ∥ = 1.5 MA/m 2 .