Direct construction of optimized stellarator shapes. Part 2. Numerical quasisymmetric solutions

Abstract

Quasisymmetric stellarators are appealing intellectually and as fusion reactor candidates since the guiding-centre particle trajectories and neoclassical transport are isomorphic to those in a tokamak, implying good confinement. Previously, quasisymmetric magnetic fields have been identified by applying black-box optimization algorithms to minimize symmetry-breaking Fourier modes of the field strength $B$. Here, instead, we directly construct magnetic fields in cylindrical coordinates that are quasisymmetric to leading order in the distance from the magnetic axis, without using optimization. The method involves solution of a one-dimensional nonlinear ordinary differential equation, originally derived by Garren & Boozer (Phys. Fluids B, vol. 3, 1991, p. 2805). We demonstrate the usefulness and accuracy of this optimization-free approach by providing the results of this construction as input to the codes VMEC and BOOZ_XFORM, confirming the purity and scaling of the magnetic spectrum. The space of magnetic fields that are quasisymmetric to this order is parameterized by the magnetic axis shape along with three other real numbers, one of which reflects the on-axis toroidal current density, and another one of which is zero for stellarator symmetry. The method here could be used to generate good initial conditions for conventional optimization, and its speed enables exhaustive searches of parameter space.

1 Introduction

Toroidal magnetic fields can possess a remarkable hidden symmetry, called quasisymmetry, in which the field strength $B=|\boldsymbol{B}|$ is independent of a particular coordinate (‘Boozer angle’) even though the magnetic field vector $\boldsymbol{B}$ is not (Boozer 1983; Nührenberg & Zille 1988; Helander 2014). Since the Lagrangian for guiding-centre particle motion in Boozer coordinates varies on magnetic surfaces only through $B$ , a symmetry direction in $B$ implies that guiding-centre trajectories behave as if the magnetic field had a true symmetry direction, and the conserved quantity that follows from Noether’s theorem implies that particle trajectories are confined. In contrast, magnetic fields without continuous symmetry generally have unconfined guiding-centre trajectories. (Quasisymmetry is sufficient but not necessary for guiding-centre confinement (Cary & Shasharina 1997).) Plasmas confined by quasisymmetric magnetic fields are also predicted to have temperature screening of impurities and to allow larger flows, which may lead to improved stability. For these reasons, quasisymmetric magnetic fields are interesting both for fusion energy and on basic physics grounds.

A number of quasisymmetric magnetic configurations have been identified to date (Nührenberg & Zille 1988; Nührenberg, Lotz & Gori 1994; Anderson et al. 1995; Garabedian 1996; Zarnstorff et al. 2001; Ku & Boozer 2011; Drevlak et al. 2013, 2018; Henneberg et al. 2018; Plunk & Helander 2018). In all of these cases, except the work of Plunk & Helander (2018), the quasisymmetric configurations have been found using optimization, by minimizing the amplitudes of symmetry-breaking Fourier modes of $B$ . The optimization algorithms used have been ‘off the shelf’ algorithms that can be applied to minimizing any function and do not exploit information about the underlying physical system. While this approach has proven successful, it does have a number of shortcomings. Little insight is provided as to the form and dimensionality of the landscape of solutions. As the results of the optimization depend on the initial guess and on manually chosen weight parameters, there is no guarantee that all interesting solutions have been found. Optimization is also computationally demanding, requiring many three-dimensional (3-D) equilibrium calculations.

A complementary approach to finding quasisymmetric geometries, developed by Garren & Boozer (1991a ), is to directly construct the geometry from the relevant equations, with no need then for optimization. Expanding in small distance $r$ from the magnetic axis (that is, large aspect ratio), Garren & Boozer derived equations for quasisymmetry to first and second order in $r$ . Their work is perhaps best known for the result that the number of equations exceeds the number of unknowns at third order, so quasisymmetry may be achieved on one surface but not throughout a volume. However the useful constructive procedure at lower order has not been fully exploited as a tool to generate quasisymmetric shapes, which can be useful as initial conditions for conventional optimization, and to understand the landscape of quasisymmetric shapes. The goal of this paper is to reinvigorate this development.

In an accompanying Part 1 (Landreman & Sengupta 2018), we derived two ways to generate a shape in standard cylindrical coordinates with prescribed $B$ using the Garren–Boozer framework, summarized in §§ 4.1 and 4.3 below. In the present paper, we develop the optimization-free approach to constructing quasisymmetric geometries in several ways. In § 3, we present a new spectrally accurate algorithm for solving the equation for quasisymmetry to first order in $r$ . Using several methods for converting the results to standard cylindrical coordinates, explained in § 4, we present in § 5 examples of quasi-axisymmetric and quasi-helically symmetric equilibria obtained without optimization. (Quasipoloidally symmetric configurations cannot be generated using this approach since this symmetry is impossible near the axis.) For each of these configurations, we use the codes VMEC (Hirshman & Whitson 1983; Hirshman, van Rij & Merkel 1986) and BOOZ_XFORM (Sanchez et al. 2000) to compute the spectrum of $B$ , confirming that the symmetry-breaking harmonics are small and that they scale as expected. One family of equilibria we consider (§ 5.3) possesses quasisymmetry but not stellarator symmetry, which may be desirable for obtaining significant intrinsic rotation. We discuss and conclude in § 6. A proof that a unique solution to the first-order quasisymmetry equation exists despite the nonlinearity of the problem is given in appendix A.

The approach here allows quasisymmetric flux surface shapes to be computed in ${<}1~\text{ms}$ on a laptop. This time scale is at least 4 orders of magnitude faster than a typical equilibrium calculation with VMEC, much less an optimization in which VMEC is iterated to find quasisymmetric equilibria. Our approach can therefore be used for extensive searches of parameter space, potentially enabling an identification of all possible quasisymmetric plasma shapes, at least in the vicinity of the magnetic axis. Also, this ‘direct construction’ approach makes clear how many degrees of freedom are available in the space of quasisymmetric magnetic fields, giving insight into the landscape of solutions.

2 System of equations

