Fast-ion-driven vertical modes in magnetically confined toroidal plasmas

It is well known that tokamak plasmas with non-circular cross sections are prone to an instability, involving axisymmetric modes with toroidal mode number n = 0, which, if left uncontrolled, would lead to vertical displacement events (VDE) and plasma current disruptions. As shown in the pioneering work by Laval et al [1], a nearby perfectly conducting wall can provide passive feedback stabilization. For the realistic case of a resistive wall, the relevant dispersion relation is cubic [2]. Two roots correspond to weakly damped oscillations, with a frequency close to the poloidal Alfvén frequency and damping rate of the order of the inverse of the resistive wall time. The third one is unstable, growing on the resistive wall time, and its complete suppression requires an active feedback stabilization system, whose principle and practical implementation have been the subject of several contributions (see, e.g., references [3–9]). Here, we focus our attention on the two oscillatory roots, and show that they can be driven unstable by the resonant interaction with energetic ions. This fast-ion-driven vertical mode (in brief, FIDVM; see reference [10] for a preliminary account) adds to the catalog of potentially dangerous fast ion instabilities in tokamak fusion plasmas [11, 12]. We suggest that the FIDVM may provide a possible explanation for the saturated n = 0 Alfvén oscillations observed in recent beam-injected JET experiments with additional ICRF heating, tentatively interpreted in reference [13] in terms of n = 0 global Alfvén eigenmodes (GAEs). We point out that the FIDVM is a mode with a discrete frequency below the minimum value of the poloidal Alfvén frequency, not affected by continuum damping. Thus, its frequency does not depend on transient details of the safety factor profile, and as such it appears to be also a robust candidate, on the same footing as GAE, for the interpretation of the JET n = 0 observations.

The FIDVM dispersion relation is derived within the framework of the hybrid kinetic-MHD model [14]. The bulk thermal plasma is described according to ideal-MHD, while the fast particles are modelled by the collisionless drift-kinetic equation [15]. We study the resonant interaction between the oscillatory branch of the n = 0 dispersion relation and fast ions. We evaluate the contributions of particles with trapped and passing orbits to the mode growth rate and the conditions under which these can overcome the damping associated with wall resistivity. We show that instability requires a fast ion distribution function with positive derivative as function of energy. The experimental conditions under which such distribution can be realized are discussed below.

The magnetic field is represented as B = e_ϕ × ∇ψ + B_ϕ e_ϕ , where e_ϕ is the unit vector along the toroidal direction, and B_ϕ is nearly constant. The plasma flow is v = e_ϕ × ∇φ. The magnetic flux function, ψ, and the stream function, φ, obey the dimensionless equations: ${\partial }_{t}\psi +\left[\varphi ,\psi \right]=0$; ${\partial }_{t}\,\nabla \cdot \left(\varrho \,\nabla \varphi \right)+\left[\varrho ,{(\nabla \varphi )}^{2}\right]/2+U\left[\varphi ,\varrho \right]+\left[\varphi ,U\right]=\left[\psi ,J\right]-{\boldsymbol{e}}_{\phi }\cdot \nabla \times \nabla \cdot {\boldsymbol{\mathcal{P}}}_{\text{h}}$, where $\left[\chi ,\eta \right]={\mathbf{e}}_{\phi }\cdot \nabla \chi \times \nabla \eta$, J = ∇² ψ is the current density, and U = ∇² φ is the flow vorticity. The fast particles pressure tensor is defined as ${\boldsymbol{\mathcal{P}}}_{\text{h}}={p}_{\perp \text{h}}\boldsymbol{\mathcal{I}}+{({p}_{{\Vert}}-{p}_{\perp })}_{\text{h}}{\boldsymbol{e}}_{{\Vert}}{\boldsymbol{e}}_{{\Vert}}$, where e _∥ = B /B, p_⊥h and p_∥h are moments of the fast ions distribution function, which obeys the collisionless drift-kinetic equation. The normalization is standard, see reference [10] for details. Time is normalized to the poloidal Alfvén time, ${\tau }_{\text{A}}={\left(4\pi \,{\varrho }_{\text{m}}\right)}^{1/2}/{B}_{\text{p}}^{\prime }$, with ϱ_m the on-axis mass density and ${B}_{\text{p}}^{\prime }$ the on-axis derivative of the equilibrium poloidal field; space is normalized to a convenient equilibrium scale length, ${r}_{0}=a\,b/{[\left({a}^{\,2}+{b}^{\,2}\right)/2]}^{1/2}$ with a and b the minor and major semi-axes of the elliptical plasma boundary.

