3.1.1. General statistical description

The scalar $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ exhibits an immediately evident trend across the database: stellarators with a lower number of field periods, $N$ , favour longer field scales (see figure 1). Stellarators with $N=1,\,2$ can exhibit values as low as $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}\sim 0.4$ , that is, coils with as low an aspect ratio as $A\sim 2$ . To put these values in context, we may compare them to the survey in Kappel et al. (Reference Kappel, Landreman and Malhotra2024, figure 7). The comparison to $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ values of other configurations must however be made with some words of caution. The near-axis $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ should be interpreted as a distance from the axis of the stellarator to coils, and not from the outermost flux surface of a finite aspect ratio configuration; that is, some sort of upper bound.Footnote ² With this in mind, and applying the necessary spatial scaling, the lowest $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ values in the database correspond to $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}^\star \sim 7$ , matching the best-performing configurations in that work. It is no coincidence that these best-performing configurations do also correspond to $N=1$ QI stellarators.

This improved behaviour at lower $N$ translates also to an increased number of configurations found above any reasonable realistic threshold. We take as a reference value $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}=0.25$ , which would amount to a distance of approximately 1.5 m from the plasma to the coils (in the ballpark of the needed space for various components in a fusion power plant; Beidler et al. Reference Beidler2001) for an aspect ratio 10, $R\sim 10\,\text{m}$ configuration. Setting a threshold for the analysis aligns with the notion of requiring a certain absolute physical plasma–coil distance to fit whichever elements are necessary (blanket, shielding, etc.). In this case, 96 % of the $N=1$ population exceeds the threshold, 56 % for $N=2$ and merely a per cent for $N=6$ (see the distributions in figure 1). This implies the need of larger aspect ratios for the practical realisation of larger $N$ stellarators.

The shrinking of the $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ scale with field period number is expected from a purely geometric perspective. All other things being equal, a larger number of field periods requires any existing variation of the field to be packed into an increasingly narrower sector of the stellarator. Following this simple argument, we expect a linear scaling $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}\propto 1/N$ . The dependence seen in the database is however weaker (especially at lower $N$ ; see figure 1) and although $N=1$ remains the best-performing, excellent candidates may also be found for a larger number of field periods. We show the configurations with the largest value of $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ for each field period number in figure 1(b).

Figure 1.

Behaviour of $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ across different field periods $N$ . (a) Dependence of the maximum and mean value of $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ as a function of the number of field periods $N$ in the database. Reference scalings are given as dotted lines and the distribution of the configurations as violin plots. (b) Rendering of the 3D finite aspect ratio flux surface of $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ maximising fields for each number of field periods in the database. As a result of the increased surface shaping with $N$ , the finite-volume representation of the configurations are shown for increasingly higher aspect ratios. The locations of $B_{\mathrm{min}}$ and $B_{\mathrm{max}}$ are indicated for the $N=1$ configuration, which has the largest value of $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ of all configurations (see table 6).

3.1.2. Geometric dependencies of $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$

It is interesting to ask what these large $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ configurations have in common and, in particular, which of their defining features (if any) consistently yield the observed behaviour. At a qualitative level, the configurations in figure 1(b) share a small off-the-plane excursion (relative to the $x$ - $y$ plane), yet an unevenly distributed torsion, particularly noticeable across the minimum of the magnetic field.

Figure 2.

Statistical summary of feature dependence for $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ . This figure is representative of the input feature dependence of derived quantities employed in the analysis and discussion of the database, here exemplified by $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ . (a) Scatter plots illustrating any univariate dependence of $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ on the features $\tau _{c1}$ and $\tau _{c2}$ for the $N=2$ subset. (b) Summary of key statistical measures describing the dependence of $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ on input features. (b i) Linear correlation between the input features and $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ . The different symbols indicate the Pearson and Spearman coefficients, while the colour distinguishes a different number of field periods. (bii) Nonlinear correlation coefficient, where the colour distinguishes a different number of field periods. (b iii) Friedman H-statistic where the size of the scatter indicates the significance of the measure compared with a null reference distribution and the colour represents a different numbers of field periods. The labels indicate which other feature they are most closely linked to. (b iv) Relative feature importance with symbols representing the PFI and cSAGE measures and the colour a different number of field periods. The values listed to the right next to the scatter indicate the coefficient of determination of a multivariate SVM model fitted to $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ as a function of the input features.

A more quantitative analysis of such observations requires us to be more systematic with the approach to the problem. In that direction we separate the features associated to each configuration in the database into: input features (i.e. the minimal set of parameters uniquely describing configurations) and derived features (i.e. derived properties describing the configuration, such as $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ ). The former constitutes the fundamental feature vector (following the standard machine learning nomenclature), which by virtue of the near-axis construction represents not just labels, but a precise physical description of the configurations: namely, their fundamental shape and field features. We will base our quantitative analysis on various statistical measures (detailed in Appendix B) focused on the role of these input features, which by virtue of the near-axis construction constitute a low-dimensional space.

Let us start the analysis of the influence of the input features on $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ by looking at linear correlations (see figure 2 b). The lack of a strong correlation (barring the special $N=1$ case) reflects the lack of a monotonic relation with any of the input features (see figure 2 a), but not necessarily a complete disconnection. We therefore need to supplement the analysis with an alternative that relaxes the monotonicity assumption. We construct a nonlinear correlation measure (see Appendix B for details) based on the goodness of fit of a univariate nonlinear model. This measure (figure 2 b) unveils a significant non-monotonic dependence, highlighting details of the axis shape, with pre-eminence of $\tau _{c2}$ ( $R^2\gtrsim 0.5$ across all different $N$ ).

Despite the significance of the univariate description of $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ , any such isolated model fails to capture the whole complexity of the database. We need to consider the dependence on multiple features simultaneously. When we consider all input features simultaneously, even a simple support vector regression (SVR) (Schölkopf & Smola 2002; Hastie, Tibshirani & Friedman Reference Hastie, Tibishirani and Friedman2009, Chapter 12.3.6.) model for $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ can account for more than 95 % of the database variability ( $R^2\gtrsim 0.95$ ) without falling into overfitting (see details in Appendix B). Clearly, a multivariate description is more predictive, but does not unveil how much each feature contributes to predictions. To assess such relative importance, we use permutation feature importance (PFI) and conditional Shapley additive global importance (cSAGE) (see full details in Appendix B). Both measures (figure 2 b) corroborate the dominant role of torsion and $\tau _{c2}$ , in particular, providing them high scores.

The $\tau _{c1}$ parameter does also appear to be important, but its role appears to swing wildly between different $N$ . This variability can in part be explained by the interaction of $\tau _{c1}$ with $\tau _{c2}$ in the $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ prediction. To describe this feature dependence, we compute the Friedman H-statistic (see Appendix B), which provides a measure of how much of the prediction really requires simultaneous knowledge of both features. In the case of $\tau _{c1}$ , figure 2(b) shows that approximately ${\sim}20\,\,\%$ of the behaviour requires the explicit involvement of both torsion features.

3.1.3. Understanding the role of torsion

Although it is unsurprising that torsion plays a role in setting $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ , the relative unimportance of curvature may come as a surprise, as the gradients in $\|\boldsymbol{B}\|$ are controlled by curvature, and thus, we naturally expect it to contribute to setting the length scales relevant to the field.

