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

This paper presents the calculation of the bounce-averaged drift of trapped particles in a near-axis framework for axisymmetric and quasisymmetric magnetic fields. 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-J property is unattainable in quasisymmetric stellarators, one may approach it through increased plasma β and triangular shaping. 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.


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 ω α .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 ω α rarely allow for analytical calculation (Velasco et al. 2023).This ends up impeding the dissection of the underlying physics and effects.
This paper carries out such calculations in a more general scenario.To make the problem tractable, we consider two main simplifications.First, we specialise to quasisymmetric (Boozer 1983; Nührenberg & Zille 1988;Rodriguez et al. 2020) and axisym- where v ∥ is the parallel velocity, the arc-length along a magnetic field-line is parametrised by ℓ, and the domain of integration is taken to be a simply connected region that satisfies v ∥ (ℓ) ⩾ 0 (which is typically referred to as a bounce-well).This is an integral at constant ψ and α, but also particle energy H and first moment µ.
It is convenient to write the bounce-average drift in Eq. (2.1) in terms of derivatives of a single scalar quantity, J ∥ (Helander 2014).This scalar quantity is the so-called second adiabatic invariant, and is an approximately conserved quantity of trapped particles.Importantly, this quantity serves as the Hamiltonian of the trapped-averaged dynamics of trapped-particles, meaning, as can be shown explicitly (Helander 2014;Calvo et al. 2017), that (2.4) Here q is the charge of the particle (which we shall take to be q = −1 for electrons), and ∆α and τ b have been defined as the total α-excursion and elapsed time in a particle bounce, respectively.Because J ∥ 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 J ∥ as a first step towards finding ω α .

