Predicting ion cyclotron emission from neutral beam heated plasmas in Wendelstein7-X stellarator

Ion cyclotron emission (ICE) is spontaneously generated strongly suprathermal radiation peaking spectrally at multiple harmonics of the cyclotron frequency of an energetic ion population within a plasma. The narrowband character of the ICE spectral peaks implies that the source is strongly spatially localised. ICE spectra are widely observed in magnetic confinement fusion (MCF) plasmas, and the peak intensities are usually a few orders of magnitude stronger than the thermal radiation.

ICE was early observed from fusion-born ion populations in pure deuterium and DT plasmas in JET [1–4] and TFTR [5, 6], and from neutral beam injected (NBI) ions in TFTR [7]. Since 2017, ICE has been detected and analysed from NBI and fusion-born ions in the KSTAR [8–10], DIII-D [11–16], ASDEX-Upgrade [17–21], NSTX-U [22, 23], JT-60U [24–26], TUMAN-3M [27–29], EAST [30–33] and HL-2A [34, 35] tokamaks, and the LHD heliotron-stellarator [36–41]. The magnetoacoustic cyclotron instability (MCI) has been identified as the excitation mechanism for ICE through analytical studies [6, 7, 42] and simulations [8–10, 12, 38, 39, 43–46]. The MCI can arise at spatial locations where the velocity space distribution of an energetic ion population has a positive slope, $\frac{\partial f}{\partial v_\perp}\gt0$ for values of $v_\perp \gtrsim V_A$, in the velocity space. Under these conditions, waves propagating nearly perpendicular to the magnetic field are excited on the fast Alfvén-cyclotron harmonic branch through wave-particle cyclotron resonance, with the fast ions transferring some of their energy to the excited waves. The spatial location of the population inversion can be inferred by matching the frequency separation between successive ICE spectral peaks to the local cyclotron frequency of the energetic ions, and hence to the local magnetic field strength and the corresponding radial location. The energetic ion species whose free energy has been found to drive ICE through the MCI include fusion born ions [1–3, 5, 8, 10, 14, 24, 25, 47–51], neutral beam injected(NBI) ions [5, 7, 9, 11–14, 19, 24, 38, 48, 49] and ions heated via ion cyclotron resonant heating [52, 53].

In this paper, we focus on the likely prospects, physics basis and spectral character for ICE associated with NBI of energetic ions in future Wendelstein 7-X stellarator plasmas. As noted above, ICE due to NBI ion populations has been previously detected from several tokamaks and the LHD heliotron-stellarator. Comparison between ICE spectra from tokamaks and stellarators sheds light on the relative importance of overall magnetic field structure compared to spatially localised physics. Edge NBI ICE has been simulated in [9, 12, 38, 39] from first principles, using particle-in-cell (PIC) kinetic codes which solve the Maxwell–Lorentz system of equations self-consistently for tens of millions of gyro-orbit resolved particles. We note that simulations [38, 39] for the LHD stellarator, unlike those presented here, were hybrid simulations with fluid electrons and did not fully resolve electron gyrodynamics. Here we report the results of PIC simulations using the EPOCH code [54] with fully gyro-resolved thermal ions, energetic ions and electrons. We use a locally uniform (slab geometry) treatment of the plasma physics of ICE emission in these simulations, as in the numerous previous PIC simulations for a wide range of MCF plasmas. This local treatment does not incorporate toroidal magnetic field geometry, whether tokamak or stellarator. Its past success in matching observed ICE signals suggests that these spectra are generated by the spatially localised physics that is captured in our code. In particular, the recent success in matching ICE simulations to ICE observations from the LHD heliotron-stellarator [38–40] provides support for adopting here our predictive approach for future ICE physics in NBI-heated plasmas in the W7-X stellarator. These simulations are computationally resource-intensive, partly due to the low ratio of the perpendicular velocity of the NBI ions to the local Alfvén velocity, $V_\textrm{NBI}/V_A = 0.14$. The sub-Alfvénic regime of ICE, though uncommon, has been previously observed in TFTR [5], ASDEX-U [19] and LHD [39]. Our simulations capture the full frequency range from ion cyclotron through lower hybrid and beyond, with high resolution. It appears that both the MCI and the lower hybrid drift instability, found in related simulations [55, 56], may operate simultaneously under the Wendelstein 7-X plasma conditions. It appears relevant that turbulence in the lower hybrid frequency range has been observed from NBI-heated plasmas in W7-X [57] and LHD [41]. Also important, and explored here, is the extent to which NBI-driven ICE phenomenology is predominantly electromagnetic or electrostatic. In the stellarator context, we refer again to [38, 39] on this topic, together with the PIC studies in [58, 59] which also bear on lower hybrid aspects. The latter papers are in the same spirit as the present work, and differ insofar as in [58] the higher-frequency excitation is electromagnetic as distinct from electrostatic, and in [59] the ion cyclotron and lower hybrid frequencies are substantially closer to each other. These two instabilities, which are taken to be well separated in frequency and wavenumber in a linear analytical description, and hence decoupled, appear here as manifestations of a single coupled plasma phenomenon at the level of first principles self-consistent nonlinear kinetics. The development of a predictive capability, in addition to interpretive, for linking the spectral structure of ICE to the velocity-space structure and spatial location of the emitting ion population is important for the diagnostic exploitation of ICE in present and future [60, 61] fusion experiments.

