Microstability of β ~ 1 tokamak equilibria

Abstract

High-power-density tokamaks offer a potential solution to design cost-effective fusion devices. One way to achieve high power density is to operate at a high $\beta$ value (the ratio of thermal to magnetic pressure), i.e. $\beta \sim 1$. However, a $\beta \sim 1$ state may be unstable to various pressure- and current-driven instabilities or have unfavourable microstability properties. To explore these possibilities, we generate $\beta \sim 1$ equilibria and investigate their stability. First, we demonstrate the generation of high-$\beta$ equilibria with the computer code VMEC. We then analyse these equilibria to determine their stability against the infinite-$n$ ideal-ballooning mode. We follow that by engaging in a detailed microstability study using the GS2 code, beginning with assessments of electrostatic ion-temperature-gradient and trapped election mode instabilities. We observe interesting behaviour for the high-$\beta$ equilibria – stabilization of these modes through two distinct mechanisms – large negative local shear and reversal of electron precession drift. Finally, we perform electromagnetic gyrokinetic simulations and observe enhanced stability in the outer core of high-$\beta$ equilibria and absence of kinetic ballooning modes in the negative-triangularity, high-$\beta$ equilibria. The enhanced outer-core stability of high-$\beta$ equilibria is different from their lower-$\beta$ counterparts and offers an alternative, potentially favourable regime of tokamak operation.

1 Introduction

The most advanced fusion reactor designs are currently based on the tokamak concept. These tokamak designs have historically led to large volume, large capital cost, plants such as ITER (Shimada et al. 2007) and EU-DEMO (Federici et al. 2014). Improvements in these designs could be achieved by operating at high power density, with reduced plasma volume (Menard et al. 2022).

For a fixed magnetic field strength, the power density of a tokamak $P$ scales as $\beta ^2$, where $\beta$ is the ratio of the plasma pressure to the magnetic pressure. Present-day tokamaks are low-$\beta$ devices. The achievable $\beta$ is typically limited by plasma instabilities. These can lead to disruptions or large turbulent transport. The higher $\beta$ is, the higher the pressure and current are, and therefore the larger the free energy available to drive these instabilities is. If these problems could be overcome, high-$\beta$ operation could be an attractive choice for future high-power-density (Menard et al. 2022) devices.

The high-beta, $\beta \sim 1$, regime has been explored previously in the context of asymptotic magnetohydrodynamic (MHD) equilibria by solving the Grad–Shafranov equation in the limit $\epsilon /(\beta q^2) \ll 1$ (Hsu, Artun & Cowley 1996) where $\epsilon$ is the aspect ratio of the tokamak and q is the safety factor defined in (??):. There have also been experimental explorations (Sykes et al. 2000; Gates & NSTX National Research Team 2003) of high-$\beta$ operation of the MAST and START tokamaks. There have only been a few studies (Chance et al. 1990; Hurricane, Chandran & Cowley 2000) that investigated the process of accessing these states while maintaining ideal-MHD stability; even fewer studies that study the microstability properties of $\beta \sim 1$ equilibria in detail (Wilson et al. 2004). Therefore, a detailed numerical analysis of these types of equilibria is required.

To that end, we generate a set of high-$\beta$ equilibria and study their susceptibility to local MHD and gyrokinetic instabilities. In context of local MHD, we study stability of these equilibria to the infinite-$n$, ideal-ballooning mode. Then we perform linear gyrokinetic analyses against various electrostatic and electromagnetic modes of instability in these equilibria. These modes are known to cause significant heat and particle transport in existing devices (White et al. 2013; Creely et al. 2017)

