Axisymmetric global Alfvén eigenmodes within the ellipticity-induced frequency gap in the Joint European Torus

Global Alfvén eigenmodes (GAEs) are normal modes of plasma oscillation that exist in a cylindrical magnetised plasma with current and density gradients producing a minimum in the shear Alfvén continuum. GAEs were first discovered numerically by Ross et al.,¹ described analytically by Appert et al.,² and observed in a tokamak by de Chambrier et al.³ Since this initial work, GAEs have been observed in a wide range of plasmas, including in stellerators with small magnetic shear⁴ and Ohmic discharges in spherical tokamaks.⁵ The study of GAEs in conventional tokamaks slowed somewhat after Fu and Van Dam found GAEs were highly continuum damped, stabilising them in a burning plasma.⁶ 

In cylindrical geometry, the structure of the GAE is dominated by a single poloidal harmonic m. This contrasts with gap modes, such as toroidicity-induced Alfvén eigenmodes (TAEs)⁷ and ellipticity-induced Alfvén eigenmodes (EAEs),^8,9 which are dominated by at least two poloidal harmonics coupled by geometrical effects— m and m—1 for TAEs, and m and m—2 for EAEs.

The axisymmetric (toroidal mode number n = 0) GAE was initially highlighted by Appert et al.² as an attractive prospect for plasma heating because it does not couple to the plasma surface. Since then, the axisymmetric mode has been observed in TFTR¹⁰—later identified as a GAE by Villard and Vaclavik¹¹—and in MAST.⁵ McClements et al. concluded that because these modes are easily excited in both conventional and spherical tokamaks, the modes are weakly damped and Alfvén heating by the n = 0 GAE may not be viable.⁵ 

We present recent experimental observations of the n = 0 mode in JET in Sec. II. The frequency of this mode coincides with the n = 0, $| m | = 1$ GAE, which falls within the ellipticity-induced gap. These observations are investigated numerically using MHD codes in Sec. III. In Sec. IV, elliptical effects are incorporated into the Hain-Lüst equation to describe the results analytically. Finally, the mode drive and context of the observations are discussed in Sec. V, and the conclusions are presented in Sec. VI.

We present two typical discharges from an experiment designed to study fusion products in JET.¹² Deuterium-helium-3 (D-³He) plasmas were heated with deuterium neutral beam injection (NBI). Ion cyclotron resonance heating (ICRH) tuned to the third harmonic of the beam cyclotron frequency on the magnetic axis, R ≈ R₀, was used to further accelerate the beam ions. Heating the beam ions from their injection energy ∼100 keV to the MeV range increased the neutron rate by a factor of 5–10 because of the higher D-D fusion cross section.¹³ 

Typical plasma parameters for this experiment were as follows: major axis R₀ ≈ 3 m, equilibrium on-axis magnetic field B(R₀) ≡ B₀ ≈ 2.24 T, minor radius a ≈ 0.9 m, and ³He concentration of puffed gas 7%–12%. Typical electron density (n_e), safety factor (q), and elongation (κ) profiles are presented in Fig. 1. Plasma elongation gives the ratio of the plasma vertical and horizontal radii, κ = b/a. q and elongation profiles are obtained from the equilibrium code EFIT^14,15 constrained by motional Stark effect (MSE) observations.¹⁶ The electron density profile was measured using LIDAR.¹⁷ Typical NBI and ICRH waveforms are shown in Fig. 2.

MHD waves excited by fast ions resulting from NBI and ICRH were detected using toroidally separated Mirnov coils.¹⁸ Toroidal mode numbers were calculated using the relative phase shifts of the magnetic perturbations at the coil toroidal positions. Positive and negative toroidal mode numbers correspond to waves propagating in the co-current and counter-current directions, respectively.

In plasma discharge #86775, a monster sawtooth with a period 2.53 s crashed at t = 11.40 s, as shown in Fig. 2. During the saturated phase of the monster sawtooth, TAEs were observed to pitchfork split and burst. This non-linear behaviour^19,20 can cause significant energetic particle transport. The crash triggered a high amplitude n = 1 neoclassical tearing mode (NTM) that persisted until the end of the pulse. After a short delay, axisymmetric and bidirectional modes with toroidal mode numbers –5 ≤ n ≤ 8 were observed in the EAE frequency range, shown in Fig. 3.