At equilibrium, we assume stationary fields, no flows, and isotropic fast ions (p_∥h,eq = p_⊥h,eq). To allow for analytic work, and since toroidal effects and current density gradients are not expected to affect the MHD response of the bulk plasma to axisymmetric perturbations, the particle and current densities are taken uniform up to an elliptical boundary, μ = μ_b, where (μ, ϑ) are elliptical coordinates: x = A sinh(μ)cos(ϑ); y = A cosh(μ)sin(ϑ); $A=\sqrt{{b}^{2}-{a}^{2}}$. The equilibrium current density is J_eq(μ) = 2H(μ_b − μ), where H(x) is the unit step function. The equilibrium magnetic flux is ψ_eq = (x²/a² + y²/b²)/2 for μ ⩽ μ_b, ${\psi }_{\text{eq}}=1/2+{{\alpha }_{0}}^{2}\left\{\mu -{\mu }_{b}-{e}_{0}\,\mathrm{sinh}\left[2(\mu -{\mu }_{\text{b}})\right]\mathrm{cos}(2\vartheta )\right\}/2$ for μ > μ_b, where ${\alpha }_{0}^{2}=ab/{r}_{0}^{2}$ and e₀ = (b² − a²)/(b² + a²) is the ellipticity parameter. The equilibrium flux exhibits X-points outside the plasma, located at μ = 2μ_b, ϑ = (π/2, 3π/2). More details on the equilibrium solution can be found in references [16, 17].

Let $\psi (\mu ,\vartheta ,t)={\psi }_{\text{eq}}(\mu ,\vartheta )+\tilde{\psi }(\mu ,\vartheta )\,{e}^{\,\gamma \,t}$ and $\varphi (\mu ,\vartheta ,t)=\tilde{\varphi }(\mu ,\vartheta )\,{e}^{\,\gamma \,t}$, where the over-tilde denotes small perturbations. To first order in perturbed quantities,

$$\begin{equation}\gamma \,\tilde{\psi }+\left[\tilde{\varphi },{\psi }_{\text{eq}}\right]=0,\end{equation} \tag{ 1 }$$

$$\begin{equation}[0.5ex]\gamma \nabla \!\!\cdot ({\rho }_{\text{eq}}\nabla \tilde{\varphi })=[\tilde{\psi },{J}_{\text{eq}}]\!+[{\psi }_{\text{eq}},\tilde{J}]-{\boldsymbol{e}}_{z}\cdot \nabla \times \nabla \cdot {\tilde{\boldsymbol{\mathcal{P}}}}_{\text{h}}.\end{equation} \tag{ 2 }$$

We multiply equation (2) by ${\tilde{\varphi }}^{\ast }/2{\gamma }^{\ast }$ and integrate over the whole volume. After straightforward algebra, we obtain a standard expression for the dispersion relation in terms of quadratic forms,

$$\begin{equation}-{\gamma }^{2}\delta I=\delta {W}_{\mathrm{M}\mathrm{H}\mathrm{D}}+\delta {W}_{\mathrm{f}\mathrm{a}\mathrm{s}\mathrm{t}},\end{equation} \tag{ 3 }$$

where $\delta I={\int }_{\mathcal{V}}{\mathrm{d}}^{3}x\,\rho \boldsymbol{\xi }\cdot {\boldsymbol{\xi }}^{\ast }/2$, $\delta {W}_{\mathrm{M}\mathrm{H}\mathrm{D}}=-{\int }_{\mathcal{V}}{\mathrm{d}}^{3}x\,{\boldsymbol{\xi }}^{\ast }\cdot \boldsymbol{F}(\boldsymbol{\xi })/2$ and $\delta {W}_{\mathrm{f}\mathrm{a}\mathrm{s}\mathrm{t}}={\int }_{\mathcal{V}}{\mathrm{d}}^{3}x{\boldsymbol{\xi }}^{\ast }\cdot \nabla \cdot {\tilde{\boldsymbol{\mathcal{P}}}}_{\text{h}}/2$. All the integrals extend over the plasma volume $\mathcal{V}$, the displacement vector is $\boldsymbol{\xi }={\boldsymbol{e}}_{z}\times \nabla \tilde{\varphi }/\gamma$ and the reduced MHD force operator is $\boldsymbol{F}(\boldsymbol{\xi })=(\tilde{\boldsymbol{J}}\times {\boldsymbol{B}}_{\text{eq}}+{\boldsymbol{J}}_{\text{eq}}\times \tilde{\boldsymbol{B}})$.

