1. 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 Reference Boozer2015; Lehnen et al. Reference Lehnen, Aleynikova, Aleynikov, Campbell, Drewelow, Eidietis, Gasparyan, Granetz, Gribov and Hartmann2015; Breizman et al. Reference Breizman, Aleynikov, Hollmann and Lehnen2019). 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. Reference Finken, Watkins, Rusbüldt, Corbett, Dippel, Goebel and Moyer1990; Jaspers et al. Reference Jaspers, Finken, Mank, Hoenen, Boedo, Cardozo and Schüller1993). 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 in 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. (Reference Tinguely, Granetz, Hoppe and Embréus2018a,Reference Tinguely, Granetz, Hoppe and Embréusb, Reference Tinguely, Hoppe, Granetz, Mumgaard and Scott2019), synchrotron spectra, images and polarization data were analysed with the help of the synthetic diagnostic Soft (Hoppe et al. Reference Hoppe, Embréus, Tinguely, Granetz, Stahl and Fülöp2018b) 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 three-dimensional magnetic fields (Carbajal & del Castillo-Negrete Reference Carbajal and del Castillo-Negrete2017; del Castillo-Negrete et al. Reference del Castillo-Negrete, Carbajal, Spong and Izzo2018).
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. Reference Pautasso, Bernert, Dibon, Duval, Dux, Fable, Fuchs, Conway, Giannone and Gude2016). 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 (Helander et al. Reference Helander, Smith, Fülöp and Eriksson2004; Smith et al. Reference Smith, Helander, Eriksson and Fülöp2005). During the subsequent current quench (CQ), the initial trace runaway population gets exponentially multiplied through large-angle collisions with the cold thermalized electrons in a runaway avalanche (Rosenbluth & Putvinski Reference Rosenbluth and Putvinski1997). 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, similar to the observation by Lvovskiy et al. (Reference Lvovskiy, Paz-Soldan, Eidietis, Aleynikov, Austin, Molin, Liu, Moyer, Nocente and Shiraki2020) on the DIII-D tokamak, using bremsstrahlung X-ray imaging.
We briefly review the relation between the RE distribution function and the observed synchrotron radiation pattern in § 2, then present the experimental set-up and parameters of the ASDEX Upgrade discharge analysed in this paper. To determine the spatiotemporal evolution of the RE 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. Reference Smith, Helander, Eriksson, Anderson, Lisak and Andersson2006; Fehér et al. Reference Fehér, Smith, Fülöp and Gál2011; Papp et al. Reference Papp, Fülöp, Fehér, de Vries, Riccardo, Reux, Lehnen, Kiptily, Plyusnin and Alper2013) which captures the radial dynamics, and the kinetic solver Code (Landreman, Stahl & Fülöp Reference Landreman, Stahl and Fülöp2014; Stahl et al. Reference Stahl, Embréus, Papp, Landreman and Fülöp2016) that models the momentum space evolution, is presented in § 3. The collision operator used in Code includes detailed models of partial screening (Hesslow et al. Reference Hesslow, Embréus, Hoppe, DuBois, Papp, Rahm and Fülöp2018), 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 § 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 § 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.
2. Synchrotron radiation from REs
When observed using a fast camera, the synchrotron radiation emitted by REs 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 RE distribution function. Disentangling such dependencies has been the subject of several studies (Pankratov Reference Pankratov1996; Zhou et al. Reference Zhou, Pankratov, Hu, Xu and Yang2014; Hoppe et al. Reference Hoppe, Embréus, Paz-Soldan, Moyer and Fülöp2018a,Reference Hoppe, Embréus, Tinguely, Granetz, Stahl and Fülöpb); here we briefly review the aspects that are most important for our present purposes.
2.1. Interpretation of the synchrotron pattern
Although REs tend to occupy a large region in momentum space, the corresponding synchrotron pattern can often be well characterized 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,\theta )$ that maximizes the quantity
\begin{equation} \hat{I} = G(r,p,\theta) f(r,p,\theta) p^2\sin\theta, \end{equation}
where 
$f(r,p,\theta )$ is the RE distribution function and 
$p^2\sin \theta$ is the momentum space Jacobian in spherical coordinates. The Green's function 
$G(r,p,\theta )$ quantifies the radiation received by a camera from a particle on the orbit labelled by radius 
$r$, with momentum 
$p$ and pitch angle 
$\theta$ (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 (Hoppe et al. Reference Hoppe, Embréus, Paz-Soldan, Moyer and Fülöp2018a), causing the fraction of emitted visible light to sometimes vary by orders of magnitude as an electron travels from the low-field side (LFS) to the high-field side (HFS) 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 1(a).
Figure 1. (a,b) Illustration of how the runaway energy affects the observed synchrotron pattern. At low energies, most radiation originates from the HFS, while at higher energies, a significant amount of radiation can also be seen on the 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 panel (h), synchrotron radiation from low energy runaways, such as in panel (a), could take on a more uniform intensity distribution.
In a fixed detector/magnetic field set-up, this effect can be thought of as an indicator of the runaway energy, as the peak wavelength of the synchrotron spectrum scales as 
$1/(\gamma ^2 B\sin \theta)$, where 
$\gamma$ is the Lorentz factor of the electron and 
$B$ the magnetic field strength. (While the peak also depends on pitch angle, the pitch angle additionally alters the vertical and toroidal extent of the pattern, thus clearly distinguishing a change in energy from a change in pitch angle.) A sketch of two typical pattern shapes at low and high runaway energy (relative to the camera spectral range) are shown in figures 1(a) and 1(b), respectively.
Figures 1(a) and 1(b) also illustrate another important consequence of changing the energy, which is related to the guiding-centre drift motion. At higher energies, the guiding-centre orbits shift significantly towards the outboard side of the tokamak, as does the corresponding synchrotron pattern. Although the guiding-centre 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. Reference Hirvijoki, Decker, Brizard and Embréus2015). The details of the recently implemented support for guiding-centre drifts in Soft are provided in appendix A.
Synchrotron patterns are also sensitive to the spatial distribution of REs, which will be utilized in § 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 leaves the minor radius as a single spatial coordinate for the parameterization of guiding-centre 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 1(c) and 1(d). In a simulation with a runaway population distributed radially in a series of rectangle functions, as in figure 1(g), the synchrotron pattern will be similar to the sketch in figure 1(e), 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 semicircular shape illustrates that, at low runaway energies, the radiation intensity is greater on the HFS (as in figure 1a) and at larger radii. As a consequence, if the radial density decreases with 
$r$, as in figure 1(h), the corresponding synchrotron pattern can appear to have uniform intensity across all radii, as in figure 1(f).
A large variety of synchrotron patterns have been reported in the literature, although circular and crescent patterns seem to be among the more common ones. In this paper, and in particular in § 4, we will analyse the transition from a circular synchrotron pattern into a crescent pattern and find, as suggested above, that the transition is due to a redistribution of the runaway density. Similar transitions have been observed before, most recently by Lvovskiy et al. (Reference Lvovskiy, Paz-Soldan, Eidietis, Aleynikov, Austin, Molin, Liu, Moyer, Nocente and Shiraki2020) who observe a similar submillisecond transition as observed here. Earlier reports also show that transitions from ellipses to crescents, and vice versa, can occur over longer time scales in both disruptions (Hollmann et al. Reference Hollmann, Austin, Boedo, Brooks, Commaux, Eidietis, Humphreys, Izzo, James and Jernigan2013) and quiescent flat-top plasmas (England et al. Reference England, Chen, Seo, Chung, Leev, Yoo, Kim, Bae, Jeonv and Kwak2013).
2.2. Experimental set-up
ASDEX Upgrade is a medium-sized tokamak (major radius 
$R = {1.65}\ \textrm {m}$, minor radius 
$a = {0.5}\ \textrm {m}$) located at the Max Planck Institute for Plasma Physics in Garching, Germany (Meyer et al. Reference Meyer, Angioni, Albert, Arden, Parra, Asunta, de Baar, Balden, Bandaru and Behler2019). 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\ \textrm {MW}$ ECRH core electron heating applied 
${100}\ \textrm {ms}$ before the disruption was deliberately triggered by injecting 
$N_\textrm {Ar}\approx 0.98\times 10^{21}$ argon atoms into the plasma at 
$t={1}\ \textrm {s}$. Due to the circular plasma shape, the plasma was vertically stable during the discharge, consistent with diagnostic camera recordings. Approximately 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 
${\approx } {200}\ \textrm {kA}$ and a duration of 
${\approx }{200}\ \textrm {ms}$. Before the disruption, the plasma current was 
$I_{p} \approx {800}\ \textrm {kA}$, the on-axis toroidal magnetic field was 
$B_{T} = {2.5}\ \textrm {T}$, the central electron temperature was 
$T_e = {4.7}\ \textrm {keV}$, and the central electron density was 
$n_e = 2.6\times 10^{19}\ \textrm {m}^{-3}$. A drop in electron temperature is observed shortly before the disruption due to internal mode activity.
Figure 2. Overview of the most relevant plasma parameters in ASDEX Upgrade discharge #35628. (a) Total plasma current, with the smaller, zoomed-in figure showing a small secondary current spike, (b) line-averaged electron density from central chord 
$\textrm {CO}_2$ interferometry, (c) electron temperature from central electron cyclotron emission (note that the temperature decreases somewhat just before the disruption to approximately 4.7 keV), (d) ex-vessel hard X-ray counts.
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. Reference Yang, Krieger, Lunt, Brochard, Briancon, Neu, Dux, Janzer, Potzel and Pütterich2013) was equipped with a narrow-band wavelength filter with central wavelength 
$\lambda _0 = {708.9}\ \textrm {nm}$ and full width at half-maximum (FWHM) of 
${8.6}\ \textrm {nm}$. The filter wavelength was chosen as to minimize background line radiation and emphasize 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 3(a), 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 predisruption magnetic equilibrium reconstructions from Cliste (McCarthy Reference McCarthy1999) for our synchrotron simulations.
Figure 3. (a) Simulated camera view of the Phantom v711 fast camera in the configuration used for discharge 35628. (b,c) Synchrotron radiation images at (b) 
$t = 1.029\ \textrm {s}$ and (c) 
$t = 1.030\ \textrm {s}$, observed using a filtered visible light camera in ASDEX Upgrade during discharge 35628. A sudden, submillisecond transition from a circular to a crescent shape is observed around 
$t = {1.030}\ \textrm {s}$.
Table 1. Parameters of the image recorded by a Phantom v711 visible light camera, which was used for synchrotron radiation imaging. Only parameters relevant to synthetic diagnostic simulation are shown.
A few milliseconds after the gas has been injected, a circular synchrotron pattern appears in the visible light camera images. During the next ~20 ms, in the runaway plateau phase of the disruption, the synchrotron pattern is observed to gradually increase in brightness while maintaining approximately the same size and shape. Eventually, the pattern attains its maximum brightness near 
$t={1.029}\ \textrm {s}$, shown in figure 3(b). In the very next frame, at 
$t={1.030}\ \textrm {s}$ the synchrotron pattern has been turned into a crescent shape, as shown in figure 3(c). Around this time, a modest 5 kA spike is observed in the total plasma current, shown in figure 2(a), corresponding roughly to a 
$2.5\,\%$ increase. During this current spike, broadband transient magnetic activity is measured by the magnetic pick-up coils in frequencies ranging from a few kHz to approximately 80 kHz (see figure 4a). Mode number analysis (figure 4b,c) was performed using the NTI Wavelet Tools (https://github.com/fusion-flap/nti-wavelet-tools) program package primarily with a method based on looking for the best fitting integer mode number on the measured cross-phases as function of relative probe positions (Horváth et al. Reference Horváth, Poloskei, Papp, Maraschek, Schuhbeck and Pokol2015). The toroidal array of ‘ballooning coils’, which measures variations to the radial magnetic field, was used for the toroidal mode number analysis, applying the phase corrections of the measured transfer functions of the probes (Horváth et al. Reference Horváth, Poloskei, Papp, Maraschek, Schuhbeck and Pokol2015). The poloidal mode numbers were determined using the C09 Mirnov coil array, using the probe positions transformed into the straight-field-line coordinate system of the 
$q=2$ surface. Measured transfer functions were not available for this set of probes, so the confidence in the results had to be improved by applying a complementary method of mode number estimation that is based on the monotonization of the phase functions (Pokol et al. Reference Pokol, Papp, Por, Zoletnik and Weller2008). The analysis clearly indicated a 
$(m,n)=(1,1)$ mode propagating in the electron diamagnetic drift direction in the frequency range of 
$8\text {--}{20}\ \textrm {kHz}$. No precursor activity is observed; the only signal components preceding the event are some low frequency (
${\sim }{100}\ \textrm {Hz}$) oscillations, which is attributed to vessel and diagnostic vibrations.
Figure 4. Time-frequency analysis of the transient MHD event in the runaway plateau stage of AUG discharge #35628. (a) Representative spectrogram of a magnetic pick-up coil signal shows wide-band activity with signal energy concentrated to below 20 kHz; (b) toroidal mode numbers fitted using the MHI-B31 toroidal ballooning coil array; (c) poloidal mode numbers fitted using the MHI-C09 poloidal Minrov coil array. Mode number plots show only the good fits and only in regions of sufficient signal energy. The mode below 20 kHz has 
$(n,m)=(1,-1)$ mode numbers in machine coordinates, which corresponds to 
$(n,m)=(1,1)$ propagating in the electron diamagnetic drift direction in plasma coordinates.
3. Numerical modelling of the RE distribution
A key objective of synchrotron radiation analysis is to validate theoretical models for RE dynamics. Given a RE 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 one-dimensional fluid code Go (Smith et al. Reference Smith, Helander, Eriksson, Anderson, Lisak and Andersson2006; Fehér et al. Reference Fehér, Smith, Fülöp and Gál2011; Papp et al. Reference Papp, Fülöp, Fehér, de Vries, Riccardo, Reux, Lehnen, Kiptily, Plyusnin and Alper2013) to the two-dimensional kinetic solver Code (Landreman et al. Reference Landreman, Stahl and Fülöp2014; Stahl et al. Reference Stahl, Embréus, Papp, Landreman and Fülöp2016), used to solve for the runaway electron distribution in an ASDEX Upgrade-like disruption.
3.1. Description of numerical model
We describe the evolution of the parallel electric field 
$E_\parallel$ in radius and time by the induction equation in a cylinder, which is solved using Go:
\begin{equation} \frac{1}{r} \frac{\partial }{\partial r}\left( r\frac{\partial E_\parallel}{\partial r} \right) = \mu_0\frac{\partial j}{\partial t}. \end{equation}
Here 
$r$ denotes the minor radius, 
$j$ the plasma current density and 
$\mu _0$ is the permeability of free space. We assume that the plasma is surrounded by a perfectly conducting wall, and that the wall and plasma are separated by a vacuum region that is 8 cm wide. The coupling to the kinetic runaway model enters in the current density, 
$j$, which is decomposed into an ohmic component 
$j_\varOmega$ and a runaway component 
$j_\textrm {RE}$,
\begin{equation} j = j_\varOmega + j_\textrm{RE} = \sigma E_\parallel + e\int v_\parallel\, f_\textrm{RE}\,\mathrm{d}^3p, \end{equation}
where 
$\sigma$ is the electrical conductivity with a neoclassical correction (Smith et al. Reference Smith, Helander, Eriksson, Anderson, Lisak and Andersson2006), 
$e$ is the elementary charge, 
$v_\parallel$ is the electron parallel velocity and 
$f_\textrm {RE}$ is the RE 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,
\begin{equation} \frac{\partial f}{\partial t} + eE_\parallel\frac{\partial f}{\partial p_\parallel} = C\left\{\, f \right\} + S_\textrm{ava}. \end{equation}
The linear collision operator 
$C\{\, f \}$ accounts for collisions between electrons and (partially screened) ions (Hesslow et al. Reference Hesslow, Embréus, Stahl, DuBois, Papp, Newton and Fülöp2017, Reference Hesslow, Embréus, Hoppe, DuBois, Papp, Rahm and Fülöp2018), and between electrons using a relativistic test-particle operator (Pike & Rose Reference Pike and Rose2014), while the source term 
$S_\textrm {ava}$ accounts for secondary REs generated through the avalanche mechanism. We use the simplified avalanche source derived in (Rosenbluth & Putvinski Reference Rosenbluth and Putvinski1997) as the difference to the fully conservative operator is small (Embréus, Stahl & Fülöp Reference Embréus, Stahl and Fülöp2018). Radiation losses were initially also considered, but were found to be negligible in this scenario.
The use of a test-particle collision operator in (3.3), given in appendix A of Hesslow et al. (Reference Hesslow, Unnerfelt, Vallhagen, Embréus, Hoppe, Papp and Fülöp2019), is what makes the combined solution of (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 approximately a factor of two (Helander & Sigmar Reference Helander and Sigmar2005), which we must account for in our coupled fluid-kinetic calculations with a self-consistent electric field evolution. The linear relation between 
$j_\varOmega$ and 
$E_\parallel$ simplifies the correction procedure. By scanning over a wide range of effective charge and temperature, it is found that the conductivity obtained with the test-particle operator is related through a multiplicative factor 
$g(Z_{\textrm {eff}})$ to the fully relativistic conductivity 
$\sigma$ obtained by Braams & Karney (Reference Braams and Karney1989):
\begin{equation} \sigma_{\textrm{CODE}, \textrm{tp}} = g\left(Z_{\textrm{eff}}\right)\sigma. \end{equation}
Hence, in order to calculate the runaway contribution to (3.2), we subtract the corrected ohmic contribution from the total current density 
$j_{\textrm {CODE}}$ in Code as follows:
\begin{equation} j_\textrm{RE} = j_\textrm{CODE} - \sigma_{\textrm{CODE}, \textrm{tp}}E_\parallel = j_\textrm{CODE} - g\left(Z_{\textrm{eff}}\right)\sigma E_\parallel. \end{equation}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. Reference Hoppe, Embréus, Tinguely, Granetz, Stahl and Fülöp2018b). 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 parameterized using guiding-centre orbits. The synthetic diagnostic tool Soft can also be used to calculate so-called radiation Green's functions 
$G(r,p,\theta )$, as introduced in (2.1), which relate phase-space densities to measured diagnostic signals. This mode of running Soft is used extensively for the backward modelling in § 4.
3.2. 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 et al. Reference Insulander Björk, Papp, Embreus, Hesslow, Fülöp, Vallhagen, Lier, Pautasso and Bock2020). The avalanche exponentiation factor is robust and depends mainly on the change in the poloidal flux profile. In the ideal theory, where the avalanche growth rate is directly proportional to the electric field, and radial transport during the CQ is assumed to be negligible, the final plateau runaway current profile is completely determined by the surviving post-TQ runaway seed and plasma current profile. In this work we assume that the loop voltage is constant across flux surfaces just before the disruption, and take the initial current profile to be the corresponding ohmic current profile, which leaves the post-TQ seed profile as the main unknown of the simulation. Due to the relatively large current drop of 
${\sim }{600}\ \textrm {kA}$, we expect a significant number of avalanche multiplications to occur throughout the plasma, and therefore a relatively weak sensitivity to the chosen runaway seed density profile.
Simulations with Go 
$+$ Code show that taking into account all the hot-tail electrons obtained from kinetic theory would overestimate the final runaway current by 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 to prescribe a radially uniform seed population such that the final runaway current is matching the experimentally observed plasma current during the plateau.
With this picture in mind, we take the following approach to modelling discharge #35628.
(i) First, we perform a purely fluid modelling of the TQ using Go, to obtain the initial electric field evolution. Then, we start the combined kinetic-fluid simulation after the TQ, thereby effectively disabling the ‘natural’ hot-tail 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_\parallel =3m_ec, p_\perp = 0$, and with standard deviation $\Delta p = 3m_ec$.(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 CQ, 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}\ \textrm {eV}$ throughout the plasma, which is largely consistent with simulated values during the CQ of ASDEX Upgrade disruptions induced by argon gas injection (Insulander et al. Reference Insulander Björk, Papp, Embreus, Hesslow, Fülöp, Vallhagen, Lier, Pautasso and Bock2020). 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_\textrm {Ar} = 0.83\times 10^{20}\ \textrm {m}^{-3}$ (corresponding to 20 % of the injected atoms Pautasso et al. Reference Pautasso, Dibon, Dunne, Dux, Fable, Lang, Linder, Mlynek, Papp and Bernert2020; Insulander et al. Reference Insulander Björk, Papp, Embreus, Hesslow, Fülöp, Vallhagen, Lier, Pautasso and Bock2020) 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 ionization and recombination; that is, $n^{i}_k$ are computed from
(3.6)where\begin{equation} \left. \begin{gathered} R^{i+1}_k n^{i+1}_k-I^i_k n^i_k=0, \quad i=0,1,\ldots, Z-1,\\ \sum_i n^{i}_k=n^\mathrm{tot}_k, \quad i=0,1,\ldots, Z, \end{gathered}\right\} \end{equation}$n^{\mathrm {tot}}_k$ is the total density of species $k$. Here $I_k^{i}$ denotes the electron impact ionization rate and $R_k^i$ the radiative recombination rate for the $i\textrm {th}$ charge state of species $k$, respectively. The ionization and recombination rates are extracted from the Atomic Data and Analysis Structure (known as ADAS) database (Summers et al. Reference Summers, O'Mullane, Whiteford, Badnell and Loch2007; Summers & O'Mullane Reference Summers and O'Mullane2011).
3.3. 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}\ \textrm {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. Reference Igochine, Dumbrajs, Zohm and Flaws2006; Papp et al. Reference Papp, Pokol, Por, Magyarkuti, Lazányi, Horváth, Igochine and Maraschek2011), in relation to an internal MHD instability – a physical mechanism beyond the modelling capabilities of the Go 
$+$ Code tool.
Figure 5. Time evolution of the electron energy spectrum (pitch-averaged distribution function) at the magnetic axis. The RE seed starts close to 
$p=5m_ec$ at 
$t=0$, and is then quickly accelerated to above 
$p=100m_ec$ within a few milliseconds. During the remainder of the runaway plateau, the initial seed sits around 
$p=100m_ec$ 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=0m_ec$ in this figure.
As shown in figure 5, the seed runaway population is quickly accelerated to a maximum energy during the CQ, which lasts for approximately 
${4}\ \textrm {ms}$. During this phase, a population of secondary runaways gradually builds up, overtaking the plasma current. The maximum energy varies across radii – from 
$p\approx 100m_ec$ in the core to only a few 
$m_ec$ 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 CQ, 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 electrons calculated with Go 
$+$ Code. The resulting synthetic camera image at 
$t={1.029}\ \textrm {s}$, just before the pattern transition, is presented in figure 6(c) (along with two preceding times in figure 6a,b). Comparing the synthetic image with the experimental image in figure 6(d), 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 REs, 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 TQ – is different. Deviations in plasma parameters such as the density and temperature could also have an effect on the avalanche gain during the CQ, 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.
Figure 6. Comparison of the synthetic synchrotron image produced with Soft, taking the distribution function calculated with Go 
$+$ Code at (a) 
$t={1.008}\ \textrm {s}$, (b) 
${1.018}\ \textrm {s}$ and (c) 
${1.029}\ \textrm {s}$ as input, with the (d) synchrotron image taken in ASDEX-U #35628, also at 
$t={1.029}\ \textrm {s}$. 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.
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 CQ and then remained at their peak energy. In figures 7(a) and 7(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 8(a) and 8(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.
Figure 7. Amount of synchrotron radiation observed from different parts of momentum space. Panels (a,b) show the contributions at 
$t={1.029}\ \textrm {ms}$ from two individual radii, indicating that the emission is dominated by the remnant hot-tail seed. The very sharp features running along almost constant pitch in panels (a,b) are physical, and are connected to the very bright edges usually seen in synchrotron images from mono-energetic and mono-pitch distribution functions. Panels (c,d) compare the radially integrated synchrotron radiation at 
$t={1.008}\ \textrm {ms}$ and 
$t={1.029}\ \textrm {ms}$, respectively.
Figure 8. Momentum space distribution functions (multiplied by the momentum-space Jacobian 
$p^2\sin \theta$) from the Go 
$+$ Code simulations at two select radii, chosen to correspond approximately to the particles contributing to figures 7(a) and 7(b). The remnant seed appears as a bump in the distribution function around (a) 
$p_\parallel = 80 m_ec$ and (b) 
$p_\parallel = 55m_ec$.
When integrating over all radii, a wider dominant region appears in momentum space, as shown in figures 7(c) and 7(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_\xi (\xi )\sim \exp (C\xi )$, with 
$C$ a time-dependent constant, and 
$\xi =p_\parallel /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
\begin{equation} C(t)\approx\frac{(p/m_ec)^2}{8n_{\textrm{Ar},20}(t-t_0)_\textrm{ms}}, \end{equation}
where 
$n_{\textrm {Ar},20}$ is the argon density, measured in units of 
$10^{20}\ \textrm {m}^{-3}$, and the times in the denominator are given in milliseconds. The emitted synchrotron power at a frequency 
$\omega$ can, furthermore, be approximated by the contribution from the strongest emitting particle of such a distribution which is proportional to
\begin{equation} \mathcal{P} = \exp\left[ -\left( \frac{\omega m_e}{3\sqrt{2}eB}\right)^{2/3}\frac{1}{\gamma^{\star 4/3}}\frac{1}{\left[\dfrac{1}{C(t_0)} + (t-t_0)\nu_D \right]^{1/3}}\right], \end{equation}
with 
$\gamma ^\star = \sqrt{1 + (p^\star/m_ec)^2}$, and 
${p^\star }$ the momentum reached by the hot tail seed. The free parameters in this expression are 
$\gamma ^\star$, 
$n_\textrm {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_\textrm {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.
Figure 9. Total detected synchrotron intensity as predicted by combined Go 
$+$ Code and Soft simulations (black, solid), and as recorded in the experiment (red, dashed). Since the experimental measurements are not absolutely calibrated, the curves have been rescaled to aid comparison of the slopes.
4. Backward modelling
The fluid-kinetic model described in § 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\approx 1.030\ \textrm {s}$. In this section, we instead analyse the synchrotron radiation images directly and extract information from the camera images using a regularized, 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 RE distribution building on the results of § 3.3, which allows us to better constrain runaway parameters.
4.1. Inversion procedure
We may capitalize on what we have learned from the fluid-kinetic simulations of § 3, in order to find constraints to regularize the inversion of the synchrotron images. We found that the synchrotron radiation pattern seen in the camera images was dominated by the remnant RE 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
\begin{equation} f_\textrm{seed}(r,p,\xi) = f_r(r)\exp\left[ -\left( \frac{p-{p^\star}}{\Delta p} \right)^2 \right]\exp(C\xi)\approx f_r(r)\delta\left( p-{p^\star} \right)\exp\left( C\xi \right), \end{equation}
where 
$f_r(r)$ is an arbitrary function describing the radial runaway density, and 
${p^\star }$, 
$\Delta 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^\star }$ and, since 
$C$ is mainly a function of 
${p^\star }$, 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 utilizes the capability of Soft to generate weight functions for a given tokamak/detector set-up, as described in § 2. The brightness of pixel 
$i$ in a synchrotron image 
$I_i$ is related to the distribution function 
$f(r,p,\xi ) = f_r(r)f_{p\xi }(p,\xi )$ through the weight function 
$G(r,p,\xi )$ as
\begin{equation} I_{i} = \int G_{i}(r,p,\xi) f(r,p,\xi)p^2\,\mathrm{d} r\,\mathrm{d} p\,\mathrm{d}\xi. \end{equation}By representing the image and the discretized radial runaway density profile as vectors, the discretized version of the equation system (4.2) can be formulated as
\begin{equation} I_{i} = \sum_k \tilde{G}_{ik}(r) f_r^k(r), \end{equation}where we have introduced the reduced weight function matrix
\begin{equation} \tilde{G}_{ik} = \Delta r_k \int G_{i}(r_k,p,\xi)f_{p\xi}(p,\xi)p^2\,\mathrm{d} p\,\mathrm{d}\xi, \end{equation}
with 
$\Delta r_k=|r_{k+1}-r_{k}|$. Given a momentum-space distribution function 
$f_{p\xi }(p,\xi )$, which can be characterized with the parameters 
${p^\star }$ and 
$C$ in (4.1), we thus seek to minimize the sum of squares of differences between pixels in the synthetic and experimental images 
$I_i$ and 
$I_i^\textrm {exp}$. Since the problem is still ill-posed, we regularize it using the Tikhonov method (Press et al. Reference Press, Teukolsky, Vetterling and Flannery2007), resulting in
\begin{align} f_r(r) &= \arg\min_{f_{r}}\left[ \left\lVert \sum_j\alpha\varGamma_{ij}\,f_r^j \right\rVert_2^2 + \sum_i\left\lVert I_i^\textrm{exp} - I_i \right\rVert_2^2 \right] \nonumber\\ &=\arg\min_{f_r} \left[ \left\lVert \sum_j \alpha\varGamma_{ij}\,f_r^j \right\rVert_2^2 + \sum_i\left\lVert I_i^\textrm{exp} - \sum_j \tilde{G}_{ij}\, f_r^j \right\rVert_2^2 \right], \end{align}
where the Tikhonov matrix 
$\varGamma _{ij}$ is taken to be the identity matrix and the Tikhonov parameter 
$\alpha$ is determined using the L-curve method (Hansen & O'Leary Reference Hansen and O'Leary1993). The ability to solve this problem using a linear least-squares method allows us to efficiently explore the space of possible combinations of 
$({p^\star }, C)$, and hence to use typical minimization methods to solve the full problem.
4.2. 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}\ \textrm {ms}$ and 
$t={1.030}\ \textrm {ms}$, corresponding to the video frames 3(b) and 3(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 
$\Delta I$ during the event, then it follows from a circuit equation for the plasma that the energy of the electrons changes according to
\begin{equation} \Delta W = \int ecE\,\mathrm{d} t = -\int \frac{ecL}{2{\rm \pi} R_0}\frac{\mathrm{d} I}{\mathrm{d} t}\,\mathrm{d} t = -\frac{ecL}{2{\rm \pi} R_0}\Delta I, \end{equation}
where 
$L$ is the self-inductance and 
$R_0$ is the major radius of the torus. Assuming 
$L\sim \mu _0 R_0$, with 
$\mu _0$ the vacuum permeability, the induced energy is found to increase by 
${60}\ \textrm {eV}$ for every ampere decrease in the plasma current. In our case, where the second current spike rises by 
$\Delta I\approx {5}\ \textrm {kA}$, the energy transferred to the electrons should be below 
${0.5}\ \textrm {MeV}$, which is small compared with typical runaway energies at 
$t={1.029}\ \textrm {ms}$ (
${\approx }{25}\ \textrm {MeV}$ at midradius, 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; 
$\Delta \theta \approx -(\Delta W_\| / W_\|) \theta$.
A more general argument against a change of the momentum-space distribution is that the synchrotron intensity is more sensitive to changes in the energy and pitch angle than to changes in the radial density profile. As discussed in § 2, the observed synchrotron intensity is exponentially sensitive to 
$p$ and 
$\xi$ in the short wavelength limit, whereas the radial density always appears as a multiplicative factor. Since the synchrotron intensity does not change significantly in the spot shape transition, the change to the momentum-space distribution should not be significant either. On the other hand, a parameter scan indicates that a significant change to the momentum-space distribution would be required for a visible spot shape transition. 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 3(a) and 3(b) simultaneously, assuming 
${p^\star }$ and 
$C$ to remain unchanged in the transition.
The sum of pixel differences squared, as a function of the fitting parameters 
${p^\star }$ and 
$C$, is shown in figure 10. In each point, the best radial density is constrained using (4.5). Optimal agreement is obtained with 
${p^\star }=57.5m_ec$ 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, 
${\theta ^\star }\approx 0.30\ \textrm {rad}$. This is in agreement with the Go 
$+$ Code and Soft simulations presented in § 3.3 which had 
${p^\star }\approx 55m_ec$, 
$C\approx 25$ (at 
$p={p^\star }$) and 
${\theta ^\star }\approx {0.36}\ \textrm {rad}$ at 
$t={1.029}\ \textrm {ms}$. Since 
$C$ is inversely proportional to the (relatively poorly diagnosed) argon density 
$n_\textrm {Ar}$, these are well within the uncertainties of the inversion and in the plasma parameters.
Figure 10. Sum of pixel differences squared for both of figures 11(a) and 11(b), given different combinations of 
${p^\star }$ and 
$C$. For each combination of 
${p^\star }$ 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^\star } = 57.5$ and 
$C = 45$. The white regions correspond to unreasonable combinations of the two parameters.
The radial density profiles obtained in the inversion for the two video frames are presented in figure 11(a). 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^\star }$ are considered, the inverted radial profile is exactly zero within the grey region 
$r/a\approx 0.1$, corresponding to the drift orbit shift for these particles. The blue and red shaded regions in figure 11(a) indicate the maximum deviation of density profiles corresponding to the optima of all combinations 
$({p^\star }, C)$ with normalized 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 evaluating the maximum deviation. The maximum deviations suggest that although the uncertainty in 
$({p^\star },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.
Figure 11. (a) Inverted radial density profiles for the video frames at 
$t={1.029}\ \textrm {s}$ (black) and 
$t=1.030\ \textrm {s}$ (red) for the best fitting values of (
${p^\star },C)$, and the corresponding inverted synthetic synchrotron radiation images at (b) 
$t={1.029}\ \textrm {s}$ and (c) 
${1.030}\ \textrm {s}$. The optimal values of 
${p^\star }$ and 
$C$ extracted from figure 10 and used to generate the images are 
${p^\star } = 57.5m_ec$ and 
$C\approx 45$. The blue and red shaded regions in panel (a) indicate the maximum variation of the radial profiles among all solutions with normalized likeness 
$\leqslant 2$ (corresponding to all points within the 2 contour of figure 10). The grey shaded region has the size of the drift orbit shift and contains no particles since Soft only counts particles in the outer midplane.
The radial density profile inversion reveals that the synchrotron pattern change is caused by a rapid redistribution of the RE density. As described in § 2, the shifted density profile causes the particles on the HFS to collectively emit relatively more radiation than those on the LFS, yielding a crescent pattern after the event. By integrating the radial density profiles in figure 11(a) 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 (Eriksson et al. Reference Eriksson, Helander, Andersson, Anderson and Lisak2004), 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 4(b).
5. Summary
We analysed visible light camera images of synchrotron emission from REs 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 spatiotemporal evolution of the runaway distribution function with a coupled kinetic-fluid code taking into account the evolution of the electric field self-consistently. 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 REs 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 RE 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.
Acknowledgements
The authors are grateful to R.A. Tinguely, K. Insulander Björk, O. Vallhagen, M. Dunne, V.G. Igochine and M. Maraschek for fruitful discussions, and to M. Gruber for operating the synchrotron camera diagnostics. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 and 2019–2020 under grant agreement no. 633053. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement no. 647121. The views and opinions expressed herein do not necessarily reflect those of the European Commission. The work was also supported by the Swedish Research Council (Dnr. 2018-03911) and the EUROfusion-Theory and Advanced Simulation Coordination (E-TASC). G.P. acknowledges the support of the National Research, Development and Innovation Office (NKFIH) Grant FK132134.
Editor Paolo Ricci thanks the referees for their advice in evaluating this article.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Accounting for guiding-centre drifts
To reduce the size of the phase-space from six to three dimensions, Soft utilizes a guiding-centre transformation which eliminates the orbit time, as well as the toroidal and gyrophase angles from the distribution function. However, the dependence on the particle's position and momentum components still remains in the radiation emission function, 
$\mathrm {d} P/\mathrm {d} \varOmega$; these must be computed within the guiding-centre framework. The RE 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(\epsilon )$ effect, where 
$\epsilon$ is the usual guiding-centre ordering parameter. An expression for the Larmor radius vector to first order was given by Hirvijoki et al. (Reference Hirvijoki, Decker, Brizard and Embréus2015), 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):
\begin{align} \boldsymbol{p} &= \boldsymbol{P} + p_\perp\,\mathrm{sign}(q)\hat{\boldsymbol{\perp}} + \epsilon\hat{\boldsymbol{b}}\left[ p_\parallel\rho(\hat{\boldsymbol{\rho}}\boldsymbol{\cdot}\boldsymbol{\kappa}) + \frac{m\mu}{q}(\boldsymbol{\mathsf{a}}_1 \boldsymbol{:} \boldsymbol{\nabla}\hat{\boldsymbol{b}} ) \right] \nonumber\\ &\quad + \epsilon\hat{\boldsymbol{\rho}}\left[ qB\frac{\rho^2}{2}( \hat{\boldsymbol{\perp}}\boldsymbol{\cdot}\boldsymbol{\nabla}\ln B) - \frac{\rho p_\parallel}{2}( \hat{\boldsymbol{b}}\boldsymbol{\cdot}\boldsymbol{\nabla}\ln B) - \rho p_\parallel (\boldsymbol{\mathsf{a}}_2 \boldsymbol{:} \boldsymbol{\nabla}\hat{\boldsymbol{b}} )\right] \nonumber\\ &\quad + \epsilon\hat{\boldsymbol{\perp}}\left[ qB\frac{\rho^2}{2}(\hat{\boldsymbol{\rho}}\boldsymbol{\cdot}\boldsymbol{\nabla} \ln B) - \frac{\rho p_\parallel}{2}(\boldsymbol{\mathsf{a}}_1 \boldsymbol{:} \boldsymbol{\nabla}\hat{\boldsymbol{b}}) \right] + O(\epsilon^2), \end{align}with the guiding-centre momentum vector
\begin{equation} \boldsymbol{P} = p_\parallel\hat{\boldsymbol{b}} + \epsilon\frac{\hat{\boldsymbol{b}}}{qB^\star_\parallel}\times (m\mu \boldsymbol{\nabla} B + p_\parallel^2\boldsymbol{\kappa}), \end{equation}
where 
$p_\parallel$ and 
$p_\perp$ are the particle momenta parallel and perpendicular to the magnetic field, respectively, 
$B$ is the magnetic field strength, 
$\hat {\boldsymbol {b}}$ and 
$\hat {\boldsymbol {\rho }}$ are mutually perpendicular unit vectors in the direction of the magnetic field and (lowest-order) Larmor radius vector, respectively, 
$\hat {\boldsymbol {\perp }} = \hat {\boldsymbol {\rho }}\times \hat {\boldsymbol {b}}$, 
$\rho = p_\perp /|q|B$ is the Larmor radius, 
$\mu =p_\perp ^2/2mB$ is the magnetic moment, 
$\boldsymbol {\kappa } = (\hat {\boldsymbol {b}}\boldsymbol {\cdot }\boldsymbol {\nabla })\hat {\boldsymbol {b}}$ is the inverse curvature vector, 
$B^\star _\parallel = B+p_\parallel \hat {\boldsymbol {b}}\boldsymbol {\cdot }(\boldsymbol {\nabla }\times \hat {\boldsymbol {b}})/q$, and the dyads 
$\boldsymbol{\mathsf{a}}_1$ and 
$\boldsymbol{\mathsf{a}}_2$ are defined as
\begin{equation} \boldsymbol{\mathsf{a}}_1 = -\tfrac{1}{2}(\hat{\boldsymbol{\rho}}\hat{\boldsymbol{\perp}} + \hat{\boldsymbol{\perp}}\hat{\boldsymbol{\rho}}), \quad \boldsymbol{\mathsf{a}}_2 = -\tfrac{1}{4}( \hat{\boldsymbol{\perp}}\hat{\boldsymbol{\perp}} - \hat{\boldsymbol{\rho}}\hat{\boldsymbol{\rho}}). \end{equation}
To lowest order, the momentum vector (A 1) describes circular gyro motion around the magnetic field line. However, when including 
$O(\epsilon )$ 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 field 
$\hat {\boldsymbol {b}}$ to the guiding-centre direction of motion 
$\boldsymbol {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 
$\theta =\arctan (p_\perp /p_\parallel )$, around the guiding-centre momentum vector 
$\boldsymbol {P}$.
A.2. Implementation and validation
The conclusions from the analysis of the electron momentum vector suggest that the cone model, implemented previously Soft (Hoppe et al. Reference Hoppe, Embréus, Tinguely, Granetz, Stahl and Fülöp2018b), can be modified to take guiding-centre drifts into account by including the first-order drift terms in the guiding-centre momentum vector 
$\boldsymbol {P}$, and consider the radiation to be emitted in an infinitesimally thin circular cone around 
$\boldsymbol {P}$. The method is, however, not exact since the 
$O(\epsilon )$ 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 full-orbit/guiding-centre model that can indicate when the cone model assumptions are violated.
The synchrotron diagnostic tool Soft computes the radiation observed from an ensemble of particles with distribution function 
$f_\textrm {gc}(r,\boldsymbol {p})$ by evaluating the integral
\begin{equation} P = \int\varTheta\left(\frac{\boldsymbol{r_{\text{cp}}}}{r_{\text{cp}}}\right) \frac{\hat{\boldsymbol{n}}\boldsymbol{\cdot}\boldsymbol{r_{\text{cp}}}}{r_{\text{cp}}^3} \frac{\mathrm{d} P(\boldsymbol{x},\boldsymbol{p},\boldsymbol{r_{\text{cp}}})}{\mathrm{d}\varOmega} f_\textrm{gc}(r,p,\mu)\,J\,\mathrm{d} r\, \mathrm{d}\tau \,\mathrm{d}\phi\,\mathrm{d}\boldsymbol{p}^{(0)}\,\mathrm{d} A. \end{equation}
Here, 
$\boldsymbol {r_{\text {cp}}} = \boldsymbol {x}_0-\boldsymbol {x}$ is the vector extending from the detector at 
$\boldsymbol {x}_0$ to the particle at 
$\boldsymbol {x}=\boldsymbol {x}(r,\tau ,\phi )$, 
$\boldsymbol {p}^{(0)}$ is the particle momentum at the outer midplane, 
$\hat {\boldsymbol {n}}$ is the detector surface normal, 
$\varTheta$ is a step function that is one whenever 
$\boldsymbol {r_{\text {cp}}}/r_{\text {cp}}$ lies in the detector field-of-view, 
$\mathrm {d} P/\mathrm {d}\varOmega$ is the angular distribution of radiation emitted by each particle, and 
$J$ is the Jacobian for the full coordinate transformation (including the guiding-centre 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 
$\tau$; and the toroidal angle 
$\phi$. The momentum 
$\boldsymbol {p}$ is therefore a function of both 
$\boldsymbol {x}$ and 
$\boldsymbol {p}^{(0)}$.
Note that, although we use a guiding-centre transformation, the integrand (A 4) is still dependent on the particle position and momentum 
$\boldsymbol {x}$ and 
$\boldsymbol {p}$. This means that we have to use (A 1) and the corresponding expression for the Larmor radius vector 
$\boldsymbol {\rho }$ to consistently evaluate (A 4) to 
$O(\epsilon )$. 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(\epsilon )$ terms is to rotate the circular cone of radiation around the guiding-centre velocity vector, thus we simply modify the cone model according to
\begin{equation} \frac{\mathrm{d} P}{\mathrm{d}\varOmega}\sim\delta(\boldsymbol{r_{\text{cp}}}\boldsymbol{\cdot}\hat{\boldsymbol{b}} - r_{\text{cp}}\cos\theta) \longrightarrow \frac{\mathrm{d} P}{\mathrm{d}\varOmega}\sim\delta( \boldsymbol{r_{\text{cp}}}\boldsymbol{\cdot}\boldsymbol{P} - r_{\text{cp}} P\cos\theta ), \end{equation}
where 
$\hat {\boldsymbol {b}}$ denotes the local magnetic field unit vector and 
$\cos \theta = \boldsymbol {p}\boldsymbol {\cdot }\hat {\boldsymbol {b}}$.
To verify that the approximation (A 5) gives sufficiently accurate results, 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 (A 1), we calculate a set of particle orbits, from which 
$\boldsymbol {r_{\text {cp}}}$ and 
$\boldsymbol {p}$ are obtained exactly. While this approach may at first seem to be more accurate, we point out that (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 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 
$\theta = {0.2}\ \textrm {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 12(b). 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.
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 (
$\theta ={0.2}\ \textrm {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 LFS, corresponding to the drift orbit shift direction.
Appendix B. Estimation of ${p^\star }$ and ${\theta ^\star }$ 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^\star }$ and 
${\theta ^\star }$, 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 electron moving in a uniform magnetic field (Bekefi Reference Bekefi1966):
\begin{equation} \left.\begin{gathered} \frac{\partial P}{\partial \omega} = \frac{m_e c r_0}{\sqrt{3}{\rm \pi}} \frac{\omega}{\gamma^2}\int_{\omega/\omega_c}^\infty K_{5/3}(x)\,\mathrm{d} x, \\ \omega_c = \frac{3}{2}\frac{eB}{m_e} \frac{\gamma^2}{\gamma_\parallel} = \frac{3}{2}\frac{eB}{m_e}\gamma\sqrt{1+(1-\xi^2)(p/m_e c)^2}, \\ \gamma_\parallel = \frac{1}{\sqrt{1-\xi^2(v/c)^2}} = \frac{\gamma}{\sqrt{1+(1-\xi^2)(p/m_e c)^2}}, \\ \omega = 2{\rm \pi} c/\lambda. \end{gathered}\right\} \end{equation}
When 
$p_{\perp } \gg m_e c$, the critical frequency – where the emission is the strongest – is approximately given by 
$\omega _c \approx 3eB\gamma ^2\theta /(2m_e)$. Considering an observed wavelength of 
$\lambda \approx 700\ \textrm {nm}$, a representative pitch angle 
$\theta \approx 0.2$ and energy 
$\gamma \approx 50$, we find that 
$\omega /\omega _c \approx 8$. Then, the emission is well described by the formula evaluated at the 
$\omega /\omega _c \gg 1$ limit
\begin{equation} \frac{\partial P}{\partial \omega} \approx \frac{m_e c r_0}{\sqrt{3}{\rm \pi}}\frac{\sqrt{\omega_c\omega}}{\gamma^2}\textrm{e}^{-\omega/\omega_c}. \end{equation}
The total emission from a distribution of runaways, which we assume can be approximated as 
$f_{\mathrm {RE}} \approx F(t,p)\exp (-C(t,p)\theta ^2/2)$, is given by the phase-space integral of
\begin{equation} \mathcal{P} = \exp\left( -\frac{2\omega m_e}{3eB}\frac{1}{\gamma^2\theta} - \frac{C}{2}\theta^2\right), \end{equation}
up to a factor depending on plasma and runaway parameters at most polynomially (in contrast to the dominant, exponential dependence captured by 
$\mathcal {P}$). At a given 
${p^\star }$ (the precise value of which depends on the energy distribution of runaways), this function is maximized by the pitch angle 
${\theta ^\star }$, and takes the maximum value 
$\mathcal {P}({\theta ^\star })$, given by
\begin{equation} \left.\begin{gathered} {\theta^\star} = \left(\frac{2\omega m_e}{3eB C\gamma^{\star2}}\right)^{1/3}, \\ \mathcal{P}({\theta^\star}) = \exp\left[ -\left( \frac{\omega m_e}{3\sqrt{2}eB}\right)^{2/3} \frac{C^{1/3}}{\gamma^{\star 4/3}}\right]. \end{gathered}\right\}\end{equation}B.1. Determination of $C$ from kinetic dynamics
At the end of the CQ, 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 Reference Breizman2014). 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
\begin{equation} \frac{\partial f}{\partial t} = \frac{\partial}{\partial \xi}\left[(1-\xi^2)\left(-\frac{E}{p}f + \frac{\nu_D(p)}{2}\frac{\partial f}{\partial \xi}\right)\right], \end{equation}
which is similar to the equation studied by Aleynikov & Breizman (Reference Aleynikov and Breizman2015) to determine the equilibrium distribution. Taking integral moments in 
$\xi$ of this equation yields
\begin{align} \langle 1-\xi \rangle &= \frac{\int_{-1}^1\mathrm{d}\xi (1-\xi) f }{\int_{-1}^1\mathrm{d}\xi f }, \nonumber\\ \frac{\partial \langle 1-\xi \rangle}{\partial t} &= \frac{\int_{-1}^1 \mathrm{d} \xi \left[ -(1-\xi^2) \dfrac{E}{p}f + \nu_D \xi f\right]}{ \int_{-1}^1\mathrm{d}\xi f} \nonumber\\ &\approx - \frac{2E}{p}\langle 1-\xi\rangle + \nu_D, \end{align}
where the last relation is valid for small pitch angles 
$1-\xi \ll 1$. Since the initial pitch angle distribution is highly anisotropic – having been generated by an acceleration in a strong electric field – the first term will initially be negligible, 
$\nu _D \gg 2E\langle 1-\xi \rangle /p$. For a distribution of the form 
$f\propto \exp (-C\theta ^2/2)$, one finds 
$\langle 1 - \xi \rangle = 1/C$, which then yields the time evolution
\begin{equation} C(t) = \frac{1}{\dfrac{1}{C(t_0)} + (t-t_0)\nu_D}\approx \frac{p^{\star 2}}{8n_{\mathrm{Ar},20}(t-t_0)_{\mathrm{ms}}}. \end{equation}In a cold post-disruption plasma where the pitch-angle scattering of electrons is dominated by collisions with argon impurities, the diffusion rate in pitch-angle is given by Hesslow et al. (Reference Hesslow, Embréus, Hoppe, DuBois, Papp, Rahm and Fülöp2018)
\begin{align} \nu_D &\approx 4{\rm \pi} n_\mathrm{Ar} c r_0^2\frac{\gamma}{p^3} \ln 90p \nonumber\\ &\approx \frac{8n_{\mathrm{Ar},20}}{p^{\star2}}\,\times (1 \ \text{ms}^{-1}), \end{align}
where the argon density is expressed in units of 
$10^{20}\ \textrm {m}^{-3}$ in the last expression.
Combining this result with (B 4), we finally obtain the following expression for the time evolution of the total emitted synchrotron power:
\begin{equation} \mathcal{P} = \exp\left[ -\left( \frac{\omega m_e}{3\sqrt{2}eB}\right)^{2/3}\frac{1}{\gamma^{\star 4/3}}\frac{1}{\left[\dfrac{1}{C(t_0)} + (t-t_0)\nu_D \right]^{1/3}}\right]. \end{equation}
Open access
