2.1. Expanding the second adiabatic invariant

Let us write $\mathcal {J}_\parallel$ explicitly as a function of the field line following coordinate $\ell$ (where we have taken the particle mass $m=1$),

(2.5)\begin{equation} \mathcal{J}_\parallel{=} \sqrt{2 H} \int \sqrt{1 - \lambda \hat{B}(\ell)} \, \mathrm{d} \ell. \end{equation}

Here we have introduced the pitch-angle $\lambda = \mu B_0 / H$ (using for the first adiabatic invariant $\mu =(2H-v_\parallel ^2)/B$), which distinguishes between different trapped particles (deeper or more shallowly trapped), and we have normalised the magnetic field by some reference field strength $\hat {B} = B/B_0$. The integral is taken between bounce points, i.e. between points at which $\hat {B}=1/\lambda$. We are assuming there is no electric field within each flux surface and, as such, no electrostatic potential appears in (2.5).

It will prove convenient to express the field line following coordinate $\ell$ in terms of Boozer angles (Boozer Reference Boozer1981). We define for that purpose the helical angle,

(2.6)\begin{equation} \chi\stackrel{\cdot}{=}\theta-N\varphi, \end{equation}

where $N$ is an integer that defines a helical angle, foreseeing the application to quasisymmetric devices with a helical symmetry. Using this angular parametrisation, the Boozer-coordinate Jacobian $\mathcal {J}=B_\alpha (\psi )/B^2$, where $B_\alpha =G+\iota I$ (in the standard Boozer notation (Boozer Reference Boozer1981; Helander Reference Helander2014)), and defining $\bar {\iota }=\iota -N$, the second adiabatic invariant in these coordinates may now be written as

(2.7)\begin{equation} \mathcal{J}_\parallel{=} \frac{\sqrt{2H}}{\bar{\iota}} \frac{B_\alpha}{B_0} \int \frac{1}{\hat{B}} \sqrt{1 - \lambda \hat{B}} \, \mathrm{d}\chi. \end{equation}

It is crucial to note that $\mathcal {J}_\parallel$ depends directly on the magnetic field magnitude along a field line, with minimal involvement of other geometric elements. This simplifies the calculation of $\mathcal {J}_\parallel$ ostensibly when considering quasi- and axisymmetric configurations, and makes them rather analogous.

To construct a representative $\hat {B}$, we resort to the inverse-coordinate near-axis framework (Garren & Boozer Reference Garren and Boozer1991b; Landreman & Sengupta Reference Landreman and Sengupta2019), in which the asymptotic form of $\hat {B}$ near the axis is rather simple. We will employ essential results from the near-axis theory of quasisymmetric equilibrium fields as needed without re-deriving them and refer to the literature for details (Garren & Boozer Reference Garren and Boozer1991b; Landreman & Sengupta Reference Landreman and Sengupta2019). That way, and to second order in the distance from the magnetic axis, $r=\sqrt {2\psi /B_0}$, we write $\hat {B}$, following Garren & Boozer (Reference Garren and Boozer1991b, (1)) or Landreman & Sengupta (Reference Landreman and Sengupta2019, (2.15)), and noting that this behaviour goes beyond the particular form of equilibrium assumed (Rodríguez, Sengupta & Bhattacharjee Reference Rodríguez, Sengupta and Bhattacharjee2022), as

(2.8)\begin{equation} \hat{B}=\bar{B}+ \delta B(\varphi,\chi), \end{equation}

where we have separated

(2.9)\begin{equation} \bar{B} = 1 - r \eta \cos \chi \end{equation}

and the second order

(2.10)\begin{equation} \delta B= r^2\left(B_{20}+B_{22}^C\cos2\chi+B_{22}^S\sin2\chi \right). \end{equation}

In the ideal quasi- or axisymmetric limit, all parameters ($\eta$, $B_{20}$, $B_{22}^C$ and $B_{22}^S$) are constants instead of functions of the toroidal angle $\varphi$, and for stellarator symmetry $B_{22}^S=0$. From here on, we shall assume that this condition is satisfied exactly. Note that in practice, this condition is only achieved approximately: the near-axis expansion generally fails to find exactly quasisymmetric solutions for equilibria at second order in $r$ (unless exactly axisymmetric fields are considered). This is commonly referred to as the Garren–Boozer overdetermination problem (Garren & Boozer Reference Garren and Boozer1991a), and results from a clash between the symmetry and the equilibrium (Rodriguez & Bhattacharjee Reference Rodriguez and Bhattacharjee2021; Rodríguez et al. Reference Rodríguez, Sengupta and Bhattacharjee2022). In that sense, the idealised framework is only an approximation, but it will prove useful and, as we shall see, practical in approximating quasisymmetric configurations. The constants that define $|\boldsymbol {B}|$ can then be interpreted as parameters that describe different configurations. In fact, with these parameters together with a magnetic axis shape, the field and flux surfaces may be constructed explicitly (Landreman & Sengupta Reference Landreman and Sengupta2019). Thus, expressing the dependence of the second adiabatic invariant on these parameters will provide a direct link to the configuration and its distinguishing features. A note of caution: although we are considering the asymptotic limit in the distance to the magnetic axis, we cannot approach it closer than the gyroradius-related ‘potato-orbit’ size (Helander & Sigmar Reference Helander and Sigmar2005, Ch. 7) without violating the ‘thin orbit limit’ of our $\mathcal {J}_\parallel$ calculation.Footnote ¹

The way that we have grouped the terms in (2.8) might be unexpected, given that we have mixed asymptotically unequal terms: the constant on-axis field ($r=0$) with the first-order variation. However, since we are interested in trapped particles, variation in $|\boldsymbol {B}|$ along the field line is necessary, otherwise trapped particles would not exist. Therefore, in our perturbative description of the problem, we must take $\bar {B}$ to constitute the leading order magnetic field magnitude. It would then appear natural to write

(2.11)\begin{equation} \mathcal{J}_\parallel{=}\sqrt{2H}\frac{B_\alpha}{\bar{\iota}B_0}\int_{\chi_b^-}^{\chi_b^+}\frac{\sqrt{1-\lambda(\bar{B}+ \delta B)}}{\bar{B}+\delta B}\,\mathrm{d}\chi, \end{equation}

where for every $\lambda$, the bounce points $\chi _b^\pm$ (left and right) of the integral are given by

(2.12)\begin{equation} \frac{1}{\lambda}=\bar{B}(\chi_b^\pm)+\delta B(\chi_b^\pm), \end{equation}

and attempt to expand it in powers of $\delta B$ (i.e. expand around smallness of $\delta B$). There are however two important sources of conflict in doing so. First, and formally, evaluating $\sqrt {1-\lambda \bar {B}}$ at the bounce points $\chi _b^\pm$ can lead to imaginary contributions near these points without additional careful consideration of the bounce points. Second, physically, there are also issues related to the behaviour of classes such as barely trapped particles, which under an infinitesimal perturbation may undergo a finite (non-infinitesimal) change. This translates into diverging expressions in the perturbative construction.

To deal with these issues consistently, we start by defining a correction field $\delta B_b (\lambda )$, defined to be the perturbed field $\delta B(\chi )$ evaluated, for a given particle class $\lambda$, at the bouncing points, see (2.12). We shall assume, for simplicity, that the device is stellarator symmetric about the bottom of the magnetic well (i.e. we ignore the $B^S_{22}$ component), so that $\delta B_b (\lambda )$ is unique (i.e. it does not have left and right values). Introducing this correction, let us rewrite $\mathcal {J}_\parallel$ for a stellarator symmetric field in the following form:

(2.13)\begin{equation} \mathcal{I}(\epsilon)=\sqrt{2H}\frac{B_\alpha}{\bar{\iota}B_0}\int_{-\chi_b}^{\chi_b}\frac{\sqrt{1-\lambda[(\bar{B}+\delta B_b)+ \epsilon(\delta B-\delta B_b)]}}{\bar{B}+\epsilon\delta B}\,\mathrm{d}\chi, \end{equation}

so that $\mathcal {J}_\parallel =\lim _{\epsilon \rightarrow 1^-} \mathcal {I}(\epsilon )$. Expressing the integral in this form ensures that the integrand upholds positive definiteness for all $\epsilon \in [0,1)$, evading the issue of the integrand becoming imaginary. This furthermore circumvents the need to expand the bounce points of the integral. Expressed in this form, the integral may now be Taylor expanded in $\epsilon$,

(2.14)\begin{equation} \mathcal{I} \approx \underbrace{\mathcal{I}(0)}_{\stackrel{\cdot}{=}\mathcal{J}_\parallel^{(0)}} + \epsilon\underbrace{\left.\frac{\partial \mathcal{I}}{\partial \epsilon}\right|_{\epsilon=0}}_{\stackrel{\cdot}{=}\mathcal{J}_\parallel^{(1)}}, \end{equation}

where we shall take $\epsilon \rightarrow 1$ in the final answer so that $\mathcal {J}_\parallel \approx \mathcal {J}_\parallel ^{(0)}+\mathcal {J}_\parallel ^{(1)}$. If each of these terms is evaluated to the right order, we shall retrieve a consistent expression for the adiabatic invariant $\mathcal {J}_\parallel$.

2.2. Leading order expression

