Trapped-particle precession and modes in quasisymmetric stellarators and tokamaks: a near-axis perspective

Abstract

This paper presents the calculation of the bounce-averaged drift of trapped particles in a near-axis framework for axisymmetric and quasisymmetric magnetic fields that possess up-down and stellarator symmetry, respectively. This analytic consideration provides important insight on the dependence of the bounce-averaged drift on the geometry and stability properties of the field. In particular, we show that although the maximum-$\mathcal {J}$ property is unattainable in quasisymmetric stellarators, one may approach it through increased plasma $\beta$ and triangular shaping, albeit going through a reduced precession scenario with potentially higher particle losses. The description of trapped particles allows us to calculate the available energy of trapped electrons analytically in two asymptotic regimes, providing insight into the behaviour of this measure of turbulence. It is shown that the available energy is intimately related to magnetohydrodynamics (MHD) stability, providing a potential synergy between this measure of gyrokinetic turbulence and MHD stability.

1. Introduction

It has long been known that the bounce-averaged drift that a trapped particle experiences is central to both linear and nonlinear stability of gyrokinetic trapped-particle modes (Kadomtsev & Pogutse 1967; Rosenbluth 1968; Helander, Proll & Plunk 2013; Helander 2017). This motion of trapped particles can serve as an energy source or sink for various instabilities, and thus their study is central to understanding their behaviour in any plasma-field scenario.

The behaviour of trapped particles depends crucially on the class of fields considered. In an effort to study stellarators, so-called omnigeneous fields (Hall & McNamara 1975; Bernardin, Moses & Tataronis 1986; Cary & Shasharina 1997; Landreman & Catto 2012; Helander 2014) are of particular interest. In such fields, composed of nested flux surfaces on which field lines live, trapped particles have, by definition, a vanishing averaged drift in the direction normal to flux surfaces (i.e. radially). This restricts the dynamics of trapped particles to an average drift within flux surfaces, often referred to as precession and denoted by $\omega _\alpha$. Many authors have investigated the behaviour of this quantity (White & Chance 1984; Roach, Connor & Janjua 1995), and analytic expressions exist for large-aspect ratio tokamaks with circular cross-sections and small non-axisymmetric perturbations (Connor, Hastie & Martin 1983; Hegna 2015). For general omnigeneous stellarators (and even more so, upon relaxing omnigeneity), expressions for $\omega _\alpha$ rarely allow for analytical calculation (Velasco et al. 2023). This ends up impeding the dissection of the underlying physics and effects.

This paper carries out such calculations in a more general scenario. To make the problem tractable, we consider two main simplifications. First, we specialise to quasisymmetric (Boozer 1983; Nührenberg & Zille 1988; Rodriguez, Helander & Bhattacharjee 2020) and axisymmetric fields with stellarator symmetry (Dewar & Hudson 1998) and up-down symmetry, respectively, two special sub-classes of the wider class of omnigenous systems. The central feature of both of these classes is the symmetry of their magnetic field magnitude, $|\boldsymbol {B}|$. This reduces the complexity of the particle dynamics significantly, especially in the region close to the magnetic axis (the centremost field line of the field, around which flux surfaces accrue). This leads to the second simplifying consideration in this paper: the asymptotic description near the axis. In this near-axis approach, the field magnitude may be directly parametrised, and the framework developed by Garren & Boozer (1991b) and Landreman & Sengupta (2019) may be employed directly. Both these simplifying considerations enable an explicit description of the trapped particle motion, whose construction and interpretation we present in §§ 2 and 3.

Once the particle precession is known, we next investigate its role on trapped-particle mode stability. We do so by studying the available energy (Æ) of trapped electrons (Helander 2020); that is, a measure of the available thermal energy that may be liberated by appropriate rearrangements of the electron distribution function. We compute Æ analytically in § 4, explicitly showing its dependence on various important parameters. This enables a direct comparison with other physically relevant properties such as MHD stability and flux surface shaping.

2. Asymptotic expression for the second adiabatic invariant

