2.1 Numerical equilibria

We generate equilibria numerically using the three-dimensional (3-D) equilibrium code VMEC (Hirshman & Whitson Reference Hirshman and Whitson1983). VMEC works by minimizing the integral

(2.7)\begin{align} W = \int\left(\frac{p}{\gamma -1} + \frac{B^2}{2 \mu_0}\right) \,{\rm d} V, \end{align}

subject to multiple constraints, which for axisymmetric equilibria is equivalent to solving the Grad–Shafranov equation (Kruskal & Kulsrud Reference Kruskal and Kulsrud1958). For our study, VMEC takes the shape of the boundary surface along with the global radial pressure $p(s)$ and safety factor $q(s)$ or enclosed toroidal current $G(s)$ as a function of the normalized toroidal flux $s$. It then creates flux surfaces to minimize the integral in (2.7) on each surface for a fixed $p$ and $q$ (or $G$) subject to various constraints. We choose the safety factor instead of toroidal current as it varies slowly (Abel & Cowley Reference Abel and Cowley2013) compared with the plasma currents in the limit of small electron-to-ion mass ratio.

In the following paragraphs, we will explain the process of generating the data for this study using VMEC. For each VMEC equilibrium, we pick two different radially local regions – equilibria that satisfy (2.6) and are localized to a flux surface. It is important to point out that for most of this paper, we will only investigate modes that have small wavelengths perpendicular to the field line, i.e. modes that are localized to a flux surface, henceforth referring to our study as a local stability analysis. Since our aim is to have ample variability in our input data, at the end of this section we justify our choices by looking at the chosen flux surfaces and their corresponding $\beta$ values.

We produce high-radial-resolution equilibria using the fixed-boundary solver in VMEC after providing it with an L-mode-like pressure profile $p = p(s)$ – a profile that does not have a sudden drop in pressure over a short radial distance, a monotonic safety factor $q = q(s)$ profile as a function of the normalized toroidal flux $s = \chi /\chi _{\mathrm {LCFS}}$ and the Last Closed Flux Surface (LCFS) shape. We choose a simple form for the profiles $p = p(s)$ and $q = q(s)$ given by

(2.8)\begin{equation} \left.\begin{array}{c@{}} p = n T, \quad p_0 = \tilde{p}_0 n_0 T_0 \\ n(s) = n_0(1 + \nu_n)(1-s^2)^{\nu_n},\quad T(s) = T_0(1 + \nu_T)(1-s^2)^{\nu_T}\\ q = q_0(1+s^{2 \nu_q})^{1/(2 \nu_q)} \end{array}\right\}.\end{equation}

The different parameters are given in table 1. The parameter $\tilde {p}_0 \in [1, 10, 70]$ for the low, intermediate and high-beta equilibria, respectively. For each triangularity value, we choose a different LCFS shape described by a Miller parameterization (Miller et al. Reference Miller, Chu, Greene, Lin-Liu and Waltz1998)

(2.9)\begin{equation} \left. \begin{array}{c@{}} R = R_0 + a \cos(t + (\sin^{{-}1}\delta)\sin t),\\ Z = a \kappa \sin(t). \end{array} \right\} \end{equation}

The parameter $t$ varies in the range $[-{\rm \pi}, {\rm \pi})$. The values of the rest of the parameters in (2.9) are given in table 2. The radial coordinate that we will use for all the stability analyses is $\rho = \sqrt {s} = \sqrt {\chi /\chi _{\mathrm {LCFS}}}$ since it is a better measure of the radial distance from the magnetic axis than the normalized poloidal flux $\psi /\psi _{\mathrm {LCFS}}$. The safety factor and pressure profiles as a function of $\rho$ are given in figure 1. For all our studies, we use the same safety factor $q$ and the normalized pressure profile $p$ with different values of $\tilde {p}_0$. In this way, we are able to create three different pressure profiles with on-axis $\beta \sim 0.01, 0.1, 1$ corresponding to $\tilde {p}_0 = 1, 10, 70$, respectively. Henceforth, we will refer to the equilibria with $\tilde {p}_0 = 1, 10, 70$ as low, intermediate and high-$\beta$ or $\beta \sim 0.01, \beta \sim 0.1$ and $\beta \sim 1$, respectively. We need to pick flux surfaces for our local stability analyses. In this study, we choose surfaces at normalized radii $\rho = 0.5$ and $0.8$. In total, there are twelve local equilibria in our study: $3 \beta$ values $\times 2$ boundary shapes $\times 2 \rho$ values. Because $\beta$ varies over a flux surface, it will be convenient to introduce a reference magnetic field for each global equilibrium and redefine

(2.10)\begin{gather} \beta(\rho) = 2 \mu_0 p(\rho)/B_{{N}}^2, \end{gather}

(2.11)\begin{gather}B_{{N}} = \chi_{\mathrm{LCFS}}/({\rm \pi} a_{{N}}^2), \end{gather}

where $B_{{N}}$ is a reference magnetic field and $a_{{N}}$ is the effective minor radius such that ${\rm \pi} a_{{N}}^2$ is equal to the area enclosed by the boundary and $\chi _{\mathrm {LCFS}}$ is the toroidal flux enclosed by the LCFS. For this study, $a_{{N}} = 0.684 \ \textrm {m}, B_{{N}} = 0.681\ \textrm {T}$. The values of $\beta$ obtained from VMEC are given in table 3.

Figure 1. This figure shows the safety factor and normalized pressure profiles used for creating the equilibria. The two red lines correspond to the values of the normalized radius $\rho$ at which the local equilibria will be analysed for their stability.

Table 1. VMEC equilibria input parameters. Throughout this study, in this table, every parameter remains fixed.

Table 2. Miller parameters for the outer boundary.

Table 3. Reference $\beta$ values for selected surfaces.

Table 4. Nominal gradient scale length values.

Each equilibrium has $512$ surfaces with each surface represented by $40$ poloidal modes. We found that the equilibria were converged with this choice of resolution. All the equilibria that we investigate in this study are up–down symmetric which is why we only show the upper half. The flux surface contours for the twelve equilibria are shown in figure 2.

Figure 2. This figure shows the flux surfaces for all the equilibria generated using VMEC. The local equilibria that will be studied in this paper are highlighted in red. The magnetic axis in each figure is the black cross.

The numerical high-$\beta$ equilibria show qualitative features like the vertical ‘core’ on the inboard side and thin boundary layer on the outboard side as shown by Hsu et al. It is interesting to see that the negative-triangularity equilibria high-$\beta$ equilibria are more strongly shaped than the positive-triangularity ones – the vertical inboard solution causes the flux surface to develop a ‘squareness’. We illustrate the ‘squareness’ in figure 3.Footnote ¹ Most importantly, these numerical equilibria do not suffer from any of the issues observed in Hsu et al. Therefore, all the resulting geometric coefficients are smooth which allows the local stability analyses in the following sections.

Figure 3. This figure shows two high-$\beta$ equilibria and their corresponding best-Miller fit. We can see that the fit for the negative triangularity is worse due to the ‘squareness’ of the flux surface on the inboard side. The agreement between gradients of various physical quantities will be even worse.

3.1 Physical and mathematical description

One of the most important MHD instability for us to investigate is the ideal-ballooning instability (Connor, Hastie & Taylor Reference Connor, Hastie and Taylor1979) – a field-aligned, pressure-driven Alfvén wave that grows when the destabilizing pressure gradient in the region of ‘bad’ curvature exceeds the stabilizing effect of field-line bending. The region of ‘bad’ curvature is a region of a flux surface where $\boldsymbol {\kappa } \boldsymbol {\cdot } \boldsymbol {\nabla } p > 0$, such that the field-line curvature is in the same direction as the plasma pressure gradient. For most tokamak equilibria, this region lies on the outboard side.

The equation governing the ideal-ballooning mode can be obtained by minimizing the ideal-MHD energy integral (Bernstein et al. Reference Bernstein, Frieman, Kruskal and Kulsrud1958) for incompressible modes in the limit of large toroidal mode number. Doing so gives us a differential equation that determines $X$, the radial displacement of said mode along a field line. To ensure that the displacement $X$ satisfies the periodicity condition on the surface of interest and nearby surfaces, one uses the ballooning transformation

(3.1)\begin{equation} X = \sum_{N ={-}\infty}^{\infty} \hat{X}(\theta - 2{\rm \pi} N)\, {\rm e}^{{\rm i}n(\phi - q (\theta- \theta_0 - 2{\rm \pi} N)))}, \quad N \in \mathbb{Z}, \end{equation}