Here we summarize the equations relevant to first-order quasisymmetry derived in Garren & Boozer (1991a ). The position vector can be written in flux coordinates as

(2.1) $$\begin{eqnarray}\boldsymbol{r}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})=\boldsymbol{r}_{0}(\unicode[STIX]{x1D711})+rX_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})\boldsymbol{n}(\unicode[STIX]{x1D711})+rY_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})\boldsymbol{b}(\unicode[STIX]{x1D711})+O(r^{2}),\end{eqnarray}$$

where $\boldsymbol{r}_{0}$ is the position of the magnetic axis, $(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$ are the poloidal and toroidal Boozer angles, $r=\sqrt{2|\unicode[STIX]{x1D713}|/B_{0}}$ is an effective minor radius that labels flux surfaces, $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D713}$ is the toroidal flux and $B_{0}$ is the magnetic field strength along the axis, which must be a constant due to quasisymmetry. The unit normal vector $\boldsymbol{n}$ and unit binormal $\boldsymbol{b}$ are defined in terms of the magnetic axis shape by the Frenet–Serret relations

(2.2) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\text{d}\boldsymbol{t}/\text{d}\ell =\unicode[STIX]{x1D705}\boldsymbol{n},\\ \text{d}\boldsymbol{n}/\text{d}\ell =-\unicode[STIX]{x1D705}\boldsymbol{t}+\unicode[STIX]{x1D70F}\boldsymbol{b},\\ \text{d}\boldsymbol{b}/\text{d}\ell =-\unicode[STIX]{x1D70F}\boldsymbol{n},\end{array}\right\}\end{eqnarray}$$

where $\boldsymbol{t}(\unicode[STIX]{x1D711})=\text{d}\boldsymbol{r}_{0}/\text{d}\ell$ is the unit tangent vector, $\boldsymbol{t}\boldsymbol{\cdot }\boldsymbol{n}\times \boldsymbol{b}=1$ , $\ell$ denotes arclength, $\unicode[STIX]{x1D705}(\unicode[STIX]{x1D711})$ is the curvature and $\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D711})$ is the torsion. (Garren and Boozer use the opposite sign convention for torsion.) The $\boldsymbol{t}$ component in (2.1) at $O(r)$ can be shown to vanish. To first order in the distance from the magnetic axis, the flux surface shape is described by

(2.3a,b ) $$\begin{eqnarray}X_{1}=X_{1s}(\unicode[STIX]{x1D711})\sin \unicode[STIX]{x1D703}+X_{1c}(\unicode[STIX]{x1D711})\cos \unicode[STIX]{x1D703},\quad Y_{1}=Y_{1s}(\unicode[STIX]{x1D711})\sin \unicode[STIX]{x1D703}+Y_{1c}(\unicode[STIX]{x1D711})\cos \unicode[STIX]{x1D703},\end{eqnarray}$$

and the magnetic field strength satisfies

(2.4) $$\begin{eqnarray}B(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})=B_{0}+r[B_{1s}(\unicode[STIX]{x1D711})\sin \unicode[STIX]{x1D703}+B_{1c}(\unicode[STIX]{x1D711})\cos \unicode[STIX]{x1D703}]+O(r^{2}),\end{eqnarray}$$

where

(2.5a,b ) $$\begin{eqnarray}X_{1s}(\unicode[STIX]{x1D711})=\frac{B_{1s}(\unicode[STIX]{x1D711})}{B_{0}\unicode[STIX]{x1D705}(\unicode[STIX]{x1D711})},\quad X_{1c}(\unicode[STIX]{x1D711})=\frac{B_{1c}(\unicode[STIX]{x1D711})}{B_{0}\unicode[STIX]{x1D705}(\unicode[STIX]{x1D711})}.\end{eqnarray}$$

In the case of quasisymmetry, we can choose the origin of the $\unicode[STIX]{x1D703}$ coordinate so $B=B_{0}+r\bar{\unicode[STIX]{x1D702}}B_{0}\cos (\unicode[STIX]{x1D703}-N\unicode[STIX]{x1D711})+O(r^{2})$ for some constant $\bar{\unicode[STIX]{x1D702}}$ and fixed integer $N$ , with quasi-axisymmetry defined by $N=0$ and quasi-helical symmetry defined by $N\neq 0$ . Then $B_{1s}=\bar{\unicode[STIX]{x1D702}}B_{0}\sin (N\unicode[STIX]{x1D711})$ and $B_{1c}=\bar{\unicode[STIX]{x1D702}}B_{0}\cos (N\unicode[STIX]{x1D711})$ . Furthermore,

(2.6a,b ) $$\begin{eqnarray}Y_{1s}=\left(\frac{s_{G}s_{\unicode[STIX]{x1D713}}\unicode[STIX]{x1D705}}{\bar{\unicode[STIX]{x1D702}}}\right)[\unicode[STIX]{x1D70E}\sin (N\unicode[STIX]{x1D711})+\cos (N\unicode[STIX]{x1D711})],\quad Y_{1c}=\left(\frac{s_{G}s_{\unicode[STIX]{x1D713}}\unicode[STIX]{x1D705}}{\bar{\unicode[STIX]{x1D702}}}\right)[\unicode[STIX]{x1D70E}\cos (N\unicode[STIX]{x1D711})-\sin (N\unicode[STIX]{x1D711})],\end{eqnarray}$$

where $s_{G}=\pm 1$ is positive (negative) if $\boldsymbol{B}$ points towards increasing (decreasing) $\unicode[STIX]{x1D711}$ , $s_{\unicode[STIX]{x1D713}}=\text{sign}(\unicode[STIX]{x1D713})=\pm 1$ and the periodic function $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D711})$ satisfies the Riccati-type equation

(2.7) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D70E}}{\text{d}\unicode[STIX]{x1D711}}+(\unicode[STIX]{x1D704}-N)\left[\frac{\bar{\unicode[STIX]{x1D702}}^{4}}{\unicode[STIX]{x1D705}^{4}}+1+\unicode[STIX]{x1D70E}^{2}\right]-\frac{2G_{0}\bar{\unicode[STIX]{x1D702}}^{2}}{B_{0}\unicode[STIX]{x1D705}^{2}}\left[\frac{I_{2}}{B_{0}}-s_{\unicode[STIX]{x1D713}}\unicode[STIX]{x1D70F}\right]=0.\end{eqnarray}$$

