Spatiotemporal analysis of the runaway distribution function from synchrotron images in an ASDEX Upgrade disruption

Synchrotron radiation images from runaway electrons (REs) in an ASDEX Upgrade discharge disrupted by argon injection are analyzed using the synchrotron diagnostic tool SOFT and coupled fluid-kinetic simulations. We show that the evolution of the runaway distribution is well described by an initial hot-tail seed population, which is accelerated to energies between 25-50 MeV during the current quench, together with an avalanche runaway tail which has an exponentially decreasing energy spectrum. We find that, although the avalanche component carries the vast majority of the current, it is the high-energy seed remnant that dominates synchrotron emission. With insights from the fluid-kinetic simulations, an analytic model for the evolution of the runaway seed component is developed and used to reconstruct the radial density profile of the RE beam. The analysis shows that the observed change of the synchrotron pattern from circular to crescent shape is caused by a rapid redistribution of the radial profile of the runaway density.


Introduction
Understanding runaway electron (RE) dynamics during tokamak disruptions is of utmost importance for the successful operation of future high-current tokamaks, such as ITER (Boozer 2015;Lehnen et al. 2015;Breizman et al. 2019). Disruptions are notoriously hard to diagnose, and existing numerical models cannot simultaneously capture all aspects of their temporally and spatially multiscale nature, including the associated runaway dynamics. Nevertheless, progress towards a reliable predictive capability requires model validation, and therefore finding ways of connecting experimental observations with theoretical predictions is essential.
A powerful and non-intrusive technique for diagnosing relativistic REs in tokamaks is to measure their synchrotron radiation (Finken et al. 1990;Jaspers et al. 1993). The toroidally asymmetric nature of the synchrotron radiation-due to being strongly biased in the direction of motion of the electrons-in addition to its continuum spectrum, help differentiating it from background line radiation using spectral filtering. Recent developments of synthetic synchrotron diagnostics have allowed detailed analysis of experimental synchrotron data. Most recently, in Tinguely et al. (2018aTinguely et al. ( ,b, 2019, synchrotron spectra, images and polarisation data were analysed with the help of the synthetic diagnostic Soft ) in a series of Alcator C-Mod discharges, providing valuable constraints on runaway energy, pitch angle and radial density. Full-orbit simulations have also recently provided deeper insight into observations of synchrotron radiation in 3D magnetic fields (Carbajal & del Castillo-Negrete 2017;del Castillo-Negrete et al. 2018).
In this paper, we examine synchrotron emission from REs in discharge #35628 of the ASDEX Upgrade tokamak, deliberately disrupted using an injection of neutral argon (Pautasso et al. 2016). The injection of argon leads to a rapid cooling-a thermal quench (TQ). As the TQ duration is shorter than the collision time at the critical velocity for runaway acceleration, a fraction of the most energetic electrons takes too long to thermalize and is left in the runaway region-a mechanism referred to as hot-tail generation Smith et al. 2005). During the subsequent current quench (CQ), the initially trace runaway population gets exponentially multiplied through large-angle collisions with the cold thermalized electrons in a runaway avalanche (Rosenbluth & Putvinski 1997). As the RE beam forms, its synchrotron emission can be observed using fast, wavelength-filtered visible light cameras. In this particular discharge, a sudden transition of the synchrotron image from circular to crescent shape was observed during the plateau phase. The probable cause of this spatial redistribution of the current is a magnetic reconnection caused by a (1,1) magnetohydrodynamic (MHD) mode, similarly to the observation by Lvovskiy et al. (2020) on the DIII-D tokamak, using bremsstrahlung X-ray imaging.
We briefly review the relation between the runaway electron distribution function and the observed synchrotron radiation pattern in section 2, then present the experimental setup and parameters of the ASDEX Upgrade discharge analyzed in this paper. To determine the spatiotemporal evolution of the runaway electron distribution, we use a coupled fluid-kinetic numerical tool, that takes into account the evolution of the electric field during the CQ self-consistently. This tool, based on coupling the fluid code Go (Smith et al. 2006;Fehér et al. 2011;Papp et al. 2013) which captures the radial dynamics, and the kinetic solver Code (Landreman et al. 2014;Stahl et al. 2016) that models the momentum space evolution, is presented in section 3. The collision operator used in Code includes detailed models of partial screening (Hesslow et al. 2018), which is particularly important in this case, due to the presence of a large amount of partially ionized argon.
Using the electron distribution function obtained by the coupled fluid-kinetic simulation, we show in section 3 that the resulting synchrotron radiation, computed with Soft, agrees well with the observed image, both regarding its shape, as well as the growth and spatial distribution of the intensity. Furthermore, inspired by the numerical simulations, we develop an analytical model, which is used in section 4, to reconstruct the radial density profile of the RE beam. The analysis shows that the change of the synchrotron pattern from circular to crescent shape is caused by a rapid redistribution of the radial profile of the RE density.

Synchrotron radiation from runaway electrons
When observed using a fast camera, the synchrotron radiation emitted by runaway electrons typically appears as a pattern on one side of the tokamak central column.
The size and shape of the synchrotron pattern is directly related to the energy, pitch and position of the runaway electrons-the runaway electron distribution function. Disentangling such dependencies has been the subject of several studies (Pankratov 1996;Zhou et al. 2014;Hoppe et al. 2018b,a); here we briefly review the aspects that are most important for our present purposes.

Interpretation of the synchrotron pattern
Although runaway electrons tend to occupy a large region in momentum space, the corresponding synchrotron pattern can often be well characterised with the energy and pitch angle of the particle that has the highest contribution to the camera image. We refer to these as the dominant or super particle of the pattern. We therefore define the super particle as the momentum space location (p, θ) that maximises the quantitŷ where f (r, p, θ) is the runaway electron distribution function and p 2 sin θ is the momentum space Jacobian in spherical coordinates. The Green's function G(r, p, θ) quantifies the radiation received by a camera from a particle on the orbit labelled by radius r, with momentum p and pitch angle θ (given at the point along the orbit of weakest magnetic field).
In present-day tokamaks the synchrotron spectra typically peak at infra-red wavelengths, and only a small fraction of the emitted intensity falls into the visible range. The visible-light intensity is, however, usually sufficient to be clearly distinguished from the background radiation, thus it is common to use visible-light cameras for synchrotron radiation imaging. The observed short wavelength tail of the spectrum is exponentially sensitive to the magnetic field strength, causing the fraction of emitted visible light to sometimes vary by orders of magnitude as an electron travels from the low-field side to the high-field side of the tokamak. As a result, whenever the synchrotron peak of the spectrum is located far from the spectral range of the camera, a crescent-like pattern emerges, as illustrated in figure 1a.
In a fixed detector/magnetic field setup, this effect can be thought of as an indicator of the runaway energy, as the peak wavelength of the synchrotron spectrum scales as 1/γ 2 B sin θ, where γ is the Lorentz factor of the electron and B the magnetic field strength. † A sketch of two typical pattern shapes at low and high runaway energy (relative to the camera spectral range) are shown in figure 1a and b, respectively.
Figures 1a and b also illustrate another important consequence of changing the energy, which is related to the guiding-center drift motion. At higher energies, the guidingcenter orbits shift significantly towards the outboard side of the tokamak, as does the corresponding synchrotron pattern. Although the guiding-center drift motion is routinely solved for in modern orbit following codes, accurately accounting for the effects of drifts in simulations of synchrotron radiation images is non-trivial and has, to our knowledge, previously only been employed in calculating the effect of synchrotron radiation-reaction (Hirvijoki et al. 2015). The details of the recently implemented support for guiding-center drifts in Soft are provided in Appendix A.
Synchrotron patterns are also sensitive to the spatial distribution of runaway electrons, which will be utilised in Section 4. In Soft, toroidal symmetry is assumed, along with that the poloidal transit time of a runaway is much shorter than the collision time. This  At low energies, most radiation originates from the high-field side (HFS), while at higher energies, a significant amount of radiation can also be seen on the low-field side (LFS).
(c/d) The runaway beam radius primarily determines the synchrotron pattern radius. (e/g) Illustration of how the runaway radial density affects the synchrotron pattern (darker colours indicate more radiation). (f) If the radial density is decreasing with r/a, as in (h), synchrotron radiation from low energy runaways, such as in (a), could take on a more uniform intensity distribution.
leaves the minor radius as a single spatial coordinate for the parametrization of guidingcentre orbits, taken here to be the minor radius r where the electron passes through the outer midplane. The larger the radius of the runaway beam, the larger the size of the corresponding synchrotron pattern, as illustrated in figures 1c and d. In a simulation with a runaway population distributed radially in a series of rectangle functions, as in figure 1g, the synchrotron pattern will be similar to the sketch in figure 1e, where darker colours correspond to higher-and white to zero-observed intensity. Thus, each radial point r contributes a thin band of radiation, weighted by the value of the distribution function in that point. The semi-circular shape illustrates that, at low runaway energies, the radiation intensity is greater on the high-field side (as in figure 1a) and at larger radii. As a consequence, if the radial density decreases with r, as in figure 1h, the corresponding synchrotron pattern can appear to have uniform intensity across all radii, as in figure 1f.

Experimental setup
ASDEX Upgrade is a medium-sized tokamak (major radius R = 1.65 m, minor radius a = 0.5 m) located at the Max Planck Institute for Plasma Physics in Garching, Germany (Meyer et al. 2019). An overview of plasma current, electron density, electron temperature and the hard x-ray count rate in ASDEX Upgrade discharge #35628 is shown in figure 2. This circular, L-mode discharge with 2.5 MW ECRH core electron heating applied 100 ms before the disruption was deliberately triggered by injecting N Ar ≈ 0.98 × 10 21 argon atoms into the plasma at t = 1 s. About 3 MW of ICRH heating was also applied for 200 ms before the disruption, as part of a different experiment, in a configuration where the power was poorly coupled. The discharge developed a subsequent runaway plateau  with a starting current of ≈ 200 kA and a duration of ≈ 200 ms. Before the disruption, the plasma current was I p ≈ 800 kA, the on-axis toroidal magnetic field was B T = 2.5 T, the central electron temperature was T e = 4.7 keV, and the central electron density was n e = 2.6 × 10 19 m −3 . A drop in electron temperature is observed shortly before the disruption due to internal mode activity. For the experiment, a Phantom V711 fast visible camera (connected to the in-vessel optics with an image guide and housed in a shielding box near the tokamak (Yang et al. 2013)) was equipped with a narrow-band wavelength filter with central wavelength λ 0 = 708.9 nm and FWHM of 8.6 nm. The filter wavelength was chosen as to minimise background line radiation and emphasise the synchrotron radiation, which is emitted in a continuous spectrum and had a higher intensity at longer wavelengths in these plasmas. A simulation of the camera view in discharge #35628 based on a CAD model is shown in figure 3a, with details of the configuration presented in table 1. We note that due to the lack of reliable post-disruption magnetic equilibrium reconstructions, we use the more accurate pre-disruption magnetic equilibrium reconstructions from Cliste (McCarthy 1999) for our synchrotron simulations.

Numerical modelling of the runaway electron distribution
A key objective of synchrotron radiation analysis is to validate theoretical models for runaway electron dynamics. Given a runaway electron distribution function from such a model, we should require that the corresponding synthetic synchrotron radiation image matches the experimental image well with regard to pattern shape, size and intensity distribution. Failure to predict the observed synchrotron radiation pattern can provide insight into which effects are missing from the model. In this section, we will discuss the coupling of the 1D fluid code Go (Smith et al. 2006;Fehér et al. 2011;Papp et al. 2013) to the 2D kinetic solver Code (Landreman et al. 2014;Stahl et al. 2016), used to solve for the runaway electron distribution in an ASDEX Upgrade-like disruption.

Description of numerical model
We describe the evolution of the parallel electric field E in radius and time by the induction equation in a cylinder, which is solved using Go: Here r denotes the minor radius, j the plasma current density and µ 0 is the permeability of free space. The coupling to the kinetic runaway model enters in the current density, j, which is decomposed into an Ohmic component j Ω and a runaway component j RE , where σ is the electrical conductivity with a neoclassical correction (Smith et al. 2006), e is the elementary charge, v is the electron parallel velocity and f RE is the runaway electron distribution function. The distribution function is in turn calculated in every time step, at each radius, using the local plasma parameters, by solving the kinetic equation The linear collision operator C {f } accounts for collisions between electrons and (partially screened) ions (Hesslow et al. 2017(Hesslow et al. , 2018, and between electrons using a relativistic testparticle operator (Pike & Rose 2014), while the source term S ava accounts for secondary runaway electrons generated through the avalanche mechanism. We use the simplified avalanche source derived in (Rosenbluth & Putvinski 1997) as the difference to the fully conservative operator is small . Radiation losses were initially also considered, but were found to be negligible in this scenario. The use of a test-particle collision operator in equation (3.3), given in Appendix A of Hesslow et al. (2019), is what makes the combined solution of equations (3.1)-(3.3) computationally feasible. Neglecting the field particle part of the collision operator, however, also means that the Ohmic current in Code is underestimated by about a factor of two (Helander & Sigmar 2005), which we must account for in our coupled fluid-kinetic calculations with a self-consistent electric field evolution. The linear relation between j Ω and E simplifies the correction procedure. By scanning over a wide range of effective charge and temperature, it is found that the conductivity obtained with the testparticle operator is related through a multiplicative factor g(Z eff ) to the fully relativistic conductivity σ obtained by Braams & Karney (1989): Hence, in order to calculate the runaway contribution to equation (3.2), we subtract the corrected Ohmic contribution from the total current density j CODE in Code: With this approach, the runaway current contribution can be calculated without arbitrarily defining a runaway region in momentum space, while providing a more accurate estimate than assuming all runaways to travel at the speed of light parallel to the magnetic field, which is otherwise usually done in Go and other fluid codes.
To compute the synchrotron radiation observed by the visible camera from the population of electrons calculated using the model above, we use the synthetic diagnostic tool Soft (Hoppe et al. 2018b). This tool calculates, for example, a synchrotron image by summing contributions from all parts of real and momentum space and weighting them with the provided distribution function. To reduce memory consumption, phase space is parameterised using guiding-centre orbits. Soft can also be used to calculate so-called radiation Green's functions G(r, p, θ), as introduced in equation (2.1), which relate phase-space densities to measured diagnostic signals. This mode of running Soft is used extensively for the backward modelling in section 4.

Plasma parameters used in the numerical simulations
Runaway electrons in ASDEX Upgrade disruptions, such as #35628, are typically generated through a combination of the hot-tail and avalanche mechanisms (Insulander Björk et al. 2020). The avalanche exponentiation factor is robust and depends mainly on the change in the poloidal flux profile, which can be estimated from measurements of the plasma current. The most important parameter to fit is therefore the density profile of the post-TQ seed runaway electrons, which determines the magnitude of the final plasma current. Simulations with Go+Code show that taking into account all the hottail electrons would overestimate the final runaway current with approximately a factor of four. The reason for this is that due to the presence of intense magnetic fluctuations during the TQ, a large part of the hot-tail runaway seed is likely to be deconfined, and the corresponding radial losses are not taken into account in the model. We therefore choose the seed population such that the final runaway current is matched.
With this picture in mind, we take the following approach to modelling discharge #35628: (i) First, we perform a purely fluid modelling of the thermal quench using Go, to obtain the initial electric field evolution. Then, we start the combined kinetic-fluid simulation after the thermal quench, thereby effectively disabling the "natural" hottail generation otherwise obtained in the suddenly cooling plasma. Thus, as an initial condition for the post TQ distribution function, we assume that, at each radius, it consists of current-carrying thermal electrons, along with a smaller electron population, representing the hot-tail that is uniform in radius and Gaussian in the momentum p, centred at p = 3m e c, p ⊥ = 0, and with standard deviation ∆p = 3m e c.
(ii) The evolution of the temperature during the TQ itself is taken from the experiment. While the uncertainties of this data are large, we have confirmed that, as a result of prescribing the runaway seed to match the final plasma current, our final results are not sensitive to the details of the temperature evolution. Furthermore, even though the temperature evolution affects the self-consistent electric field evolution during the TQ, the final runaway current is mainly sensitive to the time-integrated electric field, which is independent of temperature evolution. The post-disruption temperature is therefore taken to be T = 5 eV throughout the plasma, which is largely consistent with simulated values during the current quench of ASDEX Upgrade disruptions induced by argon gas injection (Insulander Björk et al. 2020). Although the temperature is expected to drop to a significantly lower value during the runaway plateau, it only weakly impacts the runaway dynamics and consequently we neglect this effect here.
(iii) We assume that neutral argon atoms with density n Ar = 0.83 × 10 20 m −3 are uniformly distributed in radius at the beginning of the simulation, and neglect their radial transport throughout the simulation. The densities of the various ionization states, and the corresponding electron density, are calculated assuming an equilibrium between . Time evolution of the electron energy spectrum (pitch-averaged distribution function) at the magnetic axis. The runaway electron seed starts close to p = 5mec at t = 0, and is then quickly accelerated to above p = 100mec within a few milliseconds. During the remainder of the runaway plateau, the initial seed sits around p = 100mec while new runaway production is dominated by large-angle collisions, causing the energy spectrum to slowly approach an exponential. The distribution also contains a thermal Maxwellian component, but due to its low temperature, it only appears as a vertical line at p = 0mec in this figure.
ionization and recombination; that is, where n tot k is the total density of species k. I i k denotes the electron impact ionization rate and R i k the radiative recombination rate for the i th charge state of species k, respectively. The ionization and recombination rates are extracted from the Atomic Data and Analysis Structure (ADAS) database (Summers et al. 2007;Summers & O'Mullane 2011).

Simulation results
The time evolution of the electron energy spectrum in the Go+Code simulation of discharge #35628 is shown in figure 5. We only evolve the simulation until around t = 1.029 s, before the synchrotron pattern suddenly changes. The only plausible mechanism that could cause such a rapid transition of the pattern is a relaxation of the current profile in fast magnetic reconnection (Igochine et al. 2006;Papp et al. 2011), in relation to an internal magnetohydrodynamic instability-a physical mechanism beyond the modelling capabilities of the Go+Code tool.
As shown in figure 5, the seed runaway population is quickly accelerated to a maximum energy during the current quench, which lasts for approximately 4 ms. During this phase, a population of secondary runaways gradually builds up, overtaking the plasma current. The maximum energy varies across radii-from p ≈ 100m e c in the core to only a few m e c at the edge-as it primarily depends on the magnitude of the induced electric field during the disruption, which, in turn, depends on the change in the current profile; in this scenario the maximum energy decreases monotonically with radius from its maximum value on the magnetic axis. After the current quench, during the remainder of the simulation, the seed electrons remain around this maximum energy as the electric field has dropped to the low level required to sustain the runaway current.
Using Soft, we compute the synchrotron radiation observed from the distribution of Although the synthetic synchrotron pattern is larger than the experimental pattern, the overall shape of the two patterns is the same, indicating that the overall runaway dynamics are well explained by Go+Code.
electrons calculated with Go+Code. The resulting synthetic camera image at t = 1.029 s, just before the pattern transition, is presented in figure 6c (along with two preceding times in figure 6a-b). Comparing the synthetic image to the experimental image in figure 6d, we find qualitative agreement, with both synchrotron patterns taking a round shape. The synthetic pattern is, however, significantly larger than the experimental pattern, both horizontally and vertically. The size of the synchrotron pattern is directly related to the radial density of runaway electrons, suggesting that the radial runaway density profile is more sharply peaked in the experiment than in the Go+Code simulation, which has a flat runaway density profile. An explanation for this discrepancy could be that the assumed runaway seed profile differs from that in the experiment, or that the initial current profile-which would experience flattening during the current spike in the thermal quench-is different. Deviations in plasma parameters such as the density and temperature could also have an effect on the avalanche gain during the current quench, as could the radial transport, which we do not model. All of these are assumed parameters in our model, due to the lack of low uncertainty experimental data. These parameters could be adjusted to give a better matching radial distribution of synchrotron radiation. However, improving agreement this way would be both computationally expensive and of limited value in better illuminating the underlying physics, so we leave this exercise for future studies. Instead, we turn our attention to the source of the observed radiation; a closer analysis reveals that most of it is emitted by an ensemble of particles that originally constituted the hot-tail seed. As was shown in figure 5, these electrons were rapidly accelerated during the current quench and then remained at their peak energy. In figures 7a and b, we show the synchrotron radiation observed from the particles associated with the r/a = 0.37 and r/a = 0.53 flux surfaces, respectively. By comparing the origin of the radiation at these radii with the local momentum space distributions in figures 8a and b respectively, we find that the region of momentum space that dominates synchrotron emission at each radius coincides with the location of the local seed population. Hence we conclude that it is the remnant seed runaways that dominate synchrotron radiation in these simulations.
When integrating over all radii, a wider dominant region appears in momentum space, as shown in figures 7c and d for an early and a late simulation time, respectively. A comparison of the emission at the two times reveals that the dominant region moves towards greater perpendicular momentum as time passes. This is caused by collisional pitch-angle scattering which increases the average perpendicular momentum in the distribution. As a result of the increased perpendicular momentum, the runaways emit  more synchrotron radiation, leading to a gradual increase of the total intensity in the camera images. The change in pitch-angle, however, is sufficiently small to not affect the synchrotron pattern shape significantly. Figure 9 shows the time evolution of the total intensity in the simulated (solid black line) and the experimental (dashed red) images, respectively. Although both intensities increase steadily, they do so at slightly different rates. This can be explained by a discrepancy in the argon density used for the simulations. As we show in appendix B, kinetic theory predicts that electrons with momentum p have an exponential pitch-angle dependence in the disruption plateau phase, f ξ (ξ) ∼ exp(Cξ), with C a time-dependent constant, and ξ = p /p. During the runaway plateau, C is roughly inversely proportional to time until it reaches an equilibrium value of 0.1p. In appendix B we show that the pitch parameter C evolves approximately as where n Ar,20 is the argon density, measured in units of 10 20 m −3 . The emitted synchrotron power at a frequency ω can furthermore be approximated by the contribution from the strongest emitting particle of such a distribution which is proportional to with γ = m 2 e c 2 + p 2 , and p the momentum reached by the hot tail seed. The free parameters in this expression are γ , n Ar , and the unknown prefactor, and it should therefore in principle be possible to fit this expression to the curves in figure 9, assuming a constant background plasma parameters and no radial transport. Unfortunately, however, such a fit can be rather ill-conditioned when the data is nearly linear, as is the case here. This is partly due to the relatively short time before the transition in the synchrotron patterns happens, which does not allow for significant pitch angle relaxation. Therefore, in practice, it is not possible to extract a reliable estimate of n Ar in this case. Nevertheless, the ability for (3.8) to fit both curves in figure 9 lends credibility to the physical picture obtained from Go+Code and may be used to estimate the impurity density in other experiments in the future.

Backward modelling
The fluid-kinetic model described in section 3 appears to capture the runaway evolution during the first part of the runaway plateau phase fairly well, but it does not contain the physics necessary to describe the sudden synchrotron pattern transition occurring in the experiment at t ≈ 1.030 s. In this section, we instead analyse the synchrotron radiation images directly and extract information from the camera images using a regularised, direct inversion. Without further constraints a direct inversion of the synchrotron image would be an ill-posed problem, thus we derive an analytical model for the dominant part of the runaway electron distribution building on the results of section 3.3, which allows us to better constrain runaway parameters.

Inversion procedure
We may capitalise on what we have learned from the fluid-kinetic simulations of section 3, in order to find constraints to regularise the inversion of the synchrotron images. We found that the synchrotron radiation pattern seen in the camera images was dominated by the remnant runaway electron seed population, situated at some maximum energy and slowly relaxing towards a steady-state pitch distribution. This evolution suggests an accelerated seed electron distribution function of the form where f r (r) is an arbitrary function describing the radial runaway density, and p , ∆p and C are free fitting parameters. Approximating the Gaussian with a delta function is motivated by the fact that the synchrotron pattern is usually rather insensitive to variations in the runaway energy distribution. The same argument is used to fully decouple the spatial coordinate r from the momentum parameter p and, since C is mainly a function of p , also from C. In this model we thus assume for simplicity that the momentum and pitch distributions are the same across all flux surfaces. Our inversion method utilises the capability of Soft to generate weight functions for a given tokamak/detector setup, as described in section 2. The brightness of pixel i in a synchrotron image I i is related to the distribution function f (r, p, ξ) = f r (r)f pξ (p, ξ) through the weight function G(r, p, ξ) as I i = G i (r, p, ξ)f (r, p, ξ) p 2 dr dp dξ. (4.2) By representing the image and the discretised radial runaway density profile as vectors, the discretised version of the equation system (4.2) can be formulated as where we have introduced the reduced weight function matrix (4.4) with ∆r k = |r k+1 − r k |. Given a momentum-space distribution function f pξ (p, ξ), which can be characterised with the parameters p and C in equation (4.1), we thus seek to minimise the sum of squares of differences between pixels in the synthetic and experimental images I i and I exp i . Since the problem is still ill-posed, we regularise it using the Tikhonov method (Press et al. 2007), resulting in f r (r) = arg min where the Tikhonov matrix Γ ij is taken to be the identity matrix and the Tikhonov parameter α is determined using the L-curve method (Hansen & O'Leary 1993). The ability to solve this problem using a linear least-squares method allows us to efficiently explore the space of possible combinations of (p , C), and hence to use typical minimization methods to solve the full problem.

Inverted distribution function
One of the main reasons for applying backward modelling to this ASDEX-U discharge is to better understand what gives rise to the synchrotron pattern transition between t = 1.029 ms and t = 1.030 ms, corresponding to the video frames 3b and c. To see that the transition is not due to a change of the energy distribution, we can estimate the energy gained by the runaways in a magnetic reconnection event with the following argument: if the plasma current changes by an amount ∆I during the event, then it follows from a circuit equation for the plasma that the energy of the electrons changes according to where L is the self-inductance and R 0 is the major radius of the torus. Assuming L ∼ µ 0 R 0 , with µ 0 the vacuum permeability, the induced energy is found to increase by 60 eV for every ampere decrease in the plasma current. In our case, where the second current spike rises by ∆I ≈ 5 kA, the energy transferred to the electrons should be below 0.5 MeV, which is small compared to typical runaway energies at t = 1.029 ms (≈ 25 MeV at mid-radius, see figure 8b). Furthermore, since the energy is gained exclusively in the parallel direction, and the pitch angles are initially small, the pitch angles should also be negligibly affected; ∆θ ≈ −(∆W /W )θ. Hence, we expect the observed synchrotron pattern transition to be caused by a spatial redistribution of runaways. In what follows, we will therefore seek the best fit between theory and experiment for both frames 3a and b simultaneously, assuming p and C to remain unchanged in the transition. The sum of pixel differences squared, as a function of the fitting parameters p and C, is shown in figure 10. In each point, the best radial density is constrained using equation (4.5). Optimal agreement is obtained with p = 57.5m e c and C = 45, although the region of good agreement in figure 10 is fairly large. However, most of the optimal combinations of p and C yield approximately the same value for the dominant pitch angle, θ ≈ 0.30 rad. This is in agreement with the Go+Code and Soft simulations presented in section 3.3 which had p ≈ 55m e c, C ≈ 25 (at p = p ) and θ ≈ 0.36 rad at t = 1.029 ms. Since C is inversely proportional to the (relatively poorly diagnosed) argon density n Ar , these are well within the uncertainties of the inversion and in the plasma parameters.
The radial density profiles obtained in the inversion for the two video frames are presented in figure 11a. Note that the radial coordinate denotes the particle position along the outer midplane, and that r = 0 corresponds to the magnetic axis. Since particles are counted on the outboard side, and since only particles with p = p are considered, the inverted radial profile is exactly zero within the grey region r/a ≈ 0.1, corresponding to the drift orbit shift for these particles. The blue and red shaded regions in figure 11a indicate the maximum deviation of density profiles corresponding to the optima of all combinations (p , C) with normalised likeness less than 2 in figure 10. Since the radial density profiles contain an uninteresting scaling factor in order for the inversion algorithm to match the absolute pixel values in both experimental and synthetic images, we rescale all density profiles by a scalar multiplicative factor before Normalised likeness Figure 10. Sum of pixel differences squared for both of frames 11a and b, given different combinations of p and C. For each combination of p and C, the corresponding optimal radial density is calculated and used for comparing the images. The global minimum is marked with a green cross and is located in p = 57.5 and C = 45. The white regions correspond to unreasonable combinations of the two parameters.
evaluating the maximum deviation. The maximum deviations suggest that although the uncertainty in (p , C) is relatively large, the radial density profiles are somewhat more robust. Specifically, the analysis shows that the synchrotron pattern transition must be due to a spatial redistribution of particles.
The radial density profile inversion reveals that the synchrotron pattern change is caused by a rapid redistribution of the runaway electron density. As described in section 2, the shifted density profile causes the particles on the high-field side to collectively emit relatively more radiation than those on the low-field side, yielding a crescent pattern after the event. By integrating the radial density profiles in figure 11a over the tokamak volume, one finds that approximately 14% of particles are predicted to be lost in the event.
A redistribution of the runaway density profile is also consistent with the occurrence of magnetic reconnection, which, in addition to causing the observed current spike, would redistribute the current profile. We do not have sufficient experimental information to perform an MHD stability analysis, thus it is uncertain exactly how this reconnection event is triggered. However, the fact that runaway generation is often more efficient in the plasma core , it is plausible that a corresponding increase in the peaking of the current profile and decreasing central safety factor, could destabilize an internal kink instability, consistent with the strong (1,1) mode activity evident in 4b.

Summary
We analyzed visible-light camera images of synchrotron emission from runaway electrons in the ASDEX Upgrade discharge 35628 that was deliberately disrupted by argon massive gas injection. In the runaway plateau phase of this particular discharge, a sudden transition of the synchrotron image from circular to crescent shape was observed, correlated with a spike in the plasma current as well as (m, n) = (1, 1) MHD activity.
We simulated the spatio-temporal evolution of the runaway distribution function with a coupled kinetic-fluid code taking into account the evolution of the electric field selfconsistently. This distribution function was then used as an input to the Soft synthetic diagnostic tool to calculate the shape and intensity evolution of the synchrotron images. We find that a hot-tail seed population of runaway electrons was multiplied by close collisions, resulting in a distribution with an exponentially decreasing energy spectrum. The remnant seed was accelerated to high energies and was responsible for most of the synchrotron emission. We derived an analytical expression for the time evolution of the total emitted synchrotron power and found qualitative agreement with the observed intensity evolution. By constraining the distribution function in momentum space we inverted the radial profiles of the runaway beam and found that the change from circular to crescent shape was caused by a redistribution of the radial profile of the RE density. In particular the momentum-space distribution was found to be difficult to extract from the images, and to better constrain it one could potentially combine several independent diagnostic signals, such as synchrotron spectra, bremsstrahlung spectra, and camera images that view the plasma from different angles, or are sensitive to different wavelength ranges.
A number of directions for future studies could be envisaged based on the results of this study. Radial transport is expected to significantly affect runaway dynamics, both in present-day as well as future tokamaks, however this has often been neglected in runaway models due to its complexity, including in this work. Finite-aspect-ratio effects have also been neglected in the fluid-kinetic modelling part of this study, but could play a role in the evolution of the runaway electron distribution function. Finally, the "two-component" runaway population observed in ASDEX Upgrade discharge #35628, studied here, consisting of a high-energy remnant seed and a current-carrying avalanche component, should be investigated in different devices and discharges. If the situation observed in ASDEX Upgrade is representative for many devices, the two-component model used here for fitting synchrotron radiation could be used more widely to infer runaway parameters at many experiments.

Appendix A. Accounting for guiding-centre drifts
To reduce the size of the phase-space from six to three dimensions, Soft utilises a guiding-centre transformation which eliminates the orbit time, as well as the toroidal and gyro-phase angles from the distribution function. However, the dependence on the particle's position and momentum components still remains in the radiation emission function, dP/dΩ; these must be computed within the guiding-centre framework. The runaway electron drift orbit shift can however be significant, and so the effect is of importance in many, or even most, practical scenarios. In this section we describe the incorporation of guiding-centre orbit drifts in Soft.

A.1. Electron momentum vector
Formally, the guiding-centre drift motion appears in the standard theory as an O( ) effect, where is the usual guiding-centre ordering parameter. An expression for the Larmor radius vector to first order was given by Hirvijoki et al. (2015), and using expressions found therein one can also derive the following expression for the momentum vector of a particle with mass m and charge q (negative for electrons): with the guiding-center momentum vector where p and p ⊥ are the particle momenta parallel and perpendicular to the magnetic field respectively, B is the magnetic field strength,b andρ are mutually perpendicular unit vectors in the direction of the magnetic field and (lowest order) Larmor radius vector, respectively,⊥ =ρ ×b, ρ = p ⊥ /|q|B is the Larmor radius, µ = p 2 ⊥ /2mB is the magnetic moment, κ = (b · ∇)b is the inverse curvature vector, B = B + p b · (∇ ×b)/q, and the dyads a 1 and a 2 are defined as To lowest order, the momentum vector (A 1) describes circular gyro motion around the magnetic field line. However, when including O( ) terms, the gyro motion component picks up contributions which alter this circular motion. A detailed analysis reveals that the first-order terms primarily shift the axis-of-rotation of the gyro motion from the magnetic fieldb to the guiding-centre direction of motion P /P . Hence, to first order in guiding-centre theory, synchrotron radiation is emitted in an approximately circular cone, with opening angle equal to the pitch angle θ = arctan(p ⊥ /p ), around the guiding-centre momentum vector P .

A.2. Implementation and validation
The conclusions from the analysis of the electron momentum vector suggest that the cone model, implemented previously Soft , can be modified to take guiding-centre drifts into account by including the first-order drift terms in the guidingcentre momentum vector P , and consider the radiation to be emitted in an infinitesimally thin circular cone around P . The method is however not exact since the O( ) terms in (A 1) can break gyrotropy, making the emission cone non-circular. Therefore, when we implement the modified cone model in Soft, we also implement a hybrid fullorbit/guiding-centre model that can indicate when the cone model assumptions are violated.
Soft computes the radiation observed from an ensemble of particles with distribution function f gc (r, p) by evaluating the integral Here, r cp = x 0 − x is the vector extending from the detector at x 0 to the particle at x = x(r, τ, φ), p (0) is the particle momentum at the outer mid-plane,n is the detector surface normal, Θ is a step function that is one whenever r cp /r cp lies in the detector field-of-view, dP/dΩ is the angular distribution of radiation emitted by each particle, and J is the Jacobian for the full coordinate transformation (including the guidingcentre transformation). The spatial coordinates used are the minor radius of the flux surface in the outer midplane, r; the time coordinate along a guiding-centre orbit τ ; and the toroidal angle φ. The momentum p is therefore a function of both x and p (0) . Note that although we use a guiding-centre transformation, the integrand (A 4) is still dependent on the particle position and momentum x and p. This means that we have to use equation (A 1) and the corresponding expression for the Larmor radius vector ρ to consistently evaluate (A 4) to O( ). However, since it would be prohibitively computationally expensive to introduce the drifts exactly in the cone model, we use the observation that the primary effect of the O( ) terms is to rotate the circular cone of radiation around the guiding-centre velocity vector, thus we simply modify the cone model according to whereb denotes the local magnetic field unit vector and cos θ = p ·b.
To verify that the approximation (A 5) gives sufficiently accurate result, we now also evaluate (A 4) using a hybrid approach. Instead of calculating a set of guiding-centre orbits which are used to evaluate (A 4) together with equation (A 1), we calculate a set of particle orbits, from which r cp and p are obtained exactly. While this approach may at first seem to be more accurate, we point out that equation (A 4), which is still used in this approach, was derived by making a guiding-centre transformation. Thus, the use of particle orbits to evaluate (A 4) is not expected to be quantitatively accurate and, in addition, the Jacobian J causes the amount of detected radiation to be weighed  Figure 12. Outline of the synchrotron pattern calculated using the particle orbit method (red), the modified cone model (black, solid), and the cone model neglecting guiding-centre drifts (black, dashed), at fixed pitch angle (θ = 0.2 rad) and two different momenta. In a synchrotron radiation image, every pixel enclosed by the contour would be illuminated by radiation. The main effect of the drifts is to compress the synchrotron pattern towards the low-field side, corresponding to the drift orbit shift direction.
incorrectly. The purpose of the validation is therefore to ensure that the size, shape and position of the radiation pattern is correct, and this only requires that the guiding-centre drifts and gyro-orbit can be evaluated correctly. If the assumptions of the modified cone model (A 5) are violated, we expect the synchrotron pattern shapes calculated by the two models to deviate significantly. Figure 12 compares the regular cone model, the modified cone model and the particle orbit model at two different momenta and a fixed pitch angle θ = 0.2 rad. The figures show the outline contours of a number of synchrotron images, in order to emphasize the size, shape and position of the radiation pattern. The effect of the drifts is primarily to compress the synchrotron pattern, as seen by comparing the dashed (no drifts) and solid (drifts included) curves. The modified cone model agrees mostly with the particle orbit model. As the energy is increased, the pattern resulting from the modified cone model also begins to deviate from the particle orbit model, mainly along the upper edge. This deviation is a sign of that the guiding-centre is breaking down, and gyrotropy is being violated. As an additional consequence of the theory breaking down, we observe unphysical oscillations along the lower right edge in figure 12b. From figure 12 we can however conclude that the modified cone model captures the dominant effects of first order guiding-centre corrections, allowing us to more accurately model high-energy runaway scenarios without a significant increase in numerical complexity.
Appendix B. Estimation of p and θ evolution from kinetic theory By constraining the time evolution of the runaway pitch-angle distribution using kinetic theory, the evolution of the super-particle parameters, p and θ , can be inferred from the time variation of the experimentally measured total intensity-as done in figure 9.
To obtain a simplified description of the synchrotron radiation emitted by runaways during the runaway plateau phase, we start with the expression for the emission from an When p ⊥ m e c, the critical frequency-where the emission is the strongest-is approximately given by ω c ≈ 3eBγ 2 θ/(2m e ). Considering an observed wavelength of λ ≈ 700 nm, a representative pitch angle θ ≈ 0.2, and energy γ ≈ 50, we find that ω/ω c ≈ 8. Then, the emission is well described by the formula evaluated at the ω/ω c 1 limit ∂P ∂ω ≈ m e cr 0 √ 3π √ ω c ω γ 2 e −ω/ωc .

(B 2)
The total emission from a distribution of runaways, which we assume can be approximated as f RE ≈ F (t, p) exp(−C(t, p)θ 2 /2), is given by the phase-space integral of up to a factor depending on plasma and runaway parameters at most polynomially (in contrast to the dominant, exponential dependence captured by P). At a given p (the precise value of which depends on the energy distribution of runaways), this function is maximized by the pitch angle θ , and takes the maximum value P(θ ), given by At the end of the current quench, when the plasma current is turning into its plateau phase, the electric field approaches a critical value where the runaway population is marginally sustained (Breizman 2014). In that case, the momentum flux in phase space becomes negligible, and the subsequent evolution throughout the plateau phase is mainly pitch angle relaxation, described by the equation