The kinetic ion-temperature-gradient-driven instability and its localisation

Abstract

We construct a description of ion-temperature-gradient (ITG)-driven localised linear modes which retains both wave–particle and magnetic drift resonant effects while capturing the field-line dependence of the electrostatic potential. We exploit the smallness of the magnetic drift and the strong localisation of the mode to resolve the problem with a polynomial–Gaussian expansion in the field-following coordinate. A simple semianalytical formula for the spectrum of the mode is shown to capture long wavelength Landau damping, ion-scale Larmor radius stabilisation, weakening of Larmor radius effects at short wavelengths and magnetic-drift resonant stabilisation. These elements lead to linear spectra with multiple maxima as observed in gyrokinetic simulations in stellarators. Connections to the transition to extended eigenfunctions and those localised by less unfavourable curvature regions (hopping solutions) are also made. The model provides a clear qualitative framework with which to interpret numerically simulated ITG modes’ linear spectra with realistic geometries, despite its limitations for exact quantitative predictions.

1 Introduction

The ion temperature gradient (ITG) is one of the primary modes driving turbulence in magnetic confinement fusion devices (Rudakov & Sagdeev 1961; Coppi, Rosenbluth & Sagdeev 1967; Terry, Anderson & Horton 1982). As such, much literature exists devoted to understanding this type of instability. The most basic understanding of this mode can be gained by studying the linear response of the system, as described by the linearised gyrokinetic (GK) equation (Connor, Hastie & Taylor 1980; Romanelli 1989) used to evaluate quasineutrality. Extending our understanding of the mode and its response to the geometry is particularly important in the context of stellarator physics.

Linear studies of the ITG-driven mode, either analytical or numerical, are numerous in the literature, but they are all fundamentally traceable to the review of Kadomtsev & Pogutse (1970b). Analytical progress often requires some simplifying assumptions. These occur at two levels. First, by considering separately the various destabilising mechanisms in the problem. Second, by reducing the complexity of the magnetic field along field lines (e.g. curvature, flux compression or the magnetic field magnitude), often described as constant. Our knowledge of the ITG builds on consideration of special cases that respond to different approximations of the first kind (Zocco et al. 2018). When the bad curvature, together with the temperature gradient, drives the mode to being unstable, we have a toroidal ITG (Terry et al. 1982). When destabilisation does not require curvature, but is enabled by the propagation of density along the field lines, with a specific relative phase with respect to temperature (Cowley, Kulsrud & Sudan 1991), we have a slab ITG (Rudakov & Sagdeev 1961; Kadomtsev & Pogutse 1970b). While any simplified perspective helps shed light on the physics of the ITG mode, the selective treatment of the physics can hinder their scope. A pressing case of this is the overlooking of the mode localisation along the field line in the toroidal branch (Terry et al. 1982), where the presence of bad curvature is paramount. Magnetic fields generally exhibit alternating regions of good and bad curvature every connection length.

A simple theoretical description of the inhomogeneity along the magnetic field line, and the consequent behaviour of the ITG mode, is approachable through a fluid description (Horton, Choi & Tang 1981; Hahm & Tang 1988; Romanelli 1989; Plunk et al. 2014). Assuming particle parallel streaming to be small (slow ion transit time), a strong drive and small curvature drift, one can describe the ITG through a second-order ordinary differential equation (ODE) along the field line. Such a description incorporates key physics ingredients to the problem, and crucially couples the behaviour of the instability to its structure along the field line. However, while localisation seems so important for a mathematical characterisation of the fluid branches (Wesson & Campbell 2011; Zocco et al. 2016) and interpretation of numerical results, its analytical treatment becomes increasingly intricate when Landau damping and the magnetic drift resonance are included. Without a consistent inclusion of these kinetic elements, the fluid description suffers from fundamental shortcomings in describing the behaviour of microinstabilities at small scales (either large perpendicular wavenumber or small parallel scale). This limitation can lead to wrong expectations in the mode response to changes in geometry. But the latter being of crucial importance, especially in the context of current activities such as stellarator optimisation, we need to integrate all elements together. While much progress has been made in this direction for the development of the theory of geodesic acoustic modes (Zonca, Chen & Santoro 1996; Sugama & Watanabe 2006a,b; Qiu, Chen & Zonca 2018), a manageable kinetic theory that retains resonant effects and includes geometric localisation is desirable. We address this problem in this work.

In this article we introduce a formal approach to the GK equation to describe localised ITG modes and their structure in a simple geometry with good and bad curvature alongside kinetic effects, both resonant (due to particles parallel streaming and magnetic drift) and non-resonant (due to Larmor radius effects). In § 2 we start from the linearised GK equation and transform it into a linear system of equations, introducing the key ordering in the curvature drift and localisation of the mode. In the following section, § 3, we focus on constructing the dispersion equation for a model simple Gaussian-shaped mode, and detail its derivation and numerical solution. Section 4 then studies the physical elements of the model constructed, including the roles of Landau damping, finite Larmor radius (FLR) stabilisation and the curvature drift resonance. Estimates for important features such as critical temperature gradients are also given. Section 5 then discusses the behaviour of other modes other than those with the simplest Gaussian structure, to end with some numerical GK simulations. Although the model is too simple to be quantitatively correct, we show it serves as a blueprint to interpret the results from the simulations. We use this to make connections to the presence of delocalised (slab) and hopping modes.

2 Localised description of the GK equation

2.1 Rewriting the GK equation

The starting point of our model construction is the linear GK equation for the ions, which we write as follows:

(2.1)\begin{equation} {\rm i}v_\parallel\partial_\ell g + (\omega - \tilde{\omega}_d)g = \frac{q_i}{T_i}F_{0i} {\rm J}_0(\omega-\tilde{\omega}_\star)\phi. \end{equation}

