Effects of magnetic perturbations and radiation on the runaway avalanche

The electron runaway phenomenon in plasmas depends sensitively on the momentum-space dynamics. However, efficient simulation of the global evolution of systems involving runaway electrons typically requires a reduced fluid description. This is needed for example in the design of essential runaway mitigation methods for tokamaks. In this paper, we present a method to include the effect of momentum-dependent spatial transport in the runaway avalanche growth rate. We quantify the reduction of the growth rate in the presence of electron diffusion in stochastic magnetic fields and show that the spatial transport can raise the effective critical electric field. Using a perturbative approach we derive a set of equations that allows treatment of the effect of spatial transport on runaway dynamics in the presence of radial variation in plasma parameters. This is then used to demonstrate the effect of spatial transport in current quench simulations for ITER-like plasmas with massive material injection. We find that in scenarios with sufficiently slow current quench, due to moderate impurity and deuterium injection, the presence of magnetic perturbations reduces the final runaway current considerably. Perturbations localized at the edge are not effective in suppressing the runaways, unless the runaway generation is off-axis, in which case they may lead to formation of strong current sheets at the interface of the confined and perturbed regions.


Introduction
Electron runaway is seen as one of the main threats to successful operation of magnetic confinement fusion devices with large plasma currents, such as ITER Breizman et al. 2019). The number of e-foldings in the runaway avalanche during a plasma-terminating disruption increases drastically when a tokamak is scaled up to ITER parameters from those currently in operation (Rosenbluth & Putvinski 1997). This calls for accurate models for the runaway generation and losses to ensure the design of a successful disruption mitigation system .
There is a wealth of experimental evidence that magnetic perturbations, occurring either naturally after a disruption or induced by external magnetic coils, can prevent or reduce runaway electron beam formation. In JET, a high level of magnetic fluctuations following a disruption has been seen to correlate with the absence of runaways (Gill et al. 2002). Broadband magnetic turbulence has been observed to lead to suppression of runaway current if the perturbation exceeds a certain level also in TEXTOR (Zeng et al. 2013) and in J-TEXT (Zeng et al. 2017). Kinetic instabilities driven by the runaways themselves can also induce local magnetic perturbations increasing the radial transport. Observations at DIII-D indicate that when the power in the instabilities exceeds a threshold, runaway plateau formation is absent (Lvovskiy et al. 2018). Perturbations imposed by external magnetic coils have also been shown to suppress the formation of runaway beams in several tokamaks (Yoshino & Tokuda 2000;Lehnen et al. 2008Lehnen et al. , 2009Mlynar et al. 2018).
The avalanche generation of runaway electrons is a result of momentum transfer between an existing runaway electron and a thermal one in a close collision. This leads to a growth of the runaway population that is proportional to the existing number of runaway electrons. Consequently, as the radial transport of runaways is also proportional to their number, it can reduce the growth rate of the exponentiation (Helander et al. 2000). Perturbations in the plasma confining magnetic field result in spatial transport and subsequent losses of runaway electrons (Rechester & Rosenbluth 1978). These losses reduce the number of runaway electrons participating in the avalanche mechanics and thereby have the potential to reduce the conversion of the initial plasma current to a runaway beam.
Modelling of a disrupting tokamak plasma resolved both in momentum, as needed for the runaway problem, and spatially, as needed to describe the evolution of plasma parameters, is computationally costly. Therefore, to follow the evolution of the disruption, simplified fluid models for the runaway populations are often used (Smith et al. 2006;Papp et al. 2013;Matsuyama et al. 2017;Bandaru et al. 2019;Fülöp et al. 2020). In these, the momentum space dynamics has been captured approximately, and an effective theory only dependent on spatially varying quantities is used to describe the growth and loss of the runaway population. Such simplified disruption modelling has been used to estimate the post-disruption runaway population in the presence of massive material injection (Martín-Solís et al. 2017;Vallhagen et al. 2020;Linder et al. 2020). However, these studies focused on the generation rates of the runaways and neglected the losses due to spatial transport.
The transport due to the perturbations is in general momentum dependent (Hauff & Jenko 2009), preventing a straightforward fluid description of the phenomena. The goal of this paper is to present a theory describing an effective rate of generation for runaway electrons which incorporates the effects of a momentum-dependent spatial diffusion. The diffusion considered here may originate from the motion of electrons in regions of stochastic magnetic fields as well as other perturbed magnetic field structures. In Sec. 3 we derive a self-consistent expression for the reduced avalanche growth rate of runaway electrons, including the effect of spatial transport, as well as radiation reaction forces and partially ionised impurities in a homogeneous plasma. Spatial variations in the plasma are investigated in Sec. 4 with a perturbation approach which conserves particle number to investigate the impact of radial transport in more realistic disruption simulations.
We find that, if the time-scale of the losses is comparable with that of the avalanche, spatial transport can raise the critical electric field for runaway generation significantly. The reason is that, even as runaway electrons are generated through close collisions by the avalanche dynamics, there need not be a net growth of the population if the relativistic electrons are transported out of the plasma. This may be part of the explanation of experimental observations which show strongly elevated critical electric fields (Martín-Solís et al. 2010;Hollmann et al. 2013;Granetz et al. 2014;Paz-Soldan et al. 2014;Popovic et al. 2016). The value of the critical field is also important for the dynamics in the current decay phase of disruptions, where the electric field tends to a value at which the loss and gain of runaway electrons is balanced (Breizman & Aleynikov 2017).
We demonstrate the effect of magnetic perturbations on runaway evolution in simplified disruption simulations in Sec. 5, taking into account the evolution and transport of runaways self-consistently with the electric field. We consider ITER-like plasmas with a combination of neon and deuterium injection and find that the runaway current can be suppressed, if the perturbations reach all the way to the plasma centre. The mixed magnetic topology common in disruptions is seen to have the potential to generate strong current sheets. Their stability may in turn be expected to impact the magnetic perturbation profile.

