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

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$),

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 normalized 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 Eq. (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,

where $N$ is an integer that defines a helical angle, foreseeing the application to quasi-symmetric devices with a helical symmetrry. Using this angular parameterisation, 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

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.

In order 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 Eq. (1) in Garren & Boozer (1991b) or Eq. (2.15) in Landreman & Sengupta (2019), and noting that this behaviour goes beyond the particular form of equilibrium assumed (Rodríguez et al., 2022), as

where we have separated,

and the second order,

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 quasi-symmetric solutions for equilibria at second order in $r$ (unless exactly axisymmteric fields are considered). This is commonly referred to as the Garren-Boozer oveerdetermination 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 $|\mathbf{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.¹¹1The back-of-the-envelope calculation is as follows. Simplify things by considering the tokamak notation $\eta\sim B_{0}/G_{0}\sim 1/R$ and $q=1/\iota$ for the safety factor. The adiabatic invariant calculation is correct when, along the bounce-orbit, the guiding centre approximately follows the field-line. That is, the argument $v_{\parallel}\mathrm{d}\ell$ should not change significantly over the distance the particle drifts in a bounce time. As $v_{\parallel}\sim v_{t}\sqrt{\epsilon/R}$, the argument changes $\mathrm{d}(v_{\parallel}\mathrm{d}\ell)/\mathrm{d}r\sim\mathcal{J}_{\parallel}/r$. As the transit time is $\omega_{t}^{-1}\sim qR/v_{t}$ and $\mathbf{v}_{d}\cdot\nabla r\sim v_{t}\rho/R$, the error $\Delta\mathcal{J}_{\parallel}/\mathcal{J}_{\parallel}\sim\sqrt{q^{2}\rho^{2}R/r^{3}}$. Thus, $r\gg r_{p}=(q^{2}\rho^{2}R)^{1/3}$, where $r_{p}$ is the range where potato orbits come into play (Helander & Sigmar, 2005).

The way that we have grouped the terms in Eq. (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 $|\mathbf{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,

where for every $\lambda$, the bounce-points $\chi_{b}^{\pm}$ (left and right) of the integral are given by,

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. Secondly, 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 Eq. (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 a left and right values). Introducing this correction, let us rewrite $\mathcal{J}_{\parallel}$ for a stellarator symmetric field in the following form,

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$,

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}$.

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,

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

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}=\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), Sec. 19), which we define as

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})$,

where the functions $I_{1}$ and $I_{2}$ describe the behaviour of different trapped-particle classes via $k$, {subeqnarray} I_1(k) &= 2[ (k^2 - 1) K(k) + E(k) ],
I_2(k) = 23 [(2 k^2-1) E(k)-(k^2-1) K(k)].

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{Hr\eta}$ is directly related to the reduced parallel velocity that trapped particles must have in order for them to be trapped.²²2From energy conservation, one can readily show that $mv_{\parallel}^{2}/2\approx Hr\eta\left(2k^{2}-1\right)\hat{B},$ and hence that $v_{\parallel}\sim\sqrt{Hr\eta}$. The factor $B_{\alpha 0}/\bar{\iota}_{0}B_{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}_{0}B_{0}\sim R_{0}/(\iota_{0}-N)$, where $R_{0}$ is the major radius of the device and $N$ represents the helicity of the $|\mathbf{B}|$ symmetry determined by the shape of the axis (Landreman & Sengupta, 2019; Rodriguez et al., 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 $|\mathbf{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 so 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)}$.

To evaluate the second order term we follow Eq. (14), for which we need the following integral,

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 $|\mathbf{B}|$ for a quasi- and stellarator symmetric field, Eq. (10),

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$,

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

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