subject to the condition

(3.2)\begin{equation} \lim_{\theta \rightarrow \pm \infty} \hat{X}(\theta; \psi, \theta_0) = 0, \end{equation}

and solves for $\hat {X}$. The variable $n$ is the toroidal mode number, $\theta _0$ is the ballooning parameterFootnote ² and the rest of the terms are defined in § 2. Upon minimizing the ideal-MHD energy integral and using the ballooning transformation, one obtains the ideal ballooning equation (Connor et al. Reference Connor, Hastie and Taylor1979; Dewar & Glasser Reference Dewar and Glasser1983)

(3.3)\begin{equation} \frac{1}{\mathcal{J}}\frac{\partial}{\partial \theta}\left( \frac{\lvert \boldsymbol{\nabla}\alpha \rvert^2 }{\mathcal{J}\, B^2} \frac{\partial \hat{X}}{\partial \theta}\right) + 2 \frac{{\rm d} p}{{\rm d} \psi}\left[ \boldsymbol{B}\times\boldsymbol{\nabla}\left(\mu_0 p + \frac{B^2}{2}\right)\boldsymbol{\cdot}\boldsymbol{\nabla}{\alpha}\right] \hat{X} = \rho \omega^2 \frac{\lvert\boldsymbol{\nabla}\alpha\rvert^2}{B^2} \hat{X}, \end{equation}

where $\rho$ is the plasma mass density and $\hat {X} = \hat {X}(\theta ; \psi, \theta _0)$ is the eigenfunction in ballooning space and $\omega ^2$ is the eigenvalue. The ballooning equation balances the stabilizing field-line bending term and destabilizing pressure gradient with the inertia of the resulting Alfvén wave, oscillating with a frequency $\omega$. Note that (3.3) depends on $\psi$ only as a parameter and we can compute the coefficients from the on-surface equilibrium quantities and their first derivatives. Therefore, it is possible to study the ballooning stability of the local equilibria that we chose in the previous section.

Due to the self-adjoint nature of ideal-MHD, all the eigenvalues ($\omega ^2$) of (3.3) will be real numbers. Hence, $\omega$ will either be purely real, an oscillating mode or purely imaginary, a growing mode. We refer to the oscillating modes as stable and growing modes as unstable. Since unstable modes are of more significance to us, we will plot the normalized growth rate when plotting the results in § 3.3. We define the normalized growth rate as follows:

(3.4)\begin{equation} \gamma ={-}{\rm i} \omega a_{{N}}/v_{\mathrm{th}, i}, \end{equation}

where $v_{\mathrm {th},i} = \sqrt {2T_{{i}}/m_{{i}}}$ is the ion-thermal velocity and $T_{{i}}$ and $m_{{i}}$ are the temperature and mass of the ions, respectively. We do this to establish a common normalization scheme throughout the paper.

3.2 The Greene–Chance analysis

To better understand the stability of a local equilibrium, we need the ability to vary that equilibrium. This can be done by changing the magnetic shear and pressure gradient independently about their nominal values – equivalent to varying the gradients of both the local current and plasma pressure – quantities that determine the solution to (2.6). This gives us the ability to generate multiple local equilibria satisfying the Grad–Shafranov equation and do a stability analysis without recalculating the global equilibrium. We define

(3.5)\begin{gather} \hat{s} = \frac{\rho}{q}\frac{{\rm d} q}{{\rm d} \rho}, \end{gather}

(3.6)\begin{gather}\alpha_{\mathrm{MHD}} ={-}\frac{2 \mu_0 \rho q^2}{\epsilon B_{{N}}^2}\frac{{\rm d} p}{{\rm d} \rho}, \end{gather}

as the magnetic shear and pressure gradient, respectively. This method of varying a local equilibrium through $\hat {s}$ and $\alpha _{\mathrm {MHD}}$ is known as an $\hat {s}-\alpha _{\mathrm {MHD}}$ analysis. This technique was developed by Greene & Chance (Reference Greene and Chance1981) and has been used extensively to study local stability of different axisymmetric equilibria (Connor, Hastie & Taylor Reference Connor, Hastie and Taylor1978; Connor et al. Reference Connor, Hastie and Taylor1979; Bishop Reference Bishop1985; Miller et al. Reference Miller, Chu, Greene, Lin-Liu and Waltz1998). Figure 4 illustrates the main point – we can change the gradient of the pressure and the safety factor locally by a finite amount without significantly changing their respective values. We will use this idea again in §§ 5 and 6 to vary the pressure gradient at the nominal magnetic shear when we examine the microstability of different equilibria. Details explaining the Greene–Chance analysis are given in Appendix C.

Figure 4. This figure summarizes the idea of Greene and Chance. The new pressure profile (black) with localized variation over the flux surface $\psi = \psi _0$ lies over the equilibrium profile (dashed red). Even though the variation in pressure at $\psi = \psi _0$ is small, the change in the pressure gradient can be large.

To obtain useful maximum growth rate scans it is computationally advantageous to know where the equilibrium transitions from being stable to unstable, i.e. the region of marginal stability. This is because stable modes are extended and require many more points and a wider range in $\theta$ than unstable modes, leading to a longer computation time. To that end, we first integrate the marginally stable ballooning equation

(3.7)\begin{equation} \frac{1}{\mathcal{J}}\frac{\partial}{\partial \theta}\left( \frac{\lvert \boldsymbol{\nabla}\alpha \rvert^2 }{\mathcal{J}\, B^2} \frac{\partial \hat{X}}{\partial \theta}\right) + 2 \frac{{\rm d} p}{{\rm d} \psi}\left[ \boldsymbol{B}\times\boldsymbol{\nabla}\left(\mu_0 p + \frac{B^2}{2}\right)\boldsymbol{\cdot} \boldsymbol{\nabla}{\alpha}\right] \hat{X} = 0, \end{equation}

along the field and count the zeros of the function $\hat {X}(\theta )$ – if $\hat {X}$ has at least one zero, the mode is unstable, else it is stable. This criterion was originally developed by Newcomb (Reference Newcomb1960) for a screw pinch. He used it as a method to asses the stability of a screw pinch without explicitly finding the growth rates or the eigenfunctions. It is briefly explained in Appendix D. Using this criterion, one can obtain the sign of $\gamma ^2$ and infer the stability significantly more rapidly than by exactly solving (3.3). Coupling Newcomb's criterion with the Greene–Chance analysis gives us the ability to scan the $\hat {s}-\alpha _{\mathrm {MHD}}$ space and plot the marginal stability contour ($\gamma = 0$) cheaply. For axisymmetric equilibria, the marginal stability contour is a single continuous line. Upon obtaining the contour, we choose a region around it where we solve (3.3).

To solve (3.3) we use the procedure described by Sanchez et al. (Reference Sanchez, Hirshman, Whitson and Ware2000). Our two-part codeFootnote ³ first finds the contour of marginal stability, then takes a region around the contour in the $\hat {s}-\alpha _{\mathrm {MHD}}$ space and implements the algorithm given in Sanchez et al. (Reference Sanchez, Hirshman, Whitson and Ware2000). It outputs the maximum eigenvalue and the corresponding eigenfunction for each value of $\hat {s}$ and $\alpha _{\mathrm {MHD}}$. The plots of the maximum eigenvalue along with the curve of marginal stability are shown in the next section.Footnote ⁴

3.3 Ideal-ballooning analysis results

This section contains the results of the $\hat {s}-\alpha _{\mathrm {MHD}}$ analyses of the twelve local equilibria that we chose in § 2.1. We plot a 2-D contour plot of the magnitude of the growth rate as defined by (3.4). We begin by discussing the positive-triangularity equilibria in figure 5.

Figure 5. This figure shows the normalized growth rate $\gamma a_{{N}}/v_{\mathrm {th}, i}$ contours along with the curve of marginal stability (white line) for the positive-triangularity equilibria. Columns correspond to the low, intermediate and high-$\beta$ regimes, respectively. The nominal equilibrium value is given by the green cross. The difference between the growth rates from the low and high-$\beta$ equilibria is due to our choice of normalization in (3.4).

