Nonlinear wave interactions generate high-harmonic cyclotron emission from fusion-born protons during a KSTAR ELM crash

The radio frequency detection system on the KSTAR tokamak has exceptionally high spectral and temporal resolution. This enables measurement of previously undetected fast plasma phenomena in the ion cyclotron range of frequencies. Here we report and analyse a novel spectrally structured ion cyclotron emission (ICE) feature in the range 500 MHz to 900 MHz, which exhibits chirping on sub-microsecond timescales. Its spectral peaks correspond to harmonics l of the proton cyclotron frequency fcp at the outer midplane edge, where l  =  20–36. This frequency range exceeds estimates of the local lower hybrid frequency fLH in the KSTAR deuterium plasma. The new feature is time-shifted with respect to a brighter lower-frequency chirping ICE feature in the range 200 MHz (8fcp) to 500 MHz (20fcp), which is probably driven (Chapman et al 2017 Nucl. Fusion 57 124004) by 3 MeV fusion-born protons undergoing collective relaxation by the magnetoacoustic cyclotron instability (MCI). Here we show that the new, fainter, higher-frequency chirping ICE feature is driven by nonlinear wave coupling between different neighbouring spectral peaks in the lower-frequency ICE feature. This follows from bispectral analysis of the measured KSTAR fields, and of the field amplitudes output from particle-in-cell (PIC) simulations of the KSTAR edge plasma containing fusion-born protons. This reinforces the identification of the MCI as the plasma physics process underlying proton harmonic ICE from KSTAR, while providing a novel instance of nonlinear wave coupling on very fast timescales.