The remainder of this paper is divided as follows: in § 2, we briefly describe the fundamentals of an axisymmetric equilibrium. In § 2.1, we obtain numerical equilibria using the VMEC code (Hirshman & Whitson 1983). We describe and plot the equilibria for three different beta values: low-$\beta$ ($\beta \sim 0.01$), intermediate-$\beta$ ($\beta \sim 0.1$) and high-$\beta$ ($\beta \sim 1$) and two different triangularity values. In § 3, we introduce the technique devised by Greene & Chance (1981) to vary gradients of a local equilibria. We then use this to analyse the susceptibility of the chosen local equilibria to the ideal-ballooning instability and explain our observations. After ensuring ideal-ballooning stability for our high-$\beta$ equilibria, in § 4, we briefly explain the linear, collisionless, gyrokinetic model and process for analysing kinetic stability of local equilibria using the gyrokinetic code GS2. In § 5, we solve the gyrokinetic model for electrostatic fluctuations with adiabatic electrons and use the same techniques that we used in § 3 to scan the growth rates of all the local equilibria with respect to the temperature gradient. In § 6, we solve the gyrokinetic model for electrostatic fluctuations after relaxing the assumption of adiabatic electrons and use the same techniques as § 3 to scan the growth rates of all the local equilibria with respect to the density gradient. In § 7, we study the stability of local equilibria at their nominal gradients to electromagnetic fluctuations by solving the full gyrokinetic model and explain the advantages of negative-triangularity, high-$\beta$ equilibria. Finally in § 8, we summarize our work and discuss the directions in which it can be extended.

2 Generating axisymmetric $\beta \sim 1$ equilibria

In this section, we will start with the general form of a divergence-free magnetic field, simplify it for an axisymmetric configuration and using the simplified form, obtain an explicit form of the steady-state momentum equation – the Grad–Shafranov equation. In § 2.1, we briefly explain how VMEC works and provide details for the equilibria generated for this work.

A divergence-free magnetic field $\boldsymbol {B}$ can be written in the Clebsch form (D'haeseleer et al. 2012)

(2.1)\begin{equation} \boldsymbol{B} = \boldsymbol{\nabla}\alpha \times \boldsymbol{\nabla} \psi.\end{equation}

We shall restrict our attention to solutions whose magnetic field lines lie on closed nested toroidal surfaces, known as flux surfaces. We label the flux surfaces by their enclosed poloidal flux $\psi$. On each flux surface, the line of constant $\alpha$, the field-line label, gives us the path of the magnetic field line.

We will use the right-handed, cylindrical coordinate system $(R, \phi, Z)$ where $R$ and $Z$ are the radial and vertical distance from the origin and $\phi$ is the azimuthal angle about the symmetry axis. We also define a curvilinear coordinate system $(\psi, \phi, \theta )$ with $\psi$ being the flux surface label, $\phi$ being the cylindrical azimuthal angle and $\theta$ being the ‘straight-field-line’ poloidal angle (D'haeseleer et al. 2012) such that $\alpha = \phi - q(\psi ) \theta$ where

(2.2)\begin{equation} q(\psi) \equiv \frac{{\rm d} \chi}{{\rm d} \psi} = \frac{1}{2{\rm \pi}} \oint\, {\rm d} \theta \frac{\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla}\phi}{\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla} \theta} = \frac{\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\phi}{\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla} \theta}, \end{equation}

is the safety factor, $\chi$ is the enclosed toroidal flux, and the line integral is over $\theta \in [0, 2{\rm \pi} ]$. The relation for $q(\psi )$ is consistent with the definition of $\alpha = \phi - q\theta$ as can be seen by substituting $\alpha$ into (2.1), obtaining

(2.3)\begin{equation} \boldsymbol{B} = \boldsymbol{\nabla}\phi \times \boldsymbol{\nabla}\psi - q \boldsymbol{\nabla} \theta \times \boldsymbol{\nabla}\psi, \end{equation}

and noting that (2.3) satisfies (2.2). For an axisymmetric $\boldsymbol {B}$, (2.3) can be reduced further, to

(2.4)\begin{equation} \boldsymbol{B} = \boldsymbol{\nabla}\phi \times \boldsymbol{\nabla}\psi + F \boldsymbol{\nabla}\phi. \end{equation}

Here, $F = F(\psi, \theta )$. To generate an MHD equilibrium one has to solve the steady-state, ideal-MHD force balance equation (Freidberg 2014) which, in SI units, is