Here, $\unicode[STIX]{x1D704}$ is the rotational transform on axis; $G_{0}$ is the on-axis value of $G(r)$ , the poloidal current outside the flux surface times $\unicode[STIX]{x1D707}_{0}/(2\unicode[STIX]{x03C0})$ ; and $I_{2}$ is the leading coefficient in $I(r)=r^{2}I_{2}+O(r^{4})$ , the toroidal current inside the flux surface times $\unicode[STIX]{x1D707}_{0}/(2\unicode[STIX]{x03C0})$ . The functions $G(r)$ and $I(r)$ here are those appearing in the Boozer coordinate representation

(2.8) $$\begin{eqnarray}\boldsymbol{B}=\unicode[STIX]{x1D6FD}(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})\unicode[STIX]{x1D735}r+I(r)\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+G(r)\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}.\end{eqnarray}$$

Using equation (2.20) in Part 1, equation (2.7) can be written in terms of the standard toroidal angle $\unicode[STIX]{x1D719}$ (the azimuthal angle in cylindrical coordinates $(R,\unicode[STIX]{x1D719},z)$ ) as

(2.9) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D70E}}{\text{d}\unicode[STIX]{x1D711}}=\frac{|G_{0}|}{\ell ^{\prime }B_{0}}\frac{\text{d}\unicode[STIX]{x1D70E}}{\text{d}\unicode[STIX]{x1D719}},\end{eqnarray}$$

where

(2.10) $$\begin{eqnarray}\ell ^{\prime }=\frac{\text{d}\ell }{\text{d}\unicode[STIX]{x1D719}}=\sqrt{R_{0}^{2}+(R_{0}^{\prime })^{2}+(z_{0}^{\prime })^{2}},\end{eqnarray}$$

the magnetic axis has cylindrical coordinates $R_{0}(\unicode[STIX]{x1D719})$ and $z_{0}(\unicode[STIX]{x1D719})$ , and primes denote $\text{d}/\text{d}\unicode[STIX]{x1D719}$ . Also, $G_{0}$ can be related to the magnetic axis shape by $G_{0}=s_{G}B_{0}L/(2\unicode[STIX]{x03C0})$ where $L=\int _{0}^{2\unicode[STIX]{x03C0}}\,\text{d}\unicode[STIX]{x1D719}\;\ell ^{\prime }$ is the length of the axis. The above equations all apply even if the plasma pressure is non-zero, although the pressure turns out not to appear in these expressions to this order.

As discussed in § 5.2 of Part 1, $N$ can be determined directly from the axis shape. The integer $N$ is the number of times the normal vector $\boldsymbol{n}$ rotates poloidally around the magnetic axis as the axis is traversed once toroidally. We can determine $N$ this way because the vector at each $\unicode[STIX]{x1D711}$ pointing from the magnetic axis to the maximum- $B$ contour on a flux surface $r$ is $\boldsymbol{n}r|\bar{\unicode[STIX]{x1D702}}|/\unicode[STIX]{x1D705}+\boldsymbol{b}rY_{1}$ , which has a positive projection onto $\boldsymbol{n}$ at every $\unicode[STIX]{x1D711}$ . Hence these two vectors are always within $90^{\circ }$ of each other, and so the $B$ contours loop around the magnetic axis the same number of times $\boldsymbol{n}$ does so.

Stellarators, whether quasisymmetric or not, typically are designed to possess stellarator symmetry, since this symmetry reduces the dimensionality of the parameter space for optimization, reduces the computational cost of many calculations, and typically reduces the number of unique coil shapes. Stellarator symmetry corresponds to $R(-\unicode[STIX]{x1D703},-\unicode[STIX]{x1D719})=R(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ , $z(-\unicode[STIX]{x1D703},-\unicode[STIX]{x1D719})=-z(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ and $B(-\unicode[STIX]{x1D703},-\unicode[STIX]{x1D719})=B(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ . For a magnetic field described by (2.1)–(2.7) to be stellarator symmetric, the magnetic axis shape must be stellarator symmetric, and $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D711})$ should be odd.

As proved in appendix A, even though (2.7) is nonlinear in $\unicode[STIX]{x1D70E}$ , this equation can be posed in such a way that there is guaranteed to be precisely one solution. Specifically, given well-behaved $\unicode[STIX]{x1D705}(\unicode[STIX]{x1D711})$ , $\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D711})$ , $I_{2}/B_{0}$ , $G_{0}/B_{0}$ , $\bar{\unicode[STIX]{x1D702}}$ and an initial condition $\unicode[STIX]{x1D70E}(0)$ , there is precisely one solution pair $\{\unicode[STIX]{x1D704},\unicode[STIX]{x1D70E}\}$ such that $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D711})$ is periodic. As a result, for any magnetic axis shape with non-vanishing curvature, there is an infinite number of magnetic fields in the vicinity of that axis which are consistent with quasisymmetry to first order in $r$ . The possible magnetic fields are parameterized by three numbers: $I_{2}/B_{0}$ , $\bar{\unicode[STIX]{x1D702}}$ and $\unicode[STIX]{x1D70E}(0)$ . If the current density vanishes on axis, which is common in stellarators even at finite plasma pressure since the bootstrap current density vanishes on axis, $I_{2}=0$ . Furthermore, for stellarator-symmetric fields, $\unicode[STIX]{x1D70E}(0)=0$ . Therefore, in practice usually only one of the three scalar input parameters is free. While every magnetic axis shape with non-vanishing curvature admits an infinite number of quasisymmetric fields, for many axis shapes it is found numerically that the elongation of the surrounding flux surfaces reaches enormous values (tens, hundreds or more), making the shape uninteresting.

3 Numerical method

