Resonant theory of kinetic ballooning modes in general toroidal geometry

Abstract

The linear theory of the kinetic-ballooning-mode (KBM) instability is extended to capture a weakly driven regime in general toroidal geometry where the destabilization is caused by the magnetic-drift resonance of the ions. Such resonantly destabilized KBMs are characterized by broad eigenfunctions along the magnetic-field line and near-marginal positive growth rates, even well below the normalized-plasma-pressure ($\beta$) threshold of their non-resonant counterparts. This unconventional (or sub-threshold) KBM (stKBM), when destabilized, has been shown to catalyze an enhancement of turbulent transport in the Wendelstein 7-X (W7-X) stellarator (Mulholland et al. 2023 Phys. Rev. Lett. vol. 131, 185101; 2025 Nucl. Fusion vol. 65, 016022). Simplifying the energy dependence of key resonant quantities allows for an analytical treatment of this KBM using the physics-based ordering from the more general equations of Tang et al. (1980 Nucl. Fusion vol. 20, 1439). Results are then compared with high-fidelity gyrokinetic simulations for the (st)KBM in W7-X and the conventional KBM in a circular tokamak at both high and low magnetic shear, where good agreement is obtained in all cases. This reduced KBM model provides deeper insight into (sub-threshold) KBMs and their relationship with geometry, and shows promise for aiding in transport model development and geometry-based turbulence optimization efforts going forward.

1. Introduction

Achieving high normalized plasma pressure $\beta \equiv \beta _{\textrm{ref}} = 8\pi n_{\textrm{e0}} T_{\textrm{e0}}/B_0^2$ is paramount for the successful operation of a future magnetic-confinement-fusion (MCF) power plant. Here, $n_{\textrm{e0}}$ and $T_{\textrm{e0}}$ respectively are the background electron density and temperature, and $B_0$ is the reference magnetic-field strength on axis. High $\beta$ provides the necessary conditions for efficient energy generation from fusion reactions, but such high-pressure plasmas are prone to hosting electromagnetic instabilities, such as kinetic ballooning modes (KBMs). These instabilities modify turbulent behavior and facilitate unwanted heat and particle transport, lowering the energy-confinement time of the fusion reactor. These KBMs are expected to thrive in high- $\beta$ regimes near the ideal magnetohydrodynamic (iMHD) threshold $\beta _{\textrm{crit}}^{\textrm{MHD}}$ . Both KBMs and iMHD ballooning modes are pressure-gradient-driven instabilities, but KBMs differ from their iMHD counterpart in that they are modified by kinetic effects such as magnetic-drift resonances, finite-Larmor-radius (FLR) stabilization, trapped particles, Landau resonances and collisions. These kinetic effects can have both destabilizing and stabilizing influences on the instability (Hastie & Hesketh 1981; Cheng 1982a , b ; Hirose et al. 1994, 1995; Hirose & Elia 1996).

Many theoretical studies of KBMs in toroidal geometry followed from the development of the ballooning transformation, first detailed in Connor et al. (1978, 1979), which constitutes a procedure for reconciling the long parallel wavelength characteristic of plasma instabilities with periodicity of the geometry in a sheared toroidal magnetic field. Utilizing this procedure, the theory of KBMs in a general toroidal geometry was developed in the seminal works of Antonsen & Lane (1980) and Tang, Connor & Hastie (1980). Following this, the influence of kinetic effects was studied in greater detail, including the impact of wave–particle resonances on KBMs (Hastie & Hesketh 1981; Cheng 1982a,b; Hirose et al. 1994, 1995; Hirose & Elia 1996). Further theoretical studies of KBMs in a tokamak geometry can be found in Tang et al. (1985), Kotschenreuther (1986), Hong, Horton & Choi (1989), Hong & Horton (1989), Biglari & Chen (1991), Kim, Horton & Dong (1993), Tsai & Chen (1993), Zonca et al. (1996, 1999), Zonca & Chen (2006), Aleynikova & Zocco (2017) and Zocco et al. (2018a ). To support and advance this research, extensive numerical studies of electromagnetic instabilities and turbulence have been carried out both in tokamaks (Snyder 1999, Pueschel, Kammerer & Jenko 2008; Pueschel & Jenko 2010; Bierwage et al. 2010a , b ; Belli & Candy 2010; Ishizawa et al. 2013; Maeyama et al. 2014; Aleynikova & Zocco 2017; Aleynikova et al. 2018; Kumar et al. 2021; Chen et al. 2023; Pueschel et al. 2025) and stellarators (Baumgaertel et al. 2012; Ishizawa et al. 2014, 2015, 2017; Aleynikova et al. 2018; Ishizawa et al. 2019; McKinney et al. 2021; Mishchenko et al. 2021, 2022; Aleynikova, Zocco & Geiger 2022; Mishchenko et al. 2023; Mulholland et al. 2023, 2025). Recent studies in the Wendelstein 7-X (W7-X) stellarator have revealed the deleterious manifestation of unconventional KBMs far below $\beta _{\textrm{crit}}^{\textrm{MHD}}$ (Mulholland et al. 2023, 2025). This variant has been classified as the sub-threshold KBM (stKBM), characterized by a destabilization threshold $\beta _{\textrm{crit}}^{\textrm{stKBM}}$ well below $\beta _{\textrm{crit}}^{\textrm{MHD}}$ , and a consistently low but finite growth rate until reaching the conventional KBM threshold $\beta _{\textrm{crit}}^{\textrm{KBM}}$ , at which point the stKBM transitions to a conventional KBM rapidly destabilized by further increases in $\beta$ . The stKBM may be subdominant to coexisting instabilities present at low $\beta$ . Despite its low growth rate, this instability has been shown to facilitate an enhancement of turbulent transport at low $\beta \geqslant \beta _{\textrm{crit}}^{\textrm{stKBM}}$ in ion-temperature-gradient (ITG)-driven turbulence regimes (Mulholland et al. 2023, 2025). Consequently, controlling stKBMs is of urgent interest in order to successfully achieve high performance in future MCF devices. Before stKBMs can be adequately controlled, a deeper understanding of them must first be established. This work aims to expand the knowledge basis for stKBMs by developing KBM theory that accounts for magnetic-drift-resonance effects in three-dimensional flux-tube geometry. A numerical implementation of this theory is then shown to be capable of quickly and reliably predicting the existence and behavior of both stKBMs and conventional KBMs in general toroidal geometry.

In this work, gyrokinetic theory is developed to obtain a simplified KBM equation suitable for general toroidal geometry that retains sufficient physics to capture distinctions between conventional KBMs and iMHD-like ballooning modes, as well as stKBMs. Previous studies have shown that stKBM growth rates are slightly augmented by the absence of trapped particles and parallel magnetic-field fluctuations $\delta B_\parallel$ , but these effects need not be present for stKBMs to be unstable (Mulholland et al. 2025). These findings have informed the simplifications made in the theory developed here: kinetic effects such as FLR stabilization and wave–particle resonances from passing ions subject to the magnetic-curvature drift (i.e., ion magnetic-drift resonances or toroidal resonances) are retained, while effects due to trapped particles, Landau resonances and compressional magnetic fluctuations are neglected. The absence of $\delta B_\parallel$ is partially compensated for by replacing the $\boldsymbol{\nabla }B$ -drift frequency $\omega _{\boldsymbol{\nabla }B}$ with the magnetic-curvature-drift frequency $\omega _\kappa$ (see definitions in table 2) (Zocco, Helander & Connor 2015; Aleynikova et al. 2022). A collisionless plasma is assumed. The resulting KBM equation is reformulated as a nonlinear eigenvalue problem (i.e., having nonlinear dependence on the eigenvalue), amenable to numerical evaluation to obtain self-consistent (st)KBM eigenvalues and eigenvectors in any flux-tube geometry. Applying this approach in W7-X reveals that ion magnetic-drift resonances are necessary to explain the low- $\beta$ destabilization of the stKBM. In geometries where stKBMs remain absent (such as a circular tokamak), introducing resonance effects can lead to a moderate reduction in the destabilization threshold of the conventional KBM (relative to its non-resonant counterpart), a softer onset (i.e., a less rapidly increasing growth rate with $\beta$ near marginality) and a reduced growth rate at higher $\beta$ . In all cases, the resonant-destabilization threshold differs from the non-resonant-destabilization threshold only in the presence of a finite ion-temperature gradient ( $\eta _{\textrm{i}}\neq 0$ ; see definition in table 2), and the separation between these thresholds is proportional to $\eta _{\textrm{i}}$ (compare Cheng 1982a , b ), Hirose et al. (1994, 1995), Zonca et al. (1996, 1999)). The results from this reduced KBM model are in good agreement with numerical studies using a high-fidelity gyrokinetic code, highlighting that the physics included here is sufficient to explain and predict the existence (or absence) of stKBMs. These studies indicate that the emergence of stKBMs is attributable to a combination of geometric characteristics, including local/background magnetic shear and magnetic-field-line curvature. This enhanced understanding of stKBMs marks the first step toward designing an stKBM-resilient stellarator, more suitable for high-performance scenarios.

The remainder of this article is structured as follows. In § 2, an analytical theory of KBMs is outlined. In § 2.2, the KBM equation is reformulated as a nonlinear eigenvalue problem (having nonlinear dependence on the eigenvalue), providing a means to obtain self-consistent eigenvalues and eigenvectors for (st)KBMs in general toroidal geometry. In § 3, results of the KBM model are presented and compared with high-fidelity simulations in the W7-X stellarator and a high-/low-magnetic-shear tokamak, followed by a physical interpretation and discussion. In § 4, the article is closed with a summary of key results and an outlook.

2. Kinetic-ballooning-mode theory

In this section, an eigenmode equation for the KBM is developed in general toroidal geometry; this equation is a radially local (i.e., at fixed radial position $r/a$ , with minor-radial coordinate $r$ and minor radius $a$ ) second-order ordinary differential equation for the eigenmode along the magnetic-field line. This derivation starts from the KBM theory presented in Tang et al. (1980), and leads to the general toroidal geometry equivalent of the KBM equation derived in Cheng (1982b ), where the latter was formulated for circular tokamak geometry. The resulting equation includes ion-magnetic-drift resonances and FLR effects, but neglects trapped particles, Landau resonances, collisions and compressional magnetic fluctuations. Despite these reductions, the resulting equation is found to contain sufficient physics to capture important (st)KBM behavior in both tokamaks and stellarators. Following this derivation, the equation is presented in a dimensionless form amenable to numerical implementation.

2.1. Simplified KBM equation in general toroidal geometry

Given the highly anisotropic character of plasma instabilities and turbulent fluctuations in fusion devices – where correlation lengths parallel to the magnetic field are typically 2–3 orders of magnitude greater than their perpendicular counterparts (Goerler et al. 2011) – it is convenient to employ a field-aligned coordinate system (Helander 2014). A general coordinate system is first introduced here followed by a related system utilized in flux-tube simulations (Xanthopoulos et al. 2009). To start, one can introduce $(\psi ,\theta ,\zeta )$ , where $\psi$ is the toroidal magnetic flux within a magnetic surface and acts as a radial variable, $\theta$ is the straight field-line poloidal angle and $\zeta$ is the straight field-line toroidal angle (Boozer 1981; Beer 1995; Faber et al. 2018). A field-line-labeling coordinate (i.e., fixed for a given field line) can then be defined as $\alpha = q(\psi )\theta - \zeta$ , where $q(\psi ) = {\rm d}\zeta /{\rm d}\theta = \iota (\psi )^{-1}$ is the safety factor and $\iota (\psi )$ is the rotational transform. This allows one to write the magnetic field in Clebsch form

