# Measures of quasisymmetry for stellarators

E. Rodríguez1,2<EMAIL_ADDRESS>E. Paul1,2  A. Bhattacharjee1,2  1
Department of Astrophysical Sciences, Princeton University, Princeton, NJ,
08543 2 Princeton Plasma Physics Laboratory, Princeton, NJ, 08540

###### Abstract

Quasisymmetric stellarators are an attractive class of optimised magnetic
confinement configurations. The property of quasisymmetry (QS) is in practice
limited to be approximate, and thus the construction requires measures that
quantify the deviation from the exact property. In this paper we study three
measure candidates used in the literature, placing the focus on their origin
and the comparison of their forms. The analysis shows clearly the lack of
universality in these measures. As these metrics do not directly correspond to
any physical property (except when exactly QS), optimisation should employ
additional physical metrics for guidance. Especially close to QS minima,
treating QS metrics through inequality constraints so that additional physics
metrics dominate optimisation is suggested. The impact of different
quasisymmetric measures on optimisation is presented through an example.

##  1  Introduction

The stellarator concept (Spitzer (1958)) has seen a revival in the last
decades as an option to achieve controlled thermonuclear fusion. Stellarators
are three dimensional magnetic field configurations which confine hot plasmas
within. The three-dimensionality of these devices provides the necessary
freedom to avoid some of the limiting features that axisymmetric tokamaks
suffer from: most notably, current-driven instabilities. However, this freedom
can come at a price: the three-dimensionality of the magnetic field might lead
to a loss of confinement that continuous symmetry confers upon the tokamak.

Quasisymmetry (QS) (Boozer (1983); Nührenberg & Zille (1988); Tessarotto et
al. (1996); Helander (2014); Rodríguez et al. (2020); Burby et al. (2020)) is
a property of the field that attempts to preserve in 3D stellarators the good
neoclassical properties of the tokamak. The promise of this property has led
to a search of designs that bear it. Although theoretical work suggests that
in an ideal, static equilibrium with isotropic pressure one may not construct
exact QS solutions (Garren & Boozer (1991); Landreman & Sengupta (2018);
Rodríguez & Bhattacharjee (2021b)), one may design such configurations
approximately. In practice there have been two main approaches to finding
approximately quasisymmetric equilibria. One is to directly construct
approximate solutions by employing an asymptotic expansion near the magnetic
axis (Garren & Boozer (1991); Landreman & Sengupta (2018, 2019); Rodríguez &
Bhattacharjee (2021c)) or near axisymmetry (Plunk & Helander (2018)). The
other, and most commonly employed approach, is to treat the construction of QS
configurations as part of an optimisation problem (Bader et al. (2019);
Henneberg et al. (2019)).

To proceed with optimisation it is necessary to construct cost functions that
penalise deviations from QS. This has been done in a number of contexts and
forms (Bader et al. (2019); Henneberg et al. (2019); Paul et al. (2020)).
However, at the present time there appears to be a lack of comparative studies
of the different forms of cost functions used. The purpose of this paper is to
re-evaluate carefully the origin of some of these commonly-used forms and
investigate their differences, significance and content.

The paper is organised as follows. In Section 2 three different measures of QS
are considered. From these basics, the most apparent formal aspects of the
metrics are discussed. In Section 3 the measures are formally compared on the
same footing. The connection between physical properties associated with QS
and the form of the metrics is then explored. Section IV presents a more
practical perspective on the measures by looking at the universality of QS
measures (how they compare across different near-QS designs) and the
importance of the cost functions for optimisation in practice. We close the
paper with some concluding remarks and suggestions for future optimisation.

##  2  Measures of quasisymmetry

Let us start by defining the concept of quasisymmetry from the fundamental
perspective of single-particle dynamics. The class of magnetic fields that
grants an approximate conserved momentum to the gyrocentre dynamics of charged
particles is called quasisymmetric. In particular, we focus our attention on
weakly quasisymmetric magnetic fields as recently defined (Rodríguez et al.
(2020); Constantin et al. (2021)). Formally, a magnetic field $\mathbf{B}$ is
weakly quasisymmetric iff there exists a symmetry vector field $\mathbf{u}$
such that (Rodríguez et al. (2020)),

$\displaystyle\mathbf{u}\cdot\nabla B=0,$ (1)
$\displaystyle\mathbf{B}\times\mathbf{u}=\nabla\Phi,$ (2)
$\displaystyle\nabla\cdot\mathbf{u}=0,$ (3)

where $\Phi$ labels magnetic flux surfaces that we take to be nested
(Rodríguez & Bhattacharjee (2021a)). We are assuming the electrostatic
potential in the problem to be symmetric.

Although this definition has the benefit of directly emerging from single-
particle considerations, it is difficult to use it practically to assess
whether a field is quasisymmetric or not. Doing so would require solving for
the field $\mathbf{u}$. Instead, we would like to have a condition that only
involves magnetic field quantities directly. This can be achieved in a
straightforward way by constructing $\mathbf{u}$ such that it satisfies Eqs.
(1) and (2), $\mathbf{u}=\nabla\Phi\times\nabla B/(\mathbf{B}\cdot\nabla B)$.
At the singular points where $\mathbf{B}\cdot\nabla B=0$ this construction
breaks down unless $\nabla\Phi\times\nabla B\cdot\mathbf{B}=0$. (Note that
this follows from the conditions (1) and (2).) The vanishing of
$\nabla\Phi\times\nabla B\cdot\mathbf{B}$ when $\mathbf{B}\cdot\nabla B=0$ is
referred to as pseudo-symmetry (Mikhailov et al. (2002); Skovoroda (2005)),
and it guarantees that all of the contours of constant $|\mathbf{B}|$ on a
flux surface are linked to the torus in the same way, avoiding local extrema
on the surfaces (and thus providing well-behaved streamlines of $\mathbf{u}$).
Breaking pseudo-symmetry leads to a diverging orbit width (Rodríguez et al.
(2020)) (i.e., loss of deeply-trapped and barely-trapped particles) at these
$|\mathbf{B}|$ extrema. In a QS solution, the construction of $\mathbf{u}$ is
then well defined because QS implies PS.

With the given construction of $\mathbf{u}$, we may substitute it into Eq. (3)
to obtain

$\frac{\nabla\Phi\times\nabla B\cdot\nabla(\mathbf{B}\cdot\nabla
B)}{(\mathbf{B}\cdot\nabla B)^{2}}=0.$

This yields the triple vector formulation of QS. Namely, if we define the
scalar quantity,

$f_{T}=\nabla\psi\times\nabla B\cdot\nabla(\mathbf{B}\cdot\nabla B),$ (4)

the magnetic field is quasisymmetric iff $f_{T}=0$. (By continuity, this
includes pseudo-symmetry.)

For the purpose of this paper, all magnetic configurations have well-behaved
nested flux surfaces and are non-isodynamic. By isodynamic we mean a field
whose magnitude $|\mathbf{B}|$ is a flux function (Bernardin et al. (1986)).
With these provisions, $f_{T}$ can be evaluated for non-quasisymmetric
configurations and serves as a measure of QS. We have here used $\psi$, the
toroidal magnetic flux, as a label of flux surfaces rather than the arbitrary
flux function $\Phi(\psi)$. This is one of many possible choices, and one
could, for instance, replace $\nabla\psi$ with $\hat{\mathbf{n}}$ (the normal
unit vector to the flux surface) in evaluating (4) on a magnetic surface
within a region of continuously nested surfaces.

The measure $f_{T}$ is a coordinate-independent, scalar measure of the
departure from QS. It is a local quantity, although it also has some global
properties. For instance, $\langle f_{T}\rangle=0$, where
$\langle\dots\rangle$ is the flux-surface average.

We have seen how (4) results from the fundamental definition of QS in Eqs.
(1)-(3). This grants $f_{T}$ the special feature that its application is not
tied to a particular form of equilibria. This could be relevant for
constructing cost-functions in so-called one-shot optimisation schemes
(Akçelik et al. (2006); Dekeyser et al. (2014)), i.e., a formulation in which
finding equilibria is part of the optimisation problem. Although QS may be
explored in the presence of forces more general than isotropic plasma
pressure, $p$, it is standard to adopt the magnetohydrostatic (MHS) assumption
$\mathbf{j}\times\mathbf{B}=\nabla p$ where
$\mathbf{j}=\nabla\times\mathbf{B}$ is the current density. Following standard
practice, we shall do so for the remainder of this work.

A key consequence of this form of equilibria is that
$\mathbf{j}\cdot\nabla\psi=0$. Thus, there exist Boozer coordinates (Boozer
(1981); Helander (2014)) $\\{\psi,\theta,\phi\\}$. In these straight-field
line coordinates the covariant form of the magnetic field is
$\mathbf{B}=I(\psi)\nabla\theta+G(\psi)\nabla\phi+B_{\psi}\nabla\psi$ and the
coordinate Jacobian $\mathcal{J}=(G+\iota I)/B^{2}$ where $\iota$ is the
rotational transform. Equation (4) may then be written explicitly in the form,

$f_{T}=\left[\partial_{\theta}B\partial_{\phi}-\partial_{\phi}B\partial_{\theta}\right](\partial_{\phi}+\iota\partial_{\theta})B.$

The QS condition can be rewritten in the form
$(\partial_{\theta}B)(\partial_{\phi}+\iota\partial_{\theta})(\partial_{\phi}B/\partial_{\theta}B)=0$,
from which it follows that $\partial_{\phi}B/\partial_{\theta}B$ must be a
flux function. This implies $|\mathbf{B}|$ has an explicit symmetry and the
conventional definition of QS then follows (Rodríguez et al. (2021); Boozer
(1983)). We refer to it as the Boozer formulation of QS: a magnetic field is
QS iff the magnetic field magnitude can be written as a function
$B=B(\psi,M\theta-N\phi)$, where $N,M\in\mathbb{Z}$. It is convenient to
define the symmetry helicity $\tilde{\alpha}=N/M$ and a helical angle
$\chi=\theta-\tilde{\alpha}\phi$. (We will generally consider all forms of
symmetry other than quasi-poloidal symmetry.) There is no unique way to
construct a scalar measure that quantifies how close a given magnetic field is
to QS in the presented formulation. It is however customary to define the
minimal measure,