Expanding the second adiabatic invariant
Let us write J ∥ explicitly as a function of the field line following coordinate ℓ (where we have taken the particle mass m = 1), (2.5) Here we have introduced the pitch-angle λ = µB 0 /H (using for the first adiabatic invariant µ = (2H − v 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 B = B/B 0 .The integral is taken between bounce points, i.e., between points at which B = 1/λ.We are assuming there is no electric field within each flux surface, and as such no electrostatic potential appears in Eq. (2.5).
It will prove convenient to express the field-line following coordinate ℓ in terms of Boozer angles (Boozer 1981).We define for that purpose the helical angle, where N is an integer that defines a helical angle, foreseeing the application to quasisymmetric devices with a helical symmetrry.Using this angular parameterisation, the Boozer-coordinate Jacobian J = B α (ψ)/B 2 , where B α = G + ιI (in the standard Boozer notation (Boozer 1981;Helander 2014)), and defining ῑ = ι − N , the second adiabatic invariant in these coordinates may now be written as (2.7) It is crucial to note that J ∥ depends directly on the magnetic field magnitude along a field line, with minimal involvement of other geometric elements.This simplifies the calculation of J ∥ ostensibly when considering quasi-and axisymmetric configurations, and makes them rather analogous.
In order to construct a representative B, we resort to the inverse-coordinate nearaxis framework (Garren & Boozer 1991b;Landreman & Sengupta 2019), in which the asymptotic form of 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 = 2ψ/B 0 , we write 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 B = B + δB(φ, χ), (2.8) where we have separated, B = 1 − rη cos χ, (2.9) and the second order, δB = r 2 B 20 + B C 22 cos 2χ + B S 22 sin 2χ . (2.10) In the ideal quasi-or axisymmetric limit, all parameters (η, B 20 , B C 22 and B S 22 ) are constants instead of functions of the toroidal angle φ, and for stellarator symmetry B S 22 = 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 |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 J ∥ calculation.† The way that we have grouped the terms in Eq. (2.8) might be unexpected, given that we have mixed asymptotically unequal terms: the constant on-axis field (r = 0) with the first order variation.However, since we are interested in trapped particles, variation in |B| along the field line is necessary, otherwise trapped particles would not exist.Therefore, in our perturbative description of the problem we must take B to constitute the leading order magnetic field magnitude.It would then appear natural to write, (2.11) † The back-of-the-envelope calculation is as follows.Simplify things by considering the tokamak notation η ∼ B0/G0 ∼ 1/R and q = 1/ι 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 ∥ dℓ should not change significantly over the distance the particle drifts in a bounce time.As v ∥ ∼ vt ϵ/R, the argument changes d(v ∥ dℓ)/dr ∼ J ∥ /r.As the transit time is ω −1 t ∼ qR/vt and v d • ∇r ∼ vtρ/R, the error ∆J ∥ /J ∥ ∼ q 2 ρ 2 R/r 3 .Thus, r ≫ rp = (q 2 ρ 2 R) 1/3 , where rp is the range where potato orbits come into play (Helander & Sigmar 2005).
where for every λ, the bounce-points χ ± b (left and right) of the integral are given by, (2.12) and attempt to expand it in powers of δB (i.e.expand around smallness of δB).There are however two important sources of conflict in doing so.First, and formally, evaluating √ 1 − λ B at the bounce points χ ± b 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 δB b (λ), defined to be the perturbed field δB(χ) evaluated, for a given particle class λ, at the bouncing points, see Eq. (2.12).We shall assume, for simplicity, that the device is stellarator symmetric about the bottom of the magnetic well (i.e.we ignore the B S 22 component), so that δB b (λ) is unique (i.e., it does not have a left and right values).Introducing this correction, let us rewrite J ∥ for a stellarator symmetric field in the following form, so that J ∥ = lim ϵ→1 − I(ϵ).Expressing the integral in this form ensures that the integrand upholds positive definiteness for all ϵ ∈ [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 ϵ, where we shall take ϵ → 1 in the final answer so that J ∥ ≈ J ∥ .If each of these terms is evaluated to the right order, we shall retrieve a consistent expression for the adiabatic invariant J ∥ .

Leading order expression
Let us start by investigating the leading order term of the expansion in ϵ.Setting the expansion parameter ϵ to zero results in the following integral, This integral is highly reminiscent of that occurring in the magnetic field of a large-aspect ratio tokamak with circular cross-section, which has been considered by many authors before in various asymptotic regimes (Connor et al. 1983;Roach et al. 1995;Helander & Sigmar 2005;Hegna 2015), and as such the derivation closely mirrors these calculations.
One may refer to Appendix A for a complete derivation.The main step required is to re-express the integral in terms of a trapping parameter k, which we define as where χ = 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 χ b = π, and thus correspond to k 2 = 1.The two trapped particle classes are connected monotonically, in the sense that λ and k maintain the order of trapped classes.With this definition, the integral may be expressed in terms of complete elliptic integrals of the first and second kind (also known as Legendre's Integrals, see e.g.Olver et al. (2020), Sec. 19), which we define as With these definitions, one can express J (0) ∥ in closed form expanding I(0) around the smallness of r.The result is, to order O(r 5/2 ), where the functions I 1 and I 2 describe the behaviour of different trapped-particle classes via k, One can readily interpret the leading form of J (0) ∥ by referring back to its basic form in terms of a parallel velocity and bounce distance, J ∥ ∼ v ∥ ℓ.The scaling with √ Hrη is directly related to the reduced parallel velocity that trapped particles must have in order for them to be trapped.† The factor B α0 /ῑ 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 α0 /ῑ 0 B 0 ∼ R 0 /(ι 0 − N ), where R 0 is the major radius of the device and N represents the helicity of the |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 J ∥ .Similarly, decreasing ῑ increases the distance between bounce points, as field lines become further misaligned with respect to |B| contours.Finally, we observe that I 1 , which describes the trapped-particle class dependence of J ∥ 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 J ∥ 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 J ∥ to that order.The expression is incomplete, as it is missing the contribution from

The first order correction
To evaluate the second order term we follow Eq.(2.14), for which we need the following integral, The second term in the integrand is a factor r smaller than the first one (which may be verified by writing expressions explicitly in terms of k), and hence, for our current expansion, only the first term needs to be taken into account.This term requires rewriting to involve the integration variable χ explicitly.Using the near-axis form of |B| for a quasiand stellarator symmetric field, Eq. (2.10), From this point the procedure is analogous to the steps followed in the leading order case; the details are given, once again, in Appendix A. We simply denote the central result here, which is the expression for J (1) ∥ expanded to leading order in r, where the function I C 22 (k) is defined to be, (2.23) Combining this result with Eq. (2.18), we may complete the asymptotic form of the second adiabatic invariant to order O(r 5/2 ), (2.24)

Trapped particle precession
With the second adiabatic invariant constructed, we are in a position to evaluate the bounce-average precession.We remind ourselves that we considered the exact quasisymmetric limit and stellarator symmetry † (e.g. a tokamak with up-down symmetry) when constructing J ∥ .Because of the idealised omnigeneous nature of the field, we need not worry about the field-line dependence (i.e., α dependence), as the behaviour is 'identical' in all field lines as far as the precession is concerned.This is apparent in Eq. (2.24).The formalism presented could however be extended to incorporate a description of said α 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. (2.4).We have almost all that is needed to compute ω α .The only remaining step is taking partial derivatives of Eq. (2.24) with respect to ψ 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 where the leading term scales like 1/r and ω α,1 ∼ O(r 0 ) (details of this derivation may be found in Appendix C).
The leading order term ω α,0 is 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, ω α,1 , where the functions G may be expressed in terms of elliptic integrals, G 20 (k) (3.4c)

Leading order precession: a tokamak-like behaviour
Let us start by analysing the leading order precession of trapped particles focusing on ω α,0 , Eq. (3.2).The expression includes physics in two ways: the overall scaling factors in front, and the k dependence through G(k), which describes the behaviour of different classes of trapped particles.We first investigate the former.
Precession is proportional to η, parameter defined in Eq. (2.9) as a measure of the leading order variation of |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 crosssections, the elongation along the horizontal axis (i.e., ratio of horizontal to vertical axes of the cross-section) is E = (η/κ) 2 , where κ is the curvature of the magnetic axis at the point where the cross-section is being assessed.Thus, for a fixed elongation, η ∼ κ.In the special case of an axisymmetric field, this means that ω α ∼ √ E/R 0 , where R 0 is the major radius.Going back to ω α,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 |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 Table 1: Some η values for QS-optimised configurations.This table presents the values of η 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 η correspond to quasi-helically symmetric configurations (namely HSX, QHS48 and precise QH), although there exist significant variations within these classes.
quasi-axisymmetric (QA) and quasi-helically symmetric (QH) stellarators (Rodriguez et al. (2022b), see some values of η in Table 1).† We finally take note of the divergent nature of ω α 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 v D • ∇θ ∼ v ∇B |∇θ|, where v ∇B denotes the gradient drift driven by the radial variation of B. As ∇θ ∼ 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. (3.2).A plot of G as a function of k is presented in Fig. 1b.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. 1a), an important feature which will become relevant for the discussion on trapped particle modes later.This behaviour only changes for the barely trapped particle classes which end up spending a significant fraction of their bouncetime near the maximum of |B|, where there is 'good curvature'.The transition point occurs at k 0 , where G(k 0 ) = 0, roughly at k 0 ∼ 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-J .This simply follows from the definition of maximum-J as the property of a field that guarantees ∂ ψ J ∥ < 0 for all trapped particles, which in terms of the electron precession is equivalent to ω α (k) < 0 for all k.Of course, this is not to say that the behaviour of a quasisymmetric field cannot become closer to maximum-J , as higher order shaping and equilibrium parameters modify the leading order behaviour above.However, one may not achieve it exactly everywhere, and especially, close to the axis.In practice, one may only get around this issue at a finite radius if the higher order contributions are strong enough.
Comparing this result against previous investigation, we find, as we already pointed, the behaviour to be identical to that shown in the work by Connor et al. (1983) for a large-aspect-ratio circular tokamak.That this correspondence is found in the more general quasisymmetric case should not come as a surprise, given the existing isomorphism between axisymmetric and quasisymmetric fields (Boozer 1983).We have, † In 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' η-behaviour.however, gained generality beyond circular-shaped cross-sections as η 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 |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.

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. (3.3), the order r correction on the leading order precession consists of three different terms: an "intrinsic" term (purely a consequence of consistency of the field with its elliptical shape), a term that is proportional to B 20 (and thus the radial increase of the average B), and another proportional to B C 22 .The behaviour of each one of these terms with k is illustrated in the left plot of Figure 2.
From these three contributions, that coming from B 20 (often called the magnetic well term) is the simplest: a positive B 20 pushes particles against the diamagnetic drift.That is, deeply trapped particles decrease their precession rate, while barely trapped ones precess even faster.This behaviour is a direct consequence of the influence of B 20 on the gradient ∇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 ∇B ∼ B × ∇(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 ∂ ψ (1/B) ∼ η 2 /2 − B 20 , this explains not only the B 20 contribution, but also the "intrinsic" 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 ω α,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 C 22 explicitly, and MHD stability is further influenced by pressure (as becomes clear when considering the Mercier criterion (Mercier 1962(Mercier , 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 C 22 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 χ dependence of |B| ∼ B C 22 cos 2χ.Naturally, both the most deeply and shallowly trapped particles, who live at χ = 0, π respectively, feel the same effect, marked by 1).In between these classes, particles perform an average of the gradient over their bounce, leading to a maximum at k ∼ 0.8.
The components B C 22 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 dp/dψ| ψ=0 )/2, as well as the triangularity of cross-sections (normalised by r) † as δ.These two parameters can substitute B C 22 and B 20 to write the precession explicitly in terms of p 2 and δ.
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 δ truly represents.For an updown 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 δ, 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 δ we are changing triangularity in the lab frame by the same amount, although δ 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, δ > 0, indicates a shift of the turning points in the direction of the curvature.‡ In the general quasisymmetric scenario, following this definition of δ, there appears to be a multitude of 'triangularities'.Each cross-section around the torus is generally different (but for an N -fold symmetry), but it is sufficient to describe the triangularity of any of its cross-sections to describe the field uniquely (given some choice of lower order parameters and a pressure gradient).† Once such a cross-section has been chosen, then one should interpret δ 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.

Relation to geometric and equilibrium parameters
Let us then proceed to write and analyse the precession of trapped particles ω α,1 in terms of triangularity, δ, and pressure gradient, p 2 .The details of how to do so are 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.
‡ One must be careful with how the sign of η is defined in this paper and in the definition of δ.For the derivation of the precession in this paper we assumed η > 0, but used the unusual convention of writing B = B0(1 − rη cos χ), with a minus sign.For the triangularity to have the meaning above, one must include a sign(η) to the definition presented in Rodríguez (2023), which is defined assuming η > 0 in the opposite convention to that considered here.
† This 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.
presented in Appendix D and rely heavily on Rodríguez (2023).The result reads, The function G p2 encodes the effect of the pressure gradient and G δ that of the triangularity, and their full explicit forms may be found in Appendix D. The function G is a complicated functions of lower order quantities that we do not present explicitly and shall ignore for the analysis in this paper.For a discussion on the form and meaning of the other contributions, we specialise now to a generally shaped, up-down symmetric tokamak configuration.
In the axisymmetric limit the functions become, using the definitions in Eqs.(3.4).Here the parameter ᾱ = (ηR 0 ) 4 = 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 ῑ0 = 2 √ ᾱG 0 I 2 /B 2 0 (1 + ᾱ), 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. (3.6a) that G axi p2 < −2 for all k and possible combinations of parameters, as G C 22 ⩾ −2, and therefore 1 + G C 22 /4 ⩾ 1/2.This negative sign of G axi p2 indicates that the usual peaked pressure profile (i.e., p 2 < 0 for the assumed sgn(ψ) = +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 β (at fixed triangularity) assists in making the behaviour of trapped particles more maximum-J .This is a wellknown 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-J is being approached by increasing β and this is generally regarded as a positive effect, at an intermediate stage particle precession reaches ω α ∼ 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 2 0 (1+ ᾱ) 2 /R 2 0 I 2 2 (ᾱ+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 β is sufficient to reverse their direction.We interpret the resulting as the critical plasma β for which the field becomes barely maximum-J .Formally, this involves equating two different asymptotic orders, which goes against the very nature of the asymptotic treatment.One may nevertheless interpret this as an estimate of the precession at a 'finite' radius.† This shows that one may try to increase the maximum-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-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-J ), ω ᾱ,0 = Hη/rB 0 from Eq. (3.2), and for the pressure we have ω Equating the two, we find For a parabolic pressure profile, a 2 p 2 = −p 0 , where a is the minor radius and p 0 the pressure on axis, one finds that the critical plasma β on axis is, We thus see that the most susceptible fields are those with a large-aspect ratio, vertical elongation (small E) and lower current (large f ).As expected, these finite β 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 ∼ 3, for which at the edge Reversing the precession of deeply trapped particles is thus predicted to require a significant plasma β.
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, ᾱ = 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 ᾱ beyond which G axi δ > 0 for all k.As G C 22 has a maximal value max(G C 22 ) ≈ 1.1, one can find that this critical point occurs at Thus, for tokamaks that are horizontally elongated (beyond some ∼ 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 ᾱ > ᾱ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-β.Most commonly, though, most tokamak fields are vertically elongated, and thus have and G p2 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 φ = π/N rather than φ = 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.
ᾱ < 1 < ᾱ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, G δ = 4ᾱ/(3+ ᾱ) > 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-J behaviour.As noted in the plasma β 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 β, 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δG δ (k)/2 ∼ G(k) for k = 0. Precisely as in the case of pressure, . (3.10) The interpretation of rδ 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δ ceases to be accurate upon approaching unity (especially as a significant bean-shape is developed).This limit of large rδ 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δ) 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. (3.10) indicates that a very significant triangularity (of order 1) is needed to significantly affect precession of deeply trapped particles in a tokamak.Hence we conclude that, although triangular shaping may help in achieving the maximum-J property together with finite β 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 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 δ represents the triangularity of the updown symmetric cross-section with the largest binormal elongation, we show in Fig. 3 the values of G δ and G p2 for a number of QS configurations.† Each segment consists of pairs (G δ , G p2 ) 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 β 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 G δ and G p2 may be provided by considering the critical β and rδ once more.As in the tokamak scenario, at the edge of the configuration β crit ∼ 2a|η|/G p2 (0) and (rδ) crit ∼ 2/G δ .A complete table for the configurations in Fig. 3 is included in Appendix E. We note here that in QA configurations, β crit ∼ 5%, † Many 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 F in Table 1 in Rodríguez (2023).Here we consistently consider its vertical cross-section (equivalently, the φ = 0 cross-section of the rotated configuration).All other configurations are similarly represented by their most vertically elongated cross-section for consistency.
while QH stellarators generally exhibit a more resilient behaviour with β crit ∼ 10 − 15%.Reversing the precession of deeply trapped particles via finite β effects thus appears to be a possibility most easily in QA configurations.Given the simple form of (rδ) 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).

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. (2.16), which we compute for the numerical precession calculation.Note that upon significant deviations from QS this choice ceases to be convenient, especially when there exist differently shaped wells within a flux surface.† This comparison is shown in Fig. 4, where the bounce-averaged precession is plotted as a function of k and for two radial locations, ϱ ψ/ψ edge .The correspondence is excellent for all displayed ϱ 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 ∼ 4), and significant shaping beyond the near-axis behaviour (see ϱ = 1).Although the trends and magnitude seem correct, there are clear deviations from the idealised limit.

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 ω α has, to leading order, a zero crossing.That is, there always exists a subgroup of trapped particles † In the more general case it is then more convenient to group precession through the class label λ, normalised as in Roach et al. (1995), namely k 2 = (1 − λBmin)Bmax/(Bmax − Bmin).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 (C 1b), which defines it in terms of λ. (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(Helander , 2020)).
Available energy (AE) 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 J ∥ .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 AE 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. (2022Mackenbach et al. ( , 2023c)).For a flux tube of length L and elliptical cross section ∆α and ∆ψ in (ψ, α)-space (principal axes), the available energy in an omnigeneous (ω ψ = 0) field may be written as, where Ĝ1/2 = (1 − λ B) −1/2 dℓ/L is the normalized bounce-time, 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 ω = ∆ψ ω/H.The integral is performed over z = H/T and λ, 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 AE 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 Hence the normalized AE becomes where the normalising factor in front will henceforth be succinctly referred to as V = 2 √ π dℓ/L B. To make further progress we realize that the integral over z is analytically tractable if one splits up ωT * as ωT * = ωT * ,0 /z + ωT * ,z , (4.4) where ωT * ,z = −∆ψ∂ ψ ln T and ωT ⋆,0 = −∆ψ(∂ ψ ln n − 3 2 ∂ ψ 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.∂ ψ ln n is negative for ψ > 0); we are specialising to density-driven trapped-particle instabilities.This assumption makes the diamagnetic drift ω ⋆ particle-energy independent, leading to, where, and the Θ function is a Heaviside function that vanishes for ω α (λ) < 0, which physically represents the inability of counter-rotating particles to drive the TPM.The function 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 | ≪ 1) and which are drifting too slow (i.e.|c 1 | ≫ 1) as compared to the drift wave have a vanishing contribution to the AE, as F → 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 ≈ 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 AE in Eq. (4.5) as an integral over k.

Leading order form of AE
Let us begin the asymptotic analysis by considering the first order expression of the trapped-particle precession ω α ≈ ω α,0 (k), defined in Eq. (3.2).To perform the integral in Eq. (4.5) we need a number of ingredients.First of all, we must consider the integral only over trapped-particle classes that co-rotate with ωT * ,0 (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 ω α,0 , i.e. k ∈ [0, k 0 ] (with the definition of k 0 |ω α,0 = 0 from before).
The second ingredient that naturally arises is then the need to express the integration variable λ in terms of k.The presence of c 1 = ωT * ,0 /ω α,0 (k) as a function of k inside 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 ∼ r.This appears to offer a way to proceed by expanding F in the small c 1 limit.However, that would be wrong, as it would neglect the most important contribution to AE. 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 ∆k ∼ r, close to k 0 , where c 1 ∼ 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 ω α 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 AE out, one last consideration is required.In this case, one needs to make an explicit assumption regarding the width of the flux tube ∆ψ, which the AE 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. (2022Mackenbach et al. ( , 2023c)), we estimate such a flattening length-scale to be the correlation length, and to be on the order of the Larmor radius ρ.As such, we write where C r may formally be dependent on other equilibrium parameters (e.g., the rotational transform ι or the flux expansion |∇ψ|).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 ϱ = r/a.This way, we define the density gradient ωn ≡ −a∂ r ln n = −∂ ϱ ln n, which scales like ϱ for a quadratic (in ϱ) density profile.In terms of these variables, the AE becomes where the common gyrokinetic expansion parameter is defined as ρ *

•
= ρ/a.The leading order expression for A w includes important information regarding the behaviour of the available energy near the axis.We highlight various scalings here: (i) One observes a strong scaling with the distance from the axis ϱ, whose origin may be presented in simple terms as follows (see the integral expression in Eq. (4.5) for reference).One factor of ϱ 1/2 may be accounted for due to trapping fraction of particles, which leaves one with an overall ϱ 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 ωn ∝ ϱ near the axis.Thus, two powers of ϱ can be argued to come from the kinetic drive of the diamagnetic drift.The final power of ϱ comes from the small fraction of trapped particles that contribute to the the available energy, as the band near k 0 scales with ϱ.
(ii) Another scaling of interest in Eq. (4.8) is its dependency on the stellarator shaping parameter η.Increased η leads to a more restricted access to energy and thus, a reduced TPM activity (as measured by AE).Thus, horizontally elongated shapes would seem to favour stability, and in the context of quasisymmetric stellarators, quasi-helically symmetric configurations.In tokamak terms, aη ∼ a √ E/R 0 ∼ ε √ E, where ε is the inverse aspect ratio.Thus, A w ∼ 1/(E 1/4 √ ε).Hence, the available energy increases with aspect-ratio keeping the minor radius fixed.This dependency becomes even stronger if one chooses ρ * ε = ρ/R 0 to be constant.
(iii) One finally observes a scaling with (C r ρ * ) 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 ρ ∼ 1/B 0 .Hence, in terms of AE, it is beneficial to have a strong magnetic field so that the length-scale over which energy is available decreases as A w ∼ 1/B 2 0 .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 AE as a measure of heat transport in practice (see Mackenbach et al. (2022)).†

A strongly-driven finite-radius asymptotic regime
In the derivation above it was key to consider the contribution to available energy from a narrow band of trapped particles with 'slow' precession.As such, the particular form of the expression in Eq. (4.8) is only valid in the regime where ω T ⋆,0 /ω α 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.‡ Although this might appear inconsistent, it is not, as any ordering assumption about ω n only affects † An 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.
‡ Note 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 |B| into the picture.However, we may formally proceed in this form.the evaluation of the AE integral, but not the particle precession itself.From the setup, it should be clear that this "strongly-driven" regime gains relevance away from the magnetic axis, where the precession of trapped electrons decreases and the diamagnetic drift increases.¶ In this new scenario, the integral for available energy may be recomputed (see Appendix F) using standard methods, as there no longer exists a narrow contributing band (see Fig. 7).All in all, one finds that the integral reduces to A different regime brings a different scaling with ϱ and ω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 AE, and thus introducing a 1/ϱω 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 ωn ϱ factor.Importantly, the scaling with η inverts compared to the weak regime, A w .Using the same tokamak-estimation for η as there, one finds A s ∼ ε 3/2 E 3/4 , suggesting that a large-aspect-ratio tokamak which is vertically elongated reduces AE.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 A w ≈ A s , denoting the Using some typical tokamak values, we estimate aη ∼ ε √ E ∼ 0.3, and thus ϱ crit ∼ 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 AE 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 devices †.We observe excellent agreement in the asymptotic behaviour between the found results and the numerical calculation.

Additional dependence of AE
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 A. We show how to do this to obtain the dependence of A on p 2 and δ to leading order at the end of Appendix F. ‡ After such considerations, we may write A ≈ A 0 + A 1 + . . ., where A 0 is the leading order expression and 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 where R 20 ≈ 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 ωn , is fixed.† To picture what is happening, we resort to the discussion on precession presented before.As we increase the pressure gradient, we learnt that the trapped-electron precession goes against the diamagnetic drift, which means that the trapped population is brought further away from resonance.But from this behaviour, one would expect the drive of the instability to decrease and with it AE 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 R 20 .The sign is, to a large extent, a result of the correction to the integral measure dk/dc 1 needed when writing the AE 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 ω α with k (see Fig. 1b).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 AE with plasma β.This is an important feature of available energy, which not only assigns value to the magnitude of ω α (ω T * ,0 − ω α ), 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, δ, in Eq. ( 4.12) depends critically on whether cross-sections are elongated vertically or horizontally, as we saw MHD stability to do in the preceding discussion.In the most common case of vertically elongated cross-sections, negative triangularity is seen to reduce the AE (which increases the precession of the trappedparticle 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. (F 24) for R δ directly to Eq. (4.2) in Rodríguez (2023) for T δ .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 β, which acts in a detrimental form on both MHD and TPM activity.
Performing a similar analysis in the strongly driven regime, we find A 1 / A 0 | strong ≈ −2.85A 1 / A 0 | weak , which presents the opposite sign to the weak regime.That is, an increased plasma β (in the form of pressure gradient) always decreases the AE, and in a standard tokamak with ᾱ < 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-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 β.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 β 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.A 1 / A 0 ∼ 1, we state that we have a critical scenario.At the edge, the critical β becomes β crit ∼ 2a|η|/R 20 (1 + f ), with f as defined before (with its QS generalisation, which may be found in Eq. (D 5a)).This shows that plasma β becomes effective in significantly changing AE in the strong regime for QAs in the regime of β crit ∼ 2 − 3%, while QH β crit ∼ 4 − 7%.Finite β 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. (3.10) may be used for AE with there.The values for QS configurations may be found in Appendix E, and range (rδ) crit ∼ 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 AE (only a very small one in the weak regime from R C 22 ).Given the changes in the weak and the strong regime, the critical transition gradient also changes and may be computed, 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 AE 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 AE 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 AE 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 AE, and the very character of the weak regime may make it numerically elusive (as the narrow AE band would have to be resolved in simulations).As such, it is natural to take the AE 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 AE, 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 AE 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 AE 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 AE, where A w ∼ ω3 n below criticality and A s ∼ ωn above it, one may postulate that the behaviour in Eq. (4.13) 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 (4.13): the system would reside on the weak-to-strong AE 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.

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-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-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(Helander , 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 λ-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 AE 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 AE 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 β).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-J is reached.The synergy between MHD and turbulent activity through triangularity inverts in the weakly driven regime, where it is under plasma β 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-J , omnigeneity (through the bounce-averaged radial drift) or AE, 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 AE).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.
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, where we remind the reader that χ b denotes the bounce angle.This shows, as did Eq.(A 1), 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, in which case the bounce points simply become ζ = ±π/2, independent of k.This transformation is well-defined because k ∈ [0, 1].The leading order contribution to the second adiabatic invariant is now equal to 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, †, we find where we have introduced, and define complete elliptic integrals of the first and second kind (Olver et al. 2020, † If 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α/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).
Sec. 19), (2.17b) Our final step is to expand 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 Collecting all the terms of order O(r) in We now turn to the next order correction in ϵ following Eq.(2.14), 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 χ.Using the near-axis expansion form of B for a quasi-and stellarator symmetric fields, we know that With these, introducing the trapping parameter, using double angle identities, and employing ζ as the integration parameter (like in the previous order), the integral reduces to where we have kept leading order terms in r.The function describing the behaviour of different trapped particle classes is I C 22 (k), which we have defined as Expanding to order r 2 and combining the first and second order result, Eq. (2.14), gives us our final expression for the second adiabatic invariant, (2.24) Figure 9 shows the comparison of this analytic expression to the numerical calculation of J ∥ .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 , where we assumed stellarator symmetry (i.e., B 20 (−φ) = B 20 (φ) 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 α), which meand that the Boozer toroidal angle φ = −(α − χ)/ῑ 0 becomes a function of both α and χ.That means that ∂ α f (φ) = −∂ φ f /ῑ 0 , which keeping the leading order near-axis term and using the same notation as in Appendix A gives, where the prime indicates derivative respect to the only argument of B 20 , φ.The integration variable is ζ, and one should interpret χ(ζ, k) = arcsin (k sin ζ).The integral may be readily evaluated using standard numerical methods.To find ω ψ , the bounce-averaged radial drift, though, we need to evaluate the leading order bounce time, τ b .Using the expression for J (1) ∥ , the bounce time is equal to Hence, the bounce-averaged radial drift to leading order can readily be found to be The averaged radial-drift ω ψ 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 ω ψ is however a function of both α 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 ω ψ over all field lines (i.e., α) 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 ϵ eff is insensitive to the sign of ω ψ .To find a single measure then of its magnitude, then, we attempt to find an upper bound for ω ψ .
Noting that the denominator in the integrand of I 20,α is an always positive function within the integration domain, with the equality only holding when the field line label α makes the largest gradient B ′ 20,max 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 .Any other value will necessarily be smaller.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 C 22 and B 20 can be re-expressed as linear combinations of pressure gradient, p 2 , and triangularity of an up-down cross-section, δ.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)dp/dψ.The notion of triangularity, δ, requires additional care.We will define δ 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 (κ, τ ) 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), 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. (D 1) has units of inverse length.This is so because δ 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, δ 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 δ lab = δ + Λ, where Λ is a term that depends only on the axis shape and first order near-axis shaping.† In the special case of an axisymmetric field, this geometric transformation term Λ vanishes.In general, though, varying δ lab or δ 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 crosssection 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, κ, and torsion, τ , it is natural to construct a measure of said asymmetry like F , introduced in Rodríguez † We note that the expression in Eq. (C2) in Rodríguez (2023) for Λ 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.
(2023), and defined to be, where I 2 is the toroidal current and σ is a measure of up-down asymmetry, solution to the Riccati σ equation (see Landreman & Sengupta (2019)), and thus not a degree of freedom.
The measure F must be evaluated at the location of an up-down symmetric crosssection.As we are assuming stellarator symmetry, there are at least two such distinct poisitions (in the axisymmetric limit, an infinite number of them, but 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 C 22 may be written explicitly in terms of the pressure gradient p 2 and the triangularity δ (defined to be positive in the direction of the normal curvature) of the cross-section at φ = 0, where ᾱ = (η/κ(0)) 4 and the dots denote terms independent of second order choices that we shall not focus on.With these expressions in place, we may then rewrite the effects of second-order on the trapped-particle precession (noting that we assumed η > 0 in the adiabatic invariant calculation in this paper), and write, These two expressions describe the effect of the pressure gradient and the cross-section triangularity on the trapped-particle precession as a function of the class of trapped particle.Note that G is independent of the pressure gradient and triangularity both.As such, when investigating dependencies on these parameters (at fixed η and axis), we may safely ignore contributions of G.The derivation may be checked with computer algebra, as provided in the repository associated to this paper.
Relation to the large-aspect ratio circular tokamak limit Here we relate the found results to the large-aspect ratio tokamak limit investigated by Connor et al. (1983) (D 6) Specialising to the circular tokamak case (i.e.ᾱ = 1 and F = 0), we find .
(D 7) In the circular limit then, we may write, acknowledging that we are mixing terms of different order in r, and not keeping all relevant terms consistently to the right order.This expression has a similar form to Eq. ( 9) in Connor et al. (1983).The first two terms are exactly identical, while the latter is different from, G α,Connor (k) = 2 3 We thus see that, though the functional form is similar, it is not the same, especially near deeply trapped particles.These differences correspond to the difference between the models considered and the meaning of changing a certain parameter thereof.In the approach of Connor, the field is being constructed in a way that locally, at finite radius, equilibrium is satisfied (as explicitly shown in Roach et al. (1995)).Such a method requires the shape of the flux surface and radial variations of various geometric and equilibrium quantities to fully specify the local equilibrium conditions.Furthermore, many of these parameters are treated independently, and one requires subsidiary assumptions in order to set these values (in particular, the 'artifice' (Connor et al. 1983) of ordering the pressure in a particular form).Though this does allow one to investigate the effect of such parameters independently, it is not guaranteed that there exists a global solution adhering to the set of local conditions.The near-axis approach presented, although asymptotic in nature, is valid in some region near the axis, and as such can perhaps be thought of a 'more global' solution.This translates into a larger degree of coupling between the geometry and the magnitude of the field (due to the singular nature of the axis) compared to the radially local approach, which ends up reducing the number of free parameters while it retains realism.
The difference is especially noticeable in the involvement of the magnetic shear, which appears explicitly in Connor et al. (1983), but becomes a higher order effect in the current treatment.In fact, we may formally obtain a sense of the involvement of the magnetic shear by considering the next order in the expansion for the precession, and focusing on the variation of the field line length due to the change in ι when taking the derivative of  3. The short names for the configurations follow straightforwardly from their full labels in Fig. 3.The starred triangularity corresponds to the triangularity for an aspect ratio 10 consideration (as there is difficulty in computing r c ).An aspect ratio of 10 was assumed to evaluate the β crit values.One should take the case of HSX sceptically, as the near-axis description is far from being quasisymmetric (a clearly indicated by ∆ QS ; in fact, it also has a very small r c , hence the small value of attainable triangularity).One could attempt modifying the near-axis configuration to make it more quasisymmetric and a better description, but we shall not do this here.
where a vanishing value indicates an exactly QS configuration to 2nd order.This measure shows that some near-axis configurations (especially that of HSX) lie far from ideal QS.
In the special case of HSX, the near-axis model is constructed to reproduce the leading harmonic content of |B|; small variations can however lead to significant effects especially at second order in the near-axis expansion (that is, where we asssess QS or compute triangularity).Because we are using these as illustrating examples, we do however not consider a further refinement of these configurations.The values in Table 2 show that physically relevant finite beta effects can have significant effect on the leading order behaviour of both deeply trapped particle precession and AE.This effect is stronger in quasi-axisymmetric configurations, and we should remind ourselves, on the outermost surface.As the magnetic axis is approached, the critical β needed will grow.Interestingly, AE is more susceptible than the precession itself.The case of triangularity shows that a significant degree of shaping is required to affect either precession or AE.This level of shaping is in many configurations present or exceeded, as the relative comparison of the critical rδ to r c δ shows.Once again, these effects become strongest as the outer surface is approached.To evaluate higher order corrections to the integral in the weak regime, we remind the reader first that the calculation involves the local expansion of the integrand.In that spirit, the changes due to second order on the first factor in the integral is straightforward, and may be simply read off Eq. (C 1b).It gives as a correction on the leading order integral where we have evaluated it at the original k 0 to leading order.A similarly simple correction arises from the corrections to the bounce time τ b , which may simply be read off Eq. (C 3), The contributions left come from the c 1 dependent terms, F(c 1 )dk/dc 1 .Here the change in the zero-crossing of precession due to second order corrections on |B| contribute to the available energy.We must thus assess where the new zero is.With the precession defined as ω α = ω α,0 + ω α,1 , Eq. (3.1), we rewrite it as With this "displaced" k ⋆ , we may assess the change in the available energy by simply considering the perturbation of the leading order k 0 dependence, namely k 0 K(k 0 )/G ′ (k 0 ), noting that the denominator comes from ω ′ α and thus it has an additional contribution.Thus, the third correction may be written as, Evaluating all these terms numerically, † the total correction due to B C 22 and B 20 is, where this factor should be understood to be a relative modification of the leading order available energy A = A 0 (1 + R).That is, the correction to available energy due to B 20 and B C 22 , which we denote as A 1 , is equal to where R 20 ≈ 1.37 and R C 22 ≈ 0.03.In the case of the strongly driven scenario, the integrals are over the whole k-space,

Figure 1 :
Figure 1: Precession of trapped particles near the magnetic axis.Left: the diagram shows the main drifts in a magnetic configuration with a dominant ∇B direction, such as in a tokamak.The electron particle ∇B drift, v ∇ , and the diamagnetic drifts, v * , are shown (the latter proportional to −B × ∇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.

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 ω α,1 .(Top right) Dependence of precession on triangularity of an axisymmetric tokamak as a function of k, for a number of values of α = E 2 , where 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 2 0 (1 + α) 2 /R 2 0 I 2 2 (α + 3).

Figure 3 :
Figure3: Effect of triangularity and pressure gradient in the precession of trapped electrons in some QS configurations.The plot shows the values of G δ and G p2 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 φ = π/N rather than φ = 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.

Figure 4 :
Figure4: A comparison between the analytical and numerical precession.The precession ω α for three QS fields and two radial positions ϱ = ψ/ψ edge is presented, normalized to H/r 2B 0 ψ edge , where a is the minor radius of the configuration, B 0 the B on axis, and 2πψ 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 α = 0), whereas the dashed red and green lines are the first and second order analytic precessions respectively.

Figure 5 :
Figure 5: Plot of the function weighing the contribution to AE of various trapped particle classes.Plot of F as a function of the diamagnetic to precession drifts c 1 .The dashed line indicates the value c 1 ≈ 3.9, which corresponds to the maximum of F, and in a sense represents the subset of particles that most vigorously resonate and drive the drift wave.
(a) A schematic of the AE distribution (b) A numerical calculation of AE distribution.

Figure 6 :
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 → 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 F[c 1 (k)] is shown to the right for rω T * ,0 /η = 0.1 as an example.Right: a numerical calculation showing the distribution of AE across the magnetic well, normalized by the total AE.Plotted for the precise QA device at a minor radial coordinate of r = 10 −3 .The points where ω α = 0 is denoted by a green dashed line, and ωn = 0.1.

Figure 7 :
Figure 7: Schematic of the contribution to available energy.Left: diagram showcasing contribution of F to the AE integrand.The broken line satifies ω α (k 0 ) = 0 to leading order.The plot of F[c 1 (k)] is shown to the right for rω T * ,0 /η ≫ 1 as an example.Right: a numerical calculation showing the distribution of AE across the magnetic well for the precise QA device at a minor radial coordinate of r = 10 −3 .The points where ω α = 0 is denoted by a green dashed line, and ωn = 10 4 .

Figure 8 :
Figure8: A comparison of the numerical calculation of AE against the analytical result.Left: comparison of the AE in the two asymptotic regimes in the precise QA configuration as a function of ϱ.The red dashed and dotted lines denote the analytic asymptotic results in the strongly and weakly driven regime, respectively.The critical radial coordinate ϱ crit is shown by a blue dash-dotted line.The crosses are numerical evaluations of the AE using the near-axis equilibrium of the precise QA configuration inLandreman & Jorge (2020).The plot has been generated with a gradient value of ωn /ϱ = 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 A, whereas the abscissa corresponds to the asymptotic result of Eq. (4.8),A w .For this plot, ωn /ϱ = 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 ϱ = r

Figure 9 :
Figure9: Residual between numerical J ∥ and analytic approximation.The plot shows, in log scale, the difference between the numerically computed J ∥ and the analytic expression, Eq. (2.24), to O(r 1/2 ) (solid line) and O(r 3/2 ) (broken line).The dotted line shows a reference ∼ 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 |B| profile with η = 1, B 20 = 1, B C 22 = 3 and k = 0.5.

Figure 10 :
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 φ, and (right) the bounce-average radial drift ω ψ as a function of α and k.The maximum bounceaverage radial drift ω ψ occurs for the deeply-trapped population at the field line where the largest gradient of B 20 aligns with the minimum of the well.
. Let us start comparing the pressure dependence of precession, and as such let us define α p = −4µ 0 rp 2 /|η|B 2

Table 2 :
Critical plasma β and triangularity in QS configurations.The table gathers critical plasma β, triangularity rδ and QS measure ∆ QS values for a number of near-axis QS configurations, configurations used in Fig.