(2.5)\begin{equation} \boldsymbol{\nabla} p = \frac{(\boldsymbol{\nabla}\times \boldsymbol{B})\times \boldsymbol{B}}{\mu_0} = \frac{\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{B}}{\mu_0}-\boldsymbol{\nabla}\left(\frac{B^2}{2 \mu_0}\right), \end{equation}

where p is the plasma pressure and $\mu_0$ is the vacuum magnetic permeability. For an axisymmetric system, we can substitute the form of $\boldsymbol {B}$ from (2.4) and choose the toroidal component of (2.5) to get $\boldsymbol {B}\boldsymbol {\cdot } \boldsymbol {\nabla } F = 0$ which implies $F = F(\psi )$. Upon choosing the radial component, we obtain the Grad–Shafranov (Shafranov 1957; Grad & Rubin 1958) equation

(2.6)\begin{equation} \boldsymbol{\varDelta}^{*}\psi \equiv R^2 \boldsymbol{\nabla}\boldsymbol{\cdot} \left(\frac{\boldsymbol{\nabla} \psi}{R^2}\right) ={-}\mu_0 R^2 \frac{{\rm d} p}{{\rm d} \psi} - F\frac{{\rm d} F}{{\rm d} \psi}. \end{equation}

The Grad–Shafranov equation is a nonlinear equation for the poloidal flux $\psi (R, Z)$ that depends on the pressure $p(\psi )$ and current $F(\psi )$ profiles. As discussed in the introduction, for a fixed field strength $B$, the fusion power output $P$ of a tokamak scales as $\beta ^2$ – $\beta \equiv p/(B^2/(2\mu _0))$ is the ratio of plasma pressure to magnetic pressure – which makes high-$\beta$ operation an attractive concept. Therefore, our objective is to analyse the stability of equilibria that have $\beta \sim 1$.

To create high-$\beta$ equilibria, one may start with analytical solutions of the Grad–Shafranov equation. The most general analytical solution in the $\beta \sim 1$ limit was obtained by Hsu, Artun and Cowley (Hsu et al. 1996). However, these analytical $\beta \sim 1$ equilibria are unfit for our study as the geometric quantities required for a local stability analysis can be discontinuous and deviate significantly from the exact numerical solution. We briefly explain the analytical procedure used to calculate such equilibria and their limitations in Appendix A. To avoid these issues, we will use the equilibrium solver VMEC to generate high-$\beta$, axisymmetric equilibria in the next section.

2.1 Numerical equilibria

We generate equilibria numerically using the three-dimensional (3-D) equilibrium code VMEC (Hirshman & Whitson 1983). 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 1958). 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 2013) 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. 1998)

(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.¹ 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 Infinite-$n$ ideal-ballooning stability

In this section, we will investigate the equilibria generated in the previous sections for their stability to ideal-ballooning modes. In § 3.1, we describe the physics basis and mathematical formulation of the ideal-ballooning problem. In § 3.2, we review the tools for locally varying equilibria and introduce the concept of an $\hat {s}-\alpha _{\mathrm {MHD}}$ analysis. In the final section, we present the results from our study and discuss their implications.

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 1979) – 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. 1958) 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 parameter² 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. 1979; Dewar & Glasser 1983)

(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 (1981) and has been used extensively to study local stability of different axisymmetric equilibria (Connor, Hastie & Taylor 1978; Connor et al. 1979; Bishop 1985; Miller et al. 1998). 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 (1960) 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. (2000). Our two-part code³ 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. (2000). 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.⁴

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 (1981). 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 (2022) 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. (1990). 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.⁵ 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 Microstability analysis

This section contains the general theoretical and numerical details of our microstability analysis. In § 4.1, we will explain the physical basis and theoretical details of the gyrokinetic model. In the next section, we will explain how the model is implemented numerically using the GS2 code and provide the general details of our numerical study.

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 1977)