$f_{B}={\sum_{\begin{subarray}{c}n,m\\\
n\neq\tilde{\alpha}m\end{subarray}}}|B_{nm}|^{2}.$ (5)

Here $B_{nm}$ stands for the Fourier components of the magnetic field
magnitude $B$ in Boozer coordinates, and the sum is over the non-symmetric
components. All the mode components are weighted equally in this form
(although this could be modified). In practice, only the lower-order modes
significantly contribute to $f_{B}$, a result that follows from the smoothness
of $|\mathbf{B}|$: an $s$ times differentiable function in $L^{1}$ on a torus
will have Fourier coefficients that generally scale like $\hat{f}(m)\sim
1/|m|^{s}$ (Grafakos (2008)).

Unlike $f_{T}$, the helicity of the symmetry $\tilde{\alpha}$ is required for
the evaluation of $f_{B}$. Such knowledge is needed to avoid the symmetric
modes of $|\mathbf{B}|$ in the summation of (5). In a sense, therefore,
$f_{B}$ measures deviations from a particular form of QS. The possibility to
enforce a particular helicity can be of practical importance when optimising
for QS configurations.

This measure is also different from $f_{T}$ in other two important respects.
First, $f_{B}$ is a single flux surface scalar instead of a local measure.
Second, as the Fourier mode resolution of $|\mathbf{B}|$ needed for Eq. (5) is
in Boozer coordinates, an evaluation of $f_{B}$ requires explicit knowledge of
these coordinates. This prevents application of $f_{B}$ to scenarios requiring
$\mathbf{j}\cdot\nabla\psi\neq 0$, unlike $f_{T}$. The need to compute Boozer
coordinates and Fourier resolve $B$ also imposes an additional numerical
burden (e.g., using existing codes such as BOOZXFORM) (Sanchez et al. (2000);
Paul et al. (2020)).

While $f_{T}$ does not require a Boozer coordinate transformation, it also
does not allow specification of the desired helicity, $\tilde{\alpha}$. It is
thus of interest to explore alternative measures of QS that incorporate both
aspects. We shall now obtain what we call the two-term formulation of QS (Paul
et al. (2020); Helander (2014)).

We begin from Eqs. (1)-(3) again. Defining the inner product
$C=\mathbf{B}\cdot\mathbf{u}/\Phi^{\prime}$ and using Eqs. (1) and (2), we can
write $C\mathbf{B}\cdot\nabla B=\mathbf{B}\cdot\nabla\psi\times\nabla B$. This
expression holds even when $\mathbf{B}\cdot\nabla B=0$. To involve $C$
directly in the remaining Eq. (3), we rewrite it in the form (Rodríguez et al.
(2021)),

$\mathbf{B}\cdot\nabla(\mathbf{B}\cdot\mathbf{u})=\mathbf{j}\cdot\nabla\Phi.$
(6)

Under the assumption of $\mathbf{j}\cdot\nabla\psi=0$, it follows that
$C=C(\psi)$ must be a flux function. We can then write

$f_{C}=\mathbf{B}\cdot\nabla\psi\times\nabla B-C(\psi)\mathbf{B}\cdot\nabla
B,$ (7)

which must vanish for some flux function $C$ for a QS field. This is the two-
term formulation of QS. To avoid having to perform a search for $C$ every time
(7) is evaluated, we need some form for the flux function $C$ in terms of more
recognisable quantities. To do so, we consider the limit of a QS configuration
and adopt Boozer coordinates. In that case, and for $\partial_{\chi}B\neq 0$
(otherwise $f_{C}=0$ for any value of $C$),

$C=\frac{\mathbf{B}\cdot\nabla\psi\times\nabla B}{\mathbf{B}\cdot\nabla
B}=\frac{G+\tilde{\alpha}I}{\iota-\tilde{\alpha}}.$ (8)

(This flux function $C$ is often given the symbol $F$ in the literature
(Helander (2014); Paul et al. (2020))). Now $C$ is expressed in terms of
physical quantities: the Boozer currents $I$ and $G$ and the rotational
transform $\iota$. The currents may be calculated without resorting to Boozer
coordinates directly (see Eq. (16) in (Helander (2014))). Any flux coordinate
system will do, under the assumption of well-defined, nested flux surfaces.

It is clear that the helicity of the symmetry $\tilde{\alpha}$ needs to be
specified, just as one must do to evaluate $f_{B}$. The function $C$ for which
it is needed has important physical meaning (Rodríguez et al. (2020)), as it
is a direct measure of the banana width of trapped particles, $\Lambda$, in
the QS limit. In more detail, $\Lambda|\nabla\psi|\sim\rho_{\parallel}C$,
where $\rho_{\parallel}$ is the Larmor radius associated with the parallel
velocity. The resonance at the surface $\iota=\tilde{\alpha}$ leads to
$C\rightarrow\infty$, unless $\mathbf{B}\cdot\nabla\psi\times\nabla B=0$ on
the surface. This latter condition forces the surface to be isodynamic, which
also prevents the appearance of a current singularity and the potential
opening of a magnetic island (Rodríguez & Bhattacharjee (2021a)). In the
present work we shall avoid such surfaces altogether.

In summary, we have a measure $f_{C}$ that is local like $f_{T}$, for which
the helicity of the symmetry must be prescribed as in $f_{B}$, and for which
no Boozer coordinates are needed. This formulation is also more amenable to
gradient-based optimisation techniques such as adjoint methods (Paul et al.
(2020)).

This constitutes the derivation of the three measures of QS—$f_{B}$, $f_{T}$
and $f_{C}$—that will occupy us in this paper. In Table 1 we summarise the
main formal and construction differences of the three forms.

| Helicity | Boozer | $\mathbf{j}\cdot\nabla\psi=0$ | Info
---|---|---|---|---
$f_{B}$ | I | Y | Y | G
$f_{C}$ | I | N | Y | L
$f_{T}$ | O | N | N | L
Table 1: Basic properties of QS measures. Table summarising the basic formal
features of the QS measures introduced in Section 2. These include: whether
the helicity needs to be specified as an input (I) or not (O), whether Boozer
coordinates are required to evaluate it, whether $\mathbf{j}\cdot\nabla\psi=0$
is a necessary assumption and whether the measure is local (L) or global (G).

##  3  Comparison of cost functions

All three forms introduced in the previous section can be interpreted as
measures of QS: they all share the basic feature that they vanish only for a
configuration that is exactly quasisymmetric. However, each form treats
deviations from QS differently. The focus of this section is to evaluate how
these measures compare with each other formally, and to connect them to some
of the physical properties usually associated with QS.

###  3.1 Formal relation

It is convenient to measure deviations from QS in terms of the asymmetric
components of the magnetic field magnitude. Thus, to understand the forms of
$f_{T}$ and $f_{C}$, it is informative to learn how the different Fourier
modes of $|\mathbf{B}|$ are weighted in each of these metrics. Then a
comparison to $f_{B}$ may be done on an equal footing.

Let us start with the definition of $f_{C}$ and write the magnitude of the
magnetic field as a Fourier series,

$B=\sum_{n,m}B_{nm}e^{i(m\theta-n\phi)},$

where the angular Boozer coordinates $\theta$ and $\phi$ are used for the
Fourier resolution of the field. Using the construction of $C$ in Eq. (8) we
may write after some algebra,

$\displaystyle f_{C}=$
$\displaystyle\frac{B^{2}}{\tilde{\alpha}-\iota}(\partial_{\phi}+\tilde{\alpha}\partial_{\theta})B$
(9) $\displaystyle=$ $\displaystyle
i\frac{B^{2}}{\iota-\tilde{\alpha}}\sum_{n/m\neq\tilde{\alpha}}(n-m\tilde{\alpha})B_{nm}e^{i(m\theta-n\phi)}.$

The first line shows that $f_{C}/B^{2}$, although a local measure, has some
global properties like $f_{T}$ does. In particular, the integral along closed
streamlines of $\mathbf{Q}=\mathcal{J}\nabla\psi\times\nabla\chi$ must vanish,
that is, $\oint(f_{C}/B^{2})\mathrm{d}l/|\mathbf{Q}|=0$. The function
$f_{C}/B^{2}$ is also clearly linear in the Fourier content of $B$, but
interestingly, not all modes are weighted equally. The modes ‘furthest’ from
the symmetry direction are weighted more heavily than those close to it.

Although we have successfully written $f_{C}$ in terms of the Fourier content
of $B$, its form is still far from that of $f_{B}$. To make a closer
comparison to $f_{B}$ we need to construct a single flux-surface scalar from
$f_{C}$. We take Eq. (9) and write it as,

$B_{nm}=-i\frac{\iota-\tilde{\alpha}}{n-m\tilde{\alpha}}\left(\frac{f_{C}}{B^{2}}\right)_{nm},$

where $(\dots)_{nm}$ represents the $(n,m)$ Fourier component of the function
in brackets in Boozer coordinates. Using Parseval’s theorem, we write

$\sum_{n,m}(n-m\tilde{\alpha})^{2}|B_{nm}|^{2}=\frac{(\iota-\tilde{\alpha})^{2}}{(2\pi)^{2}}\int\mathrm{d}\theta\mathrm{d}\phi\left(\frac{f_{C}}{B^{2}}\right)^{2},$

which can be rewritten as,

$\overline{\left(\frac{f_{C}}{B^{2}}\right)^{2}}=\sum_{n,m}\left(\frac{n-m\tilde{\alpha}}{\iota-\tilde{\alpha}}\right)^{2}|B_{nm}|^{2},$
(10)