All the positive-triangularity equilibria studied here are stable at their nominal values. The low-$\beta$ equilibria lie below the marginal stability contour whereas the high-$\beta$ equilibria lie above it. For low-$\beta$ equilibria the ballooning threshold is well known to be $\alpha _{\mathrm {MHD}} \sim 1$ but for the high-$\beta$ equilibria, using (3.6), $\alpha _{\mathrm {MHD}} \sim 1/\delta _{\mathrm {Hsu}}^2 \gg 1$ which pushes these equilibria into the region of ‘second’ stability, first discovered by Greene & Chance (Reference Greene and Chance1981). Note that the high-$\beta$ equilibrium in figure 5( f) is close to marginally stable. Also, since ballooning modes are Alfvén waves and we use ion-thermal speed in (3.4), we see that the maximum growth rate decreases as $v_{\mathrm {th},i}/v_{{A}} = \sqrt {\beta /2}$ –high-$\beta$ equilibria having the smallest maximum growth rates. This non-conventional choice of normalization will help us quantitatively compare ideal-ballooning growth rates with the ones obtained from various microstability studies in §§ 5–7. In the next paragraph, we discuss the results from the ideal-ballooning scans of the negative-triangularity equilibria.

In figure 6, we can see the nominal equilibria for the negative triangularity are stable for all cases except figure 6(e). The trends follow those of the positive-triangularity equilibria with one important exception. Unlike the positive-triangularity equilibria, the high-$\beta$ negative-triangularity equilibria move closer to the marginal stability boundary in the inner core region as seen in figure 6(c). This is different from the usual result – ballooning-stable equilibria with peaked pressure profiles approach marginal stability as we move towards the boundary and pressure gradient gets steeper (as seen in figures 5, 6a, 6b, 6d and 6e) – indicating that negative-triangularity, high-$\beta$ equilibria have enhanced ballooning stability in the outer core. We will observe a similar trend when we present results from the electromagnetic microstability analyses in § 7.

Figure 6. This figure shows the normalized growth rate contours along with the curve of marginal stability for the negative-triangularity equilibria. The nominal equilibrium value is denoted by the green cross; (a) $\rho = 0.5, \beta \sim 0.01$, (b) $\rho = 0.5, \beta \sim 0.1$, (c) $\rho = 0.5, \beta \sim 1$, (d) $\rho = 0.8, \beta \sim 0.01$, (e) $\rho = 0.8, \beta \sim 0.1$ and ( f) $\rho = 0.8, \beta \sim 1$.

Recently, Davies, Dickinson & Wilson (Reference Davies, Dickinson and Wilson2022) have published a study investigating access to high-$\beta$ spherical tokamak equilibria where they find that high-$\beta$, negative-triangularity equilibria are more unstable and less accessible than their positive-triangularity counterparts. At first, this may seem to contradict our results. However, looking at the plasma $\beta$ values in Davies et al., we realize that their high-$\beta$ equilibria correspond to intermediate-$\beta$ in our work. Indeed, the negative-triangularity intermediate-$\beta$ equilibrium in figure 6(e) is ballooning unstable whereas the positive-triangularity equilibrium in figure 5(e) is ballooning stable. Furthermore, Davies et al. define accessibility as the ability to reach the nominal $\alpha _{\mathrm {MHD}}$ from $\alpha _{\mathrm {MHD}} = 0$ at the nominal $\hat {s}$ – a straight line in $\hat {s}-\alpha _{\mathrm {MHD}}$ space – the actual path of an equilibrium from startup to a steady-state operation in the $\hat {s}-\alpha _{\mathrm {MHD}}$ space is a 2-D curve, similar to figures 2(a) and 2(b) in Chance et al. (Reference Chance, Jardin, Kessel, Manickam, Monticello, Peng, Holmes, Strickler, Whitson and Glasser1990). Therefore, our results corroborate rather than contradict the findings in Davies et al.

The most important takeaway from this study is that all low- and high-$\beta$ equilibria are stable to ideal-ballooning mode. This indicates that it might be possible to generate high-$\beta$ equilibria that are ballooning stable.Footnote ⁵ Even though one of the intermediate-$\beta$ equilibria is unstable, we will use it in our study as it will help us understand the behaviour of microstability with changing $\beta$.

4.1 The gyrokinetic model

The electromagnetic gyrokinetic model is a simplification of the $6D$ Vlasov–Maxwell system of equations to a $5D$ system that predicts the self-consistent evolution of a distribution of charged particles and its electromagnetic fields in the presence of low-frequency, small-scale fluctuations. We define the distribution and fields as

(4.1)\begin{equation} \left. \begin{array}{c@{}} f_s = F_{0s} + \delta f_s,\\ \boldsymbol{E} = \boldsymbol{E}_0 + \delta \boldsymbol{E},\\ \boldsymbol{B} = \boldsymbol{B}_0 + \delta \boldsymbol{B},\\ \end{array} \right\} \end{equation}

where the fields comprise their equilibrium components (with a subscript $0$) plus small fluctuations. The fluctuations $\delta f_s, \delta \boldsymbol {E}$ and $\delta \boldsymbol {B}$ are defined such that they vanish when averaged over length and time scales much larger than the particle gyroradius $\rho _s$ and turbulence frequency $\omega$, respectively. The gyrokinetic model is derived in the limit

(4.2)\begin{equation} \tilde{\epsilon} \equiv \frac{\omega}{\varOmega_{s}} \sim \frac{\rho_s}{a_{{N}}} \sim \frac{k_{{\parallel}}}{k_{{\perp}}} \sim \frac{\delta f_s}{F_{0s}} \sim \frac{Z_s e \varphi}{T_{s}} \sim \frac{|{\delta \boldsymbol{B}}|}{|\boldsymbol{B}_0|} \ll 1, \end{equation}

where $\varOmega _{s} = (Z_s e B)/(m_{s}c)$ is the cyclotron frequency and subscript $s$ denotes the species. The variable $\varphi$ is the perturbed electrostatic potential. The particle gyroradius $\rho _s$, given by

(4.3)\begin{equation} \rho_s \equiv \frac{v_{\mathrm{th},s}}{\varOmega_s}, \end{equation}

is the perpendicular length scale of the turbulent fluctuations and $v_{\mathrm {th},s} = \sqrt {(2 T_s)/m_s}$ is the thermal velocity. The length scale $a_{{N}}$ is the effective minor radius, defined in § 2.1. The wavenumbers of the mode perpendicular and parallel to the equilibrium magnetic field are denoted by $k_{\perp }$ and $k_{\parallel }$, respectively.

In the small $\tilde {\epsilon }$ limit, one can reduce the dimensionality of the problem from $6D (\boldsymbol {r}, w_{\perp }, w_{\parallel }, \vartheta )$ to $5D (\boldsymbol {r}, w_{\perp }, w_{\parallel })$ by averaging over the gyrophase $\vartheta$. For this $5D$ coordinate system, we will transform back and forth between two different coordinates, the particle position and velocity coordinates $(\boldsymbol {r}, w_{\perp }, w_{\parallel }, t)$ and the guiding-centre coordinates $(\boldsymbol {R}_s, E_s, \mu _s, t)$ where

(4.4)\begin{gather} E_s = \frac{1}{2} m_s w^2, \end{gather}

(4.5)\begin{gather}\mu_s = \frac{m_s w_{{\perp}}^2}{2B}, \end{gather}

are the total kinetic energy and the magnetic moment of the particle. The guiding centre is given in terms of the particle position by the Catto transformation (Catto & Tsang Reference Catto and Tsang1977)

(4.6)\begin{equation} \boldsymbol{R}_s = \boldsymbol{r} - \frac{\boldsymbol{b}\times \boldsymbol{w}_{{\perp}}}{\varOmega_s}. \end{equation}

The gyroaveraging operators $\langle \rangle _{\boldsymbol {R}_s}$ and $\langle \rangle _{\boldsymbol {r}}$

(4.7)\begin{gather} \langle X \rangle_{\boldsymbol{R}_s} = \frac{1}{2{\rm \pi}}\int_{0}^{2{\rm \pi}} X(\boldsymbol{r}) \, {\rm d} \vartheta, \end{gather}

(4.8)\begin{gather}\langle X \rangle_{\boldsymbol{r}} = \frac{1}{2{\rm \pi}}\int_{0}^{2{\rm \pi}} X(\boldsymbol{R}_s) \, {\rm d} \vartheta, \end{gather}

