## 1 Introduction

Quasisymmetry is a type of continuous symmetry in the strength of a toroidal magnetic field
$B=|\boldsymbol{B}|$
that does not require continuous symmetry of the magnetic field vector
$\boldsymbol{B}$
(Boozer Reference Boozer1983; Nührenberg & Zille Reference Nührenberg and Zille1988; Boozer Reference Boozer1995; Helander Reference Helander2014). As a consequence of the conservation laws associated with quasisymmetry or full axisymmetry of
$\boldsymbol{B}$
, both symmetries enable confinement of charged particles and plasma. However, confinement with axisymmetric
$\boldsymbol{B}$
requires a large electric current in the confinement region that is prone to instabilities and hard to sustain, while quasisymmetric confinement does not. Hence, non-axisymmetric toroidal magnetic fields (‘stellarators’) with quasisymmetry offer the promise of stable and efficient confinement of high-temperature plasma for fusion energy. Quasisymmetric stellarators would also enable magnetic confinement of plasmas with density that is too low to support a substantial electric current, such as electron–positron plasmas for basic physics studies (Pedersen *et al.*
Reference Pedersen, Danielson, Hugenschmidt, Marx, Sarasola, Schauer, Schweikhard, Surko and Winkler2012).

Several quasisymmetric magnetic field configurations have been found numerically, mostly by using optimization over the space of boundary magnetic surface shapes to minimize symmetry-breaking Fourier modes of
$B$
(Nührenberg & Zille Reference Nührenberg and Zille1988; Anderson *et al.*
Reference Anderson, Almagri, Anderson, Matthews, Talmadge and Shohet1995; Garabedian Reference Garabedian1996; Zarnstorff *et al.*
Reference Zarnstorff, Berry, Brooks, Fredrickson, Fu, Hirshman, Hudson, Ku, Lazarus and Mikkelsen2001; Ku & Boozer Reference Ku and Boozer2011; Drevlak *et al.*
Reference Drevlak, Brochard, Helander, Kisslinger, Mikhailov, Nührenberg, Nührenberg and Turkin2013; Henneberg *et al.*
Reference Henneberg, Drevlak, Nührenberg, Beidler, Turkin, Loizu and Helander2019). While this optimization approach is effective, it does not provide much insight into the size and character of the solution space, and it requires good initial guesses for the numerical iteration. Hence, one can never be sure that all the interesting regions of parameter space have been found, for perhaps a different initial guess would yield a new solution. The optimization approach requires significant computation time, and so it is expensive to generate parameterized families of solutions.

An alternative to the optimization approach is to construct quasisymmetric configurations directly using an analytic expansion, in the smallness of either the departure from axisymmetry of
$\boldsymbol{B}$
(Plunk & Helander Reference Plunk and Helander2018) or in the distance from the magnetic axis (equivalent to an expansion in large aspect ratio). Near-axis expansions have been explored by several authors (Mercier Reference Mercier1964; Solov’ev & Shafranov Reference Solov’ev and Shafranov1970; Lortz & Nührenberg Reference Lortz and Nührenberg1976), with the particular case of quasisymmetry examined by Garren & Boozer (Reference Garren and Boozer1991*a*
,Reference Garren and Boozer
*b*
). The near-axis expansion, although it is an approximation, is always accurate in the core of any stellarator, even stellarators for which the aspect ratio of the outermost surface is not large. In a recent series of papers (Landreman & Sengupta Reference Landreman and Sengupta2018; Landreman, Sengupta & Plunk Reference Landreman, Sengupta and Plunk2019; Plunk, Landreman & Helander Reference Plunk, Landreman and Helander2019), the near-axis expansion was developed into practical procedures for constructing fields with quasisymmetry, or the more general condition of omnigenity. It was also shown that, close to the axis, quasisymmetric configurations obtained by conventional optimization closely match configurations generated by the construction (Landreman Reference Landreman2019). The configurations presented to date from this near-axis construction have been quasisymmetric to first order in
$r/{\mathcal{R}}$
, where
$r$
is the typical distance from the axis and
${\mathcal{R}}$
denotes a scale length of the magnetic axis (such as a typical radius of curvature). At this order, due to regularity conditions at the magnetic axis, the magnetic surfaces must have an elliptical cross-section. It was found that the space of quasisymmetric configurations to this order can be parameterized by the shape of the axis together with three other numbers.

In the present paper, we extend the construction to next order in
$r/{\mathcal{R}}$
. The equations describing quasisymmetry to
$O((r/{\mathcal{R}})^{2})$
were derived in the appendix of Garren & Boozer (Reference Garren and Boozer1991*a*
), but no solutions were presented before now. At this order, several important effects appear for the first time, including triangularity and Shafranov shift. By extending the construction to
$O((r/{\mathcal{R}})^{2})$
, more complicated and realistic stellarator shapes will be generated, and quasisymmetry will be achieved to higher accuracy. The extension of the model to
$O((r/{\mathcal{R}})^{2})$
only slightly increases the computational cost of solving the equations, which remains on the level of milliseconds on one CPU. This time is far faster than a traditional three-dimensional (3-D) equilibrium calculation, which typically requires of the order of 10 seconds or more.

In Garren and Boozer’s original work, it was argued quasisymmetry can be achieved (in the absence of axisymmetry) to $O((r/{\mathcal{R}})^{2})$ but not to $O((r/{\mathcal{R}})^{3})$ , so departures from quasisymmetry should scale as the cube of the inverse aspect ratio. However, despite various numerical calculations of quasisymmetric configurations using optimization since 1991, there does not appear to have been a numerical demonstration of this predicted scaling. Using the construction here we are able to numerically demonstrate this predicted ideal scaling for the first time (figure 8). An implication of this scaling is that quasisymmetry can be achieved to arbitrary precision, in the following sense. Given any desired small level of symmetry-breaking Fourier modes in $B$ , and given any desired axis shape (constrained only by the requirement that its curvature cannot vanish), there is some aspect ratio above which quasisymmetry can be achieved to the desired precision. To emphasize this point, we will present examples of non-axisymmetric configurations in which quasisymmetry is realized to unprecedented precision, with the symmetry-breaking mode amplitudes orders of magnitude smaller than in previously reported configurations.

A primary application of the work here is to generate input data for stellarator equilibrium codes such as VMEC (Hirshman & Whitson Reference Hirshman and Whitson1983) or optimization codes such as STELLOPT (Spong *et al.*
Reference Spong, Hirshman, Whitson, Batchelor, Carreras, Lynch and Rome1998; Reiman *et al.*
Reference Reiman, Fu, Hirshman, Ku, Monticello, Mynick, Redi, Spong, Zarnstorff and Blackwell1999) or ROSE (Drevlak *et al.*
Reference Drevlak, Beidler, Geiger, Helander and Turkin2019). For these applications, the input we must generate is the shape (or initial shape) of a boundary magnetic surface. In the
$O(r/{\mathcal{R}})$
construction (Landreman *et al.*
Reference Landreman, Sengupta and Plunk2019; Plunk *et al.*
Reference Plunk, Landreman and Helander2019), it was possible to plug a finite value of
$r$
into the near-axis expansion to obtain the boundary surface. In turns out that, at
$O((r/{\mathcal{R}})^{2})$
, this substitution requires some care. We will show that, in fact, a part of the
$O((r/{\mathcal{R}})^{3})$
shape must be retained in order to generate a boundary surface that is consistent with the desired field strength to
$O((r/{\mathcal{R}})^{2})$
. Once this step is taken, we will construct boundary surfaces, then use the VMEC code to generate 3-D equilibria inside the boundaries, and show that quasisymmetry-breaking modes of
$B$
in these equilibria are small and scale with the aspect ratio, as expected.

In the following section, notation will be introduced and the near-axis expansion will be outlined. Section 3 describes the analysis of generating a finite-aspect-ratio boundary surface from the near-axis expansion, and the need for including some $O((r/{\mathcal{R}})^{3})$ terms. The numerical method for solving the equations is detailed in § 4, and several examples of constructed quasi-axisymmetric and quasi-helically symmetric configurations are presented in § 5. We discuss the results and conclude in § 6. Several detailed analytic calculations can be found in the appendices. Appendix A gives the equations for $O((r/{\mathcal{R}})^{2})$ quasisymmetry, derived using a new method that reduces the algebra required. Appendix B gives a detailed proof of results presented in § 3. Finally, one method for converting the constructed boundary shapes to cylindrical coordinates is presented in appendix C.

## 2 Near-axis expansion

Our goal is to relate the three-dimensional shapes of flux surfaces to the magnetic field strength in Boozer coordinates
$(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$
. In this section, we introduce the main features of the expansion, and many of the explicit expressions are given in appendix A. While the expansion here is equivalent to the one in Garren & Boozer (Reference Garren and Boozer1991*a*
,Reference Garren and Boozer
*b*
), our approach in appendix A provides a streamlined method to derive the equations at each order. Throughout the analysis, we assume that good nested flux surfaces exist in the region of interest near the axis.

In Boozer coordinates, the magnetic field has the forms

where $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D713}$ is the toroidal flux, $\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})$ is the rotational transform, $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D711}$ are the poloidal and toroidal Boozer angles and $I$ and $G$ are constant on $\unicode[STIX]{x1D713}$ surfaces. To consider quasi-helical symmetry later in the analysis, it is convenient to introduce a helical angle $\unicode[STIX]{x1D717}=\unicode[STIX]{x1D703}-N\unicode[STIX]{x1D711}$ , where $N$ is a constant integer. Then