An explanation of this observation necessitates a more precise consideration of the contributions to $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ at a theoretical level. In that direction, consider the leading asymptotic form of $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ (Landreman Reference Landreman2021, equation (3.11)), written in the ideal QI limitFootnote ³ as

(3.2) \begin{align} {\boldsymbol{\nabla }\kern-2.1pt\boldsymbol{B}}\approx \frac {B_0}{\ell '}\left (\frac {\bar {e}'}{2\bar {e}}-\frac {B_0'}{2B_0}\right )\hat {\boldsymbol{\kappa }}\hat {\boldsymbol{\kappa }} - \frac {B_0}{\ell '}\left (\frac {\bar {e}'}{2\bar {e}}+\frac {B_0'}{2B_0}\right )\hat {\boldsymbol{\tau }}\hat {\boldsymbol{\tau }} \nonumber \\[5pt] +\,B_0\tau (\hat {\boldsymbol{\kappa }}\hat {\boldsymbol{\tau }}-\hat {\boldsymbol{\tau }}\hat {\boldsymbol{\kappa }})+B_0\kappa (\hat {\boldsymbol{t}}\hat {\boldsymbol{\kappa }}+\hat {\boldsymbol{\kappa }}\hat {\boldsymbol{t}})+\frac {B_0'}{\ell '}\hat {\boldsymbol{t}}\hat {\boldsymbol{t}}. \end{align}

There are four physical sources of length scales in the field: the variation of the field strength along the axis ( $B_0'$ ), the variation of elongation ( $\bar {e}'$ ),Footnote ⁴ the curvature ( $\kappa$ ) and the torsion ( $\tau$ ). The main curvature term is directly related to the bending ( $\hat {\boldsymbol{\kappa }}$ ) of field lines ( $\hat {\boldsymbol{t}}$ ), while the main torsion term clearly corresponds to twisting of the frame ( $\hat {\boldsymbol{\tau }}\hat {\boldsymbol{\kappa }}-\hat {\boldsymbol{\kappa }}\hat {\boldsymbol{\tau }}$ ). In practice, in the context of the database, we may focus on these two contributions to explain the behaviour of $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ . Both $B_0'$ and $\bar {e}'$ (which depends directly on $\tau$ and $\rho$ ) can be largely treated as secondary so long as their variations are tamed, which the direct control of $B_0$ and elongation in the construction enables.

Following (3.2), $\kappa$ and $\tau$ contribute on equal footing to $\|{\boldsymbol{\nabla }\kern-2.1pt\boldsymbol{B}}\|$ . The difference in their roles can be traced to the fact that the omnigenous (or rather the pseudosymmetric; see Skovoroda (Reference Skovoroda2005)) requirement on the magnetic field constrains curvature in particular, forcing it to vanish at the maxima and minima of the magnetic field (Plunk et al. Reference Plunk, Landreman and Helander2019). The field gradient at these straight sections is thus uniquely set by torsion, whose value at the bottom we denote by $\check {\tau }$ . Away from the flattening points, curvature grows and the scale of the field is set by a combination of combining curvature and torsion that peaks somewhere half-way to the maximum.

Making a direct connection between these theoretical observations and the role played by the input features is highly non-trivial. A prime example of this is the relation of the input features to $\check {\tau }$ (or any other value of torsion around the axis), which requires knowledge of $\tau _0$ after solving the axis closure problem, besides $\tau _{c1}$ and $\tau _{c2}$ . It is then natural to set up a statistical analysis of the database with an extended feature vector that includes $\tau _0$ and $\kappa _1$ . We may go beyond and complement these with a plethora of other meaningful feature combinations (see a full account of these in Appendix B), including, for instance, $\check {\tau }$ or the position of the maximum of curvature $\ell _\kappa$ . A similar dependence question on $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ could now be posed on this extended set of features. The problem is however significantly different from the original, as the features considered are no longer strongly independent. Adapting the analysis of the problem to this reality, we find it fitting to adopt the forward sequential feature selection (FSFS) (Whitney Reference Whitney1971). The importance of features in $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ is defined by sequentially choosing (and ranking) features so as to maximise the goodness of a model for $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ fitted to a sequentially increasing number of features (see the details in Appendix B). The results are summarised in table 1.

Table 1.

Dominant extended features for $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ . The table summarises the top selected features by FSFS on the extended features detailed in Appendix B. The numerical value represents the coefficient of determination $R^2$ of the $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ model using the variable on that row alongside those above. The meaning of the various parameters in the table are presented in Appendix B: $\hat {\tau },\check {\tau }$ is the torsion at $B_{\mathrm{max}}$ and $B_{\mathrm{min}}$ , $\tau _{\mathrm{min}},\tau _{\mathrm{max}}$ are the minimum and maximum torsion, respectively, $\sigma _\tau$ is the standard deviation of torsion and $\mathcal{C}=\sqrt {\kappa _{\mathrm{max}}^2+\tau _{\mathrm{\kappa _{\mathrm{max}}}}^2}$ , with $\kappa _{\mathrm{max}}$ the maximum curvature and $\tau _{\mathrm{\kappa _{\mathrm{max}}}}$ the value of torsion there.

Figure 3.

Analysis of torsion dependence of $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ . Panel (a) shows $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ in the $N=4$ subset (light grey scatter) as a function of the value of torsion at the field minimum, $\check {\tau }$ , and the combination of torsion and curvature at the point where curvature is maximum, $\mathcal{C}$ . The black broken lines represent the two different predicted upper bounds for $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ , illustrating the two limiting contributions indicated in the diagram (panel b). The behaviour of $\|{\boldsymbol{\nabla }\kern-2.1pt\boldsymbol{B}}\|$ in these two limits is exemplified for two example configurations (i) and (ii) in panel (a). Curvature and torsion are also shown for comparison.

The most relevant features appear directly related to torsion, in agreement with the reduced input feature analysis. Remarkably, the statistical analysis points at the torsion at the field minimum to be of particular importance (describing ${\sim}80\,\%$ of the data variance), connecting directly to our preceding theoretical discussion. Figure 3 illustrates $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ (for the $N=4$ case) as a function of $\check {\tau }$ , clearly showing that the value of the torsion at the field minimum sets an upper bound $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}\sim 1/\check {\tau }$ . Having an odd order of curvature at the field minimum forces a non-zero off-plane excursion of the axis, see Rodríguez et al. (Reference Giuliani, Rodríguez and Spivak2025, Appendix A) or Plunk & Rodríguez (Reference Plunk and Rodríguez2026, Appendix A), and thus, generally an indelible finite torsion in the neighbourhood of $B_{\mathrm{min}}$ , unlike near the maximum. The contribution to $\|{\boldsymbol{\nabla }\kern-2.1pt\boldsymbol{B}}\|$ is minimised there by keeping torsion spatially separated from curvature, and thus, this finite torsion through the minimum is the least detrimental form of torsion that sets a minimal scale for $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ . This limiting torsion (i.e. local twist) around $B_{\mathrm{min}}$ aligns with the configurations of figure 1, especially clear for $N=1$ .