Let us start by investigating the leading order term of the expansion in $\epsilon$. Setting the expansion parameter $\epsilon$ to zero results in the following integral:

(2.15)\begin{equation} \mathcal{J}_\parallel^{(0)} = \sqrt{2H}\frac{B_\alpha}{\bar{\iota}B_0}\int\frac{\sqrt{1-\lambda(1 +\delta B_b - r \eta \cos \chi )}}{1 - r \eta \cos \chi}\,\mathrm{d}\chi. \end{equation}

This integral is highly reminiscent of that occurring in the magnetic field of a large-aspect-ratio tokamak with circular cross-section, which has been considered by many authors before in various asymptotic regimes (Connor et al. Reference Connor, Hastie and Martin1983; Roach et al. Reference Roach, Connor and Janjua1995; Helander & Sigmar Reference Helander and Sigmar2005; Hegna Reference Hegna2015) and, as such, the derivation closely mirrors these calculations. One may refer to Appendix A for a complete derivation. The main step required is to re-express the integral in terms of a trapping parameter $k$, which we define as

(2.16)\begin{equation} k^2 = \sin^2 \left( \frac{\chi_b}{2} \right), \end{equation}

where $\chi =0$ is defined to be the magnetic well minimum. The most deeply trapped particles reside here, and thus have $k^2 = 0$. The most shallowly trapped particles reside at the tops of the well, namely $\chi _b={\rm \pi}$, and thus correspond to $k^2 = 1$. The two trapped particle classes are connected monotonically, in the sense that $\lambda$ and $k$ maintain the order of trapped classes. With this definition, the integral may be expressed in terms of complete elliptic integrals of the first and second kind (also known as Legendre's integrals, see e.g. Olver et al. Reference Olver, Olde Daalhuis, Lozier, Schneider, Boisvert, Clark, Mille, Saunders, Cohl and McClain2020, § 19), which we define as

(2.17a)\begin{gather} K(k) \stackrel{\cdot}{=} \int_0^{{\rm \pi}/2} \frac{ \mathrm{d}\zeta}{\sqrt{1 - k^2 \sin^2 \zeta} }, \end{gather}

(2.17b)\begin{gather}E(k) \stackrel{\cdot}{=} \int_0^{{\rm \pi}/2} \sqrt{1 - k^2 \sin^2 \zeta} \, \mathrm{d}\zeta . \end{gather}

With these definitions, one can express $\mathcal {J}_\parallel ^{(0)}$ in closed form expanding $\mathcal {I}(0)$ around the smallness of $r$. The result is, to order $O(r^{5/2})$,

(2.18)\begin{equation} \mathcal{J}_\parallel^{(0)} = 4\sqrt{H r\eta }\frac{B_{\alpha0}}{\bar{\iota}_0B_0} \left[ I_1(k) +r\eta\left(I_1(k)\left(\frac{1}{2}-k^2\right) + I_2(k)\right) + O(r^2)\right], \end{equation}

where the functions $I_1$ and $I_2$ describe the behaviour of different trapped-particle classes via $k$,

(2.19a)\begin{gather} I_1(k) = 2\left[ (k^2 - 1) K(k) + E(k) \right], \end{gather}

(2.19b)\begin{gather}I_2(k) = \frac{2}{3} \left[(2 k^2-1) E(k)-(k^2-1) K(k)\right]. \end{gather}

One can readily interpret the leading form of $\mathcal {J}_\parallel ^{(0)}$ by referring back to its basic form in terms of a parallel velocity and bounce distance, $\mathcal {J}_\parallel \sim v_\parallel \ell$. The scaling with $\sqrt {H r\eta }$ is directly related to the reduced parallel velocity that trapped particles must have for them to be trapped.Footnote ² The factor $B_{\alpha 0}/\bar {\iota }_0B_0$ represents the changes in the ‘connection length’ along the magnetic field line. In terms of geometric quantities and using Ampére's law, one can estimate this length scale to be $B_{\alpha 0}/\bar {\iota }_0B_0\sim R_0/(\iota _0-N)$, where $R_0$ is the major radius of the device and $N$ represents the helicity of the $|\boldsymbol {B}|$ symmetry determined by the shape of the axis (Landreman & Sengupta Reference Landreman and Sengupta2019; Rodriguez, Sengupta & Bhattacharjee Reference Rodriguez, Sengupta and Bhattacharjee2022b). As the major radius increases, the total distance travelled by a trapped particle grows and so does $\mathcal {J}_\parallel$. Similarly, decreasing $\bar {\iota }$ increases the distance between bounce points, as field lines become further misaligned with respect to $|\boldsymbol {B}|$ contours. Finally, we observe that $I_1$, which describes the trapped-particle class dependence of $\mathcal {J}_\parallel$ to leading order, monotonically increases with $k$, as does the bounce distance. Under such a perspective, it is then clear that it must vanish for the deeply trapped particles (i.e. $I_1(k=0)=0$).

To provide a full description of $\mathcal {J}_\parallel$ to next order, it is important to note that although the expression we found is correct to order $O(r^{5/2})$, it does not correspond to the value of $\mathcal {J}_\parallel$ to that order. The expression is incomplete, as it is missing the contribution from $\mathcal {J}_\parallel ^{(1)}$.

2.3. The first-order correction

To evaluate the second-order term, we follow (2.14), for which we need the following integral:

(2.20)\begin{equation} \mathcal{J}_\parallel^{(1)}={-}\sqrt{2H}\frac{B_\alpha}{\bar{\iota}B_0}\int_{-\chi_b}^{\chi_b} \left( \frac{\lambda}{2} \frac{\delta B-\delta B_b}{\sqrt{1-\lambda(\bar{B}+\delta B_b)}} + \frac{\delta B \sqrt{1 - \lambda (\bar{B} + \delta B_b)}}{\bar{B}^2} \right)\,\mathrm{d}\chi. \end{equation}

The second term in the integrand is a factor $r$ smaller than the first one (which may be verified by writing expressions explicitly in terms of $k$) and hence, for our current expansion, only the first term needs to be taken into account. This term requires rewriting to involve the integration variable $\chi$ explicitly. Using the near-axis form of $|\boldsymbol {B}|$ for a quasi- and stellarator symmetric field, (2.10),

(2.21a)\begin{equation} \delta B-\delta B_b= r^2 B_{22}^C \left( \cos 2\chi - \cos 2 \chi_b \right). \end{equation}

From this point, the procedure is analogous to the steps followed in the leading order case; the details are given, once again, in Appendix A. We simply denote the central result here, which is the expression for $\mathcal {J}_\parallel ^{(1)}$ expanded to leading order in $r$,

(2.22)\begin{equation} \mathcal{J}_\parallel^{(1)} \approx{-} 4 \sqrt{H r \eta} \frac{B_\alpha}{\bar{\iota}B_0 } \frac{r B_{22}^C}{\eta} \mathcal{I}_{22}^C(k), \end{equation}

where the function $\mathcal {I}_{22}^C(k)$ is defined to be

(2.23)\begin{equation} \mathcal{I}_{22}^C(k)\stackrel{\cdot}{=} I_2(k)-(2k^2-1)I_1(k). \end{equation}

Combining this result with (2.18), we may complete the asymptotic form of the second adiabatic invariant to order $O(r^{5/2})$,

(2.24)\begin{equation} \mathcal{J}_\parallel\approx4\sqrt{H r\eta }\frac{B_{\alpha 0}}{\bar{\iota}_0 B_0}\left\{I_1(k) +r\eta\left[\left(\frac{1}{2}-k^2\right)I_1(k)+I_2(k)-\frac{B_{22}^C}{\eta^2}\mathcal{I}_{22}^C(k)\right]\right\}. \end{equation}

3.1. Leading order precession: a tokamak-like behaviour

Let us start by analysing the leading order precession of trapped particles focusing on $\omega _{\alpha,0}$, (3.2). The expression includes physics in two ways: the overall scaling factors in front and the $k$ dependence through $G(k)$, which describes the behaviour of different classes of trapped particles. We first investigate the former.

Precession is proportional to $\eta$, the parameter defined in (2.9) as a measure of the leading order variation of $|\boldsymbol {B}|$ over flux surfaces. This variation is however intimately linked to the near-axis elliptical shaping of cross-sections (see Garren & Boozer Reference Garren and Boozer1991a; Landreman & Sengupta Reference Landreman and Sengupta2019; Rodríguez Reference Rodríguez2023). In fact, for up-down symmetric cross-sections, the elongation along the horizontal axis (i.e. ratio of horizontal to vertical axes of the cross-section) is $\mathcal {E}=(\eta /\kappa )^2$, where $\kappa$ is the curvature of the magnetic axis at the point where the cross-section is being assessed. Thus, for a fixed elongation, $\eta \sim \kappa$. In the special case of an axisymmetric field, this means that $\omega _\alpha \sim \sqrt {\mathcal {E}}/R_0$, where $R_0$ is the major radius. Going back to $\omega _{\alpha,0}$, for a fixed cross-sectional shape, increasing the major radius reduces precession, a consequence of the field becoming more straight and the gradients in $|\boldsymbol {B}|$ (and with it the drift) becoming smaller. In quasisymmetric fields, the local curvature of the axis defines a ‘major radius’, which leads to strongly curved magnetic axis shapes having increased precession. Quasi-helically symmetric fields require more strongly shaped magnetic axes (for a fixed average major radius) and thus will tend to have a stronger precession. This provides a qualitative separation between quasi-axisymmetric (QA) and quasi-helically symmetric (QH) stellarators (Rodriguez et al. (Reference Rodriguez, Sengupta and Bhattacharjee2022b), see some values of $\eta$ in table 1).Footnote ⁴

Table 1. Some $\eta$ values for QS-optimised configurations. This table presents the values of $\eta$ for many quasisymmetric designs, normalised to have an axis whose average major radius is $R_{00}=1$. This is a form of comparing them on the same basis. The largest values of $\eta$ correspond to quasi-helically symmetric configurations (namely HSX, QHS48 and precise QH), although there exist significant variations within these classes.

We finally take note of the divergent nature of $\omega _\alpha$ with the radial coordinate. The $1/r$ behaviour may initially appear worrying, but it can be easily understood in the following terms. From the form of the poloidal drift, we estimate $\boldsymbol {v}_D \boldsymbol {\cdot } \boldsymbol {\nabla } \theta \sim v_{\boldsymbol {\nabla } B} |\boldsymbol {\nabla } \theta |$, where $v_{\boldsymbol {\nabla } B}$ denotes the gradient drift driven by the radial variation of $B$. As $\boldsymbol {\nabla } \theta \sim 1/r$ and the gradient drift does not scale with $r$ to leading order, the result is the $1/ r$ dependence (the result of a finite drift velocity over an ever shrinking surface).

We next shift our attention to the dependency on the trapped class dependence of (3.2). A plot of $G$ as a function of $k$ is presented in figure 1(b). The plot shows that for the vast majority of trapped particles, the drift of the electrons is positive. Physically, positive values imply that the trapped particles precess in the direction of the diamagnetic drift (see figure 1a), an important feature which will become relevant for the discussion on trapped particle modes later. This behaviour only changes for the barely trapped particle classes which end up spending a significant fraction of their bounce-time near the maximum of $|\boldsymbol {B}|$, where there is ‘good curvature’. The transition point occurs at $k_0$, where $G(k_0)=0$, roughly at $k_0\sim 0.9$. Such a class of particles is, to leading order, stationary. The existence of these groups of trapped particles precessing in opposite directions proves the impossibility of making quasisymmetric configurations exactly maximum-$\mathcal {J}$. This simply follows from the definition of maximum-$\mathcal {J}$ as the property of a field that guarantees $\partial _\psi \mathcal {J}_\parallel <0$ for all trapped particles, which in terms of the electron precession is equivalent to $\omega _\alpha (k)<0$ for all $k$. Of course, this is not to say that the behaviour of a quasisymmetric field cannot become closer to maximum-$\mathcal {J}$, as higher order shaping and equilibrium parameters modify the leading order behaviour above. However, one may not achieve it exactly everywhere and especially close to the axis. In practice, one may only get around this issue at a finite radius if the higher order contributions are strong enough.

Figure 1. Precession of trapped particles near the magnetic axis. (a) Diagram of the main drifts in a magnetic configuration with a dominant $\boldsymbol {\nabla } B$ direction, such as in a tokamak. The electron particle $\boldsymbol {\nabla } B$ drift, $v_{\boldsymbol {\nabla }}$, and the diamagnetic drifts, $v_*$, are shown (the latter proportional to $-\boldsymbol {B}\times \boldsymbol {\nabla } p$). Resonance between the two occurs on the outboard side (where the deeply trapped particles live), defining the bad curvature region. (b) Plot showing the function $G(k)$ as a function of the trapped-particle class $k$, and thus the behaviour of precession as a function of trapped particle class near the magnetic-axis.

Comparing this result against previous investigation, we find, as we already pointed out, the behaviour to be identical to that shown in the work by Connor et al. (Reference Connor, Hastie and Martin1983) for a large-aspect-ratio circular tokamak. That this, correspondence found in the more general quasisymmetric case should not come as a surprise, given the existing isomorphism between axisymmetric and quasisymmetric fields (Boozer Reference Boozer1983). We have, however, gained generality beyond circular-shaped cross-sections as $\eta$ allows for non-unity ellipticity. We also note that the asymptotic considerations here and in the literature are in many respects different. Many previous investigations have focused on employing radially local solutions to the MHD equation (see e.g. Mercier & Luc Reference Mercier and Luc1974; Miller et al. Reference Miller, Chu, Greene, Lin-Liu and Waltz1998; Hegna Reference Hegna2000) to discuss precession, which weakens the coupling between $|\boldsymbol {B}|$ and the geometry of the field that exists in the near-axis treatment. We take that additional coupling in the near-axis consideration to form part of a more globally consistent description of the field. This results in higher order effects showing quantitative differences (though qualitative trends are retained), as we shall now explore.

3.2. The higher order effects on precession

Let us now focus on the effect of second-order elements on precession. In the form presented in (3.3), the order $r$ correction on the leading order precession consists of three different terms: an ‘intrinsic’ term (purely a consequence of consistency of the field with its elliptical shape), a term that is proportional to $B_{20}$ (and thus the radial increase of the average $B$) and another proportional to $B_{22}^C$. The behaviour of each one of these terms with $k$ is illustrated in figure 2(a).

Figure 2. Effect of second-order quantities on precession. (a) Plot showing the dependence on the trapped particle class $k$ of the three components of $\omega _{\alpha,1}$. (b) Dependence of precession on triangularity of an axisymmetric tokamak as a function of $k$ for a number of values of $\alpha =\mathcal {E}^2$, where $\mathcal {E}$ is the elongation in the major radius direction. (c) Dependence of precession on pressure gradient in an axisymmetric tokamak as a function of $k$ for a number of values of $f=B_0^2(1+\alpha )^2/R_0^2I_2^2(\alpha +3)$.

From these three contributions, that coming from $B_{20}$ (often called the magnetic well term) is the simplest: a positive $B_{20}$ pushes particles against the diamagnetic drift. That is, deeply trapped particles decrease their precession rate, while barely trapped ones precess even faster. This behaviour is a direct consequence of the influence of $B_{20}$ on the gradient $\boldsymbol {\nabla } B$. The magnetic well term reinforces the outwards magnetic field gradient everywhere, affecting all particles equally and in the direction opposed to the diamagnetic drift. More precisely, the drift $v_{\boldsymbol {\nabla } B}\sim \boldsymbol {B}\times \boldsymbol {\nabla }(1/B)$ is driven by the gradient of $1/B$ and, thus, it is the change in the gradient of $1/B$ that most directly affects precession. As $\partial _\psi (1/B) \sim \eta ^2/2-B_{20}$, this explains not only the $B_{20}$ contribution, but also the ‘intrinsic’ $\mathcal {G}$ one.

The direct effect of $B_{20}$ relates precession naturally to MHD stability. MHD stability of interchange modes is improved by the enhancement of the so-called magnetic well of the field (Greene Reference Greene1997), which corresponds in the near-axis limit to increasing $B_{20}$ (the radial derivative of the average $B$) (Landreman & Jorge Reference Landreman and Jorge2020; Kim, Jorge & Dorland Reference Kim, Jorge and Dorland2021; Rodríguez Reference Rodríguez2023). Thus, there is a natural synergy between improving MHD stability and making particles precess in the direction opposite to the diamagnetic drift. This behaviour, obvious from the $B_{20}$ dependence of $\omega _{\alpha,1}$, aligns with the general observation made in Helander (Reference Helander2014, § 3.7) relating the ‘averaged’ behaviour of precession over a flux surface and all particle classes to the magnetic well. However, the problem of MHD stability is more subtle, as precession is also affected by the variation of the magnetic field $B_{22}^C$ explicitly, and MHD stability is further influenced by pressure (as becomes clear when considering the Mercier criterion (Mercier Reference Mercier1962, Reference Mercier1974; Greene & Johnson Reference Greene and Johnson1962; Bauer, Betancourt & Garabedian Reference Bauer, Betancourt and Garabedian2012; Freidberg Reference Freidberg2014) to assess it). We shall revisit this relation on more solid grounds later.

Setting $B_{20}$ aside for now, the $B_{22}^C$ contribution to precession introduces a richer $k$ dependence than considered so far at this order. This is so because different trapped particles experience different modifications of the magnetic field through the $\chi$ dependence of $|\boldsymbol {B}|\sim B_{22}^C\cos 2\chi$. Naturally, both the most deeply and shallowly trapped particles, who live at $\chi =0,{\rm \pi}$, respectively, feel the same effect, marked by $\mathcal {G}_{22}^C(0)=\mathcal {G}_{22}^C(1)$. In between these classes, particles perform an average of the gradient over their bounce, leading to a maximum at $k \sim 0.8$.

The components $B_{22}^C$ and $B_{20}$ have proven especially convenient to describe the higher order effects on drifts. However, they do not provide a good sense for what the field looks like in terms of its geometry or its equilibrium. A physical discussion requires us to make this connection, which we shall do within the near-axis framework. To that end, we introduce the pressure gradient supported by the field as $p_2=(B_0\,\mathrm {d}p/\mathrm {d}\psi |_{\psi =0})/2$, as well as the triangularity of cross-sections (normalised by $r$)Footnote ⁵ as $\delta$. These two parameters can substitute $B_{22}^C$ and $B_{20}$ to write the precession explicitly in terms of $p_2$ and $\delta$.

The details of this linear mapping between parameters were presented by Rodríguez (Reference Rodríguez2023). We include in this paper only the key elements necessary to make that connection. This is particularly important to make sense of what $\delta$ truly represents. For an up-down symmetric cross-section, we define triangularity as the relative displacement of the vertical turning points of a cross-section with respect to its centre-point along the symmetry line normalised to its width (Rodríguez Reference Rodríguez2023, Appendix C), and $\delta$, as its asymptotic form in $r$, normalised by $r$. For simplicity, we define this triangularity in the near-axis, Frenet–Serret frame. This makes the regular notion of triangularity in the lab frame (i.e. the triangularity of the cross-section at a constant cylindrical angle) generally different by an offset when the magnetic axis is not perpendicular to a constant cylindrical angle plane (see some details in Appendix D). However, by changing $\delta$, we are changing triangularity in the lab frame by the same amount, although $\delta$ is generally not the triangularity there. The axisymmetric case is an exception in which this correspondence is exact. We also pick the sign of triangularity to be defined with respect to the direction of curvature of the axis; i.e. positive triangularity, $\delta >0$, indicates a shift of the turning points in the direction of the curvature.Footnote ⁶

In the general quasisymmetric scenario, following this definition of $\delta$, there appears to be a multitude of ‘triangularities’. Each cross-section around the torus is generally different (but for an $N$-fold symmetry), but it is sufficient to describe the triangularity of any of its cross-sections to describe the field uniquely (given some choice of lower order parameters and a pressure gradient).Footnote ⁷ Once such a cross-section has been chosen, then one should interpret $\delta$ as a measure of its triangularity in the Frenet–Serret frame and any discussion about the effect of modifying triangularity should be interpreted as the effect of changing the triangularity of that very cross-section. We shall conveniently choose an up-down symmetric cross-section to represent the quasisymmetric stellarator and when it comes to analysing real configurations, we shall take the most vertically elongated one (often referred to as the bean-shaped cross-section Rodríguez Reference Rodríguez2023). We do so by analogy with the axisymmetric scenario. As this cross-section is changed, the remainder of the field must also change as a consequence of satisfying both equilibrium and quasisymmetry simultaneously.

In short, the pressure gradient and the triangularity as defined above provide sufficient information and a minimal second-order parametrisation for both axisymmetric and quasisymmetric configurations.

3.3. Relation to geometric and equilibrium parameters

Let us then proceed to write and analyse the precession of trapped particles $\omega _{\alpha,1}$ in terms of triangularity, $\delta$, and pressure gradient, $p_2$. The details of how to do so are presented in Appendix D and rely heavily on the work of Rodríguez (Reference Rodríguez2023). The result reads

(3.5)\begin{equation} \omega_{\alpha,1}=\frac{H \eta}{B_0 r}\left[ r\eta \tilde{\mathcal{G}}(k) - \frac{r\mu_0p_2}{|\eta| B_0^2}\mathcal{G}_{p_2}(k) + \frac{r \delta}{2}\mathcal{G}_{\delta}(k)\right]. \end{equation}

The function $\mathcal {G}_{p_2}$ encodes the effect of the pressure gradient and $\mathcal {G}_\delta$ that of the triangularity, and their full explicit forms may be found in Appendix D. The function $\tilde {\mathcal {G}}$ is a complicated function of lower order quantities that we do not present explicitly and shall ignore for the analysis in this paper. For a discussion on the form and meaning of the other contributions, we specialise now to a generally shaped, up-down symmetric tokamak configuration.

In the axisymmetric limit, the functions become

(3.6a)\begin{gather} \mathcal{G}_{p_2}^\mathrm{axi}={-}2\left[1+ \left(\frac{B_0}{R_0 I_2}\right)^2\frac{(1+\bar{\alpha})^2}{\bar{\alpha}+3}\left(1 +\frac{\mathcal{G}_{22}^C(k)}{4}\right)\right], \end{gather}

(3.6b)\begin{gather} \mathcal{G}_{\delta}^\mathrm{axi}={-}6 \left(\frac{1-\bar{\alpha}}{3+\bar{\alpha}} \right)-\frac{3-\bar{\alpha}}{3+\bar{\alpha}}\mathcal{G}_{22}^C(k), \end{gather}

using the definitions in (3.4). Here the parameter $\bar {\alpha }=(\eta R_0)^4=\mathcal {E}^2$ is the square of the elongation of the flux surfaces along the major radial direction. To arrive at such an expression, we used the relation $\bar {\iota }_0=2\sqrt {\bar {\alpha }}G_0I_2/B_0^2(1+\bar {\alpha })$, which holds true for a tokamak, whose rotational transform must be fully driven by a toroidal plasma current.

Let us analyse the behaviour of the finite pressure term first. It is clear from (3.6a) that $\mathcal {G}_{p_2}^\mathrm {axi}<-2$ for all $k$ and possible combinations of parameters, as $\mathcal {G}_{22}^C\geq -2$ and, therefore, $1+\mathcal {G}_{22}^C/4\geq 1/2$. This negative sign of $\mathcal {G}_{p_2}^\mathrm {axi}$ indicates that the usual peaked pressure profile (i.e. $p_2<0$ for the assumed $\mathrm {sgn}(\psi )=+1$) leads to an increase in precession in the direction opposite to the diamagnetic frequency; a direct consequence of the magnetic well term discussed in the previous section. A finite $\beta$ (at fixed triangularity) assists in making the behaviour of trapped particles more maximum-$\mathcal {J}$. This is a well-known effect (Rosenbluth & Sloan Reference Rosenbluth and Sloan1971; Connor et al. Reference Connor, Hastie and Martin1983), referred to as the diamagnetic effect of the toroidal field. In fact, we note that in the circular tokamak limit, the expression becomes almost identical to the result of Connor et al. (Reference Connor, Hastie and Martin1983), see the details in Appendix D. Although maximum-$\mathcal {J}$ is being approached by increasing $\beta$ and this is generally regarded as a positive effect, at an intermediate stage, particle precession reaches $\omega _\alpha \sim 0$ for many trapped particles. This slow precession scenario is generally associated with enhanced fast particle losses (especially of deeply trapped particles) whenever deviations from QS exist (Nemov et al. Reference Nemov, Kasilov, Kernbichler and Leitold2008; Bader et al. Reference Bader, Anderson, Drevlak, Faber, Hegna, Henneberg, Landreman, Schmitt, Suzuki and Ware2021; Velasco et al. Reference Velasco, Calvo, Mulas, Sánchez, Parra and Cappa2021).

Different trapped particles are affected differently by pressure as a consequence of the evolving Shafranov shift (Shafranov Reference Shafranov1963; Wesson Reference Wesson2011; Rodríguez Reference Rodríguez2023), which perturbs the field in a non-symmetric way. This underlying role of the Shafranov shift is clear from the contribution of the factor $f=B_0^2(1+\bar {\alpha })^2/R_0^2I_2^2(\bar {\alpha }+3)$, which shows an amplified effect of pressure for low currents (i.e. or low rotational transform), larger horizontal elongation or smaller major radius. All of these are known to increase the Shafranov-shift effect and will enhance the counter-precession of particles with respect to the diamagnetic drift (see figure 2).

Because, to leading order, deeply trapped particles co-rotate with the diamagnetic drift, we may estimate when the plasma $\beta$ is sufficient to reverse their direction. We interpret the resulting as the critical plasma $\beta$ for which the field becomes barely maximum-$\mathcal {J}$. Formally, this involves equating two different asymptotic orders, which goes against the very nature of the asymptotic treatment. One may nevertheless interpret this as an estimate of the precession at a ‘finite’ radius.Footnote ⁸ This shows that one may try to increase the maximum-$\mathcal {J}$ behaviour of a QS by enpowering some of the higher order contributions (pressure and other shaping). In practice, the effective radius in which the leading order is dominant may be small enough that we may refer to the field as maximum-$\mathcal {J}$. As shown in the examples of figure 4, this does not seem to be the natural tendency of QS fields and certainly is not asymptotically.

Focusing on the behaviour of $k=0$ (such deeply trapped particles are typically the least maximum-$\mathcal {J}$), $\omega _{\bar {\alpha },0}=H\eta /rB_0$ from (3.2) and for the pressure, we have $\omega _{\alpha,1}=-(H/B_0)\mu _0|p_2|(2+f)/B_0^2$. Equating the two, we find

(3.7)\begin{equation} \frac{\mu_0 |p_2|}{B_0^2} \sim \frac{\eta}{r(2+f)} \sim \frac{\sqrt{\mathcal{E}}}{R_0 r(2+f)}. \end{equation}

For a parabolic pressure profile, $a^2 p_2 = -p_0$, where $a$ is the minor radius and $p_0$ the pressure on axis, one finds that the critical plasma $\beta$ on axis is

(3.8)\begin{equation} \beta_\mathrm{crit} \stackrel{\cdot}{=} 2a^2\frac{\mu_0 p_2}{B_0^2} \sim \frac{a^2}{R_0 r} \frac{2\sqrt{\mathcal{E}}}{2 + f}. \end{equation}

We thus see that the most susceptible fields are those with a large-aspect ratio, vertical elongation (small $\mathcal {E}$) and lower current (large $f$). As expected, these finite $\beta$ effects become more pronounced as we move away from the magnetic axis. For further illustration, consider the scenario of a circular-shaped tokamak with a representative safety factor of $q=2$ and aspect ratio $\sim 3$ for which, at the edge, $\beta _\mathrm {crit}\sim (a/R_0)/(1+q^2/2)\sim 11\,\%$. Reversing the precession of deeply trapped particles is thus predicted to require a significant plasma $\beta$.

The effects of triangularity are markedly more involved than those of pressure, which prevents us from writing a parameter-insensitive bound like we did for the effect of pressure (see figure 2). Depending on the value of elongation, $\bar {\alpha }=\mathcal {E}^2$, the precession due to triangularity will tend to be in one direction or the other, a behaviour that also changes depending on the particle class considered. There exists, though, a critical value of $\bar {\alpha }$ beyond which $\mathcal {G}_\delta ^{\mathrm {axi}}>0$ for all $k$. As $\mathcal {G}_{22}^C$ has a maximal value $\max (\mathcal {G}_{22}^C) \approx 1.1$, one can find that this critical point occurs at

(3.9)\begin{equation} \bar{\alpha}_\mathrm{crit} > \frac{6 + 3 \max(\mathcal{G}_{22}^C)}{6 + \max(\mathcal{G}_{22}^C)} \approx 1.3. \end{equation}

Thus, for tokamaks that are horizontally elongated (beyond some $\sim 14\,\%$), negative triangularity tends to make all particles precess against the diamagnetic drift. This distinction regarding the effect of triangularity is reminiscent of the effects of triangularity on MHD stability. In that case and describing stability through the Mercier criterion Solov'ev & Shafranov Reference Solov'ev and Shafranov1970; Lortz & Nührenberg Reference Lortz and Nührenberg1978; Shafranov Reference Shafranov1983; Freidberg Reference Freidberg2014; Rodríguez Reference Rodríguez2023), one can show that for $\bar {\alpha }>\bar {\alpha }_\mathrm {MHD}=1$, negative triangularity contributes positively to stability. Thus, MHD stability seems to align with precession against the diamagnetic drift, at least for sufficiently horizontally elongated configurations. That is, there is some synergy, which is the opposite to the changes due to plasma-$\beta$.