The equation is to be understood as a first-order ODE in $\ell$ (Taylor 1976), the distance along a field-line, for the non-adiabatic part of the distribution function, denoted by $g$, which is a perturbation with respect to the background ion Maxwellian $F_{0i}$. To write the equation in this one-dimensional form, the ballooning transform has been considered (Taylor 1976; Connor, Hastie & Taylor 1978; Tang, Connor & Hastie 1980; Antonsen & Lane 1980) so that $\boldsymbol {k}_\perp =k_\alpha \boldsymbol {\nabla }\alpha +k_\psi \boldsymbol {\nabla }\psi$. Here $\alpha$ and $\psi$ are defined to be the straight field line Clebsch variables (D'haeseleer et al. 2012), such that $\boldsymbol {B}=\boldsymbol {\nabla }\psi \times \boldsymbol {\nabla }\alpha$ and the toroidal magnetic flux is $2\pi \psi$. The response of the system is coupled to the electrostatic potential $\phi$, and is modulated by the geometry of the magnetic field. Elements of the geometry are present in $\tilde {\omega }_d$ and ${\rm J}_0$. The former represents the ion magnetic particle drifts, which in the small $\beta$ limit may be written as $\tilde {\omega }_d=\omega _d(x_\parallel ^2+x_\perp ^2/2)$, where $\omega _d=\boldsymbol {v}_D\boldsymbol {\cdot }\boldsymbol {k}_\perp =v_{Ti}\rho _i\kappa \times \boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {k}_\perp /B$, where $x_\parallel =v_\parallel /v_{Ti}$ and $x_\perp =v_\perp /v_{Ti}$ are velocities normalised to the ion thermal speed $v_{Ti}=\sqrt {2T_i/m_i}$, and $\rho _i=m_i v_T/q_iB$ is the ion Larmor radius, where $T_i$, $q_i$ and $m_i$ are the ion temperature, charge and mass. For simplicity, we will focus on the simpler limit $k_\psi =0$. The second element of geometry is included in the FLR term ${\rm J}_0={\rm J}_0(x_\perp \sqrt {2b})$, where ${\rm J}_0$ is the Bessel function of the first kind, and $b=(k_\perp \rho _i)^2/2$. The drive of the instability is included in the diamagnetic drift, represented here by $\tilde {\omega }_\star =\omega _\star [1+\eta (x^2-3/2)]$, where $\eta =\mathrm {d}\ln T_i/\mathrm {d}\ln n_i$ is the ratio of ion temperature to ion density gradients, and $\omega _\star =(k_\alpha T_i/q_i) (\mathrm {d}\ln n_i/\mathrm {d}\psi )$.

Because the GK equation is a first-order ODE, a solution for $g$ can be written in its most general form using an integrating factor, as originally presented by Connor et al. (1980). However, the resulting integral expressions, in their generality, do not always manifest clearly the role played by the different physical elements in the problem. After imposing quasineutrality, the relation between the mode structure and the degree of instability is not clear, as an involved integral eigenvalue problem ensues (Romanelli 1989). To circumvent this complexity, we present here an approach that makes the resonant kinetic problem as close as possible to a second-order ODE, much in the way that it occurs in the fluid limit.¹

To that end, we first invoke symmetry arguments to simplify the construction of the solution as much as possible. Given the explicit involvement of $v_\parallel$ in the GK equation, it is convenient to separate $g$ into even and odd parts in $v_\parallel$, namely

(2.2a)\begin{gather} g^e(\ell,\mathbf{v}) = \frac{1}{2}\left[g(v_\parallel)+g({-}v_\parallel)\right], \end{gather}

(2.2b)\begin{gather}g^o(\ell,\mathbf{v}) = \frac{1}{2}\left[g(v_\parallel)-g({-}v_\parallel)\right]. \end{gather}

With these well-defined parity functions, the original GK equation, which had mixed parity, can be separated into two coupled first-order differential equations,

(2.3a)\begin{gather} {\rm i}v_\parallel\partial_\ell g^o+(\omega-\tilde{\omega}_d)g^e =\frac{q_i}{T_i}F_{0i}{\rm J}_0(\omega-\tilde{\omega}_\star)\phi, \end{gather}

(2.3b)\begin{gather}{\rm i}v_\parallel\partial_\ell g^e+(\omega-\tilde{\omega}_d)g^o = 0. \end{gather}

The latter is used to eliminate $g^o(\ell,\boldsymbol {v})$ from the former, to give an equation that only involves the even part of $g$ in $v_\parallel$,²

(2.4)\begin{equation} v_\parallel\partial_\ell\left[\frac{v_\parallel}{\omega-\tilde{\omega}_d}\partial_\ell g^e\right]+(\omega-\tilde{\omega}_d)g^e=\frac{q_i}{T_i}F_{0i} {\rm J}_0(\omega-\tilde{\omega}_\star)\phi. \end{equation}

We must not forget that the linearised GK equation, (2.1), does not come on its own. First, $g^e$ must satisfy vanishing boundary conditions at $\ell \rightarrow \pm \infty$ for a physically reasonable ballooning solution of the equation (Connor et al. 1978). Second, it must be complemented by the quasineutrality condition, the condition preventing charge separation from building up in the system. This imposes an additional relation between the velocity-space function $g^e$ and the real space function $\phi$, which is necessary to complete the $\omega$-eigenvalue equation. As the velocity-space average of the odd $g^o$ vanishes due to parity, the quasineutrality condition reduces to

(2.5)\begin{equation} \int {\rm J}_0 g^e \,\mathrm{d}^3\boldsymbol{v} = \bar{n}(1+\tau)\phi, \end{equation}

where $\tau =T_i/ZT_e$ and $\bar {n}$ is the equilibrium density. We are taking the electron response to be adiabatic here. Equations (2.4) and (2.5), describing our ion response, will be our starting point.

2.2 Approximations: localisation and weak curvature

So far, the problem defined by (2.4) and (2.5) is as general as the standard form of the linearised GK equation. To proceed, we introduce a number of simplifying assumptions that will make it tractable while retaining the main physical elements of the problem.

2.2.1 Simplified geometry

The first element of simplification in the problem is the geometry along the field line. We restrict all the inhomogeneity in the field to the $\ell$ dependence of the curvature drift $\omega _d(\ell )$. Doing so lets us capture the key aspect of having good and bad curvature regions along the field line. On physical grounds (Terry et al. 1982), we expect the unstable toroidal mode to have a tendency to localise around bad curvature regions, which are energetically favourable, while being repelled from good curvature ones.

To introduce a sense of this feature, we describe a single region of bad curvature, taken to be symmetric in $\ell$ and to have a finite extent, $\varLambda$, beyond which the curvature changes sign. In the usual jargon, this length scale may be deemed the connection length, and the region within can be thought as a bad curvature well. We model the well as a two-parameter simple symmetric quadratic function in $\ell$, where the parameters represent its depth and width. The magnitude of the bad curvature at the origin is defined to be $-\bar {\omega }_d<0$. The curvature remains bad, i.e. negative, in the region $|\ell |<\varLambda$. We then write $\omega _d$ as

(2.6)\begin{equation} \omega_d=\bar{\omega}_d\left[\left(\frac{\ell}{\varLambda}\right)^2-1\right]. \end{equation}

Because $\varLambda$ is the only existing length scale along the field line, as both $B$ and $k_\perp \rho _i$ are assumed to be constant, we shall use it to normalise lengths and write $\bar {\ell }=\ell /\varLambda$.

Before moving on, we should reflect on the consequences of our simplified geometry. Choosing the magnetic field magnitude to be flat along the field line eliminates any contribution from trapped ions, which do not exist in our model. The physics associated with the variation of $k_\perp \rho _i$ are also lost, and with it the possible modulation of FLR effects. In particular, this approximation erases the effects of global magnetic shear, which would have appeared as a secular term in $k_\perp$. The global shear can become important especially in localising significantly extended, slab-like modes, and thus forcing the correct behaviour of $g$ at $\ell \rightarrow \pm \infty$. Instead, in our problem, in the limit of large $|\bar {\ell }|$ we have a strongly positive good curvature. Physically, we expect this to play a stabilising role for the typical toroidal ITG that precludes modes from becoming completely delocalised. In that sense, our model of $\omega _d(\ell )$ mimics in part the action of global shear.³ Unavoidably, this simplification couples the local and global behaviour of the system. All the simplifications considered will render less accurate a quantitative comparison of the analytical results with real fields, but will serve as an important qualitative and semiquantitative tool. We prove numerically that this particular approximation gives excellent results for the real geometry of the quasisymmetric (Boozer 1983; Nührenberg & Zille 1988; Rodríguez, Helander & Bhattacharjee 2020) HSX (helically symmetric experiment) stellarator (Anderson et al. 1995).

Reducing all the field inhomogeneity to $\omega _d$ is a significant formal simplification, making the differential operator $\partial _\ell$ commute with everything in (2.4) except $g$ (we shall drop the superscript $e$ describing parity for simplicity), $\phi$ and $\omega _d$. The assumption of a constant $B$ is of particular importance here. The partial derivative respect to $\ell$ is meant to be taken keeping the velocity-space variables $\mu =m_iv_\perp ^2/2B$ and $\mathcal {E}=m_iv^2/2$ constant. Only when $B$ is constant is $v_\parallel$ independent of real space.

With this observation, (2.4) becomes, upon commutation of the differential operator

(2.7)\begin{align} \left(1-\tilde{\omega}_d/\omega\right)^{{-}1}\frac{v_\parallel^2}{\omega^2}\partial_\ell^2 g+\left(1-\frac{\tilde{\omega}_d}{\omega}\right) g+\frac{v_\parallel^2}{\omega^2}\frac{\partial_\ell\left(\tilde{\omega}_d/\omega\right)}{\left(1-\tilde{\omega}_d/\omega\right)^2}\partial_\ell g=\frac{q_i}{T_i}F_{0i}{\rm J}_0\left(1-\frac{\tilde{\omega}_\star}{\omega}\right)\phi. \end{align}

We have a second-order ODE in which all the explicit spatial dependence can be expressed in terms of polynomials in $\ell$.

2.2.2 Weak curvature drift, strong drive and strong localisation

It is clear from the above equation that the drift frequency in this scenario plays a central role in prescribing the behaviour of the instability. The mode will adapt to the geometry described by $\omega _d(\bar {\ell })$, and thus to make the treatment more manageable, it is natural to order the drift. We introduce the ordering parameter $\delta \sim \bar {\omega }_d/\omega \ll 1$. Such a restriction is not completely artificial, and it is particularly appropriate in scenarios in which the turbulence is strongly driven, namely, in cases where the $\omega _\star$ drive (either in its density or temperature gradient form) is much larger than $\omega _d$. We shall be focusing on this strongly driven scenario.

Note, however, that ordering $\bar {\omega }_d/\omega$ is not enough to simplify (2.7). Although $\bar {\omega }_d/\omega$ may be small, this is only its value at the bottom of the bad curvature well. Thus, there always exists a sufficiently large value of $\bar {\ell }$ such that $\omega _d/\omega \gtrsim 1$. This appears hopeless for an approximated approach to (2.7). However, we must not consider $\omega _d$ on its own when doing so, but rather alongside $g$ and $\phi$. If the distribution function and the potential are finite only over some finite length scale, then the effective value of $\bar {\ell }$ should reflect that finite extent. Considering a Gaussian-like envelope of the form $\sim {\rm e}^{-\lambda \bar {\ell }^2/2}$, so that $g$ and $\phi$ have length scales $\sim 1/\sqrt {\mathrm {Re}\{\lambda \}}$, where $\mathrm {Re}\{\lambda \}$ denotes the real part of $\lambda$ which is generally complex, we limit this problem with a new ordering parameter,

(2.8)\begin{equation} \epsilon=\frac{1}{\mathrm{Re}\{\lambda\}}\frac{\bar{\omega}_d}{\omega}\ll1. \end{equation}

The parameter $\epsilon$ restricts the following construction to modes that are sufficiently localised. Thus, modes may have some degree of delocalisation, but limited by the ordering parameter $\delta$. We will then consider two simultaneous ordering assumptions, namely $\delta,\epsilon \ll 1$. The precise implications of these orderings and how we proceed with them is detailed in what follows.

First, let us employ the newly introduced approximations to simplify (2.7). Expanding the equation to, and keeping, order $O(\epsilon,\delta )$,

(2.9)\begin{align} 2\left(\frac{\omega_tx_\parallel}{\omega}\right)^2\partial_{\bar{\ell}}^2 g+\left(1-\frac{\tilde{\omega}_d}{\omega}(\bar{\ell}^2-1)\right)g+4\tilde{\omega}_d\left(\frac{\omega_tx_\parallel}{\omega}\right)^2\bar{\ell}\partial_{\bar{\ell}}g=\frac{q_i}{T_i}F_{0i}{\rm J}_0\left(1-\frac{\tilde{\omega}_\star}{\omega}\right)\phi, \end{align}

where allowing some looseness in the notation, $\tilde {\omega }_d=\bar {\omega }_d(x_\parallel ^2+x_\perp ^2/2)$ here, and $\omega _t$ is the transit frequency defined as

(2.10)\begin{equation} \omega_t=\frac{v_{Ti}}{\varLambda\sqrt{2}}. \end{equation}

The consistent ordering of $\omega _t/\omega$ will be explicitly discussed later. For now, we take it to be order 1, but consider $\bar {\omega }_d/\omega$ corrections to the first term in (2.9) small.

Localised solutions to a second-order ODE like (2.9) may be approximated by considering a representation of $g$ and $\phi$ in a Taylor–Gauss basis,

(2.11)\begin{equation} g=\sum_{n=0}^\infty \frac{g_n\bar{\ell}^n}{N(n)} {\rm e}^{-\lambda\bar{\ell}^2/2}, \end{equation}

where $N(n)=\sqrt {n!/\mathrm {Re}\{\lambda \}^n}$ is a normalisation factor. The Taylor part of the basis (i.e. the expansion in powers of $\bar {\ell }$) naturally describes the mode near the bottom of the well, while the Gaussian part provides an overall envelope that localises the mode. The latter requires $\mathrm {Re}\{\lambda \}>0$, although it is consistent with having a non-zero imaginary part.⁴ The normalisation factor includes powers of $\mathrm {Re}\{\lambda \}$ to account for the scale associated with the monomial $\bar {\ell }^n$. Note that the higher the mode considered, the increasingly hollower the shape of the mode is, providing a characteristic length scale $\sqrt {n}/\mathrm {Re}\{\lambda \}^{n/2}$.

Once we have this basis representation, the ODE (2.9) becomes an infinite set of coupled algebraic linear equations on $\{\phi _n\}$ and $\{g_n\}$. The benefit of the particular basis used is the simplicity of the form that the operations in (2.9) take. The product by $\bar {\ell }^2$ in the second term of (2.9) simply upshifts the mode number by two, making the regularising role of $\epsilon$ manifest. It restricts the coupling of the different $\{g_n\}$ and $\{\phi _n\}$ modes, and thus is critical in achieving a useful truncation of the linear system of equations. Differentiation plays a similar, albeit more involved, coupling role (see Appendix A).

2.3 The Taylor–Gauss form of the GK equation

We are now in a position to resolve (2.9) in the new basis of (2.11). This can be achieved by substituting the expansion in (2.11), and applying the coupling rules in (A1) given in the Appendix A. Upon substitution, the resulting equation takes the form

(2.12a)\begin{equation} \sum_{n=0}^\infty E_{n}\sqrt{\frac{\mathrm{Re}\{\lambda\}^n}{n!}}\bar{\ell}^n{\rm e}^{-\lambda\bar{\ell}^2/2}=0, \end{equation}

where, the general expression for the nth mode equation is

(2.12b)\begin{align} E_n& =2\mathrm{Re}\{\lambda\}\left(\frac{\omega_t x_\parallel}{\omega}\right)^2\sqrt{(n+2)(n+1)}g_{n+2} \nonumber\\ & \quad +\left[1+\frac{\tilde{\omega}_d}{\omega} -2\lambda(2n+1)\left(\frac{\omega_tx_\parallel}{\omega}\right)^2+4n\frac{\tilde{\omega}_d}{\omega}\left(\frac{\omega_t x_\parallel}{\omega}\right)^2\right]g_n\nonumber\\ & \quad +\frac{2}{\mathrm{Re}\{\lambda\}}\left[\lambda^2\left(\frac{\omega_t x_\parallel}{\omega}\right)^2-\frac{\tilde{\omega}_d}{2\omega}-2\lambda\frac{\tilde{\omega}_d}{\omega}\left(\frac{\omega_t x_\parallel}{\omega}\right)^2\right]\sqrt{n(n-1)}g_{n-2}\nonumber\\ & \quad -\frac{q_i}{T_i}F_{0i}{\rm J}_0\left(1-\frac{\tilde{\omega}_\star}{\omega}\right)\phi_n. \end{align}

For (2.9) to be true, (2.12a) must hold for all $\bar {\ell }$, meaning that each of the equations $E_{n}$ must vanish separately. Hence, the original ODE becomes an infinite-dimensional system of linear equations, $\{E_n=0\}$, as promised.

An inspection of the structure of this equation shows two important features. First, even and odd modes in $\bar {\ell }$ are completely decoupled. This is a direct consequence of the original form of the equation having well-defined parity. Unless otherwise stated, we shall consider the even set, of which its most basic form is a Gaussian mode. The results follow similarly for the odd set. Second, within each of these subspaces, the system has a tridiagonal structure. That is, the equation couples every mode to the immediately adjacent ones.

3 Gaussian mode model of kinetic ITG

The resolution of the original equation into the Taylor–Gauss basis leaves us with an infinite set of algebraic equations to solve in phase space (i.e. in $\boldsymbol {v}$ and $\boldsymbol {\ell }$). It is our goal now to truncate this hierarchy down to obtain a useful model of the problem.

3.1 Constructing the dispersion function

Let us start first by understanding what the structure of the set of equations is when we truncate the system at order $N$. That is, when we set all $g_n,~\phi _n=0$ for $n>N$. In that case the finite set of equations we are left with is

(3.1)\begin{equation} \{E_0(0,2),\ldots,E_n(n-2,n,n+2),\ldots,E_{N}(N-2,N),E_{N+2}(N)\}, \end{equation}

where each element must vanish and the numbers in parenthesis denote the modes (both in $g$ and $\phi$) involved in the equations. This system of equations may be rewritten by rearranging the first $N/2+1$ equations to solve explicitly for $\{g_n:n\leq N\}$ in terms of $\{\phi _n:n\leq N\}$ by appropriate linear combinations. This is generally possible, and an explicit construction to order $O(\epsilon ^2)$ is provided in Appendix B following (2.12b). We then write

(3.2)\begin{equation} \left. \begin{aligned} g_n(\boldsymbol{v}) & =\sum_{m=0}^{N}\mathbb{D}^{(N)}_{nm}(\boldsymbol{v})\phi_m \quad (n\leq N),\\ g_N(\boldsymbol{v}) & =\bar{\mathbb{D}}_{N}(\boldsymbol{v})\phi_N, \end{aligned} \right\} \end{equation}