Radial diffusion of runaway electrons in the presence of radiation
Runaway electrons are almost collisionless and, as such, tend to follow magnetic field lines closely. Thus, in a stochastic magnetic field, the trajectories of runaway electrons generated close to one another will diverge with a rate dependent on the particle velocity along the field line and the rate of divergence of nearby field lines themselves (Rechester & Rosenbluth 1978). For a population of runaway electrons this leads to diffusive crossfield transport, the magnitude of which depends on the perturbation strength. However, at relativistic energies the electrons do not follow field lines closely, which causes the transport to decrease with increasing energy due to the effects associated with the finite orbit width (Myra & Catto 1992;Hauff & Jenko 2009;Särkimäki et al. 2020). Furthermore, it has been shown that modelling the transport as purely diffusive is insufficient in mixed magnetic topologies containing both islands and stochastic regions , but this can be addressed by including an advection term in the model (Särkimäki et al. 2016). A simplified theory to account for the momentum dependent radial diffusion in the avalanche growth rate was proposed by Helander et al. (2000), a theory which we will build on and extend to account for radiation reaction forces and the presence of partially ionised impurities. We address a case with mixed magnetic topologies in Sec. 5.
The momentum-space dynamics in the electron runaway problem is often described by the high-energy limit of the gyro-averaged kinetic equation with an accelerating electric field parallel to the magnetic field B (Hesslow et al. 2018b): Here, f is the electron distribution function, p is the momentum normalised to m e c, ξ = p·B/(pB) is the cosine of the pitch angle, E is the electric field strength normalised to the critical electric field E c = n e e 3 ln Λ c / 4πε 2 0 m e c 2 (Connor & Hastie 1975), where n e is the electron density, e the elementary charge, m e the electron mass, ε 0 the permittivity of free space, c the speed of light and ln Λ c 14.6 + 0.5 ln T eV /n e20 is the relativistic Coulomb logarithm, with T eV being the temperature measured in electronvolts and n e20 the electron density normalised to 10 20 m −3 . The relativistic collision time between electrons is τ = m e c/ (eE c ), C{f } is the relativistic collision operator and the last term on the right hand side of (2.1) represents radiation reaction forces from synchrotron radiation and bremsstrahlung. The relativistic test-particle collision operator is given by (Helander & Sigmar 2005) C{f where ν s and ν D are the slowing down and deflection frequencies, respectively. For relativistic electrons ν s and ν D take the form where for the case of a fully ionised plasmaν s = 1 andν D = 1 + Z eff (Hesslow et al. 2018b). Here γ = 1 + p 2 is the Lorentz factor, Z eff = n −1 e j n j Z 2 j is the effective charge and j is an index which runs over all ion species in the plasma, each with density n j and charge Z j .
In partially ionised plasmas, the slowing-down and deflection frequencies are influenced by the extent to which fast electrons can penetrate the bound electron cloud around the impurity ion, i.e. the effect of partial screening (Martín-Solís et al. 2015;Breizman & Aleynikov 2017). The collision frequencies ν s and ν D can be generalised to account for the differences in the collisional dynamics at different energy scales arising when screening effects are introduced, and in the relativistic limit take the form (Hesslow et al. 2017(Hesslow et al. , 2018a which we will denote asν D ≈ν D0 +ν D1 ln p andν s =ν s0 +ν s1 ln p. The collision frequencies now depend on atomic parameters of species j: ionisation degree Z 0j , charge number Z j , number of bound electrons of the nucleus for species j, N j = Z j − Z 0j , mean excitation energy of the ion I j and effective ion sizeā j determined from density functional theory calculations, given by Hesslow et al. (2017). The effect of partially ionised ions in the plasma will influence the runaway generation (Hesslow et al. 2019a,b) as well as increase the critical electric field (Hesslow et al. 2018b). Synchrotron radiation and bremsstrahlung hinder the acceleration of runaway electrons. The effective term in the kinetic equation resulting from synchrotron radiation is (Stahl et al. 2015;Hirvijoki et al. 2015a where τ syn is the synchrotron radiation-damping time scale τ syn = 6πε 0 m 3 e c 3 / e 4 B 2 . (2.6) Similarly to the treatment by Hesslow et al. (2017), the effect of bremsstrahlung is here incorporated into the kinetic equation via a mean-force model, which has been shown to capture the mean-energy accurately (Embréus et al. 2016), Here F br is approximated by and α FS is the fine structure constant. The screening and radiation effectively increase the friction at large momenta, which will prevent runaway electrons from reaching arbitrarily large energies when given a long enough time to accelerate. Equation (2.1) describes the momentum space dynamics of the runaway phenomena, however it does not include any terms allowing for spatial transport. Helander et al. (2000) amended the kinetic description by adding the radial component of the diffusion operator in cylindrical geometry, characterised by the phase-space dependent diffusion coefficient D. Including this in our formulation we obtain the full kinetic equation of interest here, (2.9) In the next section we formulate a general expression for the change in the runaway avalanche growth rate resulting from such a finite D.

Reduced avalanche growth rate and effective critical electric field
The appearance of the diffusion term in the kinetic equation adds another dimension to the problem, a radial one, compared to the standard avalanche growth rate calculation (Jayakumar et al. 1993;Rosenbluth & Putvinski 1997). We develop an approximate solution by taking advantage of a separation of time scales, following the approach outlined by Helander et al. (2000). The avalanche generates secondary electrons with momentum predominantly close to the critical momentum for the runaway process p c † and the electron is accelerated from this region up to relativistic momenta on the short timescale τ acc ∼ τ /E. The transport timescale represented by D will typically be longer than this acceleration time, so the diffusion will not be strong enough to significantly reduce or alter the generation process. However, the timescale of the avalanche growth is significantly longer, namely γ −1 r ∼ 2 ln Λ τ acc (Jayakumar et al. 1993) and thus the spatial diffusion may be expected to have a substantial impact on the avalanche.
To this end, the momentum space is divided into a low energy region, p < p * , where all the runaway generation occurs and the effects of radial diffusion are neglected, and a high energy region, p > p * , where all the radial transport takes place. After this division in momentum space, the high energy region is modelled as source free and the generation of runaway electrons is modelled as a flux through the lower boundary in momentum space at p * . Furthermore, the theory is reduced to only a single momentum-space coordinate in the high energy region.
This was done neglecting the effect of radiation in (Helander et al. 2000), by recognising that runaway electrons often have small pitch-angles, ξ ≈ 1 and so expanding the collision operator assuming p ⊥ p , where p and p ⊥ are the projections of the momentum along and perpendicular to the magnetic field line, respectively. The rate of change of the distribution function integrated over perpendicular momenta F = d 2 p ⊥ f = ∞ 0 dp ⊥ 2πp ⊥ f can then be obtained by integrating the kinetic equation (2.9), leading to ∂F ∂t where the diffusion coefficient for particles travelling purely along the magnetic field line is used to first order. This is the kinetic description of the runaway electrons given in equation (12) of (Helander et al. 2000). The synchrotron radiation reaction force is however strongly dependent on the momentum perpendicular to the magnetic field line, † The avalanche source term is 1 p 2 ∂ ∂p 1 γ−1 ∼ 1 p 5 for low momenta, and therefore does not extend far into the runaway region.
as it is a consequence of the gyration around the field line, and a treatment of radiative effects needs to account for the pitch angle distribution of particles.
Radiative effects are important close to the critical electric field where the acceleration from the electric field is close to being balanced by the radiation reaction forces, making the dynamics in the energy direction of momentum space comparatively slow. Therefore, as in Lehtinen et al. (1999); Aleynikov & Breizman (2015); Hesslow et al. (2018b), we consider the pitch-angle evolution to be a rapid process compared to the energy evolution dynamics and assume a steady-state in the pitch-angle distribution for a given p. This requires that the pitch-angle flux of particles vanishes, the condition following from the kinetic equation (2.9) as Since τ syn τ we formally neglect the effect of the synchrotron radiation on the pitchangle distribution, retaining only the balance between the diffusive effect of pitch-angle scattering and the collimating effect of the electric field. This can be used to solve for the pitch angle distribution and the distribution function may now be written as where the reduced distribution function F includes the 2πp 2 of the momentum-space Jacobian, such that the radial density of runaway electrons is and the inverse of A(p) = 2E/ (pν D τ ) determines the extent of the distribution function in ξ. In the limit of large p, the pitch angle distribution is narrow in agreement with the treatment by Helander et al. (2000), as ν D ∼ p −2 and therefore A −1 ∼ p −1 . Integrating the kinetic equation over pitch-angle, a reduced kinetic equation now accounting for radiation reaction forces and screening effects is obtained where the pitch averaged force U (p) is and the pitch-angle averaged diffusion coefficient is Note that for large p, neglecting screening effects, we recover the non-radiative result U (p) ≈ E coth A − pτ ν s → E − 1. A qualitative difference between the models with and without radiative effects is the appearance of a momentum scale p max where the pitchangle averaged advection in momentum disappears, U (p max ) = 0. This limits the energy of the relativistic particles and corresponds to the energy scale where radiation reaction forces balance the electric field acceleration. As outlined by Helander et al. (2000), we impose the boundary condition that the particle flux through p * is given by the avalanche growth, γ r n RE , where γ r is the growth rate without the impact of diffusion. Given the structure of equation (3.5) this translates to a condition on F as which is not a typical boundary condition as the value at the lower boundary in momentum space is dependent on the solution in the whole high energy region. A closedform expression for the solution can be found for the simpler problem of radially-uniform plasma parameters in a quasi-steady state in terms of a Bessel mode expansion †. The effective strength of the diffusion experienced by each Bessel mode is scaled by the square of its inverse radial length scale (Abramowitz & Stegun 1948). The solution is, The coefficients c i are determined by the initial condition, or seed profile of the avalanche process. The growth rate of the modes, γ i , will be determined by the boundary condition (3.8). Helander et al. (2000) considered only the first Bessel mode, i = 1. This gave a conservative estimate of the effect of diffusion as higher mode numbers have a smaller characteristic length scale so will experience a larger effect of diffusion, as noted above.
Here we choose to retain all the modes, which would allow the runaway distribution function to be propagated in time. As a consequence of the orthogonality of the Bessel modes, the boundary condition can be projected on each mode separately, which decouples them from one another. The equation for γ i then follows from inserting (3.9) in the boundary condition (3.8) as , (3.10) where the upper limit of the integration is p max as no particles can gain energy larger than this. The theory for the pitch-angle distribution described above is valid for large p -which we are mostly concerned with here -and not in general close to the critical momentum p c , where U (p c ) = 0. If the free parameter p * is chosen close to p c , the result will be sensitive to the choice. We can consistently minimise this impact of p * by expanding U in large p, such that the theory still retains a p max , and safely set p * = p c as a typical momentum scale for the onset of the runaway region. A large-p expansion keeping terms of order p −1 and larger gives (3.11) Finally, a model for the avalanche growth rate without magnetic perturbations is † The same type of solution is still possible with a radially dependent diffusion coefficient, as the radial part of the problem forms a Sturm-Liouville problem, however the expansion would no longer be in terms of Bessel functions but rather the eigenfunctions of the transport term in question.
needed. For continuity with our collision operator we will use the model by Hesslow et al. (2019a), which incorporates the effect of screening by the evaluation of the collision frequencies at an effective critical momentum scale p eff c , implicitly given by Furthermore the model has an increased threshold field for the avalanche generation,Ē eff c , given by Hesslow et al. (2018b). Here, the bar onĒ eff c indicates the critical electric field without the effect of radial diffusion, to distinguish from its value E eff c when the radial diffusion is taken into account in section 3.2. The expression for the avalanche generation is then where n tot e is the total number of electrons in the system (bound + free). Although we will use this model for the avalanche generation in the next section in order to consider its reduction due to radial transport, our method is agnostic to this choice and can be adapted to any avalanche description.

Reduced avalanche growth rate
As the radial transport allows for runaway electrons to be lost from the system, preventing these electrons from multiplying by the avalanche mechanism, the exponential growth rate of the Bessel modes, γ i , will be reduced compared to the uncorrected value. The magnitude of the transport coefficients considered below reduces at large momentum, and in the avalanche distribution the particle density in phase space also decreases at large energies. Therefore, the problem is to a large extent determined by the dynamics at small energies. At these energies, where p p max , the radiation reaction forces do not influence the problem significantly, and the acceleration dynamics is dominated by the electric field. For large field strength we have U ≈ E. However, at the critical electric field for runaway generation,Ē eff c , the runaway generation and the advection in momentum space fall to zero. A linear interpolation between these regions gives U ≈ E −Ē eff c . The uncorrected growth rate in equation (3.12) scales with electric field strength as E −Ē eff c , which was just noted to be the approximate dependence of U at low momentum (p p max ), such that the prefactor of the exponent in (3.10) does not significantly vary with electric field strength at low p. The impact of the diffusion at low momentum is then characterised by and the dependence of the corrected growth rate on this parameter is shown in figure 1a, which is obtained by solving equation (3.10) numerically. In the limit of small α the effect of diffusion is rather well parametrised by this normalised ratio of diffusion strength to the electric field acceleration. However, the correlation is lost when the growth rate is strongly reduced and approaches zero. Qualitatively, this can be understood as the effect of p max and the high energy particles, as a finite p max limits the number of energetic particles that can contribute to the avalanche without being significantly affected by the diffusion. The effect is demonstrated in figure 1a by the case U = E −Ē eff c where p max is formally infinite and the growth rate remains slightly above zero for a relatively wide range of diffusion strengths †. This limit is approached generally at large electric field strength. † Formally, the growth rate cannot be corrected down to zero in the theory where U = E−Ē eff c if the D decays asymptotically in momentum space faster than p −1 . The latter is the marginal case where it is impossible if α0 < 1 where α ∼ α0/p asymptotically.
In figure 1a a diffusion coefficient of the following form was used, (3.14) The motivation for this expression is to capture the expected low energy behaviour where the diffusion is proportional to the particle velocity along the magnetic field line (Rechester & Rosenbluth 1978), as well as the expected effects of orbit decorrelation at high energy, as described for example in Hauff & Jenko (2009). The latter authors made an estimate for D 0 in the small Kubo number limit, such that D v 2 B τ , where the radial velocity of the particles v B = v δB/B is due to the projection of the motion along the perturbed field line, with δB the root mean square of the magnetic perturbation amplitude. The parallel correlation time is assumed to be set by the particle motion through the perturbed poloidal magnetic field structure, τ = λ /v πqR/v , where λ is the parallel connection length and q is the safety factor. Therefore, the estimate of D 0 is The strength of the diffusion is parametrised in this paper by τ k 2 i D 0 which is thus related to the magnetic perturbation level as follows, and a[m] are the major and minor radii in meters, respectively, and n e,20 is the electron density in units of 10 20 m −3 . Consequently, for standard ITER parameters without any material injection, R = 6.2 m, a = 2 m, n e = 10 20 m −3 , ln Λ ≈ 15 and q ≈ 1, we have τ k 2 1 D 0 ≈ 2(δB/B) 2 × 10 8 for the least suppressed mode (b 1 ≈ 2.4). In absolute units, the uncorrected growth rate scales linearly with the electric field strength and in figure 1b we see that the corrected growth rate also shows a linear relation with the electric field strength for large fields. The corrected growth rate is offset from the uncorrected one because the inverse of the characteristic diffusion parameter α and the uncorrected growth rate depend similarly on the electric field. This can be seen by expanding equation (3.10) in small α yielding † for the values of p with the largest contribution to the integral. This approximation is compared to the full numerical solution in figure 1b where the offset from the uncorrected result is evident for large E and small α.

Effective critical electric field
As the radial transport reduces the growth rate, the critical electric field strength for net generation of runaway electrons may increase. In the current decay phase of the disruption the electric field stays close to the critical electric field and the current decay rate is proportional to its value (Breizman & Aleynikov 2017;Hesslow et al. 2018b). Therefore the critical electric field has direct relevance for disruption mitigation strategies. † This formula can be arrived at by using integration by parts and change the momentum variable q = p * + p p * dp E −Ē eff c /U (p ) in the intermediate steps. The transformation maps the problem to a theory without radiative effects.
The corrected growth rate as a function of electric field strength. At large electric field strength the offset of the corrected growth rate depends on the diffusion strength τ k 2 i D 0 , which is expected to be around unity in an ITER-sized machine with normalized magnetic perturbation level of δB/B 10 −4 .
The mode least suppressed by the transport is the lowest order mode, with index i = 1, and therefore can be expected to dominate the runaway profile in the late stages of the disruption. We therefore choose to define the effective critical electric field E eff c as the field strength at which the growth rate of the first mode vanishes, namely γ 1 = 0. Figure 2a shows numerical solutions for the critical electric field, obtained from equation (3.10) under the constraint γ 1 = 0 in fully ionised plasmas with different densities. For large diffusion strengths, when the effective critical electric fields are relatively large, a linear relation between diffusion strength and electric field is found. This follows naturally for large electric fields E in (3.10), U ≈ E −Ē eff c and p max is at large enough energy scales not to be relevant, as the condition γ 1 = 0 translates to a condition on α. Given this linear relation in diffusion coefficient, the critical electric field strength is expected to be quadratic in the magnetic perturbation level, δB/B. Figure 2b shows the critical effective field as a function of temperature, which introduces screening effects at the lowest temperatures, where some electrons remain bound to ions. The ionisation states are determined here by assuming equilibrium based on the ADAS coefficients of ionisation and recombination †.
In absolute units, the correction to the effective critical field strength increases with temperature, primarily due to the increase in free electrons which raises E c . Massive material injection will also raise the density and so the critical electric field for runaway generation. We see from figure 2b that this effect can be combined with the effect of spatial diffusion to further raise E eff c . However unlike massive material injection, where changes to E eff c are linked to changes in the electron density, the correction based on magnetic perturbations is only weakly dependent on the density (as long as perturbation strength is treated independent of plasma density), in absolute units. This weak dependence follows as τ γ r /E is density independent and therefore any sensitivity originates from onlyĒ eff c † http://www.adas.ac.uk and the large-p dependence of U . Further, as was shown by Hesslow et al. (2018b), the effect of partial screening on the critical electric fieldĒ eff c was to raise it to the order of E tot c , which has the same form as the usual Connor-Hastie expression E c , but with the combined density of free and bound electrons instead of only the density of free electrons. Therefore, only a weak dependence of temperature is seen in figure 2b, as the total number of electrons are kept fixed in these simulations. The two mechanisms for increasing the effective electric field (screening and magnetic perturbations) can therefore be combined. We note, however, thatĒ eff c is only weakly dependent on the temperature and ionisation state.
The theory discussed so far, in which radial variations in the plasma have been neglected, allows for self-consistent analytic solutions to the distribution function and an understanding of the dependencies of the growth rate correction due to transport. However, in a tokamak disruption, the electric field dynamics is essential and will vary spatially through its dependence on plasma properties. Therefore in the next section we take a perturbative approach to solving equation (3.5), to include effects due to radial plasma variation.

Transport in an inhomogeneous plasma
During the current quench of a disruption the flux surfaces are not completely stochastic. Instead, they often exhibit a mixed magnetic topology consisting of intact flux-surfaces, magnetic islands and stochastic regions. In such circumstances, the transport will no longer be well described by the expression given by Rechester & Rosenbluth (1978) and a more general transport model consisting of spatially dependent diffusive and advective components is often formulated, based on particle following simulations Särkimäki et al. 2016). Assuming cylindrical symmetry, the transport term on the right hand side of equation (3.5) becomes where V ξ is the pitch-angle average of the radial component of the advection coefficient defined equivalently to (3.7). Fundamentally the above transport term conserves particle number, which is a property not guaranteed by an approximate perturbative solution. Therefore, this conservation property will be imposed on the solution to prevent anomalous losses of particles. Note, that the conservation of particle number is local, and particles can be lost at the edge. The approach to solving the kinetic equation in the high energy region with a momentum space dependent diffusion coefficient in the previous section, by means of a Bessel mode expansion, breaks down when a radial dependence is introduced in the plasma parameters. To include the effects of radially varying plasma parameters we instead perform an expansion in small radial transport compared to the electric field acceleration α 1. The full form of U gives U = 0 in the vicinity of p max , so the transport would not be subdominant for such momenta. We therefore treat this with a simplified approach, taking U = E −Ē eff c which is similar to the advection in momentum space given by Helander et al. (2000), modified to include the effects of an increased critical electric field due to screening. In this model, the advection in momentum-space does not vanish, and in the previous section, we saw that the full U approaches this in the limit of large electric field. As this model does not treat the effects of radiation correctly, it may be seen as a better approximation significantly below p max . As the assumed transport coefficients and particle density decrease with momentum in avalanche dominated scenarios, these lower energy scales are anyway expected to dominate the transport here. We will also demonstrate explicitly the weak sensitivity of the results to p * .
Given this model for U , the zeroth order solution given by neglecting the transport, so the solutions at different radii are independent, is which is consistent with the incoming runaway electron flux γ r n RE and respects the definition of n RE in (3.4). Under the assumption that the radial transport is small this distribution may be used to evaluate the transport term. Integrating equation (3.5) over p then gives where the term on the right represents the source of runaway electrons and is the radial flux of the runaway electrons. Using equation (4.2) for F 0 to evaluate Γ 0 results in a form Γ 0 =Γ 0 n RE +Γ 0 ∂ r n RE and equation (4.3) can be stably solved numerically using a method based on the Crank-Nicholson scheme. This transport model can be incorporated into any suitable runaway simulation to similarly capture the effects of transport due to magnetic perturbations. In the following subsections the transport model will be integrated with the electric field evolution, but it should be noted that the momentum space shape of the distribution function, equation (4.2), is insensitive to the electric field strength, due to the appearance of the ratio γ r /(E −Ē eff c ), which is consistent with the neglect of the temporal evolution of E in its derivation.

Simplified disruption simulations including the effect of radial transport
Under normal operation of a tokamak there are radial gradients in the temperature and plasma current, both of which contribute to a radially varying electric field as the plasma is suddenly cooled in a disruption. The subsequent evolution of the electric field is described by the induction equation, which for a plasma with only radial variations is (Smith et al. 2006) where E and j are the electric field strength and the current density along the plasma cylinder.
The time evolution of the electric field is dependent on its radial profile, as the current j has both an Ohmic and a runaway component: j = σ sp E + j RE ≈ σ sp E + ecn RE where σ sp is the Spitzer conductivity (Spitzer & Härm 1953). Accordingly, the radial profile of runaway generation also plays a crucial part in understanding the electric field evolution. Furthermore, the avalanche growth rate is proportional to the electric field strength, for fields large compared to the critical one, such that the cumulative generation is highly dependent on the evolution of the electric field. Consequently, for a self-consistent treatment of both the electric field dynamics and runaway generation, it is of the utmost importance to be able to treat the runaway generation in a region of space with an electric field gradient. Such computations can be carried out within the go-framework (Smith et al. 2006;Fehér et al. 2011;Vallhagen et al. 2020) or similar codes. Introducing the effect of transport losses into such a model, including the effects of impurities and allowing for partial screening, would allow us to quantify the reduction given by transport of the total number of runaway electrons at the end of a disruption.
To achieve this, we have extended the go-framework to solve the coupled equations (4.3) and (4.5). The runaway generation and transport are described by equation (4.3), while the runaway electron dynamics couples to the electric field evolution through the current term in equation (4.5). The right hand side of the former includes primary sources of runaway electrons: Dreicer generation, tritium decay and Compton sources. When avalanching dominates, as is the case for high current devices, the flux Γ 0 derived using only the avalanche source should be valid as it describes the momentum space distribution of the majority of the population. Finite-aspect-ratio effects on the generation will be neglected here, as recent work by McDevitt & Tang (2019) indicate that their effect is negligible at the high densities and electric fields that we will consider here.
In the go-framework, the Dreicer generation is evaluated using a neural network (Hesslow et al. 2019b), trained on kinetic simulations with code (Landreman et al. 2014), using the collision operator that includes the effect of partially ionized impurities given by Hesslow et al. (2018a). β-decay of tritium will also result in a source of runaway electrons, which in the deuterium-tritium phase of operation is expected to be the dominant source of seed electrons in ITER in the absence of hot-tail electrons (Martín-Solís et al. 2017). Furthermore, neutrons produced in the fusion reactions will activate the wall which in turn will emit γ-photons. Through Compton scattering events between the γ-photons and the bulk electrons, runaway electrons can be generated (Martín-Solís et al. 2017; Vallhagen et al. 2020). In the simulations presented here, we neglect the hot-tail generation occurring in a rapidly cooling plasma. This generation occurs during the thermal quench which is typically associated with large magnetic fluctuations and corresponding transport. By neglecting the hot-tail seed we implicitly assume that the transport during the thermal quench is large enough to lead to the prompt loss of these runaways.
The go-framework has the capability to compute the plasma temperature evolution from the energy balance between heat diffusion, Ohmic heating, line radiation, bremsstrahlung losses and ionisation, as presented by Vallhagen et al. (2020). However, in the initial phase of the disruption, the energy loss is expected to be dominated by the MHD contribution, due to its strong temperature scaling ∼ T 5/2 (Ward & Wesson 1992). This phase is, for simplicity, modelled as an exponential drop in temperature until the temperature of the inner part of the plasma drops to ∼ 100 eV, with the form where t 0 is the time constant for the thermal quench and T i and T f are the initial and final temperatures, respectively. This mode of the temperature evolution uses a flat final temperature profile T f = 50 eV and is used for 6 ms with a time constant of t 0 = 1 ms. After this time the temperature is evolved based on the energy balance. The ionisation states in the background plasma are evolved in time based on the ADAS coefficients for ionisation and recombination.

Simulations of ITER-like disruptions with uniform perturbations
To investigate the large scale effect of radial transport, in tokamak disruption scenarios where the runaway generation is expected to be dominated by the avalanche mechanism, an ITER-like case with deuterium and neon injection was simulated using the goframework. The parameters considered are the same as those used by Martín-Solís et al. (2017) and Vallhagen et al. (2020): initial plasma current I p (t = 0) = 15 MA, minor and major radii a = 2 m and R = 6.2 m, respectively, initial electron, deuterium and tritium densities n e0 = 2n D0 = 2n T 0 = 10 20 m −3 . The simulation was initiated with one dimensional profiles in temperature T e = 20 1 − (r/a) 2 keV and current density j (t = 0) = j 0 1 − (r/a) 2 0.41 , where j 0 is chosen so that the current integrates to 15 MA.
At the start of the simulations we assume a rapid injection of deuterium and neon, with respective densities n D and n Ne , where the impurity is distributed uniformly throughout the plasma in the neutral state. The current evolution for three such scenarios with different combinations of injected neon and deuterium is demonstrated in figure 3, for a set of radially constant magnetic perturbation levels and a momentum space dependent diffusion coefficient of the form (3.14). The three cases considered are the same as those investigated in Vallhagen et al. (2020) and denoted both here and there as Case 1, Case 3 and Case 4 †. In each of these cases the injected material is large enough to induce a complete thermal quench: Case 1 (n Ne = 1 × 10 20 m −3 , n D = 0), Case 3 (n Ne = 8 × 10 18 m −3 , n D = 4 × 10 21 m −3 ) and Case 4 (n Ne = 8 × 10 18 m −3 , n D = 7×10 20 m −3 ). In the absence of perturbations large runaway currents were obtained in all of these three cases, even without hot-tail generation.
Interestingly, the maximum runaway current increases for small transport coefficients (δB/B 2 · 10 −4 ) compared to the baseline case of no radial transport, as the runaway electron seed is radially flattened by the transport. This agrees with previous results by Three cases of injected material are considered: A pure neon injection with density n Ne /n e0 = 1 (Case 1), and two cases with the same amount of injected neon n Ne /n e0 = 0.08, but different amount of injected deuterium n D /n e0 = 40 (Case 3) and n D /n e0 = 7 (Case 4), respectively. The coloured area corresponds to p * in the range 0.1 − 1. Fehér et al. (2011). However, for large enough perturbations we note a reduction of the maximal current carried by the runaway electrons. How large the reduction is depends on the particular scenario. Case 4, with a combination of moderate neon and deuterium injection shows the largest reduction.
The effectiveness of the radial transport in modifying the runaway evolution is closely related to the time scale of the current evolution. Generally, the longer the time scale of the current quench the more pronounced is the effect of transport, as particles have more time to be transported out of the plasma. Due to this effect Case 4 has a slower growth rate of runaway electrons than in Cases 1 and 3, diffusion therefore having a larger impact, as can be noted in figure 3. Figure 4 shows the maximum runaway current, just before the onset of the dissipation phase where the plasma current carried by the runaway electrons decays, in the three cases as a function of (δB/B) 2 . We note that almost full suppression of the runaway current can be achieved in Case 4, for a normalised perturbation δB/B 5 · 10 −4 . We also show the time it takes for the runaway current to rise from 10% to 90% of its maximum value, denoted by t 10-90 , for δB/B = 2 · 10 −4 . Clearly, Case 4 has considerably longer t 10-90 than the other two and this is the main reason for the larger suppression. The diffusion time scale can be estimated to be t diff = a 2 / D 0 a 2 /(π q Rc) (δB/B) −2 , and is 17 ms for q 1 and δB/B = 2 · 10 −4 . Here, · · · denotes a volume average. We note that the uncertainty in p * influences the result, but the effect of the magnetic perturbation is clearly dominant.
When the runaway plateau phase is reached in the final stages of the disruption the loss of plasma current is dominated by the loss of runaway electrons. If the time derivative of the Ohmic current is neglected in equation (4.5), combining with equation (4.3) yields an approximate equation for the electric field in the runaway plateau,   Based on this expression the decay of the plasma current is for a given radial profile of the runaway electron density. This approximation recovers the electric field structure in the current decay phase to a large extent, however an even simpler consideration can be made where the local growth of runaway electrons is set to zero, in equation (4.3), which in terms of equation (4.7) corresponds to neglecting the left hand side of the equation. Using the expression for the uncorrected growth rate under consideration here, equation (3.12), yields an explicit expression for the electric field given a profile n RE . Figure 5 shows the electric field obtained in the simulations corresponding to the cases shown in figure 3, together with the electric field strength which zeros the local growth of runaway electrons. The electric field profiles are seen to agree with one another in regions with small gradients in E and especially in the central part of the plasma, which is related to the neglected term from equation (4.7). Therefore, the electric field is generally higher than the threshold field for runaway generationĒ eff c in the centre of the plasma, where there is a balance between the generation and transport. However, close to the edge the current density is low enough so that the prefactor on the source term in equation (4.7) is small enough to allow a significant deviation belowĒ eff c . Furthermore, in the cases with large transport, the electric field in the central region is not as flat asĒ eff c . Instead it has a similar functional form to the runaway current profile, indicating that it is the transport which dominates the electric field in the plateau phase. In these situations, the electric field is highly dependent on the profile of runaway electrons, which in turn depends on the full temporal evolution of the system, and in particular the transport coefficients. This suggests that approaches where the current decay phase is described byĒ eff c are only valid if the transport is negligible, otherwise the coupled dynamics with the runaway electrons must be considered.

Runaway dynamics in the presence of artificial resonant perturbations
The magnetic field is expected to become fully stochastic at the end of the thermal quench, after which it begins to heal during the current quench. To mimic the conditions during the current quench in ITER, we choose the 15 MA / 5.3 T baseline scenario (Parail et al. 2013) and introduce artificial resonant magnetic perturbations at the plasma edge to create a stochastic layer. We have chosen the pre-disruption current flat top equilibrium for this exercise as obtaining realistic current quench equilibrium would require dedicated MHD modelling. The introduced perturbations are stationary and have a helical structure, where (ρ, θ, ζ) are the radial, poloidal, and toroidal Boozer coordinates, respectively, B is the unperturbed field and the phase φ nm is chosen to be random. This method is the same as used in Särkimäki et al. (2020), and also here the total perturbation consists of several modes with low mode numbers (n, m 20). The mode eigenfunctions are Gaussians, that peak at the corresponding resonance r nm and all have the same width σ = 0.03 m which is large enough for the modes to overlap and create a continuous stochastic region. The perturbation level is set to δB/B ≈ 10 −3 at which significant runaway transport is expected (Helander et al. 2000). The resulting field is illustrated in figure 6. The transport coefficients are evaluated numerically with the orbit-following code ASCOT5 (Varje et al. 2019). Markers representing guiding centers of collisionless electrons were traced in the perturbed field, and their radial position was recorded for each orbit. As particles with finite orbit width oscillate radially during their orbit, the radial position was always recorded at the same poloidal position (at the outer mid-plane) so that all changes in the radial position were due to the transport alone. No collisions, electric field Figure 6: Magnetic field Poincaré-plot at the outer mid-plane for an ITER current flattop equilibrium perturbed with artificial resonant magnetic perturbations according to (5.1). The stochastic region begins at r/a ≈ 0.6 (q = 1 surface is at r/a ≈ 0.5).
or radiation reaction force was present in the simulation to isolate the transport due to the magnetic field perturbations and to keep momentum constant in order to calculate momentum dependent coefficients. The transport coefficients are evaluated from the recorded radial positions as where the brackets denote an average over all collected data points for a marker i, ∆t is the orbit circulation time, ∆r is the change in radial position between subsequent recordings, and the sum is taken over all N markers that were traced. This scheme is similar to the one originally presented by Boozer & Kuo-Petravic (1981). At the edge markers are lost within a few orbits, making estimates (5.3) and (5.4) unreliable, and so the coefficients are evaluated from the loss-time distribution using the method described in (Särkimäki et al. 2016(Särkimäki et al. , 2020. This latter method is used if more than half of the markers are lost. Markers are simulated for 2 × 10 −5 s which corresponds to approximately 100 orbit transit times. The simulation time has to be long enough as early on the particle orbits are correlated and the motion is not diffusive. However, longer simulation time decreases the radial resolution of the evaluated coefficients as ∆r ≈ √ 2Dt. The markers have identical radial position, pitch and momentum but a random toroidal location. The transport coefficients to be used in the go simulation are found by repeating the orbit-following simulation with different values of radius, pitch, and momentum. In these simulations, the phase space is divided into 15 radial, 11 momentum, and 10 pitch slots (covering the passing particle regime). For each volume element 200 markers are simulated to calculate the coefficients at that point. The resulting advection and diffusion coefficients are shown in figure 7 for a fixed pitch, as the coefficients show no strong pitch dependence as long as the particles are passing. Radially the transport is almost uniform in the region where the field is stochastic (recall figure 6) while the momentum dependence shows decreasing transport for higher energies due to the finite orbit width effects (Hauff & Jenko 2009;Särkimäki et al. 2020).
In the inner region of the plasma (r/a < 0.58), the runaway electrons are not transported as the flux surfaces are intact, and markers initiated in the stochastic region will not be transported into this region. To properly capture this effect in a simulation with the go-framework, a reflective boundary condition was imposed at the first complete flux surface, and no transport could occur between the regions. However, the regions are still coupled through the electric field evolution.
In scenarios where the runaway generation is dominant in the central part (such as in Case 1 and Case 4) the stochastic plasma edge is not expected to affect the runaway dynamics significantly. By using the advection and diffusion coefficients shown in figure 7, and simulating an ITER-like scenario with material injection corresponding to Case 4, we find that the maximum runaway current is reduced only marginally, from 3.7 MA to 3.5 MA. To illustrate a case when the effect of a stochastic plasma edge is more pronounced, we consider Case 3, where in the absence of radial losses an off-axis final runaway profile is found. Therefore the transport should have a larger impact in this case compared to scenarios with an on-axis final current profile, where a larger part of the runaway electrons are generated in the non-transporting region. Figure 8 shows the radial profiles of the runaway current after 45 ms in the ITERlike disruption simulation of Case 3. The maximum runaway current, just before the dissipation phase, in the absence of radial transport due to magnetic perturbations is 7 MA, with a constant δB/B is 5.8 MA and with the coefficients presented in figure 7 is 4.6 MA. Without magnetic perturbations, the profile of runaway electron density has an off-axis maximum. This is due to strong radiative losses, leading to significant plasma cooling, and corresponding efficient runaway generation in the outer part of the plasma, as was pointed out by Vallhagen et al. (2020). Note, that such off-axis current profiles may become MHD unstable. The resulting magnetic activity may act to mitigate the build up of runaways.
In the presence of magnetic perturbations electrons can diffuse and this results in a final runaway current profile that is peaked on-axis. However, in the case with the stochastic edge, with transport coefficients shown in figure 7, the transport in the edge region is strong enough to prevent any significant build up of runaway electrons there. The increase of the transport at the edge results in only partial reduction of the total runaway current, as the runaway population is free to build up in the centre of the plasma. In the confined inner region strong gradients in the current density can form, which is evident in the simulation. At the transition from the confined region to the stochastic one (at r/a = 0.58) a discontinuity is formed in the current profile (however not in the electric field), as particles in the stochastic outer region are continuously transported away, but in the confined region they are free to build up, eventually forming a current sheet.
The radial profiles of the temperature, electric field and number of e-foldings (the timeintegral of the runaway growth rate) are shown in figure 9 for Case 3, at a few time slices. The vertical dashed line denotes the radial position for the transition between confined and stochastic regions. Figure 9a shows that, both with and without perturbations, the plasma is divided into two regions by a cold front, with an inner region with a temperature of about 6 eV, and an outer region with a temperature as low as about 1 eV. At such low temperatures a large fraction of the deuterium recombines, and this leads to an increased avalanche multiplication of the seed runaway electrons in the outer region, see figure 9c, quantified by the number of e-foldings †, N exp (r) = t 0 dt γ(t , r).
(5.5) † exp (Nexp) is the factor by which the avalanche mechanics amplifies the seed in a non-transporting model.  The avalanche production continues throughout the simulation, but is counteracted by the strong radial transport in this region.
The formation and strength of the current sheet seen in figure 8 is a result of the interaction between runaway transport and strong diffusion of the electric field. The location is tightly connected to that of the cold front, which propagates in from the plasma edge in the later stages of the simulation. In the scenario without perturbations, the cold front propagates inwards faster. Figure 9 compares the evolution of the electric field in the two scenarios after the cold front has crossed out of the stochastic region, at times when the cold front has reached the same position, to highlight the dynamics behind the current sheet formation. In the case with transport there is less conversion from Ohmic to runaway current in the outer parts of the plasma, so a significantly larger electric field is maintained in the outer region, as is seen in figure 9b. When the conversion starts in the inner regions, a sharp change in the gradient of the electric field develops at the interface to the stochastic region, which enhances the diffusion of the strong electric field in the inner region. This results in a larger avalanche multiplication of the seed runaway electrons in the boundary between the regions, which is demonstrated in figure 9c. Despite the strong amplification in the outer region, the transport prevents a significant number of runaway electrons building up, and so the current sheet is formed.
Another way of thinking of this phenomenon is to consider the effect of the runaway electrons on the electric field evolution. As the transport in the outer region is strong enough to prevent a runaway build up, the runaways do not affect the electric field evolution -this is akin to the assumption of a trace electron population used in the estimate by Rosenbluth & Putvinski (1997) -which leads to a very large amplification factor. However, the coupled dynamics must be considered in the presence of a large runaway population, and the amplification saturates as the current carried by the runaways in the inner region here approaches the Ohmic current. Then the current sheet is formed at the interface between the two regions. Such a current profile is likely to be very unstable, and this could affect the magnetic equilibrium and lead to magnetic perturbations penetrating deeper into the plasma core, giving further runaway mitigation. Such a study is beyond the scope of the current paper.

Conclusions
During tokamak disruptions the magnetic field lines can be severely distorted from their usual confining structure. The magnetic topology evolves in time, being almost fully stochastic during the thermal quench, and often displaying a mixed topology of intact flux-surfaces, magnetic islands and stochastic regions during the current quench. As runaway electrons travel rapidly along the tokamak magnetic field lines, their evolution during a disruption can be strongly affected by such magnetic perturbations. This introduces the possibility for radial losses of runaway electrons to offset the avalanche growth of the population, preventing the formation of a high current, potentially damaging, runaway electron beam.
In this paper, we have presented a model which generalises previous treatments of the effects of radial transport due to the interaction of runaway electrons with magnetic field perturbations on the runaway evolution. We continue to take advantage of the separation of timescales in the runaway generation dynamics, between the acceleration to relativistic energies after a knock-on collision and the characteristic avalanche population growth time. This allows us to neglect the effect of transport due to magnetic perturbations on the generation process and focus on solving the kinetic equation in the high energy limit, which simplifies the collision operator. The extension here allows for a generalised pitch-angle distribution formed by rapid pitch-angle scattering at high energy, the impact of radiation reaction and the presence of partially ionised impurity atoms. The effect of including radiation is to introduce an upper limit in momentum space, so particles are prevented from reaching very high energies where they would be well confined.
In particular, we have determined an expression which can be used to correct the growth rate of the runaway electron population by the avalanche mechanism. This takes the form of the solution of an integral equation. The introduction of radial transport raises the effective critical electric field for avalanche generation because even though particles can be kicked into the runaway region through the avalanche mechanism, they can be lost due to spatial diffusion resulting in no net gain of runaways. The increase in E eff c is weakly dependent on the plasma density, unlike the case of massive material injection.
The derivation of the integral equation for the growth rate corrections is only valid when radial variations in the plasma are neglected. To treat non-homogeneous plasmas, a perturbation approach in small radial transport has been used to estimate the radial flux of runaway electrons. By computing the radial fluxes instead of an effective (local) growth rate the formulation is particle conserving. This flux can be included in general runaway simulation frameworks. The formulation was used here in a simplified disruption simulation for ITER-like plasmas, where an induction equation was used to give the selfconsistent time evolution of the electric field in the presence of the runaways.
We find that in scenarios with a moderate amount of impurity and deuterium injection, the runaway current can be suppressed for perturbation levels of the order of δB/B ∼ 5 · 10 −4 , which is about an order of magnitude higher than the perturbation level measured in fixed magnetic field experiments where the current was scanned (Gill et al. 2002). Furthermore, it is difficult to fully dissipate the runaway electrons without significant transport in the centre where most of the runaway electrons are generated. Earlier results investigating the potential of employing edge-localized mode (ELM) coils for runaway suppression show that perturbations created at the plasma edge are generally not sufficient for runaway suppression in ITER (Papp et al. 2011). Perturbations at the edge might have an effect on scenarios with an off-axis runaway current profile, which can arise in the case of massive material injection . We investigated such a scenario, using orbit-following simulations with ASCOT5 (Särkimäki et al. 2020) to determine the diffusion coefficients for energetic electrons in an ITER-like plasma with a stochastic region at the plasma edge. Disruption simulations with these diffusion coefficients show that the final runaway current can be reduced, but not suppressed completely, in agreement with the conclusions of earlier work.
The analysis presented in this paper is valid when the distribution function is in quasisteady state, i.e. the runaway generation is balanced by transport, and the transport can be described by an advection-diffusion model. This assumption is not valid during the thermal quench phase in the disruption or in the presence of major magnetic islands in which a fraction of runaways could remain confined. Furthermore, a kinetic effect lacking in the model is the impact of diffusion on the avalanche generation dynamics at momentum scales close to the critical one. We anticipate that this effect could be small compared to the effects investigated so far, as the runaway electrons spend a comparatively short amount of time close to the critical momentum. Analytical progress in this direction would require the addition of a radial dimension in the full kinetic calculation with a source term to treat the dynamics close to the critical momentum, then the development of a solution to the kinetic equation valid for large momenta in the same calculation. Numerical progress, on the other hand, could be made directly by implementing the radial transport in kinetic frameworks. This would have the added benefit of capturing the pitch-angle dynamics in the presence of pitch-angle dependent transport coefficients, an effect which has not been considered in the present work.