High-order coupling of shear and sonic continua in JET plasmas

A recent model coupling the shear-Alfv\'{e}n and acoustic continua, which depends strongly on the equilibrium shaping and on elongation in particular, is employed to explain the properties of Alfv\'{e}nic activity observed on JET plasmas below but close to the typical frequency of toroidicity-induced Alfv\'{e}n eigenmodes (TAEs). The frequency gaps predicted by the model result from high-order harmonics of the geodesic field-line curvature caused by plasma shaping (as opposed to lower-order toroidicity) and give rise to high-order geodesic acoustic eigenmodes (HOGAEs), their frequency value being close to one-half of the TAEs one. The theoretical predictions of HOGAE frequency and radial location are found to be in fair agreement with measurements in JET experiments, including magnetic, reflectometry and soft x-ray data. The stability of the observed HOGAEs is evaluated with the linear hybrid MHD/drift-kinetic code CASTOR-K, taking into account the energetic-ion populations produced by the NBI and ICRH heating systems. Wave-particle resonances, along with drive/damping mechanisms, are also discussed in order to understand the conditions leading to HOGAEs destabilization in JET plasmas.


I. INTRODUCTION
The stability of Alfvén eigenmodes (AEs) is a key issue in magnetically confined fusion, both for currently operating machines (ASDEX-Upgrade, DIII-D, JET, etc) and for next-step devices such as JT60-SA, being one of the major concerns regarding ITER operation [1]. Once destabilized, AEs may degrade the confinement of energetic ions produced by the heating systems (NBI, ICRH) or by fusion reactions in burning plasmas and transport them radially away from the core, compromising in this way the heating process and eventually damaging the device walls or plasma facing components [2]. Unstable AEs are known to arise in frequency gaps of the MHD continuum that are produced by the coupling of Alfvén waves traveling along magnetic-field lines with a periodic refractive index [3]. The specific coupling mechanism (equilibrium shaping, geodesic curvature, pressure) and the dispersion-relation branches involved (shear-Alfvén or slow acoustic, etc) establish some of the AEs properties, like their frequency and dominant polarisation, which are fundamental to understand their interaction with the charged particles confined within the device and the nature of the anomalous transport they may eventually induce.
Toroidicity-induced AEs (TAEs) with frequencies near the reference value ω TAE = ω A /(2q), where ω A = v A /R 0 is the angular Alfvén frequency, v A is the Alfvén speed, R 0 is the torus major radius, and q is the safety factor, are produced by the coupling of two shear-Alfvén (SA) waves and are one of the most extensively studied Alfvénic instabilities in fusion devices due to their potential to cause energetic-ion redistribution and losses [1,4]. On the other hand, unstable lower frequency AEs (i.e., with ω < ω TAE ) were first observed in high-beta plasmas under intense NBI heating in DIII-D [5], in association with high energetic-ion loss levels much similar to those caused by TAEs [6]. Since then, they have been documented in many other devices [7,8]. Some properties of these AEs have been explained by a MHD model coupling a single SA wave with two acoustic waves via the lowestorder harmonic of the field-line geodesic curvature due to toroidicity [9][10][11], which results in two characteristic frequencies, respectively ω BAAE = 2β 1 2 ω TAE and ω GAM = 2β 1 2 1 + 2q 2 ω TAE . (1) The first value corresponds to the top of the beta-induced acoustic AEs (BAAEs) gap, with β = c 2 S /v 2 A the plasma beta and c S the sound speed, whilst the second corresponds to the bottom of the SA continuum branch that is upshifted by finite pressure to a frequency value that is typical of geodesic-acoustic modes (GAMs). Alone or with kinetic corrections, this model has been successfully employed to explain experimental observations of AEs inside the frequency gaps near ω BAAE [12,13] and ω GAM [8].
In this work, we report and discuss Alfvénic activity observed on JET experiments, at about one-half of the frequency ω TAE , in plasmas heated by NBI and ICRH. However, their properties (i.e., typical frequency and radial location) are not easily described by the aforementioned coupling model. In alternative, we show here that such experimental measurements can be explained by a family of previously unexplored continuum gaps that open when several acoustic waves couple with a single SA wave via higher-order harmonics of the field-line geodesic curvature (as opposed to lower-order toroidicity), which have their origin on equilibrium shaping [14]. In particular, the finite elongation of JET plasmas is shown to play an essential role in the argument. Although many families of gaps are predicted by the analytic coupling model [14], some of which have also been found by numerical approaches [11,15,16], here we focus on those gaps that open when acoustic continuum branches are able to cross a SA branch whose bottom has been lifted up to ω GAM . Therefore, the frequencies of interest lie in the range ω GAM ω ω TAE , that is, significantly above the values predicted by previous lower-order coupling models [9][10][11] and closer to the TAEs frequency. Although kinetic corrections are known to change the Alfvén-wave spectrum near the SA-continuum accumulation point [17] and produce kinetic eigenmodes with frequencies above ω GAM [18,19], this work is kept within the framework of ideal MHD, with a perturbative drift-kinetic extension to handle resonant wave-particle interactions.
The outline if this paper is as follows: In section II, some of the lower-frequency AEs most striking features are presented, along with the main characteristics of the experimental scenarios where such observations were made. The limitations of simpler low-order coupling models are discussed, as well as the improved frequency predictions provided by the more complete high-order model. In section III, a specific experimental scenario is modelled numerically with an ideal-MHD code and the AEs predicted frequency and radial location are compared with measurements. These are found to agree with each other and with the analytical estimates of the high- Example continua of SA (ξ A m ) and acoustic (ξ S m ) harmonics of the plasma displacement ξ for q = 1 + 4ρ 2 pol , β = 10 −2 , ε = 0.3, poloidal mode number m = 4 and n = 2: (a) uncoupled cylindrical-plasma limit (dots) and lowest-order coupling for circular equilibria (solid lines); (b) high-order couplings enabled by the elongation value κ = 1.3. order coupling model. In section V, the resonant interaction between low-frequency AEs and relevant plasma species (thermal ions and energetic ions heated by NBI and ICRH) is addressed in its linear stage. For the experimental scenario considered, the AEs stability is found to be mainly established by the interplay between thermal-ion Landau damping and the drive provided by the anisotropic population of ICRH-accelerated energetic ions. The conclusions are summarised in section VI.