It is widely accepted that ICE from MCF plasmas is driven by the MCI, whereby a minority non-Maxwellian population of energetic ions excites modes on the fast Alfven-cyclotron harmonic wave branches, transferring some of its energy also to the kinetic energy of the thermal ions and electrons [42, 62–66]. Here we study the MCI under conditions relevant to NBI populations (freshly ionised, or undergoing prompt loss) in the W7-X edge plasma. This location is the most likely one for excitation of NBI-driven ICE in W7-X, as previously observed in the LHD heliotron-stellarator [38, 49] and the TFTR [5, 7], KSTAR [9], DIII-D [11–14] and ASDEX-Upgrade [21, 48] tokamaks, for example. The reason is that, near the NBI injection point in the edge plasma, the velocity-space distribution of freshly ionised NBI ions, which have not yet slowed down through collisions, is closest to a delta-function; and hence most strongly unstable against the MCI which gives rise to the observed ICE signal. We adopt the first principles approach previously used to interpret observations of NBI-driven ICE from other MCF plasmas [9, 12, 38, 39]. This approach involves the use of the PIC code [67, 68] EPOCH [54], which captures the fully gyro-resolved Lorentz dynamics of all particle species, and evolves these dynamics self-consistently with the electric and magnetic fields by solving the full Maxwell's equations. We use EPOCH in 1D3V mode, that is, the simulation domain incorporates one spatial axis, which can be oriented at any angle to the direction of the ambient magnetic field, and all three vector components in velocity-space. It is known that the 1D3V approach, which demands substantial but not unrealistic computational resources, provides a good description of MCI physics for ICE simulation and interpretation [9, 12, 38, 39, 43]. Hence one can be moderately confident in now deploying this capability in predictive mode for W7-X.

EPOCH self consistently evolves the Maxwell equations using a finite-difference time-domain method for all three components of electric and magnetic fields defined on a Yee staggered grid and the gyro orbit resolved charged particle dynamics using a modified version of the Boris algorithm. The charged particles are represented by macroparticles. To reduce noise in the system, we use a fifth order weighing scheme. The electrons and majority thermal ion(protons) populations form the maxwell distributed pseudorandom thermal background. The minority fast ions representing the NBI protons are initialised as a drifting ring-beam distribution, originally modelled with a delta function for the perpendicular component in [42, 65] and later investigated with a spread in the perpendicular component in [69], given by

$$\begin{align} f\left(v_\parallel,v_\perp\right)\propto \textrm{exp}\left(-\frac{\left(v_\parallel-v_d\right)^2}{v_r^2}\right)\textrm{exp}\left(-\frac{\left(v_\perp-u_0\right)^2}{u_r^2}\right) \end{align} \tag{ 1 }$$