For a practical solution of (2.7) we consider the inputs to be $\bar{\unicode[STIX]{x1D702}}$ , $I_{2}/B_{0}$ , $\unicode[STIX]{x1D70E}(0)$ and the shape of the magnetic axis $\{R_{0}(\unicode[STIX]{x1D719}),z_{0}(\unicode[STIX]{x1D719})\}$ . The outputs are $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D719})$ and $\unicode[STIX]{x1D704}$ . Given the inputs, we solve (2.7) with (2.9) for $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D719})$ using Newton iteration with a pseudo-spectral collocation discretization. A uniform grid of $N_{\unicode[STIX]{x1D719}}$ points, $\unicode[STIX]{x1D719}_{j}=(j-1)2\unicode[STIX]{x03C0}/(N_{\unicode[STIX]{x1D719}}n_{fp})$ where $j=1\ldots N_{\unicode[STIX]{x1D719}}$ , is defined on $[0,2\unicode[STIX]{x03C0}/n_{fp})$ , where $n_{fp}$ is the number of identical field periods. The vector of unknowns is taken to be $[\unicode[STIX]{x1D704},\unicode[STIX]{x1D70E}_{2},\ldots ,\unicode[STIX]{x1D70E}_{N_{\unicode[STIX]{x1D719}}}]^{\text{T}}$ where $\unicode[STIX]{x1D70E}_{j}=\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D719}_{j})$ , so there are $N_{\unicode[STIX]{x1D719}}$ unknowns. There is no need to include $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D719}_{1})=\unicode[STIX]{x1D70E}(0)$ as an unknown since it is a prescribed input. A system of $N_{\unicode[STIX]{x1D719}}$ equations is obtained by imposing (2.7) at all $\unicode[STIX]{x1D719}_{j}$ . The $\text{d}\unicode[STIX]{x1D70E}/\text{d}\unicode[STIX]{x1D719}$ derivative is discretized using the Fourier pseudo-spectral differentiation matrix (Weideman & Reddy 2000). Newton iteration proceeds by solving linear systems involving the Jacobian matrix $[\unicode[STIX]{x2202}\boldsymbol{R}/\unicode[STIX]{x2202}\unicode[STIX]{x1D704},\unicode[STIX]{x2202}\boldsymbol{R}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{2},\ldots ,\unicode[STIX]{x2202}\boldsymbol{R}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70E}_{N_{\unicode[STIX]{x1D719}}}]$ , where $\boldsymbol{R}$ is the residual vector. It is straightforward to analytically evaluate the derivatives in the Jacobian in terms of the differentiation matrix. For the examples shown below, the residual $L^{2}$ norm is reduced by 15 orders of magnitude in ${\leqslant}5$ Newton iterations. The numerical solution is extremely robust: in parameter scans to date we have not observed any examples in which the Newton iteration fails to converge.

Figure 1 demonstrates the convergence of the rotational transform computed by this method as the number of grid points increases, for the example of § 5.1. As expected, the convergence is spectral, and $\unicode[STIX]{x1D704}$ can be computed to 15 digits of precision with $N_{\unicode[STIX]{x1D719}}\sim 50$ . In the figure, the ‘true’ value of $\unicode[STIX]{x1D704}$ is taken to be the result for $N_{\unicode[STIX]{x1D719}}=149$ .

Figure 1. The algorithm of § 3 allows the equations of § 2 or § 4.3 to be solved to machine precision with a modest number of grid points $N_{\unicode[STIX]{x1D719}}$ . The equations of these two sections yield results that are identical to machine precision (for sufficient $N_{\unicode[STIX]{x1D719}}$ ) since the equations are equivalent, as proved in Part 1.

4 Conversion to cylindrical coordinates

To take advantage of stellarator physics codes that accept VMEC equilibrium files, such as the STELLOPT optimization suite, we wish to transform the solutions from the Frenet–Serret frame to the VMEC input representation. In this representation, the plasma boundary is expressed as a Fourier expansion of the cylindrical coordinates $R(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ and $z(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ , where $\unicode[STIX]{x1D703}$ can be any poloidal angle. Here, we will continue to let $\unicode[STIX]{x1D703}$ be the poloidal Boozer angle. We can compute $R(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ and $z(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ from the asymptotic large-aspect-ratio solution in several ways, described in the following subsections. The first three approaches have the common feature that an expansion in $r$ for the surface shape is evaluated at a finite value of  $r$ .

4.1 First-order method

In one approach, the solution in the Frenet–Serret frame is transformed to cylindrical coordinates using the method detailed in Part 1, summarized here. The position vector is expressed as

(4.1) $$\begin{eqnarray}\boldsymbol{r}=\hat{\boldsymbol{r}}_{0}(\unicode[STIX]{x1D719})+r[R_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})\boldsymbol{e}_{R}(\unicode[STIX]{x1D719})+z_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})\boldsymbol{e}_{z}]+O(r^{2}),\end{eqnarray}$$

where $\hat{\boldsymbol{r}}_{0}(\unicode[STIX]{x1D719})=R_{0}(\unicode[STIX]{x1D719})\boldsymbol{e}_{R}(\unicode[STIX]{x1D719})+z_{0}(\unicode[STIX]{x1D719})\boldsymbol{e}_{z}$ , and equated to (2.1). The $\boldsymbol{n}$ and $\boldsymbol{b}$ components of the result give (to leading order in $r$ )

(4.2) $$\begin{eqnarray}\left(\begin{array}{@{}c@{}}X_{1}\\ Y_{1}\end{array}\right)=\left(\begin{array}{@{}cc@{}}n_{R} & n_{z}\\ b_{R} & b_{z}\end{array}\right)\left(\begin{array}{@{}c@{}}R_{1}\\ z_{1}\end{array}\right),\end{eqnarray}$$

where $n_{R}=\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_{R}$ , $b_{z}=\boldsymbol{b}\boldsymbol{\cdot }\boldsymbol{e}_{z}$ , etc., and this matrix equation can be inverted to give