Firstly, we neglect energetic particles. As shown in [10], the rigid-shift vertical displacement, $\tilde{\varphi }=\gamma \xi x$, with ξ = const, is the exact solution of equations (1) and (2) for the considered equilibrium. The perturbed flux in the plasma region is obtained from equation (1). In the vacuum region, the flux is modified by the presence of a nearby resistive wall, assumed to lie on an elliptical surface, μ = μ_w, confocal with the plasma elliptical boundary, as in reference [1]. Considering a thin wall, we solve analytically for $\tilde{\psi }$ and obtain the expected cubic dispersion relation for the n = 0 mode:

$$\begin{equation}{\gamma }^{3}+{\gamma }^{2}\frac{1}{{\tau }_{\eta }}-\gamma {\gamma }_{0}^{2}-{\gamma }_{0}^{2}\frac{1}{{\tau }_{\eta }}\frac{1}{1-D}=0,\end{equation} \tag{ 4 }$$

where ${\gamma }_{0}^{2}=-{\omega }_{0}^{2}$ and

$$\begin{equation}{\omega }_{0}\approx {e}_{0}^{1/2}{\tau }_{\text{A}}^{-1}\sqrt{(D-1)},\end{equation} \tag{ 5 }$$

with τ_η the resistive wall time and $D=\left[\mathrm{exp}(4{\mu }_{\text{b}})+1\right]/\left[\mathrm{exp}(2{\mu }_{\text{w}})+1\right]$ a geometrical factor that depends on the distance between the plasma boundary and the wall [1]. If μ_w = 2μ_b, the wall intercepts the X-points and D = 1. Values of D > 1 are obtained if the X-points lie outside the wall, 2μ_b > μ_w > μ_b. The criterion D > 1 corresponds to passive feedback stabilization of the n = 0 mode in the ideal wall limit (τ_η → ∞) and it is the scenario assumed in the remainder of this article. In an actual tokamak with a divertor separatrix, the plasma extends to the X-points, which lie inside the wall, and additional plasma facing components are introduced inside the vacuum chamber to ensure passive feedback stabilization. Furthermore, as shown in [18], even in the no-wall limit (μ_w → ∞), the ideal-MHD constraint imposed on the X-points leads to current sheets, localized on the magnetic separatrix, adding an important ingredient to the n = 0 stability properties. The latter scenario requires careful treatment of the ideal-MHD X-point resonance and is beyond the scope of the present article.

Solving for γ = −iω in the limit D > 1 and small resistivity γ₀ τ_η ≫ 1, we obtain the three relevant roots.

$$\begin{equation}\omega \approx \pm {\omega }_{0}-i\frac{1}{2{\tau }_{\eta }}\frac{D}{(D-1)}=\pm {\omega }_{0}-i{\gamma }_{\eta },\end{equation} \tag{ 6 }$$

$$\begin{equation}\gamma =\frac{1}{(D-1){\tau }_{\eta }}.\end{equation} \tag{ 7 }$$

The two stable roots in (6) oscillate at frequency ω = ±ω₀. As anticipated, the damping rate ${\gamma }_{\eta }=D/\left[2{\tau }_{\eta }(D-1)\right]$ is related to the resistive wall time, but the resonant interaction with fast ions can destabilize these modes, as shown in the following. The third root corresponds to the unstable n = 0 resistive wall mode, which is normally suppressed by active feedback stabilization.

Let us now consider the effects of energetic particles, represented by the last term in equation (3). We adopt a perturbative approach. We introduce the small expansion parameter, _h = β_h/β_c ≪ 1, i.e., the ratio between energetic particle and core plasma pressures. In this limit, |δW_MHD| ≫ |δW_fast|. The cubic dispersion relation, equation (4), is recovered to zeroth order in _h. Following a standard procedure [15, 19], the quadratic term related to energetic particles can be written as δW_fast = δW_h,ad + δW_h,nad. Here, the first term (the 'adiabatic' part) is fluid-like and purely real; the second term ('non-adiabatic') has both real and imaginary parts, the latter contributed by the resonant fast ions. The real part of δW_fast gives rise to $\mathcal{O}({{\epsilon}}_{\text{h}})$ corrections to the n = 0 oscillation frequency and thus can be neglected. On the other hand, the imaginary part of δW_fast contributes to the n = 0 growth rate, competing with the damping term γ_η . Accordingly, only the imaginary part of δW_fast is retained in the following.