(2.1) \begin{equation} \boldsymbol{B} = \boldsymbol{\nabla }\psi \boldsymbol{\cdot} \boldsymbol{\nabla }\alpha , \end{equation}

or alternatively

(2.2) $$\begin{align}\boldsymbol{B} &= \boldsymbol{\nabla }\psi \boldsymbol{\cdot} \boldsymbol{\nabla }\theta - \iota \boldsymbol{\nabla }\psi \boldsymbol{\cdot} \boldsymbol{\nabla }\zeta,\end{align}$$

(2.3) $$\begin{align}&= \boldsymbol{\nabla }\psi \boldsymbol{\cdot} \boldsymbol{\nabla }\theta + \boldsymbol{\nabla }\zeta \boldsymbol{\cdot} \boldsymbol{\nabla }\chi , \end{align}$$

where $\chi$ is the poloidal flux such that $\iota (\psi ) = \text{d}\chi /\text{d}\psi$ . In this coordinate system, the Jacobian is given by

(2.4) \begin{equation} J = (\boldsymbol{\nabla }\psi \boldsymbol{\cdot} \boldsymbol{\nabla }\theta \boldsymbol{\cdot }\boldsymbol{\nabla }\zeta )^{-1} = (\boldsymbol{B} \boldsymbol{\cdot }\boldsymbol{\nabla }\zeta )^{-1}. \end{equation}