The description of the trapped particle precession requires the evaluation of the bounce-averaged drift around flux-surfaces, denoted as $\omega _\alpha$. To calculate such a quantity, we must begin by appropriately defining flux surfaces and a notion of direction over them. To this end, we first introduce the Clebsch representation (D'haeseleer et al. 2012) of the magnetic field; namely, $\boldsymbol {B}=\boldsymbol {\nabla }\psi \times \boldsymbol {\nabla }\alpha$, where $\psi$ is the magnetic toroidal flux divided by $2{\rm \pi}$ and $\alpha$ is an angular potential defined as $\alpha = \theta - \iota \varphi$. Here $\theta$ and $\varphi$ are straight-field line (D'haeseleer et al. 2012) poloidal and toroidal angular coordinates, respectively, and $\iota$ is the rotational transform. The flux surfaces are assumed to be nested and correspond to constant pressure surfaces, following a magnetic field that is in equilibrium, $\boldsymbol {j}\times \boldsymbol {B}=\boldsymbol {\nabla } p$. The angle $\alpha$ can be interpreted as a field line label within flux surfaces (following $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\alpha =0$). Thus, we define the precession $\omega _\alpha$ as the rate at which trapped particles change the field line within a flux surface, formally,

(2.1)\begin{equation} \omega_\alpha=\overline{\boldsymbol{v}_D\boldsymbol{\cdot}\boldsymbol{\nabla}\alpha}. \end{equation}

This is the common bounce-averaged binormal drift (Kadomtsev & Pogutse 1967; White & Chance 1984; Helander 2014). Here, the overline operator denotes a bounce-averaging: that is, a time average over the back-and-forth motion of the trapped particle along the field line (thus assuming a ‘thin-orbit’),

(2.2)\begin{equation} \overline{\ldots}=\frac{\oint \cdots \dfrac{\mathrm{d}\ell}{v_\parallel}}{\oint \dfrac{\mathrm{d}\ell}{v_\parallel}}, \end{equation}

where $v_\parallel$ is the parallel velocity, the arc-length along a magnetic field-line is parametrised by $\ell$ and the domain of integration is taken to be a simply connected region that satisfies $v_\parallel (\ell ) \geq 0$ (which is typically referred to as a bounce well). This is an integral at constant $\psi$ and $\alpha$, but also particle energy $H$ and first moment $\mu$.

It is convenient to write the bounce-average drift in (2.1) in terms of derivatives of a single scalar quantity, $\mathcal {J}_\parallel$ (Helander 2014). This scalar quantity is the so-called second adiabatic invariant,

(2.3)\begin{equation} \mathcal{J}_\parallel{=} \int m v_\parallel \, \mathrm{d} \ell, \end{equation}

and is an approximately conserved quantity of trapped particles. Importantly, this quantity serves as the Hamiltonian of the trapped-averaged dynamics of trapped particles, meaning, as can be shown explicitly (Helander 2014; Calvo et al. 2017), that

(2.4)\begin{equation} \omega_\alpha \stackrel{\cdot}{=} - \frac{1}{q} \frac{ \left( \partial \mathcal{J}_\parallel{/} \partial \psi \right)_{\mu,H,\alpha}}{\left( \partial \mathcal{J}_\parallel{/} \partial H \right)_{\mu,\psi,\alpha}} = \frac{\Delta \alpha}{\tau_b}. \end{equation}

Here $q$ is the charge of the particle (which we shall take to be $q=-1$ for electrons), and $\Delta \alpha$ and $\tau _b$ have been defined as the total $\alpha$-excursion and elapsed time in a particle bounce, respectively.

Because $\mathcal {J}_\parallel$ encodes the dynamics of trapped particles in a single scalar expression (and more generally, also allows one to calculate the radial drift), we shall explicitly calculate the asymptotic expression for $\mathcal {J}_\parallel$ as a first step towards finding $\omega _\alpha$.

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 1981). 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 1981; Helander 2014)), 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 1991b; Landreman & Sengupta 2019), 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 1991b; Landreman & Sengupta 2019). 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 (1991b, (1)) or Landreman & Sengupta (2019, (2.15)), and noting that this behaviour goes beyond the particular form of equilibrium assumed (Rodríguez, Sengupta & Bhattacharjee 2022), 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 1991a), and results from a clash between the symmetry and the equilibrium (Rodriguez & Bhattacharjee 2021; Rodríguez et al. 2022). 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 2019). 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 2005, Ch. 7) without violating the ‘thin orbit limit’ of our $\mathcal {J}_\parallel$ calculation.¹

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. 1983; Roach et al. 1995; Helander & Sigmar 2005; Hegna 2015) 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. 2020, § 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.² 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 2019; Rodriguez, Sengupta & Bhattacharjee 2022b). 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. Trapped particle precession