where $\mathbb {D}_{nm}$ and $\bar {\mathbb {D}}_{N}$ are known functions of velocity space. The last isolated equation on $g_N$, coming from $E_{N+2}$, is a result of the truncation at a finite $N$, and leads to $g_N$ satisfying two equations simultaneously. The lack of degrees of freedom to salvage this overdetermination in velocity space is indicative of our failure to satisfy the original equation exactly with a finite number of modes. After all, we are attempting an approximate solution. We must therefore assess what is the error made by relaxing this last constraint. To do so, consider $M$ to be the dominant mode in the system, which is taken to be $O(1)$. As shown in Appendix B, since the mode coupling terms are order $\epsilon$, we expect the magnitude of the error made by dropping that last equation to be $\sim \epsilon ^{(N-M)/2}$. The error made is thus small provided we restrict ourselves to $M\ll N\ll N_\epsilon =1/\epsilon$. The upper bound $N_\epsilon$ is set to preserve the presumed smallness of $\epsilon$ as it is amplified by the mode number of the solution to our system of equations. If $N$ is chosen to be sufficiently large, the matrix elements $\mathbb {D}_{nm}^{(N)}$ should then become largely independent of $N$, and we may drop the $(N)$ superscript. Having a model that is approximately independent of the truncation is appropriate, and we shall see that a special choice of $\lambda$ enacts this truncation most snugly.

Our GK problem is now reduced to the main linear system of equations in (3.2). To complete the problem, we apply quasineutrality, (2.5). Taking the appropriate velocity moment of $\mathbb {D}_{nm}$, and defining

(3.3)\begin{equation} \mathcal{D}_{nm}=\frac{1}{\bar{n}}\frac{T_i}{q_i}\int {\rm J}_0 \mathbb{D}_{nm}(\boldsymbol{v})\,\mathrm{d}^3\boldsymbol{v}, \end{equation}

the resulting eigenvalue problem becomes

(3.4a)\begin{gather} \sum_{j=0}^N \mathbb{M}_{ij}\phi_j=0, \end{gather}

(3.4b)\begin{gather}\mathbb{M}_{ij}=(1+\tau)\delta_{ij}-\mathcal{D}_{ij}. \end{gather}

Once we solve this system of equations, we may then construct $g_n$ for $n< N$, acknowledging the order $\epsilon ^{N-M}$ error made as described above.

We now illustrate our way forward for the $M=0$ mode (the approach for a general $M$ is given in Appendix B). This is to take the Gaussian centred about the bad curvature to be our dominant mode; formally, $\phi _0\sim O(1)$. Following the principle of staying as distant as possible from the truncation point of the system, we may ask when the solution to the problem is consistent with $\phi _n=0$ for all $n>0$. We call this a pure mode, in so much as it is consistent with a ‘complete’ decoupling from higher mode numbers.

To order $\epsilon$, the system describing this ‘pure’ $n=0$ mode reduces to the following two equations:

(3.5a)\begin{gather} \mathcal{D}_{20}=0, \end{gather}

(3.5b)\begin{gather}1+\tau-\mathcal{D}_{00}=0, \end{gather}

where the expressions for $\mathcal {D}_{00}$ and $\mathcal {D}_{20}$ to order $\epsilon$ are

(3.6a)\begin{gather} \mathcal{D}_{20}=\frac{1}{\bar{n}}\dfrac{2\sqrt{2}}{\mathrm{Re}\{\lambda\}}\int\,\mathrm{d}^3\boldsymbol{v}{\rm J}_0^2 F_{0i}\left(1-\dfrac{\tilde{\omega}_\star}{\omega}\right)\left[\lambda^2\left(\dfrac{\omega_t x_\parallel}{\omega}\right)^2-\dfrac{\tilde{\omega}_d}{2\omega}\right], \end{gather}

(3.6b)\begin{gather}\mathcal{D}_{00}=\dfrac{1}{\bar{n}}\int\,\mathrm{d}^3\boldsymbol{v} {\rm J}_0^2F_{0i}\left(1-\dfrac{\tilde{\omega}_\star}{\omega}\right)\dfrac{1}{1+\dfrac{\tilde{\omega}_d}{\omega}-2\lambda\left(\dfrac{\omega_tx_\parallel}{\omega}\right)^2}. \end{gather}

Both of these expressions constitute the governing dispersion relation for our mode. It might appear that these make the problem overconstrained;however, we must bear in mind that $\omega$ is not the only unknown here. The localisation $\lambda$ is as well, and it must be chosen alongside the frequency of the mode. This is analogous to what happens in the fluid limit (Hahm & Tang (1988); R.J. Hastie, private communication, 2013; Connor, Hastie & Zocco (2013); Zocco et al. (2016)). In fact, a closed form for $\lambda$ may be obtained from satisfying (3.6a). For simplicity, considering the small $b$ limit and the leading-order form of the expression in $\delta,~\epsilon$,

(3.7)\begin{gather} \mathcal{D}_{20}\approx \frac{\sqrt{2}}{\mathrm{Re}\{\lambda\}}\left[1-\frac{\omega_{{\star}}}{\omega}(1+\eta)\right]\left[\lambda^2-\frac{\bar{\omega}_d\omega}{\omega_t^2}\right]\approx0, \nonumber\\ \therefore \lambda=\sqrt{\frac{\omega\bar{\omega}_d}{\omega_t^2}}. \end{gather}

For physically meaningful modes, and to restrict ourselves to $\mathrm {Re}\{\lambda \}>0$, we define the square root in (3.7) with its branch cut along the negative real axis in the complex plane (its conventional definition). The mode envelope $\lambda$ indicates a balance between the parallel streaming, $\omega _t^2$, and the curvature drift $\omega _d$, as is well known to be behind the localisation mechanism in the fluid limit of the ITG (Hahm & Tang 1988; Connor et al. 2013; Zocco et al. 2016). Upon approaching marginal stability, the mode will tend to become increasingly delocalised, with an oscillating mode structure when it corotates with ions.

The basic scaling of $\lambda$ together with the ordering of $\epsilon$ implies that $\xi =\omega _t/\omega$ follows $\delta ^{1/2}\xi \leq \epsilon \ll 1$. It is convenient then to take for the ordering arguments $\xi \delta ^{1/2}\sim \epsilon$. This is consistent with the assumptions made to reach this point. Although such consistency is reassuring, the precise form of $\lambda$ in (3.7) is a direct consequence of the assumption of a ‘pure’ mode. This purity assumption is, however, not whimsical, but particularly representative. Not only is it consistent with the ordering, but it also brings an elegant and efficient closure to our system of equations (see a discussion on this in Appendix B), as well as granting the correct asymptotic fluid limit of the dispersion equation, as we shall later discuss. The choice of $\lambda$ may thus be interpreted as a ‘fluid’ choice, with all the approximations that this involves. However, as one may check upon lifting the purity assumption, this remains a simplifying choice. In fact, such a relaxation eliminates in our case the $\mathcal {D}_{20}=0$ constraint, leaving $\lambda$ undetermined.⁵ Future work may be devoted to improving the model by determining $\lambda$ in some other way, perhaps by treating it as a free parameter to optimise in order to maximise the growth rate of the ITG. However, in this first work we willingly sacrifice a precise quantitative prediction for simplicity. We will illustrate this expected quantitative mismatch comparing the model with some example GK simulations. Nevertheless, it will be seen that the model remains a highly useful tool to interpret complicated linear turbulence spectra.

With this simple form for $\lambda$ as a function of $\omega$, we then interpret (3.6b) as a dispersion relation for $\omega$,

(3.8)\begin{equation} \mathcal{D}=1+\tau+\frac{\zeta}{\bar{n}}\int {\rm J}_0^2\left(1-\dfrac{\tilde{\omega}_\star}{\omega}\right)\dfrac{F_{0i}}{x_\parallel^2-\zeta\left(1+\dfrac{\bar{\omega}_d}{2\omega}x_\perp^2\right)}\,\mathrm{d}^3\boldsymbol{v}, \end{equation}

where

(3.9a)\begin{equation} \zeta=\frac{1}{2}\left[\lambda\left(\frac{\omega_t}{\omega}\right)^2-\frac{\bar{\omega}_d}{2\omega}\right]^{{-}1}, \end{equation}

or using the form of $\lambda$ in (3.7),

(3.9b)\begin{equation} \zeta=\frac{\omega^2}{2\lambda(1-\lambda/2)\omega_t^2}. \end{equation}

The parameter $\zeta$ can be interpreted as a measure of the kinetic effects, being important for $|\zeta |\sim 1.$ Equation (3.9a) includes resonant kinetics that can come in either through a Landau-type resonance involving the structure along the field line (the $\lambda$ piece), or the magnetic drift resonance. The relative importance of these terms will depend on the scale of the mode structure along the field line, i.e. the magnitude of $|\lambda |$.

3.2 Simplifying the dispersion function

To evaluate the dispersion relation in (3.8) we must perform the necessary velocity-space integrals, which includes resolving a resonant denominator. We shall deal with it by first performing the integral over $x_\parallel$ (and ignoring the resonance in $x_\perp$, as detailed in Appendix C). To that end, we need to resolve integrals of the following form:

(3.10)\begin{equation} I_{{\parallel},\beta}(\zeta)=\frac{1}{\sqrt{\pi}}\int_{-\infty}^\infty x_\parallel^{2\beta}\frac{{\rm e}^{{-}x_\parallel^2}}{x_\parallel^2-\zeta}\,\mathrm{d} x_\parallel. \end{equation}