Most commonly, however, most tokamak fields are vertically elongated and thus have $\bar {\alpha }<1<\bar {\alpha }_\mathrm {crit}$. In that usual scenario, different trapped particles respond differently, some tending to precess in one direction and others in the opposite. The most deeply and barely trapped particles are a special case, as taking $k=0,~1$, $\mathcal {G}_\delta = 4 \bar {\alpha }/(3 + \bar {\alpha })>0$ has a maximum value (see figure 2). With a sign independent of elongation, one can conclude that positive-triangularity tokamaks will always tend to make deeply and shallowly trapped particles precess in the direction of the diamagnetic drift. Thus, only through negative triangularity can this shaping be used to reverse the behaviour of deeply trapped particles. In the vertically elongated regime, negative triangularity hampers MHD stability, thus opposing the tendency to improve the maximum-$\mathcal {J}$ behaviour. As noted in the plasma $\beta$ scenario, as precession of particles is reduced, fast-ion confinement can be negatively affected in an imperfect QS stellarator. The different behaviour of each trapped particle makes an assessment of the overall effect complex.

It would be appropriate, as we did to illustrate the effect of plasma $\beta$, to introduce the idea of a critical triangularity: a value of triangularity such that one expects precession of deeply trapped particles to leading order to be significantly affected, i.e. $r \delta \mathcal {G}_\delta (k)/2 \sim G(k)$ for $k=0$. Precisely as in the case of pressure,

