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 Reference Taylor1976), 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 Reference Taylor1976; Connor, Hastie & Taylor Reference Connor, Hastie and Taylor1978; Tang, Connor & Hastie Reference Tang, Connor and Hastie1980; Antonsen & Lane Reference Antonsen and Lane1980) 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. Reference D'haeseleer, Hitchon, Callen and Shohet2012), 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. (Reference Connor, Hastie and Taylor1980). 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 Reference Romanelli1989). 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.Footnote ¹

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$,Footnote ²

(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. Reference Connor, Hastie and Taylor1978). 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. Reference Terry, Anderson and Horton1982), 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.Footnote ³ 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 Reference Boozer1983; Nührenberg & Zille Reference Nührenberg and Zille1988; Rodríguez, Helander & Bhattacharjee Reference Rodríguez, Helander and Bhattacharjee2020) HSX (helically symmetric experiment) stellarator (Anderson et al. Reference Anderson, Almagri, Anderson, Matthews, Talmadge and Shohet1995).

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.Footnote ⁴ 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.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 (Reference Hahm and Tang1988); R.J. Hastie, private communication, 2013; Connor, Hastie & Zocco (Reference Connor, Hastie and Zocco2013); Zocco et al. (Reference Zocco, Plunk, Xanthopoulos and Helander2016)). 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 Reference Hahm and Tang1988; Connor et al. Reference Connor, Hastie and Zocco2013; Zocco et al. Reference Zocco, Plunk, Xanthopoulos and Helander2016). 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.Footnote ⁵ 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 Reference Fried and Conte2015). 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. Reference Kuroda, Sugama, Kanno, Okamoto and Horton1998; Sugama Reference Sugama1999) 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 Reference Fried and Conte2015) 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 Reference Woods1926, 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 Reference Abramowitz and Stegun1968, § 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 Reference Kadomtsev and Pogutse1970a) or local, short wavelength ITG (Smolyakov, Yagi & Kishimoto Reference Smolyakov, Yagi and Kishimoto2002) 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 Reference Johnson2024) 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.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. Reference Zocco, Plunk, Xanthopoulos and Helander2016). 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. Reference Hirose, Elia, Smolyakov and Yagi2002; Smolyakov et al. Reference Smolyakov, Yagi and Kishimoto2002; Gao et al. Reference Gao, Sanuki, Itoh and Dong2003).

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. Reference Zocco, Loureiro, Dickinson, Numata and Roach2015).

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 Reference Fried and Conte2015)

(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. (Reference Connor, Hastie and Taylor1980) (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 Reference Romanelli1989, equation (19)) is possible,Footnote ⁶

(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 Reference Wesson and Campbell2011; Zocco et al. Reference Zocco, Plunk, Xanthopoulos and Helander2016). Let the mode in which the streaming terms are negligibleFootnote ⁷ 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$,Footnote ⁸ 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. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018, Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022; Podavini et al. Reference Podavini, Zocco, García-Regaña, Barnes, Parra, Mishchenko and Helander2024). 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 Reference Jenko, Dorland and Hammett2001; Roberg-Clark, Plunk & Xanthopoulos Reference Roberg-Clark, Plunk and Xanthopoulos2022a).

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.Footnote ⁹ This type of ITG activity is referred to in the literature as SWITG (Hirose et al. Reference Hirose, Elia, Smolyakov and Yagi2002; Smolyakov et al. Reference Smolyakov, Yagi and Kishimoto2002; Gao et al. Reference Gao, Sanuki, Itoh and Dong2003).

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 Reference Abramowitz and Stegun1968, § 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 (Reference Romanelli1989) and Guo & Romanelli (Reference Guo and Romanelli1993), with a visual estimate obtained from figure 2 of Biglari, Diamond & Rosenbluth (Reference Biglari, Diamond and Rosenbluth1989), with the Jenko–Dorland–Hammett formula (Jenko et al. Reference Jenko, Dorland and Hammett2001) and other critical gradient estimates (Roberg-Clark et al. Reference Roberg-Clark, Plunk and Xanthopoulos2022a).

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. (Reference Biglari, Diamond and Rosenbluth1989) 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 Reference Hahm and Tang1988; Jenko et al. Reference Jenko, Dorland and Hammett2001; Roberg-Clark et al. Reference Roberg-Clark, Plunk and Xanthopoulos2022a), 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. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018, Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022).