Using the difference of two squares for the denominator, we rewrite the integral

(3.11)\begin{equation} I_{{\parallel},\beta}(\zeta)=\frac{1}{\sqrt[*]{\zeta}}\frac{1}{\sqrt{\pi}}\int_{-\infty}^\infty\frac{x_\parallel^{2\beta} {\rm e}^{{-}x_\parallel^2}}{x_\parallel{-}\sqrt[*]{\zeta}}\,\mathrm{d} x_\parallel, \end{equation}

which will soon be related to the plasma dispersion function (Fried & Conte 2015). Here, $\sqrt [*]{\zeta }$ denotes a choice for a branch cut of the function $f(\zeta )=\sqrt {\zeta }$ that maps the portion of the $\omega$-plane $\mathrm {Im}\{\omega \}\rightarrow \infty$ to $\mathrm {Im}\{\sqrt [*]{\zeta }\}>0$. This choice is important for the problem to be consistent with the time-dependent description and the inverse Laplace transform. For large positive growth rates this avoids pole contributions to the Bromwich contour (and thus to the plasma dispersion function), making the inverse Laplace transform well defined. To evaluate the latter it is convenient to deform the complex plane contour from its original position in the positive imaginary part of the plane downwards. Thus, in addition to the choice regarding the sign of $\mathrm {Im}\{\sqrt [*]{\zeta }\}$, it is convenient to construct a Riemann sheet by placing branch-cuts southwards.

To enforce the above, we choose two branch cuts in $\omega$-space originating from the critical points $\lambda =0$ and $\lambda =1$. The latter is, in addition to a branch point, also a singular point, $\sqrt [*]{\zeta }\sim 1/\sqrt {1-\lambda /2}$. The branch points represented will lead to some secular time dependence for damped modes (Kuroda et al. 1998; Sugama 1999) which we shall not explore further in this work. With these branch points localised, the natural branch-cut choice in figure 1(a) is problematic, as the branch cut emanating from $\omega =0$ is directly linked to the definition of $\lambda$, (3.7), and the vertical choice of the cut does not guarantee the physical requirement of localised $\mathrm {Re}\{\lambda \}>0$ modes. To avoid this, one must place the branch cut along the real line, as in figure 1(b). Although this changes the continuation of the Bromwich contour to the negative $\mathrm {Im}\{\omega \}<0$ part of the plane, it should not affect the description of unstable modes, which is our main concern here.

Figure 1. Branch cuts for $\sqrt [*]{\zeta }$. Two different Riemann sheets are shown (a,b) for $\sqrt [*]{\zeta }$ in complex $\omega$-space. Panels (a i,b i) and (a ii,b ii) show the real and imaginary parts of $\sqrt [*]{\zeta }$, respectively. The plots in (a) represent the natural Laplace continuation choice for the branch cuts, while (b) is the choice that represents localised solutions everywhere in the $\omega$ plane ($\mathrm {Re}\{\lambda \}>0$). The branch cuts are denoted by the red wiggly lines across which the function is discontinuous. The function has an integrable singularity at $\omega =4\omega _t^2/\omega _d$ as indicated in the text. Frequency is normalised to $\omega _t$ and $\omega _t/\omega _d=1/2$ is chosen for illustration, with the colormaps normalised for appropriate visualisation.

With either definition in our hands, we proceed and write

(3.12)\begin{equation} I_{{\parallel},\beta}(\zeta)=\frac{({-}1)^\beta}{\sqrt[*]{\zeta}} \partial_s^\beta \left[Z\left(\sqrt{s}\sqrt[*]{\zeta}\right)\right]_{s=1}, \end{equation}

where $Z(\cdot )$ is the well-known plasma dispersion function (Fried & Conte 2015) analytical continuation of the integral

(3.13)\begin{equation} Z(a) = \frac{1}{\sqrt{\pi}}\int_{-\infty}^\infty \frac{{\rm e}^{{-}x^2}}{x-a}\,\mathrm{d} x, \end{equation}

beyond $\mathrm {Im}\{a\}>0$. The introduction of $s$ as a dummy parameter in (3.12) follows from the application of the Feynman trick (Woods 1926, Ch. VI) to compute the integral in (3.11) in terms of (3.13).

We then have, after integration over $x_\perp$ (see the details in Appendix C),

(3.14)\begin{align} \mathcal{D} & = 1+\tau + F_0(b)\left[\left(1-\frac{\omega_\star}{\omega}+\frac{3}{2}\frac{\omega_\star^T}{\omega}\right)\sqrt[*]{\zeta}Z(\sqrt[*]{\zeta})-\frac{\omega_\star^T}{\omega}\zeta Z_+\right]-\frac{\omega_\star^T}{\omega}F_2(b)\sqrt[*]{\zeta}Z(\sqrt[*]{\zeta})\nonumber\\ & \quad +\frac{\omega_\star\bar{\omega}_d}{4\omega^2}\left\{F_2(b)\left[\eta\left(\frac{3}{2}+\zeta\right)-1+\left(1+2\zeta+\eta\left(2\zeta(\zeta-2)-\frac{3}{2}\right)\right)Z_+\right]\right.\nonumber\\ & \quad \left.+\,F_4(b) \eta \left[(1+2\zeta)Z_+-1\right]\vphantom{\frac{\omega_\star\bar{\omega}_d}{4\omega^2}}\right\}, \end{align}

where the Larmor radius functions are encapsulated in the functions

(3.15)\begin{gather} F_0(b)=\varGamma_0, \end{gather}

(3.16)\begin{gather}F_2(b)=(1-b)\varGamma_0+b\varGamma_1, \end{gather}

(3.17)\begin{gather}F_4(b)=2\left[(1-b)^2\varGamma_0+\left(\frac{3}{2}-b\right)b\varGamma_1\right], \end{gather}

with $\varGamma _n=\textrm {e}^{-b}\textrm {I}_n(b)$, and $\textrm {I}_n$ is the modified Bessel function of the first kind (Abramowitz & Stegun 1968, § 9.6). We use the shorthand notation $Z_+=-Z'/2=1+\sqrt [*]{\zeta }Z(\sqrt [*]{\zeta })$. In this form of the dispersion relation, we have presented all the terms that are, at least explicitly, order 1 or larger, taking $\omega _\star /\omega \sim \omega /\bar {\omega }_d$ in ordering. Additional terms in the expansion could be included (see table 1 in Appendix C), although these would stretch the original ordering in $\delta$ and $\epsilon$ and do not include any additional physics. This dispersion relation describes the behaviour of a linear localised ITG mode responding to a quadratic magnetic curvature with both a good and bad curvature. The dispersion is obtained through ordering of the localisation and the magnetic drift, which are taken to be strong and weak, respectively. In addition, to reach such simple form, we consider what we refer to as a ‘pure’ mode, which simplifies the localisation response of the mode.

Table 1. Contributions to the dispersion relation. Different term contributions to the dispersion relation $\mathcal {D}$ which may or may not be included depending on the required approximation. The columns $\alpha$ and $\beta$ denote the powers of $\bar {\omega }_d/\omega$ and $\omega _\star /\omega$, respectively, that these terms are multiplied by, Column $\gamma$ labels the FLR function that is multiplied by these terms, related directly to the number of powers of $x_\perp$ prior to integrating over $x_\perp$. The notation $Z_+(\sqrt {\zeta })=1+\sqrt {\zeta }Z(\sqrt {\zeta })$ for brevity.

3.3 Computing the dispersion function and finding unstable modes

Some elements of the dispersion relation are reminiscent of the common local, slab (Kadomtsev & Pogutse 1970a) or local, short wavelength ITG (Smolyakov, Yagi & Kishimoto 2002) limits. We extend on these by a more careful consideration on mode localisation.

The information about the linear stability of our localised mode is encoded in the solutions to the dispersion relation $\mathcal {D} = 0$. In general, this constitutes a complicated transcendental equation whose solutions need to be found numerically. The function $\mathcal {D}(\omega )$ can be evaluated numerically in $\omega$-space using the definition of $\sqrt [*]{\zeta }$ above, and using one of the many efficient implementations of the plasma dispersion function (in this case, we use the function wofz (Johnson 2024) in scipy.special). We accept as valid roots the values of $\omega$ for which $|\mathcal {D}|=0$, which holds true when both its real and imaginary parts vanish. In figure 2 we show some examples of what the dispersion function looks like in $\omega$-space for some particularly simple cases in which no Larmor-radius effects are present.

Figure 2. Dispersion function $\mathcal {D}$ in $\omega$ space. The plots show $|\mathcal {D}|$ as a function of complex $\omega$ for different combinations of $\bar {\omega }_d$ and $\omega _\star ^T$ (all frequencies normalised to the transit frequency $\omega _t=v_{Ti}/\varLambda \sqrt {2}$). The set of three panels can interpreted as the change in the instability due to an increase in $\varLambda$, the width of the bad curvature region, where the positive $\mathrm {Im}\{\omega \}$ part of the plane denotes instability: (a) $\omega _d/\omega _t=1\times 10^{-3}$ and $\omega _\star ^T/\omega _t=-1$; (b) $\omega _d/\omega _t=1\times 10^{-2}$ and $\omega _\star ^T/\omega _t=-10$; (c) $\omega _d/\omega _t=1\times 10^{-1}$ and $\omega _\star ^T/\omega _t=-100$. In this case we have chosen $b=0$, $\tau =1$ and $\omega _\star =0$ for simplicity. The red line represents one of the branch cuts; the vertical branch cut is not present in the domain plotted, as the unstable modes live near $\omega =0$ as shown.

A single unstable mode with $\mathrm {Re}\{\omega \}<0$ is seen in the plots of figure 2. As the drift and diamagnetic frequencies are varied, the mode location evolves in $\omega$-space. In fact, if we interpret the plots in figure 2(a–c) as the change due to increasing $\varLambda$, the width of the bad curvature region along the field line, we see that decreasing $\varLambda$ below a certain threshold appears to lead to the unstable mode eventually vanishing. The evolution of this instability with $\varLambda$ and the other parameters is what we are ultimately interested in, and shall be the main focus of our study in the following section. To automate the root-find of $\mathcal {D}$ and be able to study the various interesting mode dependencies, we implement the following approach: (i) we evaluate $|\mathcal {D}|$ in a coarse-grained $\omega$-space roughly bounded by $\omega _\star ^T$; (ii) we find the regions of lowest $|\mathcal {D}|$ residual; (iii) we perform local least-squares minimisation around them; and (iv) from the multiple roots found we select the most unstable one (see diagram in figure 3).

Figure 3. Diagram sketching the procedure to find the most unstable mode. This algorithm is used when numerical roots of $\mathcal {D}$ are required.

4 Physics and features of the localised kinetic ITG mode

In this section we investigate the behaviour of the localised ITG mode, including kinetic effects, by exploring (3.14). By inspection, the modes can depend on the following parameters: the diamagnetic drive $\omega _\star$ and $\omega _\star ^T$ (density and temperature gradients, respectively); the curvature drift magnitude $\bar {\omega }_d$; the transit frequency $\omega _t=v_{Ti}/\varLambda \sqrt {2}$; the poloidal wavenumber $k_\alpha$; and the ratio of electron to ion temperature $\tau$.