(3.10)\begin{equation} (r\delta)_\mathrm{crit}\sim\frac{3+\bar{\alpha}}{2\bar{\alpha}}. \end{equation}

The interpretation of $r\delta$ as triangularity of the cross-section in the Frenet frame of the magnetic axis (see Appendix D) is an asymptotic concept. As such, this geometric interpretation of $r\delta$ ceases to be accurate upon approaching unity (especially as a significant bean-shape is developed). This limit of large $r\delta$ is nevertheless a limit of strongly shaped flux surfaces, which could even describe surfaces that self-intersect or intersect with other flux surfaces (Landreman Reference Landreman2021; Rodríguez Reference Rodríguez2023). To make sense of whether $(r\delta )_\mathrm {crit}$ is feasible in practice, we should compare this measure to the maximum triangularity achievable without incurring into these geometric defects. The critical $r_c$ was defined by Landreman (Reference Landreman2021). In any case, (3.10) indicates that a very significant triangularity (of order 1) is needed to significantly affect precession of deeply trapped particles in a tokamak. Hence, we conclude that although triangular shaping may help in achieving the maximum-$\mathcal {J}$ property together with finite $\beta$ effects, achieving it via shaping alone is more difficult in tokamaks.

Before concluding this section, let us briefly consider the full quasisymmetric case, beyond the special case of axisymmetry, through some examples. Unlike in the axisymmetric case, this general scenario preserves a measure of symmetry-breaking of the geometry. The measure $\bar {F}$ (Rodríguez Reference Rodríguez2023), for which explicit expressions are presented in Appendix D, depends on the shape of the axis and modifies the effects of pressure and triangularity. Reminding ourselves that in such a scenario, $\delta$ represents the triangularity of the up-down symmetric cross-section with the largest binormal elongation, we show in figure 3 the values of $\mathcal {G}_\delta$ and $\mathcal {G}_{p_2}$ for a number of QS configurations.Footnote ⁹ Each segment consists of pairs $(\mathcal {G}_\delta,\mathcal {G}_{p_2})$ for different $k$ for the same configuration, taking the ideal QS limit of the configurations (which are only approximately so).

