Interpreting observations of ion cyclotron emission from large helical device plasmas with beam-injected ion populations

Ion cyclotron emission (ICE) is detected from all large toroidal magnetically confined fusion (MCF) plasmas. It is a form of spontaneous suprathermal radiation, whose spectral peak frequencies correspond to sequential cyclotron harmonics of energetic ion species, evaluated at the emission location. In ICE phenomenology, an important parameter is the value of the ratio of energetic ion velocity to the local Alfvén speed . Here we focus on ICE measurements from heliotron-stellarator hydrogen plasmas, heated by energetic proton neutral beam injection (NBI) in the large helical device, for which takes values both larger (super-Alfvénic) and smaller (sub-Alfvénic) than unity. The collective relaxation of the NBI proton population, together with the thermal plasma, is studied using a particle-in-cell (PIC) code. This evolves the Maxwell–Lorentz system of equations for hundreds of millions of kinetic gyro-orbit-resolved ions and fluid electrons, self-consistently with the electric and magnetic fields. For LHD-relevant parameter sets, the spatiotemporal Fourier transforms of the fields yield, in the nonlinear saturated regime, good computational proxies for the observed ICE spectra in both the super-Alfvénic and sub-Alfvénic regimes for NBI protons. At early times in the PIC treatment, the computed growth rates correspond to analytical linear growth rates of the magnetoacoustic cyclotron instability (MCI), which was previously identified to underlie ICE from tokamak plasmas. The spatially localised PIC treatment does not include toroidal magnetic field geometry, nor background gradients in plasma parameters. Its success in simulating ICE spectra from both tokamak and, here, heliotron-stellarator plasmas suggests that the plasma parameters and ion energetic distribution at the emission location largely determine the ICE phenomenology. This is important for the future exploitation of ICE as a diagnostic for energetic ion populations in MCF plasmas. The capability to span the super-Alfvénic and sub-Alfvénic energetic ion regimes is a generic challenge in interpreting MCF plasma physics, and it is encouraging that this first principles computational treatment of ICE has now achieved this.

Ion cyclotron emission (ICE) is detected from all large toroidal magnetically confined fusion (MCF) plasmas. It is a form of spontaneous suprathermal radiation, whose spectral peak frequencies correspond to sequential cyclotron harmonics of energetic ion species, evaluated at the emission location. In ICE phenomenology, an important parameter is the value of the ratio of energetic ion velocity v Energetic to the local Alfvén speed V A . Here we focus on ICE measurements from heliotron-stellarator hydrogen plasmas, heated by energetic proton neutral beam injection (NBI) in the large helical device, for which v Energetic /V A takes values both larger (super-Alfvénic) and smaller (sub-Alfvénic) than unity. The collective relaxation of the NBI proton population, together with the thermal plasma, is studied using a particle-in-cell (PIC) code. This evolves the Maxwell-Lorentz system of equations for hundreds of millions of kinetic gyro-orbit-resolved ions and fluid electrons, self-consistently with the electric and magnetic fields. For LHD-relevant parameter sets, the spatiotemporal Fourier transforms of the fields yield, in the nonlinear saturated regime, good computational proxies for the observed ICE spectra in both the super-Alfvénic and sub-Alfvénic regimes for NBI protons. At early times in the PIC treatment, the computed growth rates correspond to analytical linear growth rates of the magnetoacoustic cyclotron instability (MCI), which was previously identified to underlie ICE from tokamak plasmas. The spatially localised PIC treatment does not include toroidal magnetic field geometry, nor background gradients in plasma parameters. Its success in simulating ICE spectra from both tokamak and, here, heliotron-stellarator plasmas suggests that the plasma parameters and ion energetic distribution at the emission location largely determine the ICE phenomenology. This is important for the future exploitation of ICE as a diagnostic for energetic ion populations in MCF plasmas. The capability to span the super-