denote the average of $X$ over a gyration period at fixed guiding centre $\boldsymbol {R}_s$ and at fixed position $\boldsymbol {r}$, respectively. It is convenient to define the gyrokinetic model in terms of the parallel component $\delta A_{\parallel }$ of the magnetic vector potential, the magnetic field strength fluctuation $\delta B_{\parallel }$

(4.9)\begin{equation} \delta B_{{\parallel}} = \boldsymbol{b}\boldsymbol{\cdot} (\boldsymbol{\nabla} \times \delta \boldsymbol{A}_{{\perp}}), \end{equation}

the electrostatic potential

(4.10)\begin{equation} \delta \boldsymbol{E} ={-}\boldsymbol{\nabla}\varphi - \frac{1}{c}\frac{\partial \boldsymbol{A}}{\partial t}, \end{equation}

and the gyrokinetic distribution function in the guiding-centre coordinate system $(\boldsymbol {R}_s, E_s, \mu _s, t)$

(4.11)\begin{equation} h_s(\boldsymbol{R}_s, E_s, \mu_s, t) = \frac{Z_s e \varphi(\boldsymbol{r}, t) F_{0s}}{T_s} + \delta f_s(\boldsymbol{R}_s, E_s, \mu_s, t). \end{equation}

Using these new fields we can now introduce the $\delta f$ gyrokinetic theory that was first derived for the linear electromagnetic case by Antonsen & Lane (Reference Antonsen and Lane1980) and nonlinear case by Frieman & Chen (Reference Frieman and Chen1982). For a collisionless, linear electromagnetic model, following the notation of Abel et al. (Reference Abel, Plunk, Wang, Barnes, Cowley, Dorland and Schekochihin2013), the governing equations are

(4.12)\begin{gather} \frac{\partial h_s}{\partial t} + (w_{{\parallel}} \boldsymbol{b}+ \boldsymbol{v}_{Ds})\boldsymbol{\cdot} \frac{\partial h_s}{\partial \boldsymbol{R}_s} = \frac{Z_s e F_{0s}}{T_s}\frac{\partial \langle \varphi - \boldsymbol{w} \boldsymbol{\cdot} \delta \boldsymbol{A}/c \rangle_{\boldsymbol{R}_s}}{\partial t} - \boldsymbol{V}_{E} \boldsymbol{\cdot} \boldsymbol{\nabla} F_{0s}, \end{gather}

(4.13)\begin{gather}\sum_{s} \frac{(Z_s e)^2 \varphi}{T_s} = \sum_s Z_s e\int d^3\boldsymbol{w}\, \langle h_{s}\rangle_{\boldsymbol{r}}, \quad \tau = \frac{T_{{e}}}{T_{{i}}}, \end{gather}

(4.14)\begin{gather}-\boldsymbol{\nabla}_{{\perp}}^2 \delta A_{{\parallel}} = \frac{4{\rm \pi}}{c} \sum_s Z_s e \int d^3\boldsymbol{w}\, w_{{\parallel}} \langle h_s\rangle_{\boldsymbol{r}}, \end{gather}

(4.15)\begin{gather}\boldsymbol{\nabla}_{{\perp}}^2 \frac{\delta B_{{\parallel}}B}{4{\rm \pi}} ={-}\boldsymbol{\nabla}_{{\perp}} \boldsymbol{\nabla}_{{\perp}} \boldsymbol{:} \sum_s \int d^3\boldsymbol{w}\, \langle m_s \boldsymbol{w}_{{\perp}} \boldsymbol{w}_{{\perp}} h_s\rangle_{\boldsymbol{r}}, \end{gather}

where the velocity integrals in (4.13) are taken at fixed $\boldsymbol {r}$. The velocities $\boldsymbol {V}_E$ and $\boldsymbol {v}_{Ds}$ are the $\boldsymbol {E}\times \boldsymbol {B}$ and the magnetic drift velocities, respectively,

(4.16)\begin{gather} \boldsymbol{V}_{E} = \frac{c}{B} \boldsymbol{b}\times \langle \boldsymbol{\nabla} \varphi \rangle_{\boldsymbol{R}_s} - \frac{1}{B} \boldsymbol{b}\times \langle \boldsymbol{\nabla}( \boldsymbol{w}\boldsymbol{\cdot}\delta \boldsymbol{A}) \rangle_{\boldsymbol{R}_s}, \end{gather}

(4.17)\begin{gather}\boldsymbol{v}_{Ds} = \frac{w_{{\parallel}}^2}{\varOmega_s} \boldsymbol{b}\times (\boldsymbol{b}\boldsymbol{\cdot} \boldsymbol{\nabla}\boldsymbol{b}) + \frac{w_{{\perp}}^2}{2\,\varOmega_s} \frac{\boldsymbol{b}\times \boldsymbol{\nabla} B}{B}. \end{gather}

This completely defines the linear gyrokinetic system. The gyrokinetic model is a good approximation for the core region of a tokamak plasma where our study is being conducted (White et al. Reference White, Howard, Greenwald, Reinke, Sung, Baek, Barnes, Candy, Dominguez and Ernst2013; Creely et al. Reference Creely, Howard, Rodriguez-Fernandez, Cao, Hubbard, Hughes, Rice, White, Candy and Staebler2017). We solve the gyrokinetic model in a 3-D flux tube – a tube with a rhombus-shaped cross-section following the field line. The appropriate length of the flux tube and the boundary conditions at the ends are determined using the ideas developed by Beer, Cowley & Hammett (Reference Beer, Cowley and Hammett1995). All the variables are assumed to be periodic perpendicular to the field line. This allows us to further simplify (4.12) by writing the fluctuating fields as a Fourier series

(4.18)\begin{equation} \left. \begin{array}{c@{}} \displaystyle h_s = \sum_{k} h_{k_{{\perp}}, s}(\theta, E, \mu, t)\exp({\rm i} \boldsymbol{k}_{{\perp}}\boldsymbol{\cdot} \boldsymbol{R}_s ),\\ \displaystyle \varphi = \sum_{k} \varphi_{k_{{\perp}}}(\theta, t)\exp({\rm i} \boldsymbol{k}_{{\perp}}\boldsymbol{\cdot} \boldsymbol{r} ),\\ \displaystyle \delta A_{{\parallel}} = \sum_{k} \delta A_{{\parallel},k_{{\perp}}}(\theta, t)\exp({\rm i} \boldsymbol{k}_{{\perp}}\boldsymbol{\cdot} \boldsymbol{r} ),\\ \displaystyle \delta B_{{\parallel}} = \sum_{k} \delta B_{{\parallel}, k_{{\perp}}}(\theta, t)\exp({\rm i} \boldsymbol{k}_{{\perp}}\boldsymbol{\cdot} \boldsymbol{r} ). \end{array} \right\} \end{equation}

Applying this ansatz to (4.12), we obtain

(4.19)\begin{align} & \left(\frac{\partial}{\partial t} - {\rm i}\omega_{Ds}\right) h_{k_{{\perp}},s} +(\boldsymbol{b}\boldsymbol{\cdot}\boldsymbol{\nabla}\theta) w_{{\parallel}} \frac{\partial h_{k_{{\perp}},s}}{\partial \theta } = \left\{\frac{\partial}{\partial t} - i{\omega}_{*,s}\left[1 + \eta_s\left(\frac{E_s}{T_s} - \frac{3}{2}\right)\right]\right\} \nonumber\\ & \quad \times \left[ J_0\left(\frac{k_{{\perp}}w_{{\perp}}}{\varOmega_s}\right) \left(\varphi_{k_{{\perp}}} -\frac{w_{{\parallel}}\delta A_{{\parallel}}}{c}\right) + J_1\left(\frac{k_{{\perp}}w_{{\perp}}}{\varOmega_s}\right) \frac{w_{{\perp}}}{k_{{\perp}}} \frac{\delta B_{{\parallel}}}{c}\right] F_{0s}, \end{align}

where

(4.20)\begin{equation} \omega_{Ds} = \boldsymbol{k}_{{\perp}}\boldsymbol{\cdot} \boldsymbol{v}_{Ds}, \end{equation}

is the magnetic drift frequency, $J_0(k_{\perp }\rho _s)$ and $J_1(k_{\perp }\rho _s)$ are the zeroth- and first-order cylindrical Bessel functions, respectively,

(4.21)\begin{equation} \frac{a_{{N}}}{L_{Ts}} ={-} \frac{{\rm d} \log(T_s)}{{\rm d} \rho},\quad \frac{a_{{N}}}{L_{{n}_s}} ={-} \frac{{\rm d} \log(n_s)}{{\rm d} \rho}, \quad \eta_s = \frac{L_{{n}_s}}{L_{{T}_s}}, \end{equation}

