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 $|\boldsymbol {B}|$ are weighted in each of these metrics. Then a comparison with $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:

(3.1)\begin{equation} B=\sum_{n,m}B_{nm}{\rm e}^{{\rm i}(m\theta-n\phi)}, \end{equation}

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

(3.2)\begin{align} f_C& =\frac{B^{2}}{\tilde{\alpha}-\iota}(\partial_\phi+\tilde{\alpha}\partial_\theta)B\nonumber\\ & ={\rm i}\frac{B^{2}}{\iota-\tilde{\alpha}}\sum_{n/m\neq\tilde{\alpha}}(n-m\tilde{\alpha})B_{nm} {\rm e}^{{\rm i}(m\theta-n\phi)}. \end{align}

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 $\boldsymbol {Q}=\mathcal {J}\boldsymbol {\nabla }\psi \times \boldsymbol {\nabla }\chi$ must vanish, that is, $\oint (f_C/B^{2})\,\mathrm {d}l/|\boldsymbol {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 with $f_B$ we need to construct a single flux-surface scalar from $f_C$. We take (3.2) and write it as

(3.3)\begin{equation} B_{nm}={-}{\rm i}\frac{\iota-\tilde{\alpha}}{n-m\tilde{\alpha}}\left(\frac{f_C}{B^{2}}\right)_{nm}, \end{equation}

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

(3.4)\begin{equation} \sum_{n,m}(n-m\tilde{\alpha})^{2}|B_{nm}|^{2}=\frac{(\iota-\tilde{\alpha})^{2}}{(2{\rm \pi})^{2}}\int\,\mathrm{d}\theta\,\mathrm{d}\phi\left(\frac{f_C}{B^{2}}\right)^{2}, \end{equation}

which can be rewritten as

(3.5)\begin{equation} \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}, \end{equation}

where $\overline {(\cdots )}$ denotes the $(0,0)$ Fourier component of the expression in parentheses. Equation (3.5) 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 (2.7) and (3.5) 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 nonlinear character of the triple product formulation is clear. We explicitly write

(3.6)\begin{equation} (\mathcal{J}^{2}f_T)_{nm}={\rm i}\sum_{k,l}(\iota l-k)(km-nl)B_{n-k,m-l}B_{k,l}, \end{equation}

where $\mathcal {J}$ is the Jacobian of Boozer coordinates and $k,l$ are integers. The nonlinearity 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 of 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

(3.7)\begin{equation} (\mathcal{J}^{2}f_T)_{nm}={\rm i}(\tilde{\alpha}m-n)(\iota m-n)\sum_llB_{\tilde{\alpha}l,l}\tilde{B}_{n-\tilde{\alpha}l,m-l}. \end{equation}

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.

3.2.1 Single-particle dynamics

The concept of QS as presented in this work is built on single-particle dynamics. In § 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 single-particle dynamics.