Interpreting observations of ion cyclotron emission from large helical device plasmas with beam-injected ion populations
Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Introduction
Suprathermal ion cyclotron emission [1,2] (ICE) is detected from all large toroidal magnetic confinement fusion (MCF) plasmas including the tokamaks TFR [3], PDX [4], JET [5], TFTR [6], JT-60U [7], ASDEX-U [8], KSTAR [9], DIII-D [10] and the stellarators LHD [11,12] and W7-AS [13]. ICE is notable as the first collective radiative instability driven by confined fusion-born ions that was observed in deuteriumtritium (D-T) plasmas in JET and TFTR [14][15][16][17]. The frequency spectrum of ICE typically exhibits narrow peaks at values which can be identified with sequential local cyclotron harmonics of a distinct energetic ion population in a spatially localised emitting region. The numerical value of the inferred ion cyclotron frequency Ω c = Z i eB/m i , where Z i is the ion charge and m i its mass, then determines the local value of the magnetic field strength in the emitting region, and hence its radial location. Typically, but not invariably, this is at the outer mid-plane edge of the toroidal plasma; ICE from the core plasma has been reported recently from DIII-D [18] and from ASDEX Upgrade [19], and earlier from JT-60U [7]. This development suggests great potential for the exploitation of ICE as a diagnostic for energetic particles in ITER [20].
Measurements of ion cyclotron emission (ICE) spectra have been obtained from heliotron-stellarator plasmas in the large helical device (LHD), both with an ICRH antenna during NBI heated plasmas [11] and by magnetic probes [12] during toroidal Alfvén eigenmodes (TAE's) [21][22][23][24]. Related numerical studies can be found in [24]. In combination with other advanced diagnostics, notably for MHD activity, these ICE measurements from LHD, can yield fresh insights into the physics of energetic ions in magnetically confined fusion (MCF) plasmas. We attribute this ICE to a neutral beam injected (NBI) proton population at energies ≈40 keV in the outer midplane edge regions of hydrogen plasmas in LHD [11,12], where the local electron temperature T e ≈ 20 eV to 150 eV, number density n e ≈ 10 19 m −3 and magnetic field strength B ≈ 0.5 T. These spectra were measured with an ICRF heating antenna in receiver mode. Importantly, these spectra span plasma regimes where the ratio of the velocity of the energetic ions V NBI to the local Alfvén speed V A in the ICE-emitting region of the LHD edge plasma takes values that can be either smaller or larger than unity. The transition between super-Alfvénic and sub-Alfvénic energetic ion phenomenology is of fundamental interest in MCF plasma physics, see section 2 below. Here, in particular, we examine LHD plasmas 79126 and 79003 where V NBI /V A = 0.872 and 1.125, respectively, in the ICE-emitting region.
Our interpretation rests on first-principles numerical solutions of the Maxwell-Lorentz system of equations. We follow the full velocity-space trajectories, including gyromotion, of tens to hundreds of millions of fully kinetic energetic and thermal ions, together with all three vector components of the evolving electric and magnetic fields, with a massless neutralising electron fluid, using a fully nonlinear 1D3V PIC-hybrid particle-in-cell code [25]. The kinetic ions, fluid electrons, and fields are coupled self-consistently through the Lorentz force and Maxwell's equations in Darwin's approximation [26]. In this hybrid scheme [25], the Debye length does not need to be resolved. It therefore requires less computational resource than the full PIC scheme implemented in EPOCH [27], which retains electron kinetics, and is also used for contemporary theor etical studies of ICE [28][29][30][31][32]. We follow these simulations through the linear phase of an instability that we show is the magnetoacoustic cyclotron instability (MCI) [25,29,[33][34][35][36], and then deeply into its nonlinear saturated phase. The Fourier transforms of the excited fields in these simulations yield frequency spectra that qualitatively match the observed ICE spectra from the LHD plasmas. These simulation results for heliotron-stellarator plasmas complement and confirm earlier interpretation of ICE driven by sub-Alfvénic NBI ions in TFTR tokamak plasmas [15,37]. In general, ICE phenomenology in MCF plasmas has been found to reflect the plasma parameters and magnetic field strength in the emitting region, together with the velocity distribution of the driving energetic ion population. Important aspects of these two key features come together in the dimensionless parameter v Energetic /V A , which is the ratio of energetic ion velocity v Energetic to the local Alfvén speed V A .

The wider experimental context and motivation
In this work we focus on the role of the parameter v Energetic /V A and, in particular, whether it exceeds, or is less than, unity. In general, one anticipates differences in how a particle interacts with the coherent oscillations that are supported by a continuous medium, depending on whether the particle is travelling faster or slower than the speed at which phase information can propagate in the medium. In a magnetised plasma, in the frequency range of interest to ICE, this speed is the Alfvén speed V A . Even if no complete experimental or theoretical information related to this frequency and velocity range in fusion plasmas is available, it would be necessary to investigate it when establishing the physics basis for the exploitation of ICE. However there exists, in addition, extensive evidence for the importance of the value of the ratio v Energetic /V A compared to unity: in relation to the ICE detected from fusionborn alpha-particles in the only DT MCF plasmas thus far, in JET and TFTR during the 1990s; and in relation to the fundamental theory of the MCI.
The fusion-born alpha-particles responsible for the ICE observed throughout the duration of JET DT H-mode plasmas with high edge density were super-Alfvénic, see the sixteenth item in table 1 of [14]; whereas in TFTR supershot DT plasmas with strong central density peaking and low edge density, they were sub-Alfvénic, as shown in figure 6 of [15]. Apparently related to this, and shown in figure 7 of [15], alpha-particle ICE from TFTR DT plasmas arose only during the first 100 ms to 200 ms, in the early stage of the density ramp, notwithstanding the greater measured production rate for fusion alpha-particles at later times with higher density. Conversely, in TFTR plasmas with edge densities such that newly born alpha-particles were super-Alfvénic, the corresponding alphaparticle ICE spectrum persisted, see [15]. Evidently the value of v Energetic /V A is central to the ICE phenomenology in the DT plasmas which, thus far, provide the most relevant signposts for ICE if it arises in ITER.
The earliest study [37] of ICE driven by NBI ions in TFTR DT plasmas reinforced the centrality of v Energetic /V A , bearing in mind also that NBI ICE was not detected from JET DT plasmas. For example, deuterium and tritium NBI ions were at v Energetic /V A = 0.1 in the plasma region where they excited ICE. The analytical MCI theory available at that time [34][35][36] was aligned to the JET super-Alfvénic alpha-particle regime, and could not immediately address the strongly sub-Alfvénic TFTR NBI regime; instead an interpretation was obtained based on a primarily electrostatic treatment [37].
The foregoing motivated the extension of the analytical theory of the MCI into the sub-Alfvénic regime, reported at length in [36], and exploited in [14] and in relation to both NBI ions and fusion-born ions in [16,38]. By this stage, the linear analytical theory of the MCI had successfully exhausted its potential at the leading edge of ICE interpretation. Only from 2010 onwards was it possible to carry out first principles kinetic calculations using PIC codes, which carry the MCI selfconsistently into its fully nonlinear regime [25,[30][31][32]39]. The saturated power spectra emerging from these computations then provide theoretical counterparts to measured ICE spectra. Hitherto, with one exception, these treatments have addressed only the super-Alfvénic regime for fusion-born ions. The exception is the set of multiple PIC simulations [31,32] of ICE driven by fusion-born protons under rapidly evolving edge plasma conditions in KSTAR. For local electron densities below n e ∼ 1.05 × 10 19 m −3 , corresponding to the lower panels of figure 4 of [31], the perpendicular velocity of the protons is sub-Alfvénic.
It is therefore timely to address the standing issue of the impact on ICE of the value of the ratio v Energetic /V A compared to unity, by using a combination of contemporary NBI ICE measurements and PIC computations. Here we exploit the fact that ICE from the outer midplane edge of LHD heliotronstellarator plasmas is generated by NBI ions in both super-Alfvénic and sub-Alfvénic regimes. We use the hybrid PIC code of [25], as previously exploited along with the EPOCH PIC code [27] in [25,[30][31][32]39], to simulate the collective relaxation of these NBI ion populations under LHD-relevant conditions.
To progress, several novel steps are necessary, which are reported in this paper. We carry out the first extensive set of 1D3V kinetic PIC-hybrid simulations of ICE in the sub-Alfvénic regime, where the minority NBI ion population has kinetic energy more than an order of magnitude lower than in previous simulations [25,[30][31][32]39] relating to super-Alfvénic fusion-born ions in JET. These are also the first ICE simulations with wavevectors inclined more than one degree from perpendicular to the magnetic field direction. The analytical theory of the linear MCI which underlies ICE is well developed, and we are able to show that its predictions regarding initial growth rates in our fully nonlinear simulations are successful, for the first time in the sub-Alfvénic regime.

Analytical theory of the magnetoacoustic cyclotron instability
The magnetoacoustic cyclotron instability (MCI) [25,29,33,34,36,37,38] is the most likely emission mechanism to account for ICE generation. The theory of the MCI was first developed analytically [33][34][35][36], and then using large PIC numerical simulation from 2010 onwards [29,30,39]. At the plasma wave-wave resonant level of description, the MCI essentially involves the resonance of a fast Alfvén wave supported by the background plasma with negative-energy ion cyclotron harmonic waves (terminology described in the next paragraph) sustained by minority fast ions whose non-Maxwellian velocity distribution incorporates a population inversion. This is evident from the structure of the left-hand side of equation (28) of [34] for example. MCI theory was originally developed by [33] for purely perpendicular propagating waves satisfying ω Ω i , the background ion cyclotron frequency, including a ring beam distribution for fast ions. The theory was revisited and extended in [34] to lower frequencies for perpendicular propagation, and MCI growth rates were further obtained for energetic ion distributions in velocity space that have the form of both a spherical and an extended-spherical shell [35], in addition to monoenergetic ring beams [33,40]. The interest in the ring beam-type distribution arose from the subset of fast ions thought to be responsible for generating ICE from DT plasmas in JET and, subsequently, TFTR. These ions, born in a very narrow range of pitch angles, undergo large drift orbit excursions from the core whose trajectories intersect the outer midplane edge. This leads to a local population inversion in velocity space, in the form of a thin cone shape which is limited by: the maximum energy of the α particles; their narrow range of pitch angles; and, at the lower bound, the strong decrease of radial excursion with decreasing energy (figure 15 of [14]).
The negative energy character of the ion cyclotron harmonic waves supported by the energetic ion population, alluded to in the preceding paragraph, is a key elementexplicit or otherwise-in all theories of the MCI. Less importantly, the ion cyclotron harmonic waves supported by the bulk plasma also enter the left-hand-side of equation (28) of [34] and, if resonant, can change the order (quadratic to cubic) of the analytical expression for the growth rate in the linear limit. Neglecting this complication, when the linear MCI occurs, there is spontaneous growth in amplitude of both the fast Alfvén wave supported primarily by the bulk plasma, and the cyclotron harmonic wave supported by the non-Maxwellian alpha-particle population. This does not violate energy conservation because the latter wave is negative-energy: its excitation lowers the overall energy of the system which supports it. Negative-energy waves can only arise where there is a distinctly non-Maxwellian population, for example beam-type, in velocity space; for an account of this, see [41].
The analytical theory of the MCI was extended to include wave-particle cyclotron resonance in [34]. In this form, it was applied to the interpretation of ICE measurements from DT plasmas in JET and TFTR, including sub-Alfvénic regimes in [36,37]. The MCI is also believed to play a central role in the excitation, by energetic ion populations in spherical tokamaks, of compressional Alfvén eigenmodes [42][43][44][45] (CAEs). CAEs can be viewed as a toroidal generalisation of the fast Alfvén wave (or of the magnetosonic wave) at frequencies comparable to, or lower than, Ω H , the energetic ion cyclotron frequency. Hence this phenomenon may in some respects be a lower-frequency continuation of ICE phenomenology. Among the most sophisticated numerical approaches in this context is the HYM code [46], which represents the energetic ion population using a full-orbit delta-f approach, and the background plasma by single-fluid resistive MHD equations. The delta-f method [47] is a particle-in-cell approach for the solution of the perturbed velocity distribution function. With this method, only the deviation of the initial velocity distribution is treated in a PIC-fashion. This has the advantage to significantly reduce the noise in the simulations since only a subset of the velocity space is represented by means of macroparticles.
Finally, we recall that with the introduction of tritium into JET in 1991 [48], ICE was detected at successive cyclotron harmonics of α particles. The intensity of this ICE extended the linear correlation with the measured neutron flux [14] to over six decades of signal intensity across all classes of JET plasma. Also, a striking correlation was obtained between the maximum linear growth rates computed via the MCI linear analytical theory and the time evolution of the ICE amplitude averaged over six TFTR supershots, see figure 6 of [38]. This calculation was carried out using a drifting Maxwellian for the parallel distribution function of the fusion-born alphas and a ring-beam for the perpendicular velocity distribution function. In [49], additional conclusions are drawn. In particular, the aforementioned distribution function gives rise to a sublinear scaling of the linear analytical growth rates (∝ n α /n e ) with increasing alpha density, in a regime where linear theory is not applicable, see figure 7 of [49]. This was taken to indicate that the linear scaling of ICE intensity with fusion alpha density in TFTR suggested that the relative alpha density at the emission location was smaller than 10 −4 . The PIC simulations of ion cyclotron emission have mostly, so far, assumed ring beam distributions or drifting ring beams; this includes the present paper. Such ring beams were used as initial alpha particles velocity distribution of hybrid PIC simulations of the MCI relevant to JET ICE [25], which showed that the linear growth rates scaled sublinearly with the alpha particles density, in agreement with the relevant linear theory (equation (36) of [34]). These simulations further suggested that the spectral power of the nonlinearly saturated fields scaled linearly with the alpha density as observed experimentally in JET.

First principles numerical simulations of ICE
Direct numerical simulations of ICE scenarios were first reported in 2013 [29]. These used a particle-in-cell (PIC) [50][51][52] code [27] (see also section 5) to evolve the full orbit kinetics of millions of thermal ions and electrons, together with the self-consistent electric and magnetic fields, all governed by the Maxwell and Lorentz equations. The distribution of energetic ions in velocity space is initialised to reflect physics considerations relevant to the observations of ICE. PIC simulations motivated by ICE measurements from JET show [25,29] that energetic minority ions relax, under Maxwell-Lorentz dynamics, in ways that replicate the linear MCI at early times and, at later times, produce power spectra capturing measured ICE features. Multiple simulations with different concentrations of energetic ions predict a linear scaling of spectral peak intensity that matches the observed linear scaling of ICE intensity with fusion reactivity [39] in JET. An ICE-related scenario relevant to α-channelling [53] has been proposed on the basis of PIC simulations [30]. It rests on a process that can arise when a radially inward propagating fast Alfvén wave, unstable against the MCI in the outer edge plasma, thereby extracts energy from a fast ion population and transfers it to the bulk plasma.

ICE measurements from NBI-heated LHD hydrogen plasmas
Highly resolved ICE signals were measured in LHD hydrogen plasmas using an ICRF antenna in receiver mode [11] during perpendicular neutral beam injection (NBI) of hydrogen ions. The large antenna loop area (≈600 cm 2 ) enhances the quality of the data, which was recorded at a maximum sampling rate of 5GSas −1 and processed via fast Fourier transform, with a rectangular window of typical duration 100 µs. Examples of time profiles of heating with four proton NBI sources, together with ICE spectra, are shown in figures 1 and 2. These spectra obtained from LHD plasmas 79126 and 79003 are not discussed in 11 and 12, but share key similarities in that they all pertain to hydrogen plasmas heated by neutral beam injection: specifically, 40 keV perpendicular NBI giving rise to ICE, as presented in [11] and [12] and simulated in the present work. These ICE signals discussed in [11,12] are detected shortly after the turn-on of the perpendicular positive-ion based NB injector #4 by an antenna located close to it, at the outer midplane of LHD.
The fundamental ICE frequency f 0 is defined by the measured interval between successive spectral peaks, and also typically corresponds to the frequency of the first measured spectral peak. A linear relation was obtained between f 0 and the magnitude of the magnetic field on axis (at a major radius of 3.6 m) across several LHD plasmas, confirming the cyclotronic character of the detected signal. It follows that the location at which this ICE signal is generated in LHD lies along a magnetic field line on which the proton cyclotron frequency f cH corresponds to the measured fundamental ICE frequency f 0 , see figures 3 and 4 of [11]; this is found to be at both the LHD plasma inner and outer edge in the LHD magnetic configuration (see figure 3 of [11]). The observation of ICE from high density (n e > 5 × 10 20 m −3 ) LHD plasmas, into which NBI cannot penetrate as far as the inner edge, further supports the interpretation that ICE originates from the outer region near the NBI #4 injection point. The ICE signal in the plasmas of [11] disappeared roughly 0.1ms after the turnoff of the perpendicular NBI (see figure 5 of [11]). This synchronization suggests that ICE is driven by the fast injected protons. Particle orbit calculations [11] for the relevant LHD plasma and magnetic field show that NBI protons are lost in a few tens of microseconds, consistent with the observed decay time of ICE.

Direct numerical simulation of LHD ICE using a kinetic PIC-hybrid code
In this paper, we will interpret the ICE spectra shown in figures 1 and 2, in terms of the collective relaxation of the NBI proton population in the LHD plasma edge. To that purpose, we first introduce the PIC-hybrid modelling approach. We then present the physical and computational parameters involved in the calculations. In particular, this study requires, as input, a representation of the distribution in velocity-space of the NBI protons which is based on the kinetic modelling of [54] that we briefly describe. We then move on to present our calculations results.

The PIC-hybrid approach
We use the one spatial dimension and three velocity space dimensions (1D3V) PIC-hybrid code [55,56] approach  implemented in [25]. This follows full gyro-orbit kinetics for ions in collisionless plasmas, where the energetic and thermal populations are represented by hundreds of thousands to hundreds of millions of macro-particles. These interact with each other, and with the self-consistent electric and magnetic fields, through the Lorentz force and Maxwell's equations. The code incorporates all three vector components of the electric and magnetic fields, and of each particle's velocity, and represents the electrons as a massless neutralising fluid. It self-consistently solves and iterates the Lorentz force equation for each particle together with Maxwell's equations in the Darwin approximation [20,55] which neglects the displacement current and alleviates the needs to resolve light waves; it is fully nonlinear. The code resolves ion gyromotion, which is necessary to simulate phenomena such as ICE where key physical length scales and time scales are of the order of the ion gyro-radius (and ion skin depth) and ion cyclotron frequency. In particular, the code captures the full spatial and gyrophase dynamics of resonant particle-field interactions close to the ion cyclotron frequency and its harmonics. We assume quasineutrality: N l=1 Z l n l = n e (1) with n e the number density of electrons and n l , Z l the number density and electric charge of each ion species l. We also have where x denotes distance along the 1D slab geometry spatial domain of our code. We solve for B using Faraday's law while Ampère's law in the Darwin approximation [26] combines with the massless electron momentum equation to give the generalized Ohm's law [57], Here the charge-weighted mean ion velocity V i is defined by where u l denotes the bulk ion velocity of species l. We assume an isothermal pressure law, p e = n e k B T e , with T e the electron temperature. A quiet start [50,58] is used to launch the majority thermal ions in phase space, so as to reduce the noise in our simulations. The code makes use of periodic boundary conditions. This level of approximation is known to sustain fluid waves such as the fast Alfvén and whistler waves, as well as kinetic waves such as electrostatic Bernstein and ion cyclotron modes.

Physical and computational parameter sets
The NBI protons are sub (super)-Alfvénic in the emitting region of LHD hydrogen plasmas 79126 (79003), whose measured ICE spectra are shown in figure 2   The increased number of particles per cell in the sub-Alfvénic calculations reduces the noise in the simulations, and facilitates the calculation of the growth rates presented in Section 6. The simulations at 85° between ω and k use fewer particles per cell and a larger grid length, which enables us to run the simulations for longer as shown in figure 4.
The time step is 0.00025 (0.0001) × τ H , where τ H = 0.143µs (τ H = 0.273µs) is the proton gyroperiod in each simulation. Cyclotron motion is thus highly resolved, in space and time, for the energetic ions whose cyclotron resonant collective relaxation underlies the observed ICE signals. The simulations run for 10-360τ H ; this is determined by the time taken for the instability driven by the NBI ions to satur ate, which for a given set of plasma parameters depends on ξ and on propagation angle. Physical and computational parameters for our simulations are summarized in tables 1 and 2 respectively. When necessary, subcycling for the electric and magn etic fields is used to satisfy the Courant-Friedrichs-Lewy (CFL) condition for the Alfvén wave [47], which states that to ensure the simulations are stable, v max × ∆t < ∆x, v max being the maximum velocity in the calculation. This means for example that a macroparticle can cross at most one cell over a timestep ∆t. In our simulations, energy is conserved to within 0.2%.

Distribution in velocity-space of the NBI protons
Kinetic modelling [54] has previously been used to obtain the steady-state distribution function of NBI fast ions in stellarators, for both TJ-II and LHD. The LHD case focused on hydrogen plasmas heated by 40 keV perpendicular NBI, relevant to our study. The orbit code ISDEP (integrator of stochastic differential equations for plasmas) [54] is a Monte Carlo orbit-following code which solves the Fokker-Planck equation in 5D phase space, namely: (x, y, z, v 2 , λ). The three spatial coordinates (x, y, z) are the guiding centre position, v 2 is a normalised kinetic energy, and λ = v · B/vB is a pitch angle. ISDEP includes collisions of fast ions with background ions and electrons, and treats re-entering particles. The initial NBI ion distribution function is calculated with HFREYA [59] which simulates the evolution of fast neutral particles by modelling their propagation, charge exchange and ionization processes. The resulting distribution in [54]'s work, which is relevant to perpendicular NBI in LHD, presents distinct features. Among them, the velocity distribution in pitch angles λ displays a peak near zero. The calculated velocity distribution function presents a ∼34 keV NB component which is close to the 36.5 keV super-Alfvénic initial fast proton population which, as we shall show, drives the ICE. The distribution is localized toroidally and the related density of fast neutrals increases with ρ, the dimensionless plasma radius defined in terms of toroidal flux surfaces. This NBI proton population is thus expected to arise close to NBI # 4, and therefore close to the ICRF receiver antenna. We shall therefore incorporate this form of minority energetic proton population in our first principles modelling of ICE in LHD.
The former motivates the initialisation of the NBI protons in velocity-space by means of a ring-beam velocity distribution function and the NBI protons are initially uniformly and randomly distributed in gyro-angle. The orientation of the velocity of each NBI proton is thus perpendicular to the background magnetic field, in order to reflect the perpendicular orientation of the relevant NBI in LHD. The inclusion of this NBI ion population in our 1D3V simulations effectively defines the spatial location to which these simulations apply: that is, the location, near the NBI injection point at the outer midplane edge of LHD, where this NBI ion population is present. Figure 4 plots the time evolution, in our simulations, of the change in energy density of the electric and magnetic fields, and of the change in the kinetic energy density of the bulk proton and NBI proton populations. The four panels of figure 4 are for sub-Alfvénic (left) and super-Alfvénic (right) NBI protons, and for propagation angles θ of k with respect to B 0 of 89.5 • (top) and 85 • (bottom). It is evident that the time taken for the NBI fast protons, which are not replenished, to relax, and for the instability which we identify below with the MCI to unfold, saturates on time scales of between 10τ H to 360τ H . This corresponds to a few microseconds to a few 100 microseconds. The relative NBI proton concentration ξ = n NBI /n e is chosen small enough to observe the unfolding of the MCI, yet not too small to reach saturation in a trac- All the simulated spectra in figure 6 show strong excitation at multiple successive proton cyclotron harmonics; and the ICE signal in the simulations is a hundred to a thousand times more intense than the thermal noise. Decreasing the angle between k and B 0 shows preferential excitation at lower harmonics, as well as longer timescales for mode excitations, see figures 4 and 5. Since the power spectra at the bottom panels of figure 6 are averaged over a longer time duration (in particular compared to the proton gyroperiod), high spectral resolution is achieved that translates into the thick curves. The similar intensities on the top right and bottom right panels are due to comparable NBI proton densities. Conversely,  figure 6 (top left) could be due to a numerical artifact during the post processing. This spike could also result from three-wave interaction between oppositely propagating waves such that ω 1 + ω 2 = 0 and k 1 + k 2 = 0. In the PIC-hybrid computations, the simulated emission is most strongly driven for propagation angles that are close to perpendicular to the magn etic field, for which the MCI is most unstable and satur ates quickly as seen from the simulations at 89.5°. Figure 7 shows good qualitative agreement between the spectra generated by relaxation of the NBI ion population in

Results of PIC-hybrid simulations
Bright spots at sequential proton cyclotron harmonics along the fast-Alfvén branch result from the MCI, driven by NBI protons, for waves propagating in the x direction, almost perpendicular to the background magnetic field. The cold plasma dispersion relation for the fast Alfvén wave, whose tangent at the origin satisfies ω = kV A , is shown by the dark dashed lines. The extension of the Alfvén wave below the proton cyclotron frequency arises from the finite parallel wavevector component on the bottom panels. The line ω = kv NBI is shown by the dark trace, and lies above (or below) ω = kV A in the super-Alfvénic (or sub-Alfvénic) NBI regimes. The wedge in (ω, k) space defined by these two lines approximately delineates the region of (ω, k) space where waves can in principle resonate with NBI protons. These plots show that, in our computations, excitation occurs for modes on the fast Alfvén branch and preferentially close to ω = kv NBI . our first principles Maxwell-Lorentz computations with the 1D3V PIC-hybrid code, and the observed ICE spectra from LHD shown in figure 2, in both the sub-Alfvénic and super-Alfvénic NBI regimes in LHD. The ICE signals from LHD studied here were obtained in the same way as the early JET ICE measurements [5,14], that is, using an ICRH antenna in receiver mode. The antenna has broad directionality, and in consequence can receive radiation across a wide range of incident angles. This is the reason for considering the emission at an angle of 85° whose corresponding simulated spectra are also shown at the bottom of figure 7. They show that ICE is still emitted away from pure perpendicular propagation under LHD edge plasma conditions and suggest the robustness of the MCI to drive ICE. In particular we infer that simulations at intermediate angles would generate ICE as well. Although the antenna cannot distinguish between incoming waves with a propagation angle of 89.5° and 85° separately, these would be summed together over the considered range of angles in forming the observed ICE spectrum. Here we compare the measured ICE spectra with our simulated spectra at 85°, which share a similar qualitative rise and then fall in the distribution of power across successive cyclotron harmonic peaks. At 89.5°, there is a monotonic increase of the power with harmonic number for the range of frequencies considered.
In figure 8, we plot the time evolution of the intensity of the spatial FFT of the magnetic field δB 2 z /B 2 0 . There is a oneto-one mapping between the excited k-values and the cyclotron harmonics Ω H , inferred from the dispersion relation of the fast Alfvén wave using information in figure 5. Thus by plotting the time evolution of the distribution of energy across wavenumbers, as in Figure 8, we can identify the time sequence in which specific cyclotron harmonics in the simulated ICE spectrum are excited. A notable feature of figure 8 is the late excitation of the spectral peak at the fundamental cyclotron frequency of the protons in the nonlinear phase. This is of interest because, across multiple simulations, the fundamental cyclotron frequency is often the most strongly linearly stable against the MCI. The observational counterpart of this was noted early on, for example the early experimental comparison of ICE intensity spanning different levels of fusion reactivity in JET plasmas, figure 5 of [14], compares measured intensity at the second harmonic rather than the fundamental. In many simulations, ICE at the fundamental cyclotron frequency is driven by nonlinear beating between pairs of neighbouring higher harmonics that are linearly unstable. Examples include figure 1 of [25] for the MCI linear and nonlinear stability aspect, and figures 3 and 4 of [32] for nonlinear beating between neighbouring cyclotron harmonics in experimental ICE data and PIC simulation output respectively.

Identification of the excitation process in the PIC-hybrid simulations for NBI proton driven ICE in LHD plasmas
The frequencies of the modes excited by ICE in our simulations typically range from ω ≈ 5Ω H to ≈45ΩH. We focus primarily on modes up to ω ≈ 15Ω H , marginally stable and unstable, to be in line with the upper frequency limit of the experimental measurements. These modes are electromagnetic and lie on the fast Alfvén branch, so that ω and k are related by ω ≈ V A k , where V A is the Alfvén speed. In all the calculations that follow, we used the numerical dispersion relation to map k to ω; this relation satisfies the cold plasma dispersion of [60] and is necessary when ω starts to deviate from kV A . In our simulations, we consider waves propagating nearly perpendicular to the local background magnetic field. Such waves can leave an MCF plasma, propagating radially, and be detected beyond it. The simulation outputs encapsulated in figure 8 enable us to infer the rate at which the energy in a given cyclotron harmonic spectral peak grows over time. It is particularly helpful to calculate this during the early phase of growth, because this enables quantitative comparison with counterpart linear growth rates obtained from analytical theory. We denote the early phase growth rate inferred from the simulations at the th harmonic by γ . A primary objective is to quantify the scaling of γ with NBI proton number density. This we shall compare with the corresponding scaling of analytical linear growth rate for the MCI, γ lin ( ), defined [14] by equations 7 and 8 below. Equation 9 shows that γ lin ∼ ξ 1/2 ; hence figures 9 and 10 compare early phase simulation outputs with linear theory by plotting γ versus ξ 1/2 for multiple simulations at a propagation angle of 89.5 • , focusing on = 11 and = 12 for sub-Alfvénic NBI LHD plasma 79126 and super-Alfvénic NBI LHD plasma 79003 respectively. These harmonics are associated with the smallest error bars in our analysis. The agreement shown is good; this further confirms the role of the MCI in our simulations and, by extension, the LHD experiments. Figures 9 and 10 are obtained as follows. The Alfvén dispersion relation provides a one-to-one mapping between the excited ω modes and the excited k modes. In addition, our simulations use an initially uniform density for both the NBI protons and the background plasma, and the domain has periodic boundary conditions. This means there is neither loss of information, nor need for windowing, when taking spatial Fourier transforms of the electric and magnetic fields. We may therefore compute the growth rates of k-modes by taking the spatial Fourier transform of B z (x,t), which we perform for a propagation angle of 89.5 • , leading to B z (k,t) as shown on the top panels in figure 8. One selects an ω-mode at ω = Ω H , 5 15, to which a unique kmode, k = k ΩH is associated through the dispersion relation as in figure 5. The time evolution of B z (k ΩH , t) is then plotted and best fits are constructed to extract the empirical growth rate γ of this mode, as described below. This approach is convenient because it does not require transformations from   the time domain. Once B z (k ΩH , t) is calculated, we identify the interval [t 0 , t 1 ] over which the initial exponential growth phase takes place in our simulations. We find the duration of this initial exponential growth phase in the simulations to be ≈1.0τH for the fastest-growing modes, while the slowestgrowing ones unfold over ≈20τH. We perform multiple fits of B z (k ΩH , t) between [t 1/2 − n∆t, t 1/2 + n∆t] where: t 1/2 is at the centre of [t 0 , t 1 ], which are the start and end times of the initial exponential growth; ∆t ≈ 0.001τ H ; and n varies between 1 and n max , such that [t 1/2 − n max ∆t, t 1/2 + n max ∆t] is the smallest interval to contain [t 0 , t 1 ]. This yields a family of n growth rates γ ,n , 1 n n max for a given mode at ω = Ω H . We take the mean of the individual best fits as the growth rate value γ =γ ,n , and define the associated error ∆γ = σ (γ ,n ), where the bar and sigma respectively represent the average and the standard deviation. The average and variance are taken between n min n n max , where n min satisfies n min ∆t 0.5τ H / . That is, computation of the average starts from ∆t corresponding to half an oscillation of the unstable th mode. This enables us to use the same value of ∆t for each cyclotron harmonic . The procedure described above enables us to calculate growth rates denoted γ for the th harmonic during the early phase of simulations. These are next compared with the scaling of the analytical expressions for the corresponding growth rate γ lin ( ) of the MCI, notably equation (36) of [34] which reads: Here ω p,NBI and ω pi are the plasma frequencies of the NBI protons and of the bulk protons respectively, and v NBI is again the initial velocity of NBI protons. We define where Π x,x , Π x,y and Π y,y are the functions of z NBI = kv NBI /Ω H as given in the appendix of [34]. Near resonance, For a given mode , if all parameters are kept fixed except for the NBI proton density ξ = n NBI /n e , equation (7) yields the scaling γ lin ( ) Ω H = α ξ (9) as in figure 10, where which depends on solely. As stated above, the function χ 0 depends on the dimensionless parameter v NBI /V A as well as on mode number . We have evaluated χ 0 in the super-Alfvénic and sub-Alfvénic cases for cyclotron harmonics between 5 and 15, and obtained roughly constant results of 0.35 and 0.15 respectively. Taking their ratios across these mode numbers range therefore yields a constant value of approximately 2.5, and translates into ratios of α of approximately 2.05. Equation 36 of [34] hence predicts higher growth in the super-Alfvénic regime at constant beam density. The values α are close to 0.30 and to 0.16 in the super-and sub-Alfvénic regimes respectively. This factor of 2 between these two regimes from analytical linear theory differs from the simulations. A possible explanation rests on the fact that analytical linear growth rates in the 1990s literature were calculated on the basis of a first order expansion about the cold plasma linear dispersion relation ω = kV A . By contrast, in our present first principles computations, the real frequency of the is the growth rate inferred from the simulations for mode number , and depends on the relative NBI density ξ at the ICE location. The translated compensated plots represent the quantity (γ /Ω H ) / √ ξ versus ξ, with γ inferred from simulations. If the simulation outputs match the linear theory of the MCI, we expect γ /Ω H = α √ ξ as in equation (9). In this case, it follows that (γ /Ω H ) / √ ξ = α , a quantity that does not depend on ξ, but on the mode value only. This outcome is reflected by the sequence of horizontal lines in the compensated plots. The values of ξ span one order of magnitude, between 10 −4 and 10 −3 . Together, these graphs show that, within modest error bars, γ ∝ √ ξ . excited modes is determined by the self-consistent Maxwell-Lorentz dynamics, and is found to deviate progressively from the linear expression at increasingly high frequencies. Hence the location in (ω, k)-space of the dominant MCI drive differs increasingly from the location assumed in linear theory, at increasingly high cyclotron harmonics. This is visible in figure 4; the extent of any deviation of early-time growth rates from linear theory can only be established reliably, for a given scenario, by running the code. Our code outputs yield continuity between early-time growth rates as one transitions across the boundary between sub-Alfvénic and super-Alfvénic regimes for the particular parameter sets considered, whereas linear analytical theory yields a difference by a factor of order two. We test the hypothesis that γ γ lin ( ) by running multiple simulations for different values of the beam proton density ratio ξ, with all other parameters kept fixed. The computed growth rates at early times γ for a given mode , obtained by means of B ξ (k ΩH , t) as appeared in figure 8, are then plotted against √ ξ , in line with the analytical scaling of γ lin ( ) [5,34] in equation (9) as shown in figure 9. The values of the growth rates obtained in our PIC hybrid calculations are similar between the sub-and super-Alfvénic regimes for comparable NBI proton densities; and the numerically computed values of α are about five for harmonic numbers between 5-15 in both regimes.
More generally, figure 9 shows congruence between the early phase of collective relaxation of the NBI ion population in our first principles PIC-hybrid simulations, and the analytical theory of the linear stage of the MCI. It tends to confirm that the resulting ICE spectra can be understood in terms of the MCI, which our simulations extend into the analytically inaccessible nonlinear regime, for both sub-Alfvénic and super-Alfvénic NBI proton populations.
The procedure is now applied accross multiple harmonic modes with the results shown figure 10.
As in figure 9, in figure 10 we plot the dependence of the growth rate γ of the fields, calculated at early times in multiple simulations, on NBI concentration ξ in each simulation. The computations are performed for parameters corresponding to LHD plasmas 79126 and 79003. In both scenarios, the propagation angle between B 0 and k is 89.5 • . The pair of panels in figure 10 instead reformulates equation (9) as This implies that for a fixed mode number , the quantity (γ lin ( ) /Ω H ) / √ ξ is independent of ξ and equals a constant that depends on solely. We translate that constant such that it equals to α = . Figure 10 strongly suggests that in general γ ∼ ξ 1/2 in our simulations, across an extended range of modes. This dependence is the same as for growth rates from linear analysis of the MCI [34], and from previous simulations e.g. figure 3 of [39]. Additional linear and cubic root scaling tests have been performed, and F-test statistics [61] applied with a 99% significance further confirm the square root scaling.
We have also investigated, through PIC-hybrid simulations, the collective relaxation of NBI proton populations that have artificially enhanced sub-and super-Alfvénic characteristics under the LHD plasma conditions already considered. We have run multiple simulations of a sub-Alfvénic fast proton population of 25 keV, for which v NBI /V A = 0.7 under LHD plasma 79126 conditions, as well as simulations with a 56 keV super-Alfvénic fast proton beam, for which v NBI /V A = 1.4 under LHD 79003 plasma edge conditions. The same conclusions are obtained.

Conclusions
Several advances are reported in this paper. First, we have presented the first 1D3V PIC-hybrid simulations of ICE where the minority energetic ion population arises from NBI, hence with kinetic energy more than an order of magnitude lower than in previous simulations [25,[30][31][32]39] relating to super-Alfvénic fusion-born ions in JET. Second, first principles kinetic simulations of ICE and MCI physics in the sub-Alfvénic regime for energetic ions have not previously been carried out, with one exception. This exception is the set of multiple PIC simulations [31,32] of ICE driven by fusion-born protons under rapidly evolving edge plasma conditions in KSTAR where, for local electron densities below n e ∼ 1.05 × 10 19 m −3 , corresponding to the lower panels of figure 4 of [31], the perpendicular velocity of the protons is sub-Alfvénic. We have found that, in the LHD-relevant context, the transition between sub-Alfvénic and super-Alfvénic ICE phenomenology is extremely smooth, both in observations and simulations.
A third novel aspect of this paper is that it is the first to report first principles kinetic ICE simulations with wavevectors inclined more than one degree from perpendicular to the magnetic field direction. While more challenging computationally, this leads to better congruence of simulation outputs with the observed ICE spectra, as in figure 6, and is helpful in a context where antenna sensitivity may not be known in detail as a function of k in the relevant range. Fourth, we have carried out the first study of early-time growth rates inferred from simulations that span sub-Alfvénic and super-Alfvénic energetic ion regimes. Of particular interest is how these growth rates depend on energetic ion concentration ξ = n NBI /n e ; figure 10 establishes a square root scaling with ξ, in line with prediction from the corresponding linear analytical MCI theory [34]. For simulations relevant to the super-Alfvénic regime, this scaling was established in [39].
In summary, the measured ion cyclotron emission (ICE) spectra (figure 2) from LHD hydrogen plasmas with both sub-Alfvénic and super-Alfvénic perpendicular proton NBI have been successfully simulated (figure 6) using a first principles approach. Direct numerical simulation of kinetic ions (bulk protons and minority energetic NBI protons) and fluid electrons using a 1D3V PIC-hybrid code captures the self-consistent Maxwell-Lorentz dynamics of the plasma and fields. It is found from the Fourier transforms and time evolution of the energy and field components in the PIC-hybrid code outputs that the dominant physical process in our first principles Maxwell-Lorentz computations is the magnetoacoustic cyclotron instability (MCI). The many correlations between our code outputs and the measured ICE spectra suggest that an emission mechanism, which corresponds essentially to the nonlinear MCI, is responsible for the main features of ICE in these LHD stellarator plasmas. In the context of the extensive prior research on the role of the MCI in ICE from tokamak plasmas, this outcome suggests a significant degree of commonality across tokamak and stellarator ICE physics. This appears to be a consequence of the strongly spatially localised character of ICE physics, and implies that overall magnetic geometry plays a relatively minor role. The spontaneously excited electric and magnetic fields in our simulations, which are carried out in local slab geometry, correspond to the fast Alfvén wave. This work helps establish a baseline for future energetic particle experiments in LHD, where magnetic geometry and toroidal eigenfunctions [62] may play a larger role. ICE links beam ion physics in LHD to fusion-born ion physics in tokamaks, and has significant diagnostic potential [20].