Thus, we are led to the following dispersion relation for the two oscillatory n = 0 modes modified by fast ions:

$$\begin{equation}{\omega }^{2}={\omega }_{0}^{2}-2i{\omega }_{0}{\gamma }_{\eta }+i{\omega }_{0}^{2}{\lambda }_{\text{h}}+\mathcal{O}({\gamma }^{2}/{\omega }_{0}^{2}),\end{equation} \tag{ 8 }$$

where ${\lambda }_{\text{h}}=\mathrm{I}\mathrm{m}(\delta {\hat{W}}_{\text{h},\text{nad}})$ is a dimensionless parameter and $\delta {\hat{W}}_{\text{h},\text{nad}}$ the properly normalized quadratic form. Leading order solutions for the real and imaginary parts of the mode frequency are ω = ±ω₀ + iγ_tot, where γ_tot = ω₀ λ_h/2 − γ_η . In order to determine the stability threshold, detailed analysis of δW_h,nad is necessary. For n = 0 modes:

$$\begin{equation}\delta {W}_{\text{h},\text{nad}}={\zeta }_{1}\enspace \sum\limits _{\sigma }\int \mathrm{d}{P}_{\phi }\,\mathrm{d}\mathcal{E}\,\mathrm{d}{\mu }_{\perp }\,{\tau }_{{\Omega}}\omega \frac{\partial F}{\partial \mathcal{E}}\sum\limits _{p\,=-\infty }^{+\infty }\frac{{\left\vert {{\Upsilon}}_{\mathit{\text{p}}}\right\vert }^{2}}{\omega +p{\omega }_{{\Omega}}},\end{equation} \tag{ 9 }$$

where ζ₁ = −(2π² c)/(Zem²) and F is the unperturbed fast ion distribution function. The σ-summation is over co- and counter-circulating particle orbits, and integration is performed over the three invariants of particle motion: toroidal canonical momentum P_ϕ , energy $\mathcal{E}$, and magnetic moment μ_⊥; ω_Ω is the bounce frequency of trapped/passing particle orbits, and τ_Ω = 2π/ω_Ω. The p-summation is the Fourier expansion of the perturbed Lagrangian over the harmonics of the unperturbed particle orbits, labelled by integer p, with Fourier coefficients

$$\begin{equation}{{\Upsilon}}_{p}(\mathcal{E},{\mu }_{\perp },{P}_{\varphi })=\oint \frac{\mathrm{d}\tau }{{\tau }_{{\Omega}}}\tilde{\mathcal{L}}\,\mathrm{exp}(ip{\omega }_{{\Omega}}\tau ),\end{equation} \tag{ 10 }$$

where the loop integral is over a closed banana/transit orbit, and the perturbed Lagrangian is $\tilde{\mathcal{L}}\approx \mathcal{E}(2-{\Lambda})\boldsymbol{\xi }\cdot \boldsymbol{\kappa }$, with r the radial variable and κ the curvature vector [15]. We consider the 'thin-orbit' approximation, where the radial excursion of both passing and trapped orbits from a reference magnetic surface is negligibly small.

To leading order in the parameter _h, the displacement vector ξ corresponds to the rigid-shift vertical displacement, which is orthogonal to the main direction of toroidal curvature. Since the perturbed Lagrangian is proportional to the scalar product ξ ⋅ κ , the interesting consequence is that toroidal curvature contributions to $\delta {\hat{W}}_{\text{h},\text{nad}}$ are negligibly small. The perturbed Lagrangian can be rewritten as $\tilde{\mathcal{L}}\approx {{\epsilon}}^{2}\mathcal{E}(2-{\Lambda})\xi \,\mathrm{sin}(\theta )/r$ with = r/R₀, R₀ the tokamak major radius and θ the standard poloidal angle. The pitch angle variable ${\Lambda}={B}_{0}{\mu }_{\perp }/\mathcal{E}$ ranges from 0 ⩽ Λ ⩽ 1 − for passing orbits, to 1 − ⩽ Λ ⩽ 1 + for trapped orbits.