The $N=1$ case is rather especial, in part because of its reduced range of input features (a single free parameter $\tau _0$ is scanned instead of $\{\tau _{c1},\tau _{c2},\kappa _1\}$ ). However, $\check {\tau }$ remains to capture ${\sim}99.5\,\%$ of the data variance, even if it does not appear in table 1. This highlights an important caveat of FSFS: it may hide features of similar importance just because a similar amount of information is given by an already included feature. This is partly what occurs at larger $N$ where the variation of the torsion along the axis, $\sigma _\tau$ , overtakes $\check {\tau }$ , even though the later remains highly important ( $R^2\gt 0.8$ ). This shift might nevertheless indicate an increased struggle to keeping the value of torsion tamed elsewhere, driven by $N$ .

This perspective, emphasising torsion at $B_{\mathrm{min}}$ , is however not complete, as it would suggest $\check {\tau }\rightarrow 0$ to maximise $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ , but figure 3 shows that this fails below some critical $\check {\tau }$ value. In this limit the configurations are no longer limited by the field gradients at $B_{\mathrm{min}}$ , but instead a combination of curvature and torsion somewhere between $B_{\mathrm{min}}$ and $B_{\mathrm{max}}$ (see $\|{\boldsymbol{\nabla }\kern-2.1pt\boldsymbol{B}}\|$ for configuration (ii) in figure 3). Such a combination mentioned in the theoretical discussion does also appear in the FSFS analysis in the form of $\mathcal{C}=\sqrt {\kappa _{\mathrm{max}}^2+\tau _{\mathrm{\kappa _{\mathrm{max}}}}^2}$ . This combination serves indeed as an upper bound $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}\sim 1/\mathcal{C}$ ; see figure 3. Despite this clear role, the FSFS analysis only gives $\mathcal{C}$ a secondary role. Here $\check {\tau }$ appears to capture its behaviour to a good degree; after all, torsion and curvature are tightly linked.

The best configurations will therefore attempt to minimise both these features as much as possible. First, they will try to minimise torsion, being left with some remanent through $B_{\mathrm{min}}$ . Second, they will limit the peak curvature to the largest extent possible and spatially separate it from the torsion contribution. It is then expected (and observed) for curves of maximal $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ to be rather flat, but retain a small off-the-plane excursion. Extreme elongation variations or a mirror ratio should also be avoided, but these issues are a lesser problem in practice.

3.2.1. General statistical description

As we did with $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ , consider the variation with $N$ of $A_c^{\mathrm{mhd}}$ across the database (figure 4 a). Configurations with larger $N$ are increasingly shaped, owing to the need to pack the field in an increasingly narrow toroidal domain. A naïve estimate of the scaling of the aspect ratio with $N$ based on the involvement of toroidal derivatives in the near-axis constructionFootnote ⁶ gives a pessimistic $A \sim N^2$ , while the database shows a scaling closer to, but slightly above, linear ( $A_c^{\mathrm{mhd}}\sim N^{5/4}$ ). A linear dependence on $N$ can be pictured as the natural change in aspect ratio for a stellarator put together stitching $N$ identical sectors toroidally. This stark discrepancy is not necessarily due to the imposition of MHD stability, but rather the minimal choice of second-order shaping. This appears clear when a similar analysis is performed through the database choosing the second-order completion of the configurations without targeting the magnetic well. In this case, $A_c\sim N^{3/4}$ , roughly a factor of ${\sim}2N^{1/2}$ smaller than the stability constrained solutions.

In practice, the deterioration in shaping is important, but not forbidding: configurations with reasonable aspect ratios can be found at larger $N$ (see figure 4). We take as a measure of ‘reasonable’ aspect ratio $A_c^{\mathrm{mhd}}\leqslant 10$ , which is similar to that of the Wendelstein 7-X (W7-X) modular stellarator and typical QI stellarator power plant designs (Beidler et al. Reference Beidler2001; Hegna et al. Reference Hegna2025; Lion et al. Reference Lion2025). Given the $N$ dependence of $A_c^{\mathrm{mhd}}$ , such an absolute reference becomes increasingly restrictive with $N$ , e.g. corresponding to 1 % of the database population for $N=6$ . We note, however, that there is no fundamental reason that large aspect ratios should be excluded from consideration absolutely.

Figure 4.

Behaviour of $A_c^{\mathrm{mhd}}$ across different field periods $N$ . (a) Dependence of the minimum and median value of $A_c^{\mathrm{mhd}}$ as a function of the number of field periods $N$ in the database. The reference scaling is given as a dotted line. (b) Rendition of the 3D finite aspect ratio flux surface of $A_c^{\mathrm{mhd}}$ maximising fields for each number of field periods in the database. The locations of $B_{\mathrm{min}}$ and $B_{\mathrm{max}}$ are indicated for the $N=2$ configuration, which has the lowest value of $A_c^{\mathrm{mhd}}$ .

We illustrate configurations at different $N$ showing those with the lowest $A_c^{\mathrm{mhd}}$ (for each $N$ ) in figure 4(b). The absolute lowest critical aspect ratio in the whole database is $A_c^{\mathrm{mhd}}\approx 1.6$ and is achieved by the figure-eight configuration at $N=2$ . Unlike the $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ maximising configurations in figure 1(b), configurations can present a strong off-the-plane excursion, epitomised by the figure-eight configuration (Plunk et al. Reference Plunk, Drevlak, Rodriguez, Babin, Goodman and Hindenlang2024). Such a trait does however not appear to be a necessity, a point that we shall consider in detail in the following sections. Interestingly, when the lowest aspect ratios are calculated for configurations constructed at a stronger magnetic well strength, take $W=-3\,\%$ , the shaping does increase overall and the absolute lowest $A_c^{\mathrm{mhd}}$ no longer resides on the $N=2$ figure-eight configuration, but instead shifts to $N=3$ .

Before moving on to a discussion on the different ways of achieving low $A_c^{\mathrm{mhd}}$ , we must make an aside on interpreting exceedingly low values of $A_c^{\mathrm{mhd}}$ . Any attempt to construct such configurations will quickly encounter significant limitations. Besides the deterioration of the near-axis model reliability (likely the main concern), unintended flux surface overlap could also occur. This is true even if $A\gt A_c^{\mathrm{mhd}}$ (Landreman Reference Landreman2021, § 4)), as is apparent near the crossing of the figure eight. In this particular case, the true bound on the aspect ratio becomes $A\approx 2.3$ .Footnote ⁷

Figure 5.

Average percentile measure for the lowest $A_c^{\mathrm{mhd}}$ . (a) Schematic diagram for the APM calculation. (b) The APM of the best $A_c^{\mathrm{mhd}}$ 150 configuration subset (roughly the lowest 1000th quantile for $N\geqslant 2$ ) with respect to the total population as a function of the number of field periods.

3.2.2. Dependencies of $A_c^{\mathrm{mhd}}$