With the second adiabatic invariant constructed, we are in a position to evaluate the bounce-average precession. We remind ourselves that we considered the exact quasisymmetric limit and stellarator symmetry³ (e.g. a tokamak with up-down symmetry) when constructing $\mathcal {J}_\parallel$. Because of the idealised omnigeneous nature of the field, we need not worry about the field-line dependence (i.e. $\alpha$ dependence), as the behaviour is ‘identical’ in all field lines as far as the precession is concerned. This is apparent in (2.24). The formalism presented could however be extended to incorporate a description of said $\alpha$ dependence upon controlled deviations from omnigeneity. We present in Appendix B an explicit estimation of the radial average drift in configurations that only achieve quasisymmetry approximately, providing a previously lacking physically meaningful measure of deviations from quasisymmetry within the near-axis framework, which could aid as an optimisation target (Landreman 2022; Rodriguez, Paul & Bhattacharjee 2022a; Rodriguez et al. 2022b).

Let us nevertheless return to the calculation of precession in our idealised scenario, (2.4). We have almost all that is needed to compute $\omega _\alpha$. The only remaining step is taking partial derivatives of (2.24) with respect to $\psi$ and the particle energy $H$. Acknowledging the dependence of the trapped particle label $k$ on both these variables, the result of this calculation may be written as

(3.1)\begin{equation} \omega_{\alpha}=\omega_{\alpha,0}+\omega_{\alpha,1}+O( r), \end{equation}

where the leading term scales like $1/r$ and $\omega _{\alpha,1}\sim O(r^0)$ (details of this derivation may be found in Appendix C).

The leading order term $\omega _{\alpha,0}$ is

(3.2)\begin{equation} \omega_{\alpha,0}=2 \frac{H \eta}{r B_0}\left(\frac{E(k)}{K(k)}-\frac{1}{2}\right)\stackrel{\cdot}{=}\frac{H \eta}{r B_0}G(k), \end{equation}

which is precisely of the form found by Connor et al. (1983) for a large-aspect-ratio tokamak, without magnetic shear or a pressure gradient. Elements of pressure and shaping are involved in the present approach (although not explicitly) through the next order correction, $\omega _{\alpha,1}$,

(3.3)\begin{equation} \omega_{\alpha,1}= \frac{H \eta}{r B_0}\left[ r \eta \mathcal{G}(k) + \frac{r B_{20}}{\eta}\mathcal{G}_{20}(k) + \frac{r B_{22}^C}{\eta}\mathcal{G}_{22}^C(k)\right], \end{equation}

where the functions $\mathcal {G}$ may be expressed in terms of elliptic integrals,

(3.4a)\begin{gather} \mathcal{G}(k)\stackrel{\cdot}{=}-4 \left(\frac{E(k)}{K(k)}\right)^2 + 2 (3 - 2 k^2) \frac{E(k)}{K(k)} - (1 - 2 k^2), \end{gather}

(3.4b)\begin{gather}\mathcal{G}_{20}(k)\stackrel{\cdot}{=}-2, \end{gather}

(3.4c)\begin{gather}\mathcal{G}_{22}^C(k)\stackrel{\cdot}{=}-4\left[\left(\frac{E(k)}{K(k)}\right)^2 - 2 k^2 \frac{E(k)}{K(k)} + \left(k^2-\frac{1}{2}\right) \right]. \end{gather}

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 1991a; Landreman & Sengupta 2019; Rodríguez 2023). 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. (2022b), see some values of $\eta$ in table 1).⁴

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. (1983) 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 1983). 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 1974; Miller et al. 1998; Hegna 2000) 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 1997), which corresponds in the near-axis limit to increasing $B_{20}$ (the radial derivative of the average $B$) (Landreman & Jorge 2020; Kim, Jorge & Dorland 2021; Rodríguez 2023). 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 (2014, § 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 1962, 1974; Greene & Johnson 1962; Bauer, Betancourt & Garabedian 2012; Freidberg 2014) 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$)⁵ 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 (2023). 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 2023, 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.⁶

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).⁷ 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 2023). 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 (2023). 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 1971; Connor et al. 1983), 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. (1983), 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. 2008; Bader et al. 2021; Velasco et al. 2021).