where $\unicode[STIX]{x1D704}_{N}=\unicode[STIX]{x1D704}-N$ .

The position vector $\boldsymbol{r}$ at a general point in a neighbourhood of the axis can be described by

where $r$ is the flux surface label defined by $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D713}=\unicode[STIX]{x03C0}r^{2}\bar{B}$ , with $\bar{B}$ a constant reference field strength, and $\boldsymbol{r}_{0}(\unicode[STIX]{x1D711})$ is the position vector along the magnetic axis. Here, the orthonormal vectors $(\boldsymbol{t},\boldsymbol{n},\boldsymbol{b})$ give the Frenet–Serret frame of the magnetic axis. These vectors satisfy

*a*-

*d*) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D711}}{\text{d}\ell }\frac{\text{d}\boldsymbol{r}_{0}}{\text{d}\unicode[STIX]{x1D711}}=\boldsymbol{t},\quad \frac{\text{d}\unicode[STIX]{x1D711}}{\text{d}\ell }\frac{\text{d}\boldsymbol{t}}{\text{d}\unicode[STIX]{x1D711}}=\unicode[STIX]{x1D705}\boldsymbol{n},\quad \frac{\text{d}\unicode[STIX]{x1D711}}{\text{d}\ell }\frac{\text{d}\boldsymbol{n}}{\text{d}\unicode[STIX]{x1D711}}=-\unicode[STIX]{x1D705}\boldsymbol{t}+\unicode[STIX]{x1D70F}\boldsymbol{b},\quad \frac{\text{d}\unicode[STIX]{x1D711}}{\text{d}\ell }\frac{\text{d}\boldsymbol{b}}{\text{d}\unicode[STIX]{x1D711}}=-\unicode[STIX]{x1D70F}\boldsymbol{n},\end{eqnarray}$$

and $\boldsymbol{t}\times \boldsymbol{n}=\boldsymbol{b}$ , where $\ell$ is the arclength along the axis, $\unicode[STIX]{x1D705}(\unicode[STIX]{x1D711})$ is the axis curvature and $\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D711})$ is the axis torsion. (Garren and Boozer use the opposite sign convention for torsion.) Using the dual relations,

where $\sqrt{g}=(\unicode[STIX]{x2202}\boldsymbol{r}/\unicode[STIX]{x2202}r)\cdot (\unicode[STIX]{x2202}\boldsymbol{r}/\unicode[STIX]{x2202}\unicode[STIX]{x1D717})\times (\unicode[STIX]{x2202}\boldsymbol{r}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711})$ is the Jacobian, then (2.2)–(2.3) can be expressed in terms of the $(\boldsymbol{t},\boldsymbol{n},\boldsymbol{b})$ vectors and derivatives of $(X,Y,Z)$ . Equating (2.2)–(2.3) then gives three scalar equations, (A 2)–(A 4). The field strength can be expressed in terms of derivatives of $(X,Y,Z)$ using the square of either (2.2) or (2.3). The former turns out to be more useful, and is given in (A 19).

These equations are supplemented by the equilibrium condition [𝜵 × (2.3) ]× (2.2) =𝜇 _{0}𝜵*p*, where
$p(r)$
is the pressure. The average of this condition over
$\unicode[STIX]{x1D717}$
and
$\unicode[STIX]{x1D711}$
gives

while the $\unicode[STIX]{x1D717}$ and $\unicode[STIX]{x1D711}$ dependence of the equilibrium condition implies

The near-axis expansion is then introduced by writing

with analogous expressions for $Y$ and $Z$ . Other than $r$ , all scale lengths in the system are ordered as ${\mathcal{R}}$ , so (2.9) represents an expansion in $r/{\mathcal{R}}$ . The field strength is expanded similarly but with an $O((r/{\mathcal{R}})^{0})$ term,

and $\unicode[STIX]{x1D6FD}(r,\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})$ is expanded in the same way. The profile functions $G(r)$ , $I(r)$ , $p(r)$ and $\unicode[STIX]{x1D704}_{N}(r)$ are analytic functions of $\unicode[STIX]{x1D713}$ , so their expansions contain only even powers of $r$ ,

Since $I(r)$ is proportional to the toroidal current inside the surface $r$ , then $I_{0}=0$ . From analyticity considerations near the axis (see appendix A of Landreman & Sengupta (Reference Landreman and Sengupta2018)), the expansion coefficients have the form

The expansions of $Y$ , $Z$ , $B$ and $\unicode[STIX]{x1D6FD}$ have the same form. The expansions (2.9)–(2.12) are then substituted into (A 2)–(A 4), (A 19) and (2.7)–(2.8), and terms at each order in $r/{\mathcal{R}}$ are collected. We can thereby relate the surface shape coefficients $(X,Y,Z)$ to the field strength through a desired order in $r/{\mathcal{R}}$ . Explicit results through $O((r/{\mathcal{R}})^{2})$ are given in (A 20)–(A 52).

Quasisymmetry is achieved through order $O((r/{\mathcal{R}})^{j})$ when $\unicode[STIX]{x2202}B_{k}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}=0$ for $k\leqslant \,j$ . At $O((r/{\mathcal{R}})^{1})$ , the analysis in appendix A shows quasisymmetric fields are described by

where $s_{\unicode[STIX]{x1D713}}=\text{sign}(\unicode[STIX]{x1D713})$ , $s_{G}=\text{sign}(G_{0})$ , $\bar{\unicode[STIX]{x1D702}}$ is a constant and $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D711})$ is a solution of

Here, the constant
$I_{2}$
is the leading term in the coefficient
$I(r)$
of (2.1), and is proportional to the on-axis toroidal current density, which is typically zero. The constant
$\bar{\unicode[STIX]{x1D702}}=B_{1c}/B_{0}$
(introduced by Garren & Boozer (Reference Garren and Boozer1991*b*
)) reflects the magnitude by which
$B$
varies on surfaces,

The surface shapes corresponding to (2.13) are ellipses in the plane perpendicular to the axis (the $\boldsymbol{n}$ – $\boldsymbol{b}$ plane). The ellipses are centred on the magnetic axis, so there is no Shafranov shift to this order.

It can be proved (Landreman *et al.*
Reference Landreman, Sengupta and Plunk2019) that quasisymmetric fields to
$O(r/{\mathcal{R}})$
can be parameterized by the shape of the axis – which can be any curve for which
$\unicode[STIX]{x1D705}$
never vanishes – along with three real numbers:
$I_{2}$
,
$\bar{\unicode[STIX]{x1D702}}$
and
$\unicode[STIX]{x1D70E}(0)$
. Varying
$\bar{\unicode[STIX]{x1D702}}$
has the effect of varying the average elongation. The parameter
$\unicode[STIX]{x1D70E}(0)$
is the value of
$\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D711})$
at
$\unicode[STIX]{x1D711}=0$
, and it reflects the angle by which the major and minor axes of the elliptical flux surfaces are oriented with respect to
$\boldsymbol{n}$
at
$\unicode[STIX]{x1D711}=0$
. Stellarators typically possess ‘stellarator symmetry’ (unrelated to quasisymmetry) which implies
$\unicode[STIX]{x1D70E}(0)=0$
. The pressure profile does not appear in the model to this order.

Proceeding to
$O((r/{\mathcal{R}})^{2})$
, nine new functions of
$\unicode[STIX]{x1D711}$
arise in the surface shapes:
$X_{20}$
,
$X_{2s}$
,
$X_{2c}$
,
$Y_{20}$
,
$Y_{2s}$
,
$Y_{2c}$
,
$Z_{20}$
,
$Z_{2s}$
and
$Z_{2c}$
. The corresponding surface shapes can now possess triangularity, and
$X_{20}$
and
$Y_{20}$
enable the centre of the surfaces at any given
$(r,\unicode[STIX]{x1D711})$
to be shifted in the
$\boldsymbol{n}$
–
$\boldsymbol{b}$
plane with respect to the axis (Shafranov shift). These nine functions are constrained by 10 new
$\unicode[STIX]{x1D711}$
-dependent equations: (A 27)–(A 29), (A 32)–(A 36) and (A 41)–(A 42) (with
$X_{1s}=\unicode[STIX]{x1D6FD}_{0}=\unicode[STIX]{x1D6FD}_{1c}=0$
as explained in A.3). Therefore, as noted by Garren & Boozer (Reference Garren and Boozer1991*b*
,Reference Garren and Boozer
*a*
), there is a net loss of one
$\unicode[STIX]{x1D711}$
-dependent degree of freedom, and so most axis shapes are not consistent with quasisymmetry through this order, in contrast to the
$O(r/{\mathcal{R}})$
case.

Also, four new scalar parameters appear at $O((r/{\mathcal{R}})^{2})$ that are independent of $\unicode[STIX]{x1D711}$ . One is $p_{2}$ , providing the first information about the pressure profile. Also appearing are the numbers $B_{20}$ , $B_{2s}$ and $B_{2c}$ , describing the variation of $B$ with $\unicode[STIX]{x1D717}$ . The global magnetic shear does not yet enter the system of equations; it first appears at $O((r/{\mathcal{R}})^{3})$ .

## 3 Generating a finite-minor-radius boundary

Our goal is ultimately to construct the shape of a boundary magnetic surface such that the magnetic field in the interior is quasisymmetric to a good approximation. The interest in constructing boundary surfaces comes from the fact that magnetohydrodynamic (MHD) equilibrium codes such as VMEC naturally take the shape of a boundary magnetic surface as an input. Conventional stellarator optimization codes such as STELLOPT and ROSE are built upon VMEC, so if we can construct a boundary surface for a quasisymmetric configuration, then we can construct a good initial condition for optimization.