Figure 3. Effect of triangularity and pressure gradient in the precession of trapped electrons in some QS configurations. The plot shows the values of $\mathcal {G}_{\delta }$ and $\mathcal {G}_{p_2}$ for a number of quasisymmetric configurations in the ideal quasisymmetric limit represented by their most vertically elongated up-down symmetric cross-section (in some configurations, this occurs at $\varphi ={\rm \pi} /N$ rather than $\varphi =0$). The scatter plot corresponds to the values of both $k=0,1$ for different configurations, while each segment represents the other $k$ values. The near-axis models can be found in the acknowledged repositories; further details, such as their QS quality, can be found in Appendix E.

From the plot, it is clear that for those quasiymmetric configurations analysed, the effect of a finite pressure gradient is the same as in the axisymmetric limit (i.e. an increase in pasma $\beta$ leads to precession in the direction opposite to the diamagnetic drift). The effect of triangularity is analogous to a tokamak that is elongated in the horizontal direction, where negative triangularity pushes particles against the diamagnetic drift. From the MHD stability analysis of these configurations by Rodríguez (Reference Rodríguez2023), one recovers the synergy of horizontally elongated tokamaks for quasisymmetric stellarators. Thus, the behaviour is quite different from that of regularly shaped tokamaks.

An interpretation of the magnitudes of $\mathcal {G}_{\delta }$ and $\mathcal {G}_{p_2}$ may be provided by considering the critical $\beta$ and $r\delta$ once more. As in the tokamak scenario, at the edge of the configuration, $\beta _\mathrm {crit}\sim 2a|\eta |/\mathcal {G}_{p_2}(0)$ and $(r\delta )_\mathrm {crit}\sim 2/\mathcal {G}_\delta$. A complete table for the configurations in figure 3 is included in Appendix E. We note here that in QA configurations, $\beta _\mathrm {crit}\sim 5\,\%$, while QH stellarators generally exhibit a more resilient behaviour with $\beta _\mathrm {crit}\sim 10\unicode{x2013}15\,\%$. Reversing the precession of deeply trapped particles via finite $\beta$ effects thus appears to be a possibility most easily in QA configurations. Given the simple form of $(r\delta )_\mathrm {crit}$, it is straightforward to see that $O(1)$ triangularity is generally required to observe a significant effect. In many of these configurations, such values are indeed achievable without incurring in forbidden flux surface shapes (see Appendix E).

3.4. Verification of the expansion

In the preceding sections, we investigated the analytical behaviour of the trapped particle precession. We derived these under two important assumptions: (i) exact quasisymmetry (or symmetry) of the fields and (ii) the near-axis expansion. It is thus natural to wonder how close trapped particle precession is to the idealised limit in realistic configurations, as these assumptions become increasingly less accurate. We check this through three numerical examples, in which we compare the bounce-averaged drift computed from a global MHD code (in this case, VMEC (Hirshman & Whitson Reference Hirshman and Whitson1983), using numerical methods detailed by Mackenbach et al. Reference Mackenbach, Duff, Gerard, Proll, Helander and Hegna2023a) against the analytical results. For this practical comparison, we employ the definition of $k$ in terms of the bounce point, (2.16), which we compute for the numerical precession calculation. Note that upon significant deviations from QS, this choice ceases to be convenient, especially when there exist differently shaped wells within a flux surface.Footnote ¹⁰