Let us now delve into the dependencies of $A_c^{\mathrm{mhd}}$ on input features (see figure 23). The correlation analysis indicates that torsion is the most important feature, dominated by $\tau _{c1}$ . A significant linear correlation with $\tau _{1c}$ ( $R^2\gtrsim 0.6$ ) indicates that the critical aspect ratio of configurations tends to decrease as $\tau _{c1}$ (related to the inclination of the configurations; see figure 6 of Plunk et al. Reference Plunk, Drevlak, Rodriguez, Babin, Goodman and Hindenlang2024) increases (see figure 23 a), pushing the lowest $A_c^{\mathrm{mhd}}$ configurations towards the sampling boundary of $\tau _{c1}$ . To assess this stark difference with the behaviour for $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ , we introduce an average percentile measure (APM). The measure calculates (for each $N$ ) the database population percentile (in the feature of interest, say $\tau _{c1}$ ) to which each configuration in a small group of interest (say, the lowest 1000th quantile, 150 configurations, in $A_c^{\mathrm{mhd}}$ ) belongs and averages it (see diagram in figure 5). Applied to $\tau _{c1}$ (see figure 5) it shows a clear skewness (above 80 %), a sign that the sampling approach to the database may be missing configurations nearby and beyond this point. Although there is some truth to this (a point that a future extension of the database shall consider), the bound on the sampling space is not arbitrary, as it reflects to an intrinsic limit encountered in the construction of a family of magnetic axes with the desired half-helicity (Plunk & Rodríguez Reference Plunk and Rodríguez2026, figure 6).

Besides the involvement of torsion, flux-surface elongation also plays a non-trivial part with $A_c^{\mathrm{mhd}}$ . We can see this in both the non-trivial importance of $\rho _1$ as well as the behaviour of $\rho$ -related features in figure 5. Unlike with $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ , field shaping beyond the axis geometry has direct important practical implications. This is also evident on the FSFS analysis summarised in table 2, where a description of $A_c^{\mathrm{mhd}}$ based only on torsion is clearly lacking. Elements of elongation are thus key, as are some considerations of curvature (especially at large $N$ ).

Table 2.

Most important combined input features for $A_c^{\mathrm{mhd}}$ . The table summarises the top selected features by FSFS on the combined input features detailed in Appendix B. The numerical value represents the coefficient of determination $R^2$ of the $A_c^{\mathrm{mhd}}$ model using the variable on that row alongside those above. The meaning of the various parameters in the table are presented in Appendix B: $\hat {\tau },\check {\tau }$ is the torsion at $B_{\mathrm{max}}$ and $B_{\mathrm{min}}$ , $\tau _0$ is the integrated torsion, $\tau _{\mathrm{rms}}$ is the root-mean-square torsion, $\check {\rho }$ and $\hat {\rho }$ are the value of $\rho$ at $B_{\mathrm{min}}$ and $B_{\mathrm{max}}$ , respectively, $\hat {\rho }''$ is the variation of elongation at $B_{\mathrm{max}}$ , $\int \kappa \mathrm{d}\ell$ is the integrated curvature over half a period, $\tau _{\kappa _{\mathrm{max}}}$ is the torsion at the point of maximum curvature, $\check {\kappa }^{(3)}$ is the third derivative of curvature at $B_{\mathrm{min}}$ and $\varDelta$ is the mirror ratio.

3.2.3. The physics and dependencies of $A_c^{\mathrm{mhd}}$

To be able to frame the statistical observations made above, we must introduce some additional theoretical notions involving MHD stability and shaping of configurations.

Let us start with the notion of MHD stability, or more precisely, that of a vacuum magnetic well. When a field has a magnetic well it means that, on average, the gradient of the magnetic field magnitude points outwards, leaving the core a magnetic minimum. Formally, $W\sim \mathrm{d}(1/B^2)/\mathrm{d}\psi$ , which in the context of a near-axis field is (Landreman & Jorge Reference Landreman and Jorge2020, equation (3.6); Rodriguez & Plunk Reference Rodríguez, Plunk and Jorge2025, equation (6.1-2))

(3.4) \begin{equation} W\propto \left [\int _0^{2\pi }\frac {\mathrm{d}\varphi }{B_0^2}\right ]^{-1}\int _0^{2\pi }\frac {\mathrm{d}\varphi }{B_0^4}\left [3\left (B_{1s}^2+B_{1c}^2\right )-4B_0B_{20}\right ], \end{equation}

where quantities have their usual near-axis meaning. The first term in the integrand involves the variation of the magnetic field strength due to the curvature of the axis, $B_1$ , and its contribution to the well is always adverse (i.e. the term is always positive, pushing towards $W\gt 0$ ).Footnote ⁸ The poloidal average of the radial $|\boldsymbol{B}|$ gradient, $B_{20}$ , can in turn contribute positively to stability if $B_{20}\gt 0$ . Thus, if we take these two contributions as independent of each other, minimising curvature and making $B_{20}$ large and positive would be seen to promote stability.

However, the problem is not so simple, as $B_{20}$ is the consequence of a complex balance between geometry and the field. Although it can be strongly altered by choice of flux-surface shaping at second order in the expansion (i.e. triangularity), it also inherits part of its behaviour from lower orders. This is particularly apparent at $B_{\mathrm{min}}$ (where unlike at $B_{\mathrm{max}}$ deviations from QI are by construction minimal), which following the analysis in Rodríguez et al. (Reference Rodríguez, Helander and Goodman2024) and Plunk et al. (Reference Plunk, Drevlak, Rodriguez, Babin, Goodman and Hindenlang2024, equation (15)):

(3.5) \begin{equation} B_{20,\mathrm{b}} = \left . B_{20}\right |_{\varphi = \pi /N} = \frac {\rho }{4}\left (\frac {\text{d}\varphi }{dl}\right )^2 \left [-\frac {B_0^{\prime \prime }}{2B_0} + \frac {1}{2} \frac {\rho ^{\prime \prime }}{\rho } -\left (\tau \frac {\mathrm{d}\ell }{\mathrm{d}\varphi }\right )^2 \right ]\!. \end{equation}

Here primes denote toroidal derivatives. Although such a measure is by definition local, it can serve as guidance, as it exhibits sufficient correlation with $A_c^{\mathrm{mhd}}$ across the database and exposes clear physical contributions. The $B_0''$ term represents the detrimental inwards bending of field lines that trace the shrinking cross-sections (at constant flux) away from $B_{\mathrm{min}}$ , i.e. breathing. By virtue of belonging to the $(2,3)$ class, this local contribution vanishes. The same cannot be said of the $\rho ''$ term, which favours configurations whose elongation grows away from the minimum, i.e. stretching. This is reflected in the database, as both the APM of $\check {\rho }''$ and $\check {\rho }$ in figure 5 and the FSFS analysis (table 2) show.Footnote ⁹ Finally, we have the always detrimental contribution of torsion, representing the inwards pointing curvature of any helical field line. This universality in its contribution presents minimising torsion everywhere as a solid principle to ameliorate stability. This is reflected in the involvement of $\tau _0$ and $\tau _{\mathrm{rms}}$ in the FSFS analysis, as well as APM in figure 5.

