Sub-microsecond temporal evolution of edge density during edge localized modes in KSTAR tokamak plasmas inferred from ion cyclotron emission

Ion cyclotron emission (ICE) is detected during edge localised modes (ELMs) in the KSTAR tokamak at harmonics of the proton cyclotron frequency in the outer plasma edge. The emission typically chirps downward (occasionally upward) during ELM crashes, and is driven by confined 3MeV fusion-born protons that have large drift excursions from the plasma core. We exploit fully kinetic simulations at multiple plasma densities to match the time-evolving features of the chirping ICE. This yields a unique, very high time resolution diagnostic of the collapsing edge pedestal density.


Introduction
Understanding the physics of edge localised modes (ELMs) [1][2][3][4] in magnetically confined fusion (MCF) plasmas is crucial for the design of future fusion power plants. The same is true of the physics of the energetic ions born at MeV energies [5][6][7] from fusion reactions between fuel ions in a multi keV thermal plasma. The crash of an ELM involves impulsive relaxation of the edge magnetic field, releasing energy and particles from the plasma at levels which may not be compatible with sustained operation of the next step fusion experiment, ITER [8,9].
In recent years the KSTAR tokamak has performed a series of experimental campaigns aimed at advancing the current understanding of ELMs A thorough account of ELM dynamics in KSTAR can be found in [10,11], in which three distinct stages of ELM filament evolution in KSTAR are detailed: (1) the initial growth of the ELM filaments near the last closed flux surface (LCFS), which grow to a saturated state in ∼300 µs; (2) the interim quasi-steady state of saturated filaments typically persisting for ∼100 µs; and (3) the collapse of the pedestal through multiple toroidally and poloidally localised filament bursts occurring over a time-scale of ∼100 µs. This third and final stage of the ELM dynamics is referred to henceforth as the 'ELM crash'. Towards the end of the ELM crash, one often observes bursts of ion cyclotron emission (ICE) whose spectral structure exhibits frequency chirping over time scales 10 µs. This phenomenology corresponds to 'stage (4)' in the characterisation of evolving radio frequency (RF) burst spectra in KSTAR in [11]. This unexpected conjunction between (i) the ELM physics associated with individual filament bursts during this terminal stage of the overall ELM cycle and (ii) the fusion-born ion physics responsible for ICE is the subject of this letter.
ICE comprises suprathermal radiation in the ion cyclotron frequency range, whose spectrum peaks at successive local cyclotron harmonics of the emitting ion population. ICE has previously been observed in all large toroidal MCF plasmas [12]. ICE is caused by a collective instability, which in its linear phase corresponds to the magnetoacoustic cyclotron instability (MCI) [13][14][15][16][17][18][19][20][21][22][23]. The MCI requires a spatially localised population inversion in velocity space for the energetic ions which drive it. This inversion can arise from, for example, large orbital drift excursions.
In this letter we perform, for the first time, a detailed quantitative comparison between ICE observed during KSTAR ELM crashes and fully nonlinear direct numerical simulations of kinetic particles together with electric and magnetic fields, evolving self-consistently under the Maxwell-Lorentz system of equations. We find good quantitative agreement between the simulated and observed spectra, to the extent that the simulations can be used to infer fast (∼μs) time scale dynamics of the local electron number density in the emitting region. Figure 1 shows an example of downward ICE chirping from a KSTAR plasma that has with toroidal magnetic field at the magnetic axis B 0 1.99 T, plasma current I p 600 kA, and total stored energy W ∼ 380 kJ. The upper panel displays the measured RF burst spectra as a function of time while the lower panel displays the Fourier transform of this signal. In the lower panel of figure 1, and all subsequent figures of this type, horizontal white lines denote successive proton cyclotron harmonics evaluated for the magnetic field strength at the outer mid-plane edge of the plasma, while t = 0 refers to the centre of a 200 µs segment of RF data. This data is obtained when the RF signal amplitude exceeds a threshold voltage during KSTAR pulse operation, with the acquisition times corresponding roughly to a spike in the D α signal [11]. For further details of the fast RF spectrometer system used on KSTAR see for example figure 1 of [11].