II. EXPERIMENTAL OBSERVATIONS AND CONTINUA-COUPLING MODELS
In order to optimise scenarios for the observation of TAEs driven unstable by fusion α-particles, experiments were conducted at JET with deuterium plasmas that are characterised by low density, large core temperature (usually associated with the presence of internal transport barriers) and elevated safety factor [20]. Injection of ICRH power (hydrogen minority heating, concentrations in the range 2% n H /n e 7%, where n e is the electron density) with the fundamental resonance located close to the magnetic axis drives a large number of AEs unstable, as those seen in figure 1.
Using reference values for these plasmas (on-axis field B 0 ≈ 3.4 T, deuterium density n D = n e ≈ 4 × 10 19 m −3 , torus major radius R 0 ≈ 2.9 m, and safety factor q ≈ 1.5, whence the values v A ≈ 8.3 × 10 6 m/s and ω A ≈ 2.9 × 10 6 rad/s) [20], TAEs are expected to be found near ω TAE /(2π) ≈ 150 kHz, this frequency being shifted slightly upwards, according to their toroidal mode number n, due to the Doppler shift caused by finite plasma rotation. These estimates closely describe the TAEs seen in figure 1 around the predicted 150 kHz value. In the same spectrogram, additional activity is also seen at roughly one half of the value corresponding to ω TAE , with well-defined frequency lines surrounded by wide (∼ 15 kHz) frequency bands. As the detailed plot in figure 1 shows, these bands seem to correspond to bursty events lasting a few milliseconds, in close resemblance to those reported in recent DIII-D observations [8]. In contrast, however, the sharp frequency lines between such bands is a distinctive feature of the AEs being reported here. This feature seems to suggest that such AEs arise inside well-defined frequency gaps and do not interact significantly with the continuum nor are substantially affected by other sources of strong damping.
The three-wave model (one SA, two acoustic) [9][10][11], induced by finite toroidicity alone and usually employed to describe continua coupling below ω TAE , is insufficient to explain the AE frequencies in figure 1. As illustrated by the frequency continua plotted in figure 2(a) along the radial variable ρ pol = (Ψ/Ψ b ) 1/2 (with Ψ the poloidalfield flux, Ψ b its value at the plasma boundary, ε = a/R 0 the inverse aspect ratio, and a the torus minor radius) near a single resonant surface, the model predicts a single frequency gap immediately below ω BAAE . However, the latter's relation with ω TAE for a typical value β ≈ 10 −2 sets the upper limit for eventual BAAE frequencies to be around 30 kHz, in clear disagreement with the AEs being measured near 80 kHz. On the other hand, the same three-wave coupling model lifts the bottom of the SA branch to ω GAM ≈ 70 kHz. Besides the difference between estimated and measured frequency values being still significant, it will be shown in the next paragraph that equation (1) indeed overestimates the value ω GAM , which becomes lower if plasma-shaping effects are taken into account. An analytically tractable equilibrium model [21] was recently employed to show that plasma shaping, and in particular the Shafranov shift and elongation, contribute with higher-order harmonics to the geodesic curvature that couples SA and acoustic continua, in addition to the lowest-order term arising from toroidicity [14]. As a consequence, the additional couplings open frequency gaps where the acoustic branches are able to cross the upshifted SA continuum, their frequencies being always in the range ω GAM ω ω TAE , as shown in figure 2(b). Located near a rational surface q = m/n (rather than directly over it), with m the poloidal mode number, the largest frequency gap of these high-order family have its width scaling as β 1/2 (κ 2 − 1), being thus closed for circular plasmas where the elongation is κ = 1. Its frequency is approximately given by [14] ω = 2β The toroidal mode number of the particular AE singled out in figure 1 is measured to be n = 2 by performing a conventional spectral analysis over the data collected by a magnetic-coil array distributed along the toroidal direction [22]. Therefore, equation (2) provides a frequency estimate of 80 kHz, which is in much better agreement with the measurements. Other than opening additional frequency gaps, plasma shaping (particularly elongation) also lowers the bottom of the SA continuum from the value in equation (1) down to the corrected estimate [14] ω GAM = ω BAAE 1 + 2q 2 1 − 3 2 as also depicted in figure 2. A reference elongation value κ ≈ 1.5 produces ω GAM ≈ 51 kHz, further distancing measurements and three-wave model predictions.
One important feature of the coupled SA-acoustic frequency gaps induced by plasma shape is the existence condition [14] q 2 involving the local safety factor and elongation. In short, a q value larger than 2 prevents the gaps from opening, this limit being pushed a bit upwards if κ > 1, which limits their existence to the inner part of the plasma where the safety factor is usually lowest. In certain circumstances, to be discussed later, high-order geodesic acoustic AEs (HOGAEs) can be found within these gaps, which may account for the properties of AEs below ω TAE as those depicted in figure 1.