The frequencies and length scales will be normalised to a reference transit frequency, $\omega _{t0}$. That is, the frequency characteristic of the motion of a thermal ion over a distance $\varLambda _0$ along the field line, $\omega _{t0}=v_{Ti}/\varLambda _0\sqrt {2}$. We take this length scale $\varLambda _0$ to normalise $\varLambda$ as well. The normalisation of the poloidal wavenumber is somewhat more complicated. Let us recall the definition of $k_\alpha$ from the covariant form of the wavevector $\boldsymbol {k}_\perp =k_\alpha \boldsymbol {\nabla }\alpha +k_\psi \boldsymbol {\nabla }\psi$. We took for simplicity $k_\psi =0$, so that $\boldsymbol {k}_\perp =k_\alpha \boldsymbol {\nabla }\alpha$. The parameter $k_\alpha$ is dimensionless, leaving us with the uncomfortable situation of $k_\alpha \rho _i$ having units of length. When we write $k_\alpha \rho _i$ in what follows we adopt the convention of meaning $\overline {k_\alpha \rho _i}=k_\alpha \rho _i|\boldsymbol {\nabla }\alpha |$ (we will often drop the overbar notation, but it should be clear when we do so). The Larmor radius parameter becomes $b=(\overline {k_\alpha \rho _i})^2/2$, which does not exactly match the convention in other works in which the minor radius scale $a$ is chosen to normalise $k_\alpha \rho _i$.

We now construct the linear spectrum of the localised ITG mode as a function of $\overline {k_\alpha \rho _i}$, as the width of the bad curvature region is changed. Some examples are presented in figures 4 and 5. For these cases we have chosen a strongly driven scenario, $\omega _\star ^T/\bar {\omega }_d=10^2$, in the presence of no density gradient.

Figure 4. Evolution of the linear spectrum of unstable ITG with the bad curvature region size, $\varLambda$. The plots show (a) the growth rate and (b) the negative real frequency of the most unstable mode as a function of $k_\alpha \rho _i$ and $\varLambda$. The frequencies are normalised to $\omega _{t0}=v_{Ti}/\varLambda _0$, where $\varLambda _0$ is some reference length. We only plot points when the mode satisfies the conditions $\gamma >0$ and $|\omega _d/\omega |<1$. The plots are constructed for the choice $\omega _\star ^T/\omega _{t0}=-10$, $\omega _d/\omega _{t0}=0.1$, $\omega _\star /\omega _{t0}=0.0$ and $\tau =1$.

Figure 5. Properties of the unstable modes as a function of the bad curvature region size, $\varLambda$, and $k_\alpha \rho _i$ for $\gamma >0$ and $|\omega _d/\omega |<1$. The plots show (a) the growth rate $\gamma$, (b) the frequency $-\omega _r$, (c) the Gaussian envelope scale $|\lambda |$, (d) the small scale $|\omega _d/\omega |$, (e) the kinetic measure $|\zeta |$ and (f) the approximation scale $\mathrm {Re}\{\lambda \}|\omega /\omega _d|$. Panels (a,b) can be interpreted as top views of figure 4. The red broken line in (d) is the estimate of the Landau threshold as detailed in § 4.3. The blue region in (f) shows where we expect our localised mode approximation to break down. Here $\omega _\star ^T/\omega _{t0}=-10$, $\omega _d/\omega _{t0}=0.1$, $\omega _\star /\omega _{t0}=0.0$ and $\tau =1$.

4.1 General features and assessment of the approximation

The linear spectra obtained have four distinctive features that are most clear in the limit of a wide bad curvature region. At small $k_\alpha \rho _i$, as both the diamagnetic drive and the bad curvature diminish (given that they are linear functions of it), so does the instability, figure 4(a). The mode remains corotating with the diamagnetic frequency, figure 4(b), but eventually, for small enough frequency, it reaches what we may call the Landau threshold. This occurs when, at some critical $k_\alpha \rho _i$, the frequency matches the parallel transit frequency leading to Landau damping. As this threshold is approached and kinetics becomes more relevant ($|\zeta |\rightarrow 1$), the ITG mode structure tends to be increasingly stretched over the field line ($|\lambda |$ decreases). As it does so, the ever-increasing good curvature of our model (which serves a role similar to that which the magnetic shear would play) dictates its behaviour, as the key approximation $\mathrm {Re}\{\lambda \}|\omega /\omega _d|\sim 1/\epsilon \gg 1$ fails. This naturally lends to the possibility of the behaviour in this limit being dominated by delocalised, slab-like ITG modes such as Floquet modes (Zocco et al. 2016). The precise value of $k_\alpha \rho _i$ at which the threshold for the localised mode occurs cannot be expected to be exactly described by the model, but we may use its trends as an informative guide.

As the poloidal wavenumber increases, the diamagnetic and curvature drifts drive the ITG more vigorously. The growth rate increases as the mode becomes more localised; the gains from becoming localised about the bad curvature energy source overcome the effects of increasing the effective $k_\parallel$ that invigorates Landau damping. A larger $k_\alpha \rho _i$ does, however, also enhance FLR effects. At the same time, the stabilisation effects of the geometry can suppress the ITG mode. This leads to the small-$k_\alpha \rho _i$ peak ($k_\alpha \rho _i\lesssim 1$) in the spectrum, figure 4(a). One can think of this as the standard ITG peak, and refer to the threshold (or dip) to its right as the FLR stabilisation threshold.

Interestingly, increasing $k_\alpha \rho _i$ does not continue to reinforce the stabilising effect of the field geometry (in our case what one can refer to as flux compression $|\boldsymbol {\nabla }\alpha |\sim 1/|\boldsymbol {\nabla }\psi |$). There is some point at which the instability drive beats the FLR effects again, and the mode starts to become more and more unstable. This threshold we refer to as the FLR weakening threshold. As the mode keeps growing, it becomes increasingly localised near the bad curvature region, with its characteristic mode frequency rising sharply to settle to a roughly constant value, figure 4(b). Eventually, the mode reaches a last critical $k_\alpha \rho _i$ threshold. This corresponds to the situation in which the mode resonates with $\omega \sim -\bar {\omega }_d$, figure 5. We call this the $\omega _d$ threshold. Once again, as the kinetic effects grow and the threshold is approached, the fidelity of the model falters, in this case primarily because $|\bar {\omega }_d/\omega |\sim 1$. This leads to a second ‘hump’ in the linear spectrum of the ITG instability. Such behaviour has been previously studied in the context of the so-called short-wavelength-ITG (SWITG) instability in tokamaks (Hirose et al. 2002; Smolyakov et al. 2002; Gao et al. 2003).

These features and their location evolve as the width of the bad curvature region changes. As the bad curvature region narrows, the two instability peaks move towards each other, merge and eventually, below a certain critical $\varLambda$, disappear. That is, the localised mode is stabilised by sufficiently shortening the bad curvature region. Unfortunately, and as it occurred near the $k_\alpha \rho _i$ threshold, as $\varLambda$ becomes smaller, the mode becomes increasingly delocalised, and our description tends to break down, We may expect delocalised modes to gain prominence and persist in this limit (Zocco et al. 2015).

4.2 Fluid limit

Let us start a more quantitative analysis by checking what we obtain in the fluid limit, that is, when all resonances can be neglected. This will prove useful in two ways. First, because it serves as a good check that the dispersion relation reproduces a correct asymptotic limit. Second, because many of the features of our linear spectra may be explained in the simplest of terms through the fluid perspective, as we shall see.

With this in mind, our fluid limit will be derived from the $|\zeta |\gg 1$ expansion of $\mathcal {D}$ in (3.14), when the kinetic Landau resonance is negligible. Since (Fried & Conte 2015)

(4.1)\begin{equation} Z(x)\approx{-}\frac{1}{x}\left[1+\frac{1}{2x^2}+\frac{3}{4x^4}+\cdots\right], \end{equation}

for large argument, then

(4.2)\begin{align} \mathcal{D}& \approx 1+\tau-F_0(b)\left[1-\frac{\omega_\star}{\omega}(1-\eta)\right]+\frac{\omega_\star^T}{\omega}F_2(b)-\frac{1}{2\zeta}\left[F_0(b)\left(1-\frac{\omega_\star}{\omega}\right)-\frac{\omega_\star^T}{\omega}F_2(b)\right]\nonumber\\ & \quad -\frac{\omega_\star\bar{\omega}_d}{2\omega^2}\left[(1-\eta)F_2(b) + \eta F_4(b)\right], \end{align}

where

(4.3)\begin{equation} \frac{1}{2\zeta}=\left(\frac{\omega_t\sqrt{\lambda}}{\omega}\right)^2-\frac{\bar{\omega}_d}{2\omega}=\frac{\bar{\omega}_d^{1/2}\omega_t}{\omega^{3/2}}-\frac{\bar{\omega}_d}{2\omega}. \end{equation}

We obtain,

(4.4)\begin{align} \mathcal{D}& \approx (1+\tau-\varGamma_0)+\frac{\omega_\star}{\omega}\left[F_0(b)(1-\eta)+\eta F_2(b)\right] + \frac{\bar{\omega}_d\omega_\star}{\omega^2}\left[-\varGamma_0+\frac{b}{2}(\varGamma_0-\varGamma_1)-\frac{\eta}{2} F_4(b)\right]\nonumber\\ & \quad +\left(\frac{\bar{\omega}_d}{\omega}\right)^{1/2}\frac{\omega_t\omega_\star}{\omega^2}\left[F_0(b)+\eta F_2(b)\right], \end{align}

where we have dropped terms that are order $\bar {\omega }_d/\omega$. This is an algebraic equation whose roots can be found straightforwardly. Presented in this form of (4.4), the dispersion function may appear obscure, however, it agrees with Connor et al. (1980) (see Appendix D). This evidences the particular significance of $\lambda$, which provides our construction with the correct fluid limit.

It is commonplace to consider the case of a vigorously temperature-gradient driven limit (namely, $\omega _\star ^T/\omega \gg 1$), with small Larmor radius effects, $b\ll 1$. Then, a direct comparison with (Romanelli 1989, equation (19)) is possible,⁶

(4.5)\begin{equation} \tau\omega^3-b\omega_\star^T\omega^2-\bar{\omega}_d\omega_\star^T\omega+\omega_\star^T\omega_t\sqrt{\omega\bar{\omega}_d}\approx 0. \end{equation}

The fluid limit of the dispersion function in (3.14) is therefore correct. Given this connection, we may extend the common nomenclature distinguishing between the slab and toroidal ITG modes interpreted as follows (Wesson & Campbell 2011; Zocco et al. 2016). Let the mode in which the streaming terms are negligible⁷ be referred to as the toroidal ITG mode; and let the slab ITG mode correspond to the reverse. Formally, this distinction is dictated by the relative importance of the last two terms in (4.5), which is simply $|\lambda |$. Thus, when the mode has significant structure within the bad curvature well ($|\lambda |>1$), we shall refer to the mode as toroidal (note that this mode is not necessarily strongly localised, which would be a statement about $\mathrm {Re}\{\lambda \}$). The slab modes will tend to be delocalised and thus care mostly about the larger $|\bar {\ell }|$ good curvature part of the problem.

4.3 The Landau threshold

We discussed qualitatively the presence of a cutoff at long wavelengths $k_\alpha \rho _i\ll 1$ in the linear ITG mode spectrum. We argued that this threshold was indicative of a stabilisation of the mode by Landau damping; that is, that the localised mode near the threshold becomes resonant with the parallel streaming frequency. Although very near the threshold the construction of the dispersion relation in (3.14) is not formally valid and delocalised modes may exist, we may nevertheless use the model to estimate where this threshold for the localised modes occurs, and how it is affected by the various properties of the field.