The magnetoacoustic cyclotron instability (MCI) is of fundamental interest because it combines three of the most characteristic and distinctive features of MCF plasma physics. These are: the cyclotron gyration of ions, which underpins their confinement; the fast Alfvén wave, which the MCI excites; and the existence of distinct minority energetic ion populations, notably fusion-born ions. The MCI can arise at frequencies such that the cyclotron frequency (and its harmonics) of the energetic ions matches the frequency of the Alfvén wave. This is a necessary, but not sufficient, condition for the MCI to occur. In addition, the velocity distribution of the energetic ion population must be sufficiently non-Maxwellian that sufficient free energy is available to be unlocked via resonant processes. The first process, identified by Belikov and Kolesnichenko in 1976 [16][17][18], involves wave-wave resonance between additional cyclotron harmonic modes supported by the energetic ion population, and the fast Alfvén wave. The second process [19,20] involves wave-particle cyclotron resonance between gyrating energetic ions and the fast Alfvén wave. In general there will be a link between the spectrum of waves excited by the MCI, and the character of the velocity distribution of the energetic ions which drive the waves. Mapping between the measured ICE signal and the theoretical characteristics of the MCI-excited waves is thus essential to interpreting the physics underlying ICE observations. The duration of the proton ICE features during KSTAR ELM crashes is brief, typically a few microseconds; see for example, figure 6 of [1] and figure 1 below. Temporal correlation between an ELM crash and the ICE phenomena observed on KSTAR [1] may be due in part to the action of the ELM in 'flushing out' particles from the edge region. Focusing on the consequences for the fusion-born proton populations in the edge region, soon after the ELM-crash, only freshly born protons close to the birth energy of 3 MeV are present; slowing-down due to collisions has not yet had time to act. The observed frequency chirping has been explained [2] in terms of rapid changes in the density of the ambient plasma in which the energetic ions are embedded, caused by ELM filament motion during the crash. The changes in density alter the spectral character of the MCI-excited waves. Hence these chirping spectral features have been used to obtain uniquely high time resolution measurements of the time-varying local plasma density. Some of the KSTAR chirping ICE features below ∼500 MHz are observed to be accompanied, after a slight time delay, by a fainter detached 'ghost' chirping feature in the higher frequency range 500 MHz to 900 MHz; see in particular figure 1. This frequency range exceeds the local lower hybrid frequency corresponding to the local plasma parameters inferred from [1]; consequently, linear cold plasma waves propagating quasi-perpendicular to the magnetic field are expected to be evanescent here, see appendix A. The question therefore arises whether the 'ghost' ICE chirping features detected in KSTAR may reflect instrumental and signal processing issues, or are a real plasma physics effect. It is important to resolve this issue, because understanding observations (if real) of radiation in this frequency range can assist understanding of the physics of energetic ion populations in magnetically confined fusion (MCF) plasmas [33].
Here we show that the 'ghost' chirping ICE feature above ∼500 MHz in figure 1 is a real physical phenomenon, which is generated by strong nonlinear wave-wave coupling between different spectral peaks within the primary chirping ICE feature below ∼500 MHz. We demonstrate this by direct bicoherence analysis (see appendix B and [34][35][36][37][38]) of: first, the KSTAR data files for field magnitudes; and second, the fields generated from direct numerical simulations using the particle-in-cell (PIC) [39] code EPOCH [40]. We solve the self-consistent Maxwell-Lorentz system of equations for fully kinetic electrons and thermal deuterons, together with the minority ring-beam distribution of confined fusion-born 3 MeV protons that drives the primary ICE. The code retains full gyro-orbit kinetics for each of the ∼1 million macroparticles in the simulation. Full gyro-orbit kinetics are essential for capturing cyclotron harmonic resonance effects including, as we shall see, coupling between modes driven at different harmonics by collective instability. The simulations detailed in [2,14,[23][24][25] are set up in slab geometry, hence they do not incorporate realistic toroidal geometry and the associated compressional Alfvén eigenmode structure [26][27][28][29][30]. Nevertheless this spatially localised physics approach is successful in capturing most of the observed features of ICE, including also recent results from the heliotron-stellarator LHD [31,32]. Identification of the physics effects underlying the simulation output is assisted by the fact that, in the linear regime, the simulation approach aligns with the original slabgeometry analytical theory of the MCI. We refer to appendix C for a summary of the EPOCH PIC approach and [2] for details of the plasma parameters used in the simulations.
Bicoherence analysis techniques [34][35][36][37][38] are designed to capture nonlinear wave-wave coupling, and we refer to appendix B for a brief account. We show that the 'ghost' spectral features are able to exist and grow in the higher-frequency, potentially evanescent, region because they are nonlinearly driven by coupled MCI-excited waves that lie within the lower-frequency, propagating (non-evanescent), region. The 'ghost' feature thus owes its existence to both a minority suprathermal ion population -here, the confined subset of fusion-born protons that relax through the MCI [2] in KSTAR deuterium plasmas -and on the capacity of the plasma to nonlinearly couple together the modes initially driven by these protons.

The role of the lower hybrid frequency in the PIC simulations
In [2], a series of PIC simulations at successive neighbouring fixed values of plasma electron number density n e were used to model the chirping of the primary ICE feature in figure 1. We first address the role of the lower hybrid frequency f LH defined by equation (A.1), and in particular the question of evanescence for cold plasma waves at frequencies above f LH , including potentially those in the 'ghost' feature shown in figure 1.
It is well known [41] that perpendicular-propagating linear waves in a cold plasma cannot exist between the lower hybrid resonance frequency f LH and a cut-off frequency f 2 > f LH . Expressions for these frequencies, along with references to relevant literature, are given in appendix A. The results reported in [2] indeed show dependence on f LH , whereby the number of modes available for excitation decreases rapidly as the electron number density n e , and hence f LH , decreases. In figure 2, we show the results of multiple simulations of the MCI for two values of magnetic field strength B and a range of densities n e . In all of these simulations, the magnetic fields, ∼1.44 T (left panel) and ∼1.52 T (right panel), are oriented entirely along the z-axis of the simulation. That is, we restrict the study to strictly perpendicular wave propagation along the x-axis, which is the spatial domain of our 1D3V PIC simulations. The magnetic fields have been chosen to be representative of the magnetic field strength in the ICE-emitting region at the outer midplane edge of KSTAR at major radial position R ∼ 2.25 m(R 0 = 1.8 m, a = 0.5 m) for central field strengths B 0 = 1.8 T and B 0 = 1.9 T respectively. The value of n e in each vertical strip in figure 2 decreases from left to right in steps of 0.2 × 10 19 m −3 , with each vertical strip corre sponding to an independent simulation which yields the spectrum of MCI-excited waves at the value of n e shown. The rest of the simulation parameters are those given in section 4 of [2]. In both panels, shading indicates the log 10 of the Fourier power in the B z component of the simulation, and horizontal white dashed lines denote successive proton cyclotron harmonics. The white crosses, joined by a solid curved white line, denote the value of f LH at the density shown. In both panels, one can see a blue region in which the spectral power falls to zero. The boundary of this region at each value of n e lies close to the corresponding value of f LH .

Bicoherence analysis
The signal used to generate the spectrogram shown in figure 1 was obtained during a KSTAR ELM crash using a fast radio frequency (RF) spectrometer sampling at 5 GHz. Thus the In addition to the main chirping feature 500 MHz ≈20 f cp discussed in [2], we also observe a second, faint ('ghost'), feature at frequencies above the lower hybrid frequency f LH ≈ 529 MHz ≈ 21 f cp . This additional, spectral feature is delayed in time by approximately 1 μs with respect to the main chirping feature. Reproduced from [2]. © 2017 Culham Centre for Fusion Energy. CC BY 3.0. maximum resolvable (Nyquist) frequency is 2.5 GHz. These data were obtained in 200 μs segments 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 [1]. In figure 1, t = 0 refers to the centre of a 200 μs segment of RF data. Further details of the fast RF spectrometer system and the experimental set-up are given in [1].
We first examine the extent of nonlinear wave-wave coupling within the experimental dataset that spans the primary and 'ghost' chirping ICE features in figure 1. The best quantitative evidence for this coupling, and characterisation of its magnitude as a function of wave frequency, is obtained from bispectral analysis [34][35][36][37][38], see also appendix B.
The bispectrum, equation (B.1), measures the extent of phase coherence due to the nonlinear coupling between three waves that satisfy the frequency and wavenumber matching criteria covered in appendix B. The bicoherence, equation (B.2), is a normalised bispectrum bounded between 0 and 1 which quantitatively measures the fraction of the Fourier power of a signal that is due to nonlinear (specifically quadratic) interaction.
Thus the bicoherence sheds light on nonlinear coupling; whereas the bispectrum yields information regarding the energy flow due to nonlinear coupling, given the wave amplitudes in the system. It is therefore useful to compute them both when diagnosing possible nonlinear wave physics. A large value of bicoherence (close to unity) may reveal waves which have significant coupling, but do not drive additional waves in practice due to their relatively low amplitudes. This becomes apparent if one supplements the information given by the bicoherence with the bispectrum, because the latter also incorporates information about relative wave amplitudes. Conversely, plotting the bispectrum alone does not necessarily yield information about the intrinsic strength of coupling between waves. Bispectral analysis has been successfully applied to the MCI [23] and experimental plasma measurements [42][43][44], including those in the KSTAR tokamak [45].

Bicoherence analysis of KSTAR ICE chirping data
The bicoherence and bispectrum corresponding to the entire KSTAR signal shown in figure 1 are plotted in the left and right panels respectively of figure 3. In the bicoherence panel, shading indicates the intrinsic strength of nonlinear coupling, 1 (dark red) being completely coupled and 0 (dark blue) completely uncoupled. The shading of the bispectrum panel is displayed on a logarithmic scale. Here the averages denoted by · in equation (B.2) are taken over a time window ∆t ∼ 0.5 μs within a signal which is 5 μs long, corresponding to the data displayed in figure 1. This choice enables us to construct ten independent realisations. In consequence, the threshold for significance is comfortably below the observed coupling strength 'b' for a wide range of relevant frequencies.
We note three distinct regions of strong intrinsic nonlinear wave coupling in the left panel of figure 3: Coupling between neighbouring modes within the main chirping feature shown below f ≈ 450 MHz in figure 1. We argue that this coupling enables formation of the faint higher frequency 'ghost' chirping feature that appears above f ≈ 450 MHz in figure 1. We are primarily concerned with point (i), which strongly suggests the 'ghost' feature is a real plasma physics phenomenon. The right panel of figure 3 indicates why it is only waves in the frequency range below f ≈ 450 MHz that can drive the observed 'ghost' features: these are the waves that are not only significantly nonlinearly coupled, but also have sufficiently large amplitude. The nonlinearly driven features that could in principle arise due to the strong coupling of waves described in points (ii) and (iii) would lie below the Nyquist frequency; however, they are never observed in practice because their amplitude is lower by several orders of magnitude. We note that the auto-bispectrum and auto-bicoherence of the KSTAR RF data, that is, bispectra computed from a single time series, cannot by themselves yield information on the direction of energy transfer. To do so would require two point measurements [34,46] which at present, are not available.

Bicoherence analysis of the PIC simulation output
Having inferred from bispectral analysis of the KSTAR data that the nonlinear wave coupling between cyclotron peaks below f ≈ 500 MHz drives the 'ghost' chirping feature, the question now arises: can the same physics be inferred from analysis of the outputs of the corresponding PIC simulations? The simulations have a propagation angle θ = 90 • , for which, as noted above, the region f LH < f < f 2 is evanescent. In order to explore the hypothesis that the observed waves in this region arise from nonlinear wave coupling, let us focus on the simulations which make up the lower panels of figure 4 in [2]. Figure 4 shows the bicoherence plots along with the corre sponding spatio-temporal Fourier transform of B z for each of three different simulations in the lower panels of figure 4 in [2]. Shading indicates the log 10 of the spectral density of the oscillatory part of the B z field component. All other plasma parameters remain identical and are specified in [2]. In the lower panels, the y-axis is plotted in units of proton cyclotron frequency f cp , while the x-axis is plotted in units of f cp /V A where V A is the Alfvén speed. The value of V A differs significantly between the simulations because it is inversely proportional to the square root of the majority ion (deuteron) mass density, and hence to n e . The horizontal black line denotes f LH , below which we see excitation of the fast Alfvén wave with resonances at consecutive proton cyclotron harmonics, characteristic of the MCI [2,14,[23][24][25] which underpins ICE. Above f LH there are several weaker but significant spectrally intense regions. The location of these regions are the locations of strong resonances on the fast Alfvén branch below f LH . This condition for wave-wave coupling is necessary for conservation of momentum and energy [34]. We also note that the most dominant nonlinear spectral features above f = f LH move to increasingly high values of normalised k as density increases.
If the spectrally dense regions with co-ordinates (k 3 , f 3 ) above f LH are indeed the result of wave-wave coupling between modes below f LH , this should be borne out by bicoherence analysis of the simulated field component B z . The corre sponding bicoherence plot for each simulation is shown in the upper panels in figure 4. These plots show clearly defined sets of (k 1 , k 2 ) pairs which have strong coupling, the most striking of which are near the k 1 = k 2 (and hence f 1 = f 2 ) boundary. These are modes close to each other in k space on the fast Alfvén branch. If we pick a region of strong coupling near the k 1 = k 2 boundary for the upper leftmost panel, say k 1 ≈ 15f cp /V A and k 2 ≈ 18f cp /V A , and read off the corresponding f 1 ≈ 12f cp and f 2 ≈ 14f cp , then we should be able to see a spectrally dense region at k 3 ≈ 33f cp /V A and f 3 ≈ 26f cp in the lower leftmost plot above the f = f LH line. This is indeed the case, and a similar correspondence is seen across all panels of figure 4.
Bicoherence analysis of both experimental data (figure 3) and simulation outputs (figure 4) thus demonstrates strong coupling between modes near the f 1 = f 2 boundary below f LH . This supports our conjecture that nonlinear coupling is responsible for the faint spectral 'ghost' feature in figure 1, since this is also captured by our simulations. This lends further credence to our interpretation in [2] that the downward ICE chirping is due to declining local plasma density, which is perhaps associated with the motion of an ELM filament.

Density dependence of downward chirping
Let us now investigate in greater depth the hypothesis that the local decline of density on submicrosecond timescales may be responsible for the downward chirping characteristics of the 'ghost' ICE feature in figure 1. Due to the abundance of waves in the simulation there are many spectrally dense regions in the f > f LH regions in figure 4. Accordingly, we adapt and extend the technique which was previously applied in [2] to ICE chirping at frequencies less than f LH in KSTAR. Key to this approach is analysis of the spectral properties of multiple PIC simulations, each of which is run into the nonlinear regime of the MCI at different, fixed, neighbouring values of n e .
(1) Using the experimental bicoherence plot (figure 3) along with the experimental spectrogram (figure 1), we identify spectral features 'f 1 ' and 'f 2 ' with f < f LH , that are able to combine to produce the faint spectral features ' The simulation with number density n e , which in [2] was found to give rise to strong spectral features with frequencies f 1 and f 2 , see the left panel of figure 5, is examined.
In cases where f 1 and f 2 are present across a range of n e values, the procedure is repeated for each simulation.  then obtained by integrating and averaging between the minimum and maximum possible values of k 3 . (6) As there is a one-to-one mapping between k 1 and f 1 , and between k 2 and f 2 , there is an approximate one-to-one mapping between k 3 and f 3 . Therefore the spectral power in k 3 corresponds to the power in the vicinity of f 3 . Figure 4 of [2] is reproduced here as the left set of panels in figure 5. In the lower panels of figure 5, the spectral power in the output of multiple simulations is plotted as a function of frequency and n e , and compared with the experimental RF spectrum (upper left panel). The mapping between nearidentical spectral features in the experimental data and the simulation outputs was used in [2] to infer the time-dependence of local density. The range of values of n e shown reflect Thomson scattering measurements in the edge pedestal (see the last paragraph of section 3 in [1]. If the faint chirping features in figure 1 are a result of wave-wave interactions between modes with f < f LH , driven by the MCI at different densities, we expect the spectral power of the newly formed modes with f > f LH to exhibit a similar dependence on frequency and electron number density. To this end, the power in these modes for each simulation has been calculated, and the results are shown in the right panels of figure 5. For comparison purposes, figure 4 of [2] is reproduced as the left panel of figure 5. The procedure is as follows: The lower left and lower right panels in figure 5 have much in common. First, in each case the dominant spectral features of the simulations chirp down in frequency as electron number density decreases. Second, the density values over which this occurs declines from the pre-crash pedestal density to much smaller values, in both cases.

Conclusions
The 'ghost' ICE feature in figure 1 is a real plasma physics phenomenon. Its existence is due to a combination of energetic particle physics with linear and nonlinear wave physics, which is so far observed only in KSTAR tokamak plasmas.
The entire phenomenology underpinning the 'ghost' unfolds on sub-microsecond timescales during an ELM crash, and the frequency chirping reflects declining local deuterium plasma density. This density evolution changes the plasma environment of the 3 MeV fusion-born protons which drive the ICE through collective relaxation by the magnetoacoustic cyclotron instability (MCI), resulting in fast evolution of the spectral distribution of energy in the excited fields. Here we have shown that the separate, fainter ('ghost') chirping ICE feature observed in the frequency range 500 MHz (20f cp ) to 900 MHz (36f cp ) is driven by nonlinear wave coupling between different neighbouring cyclotron harmonic peaks in the main ICE feature below 500 MHz. This is evident from bispectral analysis of: first, the measured KSTAR fields, where we benefit from exceptionally high (up to 20 GS s −1 ) sampling rates; and second, field amplitudes output from first principles particle-in-cell code simulations of the KSTAR fusion-born proton relaxation scenario. This reinforces the MCI interpretation of chirping proton ICE in KSTAR [2]. It also provides a novel demonstration of nonlinear wave coupling on very fast timescales in a tokamak plasma.
The successful interpretation of this unexpected phenomenon spontaneously driven by fusion-born ions, helps establish interpretive capability for ICE from future deuterium-tritium plasmas in JET and ITER. Modelling of the plasma physics underlying ICE signals yields information on two key features of the ICE-emitting energetic ion population. First, the values of key parameters, notably the ratio of the characteristic perpend icular velocity of the energetic ions to the local Alfvén speed. This needs to be of order unity [16][17][18][19][20][21][22][23][24][25][26][27][28]. Second, the structure of the distribution of the energetic ions in velocity space, which needs to be strongly non-Maxwellian in order to excite the MCI which underlies ICE. A drifting ring distribution, i.e. the product of two delta-functions in parallel and perpendicular velocity, has been found to be the best few-parameter way of capturing this structure for ICE applications. This approximation has proven fruitful across more than two decades, spanning ICE measurements from deuterium-tritium plasmas in JET [11] and TFTR [12] during the mid-1990s to the most recent measurements reported from ASDEX-Upgrade in 2014 [8] and JT-60U in 2017 [13]. The new results presented here confirm the fidelity of the output of first principles PIC simulations in relation to measured ICE signals, alongside the validity of the model for ICE that is implemented in the PIC code. The agreement between the bispectral analysis of the simulation outputs and the observations of an unexpected, strongly nonlinear, transient ICE feature provides fresh validation of the ICE model, in a challenging regime. ICE is also of interest in that, stimulated by the external application of finite amplitude ICRF waves, ICE physics could contribute [47] to 'alpha-channelling'-the rapid inward transfer of energy from fusion-born ions to the bulk plasma. to evolve the relativistic particle trajectories. In the PIC simulations reported here and in [2], we initialise the thermal electron and majority deuteron populations as Maxwellian. The perpendicular velocity component of the minority fusion-born protons on deeply passing orbits that drive the ICE, corresponding to perpendicular energy 150 keV [2], is initialised as a ring distribution in velocity space.