Interpretation of suprathermal emission at deuteron cyclotron harmonics from deuterium plasmas heated by neutral beam injection in the KSTAR tokamak

Intense bursts of suprathermal radiation, with spectral peaks at frequencies corresponding to the deuteron cyclotron frequency in the outer midplane edge region, are often detected from deuterium plasmas in the KSTAR tokamak that are heated by tangential neutral beam injection (NBI) of deuterons at 100 keV. Identifying the physical process by which this deuterium ion cyclotron emission (ICE) is generated, typically during the crash of edge localised modes, assists the understanding of collective energetic ion behaviour in tokamak plasmas. In the context of KSTAR deuterium plasmas, it is also important to distinguish deuterium ICE from the ICE at cyclotron harmonics of fusion-born protons examined by Chapman et al (2017 Nucl. Fusion 57 124004; 2018 Nucl. Fusion 58 096027). We use particle orbit studies in KSTAR-relevant magnetic field geometry, combined with a linear analytical treatment of the magnetoacoustic cyclotron instability (MCI), to identify the sub-population of freshly ionised NBI deuterons that is likely to excite deuterium ICE. These deuterons are then represented as an energetic minority, together with the majority thermal deuteron population and electrons, in first principles kinetic particle-in-cell (PIC) computational studies. By solving the Maxwell–Lorentz equations directly for hundreds of millions of interacting particles with resolved gyro-orbits, together with the self-consistent electric and magnetic fields, the PIC approach enables us to study the collective relaxation of the energetic deuterons through the linear phase and deep into the saturated regime. The Fourier transform of the excited fields displays strong spectral peaks at multiple successive deuteron cyclotron harmonics, mapping well to the observed KSTAR deuterium ICE spectra. This outcome, combined with the time-evolution of the energy densities of the different particle populations and electric and magnetic field components seen in the PIC computations, supports our identification of the driving sub-population of NBI deuterons, and the hypothesis that its relaxation through the MCI generates the observed deuterium ICE signal. We conclude that the physical origin of this signal in KSTAR is indeed distinct from that of KSTAR proton ICE, and is in the same category as the NBI-driven ICE seen notably in the TFTR tokamak and LHD heliotron–stellarator plasmas. ICE has been proposed as a potential passive diagnostic of energetic particle populations in ITER plasmas; this is assisted by clarifying and extending the physics basis of ICE in contemporary magnetically confined plasmas.

NBI-driven ICE seen notably in the TFTR tokamak and LHD heliotron-stellarator plasmas. ICE has been proposed as a potential passive diagnostic of energetic particle populations in ITER plasmas; this is assisted by clarifying and extending the physics basis of ICE in contemporary magnetically confined plasmas.
Keywords: ion cyclotron emission, magnetoacoustic cyclotron instability, ELM, tokamak, particle in cell, KSTAR, neutral beam injection (Some figures may appear in colour only in the online journal)