Different trapped particles are affected differently by pressure as a consequence of the evolving Shafranov shift (Shafranov 1963; Wesson 2011; Rodríguez 2023), 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.⁸ 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 1970; Lortz & Nührenberg 1978; Shafranov 1983; Freidberg 2014; Rodríguez 2023), 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 2021; Rodríguez 2023). 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 (2021). 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 2023), 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.⁹ 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 (2023), 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 1983), using numerical methods detailed by Mackenbach et al. 2023a) 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.¹⁰

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 (2022), 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 2022) 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. 1995), 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 2022) 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. Energy available for trapped particle modes

The preceding study of particle precession was strongly motivated by its central role in driving trapped-particle modes (TPMs). In essence, TPM turbulence arises driven by the trapped-particle precession when it resonates and destabilises the diamagnetic drift wave. Simply put, whenever the trapped particles co-drift with the diamagnetic drift wave, energy may be transferred to the drift wave, driving the instability thereof. As we have learnt, a special feature of axisymmetric and quasisymmetric fields is that $\omega _\alpha$ has, to leading order, a zero crossing. That is, there always exists a subgroup of trapped particles (which includes deeply trapped particles) that co-rotate with the drift wave and thus potentially drive TPMs. To assess how the size of this group and the magnitude of its precession feeds TPMs, a more qualitative treatment is necessary. To perform such an analysis, we delve into the stability of TPMs by studying the available energy of the trapped particles (Helander 2017, 2020).

Available energy (Æ) is an upper bound on the free energy available to the plasma after a constrained rearrangement of the distribution function (called Gardner restacking, see Kolmes, Helander & Fisch 2020), rearrangement that needs to conserve certain dynamical quantities like $\mathcal {J}_\parallel$. Restricting ourselves to the available energy contained in trapped particles, one obtains a proxy measure of nonlinear turbulent activity of TPMs (and upon specialising to electrons, TEMs Proll et al. 2012). This is, nevertheless, a simplified description of the turbulence, in the best of cases a bound, which does not include a discussion of accessibility. That is, although energy may be available, it might not be possible for a system to evolve to such lower energy state and access the available energy, thus the turbulence activity would be over-estimated. The Æ measure is nevertheless useful (Mackenbach, Proll & Helander 2022) and it provides insight into TPMs.

The calculation of available energy for trapped electrons in a flux tube was recently presented by Mackenbach et al. (2022, 2023c). For a flux tube of length $L$ and elliptical cross-section $\Delta \alpha$ and $\Delta \psi$ in $(\psi,\alpha )$-space (principal axes), the available energy in an omnigeneous ($\omega _\psi =0$) field may be written as

(4.1)\begin{equation} A = \frac{1}{2 \sqrt{\rm \pi}} \frac{{\rm \pi} \Delta \psi \Delta \alpha L}{B_0} n_0 T_0 \iint \sum_{\text{wells}(\lambda)} e^{{-}z} z^{5/2} \hat{\omega}_\alpha^2 \mathcal{R} \left[ \frac{\hat{\omega}_*^T}{\hat{\omega}_\alpha} - 1\right] \hat{G}^{1/2} \,\mathrm{d}z \,\mathrm{d}\lambda, \end{equation}

where $\hat {G}^{1/2}=\int (1 - \lambda \hat {B})^{-1/2} \,\mathrm {d} \ell /L$ is the normalised bounce time, $\mathcal {R}$ is the ramp function (indicating that only faster than the diamagnetic drift co-precessing particles contribute to $A$) and the hats denote normalisation of the frequencies $\hat {\omega }=\Delta \psi ~\omega /H$. The integral is performed over $z=H/T$ and $\lambda$, the combination of which constitute all trapped particle energies and classes (the exponential in the integrand is a consequence of the chosen distribution function of which the Æ is to be calculated, a Maxwellian). The sum over wells simply indicates that the available energy is the result of adding the contributions from every well along the flux tube, as trapped particles may be rearranged within each.

Of course, the amount of energy available depends on the size of the flux tube considered; the measure is extensive. To construct an intensive measure, we normalise it to the total thermal energy available in the tube. The total thermal energy in the flux tube for a plasma of temperature $T_0$ and density $n_0$ is (using $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\ell =B$)