where $\overline{(\dots)}$ denotes the $(0,0)$ Fourier component of the
expression in brackets. Equation (10) is now in a form close to $f_{B}$. The
scalarised version of $f_{C}$, in this case normalised to $B^{2}$ and averaged
over the surface, has a form similar to $f_{B}$ in that it involves a sum over
the non-symmetric Fourier modes of $B$. The main difference comes through the
different weighting of the modes. The measure $f_{C}$ can be thought of as a
modification of $f_{B}$ in which asymmetric components of the magnetic field
magnitude furthest from QS are more heavily penalised.

The difference in weights implies that $f_{B}$ and $f_{C}$ are actually not
monotonic with respect to one another. Changing the energy content of the
Fourier modes may lead to $f_{B}$ increasing but $f_{C}$ decreasing (and vice-
versa). Having said that, there also is a class of changes for which both
metrics respond in the same way. Given that all terms in the summation of Eqs.
(5) and (10) are positive, reducing any individual mode $|B_{nm}|$ will, for
instance, lead to the reduction of both metrics. This observation has
important implications for the universality of QS measures and optimisation,
to be explored in the next sections.

Let us now see how $f_{T}$ fares in comparison. When expressed in terms of the
Fourier content of the magnetic field strength, the non-linear character of
the triple product formulation is clear. We explicitly write,

$(\mathcal{J}^{2}f_{T})_{nm}=i\sum_{k,l}(\iota l-k)(km-nl)B_{n-k,m-l}B_{k,l},$
(11)

where $\mathcal{J}$ is the Jacobian of Boozer coordinates. The non-linearity
appears in the form of a convolution that mixes modes in a non-trivial way. We
also have one higher derivative of the field strength compared to $f_{B}$ or
$f_{C}$. Scalarising this expression by computing an average like in the case
of $f_{C}$ does not help in bringing this form any closer to $f_{B}$. The
convolution couples modes together, making it unclear whether reducing a
single asymmetric mode will reduce $f_{T}$ as it did for $f_{C}$ and $f_{B}$.

We may try to gain some additional understanding on how $f_{T}$ is affected by
changes in the mode content by considering a system that lies close to QS.
Consider a small deviation $\tilde{B}$ in the field magnitude from a
dominantly QS field. Then a linearisation of $f_{T}$ takes the form,

$(\mathcal{J}^{2}f_{T})_{nm}=i(\tilde{\alpha}m-n)(\iota
m-n)\sum_{l}lB_{\tilde{\alpha}l,l}\tilde{B}_{n-\tilde{\alpha}l,m-l}.$

It is clear from this expression that decreasing $\tilde{B}$ decreases the
magnitude of $f_{T}$. However, even upon linearisation it is not obvious that
reducing a particular Fourier harmonic will guarantee a lower value of
$f_{T}$.

The measure $f_{T}$ is, in this form, distinct from the previous metrics.
However, there is actually a natural relation between the triple vector form
and $f_{C}$ through a simple magnetic differential equation (see Appendix A).
This is reasonable given the common origin of the metrics.

###  3.2 Physical relations

We have explored three forms for QS in this work. These constructions describe
deviations from QS differently. To understand the potential physical meaning
of these differences, we now explore physical phenomena often associated with
QS configurations and their relation to the cost functions.

#### Single-particle dynamics

The concept of QS as presented in this work is built on the single-particle
dynamics. In Section 2, beginning with the most fundamental definition, we
constructed three measures of QS. Thus, it is natural to expect that all these
metrics should have some relation to the single-particle dynamics.

Let us start from the definition of the momentum
$\bar{p}=-\Phi+v_{\parallel}\mathbf{b}\cdot\mathbf{u}$, where $v_{\parallel}$
is the parallel velocity of the particle under consideration (Rodríguez et al.
(2020)). For a QS configuration where $\mathbf{u}$ satisfies Eqs. (1)-(3),
this momentum is approximately conserved to order $O(\rho/L)$, where $\rho$ is
the particle Larmor radius and $L$ a characteristic length of the system. When
the symmetry is broken, we may ask how this momentum evolves in time. To
describe that evolution we need to define the vector field $\mathbf{u}$. We
can always choose this field in a way that Eqs. (1) and (2) are satisfied by
construction. To do so we shall assume the existence of flux surfaces as well
as pseudo-symmetry (Mikhailov et al. (2002); Skovoroda (2005)). The latter
guarantees that the contours of constant $|\mathbf{B}|$ are linked to the
torus avoiding local extrema on the surfaces. These well-behaved contours are
vital to the notion of a continuous symmetry, and necessary to prevent
$\mathbf{b}\cdot\mathbf{u}\rightarrow\infty$. Once $\mathbf{u}$ has been
constructed, one evaluates the dynamics of $\bar{p}$ using the leading order
form of the equations of motion from the Littlejohn Lagrangian (Rodríguez et
al. (2020); Littlejohn (1983)). We may thus write to $O(\rho/L)$,

$\frac{\mathrm{d}}{\mathrm{d}t}\bar{p}=\frac{v_{\parallel}^{2}}{B^{2}}\left[\mathbf{B}\cdot\nabla(\mathbf{B}\cdot\mathbf{u})-\mathbf{j}\cdot\nabla\psi\right]=-\left(\frac{v_{\parallel}}{\mathbf{B}\cdot\nabla
B}\right)^{2}f_{T}.$ (12)

This shows that $f_{T}$ constitutes a measure of the conservation of
$\bar{p}$. Of course, even if $f_{T}\neq 0$ the above could result in
conservation of $\bar{p}$ in the bounce-averaged sense. This averaged
conservation could be seen as the difference between QS and omnigeneity
(Landreman & Catto (2012); Helander (2014)). It is also important to note from
(12) that the cost function $f_{T}$ appears with a factor of
$1/(\mathbf{B}\cdot\nabla B)^{2}$. This factor penalises deviations from QS
close to the extrema of $|\mathbf{B}|$ along field lines. The special
importance of the regions close to the extrema is reminiscent of the
requirements for the confinement of barely and deeply trapped particles in
omnigeneity. The presence of this resonance makes the combination
$f_{T}/(\mathbf{B}\cdot\nabla B)^{2}$ numerically ill-behaved as a modified
cost function. A practical implementation could perhaps be achieved by
regularizing the combination, as has been employed in optimisation for pseudo-
symmetry (Mikhailov et al. (2002)). Although attractive from the single-
particle perspective, and even if the resonance is amended, the measure would
still exhibit sharp gradients due to the heavy weight near extrema of
$|\mathbf{B}|$ along field lines.

The argument related to the dynamics of $\bar{p}$ might appear to some extent
artificial, as $\bar{p}$ really only gains special dynamical meaning in the
limit of QS. We may alternatively look at a more physically relevant property:
the net drift of magnetised particles off flux surfaces. It is convenient to
express the off-surface bounce-averaged drift of particles, $\Delta\psi$, as
$\partial_{\alpha}\mathcal{J}_{\parallel}$, where
$\mathcal{J}_{\parallel}=\oint v_{\parallel}\mathrm{d}l$ is the second
adiabatic invariant (calculated taking the integral along the field line and
between bouncing points) and $\alpha$ is the usual magnetic field-line label
(Helander (2014)). To evaluate this derivative of $\mathcal{J}_{\parallel}$,
consider two nearby lying magnetic field lines on the same magnetic surface at
$\alpha$ and $\alpha+\delta\alpha$ (see Fig. 1). Take the turning points that
define the integral for $\mathcal{J}_{\parallel}$ to be at a field magnitude
$|\mathbf{B}|=B_{t}$. We shall assume that these bounce points exist on both
of these field lines (as we are looking at the $\mathcal{J}_{\parallel}$
associated to a class of particles defined by its bounce points). Then taking
the line integral along the contour in Fig. 1 we compute,

Figure 1: Schematic of the loop integral for
$\partial_{\alpha}\mathcal{J}_{\parallel}$. Schematic of the path integral
followed to compute the field line dependence of the second adiabatic
invariant $\mathcal{J}_{\parallel}$. The contours match field lines and
contours of constant $|\mathbf{B}|=B_{t}$.

$\left(\int_{C_{\alpha}}+\int_{C_{\delta\alpha}}+\int_{C_{B+}}+\int_{C_{B-}}\right)v_{\parallel}\mathbf{b}\cdot\mathrm{d}\mathbf{l}=\delta\alpha\int\nabla\psi\cdot\nabla\times(v_{\parallel}\mathbf{b})\frac{\mathrm{d}l}{B}.$

To obtain the equality, we apply Stokes’ theorem and take the limit of small
$\delta\alpha$, taking the latter integral with respect to length along a
field line, $l$. Assuming that $\mathbf{j}\cdot\nabla\psi=0$ and noting that
the integrals along $C_{B+}$ and $C_{B-}$ vanish as $v_{\parallel}=0$ wherever
$|\mathbf{B}|=B_{t}$, we may rewrite,

$\displaystyle\partial_{\alpha}\mathcal{J}_{\parallel}=-\int\frac{\nabla\psi\times\nabla
B\cdot\mathbf{B}}{\mathbf{B}\cdot\nabla
B}\frac{1}{B^{2}}\left(\frac{E}{v_{\parallel}}+\frac{v_{\parallel}}{2}\right)\mathrm{d}B,$

where $B$ is being used to parametrise the integral along the field line. To
be precise, we should split the integral along the field line into
concatenated pieces joint at $|\mathbf{B}|$ extrema. The vanishing of
$\mathbf{B}\cdot\nabla B$ at the stitching points (and the consequent
divergence of the integrand) is generally not problematic as one can consider
switching back to integrating with respect to $l$ arbitrarily close to them,
with that small piece of the integral constituting a vanishingly small
contribution. We shall drop this distinction for simplicity.