Let us start by analysing the leading order precession of trapped particles focusing on $\omega_{\alpha,0}$, Eq. (25). 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$, parameter defined in Eq. (9) as a measure of the leading order variation of $|\mathbf{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}$ then, for a fixed cross-sectional shape, increasing the major radius reduces precession, consequence of the field becoming more straight and the gradients in $|\mathbf{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).⁴⁴4In the language of Rodriguez et al. (2022b), the qualitative difference between the QA and QH configurations holds for cases that live deep in their corresponding quasisymmetric phases. Close to phase transitions, i.e. axis shapes close to having vanishing curvature points (and hence stronger flux surface shaping), these differences fade, and in particular so will any ‘distinct’ $\eta$-behaviour.

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 $\mathbf{v}_{D}\cdot\nabla\theta\sim v_{\nabla B}|\nabla\theta|$, where $v_{\nabla B}$ denotes the gradient drift driven by the radial variation of $B$. As $\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 Eq. (25). A plot of $G$ as a function of $k$ is presented in Fig. 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 Fig. 1(a)), 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 $|\mathbf{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.

Comparing this result against previous investigation, we find, as we already pointed, 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 is 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. C. Mercier & N. Luc (1974); Miller et al. (1998); Hegna (2000)) to discuss precession, which weakens the coupling between $|\mathbf{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.

Let us now focus on the effect of second order elements on precession. In the form presented in Eq. (26), 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 the left plot of Figure 2.

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 $\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_{\nabla B}\sim\mathbf{B}\times\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 et al., 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 Sec. 3.7 of Helander (2014) 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 et al., 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 $|\mathbf{B}|\sim B_{22}^{C}\cos 2\chi$. Naturally, both the most deeply and shallowly trapped particles, who live at $\chi=0,\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$)⁵⁵5In fact, $\delta$ is also normalised by a factor of $r$. This is so because upon approaching the axis the triangular shaping disappears in favour of elliptical cross-sections, and thus triangularity clearly changes with $r$. Because the leading order is proportional to $r$ we normalise it out. 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 in 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 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.⁶⁶6One must be careful with how the sign of $\eta$ is defined in this paper and in the definition of $\delta$. For the derivation of the precession in this paper we assumed $\eta>0$, but used the unusual convention of writing $B=B_{0}(1-r\eta\cos\chi)$, with a minus sign. For the triangularity to have the meaning above, one must include a $\mathrm{sign}(\eta)$ to the definition presented in Rodríguez (2023), which is defined assuming $\eta>0$ in the opposite convention to that considered here.

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).⁷⁷7This statement is almost always true. There are some special cases in which, however, the triangularity is not the most appropriate geometric feature to explicitly involve in the parametrisation of the field, and instead the Shafranov shift (Shafranov, 1963; Wesson, 2011; Rodríguez, 2023) should be chosen. We do not consider this special case. 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.

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 Rodríguez (2023). The result reads,

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 functions 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,