Let us start from the definition of the momentum: $\bar {p}=-\varPhi +v_\parallel \boldsymbol {b}\boldsymbol {\cdot }\boldsymbol {u}$, where $v_\parallel$ is the parallel velocity of the particle under consideration (Rodríguez et al. Reference Rodríguez, Helander and Bhattacharjee2020). For a QS configuration where $\boldsymbol {u}$ satisfies (2.1)–(2.3), this momentum is approximately conserved to $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 $\boldsymbol {u}$. We can always choose this field in a way that (2.1) and (2.2) are satisfied by construction. To do so we assume the existence of flux surfaces as well as pseudosymmetry (Mikhailov et al. Reference Mikhailov, Shafranov, Subbotin, Isaev, Hrenberg, Zille and Cooper2002; Skovoroda Reference Skovoroda2005). The latter guarantees that the contours of constant $|\boldsymbol {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 $\boldsymbol {b}\boldsymbol {\cdot }\boldsymbol {u}\rightarrow \infty$. Once $\boldsymbol {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 (Littlejohn Reference Littlejohn1983; Rodríguez et al. Reference Rodríguez, Helander and Bhattacharjee2020). We may thus write, to $O(\rho /L)$,

(3.8)\begin{equation} \frac{\mathrm{d}}{\mathrm{d}t}\bar{p} =\frac{v_\parallel^{2}}{B^{2}}\left[\boldsymbol{B} \boldsymbol{\cdot}\boldsymbol{\nabla} (\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{u}) -\boldsymbol{j}\boldsymbol{\cdot}\boldsymbol{\nabla}\psi\right] ={-}\left(\frac{v_\parallel}{\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} B}\right)^{2}f_T. \end{equation}

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 Reference Landreman and Catto2012; Helander Reference Helander2014). It is also important to note from (3.8) that the cost function $f_T$ appears with a factor of $1/(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } B)^{2}$. This factor penalises deviations from QS close to the extrema of $|\boldsymbol {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/(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } B)^{2}$ numerically ill-behaved as a modified cost function. A practical implementation could perhaps be achieved by regularising the combination, as has been employed in optimisation for pseudosymmetry (Mikhailov et al. Reference Mikhailov, Shafranov, Subbotin, Isaev, Hrenberg, Zille and Cooper2002). 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 $|\boldsymbol {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 Reference Helander2014). 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 figure 1). Take the turning points that define the integral for $\mathcal {J}_\parallel$ to be at a field magnitude $|\boldsymbol {B}|=B_t$. We assume that these bounce points exist on both of these field lines (as we are looking at the $\mathcal {J}_\parallel$ associated with a class of particles defined by its bounce points). Then taking the line integral along the contour in figure 1, we compute

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

To obtain the equality, we apply the 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 $\boldsymbol {j}\boldsymbol {\cdot }\boldsymbol {\nabla }\psi =0$ and noting that the integrals along $C_{B+}$ and $C_{B-}$ vanish as $v_\parallel =0$ wherever $|\boldsymbol {B}|=B_t$, we may rewrite

(3.10)\begin{equation} \partial_\alpha \mathcal{J}_\parallel{=}{-}\int\frac{\boldsymbol{\nabla}\psi\times\boldsymbol{\nabla} B\boldsymbol{\cdot}\boldsymbol{B}}{\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} B}\frac{1}{B^{2}}\left(\frac{E}{v_\parallel}+\frac{v_\parallel}{2}\right)\,\mathrm{d}B, \end{equation}

where $B$ is being used to parameterise the integral along the field line. To be precise, we should split the integral along the field line into concatenated pieces joint at $|\boldsymbol {B}|$ extrema. The vanishing of $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\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 drop this distinction for simplicity.

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 $|\boldsymbol {B}|=B_t$.

One may write the kinematic expression $(E/v_\parallel +v_\parallel /2)/B^{2}=\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }(v_\parallel /B)/(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\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/(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\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 (A2):

(3.11)\begin{equation} \partial_\alpha \mathcal{J}_\parallel{=}\left.\frac{\boldsymbol{\nabla}\psi\times\boldsymbol{\nabla} B\boldsymbol{\cdot}\boldsymbol{B}}{\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} B}\frac{v_\parallel}{B}\right|_{v_\parallel{=}0}+\int\frac{f_T}{(\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} B)^{2}}v_\parallel\,\mathrm{d}l. \end{equation}

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 $(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } B)^{2}$. However, this $f_T$ term is not the only term in (3.11). For the expression in (3.11) to be finite (or vanish as we want for QS), the boundary term in (3.11) should not diverge, which requires $\boldsymbol {\nabla }\psi \times \boldsymbol {\nabla } B\boldsymbol {\cdot }\boldsymbol {B}=0$ wherever $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } B=0$. This condition is precisely the condition of pseudosymmetry introduced earlier (Mikhailov et al. Reference Mikhailov, Shafranov, Subbotin, Isaev, Hrenberg, Zille and Cooper2002; Skovoroda Reference Skovoroda2005). Note that this pseudosymmetry condition appears separate from the $f_T$ integral in (3.11). 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 § 2, requiring $f_T=0$ everywhere includes this condition. Guaranteeing that 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 natural measures of QS in the single-particle perspective. However, when the configuration fails to be pseudosymmetric, 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 § 2, the function $C$ in $f_C$ is a measure of the banana width of bouncing particles in the QS limit. The expression in (3.2) 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 Reference Helander2014).

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. Reference Bader, Drevlak, Anderson, Faber, Hegna, Likin, Schmitt and Talmadge2019) or expensive simulations are needed.

