Anisotropy and shaping effects on the stability boundaries of infernal ideal MHD modes in tokamak hybrid plasmas

The toroidally symmetric tokamak configuration is one of the most promising reactor types for achieving controlled thermonuclear fusion. An important figure of merit in fusion research is β, the ratio of kinetic over magnetic pressure, which is a measure of the plasma performance. Much effort has been aimed indeed at maximising this quantity. High plasma pressures are achieved by several heating schemes such as radio frequency heating or neutral beam injection (NBI). While both may produce parallel and perpendicular pressure anisotropy [1, 2], NBI induces strong toroidal flows [3]. Large β values are typical of spherical tokamaks (STs) (devices with an aspect ratio, i.e. the ratio of major over minor radii of the toroidal chamber, close to unity) which also exhibit strong natural shaping [4]. It has also been found that centrifugal effects have a significant impact on equilibrium and stability in STs [5]. In such devices, scenarios at high β feature optimised current profiles in which the core magnetic shear, a measure of the inclination of the magnetic field lines with respect to one another, is either reversed (advanced scenario), or small over a wide region (hybrid scenario) with the safety factor (the pitch of the magnetic field lines denoted by q) always above unity [6].

Under these circumstances, long-lived saturated magnetohydrodynamic (MHD) activity, usually with low n numbers, often occur [6, 7]. Indeed, scenarios with broad current and peaked pressure profiles are prone to develop infernal type instabilities [8–10], driven when the safety factor is close to a rational over a wide region. These instabilities exhibit a dominant mode of helicity m/n which is nearly resonating with a low q rational over a broad region. We refer to such modes as nearly-resonant, since the exact resonance of the mode with q in this broad region is not a necessary requirement for the instability to develop. We point out that spherical tokamaks are more susceptible to these ideal internal modes [6] because of the tighter aspect ratio and higher β which yield stronger toroidicity induced couplings that enhances the infernal driving mechanism. Nevertheless, such kinds of instabilities have been numerically and experimentally reported in standard tokamaks such as TFTR and DIII-D [7, 11–14]. We additionally note that reversed shear configurations with internal transport barriers (ITBs) can be unstable against pressure driven modes with infernal-like features [15].

As discussed in references [7, 11–14, 16], such instabilities are usually associated with β limits (both soft or hard). Moreover, as these modes tend to saturate in amplitude (often preserving the linear mode structure [6, 17]), fast ion losses can be enhanced [6]. Since preventing the access to the linear phase avoids the development of a saturated non-linear state [17–19], understanding the linear stability properties of such perturbations is of crucial interest.

The impact of plasma shaping in isotropic static plasmas has been extensively analysed [20–28]. A vast literature also exists on the modifications to equilibrium and stability due to toroidal flows [29–32] and pressure anisotropy (see e.g. references [33–41]). Although these two effects were generally treated separately, more recent analytical and numerical investigations provided a unified framework [2, 42–49]. The reader interested in the experimental impact of plasma anisotropy on equilibrium and stability is referred to references [41, 45, 48] and references therein. The current paper concentrates specifically on the stability analysis of ideal infernal modes in shaped tokamaks characterised by strong degrees of pressure anisotropy and toroidal flows. The problem of plasma anisotropy is tackled within the guiding centre plasma model [43, 50], in which the macroscopic plasma motion across the magnetic field is fluid-like (i.e. identical to MHD), while the dynamics parallel to the magnetic field is described by a collisionless kinetic equation. Although the main aim is to describe the dynamics of compact configurations, a large aspect ratio expansion may still be employed as long as both the equilibrium and stability analysis is carried out sufficiently near to the magnetic axis [51]. Hence, the relevance of the analysis presented in this work can be naturally extended to standard tokamaks as well.

In analogy with previous results reported in the literature [26, 31, 39], it is found that plasma shaping (elongation and triangularity), anisotropy and equilibrium flows modify the magnetic well (or hill). Contributions due to a flat toroidal rotation with a monotonically decreasing pressure profile tend to increase the magnetic well implying stabilisation. It is noted that combined guiding centre-MHD model employed here is essentially isothermal, and as such differs from the way in which flows are known to modify MHD instabilities (see e.g. [52]). Pressure anisotropy, beyond trivially affecting the averaged total pressure, has a (de)stabilising effect when $T_{||}(\lt)\gt T_\perp$ [37, 39, 40, 45]. Vertical elongation decreases the magnetic well, whose depth can be restored by allowing for corrections due to positive triangularity. The same effect is achieved in negative triangularity plasmas with an oblate cross section. Finally, elongation is found to alter the strength of the coupling between neighbouring modes, improving the global stability in vertically elongated plasmas.

The paper is organised as follows. In section 2, the anisotropic single fluid MHD model is summarised and subsequently employed for the derivation of the equilibrium equations for non-circular axisymmetric toroidal equilibria, which is outlined in section 3. Section 4 is devoted to the description of the perturbed dynamics within the infernal framework. The governing equations for the mutually coupled harmonics are derived, distinguishing between regions of high and low magnetic shear. In sections 5 and 6, by solving for the perturbed system previously derived, stability boundaries are identified and analysed by exploring several plasma parameters (e.g. degree of anisotropy, shaping) for monotonic and reversed q configurations. A discussion on the features of the physical model and the results with their implications on present and future experiments is given in sections 7 and 8 respectively. Finally, appendix A presents a brief analysis for the allowance for resistive effects on the satellite harmonics.

Under the assumption of strong heating, we regard the plasma as ideally conducting with vanishing resistivity. Thus, the plasma dynamics are described by the single fluid anisotropic ideal MHD equations [43]:

$$\begin{align}\partial_t\rho+ {\boldsymbol{\nabla}}\cdot(\rho{\boldsymbol u}) = 0, \end{align} \tag{ 1 }$$

$$\begin{align}\partial_t{\boldsymbol B} = {\boldsymbol \nabla}\times({\boldsymbol u}\times{\boldsymbol B}), \end{align} \tag{ 2 }$$

$$\begin{align}\rho\left(\partial_t{\boldsymbol u}+{\boldsymbol u}\cdot{\boldsymbol{\nabla}}{\boldsymbol u} \right) = -{\boldsymbol{\nabla}}\cdot \underline{\underline{\boldsymbol P}} +{\boldsymbol J}\times{\boldsymbol B}, \end{align} \tag{ 3 }$$

where ρ is the mass density, ${\boldsymbol u}$ the plasma fluid velocity, ${\boldsymbol B}$ the magnetic field with $|{\boldsymbol B}| = B$, ${\boldsymbol J} = {\boldsymbol{\nabla}}\times{\boldsymbol B}$ and $\underline{\underline{\boldsymbol P}} = p_\perp\underline{\underline{\boldsymbol I}}+(p_{||}-p_\perp){\boldsymbol b}{\boldsymbol b}$ with $\underline{\underline{\boldsymbol I}}$ the diagonal unit tensor and ${\boldsymbol b} = {\boldsymbol B}/B$. The parallel and perpendicular pressure are defined as moment averages in guiding centre coordinates as [43, 44, 53]

$$\begin{align*} (p_{||},p_{\perp}) = &\sum_s m_s\int d^3{\boldsymbol v}f_s[(v_{||}-{\boldsymbol b}\cdot{\boldsymbol u})^2,v_\perp^2/2], \end{align*}$$

where f_s is the particle distribution function of the species s of mass m_s (s = i, e for ions and electrons respectively), ${\boldsymbol v}$ is the microscopic particle velocity (with parallel and perpendicular projections wrt the magnetic field indicated by $v_{||}$ and $v_\perp$ respectively), and the sum is extended over all species. Hereafter for a generic vector quantity ${\boldsymbol A}$ we indicate $A = |{\boldsymbol A}|$.

For the computation of $p_{\perp,||}$ knowledge of the distribution function for each plasma species is required. In guiding centre theory, this satisfies the drift-kinetic equation [43, 44]

$$\begin{align} \nonumber\partial_tf_s+&({\boldsymbol u}_\perp+v_{||}{\boldsymbol b})\cdot{\boldsymbol{\nabla}}f_s+ \end{align}$$

$$\begin{align}&[v_{||} {\boldsymbol u}_\perp\cdot({\boldsymbol b}\cdot{\boldsymbol{\nabla}}{\boldsymbol b})-{\boldsymbol b}\cdot{\boldsymbol{\nabla}}E_s]\frac{\partial f_s}{\partial v_{||}} = 0, \end{align} \tag{ 4 }$$

with $E_s = \mu B+\frac{e_s}{m_s}\Phi-\frac{u^2_\perp}{2}$ where $\mu = v_\perp^2/2B$ is the particle magnetic moment, e_s being the particle charge and the parallel electric field given by $E_{||} = -{\boldsymbol b}\cdot{\boldsymbol{\nabla}}\Phi$. Note that corrections due to collisions have been dropped. We point out that these results are unaffected by gauge transformations of the form Φ→Φ+h(r). Although such a function, namely h, does not play a role in the equilibrium calculations, it might have an effect on perturbed quantities. Its form will then be determined in section 4. For a plasma with an equilibrium flow ${\boldsymbol U} = {\boldsymbol u}_0$ (the subscript 0 indicates the corresponding equilibrium quantity), we introduce the variable $\epsilon_s = \frac{1}{2}v_{||}^2-v_{||}U_{||}+E_{s0}$ [44] with $U_{||} = {\boldsymbol b}_0\cdot{\boldsymbol U}$, so that the parallel and perpendicular pressure are given by