We thus need to find a real frequency solution to the equation $\mathcal {D}(\omega )=0$ (i.e. the mode that as in figure 2 is just touching the $\omega$-space real axis). Typically, such as in the standard local slab ITG, the real nature of the frequency makes the arguments of the plasma dispersion functions real, and separating the real and imaginary parts of the dispersion relation becomes straightforward, without the requirement of solving any transcendental equation. In the present case, though, the argument of the plasma dispersion function, (3.14), involves the square root of $\omega$. As the mode is driven by the temperature gradient, $\omega _\star ^T$, $\omega <0$ at the threshold, and thus $Z(\sqrt {\zeta })$ has both mixed real and imaginary parts. This makes solving a transcendental equation unavoidable.

We thus proceed by employing a physically motivated approach. As mentioned above, the Landau threshold corresponds to the resonance of the mode with the transit frequency, $\mathrm {Re}\{\omega \}\sim \omega _t.$ To be more precise, we have $\zeta \approx 1$, and considering the mode to be slab-like in the region near the Landau threshold (i.e. rather elongated along the field line compared with the bad curvature region $\varLambda$), we use, following (3.9a),

(4.6)\begin{equation} \mathrm{Re}\{\omega\}\approx\omega_t\sqrt{2\lambda}=\frac{v_{Ti}}{\varLambda/\sqrt{\lambda}}, \end{equation}

where $\varLambda /\sqrt {\lambda }$ is the characteristic longitudinal scale of the mode structure. The dependence of $\lambda$ on $\mathrm {Re}\{\omega \}$ is known from (3.7). So what we ultimately need is some notion of the mode frequency, $\mathrm {Re}\{\omega \}$.

To estimate $\mathrm {Re}\{\omega \}$, we crudely assume that the kinetic effects mainly affect the stability of the mode, but leave the mode frequency to a large extent unaltered. Assuming a sufficiently large $\varLambda$ and $\omega _\star ^T/\bar {\omega }_d$, we use the slab branch of the fluid model to estimate $\mathrm {Re}\{\omega \}$, from (4.5),

(4.7)\begin{equation} \mathrm{Re}\{\omega\}\approx(k_\alpha\rho_i)^{3/5}\hat{\omega}_d^{1/5}\left(\frac{\omega_t\hat{\omega}_\star^T}{\tau}\right)^{2/5}\cos\left(\frac{4\pi}{5}\right), \end{equation}

where the root with a negative real frequency must be chosen given the branch cut choice of the square root. The hat notation is meant to indicate that we have taken the $k_\alpha \rho _i$ dependence of $\bar {\omega }_d$ and $\omega _\star ^T$ out explicitly. Although this might appear complicated, and the one-fifth powers odd, the expression just found is precisely the frequency of a standard slab ITG mode with $k_\parallel =\sqrt {\lambda }/\varLambda$,⁸ the characteristic length scale of the mode. The main difference is that in our problem we have explicit field line dependence, while the pure slab ITG does not (other than an oscillating delocalised solution).

With this expression for the real frequency, we can reconsider (4.6), to write

(4.8)\begin{equation} (k_\alpha\rho_i)_\mathrm{Landau}\approx12.5 \tau\frac{\omega_t}{\hat{\omega}_\star^T}\left(\frac{\tau\hat{\omega}_d}{\hat{\omega}_\star^T}\right)^{1/2}. \end{equation}

As shown in figure 5, this is a fair estimate of the Landau threshold for the parameters considered, and most importantly, it shows the correct $\varLambda$ scaling. The threshold value increases as the bad curvature region becomes smaller, $k_\alpha \rho _i\sim 1/\varLambda$. Narrowing down the bad curvature well produces a narrower mode structure (the mode width goes like $\Delta \ell \sim \varLambda ^{3/5}$), Landau damping becomes more effective, and thus the threshold increases. If the instability is driven more vigorously (i.e. we increase the temperature gradient drive, $|\omega _\star ^T|$), then the parallel dynamics has a harder time stabilising the mode. The dependence on the temperature gradient of (4.8) is also in agreement with the numerical solution (see figure 6). Interestingly, increasing the curvature drift, and with it the magnitude of the bad curvature, does not worsen the threshold, but rather improve it. The reason is that the effect of $\bar {\omega }_d$ in narrowing the mode is, in this limit, stronger than its direct drive of the instability. In this limit in which the mode is rather delocalised, both $\bar {\omega }_d$ and $\varLambda$ should be interpreted to represent more ‘global’ properties of the geometry, as the (global) shear would.

Figure 6. Properties of the unstable modes as a function of the temperature-gradient driven diamagnetic frequency, $\omega _\star ^T$, and $k_\alpha \rho _i$. The plots show (a) the growth rate $\gamma$, (b) the real frequency $-\omega _r$, (c) the Gaussian envelope scale $|\lambda |$, (d) the small scale $|\bar {\omega }_d/\omega |$, (e) the kinetic measure $|\zeta |$ and (f) the approximation scale $\mathrm {Re}\{\lambda \}|\omega /\bar {\omega }_d|$. The red broken line in (d) is the estimate of the Landau threshold as detailed in § 4.3. The blue region in (f) shows where we expect our localised mode approximation to break down. This means that the precise instability threshold in $\varLambda$ cannot be fully trusted. We only plot points when the mode satisfies the conditions $\gamma >0$ and $|\bar {\omega }_d/\omega |<1$. The plots are constructed for the choice $\bar {\omega }_d/\omega _{t}=1.0$, $\omega _\star /\omega _{t0}=0.0$ and $\tau =1$.

Of course, near this threshold, the possibilities of other delocalised modes dominating the dynamics increases (Zocco et al. 2018, 2022; Podavini et al. 2024). In addition, the exact occurrence of the threshold will depend on the precise form of $\lambda$, which we have already acknowledged the current model to only approximate. Thus, the particular scaling of the critical $k_\alpha \rho _i$ and its value may change, but the main physical interpretation of its origin and dependence should remain.

4.4 The FLR stabilisation

The next noteworthy feature in the $k_\alpha \rho _i$ spectrum is the stabilisation of the ITG mode as the FLR becomes relevant. The result is the appearance of a characteristic peaked linear spectrum, with typical values of $k_\alpha \rho _i\lesssim 1$. This feature is not resonant-kinetic, and its presence in figure 5 for large $\varLambda$ indicates it forms part of the fluid description of the instability.

In the fluid picture, the FLR stabilisation feature of the linear ITG spectrum results from the competition between the increased drive of the toroidal ITG and the increased stabilising efficiency of FLR effects with increased poloidal wavenumber. Ignoring the streaming term in (4.5), the threshold can be shown to occur when

(4.9)\begin{equation} \left(\frac{\omega_\star^T b}{2\tau}\right)^2+\frac{\omega_\star^T\bar{\omega}_d}{\tau}\approx0. \end{equation}

The critical poloidal wavenumber is then

(4.10)\begin{equation} (k_\alpha\rho_i)_\mathrm{FLR}\approx2\left(\frac{\tau\bar{\omega}_d}{\omega_\star^T}\right)^{1/4}. \end{equation}

From the physical picture above, we expect this stabilisation threshold to have a value close to $k_\alpha \rho _i\sim 1$. A sign of this is the weak dependence on both the drift and the temperature gradient. In addition, this mechanism is not directly linked to the mode structure, and hence independent of $\varLambda$. This latter resilience is apparent in figure 5. Only a small correction term may be observed within the fluid description of the mode by treating perturbatively the streaming contribution in (4.5), $\delta (k_\alpha \rho _i)_\mathrm {FLR}\approx \tau \omega _t/2|\omega _\star ^T|$.

This resiliency of the FLR stabilisation threshold clashes with the rather fundamental dependence of the Landau threshold on $\varLambda$. As a result, we expect the space in $k_\alpha \rho _i$ available to the localised ITG mode to narrow down as the bad curvature region is squeezed. Eventually, the first peak in the linear spectrum would be eliminated, as in figure 4(a), and thus at long wavelengths only extended ITG modes could remain present in the system. From the above, we may estimate the critical width of the bad curvature region, $\varLambda$, at which this occurs. Extrapolating the behaviour of the threshold and the FLR stabilisation, (4.8) and (4.10), the balance $(k_\alpha \rho _i)_\mathrm {Landau}\sim (k_\alpha \rho _i)_\mathrm {FLR}$ yields

(4.11)\begin{equation} \varLambda_\mathrm{crit}\sim6\tau\frac{v_{Ti}}{\hat{\omega}_\star^T}\left(\frac{\tau\hat{\omega}_d}{\hat{\omega}_\star^T}\right)^{1/4}. \end{equation}

This critical $\varLambda$ could also be rewritten in terms of a critical temperature gradient $(\omega _\star ^T)_\mathrm {crit}\sim 4\tau (\omega _t^4\bar {\omega }_d)^{1/5}$. Given the pre-eminence of the Landau damping physics, the critical temperature gradient threshold below which only extended modes are left at long wavelengths is particularly sensitive to the width of the bad curvature region. This emphasis in controlling the longitudinal spread of the mode aligns with some previous work on critical thresholds (Jenko, Dorland & Hammett 2001; Roberg-Clark, Plunk & Xanthopoulos 2022a).

4.5 The FLR weakening

The above considerations of the kinetic threshold and FLR stabilisation explain the presence of a peak in the linear spectrum of the mode in $k_\alpha \rho _i$. However, as is clear from figure 4, the ITG starts growing unstable once again at larger values of the poloidal wavenumber. This might appear surprising at first, because it is natural to think of FLR effects increasing monotonically with $k_\alpha$, and thus after the FLR threshold, its stabilising effect to continue to exceed the increase in the turbulent drive. However, this picture is not correct.

The FLR effects become weakened at large $k_\alpha \rho _i$, when the relevant perpendicular scale starts to become significantly smaller than the Larmor radius. Larmor radius effects become inefficient, as they were for small $k_\alpha \rho _i$. Formally, one can ascribe this weakening of the FLR effects to the large argument behaviour of the Bessel functions introduced in the GK equation by the gyroaveraging. The small $b$ expansion of the $F_n(b)$ functions fails in this limit, and thus so does the fluid equation written in its typical form of (4.5). Retaining the appropriate behaviour leads to a critical FLR weakening threshold beyond which the mode grows back up.⁹ This type of ITG activity is referred to in the literature as SWITG (Hirose et al. 2002; Smolyakov et al. 2002; Gao et al. 2003).

We shall then consider the large $b$ limit of the fluid equation, (4.4) using the large argument asymptotics of the Bessel functions (Abramowitz & Stegun 1968, § 9.7),

(4.12)\begin{align} \sqrt{2\pi b}(1+\tau)+\frac{\omega_\star}{\omega}\left(1-\frac{\eta}{2}\right)-\frac{3}{4}\frac{\omega_\star\bar{\omega}_d}{\omega^2}\left(1+\frac{\eta}{2}\right)+\left(\frac{\bar{\omega}_d}{\omega}\right)^{1/2}\frac{\omega_t\omega_\star}{\omega^2}\left(1+\frac{\eta}{2}\right)\approx0. \end{align}

In the limit of $|\lambda |$ being large (which is indeed the behaviour at large $k_\alpha$, see figure 5), we expect the streaming term in (4.12) to be small, and thus we are left with, once again, a quadratic in $\omega$. The threshold then occurs when its discriminant vanishes, which gives