This comparison is shown in figure 4, where the bounce-averaged precession is plotted as a function of $k$ and for two radial locations, $\varrho \stackrel {\cdot }{=} \sqrt {\psi /\psi _\mathrm {edge}}$. The correspondence is excellent for all displayed $\varrho$ in the precise quasisymmetric configurations, recently presented by Landreman & Paul (Reference Landreman and Paul2022), which are characterised for having an excellent degree of quasisymmetry. The theory works remarkably well even at larger radii, where one would expect the near-axis expansion to falter, although this near-axis nature of the configurations had been previously noticed (Rodriguez Reference Rodriguez2022) and could be a feature necessary to achieve excellent global quasisymmetry. This numerical comparison evidences that the second-order calculation is also appropriate. For HSX (Anderson et al. Reference Anderson, Almagri, Anderson, Matthews, Talmadge and Shohet1995), we see more significant deviations from the analytical result. This is mainly a consequence of the breaking of QS (for a naive fitting of its near-axis behaviour, $B_{20}$ variations are $\sim 4$), and significant shaping beyond the near-axis behaviour (see $\varrho =1$). Although the trends and magnitude seem correct, there are clear deviations from the idealised limit.

Figure 4. A comparison between the analytical and numerical precession. The precession $\omega _\alpha$ for three QS fields and two radial positions $\varrho =\sqrt {\psi /\psi _\mathrm {edge}}$ is presented, normalised to $H / r \sqrt {2 B_0 \psi _{\mathrm {edge}}}$, where $a$ is the minor radius of the configuration, $B_0$ the $B$ on axis and $2 {\rm \pi}\psi _{\mathrm {edge}}$ the total toroidal flux (and we have also normalised the charge out, which we remind the reader was set to unity). There is good correspondence for the precisely quasisymmetric configurations (Landreman & Paul Reference Landreman and Paul2022) and less so for HSX, which exhibits important deviations from QS. The black and grey lines are the numerical result of precession for different field lines (black corresponding to $\alpha =0$), whereas the dashed red and green lines are the first- and second-order analytic precessions, respectively.

4.1. Leading order form of Æ

Let us begin the asymptotic analysis by considering the first-order expression of the trapped-particle precession $\omega _\alpha \approx \omega _{\alpha,0}(k)$, defined in (3.2). To perform the integral in (4.5), we need a number of ingredients. First of all, we must consider the integral only over trapped-particle classes that co-rotate with $\hat {\omega }_{*,0}^T$ (i.e. the range allowed by the Heaviside). The domain of integration thus runs from the most deeply trapped particles to the transition point of $\omega _{\alpha,0}$, i.e. $k \in [0,k_0]$ (with the definition of $k_0|\omega _{\alpha,0}=0$ from before).

The second ingredient that naturally arises is then the need to express the integration variable $\lambda$ in terms of $k$. The presence of $c_1=\hat {\omega }_{*,0}^T/\hat {\omega }_{\alpha,0}(k)$ as a function of $k$ inside $\mathcal {F}$ makes the integral highly nested and, thus, appears to make it difficult to approach analytically. However, given the form of the precession, $c_1$ is asymptotically small in the sense that $c_1\sim r$. This appears to offer a way to proceed by expanding $\mathcal {F}$ in the small $c_1$ limit. However, that would be wrong, as it would neglect the most important contribution to Æ. Recall that the particle precession vanishes at $k_0$ and, thus, near this value of $k$, the function $c_1$ cannot be considered small. This resonant behaviour is, in the asymptotic limit, the principal contributor. The integral is significant only in a narrow region $\Delta k\sim r$, close to $k_0$, where $c_1 \sim O(1)$ (see a clarifying sketch in figure 6). This teaches us that in this asymptotic limit, the most important class of particles are those with relatively small precession, as in this limit, $\omega _\alpha$ tends to be much larger than the diamagnetic drift. The evaluation of the integral may then be carried out analytically considering a local approximation of the integrand in this contributing narrow band (correct down to a correction $O(r^{1/2})$ on the leading contribution), the details of which are presented in Appendix F.

Figure 6. Schematic of the contribution to available energy. (a) Diagram showing a narrow band of trapped particles near $k_0$ (the trapped-particle class with vanishing precession) contributing to the available energy. The broken line indicates $k_0$ in the limit of $r\rightarrow 0$, as the vertical direction denotes different classes of particles with the leading order magnetic well structure shown by the black curve. The plot of $\mathcal {F}[c_1(k)]$ is shown to the right for $r\omega _{*,0}^T/\eta =0.1$ as an example. (b) A numerical calculation showing the distribution of Æ across the magnetic well, normalisedby the total Æ. Plotted for the precise QA device at a minor radial coordinate of $r=10^{-3}$. The points where $\omega _\alpha =0$ is denoted by a green dashed line and $\hat {\omega }_n = 0.1$.

Before writing the result for the Æ out, one last consideration is required. In this case, one needs to make an explicit assumption regarding the width of the flux tube $\Delta \psi$, on which the Æ will depend. This width should be interpreted as the ‘length’ over which density gradients may be flattened by the turbulence to extract energy. Following the steps taken by Mackenbach et al. (Reference Mackenbach, Proll and Helander2022, Reference Mackenbach, Proll, Wakelkamp and Helander2023c), we estimate such a flattening length scale to be the correlation length and to be of the order of the Larmor radius $\rho$. As such, we write

(4.7)\begin{equation} \Delta \psi = B_0 r \Delta r = B_0 r C_r \rho, \end{equation}

where $C_r$ may formally be dependent on other equilibrium parameters (e.g. the rotational transform $\iota$ or the flux expansion $|\boldsymbol {\nabla }\psi |$). We shall nevertheless take $C_r$ to be a constant, assuming that the flattening length scale providing free energy to the TPMs is simply proportional to the Larmor radius.

It is customary to define a minor radius $a$ for the field configuration so that the radial coordinate may be normalised as $\varrho = r/a$. This way, we define the density gradient $\hat {\omega }_n \equiv - a \partial _r \ln n = - \partial _\varrho \ln n$, which scales like $\varrho$ for a quadratic (in $\varrho$) density profile. In terms of these variables, the Æ becomes

(4.8)\begin{equation} \hat{A}_{w} \approx \frac{2\sqrt{2}}{ 9 {\rm \pi}} C_r^2 \rho_*^2 \left(\frac{\hat{\omega}_n}{\varrho}\right)^3 \frac{\varrho^3\sqrt{\varrho}}{\sqrt{a \eta}}, \end{equation}

where the common gyrokinetic expansion parameter is defined as $\rho _* \stackrel {\cdot }{=} \rho /a$.

The leading order expression for $\hat {A}_{w}$ includes important information regarding the behaviour of the available energy near the axis. We highlight various scalings here.

1. (i) One observes a strong scaling with the distance from the axis $\varrho$, whose origin may be presented in simple terms as follows (see the integral expression in (4.5) for reference). One factor of $\varrho ^{1/2}$ may be accounted for due to the trapping fraction of particles, which leaves one with an overall $\varrho ^3$ scaling. In this scenario, the energy scale is set by the diamagnetic drift (as only precessing particles with speeds analogous to those of the diamagnetic drift contribute to the available energy), which goes like $\hat {\omega }_n\propto \varrho$ near the axis. Thus, two powers of $\varrho$ can be argued to come from the kinetic drive of the diamagnetic drift. The final power of $\varrho$ comes from the small fraction of trapped particles that contribute to the available energy, as the band near $k_0$ scales with $\varrho$.

2. (ii) Another scaling of interest in (4.8) is its dependency on the stellarator shaping parameter $\eta$. Increased $\eta$ leads to a more restricted access to energy and, thus, a reduced TPM activity (as measured by Æ). Thus, horizontally elongated shapes would seem to favour stability and, in the context of quasisymmetric stellarators, quasi-helically symmetric configurations. In tokamak terms, $a\eta \sim a\sqrt {\mathcal {E}}/R_0\sim \varepsilon \sqrt {\mathcal {E}}$, where $\varepsilon$ is the inverse aspect ratio. Thus, $\hat {A}_{w} \sim 1/(\mathcal {E}^{1/4}\sqrt {\varepsilon })$. Hence, the available energy increases with aspect ratio keeping the minor radius fixed. This dependency becomes even stronger if one chooses $\rho _* \varepsilon = \rho /R_0$ to be constant.

3. (iii) One finally observes a scaling with $(C_r \rho _*)^2$, which is the square of the assumed length scale over which energy is available. Importantly, we note that an implicit magnetic-field-strength dependency arises here (for fixed minor radius and $C_r$), as $\rho \sim 1/B_0$. Hence, in terms of Æ, it is beneficial to have a strong magnetic field so that the length scale over which energy is available decreases as $\hat {A}_{w} \sim 1/B_0^2$.

As noted before, a full interpretation of these scalings for a comparison between different configurations would have to take into account the form of $C_r$ that most closely describes the volume that is accessible to the rearrangement of energy. This may be particularly important when proceeding to a comparison between highly different configurations. The normalisation with $C_r$ being a constant appears to provide a reasonable description in Æ as a measure of heat transport in practice (see Mackenbach et al. Reference Mackenbach, Proll and Helander2022).Footnote ¹¹

4.2. A strongly driven finite-radius asymptotic regime