III. MEASUREMENTS AND MHD ANALYSIS
Unlike the spectrogram in figure 1, where multiple frequency lines are seen within the range of interest (i.e., 60; 100 kHz), JET pulse #90198 depicted in figure 3 provides a simpler scenario for analysis: it shows a clear n = 2 AE below the TAE group, with a single and very sharp frequency line that is well separated from the two bursty bands around it. At 8 seconds, its frequency value is measured to be 83.2 kHz. At the same instant, a neoclassic tearing mode (NTM) with m = 3 and n = 2 is visible at 12.2 kHz, which can be used to estimate the toroidal plasma-rotation frequency Ω φ = 6.1 kHz at the location of the q = 3/2 surface. The time trace in the top of figure 3 shows the AE to arise with increasing ICRH power and on-axis plasma beta, the latter remaining above the reference value 10 −2 while the eigenmode is visible. On-axis elongation stays approximately constant and slightly below 1.4, whereas the safety factor drops steadily in time, being nonetheless always below the critical value 2.8 produced by equation (4). From the same time trace, the AE stabilisation at 8.1 seconds cannot be related with decreasing ICRH power (which provides the drive), plasma beta, or elongation (both necessary to a non-vanishing gap width). Instead, it seems to be caused by the safety-factor drop that moves the continuum outwards, as explained later in this section The magnetic equilibrium is reconstructed with EFIT [23] and constrained by Faraday-rotation data, yielding the safety factor profile depicted in figure 4(a) that places the q = 3/2 surface at ρ pol = 0.24. This plasma rotation datum is also plotted in figure 4(a), along with Charge eXchange Recombination Spectroscopy (CXRS) measurements, thus providing a Ω φ (ρ pol ) profile that is employed to convert between AE frequencies that are computed in the plasma frame and measured Doppler-shifted frequencies. Other equilibrium parameters are the magnetic axis location at R 0 = 2.944 m and the on-axis field B 0 = 3.469 T. Electron density n e (ρ pol ) and temperature T e (ρ pol ) profiles are obtained by fitting high-resolution Thompson scattering (HRTS) and LIDAR data, as displayed in figures 4(b) and (c). Because no thermal-ion density measurements are available, one assumes n D = n e and the plasma mass density being set by the Deuterium ions alone. Likewise, T D = T e is also assumed, following CXRS measurements for ρ pol 0.4. Some variations of these assumptions will be considered, later in section V, in order to take into account Z eff measurements and the uncertainties in the T D profile.
The magnetic equilibrium, after refinement with the Grad-Shafranov solver HELENA [24], and the mass-density profile are the inputs to the compressible ideal-MHD code CASTOR [25] and its continuous-spectrum extension CSCAS [26]. The latter is employed to produce the three continua plotted in figure 5 for the toroidal mode numbers n = 1, 2, and 3, while the former is used to compute the frequency (ω/ω A = 0.173094) and radial structure of an HOGAE inside the high-order gap of the n = 2 continuum. The coupled SA-acoustic nature of the computed HOGAE is illustrated in figure 6. The electric-field perturbation normal to each flux surface is the plasma response to transverse (i.e., bi-normal) displacements and, therefore, shows a strong SA polarization dominated by the same m = 3 harmonic ξ A 3 seen in figure 5. Because this particular coupling is induced by elongation, the dominant acoustic harmonic is ξ S 3+3 [14], which is clearly visible in the poloidal-section structure of the parallel displacement.
From the plots in figure 5, it is clear that only for n = 2 can the HOGAE extend freely towards the magnetic axis without crossing any continuum branch, therefore avoiding significant continuum damping. This is so because HOGAE gaps are located at the right-hand side of the upshifted SA continuum and a line of constant frequency originating there will be blocked on its left-hand side if the corresponding rational surface is too far away from the axis, as the n = 1, 3 cases in figures 5(a,c) clearly illustrate. The same blocking effect takes place if the safety factor decreases in time, as seen in figure 3, and rational surfaces move outwards from the axis. Indeed, the HOGAE plotted in figure 5(b) is just slightly above the continuum near the axis, hinting that any further decrease of the safety factor leads to its stabilisation by continuum damping.
Conversely, if the rational surface is sufficiently close to the axis, as is the case for n = 2 in figure 5(b), the continuum left-hand side stays beyond ρ pol = 0 and an eigenmode can be found through the gap. In view of this property, experimentally observed HOGAEs with given toroidal and poloidal mode numbers n and m can be regarded as proxies for the existence of a q = m/n rational surface close to the axis, which may be used to place boundaries on the value q(0).
The peculiar shape of the HOGAEs displacement harmonics, which tend to concentrate towards the axis, makes them weakly sensitive to eventual crossings with the continuum outwards from the gap radial location, which for the n = 2 case in figure 5(b) is seen to lie at ρ pol = 0.44. Using the values for R 0 and B 0 returned by the equilibrium calculation and the on-axis Deuterium density estimated from the fit to n e data in figure 4(b), one gets the Alfvén frequency on axis ω A /(2π) = 426 kHz. After correcting for the Doppler shift due to Ω φ /(2π) = 4.5 kHz inferred at the gap radial location, the HOGAE predicted frequency becomes thus 82.7 kHz, which is fairly close to the measured value that was reported at the beginning of this section.
Many sources of experimental uncertainty may contribute to the difference between predicted and measured frequency values, the most important ones being con-cerned with the safety-factor and mass-density profiles. While the former is closely related with the reliability of equilibrium reconstruction performed by EFIT, the latter is connected not only with eventual HRTS/LIDAR uncertainties but also with the assessment of the several impurity fractions.