ELMs and ICE in KSTAR
The frequency chirping observed in KSTAR is often, as in figure 1, in discrete steps coinciding with the local proton cyclotron frequency. The only energetic protons in KSTAR plasmas are those produced in deuteron-deuteron (D-D) fusion reactions, hence it is likely that ICE at spectral peaks separated by proton cyclotron harmonics is driven by fusionborn protons. If this ICE is driven by the MCI of confined fusion-born protons with spatially localised population inversion, it is necessary to identify a candidate population; and we do so in section 3. The KSTAR experiment is not built on a scale sufficiently large to confine the majority of the energetic ions that are born in fusion reactions within pure deuterium plasmas: 3.0 MeV protons, 0.82 MeV He-3 nuclei, and 1.0 MeV tritons. To drive ICE in KSTAR via the MCI, there must nonetheless exist a collectively unstable subset of fusion-born protons whose orbits are confined within the plasma and traverse the excitation region in the outer midplane edge. The distribution in velocity space of this population must differ qualitatively from that driving ICE in other MCF plasmas such as JET and TFTR [24,25], because the fusion-born ions on the marginally trapped orbits that give rise to ICE in these two larger experiments are promptly lost from KSTAR during their first drift excursion.

Single particle orbits
We first calculate orbits of 3.0 MeV protons in equilibrium magnetic fields resembling those of a KSTAR plasma with major radius R 0 = 1.8 m, toroidal magnetic field B 0 = 2.27 T, and plasma current I p = 611 kA. Our orbit calculations show that almost all centrally-born fusion protons are lost promptly from the plasma on their first drift orbit. However, a small fraction of these protons are born onto deeply passing orbits which remain confined at the plasma edge to which they are carried by their radial drift excursion. Figure 2 shows example orbits of 3.0 MeV protons born in the midplane at initial major radii R(0) equal to (a) 1.85 m and (b) 1.90 m, and with initial velocity vectors slightly offset from the co-current toroidal direction. These orbits pass through the outer midplane plasma edge, and could give rise there to a localised population in velocity space. We show below that such a distribution can radiate collectively through the MCI and is thus capable of driving ICE at proton cyclotron harmonics characteristic of the outer midplane. We conjecture that when the ELM crash starts, confinement of all energetic ions at the edge is lost; the edge is then rapidly re-populated on a drift orbit timescale by this subset of the energetic fusion-born protons. This leads to the sharp local population inversion that drives ICE. We note that an instance of positive correlation between ELMs and ICE was also seen on JET, see figure 9 of [21].
The spatial location from which ICE is omitted is inferred from the observed frequency interval between successive spectral peak features, which corresponds to the local proton cyclotron frequency. This fixed proton cyclotron frequency implies that the radial location of the emitting region does not change. The full spectral evolution takes place on a microsecond timescale which is very short compared to any characteristic timescale of evolution of the energetic ion distribution in velocity space; thus we do not attempt to model any changes in the energetic ion distribution function.
The KSTAR magnetic field strength uniquely defines a radial location for the protons driving the ICE, to high accuracy, at the outer mid-plane.