(4.13)\begin{equation} (k_\alpha\rho_i)_\mathrm{weak}\approx\frac{1}{3\sqrt{\pi}(1+\tau)}\frac{\omega_\star}{\bar{\omega}_d}\frac{(1-\eta/2)^2}{1+\eta/2}. \end{equation}

A more refined version of this threshold keeping higher orders in $1/\zeta$ yields for a flat density $(k_\alpha \rho _i)_\mathrm {weak}\approx (0.642-0.032\omega _\star ^T/\bar {\omega }_d)/(1+\tau )$. The exact numerical value here is, however, not important, especially given the ordering in $\delta$ and $\epsilon$ considered. The key is its dependence on the relative magnitude of the diamagnetic and curvature drift, which is quite strong compared with the FLR threshold $(k_\alpha \rho _i)_\mathrm {FLR}\sim (\bar {\omega }_d/\omega _\star ^T)^{1/4}$, (4.10). Then, as the curvature of the field is increased or the driving temperature gradient reduced, the weakening FLR threshold will approach the FLR stabilisation threshold. Through a simple balance between $(k_\alpha \rho _i)_\mathrm {FLR}\sim (k_\alpha \rho _i)_\mathrm {Weak}$, we expect to find the two regions merging (see figure 4) when $|\omega _\star ^T|/\bar {\omega }_d\sim 27(1+\tau )^{4/5}\tau ^{1/5}$. This is a rather large value of $\omega _\star ^T$, meaning that having two distinct peaks in the linear spectrum requires a strongly driven regime (compared with the drift).

In figure 6 we illustrate some of these linear spectra dynamics as a function of a changing temperature gradient.

4.6 The $\omega _d$ threshold

At even larger $k_\alpha \rho _i$ we clearly have another stabilisation effect that leads to the appearance of an instability threshold. The responsible mechanism is not captured in the fluid picture, and is kinetic in nature, as the value of $|\zeta |$ in figure 5 suggests. It could be tempting from our previous discussion on the Landau threshold to suggest that a similar Landau damping mechanism is present here as well. However, as the threshold is approached, the ITG mode does not seem to exhibit a clear tendency to relax its longitudinal structure. In fact, as $|\lambda |\gg 1$ the mode retains a fine structure. What is then the mechanism in action? To answer the question it suffices to look back at the definition of $\zeta$ in (3.9a), which in this limit gives $\zeta \approx -\omega /\bar {\omega }_d$. Thus, kinetic effects at large $k_\alpha \rho _i$ are dominated by the resonance of the mode rotation with the bad curvature of the field. Note that this resonance has been retained even in our expansion in $\bar {\omega }_d/\omega \sim \delta \ll 1$, as can be recognised by inspection of our kinetic equation (3.6b).

With the dominant mechanism identified, we may estimate at what wavenumbers the mode frequency balances the drift frequency. The large $k_\alpha \rho _i$ limit of the fluid equation (4.12) yields a roughly constant real frequency for the mode, which will eventually be matched by the magnetic drift, which grows with the wavenumber $\bar {\omega }_d\propto k_\alpha \rho _i$. Using the higher-order cubic version of (4.12), the balance $\mathrm {Re}\{\omega \}\sim -\bar {\omega }_d$ gives

(4.14)\begin{equation} (k_\alpha\rho_i)_{\bar{\omega}_d}\approx\frac{1}{2}\frac{1}{1+\tau}\left|\frac{\omega_\star^T}{\bar{\omega}_d}\right|. \end{equation}

The multiplicative factors in front depend on the details of the approximation of the model, but what is key is the dependence of the threshold on the ratio $\omega _\star ^T/\bar {\omega }_d$ (which correctly describes the behaviour in figure 6). The behaviour of this threshold suggests a narrowing of the linear spectrum that shall reach a critical narrow range $(k_\alpha \rho _i)^\mathrm {crit}$ when

(4.15)\begin{equation} \left.\frac{\omega_\star^T}{\bar{\omega}_d}\right|_\mathrm{crit}\approx2(k_{\alpha}\rho_{i})^\mathrm{crit}\left(1+\tau\right), \end{equation}

or in tokamak notation $R/L_T$, which is compatible with the results of Romanelli (1989) and Guo & Romanelli (1993), with a visual estimate obtained from figure 2 of Biglari, Diamond & Rosenbluth (1989), with the Jenko–Dorland–Hammett formula (Jenko et al. 2001) and other critical gradient estimates (Roberg-Clark et al. 2022a).

In this consideration there appears not to be any direct involvement of the parallel streaming dynamics. However, as the curvature well is narrowed the two kinetic elements in the problem become mixed. Signs of this behaviour are seen in the apparent decorrelation between the $\omega /\bar {\omega }_d$ ratio and the threshold in figure 5 when $\varLambda$ starts having an effect on the mode. In fact, in the above description of a critical temperature threshold, we considered some reference critical value of $(k_\alpha \rho _i)_\mathrm {crit}$. Following figure 6 though, one can expect at some point the $\bar {\omega }_d$ threshold to become close to the Landau threshold, and not just an arbitrarily chosen wavenumber. When this occurs, we may say that there will be no more localised ITG modes. In this case the threshold would scale as

(4.16)\begin{equation} (\omega_\star^T)_\mathrm{crit}\sim(\bar{\omega}_d^3\omega_t^2)^{1/5}. \end{equation}

Thus, the threshold changes because we can affect it not only by making the resonance with $\bar {\omega }_d$ appear at longer wavelengths, but also by amplifying the effects involved in the Landau threshold. As a result, increasing the bad curvature or the parallel scale both will have a positive effect on reducing the ITG, although extended modes may persist. The particular scaling obtained by balancing the two physics ingredients goes beyond Biglari et al. (1989) and $R/L_T \sim O(1)$, where $R$ is the major radius. In our case, the balance yields $R/L_T\sim 1/q^{0.6}$, involving the safety factor $q$, in relation to the connection length. The importance of the parallel physics in determining the critical gradients has been recognised by many authors (Hahm & Tang 1988; Jenko et al. 2001; Roberg-Clark et al. 2022a), and here we see is involved quite explicitly. Often parallel dynamics are associated with the global shear, which our model does not explicitly treat. However, as discussed with the Landau threshold, the dependence on $\omega _t$ can be related to the role played by global shear when the modes become rather delocalised. The exact form of the scaling with $q$ will depend on the exact behaviour of the mode that is becoming delocalised and thus should be taken with a pinch of salt (especially given the weak spot of the model determining $\lambda$). We shall not forget that although increasing the Landau threshold helps in this endeavour of erasing localised modes, one could of course leave behind delocalised modes that could also be deleterious (Zocco et al. 2018, 2022).

5 The role of higher harmonics

The kinetic considerations above have focused on the behaviour of a localised mode whose shape is in its simplest form described by a Gaussian envelope. That is, by the ‘pure’ zeroth order in our Taylor–Gauss expansion. Shapes of the modes are seldom so simple, and in the way that we expect multiple modes as solutions to a Schrödinger equation (and in fact also the fluid equation (Hahm & Tang 1988)), we may also expect to find possible solutions to our problem in which the mode has larger $n$.

A dominant ‘pure’ $M$ mode can be described in a way rather analogous to that considered above for $n=0$. It may be shown, see Appendix B, that the resulting equations are identical to the $n=0$ case with the only modification,

(5.1)\begin{equation} \zeta_{n=M}=\frac{1}{2}\left[\lambda(2M+1)\left(\frac{\omega_t}{\omega}\right)^2-\frac{\bar{\omega}_d}{2\omega}\right]^{{-}1}. \end{equation}

The only change is the different contribution of the streaming term, formally equivalent to an increased $\omega _t\rightarrow (2M+1)\omega _t$. How can such a scaling be physically interpreted? Let us picture the change to the mode structure as one increases mode number. The mode looks like a pair of peaks pushed against the Gaussian envelope. The larger the $M$, the harder they are squeezed, which leads to the width of the peaks to roughly go like $\Delta \bar {\ell }\sim 1/\sqrt {M}$, as can be explicitly shown by computing the standard deviation. This linear scaling is what leads to the form in (5.1).

The main effect of $M$ is thus to change the parallel structure of the mode, and thus any of the phenomena that are directly linked to this feature of the instability. The Landau threshold will be most readily affected, but also the particulars of critical thresholds where $\zeta \sim 1$ and $|\lambda |$ is not exceedingly large. With the effective $\omega _t$ scaling, we may write $(k_\alpha \rho _i)_\mathrm {Landau}\sim 2M+1$, showing that Landau damping becomes more prominent for the higher modes. This means that we expect to see the lowest modes excited first, as well as those to be most sensitive to stabilisation through squeezing of the connection length, $\varLambda _\mathrm {crit}\sim 2M+1$. The change in the $\omega _d$ threshold with $M$ may appear shocking, given that this we argued is rather insensitive to parallel scales. However, as $M$ is increased, the $\omega _d$ threshold mixes with the Landau threshold, as follows directly from $\zeta$, (3.9a). Similar behaviour to what we find here, figure 7, may be interpreted from the work of Gao et al. (2005). Thus we expect close to marginality the first localised mode to go unstable to be the $n=0$ mode.

Figure 7. Growth and frequency of the ITG mode for different structure. Plots showing the growth rate and real frequency of the ITG mode for the mode numbers $n=0,1,2$ using the approximate generalisation of $\zeta$. The plots are computed using $\varLambda /\varLambda _0=10$, $\bar {\omega }_d/\omega _t=0.1$, $\omega _\star ^T/\omega _t=-30$, $\omega _\star =0$ and $\tau =1$.

To study the changes on the magnitude of the growth rate, we write the $M$th mode generalisation of (4.5), which can be shown to match other treatments of the fluid equation (Hahm & Tang 1988). In this fluid case, it is clear that increasing mode number enhances the growth rate of our ITG, which is increasingly of a more slab character, $\gamma _\mathrm {slab}\sim (2M+1)^{2/5}$. The fluid picture is unbound! Only through the regularising role played by the kinetic effects is the hierarchy regulated. There is a competition then between larger modes tending towards larger growth rates (in the fluid limit), but also becoming more effectively stabilised by Landau damping. The results of such competition are presented for some example parameters in figure 7.

These modifications that occur in the model provide us with a flavour of the kind of changes that one would expect when our ‘pure’ mode assumption is relaxed and the true value of $\lambda$ is different. While the main qualitative features shall remain, we expect quantitative differences to exist. Overall scalings can be understood as in figure 7, but other more sophisticated dependences of $\lambda$ on $k_\alpha \rho _i$ could also lead to additional features in the linear spectrum. The model should help us distinguish between these as well.

6 Qualitative comparison with simulations

