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

E. Rodríguez  1 R.J.J. Mackenbach1,2  1Max Planck Institute for Plasma
Physics, 17491 Greifswald, Germany 2Eindhoven University of Technology, 5612
AZ Eindhoven, The Netherlands

###### 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 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
et al., 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 et al., 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 et
al., 1995), and analytic expressions exist for large-aspect ratio tokamaks
with circular cross-sections and small non-axisymmetric perturbations (Connor
et al., 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 quasi-symmetric (Boozer, 1983; Nührenberg & Zille, 1988;
Rodriguez et al., 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,
$|\mathbf{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); 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 Sections
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 rearrangments of the
electron distribution function. We compute Æ analytically in Section 4,
explicitly showing its dependence on various important parameters. This
enables a direct comparison to 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 quantity, then, 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, $\mathbf{B}=\nabla\psi\times\nabla\alpha$, where $\psi$ is the
magnetic toroidal flux divided by $2\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,
$\mathbf{j}\times\mathbf{B}=\nabla p$. The angle $\alpha$ can be interpreted
as a field line label within flux surfaces (following
$\mathbf{B}\cdot\nabla\alpha=0$). Thus, we define the precession
$\omega_{\alpha}$ as the rate at which trapped particles change field line
within a flux surface, formally,

$\omega_{\alpha}=\overline{\mathbf{v}_{D}\cdot\nabla\alpha}.$ (1)

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

$\overline{\dots}=\frac{\oint\dots\frac{\mathrm{d}\ell}{v_{\parallel}}}{\oint\frac{\mathrm{d}\ell}{v_{\parallel}}},$
(2)

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

$\mathcal{J}_{\parallel}=\int mv_{\parallel}\>\mathrm{d}\ell,$ (3)

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

$\omega_{\alpha}\stackrel{{\scriptstyle\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}}.$ (4)

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

$\mathcal{J}_{\parallel}=\sqrt{2H}\int\sqrt{1-\lambda\hat{B}(\ell)}\>\mathrm{d}\ell.$
(5)

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,

$\chi\stackrel{{\scriptstyle\cdot}}{{=}}\theta-N\varphi,$ (6)

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

$\mathcal{J}_{\parallel}=\frac{\sqrt{2H}}{\overline{\iota}}\frac{B_{\alpha}}{B_{0}}\int\frac{1}{\hat{B}}\sqrt{1-\lambda\hat{B}}\>\mathrm{d}\chi.$
(7)

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

$\hat{B}=\bar{B}+\delta B(\varphi,\chi),$ (8)

where we have separated,

$\overline{B}=1-r\eta\cos\chi,$ (9)

and the second order,

$\delta B=r^{2}\left(B_{20}+B_{22}^{C}\cos 2\chi+B_{22}^{S}\sin 2\chi\right).$
(10)

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.111The 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,

$\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.$ (11)

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

$\frac{1}{\lambda}=\bar{B}(\chi_{b}^{\pm})+\delta B(\chi_{b}^{\pm}),$ (12)

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,

$\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,$ (13)

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

$\mathcal{I}\approx\underbrace{\mathcal{I}(0)}_{\stackrel{{\scriptstyle\cdot}}{{=}}\mathcal{J}_{\parallel}^{(0)}}+\epsilon\underbrace{\left.\frac{\partial\mathcal{I}}{\partial\epsilon}\right|_{\epsilon=0}}_{\stackrel{{\scriptstyle\cdot}}{{=}}\mathcal{J}_{\parallel}^{(1)}}$
(14)

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,

$\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.$ (15)

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

$k^{2}=\sin^{2}\left(\frac{\chi_{b}}{2}\right),$ (16)

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

$\displaystyle K(k)$
$\displaystyle\stackrel{{\scriptstyle\cdot}}{{=}}\int_{0}^{\pi/2}\frac{\mathrm{d}\zeta}{\sqrt{1-k^{2}\sin^{2}\zeta}},$
(17a) $\displaystyle E(k)$
$\displaystyle\stackrel{{\scriptstyle\cdot}}{{=}}\int_{0}^{\pi/2}\sqrt{1-k^{2}\sin^{2}\zeta}\>\mathrm{d}\zeta.$
(17b)

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

$\mathcal{J}_{\parallel}^{(0)}=4\sqrt{Hr\eta}\frac{B_{\alpha
0}}{\bar{\iota}_{0}B_{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],$
(18)

where the functions $I_{1}$ and $I_{2}$ describe the behaviour of different
trapped-particle classes via $k$, 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.222From 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)}$.

###  2.3 The first order correction

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

$\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.$
(19)

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

$\displaystyle\delta B-\delta B_{b}=r^{2}B_{22}^{C}\left(\cos 2\chi-\cos
2\chi_{b}\right).$ (20a)

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

$\mathcal{J}_{\parallel}^{(1)}\approx-4\sqrt{Hr\eta}\frac{B_{\alpha}}{\bar{\iota}B_{0}}\frac{rB_{22}^{C}}{\eta}\mathcal{I}_{22}^{C}(k),$
(21)

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

$\mathcal{I}_{22}^{C}(k)\stackrel{{\scriptstyle\cdot}}{{=}}I_{2}(k)-(2k^{2}-1)I_{1}(k).$
(22)

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

$\mathcal{J}_{\parallel}\approx 4\sqrt{Hr\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\\}.$
(23)

##  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 symmetry333It would in
principle be straightforward to relax this condition. However, it would
require treatment of left and right parts of the well-integrals separately.
Thus, for brevity, we do not present said calculations, and instead focus on
the relevant stellarator-symmetric case. (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
Eq. (23). 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 et al., 2022a, b).

Let us nevertheless return to the calculation of precession in our idealised
scenario, Eq. (4). We have almost all that is needed to compute
$\omega_{\alpha}$. The only remaining step is taking partial derivatives of
Eq. (23) 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

$\omega_{\alpha}=\omega_{\alpha,0}+\omega_{\alpha,1}+O(r),$ (24)

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

$\omega_{\alpha,0}=2\frac{H\eta}{rB_{0}}\left(\frac{E(k)}{K(k)}-\frac{1}{2}\right)\stackrel{{\scriptstyle\cdot}}{{=}}\frac{H\eta}{rB_{0}}G(k),$
(25)

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

$\omega_{\alpha,1}=\frac{H\eta}{rB_{0}}\left[r\eta\mathcal{G}(k)+\frac{rB_{20}}{\eta}\mathcal{G}_{20}(k)+\frac{rB_{22}^{C}}{\eta}\mathcal{G}_{22}^{C}(k)\right],$
(26)

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

$\displaystyle\mathcal{G}(k)$
$\displaystyle\stackrel{{\scriptstyle\cdot}}{{=}}-4\left(\frac{E(k)}{K(k)}\right)^{2}+2(3-2k^{2})\frac{E(k)}{K(k)}-(1-2k^{2}),$
(27a) $\displaystyle\mathcal{G}_{20}(k)$
$\displaystyle\stackrel{{\scriptstyle\cdot}}{{=}}-2,$ (27b)
$\displaystyle\mathcal{G}_{22}^{C}(k)$
$\displaystyle\stackrel{{\scriptstyle\cdot}}{{=}}-4\Bigg{[}\left(\frac{E(k)}{K(k)}\right)^{2}-2k^{2}\frac{E(k)}{K(k)}+\left(k^{2}-\frac{1}{2}\right)\Bigg{]}.$
(27c)

###  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}$, 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).444In 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.

| ARIESCS | ESTELL | GAR | HSX | NCSX | QHS48 | Precise QA | Precise QH
---|---|---|---|---|---|---|---|---
$\eta R_{00}$ | 0.59 | 0.79 | 0.73 | 1.54 | 0.55 | 1.15 | 0.68 | 1.51
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 $\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.

(a) Schematic diamagnetic and gradient drifts

(b) Plot showcasing $G(k)$

Figure 1: Precession of trapped particles near the magnetic axis. Left: the
diagram shows the main drifts in a magnetic configuration with a dominant
$\nabla B$ direction, such as in a tokamak. The electron particle $\nabla B$
drift, $v_{\nabla}$, and the diamagnetic drifts, $v_{*}$, are shown (the
latter proportional to $-\mathbf{B}\times\nabla p$). Resonance between the two
occurs on the outboard side (where the deeply trapped particles live),
defining the bad curvature region. Right: the plot shows 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, 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.

###  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 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.

Figure 2: Effect of second order quantities on precession. (Left) The plot
shows the dependence on the trapped particle class $k$ of the three components
of $\omega_{\alpha,1}$. (Top right) 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. (Bottom right) 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}^{2}I_{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
$\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$)555In 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.666One 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).777This 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.

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