where $v_\perp$ and $v_\parallel$ are the velocity components perpendicular and parallel to the background magnetic field, u₀ is initial perpendicular velocity and v_d is initial drift along the background magnetic field, u_r and v_r are the perpendicular and parallel velocity spreads respectively. The system is then allowed to relax.

We ran several simulations using the physical and computational approach described above with one higher (by a factor of 8) resolution simulation used to benchmark the outputs of simulations at lower resolution. For all simulations except the higher resolution simulation, we used 150 000 cells with the time step 0.95 times that required by the CFL condition. The higher resolution simulation used 1200 000 cells. PIC simulations typically require that the cell size is less than the electron Debye length $\lambda_{D_e}$ thereby resolving the associated plasma physics. The Nyquist condition in wavenumber space requires that dk, the resolution in wavenumber (k) space, is smaller than $\omega_c/v_A$ where ω_c is the ion cyclotron frequency and v_A is the Alfvén velocity. To pre-empt numerical heating of electrons, the cell size must resolve the smallest gyroradius in the system which is the electron gyroradius. The spatial domain should also be long enough to have at least several gyroradii of all species. The data from the simulation cannot be saved at every timestep due to memory limits, but but must nevertheless be stored at sufficiently high cadence to resolve the physics being explored here. This work requires the resolution in angular frequency, $d\omega$, to be smaller than the ion cyclotron frequencies as the MCI is a cyclotron resonant wave-particle phenomenon, hence the time duration of our simulations must exceed several ion gyroperiods. The simulation Nyquist frequency must exceed the largest frequency present in the system, in our case the electron cyclotron frequency, to avoid aliasing or adding noise in the lower frequencies. The noise is also inversely dependent on the number of macroparticles used to represent the charged particles. Here we use 1995 000 macro particles for each species. In these simulations, the spatial domain, and hence the allowable wavenumber vectors, was oriented quasiperpendicular (85–89.5 degrees) to the direction of the magnetic field. Other initial parameters are listed in table 1. In the simulations presented here, we find (see e.g. figure 6) that nonlinear three-wave coupling plays a significant role in the flow of energy in the frequency-wavenumber space, and hence in determining the structure of the simulated ICE spectra. This nonlinear coupling is quantified in section 3 using bispectral analysis techniques which are outlined briefly in appendix and has previously been used to assist interpretation of ICE simulations in, for example, figure 3 of [8].

Table 1. Physical Parameters used in the simulations.

Parameter	Value
Electron density (ne )	$6\times 10^{18}$ m−3
Thermal electron and proton temperature	30 eV
Magnetic field	2.25 T
NBI proton beam energy	55 keV
NBI proton density	$0.03 \times n_e$
Beam injection angle w.r.t. B	70°
Initial drift along B (vd )	$1.11 \times 10^6$ m s−1
Initial perpendicular velocity (u0)	$3.05 \times 10^6$ m s−1
Perpendicular velocity spread (ur )	$0.001 \times u_0$
Parallel velocity spread (vr )	$0.01 \times v_d$

Our simulation results are encapsulated in the following six figures, which we discuss in detail below. To summarise these briefly, figure 1 displays the time evolution of the change in energy densities of particles and fields. Figure 2 shows the spatiotemporal Fourier transform of the energy associated with the z-component of the excited magnetic, that is, the distribution of power across different frequencies and wavenumbers. Figure 3 plots simulated ICE power spectra, derived by summing spatiotemporal Fourier transforms over wavenumber, for a range of relevant parameters. We compare the power excited in electric versus magnetic field components in figure 4. The time evolution of energy with respect to wavenumber is shown in figure 5. Nonlinear coupling of energy across wavenumber space is quantified in the bispectral plots of figure 6, for both electric and magnetic fields, at different propagation angles. Figure 1 shows the time evolution of particle and field energy densities from a 1D3V EPOCH PIC simulation of NBI ions relaxing under the MCI for W7-X edge plasma conditions. We emphasise that this reflects spontaneous collective relaxation of millions of simulation particles—NBI ions, thermal ions and electrons—together with the self-consistent electric and magnetic field under the Lorentz-Maxwell system of equations; the MCI is not present in the simulation initial conditions and emerges as the simulation evolves.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