In the derivation above, it was key to consider the contribution to available energy from a narrow band of trapped particles with ‘slow’ precession. As such, the particular form of the expression in (4.8) is only valid in the regime where $\omega _{\star,0}^T/ \omega _\alpha$ can be considered small, that is, when the trapped particle drifts are (as a group) much faster than the diamagnetic drift. As these roles revert, either because the turbulence becomes strongly driven (i.e. large density gradient) or the precession diminishes, the ‘weak’ approach presented before will cease to provide us with a good approximation to the available energy.

As the precession becomes smaller relative to the diamagnetic drift, we however find another tractable limit, which we refer to as the ‘strongly driven’ limit. That is, we still consider a near-axis description of the magnetic field and precession, but at the same time, order the diamagnetic drift to be large, i.e. vigorously driven.Footnote ¹² Although this might appear inconsistent, it is not, as any ordering assumption about $\omega _n$ only affects the evaluation of the Æ integral, but not the particle precession itself. From the set-up, it should be clear that this ‘strongly driven’ regime gains relevance away from the magnetic axis, where the precession of trapped electrons decreases and the diamagnetic drift increases.Footnote ¹³

In this new scenario, the integral for available energy may be recomputed (see Appendix F) using standard methods, as there no longer exists a narrow contributing band (see figure 7). All in all, one finds that the integral reduces to

(4.9)\begin{equation} \hat{A}_{s} \approx 1.1605~\frac{C_r^2\rho_*^2}{{\rm \pi}\sqrt{2}} \frac{\hat{\omega}_n}{\varrho} (a \eta)^{3/2}\varrho^{3/2}. \end{equation}

A different regime brings a different scaling with $\varrho$ and $\hat {\omega }_n$, in both cases with weaker dependencies than in the narrow-band regime. These changes are a result of: (i) the particle precession that serves as energy source contributing directly to Æ, and thus introducing a $1/\varrho \hat {\omega }_n$ factor compared with the weak regime (simply because, on average, particles do not quite reach the diamagnetic drift) and (ii) the contributing trapped particle fraction corresponding to the whole population with a positive precession, which no longer is a narrow band and thus does not contribute with an additional $\hat {\omega }_n\varrho$ factor. Importantly, the scaling with $\eta$ inverts compared with the weak regime, $\hat {A}_{w}$. Using the same tokamak estimation for $\eta$ as there, one finds $\hat {A}_{s} \sim \varepsilon ^{3/2} \mathcal {E}^{3/4}$, suggesting that a large-aspect-ratio tokamak which is vertically elongated reduces Æ. Once again, this is under the assumption of keeping the minor radius $a$ fixed. The behaviour will change depending on which features of the equilibrium are kept constant.

Figure 7. Schematic of the contribution to available energy. (a) Diagram showcasing contribution of $\mathcal {F}$ to the Æ integrand. The broken line satisfies $\omega _\alpha (k_0)=0$ to leading order. The plot of $\mathcal {F}[c_1(k)]$ is shown to the right for $r\omega _{*,0}^T/\eta \gg 1$ as an example. (b) A numerical calculation showing the distribution of Æ across the magnetic well for the precise QA device at a minor radial coordinate of $r=10^{-3}$. The points where $\omega _\alpha =0$ are denoted by a green dashed line and $\hat {\omega }_n = 10^4$.

The existence of these two different regimes of available energy naturally defines a transition. One can calculate this transition point by equating $\hat {A}_{w}\approx \hat {A}_{s}$, denoting the ‘critical’ transition gradient as $-a \partial _r \ln n|_\mathrm {crit} = a / L_{n}|_\mathrm {crit}$. We find

(4.10)\begin{equation} \left. \frac{a}{L_{n}}\right|_\mathrm{crit} \approx 1.61591 a \eta. \end{equation}

For a quadratic density profile ($\hat {\omega }_n/\varrho \approx 2$), the radial position $\varrho$ at which this critical transition occurs is

(4.11)\begin{equation} \varrho_\mathrm{crit} \approx 0.80795 a \eta. \end{equation}

Using some typical tokamak values, we estimate $a \eta \sim \varepsilon \sqrt {\mathcal {E}} \sim 0.3$ and thus $\varrho _\mathrm {crit} \sim 0.2$. Hence, in a standard tokamak, one can transition to the strongly driven regime fairly close to the axis and we expect the strongly driven regime description to be most suitable for most of its volume.

We conclude this leading order Æ analysis by presenting a numerical verification in figure 8, where we compare both asymptotic regimes in one device and the weakly asymptotic regime across multiple devices.Footnote ¹⁴ We observe excellent agreement in the asymptotic behaviour between the found results and the numerical calculation.

Figure 8. A comparison of the numerical calculation of Æ against the analytical result. (a) Comparison of the Æ in the two asymptotic regimes in the precise QA configuration as a function of $\varrho$. The red dashed and dotted lines denote the analytic asymptotic results in the strongly and weakly driven regime, respectively. The critical radial coordinate $\varrho _\mathrm {crit}$ is shown by a blue dash-dotted line. The crosses are numerical evaluations of the Æ using the near-axis equilibrium of the precise QA configuration of Landreman & Jorge (Reference Landreman and Jorge2020). The plot has been generated with a gradient value of $\hat {\omega }_n/\varrho =10^3$ for visualisation purposes. (b) Correlation of the numerical and analytic result in the weakly driven regime for a wide number of quasisymmetric devices (Landreman Reference Landreman2022). The ordinates correspond to the numerically evaluated $\hat {A}$, whereas the abscissa corresponds to the asymptotic result of (4.8), $\hat {A}_{w}$. For this plot, $\hat {\omega }_n/\varrho =1$. A close correspondence can be seen across many orders of magnitude. Both plots have been generated using the pyQSc code, which does not have a notion of minor radius and, as such, they have $\varrho =r$.

4.3. Additional dependence of Æ

To learn anything about how triangularity of flux surfaces and pressure may affect this available energy, and thus how TPMs may be affected by them, we need to proceed to higher order in the calculation of $\hat {A}$. We show how to do this to obtain the dependence of $\hat {A}$ on $p_2$ and $\delta$ to leading order at the end of Appendix F.Footnote ¹⁵ After such considerations, we may write $\hat {A}\approx \hat {A}_0+\hat {A}_1+\ldots$, where $\hat {A}_0$ is the leading order expression and $\hat {A}_1$ the pressure and triangularity dependent piece. As before, for the discussion in the text, we specialise to the generally shaped up-down symmetric tokamak. Full expressions that apply to the quasisymmetric case may be found in Appendix F.

In the weakly driven regime, one finds

(4.12)\begin{equation} \left. \frac{\hat{A}_1}{\hat{A}_0} \right|_\mathrm{weak} \approx \mathcal{R}_{20} \left( - \frac{r}{\eta} \frac{\mu_0 p_2}{B_0^2} \left[1+\frac{4\sqrt{\bar{\alpha}}}{(\bar{\alpha}+3)}\left(\frac{\eta B_{\alpha0}}{B_0\bar{\iota}_0}\right)^2\right] + \frac{3}{2} \frac{1-\bar{\alpha}}{3+\bar{\alpha}}r\delta\right), \end{equation}

where $\mathcal {R}_{20}\approx 1.37$.

It follows from this that, regardless of the choice of parameters, increasing the pressure gradient always leads to an increase in the available energy. Note that this is the case even if the density gradient, here controlled by the diamagnetic frequency $\hat {\omega }_n$, is fixed.Footnote ¹⁶ To picture what is happening, we resort to the discussion on precession presented before. As we increase the pressure gradient, we learnt that the trapped-electron precession goes against the diamagnetic drift, which means that the trapped population is brought further away from resonance. However, from this behaviour, one would expect the drive of the instability to decrease and with it, Æ to do so. So, how is getting further away from the diamagnetic resonance making things worse?

To figure this out, it is illuminating to formally assess the origin of the sign of $\mathcal {R}_{20}$. The sign is, to a large extent, a result of the correction to the integral measure $\mathrm {d}k/\mathrm {d}c_1$ needed when writing the Æ integral in $c_1$ (as it is necessary for the weak regime calculation, see Appendix F). This piece of the integral measures the number of trapped particles that exist with a precession that is similar to the resonant one. The question is then how this population fraction changes as the precession slows down. The answer is that the population that has a near-vanishing precession grows, as most directly seen in the smaller gradient of $\omega _\alpha$ with $k$ (see figure 1b). Because this fraction of the population is the only one that may contribute to the total available energy, the result is the increase of Æ with plasma $\beta$. This is an important feature of available energy, which not only assigns value to the magnitude of $\omega _\alpha (\omega _{*,0}^T-\omega _\alpha )$, but also to the number of particles with a particular value for its precession. As a consequence, we expect this behaviour to change in the strongly driven regime, which we will visit later.