(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 (1980) and nonlinear case by Frieman & Chen (1982). For a collisionless, linear electromagnetic model, following the notation of Abel et al. (2013), 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. 2013; Creely et al. 2017). 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 (1995). 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

GS2⁶ (Kotschenreuther, Rewoldt & Tang 1995b; Dorland et al. 2000; Jenko & Dorland 2001; Highcock 2012) 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.⁷

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$.⁸ 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$.⁹ 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. 2019). 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 The ITG study

This section contains the results and analysis of the ITG study of the equilibria that we obtained in § 2. In § 5.1, we will present the specific details, including the values of different parameters used for the simulation and the reasoning behind our choices. In § 5.2, we will introduce the local magnetic shear, a quantity that characterizes the stability of an equilibrium to the ITG mode. In the final section, we will plot and compare the results of all the different local equilibria and explain the stability of the high-$\beta$ equilibria.

5.1 Details of the study

The most important form of electrostatic instability that arises at low wavenumbers, the ITG (Cowley, Kulsrud & Sudan 1991), 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}}$.¹⁰ 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 1998; Hirose 2000) 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 1981), 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. 2003). 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 1981; Dewar & Glasser 1983)

(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. (1996) 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 (2003a) – 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 The TEM study

In this section, we will present our analysis of the TEM instability. In § 6.1, we explain the TEM and tabulate the parameters for which we perform our study. After that, we define the maximum TEM growth rate. In § 6.2, we present the electron precession drift as a characteristic of TEM stability of a local equilibrium. In the final section, we plot the results in a fashion similar to the previous section.

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 1976). 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 1983; Roach, Connor & Janjua 1995). 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. 1995). 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 2015) 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.

7 Linear electromagnetic study

To see if the stability trend seen in the electrostatic study holds when we include electromagnetic effects, we perform an electromagnetic microstability analysis for all the nominal local equilibria. This analysis is similar to the work done by Belli & Candy (2010), albeit we are testing $\beta \sim 1$ equilibria using an initial-value solver. We solve the linear, collisionless, gyrokinetic model allowing for non-zero magnetic field perturbations $\delta A_{\parallel }$ and $\delta B_{\parallel }$ using the GS2 code. We use the nominal gradient scale lengths for this study (given in table 4). First, we plot the growth rate spectrum with $k_y \rho _{{i}}$ for the positive-triangularity cases in figure 18.

Figure 18. This figure shows comparison between the electrostatic (abbreviated ES) and electromagnetic (abbreviated EM) growth rates for all the nominal positive-triangularity equilibria. Some of the branches have been labelled by their corresponding mode names. Notice the KBM in the intermediate and high-$\beta$ equilibria in figures 21(e) and 21( f), respectively and the emergence of the collisionless-microtearing and electromagnetic-ETG modes in figure 21(c) for all values of $k_y\rho _{{i}}$ and figures 21(e) and 21( f) for $k_y\rho _{{i}} > 1$.

Since GS2 only calculates the maximum growth rate, in the inner core, we observe the finite-$\beta$ stabilization of ITG (Candy & Waltz 2003b) until the emergence of collisionless-microtearing (MTM) (Kotschenreuther et al. 1995a; Guttenfelder et al. 2011; Dickinson et al. 2013) and electromagnetic-ETG (EM-ETG) modes (Kim & Horton 1991; Joiner & Hirose 2007) in figure 18(c)¹¹ . These modes arise on small scales radially ($\textit {O}(\rho _{{e}})$) and are extended in the ballooning angle $\theta$. To accurately capture the extended eigenfunctions corresponding to these modes, we have to choose a wide range in the field-line-following coordinate ($\theta \in [-119{\rm \pi}, 119{\rm \pi} ]$) with $21$ points over a $2{\rm \pi}$ interval. The eigenfunctions for two values of $k_y\rho _{{i}}$ in figure 19(c) are shown in figure 24 in  Appendix E.