MCI involves a collective relaxation of non-Maxwellian ion population, i.e. the beam ions in our simulations. This is evidenced in figure 1 which shows the change in energy densities in the field components E_x and B_z and in the different particle populations. In this sub-Alfvénic regime for the NBI ions, it is expected from previous analyses of observations of NBI-driven ICE [7, 9, 38, 39] that the MCI may include a significant electrostatic component. This is seen here in figure 1, where the initial rise is at times corresponding to the growth of field energy at $\omega \approx \omega_{LH}$(cf figures 2 and 5). This transitions at $t = 1.25 \tau_{gp}$ in figure 1(cf figure 5, below) to the growth of spectral peaks at lower ion cyclotron harmonics. The system settles into its saturated regime, which presumably corresponds to the state determining the observable ICE spectrum, well before the simulation ends at three ion gyroperiods. In figure 1, the time traces are flat after two gyroperiods. As we shall see later, when resolved with respect to wavenumber, there is correspondingly little or no time evolution in the upper third of figure 5 below, covering the time interval between two and three gyroperiods.

Magnetic field data from the simulation is used to calculate the spatiotemporal Fourier transform (2D discrete Fast Fourier Transform), as shown in figure 2(left) for the case with wave propagation at an angle of 89.5° with respect to the magnetic field B. The B_z component of the magnetic field is chosen for the analysis, because it is perpendicular to the simulation direction and hence the k-vector, thus capturing the electromagnetic physics of the system. Figure 2(left) shows that the spontaneously excited fields are concentrated in distinct regions of (ω,k) space: spectrally peaked at the second ion cyclotron harmonic; and an extended lower hybrid feature between dimensionless wavenumbers k = 1030, and k = 50–60 (hereafter the 'ELHF') close to the lower hybrid resonance frequency given by equation (5) in [70]

$$\begin{align} \omega^2 = \omega^2_\textrm{LH}\left( 1+\frac{m_i}{m_e}\textrm{cos}^2\theta\right). \end{align} \tag{ 2 }$$

Here m_i , m_e are the ion and electron masses, θ is the angle of propagation and ω_LH is the cold plasma lower hybrid frequency given by equation (1) in [70]. Using equation (2) gives 14.96, 25.21, 34.41 and 54.19, normalised to ω_c , for propagation angles $89.5^\circ, 88^\circ, 87^\circ$ and 85° respectively. The ELHFs can be observed around these frequencies also in the spatiotemporal Fourier transforms displayed in figure 6(below), where they are paired with their corresponding bispectra (discussed below). ELHF might be related to the lower hybrid wave physics studied in [55, 71].

Integrating the squared spatiotemporal FFT of B_z over k-space and multiplying by the constant $\frac{1}{2\mu_0}$ yields the frequency power spectrum, figure 2(right). Although the exact values of the spectral peak in an experiment can vary depending on various factors such as measuring device, technique, and position of antennae, the relative heights of peaks should follow the same trend in these logarithmic plots. We note that a difference of α in the log scale (sometimes referred to as 'prominence') indicates that the peak has $10^\alpha$ times the power compared to the base level. Multiple simulations (figure 3) for quasi-perpendicular propagation at angles between 85° and 89.5° to the magnetic field show that the excitation of spectral peaks at the 2nd and 3rd harmonics is a robust feature. Their amplitude is at least three orders of magnitude above the noise floor, which is plotted in the upper two panels of figure 3, and these peaks are consistently present (while others are consistently absent) in multiple simulation cases. Comparison with the noise floor (which possesses spectral structure in consequence of the fluctuation dissipation theorem, whereby noise energy is concentrated at normal modes) also shows that the simulated ICE spectrum is strongly suprathermal. This aligns with ICE observations in all experiments, and has been a prime driver of ICE interpretation in terms of collective instability since the earliest measurements, see e.g. [1]. This is a key prediction of the present work, reinforced by the overall invariance of the phenomenology with respect to the modest variations in electron temperature shown in figure 3(b) and in velocity-space spread shown in figure 3(c). The 1D3V EPOCH PIC simulations reported here follow all three vector components of both the electric and the magnetic fields, which are evaluated at each timestep at each grid cell. This enables us to quantify the balance of the simulated ICE spectrum with respect to electrostatic and electromagnetic components, information which is important also from a diagnostic perspective. For example, the field component E_x lies along the spatial domain of the simulation, and is purely electrostatic; whereas the oscillating component of B_z is purely electromagnetic. The frequency power spectrum for the electrostatic component E_x is obtained in the same way as for the electromagnetic component B_z described above, but the multiplicative constant is replaced by $\frac{\epsilon_0}{2}$. The two panels of figure 4 shows examples of these, together with a computation of the ratio of the magnitude of the electrostatic to electromagnetic components of this simulated ICE spectrum. It can be seen that the electrostatic component is dominant across most of the frequency domain, whereas the electromagnetic component is dominant around the 2nd cyclotron harmonic peak. The relative differences in the electric and magnetic field spectra also indicate which probe might be more suitable to detect the signal at different frequency ranges on W7-X. Since the current B-dot probe ICE diagnostic [73] operates primarily inductively, it is unlikely that the predominantly electrostatic high frequency waves which comprise the ELHF will be detected by it. They should however be detectable in the ICRF antenna which will be introduced in the next experimental campaign, OP2.2. The variation in nonlinear energy coupling between the electromagnetic and electrostatic channels is explored further in figure 6 using bispectral analysis.