$$\begin{equation} (p_{||},p_{\perp}) = \sum_{s,\sigma}2\pi m_s\int_0^\infty d\mu\int_{\epsilon_m}^\infty d\epsilon_s \frac{Bf_s}{|\bar{v}_{||}|}[(v_{||}-{\boldsymbol b}\cdot{\boldsymbol U})^2,\mu B] \end{equation} \tag{ 5 }$$

where $\bar{v}_{||} = v_{||}-U_{||}$, $\sigma = sign(\bar{v}_{||})$ and $\epsilon_m = \mu B_0+\frac{e_s}{m_s}\Phi_0-U^2/2$. Finally, the density of each species is

$$\begin{equation} n_s = \sum_{\sigma}2\pi\int_0^\infty d\mu\int_{\epsilon_m}^\infty d\epsilon_s \frac{Bf_s}{|\bar{v}_{||}|}, \end{equation} \tag{ 6 }$$

and $\rho = \sum_s m_sn_s\approx m_in_i$ with $n_i = Z_in_e$ (quasineutrality for single ion species). Hereafter we take Z_i  = 1. The framework in which the following analysis is performed, is completely determined by equations (1)–(6).

We analyse a tokamak configuration of major and minor radii R₀ and a respectively, with non-circular shifted toroidal surfaces. The plasma is assumed surrounded by a perfectly conducting metallic wall. Let us introduce the coordinate system (r, θ, φ) where r is a flux labelling variable with the dimensions of a length, while θ and φ are the poloidal and toroidal angles. The covariant and contravariant components of the vector ${\boldsymbol A}$ are indicated by A_i and Aⁱ respectively. For sake of clarity, hereafter perturbed quantities will be denoted by a tilde and we will drop the subscript 0 for the equilibrium ones. We use the notation (·)⁽⁰⁾ for the leading order of $\int_0^{2\pi}(\cdot)d\theta/2\pi\equiv\langle\cdot\rangle$. The equilibrium magnetic field is written as ${\boldsymbol B} = B_\phi{\boldsymbol{\nabla}}\phi-{\boldsymbol{\nabla}}\psi\times{\boldsymbol{\nabla}}\phi$ (the magnetic field strength at the magnetic axis is denoted by B₀) with the safety factor function and magnetic shear defined by $q = \langle B^\phi/B^\theta\rangle$ and $\hat{s} = rd\text{ln}\,q/dr$. Here $B^\theta = \frac{1}{\sqrt{g}}\frac{d\psi}{dr}$ where $\sqrt{g}$ is the Jacobian associated with the coordinate system (r, θ, φ). We assume that $\hat{s}\approx0$ for $r_1\lt r\lt r_2$, while $\hat{s}\gt 0$ for 0 < r < r₁ and r₂ < r < a.

From (4), at equilibrium $f_s = f_s(r,\mu,\epsilon_s)$ [35, 43, 44]. Noting the necessary independence of f_s on θ, we choose a 'modified' bi-Maxwellian equilibrium distribution function with equal temperatures for both particle species (i.e. ions and electrons) [37, 43]

$$\begin{equation*} f_s = \frac{n_s}{(2\pi/m_s)^{3/2}T_\perp\sqrt{T_{||}}}\text{exp}\,[-\frac{m_s}{2}(\frac{\bar{v}_{||}^2}{T_{||}} +\frac{v_{\perp}^2}{T_{\perp}})]. \end{equation*}$$

Note that n_s is allowed to depend upon θ also. We shall point out that, depending on the physics that is analysed, other choices for the equilibrium distribution function are possible [1, 45]. The distribution function above allows for analytically tractable solutions to the equilibrium and stability problem while also accounting for the features of anisotropy (modification of the magnetic well, average curvature, etc). It is noted that the distribution is a limiting form of a model widely used for the anisotropic features of ICRH (in references [54, 55] with B_c chosen specifically as $B_{\textrm{max}}\equiv B_0(1+r/R_0)$). By taking the parallel gradient of equations (5) and (6) and imposing quasineutrality $\sum_se_sn_s = 0$, we have [39, 43, 44]

$$\begin{align*} {\boldsymbol{\nabla}}_{||}p_{\perp}& = \sigma_{\perp}{\boldsymbol{\nabla}}_{||}\text{ln}\,B+\frac{T_\perp}{T_{||}}\rho{\boldsymbol{\nabla}}_{||}U^2/2,\\ {\boldsymbol{\nabla}}_{||}p_{||}& = \sigma_{||}{\boldsymbol{\nabla}}_{||}\text{ln}\,B+\rho{\boldsymbol{\nabla}}_{||}U^2/2,\\ {\boldsymbol{\nabla}}_{||}\text{ln}\,\rho& = (1-\frac{T_\perp}{T_{||}}){\boldsymbol{\nabla}}_{||}\text{ln}\,B+\frac{m_i}{2T_{||}}{\boldsymbol{\nabla}}_{||}U^2/2, \end{align*}$$

where ${\boldsymbol{\nabla}}_{||} = {\boldsymbol b}\cdot{\boldsymbol{\nabla}}$, $\sigma_{\perp} = 2p_\perp(1-\frac{T_\perp}{T_{||}})$, $\sigma_{||} = p_{||}-p_\perp$.

For an equilibrium velocity of the form ${\boldsymbol U} = R^2\Omega(r){\boldsymbol{\nabla}}\phi$, the system above is solved by [43]

$$\begin{align*} \rho = \bar{\rho}(r)\frac{T_\perp}{T_{||}}\text{exp}\,[\mathcal{M}^2(R^2/R_0^2-1)],\;\;\;T_\perp = \frac{BT_{||}}{B-\Theta(r)}, \end{align*}$$

where $T_{||} = T_{||}(r)$ with $p_{\perp,||} = 2T_{\perp,||}\rho/m_i$ and $\mathcal{M}^2 = \rho R_0^2\Omega^2/(2p_{||})$. Here the function Θ(r) measures the degree of anisotropy, i.e. $\Theta(r) = B(T_\perp-T_{||})/T_\perp$. For sake of simplicity we take $\bar{\rho}$, Ω and Θ constant [31, 43] and normalise B₀ = 1. Finally, from (3) it is found that $B_\phi = F(r)/(1-\sigma_{||}/B^2)$.

We employ the large aspect ratio approximation ($\varepsilon = a/R_0\ll1$) which has been proven, via comparison with codes, to give reliable, at least qualitatively, results also for compact configurations [17, 51] providing that the analysis is performed sufficiently close to the magnetic axis. Because the analysis will then focus on radially extended perturbations, we extend the local equilibrium model presented in reference [56] by allowing elongation and triangularity corrections to be dependent upon the radial variable. By taking the covariant radial projection of (3) with $\partial_t\to 0$ and assuming $p_{\perp,||}\sim\varepsilon^2$ with a sufficiently small magnetic shear [20], a tokamak shaped equilibrium is solved to leading order by

$$\begin{equation} R = R_0+r\text{cos}\,(\theta+r\frac{\delta}{a}\text{sin}\,\theta)-\Delta,\;\;\;Z = \kappa r\text{sin}\,\theta. \end{equation} \tag{ 7 }$$