IV. EIGENMODE RADIAL LOCATION AND EXPERIMENTAL CONSTRAINTS
For the two JET pulses #90198 and #90199 considered in this work, there are no reflectometry measurements available from the edge down to the axis in order to unambiguously provide the AE location. As seen in figure 7 for the pulse #90198 case, the reflectometer cut-off layer is always restricted to the outer half of the plasma (i.e., 0.5 < ρ pol < 1). However, quite significantly, no traces of the HOGAE can be found in the spectrogram of the reflectometer signal (top) when it is compared with the one corresponding to the magnetic coil (bottom). Therefore, the eigenmode must be located in the inner half of the plasma (ρ pol < 0.5), which agrees with the numerical solution computed by CASTOR and plotted in figure 5(b).
Although the HOGAE perturbation is too weak to show up clearly in soft x-rays (SXR) signals, it is possible to find statistically significant phase coherence between the latter and magnetic-coil signals for a certain number of SXR channels as depicted in figure 8. Indeed, all lines-of-sight that cross magnetic surfaces with ρ pol < 0.4 show some trace of the perturbation, similar to the example provided (right-hand side, middle panel). On the contrary, all lines-of-sight lying entirely at ρ pol > 0.45 do not show any traces (right-hand side, bottom panel). Again, these experimental constraints support the HOGAE radial structure depicted in figure 5(b).