$\omega_{\alpha,1}=\frac{H\eta}{B_{0}r}\left[r\eta\tilde{\mathcal{G}}(k)-\frac{r\mu_{0}p_{2}}{|\eta|B_{0}^{2}}\mathcal{G}_{p_{2}}(k)+\frac{r\delta}{2}\mathcal{G}_{\delta}(k)\right].$
(28)

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,

$\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],$
(29a)
$\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),$
(29b)

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.888One 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

$\frac{\mu_{0}|p_{2}|}{B_{0}^{2}}\sim\frac{\eta}{r(2+f)}\sim\frac{\sqrt{\mathcal{E}}}{R_{0}r(2+f)}.$
(30)

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,

$\beta_{\mathrm{crit}}\stackrel{{\scriptstyle\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}.$
(31)

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

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=\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; some more details,
such as their QS quality, can be found in Appendix E.

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

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

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,

$(r\delta)_{\mathrm{crit}}\sim\frac{3+\bar{\alpha}}{2\bar{\alpha}}.$ (33)

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.999Many 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).

###  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
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.101010 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$.

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, normalized to
$H/r\sqrt{2B_{0}\psi_{\mathrm{edge}}}$, where $a$ is the minor radius of the
configuration, $B_{0}$ the $B$ on axis, and $2\pi\psi_{\mathrm{edge}}$ the
total toroidal flux. There is good correspondence for the precisely
quasisymmetric configurations (Landreman & Paul, 2022), and less so for HSX,
who 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.

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.

##  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 (TPM). 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, though, 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 et al. (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 non-linear 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 et
al., 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,

$\displaystyle A=\frac{1}{2\sqrt{\upi}}\frac{\upi\Delta\psi\Delta\alpha
L}{B_{0}}n_{0}T_{0}\iint$
$\displaystyle\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,$
(34)

where $\hat{G}^{1/2}=\int(1-\lambda\hat{B})^{-1/2}\mathrm{d}\ell/L$ is the
normalized 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
consider; 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 $\mathbf{B}\cdot\nabla\ell=B$)

$E_{t}=\int\frac{nT}{B}\mathrm{d}\psi\mathrm{d}\alpha\mathrm{d}\ell\approx
n_{0}T_{0}\frac{\upi\Delta\psi\Delta\alpha
L}{B_{0}}\int\frac{1}{\hat{B}}\frac{\mathrm{d}\ell}{L}.$ (35)

Hence the normalized Æ becomes