The 1D FFTs of magnetic field component B_z in the spatial domain at a given time show the power in different wavenumbers(k) in the system at that instant. Horizontally stacking these k-space spectra taken at successively later times in the simulation yields us figure 5. The ELHF structure can be observed growing first, before the ICE at lower harmonics is initiated, and then dying out around the time when the lower harmonics gain energy. This suggests the ELHF may drive the ICE in the lower harmonics, with the energy transferred by the nonlinear coupling, which has been tested using bispectral methods as described below.

Figure 6 shows the dispersion relations and the bicoherence (see appendix) of the magnetic field B_z and electric field E_x signals. All the dispersion plots are qualitatively similar with the ELHF located around the lower hybrid frequency and significant power at the lower ICE harmonics. The bicoherence plots have colourmaps ranging from 0 to 1 since bicoherence is a normalised bispectrum, see again in the appendix. This figure only informs about the strength of nonlinear coupling between each set of three phase-matched waves, but not the direction of energy transfer among them. The wavenumbers from the bicoherence can be matched to the dispersion plot to find the corresponding frequencies. The change in energy in different wavenumbers from figure 5 indicates which modes gain or lose energy. It can be seen that the ELHF is nonlinearly coupled to the lower and middle harmonics for all cases, sometimes indirectly through a third frequency. The bicoherence of the B_z signal is similar, but not identical, to that of E_x . For example, there is nonlinear wave coupling in a small region of dimensionless k-space, which is visible in the bicoherence of B_z but not in that of E_x , around $({\sim}58,{\sim}2)$for the 89.5° cases, around $({\sim}85,{\sim}5)$ for the 88° case, and around $({\sim}115,{\sim}5)$ for the 87° case. These differences indicate the importance of self consistently evolving and exploring both electrostatic and electromagnetic components, as here, when simulating ICE phenomenology. Faint diagonal lines which intersect the ELHF are visible in some dispersion plots. These features, seen also in [58], are a factor of order one thousand weaker than the dominant features of the spatiotemporal Fourier transform, hence are close to the thermal noise floor in these simulations. While it cannot yet be established whether they are numerical or physical in origin, it seems likely that they are an effect, rather than a cause, of the ELHF; they may repay investigation in future studies.