V. EIGENMODE STABILITY AND RESONANT INTERACTIONS
The resonant energy exchange between HOGAEs and an ion species s is evaluated by the hybrid MHD/driftkinetic code CASTOR-K [27,28] following a perturbative approach, which returns the linear growth rate as Here, L (1) and f (1) s are the linear response of the guidingcenter Lagrangian and equilibrium distribution function f (0) s to the MHD displacement eigenfunction ξ computed by CASTOR. The integrals in equation (5) are evaluated in the space of the guiding-center constants of motion: energy E, toroidal momentum P φ , and Λ = µB 0 /E, with µ the magnetic moment [29].
For such purpose, thermal-bulk Deuterium ions are described by a local Maxwellian, their density and temperature being, for the moment, those in figure 4 (i.e., T th-D = T i = T e and n th-D = n i = n e , which integrates radially to yield N e = 1.97 × 10 21 electrons within the plasma volume). The distribution of NBI-heated Deuterons is computed by the heating code ASCOT [30]. In turn, the code PION [31,32] provides the distributions of ICRH-heated H and D ions, their density and temperature being plotted in figure 9. The strong anisotropy of the hot H ion population is modelled by a separable distribution where 2R 0 δ ICRH is the Doppler broadening of the ICRH resonant layer centred around the vertical line R ICRH /R 0 = Λ ICRH . In the following, the fundamental ICRH resonance of H ions is assumed to be located on axis, whence Λ ICRH = 1. On the other hand, ICRHheated Deuterons are assumed to be already thermalised, being thus also described by a local Maxwellian. The damping rates computed by CASTOR-K for the thermal and NBI Deuterium ions are γ th-D /ω = −0.015 and γ NBI /ω = −0.002, respectively. ICRH Deuterons have a negligible damping contribution (of order 10 −4 ) due to their modest number (N ICRH-D /N e ∼ 10 −2 ) and temperature (T ICRH-D ∼ 50 keV). In turn, the drive from ICRH-heated H ions depends strongly on the gradient ICRH-H /δ ICRH , which is related with the distribution-function anisotropy, and thus with the normalised Doppler broadening δ ICRH of the resonant layer. due to a 0.5% fraction of ICRH-heated H ions (as computed by PION with respect to the electron number) is not sufficient to surpass the combined damping and drive the HOGAE unstable, except for the lowest value of δ ICRH . Besides uncertainties in the fraction of hot H ions (on which the driving rate depends linearly), other effects may contribute to a different balance of the energy exchange.
One such effect is the thermal-ion dilution due to the many impurities (in unknown concentrations) that are responsible to raise the effective charge number to the measured value Z eff ∼ 2.4. To keep the analysis simple, an impurity mix of 3% Beryllium and 0.1% of Neon and Nickel is considered, along with 4% of minority H, which keeps the quasi-neutrality and Z eff consistent with measurements, while diluting the thermal Deuterium down to 80%. The latter must be further reduced by 6%, in order to account for the 3% of Deuterons heated by ICRH (as computed by PION) and 3% by NBI (as reported by ASCOT). Overall, one finds n th-D = 0.74n e and the thermal ion damping reduces accordingly, raising the critical value of the resonance-layer width below which the HOGAE becomes unstable. Another effect is caused by the uncertainty of the profile T th-D (ρ pol ), which is constrained by CXRS data only in the outer half of the plasma (recall figure 4). The impact of reducing the ratio T th-D /T e on the combined damping rate is depicted in figure 10, where the critical resonance-layer width is seen to increase further. Additional sources of damping, like continuum or radiative damping, as well as the eventual impact of other non-ideal-MHD effects, were not taken into account. The major wave-particle resonances are plotted in figure 11 for thermal Deuterons and ICRH-heated H ions. The interactions of the former are dominated by strongly passing particles (Λ ∼ 0), whose stronger resonances are located at 5 keV and 10 keV. These values are close to the plasma temperature at the AE location but significantly below the fundamental resonance at v ∼ 3c S [14], which corresponds approximately to 50 keV in this case.
In turn, resonances with ICRH hot H ions are stronger at 340 keV, with two slightly less intense peaks near 120 keV and 60 keV, all close to the H temperature provided by PION (i.e., about 150 keV).

