3.1. Near-axis perspective on $\delta {\kern-1pt}B$

By definition, $\delta {\kern-1pt}B$ measures the discrepancy between the full- $B$ and the truncation of an asymptotic form. That is, $\delta {\kern-1pt}B$ is, in a sense, a measure of the truncation error of the asymptotic expansion. The truncation error of an asymptotic series (Bender & Orszag Reference Bender and Orszag2013, § 3.5) can be gauged by the next order in the expansion,Footnote ² and thus the calculation to second order could possibly be used to estimate $\delta {\kern-1pt}B$ . There is however an important caveat: in the near-axis expansion, moving from order $N$ to order $N+1$ in the expansion brings new degrees of freedom. From first to second, we may shape the triangularity of surfaces in an infinite number of ways.

Amongst this myriad of possibilities, we should pick that which most closely corresponds to our particular scenario. This was dealt with by Landreman & Sengupta (Reference Landreman and Sengupta2019) and Landreman (Reference Landreman2021), who used the first-order near-axis field of interest at $r=a$ to impose the shape of the outermost boundary of a global equilibrium, which they then asymptotically described. Formally, the description of the global equilibrium field $\tilde {\boldsymbol{B}}$ requires a double expansion in $r$ (a near-axis expansion) and $a$ (a large aspect ratio expansion to deal with the ‘global’ nature of the field). Asymptotically computing this field and evaluating its difference with $B_{\mathrm{nae}}^{(1)}(r=a)=B_0+a B_1$ can be shown to be (Landreman Reference Landreman2021, § 5.2)Footnote ³

(3.2) \begin{equation} \tilde {B}_2 = a^2\tilde {B}_0^{(2)} + r^2(\tilde {B}_{20}^{(0)} + \tilde {B}_{2c}^{(0)}\cos 2\chi + \tilde {B}_{2s}^{(0)}\sin 2\chi ). \end{equation}

The first is a mirror term, (C24) of Landreman (Reference Landreman2021), and arises from the finite aspect ratio construction (expansion in $a$ ) and is necessary to guarantee a constant toroidal flux $\psi _a=\bar {B}a^2/2$ . As such, it is a perturbation on the near-axis input $B_0(\varphi )$ , the magnetic field on axis that sets the leading trapping well structure. The $r$ -dependent terms in (3.2) correspond to the second-order near-axis expansion of $\tilde {B}$ , whose shaping is chosen to match the outermost surface. They may therefore be found in a form analogous to the regular second-order construction, solving (A41)–(A42) of Landreman & Sengupta (Reference Landreman and Sengupta2019) using (C9)–(C12) of Landreman (Reference Landreman2021) and then applying (A34)–(A36) of Landreman & Sengupta (Reference Landreman and Sengupta2019).

The field error $\delta {\kern-1pt}B$ may therefore be estimated as

(3.3) \begin{equation} \delta{\kern-1pt} B_{\mathrm{ar}}=a^2\sqrt {\frac {1}{2\pi }\int _0^{2\pi }\left [(\tilde {B}_{20}^{(0)}+\tilde {B}_0^{(2)})^2+\frac {(\tilde {B}_{2c}^{(0)})^2}{2}+\frac {(\tilde {B}_{2s}^{(0)})^2}{2}\right ]\,\mathrm{d}\varphi }, \end{equation}

where all fields have been evaluated at the boundary $r=a$ .

3.2. Numerical implementation

We have implemented the computation of $\delta{\kern-1pt}B_{\mathrm{ar}}$ numerically in pyQIC,Footnote ⁴ following the original work of Landreman (Reference Landreman2021), and introducing all the appropriate generalisations that apply to QI fields, especially the non-uniform $B_0$ and correct handling of half-helicity axes (Camacho Mata & Plunk Reference Camacho Mata and Plunk2023; Rodriguez et al. Reference Rodriguez, Plunk and Jorge2025). In practice, the resources devoted to finding $\tilde {B}$ quantities are analogous to those employed when solving the second-order equations of a near-axis equilibrium. In fact, the very equations involved are quite similar and more specific details may be found in the code itself.

To benchmark $\delta{\kern-1pt}B$ , we need a set of near-axis configurations and their associated finite aspect ratio equilibria. The latter are obtained running the global flux-surface equilibrium code VMEC (Hirshman & Whitson Reference Hirshman and Whitson1983), constructed using the near-axis surface at $r/R=0.07$ (roughly an aspect ratio of $A=14$ ). As for the benchmark set, we use a set of 1680 approximately QI, half-helicity near-axis fields (see some details in Appendix A and a fuller description in an upcoming publication). To compute $\delta{\kern-1pt}B$ as defined in (3.1), the field strength is computed in Boozer coordinates using BOOZXFORM (Sanchez et al. Reference Sanchez, Hirshman, Ware, Berry and Spong2000). The comparison of $\delta{\kern-1pt}B$ to the near-axis estimate $\delta{\kern-1pt}B_{\mathrm{ar}}$ is shown in figure 1.

Figure 1. Benchmark of the error field measure $\delta{\kern-1pt}B_{\mathrm{ar}}$ . The main plot compares the field error $\delta{\kern-1pt}B$ , (3.1), evaluated using the global equilibrium computed with VMEC, to the near-axis estimate $\delta{\kern-1pt}B_{\mathrm{ar}}$ for the QI near-axis configurations in the benchmark database of Appendix A.1. The black broken line represents $\delta{\kern-1pt}B=\delta{\kern-1pt}B_{\mathrm{nae}}$ and the dotted one is the moving average of the scatter. The colours for $\delta{\kern-1pt}B_{\mathrm{ar}}$ denote the number of field periods of the configurations (see legend). The inset plot shows the same data scaled by the field period number as $1/N^2$ .

The agreement is excellent (the correlation is in excess of 0.99), especially strong for the lowest $\delta{\kern-1pt}B$ values. At larger $\delta{\kern-1pt}B$ , there is a clear systematic deviation, with the near-axis $\delta{\kern-1pt}B_{\mathrm{ar}}$ overestimating the field discrepancy. This overestimation is a result of fields with larger $\delta{\kern-1pt}B$ having a tendency for increased shaping and sensitivity, and thus smaller radii of convergence. Evaluation at a finite aspect ratio thus tends to overestimate the magnitude of the field. However, $\delta{\kern-1pt}B_{\mathrm{ar}}$ retains a monotonic relation with $\delta{\kern-1pt}B$ .

The behaviour of the field error is clearly seen to depend on the number of field periods. They appear to form separate clusters, with the larger field period numbers leading to larger deviations (in agreement with the arguments before, with the exception of the special case $N=1$ ). Rescaling $\delta{\kern-1pt}B/N^2$ (see inset of figure 1) clusters all different fields together. This quadratic scaling follows the involvement of second derivatives in $\varphi$ in the construction of the second-order near-axis equilibrium expansion (see Rodriguez et al. Reference Rodriguez, Plunk and Jorge2025) and approximate scaling symmetries in axis shapes as discussed by Plunk et al. (Reference Plunk2025).

Altogether, the error field near-axis measure $\delta{\kern-1pt}B_{\mathrm{ar}}$ is an excellent predictor of $\delta{\kern-1pt}B$ . The systematic deviation may lead to a failure of a qualitative agreement at larger errors, but even this departure could be refined by using a fit of the curve in figure 1 to map $\delta{\kern-1pt}B_{\mathrm{ar}}$ to $\delta{\kern-1pt}B$ .

4.1. Identification of the expansion