NTMs can degrade plasma confinement and cause a net radial drift of fast trapped particles.^21–23 Losses of trapped D ions with a broad energy distribution up to the MeV energy range were observed using a scintillator probe²⁴ and Faraday cups.²⁵ 

Similarly, in shot #86781, a monster sawtooth with a period 1.15 s crashed at t = 10.16 s. This crash was preceded by tornado modes, TAEs localised inside the q = 1 flux surface. Tornado modes can transport energetic particles out of the q = 1 surface,²⁶ removing the stabilising influence of the ICRH-accelerated fast ions on the monster sawtooth. The crash excited a n = 2 NTM that persisted until t = 11.68 s, as shown in Fig. 4. Figure 5 shows the excitation of modes in the EAE frequency range began at t = 10.26 s and lasted until t = 11.65 s, by which time the amplitude of the NTM had decreased significantly.

The axisymmetric modes were excited with the highest amplitude of the observed modes. Additionally, pitchfork splitting (sideband formation) of the modes was observed, for example, at times t = 10.55, 10.75, 10.94, and 11.34 s in Fig. 5. This suggests that the drive/damping ratio is sufficient to drive a soft non-linear regime,¹⁹ in which the mode linear growth rate is small compared with both the drive due to energetic particles and the damping rate.

The frequency separation between adjacent toroidal mode numbers Δf ∼ 5 kHz. This frequency separation is equal to the toroidal rotation frequency at the mode location.²⁷ This suggests that these modes are excited near the plasma edge, where Landau damping $∝ exp ( − 1 / 9 β i )$ is weak, with β_i(r) representing the ion beta.

The n = 0 modes are observed in the EAE frequency range. The ellipticity-induced gap in the Alfvén continuum is caused by destructive interference between two counter-propagating modes of poloidal mode numbers m and m – 2. Hence, for the ellipticity-induced gap, the two parallel wavenumbers are equal and opposite: $k ∥ m = − k ∥ m − 2$⁠, where $k ∥ m = ( n − m / q ) / R 0$⁠. Setting n = 0 gives m = +1, so that the other dominant harmonic of the n = 0 mode must be m = –1. Furthermore, the axisymmetric mode is not associated with a particular q rational surface.

Henceforth, we use shot #86775 at time t = 11.62 s as a reference time-slice because of the availability of motional Stark effect (MSE) data for a more accurate reconstruction of the q profile. The continuous Alfvén spectrum was computed using the ideal MHD code CSCAS²⁸ based on an equilibrium constructed with the code HELENA.²⁹ The n = 0 continuum is shown in Fig. 6(a), which can be compared to the continuum for $| n | = 5$ in Fig. 6(b).

The observed frequency of the axisymmetric mode f_obs = 325 kHz corresponds to a normalised frequency $ω ̃ o b s = ω o b s R 0 / v A 0 ≈ 1.11$⁠, where v_A0 is the on-axis Alfvén velocity and R₀ is the major radius. $ω ̃ o b s$ lies within the continuum gap $1.06 ≤ ω ̃ g a p ≤ 1.12$⁠, which extends from the plasma core up to the edge, where n_e/n_e(0) → 0 (not displayed in Fig. 6). At the plasma edge, the small tail of the eigenmode crosses the Alfvén continuum and experiences continuum damping. We discuss the calculation of this damping rate in Sec. III D.

The radial structure was then computed using the ideal incompressible MHD code MISHKA-1 (Ref. 30) for the same HELENA equilibrium. A mode with a normalised frequency $ω ̃ MISH = 1.08$ was found. This corresponds to a frequency f_MISH ≈ 317 kHz, which is close to the observed value f_obs = 325 kHz. The n = 0 radial mode structure is shown in Fig. 7. The eigenvector components are related to the poloidal harmonics of the radial and poloidal components of the cross-field E × B drift velocity, V_r and V_θ.

The m = ±1 poloidal harmonics dominate the radial mode structure with comparable amplitudes as anticipated. The eigenmode is referred to as an even mode because the m = ±1 harmonics share the same sign in sV_r. This form is also known as a ballooning mode structure.