One may write the kinematic expression
$(E/v_{\parallel}+v_{\parallel}/2)/B^{2}=\mathbf{B}\cdot\nabla(v_{\parallel}/B)/(\mathbf{B}\cdot\nabla
B)$. As a result, any flux function multiplying this kinematic expression will
lead to a vanishing integral. Thus, the first fraction in the integrand can be
written as $f_{C}/(\mathbf{B}\cdot\nabla B)$. This shows that $f_{C}$ has a
direct relation to the field-line dependence of the second adiabatic
invariant. The complicated nature of the kinematic factor however obscures
this relation. Proceeding with integration by parts and using Eq. (20),

$\displaystyle\partial_{\alpha}\mathcal{J}_{\parallel}=\left.\frac{\nabla\psi\times\nabla
B\cdot\mathbf{B}}{\mathbf{B}\cdot\nabla
B}\frac{v_{\parallel}}{B}\right|_{v_{\parallel}=0}+\int\frac{f_{T}}{(\mathbf{B}\cdot\nabla
B)^{2}}v_{\parallel}\mathrm{d}l.$ (13)

Thus, we have constructed a simple relation between the bounce-averaged, off-
surface drift $\Delta\psi$ and the QS measure $f_{T}$. The measure comes once
again normalised by a factor $(\mathbf{B}\cdot\nabla B)^{2}$. However, this
$f_{T}$ term is not the only term in (13). For the expression in (13) to be
finite (or vanish as we want for QS), the boundary term in (13) should not
diverge, which requires $\nabla\psi\times\nabla B\cdot\mathbf{B}=0$ wherever
$\mathbf{B}\cdot\nabla B=0$. This condition is precisely the condition of
pseudo-symmetry introduced earlier (Mikhailov et al. (2002); Skovoroda
(2005)). Note that this pseudo-symmetry condition appears separate from the
$f_{T}$ integral in (13). Thus, generally, even when $f_{T}\neq 0$ as in
ominigeneous configurations, the condition should be satisfied if
$\partial_{\alpha}\mathcal{J}_{\parallel}$ vanishes (and deeply trapped and
barely trapped particles are confined). As discussed in Section 2, requiring
$f_{T}=0$ everywhere includes this condition. Guaranteeing the boundary term
vanishes also ensures the existence of bouncing points $B_{t}$ on nearby field
lines, in agreement with the assumptions made in evaluating
$\partial_{\alpha}\mathcal{J}_{\parallel}$.

The discussion above shows that both $f_{T}$ and $f_{C}$ can be regarded as a
natural measure of QS in the single-particle perspective. However, when the
configuration fails to be pseudo-symmetric, the interpretation of $f_{T}$ is
obscured.

The measure $f_{C}$ offers some additional particle dynamics insight. As we
mentioned when constructing $f_{C}$ in Sec. 2, the function $C$ in $f_{C}$ is
a measure of the banana width of bouncing particles in the QS limit. The
expression in (9) inherits the scaling $\propto 1/\bar{\iota}$ where
$\bar{\iota}=\iota-\tilde{\alpha}$ from $C$. Hence, $f_{C}$ also has scaling
related to the orbit width and the bootstrap current, which is proportional to
the banana width ((Helander, 2014)).

Although the metrics bear some important physical content in their relation to
the behaviour of single-particle dynamics, the mapping is not one-to-one. More
sophisticated measures are required for an accurate description of the
dynamics of single particles. Especially important are alpha particles, for
which complicated metrics (see (Bader et al. (2019))) or expensive simulations
are needed.

#### Neoclassical transport

Let us now consider the important physics of neoclassical transport (also
related to enhanced neoclassical viscosity). One of the distinguishing
features of QS is to grant a low level of neoclassical transport to the
magnetic configuration (Boozer (1983)). In particular, it prevents the
detrimental $1/\nu$ transport that tends to spoil confinement in stellarators.
Of course, as we deviate from QS, this property will degrade. We would like to
understand how this degradation is captured (if at all) by the different
metrics of QS presented in the paper.

To study this, we shall introduce first the classical stellarator concept of
helical ripple. We define helical ripple as $\epsilon_{h}\sim B_{nm}/B$, a
measure of the symmetry-breaking magnetic field magnitude. This quantity is
clearly linked to $f_{B}$ in (5). When a single, small asymmetric mode on top
of a QS field is considered, in the $1/\nu$ and $\sqrt{\nu}$ regimes, this
ripple leads to an enhancement of neoclassical losses that scales roughly as
$D\sim\epsilon_{h}^{3/2}$ (Ho & Kulsrud (1987); Nemov et al. (1999); Helander
& Simakov (2008)). Therefore, there is a direct relationship between
$\epsilon_{h}$ and transport. In this single-mode classical stellarator
scenario $f_{B}$ is the ripple (unlike in the more general case). It is then
clear that $f_{B}$ does incorporate some direct information about transport.

This classical view on transport has the drawback of treating all ‘ripples’ on
the same footing. However, a more detailed approach shows that different modes
contribute distinctly (Calvo et al. (2014)). Transport in the $1/\nu$ regime
is driven by bouncing particles near minima of $|\mathbf{B}|$ along field
lines. The enhancement of transport can then come from the introduction of
significant additional wells in which new bouncing particles will live.
However, this is not the only mechanism that can enhance transport. As
particles precess following constant $\mathcal{J}_{\parallel}$ surfaces, a
significant difference in the $|\mathbf{B}|$ profile between field lines on
the same flux surface will lead to an enhancement in the drifting of bouncing
particles off the surface. (Additional effects such as collisional detrapping
may also occur.) These phenomena are quantified in (Calvo et al. (2014)).
Small deviations from QS with helical ripple $\epsilon_{h}$ that do not affect
the $|\mathbf{B}|$ profile significantly lead to a transport enhancement
$\langle\mathbf{\Gamma}\cdot\nabla\psi\rangle\sim\epsilon_{h}^{2}$. (The same
scaling is anticipated for flow damping and generally deviations from
intrinsic ambipolarity.) Now, if the deviation is due to a mode with a large
gradient along $\mathbf{B}$, then the changes lead to an enhancement by a
factor $O(\epsilon_{h})$. When the deviations lead to significant differences
between field lines, then the transport changes significantly, not amenable to
a perturbative treatment. (Aspects of this scenario can be understood from an
analysis of the second adiabatic invariant, studied in the previous section.)

In this more detailed picture it is clear that different modes of
$|\mathbf{B}|$ contribute differently to transport. The enhancement is
significant when the gradients along the field line are large in the sense
$|\mathbf{b}_{0}\cdot\nabla B_{1}|\sim|\mathbf{b}_{0}\cdot\nabla B_{0}|$, with
$B_{1}$ the asymmetric field strength. In a more illuminating form,
$|\mathbf{b}_{0}\cdot\nabla B_{1}|/|\mathbf{b}_{0}\cdot\nabla
B_{0}|\sim\left[(n-\iota
m)/\left(M(\tilde{\alpha}-\iota)\right)\right](\epsilon_{h}/\epsilon_{\mathrm{QS}})$
where we have taken $\epsilon_{\text{QS}}$ to represent the helical ripple
associated to the symmetric part of $B$ with symmetric poloidal mode number
$M$. The weight in the numerator indicates the predominance of modes with
large gradients along the field line, while the denominator serves as a
comparative measure for the the gradient along the field line due to the
symmetric field.

This difference between modes is clearly missed by the construction of
$f_{B}$. However, we have shown that the angle-averaged $(f_{C}/B^{2})^{2}$
provides a weighted version of $f_{B}$ in which higher-order modes are
increasingly penalised $\propto(n-m\tilde{\alpha})^{2}$. Although this unequal
weight provides a distinction between modes, the weight is not of the form
$(n-\iota m)$ as in the transport criterion. It instead weights higher-order
modes according to their similarity with respect to the symmetry direction.
Thus, $f_{C}$ only qualitatively matches the high-order mode weighting of the
transport scaling. The weight that appears in the definition of $f_{C}$ can be
interpreted as a consequence of the smoothness of the function. The
quasisymmetric modes are excluded from the sum through a weight that can be
written as a differential operator: the weight must be zero at
$n/m=\tilde{\alpha}$ and smoothly increasing away from it. It thus appears
natural that the weight in $f_{C}$ is centred about the symmetry direction.
This requirement is avoided in $f_{B}$, where the symmetric discrete modes are
avoided. To more accurately treat the impact on transport, one could imagine a
modified version of $f_{B}$ in which the mode weights are given by $(n-\iota
m)$ and the symmetric components are excluded from the summation. It would,
however, still require Boozer coordinates (unlike $f_{C}$). Interestingly, the
weight $(m\iota-n)$ appears in the linearised version of $f_{T}$ in Sec. 3.1.

So far we have focused on small, isolated deviations from QS. In general,
however, transport calculations are significantly more complex. The presence
of multiple ripples will lead to intricate changes in the trapping and passing
particle behaviour. There exist more sophisticated metrics that attempt to
describe these changes in more detail. A commonly used example is the
$1/\nu$-transport regime effective ripple, $\epsilon_{\mathrm{eff}}^{3/2}$
(Nemov et al. (1999, 2014)). The complexity of this measure is not reflected
in either $f_{B}$ or $f_{C}$. (See the following sections.) Thus, although
$f_{B}$ and $f_{C}$ do contain relevant information regarding transport, they
cannot be taken as accurate measures of it.

#### Current singularities and islands

Although flux surfaces are assumed throughout our analysis, we may also relate
the QS metrics to the potential presence of current singularities and magnetic
islands (Tessarotto et al. (1995); Rodríguez & Bhattacharjee (2021a)). Fourier
modes of $|\mathbf{B}|$ that resonate with the finite rotational transform can
lead to Pfirsch-Schlüter (PS) current singularities. These finite-$\beta$
currents are associated with the potential opening of magnetic islands at the
location of the singularities (Reiman & Boozer (1984); Rodríguez &
Bhattacharjee (2021a)). The PS current singularity arises due to the
periodicity constraint on the magnetic differential equation,

$\displaystyle\mathbf{B}\cdot\nabla\left(\frac{j_{||}}{B}+\frac{p^{\prime}C}{B^{2}}\right)=-2p^{\prime}\frac{f_{C}}{B^{3}},$
(14)

