Email<EMAIL_ADDRESS>
Email<EMAIL_ADDRESS>
Email<EMAIL_ADDRESS>
# Constructing the space of quasisymmetric stellarators

E. Rodríguez  [  Max Planck Institute for Plasma Physics, Greifswald, Germany,
17491  W. Sengupta  [  Department of Astrophysical Sciences, Princeton
University, Princeton, NJ, 08543  Princeton Plasma Physics Laboratory,
Princeton, NJ, 08540  A. Bhattacharjee  [  Department of Astrophysical
Sciences, Princeton University, Princeton, NJ, 08543  Princeton Plasma Physics
Laboratory, Princeton, NJ, 08540

###### Abstract

A simplified view of the space of optimised stellarators has the potential to
guide and aid the design efforts of magnetic confinement configurations
suitable for future fusion reactors. We present one such view for the class of
quasisymmetric stellarators based on their approximate description near their
centre (magnetic axis). The result is a space that captures existing designs
and presents new ones, providing a common framework to study them. Such a
simplified construction offers a basic topological approach, guided by certain
theoretical and physical choices, which this paper presents in detail.

††preprint: AIP/123-QED

##  I  Introduction

Designing magnetic configurations suitable for holding a thermonuclear plasma
is central to fusion research. However, finding the appropriate shape of
magnetic fields with the desired properties is challenging. There are
primarily two reasons for this difficulty. First, one must select an
appropriate set of objectives, imposed both by physics and technology, which
may not be mutually compatible. Second, the potential parameter space for
three-dimensional fields is very large.

Most three-dimensional configurations are not attractive candidates, mainly
because of their poor confinement properties. In an inhomogeneous magnetic
field, charged particles generally drift away from field lines unless the
field is carefully designed. The class of stellarators where charged particles
are, on average, confined (collisionlessly), are called omnigeneous[1, 2, 3,
4, 5]. This paper focuses on a subgroup of omnigeneous stellarators known as
quasisymmetric stellarators [6, 7, 8]. In these configurations, the magnitude
of the magnetic field (but not the full vector field) has a direction of
symmetry. The conventional approach to finding this particular subset of
configurations remains an extensive search in the space of all possible
stellarators. The search attempts to ‘minimise’ the asymmetries in the
magnetic field magnitude and has successfully provided multiple designs.[9,
10, 11, 12, 13] However, the optimisation procedure remains, to a large
extent, a black box. This leaves important questions about quasisymmetry and
its implications unanswered. Moreover, and perhaps most importantly, the black
box can miss out on a significant number of undiscovered designs of
potentially high value.

This paper presents an attempt to shed some light on these questions by
considering an alternative view on the configuration space of quasisymmetric
stellarators. This alternative view consists of a model based on an
approximate description of stellarators close to their magnetic axis, where
the complexity of the stellarator is significantly reduced. In such a
framework, potential quasisymmetric stellarators are reduced to a combination
of a magnetic axis shape and the choice of two scalar parameters, in terms of
which many relevant properties can be expressed. Such a description provides
the space with a basic topological structure that enables deeper
understanding. Several choices are necessary for the axes of such a model to
represent quasisymmetric stellarators. This paper presents the formal and
physical basis underlying such choices, as well as the resulting highlights of
the approach.

Section II introduces the basics of the truncated near-axis expansion that
constitutes the basic quasisymmetric stellarator model used. Sections III-V
take the elements of such a model and explore their physical implications,
guiding the appropriate choice of parameters for the model. The focus here
will not be on the axis shapes, on which a thorough discussion may be found
elsewhere[14], but rather on choosing the remaining parameters that complete
the model. Section VI then presents an example of QS configuration space,
which we analyse to illustrate the potential of this approach. We conclude
with some final remarks and open questions.

##  II Quasisymmetry and near-axis expansions

We begin by introducing the notion of quasisymmetry(QS). This hidden symmetry
is the minimal property of a magnetic field that provides the dynamics of
charged particles with an approximate conserved dynamical quantity (Tamm’s
theorem) to leading order in the gyroradius.111We simplify the picture by
assuming that the electrostatic potential shares the QS to leading gyro-order
and do not include it in our considerations. Such a conserved momentum
prevents particles from freely escaping the magnetic field, making the concept
naturally attractive for magnetic confinement.[16, 17, 8] The condition can be
formally expressed as $\nabla\psi\times\nabla
B\cdot\nabla(\mathbf{B}\cdot\nabla B)=0$, where $2\pi\psi$ is the toroidal
flux and $\mathbf{B}$ is the magnetic field. However, this form hides the
underlying nature of QS, which is to make the contours of $|\mathbf{B}|$
symmetric. Under the assumption of ideal magnetohydrostatic equilibrium,
$\mathbf{j}\times\mathbf{B}=\nabla p$, where $\mathbf{j}$ is the current
density, in Boozer coordinates[6], QS implies[18] that
$|\mathbf{B}|=B(\psi,\chi=\theta-N\phi)$ is a function that depends on a
linear combination of Boozer angles.222For a more general form of equilibrium,
a formally analogous approach exists in terms of so-called generalised Boozer
coordinates, details of which may be found in [24]. As a result, many of the
properties of quasisymmetric stellarators in this paper are independent of a
particular form of equilibrium. Here, $N\in\mathbb{Z}$ describes the pitch of
the symmetry, which leads to the distinction between quasi-axisymmetric (QA,
$N=0$) and quasi-helically symmetric (QH, $N\neq 0$) stellarators.

Our goal is to describe stellarators with the property of QS (in the second
form discussed above) close to the magnetic axis. The magnetic axis is the
centre of the stellarator, a closed magnetic field line around which magnetic
flux surfaces accrue. Because we are considering a description of the
stellarator near its axis, it is natural to use the axis as reference for our
coordinate system. We describe flux surfaces (i.e., constant $\psi$ surfaces,
which we assume to be nested) using the Frenet-Serret basis[20, 21, 22]
$\\{\hat{b},\hat{\kappa},\hat{\tau}\\}$ and Boozer coordinates
$\\{\psi,\theta,\phi\\}$, so that with respect to the axis $\mathbf{r}_{0}$,
we write

$\mathbf{x}=\mathbf{r}_{0}+X(\psi,\theta,\phi)\hat{\kappa}+Y(\psi,\theta,\phi)\hat{\tau}+Z(\psi,\theta,\phi)\hat{b},$
(1)

where $X,~{}Y$, and $Z$ are functions of all Boozer coordinates. This paper
uses the notation in [22], including the convention on the sign of the
torsion. Not every flux surface shape described by Equation (1) is consistent
with a given divergenceless magnetic field, which must be in equilibrium and
be quasisymmetric. Let us start by imposing the condition that the magnetic
field is solenoidal ($\nabla\cdot\mathbf{B}=0$) and lives on flux surfaces
($\mathbf{B}\cdot\nabla\psi=0$) formally[23] by writing both the covariant and
contravariant forms of $\mathbf{B}$ and using Boozer coordinates as
independent coordinates, [20, 22, 24],

$\displaystyle(B_{\alpha}(\psi)-$
$\displaystyle\bar{\iota}B_{\theta})\frac{\partial\mathbf{x}}{\partial\psi}\times\frac{\partial\mathbf{x}}{\partial\chi}+B_{\theta}\frac{\partial\mathbf{x}}{\partial\phi}\times\frac{\partial\mathbf{x}}{\partial\psi}+B_{\psi}\frac{\partial\mathbf{x}}{\partial\chi}\times\frac{\partial\mathbf{x}}{\partial\phi}=$
$\displaystyle=\frac{\partial\mathbf{x}}{\partial\phi}+\bar{\iota}(\psi)\frac{\partial\mathbf{x}}{\partial\chi},$
(2)

where $\mathbf{x}$ is defined in Eq. (1). Here $\bar{\iota}=\iota-N$ and
$\iota$ is the rotational transform, and $B_{i}$ are the covariant components
of the magnetic field.

In Boozer coordinates, a quasisymmetric magnetic field satisfies

$\frac{B_{\alpha}(\psi)^{2}}{B(\psi,\chi)^{2}}=\left|\frac{\partial\mathbf{x}}{\partial\phi}+\bar{\iota}\frac{\partial\mathbf{x}}{\partial\chi}\right|^{2}.$
(3)

As given, Eqs. (2) and (3) constitute a coupled set of partial differential
equations (PDEs) describing a quasisymmetric magnetic field, not only near the
magnetic axis but everywhere. Unless we consider them close to the axis, the
system of equations is overly complicated.

The asymptotic description of the fields near the magnetic axis is known as
the near-axis expansion[25, 26], in the form here presented pioneered by [20].
This procedure entails expanding all relevant fields in the problem as power
series in the distance from the magnetic axis. A pseudo-radial coordinate
$\epsilon=\sqrt{\psi}$ is defined, where $2\pi\psi$ is the toroidal flux,
which serves as the appropriate ordering parameter.333For simplicity, we have
not normalised $\psi$ respect to the magnetic field on the axis and its
curvature as it is often customary[20, 21]. Doing so simply introduces
constant rescaling factors in various quantities involved. Because of its
radial nature, the expansion in $\epsilon$ must be carefully coupled to the
poloidal-angle, $\theta$, behaviour. To avoid a coordinate singularity on the
magnetic axis, all physical quantities must have the following asymptotic
form,

$f=\sum_{n=0}^{\infty}\epsilon^{n}{\sum_{m=0|1}^{n}}\left[f_{nm}^{c}(\phi)\cos
m\chi+f_{nm}^{s}(\phi)\sin m\chi\right].$ (4)

If the function $f$ is a flux function, the expansion reduces to a Taylor
expansion in $\epsilon^{2}$. In this paper, we will use this subscript
notation repeatedly.

Expanding all functions in Eqs. (2)-(3) as the coupled Fourier (in $\chi$)
-Taylor (in $\epsilon$) series described by Eq. (4), the PDEs are reduced into
an (a priori) infinite ordered set of ordinary differential and algebraic
equations. By order $n$, we are referring here to all the elements in the
problem that have the same power $\epsilon^{n}$. At each order, various
parameters and functions are needed as inputs to the equations to uniquely
determine the solution, as summarised in Table 1. The main elements are the
magnetic field on the axis ($B_{0}$) and its leading variations nearby
($\eta,~{}B_{22}^{C}$ and $B_{22}^{S}$), the shape of the axis, the leading
contribution of the toroidal current ($B_{\theta 20}$), the stellarator-
symmetry breaking ($\sigma(0)$) and the pressure gradient ($B_{\alpha 2}$).
Each choice represents a different stellarator; in that sense, the framework
described here serves as a reduced stellarator model. We do not present the
detailed order-by-order set of equations that constitute the near-axis
description, as these may be found elsewhere, both in the form that concerns
us here[21, 28] and in the context of more general equilibria[29, 30].
Instead, we focus on the choice of parameters that make the near-axis model
represent optimised quasisymmetric configurations.

Order | Params.
---|---
$0$ |  $B_{0}$, axis ($\kappa,~{}\tau,~{}l)$
$1$ |  $B_{\theta 20}$, $\sigma(0)$, $\eta$
$2$ |  $B_{22}^{C}$, $B_{22}^{S}$, $B_{\alpha 2}$
Table 1: Quasisymmetric configuration characterising parameters. The table
gathers the free parameters (and functions) defining the leading order form of
quasisymmetric configurations.

##  III Zeroth order: magnetic axis

Let us start with the most basic element in the model: the shape of the
magnetic axis (see Tab. 1). At a fundamental level, the near-axis model
identifies every configuration with a three-dimensional closed curve (magnetic
axis). Every configuration that shares the same magnetic axis must then also
share certain properties. A detailed discussion on the role of the magnetic
axis in the context of quasisymmetric stellarators was presented in [14]. Here
we reproduce the essential elements concerned in constructing our
quasisymmetric stellarator model.

From the set of all smooth, three-dimensional closed curves, those with
inflection points, that is, points of vanishing curvature, must be excluded.
To see why this is the case, interpret the curvature of the axis, $\kappa$, as
a measure of the scale of the magnetic field magnitude gradient near the axis
(from equilibrium $\nabla_{\perp}(B^{2}/2)=B^{2}\vec{\kappa}$). Thus, to
support any finite magnetic field variation on a flux surface around a point
where $\kappa\approx 0$, an unphysical, nearly infinitely elongated flux
surface is necessary[31, 14]. Under such conditions, the set of excluded axis
shapes becomes physically interesting only outside the rigorous realm of QS,
as in the case of quasi-isodynamic stellarators[3, 32, 33]. In QS, which is
the focus of the present paper, the requirement of specialising to regular
curves makes the Frenet-Serret frame well-defined everywhere.

With such a frame defined, we can construct a self-linking number[14, 34, 35,
36, 37] ($S_{L}$), which is the number of times the curvature vector of the
axis encircles itself in a full toroidal excursion. This number is precisely
the integer $N\in\mathbb{Z}$ that appears in the quasisymmetric form of
$|\mathbf{B}|=B(\psi,\theta-N\phi)$. Thus, the shape of the axis (a local
feature) fully determines the class of the quasisymmetric stellarator (a
global feature)[14, 31].

Besides its global implications, the strength of this association is that
$S_{L}$ is a topological invariant under regular isotopies. The property
remains unchanged if the magnetic axis is continuously deformed provided that
the curvature does not vanish anywhere along the way. Thus, the space of all
closed curves is partitioned into regions identified by an integer value of
$S_{L}$. Each quasisymmetry class thus behaves as a quasisymmetric phase,
whose nature may only be changed by crossing the separating manifolds, phase-
transitions, made up of curves with inflection points. It is then natural to
see each phase as a distinct class.

Lacking a theory that directly relates axis shapes to quasisymmetric quality
of the stellarator, we must consider all (regular) axis shapes as part of our
space of potential QS configurations. Elements from higher orders in the near-
axis expansion will be needed to tell different shapes apart and delve into
additional properties of the configurations. It is important to note the
difference between this model space and the one more traditionally employed in
optimisation: the space of toroidal surface shapes. The consideration of the
axis reduces the dimensionality of the space and, crucially, provides such a
space with topological structure [14]. It is evident from the latter that an
optimisation that starts in a particular phase will remain within it.

There is no unique way to represent this space of closed curves, and we
consider a Fourier description of the curves in cylindrical coordinates
$(R(\phi),Z(\phi))$ for simplicity. Here, $\phi$ is the cylindrical angle, and
$R=\sum_{n}R_{n}\cos nN\phi$ and $Z=\sum_{n}Z_{n}\sin nN\phi$. We specialise
to stellarator-symmetric configurations. This parametrisation of the curves
guarantees they are closed and have $N-$fold symmetry. However, the torsion,
$\tau$, and curvature, $\kappa$, which are most directly involved in the near-
axis construction become byproducts that (as is the case of torsion) can be
pretty sensitive to the choice of Fourier harmonics. Parametrising the curves
providing $\kappa$ and $\tau$, which we may call the Frenet approach, would be
more natural (and in agreement with the fundamental theorem of curves[38]) but
suffer from the issue of closing the curve. Other possible alternatives, such
as the use of control points of splines, have generated renewed interest[39,
40]. The story of representations is mixed: the advantages of one appears to
be the Achilles heel of the other. However, the foundations of this paper are
independent of the particular form of a representation.

##  IV First order: elliptic shaping

With the magnetic axis in place, let us proceed to the first order. Following
Table 1, we must choose appropriate values for three parameters: the toroidal
current $B_{\theta 20}$, the stellarator-symmetry breaking $\sigma(0)$, and
the $|\mathbf{B}|$ variation, $\eta$. Different choices will describe
different configurations, thereby lifting part of the degeneracy that results
from identifying stellarators with their magnetic axis. In this order, the
model introduces elements of flux surface shaping and the rotational
transform. We will see how the choice of parameters in this order affects the
behaviour of these features, and provide guidelines for appropriate choices.

###  IV.0.1 Toroidal current, $B_{\theta 20}$

The coefficient $B_{\theta 20}$ controls the plasma current density on the
axis and is, in more familiar notation, the leading contribution to
$B_{\theta}=I(\psi)$. In constructing most QS configurations, the trivial
assumption $B_{\theta 20}=0$ is often made. We shall do so here as well. Of
course, such a choice is exact when modelling a vacuum field, but it is not
always the appropriate limit for configurations that support a finite plasma
pressure. Plasma currents, such as bootstrap currents, may be present in that
case, even without having to drive them externally. Their evaluation generally
requires a separate kinetic consideration [41, 42]. Bearing this caveat in
mind, we specialise to vacuum magnetic fields and thus take $B_{\theta 20}=0$.
The problem does not generally show singular behaviour in this limit. Thus,
this simplifying assumption is justified for illustrating the approach in this
paper.

###  IV.0.2 Surface shaping and rotational transform: $\eta$ and $\sigma(0)$

We now consider the relevant parameters $\eta$ and $\sigma(0)$. By definition,
the parameter $\eta=-B_{11}^{C}/2B_{0}$ is a measure of the variation of the
magnetic field over flux surfaces. The other parameter, $\sigma(0)$, is
$\sigma=Y_{11}^{C}/Y_{11}^{S}$, a quantity related to the shaping of flux
surfaces at the point of stellarator symmetry $\phi=0$. Beyond these
definitions, both of these parameters may be connected to natural elements of
the geometry of flux surfaces (see [43] for more details). The shapes of the
flux surfaces at the first order are elliptical and can be characterised on
the plane orthogonal to the magnetic axis by an elongation $\mathcal{E}$, and
rotation angle, $\vartheta$,

$\displaystyle\mathcal{E}=\frac{F}{\eta^{2}/\kappa^{2}}\left[1+\sqrt{1-\frac{\eta^{4}/\kappa^{4}}{F^{2}}}\right],$
(5a) $\displaystyle\tan
2\vartheta=\frac{\sigma\eta^{2}/\kappa^{2}}{\eta^{4}/4\kappa^{4}-1-\sigma^{2}},$
(5b)

where $F=1+\sigma^{2}+\eta^{4}/4\kappa^{4}$. Here elongation is defined as the
ratio of the major to the minor radius, and $\vartheta$ is the angle between
the major radius and the positive $X$ direction (i.e., the normal vector).

From Eq. (5b), it is clear that the parameter $\sigma(0)$ serves as a measure
of the misalignment of the ellipse with the Frenet-Serret frame at the origin,
$\phi=0$. That is, it is a measure of the up-down asymmetry of this cross-
section and, generally, of stellarator asymmetry. If we specialise to the
stellarator-symmetric choice ($\sigma(0)=0$), a single parameter is then left
to choose, $\eta$.

From the form of $\mathcal{E}$, Eq. (5a), in the limit of $\sigma=0$, where
$\mathcal{E}\sim\eta^{2}/2\kappa^{2},~{}2\kappa^{2}/\eta^{2}$ (the latter for
$\eta<\kappa\sqrt{2}$), we may interpret $\eta$ as an approximate measure of
elongation. A large $\eta$ indicates a large elongation in the direction of
the curvature vector, while a small $\eta$ corresponds to large elongation
along the binormal. This correspondence is only approximate because even with
the stellarator-symmetric choice $\sigma(0)=0$, $\sigma$ is generally a non-
zero function of $\phi$, which makes Eq. (5a) highly non-trivial.

This naturally takes us to the equation that governs the form of $\sigma$, a
first-order non-linear differential equation,

$\frac{\mathrm{d}\sigma}{\mathrm{d}\phi}=-\bar{\iota}_{0}\left[1+\sigma^{2}-\frac{1}{4B_{0}}\left(\frac{\eta}{\kappa}\right)^{4}\right]+\frac{B_{\alpha
0}}{2}(2\tau+B_{\theta 20})\left(\frac{\eta}{\kappa}\right)^{2},$ (6)

referred to as the $\sigma$ Riccati equation[28, 21, 29]. The choice of $\eta$
affects $\sigma$ non-trivially, and thus the shaping of the flux surfaces.
Following Mercier[25], we then expect the rotational transform of the
configuration on the axis, $\iota_{0}$ (or more generally in quasisymmetry,
$\bar{\iota}_{0}=\iota-N$, where $N$ is the self-linking number of the axis)
to depend on $\eta$ as well. Formally, $\bar{\iota}_{0}$ forms part of the
solution to the Riccati equation, Eq. (6), as there is a unique value for
which $\sigma(\phi)$ is periodic[21]. Thus, through the lens of this equation,
we may regard the rotational transform as a function of $\eta$ for a fixed
axis shape (see Figure 1).

Figure 1: Example of $\bar{\iota}_{0}$ as a function of $\eta$. Example of the
rotational transform as a function of $\eta$ for a fixed axis shape. This
example corresponds to a quasiaxisymmetric configuration, with an axis like
the precise QA in [44]. The broken lines indicate the asymptotic behaviour
described in the text.

No closed form expression can generally be obtained for
$\bar{\iota}_{0}(\eta)$. However, the behaviour of $\bar{\iota}_{0}(\eta)$ for
large and small values of $\eta$ can be determined. In the small $\eta$ limit,
the dominant balance $\sigma\sim\eta^{2}$ and $\bar{\iota}_{0}\sim\eta^{2}$ of
Eq. (6) yields, upon integration,

$\bar{\iota}_{0}\sim\bar{\iota}_{-}\eta^{2}=\frac{B_{\alpha
0}}{4\pi}\eta^{2}\int_{0}^{2\pi}\frac{2\tau+B_{\theta
20}}{\kappa^{2}}\mathrm{d}\phi.$ (7)

Thus, the rotational transform on the axis tends to $\iota_{0}\rightarrow N$.
In the case of a quasi-axisymmetric (QA) configuration, this leads to a
vanishing rotational transform. Physically, this is a result of making flux
surfaces very elongated in the curvature direction in a way that, from
Mercier’s perspective, the rotational transform on the axis is driven only by
the rotating-ellipse contribution. Of course, in the QA case, the ellipse has
no net rotation (recall the meaning of the self-linking number).444From the
Mercier perspective on rotational transform[5],
$\bar{\iota}_{0}=-\frac{1}{2\pi}\int_{0}^{L}\frac{\cosh\bar{\eta}-1}{\cosh\bar{\eta}}(d^{\prime}+\tau)\mathrm{d}l+\frac{1}{2\pi}\int_{0}^{L}\tau\mathrm{d}l,$
in the large ellipticity $\bar{\eta}\rightarrow\infty$ limit,
$\bar{\iota}_{0}=-d/2\pi$, where $d$ is the angle of rotation of the ellipse
with respect to the curvature vector. Because $\sigma\sim\eta^{2}\rightarrow
0$ in this limit, the cross-sections align with the Frenet-Serret frame. Thus,
the net rotation $d=0$. Thus, $\iota_{0}=0$. As $\eta$ increases, the torsion
and current contributions play the same role in driving rotational transform.
In a QA configuration, this growth in $\bar{\iota}_{0}$ is equivalent to
growth in rotational transform. However, in the case of a quasi-helically
symmetric (QH) configuration, this depends on the relative sign of Eq. (7) and
the QS helicity $N$, which tends, in practice (and zero current), to be
negative. In the limit of large $\eta$, the physical scenario is similar to
that in the small $\eta$ limit. The shaping becomes large, and
$\bar{\iota}_{0}$ also tends to zero. This follows formally from the dominant
balance $\bar{\iota}_{0}=\bar{\iota}_{0}^{*}/\eta^{2}$ and
$\sigma=\eta^{2}\sigma^{*}$ in Eq. (6). The asymptotic forms of
$\\{\sigma^{*},\bar{\iota}_{0}^{*}\\}$ are solutions to a ‘universal’
$\eta$-independent Riccati equation, unique to each axis shape
($\kappa,~{}\tau$).

By the mean value theorem, it follows that $\iota_{0}(\eta)$ must have at
least one turning point (except in the marginal case of $\int(2\tau+B_{\theta
20})/\kappa^{2}=0$). In Appendix A we prove that this extremum exists, its
value $\eta=\eta^{*}$ is unique and satisfies the condition,

$\int_{0}^{2\pi}\left[\sigma^{2}+\frac{1}{4B_{0}}\left(\frac{\eta}{\kappa}\right)^{4}-1\right]E\mathrm{d}\varphi=0,$
(8)

where
$E=\exp[{2\bar{\iota}_{0}\int_{0}^{\varphi}\sigma\mathrm{d}\varphi^{\prime}}]$.
From this expression, it follows that there is a value of $\eta$ in
$(0,\sqrt{2}\kappa_{\mathrm{max}}B_{0}^{1/4}]$ that extremises the rotational
transform (see Fig. 1).

To illustrate further the meaning of $\eta^{*}$, and to further connect it to
shaping, consider the familiar limit of axisymmetry. In that case, the Riccati
equation, Eq. (6), can be solved exactly (taking $B_{\alpha 0}=1=R_{0}$) to
yield

$\iota_{0}=\frac{1}{2}\frac{B_{\theta 20}\eta^{2}}{1+\eta^{4}/4}.$ (9)

The solution in Eq. (9) exhibits all the properties of $\bar{\iota}_{0}$ we
studied in the general case, including the uniqueness of $\eta^{*}=\sqrt{2}$.
This choice of $\eta$ corresponds to having a circular cross-scetion. This
suggests that in order to maximise the use of the toroidal current to generate
rotational transform, one must simultaneously minimise the amount of shaping.

The correspondence between the behaviour of the rotational transform and the
shaping of surfaces prevails beyond axisymmetry. This can be seen by
investigating the behaviour of elongation $\mathcal{E}$, Eq. (5a), with $\eta$
in an analogous form to the case of $\bar{\iota}_{0}$. Its asymptotic
behaviour was already touched upon before, as we saw that elongation diverged
both for small $\eta$ ($\mathcal{E}\approx 2\sqrt{B_{0}}(\kappa/\eta)^{2}$)
and for large $\eta$ ($\mathcal{E}\sim
2\eta^{2}\sqrt{B_{0}}\kappa^{2}[\tilde{\sigma}^{2}+1/4B_{0}\kappa^{4}]$).
Thus, there must be some minimally elongated configuration somewhere in
between.

Figure 2: Difference in behaviour between elongation measures. Plot showing
the behaviour of different elongation measures with $\eta$ for the magnetic
axis of the ‘precise QA’ configuration[44]. The plot shows the difference in
behaviour between different global scalar measures for elongation:
$\mathcal{E}(0)$ is the value of elongation at $\varphi=0$,
$\langle\mathcal{E}\rangle$ is the $\phi$-average of elongation,
$\mathrm{max}_{\varphi}[\mathcal{E}]$ is the maximum of elongation and
$\langle\mathcal{E}/\kappa^{2}\rangle$ is the weighted average of elongation.
The broken line represents the value of $\eta$ obtained from the global
solution for comparison.

However, this notion must be qualified, as $\mathcal{E}$ is not a single
scalar but a function of $\varphi$. This allows for different ways in which to
construct a single measure for the shaping of the configuration. Figure 2
shows some possibilities. It is evident that the behaviour changes from one
definition to another. Among these measures, there is one that is particularly
insightful. Consider the weighted average
$\bar{\mathcal{E}}=(1/2\pi)\int(\mathcal{E}/\kappa^{2})\mathrm{d}\phi$ which
penalises elongation in the straighter sections of the configuration. This
particular elongation measure has the same asymptotic behaviour (up to an
$\eta-$independent factor) as $\bar{\iota}_{0}$. This formal analogy makes the
choice of $\eta^{*}$ also a choice that roughly minimises elongation (see the
comparison in Fig. 8 or Table 2 for the implications in practice).

Besides all the physical implications of $\eta^{*}$, this choice also
represents a least sensitive choice, so much so that it corresponds to a
turning point. That is, a set of near-axis constructions in some range
$(\eta^{*}-\delta,\eta^{*}+\delta)$ have roughly the same properties, making
their properties more robust. Following [29], this resilience of the
rotational transform makes the magnetic shear (a quantity that would be, in
general, a third-order quantity) independent of third-order parameter choices
(see Appendix B). Thus, our model can predict the behaviour of magnetic shear
as well.

| ARIESCS | ESTELL | GAR | HSX | NCSX | QHS48 | Precise QA | Precise QH
---|---|---|---|---|---|---|---|---
$\eta_{\mathrm{VMEC}}$ | 0.11 | 0.79 | 0.51 | 1.75 | 0.62 | 0.21 | 0.96 | 2.12
$\eta_{\mathrm{VMEC}}/\eta^{*}$ | 1.00 | 0.91 | 0.95 | 0.87 | 0.92 | 0.88 | 0.95 | 0.91
$\eta_{\mathrm{VMEC}}/\eta_{\langle\nabla\mathbf{B}\rangle}$ | 1.03 | 0.12 | 0.57 | 0.99 | 0.73 | 1.07 | 0.64 | 1.04
$\eta_{\mathrm{VMEC}}/\eta_{\bar{\mathcal{E}}}$ | 1.07 | 0.95 | 1.01 | 0.91 | 0.97 | 0.91 | 1.03 | 0.95
$N$ | 0 | 0 | 0 | 4 | 0 | 4 | 0 | 4
Table 2: Comparison of $\eta$ choices in some quasisymmetric designs.
Comparison of the $\eta$ values for many quasisymmetric designs: ARIESCS[11],
ESTELL[46], GAR[47, 48], HSX[9], NCSX[10], QHS48[12], Precise QA and QH[44].
$\eta_{\mathrm{VMEC}}$ corresponds to the value of $\eta$ obtained from the
$|\mathbf{B}|$ of the global VMEC[49] equilibrium solution of the
quasisymmetric designs. $\eta^{*}$ is the parameter value that extremises
$\bar{\iota}$ for a fixed axis obtained from VMEC.
$\eta_{\langle\nabla\mathbf{B}\rangle}$ corresponds to the maxima of the
average of $L_{\nabla\mathbf{B}}$ (see Appendix C), and finally
$\eta_{\bar{\mathcal{E}}}$ from the minimum of the weighted average elongation
$\bar{\mathcal{E}}$. The integer $N$ denotes the symmetry class.

In summary, the choice of $\eta=\eta^{*}$ is a formally convenient and
representative choice for $\eta$. It is unique and always exists, maximises
rotational transform in QA configurations, and regularises the shaping of flux
surfaces, preventing them from having extreme shaping. With such a choice,
every axis in our structured space of configurations will have an associated
natural choice of $\eta$. We do not need to keep this parameter explicitly as
an added dimension in this space, reducing its complexity. A practical design
would benefit from an additional tweaking of this parameter in the
neighbourhood of $\eta^{*}$. This refinement could be seen as a subsequent
optimisation in which one can consider higher-order properties or more
sophisticated construction specifications.

To present some evidence that backs the suitability of the choice of
$\eta^{*}$, we present in Table 2 the value of $\eta$ extracted from global
equilibrium solutions designed through other means compared to different
choices of $\eta$. This shows that $\eta^{*}$ is a reasonable representative
choice across the board. In the Table, we also include a choice
$\eta_{\nabla\mathbf{B}}$, which we have not mentioned in the discussion. This
choice of $\eta$ is taken to be the value that minimises the magnitude of
$||\nabla\mathbf{B}||$, a measure of the gradients of the magnetic field. This
measure has been used by other researchers[50, 44, 51] to good effect in
optimising near-axis QS configurations. The rationale behind this measure is
that the characteristic length scale
$L_{\nabla}\stackrel{{\scriptstyle\cdot}}{{=}}1/||\nabla\mathbf{B}||$, may be
interpreted to approximately set the distance to the field-generating coils,
and thus to indicate a measure of the range of validity of the near-axis
model.[50] Guiding the choice of $\eta$ for a fixed axis could be considered
as a guiding principle in choosing $\eta$. This is successful in many cases
(see Tab. 2), but it lacks the robustness and generality of $\eta^{*}$, and
fails in the QA phase. See Appendix C for a more detailed discussion.

##  V Second order: surface triangularity and QS breaking

So far, we have identified quasisymmetric configurations with a model
consisting of an axis and elliptical cross-sections. For every curve in our
space of configurations, the shape of the latter results from the choice of
the parameter $\eta^{*}$. However, we note that this model is exactly
quasisymmetric up to this point, and that there is no way of differentiating
which configurations will exhibit better QS globally. Consideration of the
second order in the near-axis expansion is needed for this.

When incorporating the second order, two important parameter choices must be
made: (i) the plasma pressure gradient (in the notation here, related directly
to $B_{\alpha 1}$) and (ii) the second order variation of the magnetic field
magnitude, $B_{2}$. For every choice made, the model gains a different flux-
surface shaping in the form of triangularity, Shafranov shift, and a different
degree of QS. In line with the simplifying vacuum and stellarator-symmetric
assumptions, we choose $B_{\alpha 1}=0$ and $B_{22}^{S}=0$ to reduce the
number of free parameter choices to one. Only one of the harmonics of $B_{2}$
remains. Following the approach of [21], this parameter is taken to be
$B_{22}^{C}$. This leaves the remaining component of $B$ at second order,
$B_{20}$, to be found self-consistently. This lack of freedom in $B_{20}$
follows from the necessity of satisfying force-balance and the appropriate
magnetic equations simultaneously, which generally requires $B_{20}$ to have a
toroidal angle $\varphi$ dependence.[28] Of course, such a variation violates
QS, and the variation of $B_{20}$, which we may define as $\Delta B_{20}$
bottom-to-peak, becomes a measure of QS quality at this order.

Formally, $B_{20}$ is the solution to a second-order, linear differential
equation. The $B_{20}$ equation may be found in the appendix of [30], and may
be written as,

$\mathcal{A}\left(\frac{B_{20}}{B_{0}}\right)^{\prime\prime}+\mathcal{B}\left(\frac{B_{20}}{B_{0}}\right)^{\prime}+\mathcal{C}\frac{B_{20}}{B_{0}}+\mathcal{D}=0,$
(10)

where

$\displaystyle\mathcal{A}=-\frac{B_{\alpha
0}\eta}{2\kappa^{2}\bar{\iota}_{0}l^{\prime}}\left[1+\frac{4B_{0}\kappa^{4}}{\eta^{4}}(1+\sigma^{2})\right],$
(11a) $\displaystyle\mathcal{B}=\frac{2B_{\alpha
0}\eta}{\bar{\iota}_{0}l^{\prime}}\frac{\kappa^{\prime}}{\kappa^{3}}-\frac{4l^{\prime}\sigma}{\bar{\iota}_{0}\eta}\tau,$
(11b) $\displaystyle\mathcal{C}=-\frac{l^{\prime}}{2B_{\alpha
0}\eta^{3}\kappa^{2}}\left[\bar{\iota}_{0}\left(4\kappa^{4}(1+\sigma^{2})-\frac{3\eta^{4}}{B_{0}}\right)+8B_{\alpha
0}\eta^{2}\kappa^{2}\tau\right],$ (11c)

and $\mathcal{D}$ is a complicated expression given in Appendix D.2, Eq.
(D25d) of [30] (see also Appendix D).

Relaxing QS through $B_{20}$ and not other components of $|\mathbf{B}|$ is
only a choice, not a requirement. It is, however, not a whimsical choice. For
one, it has the benefit of making the self-consistent $B_{20}$ regular, in the
sense of Fuchs criteria[52]. This follows from $\mathcal{A}\neq 0$, as Eq.
(11a) is proportional to a sum of squares. This local consideration on
$B_{20}$ is not to say that a solution to the equation that satisfies the
condition of periodicity exists. Proving that consistent periodic solutions to
Eq. (10) exist is a more challenging problem. To explore this question, let us
first consider the axisymmetric limit.

In the axisymmetric limit, everything is by definition $\varphi$-independent,
and thus Equation (10) reduces (using the axisymmetric simplification of
$[1+(1/4B_{0})\eta^{4}/\kappa^{4}]\bar{\iota}_{0}=B_{\alpha 0}B_{\theta
20}\eta^{2}/2\kappa^{2}$) to,

$-\bar{\iota}_{0}^{2}\left(\eta^{4}-12B_{0}\kappa^{4}\right)\frac{B_{20}}{B_{0}}-\frac{3}{\eta^{2}\kappa^{2}}(\eta^{4}-4B_{0}\kappa^{4})B_{22}^{C}=\dots.$
(12)

The dots denote second-order independent terms. Solving the equation for
$B_{20}$ (the approach described above[21]) clearly shows ill behaviour for
$\eta^{4}=12\kappa^{4}B_{0}$, seemingly an arbitrary choice of elliptical
shape. In that scenario, $B_{20}$ decouples from the equation, leaving only
$B_{22}^{C}$ to satisfy the equation. The roles of $B_{20}$ and $B_{22}^{C}$
can be reversed to solve Eq. (12) for $B_{22}^{C}$. However, in this case, the
construction fails again at another value of $\eta$, i.e., whenever circular
cross-sections are considered. To avoid excluding this case, solving for
$B_{20}$ appears to be the more convenient choice.

The breakdown of solutions for this ‘critical’ value of $\eta$ persists beyond
axisymmetry. Numerical evidence of this is presented in Fig. 3. We plot
$\Delta B_{20}$ for a fixed axis shape as a function of $\eta$ and
$B_{22}^{C}$, which shows a clear critical value for $\eta$.

Figure 3: Quasisymmetry breaking as a function of $\eta$ and $B_{22}^{C}$ for
an example. Variation in $B_{20}$ ($\Delta B_{20}$) as a function of the
choice of parameters $\\{\eta,B_{22}^{C}\\}$ computed using the code pyQSC.
The dashed lines correspond to the $\eta_{\mathrm{VMEC}}$ (from the global
equilibrium) and $\eta^{\star}$ values (see Tab. 2).

Formally, the appearance of such singularity may be explained through the
Fredholm alternative theorem for the existence of solutions. The singularity
occurs when there is no solution to Eq. (10), which occurs whenever a solution
exists to the adjoint problem. The adjoint problem does not generally have a
solution (and thus, a solution to Eq. (10) exists and is unique), but one may
rigorously prove that there is at least one critical value
$\eta_{\mathrm{crit}}$ for which this is the case. If this critical value is
unique (which numerical evidence suggests to be the case), it follows that
$\eta^{*}<\eta_{\mathrm{crit}}$. Thus, our choice of $\eta^{*}$ at lower order
guarantees the existence of a unique solution to Eq. (10). Details on this may
be found in Appendix E.

The existence of a solution and this critical value $\eta_{\mathrm{crit}}$
have told us little about the influence of $B_{22}^{C}$, and how to choose it.
However, the choice of $B_{22}^{C}$ follows naturally when we try to maximise
the quality of QS. That is, we choose $B_{22}^{C}$ so that it minimises the
deviation of $B_{20}(\varphi)$ from being a constant. This will make our model
represent the ‘most’ quasisymmetric configuration. Because the parameter
$B_{22}^{C}$ appears only in the inhomogeneous term of Eq. (10), the choice of
$B_{22}^{C}$ has no dramatic effect on $B_{20}$. Therefore, we expect the
search of $B_{22}^{C}$ that minimises the QS residual to be smooth. More
quantitatively, we write the dependence of $\mathcal{D}$, Eq. (10), on
$B_{22}^{C}$,

$\mathcal{D}=\frac{3\kappa^{2}}{B_{\alpha
0}B_{0}(l^{\prime})^{2}\eta^{3}}\left[\frac{\bar{\iota}_{0}}{2}\left(B_{\alpha
0}^{2}\frac{\eta^{4}}{\kappa^{4}}+4(\sigma^{2}-1)\right)-\right.\\\
\left.-4\sigma^{\prime}-8\bar{\iota}_{0}\sigma^{2}\right]B_{22}^{C}+\dots,$

which in the limit of large $|B_{22}^{C}|$, makes $B_{20}$ scale with
$B_{22}^{C}$, with a solution of the ‘universal’ form,

$[\mathcal{A}\partial_{\phi}^{2}+\mathcal{B}\partial_{\phi}+\mathcal{C}]\frac{B_{20,\mathrm{univ}}}{B_{0}}+\mathcal{D}_{C}=0,$
(13)

where $B_{20}=B_{22}^{C}B_{20,\mathrm{univ}}$. This is universal in the sense
that each first-order construction has a single solution
$B_{20,\mathrm{univ}}$. This solution provides a measure of the effect of
$B_{22}^{C}$ on symmetry- breaking (see Figure 4). The influence of
$B_{22}^{C}$ decreases as $B_{20,\mathrm{univ}}$ becomes closer to a constant
to vanish for axisymmetry. Away from the neighbourhood of axisymmetry (where
$B_{22}^{C}$ has little effect on QS), this means there must be at least one
local minimum at a finite value of $B_{22}^{C}$. This makes minimising $\Delta
B_{20}(B_{22}^{C})$ modifying $B_{22}^{C}$ a well-posed 1D search problem.
This way of choosing $B_{22}^{C}$ can be shown to be representative of
optimised stellarators in practice (see Figure 4).

Figure 4: Example of $B_{22}^{C}$ choice for precise QS designs. Plots showing
the variation of the QS residual as a function of $B_{22}^{C}$ for the NAE
models of the precise QS designs [44] (top) and their respective universal
solutions $B_{20,\mathrm{univ}}$ (bottom). The plots correspond to the QA
(left) and QH (right). The scatter points represent the values of $B_{22}^{C}$
obtained from the $|\mathbf{B}|$ spectrum of the global solutions.

Although choosing $B_{22}^{C}$ this way constitutes a well-posed problem, we
should not disregard other aspects of the stellarator that $B_{22}^{C}$ also
affects. This includes MHD stability near the axis[53, 43], the shaping of
flux surfaces[43], and as a result, the smallest effective aspect ratio of the
configuration[50]. Characterising the influence of second-order parameters on
geometry in a clear way is non-trivial. Detailed analysis on shaping was
presented as part of [54], which we refer the reader to for a detailed
discussion. However, one can show that in the large $|B_{22}^{C}|$ limit, flux
surfaces become increasingly shaped to limit their achievable aspect ratio.
For this to be lower than 10, we shall limit, quite crudely,
$|B_{22}^{C}|\lesssim 10$ (see Appendix F). This is only a rough estimate, but
it provides a useful domain to perform the search of $B_{22}^{C}$ to minimise
$\Delta B_{20}$. This will leave out any good QS configuration outside the
allowed range of $B_{22}^{C}$. To include some of those cases, one may proceed
a posteriori by relaxing the parameter choices.

##  VI Space of quasisymmetric configurations

Following the arguments in the previous sections, we have a prescription to
complete a second-order near-axis model for every regular axis shape. Every
point in our space of curves represents a stellarator-symmetric, vacuum field
stellarator that tries to be as quasisymmetric as possible while preserving
some minimal requirements on the rotational transform and shaping. To
illustrate the power of the approach, we consider the space of configurations
spanned by magnetic axes described by two Fourier harmonics. We do not keep
$\\{Z_{n}\\}$ explicitly in this space, and instead, for each set
$\\{R_{n}\\}$ we look for the most quasisymmetric solution performing a
Nelder-Mead optimisation [55, Ch. 9.5] on $\\{R_{n}\\}$ under the constraint
that $Z_{n}\sim R_{n}$. This is consistent with the evidence gathered from
optimised QS designs. This way, each point in the $\\{R_{n}\\}$ space
corresponds to an optimal QS axis shape.

This reduced space is adequate to capture common QS designs, as was explicitly
shown in [14], and may be seen in Figure 5. The spaces for a number of
different field periods are presented with the colormap representing the
quality of QS, $\Delta B_{20}$. The spaces were generated using the C++
libraries qsc555See https://github.com/landreman/qsc. The script used to
obtain the main plot in Fig. 5 can be found in the Zenodo repository
associated to this paper. The same may be achieved, albeit slower, using
pyQSC, which was how it was originally done and is also included there. and
gsl for the near-axis calculations and optimisation, respectively. On average,
the evaluation of each point in this space takes less than a second running on
a single CPU in a laptop666An 11th Gen i7-11850H core was used for this
purpose. The main space in Fig. 5 (which is 300x300) took a total of about 14
hrs. Most time is devoted to the optimisation sub-problems at each point
(search for $\eta^{*}$, $B_{22}^{C}$, and $\\{Z_{n}\\}$). Of course, the
construction of such a space is trivially parallelisable.

Figure 5: Quality of quasisymmetry for the two-harmonic quasisymmetric
configuration space. The figure shows $\Delta B_{20}$ (in logarithmic scale)
for the space of configurations spanned by two harmonic magnetic axes for
several field periods $N=2,~{}3,~{}4,~{}5$ in the top figures and $N=4$ for
the lower one. As the text describes, the $\\{Z_{n}\\}$ harmonics have been
optimised at each point of space to minimise the quasisymmetry residual. The
coloured scatter points represent typical quasisymmetric designs in our
reduced space (see legend). To represent them all in the $N=4$ space, the
magnetic axis harmonics are rescaled as
$R_{4n}=R_{nN}(1+n^{2}N^{2})/(1+16n^{2})$ following the insight in [14].
Typical designs lie close to the bands of good quasisymmetry, which show the
power of the approach. The black lines represent phase transition curves for
$R_{n}=Z_{n}$. The dark purple point represents a new QH design, construction
presented in Figure 7. The gap at $R_{4}$ is numerical, as the numerical
evaluation of the space was performed in two separate runs.

The space exhibits two remarkable features. First, the QS phase structure
studied in the context of the magnetic axis makes itself clear (see phase
transitions as solid black lines). Such features, as well as others, remain
largely unchanged as the number of field periods changes. This allows us to
represent the QS designs in Table 2 together in Fig. 5. The second important
feature of this configuration space is the appearance of what we call
quasisymmetric branches. These branches consist of well-distinguished regions
of configuration space with excellent QS. (We leave a more precise definition
for the future, a definition that will be necessary for a more systematic
study of the branches.) Importantly, these branches agree with the location of
typical QS designs.[11, 46, 47, 48, 9, 10, 12, 44]. It proves the predictive
power of the approach and the role of our model as a unifying framework.

With the QS branches identified, we have all the tools from the near-axis
framework to investigate their properties. To illustrate what can be learned
from such an analysis, let us focus on the dominant branch in the QA phase
(the phase that includes the origin) that grows from the origin in the
direction of the lowest harmonic $R_{n}$. Many (if not all) standard QA
designs belong to this branch, irrespective of $N$.

Figure 6: Properties of the configurations in the main QA branch. Properties
of configurations in the main QA branch for $N=2,~{}3,~{}4$ (blue, orange,
yellow respectively) plotted against[14] $s_{N}=R_{N}(1+N^{2})$. The
properties shown are the rotational transform on the axis, the QS residual,
the limit on aspect ratio, the $B_{22}^{C}$ parameter, the average $B_{20}$
for the magnetic well criterion, and the magnetic shear. The latter is defined
so that the change in rotational transform is $\bar{\iota}_{1}$ times the
inverse aspect ratio squared. These are obtained by following the branches in
the 2-harmonic space of Fig. 5. The shaded area indicates a region where the
construction of the NAE model is not well-behaved due to proximity to
axisymmetry (see main text).

Without delving into the origin of the branches, which we leave for future
work, we now describe the properties and trends of this class of
configurations as evident from Fig. 6. These properties should be interpreted
as representative of the class, noting that many features may be changed by
additional tweaking of the axis shape and parameter choice. However, before
looking at these, note that in the region closest to the origin (i.e., near
axisymmetry), the model presents ill behaviour, as it fails to show $\Delta
B_{20}\rightarrow 0$. Such misbehaviour aligns with the uncertainty in the
choice of $B_{22}^{C}$ that arises close to axisymmetry (see $B_{22}^{C}$ in
Fig. 5). The reason for this type of behavior may be deeper than we choose to
pursue here, see [58].

Away from this region, some physically interesting trends are observed. The QS
quality degrades as the QA-QH phase transition is approached and the
rotational transform and the magnetic shear grow. The degradation may be seen
as a result of an increasingly shaped configuration driven by an increasingly
twisted axis, a behaviour predicted in [14]. This also leads to an increase in
the minimum effective aspect ratio $\sqrt{2}/\epsilon_{\mathrm{max}}$, as it
does with the number of field periods. An aspect ratio of roughly $5-7$
appears possible for $N=2$ but increases to $15-25$ for $N=4$. Thus, although
the rotational transform grows with $N$, which is of interest, the limitation
in the compactness (and the quality of QS) restricts the configurations of
interest to the lower $N$ values, as observed in optimisation efforts[10, 44].
Compared to the rest of the QA phase space, the shear, the QS residual, and
$1/\epsilon_{\mathrm{max}}$ are small along the QA branch. This appears to
match the observation that the magnetic shear of QS configurations tends to be
small in practice[44]. The magnetic shear also presents a sign opposite to
that of $\iota_{0}$, and thus the rotational transform profiles are tokamak-
like: the rotational transform decreases towards the edge of the
configuration. The magnetic well criterion[53, 43] (via the sign of $\langle
B_{20}\rangle$) shows that the majority of the QA branch gives rise to a
magnetic hill, and thus is MHD unstable. This aligns with the conclusions
reached in [44], where an additional effort was made to reach a solution with
a magnetic well. Those configurations along the QA branch further away from
the origin will be more easily stabilised by additional tweaking of the
configuration since their magnetic hill is shallow.

This discussion of the QA branch is only a partial account of the full power
of the present approach to understand optimised QA configurations. A similar
effort could be devoted to the branches in the QH phase, but we shall not do
that here. Instead, we content ourselves with a few general observations. We
call the branch the HSX branch as it lies close to the HSX design[9]. The
other branch appears not to include any existing QH design, and thus we shall
refer to it as the new QH branch. The QH space is significantly more sensitive
to the parameter choices than the QA phase, which in particular leads to the
location of the new QH branch in phase space changing under different choices
of $\eta$. If $\eta$ were to be treated as another optimisation parameter, the
branch would become a broader region bounded by large elongation
configurations. This latter branch is particularly interesting. Configurations
along it lack the common bean-shaped cross-section and exhibit a natural
magnetic well. Nevertheless, most remarkably, no existing QS-optimised design
belongs to this class. This illustrates the power of this approach in
exploring QS-optimised configurations.

We present for completeness an example of a stellarator belonging to this
branch, which we call ‘new QH’. We construct a global equilibrium solving the
equilibrium problem inside a fixed outer surface, constructed from a finite
aspect ratio evaluation of the near-axis model (see Fig. 5) following [21].
This form of linking the near-axis and global solutions is not the best, as it
uses the worst-described feature of the near-axis model (the ‘outer’ surface)
as an intermediary. However, it suffices as a first approximation. The
particular configuration in the space of Fig. 5 was chosen to allow for a
reasonable aspect ratio (an aspect ratio of $A\sim 13.5$ given by VMEC, which
has difficulty initialising the solver in more compact scenarios). To
construct the final form of the configuration, some additional refinement of
the axis shape was made within the near-axis framework by allowing for three
small additional harmonic components. The configuration in Fig. 7 shows good
QS behaviour as expected from the near-axis model. This is especially
important, considering that no optimisation has been performed in the space of
global equilibria.

This configuration serves as an example of the approach’s potential. The
complications encountered by the numerical solvers (including the large axis
excursion) could be a reason behind the design efforts missing out
stellarators belonging to the new QH branch. However, the exact reason is hard
to pinpoint. It could have also resulted from other constraints (such as
aspect ratio) or the initial guesses

Other authors have recently found stellarator designs considering optimisation
in the space of near-axis configurations[51]. The result is a multitude of
designs, many of which have previously unseen features. The approach there may
be regarded complementary to that presented here, as this should serve to
structure, guide (e.g., seeding the search to be more exhaustive), and
understand the other approach.

Figure 7: Example of new QH branch configuration. a) Different projections of
the 3D boundary of a configuration from the new QH branch. The colormap
represents contours of constant $|\mathbf{B}|$. The global equilibria were
solved using VMEC given the surface of an NAE construction with an axis given
by the following $\\{R_{n},Z_{n}\\}$ components
$R_{n}=\\{0.426,~{}0.044,~{}-6.36\times 10^{-11},~{}2.85\times
10^{-5},~{}3.89\times 10^{-8}\\}$ and $Z_{n}=\\{0.411,~{}0.043,~{}6.53\times
10^{-5},~{}1.36\times 10^{-5},~{}1.16\times 10^{-5}\\}$. The higher harmonics
have been chosen as small last tweaking choices to minimise further $\Delta
B_{20}$, which within the QS framework, we make $\Delta B_{20}\sim 7\times
10^{-3}$. b) Magnetic field magnitude on the last flux surface of the global
solution. c) Quasisymmetric residual and $\epsilon_{\mathrm{eff}}$ (a measure
of particle transport) [59, 60] as a function of radius in the global solution
showing the quasisymmetric nature of the configuration.