using the definitions in Eqs. (27). 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_{0}I_{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 Eq. (29a) 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 to enhanced fast particle losses (especially of deeply trapped particles) whenever deviations from QS exist (Nemov et al., 2008; Velasco et al., 2021; Bader 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}^{2}I_{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 Fig. 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.⁸⁸8One needs to be careful here, as increasing the pressure gradient will also increase the Shafranov shift, and with it the second order magnetic field shaping. This will start to incur on our asymptotic ordering. Anyhow, equating the different orders still results in a convenient measure. 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, though, 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 Eq. (25), 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

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,

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

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 (Rodríguez, 2023; Freidberg, 2014; Shafranov, 1983; Lortz & Nührenberg, 1978; Solov’ev & Shafranov, 1970)), 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, though, 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, others in the opposite. The most deeply and barely trapped particles are a special case , though, as taking $k=0,~{}1$, $\mathcal{G}_{\delta}=4\bar{\alpha}/(3+\bar{\alpha})>0$ has a maximum value (see Fig. 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 effect in reversing 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,

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 in Landreman (2021). In any case, Eq. (33) 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 Fig. 3 the values of $\mathcal{G}_{\delta}$ and $\mathcal{G}_{p_{2}}$ for a number of QS configurations.⁹⁹9Many of the configurations in Figure 3 were analysed for their MHD behaviour in Rodríguez (2023), and thus make it a suitable set for discussion. We note that in Figure 5 of Rodríguez (2023), the configuration named ‘2022 QA’ was represented through its less vertically elongated cross-section. This is not wrong per se, as one may represent the configuration by any of its cross-sections (and in the framework in the paper, any up-down symmetric one). However, for a discussion on the effect of bean shapes, as in said paper, it was not the appropriate choice, and we point that out here (may it serve as erratum). That this inappropriate choice of a cross-section was made may be seen from an analysis of its cross-sections (see Landreman (2022)), or from the unusual value of $\bar{F}$ in Table 1 in Rodríguez (2023). Here we consistently consider its vertical cross-section (equivalently, the $\varphi=0$ cross-section of the rotated configuration). All other configurations are similarly represented by their most vertically elongated cross-section for consistency. 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).

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 particles against the diamagnetic drift. From the MHD stability analysis of these configurations in 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 Fig. 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-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).

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 are 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 in Mackenbach et al. (2023a)) against the analytical results. For this practical comparison, we employ the definition of $k$ in terms of the bounce point, Eq. (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.¹⁰¹⁰10 In the more general case it is then more convenient to group precession through the class label $\lambda$, normalised as in Roach et al. (1995), namely $k^{2}=(1-\lambda B_{\mathrm{min}})B_{\mathrm{max}}/(B_{\mathrm{max}}-B_{\mathrm{min}})$. This maintains the 0 and 1 values for deeply and barely trapped particles respectively. This definition is only equal to the natural near-axis definition to leading asymptotic order. Thus, in that case, a comparison would require the alternative, yet asymptotically equivalent, definition of $k$ in (68b), which defines it in terms of $\lambda$.

This comparison is shown in Fig. 4, where the bounce-averaged precession is plotted as a function of $k$ and for two radial locations, $\varrho\stackrel{{\scriptstyle\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.

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 Eq. (25). To perform the integral in Eq. (38) 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 Fig. 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.

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$, which the Æ will depend on. 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 on the order of the Larmor radius $\rho$. As such, we write

where $C_{r}$ may formally be dependent on other equilibrium parameters (e.g., the rotational transform $\iota$ or the flux expansion $|\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

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

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

1. 1.

   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 Eq. (38) for reference). One factor of $\varrho^{1/2}$ may be accounted for due to 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 the available energy, as the band near $k_{0}$ scales with $\varrho$.

2. 2.

   Another scaling of interest in Eq. (41) 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, $\widehat{A}_{\mathrm{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. 3.

   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 $\widehat{A}_{\mathrm{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 rearrangment 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, though, a reasonable description in Æ as a measure of heat transport in practice (see Mackenbach et al. (2022)).¹¹¹¹11An appealing alternative to this constant value would perhaps be to consider the poloidal Larmor radius instead of the Larmor radius as the appropriate scale of the flux tube. In that case, we would have an additional factor of aspect ratio and rotational transform, which would change the comparison of behaviour between different fields. Which form is most appropriate is not a closed question.

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 Eq. (41) 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 this 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.¹²¹²12Note that there might be some issues of independence here. We assume the density gradient to be independent from the field, which we know is not true, as the field must satisfy force balance, and the pressure gradient has a density gradient component to it. The large density gradient limit will thus to an extent bring second order $|\mathbf{B}|$ into the picture. However, we may formally proceed in this form. 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.¹³¹³13This is an important point. In the asymptotic sense, the weakly-driven regime will always hold sufficiently close to the axis. However, for a finite radius description, this may become quickly unimportant.

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 Fig. 7). All in all, one finds that the integral reduces to

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 to 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 to the weak regime, $\widehat{A}_{\mathrm{w}}$. Using the same tokamak-estimation for $\eta$ as there, one finds $\widehat{A}_{\mathrm{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.

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