for the parallel current, $j_{||}$. The periodicity condition required for
solvability is,

$\displaystyle p^{\prime}\oint\frac{f_{C}}{B^{4}}\,\mathrm{d}l=0,$ (15)

where $l$ measures length along a field line. We note that $f_{C}$ appears
explicitly. This condition will be satisfied in QS, thus precluding the
existence of PS singularities (except possibly at the conflictive
$\iota=\tilde{\alpha}$ surface, see Sec. 2). Note that solvability may still
be obtained away from QS if $p^{\prime}=0$ or the line integral vanishes.

One may generally show that in MHS, the resonant radial magnetic field due to
the singularities is $B_{\rho nm}\propto(1/B^{2})_{nm}$ (in simplified
geometry), and thus $W_{i}\propto\sqrt{(1/B^{2})_{nm}}$ (Reiman & Boozer
(1984); Bhattacharjee et al. (1995); Rodríguez & Bhattacharjee (2021a)). A
non-QS field will potentially have contributions at many rational surfaces, as
driven by the asymmetric modes. In that sense, any sum over the asymmetric
modes of $B$, such as $f_{B}$ or $f_{C}$, will include some of this
information. Resonances only occur when the rotational transform matching the
asymmetric mode helicity exists within the plasma, and thus a more faithful
metric would only account for the resonant asymmetric modes. Due to the
generally asymmetric geometry of configurations, resonant fields arising from
PS currents may appear at other rational surfaces through toroidal coupling
(Rodríguez & Bhattacharjee (2021a)).

##  4  Quasisymmetry in practice: universality and cost functions

The discussion so far has explored the origin of the three different measures
of QS—$f_{B}$, $f_{C}$ and $f_{T}$—and discussed their formal properties and
related physical content. In this section we explore some more practical
aspects of these measures. First, we will address the question ‘Which
configuration is more quasisymmetric?’ and discuss the universality of QS. At
a second stage, we explore the implications of these different formulations
for the optimisation of stellarators through numerical examples.

###  4.1 Universality of QS

The goal of the metrics $f_{B}$, $f_{C}$ and $f_{T}$ is to quantify the
deviation from QS. In this sense, all these forms would appear to be equally
valid candidates. Given the variety of choice, we would like to understand to
what extent it is meaningful to state that a certain field configuration is
more quasisymmetric than another as measured by these metrics.

As demonstrated in the comparison of the previous section, QS metrics are not
generally monotonic with respect to each other. For example, decreasing
$f_{B}$ does not imply decreasing $f_{C}$. One could, for instance, exchange
the energies of a larger low-order mode to a high-order one without changing
$f_{B}$ but modifying $f_{C}$. The non-linearity of $f_{T}$ makes the lack of
monotonicity somewhat more marked. This lack of monotonicity makes it
difficult to determine magnetic configurations which are ‘closer’ to QS, or
for that matter, if the question is meaningful. As different measures of QS
yield different answers, any comparison must be defined with respect to a
particular standard. The notion thus lacks universality: the quantification of
closeness to QS is in some sense arbitrary, and does not necessarily have a
clear physical meaning attached to it.

Let us consider some practical examples. In order to use the metrics of QS
quantitatively, we first need to normalise and scalarise them. Only then will
we be able to make fair comparisons across different devices. Perhaps
unsurprisingly, we will see that this procedure is not unique, emphasising the
lack of universality.

Let us start with the metric $f_{B}$, which measures the Fourier content of
$B$ in Boozer coordinates. It is customary to normalise it to the total
magnetic field energy. So we define,

$\hat{f}_{B}=\sqrt{\frac{1}{\overline{B^{2}}}f_{B}}.$ (16)

Doing so ensures $\hat{f}_{B}\in[0,1)$. This finite range of values makes
$f_{B}$ convenient, although the requirement $B>0$ excludes the value of 1.
Although this normalisation seems the most natural, it is not unique. For
instance, we could decide to normalise $f_{B}$ to (only) the symmetric part of
the field excluding the flux-averaged part, or we could normalise by the flux-
surface average of $B^{2}$, $\langle B^{2}\rangle$. For the purpose of this
discussion, we will take the natural and more conventional form (16).

The choice of a ‘natural’ normalisation is unclear for $f_{C}$ and $f_{T}$.
Fundamentally, our normalisation should give a dimensionless measure. This
principle is, of course, not specific enough. To constrain the construction of
the measure further, we shall try making their Fourier forms look closest to
$f_{B}$. This is a convenient choice. Let us start with $f_{C}$, which
dimensionally scales as $f_{C}\sim B^{3}$. In addition to this dimensional
factor, motivated by the form of Eq. (9), we also choose to normalise by
$1/(\tilde{\alpha}-\iota)$ to bring the measure closer to the form of $f_{B}$.
One could argue for adding additional dimensionless factors such as the
inverse aspect ratio, $\epsilon$, which characterises the magnitude of
$C\mathbf{B}\cdot\nabla B$ for $\chi=\theta$, but for simplicity we shall not
do so here. In view of the normalisation of $f_{B}$, we normalise and
scalarise $f_{C}$ as follows,

$\hat{f}_{C}=\sqrt{\frac{\langle f_{C}^{2}\rangle}{\langle
B^{2}\rangle^{3}}|\tilde{\alpha}-\iota|^{2}}.$ (17)

This leads to a form that is close to but not identical to Eq. (9).

The normalisation of $f_{T}$ is inherently more complex. The triple vector
product introduces explicit length scales which must be non-dimentionalised.
Once again, we look at the form of $f_{T}$ in (10). It is then natural to
normalise and scalarise $f_{T}$ using,

$\hat{f}_{T}=\sqrt{\langle
f_{T}^{2}\rangle\frac{R^{4}}{\langle{B^{2}}\rangle^{4}}}.$ (18)

In obtaining this normalisation we have used the scaling $\mathcal{J}\sim
R/B$, where $R$ is the major radius. Once again, one could argue that the
choice is, to some extent, arbitrary. It is interesting to note that the
combination $f_{T}/(\mathbf{B}\cdot\nabla B)^{2}$, recurrent in the
description of single-particle dynamics in the previous section is by
construction dimensionless. Thus, it appears to be a natural choice. However,
we do not consider it here due to the added complication of the resonant
denominators and the associated steep gradients.

Now that we have chosen specific normalisations for the QS metrics, we may
compare their values in stellarator equilibria. We consider two scenarios.
First, we will focus on the comparison of the metrics in configurations that
have been designed to be approximately quasisymmetric. Then we will study how
these measures differ when the system has not necessarily been optimised to be
quasisymmetric.

Figure 2: Comparison of defined cost functions in quasisymmetric designs.
Plots in logarithmic scale including a comparison of the cost functions
$\hat{f}_{B}$, $\hat{f}_{C}$ and $\hat{f}_{T}$ for a number of quasisymmetric
devices. The scatter points for each device indicate values at a number of
equally spaced magnetic flux surfaces.

We present in Figure 2 a comparison of the metrics in six quasisymmetric
configurations:111 https://github.com/landreman/vmec_equilibria.git GAR
(Garabedian (2008); Garabedian & McFadden (2009)), HSX (Anderson et al.
(1995)), NCSX (Zarnstorff et al. (2001)), ESTELL (Drevlak et al. (2013)),
QHS48 (Ku & Boozer (2010)) and ARIESCS (Najmabadi et al. (2008)). It is clear
that, within a given device, the correlation between cost functions is quite
good. All measures qualitatively show the same behaviour: for each
configuration, QS is better realised closer to the axis than away from it.
(This trend seems to be preserved except in regions which are very close to
the axis.) However, the same cost function correspondence does not seem to
hold if we make comparisons between different configurations.

To quantify this lack of ordering across devices, we compute the Spearman
correlation coefficient, which considers correlation only taking order into
account. We show this in Fig. 3. The high correlation between cost functions
within a particular device (which in all cases was $\sim 0.9-1$) clearly does
not hold across devices. Although $\hat{f}_{C}$ and $\hat{f}_{T}$ retain a
high degree of correlation ($\sim 0.9$), the conventional $\hat{f}_{B}$
measure seems largely uncorrelated to the others. (Of course, in this inter-
configuration comparison the normalisation factors chosen are particularly
important.)

Figure 3: Correlation between edge values of cost functions and
$\epsilon_{\mathrm{eff}}$ of different quasisymmetric configurations. Diagram
showing the Spearman correlation between the edge values of the different cost
functions and $\epsilon_{\mathrm{eff}}$ of the configurations in Fig. 2. The
colour represents the coefficient (as does the diameter of the coloured
circles).

The loss of correlation across machines emphasises the main point on
universality from a practical perspective. Depending on the form of cost
function taken, one could say either A or B is more quasisymmetric. The
statement has only a relative meaning at best for these optimised
configurations. Without a physical basis to the measures it is difficult to
argue in favour of one definition over the other. For a more meaningful
statement, one should instead compare these functions through some physical
property. Ultimately we are interested in finding configurations that are
close to QS because we seek some of the physical properties that QS confers
upon magnetic configurations. In Figure 3 we compare $\epsilon_{\mathrm{eff}}$
(which also contains some arbitrary normalisation factors in its definition
(Nemov et al. (1999))) to the three cost functions for the devices of Fig. 2.
Other physics measures such as alpha particle confinement could also be suited
for such comparison (Bader et al. (2019)). There exists some correlation
between the transport and the normalised values of $f_{C}$ and $f_{T}$. It is
especially remarkable the seemingly little correlation between $\hat{f}_{B}$,
which is widely used for optimisation, and $\epsilon_{\mathrm{eff}}$. Overall,
for these ‘quasisymmetric’ configurations, the QS cost functions appear not to
be good indicators of their transport levels. Similar behaviour has been
observed in previous studies of quasisymmetric configurations (Martin &
Landreman (2020)).