In this section we consider by way of an example the linear spectrum of ITG modes in a realistic but nonetheless simple stellarator geometry. We use this as way of illustration of how the lessons learnt from our model can be applied in practice to understand behaviour in more complex situations, but also its limitations. The example presented is the linear ITG mode spectrum along a flux tube of the HSX (Anderson et al. 1995) stellarator, a quasisymmetric stellarator (Boozer 1983; Nührenberg & Zille 1988; Rodríguez et al. 2020) with helical symmetry. The latter means that $|\boldsymbol {B}|$ has a direction of symmetry to prevent fast loss of trapped particles, which implies a particularly simple, quasiperiodic curvature along field lines. This simplicity, together with a reduced global magnetic shear, makes this example suitable for the comparison. The electrostatic mode in the GK simulations is driven by an ion temperature gradient $a/L_T=2.5$, keeping the density flat and treating electrons adiabatically. The linear GK simulations are conducted with the STELLA code (Barnes, Parra & Landreman 2019), of which the details are available at https://doi.org/10.1017/S0022377824001120. The resulting linear spectrum is presented in figure 8, where the poloidal wavenumber and frequencies have been normalised following the notation in this paper. For comparison with this HSX simulation, an equivalent simulation has also been carried out with a modified geometry in which all field quantities are constant along the field line except the curvature, which is modelled like a single quadratic well (modelled with parameter $\bar {\omega }_d$ corresponding to the value of bad curvature at the bottom of the central well, and $\varLambda$ as the distance from the centre of the well to the first zero curvature crossing point), truncated at a finite good curvature to prevent spurious modes as discussed later.

Figure 8. Linear mode spectrum for HSX GK simulations. The growth rates (a) and real frequencies (b) are shown for linear HSX GK simulations with $a/L_T=2.5$ and $a/L_n=0$, and adiabatic electrons. The dashed line corresponds to the simulations performed with the same gradients but a modified geometry in which the only spatial dependence in the problem is $\bar {\omega }_d$, and this is modelled as a truncated quadratic well. The coloured shade show the variation in the latter half of the simulation of the mode frequency, giving a sense of trust of the spectra (blue and red for the HSX and model geometries, respectively). The vertical lines correspond to the predicted FLR stabilisation threshold and the $\omega _d$ stabilising resonance.

The first noticeable conclusion from the comparison between these two simulations (solid and dashed lines in figure 8) is that the simplified geometry appears to capture the behaviour of the ITG mode exceptionally well. This strongly backs the approach and geometric approximations considered in this paper, and supports our view on the key role played by the curvature in localising the ITG mode. Thus, we are allowed to describe realistic linear stellarator spectra from our theoretical perspective, even though we notice that the ‘pure’ mode assumption (and perhaps others such as a vigorous drive) does not provide a good quantitative description, in particular requiring somewhat larger temperature gradients in order to match magnitudes of mode frequency and growth rate of simulation results. We show an example of this lack of agreement in figure 9, where we include the analytical model prediction for the parameters characterising HSX, and one with a $27\,\%$ larger temperature gradient.

Figure 9. Comparison of the simulation to the analytic model showing quantitative discrepancies. The plot shows a comparison between the growth rate (a) and frequency (b) of the dominant linear ITG mode in the simulation of figure 8 (solid line) and the analytic model developed in the paper. The dashed line corresponds to the model prediction for $\bar {\omega }_d/\omega _t=1.95$ and $\omega _\star ^T/\omega _t=-18$, parameters obtained from the main well of the HSX geometry. The dot–dash line corresponds to the analytical prediction with a $27\,\%$ larger temperature gradient, with the shade representing ${\pm }5\,\%$ variation. This shows that the model suffers as a quantitative predictor. This suggests consideration of the model mainly as a physical qualitative framework to interpret linear spectra behaviour.

We can still investigate the key physical principles we have explored through our model to interpret the linear spectrum observed. Let us start from the rightmost part of the spectrum, where we expect to find our $\omega _d$ threshold. In fact, and as shown in figure 8 and predicted, this stabilisation point does occur approximately at $\mathrm {Re}\{\omega \}\sim -\bar {\omega }_d$ (numerically within $6\,\%$). Our interpretation of the nature of this threshold enables us to understand how the spectrum should change as the bad curvature of the field or the temperature gradient are varied. The spectrum would narrow as the curvature is increased or the gradients decreased. Such a simple perspective can explain the narrowing of spectra observed in recent efforts to optimise stellarators for improved temperature gradient thresholds (Roberg-Clark, Xanthopoulos & Plunk 2024b). In the comparison between the full geometry simulation and that of the reduced geometry, we see that the latter seems to exhibit an additional unstable branch at smaller poloidal wavelengths. This branch corresponds to an anti-ionic-temperature-driven instability (i.e. rotating in the electronic direction, but not due to kinetic electrons, since their response is adiabatic) that exists in the region of good curvature of the modified geometry. This behaviour is an interesting case study for the future: if anti-ionic modes localise in good curvature regions, attempting to stabilise the ITG mode by increasing the good curvature in the field could be problematic.

The region beyond the $\omega _d$ threshold in the HSX case considered here does not exhibit other significant instabilities. However, less symmetric geometries such as those in quasi-isodynamic stellarators (Podavini et al. 2024), often exhibit linear spectra with additional structure beyond this point. How can this be framed within our description? That is: How can the ITG mode escape stabilisation by the $\omega _d$ resonance? It is not the localisation of the mode that is suppressing the mode, as it may occur at small $k_\alpha \rho _i$, thus delocalising the mode is not a solution.¹⁰ The alternative left is for the mode to localise in a different well, what we might call a hopping mode. Moving to a well that is a priori less unstable because it has better curvature, or milder FLR effects, may, however, be beneficial for the mode because it may be sufficient to make $\mathrm {Re}\{\omega \}\neq -\bar {\omega }_d$ and avoid stabilisation through the $\omega _d$ resonance. Such a jump of the localised mode could occur more than once and we hypothesise, could lead to additional growth rate peaks in the spectrum. These changes on localisation are reminiscent of changes that occur to modes with finite $k_\psi$ (Parisi et al. 2022). The global behaviour of the localised mode including the possibility of living at different wells could be understood by constructing one spectrum like that in figure 8 for each well, with the corresponding change in the model parameters; more precisely, the magnitude of the drift, the well width and the FLR effect (with the latter strongly affecting the scale of $k_\alpha \rho _i$). The study of this is left for future work, but is an attractive way forward to understanding the behaviour of kinetic ion instabilities.

Our physics interpretation of the occurrence of the dip in the linear spectrum of figure 8 from FLR stabilisation and then weakening appears to be also correct. The dotted line in the figure corresponds in fact to the simple expression derived above in (4.10). This relative stabilisation of Larmor radius effects can also be observed in the response of the mode eigenfunction. We present some example structures from the simulations in figure 10. We use the structure from the STELLA simulations of the simplified geometry for better clarity. The figure also presents a measure of the localisation of the mode, $\lambda _r$, and its structure, $\lambda _i$. Numerically, we compute $\lambda _r$ by fitting a Gaussian $\exp [-\lambda _r\bar {\ell }^2/2]$ to $|\phi |$. Note that the wavefunctions are not really pure Gaussians; in fact, they appear to fit closer to a Gaussian in the centre and a weaker exponential decay farther away. Of course the precise behaviour depends on the details of the geometry. Larger $\lambda _r$ denotes a more localised mode. If one takes $\mathrm {Re}\{\phi \}/|\phi |$, the mode then exhibits a clear periodicity (other than in the middle region). We define $\lambda _i$ to be the main frequency of the that oscillatory behaviour in $\bar {\ell }$. As it is clear in figure 10, as predicted by our treatment, as FLR stabilises the ITG, the mode becomes less localised, while keeping the scale of its mode structure (i.e. $\lambda _i$ unchanged). Note also that the behaviour near the $\omega _d$ resonance gives us some intuition that the ‘pure’ mode form of $\lambda$ does not fully hold there, as the mode only weakly delocalises.

Figure 10. Mode structure for HSX GK simulations. The plot shows the parameters $\lambda _r$ and $\lambda _i$ describing the localisation and oscillation of the modes in figure 8. Here $\lambda _r$ is obtained by fitting a Gaussian $\exp [-\lambda _r\bar {\ell }^2/2]$ to the absolute value of the electrostatic potential $\phi$. The mode is not a pure exponential with generally wider tails. Three examples of the modes in the simplified geometry are provided in (b). The $\lambda _i$ parameter is obtained by reading the main oscillation frequency of $\mathrm {Re}\{\phi \}/|\phi |$. Both these measures are inspired by the Gaussian basis used in this paper. The qualitative behaviour observed is fully consistent with the behaviour of our model. That is, $\lambda _i$ monotonically increases with the poloidal wavenumber while the localisation decreases near the Landau threshold ($k_\alpha \rho _i\sim 0$), the FLR stabilisation threshold and the drift resonance.

We are then left with the behaviour at the longest poloidal wavelength. As predicted by our theory, the mode becomes increasingly delocalised at lower $k_\alpha \rho _i$. In terms of our Landau threshold picture, this occurs while Landau stabilisation becomes less and less effective. A causality relation between the lack of localisation and growth rate maximisation could perhaps be established via an energetic argument. Of course, at some point the model ceases to be valid, and the linear spectrum becomes dominated by completely delocalised modes. These very extended modes we may refer to as Floquet or slab-like modes (Zocco et al. 2018, 2022; Podavini et al. 2024). The difference in growth rate between the two GK simulations in figure 8 at long wavelengths can be explained to be due to the pre-eminence of these modes. The presence of such elongated structures make the simulations rather challenging, as long flux tubes are necessary. Their behaviour lies outside the realm of the present model.

We thus understand the very distinct nature of the ITG mode for large and small wavelengths, where the kinetic stabilising mechanisms are different, and thus so is the mode response. This enables one to elucidate the meaning of the linear spectrum presented, and offers, as given, a way forward to interpreting linear spectra and their physical meaning in more complex geometries.

7 Conclusions

In this work, we have proposed a theory of the kinetic ITG-driven mode which features resonant kinetic effects, with a localised mode structure induced by the field-line-dependent geometry. We focused on the localising action of the magnetic drift, allowing for general conclusions that apply both in a tokamak and stellarator context.

The magnetic drift spatial dependence models good and bad curvature regions, introduced with a local well quadratic model. The mathematical description of the problem is based on a power series expansion of the eigenfunctions in the field-following coordinate, mitigated by a Gaussian envelope. This generates a hierarchy of coupled eigenvalue problems which for small magnetic drift frequency, strong drive and localised modes, can be truncated. A relatively simple dispersion relation is constructed by considering the simplifying assumption of a singly dominant mode, which we refer to as ‘pure’. The resulting description features long-wavelength Landau damping, arbitrary Larmor radius effects and a regularising resonant action of the magnetic drift for short wavelength modes. All these salient physical features are demonstrated to be numerically observed in realistic stellarator geometry, although the simplicity of the model limits quantitative comparisons with numerical results. Venues for improving the model are also proposed for future work.

The model is also used to provide a prediction for the resonant stabilisation of the toroidal branch of the ITG mode without tacitly assuming constant (along the field line, and thus unrealistic) eigenfunctions. The result is insight into how to tailor the field-line dependence of the magnetic drift in order to suppress the ITG instability. By enforcing that all scales are either Landau damped (at long wavelengths) or kinetically suppressed by the magnetic drift resonance (at short wavelengths) one obtains a critical gradient for ITG destabilisation that scales with $a/L_T\propto 1/(q^{0.6}R).$ This explains why large inverse aspect-ratio devices feature small critical thresholds (as is well known), but also indicates a beneficial effect in having a small safety factor $q$, or more generally short connection length. This final aspect is of particular interest, since it is synergistic with the field-line-bending stabilisation of magnetohydrodynamic instabilities (Bernstein et al. 1958; Greene & Johnson 1962; Mercier 1962; Correa-Restrepo 1978; Connor et al. 1978). These expectations, scalings, synergies and behaviour alongside magnetohydrodynamic stability will be the subject of careful analysis and simulation in a future work (Rodríguez & Zocco 2024).