2.1. The Æ in any omnigenous system

We first note that the integral over $z$ can be rewritten into a convenient form. We define two functions that are independent of $z$, namely

(2.1a,b)\begin{equation} c_0 = \frac{ \Delta \psi_t}{\hat{\omega}_\alpha(\lambda)} \frac{\mathrm{d} \ln(n)}{\mathrm{d} \psi_t} \left( 1 - \frac{3}{2} \eta \right),\quad c_1 = 1 - \frac{\Delta \psi_t}{\hat{\omega}_\alpha(\lambda)} \frac{\mathrm{d} \ln(n)}{\mathrm{d} \psi_t} \eta.\end{equation}

With these functions, the integral over the normalised energy $z$ reduces to the following form:

(2.2)\begin{equation} I_z(c_0,c_1) = \frac{8}{3 \sqrt{\rm \pi}} \int_0^\infty \exp(-z) z^{3/2} \mathcal{R} \left[ c_0 - c_1 z \right] \mathrm{d} z. \end{equation}

This integral can be solved analytically, and its functional form depends on the signs of $c_0$ and $c_1$, resulting in four different conditions. The easiest case to evaluate is the case where $c_0<0$ and $c_1>0$. In this case, the argument of the ramp function is always negative, and hence the integral reduces to zero. The second case is when the argument of the ramp function is always positive, which occurs whenever $c_0\geq 0$ and $c_1 \leq 0$. The integral then reduces to the following form:

(2.3)\begin{equation} I_z = 2c_0 - 5c_1 . \end{equation}

There are two cases left to consider. First, we inspect the case where the argument of the ramp function is positive for low $z$ but becomes negative for high $z$, that is, $c_0 \geq 0$ and $c_1 > 0$. The unique point where the argument of the ramp function vanishes is the following:

(2.4)\begin{equation} z_* = \frac{c_0}{c_1}. \end{equation}

Thus, the integral becomes

(2.5)\begin{equation} I_z = \frac{8}{3 \sqrt{\rm \pi}} \int_0^{z_*} \exp(-z) z^{3/2} \left( c_0 - c_1 z \right) \mathrm{d} z. \end{equation}

This integral can be expressed in terms of the error function, $\text {erf}(x) = 2/\sqrt {{\rm \pi} } \int _0^x \exp (-t^2) \,\mathrm {d}t$,

(2.6)\begin{equation} I_z = (2c_0 - 5c_1) \text{erf}\left({\sqrt{\frac{c_0}{c_1}}}\right) + \frac{2}{3\sqrt{\rm \pi}} (4c_0 + 15c_1) \sqrt{\frac{c_0}{c_1}} \exp\left( - \frac{c_0}{c_1} \right). \end{equation}

The final case is that where the argument of the ramp function is negative for low $z$ but becomes positive for high $z$, that is, $c_0 < 0$ and $c_1 \leq 0$. The integral then becomes

(2.7)\begin{equation} I_z =(2c_0 - 5c_1) \left[ 1 - \text{erf}\left({\sqrt{\frac{c_0}{c_1}}}\right) \right] -\frac{2}{3\sqrt{\rm \pi}} (4c_0 + 15c_1) \sqrt{\frac{c_0}{c_1}} \exp\left( -\frac{c_0}{c_1} \right). \end{equation}

Note that $I_z \geq 0,\ \forall (c_0,c_1) \in \mathbb {R}^2$, which can also be seen in figure 1. The Æ can now be found by executing the integral over the remaining coordinate

(2.8)\begin{equation} A = \frac{3}{16} \frac{\Delta \psi_t \Delta \alpha_C L}{B_0} n_0 T_0 \int_{\{ \lambda \}} \mathrm{d} \lambda \sum_{\text{wells}(\lambda)} I_z(c_0,c_1) \hat{\omega}_\alpha^2 \hat{g}^{1/2}. \end{equation}

Note that this expression is completely general; no approximations have been made in executing these integrals, aside from the preceding assumption of omnigeneity.

Figure 1. Contour plot of $I_z$ as a function of $c_0$ and $c_1$.

It is also interesting to note that from this expression one can see that there are no tokamak configurations with vanishing Æ, at least in leading order near the axis. This conclusion can most readily be drawn by investigating the expression for $\omega _\alpha$ from Connor, Hastie & Martin (Reference Connor, Hastie and Martin1983). Here, one can find that there is always a zero crossing for $\omega _\alpha$ (with no pressure gradient), which implies that $c_0$ and $c_1$ must change sign. As such, the Æ must be non-zero (as either $I(c_0,c_1)$ or $I(-c_0,-c_1)$ must be non-zero). Formally, this corresponds to the fact that such a zero crossing implies that the device does not have the so-called maximum-$\mathcal {J}$ property, which is required for the linear stability of TEMs (Proll et al. Reference Proll, Helander, Connor and Plunk2012).Footnote ¹

To make further progress in solving (2.8), one requires the function $\hat {\omega }_\alpha (\lambda )$, which in turn requires a specification of the equilibrium. In this paper, we will use the local construction of the equilibrium, employing a formalism developed by Mercier & Luc (Reference Mercier and Luc1974).

2.2. Construction of local equilibria

Equilibria are constructed by finding a radially local solution to the Grad–Shafranov equation, and this solution allows us to find $\hat {\omega }_\alpha$. We highlight the essential components of this derivation, which essentially follows the steps taken by Miller et al. (Reference Miller, Chu, Greene, Lin-Liu and Waltz1998), and a thorough overview is given by Candy (Reference Candy2009).

The Mercier–Luc formalism requires the shape of the flux surface, the poloidal field $B_p$ on that flux surface, the gradients of the pressure $p(\psi )$ and the toroidal field function $f(\psi ) = R B_\phi$ on the flux surface, where $R$ is the major radial coordinate, $B_\phi$ is the toroidal component of the magnetic field and $\psi$ is the poloidal flux. We parameterise the flux surface as $R_s = R_s(l)$ and $Z_s = Z_s(l)$, where $l$ measures the poloidal arclength along the flux surface. It is also useful to define a tangential angle $u$, which measures the angle between the unit vector in the major radial direction $\boldsymbol {e}_R$ and the vector tangential to the flux surface $\boldsymbol {e}_l$ clockwise, thus

(2.9a)\begin{gather} \frac{\mathrm{d}R_s(l)}{\mathrm{d} l} = \cos u, \end{gather}

(2.9b)\begin{gather}\frac{\mathrm{d}Z_s(l)}{\mathrm{d} l} ={-} \sin u. \end{gather}

With this definition, the angle $u$ can be calculated by $\mathrm {d} u / \mathrm {d} l = - 1 / R_c$, where $R_c(l)$ is the radius of curvature of the poloidal cross-section, and the negative sign arises because the poloidal arclength is measured clockwise. We go on to introduce a radial-like expansion variable $\rho$ which is zero on the given flux surface, in terms of which the cylindrical coordinates become

(2.10a)\begin{gather} R(\rho,l) = R_s(l) + \rho \sin u, \end{gather}

(2.10b)\begin{gather}Z(\rho,l) = Z_s(l) + \rho \cos u. \end{gather}

The metric tensor in these coordinates has non-zero components only on the diagonal (which is to be expected as we ensured orthogonality in the construction),

(2.11)\begin{equation} g_{ij} = \text{diag}\left[ \left(1 - \frac{\rho}{R_c} \right)^2, 1, R^2 \right], \end{equation}

where we use the convention $x^1 = l$, $x^2 = \rho$, $x^3 = \phi$. The local solution is now constructed by expanding in $\rho$,

(2.12a)\begin{gather} \psi \approx \psi_0 + \rho \psi_1 + \frac{\rho^2 }{2}\psi_2, \end{gather}

(2.12b)\begin{gather}p'(\psi) \approx p'(\psi_0) \end{gather}

(2.12c)\begin{gather}f'(\psi) \approx f'(\psi_0) \end{gather}

and substitute into the Grad–Shafranov equation, which to leading order reduces to

(2.13)\begin{equation} \psi_2 = \left( \sin(u) + \frac{R_s}{R_c} \right) \frac{\psi_1}{R_s} - \mu_0 R_s^2 p'(\psi_0) - f(\psi_0)f'(\psi_0). \end{equation}

This allows one to find the radial variation of the poloidal magnetic field by using (Helander & Sigmar Reference Helander and Sigmar2005)

(2.14)\begin{equation} B_p = \frac{ | \boldsymbol{\nabla} \psi |}{R}, \end{equation}

resulting in

(2.15)\begin{equation} B_p(l,\rho) = \frac{\psi_1}{R_s} \left( 1 + \rho \left[ \frac{1}{R_c} - \frac{ \mu_0 R_s^2 p'(\psi_0)}{\psi_1} - \frac{f(\psi_0) f'(\psi_0)}{\psi_1} \right] \right). \end{equation}