##  VII Conclusion

In this paper, we have presented the construction of a model for
quasisymmetric stellarators based on the near-axis expansion, in which
configurations can be identified with their axis shapes only. Doing so enables
us to represent the space of configurations in a form that inherits the
topological structure of a space of closed regular curves. To concoct such a
model requires a careful choice of parameters that form part of the near-axis
expansion.

In the case of vacuum field, stellarator-symmetric stellarators (although the
extension to more general cases should be straightforward), there are two such
parameters: $\eta$ and $B_{22}^{C}$. The former is chosen to guarantee that
flux surfaces are not extremely elongated, maximising the rotational transform
in quasiaxisymmetric configurations. The choice as presented always exists and
is unique. The $B_{22}^{C}$ parameter is then chosen so that each axis shape
is represented by the ‘most quasisymmetric’ configuration. That way, we can
construct an example space of quasisymmetric stellarators.

The model constructed is a powerful tool that reproduces stellarator designs
optimised for quasisymmetry by other approaches. It does so naturally,
grouping them into families that we call quasisymmetric branches. We show
explicitly how the model may be used to study the properties of one such
branch, namely the main QA branch. We also showcase how new configurations
appear from this approach. The structure of the space and simplicity of the
model opens the door to a fundamental study of quasisymmetric stellarators and
their properties, as well as a practical exploration of designs complementary
to recent attempts[51]. However, there remains significant room for future
work, including a systematic way of identifying branches, initialising
optimisation from these, and a more robust way of connecting the near-axis
construction to global equilibrium solvers.

