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

(2.1)\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.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:

(3.1)\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}$,

(3.2)\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,

(3.3)\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):

(3.4)\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:

(3.5)\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.

1. (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$.

2. (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.

3. (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)\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}
   where $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

(3.7)\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

(3.8)\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.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

(4.1)\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

(4.2)\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

(4.3)\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

(4.4)\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

(4.5)\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

(4.6)\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).