Self consistent simulations of the MCI
To simulate the excitation of ICE by fusion protons in KSTAR we use the EPOCH particle-in-cell (PIC) code [26], with one spatial and three velocity dimensions. This method solves the full set of Maxwell's equations and the relativistic Lorentz force law, to capture the full particle dynamics with respect to all three velocity co-ordinates and a spatial co-ordinate, while self-consistently evolving all three vector components of the electric and magnetic fields. The simulations are set up in planar geometry, hence they do not incorporate realistic toroidal geometry and the associated compressional Alfvén eigenmode structure [27][28][29][30][31]. Nevertheless this approach has been shown to be successful in capturing most of the underlying physics [32][33][34], and aligns with the original slab-geometry analytical theory of the MCI [13][14][15][16][17][18][19][20][21][22][23], against which it has been benchmarked, see figure 1 of [32] and figure 1 of [33]. The simulations capture the rapid self-consistent time-evolution of the ion velocity distributions under the MCI, but do not attempt also to capture the much slower (by orders of magnitude) effects of collisional slowing-down.
The deeply passing subset of the fusion-born proton population in the KSTAR plasma edge discussed above has a speed perpendicular to the magnetic field (v ⊥0 ) much smaller than the magnitude of their speed at birth. It is therefore justifiable to represent this by the approximate model of a deltafunction ring distribution. Previous analytical [14][15][16][17] and first principle, fully self-consistent simulation studies [32][33][34] of the MCI showed that strong driving of the MCI is enabled if the perpendicular component of the energetic ion velocity is at least comparable to the local Alfvén velocity. In all our simulations, v ⊥0 corresponds to an energy 150 keV ≃ 5% of the birth energy consistent with the orbit calculations shown in figure 2. This is comparable to the local Alfvén speed c A , hence high enough to drive the MCI. The large parallel velocities v of these passing fusion protons are not represented in the simulations since the parallel dynamics of the fast ions play no role in perpendicular wave propagation. The ratio of proton number density to deuteron number density is 10 −3 and the total simulation duration is 10 proton cyclotron periods, carrying the MCI into its saturated nonlinear regime. The spatial (x) direction is orthogonal to the uniform magnetic field B = B cyc e z , so that the propagation direction of waves excited in the system is perpendicular to B. B cyc denotes the edge magnetic field inferred from the spacing between successive proton cyclotron harmonics observed in the experimental spectrograms. The bulk plasma comprises electrons and deuterons with initial temperature 1 keV.
The distribution of collectively radiated energy between different cyclotron harmonics depends on the character of the instability (the MCI) driving the ICE. A key dimensionless parameter in the analytical theory of the MCI is the ratio of energetic ion perpendicular velocity to the Alfvén speed. Changing the density changes the Alfvén speed, which in turn changes the numerical value of this ratio. This ratio determines the threshold for instability which can be inferred from both linear theory and fully nonlinear self-consistent PIC simulations. These nonlinear simulation results have many good points of contact [32][33][34] with the linear instability theory of the MCI, see in particular, figures 1-3 in [34]. This means the character of the MCI depends strongly on the local density so that if the local density changes, the set of cyclotron harmonics which are excited by the MCI also changes. Thus the spectral character of the ICE at any time reflects the densitydependence of the MCI.
To enable comparison with observations of ICE chirping during the ELM crash in KSTAR plasmas, we run multiple simulations to determine the density-dependence of the MCI in the fully nonlinear saturated regime. These simulations are carried out with initial electron number densities in the range 0.2 × 10 19 m −3 to 2.5 × 10 19 m −3 . This range reflects Thomson scattering measurements in the edge pedestal (see the last paragraph of section 3 in [11]). Figures 3 and 4 show the spectra of saturated ICE intensity obtained from EPOCH simulations with successive values of electron number density (lower panels), along with the corresponding experimentally-measured spectrograms for downward chirping ICE during ELM crashes in plasmas with B 0 = 1.7 T and B 0 = 1.99 T (upper panels). The lower left and lower right panels correspond to the high and low frequency ranges respectively. Electron number density decreases from left to right, and each vertical strip corresponds to an independent MCI simulation at the density shown. The arrows labelled (a)-(d) denote a mapping between experimentally observed and simulation proto n cyclotron harmonics. The boundary of evanescence in the cold plasma limit, defined by the lower hybrid frequency ω LH , is clearly visible. As the electron number density, and subsequently ω LH decrease, the number of modes available for excitation also reduces. For n e 1.1 × 10 19 m −3 , Figure 2. Poloidal projection of 3.0 MeV fusion proton orbits in the model KSTAR equilibrium with initial velocity vectors slightly offset from the co-current toroidal direction. The blue boxes in both the left and right panels designate the emitting region for the ICE and are centred on the radial location that corresponds to the proton cyclotron frequency whose harmonics are excited in the chirping ICE spectrum. The horizontal extent of these boxes corresponds to the local Larmor radius of cyclotron gyration of an emitting proton. the spectrum is dominated by a single cyclotron harmonic, with the harmonic number falling monotonically as n e decreases. It is apparent that for both values of central magnetic field B 0 shown in figures 3 and 4, the trend in ICE spectral power as a function of frequency and n e is the same.

ICE intensity as a function of frequency and electron number density
The variation of simulated ICE intensity with electron number density resembles the experimentally-observed variation of ICE intensity with time. We address this in greater detail in section 5.2, meanwhile we note that the striking agreement between the highest and lowest intensity regions in the upper and lower panels of both figures 3 and 4 reinforces the choice of simulation electron number density values which reflect Thompson scattering measurements [11]. There is a 'missing harmonic' in the lower panels of figures 3 and 4 at f 323 MHz and f 375 MHz respectively. This could be due to periodic boundary conditions, as well as the limitations of a 1D3V model.
For each of the downward ICE pulse-simulation pairs shown in figures 3 and 4 we are able to construct a sequence of mappings from time in the upper panels to density in the lower panels. This is discussed in the next section.