## Acknowledgements

The authors would like to acknowledge fruitful discussions with M. Landreman,
R. Jorge, R. Nies, S. Buller, and E. Paul. This research was primarily
supported by a grant from the Simons Foundation/SFARI (560651, AB) and DoE
Contract No DE-AC02-09CH11466. ER was also partially supported by the
Charlotte Elizabeth Procter Fellowship at Princeton University.

## Data availability

The data supporting this study’s findings are available at the following
Zenodo repository https://doi.org/10.5281/zenodo.7817884.

##  Appendix A Definition and uniqueness of $\eta^{*}$

The choice of $\eta$ modifies the competition between the different
contributions to the rotational transform on the axis. We learned in the main
text that at large and small $\eta$, the rotational transform
$\bar{\iota}_{0}$ vanishes. Thus, by the mean value theorem under the
assumption of continuity, it must be that the function $\bar{\iota}_{0}(\eta)$
has at least a turning point. Let us see what we can learn about this point.

Consider linearising the Riccati equation, Eq. (6), with respect to $\eta$,

$\frac{\mathrm{d}}{\mathrm{d}\phi}\frac{\delta\sigma}{\delta\eta}=-2\bar{\iota}_{0}\sigma\frac{\delta\sigma}{\delta\eta}-\bar{\iota}_{0}\frac{\eta^{3}}{B_{0}\kappa^{4}}+\frac{\eta
B_{\alpha 0}}{\kappa^{2}}(2\tau+B_{\theta
20})-\frac{\delta\bar{\iota}_{0}}{\delta\eta}\left[1+\sigma^{2}+\frac{1}{4B_{0}}\left(\frac{\eta}{\kappa}\right)^{4}\right].$
(14)