By relating the components of the cross-field drift velocity to the wave potential ϕ, we obtain sV_r ∝ mϕ and $V θ ∝ ( 1 / m ) ∂ ( s V r ) / ∂ r$⁠. This relation shows that a factor of m = ±1 accounts for the opposite signs in V_θ. The harmonics of V_θ peak in amplitude at s ≈ 0.80. Therefore, the wave energy $∝ ( ∂ ϕ / ∂ r ) 2$ peaks at s ≈ 0.80, which coincides with the position of the minimum of the upper continuum branch.

An odd (i.e., m = ±1 harmonics have opposite signs in sV_r) mode was found below the lower branch of the continuum, shown in Fig. 8. This mode has an antiballooning structure. The frequency of this mode is $ω ̃ = 0.73$⁠, or f_MISH ≈ 216 kHz. However, a n = 0 mode at this frequency was not excited experimentally.

To delineate the effects of the ellipticity, density, and current gradients—all essential for the frequency range in which the modes were observed—a numerical investigation of the Alfvén spectrum and mode structure at varied ellipticity was performed. The plasma surfaces were constructed with the parametric form of Cartesian coordinates

where a and b are still the plasma vertical and horizontal radii, τ is the triangularity, and the parameter t is related to poloidal angle θ. With τ = 0, κ = b/a was scanned over, and HELENA, CSCAS, and MISHKA were run for each step. The size of the ellipticity-induced gap is described here by $( ω + 2 − ω − 2 ) / ( ω + 2 + ω − 2 )$⁠, where ω₊ and ω₋ are the frequencies at the top and bottom of the gap, and the ellipticity parameter for a straight field line coordinate system is defined as³¹ 

The size of the n = 0 gap increases approximately linearly with ellipticity, as shown by Fig. 9. The gap size drops to zero at e ≈ 3 × 10⁻³, this offset from e = 0 is a higher order effect of toroidal coupling, described by ϵ. For e = 0, ϵ = 0, the m = –1 and m = +1 continua are degenerate.

Additionally, for zero ellipticity (ignoring the toroidicity-induced offset), the even and odd n = 0 modes are degenerate. Figure 10 shows that the n = 0 eigenmode splits into two distinct modes for finite ellipticity. The two eigenmodes are localised at the position of the minimum of each branch. For positive ellipticity, b/a > 1, the frequency of the even mode increases with ellipticity, and the frequency of the odd mode decreases. The frequency evolution is reversed for negative ellipticity, b/a < 1, although this configuration is rarely used experimentally. For $0.05 ≲ | e | < 0.36$⁠, the mode with the higher frequency hits the continuum in the core and is heavily damped. For $0.36 ≤ | e | ≤ 2.00$⁠, the gap is wide enough for the mode to exist without hitting the top of the lower branch or bottom of the higher branch in the core. Therefore, for a sufficiently elliptical plasma cross-section, the n = 0 mode can avoid heavy continuum damping.

With increasing ellipticity, the gradient of the mode at the continuum minimum becomes more steep and the radial mode number increases from 1 to 2. Consequently, the mode energy becomes more localised.

The MHD code AEGIS was used to study the mode damping numerically. AEGIS³² is a free boundary shooting code that solves the ideal MHD equations. The code uses the same EFIT equilibrium input used for Figs. 6–8.

The eigenfrequency computed by AEGIS $ω ̃ A E G = 1.10$ is close to the frequency found by MISHKA (⁠ $ω ̃ MISH = 1.08$⁠) and the experimental frequency (⁠ $ω ̃ o b s = 1.11$⁠). The small difference between $ω ̃ MISH$ and $ω ̃ A E G$ is most likely to be due to small variations in equilibria from processing in HELENA and AEGIS. The eigenmode structure calculated by AEGIS agrees closely with that of MISHKA. A singularity is present in the radial structure of the mode at s = 0.9997, where the eigenmode hits the continuum. At s = 0.9997, the mode amplitude is ∼1% of the maximum amplitude.