Given a solution of the equations of the near-axis expansion through some order, how can a boundary surface be generated? A natural approach is to plug a small but finite value
$a$
of the expansion parameter
$r$
into the series for
$(X,Y,Z)$
, which yields a toroidal surface with average minor radius
$a$
. This approach was applied successfully to the
$O((r/{\mathcal{R}})^{1})$
equations in Landreman *et al.* (Reference Landreman, Sengupta and Plunk2019) and Plunk *et al.* (Reference Plunk, Landreman and Helander2019). However, this approach of setting
$r\rightarrow a$
turns out to require modification when applied to the
$O((r/{\mathcal{R}})^{2})$
equations. It can be seen that some care is required when setting
$r$
to a finite value, as follows. When the expansion is truncated at a finite order in
$r/{\mathcal{R}}$
, a finite value of
$r$
chosen and MHD equilibrium computed within the resulting surface, a configuration results that has a slightly different axis shape and
$(X,Y,Z)$
coefficients than the ones assumed in the original expansion. The difference turns out to be unimportant for the
$O((r/{\mathcal{R}})^{1})$
construction but critical for the
$O((r/{\mathcal{R}})^{2})$
construction. The fact that plugging in a finite value of
$r$
to obtain a boundary results in a slight change to the equilibrium can be seen in figures 3, 8 and 10 of Landreman *et al.* (Reference Landreman, Sengupta and Plunk2019). These figures, showing the
$B$
spectrum when a finite
$r$
is substituted into the
$O((r/{\mathcal{R}})^{1})$
Garren–Boozer expansion and the equilibrium is computed inside the resulting boundary, show that there is a small but finite mirror mode on the magnetic axis, even though the on-axis mirror mode amplitude is precisely zero in the original expansion.

In this section, we introduce a systematic method to examine the effect of setting the expansion parameter $r$ equal to a finite number $a$ . Technical details of the calculation are given in appendix B. Here, and in the appendix, we consider a more general problem of trying to construct a configuration with any desired field strength $B(r,\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})$ , which may or may not be quasisymmetric; therefore, the analysis also applies to more general optimizations such as omnigenity. The basic approach is to introduce a second Garren–Boozer-type expansion, denoted with tildes, that describes the configuration which results from computing an MHD equilibrium inside the boundary constructed from the original Garren–Boozer expansion. In contrast to the original expansion, the ‘tilde’ expansion is not truncated, since it describes MHD equilibrium in a finite volume. The axis shapes and $(X,Y,Z)$ coefficients for the tilde and non-tilde expansions are similar but not identical. Their differences diminish as $a\rightarrow 0$ . The non-tilde expansion represents a single idealized configuration we would like to achieve, whereas the tilde expansion represents a family of different ‘real’ configurations, (real in the sense that they are what is computed by solving for MHD equilibrium without a near-axis expansion), parameterized by the finite value $a$ used to construct their boundaries. From the fact that the (truncated) non-tilde expansion and (non-truncated) tilde expansions describe the same surface at $r=a$ , and exploiting an expansion in $r/{\mathcal{R}}\ll 1$ with the ordering $a\sim r$ , we can derive the magnetic field strength that results for the constructed configurations. We will thereby rigorously show that substituting $r\rightarrow a$ in an $O(r/{\mathcal{R}})$ Garren–Boozer solution yields a configuration that has the desired $B$ to $O(r/{\mathcal{R}})$ , but substituting $r\rightarrow a$ in an $O((r/{\mathcal{R}})^{2})$ Garren–Boozer solution yields a configuration that only has the desired $B$ through $O(r/{\mathcal{R}})$ , not $O((r/{\mathcal{R}})^{2})$ . However, the achieved $B$ can be made to match the desired one at $O((r/{\mathcal{R}})^{2})$ by a small modification of the construction, in which $X_{3}$ and $Y_{3}$ terms are included.

Beginning the formal analysis, we consider a family of equilibria parameterized by $a$ , in which the position vector is

analogous to (2.4). The Boozer angles $(\tilde{\unicode[STIX]{x1D717}},\tilde{\unicode[STIX]{x1D711}})$ for the real configuration generally differ somewhat from the angles $(\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})$ of the ideal configuration, with the differences denoted by single-valued functions $t$ and $f$ ,

*a*,

*b*) $$\begin{eqnarray}\tilde{\unicode[STIX]{x1D717}}(a,\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})=\unicode[STIX]{x1D717}+t(a,\unicode[STIX]{x1D717},\unicode[STIX]{x1D711}),\quad \tilde{\unicode[STIX]{x1D711}}(a,\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})=\unicode[STIX]{x1D711}+f(a,\unicode[STIX]{x1D717},\unicode[STIX]{x1D711}).\end{eqnarray}$$

Analogous to (2.9), we have

with similar expansions for ${\tilde{Y}}$ and $\tilde{Z}$ , and analogous to (2.9), the field strength in the real configurations is

with a similar expansion for $\tilde{\unicode[STIX]{x1D6FD}}$ . All quantities in the tilde configurations are assumed to have $a$ -dependence in the form of a Taylor series, with coefficients denoted by superscripts in parentheses,

with analogous expansions for $\tilde{\boldsymbol{n}}$ , $\tilde{\boldsymbol{b}}$ and $\tilde{\boldsymbol{t}}$ , and

with analogous expansions for ${\tilde{Y}}_{j}$ , $\tilde{Z}_{j}$ , $\tilde{B}_{j}$ and $\tilde{\unicode[STIX]{x1D6FD}}_{j}$ . We similarly assume

*a*,

*b*) $$\begin{eqnarray}t(a,\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})=\mathop{\sum }_{k=0}^{\infty }a^{k}t^{(k)}(\unicode[STIX]{x1D717},\unicode[STIX]{x1D711}),\quad f(a,\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})=\mathop{\sum }_{k=0}^{\infty }a^{k}f^{(k)}(\unicode[STIX]{x1D717},\unicode[STIX]{x1D711}).\end{eqnarray}$$

To reiterate, subscripts refer to an expansion in distance from the axis in a fixed configuration, whereas superscripts in parentheses indicate a distinct expansion in the finite value of minor radius substituted into the original near-axis expansion.

The boundary of a tilde configuration, i.e. its $r=a$ surface, by definition is the surface obtained by setting $r=a$ in the non-tilde expansion. The equation representing this fact is

and it plays a central role in the analysis.

To the order of interest, the field strength in the real configurations is

The quantities in this expression are computed in terms of the non-tilde expansion in appendix B, by systematically examining (3.8) at each order in
$a/{\mathcal{R}}\sim r/{\mathcal{R}}$
. There, assuming only that the construction is carried out through
$(X_{1},Y_{1},Z_{1})$
or higher, it is found that
$\tilde{B}_{0}^{(0)}(\unicode[STIX]{x1D711})=B_{0}(\unicode[STIX]{x1D711})$
,
$\tilde{B}_{1}^{(0)}(\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})=B_{1}(\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})$
and
$\tilde{B}_{0}^{(1)}(\unicode[STIX]{x1D711})=0$
. Hence, if a finite value
$a$
is substituted into
$r$
in a
$O((r/{\mathcal{R}})^{1})$
Garren–Boozer solution to construct a boundary surface, the real configuration inside this boundary will have the desired magnetic field strength in Boozer coordinates through
$O((r/{\mathcal{R}})^{1})$
. This finding is consistent with the results in Landreman *et al.* (Reference Landreman, Sengupta and Plunk2019). However, the results at next order are more complicated. Assuming now that the construction is carried out through
$(X_{2},Y_{2},Z_{2})$
or higher, it is found in appendix B that
$\tilde{B}_{2}^{(0)}(\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})=B_{2}(\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})$
,
$\tilde{B}_{1}^{(1)}(\unicode[STIX]{x1D717},\unicode[STIX]{x1D711})=0$
, and

where

primes denote $\text{d}/\text{d}\unicode[STIX]{x1D711}$ , and

Thus, if the Garren–Boozer solution is carried out through $(X_{2},Y_{2},Z_{2})$ but then truncated so $(X_{3},Y_{3},Z_{3})$ are all set to zero when constructing the boundary, (3.11) will be non-zero. The $a^{2}\tilde{B}_{0}^{(2)}$ term in (3.9) will then be non-zero and will cause a difference between the desired and achieved field strengths that is comparable to the desired $r^{2}B_{2}$ term. Generally (3.10) will depend on $\unicode[STIX]{x1D711}$ and so it will break quasisymmetry.

Fortunately, a workaround can be achieved that does not require a full solution of the Garren–Boozer equations for $(X_{3},Y_{3},Z_{3})$ . It can be preferable to avoid a full solution through $O((r/{\mathcal{R}})^{3})$ because the equations are extremely complicated, because one then has to choose additional parameters for the construction and because the presence of squareness that grows with $r$ at this order can limit the minimum aspect ratio. In the workaround, we take $(X_{3},Y_{3})$ to be $(X_{1},Y_{1})$ scaled by some function $\unicode[STIX]{x1D706}(\unicode[STIX]{x1D711})$ ,

*a*-

*c*) $$\begin{eqnarray}X_{3c1}=\unicode[STIX]{x1D706}X_{1c},\quad X_{3s1}=\unicode[STIX]{x1D706}X_{1s},\quad X_{3c3}=X_{3sc}=0,\end{eqnarray}$$