The APM analysis is particularly insightful here, showing a consistent low integrated torsion $\tau _{0}$ throughout the low $A_c^{\mathrm{mhd}}$ configurations across the database (below the 20 % mark for $N\geqslant 2$ ). This aligns with the configurations living near the $\tau _{c1}$ domain boundary, as noted above. In this limit the average torsion of the constructed curves sharply drops, as they are able to release their ‘tension’ through large inclination (see figure 3 in Plunk & Rodríguez Reference Plunk and Rodríguez2026). Despite these low values, it must be noted that $\tau _0\approx 0$ curves are seldom seen in the database, likely impeded by the sharpness of the transition and the granularity of the sampling. Curves with small total twist sharply increase their curvature to remain closed, which leads to curvature signs in table 2 and figure 5, and could also play a regularising role through (3.4). Although this integrated torsion picture appears to nicely tie the database and theory, we must note that such curves do not have a pointwise vanishing torsion (with the exception of the figure-eight configuration). In that sense then, is the torsion in the field really being minimised? Is minimising $\tau _{\mathrm{rms}}$ not closer to the theoretical argument? Apparently not: minimising integrated torsion appears a priority over $\tau _{\mathrm{rms}}$ , especially at large $N$ .

We must discuss shaping to explain this observation and, in particular, the sensitivity of the ellipse rotation to the mean torsion. As it follows from the so-called $\sigma$ equation (Plunk & Rodríguez Reference Plunk and Rodríguez2026, equation (B3); Landreman & Sengupta Reference Landreman and Sengupta2019, equation (A26)) in the vanishing integrated torsion scenario a solution to the equation with a small $\sigma$ is possible, minimising ‘shaping proliferation’. Therefore, $\tau _0$ and $\tau _{\mathrm{rms}}$ are not the same and their difference in value grows with $N$ for $N\geqslant 2$ , explaining the trends in figure 5 and hinting an explanation of the high- $N$ burden that stability imposes on shaping requirements ( $A_c^{\mathrm{mhd}}/A_c\sim N^{1/2}$ ).

In summary, a few strategies for reducing $A_c^{\mathrm{mhd}}$ are apparent: (i) reduce torsion as much as possible, (ii) limit the amount of shaping in the construction by minimising $\tau _0$ (and if possible stick to lower $N$ ), (iii) avoid extreme curvature values, and (iv) leverage elongation stretching.

3.2.4. Different forms of field

The broad criteria for achieving low $A_c^{\mathrm{mhd}}$ is fitting to understand the statistical behaviour of the database. Such principles are however not rigid, especially given the complexity of the measures involved. This opens the door to a larger variety of configurations compared with the $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ scenario, as was appreciated in figure 4.

Figure 6.

Clustering of the least shaped $N=2$ configurations. The figure shows the clustering of the 1000th percentile (a total of 145 configurations) of the least shaped $N=2$ configurations represented as a cloud of points in two dimensions. Here MDS was used to perform the dimensionality reduction. Examples of the two resulting clusters, matching different $\tau _{c1}$ values, are illustrated. The figure-eight-like configurations are represented by crosses (orange) and flat configurations by circles (blue). The distribution of some of their properties are shown in the histograms. The discreteness of the database precludes a definitive consideration about the inherent nature of these clusters.

To investigate this, let us look for different configurations within the group of the smallest $A_c^{\mathrm{mhd}}$ , $N=2$ configurations (the same 1000th quantile as for the APM).Footnote ¹⁰ To distinguish ‘significantly’ different configurations, we proceed systematically by applying clustering algorithms (see details on the methodology in Appendix C) and representing the input feature space in an approximately distance-preserving lower-dimensional projection using multidimensional scaling (MDS) (Borg & Groenen Reference Borg and Groenen2005); see figure 6. We observe configurations to cluster into two groups, each tagged by a different value of $\tau _{c1}$ . The group with $\tau _{c1}\approx 0$ corresponds to the markedly inclined figure-eight family; the other contains flatter configurations. This separation is likely artificial, a result of the granularity of the database. It is nevertheless clear that there is a degree of flexibility in achieving low $A_c^{\mathrm{mhd}}$ with crowns and flat configurations as limits. The figure-eight family tightly conforms to the guiding principles presented above: smallest torsion and minimal shaping (see the bottom of figure 6). The more open configurations are able to maintain a similar $A_c^{\mathrm{mhd}}$ despite requiring a moderate increase in shaping (as captured by $\hat {T}_{\mathrm{mhd}}$ )Footnote ¹¹ due to having larger torsion, although in exchange for a smaller curvature.

Figure 7.

Examples of the least shaped marginally stable configurations of the ‘flat’ and ‘crown’ classes. The figures are a 3D rendition of the representative least shaped, marginally stable configurations continuing the two $N=2$ families in figure 6. The table shows the $A_c^{\mathrm{mhd}}$ values for each of these (including the $N=2$ as reference).

Figure 8.

Family of $N=1$ configurations. Examples of configurations belonging to the low $A_c^{\mathrm{mhd}}$ , $N=1$ family of configurations. A continuum of stellarators parametrised by the integrated torsion appears to exist, their corresponding critical aspect ratios are, in order, $A_c^{\mathrm{mhd}}=2.1,\,2.0,\,2.3,\,2.3$ .

A similar analysis and conclusion applies for larger $N$ , which we represent in figure 7. We have separated configurations between flat configurations and crowns. The latter can be thought of as extensions of the figure-eight configuration to higher $N$ , showing a degree of inclination, as the $N=3$ example nicely illustrates. Note that only milder versions of such crowns are found at higher $N$ , perhaps inhibited by the increase in the curvature, or some other theoretical limit in the construction of the axis curves, although a finer resolution in the database would presumably help assess this question. Despite this trend, both ways of achieving low $A_c^{\mathrm{mhd}}$ (small and modest mean torsion) remain accessible with increasing $N$ .

The keen reader will have noticed that we have deliberately left the $N=1$ configurations out of this discussion. These configurations clearly belong to a wide continuum of configurations, as explicitly shown in figure 8. This family is naturally parametrised by their integrated torsion, which is actually the sampling parameter used to construct them as part of the database. The latter circumvents some of the limitations for the other $N$ inherited from the discreteness and sensitivity on $\tau _{1c}$ , unveiling a continuous family more clearly. The lowest $A_c^{\mathrm{mhd}}$ is still found under conditions of near-vanishing integrated torsion (see figure 8), as expected from the above discussion.

3.3.1. Two approaches to large fractions

One natural approach to achieving large $f_{\mathcal{J}}$ would be indeed to seek a dominant and positive $B_{20}$ value. The field equilibrium can be pushed in this direction leveraging large second-order shaping.Footnote ¹² The resulting fields are however hardly realistic, as extremely high aspect ratios would be generally required for their realisation. It is also true that the shaping in the database has not been chosen with this particular strategy in mind, but rather to achieve a desired value of $W$ . Unlike $f_{\mathcal{J}}$ , the magnetic well measure only cares about average behaviour and not pointwise values. Thus, even when strongly shaped, we cannot expect to achieve a fully maximum- $\mathcal{J}$ field.

Even if second-order shaping was chosen more fittingly, we know there is an intrinsic part of $B_{20}$ , near flattening points (Rodríguez et al. Reference Rodríguez, Helander and Goodman2024, Reference Rodríguez, Plunk and Jorge2025), that cannot be altered by it. In practice, the lack of correlation between $B_{20,b}$ and $f_{\mathcal{J}}$ implies that this limitation has no impact. This irrelevance derives from the flatness of the field about $B_{\mathrm{min}}$ , which dilutes the contribution of $B_{\mathrm{min}}$ with its shaped surroundings. The two examples of figure 9(a) illustrate this, showing a very sudden transition in the particle precession $\omega _\alpha$ for extreme values of $\lambda$ corresponding to trapped particles, $\lambda =1/B_{\mathrm{min}}$ . Most of the population, corresponding to smaller values of $\lambda$ , is therefore not affected.