(4.2)\begin{equation} E_t = \int \frac{n T}{B} \,\mathrm{d}\psi \,\mathrm{d}\alpha\, \mathrm{d}\ell \approx n_0 T_0 \frac{{\rm \pi} \Delta \psi \Delta \alpha L}{B_0} \int \frac{1}{\hat{B}} \frac{ \mathrm{d} \ell}{L}. \end{equation}

Hence, the normalised Æ becomes

(4.3)\begin{equation} \hat{A} \equiv \frac{A}{E_t} = \left( \int \frac{2 \sqrt{\rm \pi}}{\hat{B}} \frac{\mathrm{d}\ell}{L} \right)^{{-}1} \iint \sum_{\text{wells}(\lambda)} e^{{-}z} z^{5/2} \hat{\omega}_\alpha^2 \mathcal{R} \left[ \frac{\hat{\omega}_*^T}{\hat{\omega}_\alpha} - 1\right] \hat{G}^{1/2} \,\mathrm{d}z \,\mathrm{d}\lambda, \end{equation}

where the normalising factor in front will henceforth be succinctly referred to as $\mathcal {V}=2\sqrt {{\rm \pi} }\int \mathrm {d}\ell /L\hat {B}$. To make further progress, we realise that the integral over $z$ is analytically tractable if one splits up $\hat {\omega }_*^T$ as

(4.4)\begin{equation} \hat{\omega}_*^T=\hat{\omega}_{*,0}^T/z+\hat{\omega}_{*,z}^T, \end{equation}

where $\hat {\omega }_{*,z}^T=-\Delta \psi \partial _\psi \ln T$ and $\hat {\omega }_{\star,0}^T=-\Delta \psi (\partial _\psi \ln n-\tfrac {3}{2}\partial _\psi \ln T)$. To ease the calculation (although it may be extended to the more general case), we shall take the temperature gradient to be zero and consider the limit of a peaked density profile (i.e. $\partial _\psi \ln n$ is negative for $\psi >0$); we are specialising to density-driven trapped-particle instabilities. This assumption makes the diamagnetic drift $\omega _\star$ particle-energy independent, leading to

(4.5)\begin{equation} \hat{A}=\frac{1}{\mathcal{V}}\int\mathrm{d}\lambda \sum_\mathrm{wells(\lambda)}(\hat{\omega}^T_{*,0})^2\hat{G}^{1/2}\mathcal{F}(c_1)\varTheta(\hat{\omega}_\alpha), \end{equation}

where

(4.6)\begin{equation} \mathcal{F}\left(c_1=\frac{\hat{\omega}_{*,0}^T}{\hat{\omega}_\alpha}\right)= \frac{ 2 \sqrt{c_1} \left( 15 + 4 c_1 \right)\exp({-}c_1) + 3 \sqrt{\rm \pi} (2 c_1 - 5) \mathrm{erf}(\sqrt{c_1}) }{8 c_1^2}, \end{equation}

and the $\varTheta$ function is a Heaviside function that vanishes for $\omega _\alpha (\lambda )<0$, which physically represents the inability of counter-rotating particles to drive the TPM.

The function $\mathcal {F}$, see figure 5, may be interpreted as a measure of the coupling of different particle classes to the available energy (ignoring further contributions from the normalised bounce time). With the energy dependence averaged out after the integral over $z$, the comparison between the diamagnetic drift and precession is captured by $c_1$. Both trapped particles that are drifting too fast (i.e. $|c_1| \ll 1$) and which are drifting too slow (i.e. $|c_1| \gg 1$) as compared to the drift wave have a vanishing contribution to the Æ, as $\mathcal {F}\rightarrow 0$. This is a formal statement of the resonance requirement for an effective drive of the drift-wave instability, where the coupling is largest at $c_1 \approx 3.9$.

Figure 5. Plot of the function weighing the contribution to Æ of various trapped particle classes. Plot of $\mathcal {F}$ as a function of the diamagnetic to precession drifts $c_1$. The dashed line indicates the value $c_1 \approx 3.9$, which corresponds to the maximum of $\mathcal {F}$ and, in a sense, represents the subset of particles that most vigorously resonate and drive the drift wave.