3.2.3 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, Johnson & Zheng Reference Tessarotto, Johnson and Zheng1995; Rodríguez & Bhattacharjee Reference Rodríguez and Bhattacharjee2021a). Fourier modes of $|\boldsymbol {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 Reference Reiman and Boozer1984; Rodríguez & Bhattacharjee Reference Rodríguez and Bhattacharjee2021a). The PS current singularity arises due to the periodicity constraint on the magnetic differential equation,

(3.12)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\left(\frac{j_{\|}}{B}+\frac{p'C}{B^{2}}\right)={-}2p' \frac{f_C}{B^{3}}, \end{equation}

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

(3.13)\begin{equation} p' \oint \frac{f_C}{B^{4}} \, \mathrm{d}l = 0, \end{equation}

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 § 2). Note that solvability may still be obtained away from QS if $p' = 0$ or the line integral vanishes.

One may generally show that in magnetohydrostatics, 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 Reference Reiman and Boozer1984; Bhattacharjee et al. Reference Bhattacharjee, Hayashi, Hegna, Nakajima and Sato1995; Rodríguez & Bhattacharjee Reference Rodríguez and Bhattacharjee2021a). 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 Reference Rodríguez and Bhattacharjee2021a). This is of course only a partial account of the relevant island physics, focusing on the potentially deleterious seeding of islands. The appearance of islands will in practice depend on a plethora of additional effects beyond the simple magnetohydrostatic equilibria considerations, including, for instance, the healing effects of flows or additional kinetic effects (Fitzpatrick Reference Fitzpatrick1993; Hegna Reference Hegna2011, Reference Hegna2012).

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 nonlinearity 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

(4.1)\begin{equation} \hat{f}_B=\sqrt{\frac{1}{\overline{B^{2}}}f_B}. \end{equation}

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 take the natural and more conventional form (4.1).

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 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 (3.2), 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\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } B$ for $\chi = \theta$, but for simplicity we do not do so here. In view of the normalisation of $f_B$, we normalise and scalarise $f_C$ as follows:

(4.2)\begin{equation} \hat{f}_C=\sqrt{\frac{\langle f_C^{2}\rangle}{\langle B^{2}\rangle^{3}}|\tilde{\alpha}-\iota|^{2}}. \end{equation}

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

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

(4.3)\begin{equation} \hat{f}_T=\sqrt{\langle f_T^{2}\rangle\frac{R^{4}}{\langle{B^{2}}\rangle^{4}}}. \end{equation}

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/(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\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 focus on the comparison of the metrics in configurations that have been designed to be approximately quasisymmetric. Then we study how these measures differ when the system has not necessarily been optimised to be quasisymmetric.

We present in figure 2 a comparison of the metrics in six quasisymmetric configurations:Footnote ¹ GAR (Garabedian Reference Garabedian2008; Garabedian & McFadden Reference Garabedian and McFadden2009), HSX (Anderson et al. Reference Anderson, Almagri, Anderson, Matthews, Talmadge and Shohet1995), NCSX (Zarnstorff et al. Reference Zarnstorff, Berry, Brooks, Fredrickson, Fu, Hirshman, Hudson, Ku, Lazarus and Mikkelsen2001), ESTELL (Drevlak et al. Reference Drevlak, Brochard, Helander, Kisslinger, Mikhailov, Nührenberg, Nührenberg and Turkin2013), QHS48 (Ku & Boozer Reference Ku and Boozer2010) and ARIESCS (Najmabadi et al. Reference Najmabadi, Raffray, Abdel-Khalik, Bromberg, Crosatti, El-Guebaly, Garabedian, Grossman, Henderson and Ibrahim2008). Within a given device, the overall correlation between cost functions is good. All measures qualitatively show a similar behaviour: QS improves as we move from the boundary inwards. However, in regions closest to the axis, this correspondence is not as strong (especially in configurations such as GAR and QHS48). However, the same cost function correspondence does not seem to hold if we make comparisons between different configurations.

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 joined by lines in order from the centre to the edge of the configurations (the rightmost points represent the edge).

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 figure 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 figure 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. Reference Nemov, Kasilov, Kernbichler and Heyn1999) with the three cost functions for the devices of figure 2. Other physics measures such as alpha particle confinement could also be suited for such comparison (Bader et al. Reference Bader, Drevlak, Anderson, Faber, Hegna, Likin, Schmitt and Talmadge2019). There exists some correlation between the transport and the normalised values of $f_C$ and $f_T$. Especially remarkable is 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 Reference Martin and Landreman2020).

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 suggest 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 nonlinear $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:

(4.4)\begin{equation} \left. \begin{aligned} R & = \sum_{nm} R_{nm} \cos(m \vartheta - n N_p \phi_c), \\ Z & = \sum_{nm} Z_{nm} \sin(m \vartheta - n N_p \phi_c), \end{aligned} \right\} \end{equation}

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 quasihelicallysymmetric configurationFootnote ² with $N_p=4$ field periods. We 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. Reference Lazerson, Schmitt, Zhu and Breslau2020) with implementations of the three different forms of the cost functions. 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 (4.4)), and Figure 5 the magnitude of the associated gradient. The metrics are evaluated at six equally spaced magnetic flux surfaces, and their contributions are added together.

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 p. 28. The minima for each of these are found (clockwise) at: $(-0.031,0.031)$, $(-0.031,0.031)$, $(-0.037,0.045)$ and $(-0.037,0.037)$.

Figure 5. Gradients in 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 parameter space gradients of $\hat {f}_B$, $\hat {f}_C$, $\epsilon _\mathrm {eff}$ and $\hat {f}_T$ as a function of the $(2,1)$ surface modes. Scales of the gradients are significantly different between measures.

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 figure 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 nonlinear 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. Second, 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 § 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 nonlinearity 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 figure 6 and table 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 p. 28. The optimisation uses a BFGS gradient-based method (as part of the suite MANGOFootnote ³) 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 (4.4). 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 into the differences between the various QS metrics. The optimisation is performed for several values of the finite-difference step size to ensure convergence.

Figure 6. 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.

Table 2. Averaged metric values for optimised stellarators. Summary of the average value of the metric over the volume of the optimised stellarators with respect to the different cost functions.

Comparing cross-sections in figure 6 and the values in table 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 § 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 $|\boldsymbol {B}|$. The relation is generally complex and nonlinear, but thinking in these terms may still shed some light on 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 lift 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$ terminated at a different local minimum. The presence of nonlinearities in $f_T$ is 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 figure 4), which may make difficult reaching 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. An important additional point of distinction of the $\hat {f}_T$ optimisation with respect to the others is, as emphasised in § 2, that the form of symmetry of the configuration is not specified. Hence, the optimisation must locally decide upon the type of QS towards which to evolve. If close to a configuration with good QS, optimisation will naturally bring the system towards the corresponding helicity. However, generally it is difficult to predict given an initial configuration the form of QS towards which optimisation will tend.

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 Reference Landreman and Sengupta2019). 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 in changing the aspect ratio without sacrificing QS along the way and the choice of parameterisation (Henneberg, Helander & Drevlak Reference Henneberg, Helander and Drevlak2021) (4.4) might provide an explanation for this behaviour.

What conclusions 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 a 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 minimum it does seem appropriate to prioritise optimisation with respect to more physical measures. Optimising stellarators based solely on physical properties has been shown to be successful in the literature (Hirshman et al. Reference Hirshman, Spong, Whitson, Lynch, Batchelor, Carreras and Rome1998, Reference Hirshman, Spong, Whitson, Nelson, Batchelor, Lyon, Sanchez, Brooks, Y.-Fu and Goldston1999; Spong et al. Reference Spong, Hirshman, Whitson, Batchelor, Sanchez, Carreras, Lynch, Lyon, Valanju and Miner2000). Here we propose a hybrid approach between optimising for QS alone and optimising for physical metrics (without QS). The QS optimisation should bring the search close to a region of configurations with QS-like properties following a simpler form of cost function. Then for refinement, it will follow the more detailed physical cost function optimisation. 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 with 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. (Reference Bader, Drevlak, Anderson, Faber, Hegna, Likin, Schmitt and Talmadge2019). 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.