In summary, when comparing quasisymmetric configurations, there is no
universal measure of QS. The comparison depends on the metric used, and these
differ from each another. To make a comparison physically meaningful,
especially close to optimised configurations, we are in need of resorting to
direct physical measures of the configuration (e.g., neoclassical transport or
alpha particle confinement).

###  4.2 Optimising for QS

The differences in the formal structure of the cost functions in conjunction
with the difference in behaviour close to QS suggests that the different
measures of QS will behave quite differently as optimisation cost functions.

From a formal perspective, we would expect to find $f_{B}$ and $f_{C}$ to
perform similarly, while the non-linear $f_{T}$ to differ significantly. We
now present some examples to illustrate the differences in the optimisation.
We consider an example optimisation problem in which QS is targeted throughout
the plasma volume. The optimisation parameters correspond to the Fourier modes
describing the outermost surface,

$\displaystyle R$ $\displaystyle=\sum_{nm}R_{nm}\cos(m\vartheta-
nN_{p}\phi_{c})$ (19) $\displaystyle Z$
$\displaystyle=\sum_{nm}Z_{nm}\sin(m\vartheta-nN_{p}\phi_{c}),$

where $N_{p}$ is the number of field periods, $\phi_{c}$ is the cylindrical
angle, and $\vartheta$ is a general poloidal angle. The configuration
optimised is based on a quasi-helically symmetric configuration222The basis of
the optimisation is the 34DOF_varyBoundary_targetQuasisymmetry scenario in
https://github.com/landreman/stellopt_scenarios. with $N_{p}=4$ field periods.
We shall investigate two aspects of the problem. First, we consider slices of
parameter space with the purpose of observing the differences in QS measures
in practice beyond the region close to the minima. This should give a partial
idea of how the problem changes as a result of changing the form of the
optimisation cost function. Then, we perform various optimisations using
STELLOPT (Lazerson et al. (2020)) with implementations of the three different
forms of the cost functions.

Figure 4: Parameter space spanned by $n=2,m=1$ surface mode for $\hat{f}_{B}$,
$\hat{f}_{C}$, $\hat{f}_{T}$ and $\epsilon_{\mathrm{eff}}$. The plots show,
clockwise, the cost functions $\hat{f}_{B}$, $\hat{f}_{C}$,
$\epsilon_{\mathrm{eff}}$ and $\hat{f}_{T}$ as a function of the $(2,1)$
surface modes. The parameter scan is performed using the map option of
STELLOPT around the stellarator design in the footnote of Page 15. The minima
for each of these is found (clockwise) at: $(-0.031,0.031)$, $(-0.031,0.031)$,
$(-0.037,0.045)$ and $(-0.037,0.037)$.

Let us start with the exploration of parameter space. Figure 4 presents the
evaluation of the QS metrics and $\epsilon_{\mathrm{eff}}$ as a function of
two parameters characterising the surface (the $n=2,~{}m=1$ modes of the
boundary (19)). The metrics are evaluated at six equally-spaced magnetic flux
surfaces, and their contributions are added together.

The objective landscape varies depending on the metrics used. Qualitatively,
though, all the measures show a similar behaviour. The metrics share a region
of reduced cost, which seem to roughly coincide. In a global sense, the
metrics seem not to disagree with each other. This might appear as a surprise
given the observations of the QS configurations and some of the formal
discussion above. We should, however, remember that there exist certain
changes that affect all metrics favourably. Away from a QS configuration it
thus seems that there are directions in which one may reduce many of these
modes such that all metrics decrease. This is supported by the parameter space
representations in Fig. 4. However, note that the similarity in global trends
does not preclude differences in local trends.

First, the relative variation of the cost functions in space is significantly
different. Gradients seem to be relatively strong for the highly non-linear
measures $\hat{f}_{T}$ and $\epsilon_{\mathrm{eff}}$, while the other two
metrics behave more smoothly. This distinction in the gradients also seems to
leave a region close to the minima much shallower for $\hat{f}_{T}$ and
$\epsilon_{\mathrm{eff}}$, comparatively, which may lead to difficulties in
accessing the exact minima. Secondly, beyond these differences in magnitude,
there are also differences in the shape and location of the minima. The region
near a minimum can be more or less elongated depending on the measure, and
importantly, the minimiser is located at different values. In this particular
case, minimisers vary by up to $\sim 20-50$% in the boundary coefficients
between different measures (although $\hat{f}_{B}$ and $\hat{f}_{C}$ have, in
this case, roughly the same minimiser). In brief, although far from a minimum
the various metrics may exhibit a qualitatively similar behaviour, close to
the minimum there are significant differences. These latter differences are
consistent with the discussion on the lack of universality in Sec. 4.1.

Of course, this illustration of the optimisation space is only one slice of a
multidimensional manifold. The large number of degrees of freedom means that
the parameter space slices will change as other parameters do. In general,
this makes optimisation a highly complex process in which the interaction
between different dimensions may amplify the differences between metrics. In
addition, the algorithm employed to move around in optimisation space will
also affect the result of a given optimisation. A discussion on how the
optimisation is likely to perform is therefore extremely difficult, and we
instead present an illustrative example. With that being said, from our
discussion we expect the optimisation of $\hat{f}_{B}$ and $\hat{f}_{C}$ to be
similar, while the non-linearity of $\hat{f}_{T}$ to change the optimisation
significantly. We hypothesise here that there is a rough relation between
changes in boundary shape and the Fourier content of $B$.