Mapping between computational results and KSTAR observations
Our simulations suggest that the downward chirping of ICE observed during the KSTAR ELM crashes detailed above is likely to be a direct consequence of a rapid local decline in electron number density associated with an ELM filament. An important corollary of this result is that measurements of the ICE spectra can be used to infer the locally declining electron number density as plasma is transported out of the ICEemitting region by one of the multiple filament bursts that are associated with the ELM crash.
In figure 4 (upper panel), a given ICE spectral feature at a particular cyclotron harmonic in KSTAR can be seen to arise at a time t start and end at a time t finish . For this downward chirping case, a neighbouring spectral feature at a lower cyclotron harmonic arises at a slightly later t start and persists until a slightly later t finish . In the lower panels of figure 4, a succession of simulated ICE spectra are shown. Each vertical strip represents an ICE simulation run at a different, neighbouring, density. These simulated spectra are arranged in sequence with number density declining from left to right. A particular cyclotron harmonic spectral feature (e.g. in the lower right of figure 4) is excited by simulations in a relatively narrow range of densities between n upper and n lower . Remarkably, the pattern of spectral features in the sequence of simulated spectra, depending on number density n, has much in common with the pattern of spectral features in the experimental ICE spectra, depending on time. Hence we can construct a simple mapping from time to density by exploiting these similarities.  Essentially, for a given spectral feature, we identify n upper with the density at t start , and n lower with the density at t finish .
The foregoing describes the overall approach. It yields, we argue, sub-microsecond time resolution for the ELM filament burst. This is shown in figure 5, the left panel corresponding to figure 3, and the right panel corresponding to figure 4. Our sequence of simulations at different electron number densities, together with the chirping ICE observations, reflect a phase of the ELM crash during which individual, toroidally and poloidally localised filament bursts lasting < 10μs (figure 4, [11]) are known to be present. In greater detail, the procedure for obtaining figure 5 (left) from figure 3 is as follows: This yields an available time resolution of the density decrease of order ∼0.1 µs. Figure 5 shows that the density decreases approximately linearly on a timescale of about 2.5 μs.

Upward chirping
Upward chirping is occasionally observed during ELM crashes in KSTAR. For logical consistency with the above, we should attribute this to a rare locally rising density due to the rapid motion of an ELM filament transporting additional plasma into the ICE-emitting region during the ELM crash. We test this hypothesis as follows. The time evolution of the corresponding electron number density for the case of upward chirping shown in figure 6 (left panel) has been inferred from the saturated MCI simulations at different electron number densities ( figure 6 right panel), and is shown in figure 7.  In figure 6, the separation between successive spectral peaks in the ICE is f cp ∼ 28 MHz, which corresponds to a local magnetic field strength of 1.84 T and thus implies that the ICE source is localised to the outer midplane of this KSTAR plasma. The absolute value of each spectral peak frequency is at integer multiples of 28 MHz plus, systematically, 10 MHz. We conjecture that this systematic frequency shift may be a Doppler shift. If, as we suggest, the source plasma for the upward chirping ICE in figure 6 is a moving ELM filament, this Doppler shift would be a consequence of filament motion. Testing and exploitation of this particular aspect of the phenomenology would require information of the k spectrum of the ICE, which is not at present available.

Conclusions
We have shown that harmonic ICE with spacing equal to f cp in KSTAR deuterium plasmas can be driven by a small subset of the fusion-born proton population originating in the core of the plasma. The drift orbits of these protons have large radial excursions to the outer midplane edge. We have compared the nonlinearly saturated field spectra obtained from multiple MCI simulations at different plasma densities with exper imentally observed time evolving ICE spectra. These results suggest that downward chirping of ICE occurs when the emitting fusionborn protons ions are embedded in a local plasma with rapidly falling density. These results also show that the much rarer upward chirping of ICE occurs when the local density is rapidly rising. By combining different simulation spectra with the chirping ICE observed during KSTAR ELM crashes, we have obtained sub-microsecond time resolution measurements of the evolving electron number density in the emitting region. We conjecture that the time-evolution of the electron number density on such rapid timescales in this region results from the motion of an ELM filament; for a detailed account of the experimental set-up and the detection of chirping ICE in conjunction with ELM filaments, see [11]. This letter opens unexpected lines of research linking fusion-born ion physics with ELM crashes in medium size tokamaks such as KSTAR, which generated the ICE observations studied here. It also provides a new, very high time resolution, method of quantifying the experimental phenomenology. The discovery that one can investigate collective fusion-born ion physics in medium size tokamaks, the connection to ELM phenomenology, and the high time resolution, are all new. Because of the relatively wide access that exists to medium-size tokamaks worldwide, combined with the centrality of fusion reactivity and ELMs to tokamak plasma physics, wider exploitation of these results should be possible. Combined with the passive, non-invasive character of ICE measurements, the results of this letter suggest an attractive way forward for future energetic ion measurements in ITER [8,9].