We measure the continuum damping of the eigenmode γ_d by introducing a small, localised current δj with the frequency of the wave.³³ Adding a small imaginary component γ to the real part of the frequency ω captures the effect of friction. This allows us to measure the total dissipative power

Scanning over the frequency gap reveals a resonance at $ω ̃ = 1.10$⁠, identifying the eigenmode, given by the plasma displacement ξ. The total energy of the mode E_t is contributed to equally by the kinetic energy and potential energy of the bulk plasma (with density ρ), so we can write the damping rate as

where δ is chosen such that δ ≫ γ* to give the integration limits that capture the total dissipative power and remove the singular contribution from the mode energy. The continuum damping rate converges as γ decreases. Thus, the continuum damping of the mode with frequency $ω ̃ AEGIS = 1.10$ is γ_d/ω ≈ 8.7 × 10⁻³.

To obtain the eigenmode equations for the plasma displacement $ξ$⁠, we begin from the ideal MHD equation of motion

where

where π and Q are the perturbed plasma pressure and magnetic field, respectively, p and B are the equilibrium pressure and magnetic field, respectively, ρ is the equilibrium plasma mass density, and γ is the adiabatic index. In these equations, the constant μ₀ has been dropped for convenience in the style of Goedbloed.³⁴ It can be restored using the transformation $B → B / μ 0$⁠.

We use a Fourier representation for the plasma displacement

where c.c. denotes the complex conjugate, and k = n/R₀.

A straight magnetic field line coordinate system (r, θ, z) for an elliptical plasma,³⁵ correct to first order in ellipticity e, is employed. In the infinite aspect ratio limit, the covariant metric tensor is given by

Following the method of Goedbloed,³⁴ we obtain the n = 0 eigenmode equations for coupled m = ±1 harmonics (denoted by subscript ±)

where

where χ ≡ rξ_r, F = fB, and G = gB, where f and g are the parallel and perpendicular wavenumbers.

In the incompressible limit (γ → ∞), we obtain

This set of two coupled equations can be decoupled by introducing the even and odd combinations of χ₊ + χ₋ and χ₊ – χ₋. In the limit ellipticity e → 0, these two eigenmode equations reduce to the Hain-Lüst equation.^2,36 The term containing the ellipticity parameter e causes splitting of the eigenmode and Alfvén continuum.

Equations (11) are sufficient in conventional tokamaks, but compressibility effects are important in spherical tokamaks because of the higher plasma β. The pressure gradient can be rewritten using the ideal MHD equilibrium condition

The fourth-order differential equation resulting from the two second-order coupled eigenmode equations is singular when the determinant of the coefficients of the second derivative terms in Eq. (10) equals zero, yielding the Alfvén continuum.³⁷ The ellipticity induced continuum is plotted for n = 0 and $| n | = 5$ in Fig. 11.

Without ellipticity, the n = 0 continuum appears as one branch, as the continua branches produced by the m = ±1 are degenerate at the EAE gap frequency ω_EAE(r). For finite ellipticity, the continuum splits into two branches, producing a gap. The gap positions and the qualitative trend with ellipticity are reproduced accurately, but the gap size is slightly under-predicted compared with CSCAS. Coupling with additional poloidal harmonics and second order e effects will account for some of this difference.

We conclude that the n = 0 mode observed in our JET experiment is a GAE. The mode frequency lies below the continuum minimum and the mode energy is localised at the position of the continuum minimum. Although the mode is heavily modified by ellipticity, it exists for a circular cross-section and so is not an EAE.

The other modes with n ≠ 0 observed at the same time are not GAEs because there are no minima in the n ≠ 0 continuum for a given poloidal harmonic. Ellipticity couples two poloidal harmonics to produce continuum gaps, within which these modes exist. Therefore, these modes are EAEs.

A normal mode requires a free energy source to drive the mode unstable and a resonance to transfer the energy from the particles to the wave. The unperturbed motion of a particle can be characterised by the toroidal precessional frequency ω_ϕ and the poloidal bounce frequency ω_θ. Wave-particle resonance requires³⁸ 

where l = 0, ±1, ±2,…. Hence, trapped particles drive axisymmetric (n = 0) waves.