The mode-particle resonant condition is ω + pω_Ω = 0. Note that the p = 0 harmonic is non resonant. The resonant frequency is proportional to the particle velocity, v, as we can write ω_Ω = vh_Ω(r, Λ)/(R₀ q), where the safety factor q = rB_ϕ /(R₀ B_p) is approximately constant and the dimensionless function h_Ω involves well-known combinations of elliptic integrals. Changing integration variables, the pole contribution yields

$$\begin{equation}{\lambda }_{\text{h}}={\zeta }_{2}\enspace \sum\limits _{p=1}^{+\infty }\int r\,\mathrm{d}r\,\mathrm{d}{\Lambda}\,\frac{{({v}_{p}^{\ast })}^{3}}{{h}_{{\Omega}}}{\left.\frac{\partial F}{\partial \mathcal{E}}\right\vert }_{{v=v}_{p}^{\ast }}{\left.{\left\vert {{\Upsilon}}_{p}\right\vert }^{2}\right\vert }_{{v=v}_{p}^{\ast }},\end{equation} \tag{ 11 }$$

where ${\zeta }_{2}=(8{R}_{0}{\pi }^{4})/(\mathcal{V}{\rho }_{\text{c}}{\xi }^{2}{\omega }_{0}^{2})$, ρ_c is the core mass density and ${v}_{p}^{\ast }={\omega }_{0}{R}_{0}q/\left[p{h}_{{\Omega}}(r,{\Lambda})\right]$ is the resonant velocity.

Clearly, the sign of λ_h depends on the derivative of the fast particles equilibrium distribution function with respect to energy. Instability requires $\partial F/\partial \mathcal{E} > 0$, giving rise to a positive λ_h. Such a distribution can be obtained, for instance, from Fokker–Plank calculations including a fast particles loss term, as shown in reference [20]. We consider a simplified model with a monochromatic source term at velocity v_fast, and a velocity-independent loss frequency, ν_loss. Assuming a uniform radial distribution of fast ions up to r = r_h, and solving the relevant Fokker–Plank equation, we obtain

$$\begin{equation}F(r,\mathcal{E})=CH({r}_{\text{h}}-r)\frac{H({v}_{\mathrm{f}\mathrm{a}\mathrm{s}\mathrm{t}}-v)}{{({v}^{3}+{v}_{\text{c}}^{3})}^{1-\alpha }},\end{equation} \tag{ 12 }$$

where C is a normalization constant, α = ν_loss τ_s/3 and v_c ≪ v_fast is the critical velocity. With this model, when the fast particle slowing down time, τ_s, is larger than the loss time, τ_loss = 1/ν_loss, the usual, monotonically decreasing slowing down distribution cannot form, and F develops a positive slope as function of energy. Experimental evidence indicates that $\partial F/\partial \mathcal{E} > 0$ can also be induced by a time-dependent source, see reference [21] for the case of modulated neutral beam injection. Another relevant scenario is discussed in [22], where source modulation is caused by sawtooth oscillations. Finally, as discussed in [23] in connection with the excitation of n = 0 geodesic acoustic modes, a fast ion distribution arising from neutral beam injection, with a positive slope as function of energy, can be observed on time intervals, following neutral beam injection time, that are short compared with the fast ion slowing down time.

Next, we investigate the Fourier coefficients ϒ_p . We are going to show that passing particles contribute mostly to the p = 1 harmonic, trapped particles mostly to the p = 1 and p = 2 harmonics, while all other harmonics can be neglected. This result is obtained in the thin-orbit limit, where r ≈ const along the particle orbit and the perturbed Lagrangian depends on time only through the poloidal angle θ. Therefore, we can set ${{\Upsilon}}_{p}=\mathcal{E}(2-{\Lambda})(\xi /r){X}_{p}$, where X_p = ⟨sin(θ)exp(ipω_Ω τ)⟩ and ⟨⟩ means orbit averaging.