From this equation we can immediately see that $\psi _1/R_s = B_{p,s}$, with $B_{p,s}$ being the poloidal field on the flux-surface as indicated by the subscript. As such, the poloidal field strength can be written as

(2.16)\begin{equation} \left| B_p(l,\rho)\right| = B_{p,s}\left( 1 + \rho \left[ \frac{1}{R_c} - \frac{\mu_0 R_s p'}{B_{p,s}} - \frac{ff'}{R_s B_{p,s}} \right] \right) \equiv B_{p,s} \left( 1 + \rho \partial_\rho b_p \right). \end{equation}

The toroidal field is found from its definition $B_\phi = f(\psi )/R$, resulting in

(2.17)\begin{equation} \left|B_\phi(l,\rho) \right| = B_{\phi,s} \left( 1 + \rho \left[\frac{ f'(\psi_0)}{f(\psi_0)} R_s B_{p,s} - \frac{\sin u}{R_s}\right]\right) \equiv B_{\phi,s} \left( 1 + \rho \partial_\rho b_\phi \right), \end{equation}

where $B_{\phi,s} = f(\psi _0)/R_s$. The total magnetic field strength is also readily derived:

(2.18)\begin{equation} B = \sqrt{B_{\phi,s}^2+B_{p,s}^2} \left(1 + \rho \left[ \frac{B_{\phi,s}^2 \partial_\rho b_\phi + B_{p,s}^2 \partial_\rho b_p }{B_{\phi,s}^2+B_{p,s}^2} \right]\right) \equiv B_s \left( 1 + \rho \partial_\rho b \right). \end{equation}

Note that the derivatives $\partial _\rho b_p$, $\partial _\rho b_\phi$ and $\partial _\rho b$ are given in square brackets. The radial variation of the poloidal line element is readily found from the metric tensor,

(2.19)\begin{equation} \mathrm{d}l = \left(1 - \frac{\rho}{R_c}\right) \left[ \mathrm{d}l \right]_{\rho=0}. \end{equation}

In these equations $f'(\psi _0)$ is treated as a free parameter, but it is difficult to ascertain if the chosen value of this parameter is realistic. It is more convenient, however, to specify the magnetic shear, which is related to $f'(\psi _0)$. This can be made explicit by investigating the safety factor

(2.20)\begin{equation} q = \frac{f(\psi)}{2 {\rm \pi}} \int \frac{ \mathrm{d}l }{R_s^2 B_{p,s}}. \end{equation}

Taking the derivative of the safety factor with respect to $\psi$, one finds an equation describing this relationship:

(2.21)\begin{equation} \partial_\psi q = \frac{f'}{f} q + f \frac{1}{2{\rm \pi}} \int \frac{\mathrm{d}l}{R_s^3 B_{p,s}^2} \left(-\frac{2}{R_c} - \frac{2 \sin u}{R_s} + \frac{\mu_0 R_s p'}{B_{p,s}} + \frac{f f'}{R_s B_{p,s}} \right). \end{equation}

We also wish to relate the arclength along a magnetic field line to the poloidal arclength. These quantities are related as

(2.22)\begin{equation} \mathrm{d}\ell =\left| \frac{B}{B_p} \right| \mathrm{d}l. \end{equation}

Finally, the poloidal coordinate can be expressed in terms of the poloidal angle $\theta$ rather than the poloidal arclength by

(2.23)\begin{equation} l_\theta \equiv \frac{\mathrm{d}l}{\mathrm{d}\theta} = \sqrt{ (\partial_\theta R_s)^2 + (\partial_\theta Z_s)^2 }, \end{equation}

and the total arclength thus becomes

(2.24)\begin{equation} L = \oint l_\theta \left| \frac{B}{B_p} \right| \mathrm{d}\theta. \end{equation}

2.3. Non-dimensionalisation and Æ

We proceed to make the various functions dimensionless as in Roach, Connor & Janjua (Reference Roach, Connor and Janjua1995), and in doing so we will introduce various dimensionless constants which will be useful for the remainder of the analysis. We assume that we have been given the dependencies of the various functions in terms of the minor radial coordinate $r$, which in turn relates to the major radial coordinate $R_0$ through the inverse aspect ratio of the flux surface in question $\epsilon =r/R_0$. Furthermore, we define our reference field $B_0$ through the relation $f(\psi _0) = B_0 R_0$. Let us now define various dimensionless functions of interest:

(2.25a)\begin{gather} \hat{R}_s = R/R_0, \end{gather}

(2.25b)\begin{gather}\hat{Z}_s = Z/R_0, \end{gather}

(2.25c)\begin{gather}\hat{R}_c = R_c/r, \end{gather}

(2.25d)\begin{gather}\hat{l}_\theta = l_\theta/r, \end{gather}

(2.25e)\begin{gather}\hat{B}_\phi = B_\phi/B_0, \end{gather}

(2.25f)\begin{gather}\hat{B} = B/B_0. \end{gather}

One also needs to relate $\psi$ to $r$, which can be done by investigating the poloidal field as in (2.14)

(2.26)\begin{equation} B_p = \frac{\partial_r \psi}{R_0} \frac{ \left| \boldsymbol{\nabla} r \right|}{\hat{R}_s}. \end{equation}

We go on to identify two factors in the above expression, namely

(2.27)\begin{equation} B_{p,0} \equiv \partial_r \psi /R_0 \end{equation}

and

(2.28)\begin{equation} \hat{B}_{p,s} \equiv \left| \boldsymbol{\nabla} r \right| / \hat{R}_s. \end{equation}

Inserting these into the equation for the safety factor (2.20), one finds

(2.29a,b)\begin{equation} B_{p,0} = \frac{\gamma \epsilon}{q} B_0,\quad \gamma \equiv \frac{1}{2 {\rm \pi}} \oint \frac{\hat{l}_\theta}{\hat{R}_s^2 \hat{B}_{p,s}} \mathrm{d} \theta. \end{equation}

We proceed to define a dimensionless pressure gradient, analogous to the $\alpha$ parameter used in $s$–$\alpha$ geometry:

(2.30)\begin{equation} \alpha =- \frac{2 \mu_0 \epsilon^2 R_0^2 p'}{B_{p,0}} =-\epsilon r \frac{\mathrm{d} p}{\mathrm{d} r}/ \frac{B_{p,0}^2}{2 \mu_0}. \end{equation}

Note that this dimensionless pressure gradient is not related to the Clebsch angle. The pressure gradient can in turn be used to define a dimensionless toroidal current density:

(2.31)\begin{equation} \sigma = \left( \mu_0 p' + \frac{ff'}{R_0^2} \right) \frac{\epsilon R_0^2}{B_{p,0}} = \frac{q}{\gamma} f' R_0 - \frac{\alpha}{2 \epsilon}. \end{equation}

We go on to define the shear $s$ in the following manner:

(2.32)\begin{equation} s = \epsilon R_0^2 B_{p,0} \partial_\psi \ln q = \frac{r}{q} \frac{\partial q}{\partial r} \end{equation}

which can be substituted into (2.21) to relate the shear to $f'R_0$ as

(2.33)\begin{equation} s = \frac{\gamma \epsilon^2}{q} f'R_0 - \frac{2}{\gamma} C_1 - \frac{2 \epsilon}{\gamma} C_2 - \frac{\alpha}{2\gamma \epsilon} C_3 + \frac{q}{\gamma^2} C_4 f'R_0, \end{equation}

where we have defined the geometric constants $C_1$ to $C_4$ as

(2.34a)\begin{gather} C_1 = \frac{1}{2 {\rm \pi}} \oint \frac{\hat{l}_\theta}{\hat{R}_c \hat{R}_s^3 \hat{B}_{p,s}^2} \,\mathrm{d}\theta, \end{gather}

(2.34b)\begin{gather}C_2 = \frac{1}{2 {\rm \pi}} \oint \frac{\hat{l}_\theta \sin u}{\hat{R}_s^4 \hat{B}_{p,s}^2} \,\mathrm{d}\theta, \end{gather}

(2.34c)\begin{gather}C_3 = \frac{1}{2 {\rm \pi}} \oint \frac{\hat{l}_\theta}{\hat{R}_s^2 \hat{B}_{p,s}^3} \,\mathrm{d}\theta, \end{gather}

(2.34d)\begin{gather}C_4 = \frac{1}{2 {\rm \pi}} \oint \frac{\hat{l}_\theta}{ \hat{R}_s^4 \hat{B}_{p,s}^3} \,\mathrm{d}\theta. \end{gather}

The radial derivatives of the magnetic field become

(2.35a)\begin{gather} r \partial_\rho b_p = \left( \frac{1}{\hat{R}_c} - \frac{\alpha}{2 \epsilon \hat{B}_{p,s}} \left[\frac{1}{\hat{R}_s} - \hat{R}_s \right] - \frac{\sigma}{\hat{R}_s \hat{B}_{p,s}} \right), \end{gather}

(2.35b)\begin{gather}r \partial_\rho b_\phi = \epsilon \left( \frac{\gamma^2\epsilon}{q^2} \left[ \sigma + \frac{\alpha}{2 \epsilon} \right] \hat{R}_s \hat{B}_{p,s} - \frac{\sin u}{\hat{R}_s} \right). \end{gather}

These expressions are the same as Roach et al. (Reference Roach, Connor and Janjua1995, (14)), where the differences in sign arise because the sign convention for $\hat {R}_c$ is different here and $\rho$ has the opposite sign. Finally, we express the total magnetic field length as

(2.36a,b)\begin{equation} L = \frac{q\xi}{\gamma} R_0,\quad \xi \equiv \oint \frac{\hat{l}_\theta \hat{B}_s}{\hat{B}_{p,s}} \mathrm{d} \theta. \end{equation}

We now turn our attention to the precession frequency, which we calculate from (1.3a). To simplify the calculation slightly, we note that the operator $\Delta \psi _t \partial _{\psi _t} \approx \Delta \psi \partial _\psi$ to leading order around smallness of the radial coordinate $\rho$, as we can approximate $\Delta \psi _t \approx \Delta \psi \partial _\psi \psi _t$. Using this identity, we find the same expression as in Roach et al. (Reference Roach, Connor and Janjua1995),

(2.37)\begin{equation} \hat{\omega}_\alpha(\lambda) =- \frac{\Delta \psi}{R_0^2 B_{p,0}} \left\langle \frac{1}{\epsilon} \left( 2 \left[1 - \lambda \hat{B} \right] \left[r\partial_\rho b - r\partial_\rho b_p - \frac{1}{\hat{R}_c}\right] - \lambda \hat{B} r \partial_\rho b \right)\bigg/ \hat{B}_{p,s} \hat{R}_s \right\rangle_\lambda, \end{equation}

where we define the bounce averaging operator in angular brackets as

(2.38)\begin{equation} \langle \cdots \rangle_\lambda = \frac{\displaystyle\int \mathrm{d} \theta \cdots \hat{l}_\theta \frac{\hat{B}_s}{\hat{B}_{p,s}}/\sqrt{1 - \lambda \hat{B}}}{\displaystyle\int \mathrm{d} \theta ~ \hat{l}_\theta \frac{\hat{B}_s}{\hat{B}_{p,s}}/ \sqrt{1 - \lambda \hat{B}}}. \end{equation}

We rewrite the precession frequency as

(2.39)\begin{equation} \hat{\omega}_\alpha \equiv- \frac{\Delta \psi}{R_0^2 B_{p,0}} \hat{\omega}_\lambda. \end{equation}

Next, we investigate the Jacobian $\hat {g}^{1/2}$, which is the normalised bounce time, and find that it is equal to

(2.40)\begin{equation} \hat{g}^{1/2} = \frac{\displaystyle\int \mathrm{d} \theta \hat{l}_\theta \frac{\hat{B}_s}{\hat{B}_{p,s}}/ \sqrt{1 - \lambda \hat{B}}}{\xi}. \end{equation}

We rescale it with a factor $\epsilon ^{1/2}$ to account for the fact that in smallness of $\epsilon$ the integrand of the bounce time goes as $1/\sqrt {\epsilon }$. Therefore, we define

(2.41)\begin{equation} \hat{g}_\epsilon^{1/2} \equiv \hat{g}^{1/2} \sqrt{\epsilon}. \end{equation}

The Æ now becomes

(2.42)\begin{equation} A = \frac{3}{16} \frac{\Delta \psi_t \Delta \alpha_C L}{B_0} n_0 T_0 \sqrt{\epsilon} \left(\frac{\Delta \psi}{R_0^2 B_{p,0}} \right)^2 \left( \frac{1}{\epsilon} \int_{\{ \lambda \}} \mathrm{d}\lambda \sum_{\text{wells}(\lambda)} I_z(c_0,c_1) \hat{\omega}_\lambda^2 \hat{g}_\epsilon^{1/2} \right), \end{equation}

where the prefactor $1/\epsilon$ to the integral deliberately not cancelled against the $\sqrt {\epsilon }$, so that the integral in brackets is to lowest order independent of $\epsilon$, as the integration range scales as $\epsilon$. With the above expression, we go on to define a dimensionless Æ. We take steps in accordance with Mackenbach et al. (Reference Mackenbach, Proll and Helander2022), and calculate the fraction of the total thermal energy that is available. The thermal energy of a plasma in a flux tube can be calculated by expanding around $\psi _t = \psi _{t,0}$ and $\alpha _C = \alpha _{C,0}$ and retaining only the constant terms, resulting in

(2.43)\begin{equation} E_t = \int \frac{3}{2} \frac{n T}{B} \mathrm{d} \psi_t \,\mathrm{d} \alpha_C \,\mathrm{d} \ell = \frac{3}{2} n_0 T_0 \frac{\Delta \psi_t \Delta \alpha_C L }{B_0} \frac{1}{\xi} \oint \hat{l}_\theta \hat{B}_{p,s}^{-1} \,\mathrm{d} \theta. \end{equation}

We then define the Æ as a fraction of the thermal energy as

(2.44)\begin{equation} \hat{A} = \frac{A}{E_t}. \end{equation}

Simplifying the expression using $\Delta \psi = \Delta r \partial _r \psi$, one finds that

(2.45)\begin{equation} \hat{A} = \frac{1}{8} \left( \frac{\Delta r}{R_0} \right)^2 \frac{\xi \sqrt{\epsilon}}{\displaystyle\oint \hat{l}_\theta \hat{B}_{p,s}^{-1}\,\mathrm{d} \theta } \cdot \frac{1}{\epsilon} \int_{\{ \lambda \}} \mathrm{d}\lambda \sum_{\text{wells}(\lambda)} I_z(c_0,c_1) \hat{\omega}_\lambda^2 \hat{g}_\epsilon^{1/2}. \end{equation}

Here $\Delta r$ measures the length scale over which energy is available, i.e. a typical length scale over which gradients can be flattened. We take this to be proportional to the correlation length, typically found to be the gyroradius. Therefore, let us set

(2.46)\begin{equation} \Delta r = C_r \rho_{g}, \end{equation}

where $\rho _{g}$ is the gyroradius, and $C_r$ some function of order $O(\rho _{g}^0)$. This function is not known a priori, and may vary. For example, if there are large radial streamers present in the system $C_r$ may be significantly increased. The dimensionless Æ now becomes

(2.47)\begin{equation} \hat{A} = \frac{1}{8} \left( \frac{\rho_\text{g}}{R_0} \right)^2 \frac{ C_r^2 \xi \sqrt{\epsilon}}{\displaystyle\oint \hat{l}_\theta \hat{B}_{p,s}^{-1} \,\mathrm{d} \theta} \cdot \frac{1}{\epsilon} \int_{\{ \lambda \}} \mathrm{d}\lambda \sum_{\text{wells}(\lambda)} I_z(c_0,c_1) \hat{\omega}_\lambda^2 \hat{g}_\epsilon^{1/2}. \end{equation}

This expression has various scalings which are of interest. First, we see that reducing the aspect ratio for fixed $\rho _{g}/R_0$ is beneficial since it leads to fewer trapped particles. Note that, in the limit of a large-aspect-ratio, the trapping fraction scales as $\sqrt {\epsilon }$, which is the same dependency found here. A reduction in the expansion parameter $\rho _{g}/R_0$ (at fixed $\epsilon$) is also found to help decrease Æ.

As a final step, we introduce the dimensionless density gradient

(2.48)\begin{equation} R_0^2 B_{p,0} \partial_\psi \ln n = \frac{R_0}{n} \frac{\partial n}{\partial r} \equiv- \hat{\omega}_n, \end{equation}

with which $c_0$ and $c_1$ reduce to an especially simple form

(2.49a,b)\begin{equation} c_0 = \frac{\hat{\omega}_n}{\hat{\omega}_\lambda} \left( 1 - \frac{3}{2}\eta \right),\quad c_1 = 1 - \frac{\hat{\omega}_n}{\hat{\omega}_\lambda} \eta. \end{equation}

Importantly let us make note of the fact that the radial coordinate $r$ may have different conventions. In previous investigations (Mackenbach et al. Reference Mackenbach, Proll and Helander2022, Reference Mackenbach, Proll, Wakelkamp and Helander2023a) the radial coordinated was defined via the square root of the toroidal flux

(2.50)\begin{equation} r_{\rm eff} \propto \sqrt{\psi_t}, \end{equation}

with $\psi _t$ being the toroidal flux passing through the flux surface in question. A different choice of the radial coordinate $r$ will influence various quantities on which Æ depends, such as (2.46) and (2.48). More specifically, the length scale $\Delta r$ expressed in terms of $r_{\rm eff}$ is

(2.51)\begin{equation} \Delta r = \Delta r_{\rm eff} \frac{\partial r}{\partial r_{\rm eff}}. \end{equation}

In the aforementioned investigations $\Delta r_{\rm eff}$ was chosen as $\Delta r_{\rm eff} = \rho$, resulting in a good correlation with turbulent energy fluxes. Therefore, we choose $C_r$ such that $\Delta r_{\rm eff} = \rho$, which means that

(2.52)\begin{equation} C_r = \frac{\partial r}{\partial r_{\rm eff}}, \end{equation}

and we shall use this choice of $C_r$ from here on.

2.4. Miller geometry

Finally, we choose our equilibrium to be of the type discussed in Miller et al. (Reference Miller, Chu, Greene, Lin-Liu and Waltz1998). The key step is to parameterise the flux surface as a standard D-shaped tokamak in terms of the poloidal angle $\theta$:

(2.53a)\begin{gather} R_s(\theta) = R_0 + R_0 \epsilon \cos(\theta + \arcsin [\delta] \sin \theta), \end{gather}

(2.53b)\begin{gather}Z_s(\theta) = R_0 \kappa \epsilon \sin \theta. \end{gather}

Here, $R_0(r)$ is the centre of the flux surface, $\kappa (r)$ is the elongation and $\delta (r)$ is the triangularity. An important feature of this parameterisation is that it is up–down symmetric, which can be seen by invariance under $(Z_s,\theta )\mapsto -(Z_s,\theta )$. The poloidal magnetic field can then be calculated by (2.26), and the equilibrium is fully specified by the following set of nine parameters; $[\epsilon,\kappa,\delta, s_\kappa, s_\delta, \partial _r R_0, q, s, \alpha ]$, where $s_\kappa =r \partial _r \ln \kappa$ and $s_\delta = r \partial _r \arcsin (\delta )$. Henceforth we shall refer to this set of numbers which determines the local geometry as a ‘Miller vector’:

(2.54)\begin{equation} \boldsymbol{M} = [\epsilon,\kappa,\delta, s_\kappa, s_\delta, \partial_r R_0, q, s, \alpha]. \end{equation}

Cross-sections are plotted in figure 2, to serve as a reference for the various shapes mentioned in the following sections. Furthermore, we recognise that it is possible to calculate the toroidal flux enclosed by a poloidal cross-section, and one may retrieve analytical expressions by expanding it around the smallness of $\epsilon$,

(2.55)\begin{align} \psi_t & = B_0 {\rm \pi}r^2 \frac{2}{\rm \pi} \int_{\rm \pi}^{0} \hat{Z}_s \frac{\partial_\theta R_s}{r\hat{R}_s} \mathrm{d} \theta \nonumber\\ & \approx B_0 {\rm \pi}r^2 \kappa \left( \frac{2 J_1(\arcsin \delta)}{\arcsin \delta} + \frac{J_2(2\arcsin \delta) }{2\arcsin \delta} \epsilon + O(\epsilon^2) \right), \end{align}

with $J_n(x)$ being the $n$th Bessel function of the first kind. In terms of an effective $r$, we then have

(2.56)\begin{align} r_{\rm eff} & = r \sqrt{ \frac{2}{\rm \pi} \int_{\rm \pi}^{0} \hat{Z}_s \frac{\partial_\theta R_s}{r\hat{R}_s} \mathrm{d} \theta } \nonumber\\ & \approx r \sqrt{\kappa} \sqrt{\frac{2 J_1(\arcsin \delta)}{\arcsin \delta}} \left( 1 + \frac{1}{8} \frac{J_2 (2 \arcsin \delta)}{J_1 (\arcsin \delta)} \epsilon + O(\epsilon^2) \right), \end{align}

and we find that the factor $C_r$ becomes

(2.57)\begin{equation} C_r= \left(\sqrt{\frac{2}{\rm \pi} \int_{\rm \pi}^{0} \hat{Z}_s \frac{\partial_\theta R_s}{r\hat{R}_s} \mathrm{d} \theta} \right)^{-1}. \end{equation}

For shaped equilibria (i.e. $\kappa \neq 1$, $\delta \neq 0$, or $\epsilon \rightarrow 1$), $C_r$ will differ from unity and one should keep this important caveat in mind.

Figure 2. Cross-sections in the $(R,Z)$-plane of tokamaks parameterised via (2.53). The parameters $(\kappa,\delta )$ vary from (a–d) as $(2/3,0.9)$, $(2/3,-0.9)$, $(3/2,0.9)$ and $(3/2,-0.9)$, respectively. All plots have $R_0=1$ and $\epsilon =1/3$.

2.5. An analytical limit: large-aspect-ratio $s$–$\alpha$ tokamak

We proceed to investigate a limiting case of Miller geometries; namely that of a large-aspect-ratio tokamak with circular flux surfaces and a steep local pressure gradient, which we shall henceforth refer to as the $s$–$\alpha$ limit, and this calculation is equivalent to analyses given in Connor et al. (Reference Connor, Hastie and Martin1983) and Roach et al. (Reference Roach, Connor and Janjua1995). This will serve as a computationally efficient model in such geometries, and will, furthermore, be used as a benchmark for the more general calculation of the Æ. The algebraic details of this derivation are given in Appendix A, and here we highlight the central steps. It is convenient to express $\lambda$ as a trapping parameter $k^2$, where the deeply trapped particles have $k=0$ and the barely trapped particles have $k=1$. This mapping is given by $\lambda = 1 + \epsilon (1 - 2 k^2)$, so the magnetic field may be written as

(2.58)\begin{equation} \lambda \hat{B} = 1 + \epsilon (1 - 2 k^2 - \cos \theta). \end{equation}

One can now perform the bounce-averaging integrals required for (2.37) in the $s$–$\alpha$ limit, resulting in

(2.59)\begin{equation} \hat{\omega}_\lambda =-\frac{\alpha}{2q^2} + 2 G_1(k) + 4s G_2(k) - \alpha G_3(k), \end{equation}

where we define

(2.60a)\begin{gather} G_1 = \frac{E(k)}{K(k)} - \frac{1}{2}, \end{gather}

(2.60b)\begin{gather}G_2 = \frac{E(k)}{K(k)} + k^2 - 1, \end{gather}

(2.60c)\begin{gather}G_3 = \frac{2}{3} \left[ \frac{E(k)}{K(k)} (2k^2 - 1) + 1 - k^2 \right], \end{gather}

where $K(k)$ and $E(k)$ are complete elliptic integrals of the first and second kind, respectively. The normalised bounce time, as given in (2.41) is equal to

(2.61)\begin{equation} \hat{g}_\epsilon^{1/2} = \frac{\sqrt{2}}{\rm \pi} K(k). \end{equation}

Finally, from (2.56) we see that $C_r = 1$ in this limit. The Æ now becomes a straightforward integral of known functions over $k$

(2.62)\begin{equation} \hat{A} = \frac{1}{2{\rm \pi}\sqrt{2}} \left(\frac{\rho_\text{g}}{R_0} \right)^2 \sqrt{\epsilon} \int_0^1 \mathrm{d} k^2 I_z(c_0(k),c_1(k)) \hat{\omega}_\lambda(k)^2 K(k), \end{equation}

which can efficiently be computed numerically.

3.1. Comparison with TGLF

Our first course of action is comparing Æ with turbulent energy-flux calculations in tokamak geometries, to verify its relation to turbulent transport in such geometries. At the moment, nonlinear gyrokinetic simulations are computationally too expensive for detailed parameters scans, and therefore we instead employ the quasilinear TGLF (trapped gyro-Landau fluid) code (Staebler, Kinsey & Waltz Reference Staebler, Kinsey and Waltz2007; Staebler & Kinsey Reference Staebler and Kinsey2010). Some key differences between the two models are highlighted before any comparison is made. The TGLF code computes the linear eigenmodes of a variety of instabilities, ion and electron temperature gradient modes, electromagnetic kinetic ballooning modes, as well as trapped-ion-modes and TEMs, and then applies a quasilinear saturation rule to accurately fit the fluxes from nonlinear gyrokinetic simulations. For quasineutrality purposes, TGLF requires the inclusion of at least one ion species. These are fundamental differences to the formulation of the Æ described in this work, which only accounts for the Æ of trapped electrons. Therefore, when setting up TGLF, care was taken to ensure the modelled turbulent energy-fluxes were as much as possible due to instabilities dominated by trapped electrons, using settings analogous to those used in recent gyrokinetic simulations in a similar regime (Proll et al. Reference Proll, Plunk, Faber, Görler, Helander, McKinney, Pueschel, Smith and Xanthopoulos2022). Given the lack of collisions in this regime, the expected dominant instabilities should be of the collisionless TEM variety. However, some other instabilities can also arise from interactions with the ion population. Thus, to ensure that the dominant instabilities in the TGLF simulations were as relevant as possible for our comparison, only contributions from modes propagating in the electron-diamagnetic direction were included, which excludes, for example, the ubiquitous mode (Coppi & Pegoraro Reference Coppi and Pegoraro1977), which propagates in the ion direction. Furthermore, we find that, for the scenarios considered in this work, adding an equally large electron temperature gradient to the density gradient, i.e. taking $\eta =1$, significantly decreased the number of non-TEM modes dominant in TGLF simulations, and as such we set $\eta$ to unity for the comparison. The recent SAT2 (Staebler et al. Reference Staebler, Belli, Candy, Kinsey, Dudding and Patel2021) quasilinear saturation rule for TGLF was used, as it includes the impact of plasma shaping on the quasilinear saturation (Staebler et al. Reference Staebler, Candy, Belli, Kinsey, Bonanomi and Patel2020). Although TGLF also uses a Miller parameterisation of the local equilibrium, we note that it does not use the same normalisation as Roach et al. (Reference Roach, Connor and Janjua1995) followed in this work, and care has been taken to convert between the two. We finally stress that the current model for $C_r$ is a fairly simple model, and that prediction can be refined using a more sophisticated model. This can, for example, be done by using some fitting function for $C_r$, where one finds the best-fit parameters which minimise the error between the energy flux and the prediction of Æ.

For the comparison we use the gyro-Bohm normalised energy fluxes computed by TGLF,

(3.1)\begin{equation} \hat{Q}_e = \frac{Q_e}{Q_\mathrm{GB}}, \end{equation}

where $Q_e$ is the electron energy flux from TGLF, and $Q_\mathrm {GB}$ is the gyro-Bohm energy flux. This is compared with the estimate of the gyro-Bohm normalised energy flux from Æ (Mackenbach et al. Reference Mackenbach, Proll and Helander2022, Reference Mackenbach, Proll, Wakelkamp and Helander2023a), which is

(3.2)\begin{equation} \hat{Q}_{A} \equiv C_Q \hat{A}^{\: 3/2}, \end{equation}

where the constant of proportionality is taken from the fit presented and was found to be $C_Q \approx 1.0 \times 10^3$. With such a power law, a linear correlation between $\hat {Q}_A$ and $\hat {Q}_e$ from nonlinear gyrokinetic simulations was found for pure density gradient-driven TEMs, which is different from the current comparison in which both the electron temperature and the density gradient drive the TEM ($\eta =1$). The data points in the comparison are chosen in order to verify that TGLF reproduces some trends that will be discussed in the following sections.

A comparison in the $(\kappa,\delta )$ and $(s,\alpha )$ planes is displayed in figure 3. One can see that there is good correspondence in trends: decreasing the magnetic shear and/or increasing the pressure gradient helps in reducing the energy flux, as does increasing the elongation. However, there are also differences visible between the two models for the energy flux, which are evident in the $(s,\alpha )$-plot. A clear discrepancy can be seen at high shear values $(s \approx 3)$, where the TGLF energy flux drops and the Æ estimate does not, and the Æ, furthermore, overestimates the energy flux at high shear. In the $(\kappa,\delta )$-plots the trends are well captured by Æ, with some differences. To further investigate the relationship between the two estimates of the energy flux, all the simulation data shown in figure 3 have been combined in a scatter plot shown in figure 4. In order to check consistency with previous findings, we have, furthermore, included the data points of Mackenbach et al. (Reference Mackenbach, Proll and Helander2022), which are nonlinear simulations in general geometries. Here, we see that there is a linear relationship for most of the data (we have added a red line with the expected linear relationship), although there exist data points that deviate more significantly from the linear relationship. There are various reasons why such a discrepancy may occur.

1. (i) There may be other instabilities present (though not necessarily dominant) that are not captured by the Æ of trapped electrons, such as the ubiquitous mode, or the universal instability (Romanelli Reference Romanelli1989; Helander & Plunk Reference Helander and Plunk2015; Landreman, Antonsen & Dorland Reference Landreman, Antonsen and Dorland2015; Costello et al. Reference Costello, Proll, Plunk, Pueschel and Alcusón2023). More generally, if there are instabilities present that do not derive their energy from trapped electrons, the current Æ-model is no longer expected to be an accurate measure.

2. (ii) The Æ length scale $C_r$ may vary more significantly for certain choices of equilibrium parameters and the current choice given in (2.57) may not be accurate.

3. (iii) Recall that Æ can be interpreted as an upper bound on the amount of energy that can be released. If the portion of the Æ that resides in stabilising modes deviates markedly (see, e.g. Lang, Parker & Chen Reference Lang, Parker and Chen2008; Hatch et al. Reference Hatch, Terry, Jenko, Merz, Pueschel, Nevins and Wang2011; Pueschel et al. Reference Pueschel, Faber, Citrin, Hegna, Terry and Hatch2016; Duff et al. Reference Duff, Faber, Hegna, Pueschel and Terry2022), one can reasonably expect that the data deviate more from the found relationship.

4. (iv) The TGLF's quasilinear saturated fluxes in both the $(\kappa,\delta )$ and $(s,\alpha )$ planes show occasional extreme outliers for small changes in input. The TGLF code has been extensively verified against a wide variety of nonlinear gyrokinetic simulations (although further validation for negative triangularity is currently being pursued), but the regime explored in this work is not the typical input space and could require separate verification.

5. (v) Although not present in the current set of simulation data, the Æ of trapped electrons will certainly cease to be an accurate model in situations where the trapped electrons play no role, such as in the case of a pure ion temperature gradient, and no gradients in electron temperature or density.

Figure 3. Comparison between codes, showing the correspondence between the estimate of the energy flux from Æ and TGLF. Panels (a,c,e) display the energy flux from TGLF, and panels (b,d,f) display the corresponding estimate from Æ. One can see agreement in trends, though some details differ. Panels (a,b) have a Miller vector $[\epsilon,\kappa,\delta, s_\kappa, s_\delta, \partial _r R_0, q, s, \alpha ]$ of $[1/3,\kappa,\delta, 0, 0, 0, 2, 0, 0]$, whereas panels (c,d) have $[1/3,\kappa,\delta, 0, 0, 0, 2, 1, 0]$, and panels (e,f) have a Miller vector of $[1/3,1,1/2, 0, 0, 0, 2, s, \alpha ]$. Panels (a,b) have $\hat {\omega }_n=3$; panels (c,d) and (e,f) have $\hat {\omega }_n=6$, and $\eta =1$ in all panels. The white masked-out regions have a dominant instability which is not in the electron direction, and as such we filter them out. Finally, note that the colour bars are the same scale in each column.

Figure 4. A scatter plot showing the relation between the two estimates of the nonlinear energy flux. The red dashed line is shows the expected linear relationship, $x=y$. The plot consists of $N=1171$ points. The grey points have some transparency, and as such darker regions arise due to a high density of points. Furthermore, we have added simulation data from the gyrokinetic code GENE (gyrokinetic electromagnetic numerical experiment) (Jenko, Dorland & Hammett Reference Jenko, Dorland and Hammett2001), also presented in Mackenbach et al. (Reference Mackenbach, Proll and Helander2022), as blue markers.

We stress that the scatter does mean that predictions may be faulty if one lies within the scatter of the fit. However, seeing that general trends are well captured by Æ, it may serve as a useful estimate for transport and trends at low computational cost (Æ calculations are roughly a factor $50$ faster than the presented TGLF calculations).

3.2. The $s$–$\alpha$ geometry

We now shift our attention to the behaviour of Æ on the various free parameters found in tokamak equilibria. Recalling that we have derived two Æ expressions, one for any Miller geometry and one for the large-aspect ratio limit, let us start by investigating the latter. A plot of the Æ calculated from (2.62) is given in figure 5 as a function of magnetic shear and pressure gradient. We note that the ranges for $s$ and $\alpha$ are not meant to represent realistically attainable values here, instead, we are more interested in the general structure of the Æ over the domain. There are several interesting features visible in the figure. Even in this simplest model, the Æ exhibits rich structure over the $s$–$\alpha$ plane. More precisely, Æ is large when $s$ and $\alpha$ are comparable, $s\sim \alpha$, and is otherwise much smaller, particularly when the absolute value of one of these quantities is large. These findings are consistent with previous investigations (Rosenbluth & Sloan Reference Rosenbluth and Sloan1971; Dagazian & Paris Reference Dagazian and Paris1982; Connor et al. Reference Connor, Hastie and Martin1983; Kessel et al. Reference Kessel, Manickam, Rewoldt and Tang1994; Rettig et al. Reference Rettig, Peebles, Doyle, Burrell, Greenfield, Staebler and Rice1997; Strait et al. Reference Strait, Casper, Chu, Ferron, Garofalo, Greenfield, La Haye, Lao, Lazarus and Miller1997; Kinsey, Waltz & Candy Reference Kinsey, Waltz and Candy2006). It is also interesting to note that the precise reduction in Æ depends on the drive: for a pure electron temperature gradient, significant positive shear is more helpful in reducing Æ, while Æ driven by a pure density gradient benefits more from negative shear.

Figure 5. The Æ of a large-aspect-ratio circular tokamak, as a function of magnetic shear $s$ and pressure gradient $\alpha$. The plots have been generated using $q=2$. Panel (a) has $\hat {\omega }_n=3$ and $\eta =0$, whereas (b) has a pure electron temperature gradient, i.e. $\hat {\omega }_n=0$ and $\hat {\omega }_n \cdot \eta = 3$.

Since (2.62) can be integrated numerically to high precision, it serves as a useful benchmark for the more general Æ of (2.44). Accordingly, we have compared the Æ in the large-aspect-ratio limit with circular flux surfaces using a code that solves (2.44). This comparison is shown in Appendix B, and we find that the codes agree.

3.3. Miller geometry

We now leave the realm of the $s$–$\alpha$ limit and venture into shaped, finite-aspect-ratio equilibria. Our first step is to investigate the dependence on magnetic shear and pressure gradient for a range of different Miller vectors, and the results are shown in figure 6. Here we see similar trends as in § 3.2: negative shear and large $\alpha$ tend to be especially stabilising for a pure density gradient. However, it is also clear that the magnitude and precise contours depend strongly on the chosen Miller vector, as defined in (2.54). For example, it can be seen that lowering the safety factor is stabilising, since Æ is reduced over a large region of the $s$–$\alpha$ plane as one compares figure 6(a) with figure 6(b). In figure 6(c) the elongation has been reduced produce a ‘comet’-type configuration ($\kappa < 1$, i.e. a horizontally elongated tokamak, see figure 2), which can increase the magnitude of the Æ, and the stabilising effects of $s$ and $\alpha$ become less pronounced. Finally, in figure 6(d) the sign of the triangularity has been reversed to become negative. Although the shape of the contours remains largely unchanged, the peak in Æ is changed to higher $\alpha$ and lower $s$, indicating that negative triangularity can be particularly beneficial in high-shear discharges with a modest value for $\alpha$. In a more general sense, when changing any of the parameters significantly, one should expect that the precise shape and magnitude of the contours will change.

Figure 6. Dependence of Æ on magnetic shear $s$ and pressure gradient $\alpha$. In (a), the Miller vector $[\epsilon,\kappa,\delta, s_\kappa, s_\delta, \partial _r R_0, q, s, \alpha ]$ is set to $[ 1/3, 3/2, 1/2, 0, 0, 0, 2, s,\alpha ]$, and other panels have the same vector with one parameter changed. In (b) the safety factor $q$ is reduced from $2$ to $1$, for (c) the elongation $\kappa$ decreases from $3/2$ to $1/2$, and in (d) the sign of triangularity $\delta$ is changed from $1/2$ to $-1/2$. All panels have $\hat {\omega }_n=1$ and $\eta = 0$.

With this important caveat in mind, let us investigate the influence of geometry on the Æ. To do so, we display the dependence on $\kappa$ and $\delta$ for various Miller vectors in figure 7. Several interesting general trends can be observed. First, note that increasing the elongation beyond $\kappa = 1$ generally decreases the Æ in all Miller vectors considered here, although the precise effect depends on the triangularity. Second, we see that it is not true in general that positive or negative triangularity is always stabilising; it depends on the other Miller parameters. Third, we see that tokamaks with $\kappa <1$ and $\delta < 0$, often referred to as (negative) comet cross-sections tokamaks (Kesner, Ramos & Gang Reference Kesner, Ramos and Gang1995), show a reduction in Æ in figure 7(b,d), at least for the pure density gradient considered here. This is perhaps unsurprising, since such tokamaks are close to having the maximum-$\mathcal {J}$ property as shown by Miller et al. (Reference Miller, Chu, Dominguez and Ohkawa1989). Since Æ measures deviations from the maximum-$\mathcal {J}$ property, it is thus expected that these configurations perform well in terms of Æ.

Figure 7. The effect of the geometry on Æ. In (a) the Miller vector $[\epsilon,\kappa,\delta, s_\kappa, s_\delta, \partial _r R_0, q, s, \alpha ]$ is set to $[ 1/3,\kappa, \delta, 0, 0, 0, 2, 1,0 ]$ with $\hat {\omega }_n=1$ and $\eta = 0$, and all other panels have the same Miller vector with one parameter changed. The contours shown in (b) have a higher inverse aspect ratio as $\epsilon$ is increased from $1/3$ to $2/3$, in (c) the pressure gradient $\alpha$ is increased from $0$ to $1/2$, and for (d) we have decreased the shear from $1$ to $0$.

Investigating the plots in detail, in figure 7(a) one sees that negative triangularity is beneficial for $\kappa > 1$ as can be seen by the reduction in Æ. In the following sections, we shall see that this is a consequence of the positive magnetic shear chosen. Next, note that doubling the inverse aspect ratio, as is done when going from figure 7(a) to figure 7(b), has a stabilising effect. Naively, one would expect that doubling the inverse aspect ratio would increase the Æ by roughly a factor $\sqrt {2} \approx 1.4$, due to the factor $\sqrt {\epsilon }$ in (2.45). However, going from figure 7(a) to figure 7(b) we see a decrease of the maximum Æ by some $15\,\%$. This is likely due to the fact that, in a small-aspect-ratio device, magnetic field lines spend most of their time (or more precisely, arclength) on the inboard side of the tokamak (Helander & Sigmar Reference Helander and Sigmar2005). There, $\omega _\lambda$ tends to be opposite to the drift wave and therefore these orbits do not contribute to the Æ for a pure density gradient. It is also interesting to note that negative triangularity no longer exhibits a reduction in Æ as the aspect ratio is significantly decreased, in accordance with the findings of Balestri, Ball & Coda (Reference Balestri, Ball, Coda, Garcia-Munoz and Viezzer2023). Going from figure 7(a) to figure 7(c) the pressure gradient is increased from $\alpha = 0$ to $1/2$. With this introduction of pressure gradient, it can be seen that positive triangularity shows a decrease in Æ, where negative triangularity does not. Finally, figure 7(d) has the magnetic shear reduced from $s=1$ to $s=0$ as compared with figure 7(a), which drastically changes the picture. Most importantly, we see that the lack of this positive magnetic shear results in negative triangularity no longer being stabilising. We find that the results change somewhat if one instead imposes a pure electron temperature gradient (not shown here), though the basic trends remain intact.

All in all, we conclude from these results that the Æ is very sensitive to equilibrium parameters, including quantities not investigated here such as $q$, $s_\kappa$ and $\partial _r R_0$. This sensitivity is perhaps reassuring: gyrokinetic turbulence has long been known to be strongly dependent on equilibrium parameters and even slight nudges can drastically change the picture (a sentiment perhaps best captured by the old Dutch expression wie het kleine niet eert, is het grote niet weerd). We seem to reproduce a similar sensitivity in this Æ-model for trapped electrons. This sensitivity becomes especially clear when investigating the dependence of Æ on triangularity, which we shall discuss in the next section.

3.4. When is negative triangularity beneficial?

As hinted at in the previous section, it is not possible to make a general statement about the effect of negative triangularity on Æ; its possible benefit depends strongly on other parameters describing the equilibrium. We can, however, find trends, and in order to do so we define the following fraction:

(3.3)\begin{equation} \varDelta = \frac{\hat{A}(\delta =-0.1)}{\hat{A}(\delta =+0.1)}, \end{equation}

where $\delta = \pm 0.1$ is chosen to represent a typical experimentally realisable range of parameters. This fraction can be interpreted as the factor by which the Æ changes upon switching from positive to negative triangularity, where $\varDelta <1$ implies a reduction in Æ. We present an investigation of $\varDelta$ and its dependencies in figure 8. We see two clear trends that seem to be robust for tokamaks with $\kappa > 1$. Firstly, as noted in the previous sections, in figure 8(a,b) we see that negative triangularity tends to be especially stabilising for configurations with significant positive shear. Similar conclusions were made by Merlo & Jenko (Reference Merlo and Jenko2023), who found that the turbulent energy flux in gyrokinetic simulations follows the same trend for TEM-driven turbulence: only for sufficiently high positive shear is a decrease in energy flux found at negative triangularity. Increasing $\alpha$ tends to push the $\varDelta = 1$ line (in the plot this is the $\log \varDelta = 0$ line) to even higher values of shear, implying that a significant pressure gradient may make negative triangularity less desirable. Secondly, in figure 8(c,d) we note that negative triangularity can be beneficial in situations where the gradient is small, such as in the core. The dependence on $\eta$ is non-trivial; at small density gradients a non-zero value of $\eta$ can make negative triangularity beneficial. As in the previous sections, the results here depend on the Miller vector and are not meant to serve as a quantitative measure for core and edge transport. However, we have found that the presented trends tend to be robust as long as $\kappa > 1$ and thus do have qualitative value. We finally note that a more comprehensive model of the effect of negative triangularity should likely take collisions, impurities and global effects into account (Merlo et al. Reference Merlo, Fontana, Coda, Hatch, Janhunen, Porte and Jenko2019, Reference Merlo, Huang, Marini, Brunner, Coda, Hatch, Jarema, Jenko, Sauter and Villard2021).

Figure 8. Various plots showcasing the dependencies of $\log (\varDelta )$, where $\varDelta$ is defined as $\hat {A}(\delta =-0.1)/\hat {A}(\delta =+0.1)$, on various equilibrium parameters. Panels (a,c) and (b,d) have different Miller vectors $[\epsilon,\kappa,\delta, s_\kappa, s_\delta, \partial _r R_0, q, s, \alpha ]$, which are meant to be representative of the edge and the core. Panels (a,c) have a core-like Miller vector of $\boldsymbol {M}_\mathrm {core}=[1/100,3/2,\delta,0,0,0,1,s,\alpha ]$. Panels (b,d) have an edge-like Miller vector of $\boldsymbol {M}_\mathrm {edge}=[1/3,3/2,\delta,1/2,0,-1/2,3,s,\alpha ]$ Finally, (a) has $\hat {\omega }_n=1/2$ and $\eta =0$, (b) has $\hat {\omega }_n=2$ and $\eta =0$, (c) has $s=\alpha =0$ and (d) has $s=2$ and $\alpha =1/2$.

From these results we infer that negative triangularity is expected to be especially beneficial in the core of the plasma, where gradients are necessarily small. It is not clear if the benefit extends to the edge: only with significant positive shear does negative triangularity become beneficial here as well. One should also keep in mind that $\varDelta$ measures the effect of going to negative triangularity while keeping all other parameters fixed. A more complete investigation would, for example, compare experimental equilibria with positive and negative triangularity, or use a global magnetohydrodynamics-equilibrium code to find consistent profiles. We do not attempt such an investigation here, but we note that our mathematical framework would readily allow for such a comparison. We finally remark that the above results may seem counter-intuitive as negative triangularity is often thought to automatically imply TEM stabilisation, since the bounce points of most trapped particles reside on the inboard side of the torus, where the magnetic curvature should be favourable. Consequently, it is often argued that the bounce-averaged drift is such that TEMs are stabilised. Upon calculation of (2.37), we find no such stabilisation, however, as explained further in Appendix C.

3.5. Gradient-threshold-like behaviour

Our next step is to investigate the dependence of Æ on the gradient strength $\hat {\omega }_n$. From (2.8), one can show that there are two distinct scalings (Mackenbach et al. Reference Mackenbach, Proll, Wakelkamp and Helander2023a). In a strongly driven regime, one finds that the Æ scales linearly with the gradient strength $\hat {\omega }_n$. For a weakly driven regime one can expand around small $\hat {\omega }_n$, and one finds that the Æ scales with the gradient strength as $A \propto \hat {\omega }_n^{3}$:

(3.4)\begin{equation} \hat{A} \propto \begin{cases} \hat{\omega}_n & \text{if } |\hat{\omega}_n | \gg 1, \\ \hat{\omega}_n^3 & \text{if } |\hat{\omega}_n | \ll 1. \end{cases} \end{equation}

These scalings are reminiscent of gradient-threshold (or critical gradient) type behaviour (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000). Gradient thresholds are signified by a sudden decrease in energy flux when decreasing the gradient below some threshold value. The aforementioned scaling behaviour of the Æ is displayed in figure 9 which similarly shows a rapid decrease below some threshold value. In figure 9(b) we estimate a critical-threshold-like quantity from Æ, by fitting a straight line to the strongly driven regime, i.e. we find the best-fit parameters $a_0$ and $a_1$ in the formula

(3.5)\begin{equation} \hat{A} = a_0 + a_1 \hat{\omega}_n, \end{equation}

with $\hat {\omega }_n \gg 1$. The gradient threshold, denoted by $\hat {\omega }_c$, is then defined as the interception with the abscissa, hence

(3.6)\begin{equation} \hat{\omega}_c \equiv- \frac{a_0}{a_1}. \end{equation}

One could, of course, use different definitions for $\hat {\omega }_c$, for example, one could define the intersection point between the two straight lines on the log–log plot of figure 9 as $\hat {\omega }_c$. However, we have found that the definition of (3.6) has several benefits: it is computationally cheaper, less prone to numerical noise and seems to behave more smoothly. Other attempted definitions show the same trends.

Figure 9. Example of dependence of Æ on gradient strength. Two scalings are found in (a). In (b) we define a gradient threshold by fitting a straight line to the strongly driven regime, and finding its $\hat {\omega }_n$ interception with the abscissa.

We illustrate how $\hat {\omega }_c$ varies as a function of various equilibrium parameters in figure 10. Note that figure 10(a) has the same Miller vector as figure 6(a), and figure 10(b) has the same Miller vector as figure 7(a). Focussing on figure 10(a), an interesting trend is that increasing shear tends to increase $\hat {\omega }_c$ linearly, and $\hat {\omega }_c$ tends to plateau for low shear to some value. This is similar to the findings of Jenko et al. (Reference Jenko, Dorland and Hammett2001), though their investigation focusses on electron-temperature gradient turbulence. It is also interesting to note that, in addition to the reduction in Æ in the negative-triangularity configuration, it also benefits from a high critical gradient, which is in line with the findings of Merlo et al. (Reference Merlo, Brunner, Sauter, Camenen, Görler, Jenko, Marinoni, Told and Villard2015). This effect becomes even more pronounced as one increases the shear, which, furthermore, reduces the Æ in the negative triangularity configuration. This implies that negative triangularity may be beneficial in a different sense: since the critical gradient estimated from Æ is higher in negative triangularity geometries, the profiles may be able to sustain much higher gradients and thus higher core density/temperature.

Figure 10. The gradient threshold as a function of equilibrium parameters. In (a) the Miller vector $[\epsilon,\kappa,\delta, s_\kappa, s_\delta, \partial _r R_0, q, s, \alpha ]$ is $[1/3,3/2,1/2,0,0,0,2,s,\alpha ]$, and (b) has $[1/3,\kappa, \delta,0,0,0,2,1,0]$. In all plots $\eta = 0$.

3.6. Tokamak optimisation

In this section we aim to find Æ-optimised tokamaks for a certain set of equilibrium parameters, at fixed gradients ($\hat {\omega }_n=1$ and $\eta$=0). To this end, we choose to optimise over $\kappa$ and $\delta$ while keeping all other parameters fixed. In order to find somewhat realistic solutions, we restrict ourselves to a bounded optimisation space, namely

(3.7a,b)\begin{equation} \kappa \in (1/2,2),\quad \delta \in (-1/2,1/2). \end{equation}

The SHGO (simplicial homology global optimization) algorithm from Endres, Sandrock & Focke (Reference Endres, Sandrock and Focke2018) is ideally suited for finding the global minimum in this low-dimensional bounded parameter space and is also available in SciPy. Finally, we shall vary magnetic shear and $\alpha$, and investigate its effect on the global minimum found.

The results are displayed in figure 11, where the optimal values of Æ, $\kappa$ and $\delta$ values are displayed as a function of $s$ and $\alpha$. For a visual aid of the shape of the cross-sections, we refer to figure 2. It can be seen that both the optimal triangularity and elongation tend to be in the corners of the optimisation domain, and hence one should expect that these results are strongly dependent on this domain. Firstly, we see that vertically elongated tokamaks tend to be beneficial for all parameters considered here. It is, furthermore, interesting to note that the negative triangularity solution tends to be optimal whenever there is significant shear and the pressure gradient is not too large, which is consistent with the findings of § 3.4.

Figure 11. Global Æ-minimising solutions as a function of $s$ and $\alpha$. Panel (a) showcases the Æ of the optimal solution, (b) displays the elongation and (c) shows the triangularity. Generated with a Miller vector $[\epsilon,\kappa,\delta, s_\kappa, s_\delta, \partial _r R_0, q, s, \alpha ]$ of $[1/3,\kappa,\delta,0,0,-1/2,3,s,\alpha ]$. In all panels $\hat {\omega }_n = 1$ and $\eta = 0$.

From this plot, an important conclusion can be drawn: there is no such thing as a single ‘optimal’ solution. The global minimum depends sensitively on other equilibrium parameters, such as shear and pressure gradient, which are, in turn, determined by the profiles of the safety factor, density and temperature. Therefore, if one is interested in finding an Æ-optimised tokamak, one should take care when choosing the profiles. One could also choose to let the profiles be part of the optimisation by describing them with some number of free parameters and constraints (e.g. one could use a fixed number of Fourier modes on top of a profile and optimise for the mode amplitudes). In reality, the profiles are themselves set by equilibrium conditions, making a self-consistent optimisation highly non-trivial. A more consistent investigation could perhaps solve this by coupling the current Æ-model to a transport solver, which would calculate self-consistent profiles.

3.7. Existence of solutions with high gradients yet low Æ

In this section, we investigate how this Æ model may relate to the suppression of TEMs when the density gradient is increased. To do so, we note several interesting properties that arise as one increases this gradient. First, the normalised pressure gradient $\alpha$ scales linearly with the density gradient (assuming a constant ratio of the poloidal magnetic field pressure to the thermal pressure, which, for example, occurs if one is operating at a fixed $\beta$-limit). The shear depends on the pressure gradient, as such a gradient drives the bootstrap current, which in turn changes the rotational transform profile. The bootstrap current density has an off-axis maximum in realistic scenarios, and such an off-axis maximum can locally lower the shear. This is most readily seen by inspecting the expression for shear in a large-aspect-ratio, circular tokamak, which depends on the current density profile $j(r)$ as

(3.8a,b)\begin{equation} s(r) = 2 \left( 1 - \frac{j(r)}{\bar{\jmath}(r)} \right);\quad \bar{\jmath}(r) = \frac{2}{r^2} \int_0^r x j(x)\,\mathrm{d}x, \end{equation}

where $\bar {\jmath }$ measures the average current density inside the radius $r$. From this expression, it is clear that for current density profiles that peak at $r=0$, the shear is always positive. An off-axis maximum, supplied by the bootstrap current, can cause a locally lower shear. Hence, as one raises $\hat {\omega }_n$ one simultaneously increases $\alpha$ and decreases $s$. To estimate the magnitude of the effect of the bootstrap current on the shear, we note that the bootstrap current is proportional to the density and temperature gradients, and thus to the pressure gradient

(3.9)\begin{equation} j_b \approx \,j_{b,0}\alpha(r). \end{equation}

This is an approximation since the different transport coefficients relating the bootstrap current to the various gradients are not identical (Helander & Sigmar Reference Helander and Sigmar2005), but we ignore this minor complication. We, furthermore, write the total current density as $j = j_b + j_e$, where $j_e$ is the equilibrium current, and assume $j_b \ll j_e = j_{e,0} \hat {\jmath }(r)$. To first order in the smallness of the bootstrap current, (3.7a,b) then gives

(3.10)\begin{equation} s \approx 2 - \frac{r^2 \hat{\jmath}(r)}{\displaystyle\int_0^\rho x\hat{\jmath}(x)\,\mathrm{d}x} \left(1 + \frac{j_{b,0}}{j_{e,0}} \left[\frac{\alpha(r)}{\hat{\jmath}(r)} - \frac{\displaystyle\int_0^r x \alpha(x)\,\mathrm{d}x}{\displaystyle\int_0^r x \hat{\jmath}(x)\,\mathrm{d}x}\right]\right). \end{equation}

Finally, following Miyamoto (Reference Miyamoto2005) we estimate the ratio $j_{b,0}/j_{e,0}$ as

(3.11)\begin{equation} \frac{j_{b,0}}{j_{e,0}} \approx 0.3 \langle \beta_p \rangle \sqrt{\epsilon}, \end{equation}

where $\beta _p$ is the local ratio of the thermal pressure over the poloidal magnetic field pressure, and the angular brackets denote a volume average. We shall take $j_{b,0}/j_{e,0}$ to be of the order of $10\,\%$, implying that the shear may change as $\mathrm {d} s /\mathrm {d} \alpha \sim s/10$. Finally, one can relate the pressure gradient to $\hat {\omega }_n$ as

(3.12)\begin{equation} \alpha = \epsilon \beta_p \left( 1 + \eta + \eta_i \right)\hat{\omega}_n, \end{equation}

where $\eta _i=\partial _r \ln T_i / \partial _r \ln n$, with $T_i(r)$ being the ion temperature. We assume that the factor $\epsilon \beta _p (1 + \eta + \eta _i) \sim 0.1$, so that $\mathrm {d} \hat {\omega }_n / \mathrm {d} \alpha \sim 10$.

We illustrate the competing effects of the density gradient, pressure gradient, and shear in figure 12. In figure 12(a,b) we see various isocontours of the Æ in $(\hat {\omega }_n,s,\alpha )$-space, where figure 12(a) has positive triangularity and figure 12(b) has negative triangularity. It is especially interesting to note that in figure 12(a) there are paths in parameter space in which $\hat {\omega }_n$ increases but the Æ decreases. These paths generally require that, as the density gradient increases, the pressure gradient should also increase and the shear should decrease. As we have argued, these trends are indeed found in tokamak discharges. One such path is indicated in figure 12(a) as a blue line. Importantly, the blue line has

(3.13a,b)\begin{equation} s = 4 \left(1 - \frac{\alpha}{8} \right),\quad \hat{\omega}_n = 10 \cdot \alpha \end{equation}

which is the correct order of magnitude for both $\mathrm {d}s/\mathrm {d}\alpha$ and $\mathrm {d}\hat {\omega }_n/\mathrm {d}\alpha$. Figure 12(b) exhibits drastically different features. Planes of constant Æ tend to lie parallel to planes of constant $\hat {\omega }_n$, indicating that not much stabilisation is possible by changing the shear or the pressure gradient: the Æ rises when $\hat {\omega }_n$ is increased. In figure 12(b), we again plot a line along the direction of increasing $\alpha$ and decreasing magnetic shear in red. Finally, note that for $s_\delta$ we have used the estimate from Miller et al. (Reference Miller, Chu, Greene, Lin-Liu and Waltz1998), $s_\delta \approx \delta /\sqrt {1 - \delta ^2}$.

Figure 12. Panels (a,b) showcase isocontours of the Æ as a function of $(\hat {\omega }_n,\alpha,s)$, where (a,b) have positive and negative triangularity,respectively. In both (a,b) a straight line is drawn which has increasing $\alpha$ and decreasing $s$ with increasing $\hat {\omega }_n$, and the projection of the line onto the grid-planes is shown as a dashed line. In (c), the Æ along both the blue line of (a) and the red line of (b) is plotted as a function $\hat {\omega }_n$. These plots were generated with a Miller vector $[\epsilon,\kappa,\delta, s_\kappa, s_\delta, \partial _r R_0, q, s, \alpha ]$ of $[1/3,3/2,\delta, 1/2, \delta /\sqrt {1-\delta ^2}, -1/2, 3, s, \alpha ]$, and $\eta = 0$.

In figure 12(c) we display the Æ along the blue and red lines given in figure 12(a,b) as a function of the density gradient. Note that the positive-triangularity case exhibits a distinct maximum, with low Æ both to the left and right of the peak. One could interpret the existence of the latter as two distinct low-transport regimes: one with low gradients, and one with high gradients (which also has decreased magnetic shear and increased $\alpha$). It is, furthermore, interesting to note that the negative-triangularity tokamak rises to far higher values in terms of Æ and does not seem to drop back down to low levels along the chosen domain. Hence one could perhaps conclude that reaching a low-transport state with high gradients is not feasible in a negative-triangularity discharge. This is in line with findings of Saarelma et al. (Reference Saarelma, Austin, Knolker, Marinoni, Paz-Soldan, Schmitz and Snyder2021) and Nelson, Paz-Soldan & Saarelma (Reference Nelson, Paz-Soldan and Saarelma2022), where the H-mode was found to be inaccessible in negative-triangularity tokamaks on the basis of the ballooning instability, though the physical reason is of course different. This rise in Æ in negative triangularity is perhaps unsurprising given that we have found that negative triangularity is stabilising in cases with significant positive shear, a weak pressure gradient, and a slight density gradient, exemplified in figures 8 and 11. Since, along the chosen path shear decreases and $\alpha$ increases with increasing density gradient, which is opposite to what is stabilising for negative-triangularity tokamaks, we see a sharp increase in Æ. It may be feasible, however, to have a significant reduction in transport by tailoring the $q$-profile in such a way that negative triangularity becomes favourable, which likely implies significant positive shear. With such a reduction in Æ, one could perhaps enjoy much improved transport whilst staying in an L-mode-like regime. The parameters described in Marinoni et al. (Reference Marinoni, Austin, Hyatt, Walker, Candy, Chrystal, Lasnier, McKee, Odstrčil and Petty2019) do seem to meet such requirements, especially near the edge where the reduction in transport seems greatest as compared with the positive triangularity case.

A more comprehensive investigation, which shall be undertaken in a future publication, would self-consistently calculate the bootstrap current which would give precise paths in $(\hat {\omega }_n,\alpha,s)$-space. However, given the nature of the isocontours in this three-dimensional space, we expect the observed trends to be robust, as long as the path has the correct general dependencies (i.e. decreasing shear and increasing $\alpha$ with increasing density gradient).