with analogous expressions for $Y$ and $Z_{3}=0$ . In other words, we introduce a small correction to the $O(r/{\mathcal{R}})$ elliptical flux surface shape. Setting (3.11) $=0$ , substituting (3.14) and using (A 21), we find

Adding these $(X_{3},Y_{3})$ terms to the constructed boundary surface therefore results in $\tilde{B}_{0}^{(2)}=0$ , so the real configurations have the same Boozer spectrum as the ideal target configuration through $O((r/{\mathcal{R}})^{2})$ .

Some physical intuition for this result can be given. The leading-order field strength $B_{0}$ is approximately the toroidal flux divided by the cross-sectional area of the flux surfaces. Indeed, this interpretation of (A 21) is shown precisely in Landreman & Sengupta (Reference Landreman and Sengupta2018). The area of the surfaces is primarily determined by the $\sin \unicode[STIX]{x1D717}$ and $\cos \unicode[STIX]{x1D717}$ modes of $X$ and $Y$ , which generate ellipses, and not by $\sin 2\unicode[STIX]{x1D717}$ , $\cos 2\unicode[STIX]{x1D717}$ , and independent-of- $\unicode[STIX]{x1D717}$ modes, which distort and shift the ellipses but do not expand or contract them. The $\sin \unicode[STIX]{x1D717}$ and $\cos \unicode[STIX]{x1D717}$ modes of $X$ and $Y$ that affect the cross-sectional area arise at orders $O(r/{\mathcal{R}})$ and $O((r/{\mathcal{R}})^{3})$ but not at $O((r/{\mathcal{R}})^{2})$ , due to analyticity. Thus, if the Garren–Boozer solution is truncated by setting $X_{3}=Y_{3}=Z_{3}=0$ , there is an $O((r/{\mathcal{R}})^{2})$ error in the cross-sectional area of the surfaces, which (since $B\sim \text{flux}/\text{area}$ ) implies an $O((r/{\mathcal{R}})^{2})$ error in $B_{0}$ . This error can vary toroidally, spoiling quasisymmetry or whatever other pattern of $B$ is desired. To solve the problem, we note (A 49) is a correction to (A 21) (they both derive from (A 4)), relating the cross-sectional area and $B$ to flux. Thus, (A 49) indicates how much the surfaces should be expanded or contracted to give the correct $B_{0}$ through $O((r/{\mathcal{R}})^{2})$ . Indeed, the same result (3.15) can be obtained by setting $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D713}=\int \text{d}^{2}\boldsymbol{a}\boldsymbol{\cdot }\boldsymbol{B}$ at $O((r/{\mathcal{R}})^{2}r^{2}B)$ , where $\int \text{d}^{2}\boldsymbol{a}$ is an integral over a constant- $\unicode[STIX]{x1D711}$ cross-section, using (2.3) and (3.14).

Note that by a modified choice of
$\unicode[STIX]{x1D706}$
,
$\tilde{B}_{0}^{(2)}$
can be made to cancel the toroidal dependence of
$\tilde{B}_{20}^{(0)}$
at a non-zero value
$r_{0}$
of
$r$
. In the case of quasisymmetry, such a choice has the effect of introducing an
$O((r/{\mathcal{R}})^{2})$
mirror mode on the axis, with the mirror mode amplitude vanishing at radius
$r_{0}$
. There may be advantages in this approach, for as found in a recent numerical study (Henneberg *et al.*
Reference Henneberg, Drevlak, Nührenberg, Beidler, Turkin, Loizu and Helander2019), fast particle confinement in a quasi-axisymmetric configuration was best when the quasisymmetry was optimized off axis rather than on axis.

## 4 Numerical formulation

### 4.1 Inputs and outputs

We now describe a practical numerical implementation of the equations derived in the preceding sections and associated appendices. The parameters of the algorithm here are a superset of the inputs to the
$O((r/{\mathcal{R}})^{1})$
construction detailed in Landreman *et al.* (Reference Landreman, Sengupta and Plunk2019). The latter are the shape of the magnetic axis, which must have non-vanishing curvature everywhere, and the numbers
$\bar{\unicode[STIX]{x1D702}}$
,
$I_{2}$
and
$\unicode[STIX]{x1D70E}(0)$
. The parameter
$\bar{\unicode[STIX]{x1D702}}$
effectively controls the elongation;
$I_{2}$
indicates the on-axis toroidal current and is typically zero for stellarators;
$\unicode[STIX]{x1D70E}(0)$
controls the angle of elongation at
$\unicode[STIX]{x1D711}=0$
and is zero for stellarator-symmetric configurations. In the
$O((r/{\mathcal{R}})^{2})$
construction, three new constant input parameters are needed:
$p_{2}$
,
$B_{2c}$
and
$B_{2s}$
. We will not take
$B_{20}$
as an input parameter for reasons explained shortly. The parameter
$p_{2}$
defines the pressure profile to this order. The parameters
$B_{2c}$
and
$B_{2s}$
set the desired
$\cos 2\unicode[STIX]{x1D717}$
and
$\sin 2\unicode[STIX]{x1D717}$
modes in the field strength. These two parameters have the effect of controlling the stellarator-symmetric and stellarator-asymmetric parts of the triangularity. For stellarator-symmetric configurations,
$B_{2s}=0$
.

The outputs of the calculation include the shapes of the magnetic surfaces and the rotational transform on axis,
$\unicode[STIX]{x1D704}_{0}$
. As noted by Garren & Boozer (Reference Garren and Boozer1991*a*
,Reference Garren and Boozer
*b*
) and discussed above, the number of scalar
$\unicode[STIX]{x1D711}$
-dependent unknowns exceeds the number of
$\unicode[STIX]{x1D711}$
-dependent equations by only one if quasisymmetry is imposed through
$O((r/{\mathcal{R}})^{2})$
, whereas an axis shape represents two
$\unicode[STIX]{x1D711}$
-dependent quantities (e.g.
$\unicode[STIX]{x1D705}$
and
$\unicode[STIX]{x1D70F}$
, or
$R(\unicode[STIX]{x1D719})$
and
$z(\unicode[STIX]{x1D719})$
, where
$(R,\unicode[STIX]{x1D719},z)$
are cylindrical coordinates.) Therefore, to achieve quasisymmetry through
$O((r/{\mathcal{R}})^{2})$
, one needs to solve for part of the axis shape. We proceed by temporarily relaxing the requirement that
$B_{20}$
must be independent of
$\unicode[STIX]{x1D711}$
. By reducing the number of equations by one in this way, any axis shape becomes allowed. One can still make
$B_{2s}$
and
$B_{2c}$
independent of
$\unicode[STIX]{x1D711}$
, achieving quasisymmetry partially through
$O((r/{\mathcal{R}})^{2})$
. Then
$B_{20}(\unicode[STIX]{x1D711})$
is an output of the calculation, and it generally has some toroidal variation. We can then numerically optimize the input parameters (including not only the axis shape but also
$\{\bar{\unicode[STIX]{x1D702}},\unicode[STIX]{x1D70E}(0),B_{2c},B_{2s}\}$
) such that the toroidal variation of
$B_{20}$
is minimal. While we have not proved rigorously that solutions exist in which
$B_{20}$
is exactly constant, experience with the numerical solutions described in the next section suggests
$B_{20}^{\prime }$
can be made arbitrarily small as the number of degrees of freedom in the axis shape is increased.

As explained in detail in § 5.2 of Landreman & Sengupta (Reference Landreman and Sengupta2018), quasi-axisymmetry versus quasi-helical symmetry is determined by the choice of axis shape. In particular, the integer $N$ is the number of times the axis normal vector rotates poloidally around the axis as the axis is traversed toroidally.

### 4.2 Numerical solution of the equations

Given the input parameters described in the previous subsection, the
$O(r/{\mathcal{R}})$
equations are solved as described in § 3 of Landreman *et al.* (Reference Landreman, Sengupta and Plunk2019). As a result,
$X_{1}$
and
$Y_{1}$
are computed on a uniform grid in the standard toroidal angle
$\unicode[STIX]{x1D719}$
covering one field period with
$N_{\unicode[STIX]{x1D719}}$
points, and
$\unicode[STIX]{x1D704}_{0}$
is obtained. Then,
$Z_{2}$
is computed from (A 27)–(A 29),
$X_{2s}$
is computed from (A 35) and
$X_{2c}$
is computed from (A 36).

Next, a
$(2N_{\unicode[STIX]{x1D719}})\times (2N_{\unicode[STIX]{x1D719}})$
linear system is solved. The unknowns for this system are the values of
$X_{20}$
and
$Y_{20}$
on the
$N_{\unicode[STIX]{x1D719}}$
grid points. The rows of the linear system represent (A 41) and (A 42) imposed at the
$N_{\unicode[STIX]{x1D719}}$
grid points. In these equations,
$\text{d}/\text{d}\unicode[STIX]{x1D711}$
derivatives are discretized using the same pseudospectral differentiation matrix described in Landreman *et al.* (Reference Landreman, Sengupta and Plunk2019). In this system,
$Y_{2s}$
and
$Y_{2c}$
are eliminated using (A 32)–(A 33). The dense linear system is solved with direct factorization (LAPACK). With
$X_{20}$
and
$Y_{20}$
thereby determined,
$Y_{2s}$
and
$Y_{2c}$
are computed from (A 32)–(A 33), and then (A 34) gives
$B_{20}$
. Finally, (3.14)–(3.15) give
$X_{3}$
and
$Y_{3}$
.