Looking at the extremum
$\delta\bar{\iota}_{0}/\delta\eta\stackrel{{\scriptstyle!}}{{=}}0$ makes Eq.
(14) an ODE on $\delta\sigma/\delta\eta$, the change in the solution $\sigma$
upon infinitesimal change of the parameter $\eta$. Thus,
$\delta\sigma/\delta\eta$ must, as $\sigma$, be periodic for
$\sigma+\delta\sigma$ to be so. Because $\sigma(0)=0$ is kept fixed,
$\delta\sigma/\delta\eta=0$ at $\phi=0$. The linearised ODE can then be solved
by an integration factor,

$\frac{\delta\sigma}{\delta\eta}=e^{-2\bar{\iota}_{0}\int_{0}^{\phi}\sigma\mathrm{d}\phi^{\prime}}\int_{0}^{\phi}e^{2\bar{\iota}_{0}\int_{0}^{\phi^{\prime}}\sigma\mathrm{d}\phi^{\prime\prime}}\frac{2}{\eta}\left[\frac{B_{\alpha
0}}{2}(2\tau+B_{\theta
20})\left(\frac{\eta}{\kappa}\right)^{2}-\frac{\bar{\iota}_{0}}{2B_{0}}\left(\frac{\eta}{\kappa}\right)^{4}\right]\mathrm{d}\phi^{\prime}.$
(15)

Imposing periodicity, using the Riccati equation, Eq. (6), and integrating by
parts (assuming $\bar{\iota}_{0}\neq 0$ and stellarator symmetry to drop the
boundary term), we obtain

$\int_{0}^{2\pi}\left[\sigma^{2}+\frac{1}{4B_{0}}\left(\frac{\eta}{\kappa}\right)^{4}-1\right]E\mathrm{d}\phi=0,$
(16)

where
$E=\exp[{2\bar{\iota}_{0}\int_{0}^{\phi}\sigma\mathrm{d}\phi^{\prime}}]$.

Definition of $\eta^{*}$: the first-order parameter $\eta$ extremises the
rotational transform on the axis of a stellarator symmetric construction iff,

$\int_{0}^{2\pi}\left[\sigma^{2}+\frac{1}{4B_{0}}\left(\frac{\eta}{\kappa}\right)^{4}-1\right]E\mathrm{d}\phi=0,$
(17)

where
$E=\exp[{2\bar{\iota}_{0}\int_{0}^{\phi}\sigma\mathrm{d}\phi^{\prime}}]$.

Although exact, the condition in Eq. (8) is implicit through $\sigma$. This
prevents a closed form of $\eta^{*}$, which has to be found numerically. To
this end, finding bounds on $\eta^{*}$ from Eq. (8) is helpful. Because
$\sigma^{2}$ and $E$ are both positive quantities, in order for Eq. (8) to
hold, the integrand in square brackets must cross zero somewhere. However,
when $\eta\geq(4B_{0})^{1/4}\kappa_{\mathrm{max}}$, the integrand is always
positive, and thus there cannot be a solution. This serves as an upper bound
on $\eta^{*}$. In the small $\eta$ limit (with $\sigma\sim\eta^{2}$), the $-1$
piece dominates, and the equation cannot be satisfied. We may then rigorously
give the interval $\eta^{*}\in(0,\sqrt{2}\kappa_{\mathrm{max}}B_{0}^{1/4}]$,
or with less rigour, change the lower bound to
$\eta^{*}>\sqrt{2}\kappa_{\mathrm{min}}B_{0}^{1/4}$, condition below which
$\sigma^{2}$ in Eq. (8) is the only term that may balance the other two. Given
the differences in curvature between the QA and QH phases (see [14]),
$\eta^{*}$ will tend to be larger in QH configurations.