Let us consider the transit frequency of passing particles, ${\omega }_{{\Omega}}\to {\omega }_{t}=\pi v\kappa \sqrt{{\epsilon}/2}/\left[{R}_{0}qK(1/{\kappa }^{2})\right]$, where K(x) is the complete elliptic integral of the first kind and κ² = 1/2 + (1 − Λ)/(2) (for passing orbits, κ ⩾ 1). We find that ω_t τ = πF(θ/2|1/κ²)/K(1/κ²) ≈ θ, where F(ϕ|x²) is the incomplete elliptic integral of the first kind and the last equality is valid for most values of κ, with the exception of κ very close to unity, which corresponds to the barely passing limit. Since for barely passing particles ω_t → 0, these particles are non-resonant, and so their contribution to λ_h can be safely neglected. In this way, X_p reduces to $i(1/4K(1/{\kappa }^{2})){\int }_{0}^{2\pi }\mathrm{d}\theta \left[\mathrm{sin}(\theta )\mathrm{sin}(p\theta )\right]/\sqrt{1-{\mathrm{sin}}^{2}(\theta /2)/{\kappa }^{2}}$, and in the asymptotic limit κ ≫ 1, ${X}_{p}\approx i(1/2\pi ){\int }_{0}^{2\pi }\mathrm{d}\theta \,\mathrm{sin}(\theta )\mathrm{sin}(p\theta )$. It is clear that the only harmonic that contributes in this limit is p = 1, while the average is zero for all other p values. Figure 1 displays |X_p |² for the first two harmonics as functions of κ, showing that even considering all κ values the contribution of passing particles for harmonics with p > 1 is indeed negligible.

Zoom In Zoom Out Reset image size

For trapped particles, 0 ⩽ κ ⩽ 1, the bounce frequency is ${\omega }_{{\Omega}}\to {\omega }_{\text{b}}=\pi v\sqrt{{\epsilon}/2}/\left[2{R}_{0}qK({\kappa }^{2})\right]$. Then, ω_b τ = π/2F(ζ|κ²)/K(κ²), where sin(θ/2) = κ sin(ζ). For κ sufficiently far from the barely trapped limit κ = 1, where the bounce frequency drops to 0, it is possible to approximate ω_b τ ≈ ζ. Therefore, one obtains ${X}_{p}=\kappa /K({\kappa }^{2}){\int }_{-\pi /2}^{\pi /2}\mathrm{d}\zeta \,\mathrm{sin}(\zeta )\mathrm{sin}(p\zeta )$. As shown in figure 1, this is zero for both barely trapped (κ = 1) and deeply trapped (κ = 0) particles. It can be shown that all even-p harmonics yield a non-zero X_p ; however, values of X_p with p ⩾ 4 are negligibly small. For odd-p values, X_p = 0, with the only exception of p = 1. Thus, we conclude that only the p = 1 and p = 2 harmonics are of interest. In figure 1, values of |X_p |² for the p = 1 and p = 2 harmonics are plotted as functions of κ.

The reason why deeply trapped (κ = 0) as well as barely trapped/passing (κ = 1) particles do not contribute to ϒ_p , and thus to λ_h, is very simple. For n = 0 vertical displacements, ξ ⋅ κ ∼ sin θ. Deeply trapped particles are localized near θ = 0; barely trapped particles spend most of their time near θ = π; hence, the orbit-averaged value of the scalar product ξ ⋅ κ vanishes in both these limits.

In order to determine stability thresholds, we consider the model distribution function in equation (12), where the normalization constant is determined by ∫d³ vF = n_h, with n_h the fast particle density. After straightforward algebra, we obtain

$$\begin{align}\hfill & {\lambda }_{\text{h}}=\frac{9\pi \alpha (\alpha -1){q}^{2}}{{a}^{2}{R}_{0}^{2}}\frac{{n}_{\text{h}}{m}_{\text{h}}}{{n}_{\text{c}}{m}_{\text{c}}}\sum\limits _{p=1,2}{\left(\frac{{v}_{p0}^{\ast }}{{v}_{\mathrm{f}\mathrm{a}\mathrm{s}\mathrm{t}}}\right)}^{3\alpha }\frac{1}{{p}^{2}}\hfill \\ \hfill & \qquad \times \int \mathrm{d}r\,\mathrm{d}{\Lambda}\,\frac{{r}^{3}{(2-{\Lambda})}^{2}\vert {X}_{p}{\vert }^{2}}{{h}_{{\Omega}}{(r,{\Lambda})}^{3\alpha +3}}H({r}_{\text{h}}-r)H({v}_{\mathrm{f}\mathrm{a}\mathrm{s}\mathrm{t}}-{v}_{p}^{\ast }),\hfill \end{align} \tag{ 13 }$$