The second contribution inside the square bracket in (3.8) is

(3.9) \begin{equation} T=\left (B_0^2d^2/B_0'\right )'/4. \end{equation}

The origin of this term is understood as follows. In a QI field, the leading (first-order) poloidal magnetic drift does not contribute to precession due to its parity, forced by omnigeneity. However, if one takes into account the change in the parallel velocity of the particles due to the magnetic field variation in the neighbourhood of the axis, i.e. $B_1$ contribution, a non-zero, second-order net poloidal precession due to first-order variations is obtained. The result is $T$ .Footnote ¹³ Stellarator symmetry forces this term to be even in $\varphi$ , necessarily vanishing at both $B_{\mathrm{min}}$ and $B_{\mathrm{max}}$ ,Footnote ¹⁴ and satisfying $T\geqslant 0$ ( $T\leqslant 0$ ) in their neighbourhood, respectively. The sign of $T$ changes roughly at the point of maximum $d$ (also the point of maximum curvature when we control elongation), making a majority of trapped particles contribute positively to $f_{\mathcal{J}}$ (as bounce times are biased towards barely trapped particles). Although this approach cannot attain $f_{\mathcal{J}}=1$ on its own, maximal profit is achieved if the bottom of the well is as wide as possible and the curvature maximum as close as possible to it. That way, the $T\geqslant 0$ contribution is restricted to a smaller fraction of deeply trapped particles and the positive contribution maximised. For the first time in the context of this paper, curvature does not necessarily act to the detriment of the desired behaviour.

To investigate how these two extreme approaches materialise in the database, we investigate the subset of configurations with the largest values of $f_{\mathcal{J}}$ for each field period. As shown explicitly in figure 9(b), the distributions of these configurations exhibit clear traces of both the large and small shaping approaches, especially for $N=2,\,3$ . The $N$ scaling of the critical aspect ratio is indicative of this essential difference, with a scaling $A_c^{\mathrm{mhd}}\sim N^2$ matching strong shaping and $A_c^{\mathrm{mhd}}\sim N$ the lower one. At larger $N$ it appears that the distribution of configurations with high $f_{\mathcal{J}}$ is biased towards the larger shaping configurations, as the natural magnitude of second-order quantities grow in this limit.

Figure 9.

High and low shaping approaches to $f_{\mathcal{J}}$ . (a) Near-axis detail of the high (top) and low (bottom) shaping approaches to $f_{\mathcal{J}}$ . Configurations are $N=2$ corresponding to the red scatter from the shaped cluster and (ii) in panel (c). The left plots show the contributions to the integrand of $\omega _{\alpha ,\mathrm{vac}}$ in (3.8), showing the $B_{20}$ and $T$ -dominated limits. The plots on the right show the particle precession at $r=0.1$ as a function of $\lambda$ , with the grey lines representing the precession for different $\alpha$ field lines. (b) Behaviour of the distribution of $A_c^{\mathrm{mhd}}$ across different $N$ , both of the whole database (left half of the violin plots) and the top $f_{\mathcal{J}}$ 1000th percentile subset (right half). The presence of separate clusters is apparent especially for $N=2{-}4$ and is highlighted by representing the mean and standard deviation of fitted Gaussian mixture models through error bars (Wierzchoń & Kłopotek Reference Wierzchoń and Kłopotek2018, Chapter 3.2). The line width represents the weight of the Gaussian. Reference scalings are given, following each of the cluster evolutions. (c) Illustration of the $N=2$ , top $f_{\mathcal{J}}$ 1000th percentile in a MDS representation (left). The clustering distinguishes the high and low shaping scenarios, and three different configurations in the latter are shown labelled (i), (ii) and (iii). The distribution of some features distinguishing the two clusters are shown as histograms (the orange hashed bars corresponding to the moderately shaped fields).

The presence of these two groups can also be observed by application of clustering algorithms (see Appendix C), which naturally split the population according to this criterion (see figure 9 c). Some representative shapes of the low shaping (i.e. $T$ -dominated) cluster are given, recovering shapes seen in the analysis of $A_c^{\mathrm{mhd}}$ : configurations with a limited torsion that limit the amount of shaping and promote $B_{20}\gt 0$ . In the $T$ -dominated regime the value of $B_{20}$ can still play a relevant role, in particular near points where $T=0$ . More generally, a compromise between $B_{20}$ and $T$ contributions may be needed, especially if close-to-unity values of $f_{\mathcal{J}}$ are sought; the contribution from $B_{20}$ must amend the mishaps near $B_{\mathrm{min}}$ through $T$ .

Unlike in the $A_c^{\mathrm{mhd}}$ analysis, here curvature is not an impediment, allowing highly exotic configurations to come to the forefront, with especially large inclinations (and even knotted forms), such as those illustrated in figure 10

Figure 10.

Large curvature fields in the large $f_{\mathcal{J}}$ subsets. The figure presents three examples of large curvature fields for (a) $N=3$ , (b) $N=4$ and (c) $N=5$ , as part of the low shaping clusters. These show a knotted configuration, a crown and a curled alternative. Both (a) and (c) are more a curiosity than practical alternatives, but emphasise the breadth of the database, as well as the large curvature solutions. All these configurations remain half-helicity stellarators, following the generalised definition in Appendix E.

3.3.2. Statistical description of $f_{\mathcal{J}}$

Armed with the perspective above, we now proceed to the statistical analysis of $f_{\mathcal{J}}$ across the database (see table 6). The database does not present a clear dependence of $f_{\mathcal{J}}$ on $N$ , similar best values $f_{\mathcal{J}}\gtrsim 80$ % being achievable for any $N$ . This seeming immutability is in part due to $f_{\mathcal{J}}$ not depending on the magnitude of the precession of particles, but rather their directions. That is, even if precession may scale with $N$ , the fraction does not necessarily do so. Having said that, at larger $N$ we expect the role of $B_{20}$ to become more prominent in the fight with $T$ in (3.8), given its natural ${\sim}N^2$ scaling. We see this to be indeed the case in the two-distribution analysis of figure 9(b).

The correlation analysis (see figure 24) reflects the complexity of $f_{\mathcal{J}}$ , indicating a distributed importance amongst elements of all curvature, torsion, magnetic field and elongation. It is particularly remarkable how elements of the magnetic field shape, $\lambda _B$ , and curvature, $\kappa _2$ , have gained relevance here (see also table 3). The former controls the location (in the toroidal domain) at which the magnetic field transitions from $B_{\mathrm{min}}$ to $B_{\mathrm{max}}$ , while $\kappa _2$ determines the position of peak curvature. Both these features have direct control on the trapped particle population and their drifts.

Table 3.

Most important combined input features for $f_{\mathcal{J}}$ . The table summarises the top selected features by FSFS on the combined input features detailed in Appendix B. The numerical value represents the coefficient of determination $R^2$ of the $f_{\mathcal{J}}$ model using the variable on that row alongside those above. The meaning of the various parameters in the table are presented in Appendix B: $\ell _{B_0}$ is the position of the $B_{\mathrm{min}}-B_{\mathrm{max}}$ transition, $\varDelta$ is the mirror ratio, $w_{\check {B}},w_{\hat {B}}$ are the toroidal extent of $B_{\mathrm{min}}$ and $B_{\mathrm{max}}$ , $B_0'$ is the sharpness of the $B_{\mathrm{min}}-B_{\mathrm{max}}$ transition, $\ell _{\kappa _{\mathrm{max}}}$ is the position of the curvature maximum, $\hat {\kappa }''$ is the second derivative of the curvature at $B_{\mathrm{max}}$ , $\check {\kappa }^{(3)}$ is the third derivative of the curvature at $B_{\mathrm{min}}$ , $\hat {\tau },\check {\tau }$ are the torsion at $B_{\mathrm{max}}$ and $B_{\mathrm{min}}$ , $\tau _0$ is the integrated torsion, $\check {\rho }$ is the value of $\rho$ at $B_{\mathrm{min}}$ and $\hat {\rho }''$ is the variation of elongation at $B_{\mathrm{max}}$ .

Besides curvature and $B_0$ , FSFS does also indicate the necessary involvement of torsion (particularly strong at higher $N$ ) and elongation. Both of these contributions should not come as a surprise given the involvement of $B_{20}$ ; see the analysis of $A_c^{\mathrm{mhd}}$ . Exploiting this, $f_{\mathcal{J}}$ can be improved (especially for large $N$ ) by minimising $\tau _{\mathrm{rms}}$ , and pushing variation of elongation near $B_{\mathrm{max}}$ (and not $B_{\mathrm{min}}$ ). Although all these features are clearly key and necessary, they hardly amount to a 50 %–70 % description of the total variance in the data, pointing at the complexity of the measure and the existing competing approaches to $f_{\mathcal{J}}$ .

3.3.3. Improving maximum $\mathcal{J}$ through additional shaping

The database behaviour above considered near-axis constructions achieving marginal stability minimally (Rodriguez & Plunk Reference Rodríguez, Plunk and Jorge2025). We noted that such a choice was suboptimal for $f_{\mathcal{J}}$ , and given the synergy between the magnetic well and precession, we expect shaping for larger $W$ values would improve $f_{\mathcal{J}}$ . As a matter of fact, a practical design of a stellarator would require a finite value of the magnetic well, and thus, such constructions are highly relevant. We present two examples of the behaviour of $f_{\mathcal{J}}$ for different target- $W$ values in figure 11.

Figure 11.

Evolution of $f_{\mathcal{J}}$ with different amounts of MHD stabilising shaping. The plots show the value of both the fraction of maximum $\mathcal{J}$ for different amounts of shaping tuned to achieve a magnetic well (negative values) desired, $W$ . Panel (a) corresponds to the highly shaped configuration in figure 9(c), while (b) is configuration (ii). The plots also show the explicit evolution of the shaping through $A_c^{\mathrm{mhd}}$ .

The figures show a clear improvement in $f_{\mathcal{J}}$ with increasing magnetic well (more negative $W$ ). The slope is however highly variable, opening the door to future considerations regarding robustness of configurations. As $f_{\mathcal{J}}$ is improved, the amount of shaping needed to achieve this does as well. Other physically relevant features, such as the omnigenous behaviour of the field, will also change.

Before closing this section, we note that these considerations of precession have been taken in the context of a vacuum. It is well known that finite plasma $\beta$ increases precession of trapped particles in the direction of maximum $\mathcal{J}$ . Thus, $f_{\mathcal{J}}$ will only increase.

3.4.1. Theoretical perspective on small $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$

By definition, the effective ripple is a weighted average of the net radial drift squared of trapped particles over the whole population. Configurations that minimise this drift will also tend to minimise $\epsilon _{\mathrm{eff}}$ , although contributions also depend in detail on the number of different trapped particles, which depends on how the magnetic field $B_0$ varies with $\varphi$ , its mirror ratio $\varDelta$ , etc.

Looking at the formal contributors to $\epsilon _{\mathrm{eff}}^{3/2,(2)}$ (Rodriguez & Plunk Reference Rodríguez, Plunk and Jorge2025, equation (4.11)),

(3.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

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

(3.12b) \begin{align} \Delta 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 ) \end{align}

and $B_{2c}^{\mathrm{QI}}=-(B_{2c}\cos 2\bar {\iota }_0\varphi +B_{2s}\sin 2\bar {\iota }_0\varphi )$ , where $\bar {\iota }_0 = \iota _0 - N/2$ .

The function $B_{2c}(\varphi )$ is central to $\epsilon _{\mathrm{eff}}$ , its pointwise value having to be of a very particular form, $\Delta B_{2c}^{\mathrm{QI}}=0$ , so as to ensure $\epsilon _{\mathrm{eff}}^{3/2,(2)}=0$ . Unlike many of the previously discussed measures, this is a local equality condition for all $\varphi$ . We cannot therefore easily come up with a line of attack at achieving low effective ripple by considering extreme limits. Instead we must take into account every contribution to the $T$ terms and $B_{2c}^{\mathrm{QI}}$ . Elements of curvature clearly arise in $T$ , with shaping playing a direct role in $B_{2c}$ . Generally, the choice of stabilising shaping does not minimise $\epsilon _{\mathrm{eff}}$ (Rodriguez & Plunk Reference Rodríguez, Plunk and Jorge2025), whose intrinsic behaviour we may also study as we did with $B_{20}$ . The value of $B_{2c}$ at $B_{\mathrm{min}}$ , where $B_{2c,\mathrm{b}}$  is partially correlated (with $R^2\sim 0.7$ ) to the ripple, can be written as, assuming perfect QI at first order (see (3.10) of Rodríguez et al. (Reference Rodríguez, Plunk and Jorge2025)), as

(3.13) \begin{equation} \left . \frac {B_{2c}}{B_0}\right |_{\varphi = \pi /N} =-\frac {1}{4(\ell ')^2}\frac {\bar {B}}{B_0}\left [\bar {e}''\left (1+\frac {1}{\bar {e}^2}\right )+\bar {e}(\tau \ell ')^2\left (3+\frac {1}{\bar {e}^2}\right )\right ]\!. \end{equation}

This shows that elongation and torsion are both important; after all, they control the bending of field lines over flux surfaces, and thus, the variation of $|\boldsymbol{B}|$ over them. From (3.13) and the QI condition, it is necessary that $\bar {e}''\leqslant 0$ and $|\tau |\geqslant 0$ for $B_{2c}$ to vanish, so if elongation grows away from the bottom of the well to aid with the magnetic well, then a finite amount of torsion is necessary. There is a clear tension between this and the behaviour of the magnetic well, as pointed out in Rodríguez et al. (Reference Rodríguez, Helander and Goodman2024) and Plunk et al. (Reference Plunk, Drevlak, Rodriguez, Babin, Goodman and Hindenlang2024).

3.4.2. Statistics of $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$

With the theoretical grounds laid, we consider the statistics of $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ . Aligning with the trends of previous features, the effective ripple of the bulk of the configurations becomes poorer the larger the number of field periods. In fact, the value of the population median is not far from the natural ${\sim}N^{8/3}$ scaling, making the bulk of the configurations have a forbiddingly large value of $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ .Footnote ¹⁵

Figure 12.

Behaviour of $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ across different field periods $N$ . (a) Dependence of the minimum and median value of $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ as a function of the number of field periods $N$ in the database. Reference scaling is given as a dotted line and distributions are represented as vertical violin plots. (b) Rendition of the 3D finite aspect ratio flux surface of $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ minimising fields for each number of field periods in the database. The locations of $B_{\mathrm{min}}$ and $B_{\mathrm{max}}$ are indicated for the $N=2$ configuration, which has the lowest value of $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ .

Unlike the bulk of the population, though, the best configurations (lowest values of $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ ) do not exhibit such a strong deterioration with $N$ . One is able to find configurations with values of ripple well below a percent for all $N$ , with the best candidate corresponding to $N=2$ . Depictions of the best-performing configurations are shown in figure 12(b). This stark difference between the bulk of the population and the extremes can be attributed to the necessity to carefully balance the different geometric elements of the field to achieve an omnigenous form of $|\boldsymbol{B}|$ , as we have discussed. If such balance is not struck then the mismatch will grow by an amount directly related to the magnitude of any of the ingredients involved, in this case the $N^2$ of the second-order field. Because of the need to strike such a balance, and the restricted forms used for the different inputs of the problem ( $\kappa$ $\tau$ , $\rho$ , etc.), it is likely that configurations exist beyond the database that achieve lower $\epsilon _{\mathrm{eff}}$ . Optimisation within a larger input feature near-axis space is left for future work.

The correlation analysis provides evidence of the complexity of the measure and the underlying careful balance. Although torsion of the axis appears to be most important, and in particular $\tau _{c1}$ (see figure 25), the relatively low value of the nonlinear correlation shows that there is no clear simple dominating trend. The importance analysis shows that other aspects such as elongation and magnetic shaping do also play an important role. We proceed with FSFS to further investigate this dependence, summarising the results in figure 13. The dependence on torsion is the most apparent, with a particular preference for its value at the bottom of the magnetic well, where its role was explicitly described by (3.13). Unsurprisingly, elongation also appears, although in no clear direction. The APM analysis of figure 13 reflects the requirement of a finite amount of torsion argued theoretically (see low $N$ ), but an emphasis to avoid large values at higher $N$ to avoid overwhelming shaping.

Figure 13.

The FSFS summary for $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ and the average percentile of the lowest $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ subspace. Panel (a) shows a table that summarises the top selected features by FSFS on the combined input features detailed in Appendix B. The numerical value represents the coefficient of determination $R^2$ of the $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ model using the variable on that row alongside those above. We considered $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ in log scale. (b) The average percentile (of the whole population across different field periods) to which the subset of the lowest 1000th percentile of $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ belongs. The meaning of the various parameters in the table are presented in Appendix B: $\hat {\tau },\check {\tau }$ is the torsion at $B_{\mathrm{max}}$ and $B_{\mathrm{min}}$ , $\check {\rho }$ is the value of $\rho$ at $B_{\mathrm{min}}$ , $\rho _{\kappa _{\mathrm{max}}}$ is the value of $\rho$ at the point of maximum curvature and $\int \kappa \mathrm{d}\ell$ is the integrated curvature over half a period.

As part of the FSFS analysis, we also obtain two other expected ingredients, as are the integrated curvature and the mirror ratio. We already touched on the physical grounds of both these elements, which we can connect to the radial drifts and trapped population size, respectively. Thus, it is reassuring to see them come out of the analysis. Despite this, we cannot say that configurations that achieve low ripple minimise their integrated curvature; see figure 13. Although the integrated curvature is kept small (especially at lower $N$ ), there is a clear shift toward prioritising low torsion as $N$ increases, reflecting the growing necessity to control shaping. Finally, the discreteness in the mirror ratio precludes any detailed conclusion regarding its role.

3.4.3. ‘Good’ configurations

The analysis above remains somewhat unsatisfactory, as no simple recipe to achieving low effective ripple has really been found. This is in part because of the many elements that go into the complex balance in $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ .

We could proceed like previously and study the behaviour of the lowest $\epsilon _{\mathrm{eff}}$ subset. However, we are in a position to pose a more interesting problem now. Instead of considering $\epsilon _{\mathrm{eff}}$ in isolation, we consider it alongside $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ and $A_c^{\mathrm{mhd}}$ , defining good configurations to be those that simultaneously satisfy

1. (i) $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}\gt 0.25$ ,