Introduction
In magnetically confined fusion (MCF) plasmas, the measured power spectrum of ion cyclotron emission (ICE) typically exhibits multiple narrow peaks whose frequencies are identified with sequential cyclotron harmonics of ions in the emitting region. The intensity of ICE is typically orders of magnitude greater than that of black-body radiation from thermal ions. ICE is therefore interpreted in terms of the collective relaxation of a non-Maxwellian energetic ion population which is spatially localised to the emitting region of the plasma. ICE from the outer midplane edge region was the first collective radiative instability driven by confined fusion-born ions that was observed from deuterium-tritium (D-T) plasmas in JET [1][2][3] and TFTR [4,5]. The measured intensity of the ICE signal scaled linearly with the measured neutron flux and with the inferred alpha-particle concentration, and hence with fusion reactivity [6]. ICE from the outer edge of the vessel has been detected from the largest contemporary toroidal plasmas including DIII-D [7], ASDEX-Upgrade [8], JT-60U [9], KSTAR [10,11], and LHD [12]. Recently, ICE due to both fusion-born and neutral beam injected (NBI) ions in the core of the tokamak has been detected on ASDEX-Upgrade [13,14], DIII-D [15,16], and TUMAN-3M [17]. The use of ICE as a diagnostic for confined energetic ion populations in ITER is under active consideration [18], and also for future D-T plasmas in JET. ICE is an ubiquitous plasma phenomenon, and has also been detected in solar-terrestrial plasmas [19][20][21][22], including the Earth's Van Allen belts [23], and is possibly present downstream of supernova remnant shocks [24].
The excitation mechanism for ICE is the magnetoacoustic cyclotron instability (MCI) [18,[25][26][27][28][29][30][31][32][33][34][35][36][37][38][39][40]. This arises when a minority energetic ion population enters into cyclotron resonance with a fast Alfvén wave propagating nearly perpendicular to the local magnetic field. This resonance can either be wave-wave [26][27][28], between cyclotron harmonic waves supported by the minority ions and fast Alfvén waves supported by the bulk plasma, or wave-particle [30,31], at Dopplershifted cyclotron resonance between the minority ions and the Alfvén wave. A necessary condition for the MCI to occur is that the energetic ions which drive it have a velocity distribution function that is highly non-Maxwellian. This typically arises in the ICE emitting region from drift orbit effects, for ICE driven by fusion-born ions, and from local ionisation and prompt loss effects, for ICE driven by NBI ions.
In recent years, electromagnetic emission in the radio frequency (RF) range has been observed from ELMy H-mode plasmas in the KSTAR tokamak [10,11,41,42]. The emission is highly dynamic, evolving throughout the duration of the edge localised mode (ELM) crash. There are two main distinct phases: spectral lines at multiple harmonics of the local ion cyclotron frequency at the outer midplane edge pedestal of the KSTAR plasma, at times ∼10 µs to 100 µs before the ELM crash; and broadband RF emission, sometimes chirping in frequency, during the pedestal collapse following the filament burst [41]. These two phases can be seen in figure 1, which has been adapted from figure 2 of [11]. Time intervals A and B broadly correspond to 'phase 1', and C and D to 'phase 2'. Data such as that displayed in figure 1 is obtained using a fast radio frequency (RF) spectrometer consisting of a wide band antenna (either Bowtie or Spiral) and a high speed digitiser. The operating frequency of the antennae spans a wide range from 100 MHz to 2 GHz, and the high speed digitiser records data at up to 20 giga-samples per second. For a more thorough account of the detection system, the differences between regions A to D, and their relationship to ELM dynamics, we refer to [10] and [11].
In this paper we address steady state ICE characterised by deuteron, as distinct from proton, cyclotron harmonic structure, which typically occurs during 'phase 1'. Before doing so, let us contextualise our study by offering a brief account of the proton chirping ICE observed during 'phase 2'. The proton chirping ICE observed during 'phase 2' involves a fast redistribution of intensity across the fixed multiple cyclotron harmonics at which the ion cyclotron emission (ICE) is observed. Using particle-in-cell (PIC) [43] simulations [44], this chirping has been attributed to the rapidly changing spectral characteristics of the MCI which are caused by the rapidly decreasing local plasma density. This is turn is due to the motion of ELM filaments during the pedestal collapse. Chirping ICE spectral peaks in KSTAR plasmas are observed at harmonics of either the deuteron cyclotron frequency f cD , or the proton cyclotron frequency f cp = 2 × f cD . The only protons in these KSTAR deuterium plasmas heated by deuterium neutral beam injection (NBI) are those born from deuteron-deuteron fusion reactions in the core plasma. These have initial energy 3MeV, and were previously thought to be unconfined. However the ICE observations at harmonics of f cp prompted the test particle calculations of [44] which show that a subset of the fusionborn protons remain confined on deeply passing orbits whose drift excursions traverse the ICE emitting region. In addition, some KSTAR pulses exhibit much fainter, higher frequency ICE spectral features at the proton cyclotron frequency which are time delayed with respect to the main chirping feature. One such feature, also displaying downward proton chirping, was explained using PIC simulations combined with spectral analysis of the data, as being driven by nonlinear wave-wave interactions within the lower frequency chirping ICE feature [45].
Turning to the steady state deuteron ICE in 'phase 1', we consider two KSTAR deuterium plasmas, in which beams of 80-100 keV deuterons were injected tangentially with NBI heating power in the range of 2.8-3.8 MW. Detailed calculations of NBI ion losses suggest that some of these NBI deuterons can remain confined on trapped orbits [46]. We therefore investigate whether this small subset of NBI deuterons can give rise to ICE via the MCI in the outer midplane edge plasma. In particular we study the fast collective relaxation of these NBI deuterons by means of a direct numerical kinetic treatment. We self-consistently solve the Maxwell-Lorentz system of equations for more than 100 million computational macro-particles, using the EPOCH PIC code [43] with one periodic spatial dimension and all three velocity dimensions (1D3V). We follow the full gyro-orbit kinetics of the electrons, background deuterons, and the minority energetic NBI deuteron population. The velocity-space distribution of the energetic ions is initialised as a ring-beam dis- here v and v ⊥ are the magnitudes of the velocity components parallel and perpendicular to the magnetic field respectively, and u is the magnitude of the initial perpendicular component of NBI ion velocity. In these computations, we consider only perpendicular spatial propagation k = 0 , so that the distribution of v is immaterial and can be simplified to a delta function. A more extensive summary of the EPOCH code is given in appendix. Our simulations are set up in slab geometry corresponding to the locally uniform approximation with '1D3V' Cartesian phase space (x, v x , v y , v z ). This represents all three directions in velocity space, but only one direction in configuration space. Our z-direction is the direction of the magn etic field vector, corresponding approximately to the toroidal direction; and our x-direction would correspond to the radial direction. Importantly, though, we do not model tokamak geometry and the associated eigenmode structure [33][34][35][47][48][49][50][51][52][53][54], drifts or any other non-uniformity or transport effects; so any mapping from our co-ordinates to toroidal co-ordinates applies only locally.It seems likely that the latter effects may have only marginal consequences for the ICE phenomenology addressed here, because of the success of the locally uniform approximation in explaining recent results from KSTAR [44,45] and LHD [39,40], as well as the ICE observations from JET and TFTR [1-5, 36, 37].
In section 2 we carry out the preliminary studies that enable us to identify the subset of the NBI deuterons that is likely to drive deuterium cyclotron harmonic ICE from KSTAR. These studies combine particle orbit calculations with analysis of linear MCI growth rates for deuterons with different energies and pitch angles. In particular, we identify the types of NBI deuteron orbits that can give rise to a population inversion in velocity space, spatially localised to the outer midplane edge, that is substantial enough to drive the MCI strongly. In section 3 we present the main results of this paper: direct numerical simulations of the fast collective relaxation of the NBI deuteron population, evolving self-consistently with the electric and magnetic fields which they are found to excite under the Maxwell-Lorentz system of equations. The outputs of these simulations are then compared with observed ICE power spectra from KSTAR plasmas. Section 4 offers a summary, along with comments on the potential benefits of possible future ICE measurements on ITER [25].

Identifying NBI deuterons in KSTAR that could relax via the MCI
In order to carry out direct numerical simulations using the PIC self-consistent kinetic approach, it is first necessary to identify a population of NBI deuterons that could give rise to ICE at the outer midplane edge. To inform this search, we first calculate analytical linear growth rates of the MCI given by equation (29) of [30] (see also [31]) across a range of NBI injection energies and pitch angles φ = arcsin (v ⊥ /|v|). Time t = 0 denotes the time at which the first derivative of the RF signal is almost discontinuous, which coincides with the peak in the RF signal, the bursting phase [11]. During windows A and B, the ICE signal shows spectral peaks at successive harmonics of the deuterium cyclotron frequency. The ICE signal during windows C and D shows more complex burst phenomena, which are discussed in [10,11,44,45].
The expression for the analytical linear growth rate in [30] is derived from dispersion relations that were obtained by multiplying together dielectric tensor elements within the overall framework of Maxwell's equations. Textbook kinetic integral expressions for each dielectric tensor element were used for all three species-electrons, thermal ions, and minority energetic ions-and summed together. The non-Maxwellian velocity distribution of the energetic ions, when integrated over, yields distinctive contributions to the dielectric tensor element, which fundamentally alter the dispersion relation that would otherwise govern waves in the background plasma. These give rise to the linear MCI of waves on the fast Alfvén-cyclotron harmonic wave branches. When calculating these analytical linear growth rates, it is necessary to use an angle of wave propagation that is not strictly perpendicular to the magn etic field, and we select θ = 89 • . This very closely resembles our simulation set-up for the v v ⊥ cases considered in this paper, as the Doppler shift due to the finite parallel wave-number k is minimal. The results are shown as a contour plot in figure 2. The colour scale represents the log 10 of the growth rate of the fastest growing mode in both wavevector and frequency space. We see that NBI deuterons with higher energy and pitch angle have the stronger growth rates. This is expected, since previous studies show [36][37][38][39][40] that the MCI is more readily excited where v A is the local Alfvén speed; and, of course, v ⊥ increases with both energy and pitch angle.
We now calculate orbits for deuterons with an energy of 100 keV and large initial pitch angles. The choice of 100 keV reflects the initial energy of KSTAR NBI deuterons. If this energy is primarily perpendicular, it corresponds to v ⊥ /v A ∼ 0.68 in the KSTAR outer midpane edge plasma. This choice is also helpful from a computational physics perspective: previous PIC and PIC-hybrid simulations [36][37][38][39][40] have shown that a value of v ⊥ /v A close to unity allows excitation of the MCI in a feasible amount of computational time, while maintaining high signal-to-noise ratios [36,37,39,40]. The deuteron orbit calculations are carried out using the CUEBIT test particle code [55]. Defining poloidal flux ψ such that the poloidal magnetic field is equal to ∇ψ × ∇ϕ, where ϕ is toroidal angle in (R, ϕ, Z) cylindrical coordinates, we set Here R 0 , R b , γ and ψ 0 are constants which determine the plasma major and minor radii, the plasma elongation, and the plasma current. Equation (1) belongs to a class of solutions of the Grad-Shafranov equation obtained by Solov'ev [56], and is written in a convenient form proposed by Freidberg [57].
It is applicable to a plasma in which the pressure depends linearly on ψ, and RB ϕ is uniform, B ϕ being the toroidal magnetic field. We use the following parameters to approximate typical KSTAR equilibria: R 0 = 1.8 m, R b = 1.31 m, γ = 0.7, and ψ 0 = 0.36 T m −2 . Motivated by experiment [11], we first set B 0 = 1.8 T. We initialise particles near the core with a range of pitch angles, and find that some of these lead to orbits that traverse the edge region while remaining confined. The results of some of these calculations are shown in figure 3. Each panel represents a different initial pitch angle, and is accompanied by an enlarged plot of the edge region. By matching deuteron cyclotron harmonic ICE spectral peak frequencies to the known spatial dependence of magnetic field strength in KSTAR, we know that the emitting region is at R = 2.21 ± 0.023 m [11]. We therefore discard orbits that do not reach this far out in radius, in addition to orbits that cross the plasma boundary. In principle not all of the fast ions that cross the plasma boundary will impact on the first wall, and some will remain ionised in the plasma. While such particles could contribute in part to the ICE drive, they are likely to be much less abundant than particles within the plasma boundary, and as such are excluded in this initial study. We find that an initial pitch angle of 80 degrees results in an orbit that crosses the plasma boundary, while an initial pitch angle of 84 degrees results in an orbit which does not traverse the location of ICE emission. Initial pitch angles between these two values result in orbits that are both within the plasma boundary, and traverse the ICE emitting region for the given values of NBI injection energy and magnetic field. We also consider B 0 = 2.27 T, for which KSTAR deuterium ICE was also observed [10], and arrive at similar conclusions: only NBI deuterons with energy close to the injection energy and within a narrow range of high pitch angles, are capable of driving the MCI. Now that we have isolated the subset of NBI deuterons that is in principle capable of driving the observed ICE at deuterium harmonics, in section 3 we simulate their relaxation using a PIC code and analyse the outputs. Contour plot displaying the analytical linear MCI growth rate of the fastest growing mode as a function of pitch angle and particle energy, using a log 10 colour scale. Motivated by experimental observations, these calculations are restricted to frequencies below the 30th deuteron cyclotron harmonic. The linear MCI growth rate is exponentially strongest for pitch angles in the range 78 • < φ < 85 • . For a given pitch angle, the linear MCI is strongest at higher NBI deuteron energy; the strength of this dependence increases with pitch angle. The range of pitch angles displayed reflects the range for which the CUEBIT test particle code predicts orbits which are within the plasma boundary and traverse the ICE emitting region, see also figure 3.

Comparison between kinetic simulations and experiment
We have carried out three PIC simulations with parameters corresponding to two KSTAR plasmas that exhibit deuterium ICE leading up to the ELM crash. These simulations use an initially uniform electron number density n e = 2.5 × 10 19 m −3 , corresponding to the density at the top of the edge pedestal before the ELM crash, as inferred from Thompson scattering measurements [10]. This value of n e was the upper bound of the multiple simulations at different fixed density that were used previously [44,45] to show that the frequency chirping of KSTAR proton ICE is due to the density collapse. Our simulations use 2000 particles per cell, with over 50 000 cells so as to adequately resolve the Debye length λ D . This means the simulations evolve the dynamics of over 10 8 computational particles with physically correct electron-to-ion mass ratios. Each simulation lasts ten deuteron gyro-periods τ cD , by which time the instability is well into its nonlinear saturated regime, which is crucial for comparison with experiment [36,37,45]. The temperatures of the initially Maxwellian background thermal dueterons and electrons are set to 1 keV. The beam deuterons are initialised with a pitch angle of 72.4 • ; this is a value in the edge region which corresponds to a core value within the range 80 • < φ < 84 • identified by the CUEBIT orbit calculations in the preceding section. We denote the bulk and beam deuteron number densities by n D and n NBI respectively, and the fast ion concentration by ξ = n NBI /n D . Figure 4 shows the time evolution of the different particle and field contributions to the energy density in a PIC simulation with ξ = 10 −3 and B 0z = 1.44 T; this value corresponds to the magnetic field in the ICE emitting region of a KSTAR plasma which has central B 0 = 1.8 T. The energy transfer between particles and fields qualitatively resembles that of previous work [36,37], with the minority NBI deuterons transferring their energy to the bulk plasma and to the fields on timescales of the order of several ion gyro-periods. It is interesting to note that the energy density of the electrostatic component of the fields excited by the simulation is slightly larger than that of the electromagnetic component. This approximate equipartition of energy is common when the speed of the minority ions is significantly less than v A [39,40]. Confident that the energy transfer resembles that which is characteristic of the MCI, we now run two more simulations with ξ = 10 −2 . The two values of ξ used in our simulations are unrealistically large compared to that expected in KSTAR (for example, ξ is typically ∼10 −4 -10 −5 in JET [27]), and leads to an unrealistically short saturation time, but is necessary computationally in order to obtain reasonable signal-to-noise ratios. This choice is not expected to affect our conclusions, because the simulated ICE power due to the MCI has been found to scale linearly with fast particle concentration [38]. Increasing ξ in the simulations thus acts to shift the simulated signal above the noise, but with no significant consequences for the underlying physics. Even with this choice of ξ, the noise level for the higher, experimentally relevant, cyclotron harmonic spectral peaks that are excited in the PIC simulations still poses a challenge; this we address by using an unusually high number of computational particles. To provide a baseline which quantifies the effect of the noise, we have run two additional simulations which correspond to a thermal background plasma without the minority energetic NBI deuterons. These 'background' simulations have parameter sets which are otherwise identical to their MCI counterparts, and give rise to the green traces in figure 6 below. The residual spectral structure in the green traces reflects the concentration of noise energy at normal modes in line with the fluctuation-dissipation theorem [58]. Figure 5 shows the distribution of energy in the fluctuating z-component of the magnetic field, ∆B z , across frequency-wavenumber space, obtained from a simulation identical to that of figure 4, except for the larger value of ξ = 10 −2 . This plot is a spatio-temporal Fourier transform of ∆B z over the intervals 0 x 50 000λ D and 0 t 5τ cD . The vertical axis of figure 5 is normalised to the deuteron cyclotron frequency f cD , while the horizontal axis is normalised to f cD /V A , where V A is the Alfvén speed. The solid blue diagonal line denotes the initial velocity of the minority NBI deuterons, and the ratio of the perpendicular component of the NBI deuteron velocity to the Alfvén speed is ∼0.66. This sub-Alfvénic value highlights a computational challenge associated with resolving the high harmonics that are excited by the MCI.
The phase velocity of the fast Alfvén wave is significantly greater than the perpendicular velocity of the minority NBI deuterons v NBI , which is plotted as a blue diagonal line in figure 5. For a given perpendicular wave-number k ⊥ , this causes a divergence between the value of ω = l (2πf cD ) ≈ k ⊥ V A and the closest value of k ⊥ v NBI . The size of this divergence increases with harmonic number l. It follows that the phase velocity of MCI-excited waves with high wavenumbers, and hence high frequency, is further from that of the fast Alfvén wave than it is for low frequency MCI excited waves; in consequence, they have much lower growth rates. The spectral third (blue) the energy density of the electrostatic field E x ; fourth (green) the energy density of the magnetic field perturbation ∆B z ; fifth (cyan) the change in kinetic energy density of the minority energetic NBI deuterons. Time is normalised to the deuteron gyro period. The primary energy flow from the NBI deuterons is to the thermal deuterons, whose kinetic oscillation helps support the field oscillations excited by the MCI. These field oscillations include, with comparable magnitude, an electromagnetic component (∆B z ) 2 and an electrostatic component E 2 x . The electrostatic component involves electron kinetics which are fully captured in our PIC model. The MCI saturates within five deuteron gyroperiods. This plot is a spatio-temporal Fourier transform of the B z field over the intervals spanning 0 x 50 000λ D and 0 t 5τ cD . Shading indicates the log 10 of the spectral density of the oscillatory part ∆B z of the B z field component in frequencywavenumber space. The sweep of the fast Alfvén wave from bottom left to top right is intersected by cyclotron harmonic waves at successive deuteron harmonics. The phase velocity of the fast Alfvén wave v A , and this exceeds the speed v NBI of the NBI deuterons which is plotted as a blue diagonal line. Wave excitation is strongest in the wedge between v A and v NBI , and in particular where cyclotron harmonic waves intersect the fast Alfvén wave. Simulated ICE frequency spectra, such as the lower panels of figure 6, are obtained by integrating plots such as figure 5 over wavenumber.
power is greatest at locations in (ω, k) space where the deuteron cyclotron harmonic waves intersect the fast Alfvén wave, and these are visible at lower harmonics in figure 5. At higher wavenumbers and frequency, this effect is less striking, but can be distinguished from the background noise.
The spatiotemporal Fourier transform shown in figure 5 has been integrated across the domain k > 0 to obtain the power spectrum as a function of frequency. This spectrum is show as the blue trace in figure 6(c). Its counterpart from a simulation with B z = 1.84 T which corresponds to a KSTAR plasma with central magnetic field B 0 = 2.27 T, is show in figure 6(d). These power spectra constitute the PIC simulation counterparts of the experimentally measured ICE spectra. Figure 6(a) displays the power spectrum of the experimental counterpart to the PIC simulation which gives rise to (c): KSTAR plasma 16176, in which the spacing between deuteron cyclotron harmonics f cD ∼ 11.1 MHz. In all panels, the vertical axis is plotted on a log 10 scale while the horizontal axis is normalised to the deuteron cyclotron frequency. We see that the deuteron cyclotron harmonic ICE spectral peaks observed in the KSTAR plasmas are all excited by the collective relaxation of the energetic deuteron population in the PIC simulations. Across all deuteron harmonics the intensity of the MCI-excited waves (blue traces) is one or two orders of magnitude higher than that of the thermal plasma noise (green traces). We note that any apparent spectral structure in the thermal plasma noise arises from the fluctuation-dissipation theorem and identifies normal modes. It can be seen that the spectral peaks in the simulations are less well resolved in frequency than in the experimentally measured signal. This is because the instability in the simulation reaches saturation so rapidly, owing to the unrealistically large value of ξ which is necessary to achieve a sufficient signal-to-noise ratio. It is encouraging to see the broadening and diminishing amplitude, of spectral peaks at l 21 in both experiment and simulation.  Green traces provide a noise baseline for the blue traces. They are obtained from the thermal plasma without a ring-beam, so that any spectral structure arises from the fluctuation-dissipation theorem and identifies normal modes. l 11. However several spectral peaks with l < 11 are present in the simulation but are not observed in the experiment. We consider that these peaks may well have been excited, but were not detectable because the antenna used to measure the RF data whose power spectrum is plotted in figure 6(a) was a spiral antenna, optimised for circularly polarised waves. A different, Bowtie antenna, optimised for linearly polarised waves, was used for KSTAR plasma 11474, corresponding to figure 6(b). The S 11 return loss [59] of these two antennas is plotted in figure 7. For frequencies less than 150 MHz (l < 11), the return loss of the Bowtie antenna is close to 0 dB. This implies a high degree of reflectivity for an incoming signal, which would therefore be undetectable. Combined with the linear or circular optimisations of the two antenna polarisations, this offers a likely explanation of the mismatch between experimental and simulation spectral peaks at low harmonic numbers seen in figures 6(b) and (d).

Conclusions
In this paper we have combined the linear analytical theory of the MCI, energetic particle orbit studies, and first principles PIC simulations which self-consistently solve the Maxwell-Lorentz equations for fully kinetic ion and electron populations. This has enabled us to provide an explanation for the origin of ICE at multiple deuterium cyclotron harmonics, which is observed in KSTAR deuterium plasmas heated by deuteron NBI. We first identified a small subset of the NBI deuteron population that could be responsible for the emission. This was done by performing test particle calculations for deuterons with a range of NBI-relevant energies and pitch angles. We identified those that are confined, pass through the outer midplane edge region where the ICE originates, and have the largest analytical MCI growth rates. We then carried out two high resolution PIC simulations with parameters corresponding to the ICE emitting region at the outer mid-plane edge of two KSTAR plasmas. The initial velocity distributions of the kinetic thermal deuterons and electrons in these PIC simulations is Maxwellian. In addition, there is an initial minority energetic deuteron population, whose velocity distribution reflects NBI parameters and our orbit studies. The collective relaxation of the NBI deuteron population in these two PIC simulations generates electric and magnetic field oscillations whose power spectra substantially resemble the measured ICE spectra. Some low harmonic peaks in one simulation frequency spectrum were not detected in its experimental counterpart. A probable explanation for this is that the S 11 return loss of the Bowtie antenna used to measure this RF signal was close to 0 dB, implying very high reflectivity in this low frequency range.
The present study complements the analysis in [44,45] of ICE from KSTAR deuterium plasmas, in which the observed spectral peaks correspond to harmonics of the proton cyclotron frequency. This was explained in terms of regions of the collective relaxation of a subset of confined fusion-born protons through the MCI, and occurs around 100 µs after the stage in the ELM dynamics considered here. We have now explained the main features of two distinct types of ICE from KSTAR plasmas: here, that due to NBI deuterons; and in [44,45], that due to fusion-born protons. We note that upward and downward chirping ICE with spacing f cD has been observed in [10,11], and this is a phenomenon under investigation.
The ICE physics addressed in this paper and in [10,11,44,45], would have gone unnoticed had it not been for KSTAR's sophisticated RF system and high speed digitizer [42]. Their uniquely high time resolution has already yielded new insights into the dynamics of ELMs [44,45], energetic ions, and wave phenomena [45] in tokamak plasmas. ICE is typically driven by two classes of energetic ion population, both of which are potentially of interest for ITER [60,61]. The first comprises energetic ions that are marginally confined in the outer edge region, whether fusion-born or NBI. The second, following recent observations of transient ICE from the centre of ASDEX-Upgrade [8,13,14] and DIII-D [15,16], probably comprises ions born in the core plasma while the fusion reactivity is rising. These populations spontaneously communicate with the observer through the medium of ICE. Provided that the physics underlying ICE is understood, information about such energetic particle populations and their plasma environment could be extracted from future ICE measurements in Figure 7. The return loss (S 11 ) as a function of frequency of the Spiral (blue) and Bowtie (red) antennas used for the detection of signals in (a) and (b) respectively. The y -axis denotes the return loss and is on a log scale. For frequencies less than 150 MHz, the return loss of the Bowtie antenna is close to 0 dB, offering a likely explanation of the mismatch between experimental and simulation spectral peaks at low harmonic number (l 11) seen in figures 6(b) and (d).
ITER. To develop this understanding in present-day large MCF plasmas, an extensive effort is under way, as outlined in section 1, to which the present paper contributes.