Although the existence of $\eta^{*}$ is guaranteed by the form of the $\eta$
asymptotes and the mean value theorem, we do not know whether such an extremum
is unique. To prove so, we will investigate the second variation of the
Riccati equation, Eq. (6), and its sign at $\eta=\eta^{*}$. Taking the
variation of Eq. (14) (and looking at the extrema),

$\frac{\mathrm{d}}{\mathrm{d}\phi}\frac{\delta^{2}\sigma}{\delta\eta^{2}}=-2\bar{\iota}_{0}\left[\left(\frac{\delta\sigma}{\delta\eta}\right)^{2}+\sigma\frac{\delta^{2}\sigma}{\delta\eta^{2}}\right]-3\bar{\iota}_{0}\frac{\eta^{2}}{B_{0}\kappa^{4}}+\\\
+\frac{B_{\alpha 0}}{\kappa^{2}}(2\tau+B_{\theta
20})-\frac{\delta^{2}\bar{\iota}_{0}}{\delta\eta^{2}}\left[1+\sigma^{2}+\frac{1}{4B_{0}}\left(\frac{\eta}{\kappa}\right)^{4}\right].$
(18)

The periodicity requirement in $\delta^{2}\sigma/\delta\eta^{2}$ may be
written in the following solvability form: eliminating $\tau$ using the
Riccati $\sigma$-equation, Eq. (6), integrating by parts, using the $\eta^{*}$
condition, Eq. (8), and using the odd parity of $\sigma$ in stellarator
symmetry, we obtain

$\int_{0}^{2\pi}E\left\\{-2\bar{\iota}_{0}\left[\left(\frac{\delta\sigma}{\delta\eta}\right)^{2}+\frac{\eta^{2}}{B_{0}\kappa^{4}}\right]-P\frac{\delta^{2}\bar{\iota}_{0}}{\delta\eta^{2}}\right\\}\mathrm{d}\phi=0,$
(19)

where $P=1+\sigma^{2}+(\eta/\kappa)^{4}/4B_{0}$. The term with the square
bracket has a sign of $-\mathrm{sgn}(\bar{\iota}_{0})$. As $P,E>0$, for the
integral to vanish, it must be the case that the sign of the second variation
of $\bar{\iota}_{0}$ (which is not a function of $\phi$) satisfies,

$\mathrm{sgn}\left(\frac{\delta^{2}\bar{\iota}_{0}}{\delta\eta^{2}}\right)=-\mathrm{sgn}(\bar{\iota}_{0}).$
(20)

As the sign of $\bar{\iota}_{0}$ is set by the sign of the combination
$\int(B_{\alpha 0}/\kappa^{2})(2\tau+B_{\theta 20})\mathrm{d}\phi$ (and thus
cannot change with $\eta$), the extrema of $\bar{\iota}_{0}$ can only be
either maxima or minima (but only one of these), as the sign of
$\delta^{2}\bar{\iota}_{0}/\delta\eta^{2}$ is fixed. Thus, $\eta^{*}$ is
unique.

##  Appendix B Magnetic shear decoupling for $\eta^{*}$

It was shown in [29] that, within the near-axis framework, one could evaluate
magnetic shear (and higher derivatives of rotational transform) as solvability
conditions of first-order periodic ODEs, so-called generalised $\sigma$
equations. Physically, this is reasonable since the rotational transform must
be self-consistently chosen given an average toroidal current profile. Here we
shall focus on magnetic shear, writing
$\bar{\iota}=\bar{\iota}_{0}+\psi\bar{\iota}_{1}+\dots$,
$\mathrm{d}\bar{\iota}/\mathrm{d}\psi=\bar{\iota}_{1}$ [5, Sec. 2.8]. Having
an understanding of this quantity is important, as it affects important
properties of the stellarator, such as ballooning stability [61, Ch. 6.14].

From [29], it is clear that magnetic shear is a third-order quantity in the
near-axis framework (or rather 2.5 order) in a vacuum given by

$\bar{\iota}_{1}=\frac{\int_{0}^{2\pi}\mathrm{d}\phi^{\prime}e^{2\bar{\iota}_{0}\int_{0}^{\phi^{\prime}}\sigma\mathrm{d}\phi^{\prime\prime}}\tilde{\Lambda}_{3}}{\int_{0}^{2\pi}\mathrm{d}\phi^{\prime}e^{2\bar{\iota}_{0}\int_{0}^{\phi^{\prime}}\sigma\mathrm{d}\phi^{\prime\prime}}\left[1+\sigma^{2}+\frac{1}{4B_{0}}\left(\frac{\eta}{\kappa}\right)^{4}\right]}.$
(21)

where,

$\displaystyle\tilde{\Lambda}_{3}=$
$\displaystyle\frac{1}{(Y^{S}_{1,1}){}^{2}}\left\\{2l^{\prime}\left[\tau\left(X^{C}_{3,1}Y^{S}_{1,1}+2X^{C}_{2,2}Y^{S}_{2,2}+X^{C}_{1,1}Y^{S}_{3,1}-\right.\right.\right.$
$\displaystyle\left.\left.\left.X^{S}_{3,1}Y^{C}_{1,1}-2X^{S}_{2,2}Y^{C}_{2,2}\right)+\kappa\left(2X^{C}_{2,2}Z^{S}_{2,2}+X^{C}_{1,1}Z^{S}_{3,1}-\right.\right.\right.$
$\displaystyle\left.\left.\left.-2X^{S}_{2,2}Z^{C}_{2,2}\right)\right]-2\iota_{0}\left(2X^{C}_{2,2}{}^{2}+X^{C}_{1,1}X^{C}_{3,1}+2X^{S}_{2,2}{}^{2}+\right.\right.$
$\displaystyle\left.\left.+2Y^{C}_{2,2}{}^{2}+2Y^{S}_{2,2}{}^{2}+Y^{S}_{1,1}Y^{S}_{3,1}+2Z^{C}_{2,2}{}^{2}+2Z^{S}_{2,2}{}^{2}\right)-\right.$
$\displaystyle\left.-X^{S}_{3,1}X^{C}_{1,1}{}^{\prime}-2X^{S}_{2,2}X^{C}_{2,2}{}^{\prime}+2X^{C}_{2,2}X^{S}_{2,2}{}^{\prime}+X^{C}_{1,1}X^{S}_{3,1}{}^{\prime}\right.$
$\displaystyle\left.-Y^{S}_{3,1}Y^{C}_{1,1}{}^{\prime}-2Y^{S}_{2,2}Y^{C}_{2,2}{}^{\prime}+2Y^{C}_{2,2}Y^{S}_{2,2}{}^{\prime}+Y^{C}_{1,1}Y^{S}_{3,1}{}^{\prime}\right.$
$\displaystyle\left.-2Z^{S}_{2,2}Z^{C}_{2,2}{}^{\prime}+2Z^{C}_{2,2}Z^{S}_{2,2}{}^{\prime}\right\\}.$
(22)

The weighted integral of this quantity $\tilde{\Lambda}_{3}$ drives the
magnetic shear, Eq. (21).

This quantity depends on third-order quantities, for which closed forms may be
found in the literature[22, 21, 28]. Although Eq. (22) is of the third order,
most elements in it may be written exclusively in terms of lower order
quantities, except for $X_{3}$, which explicitly introduces parameters
$B_{31}^{C}$ and $B_{31}^{S}$. Therefore, let us focus on the elements in the
shear that depend on $B_{31}^{C}$.

This component of the magnetic field may be seen as the variation of the
parameter $\eta$ with radius. An approximate variation of ellipticity of flux
surfaces with radius. It is then unsurprising that, given the central role of
$\eta$ in determining the rotational transform on the axis, $B_{31}^{C}$
directly affects the shear. Explicitly, the dependence of
$\tilde{\Lambda}_{3}$ on $B_{31}^{C}$ is,

$\tilde{\Lambda}_{3}=\left[\frac{\bar{\iota}_{0}}{B_{0}\eta}\left(\frac{1}{4B_{0}}\left(\frac{\eta}{\kappa}\right)^{4}-1\right)-\frac{\sigma^{\prime}}{2B_{0}\eta}\right]B_{31}^{C}+\dots,$
(23)

from the expressions for $X_{31}$ (see Appendix F in [29]). Define the factor
$\mathcal{M}$ as the multiplicative factor modulating the contribution of
$B_{31}^{C}$ to the magnetic shear. Using Eq. (23) in Eq. (21), integrating by
parts and assuming stellarator symmetry, we obtain

$\mathcal{M}=\frac{\bar{\iota}_{0}}{B_{0}\eta}\left[1-\frac{2\int_{0}^{2\pi}e^{2\bar{\iota}_{0}\int\sigma\mathrm{d}\phi^{\prime}}\mathrm{d}\phi}{\int_{0}^{2\pi}\left[1+\sigma^{2}+\frac{1}{4B_{0}}\left(\frac{\eta}{\kappa}\right)^{4}\right]e^{2\bar{\iota}_{0}\int\sigma\mathrm{d}\phi^{\prime}}\mathrm{d}\phi}\right].$
(24)

The expression in the large square brackets is a number between $(-1,1)$ sets
a hard upper bound to the effect of the third-order modulation on the shear.
Furthermore, it restricts the rotational transform on the axis not to exceed
$(B_{31}^{C}/B_{0})/\eta$. The minimum of $\mathcal{M}$ with finite (non-zero)
$\eta$ is achieved when the expression in brackets vanishes. Rewriting it into
a single fraction, we get

$\mathcal{M}=\frac{\bar{\iota}_{0}}{B_{0}\eta}\frac{\int_{0}^{2\pi}\left[\sigma^{2}+\frac{1}{4B_{0}}\left(\frac{\eta}{\kappa}\right)^{4}-1\right]E\mathrm{d}\phi}{\int_{0}^{2\pi}\left[1+\sigma^{2}+\frac{1}{4B_{0}}\left(\frac{\eta}{\kappa}\right)^{4}\right]E\mathrm{d}\phi},$
(25)

where $E=\exp[2\bar{\iota}_{0}\int\sigma\mathrm{d}\phi^{\prime}]$. The
numerator is precisely the extremum condition for $\eta^{*}$ in Eq. (8). Thus,
the resilience to $\eta$ presents itself by making the magnetic shear
independent of third-order choices. Thus, $\iota_{1}$ becomes a second-order
quantity upon choosing $\eta^{*}$. This makes the truncated near-axis model,
as constructed in this paper, more complete. Note that $\mathcal{M}$ vanishing
does not equal vanishing of the magnetic shear. There remains a generally non-
zero ‘intrinsic’ contribution from lower-order pieces. However, if
configurations with a reasonably low aspect ratio are sought, then the second-
order shaping should remain small, and so will the magnetic shear.

##  Appendix C Choosing $\eta$ to maximise $L_{\nabla}$

Let us consider briefly in this Appendix some of the intricacies of the
measure proposed in [50] as a guiding principle in the choice of the parameter
$\eta$. The measure is the magnitude of $||\nabla\mathbf{B}||$, using the
Frobenius norm $||M_{ij}||^{2}=\sum_{i,j}(M_{ij})^{2}$ of the
$\nabla\mathbf{B}$ tensor. This provides a characteristic length scale of the
field, approximately that corresponding to the maximum distance at which a
coil may be placed[50]. This measure has been used as a proxy of the radius of
applicability of the near-axis construction, which appears to work when
applied to optimisation within the near-axis framework in practice[51, 62].

From this perspective, it appears to be a natural choice of $\eta$ that
maximises the length scale
$L_{\nabla}\stackrel{{\scriptstyle\cdot}}{{=}}1/||\nabla\mathbf{B}||$. To
investigate this, we write the gradient of $\mathbf{B}$ at first order (in the
notation used in this paper),

$\nabla\mathbf{B}=\nabla\psi\partial_{\psi}\mathbf{B}+\nabla\chi\partial_{\chi}\mathbf{B}+\nabla\phi\partial_{\phi}\mathbf{B}.$
(26)

Using
$\mathbf{B}=\mathcal{J}^{-1}(\partial_{\phi}+\bar{\iota}\partial_{\chi})\mathbf{x}$
and the dual relations[20], we may then expand the expression in the near-axis
fashion,

$\nabla\mathbf{B}\approx\mathcal{J}_{0}^{-1}\partial_{\chi}\mathbf{x}_{1}\times\partial_{\phi}\mathbf{x}_{0}\partial_{\psi}\mathbf{B}_{1}+\mathcal{J}_{0}^{-1}\partial_{\phi}\mathbf{x}_{0}\times\partial_{\psi}\mathbf{x}_{1}\partial_{\chi}\mathbf{B}_{1}+\\\
\mathcal{J}_{0}^{-1}\partial_{\psi}\mathbf{x}_{1}\times\partial_{\chi}\mathbf{x}_{1}\partial_{\phi}\mathbf{B}_{0}.$
(27)

For $\mathbf{x}$ we take Eq. (1), to construct expressions like,

$\partial_{\chi}\mathbf{x}_{1}\times\partial_{\phi}\mathbf{x}_{0}=\frac{\mathrm{d}l}{\mathrm{d}\phi}(\partial_{\chi}Y_{1}\hat{\kappa}-\partial_{\chi}X_{1}\hat{\tau}),$

and

$\partial_{\chi}\mathbf{B}_{1}=\mathcal{J}^{-1}\left[-\kappa
l^{\prime}\partial_{\chi}X_{1}\hat{b}+(\partial_{\chi}Y_{1}^{\prime}-\bar{\iota}Y_{1}-\tau
l^{\prime}\partial_{\chi}X_{1})\hat{\tau}+\right.\\\
\left.(\partial_{\chi}X_{1}^{\prime}-\bar{\iota}X_{1}+\tau
l^{\prime}\partial_{\chi}Y_{1})\right].$

With this,

$\nabla\mathbf{B}=\frac{1}{B_{\alpha
0}B_{0}}\left[\frac{1}{2\sqrt{B_{0}}}\left(Y_{11}^{S}X_{11}^{C}{}^{\prime}+\bar{\iota}_{0}Y_{11}^{C}X_{11}^{C}\right)\hat{\kappa}\hat{\kappa}+\right.\\\
\left.+\left(l^{\prime}\tau-\frac{\bar{\iota}_{0}}{2l^{\prime}}(X_{11}^{C})^{2}\right)\hat{\tau}\hat{\kappa}+\right.\\\
\left.+\left(\frac{Y_{11}^{C}{}^{\prime}Y_{11}^{S}-Y_{11}^{S}{}^{\prime}Y_{11}^{C}+\bar{\iota}_{0}[(Y_{11}^{S})^{2}+(Y_{11}^{C})^{2}]}{2\sqrt{B_{0}}}-l^{\prime}\tau\right)\hat{\kappa}\hat{\tau}+\right.\\\
+\left.\frac{1}{2\sqrt{B_{0}}}(X_{11}^{C}Y_{11}^{S}{}^{\prime}-\bar{\iota}_{0}X_{11}^{C}Y_{11}^{C})\hat{\tau}\hat{\tau}\right]+\frac{\kappa}{\sqrt{B_{0}}}(\hat{b}\hat{\kappa}+\hat{\kappa}\hat{b}).$
(28)