Here κ ∼ 1 and δ ∼ ɛ are real numbers describing plasma elongation and triangularity respectively, with the Shafranov shift Δ∼ɛa fulfilling ($'\equiv d/dr$)

$$\begin{equation} \Delta^{\prime\prime}+\frac{3}{r}\Delta^{\prime} +\frac{4}{1+3\kappa^2}\left( \frac{2q^2R_0\bar{p}^{\prime}}{r} +\frac{2\delta}{a}-\frac{\kappa^2}{R_0} \right) = 0 \end{equation} \tag{ 8 }$$

with $\bar{p} = (p_{||}^{(0)}+p_\perp^{(0)})/2$. Finally, $\psi' = r\kappa/q$ and $F = R_0(1-\frac{r^2(1+\kappa^2)}{2q^2R_0^2}-p_{||}^{(0)})$. Equation (8) is valid from $T_{||}/T_\perp\gg1$ to $T_{||}\sim T_\perp$. However for $T_\perp/T_{||}\sim\varepsilon^{-1}$ additional harmonics in $p_{||,\perp}$ with respect to θ are generated [41, 57]. In regions with $\hat{s}\sim1$, where it is not necessary to resolve the equilibrium to such an order, we take R = R₀+rcosθ and Z = κrsinθ. We point out that for describing equilibria which have higher β values and extremely pronounced shaping, different techniques must be employed. Nevertheless, the analysis presented above proves to be adequate in dealing with most of the experimentally relevant configurations.

We transform θ → ϑ, where ϑ is the rectified angle for which the field lines are straight [58] and the ratio $B^\phi/B^\vartheta$ is a function of r only. It is easily verified that θ = ϑ + λ(r, ϑ) with $\lambda = -(r/R_0+\Delta')\text{sin}\,\vartheta$. Hence, by means of (7), the metric tensor elements $g_{ij} = \partial_iR\partial_jR+\partial_iZ\partial_jZ$ in the coordinate system (r, ϑ, φ) can be straightforwardly derived to order ɛ. This is the required accuracy needed for the stability calculations outlined in the next section.

In order to investigate the stability of the system, we employ the energy method [44, 59]. Introducing the Lagrangian displacement ${\boldsymbol\eta}$ [60], according to reference [44] the momentum equation (3) is written as $\rho(\partial_t+{\boldsymbol U}\cdot{\boldsymbol{\nabla}})^2{\boldsymbol\eta} = {\boldsymbol F}({\boldsymbol\eta})$ where ${\boldsymbol F}(\boldsymbol\eta) = \tilde{\boldsymbol J}\times{\boldsymbol B}+{\boldsymbol J}\times\tilde{\boldsymbol B}-{\boldsymbol{\nabla}\times\tilde{\underline{\underline{\boldsymbol P}}}}+ {\boldsymbol{\nabla}}\cdot[\rho{\boldsymbol\eta}{\boldsymbol U}\cdot{\boldsymbol{\nabla}}{\boldsymbol U}]$ is self-adjoint (see also [50, 61] and references therein). The equation above can be cast in the form [60]

$$\begin{equation} \rho(\partial_t^2{\boldsymbol\eta}+2{\boldsymbol U}\cdot{\boldsymbol{\nabla}}\partial_t{\boldsymbol\eta}) = {\boldsymbol F}({\boldsymbol\eta})- {\boldsymbol{\nabla}}\cdot(\rho{\boldsymbol U}{\boldsymbol U}\cdot{\boldsymbol{\nabla}}{\boldsymbol\eta}), \end{equation} \tag{ 9 }$$

where it can be easily shown that the last term on the rhs is self-adjoint. It follows that a necessary and sufficient criterion for stability can be derived [60, 62]. Note that $({\boldsymbol U}\cdot{\boldsymbol{\nabla}}{\boldsymbol\eta})_i = \Omega[\partial_\phi\eta_i-({\boldsymbol z}\times{\boldsymbol\eta})_i]$ [32] with $\sqrt{g}z^r = -R\partial_\vartheta R$, $\sqrt{g}z^\vartheta = RR'$ and $z^\phi = 0$. We express the perturbation in the form $\eta^j = \sum_{m,n}\eta^j_{m,n}e^{i(m\vartheta-n\phi)+(\gamma_r+i\omega)t}+c.c.$ (j = r, ϑ, φ) where γ_r and ω are real and c.c. stands for complex conjugate. Thus, we dot (9) with ${\boldsymbol\eta}$ and divide it by $e^{2\gamma_r t}$, with ω = nΩ where Ω is constant. Integrating the resulting equation over the plasma volume and averaging over an oscillation period 2π/ω, yields an equation for $\gamma_r^2$ showing that this quantity is real, which therefore indicates that unstable modes rotate with frequency nΩ [31]. It follows that the stability boundaries are identified by γ_r  = 0.

Having established the mode frequency characteristics, we shall now proceed with the derivation of the stability equations. Assuming that perturbed quantities have a time dependence of the form exp(γt) and vary along the poloidal and toroidal angles proportionally to $\text{exp}\,(i\ell\vartheta-in\phi)$, we introduce the perturbed fluid displacement ${\boldsymbol\xi} = \tilde{\boldsymbol u}/(\gamma-in\Omega)$. Since we assume Ω constant, it follows that γ − iΩ is constant as well, and hence is not affected by the differential operators. Let us denote $(\cdot)_\ell = \int_0^{2\pi}(\cdot)e^{-i\ell\vartheta+in\phi}d\vartheta d\phi/4\pi^2$. Within the infernal framework, we impose a wide region of flat q for $r_1\lt r\lt r_2$ (we shall specify r₁ later) which nearly resonates with a dominant mode of helicity m/n accompanied by its neighbouring sidebands with poloidal mode numbers m ± 1. Hereafter, m will always denote the poloidal mode number of the dominant harmonic. Because of the coupling induced by elongation, a larger number of harmonics, namely m ± 2, m ± 3, ... all of order ɛ, should be in principle taken into account [24, 27]. However, it will be shown later that under appropriate conditions, retaining all these harmonics is not necessary and the system can be described by linear coupling of three modes. Nevertheless, for illustrative purposes, we formally allow for higher order harmonics. We write $q = m/n+\delta q$ and we adopt the ordering m ∼ n ∼ q ∼ 1 with δq ∼ ɛ, ${\boldsymbol\xi}_{m\pm\ell}\sim\varepsilon{\boldsymbol\xi}_m$ $\ell = 1,2,\ldots$, and γ ∼ Ω∼ɛω_A with $\omega_A = 1/\sqrt{R_0\rho|_{r = 0}}$. Therefore, the contravariant radial and poloidal projections of the perturbed (2) yield respectively at leading order $(\sqrt{g}\tilde{B}^r)_\ell = i\kappa r(\ell/q-n)\xi^r_\ell$ and $\frac{1}{r}(r\xi_m^r)^{\prime}+im\xi_m^\vartheta-in\xi_m^\phi = 0$. Hereafter we indicate with $\sqrt{g}$ the Jacobian associated with the coordinate system (r, ϑ, φ). Since Ω is constant and considering $\sqrt{g}\tilde{B}^\phi$ small enough, which is verified a posteriori, from the contravariant φ projection of (2) it follows that $\xi_m^\phi\approx0$. Note that although sideband harmonics have ɛ times smaller fluid displacements compared to the dominant mode, their associated magnetic perturbations are of the same order.

In the covariant basis identified by vectors ${\boldsymbol e}_{r,\vartheta,\phi}$, we shall note that ${\boldsymbol U}_\perp = {\boldsymbol b}\times({\boldsymbol U}\times{\boldsymbol b})\approx(0,-\Omega/q,0)$, and $\widetilde{(\boldsymbol u_\perp)}\approx(\tilde{u}^r-\Omega\tilde{B}^r/B^\phi,\tilde{u}^\vartheta-\tilde{u}^\phi/q-\Omega\tilde{B}^\vartheta/B^\phi,0)$. By taking the covariant φ projection of the linearised equation (3), it can be shown that at leading order

$$\begin{equation} \tilde{B}_\phi\approx-R_0\tilde{p}_\perp. \end{equation} \tag{ 10 }$$

An expression for the perturbed distribution function is required in order to obtain the fluctuation of the mass density and the parallel/perpendicular pressure. We turn to (4) and we change variables, replacing the variable $\bar{v}_{||}$ with _s (this turns out to be more convenient when working in toroidal geometry [44, 53]). By employing the expression for the perturbed velocity given above, after some algebra it can be shown that

$$\begin{align} \nonumber[\hat{\gamma}-&(v_{||}-U_{||}){\boldsymbol b}\cdot{\boldsymbol{\nabla}}](\tilde{f}_s+{\boldsymbol\xi}\cdot{\boldsymbol{\nabla}}f_s)-(v_{||}-U_{||})\times \end{align}$$

$$\begin{align}&\times[{\boldsymbol b}\cdot{\boldsymbol{\nabla}}\tilde{E}_s+{\boldsymbol u}_\star\cdot{\boldsymbol{\nabla}}(v_{||}+R/R_0)]\frac{\partial f_s}{\partial\epsilon_s} = 0, \end{align} \tag{ 11 }$$

where $\hat{\gamma} = \gamma-in\Omega$, $\tilde{E}_s = \mu\tilde{B}+\frac{e_s}{m_s}\tilde{\Phi}-{\boldsymbol U}_\perp\cdot\widetilde{({\boldsymbol u}_\perp)}$ and ${\boldsymbol u}_\star = (\tilde{u}^r,\tilde{u}^\vartheta-\tilde{u}^\phi/q,0)$ in the basis ${\boldsymbol e}_i$, having approximated ${\boldsymbol b}\cdot{\boldsymbol{\nabla}}{\boldsymbol b}\approx-{\boldsymbol{\nabla}}_\perp R/R_0$. For the calculation of the stability boundaries, we let $\hat{\gamma}\to 0$ so that terms involving ${\boldsymbol u}_\star\propto\hat{\gamma}{\boldsymbol\xi}$ vanish. Let us first note that in choosing the gauge function h defined in section 2 to be vanishing, we ensure that $\langle\tilde{\Phi}\rangle = 0$. The term proportional to $\tilde{E}_s$ in equation (11) produces small corrections to $\tilde{p}_{||,\perp}$ and n_s , hence it can be safely neglected. Thus, the perturbed distribution function can be written as [39, 40, 53]

$$\begin{equation*} \tilde{f}_s = -{\boldsymbol\xi}\cdot{\boldsymbol{\nabla}}f_s, \end{equation*}$$

where we dropped trapped particles contributions [53]. Note that neglecting compressibility effects is allowed only at marginal stability. By using the equation above for the evaluation of $\tilde{p}_{||,\perp}$ and $\tilde{\rho}$ (obtained at leading order from (5)–(6) with the substitution $f_s\to\tilde{f}_s$) we get [53]

$$\begin{align} &\tilde{p}_{||,\perp} = -\xi^r[(p_{||,\perp}^{(0)})^{\prime}-\sigma_{||,\perp}^{(0)}\frac{r}{R_0}\text{cos}\,\vartheta], \end{align} \tag{ 12 }$$

$$\begin{align} &\tilde{\rho} = -2\xi^r\rho^{(0)}(\mathcal{M}^2)^{\prime}\frac{r}{R_0}\text{cos}\,\vartheta, \end{align} \tag{ 13 }$$

where small corrections due to $\tilde{B}$ have been dropped by virtue of (10). We recall that $\bar{\rho}$ and Θ are constant and so is ρ⁽⁰⁾. Note that $\tilde{p}_{||,\perp}$ satisfies to leading order the parallel projection of (3), i.e.

$$\begin{equation*} 0\approx{\boldsymbol B}\cdot{\boldsymbol{\nabla}}\tilde{p}_{||}+\tilde{\boldsymbol B}\cdot{\boldsymbol{\nabla}}p_{||}+ \psi'\tilde{\sigma}_{||}\partial_\vartheta R/R_0. \end{equation*}$$

All the expressions for the perturbed quantities entering the linearised equations (1)-(3) have been obtained.

We now apply the operator ${\bf{D}}\equiv\sqrt{g}{\boldsymbol{\nabla}}\phi\cdot{\boldsymbol{\nabla}}\times(1/B^\phi)$ to the linearised (3) [53, 63]. By taking the m and m ± 1 moments of the resulting expression, three equations for the corresponding Fourier modes are obtained. Employing the usual techniques for infernal modes [10], we distinguish between low and high shear regions.

4.1. Low shear region

Let us first introduce the elongation factor, or ellipticity, $e = (\kappa^2-1)/(\kappa^2+1)$. It is worth noting that with an elongation of order unity, a mode $\bar{m}$ couples to all harmonics $\bar{m}\pm2\ell$ ($\ell = 1,2,\ldots$) due to elongation induced metric oscillations. An analytic treatment which retains the whole spectrum is clearly hopeless. In order to make the algebra analytically tractable, we follow reference [26] in which a double expansion, first in ɛ and then independently to first order in e, is performed (such an approximation proves to hold for $e\lesssim\frac{1}{2}$ [26] or $\kappa\lesssim2$). It can be shown that harmonics with $\ell\geq(\leq) m+(-)2$ can be neglected, as they enter the analysis to second order in e.

Focussing the analysis in the low shear region, we note that within the above mentioned approximations, shaping and flow effects are contained in the field line bending and perturbed toroidal field contributions. The latter, along with the perturbed pressure tensor, generates additional anisotropy corrections. After some algebraic manipulations, it is possible to show that

$$\begin{align}{\bf{D}}({\boldsymbol{\nabla}}\cdot\underline{\underline{\tilde{\boldsymbol P}}})\approx r\text{sin}\,\vartheta(\tilde{p}_{||}+\tilde{p}_\perp)'+\text{cos}\,\vartheta\partial_\vartheta(\tilde{p}_{||}+\tilde{p}_\perp) \end{align}$$

$$\begin{align}&-\partial_\vartheta[\{R_0(p_\perp^{(0)})'+\langle \frac{1}{B^\phi}\frac{d\text{ln}\,B}{dr}\rangle\} \tilde{\sigma}_{||} -\langle1/B^\phi\rangle'\tilde{p}_\perp]. \end{align} \tag{ 14 }$$

The field line bending term contains coupling with the $\ell\pm2$ sidebands due to the elongation induced metric tensor oscillations. It can be shown that at leading order the equations for the sideband harmonics are [24] (for sake of clarity we drop the superscript r in denoting the radial fluid displacement, viz. $\xi^r_\ell\to\xi_\ell$)

$$\begin{align*} [r^{-1\mp2m}(r^{2\pm m}\xi_{m\pm1})']'-&\frac{e}{2}(\frac{m\pm1}{m\mp1})[\frac{1}{r}(r^{2\mp m}\xi_{m\mp1})']' = \\ &\frac{1\pm m}{1+\kappa^2}[r^{\mp m}\alpha\xi_m]', \end{align*}$$

with $\alpha = -2R_0\bar{p}'q^2$. These two equations can be integrated yielding to leading order in e

$$\begin{align} \nonumber r^{-1\mp2m}(r^{2\pm m}\xi_{m\pm1})'& = \frac{(1\pm m)}{1+\kappa^2}(1-\frac{e}{2})r^{\mp m}\alpha\xi_m+ \end{align}$$

$$\begin{align}&L_\pm -\frac{e}{2}(\frac{1\pm m}{1\mp m})r^{\mp2m}L_\mp, \end{align} \tag{ 15 }$$

where $L_\pm$ are constants, which in general depend on e, to be determined later through matching with the sheared region solutions [10]. If the pressure profile has a step at $\bar{r}$ so that $\alpha\propto\delta(r-\bar{r})$, expanding $\xi_{m\pm1}$ and $L_\pm$ in e and integrating across $\bar{r}$ shows that the constants $L_\pm$ are the same on either side of $\bar{r}$ [64] (this will turn to be useful in the next sections).

Let us introduce the Newcomb's operator [65]

$$\begin{equation*} \mathcal{L}_\ell = \frac{1}{r}\frac{d}{dr}\left[r^3(\frac{\ell}{q}-n)^2\frac{d}{dr}\right]-(\ell^2-1)(\frac{\ell}{q}-n)^2. \end{equation*}$$

Elongation driven coupling can be neglected in the analysis of the mode m as it enters the equations proportionally to e². Hence, the mth component of the perturbed field line bending term can be written as $i(\sqrt{g}{\boldsymbol B}\cdot{\nabla}\tilde{J}^\phi/B^\phi)_m\approx\frac{1+\kappa^2}{2 {\rm m} R_0} \mathcal{L}_m(\xi_m)$. Thus, employing standard techniques and taking the limit $\hat{\gamma}\to 0$, by means of (10), (12), (13) and (14) after rather lengthy algebra we obtain

$$\begin{align} \nonumber[r^3k_{||}^2\xi_m']'-r\{ (m^2&-1)k_{||}^2+rw'\}\xi_m+ \end{align}$$

$$\begin{align}&\frac{\alpha}{2}(1-\frac{3e}{2})\sum_\pm\frac{r^{1\pm m}L_\pm}{1\pm m} = 0, \end{align} \tag{ 16 }$$

where $k_{||} = 1-\frac{n}{m}q$ and the function $w'$ associated with the plasma magnetic well is given by

$$\begin{equation} w' = \frac{\alpha}{R_0}[1-\frac{1}{q^2}+(1-e)\tau -\frac{3e}{4}+\frac{3}{2}\frac{e\delta}{\varepsilon}]+(1-e)(\frac{q\Omega}{\omega_A})^2\frac{d\mathcal{M}^2}{dr}, \end{equation} \tag{ 17 }$$

with $\tau = (\frac{1}{2}+\frac{T_\perp}{T_{||}})\frac{T_{||}-T_\perp}{T_{||}+T_\perp}$ [39] where it is understood that $T_{||,\perp}$ is taken on the magnetic axis. Note that, regardless of either plasma shaping or pressure anisotropy, the term associated with plasma rotation is positive for monotonically decreasing pressure profiles. This indicates a deepening of the magnetic well, and therefore a stabilising effect [31]. The dynamics in the low shear region is completely determined by equations (15) and (16). The derivation of the harmonics governing equations in the high shear region is the aim of the next subsection.

4.2. Sheared region

In the regions of large magnetic shear, i.e. for 0 < r < r₁ and r₂ < r < a, the parallel wave vector associated with the dominant mth mode is large enough so that inertial (viz. the lhs of (3)) and coupling terms can be neglected. We recall that m denotes always the mode number of the main harmonic. Thus, to leading order (i.e. to e²) the fluid disturbance of helicity m/n obeys the Newcomb's equation [10]

$$\begin{equation} \mathcal{L}_m(\xi_m) = 0. \end{equation} \tag{ 18 }$$

Multiplying the equation above by ξ_m and integrating from r₂ to a, yields ξ_m  = 0. In the region 0 < r < r₁, the same procedure gives ξ_m  = 0 for m > 1 and ξ_m  = const for m = 1.

Thus, since the main harmonic vanishes for r₂ < r < a it follows that in this region the sidebands behave according to

$$\begin{equation} [\sqrt{g}{\boldsymbol B}\cdot{\boldsymbol{\nabla}}(\tilde{J}^\phi/B^\phi)+ \sqrt{g}\tilde{\boldsymbol B}\cdot{\boldsymbol{\nabla}}(J^\phi/B^\phi) ]_{m\pm1} = 0. \end{equation} \tag{ 19 }$$

Note that elongation driven coupling between satellite harmonics of the type $m\pm1\to m\mp1$, which are of the same order, are allowed regardless of the weakness (or strength) of the magnetic shear. For monotonic q profiles, it is sufficient for to require $\xi_m(r_2) = 0$ and smooth matching of the sidebands across r₂ for any m ≥ 1. Let us consider inverted safety factor profiles with r₁≠ 0. Focussing on the internal high shear region 0 < r < r₁ we distinguish between m > 1 and m = 1 dominant modes. For m > 1, the same logic adopted above implies that the satellite harmonics fulfil (19) in this region also. Thus, the boundary conditions at r₁ are $\xi_m(r_1) = 0$ and smooth sidebands at this point.

The analysis of the m = 1 mode is more complicated. Since ξ_m does not vanish for 0 < r < r₁, additional terms must be retained in the sideband equations. Moreover, a more careful computation of (18) (up to order ɛ² even in low shear circular plasmas [44, 53, 66, 67]) is required to obtain the correct boundary condition (viz. $\xi_m'(r_1)$) at r₁ [67]. We point out that infernal modes with m = 1 occur when q is very close to unity over a broad region. Usually, in inverted q configurations such a region is not particularly wide. This induces us to conjecture that the dynamics of the m = 1 mode with an inverted q profile are more kink-like than infernal-like, and, as such, are better described by computing δW in the region 0 < r < r₁ [68] and matching the solution across r₁ and r₂ allowing for second variations in q, i.e. $q''$ [10]. Nonetheless, we may argue that the m = 1 mode is not strictly relevant for sufficiently inverted q scenarios with the minimum of the safety factor well above unity [7, 11, 69–74]. Hence, with a reversed magnetic shear we consider m > 1 dominant modes only, i.e. infernal modes with $q_{\textrm{min}}>1$.

The mode stability is determined by equations (16) and (19) with the boundary conditions given above. This will be analysed in the next section.

Let us consider a broad region of weak magnetic shear extending from the magnetic axis to r₂ (i.e. we let r₁→ 0). We denote the value of q in this region with q_m . For r₂ < r < a we take $q = q_m(r/r_2)^2$ so that at leading order the flux surface averaged toroidal current density is vanishing. By imposing the mode m + 1 having a resonance within the plasma, the maximum allowed extension of the current channel is $r_2/a<\sqrt{m/(1+m)}$ with q_m ≈m/n.

Because of the presence of a perfectly conducting metallic wall at the plasma edge, we have $\xi_{m,m\pm1}(a) = 0$. As discussed in section 4.2, the appropriate boundary condition for ξ_m at r₂ is $\xi_m(r_2) = 0$ (we recall that this expression is exact to order e²). We cast (19) in the form

$$\begin{equation} \mathcal{L}_{m\pm1}(\xi_{m\pm1}) = \mathcal{N}_{m\mp1}(\xi_{m\mp1}), \end{equation} \tag{ 20 }$$

where $\mathcal{N}_{m\pm1}\sim e$ is a linear functional. After setting L₋ = 0 for regularity of the lower sideband on the magnetic axis [75, 76] (this will be proven also in the next section where a more general case is addressed), matching ξ_m − 1 smoothly (to leading order in e) at r₂ yields ξ_m − 1≈ 0. Hence, by joining ξ_m + 1 across r₂ we have [10, 75]

$$\begin{equation*} \frac{L_+}{1+m} = \frac{(1+m)(1-3e/2)}{r_2^{2+2 {\rm m}}} \left(\frac{2+m+\mathcal{B}_+}{m-\mathcal{B}_+}\right)\int_0^{r_2}r^{1+m}\alpha\xi_mdr, \end{equation*}$$

where $\mathcal{B}_+ = r_2\xi_{m+1}'(r_2)/\xi_{m+1}(r_2)$ is obtained through solving $\mathcal{L}_{m\pm1}(\xi_{m\pm1})\approx0$. It easily follows that

$$\begin{equation} \mathcal{B}_+ = m+\frac{2(m+1)}{(m+1)-nq_m}-\frac{2(m+1)}{1+d(r_2/r_s)^{2 {\rm m}+2}}, \end{equation} \tag{ 21 }$$

where $r_2/r_s = \sqrt{nq_m/(1 + m)}$ and d determines the behaviour of ξ_m + 1 at r_s . For an ideal mode, we choose d =−1 so that ξ_m + 1 is finite at r_s .

With a parabolic $p_{||}^{(0)}$ and $\bar{p}$, by taking $(\mathcal{M}^2)'\propto r$ it is possible to find an exact solution of (16) [28, 67, 76]. The result, however, is rather convoluted. A great deal of simplification is achieved by approximating the plasma pressure within the shear-free region with a Heaviside step function [64, 77], namely $p_{||}^{(0)}\propto H(r_p-r)$ with $0\lt r_p\lt r_2$ (this is shown in figure 1). Despite the crudeness of this approximation, all the important physical ingredients are retained. Due to the rotation being constant in r, it then follows that the Mach number $\mathcal{M}(r)$ is also distributed with a Heaviside step, so that $(\mathcal{M}^2)'\approx\delta(r-r_p)[\mathcal{M}^2(r_p^+)-\mathcal{M}^2(r_p^-)] = \delta(r-r_p)\Delta_M$. Similarly, we write $\alpha\approx-2R_0[\bar{p}(r_p^+)-\bar{p}(r_p^-)]q^2\delta(r-r_p) = r_p\delta(r-r_p)\bar{\alpha}$. We choose $\bar{p}(r_p^+) = \bar{p}(r_2)$ and $\bar{p}(r_p^-) = \bar{p}(0)$ (cf figure 1). The amplitude of the Heaviside step is determined by arguing that the volume averaged $\beta = \frac{2\mu_0}{B^2}\frac{\int_0^a\sqrt{g}\bar{p}dr}{\int_0^a\sqrt{g}dr}$ with $\sqrt{g}\approx \kappa rR_0$ is identical to the one having parabolic $\bar{p}$ everywhere (see figure 1). Thus $r_p = r_2/\sqrt{2}$ and $\bar{\alpha} = 2\beta(r_2/a)^2q_m^2/\bar{\varepsilon}$ with $\bar{\varepsilon} = r_p/R_0$. Finally, it is worth pointing out that within our model equations the ratio $T_\perp/T_{||}$ can be varied independently of β, which is kept constant throughout the following analysis as well as the current channel width r₂, by allowing variations in $T_{||}$. We nevertheless point out that other constraints, e.g. on the cyclotron frequencies Ω_s of the species s (related to the magnetic field strength) which have to fulfil $\gamma/\Omega_s\ll1$, may impose further limitations on how specific equilibrium quantities can be varied independently of each other. However, we argue that the amount of the independent variation of the physical parameters associated with the equilibrium configurations of interest is consistent with the model assumptions

Zoom In Zoom Out Reset image size

Note that, alternatively, we could have had chosen to work with a constant plasma current I_p , or equivalently q(a), which corresponds to having $\beta_N\propto\beta aB/I_p$ constant (this may be more suitable for current tailoring studies). Thus, exploiting the δ-function behaviour of α and $(\mathcal{M}^2)'$ and assuming q_m ≈m/n, integration of equation (16) across r_p [64, 77] yields (note that integrating from 0 to r₂ gives the same result)

$$\begin{align} \nonumber\frac{\bar{\alpha}^2}{2\bar{\varepsilon}}\Lambda+ &\bar{\alpha}\left(\frac{3e}{4}-\frac{3}{2}\frac{e\delta}{\varepsilon}\right)+ \frac{(1+e)}{\bar{\varepsilon}}\big[(\frac{\delta q}{m/n})^2\unicode{x27e6} r\xi_m'\unicode{x27e7}_{r_p}- \end{align}$$

$$\begin{align}&\frac{r_p\bar{\alpha}}{R_0}(1-\frac{n^2}{m^2})\big]-\frac{m^2\Omega^2}{\bar{\varepsilon}n^2\omega_A^2}\Delta_M = \bar{\alpha}\tau \end{align} \tag{ 22 }$$

where $\unicode{x27e6} A\unicode{x27e7}_{r} = A(r^+)-A(r^-)$ having normalised $\xi_m(r_p) = 1$. Note that in references [27, 28] where the m = 1 mode is studied and the quantity 1 − 1/q² vanishes, the Mercier term refers to the second term on the lhs in the equation above, i.e. the geometrical shaping contribution to the magnetic well. Nevertheless, as a matter of terminology, hereafter we call the Mercier contribution the term $1-n^2/m^2$. The constant Λ accounts for the sideband coupling, i.e. L₊, so that employing (21) gives

$$\begin{equation} \Lambda = (1+m)(1-2e) \left[(\frac{1+m}{m})^{1+m}-2\right](\frac{r_p}{r_2})^{2 {\rm m}+2}. \end{equation} \tag{ 23 }$$

Finally, the solution of (16) that is continuous at r_p and fulfilling the boundary conditions at r₂ is

$$\begin{equation} \xi_m = \left\lbrace\begin{aligned} &(r/r_p)^{m-1},& r \lt r_p, \\ &\frac{(r/r_2)^{m-1}-(r/r_2)^{-m-1}} {(r_p/r_2)^{m-1}-(r_p/r_2)^{-m-1}}, &r\gt r_p, \end{aligned} \right. \end{equation} \tag{ 24 }$$

from which $\unicode{x27e6} r\xi_m'\unicode{x27e7}_{r_p} = -2 {\rm m}/[1-(r_p/r_2)^{2 {\rm m}}]$. Hence, collating these results, equation (22) can be solved for the degree of anisotropy that yields marginal stability, i.e. we solve for marginal τ and hence marginal $T_{||}/T_\perp$ for e.g. a given q_m (a simplified expression for $T_{||}/T_\perp$ in a weak anisotropic case can be obtained by taking $\tau\approx\frac{3}{4}(\frac{T_{||}}{T_{\perp}}-1)$). An example of the stability regions in the $(q_m,T_{||}/T_\perp)$ plane is shown in figure 2.

Zoom In Zoom Out Reset image size

Let us first note that, as pointed out in section 4.1, the last term in the lhs of (22) is always negative for monotonically decreasing pressure profiles, so that its effect is stabilising [31]. It is worth noting that modes with a sufficiently large m are stable since the field line bending contribution eventually dominates over the destabilising contribution due to the pressure driving terms. We shall also note that large m, i.e. short wavelength, perturbations of infernal type might be suppressed by higher order effects, such as e.g. diamagnetic corrections [78], which have not been included in our analysis. Finally, for given q_m , the most unstable mode is expected to be the one with m and n coprime with q_m ≈m/n (i.e. the lowest m resonating mode).

It is immediate to verify that, for a fixed $\bar{\alpha}$ and m > 1, smaller aspect ratio configurations are prone to exhibit enhanced stability by having a stronger stabilising contribution from the Mercier term (always having q_m  > 1). This, however, does not hold for modes with m/n = 1, for which such term vanishes. In agreement with previous studies, we find that with anisotropic temperatures, longitudinal injection ($T_{||}\gt T_\perp$) improves stability, while transverse injection ($T_\perp\gt T_{||}$) degrades it [36–40, 45]. Assuming e > 0, we recover the elongation stabilising effect at high plasma pressure, with the stabilising influence of a positive plasma triangularity (this was previously noted in reference [28] for the standard MHD isotropic case $T_\perp = T_{||}$). Pressure destabilisation at low β and sufficiently large vertical elongation is more easily achieved by the m = 1 mode because of the vanishing Mercier contribution (which is eventually dominant for high m modes). An example of the effect of plasma triangularity on the stability boundaries of a MAST-like configuration is shown in figure 3.

Zoom In Zoom Out Reset image size

Note that for negative triangularity plasmas, oblate cross sections might be beneficial in keeping the product eδ positive, and therefore deepening the magnetic well [79]. However, one might have the side effect of enhancing the coupling with the sidebands, and therefore worsening the stability. Thus, a careful optimisation of this effect by considering the global stability against a broader spectrum of perturbations may be required. The case of inverted q profiles with a core localised high shear region is the aim of the analysis in the next section.

Let us allow a magnetic shear reversal in the region 0 < r < r₁ with $r_1\lt r_2$. We recall that in the following analysis dominant harmonics with m > 1 only are considered. This is indeed appropriate for describing experimental configurations with negative central shear which have the minimum of the safety factor well above unity [11, 74]. Since plasma anisotropy and toroidal flow, and also triangularity, enter through a modification of the magnetic well, the conclusions drawn in the previous section on their effect on the stability boundaries hold in inverted q configurations as well. Elongation, on the other hand, affects the dynamics in a more subtle manner.

Following the same analysis presented in the previous section, in approximating the pressure within $r_1\lt r\lt r_2$ with a step function (cf figure 1) and retaining the same β of the parabolic profile, we choose $r_p = \sqrt{(r_1^2+r_2^2)/2}$ and $\bar{p}(r_p^-) = \bar{p}(r_1)$. Therefore, we have $\bar{\alpha} = 2\beta[(r_2/a)^2-(r_1/a)^2]q_m^2/\bar{\varepsilon}$. As shown in section 4.2, the main harmonic is vanishing in the high shear regions. Therefore, maintaining the normalisation $\xi_m(r_p) = 1$, the solution of (16) for r > r_p is identical to (24), while for r < r_p reads

$$\begin{equation*} \xi_m = \frac{(r/r_1)^{m-1}-(r/r_1)^{-m-1}}{(r_p/r_1)^{m-1}-(r_p/r_1)^{-m-1}}, \end{equation*}$$

yielding

$$\begin{equation} \unicode{x27e6} r\xi_m'\unicode{x27e7}_{r_p} = -\frac{2 {\rm m}\{(r_p/r_2)^{2 {\rm m}}-(r_p/r_1)^{2 {\rm m}}\}}{\{1-(r_p/r_1)^{2 {\rm m}}\}\{1-(r_p/r_2)^{2 {\rm m}}\}}. \end{equation} \tag{ 25 }$$

By repeating the logical steps employed in the previous section, we arrive at (22) where the constant Λ is now defined by $\Lambda = \frac{1-e/2}{\bar{\alpha}}\sum_\pm\frac{r_p^{\pm m}L_\pm}{1\pm m}$. In an inverted shear configuration the computation of the constants $L_\pm$ requires a more careful treatment compared to the monotonic q case. The difficulty arises because of the coupling due to elongation which yields non-trivial expressions for the satellite harmonics. Let us start by expanding the sideband perturbations according to $\xi_\ell = X_\ell+eY_\ell$ with $Y_\ell/X_\ell\sim 1$ and $\mathcal{L}_{m\pm1}(X_{m\pm1}) = 0$ (cf (20)). It is helpful to introduce the constant

$$\begin{align*} & \mathcal{H}_\pm(C,B,\ell) = \nonumber\\ & \frac{\bar{\alpha}(1\pm m)(2\pm m+C)(2\pm m+B)(r_p/r_2)^{2\pm\ell}}\nonumber\\ & \times\! {(\mp m\pm\ell-B)(2\pm m+C)\!-\!(\frac{r_1}{r_2})^{2\!\pm\!\ell}(\mp\! m\pm\ell-C)(2\pm m+B)}. \end{align*}$$

Writing $L_\pm = \bar{L}_\pm+e\hat{L}_\pm$, and following the standard procedure outlined in references [10, 64], equation (15) yields $\frac{r_p^{\pm m}\bar{L}_\pm}{1\pm m} = \mathcal{H}_\pm(\bar{\mathcal{C}}_{\pm},\bar{\mathcal{B}}_{\pm},2 {\rm m})$ and $\frac{r_p^{\pm m}\hat{L}_\pm}{1\pm m} = - \mathcal{H}_\pm(\hat{\mathcal{C}}_{\pm},\hat{\mathcal{B}}_{\pm},2 {\rm m}) [\frac{3}{2}-(\frac{1\pm m}{1\mp m})\frac{r_p^{\mp m}\bar{L}_\mp/4}{\mathcal{H}_\pm(\hat{\mathcal{C}}_{\pm},\hat{\mathcal{B}}_{\pm},0)} ]$ where $\bar{\mathcal{C}}_\pm = rd\text{ln}\,X_{m\pm1}/dr|_{r_1}$, $\bar{\mathcal{B}}_\pm = rd\text{ln}\,X_{m\pm1}/dr|_{r_2}$, $\hat{\mathcal{C}}_\pm = rd\text{ln}\,Y_{m\pm1}/dr|_{r_1}$ and $\hat{\mathcal{B}}_\pm = rd\text{ln}\,Y_{m\pm1}/dr|_{r_2}$. It is worth stressing that, since the dependence in e of ξ_m is missing, there is no need to expand this quantity in the elongation parameter.

In order to compute $X_{m\pm1}$ and $Y_{m\pm1}$, we turn to equation (19) which is more easily manipulated when it is expressed in terms of the perturbed poloidal flux $\tilde{\psi}_\ell = -\kappa r(1/q-n/\ell)\xi_\ell$. This also is expanded in e yielding $\tilde{\psi}_\ell = \Psi_\ell+e\chi_\ell$ with $\chi_\ell/\Psi_\ell\sim 1$. By linking $X_\ell$ and $Y_\ell$ to $\tilde{\psi}_\ell$, it follows that $X_\ell = -\Psi_\ell/[r(1/q-n/\ell)]$ and $Y_\ell = -X_\ell-\chi_\ell/[r(1/q-n/\ell)]$. Hence, to leading order equation (19) gives ($\ell = m\pm1$)

$$\begin{equation*} {\boldsymbol{\nabla}}^2_\star\Psi_\ell\equiv (r\Psi_\ell^{\prime})^{\prime}-\Big[\frac{\ell^2}{r}+\frac{(\frac{2}{q}-\frac{rq^{\prime}}{q^2})^{\prime}}{1/q-n/\ell}\Big]\Psi_\ell = 0, \end{equation*}$$

which is equivalent to $\mathcal{L}_\ell(X_\ell) = 0$. To the next order in e, we obtain an equation for $\chi_{m\pm1}$

$$\begin{align} \nonumber{\boldsymbol{\nabla}}^2_\star&\chi_{m\pm1}+\frac{1}{2}\Big[(r\Psi_{m\mp1}^{\prime})^{\prime}\mp 2 {\rm m}\Psi_{m\mp1}^{\prime}+\frac{m^2-1}{r}\Psi_{m\mp1} \end{align}$$

$$\begin{align}&+\frac{ (\frac{rq^{\prime}}{q^2})^{\prime}(m\mp1)\Psi_{m\mp1} \mp\frac{2rq^{\prime}}{q^2}\Psi^{\prime}_{m\mp1}}{\frac{m\pm1}{q}-n} \Big] = 0. \end{align} \tag{ 26 }$$

Focussing on the region r₂ < r < a first, we employ the safety factor profile used in the previous section which is chosen such that the mode m + 1 is resonant at r_s  < a. By requiring the upper satellite harmonic to be finite at r_s with an ideal metallic wall at the boundary, we have $\Psi_{m-1} = a_-[(r/a)^{m-1}-(r/a)^{-m+1}]$ whereas Ψ_m + 1 = 0 for r_s  < r < a and $\Psi_{m+1} = a_+[(r/r_s)^{m+1}-(r/r_s)^{-m-1}]$ for $r_2\lt r\lt r_s$ (the computation of $\bar{\mathcal{B}}_\pm$ is straightforward). Proceeding further, let us denote the particular solution of (26) associated with the m + 1 mode with $\mathfrak{g}$ (which is proportional to a₋). By imposing $\chi_{m+1}(r_s) = 0$, we eventually have

$$\begin{equation*} \chi_{m+1} = b_+[(\frac{r}{r_s})^{m+1}-(\frac{r}{r_s})^{-m-1}]-\mathfrak{g}(r_s)(\frac{r}{r_s})^{-m-1}+\mathfrak{g}, \end{equation*}$$

In the region 0 < r < r₁ we use a safety factor of the form [53]

$$\begin{equation*} \frac{1}{q} = \frac{1}{q_0}+(\frac{1}{q_m}-\frac{1}{q_0})(r/r_1)^2. \end{equation*}$$

With such a profile, the sideband equations can be expressed exactly in terms of hypergeometric functions [53, 80]. Retaining the full solutions might be helpful when strongly reversed safety factor profiles with $q_0\gg1$, i.e. advanced scenarios, are considered which, however, are beyond the scope of this work. The analysis of the problem tackled here, can be significantly simplified by assuming $q_0/q_m$ not too large. Hence, within this approximation, we have $\Psi_{m-1}\approx A_-(r/r_1)^{m-1}$ and $X_{m-1}\propto r^{m-2}$. This yields $\bar{\mathcal{C}}_- = m-2$, and therefore $\bar{L}_-/(1-m) = 0$. By using (15), it follows that

$$\begin{equation} [X_{m\pm1}(2\pm m+\bar{\mathcal{B}}_\pm)]_{r_2} = r_2^{\pm m}\bar{L}_\pm, \end{equation} \tag{ 27 }$$

from which a₋ = 0 giving $\mathfrak{g} = 0$ and $\hat{\mathcal{B}}_+ = \bar{\mathcal{B}}_+$. Plugging the function Ψ_m − 1 determined above into (26) and neglecting the terms proportional to the magnetic shear, it is promptly verified that ${\boldsymbol{\nabla}}_\star^2\chi_{m+1} = 0$ which yields $\hat{\mathcal{C}}_+ = \bar{\mathcal{C}}_+$ and, as an immediate consequence, $\hat{L}_+ = -{3}/{2}\bar{L}_+$. It is worth pointing out that in the limit r₁→ 0 the expressions derived in the previous section are recovered, also for the m = 1 case.

In order to have the driving term Λ fully determined, it only remains to compute the constant $\hat{L}_-$. Noticing that, depending on the value of q₀ the mode m + 1 might have a resonance at $\bar{r}_s\lt r_1$, two cases are then examined in the next subsections.

6.1. Weak reversal (nq₀ < m + 1)

Here we employ the weakly reversed q approximation, viz. $\frac{1/q_m-1/q_0}{1/q_0}\ll1$, so that $X_\ell\propto\Psi_\ell/r$. Let us denote $\Delta\mu = 1/q_m-1/q_0$, $\bar{m} = \sqrt{(m+1)^2+8}$ and $\zeta = n/[m(m+1)]$. By introducing the variable $z = (r/r_1)^2\Delta\mu/[\Delta\mu-\zeta]$, the upper harmonic can be written in terms of the hypergeometric function F

$$\begin{equation*} X_{m+1}\propto \frac{y^{m/2}}{(1-y)^{(\bar{m}-3)/2}} F(A,A+1;2+m;y), \end{equation*}$$

with $y = z/(z-1)$ and $A = (m+1-\bar{m})/2$. Since y < 1, this expression yields approximately $\Psi_{m+1}\approx A_+(r/r_1)^{m+1}$ and $X_{m+1}\propto r^m$ so that $\bar{\mathcal{C}}_+ = m$. By employing equation (27) with the replacements $r_2\to r_1$ and $\bar{\mathcal{B}}_\pm\to \bar{\mathcal{C}}_\pm$, it readily follows that $A_+ = \frac{r_1\bar{\alpha}}{2}[n-(m+1)/q_0)(r_1/r_p)^{m}(r_p/r_2)^{2+m}[(r_s/r_2)^{2+2 {\rm m}}-2]$.

Because of the absence of internal resonances, in (26) we drop effects linear with respect to the magnetic shear. Hence it is straightforward to show that

$$\begin{equation*} \chi_{m-1} = C(r/r_1)^{m -1}-A_+ \frac{1+m}{2} (r/r_1)^{1+m}, \end{equation*}$$

where C is a constant. We point out that a similar expression could have been derived by using (15) with α → 0 and the obvious replacement $L_\pm\to A_\pm$.

Similarly to the derivation of (27), expanding equation (15) further in e yields

$$\begin{equation} [(2-m)Y_{m-1}+rY'_{m-1}]_{r_1} = r_1^{-m}\hat{L}_-+\frac{1}{2}(\frac{m-1}{m+1})r_1^m\bar{L}_+. \end{equation} \tag{ 28 }$$

By plugging the expression for χ_m − 1 into the equation above and taking q₀≈n/m, we obtain

$$\begin{equation*} \frac{r_p^{-m}\hat{L}_-}{1-m} = \frac{m(1+m)}{nr_1}A_+(\frac{r_1}{r_p})^m+\frac{1}{2}(\frac{r_1}{r_p})^{2 {\rm m}}\frac{r_p^m\bar{L}_+}{1+m}. \end{equation*}$$

It can be shown that $\hat{L}_- = 0$ so that the constant Λ is given by (23). Therefore, the stability boundaries are determined by the set of equations (22), (23) and (25). It is immediate to note that, with r₁ > 0, the stability is improved by having a stronger field line bending stabilising contribution. As previously noted, the expressions for the monotonic q case are recovered by letting r₁→ 0. In the next subsection we allow for a stronger shear reversal which produces an internal resonance of the m + 1 mode.

6.2. Moderate reversal (nq₀ > m + 1)

The analysis of the moderately reversed q configuration essentially repeats the steps outlined in the paragraph above. The only difference is a more complicated dependence upon the radial variable of the upper sideband. For nq₀ > m + 1 the harmonic X_m + 1 has a resonance at $\bar{r}_s = r_1\sqrt{\frac{1/q_0 - n/(1+m)}{1/q_0-n/m}}\lt r_1$, hence for $0\lt r\lt\bar{r}_s$ we have X_m + 1 = 0, whereas for $\bar{r}_s\lt r\lt r_1$ the eigenfunction is (F is the hypergeometric function) [53]

$$\begin{align*} \nonumber X_{m+1}\propto&\frac{z^{-(\bar{m}+3)/2}}{1-1/z}[F(B,B-C+1;B-A+1;1/z)+\\\label{FULLXIMPLUS1}&Dz^{\bar{m}}F(A,A-C+1;A-B+1;1/z)], \end{align*}$$

where $B = (m+1+\bar{m})/2$, C = 2 + m and D is such that X_m + 1 is finite at $\bar{r}_s$.

We find that, far from $\bar{r}_s$, the leading order of the upper sideband displacement is well approximated by $X_{m+1}\propto r^{\bar{m}-1}/[1/q-n/(m+1)]$. Since $|D|\gt 1$, we may take $\Psi_{m+1} = A_+(r/r_1)^{\bar{m}}$ for $\bar{r}_s\lt r\lt r_1$ and Ψ_m + 1 = 0 for $0<r<\bar{r}_s$. In analogy with the case of weak reversal, we obtain $A_+ = -r_1\frac{n}{m}(r_1/r_p)^mr_p^m\bar{L}_+/[(1+m)(2+m+\mathcal{C}_+)]$ where

$$\begin{align} & \frac{r_p^m\bar{L}_+}{1+m} \nonumber\\ & = \frac{\bar{\alpha}(1+m) [\bar{m}+(1+m)(2 {\rm m}/nq_0-1)]Z}{ \bar{m}+(1+m)(\frac{2 {\rm m}}{nq_0}-1)+(\frac{r_1}{r_2})^{2+2 {\rm m}}[\bar{m}+(1+m)(\frac{2 {\rm m}}{nq_0}-3)]Z } \end{align} \tag{ 29 }$$

with $Z = (r_s/r_2)^{2+2 {\rm m}}-2$.

It can be shown that the particular solution of the m − 1 mode is not significantly affected by dropping the last term in square brackets in (26). Hence, for the sake of simplicity in the resulting expressions, this term will be neglected. It follows that in the region $\bar{r}_s\lt r\lt r_1$ the solution of (26) for the m − 1 harmonic reads

$$\begin{equation*} \chi_{m-1} = C_0(\frac{r}{r_1})^{m -1}+C_1(\frac{r}{r_1})^{-m +1}- \frac{A_+}{2}(\frac{\bar{m}+m+1}{\bar{m}-m+1})(\frac{r}{r_1})^{\bar{m}}. \end{equation*}$$

Under the assumption $q_0/q_m$ is not too large (see section 6) we take $(m-1)/q-n\approx (m-1)/\bar{q}-n$ with $\bar{q} = 2/(1/q_0 + n/m)$. In order to find $\hat{L}_-$ it is sufficient to specify the constant C₁. Some simple integrations of (26) across $\bar{r}_s$ show that $\chi_{m-1}'$ has a jump at $\bar{r}_s$, whereas χ_m − 1 remains continuous. These conditions read $\unicode{x27e6}\chi_{m-1}\unicode{x27e7}_{\bar{r}_s} = 0$ and $\unicode{x27e6} r\chi_{m-1}'\unicode{x27e7}_{\bar{r}_s} = -{1}/{2}\bar{m}(\bar{r}_s/r_1)^{\bar{m}}A_+$. By applying such constraints we obtain

$$\begin{equation*} C_1 = \frac{(1+m+2\bar{m})}{4(1-m)} (\frac{\bar{r}_s}{r_1})^{\bar{m}+m-1}A_+. \end{equation*}$$

If $\bar{r}_s/r_1$ is sufficiently small, we let C₁→ 0. Hence, by means of (28), after some straightforward manipulations we finally have $\frac{r_p^{-m}\hat{L}_-}{1-m} = K\times\frac{r_p^m\bar{L}_+}{1+m}$ with

$$\begin{align} K \!=\! \frac{1}{2}(\frac{r_1}{r_p})^{2 {\rm m}}\big\{ 1\!+\! \frac{(1+m+\bar{m})n/m}{[(m-1)/\bar{q}-n][\bar{m}+(1+m)(2 {\rm m}/nq_0-1)]} \big\}. \end{align} \tag{ 30 }$$

Therefore, the parameter Λ expanded to leading order in e reads

$$\begin{equation*} \Lambda = \frac{[1+e(K-2)]}{\bar{\alpha}}\frac{r_p^m\bar{L}_+}{1+m}, \end{equation*}$$

which is computed by means of (29) and (30). It is found that the quantity K is rather small, and in the range of experimental relevant parameters we may approximate Λ with (23). Thus, in analogy with the weakly reversed q case, the stability boundaries are identified by the set of equations (22), (23) and (25). Thus, the same conclusions drawn in section 5 hold also for weakly and moderately reversed q configuration. We note in particular that the marginal β has a very weak dependence upon q₀ and for sufficiently small $r_1/r_p$ the stability boundaries of monotonic and weakly reversed configurations coincide [7]. This might be somehow expected, since according to references [7, 11, 12] the eigenfunctions are relatively small in the region of the shear reversal indicating weak contributions associated with this region. Indeed, the fluid displacements, viz. the eigenfunctions, tend to be more localised where either the shear is small or in the external sheared region.

The aim of this section is to point out the robustness of the anisotropic model and to discuss to what extent our results, in particular flow modifications, differ from previous findings found in the literature (see e.g. [52]). Such differences can be attributed essentially to (i) the absence of rotation and averaged density gradients (this has been chosen for mathematical ease) and (ii) the closure model.

Focussing on the first effect, it is well known that the allowance of rotation and density gradients might have a profound impact on the stability properties. Indeed, in reference [52] it was found that the additional contributions to the magnetic well (the dA₂/dr term given by equation (30) in the reference above) due to realistic radial variation of Ω and $\bar{\rho}$ can be considerably larger than the last term in (17). Radial variations of the equilibrium averaged density can be included (e.g. in the model of the current paper) rather simply, yielding

$$\begin{align*} w^{\prime} = &\frac{\hat{\alpha}}{R_0}[1-\frac{1}{q^2}+(1-e)\tau -\frac{3e}{4}+\frac{3}{2}\frac{e\delta}{\varepsilon}]+\\ &q^2\{\frac{\Omega^2}{\omega_A^2}[2-\frac{1}{q^2}-\frac{3}{2}e+(1-e)\mathcal{M}^2]\}^{\prime}, \end{align*}$$

where $\omega_A(r) = R_0^2\bar{\rho}(r)$ and $\hat{\alpha} = -2R_0[\bar{p}(1+\mathcal{M}^2)]^{\prime}q^2$ ( (we also must replace α with $\hat{\alpha}$ in equations (16) and (19)). On the contrary, dealing with rotation gradients is much more difficult. This is because, if rotation is radially dependent, we cannot set to zero the Doppler shifted eigenvalue across the whole low-shear region. This, therefore, requires a more elaborate calculation of the perturbed distribution function (cf equation (11)), which yields a complicated expression of the plasma response written in terms of the plasma dispersion function [44]. We nevertheless stress that within the infernal framework and under appropriate approximations, effects of sheared flows of the order $a\Omega^{\prime}/\Omega\sim1$ can be included in the analysis, thus showing off the power of the approach. The inclusion of these particular flows demonstrates that a generalisation beyond simple anisotropy is possible, and that future work will attempt to take into account flows that resemble more realistic experimental situations. We finally point out that a large flow shear may break the assumptions employed in our model, in particular the gyrotropic assumption [81, 82]. However, as long as we do not deal with instabilities localised in extremely pronounced ITBs or at the pedestal where the flow shear can be very large, the assumptions within our framework is developed should be reasonably fulfilled.

The other differences arise from the closure model. As is well know, in standard MHD we employ the adiabatic closure for which $d(p\rho^{-\Gamma})/dt = 0$ where Γ = 5/3. In the guiding centre plasma (GCP) approximation, the closure is provided by solving the one dimensional collisionless kinetic equation, i.e. (4), which describes the particle motion parallel to the magnetic field [42, 43]. It was pointed out in reference [42] that MHD marginal stability boundaries are recovered by the isotropic GCP model (in the case of shearless toroidal flows in cylindrical geometry) if Γ = 1, i.e with the plasma being isothermal. This is indeed apparent by inspecting equations (12) and (13). In standard MHD in toroidal geometry, the adiabatic index enters inertia and magnetic well terms [77]. In toroidal geometry, these contributions reduce to the ones obtained in the isotropic GCP limit for Γ = 1 at marginal stability with uniform flows. Thus, it follows that within the GCP model the Brunt-Väisälä stabilisation mechanism, coming from the centrifugal force which gives a stable density or entropy distribution within each magnetic surface, is lost (see e.g. [30, 83, 84] and references therein). We shall point out that other types of closure may be used [85, 86], and different results might be expected depending on the closure model.

It is worth pointing out the robustness of the anisotropic fluid model used here, for cases where strong parallel anisotropy is applicable to the experimental regime of interest. It is a more complete model in the deep parallel anisotropic limit, for low or high performance plasmas, than the isotropic MHD limit for cases where isotropy applies to a fusion grade plasma. Indeed, in the limit of strong parallel anisotropy, i.e. $T_{||}/T_\perp\to\infty$, the trapped particle fraction vanishes so that kinetic corrections to kinetic-MHD models are negligible [45]. Although kinetic models are dependent upon several parameters (collisional regime, temperature, bounce frequency resonance effects, etc), in the parallel anisotropic limit kinetic models yield weak corrections to a fluid model. What is left are anisotropic fluid effects which are robustly stabilising since pressure perturbation associated with passing particles are larger on the high field side of the flux surface (contrarily to an isotropic plasma in which passing particles compensate for trapped particles, where the latter tend to spend more time in the low field side).

It is also possible to see that combined elongation and anisotropy effects are recovered by inspecting equation (3) in reference [45], noting that $p_\perp$, $p_{||}$ and the factor $C = \sum_s m_s\int d\mu d\epsilon\frac{B}{v_{||}}(\mu B)^2\frac{\partial f_s}{\partial\epsilon}$ [34] can be written easily in terms of cosθ, $T_{||}$ and $T_\perp$ for the model used. This equation indicates also that the effect of elongation would yield exactly the elongation correction seen in the Mercier index (−eτ) in equation (17). Note that for parallel anisotropy, changes to the metric are weak, allowing the toroidal and shaping expansions presented in section 3. Toroidicity and elongation effects combine with anisotropy in such a simple way because anisotropy introduces a low order poloidal dependence, which combines with the poloidal dependence of 1/R. The other shaping effects are additive, giving the well known shaping effects on interchange modes [26]. Thus, as the physical effects are additive in nature, and fairly intuitive, more realistic, though perhaps ad-hoc, approaches to including combined fluid anisotropy and flows can be sought.

Finally, it is worth pointing out that parallel anisotropy provides stabilisation also in stellarators, whose fields decay proportionally to 1/R [87].

The problem of plasma anisotropy and certain limited equilibrium flow effects in shaped hybrid plasmas has been addressed. As discussed in the introduction, the relevance of the results presented in this work can be extended both to spherical and standard tokamaks. Using simplified class of profiles allows to have a manageable analytical treatment, nevertheless retaining relevant physical effects. The main result of this work is equation (22) which describes the stability boundaries for ideal infernal modes. Such an equation holds both for flat and q profile reversed configurations, whose specific case is identified by selecting appropriately the field line bending contribution. In line with previous results, beyond a trivial redefinition of the ballooning parameter α due to a modified averaged total pressure, it is found that plasma anisotropy effects enter through a modification of the magnetic well, yielding better stability properties with tangential injection. It is found that a uniform toroidal flow improves stability as well. Positive triangularity effects in vertically elongated plasmas are stabilising, whereas negative triangularity tends to be destabilising. We notice from figure 3 that even with a relative modest parallel pressure anisotropy ($T_{||}\gt T_{\perp}$), that should be easily realisable with the MAST neutral beam system, the infernal mode can be completely stabilised even when the minimum of q approaches unity. Oblate cross sections might be beneficial for restoring the stabilising contribution to the magnetic well in negative triangularity plasmas.

Depending on β, elongation might be either stabilising or destabilising. A positive elongation parameter e at low pressure and negligible triangularity tends to be destabilising with a sufficiently small filed line bending contribution, while its effect turns out to be stabilising at sufficiently high β in accordance with [28]. Note that equation (22) seems to suggest that mode stabilisation could be achieved not only with current profile optimisation, i.e. modifying the safety factor, but also with a careful tailoring of the plasma shaping. We point out that the stability eigen-equations for the Fourier harmonics derived in this work, and in particular the equation for the dominant mode (16) with the allowance for residual magnetic shear effects in the field line bending contribution, might be helpful for analysing a broader class of instabilities, e.g. ballooning modes.

Further work is nevertheless needed to assess more rigorously finite aspect ratio effects, highly relevant for compact machines when ɛ ∼ 1, and the coupling effects resulting from the inclusion of a larger number of satellite Fourier harmonics, viz. allowing the perturbation to be more ballooned. Moreover, additional analysis is required for the description of cases in which the safety factor is strongly reversed, such as either advanced scenarios or current hole configurations. Finally, we point out that the penetration of beams in a ST reactor might not be sufficient to provide the required levels of anisotropy envisaged to stabilise the MHD activity analysed in this work. Further work is therefore needed to assess the overall effect of NBI heating in future experiments, and particularly the role of shaping as a primary player in mode stabilisation. It is envisaged that such challenging analysis, might be more suitably addressed numerically.

This work has been funded by the RCUK Energy Programme [grant number EP/T012250/1]. To obtain further information on the data and models underlying this paper please contact PublicationsManager@ukaea.uk.

The authors gratefully acknowledge Dr J Connor for his valuable comments.

We consider a safety factor of the form shown in figure 1, i.e. monotonic, and a uniformly flat rotation profile. In the eventuality that some residual resistivity is allowed at r_s , the energy method is not valid any longer to obtain the stability boundaries. These are instead heuristically identified by the condition $\Delta^{\prime} = \frac{r\tilde{\psi}_{m+1}^{\prime}}{\tilde{\psi}_{m+1}}|_{r_s} = 0$ [88] where $\sqrt{g}\tilde{B}^r\approx-\partial_\vartheta\tilde{\psi}$ (see section 6). Hence at marginal stability, the constant d appearing in (21) is set to $d = -(r_s/a)^{2 {\rm m}+2}$. It immediately follows that we must replace (23) with

$$\begin{equation*} (1+m)\frac{\bar{\alpha}^2}{2\bar{\varepsilon}}(1-2e) [(\frac{a}{r_2})^{2+2 {\rm m}}-2](\frac{r_p}{r_2})^{2 {\rm m}+2}. \end{equation*}$$

Noting that $(r_s/r_2)^{2+2 {\rm m}} = (\frac{1+m}{m})^{1+m}$, the destabilising effect of plasma resistivity is clearly apparent. Further work is nonetheless needed to assess the effects of plasma compressibility and sheared flows in a more rigorous manner and the possibility of a second resonance in case of strongly inverted q profiles.