(4.3) $$\begin{eqnarray}\left(\begin{array}{@{}c@{}}R_{1}\\ z_{1}\end{array}\right)=\frac{\ell ^{\prime }}{R_{0}}\left(\begin{array}{@{}cc@{}}-b_{z} & n_{z}\\ b_{R} & -n_{R}\end{array}\right)\left(\begin{array}{@{}c@{}}X_{1}\\ Y_{1}\end{array}\right).\end{eqnarray}$$

Since $X_{1}$ and $Y_{1}$ each have $\sin \unicode[STIX]{x1D703}$ and $\cos \unicode[STIX]{x1D703}$ components, the same is true of $R_{1}$ and $z_{1}$ , so the flux surfaces are ellipses in the $R$ – $z$ plane. Equation (4.3) can be applied at each $\unicode[STIX]{x1D719}$ to both the $\sin \unicode[STIX]{x1D703}$ and $\cos \unicode[STIX]{x1D703}$ components. Then given any choice for $r$ , a finite-aspect-ratio magnetic surface can be formed in cylindrical coordinates from $R(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=R_{0}(\unicode[STIX]{x1D719})+rR_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ and $z(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})=z_{0}(\unicode[STIX]{x1D719})+rz_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ . Note that the transformation (4.3) represents only the leading-order behaviour in an expansion in $r$ , so for finite $r$ the flux surface geometry will depart somewhat from (2.1), meaning cross-sections of the boundary surface normal to the magnetic axis will no longer be perfect ellipses.

4.2 Alternative method