The effect of triangularity, $\delta$, in (4.12) depends critically on whether cross-sections are elongated vertically or horizontally, as we saw MHD stability to do in the preceding discussion. In the most common case of vertically elongated cross-sections, negative triangularity is seen to reduce the Æ (which increases the precession of the trapped-particle class at $k_0$). Thus, the effect is not synergistic with MHD stability, as triangularity has precisely the opposite effect there. This anti-correlation holds also in the more general case of quasisymmetric configurations, which is readily seen by comparing (F26) for $\mathcal {R}_\delta$ directly to Rodríguez (Reference Rodríguez2023, (4.2)) for $\mathcal {T}_\delta$. This intimate relation between MHD stability and what may be interpreted as TPM activity has been observed on many occasions (in fact, could be interpreted as the driver for many reverse triangularity studies in advance tokamak scenarios). Here we have formally proven in the weak asymptotic regime that a compromise between the two properties is needed in this regime. This opposed behaviour is not shared by plasma $\beta$, which acts in a detrimental form on both MHD and TPM activity.

Performing a similar analysis in the strongly driven regime, we find $\hat {A}_1/\hat {A}_0|_{\mathrm {strong}}\approx -2.85\hat {A}_1/\hat {A}_0|_{\mathrm {weak}}$, which presents the opposite sign to the weak regime. That is, an increased plasma $\beta$ (in the form of pressure gradient) always decreases the Æ and in a standard tokamak with $\bar {\alpha }<1$, positive triangularity becomes TEM-stabilising. In the strongly driven regime, reducing precession brings the zero-crossing point closer to $k=0$, thus reducing the total $k$-space available to drive TEMs. In that limit, with both precession and accessible population decreasing with increased pressure gradient and positive triangularity, we expect available energy to decrease and, regarding fast particle confinement due to non-QS behaviour, to momentarily grow before eventually decreasing as precession grows in the direction of maximum-$\mathcal {J}$. The details will depend on how different trapped particles are affected and how important is their contribution to confinement. Unlike in the weak regime, the whole trapped population becomes important and not just a narrow band, figure 6. In the strongly driven regime, there is a synergy between MHD stability and TEM activity with respect to the triangular shaping of cross-sections, but opposed in the effect of plasma $\beta$. The expected behaviour of a stellarator will thus depend critically on the particular regime considered.

In addition to the sign, there is also a difference in magnitude of roughly a factor $3$ between the relative effects of triangularity and plasma $\beta$ in the strong regime compared with the weak regime. For the discussion following, we focus on the strongly driven regime. This effect can be quantified as we did in the discussion of precession, which we do as follows. When the first-order correction significantly affects the available energy, i.e. $\hat {A}_1/\hat {A}_0\sim 1$, we state that we have a critical scenario. At the edge, the critical $\beta$ becomes $\beta _\mathrm {crit}^{\unicode{x00C6}}\sim 2a|\eta |/\mathcal {R}_{20}(1+f)$, with $f$ as defined before (with its QS generalisation, which may be found in (D5a)). This shows that plasma $\beta$ becomes effective in significantly changing Æ in the strong regime for QAs in the regime of $\beta _\mathrm {crit}^{\unicode{x00C6}}\sim 2\unicode{x2013}3\,\%$, while QH $\beta _\mathrm {crit}^{{\unicode{x00C6}} }\sim 4\unicode{x2013}7\,\%$. Finite $\beta$ QA equilibria appear, therefore, to significantly affect the behaviour of TPM-like behaviour, while QHs remain more resilient, as expected from the behaviour of precession. As far as triangularity is concerned, the expression in (3.10) may be used for Æ with $\mathcal {G}_\delta =3[(3+\bar {\alpha })-(\bar {\alpha }+1)\bar {F}]/[(1-\bar {\alpha })+(1+\bar {\alpha })\bar {F}]$ there. The values for QS configurations may be found in Appendix E in the range $(r\delta )_\mathrm {crit}\sim 0.4\unicode{x2013}1.5$, which is a significant triangularity, nevertheless consistently achievable in many configurations (see Appendix E). Note that in the circular tokamak limit, triangularity has no effect on Æ (only a very small one in the weak regime from $\mathcal {R}_{22}^C$).

Given the changes in the weak and the strong regime, the critical transition gradient also changes and may be computed by

(4.13)\begin{gather} \left. \frac{a}{L_{n}}\right|_\mathrm{crit} \approx \left.\frac{a}{L_{n}}\right|_{\mathrm{crit},0}\left[ 1-2.6389 \times \left(- \frac{r}{\eta} \frac{\mu_0 p_2}{B_0^2} \left[1+\frac{4\sqrt{\bar{\alpha}}}{(\bar{\alpha}+3)}\left(\frac{\eta B_{\bar{\alpha}0}}{B_0\bar{\iota}_0}\right)^2\right] + \frac{3}{2} \frac{1-\bar{\alpha}}{3+\bar{\alpha}} r \delta\right)\right]. \end{gather}

This means that the critical gradient decreases for an increased pressure gradient (as the factor multiplying the pressure is always positive) and increases for negative triangularity tokamaks which are vertically elongated.

To close the paper, we make a comparative study of the insight and results presented in this paper with the literature. The comparison to a recent numerical analysis of the Æ for tokamak equilibria described by a Miller (Miller et al. Reference Miller, Chu, Greene, Lin-Liu and Waltz1998) geometry (Mackenbach et al. Reference Mackenbach, Proll, Snoep and Helander2023b) is most straightforward. One can see that many of the trends found there are recovered in the current paper on an analytical basis; in particular, negative triangularity decreases the Æ for vertically elongated tokamaks only if the gradient is sufficiently small and the gradient threshold has the same dependencies as derived here. Of course, our work sheds light on the origin of such behaviour and can be applied beyond axisymmetric configurations.

The comparison to other turbulent study results and experiments requires us to carefully discern between the two distinct regimes where we have shown the Æ exhibits. Depending on which regime a given scenario is in, the beneficial or detrimental nature of various equilibrium shaping parameters may change. In general and bearing this important caveat in mind, the strongly driven regime is most often entered fairly close to the magnetic axis (as argued before). It is also, in magnitude, the principal contributor to Æ and the very character of the weak regime may make it numerically elusive (as the narrow Æ band would have to be resolved in simulations). As such, it is natural to take the Æ in the strongly driven regime as indicative of overall behaviour of a configuration to TPM mediated transport. In terms of leading order effects, the prediction that increasing the vertical elongation in tokamaks improves transport agrees with existing knowledge (Chu, Ott & Manheimer Reference Chu, Ott and Manheimer1978; Qi et al. Reference Qi, Kwon, Hahm, Yi and Choi2019). Concerning higher order effects, the beneficial nature of a pressure gradient on the trapped electron mode has been noted by many authors (Rosenbluth Reference Rosenbluth1968; Connor et al. Reference Connor, Hastie and Martin1983; Li & Kishimoto Reference Li and Kishimoto2002). The effect of triangularity on Æ, however, is more paradoxical. In experiments, it has been shown that negative triangularity exhibits improved confinement whilst remaining in L-mode (Marinoni et al. Reference Marinoni, Austin, Hyatt, Walker, Candy, Chrystal, Lasnier, McKee, Odstrčil and Petty2019). The current model, however, would predict an increase in the Æ in such a scenario and hence more unfavourable transport. Part of this discrepancy may be explained by the role of magnetic shear, which is positive and significant near the edge of the tokamak, but we have not explicitly considered it here. The trapped particle precession increases with increasing magnetic shear, which may push one to a more weakly driven regime in which negative triangularity is beneficial.

However, the discrepancy may also come from the consideration that behaviour within the strongly driven regime may not be the most adequate to describe the turbulent performance of a configuration. To illustrate this, take a scenario in which Æ is large. Turbulence is expected to be virulent in such a scenario, which will enhance transport and ultimately limit the attainable density gradient (as a maximum transport may be supported). This limiting factor to the gradient naturally leads to consider profile stiffness (Garbet et al. Reference Garbet, Mantica, Ryter, Cordey, Imbeaux, Sozzi, Manini, Asp, Parail and Wolf2004); profiles are stuck to the maximal sustainable gradient values, which are fixed by the point at which a regime of increased turbulence is accessed. Under such a prism, it is not that important what occurs within the strong and weak turbulent regimes, but rather what happens to the transition point. In such a context, support of larger gradients is seen as beneficial, as higher central densities and higher confinement times can be achieved. The key is then to understand the behaviour of this threshold. In practice, this transition point is found through gyrokinetic simulations (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther and Kritz2000) to find when a steep decrease in the nonlinear heat flux is seen when the gradient value is decreased below some threshold value. Recognising an analogous behaviour in Æ, where $\hat {A}_{w} \sim \hat {\omega }_n^3$ below criticality and $\hat {A}_{s} \sim \hat {\omega }_n$ above it, one may postulate that the behaviour in (4.13) is similar to that which one would observe for the common conception of critical gradient. As a consequence of this threshold behaviour, one would then expect that the attained gradient in experiments is the one which we have calculated in (4.13): the system would reside on the weak-to-strong Æ boundary. We do not attempt to verify this relationship here, which for a consistent consideration would require consistent plasma profiles based on (estimated) heat fluxes, which would feedback onto the geometry (e.g. Shafranov shifts). We make no attempt of solving this problem here and leave it to a future investigation, but merely note similarities in the behaviour of our transition threshold and Merlo et al. (Reference Merlo, Brunner, Sauter, Camenen, Görler, Jenko, Marinoni, Told and Villard2015) and Merlo & Jenko (Reference Merlo and Jenko2023) in the increase of gradient thesholds with negative triangularity tokamaks.