Note that although the
$O((r/{\mathcal{R}})^{1})$
equations are nonlinear in the unknowns, a unique solution always exists, as proved in Landreman *et al.* (Reference Landreman, Sengupta and Plunk2019), so the solution by Newton’s method is extremely robust. Furthermore, once the
$O((r/{\mathcal{R}})^{1})$
solution is determined, the equations of § A.2 are linear in the higher-order unknowns, so the
$O((r/{\mathcal{R}})^{2})$
construction is equally robust. At a typical resolution
$(N_{\unicode[STIX]{x1D719}}\sim 30)$
, solution of the equations for the
$O((r/{\mathcal{R}})^{2})$
construction takes
${<}2~\text{ms}$
on one core of a modern laptop, many orders of magnitude lower computational cost than a general 3-D equilibrium calculation used in each iteration of traditional stellarator optimization.

### 4.3 Optimization of input parameters

For many sets of input parameters, the model results in configurations that are not of practical interest because they are limited to extremely high aspect ratio. This limitation arises because for any solution of the model equations, beyond a certain value of $r$ , the constant- $r$ surfaces will begin to self-intersect. If $X_{2}$ , $Y_{2}$ , $X_{3}$ or $Y_{3}$ is large, this critical $r$ will be small. From another perspective, the near-axis expansion is only accurate at values of $r$ sufficiently small that the terms of successive orders in $r$ in the expansion are decaying. If $X_{2}$ , $Y_{2}$ , $X_{3}$ or $Y_{3}$ is large, the accuracy of the expansion is then limited to smaller values of $r$ .

Therefore, for some of the examples below, we use optimization – either numerical or by hand – to find solutions with small
$X_{2}$
,
$Y_{2}$
,
$X_{3}$
or
$Y_{3}$
. We also minimize the toroidal derivatives of these quantities, anticipating that large derivatives would drive large symmetry-breaking terms at next order. All the quantities targeted for minimization are squared, averaged over
$\unicode[STIX]{x1D719}$
and combined in a weighted sum to form a single objective function. For some examples, to improve quasisymmetry, we also include in the sum a term minimizing the toroidal variation of
$B_{20}$
. While doing these optimizations, it may be necessary to penalize parameters for which
$\unicode[STIX]{x1D704}_{0}$
becomes too small. Note that when optimization is applied to the
$O((r/{\mathcal{R}})^{2})$
quasisymmetry model, the objective function can be evaluated in milliseconds, approximately four orders of magnitude or more faster than the objective function evaluations in conventional stellarator optimization. Also, in principle, analytic derivatives are available for the
$O((r/{\mathcal{R}})^{2})$
model, although we will not exploit them here to accelerate optimization. For the examples below, we use Matlab’s derivative-free algorithm ‘fminsearch’, a variant of the Nelder–Mead simplex algorithm by Lagarias *et al.* (Reference Lagarias, Reeds, Wright and Wright1998).

### 4.4 Conversion to cylindrical coordinates

A principal aim of the construction is to generate boundary shapes that can be provided as input to an MHD equilibrium code such as VMEC. VMEC requires as input the boundary surface shape defined by its cylindrical coordinates
$(R,z)$
expressed as a double Fourier series in the toroidal angle
$\unicode[STIX]{x1D719}$
and any poloidal angle. As discussed in § 4 of Landreman *et al.* (Reference Landreman, Sengupta and Plunk2019), there are several ways this boundary description can be obtained from our representation (2.4). One approach is to derive the transformation between the two representations order by order in
$a$
. This approach was developed in § 4.1 of Landreman & Sengupta (Reference Landreman and Sengupta2018) to
$O(r/{\mathcal{R}})$
. In appendix C, the transformation is extended to
$O((r/{\mathcal{R}})^{3})$
. The advantage of this approach is that it results in explicit expressions for
$R(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$
and
$z(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$
that can be evaluated extremely rapidly.

A second approach to transforming from the Frenet representation to the representation required by VMEC is described in § 4.2 of Landreman *et al.* (Reference Landreman, Sengupta and Plunk2019). In this approach, nonlinear root finding is applied to (2.4). At a grid of points in
$\unicode[STIX]{x1D703}$
and
$\unicode[STIX]{x1D719}$
, one solves for the value of
$\unicode[STIX]{x1D711}$
such that the position vector has a standard toroidal angle
$\unicode[STIX]{x1D719}$
. The nonlinear root finding requires additional computation time. This approach tends to result in slightly lower magnitude of symmetry breaking, so we use it for the numerical examples that follow.

## 5 Numerical results

Several examples of constructed quasi-axisymmetric and quasi-helically symmetric configurations will now be presented. The examples are all generated ‘from scratch’, in that no fitting was done to previously optimized equilibria. All the examples are scaled such that the zero-frequency component of the axis major radius $R_{00}=(2\unicode[STIX]{x03C0})^{-1}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D719}\,R_{0}$ is 1 metre, and the on-axis field strength $B_{0}$ is 1 Tesla. In each VMEC calculation shown, the pressure profile specified was $p(r)=(1-r^{2}/a^{2})p_{2}$ for the same constant $p_{2}$ used in the construction. Also, the current profile for VMEC calculations was specified as ${\mathcal{I}}(s)=2\unicode[STIX]{x03C0}sa^{2}I_{2}/\unicode[STIX]{x1D707}_{0}$ , where ${\mathcal{I}}(s)$ is the toroidal current inside the surface with normalized toroidal flux $s=(r/a)^{2}$ , and $I_{2}$ is the constant used in the construction.

We will describe the configurations using two different measures of effective aspect ratio. The measure that is most convenient for the construction is $A=R_{00}/a$ , where again $\unicode[STIX]{x03C0}a^{2}$ is the toroidal flux. We will also quote the effective aspect ratio $A_{vmec}$ used in the VMEC code, since this measure is often reported in the literature. Its definition is $A_{vmec}=(\mathtt{Aminor\_p})/(\mathtt{Rmajor\_p})$ , where the effective minor radius $\mathtt{Aminor\_p}$ is defined by $\unicode[STIX]{x03C0}(\mathtt{Aminor\_p})^{2}=\bar{S}$ , with $\bar{S}=(2\unicode[STIX]{x03C0})^{-1}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D719}\,S(\unicode[STIX]{x1D719})$ the toroidal average of the area $S(\unicode[STIX]{x1D719})$ of the outer surface’s cross-section in the $R$ – $z$ plane, and the effective major radius $\mathtt{Rmajor\_p}$ is defined by $[2\unicode[STIX]{x03C0}(\mathtt{Rmajor\_p})][\unicode[STIX]{x03C0}(\mathtt{Aminor\_p})^{2}]=V$ with $V$ the volume of the outer surface.

### 5.1 Quasi-axisymmetry partially through $O((r/{\mathcal{R}})^{2})$

The first set of input parameters we consider includes the axis shape

corresponding to two field periods. We also choose
$\bar{\unicode[STIX]{x1D702}}=0.640~\text{m}^{-1}$
and
$B_{2c}=-0.00322~\text{T}~\text{m}^{-2}$
. These values and axis shape were obtained by minimizing
$X_{2}$
,
$Y_{2}$
,
$X_{3}$
and
$Y_{3}$
, as described in § 4.3, subject to a lower bound
$\unicode[STIX]{x1D704}_{0}\geqslant 0.42$
. The parameters
$\unicode[STIX]{x1D70E}(0)$
and
$B_{2s}$
were set to zero so the configuration is stellarator symmetric. We also choose
$I_{2}=0$
and
$p_{2}=0$
so the configuration is a vacuum field. For this first configuration, no attempt was made to make
$B_{20}$
independent of
$\unicode[STIX]{x1D711}$
, so the configuration is only partially quasisymmetric at
$O((r/{\mathcal{R}})^{2})$
. The resulting configuration has
$\unicode[STIX]{x1D704}_{0}=0.420$
as desired, and the boundary shape for
$A=10$
is shown in figure 1. VMEC is then run to compute the equilibrium inside the finite-thickness boundary without making any near-axis approximation, and then the BOOZ_XFORM code (Sanchez *et al.*
Reference Sanchez, Hirshman, Ware, Berry and Spong2000) is run to transform the VMEC result to Boozer coordinates. Figure 2 shows the resulting Fourier coefficients
$B_{m,n}(r)$
defined by
$B(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})=\sum _{m,n}B_{m,n}(r)\cos (m\unicode[STIX]{x1D703}-n\unicode[STIX]{x1D711})$
. It can be seen that the dominant
$B_{m,n}$
mode is the quasi-axisymmetric term
$B_{1,0}$
as desired. The magnitude of this mode predicted by the construction,
$B_{1,0}=r\bar{\unicode[STIX]{x1D702}}B_{0}$
, is displayed for comparison, and it is nearly identical to the VMEC result.

Keeping all input parameters fixed except for the boundary aspect ratio, we then construct boundary surfaces at a sequence of increasing aspect ratios, and repeat the VMEC and BOOZ_XFORM calculations for each case. In figure 3, the quantity $[B_{m=0}(\unicode[STIX]{x1D711},r=a)-B(\unicode[STIX]{x1D711},r=0)]/a^{2}$ is displayed for this sequence of numerical calculations. According to the construction, this quantity should be $B_{20}(\unicode[STIX]{x1D711})$ . Indeed, it can be seen that the VMEC/BOOZ_XFORM results converge to the Garren–Boozer prediction. Thus, in the limit $A\gg 1$ , the full 3-D equilibrium calculations achieve the expected field strength in Boozer coordinates.