The method of the previous subsection results in a flux surface shape that is consistent with (2.1) to $O(r)$ . Alternatively, one can compute the surface defined by the terms through $O(r)$ in (2.1) as follows. First a positive value of $r$ is chosen. Given uniform grids in $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D719}$ , a tensor product grid is formed. For each point $(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ in this tensor product grid, a 1-D root finding problem is solved to find $\unicode[STIX]{x1D711}$ such that the position vector (2.1) has toroidal angle $\unicode[STIX]{x1D719}$ . In this way, $R$ and $z$ are obtained on the $(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ grid, and so they can be Fourier transformed to provide input to VMEC. To $O(r)$ the resulting surface is identical to the surface constructed in the previous subsection, but $O(r^{2})$ differences are present. This second method ensures cross-sections of the flux surfaces perpendicular to the magnetic axis are elliptical, while cross-sections in the $R$ – $z$ plane will generally not be elliptical.

For the quasi-axisymmetric examples below, we find the two methods for converting to cylindrical coordinates yield nearly indistinguishable results, and so there is no need for the extra complexity of the second method. However, for the quasi-helically symmetric example below, we find the second method yields smaller symmetry-breaking harmonics by a factor ${\sim}2$ , and so we will use it in § 5.2.

4.3 Direct solution in cylindrical coordinates

There is also a third approach to computing $R(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ and $z(\unicode[STIX]{x1D703},\unicode[STIX]{x1D719})$ for first-order quasisymmetric magnetic surface shapes: directly solving the first-order quasisymmetry equations in cylindrical coordinates rather than in the Frenet–Serret frame. The first-order quasisymmetry equations in cylindrical coordinates were derived in Part 1 and are

(4.4) $$\begin{eqnarray}\displaystyle & \displaystyle R_{1s}z_{1c}-R_{1c}z_{1s}-\frac{s_{G}\ell ^{\prime }}{R_{0}}=0, & \displaystyle\end{eqnarray}$$

(4.5) $$\begin{eqnarray}\displaystyle & \displaystyle K_{R}R_{1c}+K_{z}z_{1c}-\frac{B_{1c}}{B_{0}}=0, & \displaystyle\end{eqnarray}$$

(4.6) $$\begin{eqnarray}\displaystyle & \displaystyle K_{R}R_{1s}+K_{z}z_{1s}-\frac{B_{1s}}{B_{0}}=0, & \displaystyle\end{eqnarray}$$

(4.7) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D704}V-T=0, & \displaystyle\end{eqnarray}$$

where

(4.8a,b ) $$\begin{eqnarray}K_{R}=\unicode[STIX]{x1D705}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_{R},\quad K_{R}=\unicode[STIX]{x1D705}\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{e}_{z}\end{eqnarray}$$

(4.9) $$\begin{eqnarray}\displaystyle T & = & \displaystyle \frac{|G_{0}|}{(\ell ^{\prime })^{3}B_{0}} [R_{0}^{2}(R_{1c}R_{1s}^{\prime }-R_{1s}R_{1c}^{\prime }+z_{1c}z_{1s}^{\prime }-z_{1s}z_{1c}^{\prime })\nonumber\\ \displaystyle & & \displaystyle +\,(R_{1c}z_{1s}-R_{1s}z_{1c})(R_{0}^{\prime }z_{0}^{\prime \prime }+2R_{0}z_{0}^{\prime }-z_{0}^{\prime }R_{0}^{\prime \prime })\nonumber\\ \displaystyle & & \displaystyle +\,(z_{1c}z_{1s}^{\prime }-z_{1s}z_{1c}^{\prime })(R_{0}^{\prime })^{2}+(R_{1c}R_{1s}^{\prime }-R_{1s}R_{1c}^{\prime })(z_{0}^{\prime })^{2}\nonumber\\ \displaystyle & & \displaystyle +\,(R_{1s}z_{1c}^{\prime }-z_{1c}R_{1s}^{\prime }+z_{1s}R_{1c}^{\prime }-R_{1c}z_{1s}^{\prime })R_{0}^{\prime }z_{0}^{\prime }]+\frac{2G_{0}I_{2}}{B_{0}^{2}}\end{eqnarray}$$

and

(4.10) $$\begin{eqnarray}\displaystyle V & = & \displaystyle \frac{1}{(\ell ^{\prime })^{2}} [R_{0}^{2}(R_{1c}^{2}+R_{1s}^{2}+z_{1c}^{2}+z_{1s}^{2})+(R_{0}^{\prime })^{2}(z_{1c}^{2}+z_{1s}^{2})\nonumber\\ \displaystyle & & \displaystyle -\,2R_{0}^{\prime }z_{0}^{\prime }(R_{1c}z_{1c}+R_{1s}z_{1s})+(z_{0}^{\prime })^{2}(R_{1c}^{2}+R_{1s}^{2})],\end{eqnarray}$$

and primes denote $\text{d}/\text{d}\unicode[STIX]{x1D719}$ . As proved in Part 1, these equations are exactly equivalent to (2.7) under the transformation (4.2)–(4.3).

The system (4.4)–(4.7) can be solved with Newton’s method using a procedure similar to the one of § 3. The vector of unknowns consists of $R_{1c}$ , $R_{1s}$ , $z_{1c}$ and $z_{1s}$ , each evaluated at $\unicode[STIX]{x1D719}_{j}$ , along with $\unicode[STIX]{x1D704}$ . The same number of equations is obtained by imposing (4.4)–(4.7) at each of the $\unicode[STIX]{x1D719}_{j}$ , along with one additional equation corresponding to the initial condition for $\unicode[STIX]{x1D70E}$ . We verified that this direct solution in cylindrical coordinates (4.4)–(4.7) indeed yields identical results to the method of § 4.1, within discretization error that can be made as small as machine precision, as shown in figure 1.

4.4 Outward extrapolation using specific coils

A fourth method for generating finite-size quasisymmetric plasma shapes from the high-aspect-ratio theory, which we now describe, can potentially generate shapes with relatively low aspect ratio that are realizable with reasonable coils, at least for the limited case of vacuum fields. In this method, first one of the methods of §§ 4.1–4.3 is used to generate a flux surface shape at a high aspect ratio. Next, a coil design code such as REGCOIL (Landreman 2017) or FOCUS (Zhu et al. 2018) is used to find coil shapes that produce this high-aspect-ratio surface, by minimizing the (squared) magnetic field normal to the surface. Typically, good flux surfaces will be produced by these coils well outside of the original target surface, and field line tracing can be used to identify a large region filled with good surfaces. If desired, VMEC can be run in free-boundary mode to obtain a representation of the field in this larger region. This fourth approach results in boundary surface shapes that are not strictly ellipses. There is substantial flexibility in this method, as the designer can choose the number of coils, the regularity of the coil shapes and any other input parameters to the coil design code. There is no particular reason the magnetic field will be quasisymmetric outside of the smaller high-aspect-ratio target volume, so this procedure tends to produce better quasisymmetry on axis than at the edge. Since this method requires a coil design code, which takes at least ${\sim}10~\text{s}$ to run in the case of REGCOIL, as well as field line tracing to find the resulting surfaces, the computational cost is higher than that of the previous methods, although still very small compared to conventional stellarator optimization. An example of this method will be shown at the end of § 5.1.

5 Examples

5.1 Quasi-axisymmetry

We now demonstrate the procedures of the previous sections to construct a variety of quasisymmetric stellarator shapes. We begin with an example of quasi-axisymmetry, considering the magnetic axis shape

(5.1a,b ) $$\begin{eqnarray}R_{0}(\unicode[STIX]{x1D719})=1+0.045\cos (3\unicode[STIX]{x1D719}),\quad z_{0}(\unicode[STIX]{x1D719})=-0.045\sin (3\unicode[STIX]{x1D719}),\end{eqnarray}$$

with $\bar{\unicode[STIX]{x1D702}}=-0.9$ . We take $\unicode[STIX]{x1D70E}(0)=0$ (stellarator symmetry). For this and the later examples, we consider $I_{2}=0$ . For these parameters, the numerical procedure above yields a rotational transform $\unicode[STIX]{x1D704}=0.418$ , and the maximum flux surface elongation in the $R$ – $z$ plane is found to be 2.40. Hereafter we call $1/r$ the aspect ratio, since the average major radius is 1. The flux surfaces for aspect ratio 10 are shown in figure 2. Supplying this surface as an input to VMEC, the resulting magnetic field strength on the boundary is shown in figure 2(a). The Fourier spectrum of $B$ in Boozer coordinates at each flux surface is then computed using the BOOZ_XFORM code (Sanchez et al. 2000). The resulting spectra for aspect ratios 10 and 80 are shown in figure 3. At aspect ratio 10, the $(m,n)=(1,0)$ harmonic is dominant across all surfaces, as desired, and the quality of the quasisymmetry increases as the aspect ratio is increased. For both aspect ratios shown, the largest symmetry-breaking mode at the edge is the mode $(m,n)=(2,-3)$ . Since modes of $B$ with given poloidal mode number $m$ have amplitude $\propto r^{m}$ near the axis (Garren & Boozer 1991a ), the symmetry breaking on axis is dominated by modes with $m=0$ (shown in brown in figure 3).

Figure 2. Quasi-axisymmetry example. (a) Flux surface shape computed by the procedure of §§ 3 and 4.1, taking aspect ratio $=$  10, showing $|B|$ computed by VMEC. (b) Cross-sections of the flux surfaces at equally spaced values of $\unicode[STIX]{x1D719}$ , with $+$ signs denoting the magnetic axis.

Figure 3. Fourier amplitudes $B_{m,n}(r)$ of the magnetic field magnitude $B(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$ computed by BOOZ_XFORM, for the quasi-axisymmetric configuration of § 5.1.

The theory here generates flux surface shapes that give quasisymmetry to first order in the distance from the magnetic axis, and at next order in this distance there will be breaking of the symmetry. Therefore, the symmetry-breaking Fourier harmonics should scale as $1/A^{2}$ where $A$ is the aspect ratio. This scaling is verified in figure 4. In this figure the amount of symmetry breaking is measured by the quantity

(5.2) $$\begin{eqnarray}S=\frac{1}{B_{0,0}}\sqrt{\displaystyle \mathop{\sum }_{n/m\neq N/M}B_{m,n}^{2}}.\end{eqnarray}$$

As expected, the symmetric modes $B_{m,n}$ are found to scale as $1/A$ (not shown). Similarly, figure 5 shows that the rotational transform computed by VMEC converges to the value predicted by (2.7) as the aspect ratio increases. For $A\geqslant 160$ , the agreement extends to at least 5 digits.

Figure 4 includes a point for the quasi-axisymmetric design NCSX, which was obtained using conventional optimization. The NCSX point falls below the trend line, so it evidently has a somewhat better quality of quasisymmetry for its aspect ratio than the configurations constructed here.

Figure 4. For all three examples presented in § 5, the symmetry-breaking Fourier components scale as $A^{-2}$ as predicted by theory.

Figure 5. For all three examples presented in § 5, the rotational transform computed by VMEC converges to the value predicted by (2.7) as the aspect ratio increases.

Figure 6. The aspect ratio 5 quasi-axisymmetic stellarator constructed by the procedure of § 4.4, using no optimization (aside from the REGCOIL linear least-squares problem). (a,b) Colour indicates $B$ on the outermost flux surface, and the four unique coil shapes are shown with four shades of grey. (c–e) Poincaré plots computed from the vacuum field of the coils, demonstrating good flux surfaces out to aspect ratio 5, at three toroidal angles. (f) Boozer spectrum, demonstrating the quasi-axisymmetric mode is dominant. (g) Profile of $\unicode[STIX]{x1D704}$ .

Figure 7. Quasi-helical symmetry example. (a) Flux surface shape computed by the procedure of §§ 3 and 4.2, taking aspect ratio $=$  40, showing $|B|$ computed by VMEC. (b) Cross-sections of the flux surfaces at equally spaced values of $\unicode[STIX]{x1D719}$ , with $+$ signs denoting the magnetic axis.

For a different approach to constructing a finite-aspect-ratio geometry from the high-aspect-ratio theory, an example of the outward extrapolation method of § 4.4 is shown in figure 6, again using the input axis shape (5.1). First, an aspect ratio 160 shape is generated by the method of § 4.1. Coil shapes to produce this magnetic surface shape are then calculated using the REGCOIL method (Landreman 2017). For this method, a coil winding surface is chosen by taking an aspect ratio 5 surface constructed using the method of § 4.1, and expanding uniformly outward by one quarter of the average major radius. REGCOIL’s regularization parameter is chosen to be the smallest value for which there are no saddle coils, i.e. there are no local maxima or minima in the current potential. Next, 24 coil shapes (four unique shapes, each repeated six times) are identified from uniformly spaced contours of the current potential. A Poincaré plot of the vacuum field produced by these coils (figure 6 c–e) shows that good flux surfaces exist out to an aspect ratio of 5.0 (using VMEC’s definition of the major and minor radius). The Fourier amplitudes of $B$ in Boozer coordinates are shown in figure 6(f), showing the quasi-axisymmetric term is dominant, as desired. Again the largest symmetry-breaking mode at the edge is the mode $(m,n)=(2,-3)$ . The symmetry-breaking harmonics reach a rather sizeable amplitude at the last closed flux surface, and no effort has been made to achieve other desirable physics properties such as a high magnetohydrodynamic $\unicode[STIX]{x1D6FD}$ limit. However, this configuration required very little computational effort to compute, compared to the hundreds or thousands of VMEC computations required for conventional optimization, and it could serve as a useful initial condition for conventional optimization.

5.2 Quasi-helical symmetry

For an example of quasi-helical symmetry, we consider the magnetic axis shape

(5.3a,b ) $$\begin{eqnarray}R_{0}(\unicode[STIX]{x1D719})=1+0.265\cos (4\unicode[STIX]{x1D719}),\quad z_{0}(\unicode[STIX]{x1D719})=-0.21\sin (4\unicode[STIX]{x1D719}).\end{eqnarray}$$

For this curve, the normal vector rotates poloidally in each field period, so solutions have quasi-helical symmetry rather than quasi-axisymmetry. We also choose $\bar{\unicode[STIX]{x1D702}}=-2.25$ and $\unicode[STIX]{x1D70E}(0)=0$ . For these parameters, the numerical procedure of §§ 3–4.1 yields a rotational transform $\unicode[STIX]{x1D704}=1.93$ , and the maximum flux surface elongation in the $R$ – $z$ plane is found to be 2.52.

The flux surfaces for aspect ratio 40, computed using the method of § 4.2, are shown in figure 7. Due to the strongly shaped axis in this example, the flux surface cross-sections in the $R$ – $z$ plane become visibly different from ellipses even at this high aspect ratio. Note that the cross-sections in the plane perpendicular to the magnetic axis are perfectly elliptical, and the cross-sections in the $R$ – $z$ plane approach ellipses as the aspect ratio is raised. The spectra of $B$ in Boozer coordinates for aspect ratios 40 and 160 are shown in figure 8. The largest symmetry-breaking mode at the edge is the mode $(m,n)=(2,-16)$ .

As pointed out by Garren & Boozer (1991a ), the relevant ratio for breaking of quasisymmetry is the minor radius divided by the scale length of the magnetic axis’ Frenet–Serret frame (e.g. $1/\unicode[STIX]{x1D705}$ , $1/\unicode[STIX]{x1D70F}$ ), not the conventional aspect ratio. For axis shapes consistent with quasi-helical symmetry, where the normal vector rotates about the axis, the scale lengths of the axis Frenet–Serret frame are smaller than for axes consistent with quasi-axisymmetry at comparable major radius. Therefore, quasi-helical symmetry is limited to higher conventional aspect ratios than quasi-axisymmetry. This trend is apparent in a comparison of the examples of §§ 5.1–5.2. The peak axis curvature and torsion are approximately twice as large in the latter compared to the former, and this ratio is squared in the symmetry breaking. Indeed, figure 3(a) for quasi-axisymmetry at aspect ratio 10 has comparable symmetry breaking to figure 8(a) for quasi-helical symmetry at aspect ratio 40.

Figure 8. Fourier amplitudes $B_{m,n}(r)$ of the magnetic field magnitude $B(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$ computed by BOOZ_XFORM, for the quasi-helically symmetric configuration of § 5.2.

Figure 4 includes a point for the quasi-helically symmetric experiment HSX, which was designed using conventional optimization. (Coil ripple is not included for the HSX and NCSX configurations in the figure; the values of $S$ for HSX and NCSX are nearly unchanged on the scale of the figure if coil ripple is included.) HSX has symmetry breaking that is an order of magnitude smaller than the configuration generated here at comparable aspect ratio. The fact that conventional optimization results in lower symmetry breaking than the construction here is not surprising, given that the construction is limited to producing shapes with elliptical cross-section.

Figure 9. Quasi-axisymmetric stellarator without stellarator symmetry. (a) Flux surface shape computed by the procedure of §§ 3 and 4.1, taking aspect ratio $=$  10, showing $|B|$ computed by VMEC. (b) Cross-sections of the flux surfaces at equally spaced values of $\unicode[STIX]{x1D719}$ , with $+$ signs denoting the magnetic axis.

Figure 10. Fourier amplitudes $B_{m,n}(r)$ of the magnetic field magnitude $B(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$ computed by BOOZ_XFORM, for the non-stellarator-symmetric quasi-axisymmetric configuration of § 5.3.

5.3 Case without stellarator symmetry

There is no reason a stellarator with quasisymmetry must also possess stellarator symmetry. For instance, a tokamak with a single null is quasi-axisymmetric but not stellarator symmetric. Plasma shapes that lack stellarator symmetry are of interest since the turbulent momentum flux is predicted to be larger by a factor ${\sim}1/\unicode[STIX]{x1D70C}_{\ast }$ than in stellarator-symmetric shapes, meaning the intrinsic rotation is larger (Peeters & Angioni 2005; Parra, Barnes & Peeters 2011; Sugama et al. 2011). The resulting rotation and/or rotation shear may improve plasma stability. While quasisymmetry reduces the strong damping of flows otherwise typical of stellarators, significant flow still requires a drive, and turbulent momentum transport associated with broken stellarator symmetry could provide such a drive. In the model considered here, stellarator symmetry can be broken by specifying a non-stellarator-symmetric axis shape, or by specifying a non-zero $\unicode[STIX]{x1D70E}(0)$ , or both. Here we present an example with both sources of symmetry breaking. We take the magnetic axis shape to be

(5.4a,b ) $$\begin{eqnarray}R_{0}(\unicode[STIX]{x1D719})=1+0.042\cos (3\unicode[STIX]{x1D719}),\quad z_{0}(\unicode[STIX]{x1D719})=-0.042\sin (3\unicode[STIX]{x1D719})-0.025\cos (3\unicode[STIX]{x1D719}),\end{eqnarray}$$

with $\bar{\unicode[STIX]{x1D702}}=-1.1$ and $\unicode[STIX]{x1D70E}(0)=-0.6$ . For these parameters, the numerical procedure above yields a rotational transform $\unicode[STIX]{x1D704}=0.311$ , and the maximum flux surface elongation in the $R$ – $z$ plane is found to be 3.29. The flux surface shape for $A=10$ is displayed in figure 9, and the Boozer spectra for $A=10$ and $A=80$ are shown in figure 10. In figure 10, it can be seen that $B$ has a significant $\sin \unicode[STIX]{x1D703}$ component (red dotted line) which is not stellarator symmetric but which preserves quasi-axisymmetry. As with the stellarator-symmetric quasi-axisymmetric example, the largest symmetry-breaking mode at the edge is the mode $(m,n)=(2,-3)$ . The $1/A^{2}$ scaling of the quasisymmetry-breaking harmonics is again plotted in figure 4, and the convergence of the VMEC rotational transform to the predicted value as $A\rightarrow \infty$ is shown in figure 5. Generally the properties of this family of configurations are quite similar to the stellarator-symmetric and quasi-axisymmetric configurations of § 5.1.

6 Discussion and conclusions

While quasisymmetric stellarator shapes have been found previously by applying black-box optimization methods to minimize the departure from quasisymmetry, such methods are computationally demanding and do not provide comprehensive information about the landscape of all possible solutions. Here we have demonstrated a complementary approach in which quasisymmetric stellarator shapes can be parameterized and computed extremely rapidly ( ${<}1~\text{ms}$ on a laptop), enabling exhaustive high-resolution parameter scans and insight into the size of the solution space. We have demonstrated that this approach is a practical method to generate configurations that can be examined numerically with VMEC and other physics codes. As part of this demonstration, we have shown using BOOZ_XFORM that the Boozer coordinate Fourier spectra of the resulting equilibria are indeed dominated by a single helicity. We have further demonstrated that the symmetry-breaking harmonics scale as $r^{2}$ as expected.

Although the ‘optimization-free’ approach here requires solving a nonlinear equation (2.7), the numerical solution is extremely robust since the problem can be formulated so a unique solution is guaranteed to exist. As proved in appendix A, every magnetic axis shape with non-vanishing curvature admits an infinity of (first-order) quasisymmetric flux surface shapes surrounding it, each labelled by the three numbers $\bar{\unicode[STIX]{x1D702}}$ , $I_{2}$ and $\unicode[STIX]{x1D70E}(0)$ . Given any values for these three numbers, as long as $\bar{\unicode[STIX]{x1D702}}\neq 0$ , and given any axis shape for which the curvature does not vanish, there is exactly one first-order quasisymmetric shape (as a function of $r$ ). However, much of this solution space is not interesting since the elongation of the surfaces is impractically high. The space of stellarator-symmetric solutions is significantly smaller than the space of all solutions both because the space of axis shapes is restricted and also since $\unicode[STIX]{x1D70E}(0)$ must be 0.

While the calculations here are limited to high aspect ratio, any stellarator with a low-aspect-ratio boundary will have a region close to the magnetic axis in which the local aspect ratio is high. Hence the results here describe the core of any quasisymmetric stellarator, even those with low aspect ratio at the plasma boundary. The observed accuracy of the solutions constructed here is consistent with the popular wisdom that good quasi-axisymmetry can be achieved at much lower aspect ratio than quasi-helical symmetry. It is likely that by extending the method here to second order in $r$ , using equations in the appendix of Garren & Boozer (1991b ), the accuracy of these parameterized solutions could be extended to lower aspect ratio.