2. (ii) $A_c^{\mathrm{mhd}}\lt 10$ ,

3. (iii) $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}\lt 0.01$ .

A configuration meeting these criteria achieves a marginal magnetic well with minimal shaping and sufficiently small effective ripple (see previous sections for the informed choice of such thresholds). By studying the subclass of such configurations, we can see how constraining each of these requirements is in practice. An exhaustive study of physics trade-offs across the database is not presented in this paper of work.

Figure 14.

‘Good’ configurations in the database. (a) A 3D rendition of the largest $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ configurations in the ‘good’ configuration subset per field period. (b) Example of the variety of good $N=2$ configurations in a lower-dimensional MDS representation showing larger and lower torsion clusters. Three examples are illustrated numbered (i)–(iii), bearing features of the configurations excelling in the various physics features discussed before. (c) Average percentile measure across different $N$ for ‘good’ configurations. Note that there is no ‘good’ configurations for $N=6$ . We also indicate as reference the marginal APM values, meaning the APM of the subset satisfying each of the thresholds separately, to illustrate the influence of each threshold on the others.

‘Good’ configurations with the largest $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ for each $N$ are shown figure 14. These are prime candidates as initial conditions of further stellarator optimisation pursues. But that is not our main focus here; we are interested in learning how restrictive each of the requirements is, i.e. how strongly they shape the overall population. We thus apply our APM to the good-configuration subsets (see figure 14 c). The shaping and ripple measures behave quite differently through $N$ . The variation of the $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ APM with $N$ appears weakest, but it proves to be the single most limiting feature: the measure stays below ${\sim}5\,\%$ for all $N$ . This is particularly noticeable at the lowest $N$ , where the APM of both $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ and $A_c^{\mathrm{mhd}}$ are significantly larger (even though the $\epsilon _{\mathrm{eff}}$ constraint significantly hampers $A_c^{\mathrm{mhd}}$ , as shown in figure 14). This is consistent with the loss of standing of the figure-eight-like configurations over the less inclined configurations at low $N$ (even though there still remains an ample variety of configurations, as shown in figure 14 b). The implications of shaping and coil compatibility requirements ( $L_{{\boldsymbol{\nabla }\kern-1.2pt\boldsymbol{B}}}$ and $A_c^{\mathrm{mhd}}$ ) significantly worsen at larger $N$ , becoming similarly constraining as $\epsilon _{\mathrm{eff}}^{\mathrm{edge}}$ . Ultimately, no $N=6$ configuration exists within the database that satisfies all requirements to be considered ‘good’. It is worth noting that additional low-ripple configurations might be found by extending the input-model flexibility, which could alleviate some of the burden imposed by the requirement on $\epsilon _{\mathrm{eff}}$ .