Many of the complexities involved in evaluating $\epsilon _{\mathrm{eff}}$ , (4.1), disappear when its asymptotic near-axis form is considered. To see this, let us start by introducing the asymptotic form of the magnetic field strength, which defines trapped particle classes and determines the radial drift. Concerned with the description of a magnetic field that is approximately QI, we can write $B\approx B_0(\varphi )+rB_1(\chi ,\varphi )+r^2B_2(\chi ,\varphi )$ , where (Rodriguez et al. Reference Rodriguez, Plunk and Jorge2025)

(4.4a) \begin{align} B_1 &= -B_0(\varphi )d(\varphi )\sin \left [\alpha -\alpha _{\mathrm{buf}}(\varphi )\right ]\!, \end{align}

(4.4b) \begin{align} B_2 &= B_{20}(\varphi ) + B_{2c}(\varphi )\cos 2\chi + B_{2s}(\varphi )\sin 2\chi .\\[10pt]\nonumber \end{align}

To zeroth order in $r$ , $B_0(\varphi )$ defines (by assumption) a single toroidal trapping well per field period. The higher order contributions in $r=\sqrt {2\psi /\bar {B}}$ then deform this underlying well to make distinct ones for different field lines, labelled by $\alpha =\theta -\iota \varphi$ . Asymptotically, they do not introduce additional wells. The particular form of $B_1$ in (4.4a ) corresponds to that of an omnigeneous field (when taking $d$ to be an odd function respect to $B_0$ , and $d=0$ at the bottom and top of the trapping well), except where $\alpha _{\mathrm{buf}}\neq 0$ , i.e. the buffer regions. These are for the most part unavoidable near the edges of the toroidal domain (Cary & Shasharina Reference Cary and Shasharina1997; Plunk et al. Reference Plunk, Landreman and Helander2019; Rodríguez & Plunk Reference Rodríguez and Plunk2023). The second-order field contribution is for now chosen simply to satisfy stellarator symmetry, so $B_{20}$ and $B_{2c}$ are even, and $B_{2s}$ odd.

As a result of the trapping wells becoming labelled by $\alpha$ , the infinite sum over trapping wells in the effective ripple, (4.1), becomes by Weyl’s lemma of equidistribution (Weyl Reference Weyl1916) (Stein & Shakarchi Reference Stein and Shakarchi2011, Lemma 2.2) (assuming an irrational rotational transform),

(4.5) \begin{equation} \sum _{j=1}^{N_{\mathrm{well}}}f_j\approx \frac {N_{\mathrm{well}}}{2\pi }\int _0^{2\pi }f(\alpha )\,\mathrm{d}\alpha , \end{equation}

where $N_{\mathrm{well}}$ is the number of wells within a field line portion of length $L$ . That is, the sum over wells simply corresponds to an average over $\alpha$ .

Changing variables from $\ell$ to $\varphi$ (with the appropriate Boozer Jacobian $\mathcal{J}=(G+\iota I)/B^2$ , D’haeseleer et al. Reference D’haeseleer, Hitchon, Callen and Shohet2012, (6.6.8b)) and writing flux surface averages in terms of $(\alpha ,\varphi )$ ,

(4.6) \begin{equation} \epsilon _{\mathrm{eff}}^{3/2}=\frac {\pi }{8\sqrt {2}}\frac {(\bar {R}\bar {B})^2}{\mathcal{G}^2}\,\int _{1/B_{\mathrm{max}}}^{1/B_{\mathrm{min}}}\lambda \hat {\mathcal{E}}(\lambda )\,\mathrm{d}\lambda , \end{equation}

where

(4.7) \begin{equation} \hat {\mathcal{E}}(\lambda ) = \frac {1}{\pi }\int _0^{2\pi }\frac {\hat {H}(\lambda ,\alpha )^2}{\hat {I}(\lambda ,\alpha )}\,\mathrm{d}\alpha , \end{equation}

and all remaining definitions can be found explicitly at the beginning of Appendix B. The functions $\hat {H}$ and $\hat {I}$ (directly linked to the previous $H_j$ and $I_j$ of Nemov et al. Reference Nemov, Kasilov, Kernbichler and Heyn1999), (B1a )–(B1b ), have been defined to be dimensionless, and $\mathcal{G}^2$ (directly related to $F$ ) has dimensions of $[L]^2$ , which matches the presence of the normalising $\bar {R}$ .

In this form, the asymptotic expansion in powers of $r$ can be carried out explicitly and we do so in Appendix B. The details may be found there, and we emphasise the importance of carefully considering bounce integrals and $1/\sqrt {1-\lambda B}$ factors (as in Rodríguez & Mackenbach Reference Rodríguez and Mackenbach2023 and Rodríguez et al. Reference Rodríguez, Helander and Goodman2024). In the main text, we summarise the resulting expressions in $\epsilon _{\mathrm{eff}}^{3/2}\approx \epsilon _{\mathrm{eff}}^{3/2,(0)}+r^2\epsilon _{\mathrm{eff}}^{3/2,(2)}+O(r^4)$ , and focus on their significance.

4.2. Contribution from buffer regions: $\epsilon _{\mathrm{eff}}^{3/2,(0)}$

The leading order contribution in the distance from the magnetic axis $r$ to $\epsilon _{\mathrm{eff}}^{3/2}$ is a constant offset that sets a lower bound to the effective ripple value in the limit of $r\rightarrow 0$ . The leading order asymptotic form can be written purely in terms of first- (and zeroth-) order quantities (see the details in Appendix B),

(4.8) \begin{equation} \epsilon _{\mathrm{eff}}^{3/2,(0)}=\frac {\pi }{8\sqrt {2}}\frac {(\bar {R}\bar {B})^2}{(\mathcal{G}^{(1)})^2}\int _{1/B_0^{\mathrm{max}}}^{1/B_0^{\mathrm{min}}}\lambda \frac {(h^{(1)})^2}{I^{(0)}}\,\mathrm{d}\lambda , \end{equation}

where

(4.9a) \begin{align} h^{(1)} &=\int _{\varphi _-}^{\varphi _+}\frac {B_0}{\bar {B}}\mathcal{H}(\lambda ,B_0)d\sin \alpha _{\mathrm{buf}}\,\mathrm{d}\varphi , \end{align}

(4.9b) \begin{align} I^{(0)} &= \int _{\varphi _-}^{\varphi _+}\frac {\sqrt {1-\lambda B_0}}{(B_0/\bar {B})^2}\,\mathrm{d}\varphi , \end{align}

(4.9c) \begin{align} (\mathcal{G}^{(1)})^2 &= 2\left (\frac {1}{2\pi }\int _0^{2\pi }\mathrm{d}\alpha \int _0^{2\pi /N}\mathrm{d}\varphi \frac {\varPsi ^{(1)}}{B_0}\right )^2\Bigg /\left (\int _0^{2\pi /N}\frac {\mathrm{d}\varphi }{B_0^2}\right ), \\[10pt]\nonumber\end{align}

where

(4.10) \begin{equation} \varPsi ^{(1)}=\left [(X_{1c}\sin \chi -X_{1s}\cos \chi )^2+(Y_{1c}\sin \chi -Y_{1s}\cos \chi )^2\right ]^{1/2}. \end{equation}