As another verification of the construction, figure 4 displays three symmetry-breaking measures

computed from the VMEC and BOOZ_XFORM results for the aspect ratio scan. (‘Config 1’ in the figure refers to the present section, while ‘Config 2’ will be described in § 5.2.) It can be seen that $S_{m>0}^{r=a}$ scales as $1/A^{3}$ , consistent with the fact that the corresponding modes were constructed to be zero through $O((r/{\mathcal{R}})^{2})$ , so symmetry breaking generally arises at next order. The on-axis mirror modes, measured by $S_{m=0}^{r=0}$ , are found to have an even stronger scaling, $\propto 1/A^{4}$ . This scaling arises because the $m=0$ modes were constructed to be zero through $O((r/{\mathcal{R}})^{2})$ , and they are automatically zero at $O((r/{\mathcal{R}})^{3})$ (only $m=1$ and $m=3$ modes exist at this order), so the first non-vanishing contribution occurs at $O((r/{\mathcal{R}})^{4})$ . Finally, $S_{m=0}^{r=a}$ shows an asymptotic scaling $\propto 1/A^{2}$ , associated with the fact that $B_{20}$ is not independent of $\unicode[STIX]{x1D711}$ . Thus, all three symmetry-breaking measures scale as predicted by the construction.

### 5.2 Quasi-axisymmetry fully through $O((r/{\mathcal{R}})^{2})$

We next consider a configuration that is similar to the one of § 5.1, but with slightly different parameters such that $B_{20}$ has significantly reduced toroidal variation, resulting in improved quasi-axisymmetry. We again consider a two-field-period device, with axis shape

The other non-zero input parameters are $\bar{\unicode[STIX]{x1D702}}=0.632~\text{m}^{-1}$ and $B_{2c}=-0.158~\text{T}~\text{m}^{-2}$ . These values and axis shape were obtained by the optimization procedure of § 4.3, again minimizing $X_{2}$ , $Y_{2}$ , $X_{3}$ and $Y_{3}$ , but now also minimizing the toroidal variation of $B_{20}$ . The parameters $\unicode[STIX]{x1D70E}(0)$ and $B_{2s}$ were again set to zero so the configuration is stellarator symmetric, and $I_{2}$ and $p_{2}$ were set to zero so the configuration is a vacuum field. The resulting configuration has $\unicode[STIX]{x1D704}_{0}=0.424$ . The function $B_{20}(\unicode[STIX]{x1D711})$ for these parameters is shown as the black dotted curve in figure 7, and it can be seen that the toroidal variation is greatly reduced compared to figure 3. The small remaining toroidal variation of $B_{20}$ could presumably be further reduced if additional Fourier modes were included in the axis shape. The constructed boundary shape for $A=10$ is shown in figure 5, and it is only slightly different from the previous example (figure 1.) Running VMEC and BOOZ_XFORM inside this boundary results in the magnetic spectrum of figure 6. Again, the desired mode $B_{1,0}$ dominates, and its magnitude is nearly identical to the prediction. Figure 7 shows that $[B_{m=0}(\unicode[STIX]{x1D711},r=a)-B(\unicode[STIX]{x1D711},r=0)]/a^{2}$ again converges to the predicted function, $B_{20}(\unicode[STIX]{x1D711})$ .

The three symmetry-breaking measures for this second configuration are displayed in figure 4, labelled as ‘Config 2’. It can be seen that $S_{m>0}^{r=a}$ and $S_{m=0}^{r=0}$ are not much changed from the first configuration, scaling as $1/A^{3}$ and $1/A^{4}$ , as before. However, $S_{m=0}^{r=a}$ is significantly changed, still scaling as $1/A^{2}$ for sufficiently large $A$ , but with the leading constant reduced by over an order of magnitude. This reduction reflects the reduced variation of $B_{20}$ . For $A<40$ , $S_{m=0}^{r=a}$ now scales as $1/A^{4}$ since it is dominated by the on-axis variation of $B$ measured by $S_{m=0}^{r=0}$ . Thus, $B_{20}$ is constant enough that it is not the dominant source of symmetry breaking for the entire range of aspect ratios shown, $A\in [5,320]$ . This point is apparent also in figure 8, in which the quantity plotted

includes all quasisymmetry-breaking modes. This figure more clearly shows the scaling predicted by Garren & Boozer (Reference Garren and Boozer1991*a*
) that the total deviation from quasisymmetry can be made to scale as
$\propto 1/A^{3}$
. The quasisymmetry of this quasi-axisymmetric configuration is sufficiently good that, for
$A\geqslant 60$
, the deviation from quasisymmetry
$S_{\text{tot}}B_{0}$
is smaller than the Earth’s magnetic field of
${\sim}0.5$
Gauss. At the right-most point (
$A=320$
),
$S_{\text{tot}}$
is
${<}4\times 10^{-7}$
, and the largest single symmetry-breaking Fourier mode has an amplitude
${<}2\times 10^{-7}~\text{T}$
.

### 5.3 Tokamak–stellarator hybrid

To verify the construction for a case in which the plasma pressure and on-axis current are non-zero, we next consider a tokamak–stellarator hybrid configuration, in which both non-axisymmetric shaping and toroidal current contribute to the rotational transform. We again consider a two-field-period geometry, with axis shape

The parameters $\unicode[STIX]{x1D70E}(0)$ and $B_{2s}$ were again set to zero so the configuration is stellarator symmetric. The other input parameters were $\bar{\unicode[STIX]{x1D702}}=0.95~\text{m}^{-1}$ , $I_{2}=0.9~\text{T}~\text{m}^{-1}$ , $p_{2}=-6\times 10^{5}~\text{Pa}~\text{m}^{-2}$ and $B_{2c}=-0.7~\text{T}~\text{m}^{-2}$ . For this value of $p_{2}$ , the volume-averaged $\unicode[STIX]{x1D6FD}$ (plasma pressure/magnetic pressure) for the configuration at $A=5$ is 2.9 %. The resulting configuration has $\unicode[STIX]{x1D704}_{0}=0.960$ . For comparison, a vacuum field inside the constructed $A=5$ boundary has an on-axis transform $\unicode[STIX]{x1D704}_{0}=0.214$ . This level of vacuum transform might be sufficient to provide stellarator-like stability. The boundary shape for $A=5$ is shown in figure 9. Figure 10 shows the Boozer spectrum of the finite- $\unicode[STIX]{x1D6FD}$ finite-current configuration inside this boundary. While the desired mode $B_{1,0}$ dominates, and it has a magnitude close to that predicted by the construction, the departures from quasisymmetry are larger than in the previous examples, associated with the smaller value of $A$ here. Figure 11 shows that as $A$ is increased, $[B_{m=0}(\unicode[STIX]{x1D711},r=a)-B(\unicode[STIX]{x1D711},r=0)]/a^{2}$ again converges to the predicted function, $B_{20}(\unicode[STIX]{x1D711})$ . The scaling of the three symmetry-breaking measures with $A$ is plotted in figure 12, and again, they scale as the expected power of $A$ . Together, figures 11–12 verify the $O((r/{\mathcal{R}})^{2})$ construction behaves correctly when $I_{2}$ and $p_{2}$ terms are included.

### 5.4 Quasi-helical symmetry

We next consider a quasi-helically symmetric configuration. The axis shape is taken to be

with $\bar{\unicode[STIX]{x1D702}}=1.569~\text{m}^{-1}$ and $B_{2c}=0.1348~\text{T}~\text{m}^{-2}$ . These values were obtained using the optimization procedure of § 4.3 to minimize $X_{2}$ , $Y_{2}$ , $X_{3}$ and $Y_{3}$ . For this axis shape, the normal vector rotates around the axis poloidally four times as the axis is traversed toroidally, so the construction yields quasi-helical symmetry rather than quasi-axisymmetry. The parameters $\unicode[STIX]{x1D70E}(0)$ and $B_{2s}$ were set to zero so the configuration is stellarator symmetric. The other input parameters were $I_{2}=0$ and $p_{2}=0$ . The resulting configuration has $\unicode[STIX]{x1D704}_{0}=1.14$ . The constructed boundary shape for $A=8$ is shown in figure 13.

Compared to the case of quasi-axisymmetry, for quasi-helical symmetry it seems relatively hard to find sets of input parameters for which $X_{2}$ , $Y_{2}$ , $X_{3}$ and $Y_{3}$ are acceptably small. If these quantities are not small, the boundary aspect ratio must be large, or else the symmetry-breaking errors tend to be large and the boundary surface may self-intersect. This challenge for finding good quasi-helically symmetric configurations likely arises from the fact that they require significant helical excursion of the axis, implying larger $\unicode[STIX]{x1D70F}$ and $\unicode[STIX]{x1D705}$ compared to quasi-axisymmetric configurations, which act to drive larger $X_{2}$ and $Y_{2}$ . The configuration in this section manages to have small values of $\{X_{2},Y_{2},X_{3},Y_{3}\}$ due to some delicate balances in the equations of appendix A. For instance, merely rounding the coefficients in the axis shape (5.6) to 3 digits of precision rather than 4 causes significant increases in $X_{3}$ and $Y_{3}$ that result in visible changes to the boundary shape.