which is consistent with [50, Eq. (3.12)], differences in notation provided.

Eq. (28) may be interpreted as a function of $\eta$ for a fixed axis shape,
much in the same way as we did for $\bar{\iota}_{0}$. As in that case, a basic
understanding of the behaviour of $||\nabla\mathbf{B}||$ can be gained from
the large and small $\eta$ asymptotics. For stellarator symmetry and
$\eta\rightarrow 0$, (taking for ease of notation $B_{0}\sim 1$),

$||\nabla\mathbf{B}||=\sqrt{2\kappa^{2}+\frac{2}{(l^{\prime})^{2}}\left(\frac{\kappa^{\prime}}{\kappa}\right)^{2}+\tau^{2}+(\tau+B_{\theta
20})^{2}},$ (29)

while for large $\eta$,

$||\nabla\mathbf{B}||=\sqrt{2\kappa^{2}+\frac{2}{(l^{\prime})^{2}}\left(\frac{\kappa^{\prime}}{\kappa}-\sigma\bar{\iota}_{0}\right)^{2}+\left(\tau-\frac{\bar{\iota}_{0}}{2l^{\prime}}\left(\frac{\eta}{\kappa}\right)^{2}\right)^{2}+\left(\tau-\frac{\bar{\iota}_{0}}{2l^{\prime}}\left(\frac{\eta}{\kappa}\right)^{2}+B_{\theta
20}\right)^{2}}.$ (30)

These are both generally different constant values, provided $\sigma\neq 0$,
which does not guarantee the existence of an extremum. An example of the
special case is that of axisymmetry, for which the asymptotes are
$||\nabla\mathbf{B}||\sim\sqrt{2+B_{\theta 20}^{2}}$, and the extremum of
$||\nabla\mathbf{B}||$ coincides with $\eta^{*}$, the circular cross-section.

The consequence of this uneven asymptotic behaviour is generally a small
prominence of the extrema, if there is one. Examples of this behaviour are
shown in Figure 8, where the $L_{\nabla}$ metric is compared to other first-
order measures such as the rotational transform and the weighted averaged
elongation for a selection of quasisymmetric designs. This uncertain behaviour
of $||\nabla\mathbf{B}||$, which shows excellent behaviour in the QH cases,
makes this gradient length not as robust a measure as the other choices
proposed (see Table 2). This excellent behaviour in the case of QH
configurations may result from the tendency of such configurations to develop
small-scale features. Like $\mathcal{E}$, $L_{\nabla}$ is also a function of
$\phi$, and differences arise from different considerations.

Figure 8: Comparison of $\eta$ measures. The main measures guiding the choice
of $\eta$ discussed in the text are compared for the same configurations as
Table 2. These include $|\bar{\iota}_{0}|$ whose extremum is $\eta^{*}$
(blue), the weighted elongation whose extremum is $\eta_{\bar{\mathcal{E}}}$
(orange), the mean $L_{\nabla}$ whose maximum is $\eta_{\nabla\mathbf{B}}$
(green), and for comparison the minimum of $L_{\nabla}$ (red). Green and red
correspond to the right-axis scale. The blue point represents $\eta$ from the
global optimised designs, namely $\eta_{\mathrm{VMEC}}$. In many designs, all
considerations provide similar $\eta$ values (see Table 2), but the
shortcoming of $L_{\nabla}$ is also patent (especially in QA examples such as
ESTELL [46]). These also show the similarity in behaviour of
$\langle\mathcal{E}/\kappa^{2}\rangle^{-1}$ and $|\bar{\iota}_{0}|$.

##  Appendix D Second order ODE on $B_{20}$

The seed of overdetermination at the second order comes from the magnetic
field having to satisfy both the QS and equilibrium conditions. The common way
to avoid this overdetermination at the second order is to relax the condition
of quasisymmetry partially. Doing so requires solving a second-order, regular
ODE for $B_{20}$, the second-order, $\theta$ independent change in
$1/|\mathbf{B}|^{2}$. We write this equation explicitly here for a
stellarator-symmetric vacuum field:

$\mathcal{A}\frac{\mathrm{d}^{2}}{\mathrm{d}\phi^{2}}\left(\frac{B_{20}}{B_{0}}\right)+\mathcal{B}\frac{\mathrm{d}}{\mathrm{d}\phi}\left(\frac{B_{20}}{B_{0}}\right)+\mathcal{C}\frac{B_{20}}{B_{0}}+\mathcal{D}=0,$
(31)

where,

$\displaystyle\mathcal{A}=-\frac{B_{\alpha
0}\eta^{2}}{2\kappa^{2}\bar{\iota}_{0}l^{\prime}}\left[1+\frac{4B_{0}\kappa^{4}}{\eta^{4}}(1+\sigma^{2})\right],$
(32a) $\displaystyle\mathcal{B}=\frac{2B_{\alpha
0}\eta^{2}}{\bar{\iota}_{0}l^{\prime}}\frac{\kappa^{\prime}}{\kappa^{3}}-\frac{4l^{\prime}\sigma}{\bar{\iota}_{0}}\tau,$
(32b) $\displaystyle\mathcal{C}=-\frac{l^{\prime}}{2B_{\alpha
0}\eta^{2}\kappa^{2}}\left[\bar{\iota}_{0}\left(4\kappa^{4}(1+\sigma^{2})-\frac{3\eta^{4}}{B_{0}}\right)+8B_{\alpha
0}\eta^{2}\kappa^{2}\tau\right],$ (32c)

and $B_{\alpha 0}^{2}B_{0}=(\mathrm{d}l/\mathrm{d}\phi)^{2}$. The homogeneous
operator only depends on zeroth and first-order considerations. However, the
inhomogeneous term $D$ is given by,

$\mathcal{D}=\frac{\mathcal{Y}_{0}^{S}}{\mathcal{Y}_{0}^{C}}\beta_{R}^{C}+\mathcal{Y}_{1}^{S}\left(\frac{\beta_{R}^{C}}{\mathcal{Y}_{0}^{C}}\right)^{\prime}+\beta_{R}^{S},$
(33)

where,

$\displaystyle\mathcal{Y}_{0}^{C}=\frac{4\bar{\iota}_{0}\kappa}{\eta},$ (34a)
$\displaystyle\mathcal{Y}_{0}^{S}=\frac{4\kappa^{\prime}-8\bar{\iota}_{0}\kappa\sigma}{\eta},$
(34b) $\displaystyle\mathcal{Y}_{1}^{S}=-\frac{4\kappa}{\eta}$ (34c)
$\beta_{R}^{C}=\frac{1}{2l^{\prime}\kappa}\left[4\eta^{2}B_{\alpha
0}^{2}\left(Y_{22}^{C}Z_{22}^{S}+Y_{22}^{S}\left(\tilde{Z}_{20}-Z_{22}^{C}\right)\right)+2B_{\alpha
0}\left(\eta^{2}Y_{22}^{C}\tau l^{\prime}+\right.\right.\\\
\left.\left.+\eta^{2}{X_{22}^{C}}^{\prime}-\kappa
l^{\prime}\tilde{Z}_{20}\left(\eta^{2}-4\kappa\left(X_{22}^{C}-\sigma
X_{22}^{S}\right)\right)-4l^{\prime}\kappa^{2}\sigma
X_{22}^{C}Z_{22}^{S}-4l^{\prime}\kappa^{2}\tilde{X}_{20}Z_{22}^{C}+\right.\right.\\\
\left.\left.+4l^{\prime}\kappa^{2}\sigma\tilde{X}_{20}Z_{22}^{S}+2\eta^{2}\iota_{0}X_{22}^{S}+4l^{\prime}\kappa^{2}\sigma
X_{22}^{S}Z_{22}^{C}+\eta^{2}\kappa
l^{\prime}Z_{22}^{C}-\eta^{2}\tilde{X}_{20}{}^{\prime}\right)+\right.\\\
\left.+4\kappa^{2}l^{\prime}\left(\sigma\tilde{X}_{20}\tau
l^{\prime}-\sigma\tau l^{\prime}X_{22}^{C}-\tau
l^{\prime}X_{22}^{S}+\sigma\tilde{Y}_{22}^{C}{}^{\prime}-2\iota_{0}(Y_{22}^{C}-\sigma
Y_{22}^{S})+\tilde{Y}_{22}^{S}{}^{\prime}\right)\right]$ (34d)
$\beta_{R}^{S}=\frac{1}{2\kappa l^{\prime}}\left[2B_{\alpha
0}\left(\eta^{2}Y_{22}^{S}\tau
l^{\prime}-2\eta^{2}\iota_{0}X_{22}^{C}+4\kappa^{2}\sigma\tilde{Z}_{20}X_{22}^{C}l^{\prime}-\right.\right.\\\
\left.\left.-4\kappa^{2}X_{22}^{C}Z_{22}^{S}-4\kappa^{2}\sigma\tilde{X}_{20}Z_{22}^{C}l^{\prime}-4\kappa^{2}\tilde{X}_{20}Z_{22}^{S}l^{\prime}+\eta^{2}{X_{22}^{S}}^{\prime}+4\kappa^{2}X_{22}^{S}Z_{22}^{C}l^{\prime}+\right.\right.\\\
\left.\left.+4\kappa^{2}\tilde{Z}_{20}X_{22}^{S}l^{\prime}+\eta^{2}\kappa
Z_{22}^{S}l^{\prime}\right)+\kappa l^{\prime}\left(4\kappa\tilde{X}_{20}\tau
l^{\prime}+4\kappa l^{\prime}\tau X_{22}^{C}-4\kappa l^{\prime}\sigma\tau
X_{22}^{S}-4\kappa\tilde{Y}_{22}^{C}{}^{\prime}-8\iota_{0}\kappa(\sigma\tilde{Y}_{22}^{C}+Y_{22}^{S})+\right.\right.\\\
\left.\left.+4\kappa\sigma\tilde{Y}_{22}^{S}{}^{\prime}\right)-4B_{\alpha
0}^{2}\eta^{2}\tilde{Y}_{22}^{C}\tilde{Z}_{20}\right],$ (34e)

The various forms of $X$, $Y$, and $Z$ are needed to complete the expression
for $D$ (remembering that for the expressions we should make $B_{20}=0$ as we
have dealt with these explicitly already). can be found explicitly in [22] in
the simplifying vacuum limit. More explicitly, the forms of $X$ are found in
Eqs. (14)-(15) and Appendix C; the expressions for $Y$ in Eqs. (27)-(28); and
the expressions for $Z$ in Eq. (24) and following. For a systematic way to
obtain these expressions, see [29]. This leaves a contribution to $D$
proportional to $B_{22}^{C}$.

##  Appendix E Solution existence for $B_{20}$

In this Appendix, we consider some existence and uniqueness properties of the
$B_{20}$ equation, Eq. (10). We do so by assessing the second-order
differential equation by applying the Fredholm Alternative theorem.

Fredholm Alternative Theorem Let $\mathbb{L}$ be a linear operator with
adjoint $\mathbb{L}^{\dagger}$. Then exactly one of the following is true:

  * •

The inhomogeneous problem $\mathbb{L}y=f$ has a unique solution $y$.

  * •

The homogeneous adjoint problem $\mathbb{L}^{\dagger}y=0$ has a non-trivial
solution.

In the event of the latter, the inhomogeneous equation has either no solution
or infinitely many. If the solvability condition $\langle y_{0},f\rangle=0$,
then there are an infinite number of solutions.

For application of the Fredholm Alternative theorem, write Eq. (10) as,

$\mathbb{L}\frac{B_{20}}{B_{0}}+\mathcal{D}=[\mathcal{A}\partial_{\phi}^{2}+\mathcal{B}\partial_{\phi}+\mathcal{C}]\frac{B_{20}}{B_{0}}+\mathcal{D}=0,$
(35)

and complete the problem with periodic boundary conditions on $B_{20}$.

It is thus clear that we first need to study the homogeneous operator
$\mathbb{L}$. To construct the adjoint of $\mathbb{L}$, it is convenient to
rewrite Eq. (35) in a self-adjoint form. This can be achieved by writing the
second-order ODE in the form of a Hill equation,

$\psi^{\prime\prime}+Q\psi=-\frac{D}{\mathcal{A}}e^{\int\frac{\mathcal{B}}{2\mathcal{A}}\mathrm{d}\phi},$
(36)

where,

$\displaystyle
Q=-\frac{1}{2}\left(\frac{\mathcal{B}}{\mathcal{A}}\right)^{\prime}-\frac{1}{4}\left(\frac{\mathcal{B}}{\mathcal{A}}\right)^{2}+\frac{C}{\mathcal{A}},$
(37a)
$\displaystyle\psi=B_{20}e^{-\int\left(\frac{\mathcal{B}}{2\mathcal{A}}\right)\mathrm{d}\phi}.$
(37b)

Using the coefficients $\mathcal{A}$ and $\mathcal{B}$, Eqs. (11a)-(11b),

$\frac{\mathcal{B}}{\mathcal{A}}=2\bar{\iota}_{0}\sigma+\left[\ln\left(-\frac{\bar{\beta}_{1}^{C}}{\mathcal{Y}_{1}^{S}}\right)\right]^{\prime},$

where $\mathcal{Y}_{1}^{S}=-4\kappa/\eta$ and $\bar{\beta}_{1}^{C}=B_{\alpha
0}\eta^{2}l^{2}\left[1+4\kappa^{4}(1+\sigma^{2})/\eta^{4}B_{\alpha
0}^{2}\right]/8\kappa^{3}\bar{\iota}_{0}$. The exponential factor that maps
$B_{20}\rightarrow\psi$ is then by construction periodic. Therefore, solving
the $B_{20}$ equation is equivalent to solving Eq. (36) with periodic boundary
conditions on $\psi$. This mapping preserves the even parity of $B_{20}$ in
stellarator symmetry, and thus we may focus on Eq. (36).

As the homogeneous operator of Eq. (36) is self-adjoint, the adjoint problem
for the Fredholm Alternative is $\psi^{\prime\prime}+Q\psi=0$. The problem is
nothing but a periodic, time-independent Schrödinger equation. With this
interpretation, the ‘quantum potential’ is $V=-Q/2$ (formally taking
$\hbar,m=1$). Due to periodicity, the potential can be thought to represent a
‘crystal’ (or periodic lattice) of period $\Delta\phi=2\pi/N$. The structure
of this potential is determined by the axis shape and choice of $\eta$.
Whenever this adjoint periodic Schrödinger equation supports a zero energy
state, then by the Fredholm Alternative, there will be no unique solution to
the original $B_{20}$ equation. This may mean that no solution or an infinite
number of them exists, a distinction in which the second-order choices will
intervene through $\mathcal{D}$.