To proceed further, and since the expressions for the precession derived in the preceding section depend on the trapping parameter $k$ explicitly, it will be natural to write Æ in (4.5) as an integral over $k$.

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. (2022, 2023c), 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. 2022).¹¹

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.¹² 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.¹³

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.¹⁴ 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 (2020). 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 2022). 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.¹⁵ 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.¹⁶ 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 (2023, (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. 1998) geometry (Mackenbach et al. 2023b) 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 1978; Qi et al. 2019). Concerning higher order effects, the beneficial nature of a pressure gradient on the trapped electron mode has been noted by many authors (Rosenbluth 1968; Connor et al. 1983; Li & Kishimoto 2002). 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. 2019). 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. 2004); 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. 2000) 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. (2015) and Merlo & Jenko (2023) in the increase of gradient thesholds with negative triangularity tokamaks.

5. Conclusions and outlook

In this paper, we explored the trapped-particle precession and its effects on the available energy to trapped particle modes in axisymmetric and quasisymmetric fields. We do so by considering a near-axis description of the fields, in which the magnitude of the magnetic field is directly prescribed and interlinked to other geometric and equilibrium features. As such, this may be regarded as an extension or alternative to previous local considerations (Connor et al. 1983; Roach et al. 1995; Hegna 2015). The precession of trapped particles is constructed explicitly, and analytic expressions are given to leading order and the first-order correction. This allows us to prove the impossibility of the maximum-$\mathcal {J}$ property in such fields to leading order, as follows from Boozer (1983). The role of a pressure gradient in helping to attain this property at a finite radius is presented. Investigating the effect of triangularity in axisymmetric devices, we find that negative triangularity may aid in attaining the maximum-$\mathcal {J}$ property, but the influence on different classes of trapped particles is generally different. In the full quasisymmetric case, closed forms are also provided and some practical examples discussed.

The influence of such precession on density-gradient driven trapped particle modes in the system is then analysed by assessing its effect on the available energy (Helander 2017, 2020; Mackenbach et al. 2023c). Two different asymptotic regimes naturally arise in the form of ‘weakly’ and ‘strongly’ driven regimes, each with a different behaviour and physical mechanism, which furthermore allows one to define a critical transition density gradient. In the weakly driven regime, a narrow band (in $\lambda$-space) of the trapped particle population is responsible for the available energy, whilst in the strongly driven regime, all resonating trapped particles contribute. This physical difference between the two asymptotic regimes leads to a difference in the dependencies of Æ on the field.

This is certainly true for the effects of pressure gradient and triangularity: its beneficial nature depends on the asymptotic regime considered. Interestingly, we find that the dependence of Æ on triangularity is of the same form as that in Mercier's criterion for MHD stability, up to a sign that changes depending on the driving regime. In the strongly driven regime, a synergy is found between improving MHD stability and lowering energy available to the trapped particles through triangularity, meaning that improving one will improve the other (the opposite being true of plasma $\beta$). The reduction in precession of deeply trapped particles behind this behaviour will, whenever deviations from ideal QS exist, lead to increased fast particle losses, at least until a regime of significant precession in the direction of maximum-$\mathcal {J}$ is reached. The synergy between MHD and turbulent activity through triangularity inverts in the weakly driven regime, where it is under plasma $\beta$ that this synergy is observed. This difference in behaviour affects the critical-gradient estimate for the transition between the two regimes. This gradient grows with negative triangularity, which could align with some of the existing observations in advanced tokamak scenarios.

All in all, one finds that the near-axis framework is a convenient model to assess properties of trapped particles in quasisymmetric magnetic fields. The notions of maximum-$\mathcal {J}$, omnigeneity (through the bounce-averaged radial drift) or Æ, for which analytical expression may be found, allows one to readily evaluate several aspects of performance of different stellarator configurations at negligible computational cost (as measured by Æ). This may be helpful as a proxy for turbulence optimisation. Finally, though we have specialised to quasisymmetric configurations, many of the techniques presented may be applied to other contexts (such as quasi-isodynamic fields or Miller geometry models) allowing one to make similar statements. For the quasi-isodynamic case, such an investigation is currently being undertaken and an even more distinct case in which the bounce-averaged radial drift may play a significant role could also benefit from the current framework.