where n_c is the core plasma density, m_c and m_h are the core and fast particle mass, respectively, and ${v}_{p0}^{\ast }={\omega }_{0}{R}_{0}q/p$ is the resonant velocity of passing particles with Λ = 0. We can write ${\lambda }_{\text{h}}=({n}_{\text{h}}/{n}_{\text{c}})({m}_{\text{h}}/{m}_{\text{c}})({q}^{2}{a}^{2}/{R}_{0}^{2})\lambda$, where $\lambda (\alpha ,{r}_{\text{h}}/a,{v}_{p0}^{\ast }/{v}_{\mathrm{f}\mathrm{a}\mathrm{s}\mathrm{t}})$ is a dimensionless factor that depends on three parameters. Inserting this expression for λ_h in equation (8), the instability threshold can be cast in the form

$$\begin{equation}\frac{{n}_{\text{h}}}{{n}_{\text{c}}}\geqslant {\left(\frac{{n}_{\text{h}}}{{n}_{\text{c}}}\right)}_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}=\frac{{m}_{\text{c}}}{{m}_{\text{h}}}\frac{{R}_{0}^{2}}{{q}^{2}{a}^{2}}\frac{2{\gamma }_{\eta }}{\lambda {\omega }_{0}}.\end{equation} \tag{ 14 }$$

To gain further insight on realistic numerical values for the critical density threshold, we consider the parameters of the JET experiments discussed in reference [13], where saturated n = 0 oscillations were observed, and MeV fast ions were produced by a combination of ICRH and NBI. Some parameters are experimentally well known, e.g., R₀ = 3 m, a = 0.9 m, elongation b/a ≈ 1.3, which corresponds to e₀ ≈ 0.3, toroidal magnetic field B_T = 2.2 T, q ≈ 1, and deuterium main ion species, yielding an inverse poloidal Alfvén time ${\tau }_{\text{A}}^{-1}\approx 2\times 1{0}^{-6}$ s⁻¹. The experimentally observed n = 0 oscillation frequency is in the range f₀ = ω₀/2π ≈ 300 kHz, which compares well with the theoretical predicted value of the frequency ω₀, defined below equation (4), if we assume realistic values for the geometrical factor, D ≈ 3. Fast ions in the mentioned JET experiments are also mostly deuterium, thus we take m_h/m_c = 1. We assume r_h/a = 0.5 and α ≈ 2. The remaining parameter to be evaluated, ${v}_{p0}^{\ast }/{v}_{\mathrm{f}\mathrm{a}\mathrm{s}\mathrm{t}}={\omega }_{0}{R}_{0}q/(p{v}_{\mathrm{f}\mathrm{a}\mathrm{s}\mathrm{t}})$, depends on the fast ion cut-off energy, ${E}_{\mathrm{f}\mathrm{a}\mathrm{s}\mathrm{t}}={m}_{\text{h}}{v}_{\mathrm{f}\mathrm{a}\mathrm{s}\mathrm{t}}^{2}/2$. Realistic JET estimates indicate values of E_fast in the MeV range. We can see that the p = 1 harmonic, contributed by the passing particles, is the leading contribution for fast ion energies up to ∼1.5 MeV. At higher energies, the p = 2 harmonic contributed by trapped particles can also give a significant contribution. On the other hand, for the considered parameters, the p = 1 harmonic contributed by trapped particles gives a non-zero contribution only for cut-off energy E_fast ⩾ 5 MeV. Let us write λ = λ_p=1,pass + λ_p=2,trap. In figure 2, we present plots of these quantities as functions of E_fast, having fixed all the other parameters to the JET-relevant values indicated above. The overall value of the parameter λ is estimated of order 3 × 10⁻¹ for the JET experiments of interest.

Zoom In Zoom Out Reset image size

The damping rate due to the resistive wall can be estimated of the order of γ_η ∼ 3 × 10² s⁻¹. As a consequence, the critical instability threshold for the relevant JET experiments can be estimated to be ${({n}_{\text{h}}/{n}_{\text{c}})}_{\mathrm{c}\mathrm{r}\mathrm{i}\mathrm{t}}\sim 1\times 1{0}^{-2}$.