For a generally shaped potential (see Fig. 9), there is no closed form for the
energy eigenstates of the adjoint. However, the discrete nature of energy
eigenvalues and the very special requirement for a vanishing energy state
suggests that, in general, such an eigenstate will not exist. Thus, a solution
to the direct $B_{20}$ problem will. Only for very particular choices will the
zero energy eigenstate exist, leading to the divergence observed in Fig. 3.
Such special value for a fixed axis shape corresponds to
$\eta_{\mathrm{crit}}$. We shall illustrate these abstract statements with an
example in which the existence or not of a solution is considered.

Figure 9: Example of effective quantum potential for the adjoint problem. This
shows examples of the potential $-Q/2=V$ for (a) a QA with $R=1+0.0144\cos
4\phi$, $Z=0.0144\sin 4\phi$, $N=4$ and $\eta=2.1$ and (b) a QH with
$R=1+0.1\cos 4\phi$, $Z=0.1\sin 4\phi$, $N=4$ and $\eta=2$. The broken line in
(a) shows a quadratic potential centred around the minimum. These are two
examples of what the effective potential wells of the adjoint problem look
like.

For example, the near-axis construction corresponding to Fig. 9, in which the
potential has a simple finite depth and width well in each lattice cell. In
general, more complex features will be present, as also shown in the figure.
Naturally, the first guess to the form of the eigenstates are states
corresponding to single finite depth well, a standard textbook problem [63,
Chap. 2.6]. Of course, the solution will be significantly more involved than
this. For once, in a lattice, wavefunctions are generally not isolated to each
cell, but there is a hopping energy and an overlap between wavefunctions in
nearby lattice sites that lead to the splitting of the energy states.
Secondly, the shapes of the potential wells change the eigenstructure of the
problem. Of course, not having an exact solution, we would like to estimate
the energy spectrum of the problem without depending excessively on the model
used. This is achieved through a variational approach in quantum mechanics.
From the orthogonality of the energy eigenstates and focusing on the ground-
state energy, this can be estimated by considering

$E_{0}\leq\frac{\int\left[(\psi^{\prime})^{2}-Q\psi^{2}\right]\mathrm{d}\phi}{\int\psi^{2}\mathrm{d}\phi},$
(38)

where $\psi$ can be any (here real) function. To get a tighter bound on the
energy, an informed guess of $\psi$ is needed. Leaving some parameters of
$\psi$ as unknowns (say concatenated single-well wavefunctions with free depth
and width parameters), a minimum to the expression may be sought, with the
resulting value being an approximation to the energy state.

Figure 10: Ground-state energy of adjoint problem as a function of $\eta$. The
plot shows the evaluation of the adjoint ground-state energies ($E_{0}$
normalised to $Q$ at the bottom of the potential troughs, $Q_{0}$) as a
function of the $\eta$ parameter for the QA example in Fig. 9. The blue curve
shows the estimate of the ground-state using the wavefunction guess from
$B_{20}$, and the orange shows the finite well variational approach. The
estimated ground-state can be considered the smallest of these two curves. The
crossing with the 0 for the $B_{20}$ estimate is precisely the point at which
the resonance occurs in Fig. 9.

For most values of $\eta$ (see the domain in which $E_{0}<0$ in Fig. 10), the
ground-state lies below the zero energy level, and the $B_{20}$ solution is
unique. This assumes that none of the higher energy states resonate, which is
true given the proximity of the energy $E_{0}$ to 0 (compared to the well
depth) and the parity requirement on the solution, which makes the next energy
level solution lie significantly higher. Although we do not prove it here, we
shall assume that only the ground-state is relevant here. The precise location
of the 0-energy crossing will vary with $\psi_{\mathrm{guess}}$, and thus only
a lower bound may be given for $\eta_{\mathrm{crit}}$. A more precise crossing
value can be obtained by constructing $\psi$ from the solution of $B_{20}$.
This is generally not a great guess for the adjoint problem, except in the
region close to where the existence of the solution starts to break down: the
critical $\eta$ value. If this critical value indicates the second alternative
of the Fredholm Alternative, then in the neighbourhood of this critical
resonance, the solution $B_{20}$ will be dominated by the ground-state
solution of the adjoint problem, and thus the zero crossing will be exact (see
Fig. 10). At this point, the choice of $B_{22}^{C}$ is important, as it
determines whether the inhomogeneous term in Eq. (36) is orthogonal to the
ground-state solution or not, controlling whether there exist none or an
infinite number of solutions, where the infinite family would correspond to
the addition of an arbitrary multiple of the ground-state $\psi_{0}$.

Solutions do exist in general and are unique, except for particular critical
values of $\eta$ for which the ground-state of the adjoint problem has
vanishing energy.

Figure 11: Diagram describing the idea behind the existence proof.
Illustration of the rationale behind the proof of existence of a singular
value $\eta_{\mathrm{crit}}$ to the $B_{20}$ equation.

The existence of such a value has not been proven yet. We do so now by showing
that the energy state of the adjoint homogeneous equation necessarily crosses
the zero. Assuming continuity with $\eta$, it will suffice to show that the
energy level is positive (negative) for small (large) $\eta$ (see Fig. 11).

Let us start with the small $\eta$ limit, for which $\mathcal{A}\sim
O(1/\eta^{5})$, $\mathcal{B}\sim O(1/\eta)$ and $\mathcal{C}\sim O(1/\eta)$,
from Eqs. (11a)-(LABEL:eqn:CB20). To leading order $O(\eta^{4})$ then
$Q\sim-(\mathcal{B}/\mathcal{A})^{\prime}/2+\mathcal{C}/\mathcal{A}$. In this
limit, it is sufficient to show that the estimate of the ground-state energy
is negative. Choose the wavefunction to be constant, $\psi=\psi_{0}$. Then,

$E_{0}\leq\frac{\int_{0}^{2\pi}\left[(\psi^{\prime})^{2}-Q\psi^{2}\right]\mathrm{d}\phi}{\int_{0}^{2\pi}\psi^{2}\mathrm{d}\phi}=-\frac{1}{2\pi}\int_{0}^{2\pi}Q\mathrm{d}\phi.$

Given the asymptotic form of $Q$,

$\int_{0}^{2\pi}Q\mathrm{d}\phi=\int_{0}^{2\pi}\frac{\mathcal{C}}{\mathcal{A}}\mathrm{d}\phi\sim\frac{3\bar{\iota}_{0}^{2}}{B_{\alpha
0}(l^{\prime})^{2}}2\pi\geq 0.$

Thus, for a small (but finite $\eta$), the energy state of the adjoint
homogeneous equation has ground-state energy $E_{0}<0$.

In the large $\eta$ limit, $\mathcal{A}\sim O(\eta^{3})$, $\mathcal{B}\sim
O(\eta^{3})$ and $\mathcal{C}\sim O(1/\eta)$. Thus, the potential
$Q\sim-(\mathcal{B}/\mathcal{A})^{\prime}/2-(\mathcal{B}/\mathcal{A})^{2}/4$
becomes independent of $\mathcal{C}$ allowing one to rewrite the adjoint
problem defining $y=\psi\exp[-\int(\mathcal{B}/2\mathcal{A})\mathrm{d}\phi]$,

$y^{\prime\prime}+\frac{\mathcal{B}}{\mathcal{A}}y^{\prime}\sim 0,$ (39)

which can be solved exactly by
$y\sim\mathbb{C}\int\exp[-\int(\mathcal{B}/\mathcal{A})\mathrm{d}\phi^{\prime}]\mathrm{d}\phi$.
The exponential must be positive; thus, its integral is non-periodic, making
$y$ non-periodic as well. Thus, the only solution is the trivial $y=0$. Due to
the adjoint problem only having a trivial solution, in the large $\eta$ limit,
the solution to the $B_{20}$ equation is once again unique. However, this does
not provide us with the ground-state energy of the adjoint equation. To learn
what that is, write $y^{\prime\prime}+(\mathcal{B}/\mathcal{A})y^{\prime}\sim-
E_{0}y$, where we are explicitly including the energy eigenvalue associated
with the ground-state. We multiply the equation by $y$ and integrate over
$\phi$, so that

$y\left[Ey^{\prime}\right]^{\prime}\sim-E_{0}Ey^{2}\rightarrow
E_{0}\sim\frac{\int(y^{\prime})^{2}E\mathrm{d}\phi}{\int
y^{2}E\mathrm{d}\phi}>0.$

Here $E=\exp[2\bar{\iota}_{0}\int\sigma\mathrm{d}\phi]$. Regardless of the
choice of $y$, the ground-state is necessarily positive. Thus, as the ground-
state energy of the adjoint problem in these two limits has opposite signs, by
the mean-value theorem, it must cross the zero at some value of
$\eta=\eta_{\mathrm{crit}}$. This proof does not guarantee the uniqueness of
this singularity (see Fig. 3), but it guarantees its presence. Given the
numerical evidence, we shall, for the remainder of this Appendix, assume it to
be unique.

All of this has introduced a new special $\eta$ value, $\eta_{\mathrm{crit}}$,
which appears not to bear any relation to the choice of $\eta^{*}$ introduced
in the main text. We would like to relate the two. For a unique
$\eta_{\mathrm{crit}}$, the sign of $E_{0}$ at a given $\eta$ determines its
position relative to $\eta_{\mathrm{crit}}$: $E_{0}<0$ for
$\eta<\eta_{\mathrm{crit}}$ and $E_{0}>0$ for $\eta>\eta_{\mathrm{crit}}$.
Thus, we would like to assess the sign of the ground-state energy associated
with $\eta^{*}$ to situate it to the right or left of $\eta_{\mathrm{crit}}$.

Let us again exploit the Ritz variational form of energy, Eq. (38), in this
case using a trial wavefunction
$\psi\sim\exp\left[\int(\mathcal{B}/2\mathcal{A})\mathrm{d}\phi\right]$. A
wavefunction of this form is suggested by the defining condition of
$\eta^{*}$, Eq. (8). For the ground-state energy threshold, the following are
required,

$\displaystyle\int\left(\frac{\mathcal{B}}{\mathcal{A}}\right)^{2}\psi^{2}\mathrm{d}\phi=4\int(\psi^{\prime})^{2}\mathrm{d}\phi,$
$\displaystyle\int\left(\frac{\mathcal{B}}{\mathcal{A}}\right)^{\prime}\psi^{2}\mathrm{d}\phi=-4\int(\psi^{\prime})^{2}\mathrm{d}\phi,$
$\displaystyle\int\frac{\mathcal{C}}{\mathcal{A}}\psi^{2}\mathrm{d}\phi=\frac{|\bar{\iota}_{0}|}{4B_{\alpha
0}\eta^{*}}\int
e^{2\bar{\iota}_{0}\int\sigma\mathrm{d}\phi^{\prime}}\mathrm{d}\phi.$

For the latter, the forms of $\mathcal{A}$ and $\mathcal{B}$ in Eqs.
(11a)-(11b) were used, as well as integration by parts using Eq. (8). Putting
them together,

$\int\left[(\psi^{\prime})^{2}-Q\psi^{2}\right]\mathrm{d}\phi=-\int\frac{\mathcal{C}}{\mathcal{A}}\psi^{2}\mathrm{d}\phi<0,$
(40)

which thus yields a ground-state energy $E_{0}(\eta^{*})<0$. Thus, it follows
that $\eta^{*}<\eta_{\mathrm{crit}}$. The choice of $\eta^{*}$ thus naturally
avoids the divergence of the $B_{20}$ equation, as will a search in
$0<\eta<\eta^{*}$. This observation may be taken as additional evidence in
favour of the choice of $\eta^{*}$ first order.

##  Appendix F Bounding $B_{22}^{C}$ search

We argued in the main text that to avoid extreme second-order shaping, we
should put some bounds on the allowable $B_{22}^{C}$. Following [43], it is
straightforward to see that both the Shafranov shift as well as the standard
measure of triangularity increase linearly with $B_{22}^{C}$ in the large
$|B_{22}^{C}|$ limit. Furthermore, with this dependence, the construction
becomes unphysical if $|B_{22}^{C}|$ is too large.

To make things more quantitative, consider an up-down symmetric cross-section
in the stellarator as representative of the behaviour of shaping in the
configuration, and focus on the behaviour of its Shafranov shift. This
describes the relative displacement of the centres of cross-sections when
going from one flux surface to the next. If the relative displacement of the
centres of the cross-sections at different $\psi$ is too large, then
eventually flux surfaces will intersect each other (this phenomenon of flux
surface intersection gives rise to the measure $r_{c}$, the largest value of
$\epsilon$ in which the near-axis description is sensible, presented in [50]).
This is unacceptable and thus can be leveraged to estimate a bound on
$B_{22}^{C}$. We reproduce from [54] a simple estimate of when this situation
is reached.

Start for simplicity with a second-order construction in which, besides the
elliptic shape, the second-order shaping has $X_{2}=X_{20}+X_{22}^{C}\cos
2\chi$. Focus then on the centre-line of the cross-section, about which it is
up-down symmetric. Looking at $\chi=0$, $X=\epsilon
X_{11}^{C}+\epsilon^{2}(X_{20}+X_{22}^{C})$. Now the intersection occurs
whenever $\partial_{\epsilon}X=0$, which can be solved for the critical
Shafranov shift,

$X_{20}+X_{22}^{C}=\frac{X_{11}^{C}}{2\epsilon}.$ (41)

Whenever the Shafranov shift exceeds the value on the RHS, then the cross-
sections at the stellarator symmetric point will intersect.

Let us see how a large $|B_{22}^{C}|$ affects the Shafranov shift. To do that,
we need to learn about the behaviour of $X_{20}$ and $X_{22}^{C}$ at $\phi=0$.
These expressions may be found by careful consideration of the symmetry
properties of the various functions involved. As a result,

$X_{20}+X_{22}^{C}\sim-\frac{1}{2\kappa}\frac{B_{22}^{C}}{B_{0}}\left(1+\frac{B_{20}^{\mathrm{univ}}}{B_{0}}\right).$
(42)

Thus, a critical $B_{22}^{C}$,

$\frac{|B_{22}^{C}|}{B_{0}}\sim\frac{\kappa
X_{11}^{C}}{\epsilon}\left(1+\frac{B_{20}^{\mathrm{univ}}}{B_{0}}\right)^{-1}.$
(43)

Crudely, taking $\eta=\kappa X_{11}^{C}\sim 1$ and ignoring the expression in
the bracket, we get $|B_{22}^{C}|/B_{0}\sim 1/\epsilon$. That is, we may
roughly consider it proportional to the ‘aspect ratio,’ $1/\epsilon$, of the
configuration. With a reasonable aspect ratio of $1/\epsilon\sim 10$, we
obtain the $B_{22}^{C}$ limit in the text (which could be relaxed by allowing
for other larger values like, e.g., 20).