ICRH couples most efficiently to trapped particles with mirror points lying on the ICRH resonant surface, so most RF heated ions in this experiment are deeply trapped particles. This is consistent with the observation that the axisymmetric modes are excited at the same time as NTMs, which are known to transport fast trapped ions radially.^21–23 This suggests that the mode is driven by ICRH-accelerated trapped beam ions ejected from the core by the NTMs.

The energy of the even mode peaks on the low field side of the tokamak, where most of the trapped particles are located. The energy of the odd mode, which has an anti-ballooning radial structure, peaks on the high field side. Therefore, the odd mode has lower drive than the even mode. This may explain why the odd mode was not observed in our experiment.

The drive due to fast particles on a shear Alfvén wave of frequency ω is proportional to^39,40

where f₀(E, P_ϕ, μ) is the distribution function of particles with energy E, magnetic moment μ, and canonical toroidal angular momentum P_ϕ = mRv_ϕ – Zeψ_P, with m, Ze, and v_ϕ giving the particle mass, charge, and toroidal velocity, respectively. E, μ, and P_ϕ are the adiabatic invariants for a particle, so are approximately constant for an unperturbed motion. μ is conserved for particles in a shear Alfvén wave field because the wave frequency is much lower than the ion cyclotron frequency. For an axisymmetric mode, P_ϕ is also conserved because ΔP_ϕ = (n/ω)ΔE.⁴¹ Therefore, the free energy source of the axisymmetric wave is provided by positive gradients in the particle energy distribution function, referred to as a bump-on-tail distribution. Consequently, n = 0 mode excitation can be used to assess the distribution function of energetic particles, as for ion cyclotron emission (ICE).⁴² Additionally, axisymmetric modes present a good opportunity to test models of non-linear drive using experimental data. We will publish a study on this in the near future.

The modes with toroidal mode number n = 0 were observed in JET D-³He plasmas heated with ICRH accelerated deuterium beam ions. The axisymmetric modes were observed at the same time as neoclassical tearing modes (NTMs) triggered by monster sawtooth crashes.

There is a gap in the n = 0 Alfvén continuum, computed with the MHD code CSCAS, which widens linearly with the ellipticity of the plasma surfaces. Within this gap, an even GAE is identified at a frequency close to the experimentally observed frequency using MHD codes MISHKA and AEGIS. The mode energy is localised at the position of the continuum minimum, and the mode structure spans the entire plasma radius.

Below the gap, an odd GAE that was not excited experimentally is found numerically. For zero ellipticity, the two branches of the continuum gap combine, and the even and odd modes are degenerate, existing as an even mode below the single continuum branch. Hence, the mode is not induced by ellipticity, but is significantly modified by ellipticity. For small ellipticity, the gap GAE is heavily continuum damped. For moderate values of ellipticity, the ellipticity-induced gap is sufficiently wide to prevent a strong continuum damping of the mode. The MHD code AEGIS also finds the even GAE with a similar frequency. By introducing a small localised current with the frequency of the wave, the continuum damping rate of the eigenmode is estimated to be γ_d/ω ≈ 0.87%.

The coupled eigenmode equations for shear Alfvén waves in a cylindrical plasma with an elliptical cross-section are derived from the ideal MHD equation of motion. From these equations, the analytical n = 0 Alfvén continuum is found.

We speculate that the mode is driven by the ICRH-accelerated trapped beam ions ejected from the core by NTMs. The axisymmetric mode is of interest because it can be used to diagnose fast particle energy distributions at the mode location. Furthermore, pitchfork splitting was observed in the excitation of the n = 0 mode. This will be the subject of a future study of non-linear drive.

This paper benefited greatly from comments by K. G. McClements and M. Baruzzo. We also thank M. Li for advice on modifications to AEGIS. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 under Grant Agreement No. 633053 and from the RCUK Energy Programme [Grant No. EP/P012450/1]. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

This research was supported by the Office of Fusion Energy Science of the U.S. Department of Energy under Grant No. DE-FG02-04ER54742. AEGIS calculations used resources of the National Energy Research Scientific Computing Center, a Department of Energy Office of Science User Facility supported under Contract No. DE-AC02-05CH11231.