VI. DISCUSSION
In summary, Alfvénic activity observed on JET experiments below but close to the typical TAE frequency was reported. Its properties were explained using a model coupling the SA and acoustic continua that depends strongly on equilibrium shaping and, in particular, on the elongation of JET plasmas. Predicted frequencies and radial locations of the HOGAEs arising from such coupling in a specific JET discharge were found to be in fair agreement with experimental measurements, which include magnetic, reflectometry, and SXR data. The stability assessment confirms HOGAEs instability in the presence of ICRH-heated ions with energies around 340 keV, if the resonance layer is sufficiently narrow or, correspondingly, if their distribution-function anisotropy is sufficiently high. The effects of thermal-ion dilution by impurities and different T i /T e ratios on the stability threshold were also addressed.
Numerical results indicate that anisotropy is a good candidate to drive HOGAEs unstable in the JET scenar-ios under analysis. Indeed, figure 10 hints that isotropic ICRH distributions (i.e., in the limit δ ICRH → ∞) are unable to surpass thermal-ion Landau damping unless T i /T e is unrealistically low. However, this conclusion applies to the specific ICRH population considered in this work, which is characterised by temperatures below 200 keV at the eigenmode location. The effect of isotropic populations with much higher energy (e.g., fusion alphas) on HOGAEs stability is still an open question.