It should be clear from this approach that the estimate is but a crude one,
yet nevertheless useful. For instance, we are ignoring
$B_{20}^{\mathrm{univ}}$, which could relax this bound significantly depending
on the situation. In practice, the construction of this crude bound allows us
to perform bounded optimisation. In most of the relevant space (where the most
reasonably shaped, quasisymmetric configurations lie), it has not had much of
an effect on the result (the $B_{22}^{C}$ minimum is well within the interval,
as an example, see Fig. 4).

## References

  * Bernardin, Moses, and Tataronis [1986] M. P. Bernardin, R. W. Moses, and J. A. Tataronis, “Isodynamical (omnigenous) equilibrium in symmetrically confined plasma configurations,” The Physics of Fluids 29, 2605–2611 (1986).
  * Cary and Shasharina [1997] J. R. Cary and S. G. Shasharina, “Omnigenity and quasihelicity in helical plasma confinement systems,” Physics of Plasmas 4, 3323–3333 (1997).
  * Hall and McNamara [1975] L. S. Hall and B. McNamara, “Three‐dimensional equilibrium of the anisotropic, finite‐pressure guiding‐center plasma: Theory of the magnetic plasma,” The Physics of Fluids 18, 552–565 (1975).
  * Landreman and Catto [2012] M. Landreman and P. J. Catto, “Omnigenity as generalized quasisymmetry,” Physics of Plasmas 19, 056103 (2012).
  * Helander [2014] P. Helander, “Theory of plasma confinement in non-axisymmetric magnetic fields,” Reports on Progress in Physics 77, 087001 (2014).
  * Boozer [1983] A. H. Boozer, “Transport and isomorphic equilibria,” The Physics of Fluids 26, 496–499 (1983).
  * Nührenberg and Zille [1988] J. Nührenberg and R. Zille, “Quasi-helically symmetric toroidal stellarators,” Physics Letters A 129, 113 – 117 (1988).
  * Rodríguez, Helander, and Bhattacharjee [2020] E. Rodríguez, P. Helander, and A. Bhattacharjee, “Necessary and sufficient conditions for quasisymmetry,” Physics of Plasmas 27, 062501 (2020).
  * Anderson  _et al._ [1995] F. S. B. Anderson, A. F. Almagri, D. T. Anderson, P. G. Matthews, J. N. Talmadge, and J. L. Shohet, “The Helically Symmetric Experiment (HSX) goals, design and status,” Fusion Technology 27, 273–277 (1995).
  * Zarnstorff  _et al._ [2001] M. Zarnstorff, L. Berry, A. Brooks, E. Fredrickson, G. Fu, S. Hirshman, S. Hudson, L. Ku, E. Lazarus, D. Mikkelsen, _et al._ , “Physics of the compact advanced stellarator ncsx,” Plasma Physics and Controlled Fusion 43, A237 (2001).
  * Najmabadi  _et al._ [2008] F. Najmabadi, A. Raffray, S. Abdel-Khalik, L. Bromberg, L. Crosatti, L. El-Guebaly, P. Garabedian, A. Grossman, D. Henderson, A. Ibrahim, _et al._ , “The ARIES-CS compact stellarator fusion power plant,” Fusion Science and Technology 54, 655–672 (2008).
  * Ku and Boozer [2010] L. Ku and A. Boozer, “New classes of quasi-helically symmetric stellarators,” Nuclear Fusion 51, 013004 (2010).
  * Bader  _et al._ [2019] A. Bader, M. Drevlak, D. T. Anderson, B. J. Faber, C. C. Hegna, K. M. Likin, J. C. Schmitt, and J. N. Talmadge, “Stellarator equilibria with reactor relevant energetic particle losses,” Journal of Plasma Physics 85, 905850508 (2019).
  * Rodríguez, Sengupta, and Bhattacharjee [2022] E. Rodríguez, W. Sengupta, and A. Bhattacharjee, “Phases and phase-transitions in quasisymmetric configuration space,” Plasma Physics and Controlled Fusion 64, 105006 (2022).
  * Note [1] We simplify the picture by assuming that the electrostatic potential shares the QS to leading gyro-order and do not include it in our considerations.
  * Tessarotto  _et al._ [1996] M. Tessarotto, J. L. Johnson, R. B. White, and L. Zheng, “Quasi‐helical magnetohydrodynamic equilibria in the presence of flow,” Physics of Plasmas 3, 2653–2663 (1996).
  * Burby, Kallinikos, and MacKay [2021] J. W. Burby, N. Kallinikos, and R. S. MacKay, “Approximate symmetries of guiding-centre motion,” Journal of Physics A: Mathematical and Theoretical 54, 125202 (2021).
  * Rodríguez, Paul, and Bhattacharjee [2022] E. Rodríguez, E. Paul, and A. Bhattacharjee, “Measures of quasisymmetry for stellarators,” Journal of Plasma Physics 88, 905880109 (2022).
  * Note [2] For a more general form of equilibrium, a formally analogous approach exists in terms of so-called generalised Boozer coordinates, details of which may be found in [rodrigGBC]. As a result, many of the properties of quasisymmetric stellarators in this paper are independent of a particular form of equilibrium.
  * Garren and Boozer [1991a] D. A. Garren and A. H. Boozer, “Magnetic field strength of toroidal plasma equilibria,” Physics of Fluids B: Plasma Physics 3, 2805–2821 (1991a).
  * Landreman, Sengupta, and Plunk [2019] M. Landreman, W. Sengupta, and G. G. Plunk, “Direct construction of optimized stellarator shapes. part 2. numerical quasisymmetric solutions,” Journal of Plasma Physics 85, 905850103 (2019).
  * Rodríguez and Bhattacharjee [2021] E. Rodríguez and A. Bhattacharjee, “Solving the problem of overdetermination of quasisymmetric equilibrium solutions by near-axis expansions. i. generalized force balance,” Physics of Plasmas 28, 012508 (2021).
  * Kruskal and Kulsrud [1958] M. D. Kruskal and R. M. Kulsrud, “Equilibrium of a magnetically confined plasma in a toroid,” The Physics of Fluids 1, 265–274 (1958).
  * Rodríguez, Sengupta, and Bhattacharjee [2021] E. Rodríguez, W. Sengupta, and A. Bhattacharjee, “Generalized boozer coordinates: A natural coordinate system for quasisymmetry,” Physics of Plasmas 28, 092510 (2021).
  * Mercier [1964] C. Mercier, “Equilibrium and stability of a toroidal magnetohydrodynamic system in the neighbourhood of a magnetic axis,” Nuclear Fusion 4, 213 (1964).
  * Solov’ev and Shafranov [1970] L. S. Solov’ev and V. D. Shafranov,  _Reviews of Plasma Physics 5_ (Consultants Bureau, New York - London, 1970).
  * Note [3] For simplicity, we have not normalised $\psi$ respect to the magnetic field on the axis and its curvature as it is often customary[20, 21]. Doing so simply introduces constant rescaling factors in various quantities involved.
  * Garren and Boozer [1991b] D. A. Garren and A. H. Boozer, “Existence of quasihelically symmetric stellarators,” Physics of Fluids B: Plasma Physics 3, 2822–2834 (1991b).
  * Rodríguez, Sengupta, and Bhattacharjee [2022] E. Rodríguez, W. Sengupta, and A. Bhattacharjee, “Weakly quasisymmetric near-axis solutions to all orders,” Physics of Plasmas 29, 012507 (2022).
  * Rodriguez [2022] E. Rodriguez,  _Quasisymmetry_ , Ph.D. thesis (2022).
  * Landreman and Sengupta [2018] M. Landreman and W. Sengupta, “Direct construction of optimized stellarator shapes. part 1. theory in cylindrical coordinates,” Journal of Plasma Physics 84, 905840616 (2018).
  * Plunk, Landreman, and Helander [2019] G. G. Plunk, M. Landreman, and P. Helander, “Direct construction of optimized stellarator shapes. part 3. omnigenity near the magnetic axis,” Journal of Plasma Physics 85, 905850602 (2019).
  * Rodriguez and Plunk [2023] E. Rodriguez and G. G. Plunk, “Higher order theory of quasi-isodynamicity near the magnetic axis of stellarators,” arXiv preprint arXiv:2303.06038 (2023).
  * Oberti and Ricca [2016] C. Oberti and R. L. Ricca, “On torus knots and unknots,” Journal of Knot Theory and Its Ramifications 25, 1650036 (2016).
  * Aicardi [2000] F. Aicardi, “Self-linking of spatial curves without inflections and its applications,” Funct Anal Its Appl 34, 79–85 (2000).
  * Fuller [1999] J. Fuller, Edgar J.,  _The geometric and topological structure of holonomic knots_ (1999) p. 112.
  * Moffatt and Ricca [1992] H. K. Moffatt and R. L. Ricca, “Helicity and the cälugäreanu invariant,” Proceedings of the Royal Society of London. Series A: Mathematical and Physical Sciences 439, 411–429 (1992).
  * Fenchel [1951] W. Fenchel, “On the differential geometry of closed space curves,” Bulletin of the American Mathematical Society 57, 44–54 (1951).
  * Paz-Soldan [2020] C. Paz-Soldan, “Non-planar coil winding angle optimization for compatibility with non-insulated high-temperature superconducting magnets,” Journal of Plasma Physics 86, 815860501 (2020).
  * Lonigro and Zhu [2022] N. Lonigro and C. Zhu, “Stellarator coil design using cubic splines for improved access on the outboard side,” Nuclear Fusion 62, 066009 (2022).
  * Helander and Sigmar [2005] P. Helander and D. Sigmar,  _Collisional Transport in Magnetized Plasmas_ , Cambridge Monographs on Plasma Physics (Cambridge University Press, 2005).
  * Ware  _et al._ [2006] A. S. Ware, D. A. Spong, L. A. Berry, S. P. Hirshman, and J. F. Lyon, “Bootstrap current in quasi-symmetric stellarators,” Fusion Science and Technology 50, 236–244 (2006).
  * Rodriguez [2023a] E. Rodriguez, “Mhd stability and the effects of shaping: a near-axis view for tokamaks and quasisymmetric stellarators,” arXiv preprint arXiv:2302.03359 (2023a).
  * Landreman and Paul [2022] M. Landreman and E. Paul, “Magnetic fields with precise quasisymmetry for plasma confinement,” Physical Review Letters 128, 035001 (2022).
  * Note [4] From the Mercier perspective on rotational transform[5],  $\bar{\iota}_{0}=-\frac{1}{2\pi}\intop\ilimits@_{0}^{L}\frac{\mathop{cosh}\nolimits\bar{\eta}-1}{\mathop{cosh}\nolimits\bar{\eta}}(d^{\prime}+\tau)\mathrm{d}l+\frac{1}{2\pi}\intop\ilimits@_{0}^{L}\tau\mathrm{d}l,$ in the large ellipticity $\bar{\eta}\rightarrow\infty$ limit, $\bar{\iota}_{0}=-d/2\pi$, where $d$ is the angle of rotation of the ellipse with respect to the curvature vector. Because $\sigma\sim\eta^{2}\rightarrow 0$ in this limit, the cross-sections align with the Frenet-Serret frame. Thus, the net rotation $d=0$.
  * Drevlak  _et al._ [2013] M. Drevlak, F. Brochard, P. Helander, J. Kisslinger, M. Mikhailov, C. Nührenberg, J. Nührenberg, and Y. Turkin, “ESTELL: A quasi-toroidally symmetric stellarator,” Contributions to Plasma Physics 53, 459–468 (2013).
  * Garabedian [2008] P. R. Garabedian, “Three-dimensional analysis of tokamaks and stellarators,” Proceedings of the National Academy of Sciences 105, 13716–13719 (2008).
  * Garabedian and McFadden [2009] P. R. Garabedian and G. B. McFadden, “Design of the DEMO fusion reactor following ITER,” Journal of research of the National Institute of Standards and Technology 114, 229 (2009).
  * Hirshman and Whitson [1983] S. P. Hirshman and J. C. Whitson, “Steepest‐descent moment method for three‐dimensional magnetohydrodynamic equilibria,” The Physics of Fluids 26, 3553–3568 (1983).
  * Landreman [2021] M. Landreman, “Figures of merit for stellarators near the magnetic axis,” Journal of Plasma Physics 87, 905870112 (2021).
  * Landreman [2022] M. Landreman, “Mapping the space of quasisymmetric stellarators using optimized near-axis expansion,” Journal of Plasma Physics 88, 905880616 (2022).
  * Bender and Orszag [1999] C. M. Bender and S. A. Orszag,  _Advanced mathematical methods for scientists and engineers: I: Asymptotic methods and perturbation theory_ (Springer, 1999).
  * Landreman and Jorge [2020] M. Landreman and R. Jorge, “Magnetic well and mercier stability of stellarators near the magnetic axis,” Journal of Plasma Physics 86, 905860510 (2020).
  * Rodriguez [2023b] E. Rodriguez, “Mhd stability and the effects of shaping: a near-axis view for tokamaks and quasisymmetric stellarators,” arXiv preprint arXiv:2302.03359 (2023b).
  * Wright, Nocedal  _et al._ [1999] S. Wright, J. Nocedal, _et al._ , “Numerical optimization,” Springer Science 35, 7 (1999).
  * Note [5] See https://github.com/landreman/qsc. The script used to obtain the main plot in Fig. 5 can be found in the Zenodo repository associated to this paper. The same may be achieved, albeit slower, using pyQSC, which was how it was originally done and is also included there.
  * Note [6] An 11th Gen i7-11850H core was used for this purpose. The main space in Fig. 5 (which is 300x300) took a total of about 14 hrs. Most time is devoted to the optimisation sub-problems at each point (search for $\eta^{*}$, $B_{22}^{C}$, and $\\{Z_{n}\\}$). Of course, the construction of such a space is trivially parallelisable.
  * Plunk and Helander [2018] G. Plunk and P. Helander, “Quasi-axisymmetric magnetic fields: weakly non-axisymmetric case in a vacuum,” Journal of Plasma Physics 84, 905840205 (2018).
  * Nemov  _et al._ [1999] V. V. Nemov, S. V. Kasilov, W. Kernbichler, and M. F. Heyn, “Evaluation of $1/\nu$ neoclassical transport in stellarators,” Physics of Plasmas 6, 4622–4632 (1999).
  * Nemov, Kasilov, and Kernbichler [2014] V. V. Nemov, S. V. Kasilov, and W. Kernbichler, “Collisionless high energy particle losses in optimized stellarators calculated in real-space coordinates,” Physics of Plasmas 21, 062501 (2014).
  * Wesson [2011] J. Wesson,  _Tokamaks; 4th ed._ , International series of monographs on physics (Oxford Univ. Press, Oxford, 2011).
  * Jorge  _et al._ [2022] R. Jorge, G. Plunk, M. Drevlak, M. Landreman, J.-F. Lobsien, K. Camacho Mata, and P. Helander, “A single-field-period quasi-isodynamic stellarator,” Journal of Plasma Physics 88, 175880504 (2022).
  * Griffiths [2004] D. J. Griffiths,  _Introduction to Quantum Mechanics (2nd Edition)_ , 2nd ed. (Pearson Prentice Hall, 2004).