and

(4.22)\begin{equation} \omega_{*,s} = \frac{T_s}{Z_s e B}[(\boldsymbol{b} \times \boldsymbol{k}_{{\perp}})\boldsymbol{\cdot} \boldsymbol{\nabla} \log{n_{s}}]. \end{equation}

For this study, we choose a hydrogen plasma, i.e. $Z_{{i}} = 1, Z_{{e}}=-1$. The variables $L_{{n}s}$ and $L_{{T}s}$ are the density and temperature-gradient scale lengths. The quantity $\omega _{*,s}$, the diamagnetic drift frequency, is the typical frequency of a drift wave – waves caused due to gradients in temperature or density. For these modes

(4.23)\begin{equation} \frac{\partial}{\partial t} \sim \omega \sim \omega_{*,s}. \end{equation}

In §§ 5 and 6, we will investigate two types of electrostatic drift wave instabilities, the ion-temperature-gradient (ITG) mode and the trapped electron mode (TEM). For purely electrostatic modes, we solve (4.12)–(4.13) only and assume the magnetic fluctuations, $\delta A_{\parallel }$ and $\delta B_{\parallel }$ to be absent. In § 7, we investigate the effect of magnetic fluctuations by solving the full set of (4.12)–(4.15). The procedure for numerically solving the gyrokinetic model is explained in the following section.

4.2 Using the GS2 code

GS2Footnote ⁶ (Kotschenreuther, Rewoldt & Tang Reference Kotschenreuther, Rewoldt and Tang1995b; Dorland et al. Reference Dorland, Jenko, Kotschenreuther and Rogers2000; Jenko & Dorland Reference Jenko and Dorland2001; Highcock Reference Highcock2012) is a parallel code that solves the gyrokinetic model as an initial-value problem. It solves (4.12)–(4.15) numerically by calculating the evolution of an initial perturbation inside a flux tube.

Before each run, one has to specify the value of the gradient scale lengths $a_{{N}}/L_{{n_{{s}}}}, a_{{N}}/L_{{T_{{s}}}}$, the range of normalized wavenumbers $k_{\perp } \rho _{{i}}$ and various geometric coefficients as a function of $\theta$. The values of these inputs are specific to the linear mode under consideration and will be provided in the following sections. Since we are only studying stability against fluctuations that vary on a small perpendicular scale $k_{\perp } \rho _{{i}} \sim 1$, we can get the geometric coefficients from the local equilibria just like we did for the ideal-ballooning stability analysis.Footnote ⁷

The perpendicular structure of different fluctuations is defined by defining the wavevector $\boldsymbol {k}_{\perp }$ that can be written as

(4.24)\begin{equation} \boldsymbol{k}_{{\perp}} = k_{y} \boldsymbol{\nabla} y + k_{x}\boldsymbol{\nabla} x, \end{equation}

where $x$ and $y$ are normalized forms of the coordinates $\psi$ and $\alpha$, respectively. For our study, we assume $k_{x} = 0$, i.e. modes with no variation in the radial direction at $\theta = 0$.Footnote ⁸ We choose around $15{-}25$ values of $k_y\rho _{{i}}$ in the range $k_y\rho _{{i}} = 0.05{-}6$. All of our simulations are well resolved in $\theta$ and well converged as can be seen in figure 7. For this study, $\theta \in [-19{\rm \pi}, 19{\rm \pi} ]$ and more than $450$ points along the $\theta$ grid unless stated otherwise.

Figure 7. This figure shows (a) the output from a typical electrostatic GS2 run showing the normalized frequency and growth rate spectrum with both electron-gradient- ($\omega a_{{N}}/v_{\mathrm {th}, i} < 0$) and ion-gradient-driven ($\omega a_{{N}}/v_{\mathrm {th}, i} > 0$) instabilities and (b) showing the variation along the field line of the square of the normalized electrostatic potential $\lvert \varphi \rvert ^2$. We can see that the potential is well resolved and decays sufficiently before reaching the boundaries.

For the velocity space structure GS2 uses an $(E, \lambda )$ grid instead of the $(w_{\parallel }, w_{\perp })$ grid. Defining the pitch angle as

(4.25)\begin{equation} \lambda(\theta_{{b}}) = \frac{\mu}{E}, \end{equation}

where $\theta _{{b}}$ is the bounce angle – the value of $\theta$ at which a trapped particle with a pitch angle $\lambda$ reflects back from a region of high magnetic field. For a given pitch angle $\lambda$, the bounce angle is defined such that $B(\theta _{{b}}) = 1/\lambda$. In GS2, resolution of the passing particle distribution function in the coordinate $\lambda$ is set by the variable $\mathrm {nlambda}$. We set $\mathrm {nlambda} = 12$. For the trapped particle distribution, we choose $11$ bounce points. Similarly, for the energy-space resolution, we set the value of the GS2 variable $\mathrm {negrid} = 10$.Footnote ⁹ We choose $27$ points along the flux tube for every $2{\rm \pi}$ interval to ensure sufficient resolution along the field line. This completely defines the resolution in GS2.

In figure 7, we show the results from a typical GS2 run. After each run, one obtains the normalized growth rate $\gamma a_{{N}}/v_{\mathrm {th}, i}$, the wave frequency $\omega a_{{N}}/v_{\mathrm {th}, i}$, the electrostatic potential eigenfunction $\varphi (\theta, t)$ for electrostatic runs and $\varphi (\theta, t), \delta A_{\parallel }(\theta, t)$ and $\delta B_{\parallel }(\theta, t)$ for electromagnetic runs. We also obtain the quasilinear particle and heat fluxes for each mode

(4.26)\begin{gather} \varGamma_{s, k_y} = \int \frac{{\rm d} \theta}{\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\theta} \int {\rm d}^3\boldsymbol{w} (\boldsymbol{V}_{E, k_y}\boldsymbol{\cdot}\boldsymbol{\nabla}\psi) h_{s, k_y}, \end{gather}

(4.27)\begin{gather}Q_{s, k_y} = \int \frac{{\rm d}\theta}{\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\theta} \int {\rm d}^3\boldsymbol{w} (\boldsymbol{V}_{E, k_y}\boldsymbol{\cdot}\boldsymbol{\nabla}\psi) E_s h_{s, k_y}, \end{gather}

where the subscript $k_y$ denotes the mode $k_y$ in the Fourier space. These quantities can be used to extract information about an unstable mode (Kotschenreuther et al. Reference Kotschenreuther, Liu, Hatch, Mahajan, Zheng, Diallo, Groebner, Hillesheim, Maggi and Giroud2019). Since this is a linear study, the absolute values of the fluxes do not contain any useful information but their ratio $\varGamma _{s, k_y}/Q_{s, k_y}$ can still be used to characterize the type of instability. We will use this ratio for the TEM study in § 6.

5.1 Details of the study

The most important form of electrostatic instability that arises at low wavenumbers, the ITG (Cowley, Kulsrud & Sudan Reference Cowley, Kulsrud and Sudan1991), occurs when a drift wave becomes unstable due to a large ion temperature gradient, i.e. large $a_{{N}}/L_{{T_i}}$. Therefore, our objective is to understand this mode by doing a scan in the temperature-gradient scale length, $a_{{N}}/L_{{T_i}}$. Using the definition of the pressure $p_s = n_s T_s$, we can write

(5.1)\begin{equation} \frac{a_{{N}}}{L_{p_s}} \equiv{-}\frac{{\rm d} \log(p)}{{\rm d} \rho} ={-}\frac{{\rm d} \log(T_s)}{{\rm d} \rho} - \frac{{\rm d} \log(n_s)}{{\rm d} \rho} = a_{{N}} \left(\frac{1}{L_{n_s}} + \frac{1}{L_{T_s}}\right).\end{equation}

Using the equation above along with (3.6), we can write

(5.2)\begin{equation} \alpha_{\mathrm{MHD}} = \frac{\beta}{2} \sum_{s} \frac{a_{{N}}}{L_{{p}_s}} \frac{\rho q^2}{2 \mu_0 \epsilon}. \end{equation}