Figure 19. This figure shows comparison between the electrostatic (ES) electromagnetic (EM) growth rates for all the nominal negative triangularity equilibria. The sudden jump in figure 19( f) around $k_y\rho _{{i}} = 4.5$ is a different branch of the collisionless MTM. Note also that the growth rate around $k_y \rho _{{i}} = 0$ in figure 19(e) goes to fixed value of $\gamma a_{{N}}/v_{\mathrm {th},i} = 0.152$ since the equilibrium is unstable to the ideal-ballooning mode.

In the outer core, we observe finite-$\beta$ stabilization only for the low-$\beta$ equilibrium. For the intermediate and high-$\beta$ cases, the electrostatic modes are replaced by kinetic ballooning modes (KBMs) at low wavenumbers and collisionless-MTM and EM-ETG at high wavenumbers. Overall, positive-triangularity high-$\beta$ equilibria are more unstable than the low or intermediate ones in the inner core due to the collisionless-MTM and the EM-ETG mode. As we move towards the outer core, high-$\beta$ equilibria become much more stable – exactly the opposite trend compared with the inner core. This means that the outer core is more stable for the high-$\beta$ equilibria.

Next, we plot the growth rates for the negative-triangularity equilibria in figure 19. The negative-triangularity, inner-core equilibria are also stabilized due to finite-$\beta$ effects. This effect is also visible for the intermediate-$\beta$ case for $k_y\rho _{{i}} > 0.5$ but absent for the high-$\beta$ cases as the ITG and TEM are superseded by electromagnetic modes: collisionless-MTM and the EM-ETG mode. For the outer-core cases, we see the exact same pattern as the inner core – finite-$\beta$ stabilization for the low-$\beta$, KBM and collisionless-MTMs for the intermediate-$\beta$ and collisionless-MTMs and EM-ETG modes for the high-$\beta$. Note that unlike the positive-triangularity case, we do not observe the KBM in the negative-triangularity, high-$\beta$ equilibria. Moreover, the growth rates are also smaller than the positive triangularity cases for a wide range of wavenumbers ($k_y \rho _{{i}} \in [0.01, 4.5]$). Since turbulence is most likely to peak at lower wavenumbers, the growth rate characteristics of negative-triangularity, high-$\beta$ are the most favourable.¹²

In summary, we find that turning on the electromagnetic effects destabilizes the high-$\beta$ equilibria. Figures 18(c), 19(c) and 19( f) show the emergence of collisionless-MTMs and EM-ETG modes whereas figure 18(c) also shows instability to the KBM in the range $k_y\rho _{{i}} \in [0.01, 0.5]$. These equilibria are much more unstable than the low-$\beta$ ones in the inner core but they are much more stable as we move towards the outer core, with negative-triangularity high-$\beta$ equilibria showing the best characteristics. We believe that the outer-core stability is due to a large Shafranov shift and strong shaping. We also argue that negative triangularity has better characteristics than positive triangularity due to stronger shaping, i.e. ‘squareness’ that we discussed in figure 3. It is also interesting to note that for the negative-triangularity equilibria, the growth rate trend matches that of the ideal-ballooning stability in figure 6 – less stable towards the core, more stable towards the edge. Since outer-core or edge transport is usually a limiting factor in experimental low-$\beta$ equilibria, these equilibria may be a novel alternative to realize higher-power devices.

In this section, the lowest wavenumber that we have scanned is $k_y \rho _{{i}} = 0.01$. Using the same analyses as presented at the end of § 6, the longest wavelength $\lambda \approx 3.1 \textrm {m}$. Assuming $n = 1$ corresponds to $\lambda = a_{{N}} = 0.68\ \textrm {m}$, the longest mode would require $n < 1$. Modes with such long wavelengths violate the local assumption of flux tube codes. However, the growth rates go to zero at low wavenumbers (with the exception of figure 19e) and peak at $k_{y, \mathrm {peak}} \rho _{{i}} \gg 0.01$. Calculating the value of $n$ at the wavenumber $k_{\mathrm {peak}}\rho _{{i}}$, we obtain the lowest $n$ in figures 19(e) and 18( f) as $n = 0$ and $n=5$, respectively. Therefore, the reader must take the low-wavenumber results shown in figures 19(e) and 18( f) with a grain of salt.