$\widehat{A}\equiv\frac{A}{E_{t}}=\left(\int\frac{2\sqrt{\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,$
(36)

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

$\hat{\omega}_{*}^{T}=\hat{\omega}_{*,0}^{T}/z+\hat{\omega}_{*,z}^{T},$ (37)

where $\hat{\omega}_{*,z}^{T}=-\Delta\psi\partial_{\psi}\ln T$ and
$\hat{\omega}_{\star,0}^{T}=-\Delta\psi(\partial_{\psi}\ln
n-\frac{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,

$\widehat{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})\Theta(\hat{\omega}_{\alpha}),$
(38)

where,

$\mathcal{F}\left(c_{1}=\frac{\hat{\omega}_{*,0}^{T}}{\hat{\omega}_{\alpha}}\right)=\frac{2\sqrt{c_{1}}\left(15+4c_{1}\right)\exp(-c_{1})+3\sqrt{\pi}(2c_{1}-5)\mathrm{erf}(\sqrt{c_{1}})}{8c_{1}^{2}},$
(39)

and the $\Theta$ 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.

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.

The function $\mathcal{F}$, see Fig. 5, may be interpreted as a measure of the
coupling of different particle classes to the available energy (ignoring
further contributions from the normalized 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$.

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 Eq. (38) 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 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).

(a) A schematic of the Æ distribution

(b) A numerical calculation of Æ distribution.

Figure 6: Schematic of the contribution to available energy. Left: 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. Right: a numerical calculation showing the distribution of Æ
across the magnetic well, normalized by 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$.

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

$\Delta\psi=B_{0}r\Delta r=B_{0}rC_{r}\rho$ (40)

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

$\widehat{A}_{\mathrm{w}}\approx\frac{2\sqrt{2}}{9\pi}C_{r}^{2}\rho_{*}^{2}\left(\frac{\hat{\omega}_{n}}{\varrho}\right)^{3}\frac{\varrho^{3}\sqrt{\varrho}}{\sqrt{a\eta}}$
(41)

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)).111111An 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.

###  4.2 A strongly-driven finite-radius asymptotic regime

(a) A schematic of the Æ distribution

(b) A numerical calculation of Æ distribution.

Figure 7: Schematic of the contribution to available energy. Left: diagram
showcasing contribution of $\mathcal{F}$ to the Æ integrand. The broken line
satifies $\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. Right: 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$ is denoted by
a green dashed line, and $\hat{\omega}_{n}=10^{4}$.

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.121212Note 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.131313This 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

$\widehat{A}_{\mathrm{s}}\approx
1.1605~{}\frac{C_{r}^{2}\rho_{*}^{2}}{\pi\sqrt{2}}\frac{\hat{\omega}_{n}}{\varrho}(a\eta)^{3/2}\varrho^{3/2}.$
(42)

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,

$\left.\frac{a}{L_{n}}\right|_{\mathrm{crit}}\approx 1.61591~{}a\eta.$ (43)

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

$\varrho_{\mathrm{crit}}\approx 0.80795~{}a\eta.$ (44)

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 Fig. 8, where we compare both asymptotic regimes in one
device, and the weakly asymptotic regime across multiple devices141414The code
is freely available on https://github.com/RalfMackenbach/AEpy. We observe
excellent agreement in the asymptotic behaviour between the found results and
the numerical calculation.

(a) Precise QA

(b) Multiple configurations

Figure 8: A comparison of the numerical calculation of Æ against the
analytical result. Left: 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 in Landreman & Jorge (2020). The plot has been generated with
a gradient value of $\hat{\omega}_{n}/\varrho=10^{3}$ for visualisation
purposes. Right: 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 $\widehat{A}$,
whereas the abscissa corresponds to the asymptotic result of Eq. (41),
$\widehat{A}_{\mathrm{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$

### 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 $\widehat{A}$. We show
how to do this to obtain the dependence of $\widehat{A}$ on $p_{2}$ and
$\delta$ to leading order at the end of Appendix F.151515Note that to do so it
is not necessary to compute the whole next order correction to $\widehat{A}$,
but only the pieces concerned with the second order parameters directly. After
such considerations, we may write
$\widehat{A}\approx\widehat{A}_{0}+\widehat{A}_{1}+\dots$, where
$\widehat{A}_{0}$ is the leading order expression and $\widehat{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

$\left.\frac{\widehat{A}_{1}}{\widehat{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_{\alpha
0}}{B_{0}\bar{\iota}_{0}}\right)^{2}\right]+\frac{3}{2}\frac{1-\bar{\alpha}}{3+\bar{\alpha}}r\delta\right),$
(45)

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.161616The change in
pressure without a change in the density gradient may appear impossible, but
it may be achieved straightforwardly by keeping the density gradient fixed and
increasing the ion temperature. 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. But 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 assess, formally, 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. And 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 Fig. 1(b)). And 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 Eq. (45) 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 Eq. (110) for $\mathcal{R}_{\delta}$ directly to Eq.
(4.2) in Rodríguez (2023) for $\mathcal{T}_{\delta}$. This intimate relation
between MHD stability and what may be interpreted as TPM activity has been
observed in many occasions (in fact, could be interepreted 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
$\widehat{A}_{1}/\widehat{A}_{0}|_{\mathrm{strong}}\approx-2.85\widehat{A}_{1}/\widehat{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, then, 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 their contribution to confinement
is. Unlike in the weak regime, the whole trapped population becomes important,
and not just a narrow-band, Fig. 6. In the strongly-driven regime then 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.

Besides 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 to 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.
$\widehat{A}_{1}/\widehat{A}_{0}\sim 1$, we state that we have a critical
scenario. At the edge, the critical $\beta$ becomes
$\beta_{\mathrm{crit}}^{\AE{}}\sim 2a|\eta|/\mathcal{R}_{20}(1+f)$, with $f$
as defined before (with its QS generalisation, which may be found in Eq.
(76a)). This shows that plasma $\beta$ becomes effective in significantly
changing Æ in the strong regime for QAs in the regime of
$\beta_{\mathrm{crit}}^{\AE{}}\sim 2-3\%$, while QH
$\beta_{\mathrm{crit}}^{\AE{}}\sim 4-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 Eq. (33)
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, and range
$(r\delta)_{\mathrm{crit}}\sim 0.4-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,

$\displaystyle\left.\frac{a}{L_{n}}\right|_{\mathrm{crit}}\approx\left.\frac{a}{L_{n}}\right|_{\mathrm{crit},0}\Bigg{[}$
$\displaystyle 1-2.6389\times$ (46)
$\displaystyle\left.\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].$

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 tokamak 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 that 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, though, 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, then, the prediction that increasing the vertical elongation in
tokamaks improves transport agrees with existing knowledge (Chu et al., 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 on 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 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 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 non-linear heat-flux is seen when
the gradient value is decreased below some threshold value. Recognizing an
analogous behaviour in Æ, where
$\widehat{A}_{\mathrm{w}}\sim\hat{\omega}_{n}^{3}$ below criticality and
$\widehat{A}_{\mathrm{s}}\sim\hat{\omega}_{n}$ above it, one may postulate
that the behaviour in Eq. (46) is similar to that 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 (46): 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); 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; Connor et al., 1983; 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 a “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 of 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 quasi-symmetric 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 quasi-symmetric configurations, many of the techniques
presented may be applied in 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.

## Data availability

The data that support the findings of this study are openly available at the
Zenodo repositories with DOI/URL 10.5281/zenodo.8185571 and
10.5281/zenodo.8199904.

## Acknowledgements

The authors would like to acknowledge fruitful discussion with Per Helander,
Jaime Caballero and Iván Calvo. E. R. was supported by a grant by Alexander-
von-Humboldt-Stiftung, Bonn, Germany, through a postdoctoral research
fellowship. R.M. is partly supported by a grant from the Simons Foundation
(560651, PH), and through the project “Shaping turbulence—building a framework
for turbulence optimisation of fusion reactors,” with Project No.
OCENW.KLEIN.013 of the research program “NWO Open Competition Domain Science”
which is financed by the Dutch Research Council (NWO).

##  Appendix A Calculation of the second adiabatic invariant

In this appendix we write out the full derivation of the second adiabatic
invariant $\mathcal{J}_{\parallel}$. Our starting point is the leading order
integral given in Eq. (15) after a straightforward expansion in the ordering
parameter $\epsilon$, defined in the main text.

To make progress with this integral, start by defining a so-called trapping
parameter $k$, which relates to the pitch-angle parameter $\lambda$ via

$\lambda=\frac{1}{1+r\eta(2k^{2}-1)+\delta B_{b}}.$ (47)

Such a definition gives $k\in[0,1]$, with the limits corresponding to deeply
and barely trapped particles respectively. It must be noted that despite its
simple appearence, this definition hides complexity in the trapped particle
class dependence of $\delta B_{b}$, Eq. (12, the field perturbation at the
bouncing point. With this definition, we rewrite the integral

$\mathcal{J}_{\parallel}^{(0)}=\underbrace{\sqrt{2H}\frac{B_{\alpha}}{\bar{\iota}B_{0}}\sqrt{r\eta\lambda}}_{\stackrel{{\scriptstyle\cdot}}{{=}}\bar{\mathcal{J}}}\int\frac{\sqrt{2k^{2}-1+\cos\chi}}{1-r\eta\cos\chi}\mathrm{d}\chi.$
(48)

This integral can be cast into a form close to the one required for elliptic
integrals by employing the double-angle identity
$\cos\chi=1-2\sin^{2}\left(\chi/2\right)$, with which the integral reduces to

$\mathcal{J}_{\parallel}^{(0)}=2\sqrt{2}\bar{\mathcal{J}}\int\frac{\sqrt{k^{2}-\sin^{2}\overline{\chi}}}{1-r\eta(1-2\sin^{2}\overline{\chi})}\mathrm{d}\overline{\chi},$
(49)

where the integration variable has been set to $\overline{\chi}=\chi/2$. The
final subtlety that one needs to take into account is that the limits of
integration are currently set by the zeros of the argument of the numerator in
the integrand (namely, the bouncing points), which also gives an intuitive
(and equivalent) definition of the trapping parameter $k$,

$k^{2}=\sin^{2}\left(\frac{\chi_{b}}{2}\right),$ (50)

where we remind the reader that $\chi_{b}$ denotes the bounce angle. This
shows, as did Eq. (47), that $k$ includes some of the higher order corrections
to $B$. This arises naturally from the integration, and importantly, preserves
the meaning of deeply and barely trapped particles at $k=0,~{}1$, regardless
of the perturbation.

A final substitution puts the integral into the standard form required for
elliptic integrals. Define,

$k\sin\zeta=\sin\overline{\chi}\implies\mathrm{d}\overline{\chi}=\frac{\sqrt{k^{2}-\sin^{2}\overline{\chi}}}{\sqrt{1-\sin^{2}\overline{\chi}}}\mathrm{d}\zeta,$
(51)

in which case the bounce points simply become $\zeta=\pm\pi/2$, independent of
$k$. This transformation is well-defined because $k\in[0,1]$. The leading
order contribution to the second adiabatic invariant is now equal to

$\mathcal{J}_{\parallel}^{(0)}=2\sqrt{2}\bar{\mathcal{J}}\int_{-\pi/2}^{\pi/2}\frac{1}{1-r\eta(1-2k^{2}\sin^{2}\zeta)}\frac{k^{2}\cos^{2}\zeta}{\sqrt{1-k^{2}\sin^{2}\zeta}}\mathrm{d}\zeta.$
(52)

As part of the asymptotic near-axis treatment, $r$ is to be considered small,
and we may thus expand the non-singular denominator of the integrand in powers
of $r$. Including terms up to the first order,171717If one wants to include
higher order effects into the presented calculation, the Boozer Jacobian must
be treated with some care. The reason is that, after all, the Jacobian
represents the geometry of flux surfaces, and thus, the form of the Jacobian
itself depends on how the near-axis expansion is interpreted. If the near-axis
expansion is treated as a model that is truncated at second order to construct
flux surfaces, then the Jacobian is actually not equal to $B_{\alpha}/B^{2}$
beyond the first order. The integral at second order would need to be
reconsidered. A comment in this regard may be found in Landreman (2021)., we
find

$\mathcal{J}_{\parallel}^{(0)}=2\sqrt{2}\bar{\mathcal{J}}\left[I_{1}(k)+I_{2}(k)r\eta+O(r^{2})\right],$
(53)

where we have introduced,

$\displaystyle I_{1}$
$\displaystyle=\int_{-\pi/2}^{\pi/2}\left(\frac{k^{2}-1}{\sqrt{1-k^{2}\sin^{2}\zeta}}+\sqrt{1-k^{2}\sin^{2}\zeta}\right)\mathrm{d}\zeta$
(54) $\displaystyle=2\left[(k^{2}-1)K(k)+E(k)\right]$ $\displaystyle I_{2}$
$\displaystyle=\int_{-\pi/2}^{\pi/2}\frac{k^{2}\cos^{2}\zeta\left(1-2k^{2}\sin^{2}\zeta\right)}{\sqrt{1-k^{2}\sin^{2}\zeta}}\mathrm{d}\zeta$
$\displaystyle=\frac{2}{3}\left[(2k^{2}-1)E(k)-(k^{2}-1)K(k)\right]$

and define complete elliptic integrals of the first and second kind (Olver et
al., 2020, Sec. 19),

$\displaystyle K(k)$
$\displaystyle\stackrel{{\scriptstyle\cdot}}{{=}}\int_{0}^{\pi/2}\left(\sqrt{1-k^{2}\sin^{2}\zeta}\right)^{-1}\mathrm{d}\zeta,$
(17a) $\displaystyle E(k)$
$\displaystyle\stackrel{{\scriptstyle\cdot}}{{=}}\int_{0}^{\pi/2}\sqrt{1-k^{2}\sin^{2}\zeta}\>\mathrm{d}\zeta.$
(17b)

Our final step is to expand $\bar{\mathcal{J}}$ to the required order in $r$.
This expansion is readily done and one can show that, neglecting terms of
order $O(r^{2})$, this reduces to

$\displaystyle\bar{\mathcal{J}}\approx\sqrt{2Hr\eta}\frac{B_{\alpha
0}}{\bar{\iota}_{0}B_{0}}\left[1+r\eta\left(\frac{1}{2}-k^{2}\right)+O(r^{2})\right].$
(55)

Collecting all the terms of order $O(r)$ in $\mathcal{J}_{\parallel}^{(0)}$,

$\mathcal{J}_{\parallel}^{(0)}=4\sqrt{Hr\eta}\frac{B_{\alpha
0}}{\bar{\iota}_{0}B_{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].$
(56)

We now turn to the next order correction in $\epsilon$ following Eq. (14),

$\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.$
(57)

The second term in the integrand is a factor $r$ smaller than the first, and
hence, for our current expansion, only the first term needs to be taken into
account. To perform that remaining integral, we ought to express the integrand
explicitly as a function of the integration variable $\chi$. Using the near-
axis expansion form of $B$ for a quasi- and stellarator symmetric fields, we
know that

$\displaystyle\delta B-\delta B_{b}=r^{2}B_{22}^{C}\left(\cos 2\chi-\cos
2\chi_{b}\right),$ (58a) $\displaystyle 1-\lambda(\bar{B}+\delta
B_{b})=\frac{r\eta(2k^{2}-1+\cos\chi)}{1+r\eta(2k^{2}-1)+\delta B_{b}}.$ (58b)

With these, introducing the trapping parameter, using double angle identities,
and employing $\zeta$ as the integration parameter (like in the previous
order), the integral reduces to

$\displaystyle\mathcal{J}_{\parallel}^{(1)}$
$\displaystyle\approx\bar{\mathcal{J}}\frac{rB_{22}^{C}}{\eta}\times\sqrt{2}\int_{0}^{\pi/2}\frac{\cos
4\bar{\chi}_{b}-\cos 4\bar{\chi}}{\sqrt{1-k^{2}\sin^{2}\zeta}}\mathrm{d}\zeta$
(59)
$\displaystyle\approx-4\sqrt{Hr\eta}\frac{B_{\alpha}}{\bar{\iota}B_{0}}\frac{rB_{22}^{C}}{\eta}\mathcal{I}_{22}^{C}(k),$

where we have kept leading order terms in $r$. The function describing the
behaviour of different trapped particle classes is $\mathcal{I}_{22}^{C}(k)$,
which we have defined as

$\mathcal{I}_{22}^{C}(k)\stackrel{{\scriptstyle\cdot}}{{=}}I_{2}(k)-(2k^{2}-1)I_{1}(k).$
(60)

Expanding to order $r^{2}$ and combining the first and second order result,
Eq. (14), gives us our final expression for the second adiabatic invariant,

$\mathcal{J}_{\parallel}\approx 4\sqrt{Hr\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\\}.$
(23)

Figure 9 shows the comparison of this analytic expression to the numerical
calculation of $\mathcal{J}_{\parallel}$.

Figure 9: Residual between numerical $\mathcal{J}_{\parallel}$ and analytic
approximation. The plot shows, in log scale, the difference between the
numerically computed $\mathcal{J}_{\parallel}$ and the analytic expression,
Eq. (23), to $O(r^{1/2})$ (solid line) and $O(r^{3/2})$ (broken line). The
dotted line shows a reference $\sim r^{5/2}$ scaling, which agrees with the
broken line as predicted by the theory. This particular case was run for an
artificial ideal second-order near-axis $|\mathbf{B}|$ profile with $\eta=1$,
$B_{20}=1$, $B_{22}^{C}=3$ and $k=0.5$.

##  Appendix B Integral expressions for the radial drift

In this section, we use the near-axis framework to find expressions for the
bounce-averaged radial drifts to leading order in the near-axis expansion. Of
course, only if the system is not exactly quasisymmetric will the radial drift
be non-vanishing. We will assume that quasisymmetry is broken at second order
in the near-axis expansion.

As is well-known, it is generally not possible to guarantee the symmetry of
$|\mathbf{B}|$ to second-order in the expansion. Thus, generally, to form a
consistent description to second-order one formally allows $B_{20}$ to be a
function of $\varphi$, rather than a constant. This is the conventional
choice, and keeps the other pieces of $|\mathbf{B}|$ constant. Let us then
consider here a field that is quasisymmetric to first order, but, at second
order, has $B_{2}=B_{20}(\varphi)+B_{22}^{C}\cos 2\chi$, which retains
stellarator symmetry for an even $B_{20}(\varphi)$. Including the
contributions from the other terms if their $\varphi$-independence was relaxed
would also be straightforward.

Let us then consider the $\alpha$ derivative of the second adiabatic invariant
needed for assessing the averaged radial drift explicitly,

$\partial_{\alpha}\mathcal{J}_{\parallel}=\sqrt{2H}\frac{B_{\alpha}}{B_{0}\bar{\iota}}\partial_{\alpha}\left[\int\frac{\sqrt{1-\lambda\hat{B}}}{\hat{B}}\mathrm{d}\chi\right].$
(61)

The integrand vanishes at the bouncing points by construction, and therefore,
by Leibniz’s rule, the boundary terms may be dropped when enacting the
derivative of the integral. This is not completely true for barely and deeply
trapped particles, as their bounce points may change in a non-smooth way as we
move from field line to field line. To picture this, think of a top of a
magnetic well coming down on one side of the well as we move to a different
field line. The original barely-trapped particle ‘leaks’, undergoing a sudden
change in behaviour, leading to a new class of trapped particle. Such a
behaviour cannot be appropriately captured in a perturbative sense, but the
importance of this particle ‘leak’ may be assessed by constructing a measure
of the leaked particle fraction
$f_{\mathrm{leak}}=\max_{\alpha}[|B_{20}((\pi-\alpha)/\bar{\iota})-B_{20}((\pi+\alpha)/\bar{\iota})|]/(B_{\mathrm{max}}-B_{\mathrm{min}})$,
where we assumed stellarator symmetry (i.e.,
$B_{20}(-\varphi)=B_{20}(\varphi)$ at second order).

With this in mind, and ignoring this fringing case, we may pull the derivative
through inside the integral. The integral is taken along field lines (i.e., at
constant $\alpha$), which meand that the Boozer toroidal angle
$\varphi=-(\alpha-\chi)/\bar{\iota}_{0}$ becomes a function of both $\alpha$
and $\chi$. That means that
$\partial_{\alpha}f(\varphi)=-\partial_{\varphi}f/\bar{\iota}_{0}$, which
keeping the leading order near-axis term and using the same notation as in
Appendix A gives,

$\displaystyle\partial_{\alpha}\mathcal{J}_{\parallel}=$
$\displaystyle-\sqrt{2H}\frac{B_{\alpha}}{B_{0}\bar{\iota}}\frac{\lambda}{2}\int\frac{\partial_{\alpha}\hat{B}}{\sqrt{1-\lambda\hat{B}}}\mathrm{d}\chi$
(62) $\displaystyle=$
$\displaystyle\sqrt{2H}\frac{B_{\alpha}}{B_{0}\bar{\iota}_{0}^{2}}\frac{r}{\eta}\sqrt{\frac{\eta
r}{2}}\overbrace{\int_{-\pi/2}^{\pi/2}\frac{B_{20}^{\prime}[\varphi(\bar{\chi},\alpha)]}{\sqrt{1-k^{2}\sin^{2}\zeta}}\mathrm{d}\zeta}^{\stackrel{{\scriptstyle\cdot}}{{=}}\mathcal{I}_{20,\alpha}}$
(63)

where the prime indicates derivative respect to the only argument of $B_{20}$,
$\varphi$. The integration variable is $\zeta$, and one should interpret
$\overline{\chi}(\zeta,k)=\arcsin\left(k\sin\zeta\right)$.

The integral may be readily evaluated using standard numerical methods. To
find $\omega_{\psi}$, the bounce-averaged radial drift, though, we need to
evaluate the leading order bounce time, $\tau_{b}$. Using the expression for
$\mathcal{J}_{\parallel}^{(1)}$, the bounce time is equal to

$\tau_{b}=\partial_{H}\mathcal{J}_{\parallel}\approx\frac{B_{\alpha
0}}{B_{0}\bar{\iota}_{0}}\frac{2K(k)}{\sqrt{Hr\eta}}$ (64)

Hence, the bounce-averaged radial drift to leading order can readily be found
to be

$\omega_{\psi}=-\frac{Hr^{2}}{\bar{\iota}_{0}}\frac{\mathcal{I}_{20,\alpha}(k,\alpha)}{2K(k)}.$
(65)

The averaged radial-drift $\omega_{\psi}$ serves as a physically meaningful
measure of the quasisymmetry quality of a configuration in this near-axis
construction, vanishing when it is omnigeneous. The expression for
$\omega_{\psi}$ is however a function of both $\alpha$ and $k$, and thus, for
a single scalar measure that characterises the radial drift performance of a
field at a given flux surface we must reduce it. Note that given the
periodicity of $B_{20}$, the average of $\omega_{\psi}$ over all field lines
(i.e., $\alpha$) vanishes. This might suggest that there is no net detrimental
effect to having this radial drift, as on average, there is the ‘same’ amount
of particles going one way or the other. However, the neoclassical transport
in the low collisionality limit as measured by $\epsilon_{\mathrm{eff}}$ is
insensitive to the sign of $\omega_{\psi}$. To find a single measure then of
its magnitude, then, we attempt to find an upper bound for $\omega_{\psi}$.

Noting that the denominator in the integrand of $\mathcal{I}_{20,\alpha}$ is
an always positive function within the integration domain,

$\mathcal{I}_{20,\alpha}\leq 2B_{20,\mathrm{max}}^{\prime}K(k).$ (66)

Thus,

$\omega_{\psi}\leq\frac{r^{2}H}{\bar{\iota}_{0}}B_{20,\mathrm{max}}^{\prime},$
(67)

with the equality only holding when the field line label $\alpha$ makes the
largest gradient $B_{20,\mathrm{max}}^{\prime}$ match the bottom of the well,
and deeply trapped particles are considered (see Fig. 10). Only in this limit
the particle samples the largest non-QS value of $B_{20}^{\prime}$. Any other
value will necessarily be smaller.

This bound provides a simple relation between the derivative of $B_{20}$ and
the radial drift of particles. The derivative of $B_{20}$ is thus a more
physical form of measuring quasisymmetry breaking within the near-axis
framework compared to simply using the peak-to-peak $B_{20}$ variation as it
is customary (Landreman & Sengupta, 2019; Landreman, 2022; Rodriguez et al.,
2022b).

Figure 10: Radial drift of precise QH in the near-axis description. The plot
shows the $B_{20}$ function of the precise QH near-axis construction as a
function of $\varphi$, and (right) the bounce-average radial drift
$\omega_{\psi}$ as a function of $\alpha$ and $k$. The maximum bounce-average
radial drift $\omega_{\psi}$ occurs for the deeply-trapped population at the
field line where the largest gradient of $B_{20}$ aligns with the minimum of
the well.

##  Appendix C Full expressions for the particle precession

In this appendix we present some key elements and expression for the
calculation of the precession of trapped particles,
$\omega_{\alpha}=\partial_{\psi}\mathcal{J}_{\parallel}/\partial_{H}\mathcal{J}_{\parallel}$.
Obtaining such expressions from the asymptotic form of
$\mathcal{J}_{\parallel}$ is straightforward, but it requires taking care of
partial derivatives appropriately.

To that end, let us remind ourselves of the set of independent variables:
$\psi$, $\alpha$, $H$ and $\mu$. Because we are using the toroidal flux over
$2\pi$ as our flux surface label,

$\displaystyle\frac{\partial
r}{\partial\psi}=\partial_{\psi}\left(\sqrt{\frac{2\psi}{B_{0}}}\right)=\frac{1}{B_{0}r},.$
(68a) Furthermore, as shown in Appendix A, the pitch angle is related to the
trapping parameter via $\lambda=\frac{1}{1+r\eta(2k^{2}-1)+\delta B_{b}}.$
(68b) We thus require expressions for $\delta B_{b}$, which may readily be
computed as $\displaystyle\delta B_{b}$
$\displaystyle=r^{2}\left(B_{20}+B_{22}^{C}\cos 2\chi_{b}\right)$ (68c)
$\displaystyle=r^{2}\left[B_{20}+B_{22}^{C}\left(1-8k^{2}+8k^{4}\right)\right].$
With this expression at hand, it is straightforward to compute the radial
derivative of $k$ (at constant $\mu$ and $H$, and thus constant $\lambda$),
$\displaystyle\frac{\partial k}{\partial
r}\approx\frac{1-2k^{2}}{4kr}+\frac{B_{22}^{C}-B_{20}}{2k\eta}.$ (68d)
Similarly the derivative with respect $H$ can be found as
$\displaystyle\frac{\partial k}{\partial
H}\approx\frac{1}{4kHr\eta}+\frac{(2k^{2}-1)(\eta^{2}-4B_{22}^{C})}{4kH\eta^{2}}.$
(68e)

The particle class label $k$ changes both with radius and particle energy, as
it follows from the change in $|\mathbf{B}|$.

We first use the above expressions to find the total bounce-averaged excursion
in $\alpha$, $\Delta\alpha=\partial_{\psi}\mathcal{J}_{\parallel}$. Using the
chain rule and keeping the leading two non-zero orders,

$\displaystyle\Delta\alpha=$
$\displaystyle\sqrt{\frac{H\eta}{r}}\frac{B_{\alpha
0}}{\bar{\iota}_{0}B_{0}}\frac{1}{B_{0}r}\Bigg{\\{}2\left[2E(k)-K(k)\right]+$
(69) $\displaystyle+r\eta\left[(2-4k^{2})E(k)+(1+2k^{2})K(k)\right]-$
$\displaystyle-\frac{4rB_{20}}{\eta}K(k)+\frac{4rB_{22}^{C}}{\eta}\left[(2k^{2}-1)E(k)+(1-k^{2})K(k)\right]\Bigg{\\}}.$

Next, the bounce time $\tau_{b}=\partial_{H}\mathcal{J}_{\parallel}$ is

$\displaystyle\tau_{b}=$ $\displaystyle\frac{B_{\alpha
0}}{B_{0}\bar{\iota}_{0}}\frac{1}{\sqrt{Hr\eta}}\Bigg{\\{}2K(k)+$ (70)
$\displaystyle+r\eta\left[4E(k)+(2k^{2}-3)K(k)\right]+\frac{4rB_{22}^{C}}{\eta}\left[E(k)-k^{2}K(k)\right]\Bigg{\\}}.$

The trapped-particle precession is now readily found by taking the ratio of
$\Delta\alpha/\tau_{b}$ and expanding to include the leading-order and the
first correction terms. Doing the expansion, one finds this to be equal to

$\displaystyle\omega_{\alpha}=$
$\displaystyle\frac{H\eta}{B_{0}r}\Bigg{\\{}2\frac{E(k)}{K(k)}-1+$ (71)
$\displaystyle+r\eta\left[-4\left(\frac{E(k)}{K(k)}\right)^{2}+2(3-2k^{2})\frac{E(k)}{K(k)}+(2k^{2}-1)\right]-$
$\displaystyle-\frac{2rB_{20}}{\eta}+\frac{2rB_{22}^{C}}{\eta}\Bigg{[}-2\left(\frac{E(k)}{K(k)}\right)^{2}+4k^{2}\frac{E(k)}{K(k)}+(1-2k^{2})\Bigg{]}+O(r^{2})\Bigg{\\}},$

which is equivalent to Eq. (24).

##  Appendix D Precession in terms of triangularity and pressure gradient

In the main text and Appendix C we showed how to get an expression for the
precession of trapped particles at second order in the near-axis expansion.
This involved the parameters $B_{22}^{C}$ and $B_{20}$, natural parameters in
the near-axis expansion which directly describe the shape of $|\mathbf{B}|$.
Although most natural in the calculations of trapped-particle precession, they
do however not offer a clear connection to the physical nor geometric features
of the field.

Within the near-axis framework, though, such a connection can be made. At
second order in the near-axis expansion of an exactly quasisymmetric,
stellarator-symmetric configuration, $B_{22}^{C}$ and $B_{20}$ can be re-
expressed as linear combinations of pressure gradient, $p_{2}$, and
triangularity of an up-down cross-section, $\delta$. The latter two can then
be used as the independent parameter choices at second order in the expansion
of the field.

The definition of pressure gradient is straightforward as the flux derivative
of the pressure on axis, $p_{2}=(B_{0}/2)\mathrm{d}p/\mathrm{d}\psi$. The
notion of triangularity, $\delta$, requires additional care. We will define
$\delta$ as the relative displacement of the vertical turning points of the
cross-section to the width of the cross section along the up-down symmetry
line (normalised by $r$), and do so in the $(\kappa,~{}\tau)$ Frenet frame.
That is, in the plane orthogonal to the magnetic axis where the cross-section
is up-down symmetric. Asymptotically, and in terms of the near-axis expansion
quantities (Landreman & Sengupta, 2019; Rodríguez, 2023),

$\delta=2~{}\mathrm{sign}(\eta)\left(\frac{Y_{22}^{S}}{Y_{11}^{S}}-\frac{X_{22}^{C}}{X_{11}^{C}}\right).$
(72)

Here the $X$ and $Y$ expansion parameters describe the flux surface shape in
the Frenet-Serret frame, and details on them may be found in Landreman &
Sengupta (2019). Note that dimensionally, the definition in Eq. (72) has units
of inverse length. This is so because $\delta$ is defined not as
triangularity, but rather as $1/r$ times triangularity, which accounts for the
fact that close to the magnetic axis, where cross-sections are elliptical, the
triangularity vanishes. This notion of normalised triangularity in the plane
normal to the axis, as explained in Rodríguez (2023), is generally different
to the common definition of triangularity in the lab-frame,
$\delta_{\mathrm{lab}}$. That is, it is not equivalent to the triangularity of
the cross-section that results from making a cut of the configuration at a
constant cylindrical angle. If the magnetic axis has a relative tilt respect
to the cylindrical coordinate system, then
$\delta_{\mathrm{lab}}=\delta+\Lambda$, where $\Lambda$ is a term that depends
only on the axis shape and first order near-axis shaping.181818We note that
the expression in Eq. (C2) in Rodríguez (2023) for $\Lambda$ is incorrectly
simplified, as it assumes certain symmetry that does generally not exist. The
correct expression will be presented in a future publication, but makes no
difference to the discussion in the present paper. In the special case of an
axisymmetric field, this geometric transformation term $\Lambda$ vanishes. In
general, though, varying $\delta_{\mathrm{lab}}$ or $\delta$ is equivalent
other near-axis features kept constant.

With the notions of triangularity and pressure gradient in place, the
equilibrium field of a quasisymmetric configuration can be uniquely defined at
second order Rodríguez (2023). In the axisymmetric limit this is evident
following the behaviour of the Grad-Shafranov equation and the set-up of its
solutions. In the general quasisymmetric configuration this is more subtle, as
the cross-section shapes change around the torus driven by the asymmetry of
the magnetic axis. Specifying the triangularity of a single cross-section
might then appear insufficient to describe uniquely the whole field, but the
conditions of quasisymmetry and equilibrium are sufficient to grant this
uniqueness. Schematically, one may picture the situation as a kind of initial
value problem in which the single cross-section is the initial value, and the
axis-shape describes the flow of evolution. There are therefore different ways
of describing the very same field, as different cross-sections may be chosen.

It should appear clear that the shape of the axis thus plays a crucial role in
the problem, especially through its asymmetry. Described by its curvature,
$\kappa$, and torsion, $\tau$, it is natural to construct a measure of said
asymmetry like $\bar{F}$, introduced in Rodríguez (2023), and defined to be,

$\bar{F}=2\left[\frac{(I_{2}-B_{0}\tau)/\kappa^{2}}{\int_{0}^{2\pi}\mathrm{d}\varphi(I_{2}-B_{0}\tau)/\kappa^{2}}\frac{\int_{0}^{2\pi}\mathrm{d}\varphi(1+\sigma^{2}+\eta^{4}/\kappa^{4})}{1+\eta^{4}/\kappa^{4}}-1\right],$
(73)

where $I_{2}$ is the toroidal current and $\sigma$ is a measure of up-down
asymmetry, solution to the Riccati $\sigma$ equation (see Landreman & Sengupta
(2019)), and thus not a degree of freedom. The measure $\bar{F}$ must be
evaluated at the location of an up-down symmetric cross-section. As we are
assuming stellarator symmetry, there are at least two such distinct poisitions
(in the axisymmetric limit, an infinite number of them, but $\bar{F}=0$ in
that case). This emphasises the freedom mentioned before about the description
of stellarator fields; there are two naturally simple ways of identifying the
very same field, depending on which up-down symmetric cross-section is chosen
to identify the configuration with. Depending on this choice, the meaning of
the ‘effects of changing the cross-section’ (and in particular triangularity)
will of course change, and with it the conclusions derived. One must interpret
it as the effects of changing the shape of that very cross-section, modifying
the remaining of the configuration in a consistent way, keeping the axis and
profiles fixed. For consistency and analogy with the typical axisymmetric
case, we shall choose to identify configurations with their most vertically
elongated, up-down symmetric cross-section, which often exhibit a bean-shape.

With the definitions above, the magnetic parameters $B_{20}$ and $B_{22}^{C}$
may be written explicitly in terms of the pressure gradient $p_{2}$ and the
triangularity $\delta$ (defined to be positive in the direction of the normal
curvature) of the cross-section at $\varphi=0$,

$B_{20}=-\frac{\mu_{0}p_{2}}{B_{0}^{2}}\left[1+\frac{4\sqrt{\bar{\alpha}}}{(\bar{\alpha}+3)-(\bar{\alpha}+1)\bar{F}}\left(\frac{\eta
B_{\alpha