Furthermore, recall that we can vary the normalized pressure gradient $\alpha _{\mathrm {MHD}}$ for a local equilibrium using the idea of Greene and Chance without recalculating the global equilibrium. This gives us the ability to self-consistently vary the temperature and density gradient scale lengths for a fixed $\beta$ as long as we recalculate the local equilibrium for the resulting value of $\alpha _{\mathrm {MHD}}$.Footnote ¹⁰ Table 3 contains the nominal density, temperature and pressure gradient scale lengths, denoted $a_{{N}}/L^{\mathrm {nom}}_{{n}_{{i}}}$, $a_{{N}}/L^{\mathrm {nom}}_{{T}_{{i}}}$ and $a_{{N}}/L^{\mathrm {nom}}_{p_{{i}}}$, respectively.

These are the values obtained from the original local equilibrium generated by VMEC and are exactly the same for all the different beta and triangularity values. For the ITG mode study, we define

(5.3)\begin{equation} \mathrm{fac} = \left.\frac{{\rm d} P}{{\rm d} \rho}\right / \left(\frac{{\rm d} P}{{\rm d} \rho}\right)_{\mathrm{nom}}, \end{equation}

as the ratio of actual pressure to the nominal pressure. We choose $\mathrm {fac} = (0.5, 1, 2, 4, 8)$ times the nominal pressure gradient for $\rho = 0.5$ and $\mathrm {fac} = (0.5, 1, 2, 4)$ times the nominal pressure gradient for $\rho = 0.8$. For each pressure gradient, we choose two density gradient scale lengths – the nominal and half of the nominal value from the local VMEC equilibrium while varying the temperature gradient scale length consistently for each gradient scale length. Tables 5 and 6 contain the resulting values.

Table 5. Values of $a/L_{{T_i}}$ at $\rho = 0.5$ for the ITG study.

Table 6. Values of $a/L_{{T_i}}$ at $\rho = 0.8$ for the ITG study.

These $18$ values of various scale lengths are exactly the same for all the triangularities as well as for all the different beta values due to the way we have defined $\rho$. From previous observations and studies, we know that the typical peak ITG growth rate lies around $k_y \rho _{{i}} = 1$. To capture the maximum growth rate, we calculate the growth rates in the range $k_y \rho _{{i}} \in [0.05, 2]$. For ITG, we have made the common assumption of adiabatic electrons to exclude the effect of kinetic electrons on the ITG mode and avoid other modes like the TEM. Mathematically, this means one assumes $h_{{e},k_y} = 0$ when solving (4.19) for the electrons.

Using the values in tables 5 and 6, we run GS2 in the electrostatic limit ($\delta A_{\parallel } = 0, \delta B_{\parallel } = 0$) and obtain the maximum normalized growth rate $\gamma a_{{N}}/v_{\mathrm {th},i}$ for each of the $108$ cases, $18$ for each beta and each triangularity value. The results showing comparison between different beta values, boundary shapes and normalized radii will be shown in § 5.3. We expect the equilibria to become more stable to the ITG mode as we increase $\beta$. This behaviour is well-known (Jarmén, Malinov & Nordman Reference Jarmén, Malinov and Nordman1998; Hirose Reference Hirose2000) in the literature for low and intermediate-$\beta$ equilibria. To try and explain this trend, in the next section we look as the local shear as a characteristic quantity that explains stabilization of the ITG mode with increasing $\beta$.

5.2 Characterizing stability to the ITG mode

In this section, we define and plot an important quantity, the local magnetic shear (Greene & Chance Reference Greene and Chance1981), which will help us understand the response of a local equilibrium to the ITG mode. We use the local shear since negative global shear $\hat {s}$ is known to stabilize ITG (Uchida et al. Reference Uchida, Sen, Fukuyama and McCarthy2003). In the following section, we will plot the local shear as a function of $\theta _{\mathrm {geo}}$ at the nominal $\hat {s}$ and $\alpha _{\mathrm {MHD}}$ for the high-$\beta$ equilibria and compare it with the low, intermediate and a low-$\beta$ shifted-circle equilibrium (abbreviated as SC in the plots). We will show that the behaviour of ITG is directly related to the local shear and argue that a large negative local shear over a wide range along the field lines stabilizes the ITG mode. Mathematically, the local shear $\nu$ is given by (Greene & Chance Reference Greene and Chance1981; Dewar & Glasser Reference Dewar and Glasser1983)

(5.4)\begin{equation} \nu ={-}\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla}_{{N}}\left(\frac{\boldsymbol{\nabla}_{{N}}\alpha\boldsymbol{\cdot} {\boldsymbol{\nabla}}_{{N}}\psi}{\lvert\boldsymbol{\nabla}_{N}\psi\rvert^2}\right),\end{equation}

where $\boldsymbol {\nabla }_{{N}} = a_{{N}} \boldsymbol {\nabla }$ is used to non-dimensionalize $\nu$. We will plot the local shear as a function of the geometric poloidal angle defined in the figure 8. The coordinate $\theta _{\mathrm {geo}}$ is advantageous as it, unlike $\theta$, is a physically intuitive poloidal angle. Mathematically, $\theta _{\mathrm {geo}}$ is a monotonic function of $\theta$ – $\nu$ can always be transformed to $\hat {\kappa }(\theta )$ and vice versa. To plot $\nu$ with respect to $\theta _{\mathrm {geo}}$, we further simplify (5.4)