Note that this is a collisionless electromagnetic analysis. It is possible that in the presence of collisions, more modes like collisional-MTMs (Applegate et al. 2007; Patel et al. 2021), collisional EM-ETG etc. arise and supersede their collisionless variants (Dickinson et al. 2013) as the most dominant mode. On the other hand, it is also possible to achieve lower growth rates if we include the stabilizing effect of velocity shear (Patel et al. 2021).

8 Summary and conclusions

We began this work by deriving the Grad–Shafranov equation and briefly discussed existing approaches used to obtain $\beta \sim 1$ axisymmetric equilibria in § 2. In § 2.1, we generated numerical high-$\beta$ equilibria using the 3-D equilibrium code VMEC. We produced different equilibria based for three beta values: low, intermediate and high ($\beta \sim 1$), each differing from the rest by at least an order or magnitude. We also chose two different boundary shapes with negative and positive triangularity. For each global equilibrium, we picked two local equilibria: one in the inner core and the other in the outer core. This gave us a total of $12$ local equilibria.

Upon ensuring that these equilibria are smooth and well resolved, we first tested their stability against the infinite-$n$ ideal-ballooning mode in § 3. Since the ideal-ballooning mode is an MHD instability, we chose the local equilibria such that all the high and low-$\beta$ equilibria are ideal-ballooning stable. We also elucidated the idea of varying the magnetic shear and pressure gradient self-consistently, the Greene–Chance analysis in § 3.2. This powerful technique enabled us to know how far the local equilibria are from the region of marginal stability.

In § 4, we briefly explained the linear, collisionless, gyrokinetic model. In §§ 5 and 6, we use this to study the stability of all the local equilibria to the two most virulent electrostatic modes of turbulence: ITG and TEM. We found a clear inverse relationship between the beta value and the growth rates – increasing the beta value stabilized both the ITG mode and the TEM. Using a Greene–Chance analysis, we also scanned the maximum ITG, TEM growth rates vs the temperature and density gradient scale lengths, respectively. This was important to ensure that these equilibria are not ‘stiff’, i.e. the growth rates do not increase sharply as the gradients exceed some threshold. In § 5.2, we explained how a large negative local shear resulting from a large Shafranov shift stabilizes the ITG mode and in § 6.2 we showed how the reversal of the precession of the trapped electrons stabilizes the TEM throughout the whole range of gradient scale lengths.

The effect of electromagnetic fluctuations can be important for all equilibria, especially the intermediate and high-$\beta$ ones. To that end, in § 7 we performed an electromagnetic study of all the nominal equilibria. We found that the stability trend seen for the electrostatic case did not hold after turning on electromagnetic effects. However, even though the high-$\beta$ equilibria were more unstable than the low-$\beta$ ones due to collisionless-MTMs and EM-ETGs in the inner core, they were much more stable than low-$\beta$ in the outer core. We found that negative-triangularity, high-$\beta$ equilibria were stable to the KBM. We believe that this is due to the strong shaping (‘squareness’) of the negative-triangularity high-$\beta$ equilibria. This indicates that turbulent transport may flatten the pressure gradient in the core but may not significantly affect the pressure gradient toward the edge for high-$\beta$ equilibria.

This work suggests many promising avenues for future research. One could repeat the same study with a wider range of input parameters to determine whether the microstability trends we have found are not strongly dependent on the input parameters. An important path to explore would be to map the trajectory of these high-$\beta$ equilibria, starting from a low-beta equilibrium in the $\hat {s}-\alpha _{\mathrm {MHD}}$ space such that it is always stable to major sources of disruption. This would guarantee a high-$\beta$ operation free from any ideal-MHD instability. Once such a path is found, the final step would be to do a microstability calculation scan for different values of collisionalities followed by a transport calculation to obtain the evolution of the density and temperature profiles.