We present the resulting cross sections and values associated with the various
metrics from optimisations in which QS is targeted in a volume using the
different forms of QS metrics as cost functions in Fig. 5 and Tab. 2. We start
from a configuration that is deliberately chosen to be some distance away in
parameter space from the QS configuration found in the footnote of Page 15.
The optimisation uses a BFGS gradient-based method (as part of the suite
MANGO333 https://github.com/hiddenSymmetries/mango.) in which there is a
single scalarised cost function with 32 degrees of freedom in the form of
Fourier components (from 0 through 3, keeping the major radius $R_{00}$ fixed)
defining the surface as in (19). The process does not penalise any additional
property such as the commonly-targeted rotational transform or aspect ratio.
By studying a simplified problem, we hope to gain more insight on the
differences between the various QS metrics. The optimisation is performed for
several values of the finite-difference step size to ensure convergence.

Comparing cross-sections in Fig. 5 and the values in Tab. 2 we observe that
the resulting optima are quite different from one another. This emphasises the
points discussed in this paper: that the precise form of the cost function has
important qualitative implications. The optimisation using the standard form
of $\hat{f}_{B}$ seems, in this case, to perform best in allowing the
configuration to reach a state of minimal transport and QS measures. The
optimisation by $\hat{f}_{C}$ is not far away, and a comparison of the cross
sections shows similar shaping features. Interestingly, optimising for
$\hat{f}_{C}$ does not seem to yield the smallest value for $\hat{f}_{C}$
(which $\hat{f}_{B}$ optimisation does). This suggests that the optimisation
for $\hat{f}_{C}$ gets stuck in some effective local minimum.

In light of the comparison made in Sec. 3, and with the hypothesis that the
surface variations are roughly related to variations in $B$, we can think of
optimisation as a process that rearranges energy in the space of Fourier modes
of $|\mathbf{B}|$. Thinking in these terms may shed some light on
understanding the differences in optimisation. In this rearrangement of energy
in Fourier space, $\hat{f}_{B}$ contains a large space of degenerate
interchanges. Energy can be freely exchanged between asymmetric modes of $B$,
regardless of their mode numbers. In that sense, reducing $\hat{f}_{B}$ means,
necessarily, reducing some asymmetric mode. On the other hand, the cost
function $\hat{f}_{C}$ distinguishes between modes. Different weights at
different mode numbers lifts part of the degeneracy in the exchange of mode-
energy that occurs for $\hat{f}_{B}$. This distinction effectively excludes
certain regions of the mode space from being used for the restacking of the
mode-energy. In addition to decreasing the magnitude of an asymmetric mode,
moving energy to modes closer to the symmetry direction will decrease
$\hat{f}_{C}$. This could be a mixed blessing. Having this additional
mechanism could be a way to avoid certain local minima by adding a new
direction of descent. However, it could also present additional local minima
if the optimisation prematurely reduces modes very close to the symmetry
direction.

Relating changes in the mode content to changes in the boundary, we expect
these two metrics to facilitate different ways of modifying the boundary. It
is not clear which metric will lead to improved optimisation. We have here
presented a single numerical exercise in which $\hat{f}_{B}$ appeared to
perform better, but the performance of these metrics will likely depend on the
particulars of the optimisation problem as well as other aspects such as the
optimisation algorithm used.

The optimisation of $\hat{f}_{T}$ in the presented example exhibits the worst
performance of all. The cross sections show large differences with respect to
the other configurations, suggesting that optimisation has stopped at a
distinct local minimum. That the measure $\hat{f}_{T}$ is lower in the
$\hat{f}_{B}$-optimised case serves as evidence that the optimisation of
$\hat{f}_{T}$ teriminated at a different local minimum. The presence of
nonlinearities in $f_{T}$ are perhaps responsible for this additional minimum.
In the light of the Fourier space energy picture, $\hat{f}_{T}$ not only
depends on the magnitude and the location of the mode energy, but also on the
relative location of modes. This nonlinear coupling complicates the problem,
potentially leading to additional local minima. Furthermore, the nonlinearity
of $\hat{f}_{T}$ results in an objective landscape which appears flatter near
the minimum (see Fig. 4) and may difficult attaining an accurate optimum.
Although this qualitative picture fits with observations made for this
optimisation example, we remark that the behaviour is likely to vary between
problems.

In the example presented, the aspect ratio of the configuration does not seem
to change significantly between and along optimisations (values lie within
$\sim 15\%$ of the initial point). This seems counterintuitive, as from a
near-axis perspective one would expect a larger aspect ratio to allow for a
more quasisymmetric solution. This follows from the possibility of satisfying
QS conditions sufficiently close to the magnetic axis, as the relevant
expansion parameter for near-axis scales as $\sqrt{\psi/B}/R$ (Landreman &
Sengupta (2019)). Since we are fixing the total flux, decreasing the minor
radius of the toroidal configuration can increase the differential flux over
which the QS error is small. In practice, however, the optimiser has
difficulty in finding these large-aspect-ratio minima. The difficulty to
change the aspect ratio without sacrificing QS along the way and the choice of
parameterisation (Henneberg et al. (2021)) (19) might provide an explanation
for this behaviour.

Figure 5: Comparison of cross sections of optimisation results for different cost functions. Plot showing the cross section of the stellarators as obtained from the optimisation using different forms of the cost function. The optimisation stops when the optimiser is unable to make further improvement (the characteristic change in the parameters during an iteration is $10^{-4}-10^{-3}$). The broken line represents the starting point, while the black, blue and red lines represent the $f_{B}$, $f_{T}$ and $f_{C}$ optimisation results respectively. The cross sections are only presented in one half of the plane but the missing plane can be reconstructed symmetrically.  Opt. | $\hat{f}_{B}$ | $\hat{f}_{C}$ | $\hat{f}_{T}$ | $\epsilon_{\mathrm{eff}}$ | Asp.
---|---|---|---|---|---
Start | 0.54 | 5.3 | 24 | 66 | 6.3
$f_{B}$ | 3.3e-3 | 2.4e-2 | 3.0e-2 | 5.1e-5 | 7.1
$f_{C}$ | 2.3e-2 | 7.6e-2 | 0.50 | 3.9e-3 | 6.8
$f_{T}$ | 6.2e-2 | 1.9e-1 | 4.8e-2 | 3.6e-2 | 6.0
Table 2: Averaged metric values for optimised stellarators. The table
summarises the average value of the metric over the volume of the optimised
stellarators with respect to the different cost functions.

What can we thus draw from this illustrative example? First, the form of the
cost function must be carefully considered. Here we only examined QS and
observed that the optimisation result varies significantly depending on the
form of the cost function. Similar behaviour will likely affect the
optimisation of other physical measures. Second, we see in this example that
the ‘simplest’ forms (in this case, $f_{B}$ and $f_{C}$) seemed to provide the
better global performance. (By simpler we mean a lower number of field
derivatives and a simpler dependence on the Fourier modes.) Increased
complexity seems to make optimisation harder. However, it is difficult to make
general statements based only on this example, and there might be exceptions
not treated in this paper.

As discussed in the previous section, close to local QS minima, improvements
in QS metrics do not always lead to improvements in physical properties. Thus
insisting on making marginal improvements on $\hat{f}_{B}$ (and for that
matter, any of the other QS measures) does not seem to be a useful exercise.
Close to a QS minima it does seem appropriate to prioritise optimisation with
respect to more physical measures. Perhaps insisting on QS at all times may
stand its ground if an absolute QS zero could be reached (a point with clear
physical meaning), such as in the optimisation for QS on a single magnetic
surface. But even then, stellarator design usually requires additional
physical considerations.

We seem then to need to alternate between an optimisation dominated by the QS
metrics and an optimisation guided by other physics metrics. A way forward is
to impose QS as an inequality constraint, so that physical metrics can
dominate optimisation when QS is ‘small enough.’ Of course there is the
unavoidable difficulty of choosing the value of a QS metric to indicate
‘sufficient QS.’ Picking a fixed value will modify the shape of optimisation
space in a non-trivial way. Instead of choosing a fixed value, a more dynamic
alternative would be to turn QS optimisation on and off guided by the relative
changes or gradient norms of the QS metrics. One might choose to not optimise
for QS at all, instead initialising the search from an already fairly
quasisymmetric configuration. This would place the optimisation in the right
basin, but lead to refinements in it based on physically-relevant features.
Possible initial seeds could be constructed from near-axis solutions.

As explored, the non-universality of the QS metrics close to minima seems to
require amendment by introducing additional physics metrics. At the same time,
detailed physical calculations are likely to lead to complicated parameter
spaces including multiple local minima compared to those of the simpler QS
metrics. Using QS metrics in conjunction with a more complex physical metric
may provide the best of both worlds. This could shed some light on recent
results in (Bader et al. (2019)). In this work, optimisation of QS through
$f_{B}$ together with a more complex physical proxy for $\alpha$ particle
losses seems to perform better than each of the pieces separately. We
hypothesise that $f_{B}$ allows for a smoother global navigation (avoiding
some potential local minima), while the physics metric helps in finding a more
physically meaningful final minimum. Using both metrics in an appropriately
balanced way seems to lead to a convenient compromise in which benefits of
both aspects are brought together.

Stellarator design often requires multi-objective optimisation. Given the
nonlinearity of the objectives of interest, adding additional terms to the
cost function can change the problem dramatically. Thus, the differences
observed above for the QS metrics will only be more complex in the presence of
other physics objectives. Those differences may be amplified or reduced
depending on their relation to the behaviour of the other metrics. Such an
analysis will be necessary to reason how to best compose QS with these other
metrics. In these complex scenarios, the idea of initially seeding
optimisation with an approximately QS configuration gains strength, as it
would make the problem simplest. Further study is required to verify these
proposed optimisation strategies.

##  5  Conclusion

In this paper we have explored three main measures of quasisymmetry: the
standard Fourier form, the two-term formulation and the triple product form.
We performed an analytical comparison of their origin and mathematical
structure. Particularly important are their differences in weighting the
asymmetric modes of $B$ in Boozer coordinates. While the standard Fourier form
treats all modes equally, the two-term form weights them differently, and the
triple vector formulation involves a mode convolution.

The differences in these valid QS metrics shows that there is no universal
measure of QS. While the measures are all rigorously equivalent when QS is
satisfied exactly (assuming $\textbf{j}\cdot\nabla\psi=0$), there is no unique
measure that quantifies deviations from QS. This lack of universality makes
the comparison between configurations particularly difficult near
approximately QS solutions, and some physically relevant property (such as
neoclassical transport through $\epsilon_{\mathrm{eff}}$ or $\alpha$ particle
confinement) would be needed for a more meaningful comparison. Near those
local minima that are not exactly QS, optimising QS cost functions for
marginal improvements generally lacks physical meaning. Away from the minima,
the global structures of the QS metrics are qualitatively similar, providing
more physical meaning to the metrics.

We also show that the mathematical form of the measure affects the
optimisation significantly via an example. In the particular example
presented, the complexity of a cost function seems to lead to ending in
alternative local minima.

In light of the above, we suggest either enforcing approximate QS through an
inequality constraint (so that optimisation close to minima is dominated by
some direct physics measure) or through an appropriate initial configuration.
Additional attention should also be paid to the form of physics or engineering
cost functions.

## Acknowledgements

The authors would like to acknowledge Daniel Dudt for fruitful discussion.
This research is primarily supported by a grant from the Simons
Foundation/SFARI (560651, AB) and DoE Contract No DE-AC02-09CH11466.

## Data availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Magnetic differential relation between $f_{C}$ and $f_{T}$

Although the triple vector product form of QS seems distinct from the other
forms given its non-linear character, we can establish a simple relation
between $f_{T}$ and $f_{C}$. To obtain this relation it is most convenient to
use the same tools we used in the derivations of the different forms of QS in
Sec. 2. Away from QS, one may still construct the symmetry field as
$\tilde{\mathbf{u}}=\nabla\Phi\times\nabla B/(\mathbf{B}\cdot\nabla B)$, where
we assume $\mathbf{B}\cdot\nabla B\neq 0$. At extremal points of the magnetic
field magnitude along $\mathbf{B}$, this construction fails. We have seen the
role of the extrema, for example, in the discussion of the single-particle
dynamics (Sec. 3). Requiring pseudo-symmetry takes care of these, as it
guarantees $\mathbf{B}\cdot\tilde{\mathbf{u}}$ to be well-behaved.
Constructing this vector field $\tilde{\mathbf{u}}$ in this way is equivalent
to enforcing (1) and (2). With these as well as $\nabla\cdot\mathbf{B}=0$,

$\mathbf{B}\cdot\nabla(\mathbf{B}\cdot\tilde{\mathbf{u}})=B^{2}\nabla\cdot\tilde{\mathbf{u}},$

using $\mathbf{j}\cdot\nabla\psi=0$ and appropriate vector identities to
rewrite the LHS. Now, from
$\mathbf{B}\cdot\tilde{\mathbf{u}}/\Phi^{\prime}=C+f_{C}/(\mathbf{B}\cdot\nabla
B)$, it then follows that $f_{C}$ and $f_{T}$ are related to each other
through a magnetic differential equation (MDE),

$\mathbf{B}\cdot\nabla\left(\frac{f_{C}}{\mathbf{B}\cdot\nabla
B}\right)=-\frac{f_{T}}{(\mathbf{B}\cdot\nabla B)^{2}}.$ (20)

This MDE shows that $f_{T}$ does, in fact, contain a higher order derivative
in comparison to $f_{C}$.

## References

  * Akçelik et al. (2006) Akçelik, V., Biros, G., Ghattas, O., Hill, J., Keyes, D. & van Bloemen Waanders, B. 2006 Parallel algorithms for pde-constrained optimization. In Parallel processing for scientific computing, pp. 291–322. SIAM.
  * Anderson et al. (1995) Anderson, F. S. B., Almagri, A. F., Anderson, D. T., Matthews, P. G., Talmadge, J. N. & Shohet, J. L. 1995 The Helically Symmetric Experiment (HSX) goals, design and status. Fusion Technology 27 (3T), 273–277.
  * Bader et al. (2019) Bader, A., Drevlak, M., Anderson, D. T., Faber, B. J., Hegna, C. C., Likin, K. M., Schmitt, J. C. & Talmadge, J. N. 2019 Stellarator equilibria with reactor relevant energetic particle losses. Journal of Plasma Physics 85 (5), 905850508\.
  * Bernardin et al. (1986) Bernardin, M. P., Moses, R. W. & Tataronis, J. A. 1986 Isodynamical (omnigenous) equilibrium in symmetrically confined plasma configurations. The Physics of Fluids 29 (8), 2605–2611.
  * Bhattacharjee et al. (1995) Bhattacharjee, A., Hayashi, T., Hegna, C. C., Nakajima, N. & Sato, T. 1995 Theory of pressure‐induced islands and self‐healing in three‐dimensional toroidal magnetohydrodynamic equilibria. Physics of Plasmas 2 (3), 883–888.
  * Boozer (1981) Boozer, A. H. 1981 Plasma equilibrium with rational magnetic surfaces. The Physics of Fluids 24 (11), 1999–2003.
  * Boozer (1983) Boozer, A. H. 1983 Transport and isomorphic equilibria. The Physics of Fluids 26 (2), 496–499.
  * Burby et al. (2020) Burby, J. W., Kallinikos, N. & MacKay, R. S. 2020 Some mathematics for quasi-symmetry. Journal of Mathematical Physics 61 (9), 093503.
  * Calvo et al. (2014) Calvo, I., Parra, F. I., Velasco, J. L. & Alonso, J. A. 2014 Flow damping in stellarators close to quasisymmetry. Plasma Physics and Controlled Fusion 57 (1), 014014.
  * Constantin et al. (2021) Constantin, P., Drivas, T. D. & Ginsberg, D. 2021 On quasisymmetric plasma equilibria sustained by small force. Journal of Plasma Physics 87 (1), 905870111.
  * Dekeyser et al. (2014) Dekeyser, W., Reiter, D. & Baelmans, M. 2014 Automated divertor target design by adjoint shape sensitivity analysis and a one-shot method. Journal of Computational Physics 278, 117–132.
  * Drevlak et al. (2013) Drevlak, M., Brochard, F., Helander, P., Kisslinger, J., Mikhailov, M., Nührenberg, C., Nührenberg, J. & Turkin, Y. 2013 Estell: A quasi-toroidally symmetric stellarator. Contributions to Plasma Physics 53 (6), 459–468.
  * Garabedian (2008) Garabedian, P. R. 2008 Three-dimensional analysis of tokamaks and stellarators. Proceedings of the National Academy of Sciences 105 (37), 13716–13719.
  * Garabedian & McFadden (2009) Garabedian, P. R. & McFadden, G. B. 2009 Design of the demo fusion reactor following iter. Journal of research of the National Institute of Standards and Technology 114 (4), 229.
  * Garren & Boozer (1991) Garren, D. A. & Boozer, A. H. 1991 Existence of quasihelically symmetric stellarators. Physics of Fluids B: Plasma Physics 3 (10), 2822–2834.
  * Grafakos (2008) Grafakos, L. 2008 Fourier Analysis on the Torus, pp. 1–88. New York, NY: Springer New York.
  * Helander (2014) Helander, P. 2014 Theory of plasma confinement in non-axisymmetric magnetic fields. Reports on Progress in Physics 77 (8), 087001.
  * Helander & Simakov (2008) Helander, P. & Simakov, A. N. 2008 Intrinsic ambipolarity and rotation in stellarators. Phys. Rev. Lett. 101, 145003.
  * Henneberg et al. (2019) Henneberg, S., Drevlak, M., Nuhrenberg, C., Beidler, C., Turkin, Y., Loizu, J. & Helander, P. 2019 Properties of a new quasi-axisymmetric configuration. Nuclear Fusion 59 (2), 026014.
  * Henneberg et al. (2021) Henneberg, S., Helander, P. & Drevlak, M. 2021 Representing the boundary of stellarator plasmas. Journal of Plasma Physics 87 (5), 905870503.
  * Ho & Kulsrud (1987) Ho, D. D. & Kulsrud, R. M. 1987 Neoclassical transport in stellarators. The Physics of Fluids 30 (2), 442–461.
  * Ku & Boozer (2010) Ku, L. & Boozer, A. 2010 New classes of quasi-helically symmetric stellarators. Nuclear Fusion 51 (1), 013004\.
  * Landreman & Catto (2012) Landreman, M. & Catto, P. J. 2012 Omnigenity as generalized quasisymmetry. Physics of Plasmas 19 (5), 056103\.
  * Landreman & Sengupta (2018) Landreman, M. & Sengupta, W. 2018 Direct construction of optimized stellarator shapes. part 1. theory in cylindrical coordinates. Journal of Plasma Physics 84 (6), 905840616\.
  * Landreman & Sengupta (2019) Landreman, M. & Sengupta, W. 2019 Constructing stellarators with quasisymmetry to high order. Journal of Plasma Physics 85 (6), 815850601.
  * Lazerson et al. (2020) Lazerson, S., Schmitt, J., Zhu, C., Breslau, J., Developers, S. & others 2020 STELLOPT. Tech. Rep.. Princeton Plasma Physics Lab.(PPPL), Princeton, NJ (United States).
  * Littlejohn (1983) Littlejohn, R. G. 1983 Variational principles of guiding centre motion. Journal of Plasma Physics 29 (1), 111–125.
  * Martin & Landreman (2020) Martin, M. F. & Landreman, M. 2020 Impurity temperature screening in stellarators close to quasisymmetry. Journal of Plasma Physics 86 (3).
  * Mikhailov et al. (2002) Mikhailov, M., Shafranov, V., Subbotin, A., Isaev, M., hrenberg, J. N., Zille, R. & Cooper, W. 2002 42 (11), L23–L26.
  * Najmabadi et al. (2008) Najmabadi, F., Raffray, A., Abdel-Khalik, S., Bromberg, L., Crosatti, L., El-Guebaly, L., Garabedian, P., Grossman, A., Henderson, D., Ibrahim, A. & others 2008 The ARIES-CS compact stellarator fusion power plant. Fusion Science and Technology 54 (3), 655–672.
  * Nemov et al. (2014) Nemov, V. V., Kasilov, S. V. & Kernbichler, W. 2014 Collisionless high energy particle losses in optimized stellarators calculated in real-space coordinates. Physics of Plasmas 21 (6), 062501.
  * Nemov et al. (1999) Nemov, V. V., Kasilov, S. V., Kernbichler, W. & Heyn, M. F. 1999 Evaluation of $1/\nu$ neoclassical transport in stellarators. Physics of Plasmas 6 (12), 4622–4632.
  * Nührenberg & Zille (1988) Nührenberg, J. & Zille, R. 1988 Quasi-helically symmetric toroidal stellarators. Physics Letters A 129 (2), 113 – 117.
  * Paul et al. (2020) Paul, E. J., Antonsen, T., Landreman, M. & Cooper, W. A. 2020 Adjoint approach to calculating shape gradients for three-dimensional magnetic confinement equilibria. part 2. applications. Journal of Plasma Physics 86 (1), 905860103.
  * Plunk & Helander (2018) Plunk, G. & Helander, P. 2018 Quasi-axisymmetric magnetic fields: weakly non-axisymmetric case in a vacuum. Journal of Plasma Physics 84 (2), 905840205.
  * Reiman & Boozer (1984) Reiman, A. & Boozer, A. H. 1984 Island formation and destruction of flux surfaces in three‐dimensional mhd equilibria. The Physics of Fluids 27 (10), 2446–2454.
  * Rodríguez et al. (2020) Rodríguez, E., Helander, P. & Bhattacharjee, A. 2020 Necessary and sufficient conditions for quasisymmetry. Physics of Plasmas 27 (6), 062501.
  * Rodríguez & Bhattacharjee (2021a) Rodríguez, E. & Bhattacharjee, A. 2021a Islands and current singularities in quasisymmetric toroidal plasmas. Physics of Plasmas 28 (9), 092506.
  * Rodríguez & Bhattacharjee (2021b) Rodríguez, E. & Bhattacharjee, A. 2021b Solving the problem of overdetermination of quasisymmetric equilibrium solutions by near-axis expansions. i. generalized force balance. Physics of Plasmas 28 (1), 012508.
  * Rodríguez & Bhattacharjee (2021c) Rodríguez, E. & Bhattacharjee, A. 2021c Solving the problem of overdetermination of quasisymmetric equilibrium solutions by near-axis expansions. ii. circular axis stellarator solutions. Physics of Plasmas 28 (1), 012509.
  * Rodríguez et al. (2021) Rodríguez, E., Sengupta, W. & Bhattacharjee, A. 2021 Generalized boozer coordinates: A natural coordinate system for quasisymmetry. Physics of Plasmas 28 (9), 092510.
  * Sanchez et al. (2000) Sanchez, R., Hirshman, S. P., Ware, A. S., Berry, L. A. & Spong, D. A. 2000 Ballooning stability optimization of low-aspect-ratio stellarators. Plasma Physics and Controlled Fusion 42 (6), 641–653.
  * Skovoroda (2005) Skovoroda, A. A. 2005 3d toroidal geometry of currentless magnetic configurations with improved confinement. Plasma Physics and Controlled Fusion 47 (11), 1911–1924.
  * Spitzer (1958) Spitzer, L. 1958 The stellarator concept. The Physics of Fluids 1 (4), 253–264.
  * Tessarotto et al. (1996) Tessarotto, M., Johnson, J. L., White, R. B. & Zheng, L. 1996 Quasi‐helical magnetohydrodynamic equilibria in the presence of flow. Physics of Plasmas 3 (7), 2653–2663.
  * Tessarotto et al. (1995) Tessarotto, M., Johnson, J. L. & Zheng, L. J. 1995 Hamiltonian approach to the magnetostatic equilibrium problem. Physics of Plasmas 2 (12), 4499–4512.
  * Zarnstorff et al. (2001) Zarnstorff, M., Berry, L., Brooks, A., Fredrickson, E., Fu, G., Hirshman, S., Hudson, S., Ku, L., Lazarus, E., Mikkelsen, D. & others 2001 Physics of the compact advanced stellarator ncsx. Plasma Physics and Controlled Fusion 43 (12A), A237.