(5.5)\begin{equation} \nu ={-}\frac{F}{q R^2} \left[\left(\frac{q F' }{F} + F{'} F \frac{q}{(R_s B_{{{p}}s})^2}\right) + \frac{q\, p{'}}{ B_{{{p}}s}^2} + \frac{2 q}{ R_s B_{{{p}}s}} \left(\frac{\sin(u)}{R_s} - \frac{1}{R_{{c}}} \right)\right]. \end{equation}

The symbols used and formalism needed to derive (5.5) is given in Appendix C. It is important to point out that

(5.6)\begin{equation} \hat{s} = \frac{\rho |\boldsymbol{\nabla}_{{N}} \psi|}{2{\rm \pi} q a^2_{{N}}} \int_{0}^{2{\rm \pi}} \frac{{\rm d}\theta\, {\rm d}\phi}{\boldsymbol{B}\boldsymbol{\cdot}{\boldsymbol{\nabla}}_{{N}}\theta}\nu. \end{equation}

Thus, the appropriately weighted average of the local shear over a flux surface gives us the global magnetic shear $\hat {s}$. The global shear is held fixed for all the local equilibria at a given $\rho$. This relation implies that the local shear can be negative for a given range in $\theta$ for a positive global shear $\hat {s}$. Finally, we plot local shear at nominal $\hat {s}$ and $\alpha _{\mathrm {MHD}}$ in figure 9 for different beta values.

Figure 8. This figure illustrates the definition of $\theta _{\mathrm {geo}}$. The coordinates of the magnetic axis, marked with a cross, are $(R_{\mathrm {ax}}, Z_{\mathrm {ax}})$.

Figure 9. This figure explains the physical meaning of the local magnetic shear with (a) showing a typical local shear plot of negative triangularity equilibria at the nominal $\hat {s} = 0.45$ and nominal $\alpha _{\mathrm {MHD}}$ values at $\rho = 0.8$ for different beta values and a low-$\beta$ shifted-circle model (abbreviated SC). (b) Illustrates a modified interpretation from Antonsen et al. (Reference Antonsen, Drake, Guzdar, Hassam, Lau, Liu and Novakovskii1996) explaining the concept of local magnetic shear. Negative local shear twists the turbulent eddies more than positive or zero local shear in the region of bad curvature, stabilizing the ITG mode.

The local shear depends on the pressure gradient $\alpha _{\mathrm {MHD}}$ which further depends on $\beta$ as well as the gradient scale lengths $L_{{n}}$ and $L_{{T}}$. The beta value increases the local negative shear through the Shafranov shift. The gradient scale length does so by increasing the poloidal current gradient $\textrm {d} F/\textrm {d} \rho$ required to balance the pressure gradient which consequently increases the toroidal magnetic field. A plot showing the effect of pressure gradient on the local shear is shown in figure 13(a). Note that in our scans decreasing the pressure gradient scale length increases the driving term $\eta = L_{{n}}/L_{{T}}$. In fact, in this study, for a given $\beta$, $\eta$ increases faster than $a_{{N}}/L_{{p}}$. A plot of $\eta$ versus the pressure gradient scaling factor is shown in figure 10.

Figure 10. This shows (a) the change in the local shear for $\mathrm {fac} = 4$ (increased pressure gradient). The local shear is much more negative with an approximately linear dependence with $\mathrm {fac}$ on outboard side. On the other hand, (b) shows the comparison between the ITG driving term $\eta$ and the pressure gradient scaling factor $\mathrm {fac}$. The term $\eta$ is calculated using the values given in table 5. We can see that $\eta$ grows linearly, but with a larger pre-factor. These figures illustrate how the ITG driving term grows more rapidly than the stabilizing local shear as we increase the pressure gradient.

The dominant mechanism for generating negative local shear is the $\beta$-induced Shafranov shift. The local shear may also depend on the shaping, especially for negative triangularities. However, we find that shaping does not have a significant effect on the ITG stability of the high-$\beta$ equilibria.

5.3 The ITG results

Since GS2 is an initial-value code, resolving a growth rate $\gamma \leq \gamma _{\mathrm {thresh}}$ would require us to run the simulation for at least $t \sim 1/\gamma _{\mathrm {thresh}}$. For large wavenumbers and a semi-implicit time-stepping scheme, it leads to a relatively high runtime cost. Hence, we define a small positive ‘threshold’ growth rate $\gamma _{\mathrm {thresh}}$ which we will use to separate stable and unstable modes. If the maximum ITG growth $\max (\gamma a_{{N}}/v_{\textrm {th}, i}) < \gamma _{\mathrm {thresh}}$, we classify it as stable. For this study, we choose $\gamma _{\mathrm {thresh}} = 0.005$.

We find that all the nominal equilibria stabilize as we increase $\beta$ as first noted by Candy & Waltz (Reference Candy and Waltz2003a) – the high-$\beta$ equilibria are stable to the ITG mode, as shown in figure 11 below. To better understand what causes this effect, we do a scan in gradient scale length scales given in tables 5 and 6. Each group of plots contains maximum growth rate for a range of $a/L_{{T}_{i}}$ at a fixed $a/L_{{n}_{i}}$ for the three beta values. Plots are grouped by the triangularity of the equilibria. For each group, there are subgroups based on the normalized radius $\rho$. The results corresponding to the positive-triangularity boundary shape are shown in figure 12.

Figure 11. These plots show the ITG $\max (\gamma a_{{N}}/v_{\mathrm {th}, i})$ (over $k_y \rho _{{i}} \in [0.05, 2]$) vs the typical $\beta$ for nominal equilibria at different radial locations.

Figure 12. This figure shows the ITG $\max (\gamma a_{{N}}/v_{\mathrm {th}, i})$ plots for positive-triangularity equilibria against the temperature gradient length scale. For the high-$\beta$ equilibria, ITG is stabilized at both $\rho = 0.5$ and $\rho = 0.8$. The rightmost figures in each row are the local magnetic shear vs the geometric theta $\theta _{\mathrm {geo}}$ at the nominal $\textrm {d} p/\textrm {d} \rho$ and $\hat {s}$. The grey line corresponds to the local shear for a low-$\beta$ shifted-circle equilibrium (abbreviated SC). The magnetic shear $\hat {s}$ is the same for all the equilibria at every $\rho$.

In figure 11, one observes that, for the high-$\beta$ cases, the ITG mode is stable (that is, $\gamma < \gamma _\mathrm {thresh})$. For the intermediate and low-$\beta$ cases, figure 12 shows destabilization with increasing temperature gradient. We believe that the stabilization of high-$\beta$ equilibria is a result of large local negative shear (rightmost panels) that spans over a wide range in $\theta$. The local shear becomes positive only after $\theta _{\mathrm {geo}} > {\rm \pi}/2$ – the whole outboard side has a large local negative shear. These large negative values are predominantly due to a large Shafranov shift but could have sub-dominant effects resulting from strong shaping, especially for the negative triangularity high-$\beta$ equilibria. For the shifted-circle equilibria, the local shear is small as compared with the rest of the equilibria even though it is negative over the outboard side. This is because the low-$\beta$ shifted-circle equilibria are neither strongly shaped nor have a large beta. The trend remains the same even for the equilibria with half-nominal density gradients. Next, we plot results for negative triangularity equilibria in figure 13.

Figure 13. Shows the ITG $\max (\gamma a_{{N}}/v_{\mathrm {th}, i})$ plots for negative-triangularity equilibria. For the high-$\beta$ equilibria, ITG is stabilized at both $\rho = 0.5$ and $\rho = 0.8$. The local shear for the high-$\beta$ equilibria is negative over the whole outboard side. The grey line corresponds to the local shear for a low-$\beta$ shifted-circle equilibrium (abbreviated SC).

We see that the high-$\beta$ equilibria with a negative-triangularity boundary shape are stabilized as well. The local magnetic shear shows the same trend with beta but is even more negative and spans an even wider range in $\theta _{\mathrm {geo}}$ as compared with the positive-triangularity equilibria. We believe this extended range in negative local shear is due to the ‘squareness’ of the high-$\beta$ profiles, which is a stronger shaping than the positive triangularity cases. Since all the high-$\beta$ equilibria are stabilized and the growth rates are small, it is hard to find a clear difference in growth rates based on triangularity.

In the next section, we will study the stability of different equilibria to an electron-driven electrostatic mode, the TEM.

6.1 Details of the study

The second type of electrostatic instability that we investigate, the collisionless TEM, becomes unstable when drift waves resonate with the precession of the electrons. This can cause significant transport loss through the electron channel that degrades plasma confinement (Adam, Tang & Rutherford Reference Adam, Tang and Rutherford1976). For the TEM, we choose five values of pressure gradient, corresponding to $\mathrm {fac} = (0.5, 1, 2, 4, 8)$ for $\rho = 0.5$ and five values corresponding to $\mathrm {fac} = (0.5, 0.75, 1, 2, 4)$ for $\rho = 0.8$. For each pressure gradient, we choose two temperature-gradient scale lengths – nominal and $30\,\%$ of the nominal, while scanning the growth rates in the density gradient scale lengths. We do this since TEM, unlike ITG, is primarily a density gradient-driven instability. Tables 7 and 8 contain the values for which we have solved the gyrokinetic model. Note that the temperature and density scale lengths are the same for the ions and electrons.

Table 7. Values of gradient scale lengths at $\rho = 0.5$ used for the TEM study.

Table 8. Values of gradient scale lengths at $\rho = 0.8$ used for the TEM study.

These $20$ values of gradients are the same for ions and electrons as well as both positive- and negative-triangularity equilibria at all the different beta values. Unlike the ITG study, we turn on the kinetic effects of electrons since TEM is an electron-driven instability. The TEM growth rate peak occurs over a wide range $k_y\rho _{{i}} \in [0.5, 6]$. Since there is an overlap with the ITG and the electron-temperature-gradient (ETG) mode, having two species makes it difficult to separate modes with purely ITG and ETG-related effects from modes with purely TEM-related effects. Therefore, to calculate the TEM growth rate, we choose the growth rate corresponding to the wavenumber $k_y\rho _{i}$ at which the ratio of the quasilinear electron flux to electron heat flux is the maximum, i.e.

(6.1)\begin{equation} \gamma_{\mathrm{TEM}}(k_y\rho_{{i}}) = \gamma\Big \lvert_{\max (\varGamma_{{e}, ky}/Q_{{e}, ky})}, \end{equation}

where the definition of quasilinear fluxes is given in equations (4.26) and (4.27).

We run GS2 for the wavenumbers in the range $k_y\rho _{{i}} \in [0.2, 6.5]$. Just like the ITG study, we run GS2 and obtain the maximum growth rate $\gamma a_{{N}}/v_{\mathrm {th}, i}$ for each of the $120$ cases: $10$ values of $(a_{{N}}/L_{{T}}, a_{N}/L_{{n}}) \times 2$ surfaces $\times 3 \, \beta$ values $\times 2$ boundary shapes. The results showing the comparison between different beta values, triangularities, and normalized radius will be shown in § 6.3.

Curvature-driven TEMs are associated with the precession of trapped electrons in a flux surface. To that end, we elucidate the definition and role of the electron precession drift frequency in the next section and how it characterizes the stability of an equilibrium to the TEM.

6.2 Characterizing stability to the TEM

The collisionless curvature-driven TEM is a drift wave that becomes unstable when it resonates with the bounce precession of trapped electrons. The precession of electrons is characterized by their precession frequency

(6.2)\begin{equation} \langle\omega_{\mathrm{De}}\rangle = \left(\int_{-\theta_b}^{\theta_{{b}}} \dfrac{{\rm d} \theta}{w_{{\parallel}}} \frac{B}{(\boldsymbol{B}\boldsymbol{\cdot} {\boldsymbol{\nabla}}\theta)}\right)^{{-}1} \int_{-\theta_b}^{\theta_{{b}}} \frac{{\rm d} \theta}{w_{{\parallel}}} \dfrac{B}{(\boldsymbol{B}\boldsymbol{\cdot} {\boldsymbol{\nabla}}\theta)} \omega_{\mathrm{De}}, \end{equation}

where the integral operator is the bounce-average operator. The precession frequency is a function of the bounce angle $\theta _{{b}}$. Depending on the convention, one usually takes $\mathrm {sign}(\omega ) = \mathrm {sign}(\omega _{*,{e}})$ for trapped electron modes since the TEM is a drift wave. Therefore, if

(6.3)\begin{equation} \mathrm{sign}(\langle \omega_{\mathrm{De}} \rangle) \mathrm{sign}(\omega_{*,{e}}) < 0,\end{equation}

the drift wave will not be able to resonate with the precession of electrons. If the precession drift satisfies (6.3) at all the different pitch angles, the curvature-driven TEM will be stabilized (Connor, Hastie & Martin Reference Connor, Hastie and Martin1983; Roach, Connor & Janjua Reference Roach, Connor and Janjua1995). The expression for $\omega _{{D}}$, given by (4.20) can be alternatively written as

(6.4)\begin{equation} \omega_{\mathrm{De}} = \frac{k_y \rho_{{e}}}{2} \frac{v_{{\rm th},e}}{a_{{N}}} E_{{e}} [2(1-\lambda B) \omega_{\kappa} + \lambda B \omega_B ], \end{equation}

where $\omega _{\kappa }(\theta )$ and $\omega _B(\theta )$ are geometric factors independent of the electron energy $E_{{e}}$ and the pitch angle $\lambda$. A semi-analytical formula for calculating $\omega _{\kappa }$ is given in (C29). As a characteristic of TEM stability, we define the quantity

(6.5)\begin{equation} \langle \bar{\omega}_{\mathrm{De}}\rangle = \langle \omega_{{\rm De}} \rangle \mathrm{sign}(\omega_{*,{e}})/E_{{e}}, \end{equation}

as the precession drift per energy in the electron diamagnetic direction. A typical plot of the $\langle \bar {\omega }_{\textrm {De}} \rangle$ is shown at nominal $\hat {s}$ and $\alpha _{\mathrm {MHD}}$ for different $\beta$ values in figure 14. We see that the precession drift is negative everywhere only for high-$\beta$ case. This will form the basis for our understanding of the growth rate trends in the following section.

Figure 14. This figure shows the precession drift in (a) and the corresponding magnetic field magnitude in (b) for negative-triangularity equilibria at $\hat {s} = 0.45$ and nominal $\alpha _{\mathrm {MHD}}$ values at $\rho =0.8$ for different beta values. Note the atypical magnetic field for the high-$\beta$ equilibria where $\min (B)$ is located at a finite $\theta$.

6.3 The TEM results

Just like the ITG study, we choose a value of a growth rate $\gamma _{\mathrm {thresh}}$ such that if $\gamma _{\mathrm {TEM}} (a_{{N}}/v_{\mathrm {th}, i}) < \gamma _{\mathrm {thresh}}$, we classify an equilibrium as stable. The reasoning behind setting a threshold is described in the first paragraph of § 5.3. For this study, we choose $\gamma _{\mathrm {thresh}} = 0.005$. First, we plot the maximum TEM growth rates for the nominal equilibria in figure 15. We find that increasing the beta stabilizes the TEM and the high-$\beta$ equilibria are stable to the TEM.

Figure 15. This figure shows the TEM $\gamma _{\mathrm {TEM}} (a_{{N}}/v_{\mathrm {th}, i})$ vs the typical $\beta$ plots for nominal equilibria at different radial locations.

To understand this effect further, we plot the result from scans in the density gradients in two groups, each group containing the maximum TEM growth rate for a range of $a/L_{{n}}$ at a fixed $a/L_{{T}}$, given in tables 7 and 8 for different equilibria. Just like the ITG study, we group the plots by triangularity and arrange them in rows based on the normalized radius $\rho$. For positive triangularity equilibria, the results are shown in figure 16

Figure 16. This figure presents $\gamma _{\mathrm {TEM}} (a_{{N}}/v_{\mathrm {th}, i})$ plots for positive-triangularity equilibria. For the high-$\beta$ equilibria, TEM is stabilized at both $\rho = 0.5$ and $\rho = 0.8$. As you can see in figures 19(c) and 19( f), the precession drift for the high-$\beta$ equilibria is negative for all values of the bounce angle $\theta _{{b}}$.

We observe that the TEM is completely suppressed for the high-$\beta$, positive-triangularity equilibria. The frequency $\langle \bar {\omega }_{\mathrm {De}}\rangle$ is negative everywhere which means that all the trapped electrons precess in a direction opposite to the electron diamagnetic direction. They cannot destabilize the drift wave by exchanging energy with them.

The TEMs we seek are curvature driven. A large pressure gradient causes the electron precession drift to become negative for all pitch angles (Roach et al. Reference Roach, Connor and Janjua1995). This suppresses the curvature-driven mode and causes the slab-like branch of the TEM called the universal mode to appear. We can see that the universal mode dominates in figures 16(b) and 17(b) for the intermediate-$\beta$ equilibria at large density gradients. However, for high-$\beta$ equilibria, the universal mode is suppressed as well since the large local shear

(6.6)\begin{equation} L_{\nu} = \frac{a_{{N}}}{\nu}, \end{equation}

combined with strong shaping reduces the shearing length scale (Landreman, Antonsen & Dorland Reference Landreman, Antonsen and Dorland2015) which makes it harder for fluctuations to grow and persist along the field line.

Figure 17. This figure presents $\gamma _{\mathrm {TEM}} (a_{{N}}/v_{\mathrm {th}, i})$ plots for negative-triangularity equilibria. For the high-$\beta$ equilibria, TEM is stabilized at both $\rho = 0.5$ and $\rho = 0.8$. The rightmost plot on each row is the electron precession drift frequency.

Next, now look at the TEM growth rate trends for the negative triangularity equilibria, shown in figure 17. The negative-triangularity TEM growth rates follow the same trend as the positive-triangularity ones. The TEM is suppressed for the high-beta equilibria due to the negative precession drift. The intermediate and low beta are more unstable for thenegative-triangularity equilibria at $\rho =0.5$ and as unstable as the positive-triangularity ones at $\rho =0.8$.

The result of our local analyses are only strictly valid for perturbations localized to a field line on a flux surface. Alternatively, one can say that for all the fluctuations, the toroidal mode number $n \gg 1$. Therefore, we must ensure that all modes in this study meet the local approximation in order to ensure that they are self-consistent. In both § 5 and § 6, the lowest wavenumber in our analyses is $k_y \rho _{{i}} = 0.05$. The largest ion gyroradius arises in the outer-core high-$\beta$ equilibria where $\rho _{{i}} \approx 0.005\ \textrm {m}$. This corresponds to a wavelength $\lambda = 0.628\ \textrm {m}$. Assuming an $n = 1$ mode has a wavelength equal to the normalized minor radius $a_{{N}} = 0.68\ \textrm {m}$, $k_y \rho _{{i}} = 0.05$ corresponds to $n \approx 1$ – the longest modes are not localized. Hence, to accurately capture these modes at low wavenumbers, we must include non-local effects. However, for all the electrostatic studies, the growth rates $\gamma a_{{N}}/v_{\mathrm {th}, i} \rightarrow 0$ as $k_y \rho _{{i}} \rightarrow 0$ and the most unstable modes always arise at $k_{y, \textrm {peak}} \rho _{{i}} \geq 0.5$ which corresponds to $n \geq 10$ implying that the local approximation is a fair assumption since we are only concerned with peak growth rates in §§ 5 and 6.

We have demonstrated the stability of high-$\beta$ equilibria against two major sources of electrostatic instabilities. However, when $\beta \sim 1$, magnetic fluctuations may play an important role in deciding the stability of an equilibrium. Therefore, in the next section, we study the effect of electromagnetic modes on the high-$\beta$ equilibria.