In this analysis, we have restricted ourselves to fast ion distribution functions that are isotropic in velocity space. In principle, pitch-angle anisotropy may also contribute to the instability drive. This should make the instability more ubiquitous in many tokamak experiments, where fast ions are produced by neutral beam injection and/or by ion cyclotron radio frequency heating. The effect of velocity space anisotropy will be investigated in a future article.

As we have already pointed out, the theory presented in this article is motivated in part by the observation of saturated n = 0 fluctuations, with a frequency of the order of the poloidal Alfvén frequency, in recent JET experiments where fast ions are produced by auxiliary heating, see reference [13]. At that time, the FIDVM instability was not a known concept, and the observations were tentatively interpreted in [13] in terms of a saturated n = 0 GAE. It is early for us to conclude whether, in fact, the mode observed at JET is a FIDVM rather than a GAE: more experiments are required, but also, the theory developed here ought to be developed further. Nevertheless, we can indicate what are the main points of distinction between GAE and FIDVM that may facilitate the experimental identification. These are basically three. First, the GAE mode frequency for n = 0 is given by [24] ${\omega }_{\mathrm{G}\mathrm{A}\mathrm{E}}={v}_{\text{A}}/qR={\tau }_{\text{A}}^{-1}$, where τ_A is the poloidal Alfvén as defined in the first page of this article. On the other hand, the FIDVM mode frequency is given by equation (5), and, taking the parameter D of order unity, it falls below the GAE mode frequency, since the ellipticity parameter e₀ is typically below unity. Secondly, the FIDVM scales as the square root of e₀, while the n = 0 GAE mode frequency is independent of elongation. Indeed, the GAE would survive in the circular limit, while the FIDVM would not. Thirdly, and perhaps most importantly, the FIDVM mode structure is different from that of the GAE. The FIDVM is a vertical mode, corresponding to a vertical shift of the plasma cross section. This signature would be easily detected by magnetic perturbation coils placed top and bottom of the tokamak midplane. The GAE mode structure favours instead a ballooning type of parity, with the relevant perturbed flux being an even function of the poloidal angle. It also happened at JET, as reported in reference [13], that the GAE mode frequency was close to the frequency of ellipticity-induced Alfvén eigenmodes (EAE) [25]. In this case, however, the distinction between the FIDVM and the EAE is easy, because EAEs are not axisymmetric (their toroidal mode number is n ≠ 0).

In conclusion, we have presented a new type of fast ion instability, dubbed FIDVM, involving n = 0 axisymmetric modes in a tokamak plasma, with a well defined, discrete frequency, ${\omega }_{0}={e}_{0}^{1/2}{\tau }_{\text{A}}^{-1}\sqrt{D-1}$, determined by passive feedback stabilization of vertical displacements in plasma with elongated cross-section. The mode is driven unstable by a mode-particle resonant interaction involving fast ions on both trapped and passing orbits. Instability requires a fast ion distribution with a positive slope as function of energy. The critical fast ion density threshold for instability is determined by the competition between the mode-particle resonance drive and damping introduced by the resistive wall.

Even though the theory presented in this article is linear, it is appropriate to speculate on what might be the nonlinear consequences of this instability. Typically, n = 0 perturbations do not cause a direct radial transport of resonant particles; at most, we might expect that they induce a faster relaxation of the fast ion distribution function, by eliminating 'bumps' in the tail of the distribution, thus relaxing it faster to a standard slowing down distribution. On the other hand, vertical displacements are potentially dangerous, because of their connection with VDEs. One aspect of axisymmetric MHD modes that deserves further attention is the fact that these modes are resonant at the X-point(s) of a tokamak divertor [18]. The resonance is likely to drive perturbed current sheets in the vicinity of the X-points and extending along the separatrix, as numerical simulations reported, e.g., in reference [5] also seem to indicate. These current sheets may in turn affect the stability of edge-localized modes and, more in general, the dynamics of plasma flows in the divertor region.

The theoretical model presented in this article is a proof of concept, indicating that the FIDVM is a potential instability in fusion burning tokamak plasmas. Preliminary estimates indicate that the theoretically predicted mode frequency and stability threshold compare favourably with JET experimental results [13, 22]. However, a more detailed comparison between theory and experiments requires numerical work using realistic equilibria and a careful assessment of the fast ion distribution function.

Acknowledgments