We predict the excitation of ICE spectral peaks at the 2nd and 3rd proton cyclotron harmonics in the edge region of hydrogen plasmas heated by proton NBI in the W7-X stellarator. The predominantly electrostatic character that we predict for these ICE peaks suggests that they could be detectable with probes, as previously in TFTR [5, 7], as well as antennas. The lower cyclotron harmonics in our simulation may be driven by nonlinear interactions within the ELHF structure, which is excited at frequencies closer to the lower hybrid frequency. The excitation of the ELHF, involving radiation at frequencies comparable to the lower hybrid frequency, is our second main prediction. Based on the change in power in the wavenumbers with time and the nonlinear coupling revealed by the bispectral analysis of the fields excited in our simulations, cross-scale energy transfer is expected between the ELHF and the lower frequency ion cyclotron harmonic peaks characteristic of ICE. Whilst these predictions are specialised to the W7-X NBI edge plasma scenario specified, we note that the nonlinear coupling between lower hybrid waves, identified here, could also arise in other MCF plasmas irrespective of their overall confining magnetic field geometry, and indeed in solar-terrestrial and astrophysical plasmas where ICE-type phenomena have been observed [74–77] or are anticipated [78, 79]. In summary, the phenomenology predicted by these simulations is not narrowly concentrated in frequency and wavenumber. Instead, energy flows initially into wavenumbers corresponding to frequencies near the lower hybrid frequency, the broadband ELHF. It appears from figures 5 and 6 that, probably, energy then flows by nonlinear self-interactions into: wavenumbers corresponding to the intermediate frequency range, which would also be expected to be excited by standard MCI and ICE physics; and to the lowest ion cyclotron harmonic peaks.

Calculations were performed using the Sulis Tier 2 HPC platform hosted by the Scientific Computing Research Technology Platform at the University of Warwick. Sulis is funded by EPSRC Grant EP/T022108/1 and the HPC Midlands+ consortium.

This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200—EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.

S C C acknowledges support via the Johannes Geiss Fellowship of the ISSI, and AFOSR Grant FA8655-22-1-7056.

We use bispectral analysis to identify nonlinear phase coupling, which is indicative of energy being transferred between waves. In particular, we indentify where three waves interact nonlinearly while satisfying the wavenumber and (angular) frequency matching conditions $k_3 = k_1+k_2$ and $\omega_3 = \omega_1+\omega_2$; i.e. where waves '1' and '2' interact to produce '3'. We use a Fourier transform based bispectrum b_s , defined as [80]

$$\begin{align} b^2_s\left(k_1,k_2\right) = |B\left(k_1,k_2\right)|^2 = |\langle F\left(k_1\right)F\left(k_2\right)F^*\left(k_1+k_2\right)\rangle |^2 \end{align} \tag{ 3 }$$

where $F(k_1)$ is the complex Fourier transform of the quantity being analysed, in our case the electric field E_x or magnetic field B_z . The angled brackets $\langle.\rangle$ can represent averaging over time or ensemble. In our case, we implement ensemble averaging by creating multiple records of the signal by using a sliding window of the signal with respect to time. Computing this for the entire space is expensive, in terms of computation time and memory. Therefore, we calculate the bicoherence around the regions with maximum power in the spatiotemporal Fourier transform. The bispectrum can be normalised in a number of ways [81, 82] to obtain the bicoherence, b_c , and we choose to do so using Schwartz's inequality:

$$\begin{align} b^2_c\left(k_1,k_2\right) = \frac{|\langle F\left(k_1\right)F\left(k_2\right)F^*\left(k_1+k_2\right)\rangle|^2}{\langle |F\left(k_1\right)F\left(k_2\right)|^2 \rangle \langle|F^*\left(k_1+k_2\right)|^2\rangle} \end{align} \tag{ 4 }$$

such that the bicoherence is bounded by $0\unicode{x2A7D} b_c \unicode{x2A7D} 1$. A value closer to 1 indicates there is coupling between the two frequencies and their sum while a value close to 0 indicates its absence. For example, if there is an intense region around (a, b) then waves with wavenumbers a, b and a + b are coupled. The value of $b_c^2$ measures the fraction of power at a + b that is coupled in the three wave coupling. The preceding definitions are in terms of wavenumber(k) which is relevant to our case, but the same definitions and relations hold in (angular) frequency space and k can be substituted by ω. We note that equation (4) posses the following symmetries: $B(k_1,k_2) = B(k_2,k_1) = B^*(-k_1,-k_2)$ and $B(k_1,k_2) = B(-k_1-k_2,k_2) = B(k_1,-k_1-k_2)$ [80]. These symmetries apply equally to the bicoherence in three wave interactions [83].