As with the earlier configurations, VMEC and BOOZ_XFORM calculations for this quasi-helically symmetric configuration confirm that the desired field strength is produced. One aspect of this verification is shown in figure 14, which displays the Boozer spectrum inside the constructed $A=8$ boundary. This time the dominant mode is $B_{1,4}$ , and the magnitude of this mode matches the prediction $r\bar{\unicode[STIX]{x1D702}}B_{0}$ . Figure 15 shows that as $A\rightarrow \infty$ , $[B_{m=0}(\unicode[STIX]{x1D711},r=a)-B(\unicode[STIX]{x1D711},r=0)]/a^{2}$ again converges to the predicted function, $B_{20}(\unicode[STIX]{x1D711})$ . Figure 16 shows that $S_{m>0}^{r=a}$ , $S_{m=0}^{r=a}$ and $S_{m=0}^{r=0}$ scale approximately as expected ( $1/A^{3}$ , $1/A^{4}$ transitioning to $1/A^{2}$ at large $A$ , and $1/A^{4}$ ), as for the previous configurations. The total deviation from quasisymmetry $S_{\text{tot}}$ is also displayed in figure 8, demonstrating again Garren and Boozer’s predicted scaling. (The range of aspect ratios plotted differs from that for the quasi-axisymmetric configuration since at the highest $A$ , it is difficult to obtain converged values from VMEC for the very small symmetry-breaking modes.)

### 5.5 Testing all terms

For a final example, we present an example in which all the parameters of the near-axis model are non-zero. This example is limited to quite a large aspect ratio due to the large $X_{2}$ and $Y_{2}$ terms, and so is not interesting as an experimental design, but it is useful here as a challenging verification test. We choose the axis shape

which yields quasi-helical symmetry with $N=5$ . The other input parameters are chosen to be $\bar{\unicode[STIX]{x1D702}}=2.5~\text{m}^{-1}$ , $\unicode[STIX]{x1D70E}(0)=0.3$ , $I_{2}=1.6~\text{T}~\text{m}^{-1}$ , $p_{2}=-5\times 10^{6}~\text{Pa}~\text{m}^{-2}$ , $B_{2c}=1~\text{T}~\text{m}^{-2}$ and $B_{2s}=3~\text{T}~\text{m}^{-2}$ . Note that stellarator symmetry is broken both by the non-zero value of $\unicode[STIX]{x1D70E}(0)$ and of $B_{2s}$ . This resulting configuration has $\unicode[STIX]{x1D704}_{0}=0.829$ . The constructed boundary shape is shown in figure 17 for $A=40$ ( $A_{vmec}=28.5$ ), and it can be seen that the surface cross-sections are not stellarator symmetric. The amplitudes of the $\cos (m\unicode[STIX]{x1D703}-n\unicode[STIX]{x1D711})$ and $\sin (m\unicode[STIX]{x1D703}-n\unicode[STIX]{x1D711})$ modes of $B$ inside this boundary, as computed by VMEC and BOOZ_XFORM, are shown in figure 18. The $\cos (\unicode[STIX]{x1D703}-5\unicode[STIX]{x1D711})$ term dominates, as desired, and its amplitude agrees with the prediction of the near-axis equations. Repeating the VMEC and BOOZ_XFORM computations for this solution of the near-axis equations for a range of aspect ratios, $[B_{m=0}(\unicode[STIX]{x1D711},r=a)-B(\unicode[STIX]{x1D711},r=0)]/a^{2}$ again converges to the predicted function $B_{20}(\unicode[STIX]{x1D711})$ , as shown in figure 19. Figure 20 shows that the symmetry-breaking modes scale as $1/A^{3}$ or better, as desired, except for the expected $1/A^{2}$ scaling of $S_{m=0}^{r=a}$ associated with the toroidal variation of $B_{20}$ . (The right-most blue points are missing since it did not seem possible to obtain values that were converged with respect to VMEC resolution parameters.) Thus, finite-aspect-ratio VMEC calculations successfully match the near-axis solution even when all parameters of the latter are non-zero.

## 6 Discussion and conclusions

In this work, we have developed a new and fast method to generate quasisymmetric magnetic fields with sophisticated shaping. In contrast to the traditional approach based on numerical optimization, the approach here uses a reduced set of equations relating the field strength in Boozer coordinates
$B(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$
to the three-dimensional shapes of the magnetic surfaces. The shapes that are describable by the
$O((r/{\mathcal{R}})^{2})$
near-axis model here are sufficiently general that they can be quite reminiscent of stellarators that have been designed previously using numerical optimization. For instance, the configuration of § 5.2 (figure 5) resembles CFQS (Liu *et al.*
Reference Liu, Shimizu, Isobe, Okamura, Nishimura, Suzuki, Xu, Zhang, Liu and Huang2018; Shimizu *et al.*
Reference Shimizu, Liu, Isobe, Okamura, Nishimura, Suzuki, Xu, Zhang, Liu and Huang2018) and the configuration of Henneberg *et al.* (Reference Henneberg, Drevlak, Nührenberg, Beidler, Turkin, Loizu and Helander2019). Also the configuration of § 5.4 (figure 13) resembles the HSX experiment (Anderson *et al.*
Reference Anderson, Almagri, Anderson, Matthews, Talmadge and Shohet1995). Despite these similarities, the examples here were generated independently of any previously known configurations. Since these shapes computed by our model are described analytically, they can be parameterized, evaluated rapidly and differentiated. As analytic expressions for the position vector in terms of both Boozer coordinates and cylindrical coordinates (appendix C) are available, one can evaluate virtually any quantity of interest, such as the geometric quantities appearing in the gyrokinetic model of turbulence. Since the system of equations involves only one independent variable (
$\unicode[STIX]{x1D711}$
), compared to three for general MHD equilibrium, the equations here are orders of magnitude faster to solve.

Through the examples in § 5, we have demonstrated that the approach here is a practical way to generate and parameterize both quasi-axisymmetric and quasi-helically symmetric configurations. For each of the examples, we showed that the departures from quasisymmetry computed by conventional codes scale with the aspect ratio as expected. In particular, we have demonstrated that quasisymmetry can be achieved (without axisymmetry) to any desired precision, at sufficiently high aspect ratio. Due to the high-order accuracy of the equations in our model, the quality of quasisymmetry can be extremely good. For example, the symmetry-breaking measures (5.2) are smaller than
$4\times 10^{-7}$
for the rightmost ‘Config 1’ point in figure 4. While at
$A=320$
this configuration is not of great experimental interest, it does represent the most accurate realization of quasisymmetry in a 3-D equilibrium ever reported. While arbitrarily small departures from quasisymmetry can also be obtained with the
$O((r/{\mathcal{R}})^{1})$
construction, as demonstrated in figure 4 of Landreman *et al.* (Reference Landreman, Sengupta and Plunk2019), higher aspect ratios would be required to obtain quasisymmetry to the same precision, due to the weaker scaling of symmetry breaking with
$1/A^{2}$
in that case. As an additional demonstration of the high accuracy to which quasisymmetry can be achieved by the
$O((r/{\mathcal{R}})^{2})$
construction, figure 21 shows the contours of
$B$
on the boundary surfaces of the configurations of §§ 5.2 and 5.4 with aspect ratios chosen by the following criterion: at a reactor-relevant mean field
$B_{0,0}=5$
T, the largest quasisymmetry-breaking Fourier mode amplitudes are only 0.5 Gauss, the approximate magnitude of Earth’s magnetic field. This condition results in aspect ratios
$A_{vmec}\sim$
80. There may be other configurations for which this condition can be met at lower aspect ratio; our goal here is merely to demonstrate that the condition can indeed be achieved. The deviation from symmetry is nearly invisible in figure 21, and the
$B$
contours are far more quasisymmetric than in other nominally quasisymmetric configurations reported previously, e.g. figures 5–7 of Beidler *et al.* (Reference Beidler, Allmaier, Isaev, Kasilov, Kernbichler, Leitold, Maassberg, Mikkelsen, Murakami and Schmidt2011).

Most of the solutions exhibited in this paper have a relatively high aspect ratio, which is not surprising since the method is based on an expansion in aspect ratio. Several methods are likely to enable configurations of lower aspect ratio to be generated. First, by including more Fourier modes in the axis shape,
$X_{2}$
and
$Y_{2}$
could perhaps be further reduced, resulting in configurations for which the
$O((r/{\mathcal{R}})^{3})$
symmetry-breaking terms have a smaller leading constant. Second, the method of § 4.4 of Landreman *et al.* (Reference Landreman, Sengupta and Plunk2019) could be used to extrapolate outward from a high-
$A$
configuration while preserving good quasisymmetry in the core. Third, a large value of
$a$
could be used in the present approach, resulting in moderate deviations from quasisymmetry, which could then be reduced by conventional optimization. Lastly, a configuration with smaller
$a$
and good quasisymmetry generated by the construction here could be used to initialize conventional optimization, in which the aspect ratio is included in the objective function for minimization.

The work here suggests many avenues for future study, some of which are enumerated here. (i) It was shown previously that practical quasisymmetric configurations obtained by optimization closely match the
$O(r/{\mathcal{R}})$
near-axis construction (Landreman Reference Landreman2019), and the comparison should be repeated for the
$O((r/{\mathcal{R}})^{2})$
model. (ii) An efficient procedure should be found to solve for model parameters such that
$B_{20}$
is independent of
$\unicode[STIX]{x1D711}$
. (iii) While a precise understanding exists of the solution space for
$O(r/{\mathcal{R}})$
quasisymmetry (Landreman *et al.*
Reference Landreman, Sengupta and Plunk2019), the same insight has yet to be developed for the
$O((r/{\mathcal{R}})^{2})$
model. It would be valuable to understand the space of solutions to the
$O((r/{\mathcal{R}})^{2})$
model to be sure all the interesting regions of parameter space have been identified. (iv) It was seen here that some toroidal variation of
$B_{20}$
could be allowed without
$B_{20}$
becoming the dominant quasisymmetry-breaking mode, so the effect of allowing small toroidal variation of
$B_{2c}$
or
$B_{2s}$
should be examined. (v) It should be investigated whether quasisymmetry could be optimized off-axis, by introducing small toroidal variation in
$B_{0}$
that is cancelled by
$B_{20}$
at a certain radius. (vi) The space of configurations that are omnigenous to
$O(r/{\mathcal{R}})$
was recently examined (Plunk *et al.*
Reference Plunk, Landreman and Helander2019), and the analysis could possibly be extended to
$O((r/{\mathcal{R}})^{2})$
omnigenity using results derived here.

Finally, extensions of the quasisymmetry model here to even higher order in
$r$
could be pursued. One motivation for such an extension is that global magnetic shear first appears at
$O((r/{\mathcal{R}})^{3})$
. Although quasisymmetry cannot generally be achieved through
$O((r/{\mathcal{R}})^{3})$
(Garren & Boozer Reference Garren and Boozer1991*a*
), the size of the
$O((r/{\mathcal{R}})^{3})$
terms informs how rapidly quasisymmetry degrades with
$r$
, so solutions could be sought in which these terms were minimized.

## Acknowledgements

This work was supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Science, under award nos DE-FG02-93ER54197 and DE-FG02-86ER53223. This work was also supported by a grant from the Simons Foundation (560651, M.L.).

## Appendix A. Derivation of the equations at each order

In this section we elaborate on § 2, showing a streamlined method to derive the required equations at each order in $r/{\mathcal{R}}$ . It is possible to obtain the same final equations without the manipulations of § A.1, but at the cost of substantial additional algebra. In particular, to derive the equations for $\{X_{j},Y_{j},Z_{j}\}$ at a given order $j$ without these manipulations, it would be necessary to first derive equations involving $Z_{j+1}$ , and then form linear combinations to eliminate this higher-order quantity. The method of § A.1 enables the equations for $\{X_{j},Y_{j},Z_{j}\}$ to be obtained directly without introducing $Z_{j+1}$ .

### A.1 Fundamental equations

To begin, note the product of (2.2) and (2.3) gives the inverse Jacobian

Then, from (2.2) and (2.3), applying (2.5)–(2.6) and using $\unicode[STIX]{x2202}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}=(r\bar{B})^{-1}\unicode[STIX]{x2202}/\unicode[STIX]{x2202}r$ , we obtain the following three scalar equations:

where $\{\ldots ,\ldots \}$ denotes a Poisson bracket in the $(r,\unicode[STIX]{x1D717})$ coordinates,

The right-hand sides of (A 2)–(A 4) are

and

where

and $\ell ^{\prime }=\sqrt{(\text{d}\boldsymbol{ r}_{0}/\text{d}\unicode[STIX]{x1D711})^{2}}$ .

As alluded to above, it turns out to be inconvenient to solve (A 2) and (A 3) as written, since these equations involve $Z$ to one higher order than $X$ or $Y$ . The reason is that $X_{1}$ and $Y_{1}$ are non-zero while $Z_{1}$ turns out to vanish (shown in the next subsection), so the left-hand sides at $O((r/{\mathcal{R}})^{j}r)$ include terms $\{rY_{1},r^{j+1}Z_{j+1}\}$ and $\{r^{j+1}Z_{j+1},rX_{1}\}$ , while the highest orders of $X$ and $Y$ appear through $\{r^{\,j}Y_{j},r^{2}Z_{2}\}$ and $\{r^{2}Z_{2},r^{\,j}X_{j}\}$ . Therefore it turns out to be convenient to form two combinations of (A 2) and (A 3), one in which $Z$ is given explicitly in terms of lower-order quantities, and the other in which the higher-order $Z$ terms are eliminated to give a constraint on the lower-order quantities. To form the first desired combination, we start by writing (A 2)–(A 3) as

This linear system can be solved to give

where the determinant can be recognized from (A 4) as $T_{Z}$ . The top row of (A 13) then gives an equation that will tell us $Z$ at each order in terms of lower-order quantities,

The equality of mixed partial derivatives $\unicode[STIX]{x2202}^{2}Z/\unicode[STIX]{x2202}r\unicode[STIX]{x2202}\unicode[STIX]{x1D717}=\unicode[STIX]{x2202}^{2}Z/\unicode[STIX]{x2202}\unicode[STIX]{x1D717}\unicode[STIX]{x2202}r$ can be used with (A 13) to obtain

This latter equation is the second desired combination of (A 2) and (A 3), giving a constraint on $X$ and $Y$ at each order without introducing $Z$ at the next order.

We can also derive a different combination of (A 2) and (A 3) with the same property, as an equivalent alternative to (A 15), which more closely corresponds to the equations of Garren and Boozer at $O((r/{\mathcal{R}})^{2})$ . This second approach begins with the observation that the problematic terms that introduce $Z$ at higher order than $X$ and $Y$ are $\{rY_{1},r^{j+1}Z_{j+1}\}$ and $\{r^{j+1}Z_{j+1},rX_{1}\}$ . To separate out these terms, we introduce $X_{{>}1}=X-rX_{1}$ and $Y_{{>}1}=Y-rY_{1}$ , so (A 2)–(A 3) give

*a*,

*b*) $$\begin{eqnarray}T_{X}-\{rY_{1},Z\}-\{Y_{{>}1},Z\}=0,\quad T_{Y}-\{Z,rX_{1}\}-\{Z,X_{{>}1}\}=0.\end{eqnarray}$$

We look for a combination of these equations in which the problematic terms $\{rY_{1},Z\}$ and $\{Z,rX_{1}\}$ are annihilated. To this end, it can be verified that

Forming the analogous combination of (A 16) then gives the desired relation, in which $Z$ appears at no higher order than $X$ or $Y$ ,

Finally, we obtain an expression for the magnetic field strength by squaring (2.2), and using (A 1),

Equations (A 4), (A 14), (A 18) and (A 19) are the four equations we will solve at each order for the corresponding unknowns $X$ , $Y$ , $Z$ and $B$ .

### A.2 Equations through $O((r/{\mathcal{R}})^{2})$

We now evaluate the first few orders of the $r/{\mathcal{R}}$ expansion, without assuming quasisymmetry. At $O((r/{\mathcal{R}})^{0})$ , (A 19) gives

where $s_{G}=\pm 1=\text{sign}(G_{0})$ , and (A 14) gives $Z_{1}=0$ . Equations (A 4) and (A 18) have no terms of this order. At $O((r/{\mathcal{R}})^{1})$ , (A 4) gives

the $\sin \unicode[STIX]{x1D717}$ and $\cos \unicode[STIX]{x1D717}$ modes of (A 19) give

*a*,

*b*) $$\begin{eqnarray}B_{1s}=\unicode[STIX]{x1D705}X_{1s}B_{0},\quad B_{1c}=\unicode[STIX]{x1D705}X_{1c}B_{0},\end{eqnarray}$$

and (A 18) gives

where primes denote $\text{d}/\text{d}\unicode[STIX]{x1D711}$ and

It is convenient to introduce $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D711})=(B_{1s}Y_{1s}+B_{1c}Y_{1c})/(s_{G}\bar{B}\unicode[STIX]{x1D705})$ , in which case (A 21)–(A 22) imply

*a*,

*b*) $$\begin{eqnarray}Y_{1s}=\frac{(B_{1c}+B_{1s}\unicode[STIX]{x1D70E})s_{G}\bar{B}\unicode[STIX]{x1D705}}{B_{1s}^{2}+B_{1c}^{2}},\quad Y_{1c}=\frac{(-B_{1s}+B_{1c}\unicode[STIX]{x1D70E})s_{G}\bar{B}\unicode[STIX]{x1D705}}{B_{1s}^{2}+B_{1c}^{2}},\end{eqnarray}$$

and (A 23)–(A 24) can be written

Next, the $\unicode[STIX]{x1D717}$ -independent, $\sin 2\unicode[STIX]{x1D717}$ and $\cos 2\unicode[STIX]{x1D717}$ modes of (A 14) give

where

At $O((r/{\mathcal{R}})^{2})$ , the $\sin \unicode[STIX]{x1D717}$ and $\cos \unicode[STIX]{x1D717}$ terms of (A 4) are

The $\unicode[STIX]{x1D717}$ -independent, $\sin 2\unicode[STIX]{x1D717}$ and $\cos 2\unicode[STIX]{x1D717}$ modes of (A 19) at $O((r/{\mathcal{R}})^{2})$ give

where

The $\cos \unicode[STIX]{x1D717}$ and $\sin \unicode[STIX]{x1D717}$ terms of (A 18) at $O((r/{\mathcal{R}})^{2})$ are

and

where

We will not need the $O((r/{\mathcal{R}})^{2})$ terms of (A 14), which give $Z_{3}$ .

Finally, for the analysis in § 3, we need the independent-of- $\unicode[STIX]{x1D717}$ mode of (A 4) at $O((r/{\mathcal{R}})^{3})$ , which gives

with $Q$ given by (