In order for the following discussion to be consistent with the numerical flux-tube approach used by Gene (Jenko et al. 2000) (see http://www.genecode.org for code details and access), the following field-aligned coordinate system $(x,y,z)$ is used (Xanthopoulos et al. 2009):

(2.5) $$\begin{align}x &= L_z\sqrt {s},\end{align}$$

(2.6) $$\begin{align}y &= L_z \frac {\sqrt {s_0}}{q_0} \alpha,\end{align}$$

(2.7) $$\begin{align}z &= \theta,\end{align}$$

where $s=\psi /\psi _{\textrm{edge}}$ is the normalized toroidal flux indicating radial position, $s_0$ is the radial position of a chosen flux surface (centre of the flux-tube domain), $L_z$ is the macroscopic reference length (typically chosen to be the minor radius $a$ , the major radius $R$ or the connection length $L_{\textrm{c}} = q_0 R$ , where $q_0 = q(s_0)$ is the safety factor on a given flux surface) and $\theta$ is a field-line-following poloidal angle coordinate (analogous to the ballooning angle (Beer 1995; Faber et al. 2018)). In this coordinate system, the normalized magnetic field and Jacobian are

(2.8) $$\begin{align}\hat {\boldsymbol{B}} &\equiv \frac {\boldsymbol{B}}{B_0}= \boldsymbol{\nabla }x \boldsymbol{\cdot} \boldsymbol{\nabla }y,\end{align}$$

(2.9) $$\begin{align}J &= (\boldsymbol{\nabla }x \boldsymbol{\cdot} \boldsymbol{\nabla }y \boldsymbol{\cdot }\boldsymbol{\nabla }z)^{-1} = (\hat {\boldsymbol{B}} \boldsymbol{\cdot }\boldsymbol{\nabla }z)^{-1}.\end{align}$$

With the appropriate coordinate system established, the simplified KBM equation can now be derived. For the theory presented here, focus is given to the intermediate-frequency regime outlined in Tang et al. (1980)

(2.10) \begin{equation} \omega _{\textrm{b,i}}, \omega _{\textrm{t,i}} \ll \omega \ll \omega _{\textrm{b,e}}, \omega _{\textrm{t,e}} \ , \end{equation}

where $\omega = \omega _{\textrm{r}} + i\gamma$ is the complex eigenfrequency of the KBM with real frequency $\omega _{\textrm{r}}$ and growth rate $\gamma$ , and $\omega _{\textrm{b},j}$ and $\omega _{\textrm{t},j}$ , respectively, are the bounce and transit frequencies of particle species $j$ . The relevant expansion parameter is $\epsilon = [v_{\textrm{th,i}}/(\omega L_z)]^2 \ll 1$ , where $v_{\textrm{th,i}} = \sqrt {2 T_{\textrm{i0}}/m_{\textrm{i}}}$ is the ion thermal velocity with ion temperature $T_{\textrm{i0}}$ and ion mass $m_{\textrm{i}}$ . Specific expressions required for the following derivation are presented in table 1, and further definitions are provided in table 2.

Table 1. Definitions of quantities in gyrokinetic KBM theory (Howes et al. 2006).

Table 2. Definitions of quantities in (2.11) and terms therein.

The linearized collisionless gyrokinetic equation (Rutherford & Frieman 1968) governing the evolution of the non-adiabatic part of the perturbed distribution function $h_{\!j}$ for particle species $j$ (i.e., the distribution of charged rings responding to the perturbed fields) is

(2.11) \begin{equation} v_\parallel \boldsymbol{\nabla} _\parallel \hat {h}_{\!j} - i(\omega - \omega _{\textrm{D}j}) \hat {h}_{\!j} = -i\frac {q_{\!j} F_{\textrm{M}j}}{T_{0j}} \left(\omega - \omega _{*j}^T\right) \langle \hat {\chi } \rangle _{\boldsymbol{R}_{\!j}}, \end{equation}

where the gyro-averaged potential is

(2.12) \begin{align} \langle \chi \rangle _{\boldsymbol{R}_{\!j}} &= \langle \hat {\chi } \rangle _{\boldsymbol{R}_{\!j}} \textrm{e}^{i(\boldsymbol{k} \boldsymbol{\cdot }\boldsymbol{R}_{\!j} - \omega t)} \nonumber \\ &= \left [ J_0 \left (\frac {k_\perp v_\perp }{\omega _{\textrm{c}j}} \right ) \left (\hat {\phi } - \frac {v_\parallel \hat {A}_\parallel }{c} \right ) + \frac {J_1\left ({k_\perp v_\perp }/{\omega _{\textrm{c}j}} \right )}{\left ({k_\perp v_\perp }/{\omega _{\textrm{c}j}} \right )} \frac {m v_\perp ^2}{q_{\!j}} \frac {\delta \hat {B}_\parallel }{B} \right ] \textrm{e}^{i(\boldsymbol{k} \boldsymbol{\cdot }\boldsymbol{R}_{\!j} - \omega t)}, \end{align}

the Maxwellian background distribution function takes the usual form

(2.13) \begin{equation} F_{\textrm{M}j} = n_0 \left (\frac {m_{\!j}}{2T_{0j}}\frac {1}{\pi } \right )^{3/2} e^{-E_{j}/T_{0j}} = n_0 \left (\frac {1}{v_{{\textrm{th},j}}^2 \pi } \right )^{3/2} e^{-\hat {v}^2}, \end{equation}

where hat quantities (except for velocity terms) denote the slowly varying amplitude of the fluctuating quantity in question, and $J_0$ and $J_1$ respectively are the zeroth- and first-order Bessel functions. The remaining quantities are defined in table 2. Using $v_\parallel \boldsymbol{\nabla} _\parallel h_{\!j} \to i k_\parallel v_\parallel h_{\!j}$ , (2.11) can be rearranged to obtain an expression for $\hat {h}_{\!j}$

(2.14) \begin{equation} \hat {h}_{\!j} = \frac {q_{\!j} F_{\textrm{M}j}}{T_{0j}} \left (\frac {\omega - \omega _{*j}^T}{\omega - \omega _{\textrm{D}j} - k_\parallel v_\parallel }\right ) \langle \hat {\chi } \rangle _{\boldsymbol{R}_{\!j}}. \end{equation}

The aforementioned simplifications are made from the outset, where the effects of $\delta B_\parallel$ and trapped particles are neglected. Given the intermediate-frequency regime (2.10), the mode frequency is assumed to be much larger than the ion transit frequency $|\omega | \gg k_\parallel v_{\textrm{th,i}}$ and much smaller than the electron transit frequency $|\omega | \ll k_\parallel v_{\textrm{th,e}}$ . The parallel current from the ions is found to be ignorable, and is largely carried by the electrons (Hirose, Zhang & Elia 1995). For the electrons, FLR effects can be ignored. The perturbed distribution functions for the ions and electrons respectively are then

(2.15) $$\begin{align}\delta f_{\textrm{i}} &= - \frac {q_{\textrm{i}}}{T_{\textrm{i0}}} \left [ F_{\textrm{Mi}} \phi - \hat {h}_{\textrm{i}} \textrm{e}^{i(\boldsymbol{k} \boldsymbol{\cdot }\boldsymbol{R}_{\textrm{i}} - \omega t)} \right ] \nonumber \\ &= - \frac {e}{T_{\textrm{i0}}} \left [ F_{\textrm{Mi}} \phi - F_{\textrm{Mi}} \left (\frac {\omega - \omega _{\textrm{*i}}^T}{\omega - \omega _{\textrm{Di}}}\right ) J_0 \left (\frac {k_\perp v_\perp }{\omega _{\textrm{ci}}} \right ) \left (\hat {\phi } - \frac {v_\parallel \hat {A}_\parallel }{c} \right ) \textrm{e}^{i(\boldsymbol{k} \boldsymbol{\cdot }\boldsymbol{R}_{\textrm{i}} - \omega t)} \right ],\end{align}$$

(2.16) $$\begin{align}\delta f_{\textrm{e}} &= - \frac {q_{\textrm{e}}}{T_{\textrm{e0}}} \left [ F_{\textrm{Me}} \phi - \hat {h}_{\textrm{e}} \textrm{e}^{i(\boldsymbol{k} \boldsymbol{\cdot }\boldsymbol{R}_{\textrm{e}} - \omega t)} \right ] \nonumber \\ &= \frac {e}{T_{\textrm{e0}}} \left [ F_{\textrm{Me}} \phi - F_{\textrm{Me}} \left (\frac {\omega - \omega _{\textrm{*e}}^T}{- k_\parallel v_\parallel }\right ) \left (\hat {\phi } - \frac {v_\parallel \hat {A}_\parallel }{c} \right ) \textrm{e}^{i(\boldsymbol{k} \boldsymbol{\cdot }\boldsymbol{R}_{\textrm{e}} - \omega t)} \right ].\end{align}$$

The perturbed densities for ions and electrons are obtained from $\delta n_{\!j} = \int \text{d}^3 v \, \delta f_{\!j}$ (where integration is performed at a fixed $\boldsymbol{r}$ , where the charges reside). Using

(2.17) \begin{equation} \int \text{d}^3 v \, h_{\!j} = \int \text{d}^3 v \, \langle h_{\!j} \rangle _{\boldsymbol{r}} = \int \text{d}^3 v \, \hat {h}_{\!j} \textrm{e}^{i(\boldsymbol{k \boldsymbol{\cdot }\boldsymbol{r}} - \omega t)} J_0 \left (\frac {k_\perp v_\perp }{\omega _{\textrm{ci}}} \right ) \end{equation}

leads to the perturbed ion density

(2.18) \begin{align} \delta n_{\textrm{i}} &= - \frac {e}{T_{\textrm{i0}}} \left [ \int \text{d}^3 v \, F_{\textrm{Mi}} \phi - \int \text{d}^3 v \, h_{\textrm{i}} \right ]\nonumber \\ &= - \frac {e}{T_{\textrm{i0}}} \left [ n_0 \phi - \hat {\phi } \textrm{e}^{i(\boldsymbol{k \boldsymbol{\cdot }\boldsymbol{r}} - \omega t)} \int \text{d}^3 v \,F_{\textrm{Mi}} J_0^2 \left (\frac {k_\perp v_\perp }{\omega _{\textrm{ci}}} \right ) \left (\frac {\omega - \omega _{\textrm{*i}}^T}{\omega - \omega _{\textrm{Di}}}\right ) \right ]\nonumber \\ &= - \frac {e n_0}{T_{\textrm{i0}}} (1-Q) \phi , \end{align}

where

(2.19) \begin{equation} Q \equiv \int \text{d}^3v \frac {F_{\textrm{Mi}}}{n_0} J_0^2(\hat {v}_\perp \sqrt {2b}) \left (\frac {\omega - \omega _{\textrm{*i}}^T}{\omega - \omega _{\textrm{Di}}} \right ), \end{equation}

with $b=(k_\perp \rho _{\textrm{th,i}})^2/2$ , and terms that are odd with respect to $v_\parallel$ vanish on integration.

Similarly, for the electrons

(2.20) \begin{align} \delta n_{\textrm{e}} &= \frac {e}{T_{\textrm{e0}}} \left [ \int \text{d}^3 v \, F_{\textrm{Me}} \phi - \int \text{d}^3 v \, h_{\textrm{e}} \right ] \nonumber \\ &= \frac {e}{T_{\textrm{e0}}} \left [ n_0 \phi - \left (\frac {\omega \hat {A}_\parallel }{k_\parallel c} \right )\textrm{e}^{i(\boldsymbol{k \boldsymbol{\cdot }\boldsymbol{r}} - \omega t)} \int \text{d}^3 v \, F_{\textrm{Me}} \left (1 - \frac {\omega _{\textrm{*e}}^T}{\omega }\right ) \right ] \nonumber \\ &\simeq \frac {e n_0}{T_{\textrm{e0}}} \left [ \phi - \left (1 - \frac {\omega _{\textrm{*e}}}{\omega }\right ) \left (\frac {\omega A_\parallel }{k_\parallel c} \right ) \right ]. \end{align}

Utilizing the ballooning transformation (Connor, Hastie & Taylor 1978) in a similar fashion as Hirose et al. (1995), one can relate

(2.21) \begin{equation} \boldsymbol{\nabla} _\parallel = i k_\parallel = \frac {1}{J\hat {B}} \frac {\partial }{\partial \theta }, \end{equation}

such that one can define (Tang et al. 1980)

(2.22) \begin{equation} \psi _\parallel \equiv \frac {\omega A_\parallel }{k_\parallel c} \implies \frac {1}{J\hat {B}} \frac {\partial }{\partial \theta } \psi _\parallel = \frac {i\omega }{c}A_\parallel , \end{equation}

where $\psi _\parallel$ is the opposite-parity representation of $A_\parallel$ .

With the perturbed densities of both particle species now obtained, they can be inserted into the quasi-neutrality condition

(2.23) \begin{equation} \sum \limits _{\!j} q_{\!j} \delta n_{\!j} = e \left (\delta n_{\textrm{i}} - \delta n_{\textrm{e}} \right ) = 0, \end{equation}

giving

(2.24) \begin{equation} \left ( 1 - \frac {\omega _{\textrm{*e}}}{\omega } \right ) \hat {\psi }_\parallel = (1 + \tau - \tau Q)\hat {\phi }, \end{equation}

where the potentials have been expanded as plane waves, leaving only the amplitudes (hat terms).

The parallel component of Ampère’s law involves taking the parallel-velocity moment of $h_{\!j}$

(2.25) \begin{equation} - \boldsymbol{\nabla} _\perp ^2 A_\parallel = \frac {4\pi }{c} \delta j_\parallel = \sum \limits _{\!j} \frac {4\pi }{c} q_{\!j} \int \text{d}^3 v \, v_\parallel \langle h_{\!j} \rangle _{\boldsymbol{r}}, \end{equation}

where the adiabatic part vanishes due to odd integration with respect to $v_\parallel$ . By considering only contributions from the electron current, this law reads (Tang et al. 1980)

(2.26) \begin{align} \frac {L_z^2}{J\hat {B}} \frac {\partial }{\partial \theta } \left ( \frac {b}{J\hat {B}} \frac {\partial }{\partial \theta } \hat {\psi }_\parallel \right ) &= \left (\frac {\omega }{\omega _{\textrm{A}}}\right )^2 \biggl \{ \left (Q - 1 - \frac {\omega _{\textrm{*e}}}{\omega \tau } \right ) \hat {\phi } \nonumber \\ &\quad - \left [ 1 - \frac {\omega _{\textrm{*e}}}{\omega } (1+\eta _{\textrm{e}}) \right ] \left ( \frac {\omega _\kappa + \omega _{\boldsymbol{\nabla }B}}{\omega } \right ) \hat {\psi }_\parallel \biggr \}, \end{align}

where $\omega _{\textrm{A}}^2 = v_{\textrm{th,i}}^2/(L_z^2 \beta _{\textrm{i}})$ is the Alfvén frequency, and $\beta _{\textrm{i}} = \beta _{\textrm{e}}/\tau$ for equal particle-species density $n_{\textrm{i0}} = n_{\textrm{e0}}$ .

Combining (2.24) and (2.26) yields

(2.27) \begin{align} \frac {L_z^2}{J\hat {B}} \frac {\partial }{\partial \theta } \left ( \frac {b}{J\hat {B}} \frac {\partial }{\partial \theta } \hat {\psi }_\parallel \right ) &= \left (\frac {\omega }{\omega _{\textrm{A}}}\right )^2 \biggl \{ \left (Q - 1 - \frac {\omega _{\textrm{*e}}}{\omega \tau } \right )(1 + \tau - \tau Q)^{-1}\left ( 1 - \frac {\omega _{\textrm{*e}}}{\omega } \right ) \nonumber \\ &\quad - \left [ 1 - \frac {\omega _{\textrm{*e}}}{\omega } (1+\eta _{\textrm{e}}) \right ] \left ( \frac {\omega _\kappa + \omega _{\boldsymbol{\nabla }B}}{\omega } \right ) \biggr \} \hat {\psi }_\parallel \nonumber \\ &= \left (\frac {\omega }{\omega _{\textrm{A}}}\right )^2 \left \{ \frac {Q - \alpha _{\textrm{0,i}}}{1 + \tau - \tau Q} \alpha _{\textrm{0,e}} - \alpha _{\textrm{1,e}} \left ( \frac {\omega _\kappa + \omega _{\boldsymbol{\nabla }B}}{\omega } \right ) \right \}\hat {\psi }_\parallel , \end{align}

where the condensed notation

(2.28) \begin{equation} \alpha _{m,j} = 1 - \frac {\omega _{*j}}{\omega }(1 + m \eta _{\!j}) \end{equation}

has been used.

The KBM eigenfunction $\hat {\psi }_\parallel$ can be rewritten in terms of a reduced scalar $\varPhi$ (Hirose et al. 1995) in the following manner. Rearranging (2.24) gives

(2.29) \begin{equation} \hat {\psi }_\parallel = \frac {(1 + \tau - \tau Q)}{\left ( 1 - \frac {\omega _{\textrm{*e}}}{\omega } \right )} \hat {\phi } = \frac {(1 + \tau - \tau Q)}{\alpha _{\textrm{0,e}}} \hat {\phi } \equiv \varPhi . \end{equation}

By making use of an approximative form (Aleynikova & Zocco 2017) of the ion integral $Q$ (for $\omega _{\textrm{Di}}/|\omega | \sim \epsilon \ll 1$ )

(2.30) \begin{align} Q &\approx 1 - \frac {\omega _{\textrm{*i}}}{\omega } + \left (\frac {\omega _{\boldsymbol{\nabla }B} + \omega _\kappa }{\omega } - b \right ) \left [1 - \frac {\omega _{\textrm{*i}}}{\omega }(1+\eta _i)\right ] + \mathcal{O}(\epsilon ^2) \nonumber \\ &\approx \alpha _{\textrm{0,i}} + \left (\frac {\omega _{\boldsymbol{\nabla }B} + \omega _\kappa }{\omega } - b \right ) \alpha _{\textrm{1,i}} \nonumber \\ &\approx \alpha _{\textrm{0,i}} + \mathcal{O}(\epsilon ), \end{align}

one obtains to zeroth order

(2.31) \begin{align} \varPhi &= \frac {(1 + \tau - \tau Q)}{\alpha _{\textrm{0,e}}} \hat {\phi } \approx \frac {(1 + \tau - \tau \alpha _{\textrm{0,i}})}{\alpha _{\textrm{0,e}}}\hat {\phi } = \hat {\phi }, \end{align}

such that $\varPhi$ is equivalent to $\hat {\phi }$ with corrections of order $\epsilon \ll 1$ .

The KBM equation is now rewritten in a dimensionless form that is suitable for numerical implementation, where changes in notation are detailed in table 3. Additional simplifications are made here: only eigenmodes centered at $k_x=0$ are considered, and due to the absence of parallel magnetic-field fluctuations, a common approximation of setting $\omega _{\boldsymbol{\nabla }B} \to \omega _\kappa$ is made, which can act to partially compensate for the lack of destabilization normally induced by $\delta B_\parallel$ (Zocco et al. 2015; Aleynikova et al. 2018).

Table 3. Definitions and normalizations for the dimensionless KBM equation (2.32).

Applying these changes, the KBM equation becomes

(2.32) \begin{equation} \frac {1}{\hat {J} \hat {B}} \frac {\partial }{\partial \theta } \left ( \frac {g^{yy}}{\hat {J}\hat {B}^3} \frac {\partial }{\partial \theta } \varPhi \right ) = \varOmega ^2 \frac {\beta _{\textrm{i}}\tau ^2}{2\epsilon _n^2} \left ( \frac {Q - \alpha _{\textrm{0,i}}}{1 + \tau - \tau Q} \alpha _{\textrm{0,e}} - \alpha _{\textrm{1,e}} \frac {2\varOmega _\kappa }{\varOmega }\right ) \varPhi . \end{equation}

This KBM equation differs from those presented in Tang et al. (1980) and Hirose, Zhang & Elia (1994) in several respects. Unlike the result of Tang et al. (1980), (2.32) omits the effect of parallel-magnetic-field fluctuations (and uses the aforementioned approximation $\omega _{\boldsymbol{\nabla }B} \to \omega _\kappa$ ), and the influence of trapped particles is neglected. The result of Hirose et al. (1994) includes similar physics as (2.32), but is only applicable to a circular tokamak geometry; (2.32) is suitable for general toroidal geometry, i.e., shaped tokamaks and stellarators. The physics contained in (2.32) can be understood as follows. The left-hand side captures stabilizing effects such as shear stabilization via $\partial g^{yy}(\theta )/\partial \theta$ ( $g^{yy}$ also appears as an argument of the Bessel functions in $Q$ , where it contributes to FLR damping) and magnetic-field-line-bending stabilization manifest via the degree of eigenmode localization along the field line, measured by $|\partial ^2 \varPhi /\partial \theta ^2|$ . When this equation is solved, these stabilizing effects are balanced by the drives on the right-hand side: the normalized plasma pressure $\beta _{\textrm{i}}$ ; normalized pressure gradients are contained in $\epsilon _n^{-1}$ , $Q$ and $\alpha _{m,j}$ ; the geometric drive of bad curvature is contained in $\varOmega _\kappa$ and is also manifest in $Q$ ; and lastly, wave–particle resonances associated with passing ions subject to the magnetic-curvature drift is provided by $Q$ . It is worth mentioning that the right-hand side of this equation is not purely destabilizing; for instance, kinetic resonances can act to dampen KBM growth rates away from marginality (shown in § 3) and $Q$ includes FLR-stabilization terms, see (2.30) and (2.48).

This equation can be evaluated using either a resonant framework, where the ion integral $Q$ is treated resonantly, thus retaining ion-magnetic-drift resonances (i.e., toroidal resonances), or by using a non-resonant framework where $Q$ is evaluated assuming $\omega _{\textrm{Di}}/|\omega | \ll 1$ , shown in (2.30). The resonant treatment is preferable if one wishes to accurately resolve the KBM’s behavior near marginality, whilst a non-resonant approximation is suitable when the KBM is strongly unstable (Aleynikova & Zocco 2017). The non-resonant KBM equation reduces to an iMHD ballooning-mode equation in the limit $k_y \to 0$ . For the resonant treatment, a simplified form of the resonant ion integral is derived here, in agreement with the result reported in Cheng (1982b ).

Utilizing notation from Cheng (1982b ), and in keeping with their convention of normalizing frequencies by the electron diamagnetic frequency $\omega _{\textrm{*e}}$

\begin{align*} F_0 &= F_{\textrm{Mi}}\frac {v_{\textrm{th,i}}^3}{n_0} = \frac {1}{\pi ^{3/2}}e^{-\hat {v}^2},\nonumber \\[5pt] \varOmega &\equiv \frac {\omega }{\omega _{\textrm{*e}}},\nonumber \end{align*}

(2.33) \begin{align} \varOmega _{\textrm{*i}}^T &\equiv \frac {\omega _{\textrm{*i}}^T}{\omega _{\textrm{*e}}} = \frac {\omega _{\textrm{*i}}}{\omega _{\textrm{*e}}} [1 + \eta _{\textrm{i}}(\hat {v}^2 - 3/2)] = -\frac {1}{\tau }[1 + \eta _{\textrm{i}}(\hat {v}^2 - 3/2)],\nonumber\\[5pt]\varOmega _{\textrm{si}} &\equiv \tau \big(\varOmega - \varOmega _{\textrm{*i}}^T\big) = \tau \varOmega + [1 + \eta _{\textrm{i}}(\hat {v}^2 - 3/2)],\nonumber \\[5pt] \varOmega _{\textrm{Di}} &\equiv \frac {\omega _{\textrm{Di}}}{\omega _{\textrm{*e}}}, \end{align}

such that $Q$ can be obtained from

(2.34) \begin{align} \tau Q &= \tau \int \text{d}^3v \frac {F_{\textrm{Mi}}}{n_0} J_0^2 \big(\sqrt {2b} \hat {v}_\perp \big) \frac {\omega _{\textrm{*e}}}{\omega _{\textrm{*e}}} \left (\frac {\omega - \omega _{\textrm{*i}}^T}{\omega - \omega _{\textrm{Di}}} \right ) \nonumber \\[4pt] &= \int \text{d}^3\hat {v} F_0 J_0^2 \big(\sqrt {2b} \hat {v}_\perp \big) \left (\frac {\varOmega _{\textrm{si}}}{\varOmega - \varOmega _{\textrm{Di}}} \right ), \end{align}

which aside from a different sign convention for $\varOmega _{\textrm{Di}}$ is equivalent to (18) in Cheng (1982b ).

To make this integral solvable analytically, the approximation is made of setting $\hat {v}_\perp ^2 = 2 \hat {v}_\parallel ^2$ in the denominator, which corresponds to a common curvature model in ITG theory (Terry, Anderson & Horton 1982); this simplification is made for KBM theory in Cheng (1982b ), which removes the resonance present in the $\hat {v}_\perp$ integral but retains it in the $\hat {v}_\parallel$ integral. As explained in Cheng (1982b ), this approximation deforms the resonant surface in velocity space into a plane, which provides analytical tractability without significantly affecting the underlying magnetic-curvature-drift-resonance effects.

The velocity-space integrals are expanded as

(2.35) \begin{equation} \int \text{d}^3 \hat {v} = 2\pi \int \limits _{-\infty }^{+\infty } \text{d} \hat {v}_\parallel \int \limits _{0}^{+\infty } \text{d} \hat {v}_\perp \hat {v}_\perp , \end{equation}

and the following relations will be utilized:

(2.36) \begin{align} \frac {1}{\sqrt {\pi }} \int \limits _{-\infty }^{+\infty } \text{d}x \frac {e^{-x^2}}{x-\zeta } \equiv Z(\zeta ) , \end{align}

(2.37) \begin{align} \frac {1}{i \pi } \int \limits _{-\infty }^{+\infty } \text{d}x \frac {e^{-x^2}}{x-\zeta } = \frac {1}{i\sqrt {\pi }} Z(\zeta ) \equiv W(\zeta ) , \end{align}

(2.38) \begin{align} \frac {2 \zeta }{i \pi } \int \limits _{0}^{+\infty } \text{d}x \frac {e^{-x^2}}{x^2-\zeta ^2} = \frac {\zeta }{i \pi } \int \limits _{-\infty }^{+\infty } \text{d}x \frac {e^{-x^2}}{x^2-\zeta ^2} \equiv W(\zeta ) , \end{align}

(2.39) \begin{align} \int \limits _{0}^{+\infty } \text{d} \hat {v}_\perp \hat {v}_\perp e^{-\hat {v}_\perp ^2} J_0^2 \big(\sqrt {2b} \hat {v}_\perp \big) = \frac {1}{2} \varGamma _0(b), \end{align}

(2.40) \begin{align} \int \limits _{0}^{+\infty } \text{d} \hat {v}_\perp \hat {v}_\perp ^3 e^{-\hat {v}_\perp ^2} J_0^2 \big(\sqrt {2b} \hat {v}_\perp \big) = \frac {1}{2} \left \{ \varGamma _0(b) + b[\varGamma _1(b)-\varGamma _0(b)] \right \}, \end{align}

where (2.36) is the plasma dispersion function, (2.37) and (2.38) define the Faddeeva function and $\varGamma _n(b) \equiv I_n(b) e^{-b}$ , where $I_n(b)$ is the modified Bessel function of the first kind. Note that the integral representations for $Z(\zeta )$ and $W(\zeta )$ are evaluated only for $\textrm {Im}(\zeta )\gt 0$ .

Using these definitions, one may write

(2.41) \begin{align} \tau Q &= \frac {2}{\sqrt {\pi }} \int \limits _{-\infty }^{+\infty } \text{d} \hat {v}_\parallel \left (\frac {e^{-\hat {v}_\parallel ^2}}{\varOmega - 4 \varOmega _\kappa \hat {v}_\parallel ^2} \right )\nonumber \\[3pt] &\quad \times \int \limits _{0}^{+\infty } \text{d} \hat {v}_\perp \hat {v}_\perp \left \{ e^{-\hat {v}_\perp ^2} J_0^2 \big(\sqrt {2b} \hat {v}_\perp \big) \left [\left(\tau \varOmega + 1 + \eta _{\textrm{i}}\hat {v}_\parallel ^2 - \eta _{\textrm{i}}3/2\right) + \eta _{\textrm{i}}\hat {v}_\perp ^2 \right ] \right \}. \end{align}

Focusing only on the $\hat {v}_\perp$ integral

(2.42) \begin{align} &\int \limits _{0}^{+\infty } \text{d} \hat {v}_\perp \hat {v}_\perp \left \{ e^{-\hat {v}_\perp ^2} J_0^2 \big(\sqrt {2b} \hat {v}_\perp \big) \left [(\tau \varOmega + 1 + \eta _{\textrm{i}}\hat {v}_\parallel ^2 - \eta _{\textrm{i}}3/2) + \eta _{\textrm{i}}\hat {v}_\perp ^2 \right ] \right\} \nonumber \\[3pt] &\quad {} = \frac {1}{2}\left\{\eta _{\textrm{i}}\hat {v}_\parallel ^2 \varGamma _0(b) + (\tau \varOmega + 1 - \eta _{\textrm{i}}/2) \varGamma _0(b) + \eta _{\textrm{i}} b[\varGamma _1(b)-\varGamma _0(b)]\right\}. \end{align}

Reinserting this result into the full expression, one arrives at

(2.43) \begin{align} \tau Q &= \frac {2}{\sqrt {\pi }} \int \limits _{-\infty }^{+\infty } \text{d} \hat {v}_\parallel \left (\frac {e^{-\hat {v}_\parallel ^2}}{\varOmega - 4 \varOmega _\kappa \hat {v}_\parallel ^2} \right )\nonumber \\[3pt] &\quad \times \frac {1}{2}\left \{\eta _{\textrm{i}}\hat {v}_\parallel ^2 \varGamma _0(b) + (\tau \varOmega + 1 - \eta _{\textrm{i}}/2) \varGamma _0(b) + \eta _{\textrm{i}} b[\varGamma _1(b)-\varGamma _0(b)]\right \} \nonumber \\[3pt] &= \frac {1}{\sqrt {\pi }} \eta _{\textrm{i}} \varGamma _0(b) \int \limits _{-\infty }^{+\infty } \text{d} \hat {v}_\parallel \hat {v}_\parallel ^2 \left (\frac {1}{-4 \varOmega _\kappa } \right ) \left (\frac {e^{-\hat {v}_\parallel ^2}}{\hat {v}_\parallel ^2 -\varOmega /4\varOmega _\kappa } \right ) \nonumber \\[3pt] &\quad + \frac {1}{\sqrt {\pi }} G \int \limits _{-\infty }^{+\infty } \text{d} \hat {v}_\parallel \left (\frac {1}{-4 \varOmega _\kappa } \right ) \left (\frac {e^{-\hat {v}_\parallel ^2}}{\hat {v}_\parallel ^2 -\varOmega /4\varOmega _\kappa } \right ), \end{align}

where

(2.44) \begin{equation} G \equiv (\tau \varOmega + 1 - \eta _{\textrm{i}}/2) \varGamma _0(b) + \eta _{\textrm{i}} b[\varGamma _1(b)-\varGamma _0(b)]. \end{equation}

By introducing the relabeling $\xi ^2 \equiv \varOmega /(4\varOmega _\kappa )$ , (2.43) can be written as

(2.45) \begin{align} \tau Q &= \frac {1}{\sqrt {\pi }} \eta _{\textrm{i}} \varGamma _0(b) \left (\frac {1}{-4 \varOmega _\kappa } \right ) \int \limits _{-\infty }^{+\infty } \text{d} \hat {v}_\parallel \hat {v}_\parallel ^2 \left (\frac {e^{-\hat {v}_\parallel ^2}}{\hat {v}_\parallel ^2 - \xi ^2} \right ) \nonumber \\ &\quad + \frac {1}{\sqrt {\pi }} G \left (\frac {1}{-4 \varOmega _\kappa } \right ) \int _{-\infty }^{+\infty } \text{d} \hat {v}_\parallel \left (\frac {e^{-\hat {v}_\parallel ^2}}{\hat {v}_\parallel ^2 - \xi ^2} \right ). \end{align}

Now utilizing the earlier defined Faddeeva function $W(\xi ) \equiv ({1}/{i\sqrt {\pi }}) Z(\xi )$

(2.46) \begin{align} \frac {Z(\xi )}{\xi } &= \frac {1}{\sqrt {\pi }} \int \limits _{-\infty }^{+\infty } \text{d}\hat {v}_\parallel \frac {e^{-\hat {v}_\parallel ^2}}{\hat {v}_\parallel ^2-\xi ^2}, \end{align}

and similarly

(2.47) \begin{align} \frac {1}{\sqrt {\pi }} \int \limits _{-\infty }^{+\infty } \text{d} \hat {v}_\parallel \hat {v}_\parallel ^2 \left (\frac {e^{-\hat {v}_\parallel ^2}}{\hat {v}_\parallel ^2 - \xi ^2} \right ) &= \frac {1}{\sqrt {\pi }} \int \limits _{-\infty }^{+\infty } \text{d} \hat {v}_\parallel (\hat {v}_\parallel ^2 - \xi ^2 + \xi ^2) \left (\frac {e^{-\hat {v}_\parallel ^2}}{\hat {v}_\parallel ^2 - \xi ^2} \right ) \nonumber \\ &= 1 + \xi Z(\xi ). \end{align}

Reinserting (2.47) into the full expression, the ion integral reduces to

(2.48) \begin{align} \tau Q &= \left \{\left (\tau \varOmega + 1 - \frac {\eta _{\textrm{i}}}{2} \right ) \varGamma _0(b) + \eta _{\textrm{i}} b [\varGamma _1(b) - \varGamma _0(b)] \right \} D_0 + \eta _{\textrm{i}} \varGamma _0 D_1, \end{align}

with

(2.49) $$\begin{align}D_0 &= -\frac {1}{4 \varOmega _\kappa } \frac {Z(\xi )}{\xi },\end{align}$$

(2.50) $$\begin{align}D_1 &= -\frac {1}{4 \varOmega _\kappa } [1 + \xi Z(\xi )].\end{align}$$

This final result is identical to (18) of Cheng (1982b ) and can be directly inserted into the KBM equation (2.32) for evaluation. Alternatively, the full unapproximated integral, (2.19), can be numerically evaluated from the outset by utilizing Gauss–Laguerre and Gauss–Hermite quadrature schemes for the velocity-space integration (Hirose et al. 1995); results from both analytical and numerical evaluations are presented in § 3. Another analytical and higher-fidelity form of $Q$ (without the simplifying approximation $\hat {v}_\perp ^2 = 2 \hat {v}_\parallel ^2$ in the denominator of $Q$ ) has been derived for ITG theory in Zocco et al. (2018b ); this result has been found to agree exactly with the full unapproximated integral (2.19) when applied to a circular tokamak (not shown here). Further refinement is required regarding the numerical treatment of this latter analytical approach in order to reduce its computation time and expand its usage to stellarator geometry.

2.2. The KBM eigenvalue problem

The dimensionless KBM eigenmode equation (2.32) is now reformulated as an eigenvalue problem amenable to straightforward numerical evaluation. The left-hand side of (2.32) is rewritten as follows:

(2.51) \begin{equation} \frac {1}{\hat {J} \hat {B}} \frac {\partial }{\partial \theta } \left ( \frac {g^{yy}}{\hat {J}\hat {B}^3} \frac {\partial }{\partial \theta } \varPhi \right ) = \frac {1}{\hat {J} \hat {B}} \left [\frac {\partial }{\partial \theta }\left (\frac {g^{yy}}{\hat {J}\hat {B}^3} \right )\frac {\partial }{\partial \theta } + \frac {g^{yy}}{\hat {J}\hat {B}^3} \frac {\partial ^2}{\partial \theta ^2} \right ] \varPhi \equiv \mathcal{L} \varPhi , \end{equation}

where $\mathcal{L}$ is an operator acting on $\varPhi$ . Similarly, for the right-hand side

(2.52) \begin{equation} \varOmega ^2 \frac {\beta _{\textrm{i}}\tau ^2}{2\epsilon _n^2} \left ( \frac {Q - \alpha _{\textrm{0,i}}}{1 + \tau - \tau Q} \alpha _{\textrm{0,e}} - \alpha _{\textrm{1,e}} \frac {2\varOmega _\kappa }{\varOmega } \right ) \varPhi \equiv \mathcal{R} \varPhi , \end{equation}

such that the KBM equation can be written compactly as

(2.53) \begin{equation} (\mathcal{L} - \mathcal{R}) \varPhi \equiv \mathcal{M} \varPhi = 0. \end{equation}

Given the nonlinear dependence of $\mathcal{M}$ on $\varOmega$ , this system can be treated as a nonlinear eigenvalue problem $\mathcal{M}(\varOmega )\varPhi = 0$ , where one seeks to calculate the eigenvalue–eigenvector pairs $(\varOmega , \varPhi )$ . As will be shown, the eigenvalues $\varOmega$ are obtained first, and are then utilized to obtain the eigenvectors $\varPhi$ .

To solve the system of linear equations represented by $\mathcal{M}(\varOmega )\varPhi = 0$ , non-trivial solutions for the eigenvector $\varPhi$ can be obtained if

(2.54) \begin{equation} \textrm{Det}[\mathcal{M}(\varOmega )] = 0. \end{equation}

Thus, eigenvalues $\varOmega$ are sought that make $\mathcal{M}(\varOmega )$ singular. In practice, determinant calculations can be susceptible to significant numerical instability, producing large-magnitude determinants even for singular matrices. As such, a more reliable method is to use condition numbers, where a singular matrix $\mathcal{M}(\varOmega )$ has condition number $\textrm{Cond}[\mathcal{M}(\varOmega )]=\infty$ . This is the approach applied here. Note that, although $\textrm{Cond}[\mathcal{M}(\varOmega )]=\infty$ is not strictly achieved, acceptable solutions – justified by the quality of the results obtained a posteriori – can be obtained by identifying $\varOmega$ that maximize $|\textrm{Cond}[\mathcal{M}(\varOmega )]|$ .

The procedure employed here is as follows. To obtain eigenvalues $\varOmega$ that satisfy (2.54), $|\textrm{Cond}[\mathcal{M}(\varOmega )]|$ is evaluated at each point on a chosen eigenvalue grid (Re $\varOmega$ , Im $\varOmega$ ), where the solutions $\varOmega$ that maximize $|\textrm{Cond}[\mathcal{M}(\varOmega )]|$ satisfy

(2.55) $$\begin{align}\frac{\text{d}}{{\rm d}\varOmega } |\textrm{Cond}[\mathcal{M}(\varOmega )]| &= 0,\end{align}$$

(2.56) $$\begin{align}\frac{\text{d}^2}{{\rm d}\varOmega ^2} |\textrm{Cond}[\mathcal{M}(\varOmega )]| &\lt 0,\end{align}$$

where the second condition ensures $\varOmega$ is at a local maximum. A new grid is then constructed in the vicinity of this previously found $\varOmega$ , where the mesh is made increasingly fine grained. The maximum value of $|\textrm{Cond}[\mathcal{M}(\varOmega )]|$ near the eigenvalue solution increases with successive grid refinement, as the calculation is able to converge increasingly to the true solution. These iterations continue until a convergence criterion is met between consecutive eigenvalue results $\varOmega$ (e.g., until $|\varOmega _n - \varOmega _{n-1}| \lt 10^{-10}$ ). Thus, although $\textrm{Cond}[\mathcal{M}(\varOmega )]=\infty$ is not reached in practice, one can get arbitrarily close to the true solution of $\varOmega$ with continuous grid refinement, where each iteration allows for higher precision of the resulting $\varOmega$ .

The converged results for $\varOmega$ are the eigenvalue solutions to the simplified KBM equation. Once these eigenvalues are obtained, all terms in $\mathcal{M}$ are known. To obtain the corresponding eigenvector $\varPhi$ , one can exploit the fact that the nonlinear eigenvalue problem $\mathcal{M}(\varOmega ) \varPhi = 0$ is a special case of a simpler linear eigenvalue problem, which has eigenvector solutions $\varPhi$ and eigenvalues $\lambda _{\mathcal{M}(\varOmega )}$

(2.57) \begin{equation} \mathcal{M}(\varOmega ) \varPhi = \lambda _{\mathcal{M}(\varOmega )} \varPhi \hspace {0.5cm} \text{(for $\lambda _{\mathcal{M}(\varOmega )}=0$)}, \end{equation}

such that the solution $\varPhi$ of the nonlinear eigenvalue problem is identical to the solution $\varPhi$ of the linear eigenvalue problem with eigenvalue $\lambda _{\mathcal{M}(\varOmega )}=0$ . This linear eigenvalue problem can be solved numerically and involves calculating the eigenvectors of $\mathcal{M}(\varOmega )$ , and selecting the solution with corresponding eigenvalue $\lambda _{\mathcal{M}(\varOmega )} \approx 0$ ; this provides one with the KBM eigenfunction $\varPhi ^{\textrm{KBM}}$ .

3. Numerical evaluation of KBM theory

In this section, the simplified KBM theory outlined in previous sections is numerically evaluated and compared with gyrokinetic flux-tube simulations performed with the Gene code (Jenko et al. 2000), including cases that were first presented in Mulholland et al. (2023, 2025). The numerical implementation of this theory is referred to as the KBM-eigenvalue-yielder (Key) code. The Key code requires only flux-tube geometry and plasma parameters as inputs, and computes self-consistent KBM eigenvalues and eigenvectors for the eigenmode along the field line. Eigenvalues (eigenvectors) are calculated in $\mathcal{O}(1)$ ( $\mathcal{O}(10)$ ) CPU seconds, respectively, where the computation time depends mostly on the field-line resolution.

In the following subsections, Key and Gene results are compared for a sub-threshold KBM study in W7-X (Mulholland et al. 2023, 2025), followed by both high and low shear ( $\hat {s} \approx 0.8$ and $\hat {s} = 0.1$ , where $\hat {s} = (-r/\iota)\textrm{d}\iota/\textrm{d}r$ ) tokamak scenarios. In the stellarator case, results from Gene include the effects of trapped particles and parallel magnetic-field fluctuations $\delta B_\parallel$ (where the $\boldsymbol{\nabla }B$ -drift and curvature-drift frequencies are treated distinctly in conjunction with a finite pressure gradient (Aleynikova et al. 2022)); these effects are absent in Key, and the $\boldsymbol{\nabla }B$ -drift frequency is replaced by the curvature-drift frequency ( $\omega _{\boldsymbol{\nabla }B} \to \omega _\kappa$ ) to partially compensate for the absence of $\delta B_\parallel$ (Zocco et al. 2015; Aleynikova et al. (2022)). In addition, the Gene results for W7-X utilize a flux-tube geometry corresponding to one poloidal turn (for which numerical convergence was achieved), where the field-line domain is extended artificially by exploiting the flux-tube parallel boundary condition (employing a finite number of $k_x$ modes) (Merz 2008); this approach is more numerically efficient but can lead to notable deviations in stellarator geometry between the artificially extended domain and a flux tube with a higher number of poloidal turns (Faber et al. 2018). Results from Key are obtained using flux-tube geometry with 20 poloidal turns and thus do not use artificially extended flux-tube domains. In the tokamak cases, results from Gene are obtained utilizing mostly the same assumptions as Key of neglecting effects from trapped particles and $\delta B_\parallel$ (using $\omega _{\boldsymbol{\nabla }B} \to \omega _\kappa$ ). In the tokamak, artificially extending the flux-tube domain is equivalent to physical extensions using more poloidal turns. In all cases, Landau resonances are absent in Key, but remain in Gene.

Results from Key include non-resonant (obtained using (2.30)), analytic-resonant (partial toroidal resonance – obtained using (2.48)) and numeric-resonant (full toroidal resonance – obtained using (2.19)) eigenvalue (EV) solutions.

3.1. Sub-threshold KBMs in Wendelstein 7-X

In this section, Key is applied to an stKBM study carried out in stellarator geometry (Mulholland et al. 2023, 2025), using an MHD-optimized high-mirror configuration of W7-X denoted as KJM (Aleynikova et al. 2018). For all $\beta$ , the equilibria are varied consistently with the pressure gradient. The normalized gradients are $a/L_{T\textrm{i}} = 3.5$ , $a/L_{T\textrm{e}} = 0$ , $a/L_{n\textrm{i}} = a/L_{n\textrm{e}} = 1$ , where $a$ is the minor radius and $L_{Tj} = -({\rm d}\ln T_{\!j}/{\rm d}r)^{-1}$ and $L_{nj} = -({\rm d}\ln n_{\!j}/{\rm d}r)^{-1}$ are the respective length scales of the temperature and density for particle species $j$ .

Figure 1. Eigenvalues (growth rates $\gamma$ and real frequencies $\omega _{\textrm{r}}$ ) for (st)KBMs with increasing $\beta$ in W7-X KJM geometry at wavenumber $k_y \rho _{\textrm{s}} = 0.2$ , comparing Gene and Key. Key curves include results from non-resonant (green triangles, dotted lines), analytic-resonant (partial toroidal resonance; red squares, dashed lines) and numeric-resonant (full toroidal resonance; magenta diamonds, dashed lines) EV solutions. Results from Gene (blue circles, solid lines) include effects of trapped particles and $\delta B_\parallel$ , which are absent in Key. For the non-resonant approach, stable KBMs are obtained until reaching the predicted conventional-KBM regime at $\beta _{\textrm{crit}}^{\textrm{KBM}} \approx 3.3\,\%$ , beyond which the mode grows rapidly; both resonant approaches detect low- $\beta$ destabilization of the stKBM at $\beta _{\textrm{crit}}^{\textrm{stKBM}} \approx 1\,\%$ and the mode’s approximately constant growth rate until reaching $\beta _{\textrm{crit}}^{\textrm{KBM}} \approx 3.3\,\%$ . All solutions return similar frequencies with $\beta$ ; beyond $\beta _{\textrm{crit}}^{\textrm{KBM}} \approx 3.3\,\%$ , the non-resonant approach recovers the theoretically predicted frequency for strongly driven KBMs of $\omega _{\textrm{r}} = \omega _{\textrm{*pi}}/2$ .

Figure 1 shows growth rates $\gamma$ and real frequencies $\omega _{\textrm{r}}$ with increasing $\beta$ for the (st)KBM from Key and Gene at normalized binormal wavenumber $k_y \rho _{\textrm{s}} = 0.2$ . Focus is given to low binormal wavenumbers as KBM studies in W7-X have shown that growth rates peak at $k_y \rho _{\textrm{s}} \to 0$ (Aleynikova et al. 2018, 2022; Mulholland et al. 2025). Results from Key include non-resonant, analytic-resonant and numeric-resonant EV solutions. In both Gene and Key, two distinct eigenmode clusters are obtained corresponding to KBMs and electromagnetic ITGs; focus is given here only to the KBM cluster. The stKBM, once rendered unstable at $\beta _{\textrm{crit}}^{\textrm{stKBM}} \approx 1\,\%$ , is a single mode within the KBM cluster that remains subdominant ( $\gamma \lt \gamma _{\textrm{max}}$ ) until $\beta _{\textrm{crit}}^{\textrm{KBM}} \approx 3.3\,\%$ , from which point it grows steeply with $\beta$ as a conventional KBM (Mulholland et al. 2023, 2025). In Key, the non-resonant approach finds that the mode is stable until surpassing $\beta = 3\,\%$ , whereas the resonant approaches show instability from as low as $\beta \approx 1\,\%$ . Thus, when accounting for ion-magnetic-drift resonances, Key is shown to capture both the unconventional characteristics of the stKBM at low $\beta$ as well as the more typical behavior expected from the conventional KBM at higher $\beta$ . Specifically, Key reproduces the low- $\beta$ destabilization at $\beta _{\textrm{crit}}^{\textrm{stKBM}} \approx 1\,\%$ and the mode’s approximately constant growth rate until reaching the conventional-KBM regime at $\beta _{\textrm{crit}}^{\textrm{KBM}} \approx 3.3\,\%$ , where the iMHD limit is reached and the mode grows rapidly with $\beta$ . All approaches return similar but distinct frequency trends with $\beta$ ; beyond $\beta _{\textrm{crit}}^{\textrm{KBM}} \approx 3.3\,\%$ , the non-resonant results closely correspond to the theoretically predicted value for strongly driven KBMs (Tang et al. 1980) of $\omega _{\textrm{r}} = \omega _{\textrm{*pi}}/2$ , where $\omega _{\textrm{*pi}} = \omega _{\textrm{*i}}(1+\eta _{\textrm{i}})$ is the ion diamagnetic frequency associated with the ion pressure gradient. A noteworthy feature is the common frequency trend for all curves (Gene, non-resonant and resonant) with increasing $\beta$ ; the non-resonant results aid in explaining this phenomenon. For $\beta \lt 3\,\%$ , where the KBM is stable in the non-resonant picture, Key returns a pair of real solutions that are symmetric about $\omega _{\textrm{r}} = \omega _{\textrm{*pi}}/2$ (only one solution is shown here); once criticality is reached at $\beta _{\textrm{crit}}^{\textrm{KBM}} \approx 3.3\,\%$ , these two real solutions transition to a complex-conjugate pair with a common real part of $\omega _{\textrm{r}} = \omega _{\textrm{*pi}}/2$ and imaginary parts $\pm i\gamma$ , corresponding to one unstable and one stable branch, respectively (only the unstable branch is shown here). Therefore, the similarities seen between the non-resonant and resonant results may indicate that the introduction of resonances mainly acts to augment the non-resonant picture, but does not drastically alter it. Regardless, accounting for the differences brought about by ion-magnetic-drift resonances – although subtle – is vital, given that their inclusion explains the deleterious manifestation of the stKBM and its secondary effect on ITG heat fluxes.

Figure 2. Eigenfunctions $\textrm{Re}\,\varPhi$ for the (st)KBM versus ballooning angle $\theta$ at different $\beta$ in W7-X KJM geometry at wavenumber $k_y \rho _{\textrm{s}} = 0.2$ , comparing Gene (blue solid lines) and key (red dashed lines), where Key eigenfunctions correspond to the analytic-resonant eigenvalues in figure 1. Results from Gene include effects of trapped particles and $\delta B_\parallel$ , which are absent in Key. Key captures the broad mode structure of the stKBM for $\beta \approx 1\,\%{-}3\,\%$ in the complex geometry of W7-X, as well as its progressive narrowing with increasing drive as it transitions to the conventional KBM for $\beta \gt 3\,\%$ .

Figure 2 shows eigenfunctions of the real part of the perturbed electrostatic potential $\textrm{Re}\,\varPhi$ versus ballooning angle $\theta$ from Gene and Key for the stKBM ( $\beta \approx 1\,\%{-}3\,\%$ ) and conventional KBM ( $\beta \gt 3\,\%$ ) at normalized binormal wavenumber $k_y \rho _{\textrm{s}} = 0.2$ . Eigenfunctions from Key correspond to the analytic-resonant EV approach; other EV approaches (non-resonant and numeric resonant) return similar eigenfunctions to those shown here. Eigenfunctions from both resonant approaches in Key have $\textrm{Re}\,\varPhi \gg \textrm{Im}\,\varPhi$ , while in the non-resonant approach, $\textrm{Im}\,\varPhi = 0$ . Only $\textrm{Re}\,\varPhi$ is shown here to highlight the positive–negative oscillations of the stKBM eigenfunction at $|\theta | \neq 0$ for low $\beta$ observed in both Gene and Key. Key is shown to capture the broad mode structure of the stKBM in the complex geometry of W7-X, as well as its progressive narrowing with increasing drive as it transitions to the conventional KBM. The shape of the stKBM eigenfunction is primarily determined by the field-line curvature $\mathcal{K}_y$ in W7-X; specifically, the high spatial variation of bad-curvature wells in W7-X gives the broad stKBM eigenfunction its modulated amplitude along the field line, where the eigenmode peaks at locations where bad-curvature wells coincide with regions of lower $g^{yy}$ , where the latter provides the stabilizing influence of magnetic shear and FLR suppression. Discrepancies between the results from Gene and Key are likely attributed to differences in the underlying assumptions of each code combined with deviations in flux-tube geometry for $\theta /\pi \gt 1$ , i.e., beyond one poloidal turn, Gene uses artificial extensions of the geometry whereas Key uses a longer flux tube of 20 poloidal turns. Improved agreement is expected if the same underlying assumptions are made and if identical flux-tube domains are employed by both codes.

3.2. High-/low-shear tokamak

Key is now applied to a tokamak scenario adapted from the Cyclone Base Case (Dimits et al. 2000), where both nominal background magnetic shear ( $\hat {s} \approx 0.8$ ) and a reduced shear ( $\hat {s} = 0.1$ ) are considered. In both cases, the equilibria are varied consistently with the pressure gradient, trapped particles and $\delta B_\parallel$ are neglected, the $\boldsymbol{\nabla }B$ -drift frequency is replaced by the curvature-drift frequency ( $\omega _{\boldsymbol{\nabla }B} \to \omega _\kappa$ ) and the normalized gradients are $R/L_{T\textrm{i}} = R/L_{T\textrm{e}} = R/L_{n\textrm{i}} = R/L_{n\textrm{e}} = 10$ , where $R$ is the major radius.

Figure 3. Eigenvalues (growth rates $\gamma$ and real frequencies $\omega _{\textrm{r}}$ ) for the KBM as functions of $\beta$ in tokamak geometry with $\hat {s} \approx 0.8$ at wavenumber $k_y \rho _{\textrm{s}} = 0.1$ , comparing Gene and Key. Key curves include results from non-resonant (green triangles, dotted lines), analytic-resonant (partial toroidal resonance; red squares, dashed lines) and numeric-resonant (full toroidal resonance; magenta diamonds, dashed lines) EV solutions. Growth rates from Key show fair but not perfect agreement with Gene (blue circles, solid lines); resonant results follow Gene more consistently than non-resonant results. For the frequencies, the resonant approaches follow the trend of Gene, while the non-resonant approach recovers the theoretically predicted value for strongly driven KBMs of $\omega _{\textrm{r}} = \omega _{\textrm{*pi}}/2$ . Resonant results from Key recover $\beta _{\textrm{crit}}^{\textrm{KBM}} \approx 0.5\,\%$ , in agreement with Gene. The increasing stabilization seen for $\beta \gt 2\,\%$ is well captured in Key.

Figure 3 shows growth rates $\gamma$ and real frequencies $\omega _{\textrm{r}}$ with increasing $\beta$ from Key and Gene for the KBM with $\hat {s} \approx 0.8$ at normalized binormal wavenumber $k_y \rho _{\textrm{s}} = 0.1$ . Results from Key include non-resonant, analytic-resonant and numeric-resonant EV solutions. In all approaches, growth rates from Key show qualitative agreement with Gene (blue circles, solid lines) in terms of capturing the overall trend with increasing $\beta$ . For the frequencies, the resonant approaches approximately follow the trend of Gene, while the non-resonant approach recovers the theoretically predicted value for strongly driven KBMs of $\omega _{\textrm{r}} = \omega _{\textrm{*pi}}/2$ . In alignment with this, frequencies from the resonant approaches and Gene tend toward this strong-drive frequency $\omega _{\textrm{r}} \to \omega _{\textrm{*pi}}/2$ as the growth rate peaks at $\beta \approx 1.5\,\%{-}2\,\%$ . Both resonant approaches in Key obtain the same destabilization threshold as Gene of $\beta _{\textrm{crit}}^{\textrm{KBM}} \approx 0.5\,\%$ , whereas the non-resonant approach overestimates this threshold to be $\beta _{\textrm{crit}}^{\textrm{KBM}} \approx 0.6\,\%$ . The distinct feature of increasing stabilization of the KBM at high $\beta \gt 2\,\%$ is well captured in all approaches.

Figure 4. Eigenfunctions $|\varPhi |$ for the KBM versus ballooning angle $\theta$ at different $\beta$ in tokamak geometry with $\hat {s} \approx 0.8$ at wavenumber $k_y \rho _{\textrm{s}} = 0.1$ , comparing Gene (blue solid lines) with Key (red dashed lines), where Key eigenfunctions correspond to the numeric-resonant eigenvalues in figure 3. Good agreement is found for all $\beta$ .

Figure 4 shows eigenfunctions of the perturbed electrostatic potential $\varPhi$ with ballooning angle $\theta$ from Gene and Key for the KBM with $\hat {s} \approx 0.8$ at normalized binormal wavenumber $k_y \rho _{\textrm{s}} = 0.1$ at different $\beta$ ; these eigenfunctions correspond to the numeric-resonant eigenvalues shown in figure 3. Near marginality, both codes obtain a centrally peaked eigenfunction with decreasing amplitude at higher $\theta$ . With increasing $\beta \gt 1\%$ , the mode narrows slightly as the growth rate increases sharply. At higher $\beta \gt 2\,\%$ , the growth rate reduces and the eigenfunction broadens into a bell shape, present in both Gene and Key. Good agreement is found between both codes for all $\beta$ .

Figure 5. Eigenvalues (growth rates $\gamma$ and real frequencies $\omega _{\textrm{r}}$ ) for the KBM as functions of $\beta$ in tokamak geometry with $\hat {s} \approx 0.1$ at wavenumber $k_y \rho _{\textrm{s}} = 0.1$ , comparing Gene and Key. Key curves include results from non-resonant (green triangles, dotted lines), analytic-resonant (partial toroidal resonance; red squares, dashed lines) and numeric-resonant (full toroidal resonance; magenta diamonds, dashed lines) EV solutions. Growth rates from Key show qualitative agreement with Gene (blue circles, solid lines); resonant results follow Gene more consistently than non-resonant results. For the frequencies, the resonant approaches follow the trend of Gene, while the non-resonant approach recovers the theoretically predicted value for strongly driven KBMs of $\omega _{\textrm{r}} = \omega _{\textrm{*pi}}/2$ . Key slightly overestimates $\beta _{\textrm{crit}}^{\textrm{KBM}}$ at $0.1\,\%{-}0.3\,\%$ , compared with Gene’s $\beta _{\textrm{crit}}^{\textrm{KBM}} \approx 0.05\,\%$ . The modulating $\gamma$ with increasing $\beta$ is well captured in Key.

Figure 5 shows growth rates $\gamma$ and real frequencies $\omega _{\textrm{r}}$ with increasing $\beta$ from Key and Gene for the KBM with $\hat {s} \approx 0.1$ at normalized binormal wavenumber $k_y \rho _{\textrm{s}} = 0.1$ . Results from Key include non-resonant, analytic-resonant and numeric-resonant EV solutions. In all approaches, Key captures the modulating but generally increasing $\gamma$ with $\beta$ ; this modulation is caused by the interchanging of eigenmodes for the dominant position (largest $\gamma$ ) as $\beta$ increases, where the results shown here correspond only to the dominant mode at any given $\beta$ . It is noteworthy to mention that Key, like Gene, can recover subdominant modes. The modulating behavior observed here resembles a previous study of KBMs in a low-shear tokamak, reported in Hirose et al. (1995). For the frequencies, the resonant approaches approximately follow the trend of Gene, while the non-resonant approach recovers the theoretically predicted value for strongly driven KBMs of $\omega _{\textrm{r}} = \omega _{\textrm{*pi}}/2$ . Agreement between Gene and Key improves with the increasing physics fidelity of each approach (i.e., non-resonant $\to$ partial toroidal resonance $\to$ full toroidal resonance). The non-resonant approach overestimates the growth rate for $\beta \gt 0.4\,\%$ , and deviations increase with $\beta$ . The analytic-resonant approach improves upon the former, but consistently underestimates the growth rate for all $\beta$ . Lastly, the numeric-resonant approach obtains the best agreement with Gene, although it slightly underestimates the growth rate for $\beta \gt 0.5\,\%$ . All approaches in Key slightly overestimate the destabilization threshold to be $\beta _{\textrm{crit}}^{\textrm{KBM}} \approx 0.1\,\%{-}0.3\,\%$ , compared with Gene’s estimated extrapolated threshold of $\beta _{\textrm{crit}}^{\textrm{KBM}} \approx 0.05\,\%$ . The relative difference between these thresholds is substantial, but is amplified by the very low $\beta$ here; in terms of absolute differences, these deviations are minor. In terms of comparing machine performance, Key accurately detects the notable differences in destabilization thresholds and growth-rate trends observed in the high- and low-shear tokamaks presented here.

Figure 6. Eigenfunctions  $|\varPhi |$ for the KBM versus ballooning angle $\theta$ at different $\beta$ in tokamak geometry with $\hat {s} \approx 0.1$ at wavenumber $k_y \rho _{\textrm{s}} = 0.1$ , comparing Gene (blue solid lines) and Key (red dashed lines), where Key eigenfunctions correspond to the numeric-resonant eigenvalues in figure 5. Good agreement is found up to $\beta \approx 1\,\%$ , beyond which some discrepancies become manifest, but substantial agreement persists. At $\beta \gt 3\,\%$ , Gene produces a symmetric eigenfunction while Key obtains two identical modes mirrored across $\theta =0$ (dashed and dotted red lines); the combination of these mirrored modes agrees with the single mode from Gene. Results obtained for $\beta \gt 1.7\,\%$ likely correspond to higher-excitation states of the KBM whose substantial amplitude at high $|\theta |$ is owed to the low background magnetic shear of this configuration and the influence of increasing $\alpha _{\textrm{MHD}}$ (see text for details).

Figure 6 shows eigenfunctions of the perturbed electrostatic potential $\varPhi$ with ballooning angle $\theta$ from Gene and Key for the KBM with $\hat {s} = 0.1$ at normalized binormal wavenumber $k_y \rho _{\textrm{s}} = 0.1$ at different $\beta$ ; these eigenfunctions correspond to the numeric-resonant eigenvalues shown in figure 5. Near marginality, both codes obtain very broad eigenfunctions that gradually decrease in amplitude up to $|\theta | \approx 20 \pi$ , where the mode has an oscillatory structure. With increasing $\beta \approx 1\,\%$ , the amplitude of the mode reduces strongly for $|\theta | \gt 5 \pi$ , and some oscillatory structure remains for $|\theta | \lt 5\pi$ . At $\beta \approx 1.6\,\%$ , the shoulders of the mode increase in amplitude. For higher $\beta$ , peaks become more prominent away from $\theta =0$ , far exceeding the central-peak value. At $\beta \approx 2.2\,\%$ , the central peak has almost fully disappeared while substantial peaks grow at $|\theta | \approx 5\pi - 10\pi$ . This trend continues for higher $\beta \gt 3\,\%$ , where the only remaining peaks occur away from $\theta =0$ , and which move to higher $|\theta |$ with increasing $\beta$ . For $\beta \gt 3\,\%$ , Gene produces a single symmetric eigenfunction, while Key obtains two identical modes which are mirrored across $\theta =0$ (shown in dashed and dotted red lines); the combination of these mirrored modes resembles the single mode from Gene. Such mirror-mode pairs with a shared eigenvalue have been found in previous KBM studies using high-fidelity eigenvalue solver approaches (McKinney et al. 2021; Mulholland et al. 2023, 2025). These mirror-mode pairs may correspond to two finite-ballooning angle modes having radial wavenumbers $+k_x$ and $-k_x$ , where the ballooning angle can be expressed as $\theta _b = -k_x/(\hat {s}k_y)$ . The results obtained here for higher $\beta \gt 1.6\,\%$ are atypical, as KBMs are generally expected to have centrally peaked eigenfunctions with moderate-to-low amplitude at $\theta \neq 0$ ; these modes at high $\beta$ are presumably higher-excitation states (Pueschel et al. 2019) of the KBM. The shape of these eigenfunctions is primarily determined by the field-line curvature $\mathcal{K}_y$ and the influence of magnetic shear manifest in $g^{yy}$ (where the latter is also related to FLR damping). In particular, in this circular tokamak geometry, the bad-curvature regions that drive the eigenmode are located at $\theta = 2\pi m$ , where $m$ corresponds to positive and negative integers; the eigenfunction amplitude peaks at these locations if the stabilizing influence of $g^{yy}$ is not too strong there. The variation of $g^{yy} = 1 + [\hat {s}\theta - \alpha _{\textrm{MHD}} \sin (\theta )]^2$ with increasing $\alpha _{\textrm{MHD}}$ causes the relocation of the eigenfunction peaks to higher $|\theta|$ with increasing $\beta$ . It is noteworthy that KBM eigenfunctions from Key and Gene can also have tearing parity (odd parity in $\varPhi$ , even parity in $A_\parallel )$ , which is less common but sometimes manifest (see McKinney et al. 2021). Despite the unconventional nature of these highly excited eigenmodes, Key manages to obtain good agreement with Gene.

3.3. Analysis and physical interpretation of (st)KBMs

The results in the previous section show that Key accurately captures important features of KBMs in both the stellarator and the tokamak; notably, Key reveals that wave–particle resonances from passing ions subject to the magnetic-curvature drift (i.e., toroidal resonances) must be accounted for in order to detect the presence of stKBMs. In the tokamak, Key qualitatively captures the behavior of conventional KBMs in both high- and low-magnetic-shear scenarios, approximately reproducing the eigenvalue trends and eigenfunctions in each case.

In this section, the results from Key are discussed further with a focus on the distinctions between KBMs in the stellarator and the tokamak; specifically, the presence of stKBMs in W7-X and their absence in the tokamak cases studied here. A physical interpretation is provided for the mechanism of the stKBM and the conditions required for its existence.

Certain conditions are known to be necessary in order for a device to host an stKBM, while other conditions remain speculative and require further scrutiny; assessing the latter is planned for future work. The known conditions pertain both to the plasma itself and to the magnetic geometry in question. Regarding the plasma, a finite ion-temperature gradient ( $\eta _{\textrm{i}} \neq 0$ ) is essential in order for the resonant threshold of destabilization to be below the non-resonant threshold; $\eta _{\textrm{i}} = 0$ leads to both approaches sharing a common $\beta _{\textrm{crit}}^{\textrm{KBM}}$ at the non-resonant threshold, and increasing $\eta _{\textrm{i}}$ leads to an incremental reduction in the resonant threshold relative to the non-resonant counterpart (Cheng 1982b ). It follows that another requirement for stKBM manifestation is the presence of ion-magnetic-drift resonances; a purely non-resonant model is insufficient to capture important details of the KBM near marginality, such that the low- $\beta$ destabilization of the stKBM would be missed. However, it should be noted that the non-resonant approach used here is capable of producing KBM eigenfunctions (not shown here) that are very similar to those obtained in the resonant approaches in figure 2, even in the stKBM region of $\beta \approx 1\,\%{-}3\,\%$ . This indicates that stKBMs arise from a combination of non-resonant and kinetic physics: the non-resonant picture largely captures the eigenfunction shape and approximate eigenvalue trends, while the resonant picture allows access to more accurate eigenvalues at lower drive near marginality, and thus determines the actual stability threshold of the (st)KBM, i.e., whether unstable stKBMs are present, where $\beta _{\textrm{crit}}^{\textrm{stKBM}} \lt \beta _{\textrm{crit}}^{\textrm{KBM}} \lt \beta _{\textrm{crit}}^{\textrm{MHD}}$ .

Regarding geometric conditions needed for the emergence of stKBMs, low background magnetic shear $\hat {s}$ is essential for allowing very broad eigenmodes to exist along the field line – a characteristic feature of the stKBM – but by itself is insufficient to fully explain their presence. This is evidenced by the notable absence of stKBMs in the low-shear tokamak considered here (see also Mulholland et al. 2025, where a similar low-shear tokamak study is presented using a lower pressure gradient), which shares a similar background magnetic shear as W7-X of $|\hat {s}| \approx 0.1$ . Despite this ostensible similarity, it is important to highlight that $\hat {s}$ is an averaged quantity and thus does not provide local information of the field-line geometry; the local differences that persist between these configurations must be factored in when aiming to explain stKBM emergence. These local quantities are those pertaining to FLR stabilization, the magnetic-field-line curvature and local shear, each of which are discussed here.

In the simplified KBM theory outlined here, $\hat {s}$ does not appear explicitly, but plays a role indirectly via magnetic-shear stabilization (a geometric effect) and FLR stabilization (arising from the Bessel functions; geometric and plasma effect), both of which depend on the metric-tensor quantity $g^{yy}(\theta )$ . As an instructive example, consider a circular tokamak equilibrium (with $\beta = \alpha _{\textrm{MHD}} = 0$ , for simplicity) where this quantity is described by $g^{yy}(\theta ) = 1 + \hat {s}^2 \theta ^2$ , such that the amount of suppression from magnetic-shear and FLR effects increases quadratically with ballooning angle $\theta$ ; for high values of $\hat {s}$ , this stabilizing effect quickly becomes substantial such that eigenmodes are restricted to exist near the outboard midplane $\theta \approx 0$ , and no substantial amplitude is expected at higher $|\theta |$ . In contrast, low values of $\hat {s}$ lead to weak suppression at $\theta \neq 0$ , such that modes can have substantial amplitude farther along the field line (higher $|\theta |$ ) and are not restricted to the outboard midplane.

The magnetic-field-line curvature – manifest as $\varOmega _\kappa$ – is a prominent geometric term included in the KBM theory developed here. Bad-curvature regions ( $\mathcal{K}_y \lt 0$ ) act to drive the mode, and so play an important role in determining KBM stability. A notable distinction between the curvature observed in W7-X and the tokamak is its spatial variation along the field line. For example, in a single poloidal turn, W7-X (KJM) has five significant bad-curvature wells, while in the same field-line domain, the tokamak has one. The higher spatial variation of curvature in W7-X ensures the mode (near marginality) is driven consistently along the field line by densely spaced neighboring bad-curvature wells; this feature, combined with the lack of strong FLR suppression and low shear, likely aids the mode in remaining correlated over large portions of the field-line domain. Such non-localized modes may be more prone to early (low- $\beta$ ) destabilization due to their ability to have significant amplitude over substantial portions of the field-line domain without significantly bending the field line in any particular location; a lack of field-line bending from the mode implies a lack of significant stabilizing field-line tension in response. For $\beta _{\textrm{crit}}^{\textrm{stKBM}} \lt \beta \lt \beta _{\textrm{crit}}^{\textrm{KBM}}$ , the stKBM eigenfunction is shown to become increasingly localized with increasing pressure, such that an increase in field-line tension should be experienced by the mode. This balance of increasing drive with increasing restorative stabilization may explain the observed trend of the stKBM’s approximately constant growth rate for $\beta _{\textrm{crit}}^{\textrm{stKBM}} \lt \beta \lt \beta _{\textrm{crit}}^{\textrm{KBM}}$ . Beyond this region, mode narrowing can no longer significantly increase, and the destabilizing drive largely overwhelms the field-line tension; this marks the transition to the fast-growing conventional KBM. These observations support the notion that important characteristics of the stKBM can be explained by employing a simplified non-resonant picture (e.g., the eigenfunction shape and eigenvalue trend), but it should be reiterated that a comprehensive understanding of the stability of this mode requires accounting for kinetic-resonance effects.

The low-shear tokamak also hosts very broad KBM eigenfunctions, and yet does not produce an stKBM; instead, a fast-growing KBM is immediately present at the lowest simulated $\beta$ . This absence of an stKBM is potentially rooted in a combination of geometric characteristics, including – but not limited to – the aforementioned background magnetic shear and magnetic-field-line curvature. Another geometric quantity that may factor in is the local magnetic shear along the field line, defined by $\hat {s}_{\textrm{loc}} = \partial (g^{xy}/g^{xx})/\partial \theta$ . Although this term is not directly present in the simplified KBM model, a related term is: $\partial g^{yy}/\partial \theta$ . Once again, an instructive example to connect these quantities is the circular tokamak equilibrium, for which $g^{yy} = 1 + (g^{xy})^2$ , and $g^{xx}=1$ . In this case, $\hat {s}_{\textrm{loc}} \propto \partial g^{yy}/\partial \theta$ . This fluctuating quantity (along $\theta$ ) has notably higher magnitude (local shear spikes) for a large portion of the field-line domain in W7-X compared with the low-shear tokamak, despite these configurations having comparable background magnetic shear. Such a feature in W7-X – possibly in combination with the field-line tension stabilization described above – may aid in suppressing the stKBM after it has first become unstable at $\beta _{\textrm{crit}}^{\textrm{stKBM}}$ such that its growth rate remains approximately constant before reaching $\beta _{\textrm{crit}}^{\textrm{KBM}}$ , and would at least in part explain the large separation between $\beta _{\textrm{crit}}^{\textrm{stKBM}}$ and $\beta _{\textrm{crit}}^{\textrm{KBM}}$ observed in W7-X. The absence of such substantial local shear spikes in the low-shear tokamak may explain the immediate manifestation of a strongly driven conventional KBM at low $\beta$ rather than an stKBM.

The effects of these geometric characteristics on (st)KBM behavior in stellarator geometry will be studied in greater detail in future work; specifically, the influence of local/background magnetic shear and Shafranov shift.

4. Conclusions

One focus of MCF research is to improve on the energy-confinement time which is largely limited by energy losses from turbulent transport. As fusion devices strive to achieve high performance (high pressure), regimes with substantial $\beta$ must be studied and understood. The pressure-gradient-driven instabilities relevant in such scenarios are then of particular interest to address, such that they can be controlled or mitigated. This work aims to expand the current understanding of one such instability – the KBM – in general toroidal geometry (i.e., in tokamaks and stellarators) that includes kinetic effects such as magnetic-drift resonances and FLR stabilization.

The KBM theory developed here (see § 2) relies on several simplifications: certain physics effects are neglected, such as particle trapping, parallel-magnetic-field fluctuations $\delta B_\parallel$ , Landau resonances and collisions; the absence of $\delta B_\parallel$ prompts the use of a magnetic-curvature-drift approximation of $\omega _{\boldsymbol{\nabla }B} \to \omega _\kappa$ ; and only eigenmodes centered at $k_x = 0$ are considered. The resulting KBM eigenmode equation (2.32) depends nonlinearly on the KBM eigenvalue, and thus is reformulated as a nonlinear eigenvalue problem (see § 2.2); the KBM eigenvalues are first obtained, followed by treating the system as a special case of a linear eigenvalue problem, which readily yields the KBM eigenvectors. Despite the aforementioned simplifications, this KBM model performs well in rapidly predicting the correct eigenvalue trends and KBM eigenfunctions in a variety of scenarios (W7-X stellarator; high- and low-shear circular tokamak), as compared with Gene simulations (see § 3).

The numerical implementation of this theory has resulted in the Key code, capable of calculating self-consistent eigenvalues and eigenvectors for (st)KBMs along the field line in any flux-tube geometry, where the only required inputs are flux-tube geometry and plasma parameters. Applying the Key code to both stellarator and tokamak geometry elucidates the required conditions for the manifestation of stKBMs (discussed in § 3.3): wave–particle resonances involving passing ions subject to the magnetic-curvature drift must be accounted for, which can heavily reduce the destabilization threshold of the KBM relative to the iMHD threshold for $\eta _{\textrm{i}} \neq 0$ ; low background magnetic shear $\hat {s}$ is necessary to allow for the existence of broad eigenmodes that remain correlated for large extents along the field line; other factors that may play an important role in determining stKBM emergence are the magnetic-field-line curvature and the local magnetic shear along the field-line domain (see discussion in § 3.3). A closer examination of these geometric influences on (st)KBM behavior is reserved for future work.

Despite its simplicity, this model shows great promise in a variety of aspects: it can be used as a predictive tool for determining stKBM manifestation in general toroidal geometry, which may aid in the development of reduced models for predicting turbulent transport at high $\beta$ ; the simplified theory in Key provides a basic framework to help elucidate the essential physics of (st)KBMs, enhancing the current understanding of these instabilities and their relationship with geometry; finally, the speed, affordability and fidelity of Key allow for its utilization in turbulence-optimization routines for both shaped tokamaks and stellarators aiming to achieve and sustain high- $\beta$ plasma. Thus, the results of this work may help inform the design of high- $\beta$ -turbulence-optimized fusion reactors.