The expression $\epsilon _{\mathrm{eff}}^{3/2,(0)}$ , as can be seen by inspecting $h^{(1)}$ , measures the impact that buffer regions have on omnigeneity. If such buffers were not present (which may only happen if the rotational transform matches the helicity of the axis, $\iota _0=M$ (Plunk et al. Reference Plunk, Landreman and Helander2019; Camacho Mata & Plunk Reference Camacho Mata and Plunk2023; Rodriguez et al. Reference Rodriguez, Plunk and Jorge2025), then $\alpha _{\mathrm{buf}}=0$ and $h^{(1)}=0$ , implying $\epsilon _{\mathrm{eff}}\rightarrow 0$ as $r\rightarrow 0$ . In general though, they will make a finite contribution, and the expression above allows us to gauge its magnitude. It is a measure of the level of non-omnigeneity of the first-order construction. Only if this is sufficiently small is it reasonable to proceed to the next order in the expansion (second order) to evaluate the omnigeneity at that order (Rodríguez & Plunk Reference Rodríguez and Plunk2023; Rodriguez et al. Reference Rodriguez, Plunk and Jorge2025).

Figure 2. Benchmark of effective ripple calculation. The plots show, in log scale, a comparison of the effective ripple as calculated using the code NEO on global equilibria (scatter points) and calculated with the near-axis estimate (solid lines), for a number of different benchmark configurations (see Appendix A). The two detail plots on the right show the individual comparison for two of the cases, including the zeroth-order ripple offset $\epsilon _{\mathrm{eff}}^{(0)}$ as reference (dotted line). The ripple is normalised to a reference $\bar {B}=1\,\rm T$ and $\bar {R}=1\,\rm m$ .

The expression for $h^{(1)}$ , (4.9a ), is then a local measure: trapped particles that do not venture into buffer regions near the trapping well tops, such as deeply trapped ones, will not contribute by any amount. Thus, a smaller buffer results in a smaller fraction of contributing particles. Those trapped particles that do, contribute to $h^{(1)}$ only through their radial magnetic drift (proportional to $d$ ) in the buffer region. By construction, and to abide by the condition of pseudosymmetry that avoids losses at turning points of the magnetic field, $d=0$ is chosen to vanish at both the bottom and the top of the trapping well. Hence, we expect $d\approx 0$ near the trapping well tops and thus in the buffer regions, automatically limiting the magnitude of the effective ripple. To make this somewhat more quantitative, consider the buffer to have a toroidal extent $\delta$ about the well top. For a magnetic axis with a curvature $\kappa$ that has a zero of order $v$ (fixing the form of $d=\bar {d}\kappa$ , Plunk et al. Reference Plunk, Landreman and Helander2019; Rodríguez & Plunk Reference Rodríguez and Plunk2023) and a field $B_0$ with a first derivative of order $u-1$ , following (4.9a ), we expect a scaling $h^{(1)}\propto \delta ^{u/2}\times \delta ^v \times \delta \sim \delta ^{u/2+v+1}$ . The first factor comes from the slowing down of the particle speed near the tops; the second comes from the low radial drift linked to $d$ itself; and the latter is the size of the buffer domain. The natural impulse to take $\delta \rightarrow 0$ has, though, dire consequences on the shaping of the resulting field (Plunk et al. Reference Plunk, Landreman and Helander2019; Camacho Mata et al. Reference Camacho Mata, Plunk and Jorge2022). One could alternatively increase the order of the power of $\delta$ , for instance, by making the top of the trapping well flatter. Without the need to resort to any of these, the behaviour of $I^{(0)}$ helps reduce the buffer contribution even further. This second-adiabatic-invariant-like factor is largest for barely trapped particles, reducing the contribution of the buffer to the effective ripple. This physically corresponds to these particles being more strongly diffused by collisions in $\lambda$ space, to a large extent because they travel for a longer time, and thus their drift contribution is offset to some degree. All in all, the contribution of the buffer tends to be small (see figure 2), setting a negligible lower bound on $\epsilon _{\mathrm{eff}}$ in the limit of $r\rightarrow 0$ .

4.3. Omnigeneity to second order: $\epsilon _{\mathrm{eff}}^{3/2,(2)}$

Given the rather benign contribution to the effective ripple from first order, it is important to assess the next order contribution to $\epsilon _{\mathrm{eff}}$ . The evaluation of the order $r^2$ correction to the effective ripple is presented in detail in Appendix B. The calculation is significantly more involved than first order, requiring the next two order corrections of many terms and even involving third-order quantities formally. Many of those terms do not, in practice, contribute significantly, and it is a great convenience to retain only second-order contributions.Footnote ⁵ For this reason, we write a reduced, simplified form of these contributions (see Appendix B for a more detailed account of this) that focuses on the next order contribution to $H$ ,

(4.11) \begin{equation} \epsilon _{\mathrm{eff}}^{3/2,(2)}=\frac {\pi }{8\sqrt {2}}\frac {(\bar {R}\bar {B})^2}{(\mathcal{G}^{(1)})^2}\int _{1/B_0^{\mathrm{max}}}^{1/B_0^{\mathrm{min}}}\lambda \frac {(h^{(2)})^2}{I^{(0)}}\,\mathrm{d}\lambda , \end{equation}

where

(4.12a) \begin{align} h^{(2)} &= -2\int _{\varphi _-}^{\varphi _+}\mathcal{H}(\lambda ,B_0)\frac {\varDelta B_{2c}^{\mathrm{QI}}}{\bar {B}}\,\mathrm{d}\varphi ,\\[6pt]\nonumber \end{align}

(4.12b) \begin{align} \varDelta B_{2c}^{\mathrm{QI}} &= B_{2c}^{\mathrm{QI}}-\frac {1}{4}\partial _\varphi \left (\frac {B_0^2d^2}{B_0'}\cos 2\alpha _{\mathrm{buf}}\right )\!,\\[10pt]\nonumber \end{align}

and $B_{2c}^{\mathrm{QI}}=-(B_{2c}\cos 2\bar {\iota }\varphi +B_{2s}\sin 2\bar {\iota }\varphi )$ . The expression in (4.12b ) is equivalent to the omnigeneity condition at second order when it vanishes, as derived explicitly by Rodríguez & Plunk (Reference Rodríguez and Plunk2023 (32c)). Thus, this approximation to the asymptotic $O(r^2)$ behaviour of the effective ripple is driven by the second-order non-omnigeneous behaviour and will generally dominate $\epsilon _{\mathrm{eff}}^{3/2,(0)}$ at a finite aspect ratio.

Unlike to leading order, in this case, the non-omnigeneous contribution is spread over the entirety of the trapping well, wherever $\varDelta B_{2c}^{\mathrm{QI}}\neq 0$ , a local measure of radial drift imbalance. An imbalance at some point $\varphi$ will affect all trapped particles with $\lambda \lt 1/B_0(\varphi )$ that pass over this point. Although they feel such deviation from omnigeneity, this is not to say that they will necessarily behave in a non-omnigeneous fashion, as the bounce-integral in (4.12a ) can lead to partial cancellations of $\varDelta B_{2c}^{\mathrm{QI}}$ , which can have either sign. However, it remains true that if $\varDelta B_{2c}^{\mathrm{QI}}\neq 0$ , then $\epsilon _{\mathrm{eff}}^{3/2,(2)}\neq 0$ .

The abovementioned scenario then describes the significance of the omnigeneity condition of Rodríguez & Plunk (Reference Rodríguez and Plunk2023), which is a property of second-order fields. As such, different choices of second-order shaping (namely, triangularity and Shafranov shift) will affect the omnigeneous behaviour of the field and their choice becomes key in constructing the stellarator field. This is not to say that the first-order near-axis construction is no longer important, as it will remain to strongly affect the amount and form of second-order shaping needed to make a field more omnigeneous. A clear example of this is the radial drift involved in (4.12b ), which by increasing axis curvature (and thus $d$ for a controlled elongation of flux surfaces) will generally require stronger shaping.

One may ask how can the shape at second order be chosen to omnigenise a field and if it is an example of an optimisation problem. As shown by Rodriguez et al. (Reference Rodriguez, Plunk and Jorge2025, (3.6)–(3.7)), there is however a unique, closed form choice of shaping that exactly achieves $\varDelta B_{2c}^{\mathrm{QI}}=0$ . This excludes, though, the neighbourhoods (in $\varphi$ ) of flattening points of the magnetic axis, where the shaping necessary to achieve this ideal omnigeneous behaviour would diverge unless the first-order construction was specially chosen. Excluding regions around flattening points (see further discussion later) and applying this choice of shaping, we say we have ‘omnigenised’ the construction.

4.4. Numerical implementation and benchmark

Let us first benchmark the near-axis estimate of the abovementioned effective ripple. To that end, we take an extended set of second-order near-axis equilibria previously used by Rodriguez et al. (Reference Rodriguez, Plunk and Jorge2025), described in Appendix A, and their associated global equilibria at a number of different aspect ratios. We then compare the near-axis estimate of $\epsilon _{\mathrm{eff}}$ (as a function of $r$ ), implemented in the pyQIC framework,Footnote ⁶ with the effective ripple computed by the neoclassical code NEO of the finite aspect ratio equilibria. The comparison is presented in figure 2 as a function of aspect ratio $A=R/r$ .

The agreement between the predicted near-axis $\epsilon _{\mathrm{eff}}$ and the finite volume equilibria calculation, although not exact, is excellent over quite a large range of $A$ . The buffer region contribution in all cases sets a rather benign lower bound to the value of $\epsilon _{\mathrm{eff}}$ at large aspect ratios and seems to be in agreement with the global calculation. This shows the critical role played by the second order in the expansion, which controls the dominant $O(r^{4/3})$ behaviour. At lower aspect ratios, $A\lesssim 10$ , the departures grow, as the deformation of the field increases and departs from its asymptotic description (see figure 3). The latter is unable to capture local ripplesFootnote ⁷ and the appearance of new trapped particle classes (as in figure 3 b), or misalignment of maxima (as in figure 3 c) with the associated abrupt modification of trapped particle behaviour. These effects become increasingly prominent at lower aspect ratios, and thus additional deviations between the near-axis and global calculation are to be expected.

Figure 3. Diagram with different contributions to $\epsilon _{\mathrm{eff}}$ . (a) Deformation of $|\boldsymbol{B}|$ within the principal well violating the equal radial drift condition of omnigeneity (depicted by broken line). (b) Appearance of local ripple or secondary wells. (c) Misalignment of field maxima, leading to multiple-well trapped particles.

4.5. Measures of omnigeneity

With the benchmark in place, we now introduce different measures based on the effective ripple to assess near-axis fields. We propose measures to characterise the effective ripple of second-order constructions, first-order constructions and the shaping involved when the latter are omnigenised at second order.

We start with the fundamental form of the effective ripple.

1. (a) Effective ripple from buffer regions, $\epsilon _{\mathrm{eff}}^{3/2,(0)}$ .

   Leading order contribution to the effective ripple from the buffer regions, see (4.8). Its calculation requires information of first order and the evaluation of bounce integrals between bounce points defined by $B_0(\varphi )$ .

2. (b) Second-order effective ripple, $\epsilon _{\mathrm{eff}}^{3/2,(2)}$ .

   Approximate, order $O(r^{4/3})$ contribution to the effective ripple due to omnigeneity breaking at second order, see (4.11). It requires information on the second-order construction and once again the evaluation of bounce integrals with bounce points defined by $B_0(\varphi )$ .

These combine into a single measure in item (c).

1. (c) Effective ripple at the edge, $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ .

   Estimation of the effective ripple value at the flux surface corresponding to a reference aspect ratio of $A_{\mathrm{ref}} = 10$ evaluated using the asymptotic expressions

   (4.13) \begin{equation} \epsilon _{\mathrm{eff}}^{\mathrm{edge}}=\left (\epsilon _{\mathrm{eff}}^{3/2,(0)}+r_{\mathrm{ref}}^2\epsilon _{\mathrm{eff}}^{3/2,(2)}\right )^{2/3}, \end{equation}
   with $r_{\mathrm{ref}}=R/A_{\mathrm{ref}}$ and $R$ is the major radius calculated using the magnetic axis length, $L$ , as $R=L/2\pi$ .

The $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ is a dimensionless measure that can be directly interpreted as one would the effective ripple of an optimised stellarator. This is a measure that diagnoses the effective ripple of a fully consistent second-order near-axis construction. However, as discussed previously in this section and by Rodriguez et al. (Reference Rodriguez, Plunk and Jorge2025), there are certain aspects of the field that are a result of lower-order choices in the near-axis construction. Thus, we would like to have some way to assess this ‘intrinsic’ behaviour. The simplest measure is item (d).

1. (d) Non-omnigeneous mismatch at the well bottom,Footnote ⁸ $\varDelta B_{\mathrm{min}}^{\mathrm{QI}}$ .

   Define

   (4.14) \begin{equation} \varDelta B_{\mathrm{min}}^{\mathrm{QI}}=r_{\mathrm{ref}}^2\left .\frac {\varDelta B_{2c}^{\mathrm{QI}}}{\bar {B}}\right |_{\mathrm{min}} \end{equation}
   evaluated at the bottom of the trapping well. Here, $\bar {B}$ is the average value of $B_0$ and $r_{\mathrm{ref}}$ a reference radial value (which we compute for $A_{\mathrm{ref}}=10$ ). The measure provides a sense of omnigeneity breaking near the bottom of the trapping well normalised to the average field strength. It only requires first -order quantities, following (3.8) of Rodriguez et al. (Reference Rodriguez, Plunk and Jorge2025).

To translate this relative breaking of omnigeneity into the language of the effective ripple, we can evaluate $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ for a near-axis field constructed using the omnigeneising shaping at second order. That is, the shaping that forces $\varDelta B_{2c}^{\mathrm{QI}}=0$ for all $\varphi$ , barring some masked region near the flattening points. The shaping needed to achieve this is given by Rodriguez et al. (Reference Rodriguez, Plunk and Jorge2025, (3.6)–(3.7)). The remaining ripple value is the result of the contribution from the straight sections, we denote as item (e).

1. (e) Effective ripple of omnigenised field, $\epsilon _{\mathrm{eff}}^{\mathrm{shap}}$ .

   Estimate of the effective ripple of an omnigenised near-axis field at $r_{\mathrm{ref}}$ . It is defined as $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ for a field in which the second-order shaping has been chosen to enforce $\varDelta B_{2c}^{\mathrm{QI}}=0$ everywhere except a 15 % of the toroidal domain around the inflexion points. This evaluation requires solving the second-order near-axis equations with the appropriate omnigenising shaping and computing the appropriate bounce integrals.

To quantify how much shaping is required to construct this omnigenised field, we define item (f).

1. (f) Average omnigenising shaping, $\hat {T}^{\mathrm{shap}}$ .

   Averaged measure of the amount of shaping introduced at second order. We define

   (4.15) \begin{equation} \hat {T}^2=\frac {r_{\mathrm{ref}}^2}{2\pi }\int _0^{2\pi }\left (X_{20}^2+Y_{20}^2+Z_{20}^2+\frac {X_{2c}^2+X_{2s}^2+Y_{2c}^2+Y_{2s}^2+Z_{2c}^2+Z_{2s}^2}{2}\right )\,\mathrm{d}\varphi . \end{equation}
   which we have non-dimensionalised with respect to the minor radius $r_{\mathrm{ref}}$ , indicating the relative deformation of the elliptical cross-sections needed to omnigenise a field. We define $\hat {T}^{\mathrm{shap}}$ as the shaping value of the omnigenised field at an aspect ratio $A_{\mathrm{ref}}=10=R/r_{\mathrm{ref}}$ . Note that $\hat {T} = \hat {T}[X_{2c}, X_{2s}]$ may be considered a functional of the two free shaping functions at second order.

The average measure $\hat {T}$ quantifies second-order shaping in absolute terms, but makes no reference to the first-order shape that it modifies. To elucidate how shaped a configuration is, it is important to consider the interaction with first-order shaping at finite radius $r$ . To capture this, we introduce the critical radius $r_c$ (or its reciprocal, the critical aspect ratio $A_c\sim R/r_c$ ) defined by Landreman (Reference Landreman2021). This measure captures the complexity of flux surfaces by reflecting the maximal radial extent at which flux surfaces continue to exist without any unphysical intersection, i.e. before the Jacobian of the coordinate system vanishes. That way, low $A_c$ indicates that the second-order shaping is compatible with lower first-order shaping for relatively compact fields.

1. (g) Minimal aspect ratio of omnigenised field, $A_c^{\mathrm{shap}}$ .

   Minimum value of the aspect ratio $A$ for which the near-axis construction of the omnigenised field does not present locally self-intersecting flux surfaces. This is based on the definition of Landreman (Reference Landreman2021),

   (4.16) \begin{equation} A_c=R/\mathrm{min}\left [r\,|\,\exists \,\theta ,\varphi : \mathcal{J}(r,\theta ,\varphi )=0\right ], \end{equation}
   where $\mathcal{J}$ is the Jacobian of the $\{r,\theta ,\varphi \}$ coordinate system.

We gather the abovementioned measures evaluated on the configurations in the benchmark in table 1. The first figure to look at is the effective ripple at the edge, which provides a single measure to contextualise these near-axis constructions amongst the space of optimised stellarators. Configurations with values below a per cent are often regarded as adequate with regards to neoclassical behaviour. It is clear that the second-order ( $O(r^{4/3})$ ) contribution dominates the behaviour of $\epsilon _{\mathrm{eff}}$ . To understand how much of the behaviour comes from the particular choice at second order and how much is intrinsic to the lower orders, we then look at the properties of the omnigenised field.

Table 1. Shaping configurations for effective ripple. The table shows information regarding the omnigeneous nature of the near-axis constructions in the benchmark. The table is separated into two main parts. Top rows show the measures of omnigeneity, in particular, the effective ripple for the original second order field in the benchmark (top), and the omnigenised form of the field (bottom). The lower rows present information regarding the shaping of the original field (top) and the omnigenised field. All the measures are defined in the main text.

The lack of omnigeneity near the inflexion points is manifest in $\epsilon _{\mathrm{eff}}^{\mathrm{shap}}$ , which does not approach the baseline value $\epsilon _{\mathrm{eff}}^{(0)}$ for any of the configurations. This is also indicated by the non-zero values of $\varDelta B_{\mathrm{min}}^{\mathrm{QI}}$ . However, the optimised version of config. 4.4 (by construction) and # 76, both show reduced ripple. What might appear most surprising, though, is that not all omnigenised configurations exhibit smaller ripple than their original forms, that is, $\epsilon _{\mathrm{eff}}^{\mathrm{shap}}\gt \epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ in some cases. This is a result of a strong non-omnigeneous drive near the bottom of the well, which, in the original configuration, partially cancelled the $\varDelta B_{2c}^{\mathrm{QI}}$ away from it.

Figure 4. Evolution of effective ripple with omnigenising second-order shaping. The plots show the evolution of the half-helicity benchmark configuration #76 with changing second-order shaping. (a) Evolution of $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ (solid) and $A_c$ (broken) as a function of shaping. A value of 0 for the scaled shape corresponds to the original second-order benchmark configuration, while a value of 1 indicates the omnigenised construction. The shaping is scaled linearly in between. (b) Examples of cross-sections in cylindrical coordinates at three values of shaping, indicated as vertical lines on the left plot.

Figure 4 shows in more detail what the act of omnigenising is, presenting the case of configuration # 76. We choose this case because it has an intrinsic near-omnigeneous behaviour (small $\varDelta B_{\mathrm{min}}^{\mathrm{QI}}$ ) at the bottom of its trapping well, and thus nicely illustrates the omnigenising effect of shaping. We observe that moderate variations in the cross-sections lead to significant changes in the ripple. The plot in panel (a) shows how gradually increasing the shaping to its perfect omnigenising value (0 meaning unshaped and 1 the completely shaped case) reduces the ripple, while the changes in the surfaces lead to a significant increase in $A_c$ . This reflects a general tendency observed with near-axis QI constructions: second-order shaping to improve omnigeneity is often costly in terms of shaping. Freedom at first order must also be used to effectively construct omnigeneous fields.

5.1. Construction of a ripple well measure, $A_w$

How do these secondary wells appear in the approximately QI near-axis scenario? By construction choice, in the limit of an infinite aspect ratio equilibrium ( $r\rightarrow 0$ ), there is a single well defined by $B_0$ and thus no trace of any secondary wells. The ripple well question is then a finite aspect ratio one, which will require evaluating the asymptotic description at a finite $A$ . We define $A_w$ as the largest aspect ratio $A$ for which a ripple well first appears along a magnetic field line within our asymptotic near-axis field model. Trapping wells are identified seeking points at which $\partial _\varphi |_\alpha B=0$ (i.e. partial toroidal derivative keeping the field line label constant, namely $\partial _\varphi |_\alpha B=(\partial _\varphi +\iota \partial _\theta )B$ ). To distinguish ripple wells from the main trapping well defined by $B_0(\varphi )$ on axis, we shall picture the appearance of secondary wells as a ‘dynamic’ action in radius. As we look at surfaces of decreasing aspect ratio, ripples may start forming but must do so by first forming inflexion points, $\partial _\varphi |_\alpha ^2 B=0$ . We therefore define the ripple-well aspect ratio,

(5.1) \begin{equation} A_w=R\big /\min \left \{r \,|\, \exists \,\theta ,\varphi : \partial _\varphi |_\alpha B(r,\theta ,\varphi )=0, \partial _\varphi ^2|_\alpha B(r,\theta ,\varphi )=0\right \}. \end{equation}

In practice, finding $A_w$ requires a procedure (detailed in Appendix C) similar to the calculation of the critical radius $r_c$ introduced previously and originally by Landreman (Reference Landreman2021). Briefly put, the two conditions in the definition of (5.1) written in terms of $B(r,\theta ,\varphi )=\sum _{n=0}^2 r^nB_n(\theta ,\varphi )$ may be interpreted as simultaneous algebraic equations on $r$ and $\theta$ , which may be solved for every $\varphi$ . The resulting multiple roots may then be written as a multi-valued function $\hat {A}_w(\alpha )$ (eliminating $\varphi$ ), of which the maximum is $A_w$ .

5.2. Assessment and application

Let us consider the behaviour of ripple wells in two of the example configurations in the benchmark (see Appendix A) and plot in figure 5 the multi-valued function $1/\hat {A}_w$ . For illustration purposes, we accompany these plots with the corresponding $B$ contours.

Figure 5. Ripple well diagnostic measure $\hat {A}_w$ for configurations 4.1 and 4.3. The plots show (left) the function $1/\hat {A}_w$ as a function of $\alpha$ , the field line label, for configs. (a) 4.1 and (b) 4.3, the spaghetti diagrams. The hatched area represents $\hat {A}_w\gt A_w$ , and the shaded regions the interval of $\alpha$ that have a ripple well. The right plots show the $|\boldsymbol{B}|$ contours corresponding to the $r$ values indicated by the horizontal silver lines on the left plot. The solid black line in the contour plots shows the direction of a magnetic field line, the broken white line contours of $\partial _\varphi |_\alpha B=0$ and the scatter, inflexions captured by the spaghetti diagram.

The plots of $1/\hat {A}_w$ , which we refer to as spaghetti diagrams (borrowing from the condensed matter lingo), are best interpreted from bottom to top. In the interval $0\lt 1/\hat {A}_w\lt 1/A_w$ , the field does not present any secondary well, even though it may present topological defects in the contours of $|\boldsymbol{B}|$ , i.e. puddles (Rodríguez & Plunk Reference Rodríguez and Plunk2023). At $\hat {A}_w=A_w$ , a band is reached; that is, an inflexion point appears along a field line for the first time. By symmetry, these inflexions occur in pairs. Pairs of bifurcating branches appear and separate with increasing $1/\hat {A}_w$ . The shaded areas these curves bound represent every field line (i.e. intervals of $\alpha$ ) containing secondary wells. This area grows at lower aspect ratios and additional bifurcations do appear. We must mention the unpaired bifurcations that appear near $\alpha /\pi \approx 1$ of panel (a), which appear to contravene the general behaviour presented previously, as they do not bound field lines with ripples. In fact, field lines in between have extended trapping wells. At least, until these branches cross.

With all this into account, at a finite $A\lt A_w$ , we expect more ripples in config. 4.1 as compared with 4.3. What the spaghetti diagram fails to indicate is the location (in $\varphi$ ) of such ripples and thus how strong is the radial drift experienced by ripple-trapped particles. To fully assess the relevance of ripples, some of that information should be taken into account, but we shall not explore this further here, and leave the relationship between $A_w$ and the level of departure of the near-axis approximation to $\epsilon _{\mathrm{eff}}$ (see previous section) for future investigation.

The comparison of different configurations can be made more quantitative by looking at $A_w$ , which we present in table 2. The ripple measure $A_w$ exhibits some correlation with the shaping measure $A_c$ (they have a Pearson correlation of $\rho =0.75$ , $p-\text{val}=0.05$ ). The more shaping means the more variation in the field, the larger the second-order field terms and thus the larger $A_w$ . However, this is far from being a one-to-one relationship and $A_w$ merits its own use as a field diagnostic. Minimising $A_w$ leads to ‘cleaner’ $|\boldsymbol{B}|$ contours at finite aspect ratios. From this brief benchmark, it appears that half-helicity configurations (as well as the second-order QI optimised one) are more resilient to the appearance of wells. There are other elements in the near-axis equilibrium that have the potential to significantly influence the appearance of ripple wells, one being the shape of $B_0$ . We would expect configurations with flatter minima, as those favouring MHD stability and maximum- $\mathcal{J}$ (Plunk et al. Reference Plunk, Drevlak, Rodriguez, Babin, Goodman and Hindenlang2024; Rodríguez et al. Reference Rodríguez, Helander and Goodman2024), to be more susceptible to the appearance of ripples.

Table 2. Ripple well distance in the benchmark configurations. The table shows the values of the aspect ratio when the first ripple-well appears, $A_w$ , and that at which the near-axis construction breaks down, $A_c$ , for the configurations used as a benchmark in the paper.

6.1. Sensitivity of vacuum well

To determine the sensitivity of the magnetic well to the shaping at second order, we compute the variation of the magnetic well defined in (6.1) with respect to the input functions $X_{2c}$ and $X_{2s}$ on which the integrand depends. Following the general prescription detailed in Appendix E for computing the shape gradient $\mathcal{G}^S=(\delta{\kern-1pt}S/\delta{\kern-1pt}X_{2c},\,\delta{\kern-1pt}S/\delta{\kern-1pt}X_{2s})$ of a general functional of the form

(6.4) \begin{equation} S=\int _0^{2\pi } f(B_{20},B_{2c},B_{2s},\varphi )\,\mathrm{d}\varphi , \end{equation}

we may compute the shape gradient associated with the magnetic well, $\mathcal{G}^{W}$ . For the magnetic well, we need

(6.5) \begin{equation} \frac {\partial f}{\partial B_{20}}=-\frac {2r_{\mathrm{ref}}^2}{B_0^3}\left [\int _0^{2\pi }\frac {\mathrm{d}\varphi }{B_0^2}\right ]^{-1}, \quad \frac {\partial f}{\partial B_{2c}}=0, \quad \frac {\partial f}{\partial B_{2s}}=0, \end{equation}

which may be used to compute explicitly the gradient, as in Appendix E.

Although explicit, the resulting gradient involves the inversion of a differential linear operator (see Appendix D). This is a result of $B_{20}$ being part of a self-consistent equilibrium solution. The necessary calculation to find the gradient is then computationally equivalent to solving a modified, second-order near-axis equilibrium. The strength of this gradient calculation is that, because the magnetic well is linear on $B_{20}$ , and thus also on the shaping, the gradient is independent of second order. Once we have computed the gradient, we know precisely how the magnetic well will change under changes in the shaping.

Before exploiting this knowledge of the gradient to concoct some diagnosing measures of the fields, we can gain some perspective on $\mathcal{G}^{W}$ by considering the simplest limit we can: that of an up-down symmetric tokamak field.Footnote ⁹ In this limit, any toroidal derivative drops out and (see details in Appendix E.1)

(6.6) \begin{equation} \mathcal{G}_{2c}^W = \frac {3r_{\mathrm{ref}}^2}{\pi R}\frac {\bar {d}^4-1}{\bar {d}^4-3}, \quad \mathcal{G}_{2s}^{W} =0, \end{equation}

where $\bar {d}=1$ corresponds to a tokamak with circular cross-sections.Footnote ¹⁰ For the case of circular cross-sections, the magnetic well is insensitive to changes in shaping (triangularity) (Rodríguez Reference Rodríguez2023). For a vertically elongated cross-section ( $\bar {d}\lt 1$ ), an increase in $X_{2c}$ (i.e. of positive triangularity) favours the vacuum well, matching the well-known tokamak intuition (Freidberg Reference Freidberg2014; Rodríguez Reference Rodríguez2023). The divergence as $\bar {d}^4\rightarrow 3$ is worrisome, as it leads to ill-behaviour for a very particular horizontally elongated elliptically shaped tokamak. Such behaviour was previously described and analysed by Rodríguez et al. (Reference Rodríguez2023) in the context of quasisymmetric stellarators, where the generalisation of the above observations is true. The divergence is indicative of only one particular form of shaping $X_{2c}$ being physical and consistent with how the construction is being carried out. This unphysicality can be seen on a similar divergence of triangularity, (D11). Fortunately, configurations of interest tend to live away from this singularity.

Having a knowledge of the gradient $\mathcal{G}^{W}$ provides complete insight into the sensitivity of the magnetic well of a configuration. In particular, we are interested in knowing what is the minimal amount of shaping that would make a given first-order near-axis configuration achieve a magnetic well (mirroring the practical approach of Plunk et al. (Reference Plunk, Drevlak, Rodriguez, Babin, Goodman and Hindenlang2024). Because we have computed the shape gradient, we know precisely which shaping combinations will make a configuration stable: we simply need to choose $X_{2c}$ and $X_{2s}$ such that

(6.7) \begin{equation} \mathcal{C}=\int _0^{2\pi }\left (\mathcal{G}_{2c}^{W}X_{2c}+\mathcal{G}^{W}_{2s}X_{2s}\right )\,\mathrm{d}\varphi +W_{\mathrm{ref}}=0, \end{equation}

where $W_{\mathrm{ref}}$ is the magnetic well value of the second-order near-axis construction for $X_{2c}=0=X_{2s}$ . This imposes marginal stability, and in the rare case of $W_{\mathrm{ref}}\gt 0$ , it would require relaxing the stability of the field. The ‘problem’, though, is that there is no unique way of choosing the shaping to satisfy this constraint. We must then impose some regularising choice, which we take to be one that minimises shaping in the form of an integral over the sum of squares of all the second-order shaping functions, as defined in (4.15).Footnote ¹¹ The problem of finding the minimal stabilising shaping can then be formulated as one of finding $X_{2c}$ and $X_{2s}$ that minimise the following constrained variational problem:

(6.8) \begin{equation} T^2[X_{2c},X_{2s}]=\hat {T}^2[X_{2c},X_{2s}]-\lambda \,\mathcal{C}[X_{2c},X_{2s}], \end{equation}

where $\lambda$ is a Lagrange multiplier. We denote the resulting shaping as $X_{2c}^{\mathrm{mhd}}$ and $X_{2s}^{\mathrm{mhd}}$ , for which explicit expressions may be found in Appendix F. These expressions involve the shape gradients computed previously, as well as some additional matrix multiplications, and could in principle be studied analytically. The main power of the approach is however providing a prescription to construct marginally stable approximately QI near-axis fields given a choice of axis, magnetic field strength on axis and first-order information.

To assess such stabilised constructions, we introduce some simple measures. One of the simplest is the following.

1. (b) Average second-order shaping of stabilised field, $\hat {T}^{\mathrm{mhd}}$ .

   The root-mean-square value of the consistent stabilised second-order shaping, evaluated using (4.15) using the minimal, stabilising choice of shaping, $\hat {T}^{\mathrm{mhd}} = \hat {T}[X_{2c}^{\mathrm{mhd}},X_{2s}^{\mathrm{mhd}}]$ . It is therefore a simple measure of the amount of shaping necessary to stabilise a field, normalised to some reference $r_{\mathrm{ref}}$ .

Table 3. Magnetic well sensitivity of benchmark configurations. The table shows the values of the magnetic well and re-shaped fields for the benchmark configurations, the shape $\hat {T}$ and $r_c^{\mathrm{mhd}}$ second-order shaping measures.

A more accurate measure of the shaping that takes toroidal variation into consideration may also be introduced as follows.

1. (c) Critical aspect ratio of stabilised field, $A_c^{\mathrm{mhd}}$ .

   The minimum aspect ratio for which the near-axis construction of the stabilised configuration is physical (non-intersecting flux-surfaces). That is, $A_c$ for the near-axis field evaluated with $X_{2c}^{\mathrm{mhd}}$ and $X_{2s}^{\mathrm{mhd}}$ . It is a dimensionless quantity.

Finally, to get a sense of the sensitivity of the shaping and the magnetic well, we also define $ \varDelta W/\varDelta A_c$ as follows.

1. (d) Magnetic well sensitivity, $ \varDelta W/\varDelta A_c$ .

   Change in the magnetic well $W$ with the change of shaping as measured by $A_c$ around the marginally stable point. Here, we define it by finite differencing as

   (6.9) \begin{equation} \frac {\varDelta W}{\varDelta A_c} = \frac {W(V''=1\,\text{T}^{-2}\,\text{m}^{-1})}{A_c(V''=1\,\text{T}^{-2}\,\text{m}^{-1})-A_c(V''=0)}, \end{equation}
   a dimensionless quantity. A larger value indicates increased sensitivity to shaping. We leave a more precise (perhaps analytic) form of this quantity to future work.

6.2. Numerical implementation and benchmark

We apply the abovementioned measures to the configurations that constitute the benchmark for this paper. We compute the shape gradient for all such configurations, compute their stabilised versions and evaluate the scalar measures presented previously. A summary of the relevant vacuum-well-related properties is gathered in table 3. We start by noting that all configurations in the benchmark are MHD unstable, which appears to be the overwhelming tendency with constructed QI fields, as is the case with the other classes of omnigeneous fields (Landreman & Jorge Reference Landreman and Jorge2020; Landreman Reference Landreman2022; Rodríguez et al. Reference Rodríguez2023). All examples encountered require some concerted shaping to push them towards MHD stability, which we illustrate in figure 6. The correctness of the approach is evidenced in the vanishing of the magnetic well machine precision for the reshaped configurations, $W_{{mhd}}$ , in table 3.

Figure 6. Stabilised configurations and sensitivity of vacuum magnetic well to shaping. (a,b) Plots show (left) the shaping input for stabilising the near-axis fields and (right) the resulting re-shaped cross-sections for the first three configurations of the benchmark. The black, solid lines represent the re-shaped configuration, while the dotted one denotes the starting point. (c) Change in the magnetic well as a function of the shaping measure $A_c$ , illustrating the sensitivity of the fields to shaping.

Beyond the qualitative assessment of the shaping needed from the cross-sections, both the measures $\hat {T}$ and $A_c^{\mathrm{mhd}}$ quantitatively describe the amount of shaping, with the simpler form correlating with the latter. Whenever the re-shaped configuration has a small value for the critical radius, it means that significant shaping is present and if there is a large difference from the original benchmark field construction, $A_c$ in table 3, we can conclude that the magnetic well is quite insensitive to shaping. Such sensitivity is more precisely measured by $\varDelta W/\varDelta A_c$ . Note its large value in the case of config. 4.2, figure 6(b), where minor changes in the shape of the cross-sections are enough to stabilise the field (see figure 6 c). This configuration is also special in that $B_{20,\mathrm{min}}\gt 0$ , as we may expect from the help of having a finite plasma beta in that case. The high sensitivity is a potentially problematic feature of the field if small changes in its shaping can change the response significantly. Note however that this sensitivity does change with $W$ (see figure 6 c) so that a given configuration could perhaps be made less sensitive if sufficiently stabilised. This merits further study.

If one was interested in finding the most MHD stable configurations, these measures could be used as a target of optimisation, which could be helped by the exploitation of ideas from Plunk et al. (Reference Plunk, Drevlak, Rodriguez, Babin, Goodman and Hindenlang2024) and Rodríguez et al. (Reference Rodríguez, Helander and Goodman2024). Before moving on, we note that this notion of MHD stability is only a minimal one, but not necessarily sufficient to achieve true MHD stability. Behaviour such as localised ballooning modes will react differently to the field geometry. In addition, this stabilising reshaping necessarily causes other second-order behaviours to change, and the ensuing trade-offs are an important venue of future work.

7.1. Shafranov shift shape gradient, $\mathcal{S}$

The variations problem is now one of $\delta{\kern-1pt}\boldsymbol{\varDelta }$ , a change in the shift, with $\delta{\kern-1pt}p_2$ , a change in the pressure gradient. As a result of the change in the plasma pressure, from the definition in (7.1),

(7.2) \begin{equation} \delta{\kern-1pt}\boldsymbol{\varDelta }=\underbrace {\begin{pmatrix} \delta{\kern-1pt}X_{20} \\ \delta{\kern-1pt}Y_{20} \end{pmatrix}}_{\text{Rigid shift}} + \underbrace {\begin{pmatrix} 0 \\ \delta{\kern-1pt}Y_{2c} \end{pmatrix}}_{\text{Triangular shaping}}, \end{equation}

where we distinguish between the shift of the ‘cross-section’ due to a rigid displacement and due to triangular shaping. Although both contribute to the Shafranov shift, we shall here focus on the $\theta$ -average, rigid shift, most useful when thinking of the problem in terms of shifted ellipses (see figure 9).

The variation required then simply involves the second-order equations in Appendix D, explicitly in (D3). In this case, unlike for the magnetic well sensitivity, we need the variation respect to $\delta{\kern-1pt}p_2$ , a scalar,

(7.3) \begin{equation} \hat {\mathbb{L}}\begin{pmatrix} \delta{\kern-1pt}X_{20} \\ \delta{\kern-1pt}Y_{20} \end{pmatrix} = \frac {\delta{\kern-1pt}\boldsymbol{f}}{\delta{\kern-1pt}p_2}\delta{\kern-1pt}p_2, \end{equation}

the variation of $\hat {\mathbb{L}}$ vanishing exactly. Solving this linear system, we may define

(7.4) \begin{equation} \boldsymbol{\mathcal{S}}=\hat {\mathbb{L}}^{-1}\frac {\delta{\kern-1pt}\boldsymbol{f}}{\delta{\kern-1pt}p_2}, \end{equation}

which is a vector function of $\varphi$ . This gradient holds for all second-order choices, as the equation is linear on $p_2$ . That is, $\boldsymbol{\mathcal{S}}$ is a property of the first-order field.

It is natural to introduce some normalisation for the gradient now. With $[\varDelta ]=[L]^{-1}$ and $[p_2]=[p][L]^{-2}$ , introducing the dimensionless plasma beta parameter $\beta _{p,2}=\mu _0p_2/(\bar {B}^2/2)$ for some reference magnetic field and the reference $r_{\mathrm{ref}}=R/A_{\mathrm{ref}}$ , we define $\boldsymbol{\hat{\mathcal{S}}}=(\bar{B}^{2}/2\mu_{0} r_{\mathrm{ref}})\boldsymbol{\mathcal{S}}$ . We should therefore interpret this gradient as a fractional displacement of the surface at a reference aspect ratio $A_{\mathrm{ref}}$ due to a change in plasma beta. We define the following.

1. (a) Maximum sensitivity of Shafranov shift, $\hat {\mathcal{S}}_{\mathrm{max}}$ .

   Maximum relative value of the Shafranov shift gradient with respect to changes of plasma $\beta$ at a reference aspect ratio $A_{\mathrm{ref}}$ ,

   (7.5) \begin{equation} \hat{\mathcal{S}}_{\mathrm{max}}= \max_{\varphi} \left [\frac {\bar {B}^{2}}{2\mu_{0} r_{\mathrm{ref}}}|{\boldsymbol{\mathcal{S}}}|\right ]. \end{equation}

It is illuminating to consider the familiar tokamak limit. Simplifying the system (see details in Appendix D.3), and considering for simplicity an up-down symmetric configuration,

(7.6) \begin{equation} \hat {\mathcal{S}}_{\mathrm{max}}=\frac {A_{\mathrm{ref}}}{\iota _0^2}\frac {\bar {d}^4}{\bar {d}^4-3}, \end{equation}

where all quantities have their usual meaning. The sensitivity grows with aspect ratio, as the same displacement becomes relatively larger compared with the minor radius. The scaling with $\propto \iota _0^{-2}$ responds to the physical picture of the poloidal field balancing the pressure gradient. A stronger poloidal component means the more resilient the field is to the push of the pressure. The role of elongation is more involved, although it guarantees for vertically elongated configurations, $\bar {d}\lt 1$ , $\hat {\mathcal{S}}_{x}\lt 0$ . That is, there will be a bunching of surfaces on the outboard side. It also shows the artificial divergence discussed in the previous section.

So far, the gradient measure defined in (7.5), $\hat {S}_{\mathrm{max}}$ , provides a scalar measure of sensitivity without distinguishing direction. The significance of any given shift does however depend on the direction in which it occurs. It is not the same to rigidly shift nested ellipses in the direction of the minor axis or the major axis. For directions in which surfaces are closer to each other, the bunching of surfaces is easier and so, potentially, is their intersection. We define an associated critical $\beta$ as that value for which, due to the rigid Shafranov shift, the underlying first-order ellipses just touch at a radius of $r_{\mathrm{ref}}$ . This definition directly uses $\hat {\boldsymbol{\mathcal{S}}}$ , as described in detail in Appendix G and illustrated in figure 9.

1. (b) Estimate of critical plasma $\beta$ , $\beta _\varDelta$ .

   Estimate of the critical plasma $\beta _p$ at which, due to the rigid part of the Shafranov shift, the first-order near-axis elliptical flux surfaces of an equilibrium just touch at $r=r_{\mathrm{ref}}$ . The value of the plasma beta is defined as a dimensionless scalar $\beta _\varDelta$ , for which a full definition is provided in (G12). Large values of $\beta _\varDelta$ indicate the resilience of the Shafranov shift to changes in plasma $\beta$ .

The metric $\beta _\varDelta$ condenses the information of the gradient $\boldsymbol{\mathcal{S}}$ into a single physically meaningful measure. However, it does so under certain simplifying assumptions. In particular, when it comes to describing when flux surfaces touch each other, it ignores the triangular shaping of flux surfaces, approximating them as ellipses. Thus, $\beta _\varDelta$ is only an estimate of the true critical plasma $\beta$ , which we may compute numerically.

1. (c) Numerical critical plasma $\beta$ , $\beta _c$ .

   Critical value of $\beta _p=r_{\mathrm{ref}}^2\beta _{p,2}$ , for an assumed aspect ratio $A_{\mathrm{ref}}$ , at which the flux surfaces of the near-axis construction just touch. That is, for $\beta _p = \beta _c$ , we have $A_c = A_{\mathrm{ref}}$ .

In other words, fixing $\beta _p = \beta _c$ , the near-axis construction will not succeed for aspect ratios below $A_{\mathrm{ref}}$ . The value of $\beta _c$ may be found by a root search (and thus is more numerically costly) and includes all elements of shaping up to (and including) second order. We emphasise that this quantity, although defined in exact terms, may still not be the most realistic $\beta$ limit, because it is calculated taking the axis fixed, allowing flux surfaces to move around it as the plasma beta is changed.

7.2. Implementation and examples

We assess the Shafranov shift sensitivity in the configurations of the benchmark. Figure 7 shows the shape gradient of the Shafranov shift for two different examples and table 4 summarises the scalar features of the configurations.Footnote ¹²

Figure 7. Shafranov shift sensitivity to plasma $\beta$ . The plots show the shape gradient $\hat {\mathcal{S}}_{x}$ and $\hat {\mathcal{S}}_{y}$ for configurations (a) 4.1 and (b) 4.3 in the benchmark set.

The sensitivity $\mathcal{S}_{\mathrm{max}}$ shows that there exists a wide range of sensitivities to changes in plasma $\beta$ . The gradient as a function of $\varphi$ delves deeper into these differences showing the stark difference in the binormal shift of the surfaces (which the half-helicity field 4.3 appears to avoid to a large extent). The gradients also vanish at certain points in $\varphi$ , corresponding to the directions that would break the up-down symmetry of the cross-sections at stellarator symmetric points. From this benchmark, the lowest number of field period configuration (config. 4.1) is the most sensitive, which also has the lowest value of rotational transform. This is true of $\beta _\varDelta$ as well, showing that half-helicity configurations (namely configs. 4.3, # 50 and # 76) consistently exhibit larger beta limits.

The tools introduced here make it possible, in principle, to explore the generality of such observations, but an in-depth analysis is left for future work.

Table 4. Shafranov shift sensitivity to plasma pressure. The table shows the values of the maximum sensitivity of the Shafranov shift, the estimate and numerical critical $\beta$ as well as the reference rotational transform on axis. The half-helicity fields appear to be the most resilient in the benchmark set. The comparison of $\beta _c$ to $\beta _\varDelta$ shows that although $\beta _c$ includes some key elements of the full phenomena, it can deviate significantly.