2 Near-axis expansion

The NAE, as described by Garren & Boozer (Reference Garren and Boozer1991), solves the equilibrium MHD equations by performing a first-order Taylor expansion about the magnetic axis in pseudo-Cartesian coordinates (Kuo-Petravic & Boozer Reference Kuo-Petravic and Boozer1987). When using Boozer coordinates $(\psi,\theta,\varphi )$ the general form of the magnetic field at first order in the distance from the axis is given by

(2.1)\begin{equation} B(\epsilon, \theta, \varphi) \approx B_{0}(\varphi) + \epsilon B_{1}(\theta,\varphi)= B_{0}(\varphi) \left(1 + \epsilon d(\varphi) \cos[\theta - \alpha(\varphi) ] \right), \end{equation}

with the expansion parameter

(2.2)\begin{equation} \epsilon=\sqrt{\psi}\ll 1, \end{equation}

being a measure of the distance from the axis. As a consequence, these solutions will only be valid in the large aspect ratio limit. The first-order correction to the magnetic field $B_{1}$ has the general form of an analytical function close to the axis. At this order, the solutions correspond to the vacuum case and magnetic surfaces are always guaranteed, as discussed in Jorge, Sengupta & Landreman (Reference Jorge, Sengupta and Landreman2020).

At first order, the spatial position is given by

(2.3)\begin{equation} \boldsymbol{x} \approx \boldsymbol{r}_{0} + \epsilon \left( X_{1}\boldsymbol{n}^{s} + Y_{1}\boldsymbol{b}^{s} \right), \end{equation}

where $\boldsymbol {r}_{0}$ describes the position of the magnetic axis, a space curve characterised by its torsion $\tau$ and signed curvature $\kappa ^{s}$. The signed Frenet–Serret frame, $(\boldsymbol {t}, \boldsymbol {n}^{s}, \boldsymbol {b}^{s})$, of this curve is used as an orthonormal basis to describe the coordinate mapping. The tangent vector $\boldsymbol {t}$ is identical to that of the traditional Frenet–Serret apparatus, and the normal and binormal signed vectors, $\boldsymbol {n}^{s}, \boldsymbol {b}^{s}$, as well as the signed curvature change signs at points where $\kappa = 0$.

When imposing the conditions for omnigenity in the NAE, as done by Plunk et al. (Reference Plunk, Landreman and Helander2019) and Plunk et al. (Reference Plunk, Landreman and Helander2021), an approximate plasma boundary corresponding to a QI magnetic equilibria can be described by

(2.4)\begin{gather} X_{1} = \frac{d(\varphi)}{\kappa^{s}} \cos{[\theta - \alpha (\varphi) ]} \end{gather}

(2.5)\begin{gather}Y_{1} = \frac{2\kappa^{s}}{B_{0}(\varphi)d(\varphi)} \left( \sin{[\theta - \alpha (\varphi)]} + \sigma(\varphi) \cos{[\theta - \alpha (\varphi) ]} \right), \end{gather}

where $B_{0}$ is the on-axis magnetic field strength, which must vary with $\varphi$ because, as also shown by Plunk et al. (Reference Plunk, Landreman and Helander2019), at first order, constant on-axis magnetic field is incompatible with the omnigenity conditions. The functions $\alpha (\varphi )$ and $d(\varphi )$ need to be prescribed and are required to satisfy specific conditions to correspond to omnigenous solutions as will be discussed in more detail later in this work. The function $\sigma (\varphi )$ appearing in (2.5) can be seen to relate to the elongation of the boundary cross-section.

The problem of constructing QI magnetic equilibria is then transformed into the problem of specifying an on-axis magnetic field strength $B_{0}$, an axis shape, two functions $\alpha (\varphi )$ and $d(\varphi )$, and solving a differential equation for the quantity $\sigma (\varphi )$

(2.6)\begin{equation} \sigma^{\prime} + ({\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} -\alpha^{\prime})\left(\sigma^{2}+ 1 +\frac{B_{0}^{2}{\bar{d}}^{4}}{4} \right)-G_{0}\bar{d}^{2}\left(\tau + I_{2}/2\right) = 0, \end{equation}

which must be solved self-consistently with the rotational transform on axis ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$. Here $G_{0}$ and $I_2$ are related to the first non-zero terms of the poloidal and toroidal current functions on axis, respectively, as expanded in Garren & Boozer (Reference Garren and Boozer1991), and $\bar {d}(\varphi ) = d(\varphi ) /\kappa ^{s}(\varphi )$. Primes represent derivatives with respect to the toroidal angular coordinate $\varphi$. In the case of stellarator symmetry, $\sigma (0) = 0$.

In order to mathematically formulate the omnigenity conditions required on $\alpha (\varphi )$ and $d(\varphi )$, a certain coordinate mapping described in Cary & Shasharina (Reference Cary and Shasharina1997) was used in Plunk et al. (Reference Plunk, Landreman and Helander2019). For each magnetic well in the on-axis magnetic field, trapping domains are delimited by the maxima defining the well. This is divided into a right-hand domain $D_{R}$ and a left-hand domain $D_{L}$, depending on which side of the minimum of the well the points lie. For each point $\varphi$ in a right-hand domain we identify its corresponding bounce point $\varphi _{b}(\varphi )$ in the left-hand domain by the condition

(2.7)\begin{equation} B_{0}(\varphi) = B_{0}(\varphi_{b}(\varphi)), \end{equation}

and vice versa. It is clear from this construction that

(2.8)\begin{equation} \varphi_{b}(\varphi_{b}(\varphi)) = \varphi. \end{equation}

The angular distance between an arbitrary point $\varphi$ and its corresponding bounce point is given by

(2.9)\begin{equation} \Delta \varphi (\varphi) \equiv \varphi - \varphi_{b}(\varphi). \end{equation}

Plunk et al. (Reference Plunk, Landreman and Helander2019) provided a way to construct the functions $d(\varphi )$ and $\alpha (\varphi )$ in $D_{R}$, giving their dependency in $D_{L}$, to guarantee omnigenity. This is done by equating the near-axis and the omnigenous forms of $B_{1}$

(2.10)\begin{equation} B_{0}(\varphi) d(\varphi) \cos{(\theta - \alpha(\varphi))} ={-} F_{1}(\theta,\varphi) B_{0}^{\prime}(\varphi), \end{equation}

where $F_{1}(\theta,\varphi )$, that comes from the expansion of the Cary–Shasharina mapping, is a small perturbation added to a zeroth-order quasi-symmetric magnetic field to satisfy omnigenity, and hence preserves the following symmetry:

(2.11)\begin{equation} F_{1}(\theta,\varphi) = F_{1}(\theta-{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \Delta\varphi(\varphi), \varphi_{b}(\varphi)),\quad \text{for } \varphi \in D_{R}. \end{equation}

Solving (2.10) for $F_{1}$ and plugging it into (2.11) we obtain the following set of conditions:

(2.12)\begin{gather} d(\varphi) = \varphi^{\prime}_{b}(\varphi) d(\varphi_{b}(\varphi)),\quad\text{for }\varphi\in D_{R} \end{gather}

(2.13)\begin{gather}\alpha(\varphi) = \alpha(\varphi_{b}(\varphi)) + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \Delta \varphi(\varphi),\quad\text{for }\varphi\in D_{R}. \end{gather}

In addition, $d$ must vanish at all extrema of the on-axis magnetic field $B_{0}$, as seen from (2.10). The curvature of the magnetic axis $\kappa ^{s}$ needs to have zeros of the same order as $d$ at these points for the plasma boundary to be well described, i.e. so that $X_{1}$ and $Y_{1}$, which are proportional to $d/\kappa ^{s}$ and $\kappa ^{s}/d$, respectively, remain non-zero and bounded.

It is important to note that periodicity cannot be enforced if the condition (2.13) on $\alpha (\varphi )$ is satisfied and the rotational transform ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$ is irrational. We can see this from evaluating (2.13) at the maximum $\varphi = 2{\rm \pi}$,

(2.14)\begin{equation} \alpha(2{\rm \pi}) -\alpha(0) = 2{\rm \pi} {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}. \end{equation}

However, continuity of $B_{1}$, $X_{1}$ and $Y_{1}$ requires

(2.15)\begin{equation} \alpha(2{\rm \pi}) - \alpha(0) = 2{\rm \pi} M, \end{equation}

with $M$ defined as the number of times the axis curvature vector $\boldsymbol {n}^{s}$ rotates during one toroidal transit. Thus, (2.14) is generally in conflict with (2.15) and omnigenity is only consistent with continuity for integer values of ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$. Plunk et al. (Reference Plunk, Landreman and Helander2019) resolved this conflict by introducing small matching regions around $\varphi =0$ and $\varphi = 2{\rm \pi}$, where omnigenity is abandoned and $\alpha$ is defined to guarantee periodicity. A different approach to solving this problem, as well as the form these conditions take for a single-well, $N$-field-period, stellarator-symmetric configuration is discussed in the following sections.

3 $N$-field periods and stellarator symmetry

A magnetic field is said to be stellarator-symmetric if it is invariant, except for sign change, under a $180^\circ$ rotation (the operations $I_n$ defined in the following) around a set of axes. It is conventional to take these axes to lie within the $x\unicode{x2013}y$ plane, passing through the origin, and regularly spaced in the toroidal angle. This type of symmetry was defined formally by Dewar & Hudson (Reference Dewar and Hudson1998), who introduced a symmetry operator $I_{0}$ by

(3.1)\begin{equation} I_{0}\boldsymbol{F}(\psi,\theta,\varphi) = \boldsymbol{F}(\psi,-\theta,-\varphi), \end{equation}

with $\theta$ and $\varphi$ being angular coordinates. We say that a vector $\boldsymbol {F}$ possesses stellarator symmetry if

(3.2)\begin{equation} I_{0}[F_{\psi},F_{\theta},F_{\varphi}] = [{-}F_{\psi},F_{\theta},F_{\varphi}], \end{equation}

where $F_{j}={\partial {\boldsymbol {x}}}/{\partial j} \boldsymbol {\cdot } \boldsymbol {F}$, with $j=\{\psi,\theta,\varphi \}$, are the covariant components of $\boldsymbol {F}$. For the case of a scalar quantity, for instance $|F|$, stellarator symmetry implies

(3.3)\begin{equation} I_{0}|F| = |F|. \end{equation}

If the vector field $\boldsymbol {F}$ possesses $N$-fold discrete symmetry about the $z$-axis, then stellarator symmetry also exists about the cylindrical inversion symmetry operation with respect to the half-line $\{\phi = 2i{\rm \pi} /N, Z= 0\}$, with $i= 1,2,\ldots,(N-1)$, and with respect to the half-line in the middle of each field period $\{\phi = (2i-1){\rm \pi} /N, Z= 0\}$, with $i= 1,2,\ldots,N$. These symmetry operators will be referred to as $I_{2i}$ and $I_{2i-1}$, respectively.

Let us now consider the case of a stellarator-symmetric $N$-field-period QI configuration, and focus on the case of one magnetic well per field period.Footnote ³ This symmetry is not a necessary condition but is used to reduce the complexity of the problem and is conventionally assumed in QI designs.

In order to be stellarator symmetric, the magnetic field strength ${B}$ needs to fulfil

(3.4)\begin{equation} B(\theta,\varphi_{\text{min}}^{(i)}+ \delta \varphi) = B(-\theta,\varphi_{\text{min}}^{(i)}- \delta \varphi), \end{equation}

which is obtained by invoking the symmetry operator $I_{2i-1}$. The angular position of the $i$th minimum, $\varphi _{\mathrm {min}}^{(i)}$ is located on the rotation axis, and is given by

(3.5)\begin{equation} \varphi_{\mathrm{min}}^{(i)} = \frac{(2i-1){\rm \pi}}{N},\quad\text{for }i = 1,\ldots,N \end{equation}

whereas the $i$th maximum is

(3.6)\begin{equation} \varphi_{\mathrm{max}}^{(i)} = \frac{2(i-1){\rm \pi}}{N},\quad\text{for }i = 1,\ldots,N. \end{equation}

The trapping domain in each period is labelled with the subscript $i$, shown in figure 1, and the left- and right-hand domains are defined as

(3.7)\begin{equation} \left.\begin{array}{ll@{}} D_{iL}(\varphi), & \text{for }\dfrac{2{\rm \pi}(i-1)}{N}\leqslant \varphi \leqslant \dfrac{(2i-1){\rm \pi}}{N} \\ D_{iR}(\varphi), & \text{for } \dfrac{(2i-1){\rm \pi}}{N}\leqslant \varphi \leqslant \dfrac{2{\rm \pi} i}{N} \end{array}\right\} \end{equation}

Now, (2.9), which gives the distance between bounce points, can be written as

(3.8)\begin{equation} \Delta\varphi(\varphi) = 2(\varphi - \varphi_{\mathrm{min}}^{(i)}), \end{equation}

and the bounce points can be found using

(3.9)\begin{equation} \varphi_{b}(\varphi) = \varphi - \Delta \varphi(\varphi) = 2\varphi_{\mathrm{min}}^{(i)} - \varphi. \end{equation}

Using this new notation, we can write the symmetry operator $I_{2i-1}$ as

(3.10)\begin{equation} I_{\varphi_{\mathrm{min}}^{i}} f(\psi,\theta,\varphi)= f(\psi,-\theta,\varphi_{b}). \end{equation}

Figure 1. One-well magnetic field strength example, plotted over a single toroidal field period. The right- and left-hand domains, $D_{iR}$ and $D_{iL}$, are indicated. The angular position of the minimum and maxima of the well are labeled by $\varphi _{\mathrm {min}}^{i}$ and $\varphi _{\mathrm {max}}^{i}$, respectively.

Let us now focus on the condition for omnigenity on the function $d(\varphi )$, which is shown in (2.12). We substitute $\varphi _{b}$ from expression (3.9) and note that $\varphi _{b}^{\prime }=-1$, to find

(3.11)\begin{equation} d(\varphi) ={-}d(\varphi_{b})={-} d(2\varphi_{\mathrm{min}}^{(i)} - \varphi),\quad\text{for }\varphi\in D_{iR}. \end{equation}

As a consequence, $d(\varphi )$ must necessarily be an odd function about $\varphi _{\mathrm {min}}^{(i)}$, the bottom of the well. Note that $\varphi _{b}^{\prime }(\varphi )$ is generally negative near bounce points for $\varphi \in D_{iL}$, even without requiring stellarator symmetry, so the order of the zeros of $d$ must always be odd, and therefore the same is true for those of $\kappa ^{s}$.

For the case of $\alpha (\varphi )$, we first define a quantity $\Delta \alpha$ in the same fashion as (2.9)

(3.12)\begin{equation} \Delta \alpha(\varphi) \equiv \alpha(\varphi) - \alpha(\varphi_{b}(\varphi)). \end{equation}

By replacing $\alpha$ by expression (2.13) in the previous definition, we obtain

(3.13)\begin{equation} \Delta \alpha(\varphi) = {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}\Delta\varphi(\varphi), \end{equation}

and using expression (3.8) we get

(3.14)\begin{equation} \Delta \alpha(\varphi) = 2{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} (\varphi - \varphi_{\mathrm{min}}^{(i)}). \end{equation}

We see that $\alpha$ needs to have a part proportional to ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \Delta \varphi (\varphi )$, to guarantee the right form of $\Delta \alpha$, and an even part $\alpha _{e}$ such that

(3.15)\begin{equation} \alpha_{e}(\varphi) = \alpha_{e}(\varphi_{b}(\varphi)). \end{equation}

Accordingly, we write $\alpha$ as

(3.16)\begin{equation} \alpha(\varphi) = \tfrac{1}{2}{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}\Delta\varphi(\varphi) + \alpha_{e}(\varphi). \end{equation}

In order to find $\alpha _{e}$, we note that the first-order correction to the magnetic field, given in (2.1),

(3.17)\begin{equation} B_{1} = B_{0}(\varphi) d(\varphi)\cos{(\theta-\alpha(\varphi))} \end{equation}

needs to possess stellarator symmetry. We know from (3.11) that $d(\varphi )$ is an odd function and, by construction, that $B_{0}$ is an even function. Therefore, $\cos {(\theta -\alpha (\varphi ))}$ must be odd under the symmetry operation defined in (3.10). In particular, this requires

(3.18)\begin{equation} \cos{(\theta_{0} - \alpha(\varphi_{\mathrm{min}}^{i}))}=0, \end{equation}

at the stellarator symmetry point $(\theta _0 = 0, \varphi _{\mathrm {min}}^i)$, which is valid for any integer $n_{i}$ such that

(3.19)\begin{equation} \alpha(\varphi_{\mathrm{min}}^{i}) = {\rm \pi}(n_{i}+ \tfrac{1}{2}). \end{equation}

Evaluating (3.16) at the minimum, we find

(3.20)\begin{align} \alpha_{e} = {\rm \pi}(n_{i}+ \tfrac{1}{2}), \end{align}

so that $\alpha$ takes the form

(3.21)\begin{equation} \alpha(\varphi) = {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} (\varphi - \varphi_{\mathrm{min}}^{(i)}) + {\rm \pi}(n_{i}+ \tfrac{1}{2}). \end{equation}

Now let us impose a necessary (but not sufficient) condition for periodicity of $B$ from one field period to the next, i.e.

(3.22)\begin{equation} \alpha(\varphi_{\rm max}^{i+1})-\alpha(\varphi_{\rm max}^{i}) = 2{\rm \pi} m, \end{equation}

where $m$ is the number of times the signed curvature vector ${\boldsymbol {n}}^{s}$ rotates per field period. This ensures that the poloidal angle $\theta$ is periodic and increases by $2{\rm \pi}$ after a full toroidal rotation. The addition of the term $2{\rm \pi} m$ to $\alpha$ compensates the poloidal rotation of the axis (measured by $m$) because $\alpha$ effectively behaves as a phase shift on $\theta$. The previous expression translates into a relation for $n_{i}$

(3.23)\begin{equation} {\rm \pi}(n_{i+1}+1/2)-{\rm \pi}(n_{i}+1/2) = 2{\rm \pi} m, \end{equation}

which can also be written as

(3.24)\begin{equation} n_{i+1} = 2im+n_{1}. \end{equation}

Without loss of generality we can choose $n_{1}=0$, and inserting this result in (3.21) we obtain an expression for $\alpha (\varphi )$ satisfying omnigenity and $N$-fold periodicity

(3.25)\begin{equation} \alpha(\varphi) = {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} (\varphi-\varphi_{\mathrm{min}}^{i})+{\rm \pi} (2 i m + \tfrac{1}{2}). \end{equation}

The shape of the boundary must also be made periodic, which can be achieved in the same way as done in Plunk et al. (Reference Plunk, Landreman and Helander2019), by introducing so-called buffer regions around the maxima of $B_{0}$. In these regions, $\alpha$ is not calculated using (3.25), but is instead chosen to give a periodic plasma boundary. The function in the buffer region needs to be constructed carefully to make sure the function itself and its derivatives are continuous and smooth in order to avoid numerical difficulties in equilibrium solvers such as VMEC (Hirshman & Whitson Reference Hirshman and Whitson1983).

We note that although we have some freedom in the choice of $n_{1}$ and, hence, in the value of $\alpha$, the term entering (2.6) and defining the equilibria is $\alpha ^{\prime }(\varphi )$, which is independent of the choice of $n_{1}$. The only freedom left in the choice of $\alpha (\varphi )$ is on how omnigenity is broken to impose continuity of the solutions.

4 Smoother $\alpha (\varphi )$

In order to avoid the problems derived from defining $\alpha (\varphi )$ as a piecewise function, we propose a different approach, namely by adding an omnigenity-breaking term to (3.25) that allows periodic solutions

(4.1)\begin{equation} \alpha(\varphi) = {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} (\varphi-\varphi_{\mathrm{min}}^{i})+{\rm \pi} (2\,m i + \tfrac{1}{2}) + a(\varphi-\varphi_{\mathrm{min}}^{i})^{2k+1}. \end{equation}

The last term in this expression goes to zero at the bottom of the well as long as $k\geqslant -1/2$ and, hence, makes the solution omnigenous for deeply trapped particles. The parameter $a$ can be chosen in such a way that periodicity of $\alpha (\varphi )$ is guaranteed. A first requirement is continuity at $\varphi _\textrm {max}^{(i)}$, i.e.

(4.2)\begin{equation} \lim_{\varphi\to\varphi_{\mathrm{max}}^{(i+1)-}} \alpha(\varphi) = \lim_{\varphi\to\varphi_{\mathrm{max}}^{(i+1)+}} \alpha(\varphi), \end{equation}

where we are considering the $i$th well when approaching from the left and the well labelled $i+1$ when approaching from the right. Hence, we obtain

(4.3)\begin{align} & {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} (\varphi_{\mathrm{max} }^{(i+1)}-\varphi_{ \mathrm{min}}^{(i)}) + a(\varphi_{\mathrm{max}}^{(i+1)}-\varphi_{\mathrm{min}}^{(i)})^{2k+1}\nonumber\\ & \quad = {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} (\varphi_{\mathrm{max}}^{(i+1)}-\varphi_{\mathrm{min}}^{(i+1)}) + a(\varphi_{\mathrm{max}}^{(i+1)}-\varphi_{\mathrm{min}}^{(i+1)})^{2k+1} + 2{\rm \pi} mi, \end{align}

and when inserting the values of $\varphi _{\mathrm {max} }^{(i)}$ and $\varphi _{\mathrm {min} }^{(i)}$

(4.4)\begin{equation} {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \left( - {\rm \pi}/ N\right) + a\left(-{\rm \pi} / N\right)^{2k+1} + 2{\rm \pi} m = {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \left({\rm \pi} / N\right) + a\left({\rm \pi} / N\right)^{2k+1}, \end{equation}

we find an expression for $a$

(4.5)\begin{equation} a = \frac{{\rm \pi} \left( m - {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} /N \right)}{({\rm \pi}/N )^{2k+1}}, \end{equation}

which finally gives us a form of $\alpha$ that breaks omnigenity in a smooth and controlled way

(4.6)\begin{equation} \alpha(\varphi) = {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} (\varphi-\varphi_{\mathrm{min}}^{i})+{\rm \pi} \left(2\,m i + \frac{1}{2}\right) + {\rm \pi}\left( m - {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} /N \right) \left(\frac{\varphi-\varphi_{\mathrm{min} }^{i}}{{\rm \pi}/N} \right)^{2k+1}. \end{equation}

In figure 2, we can see the effect the choice of the parameter $k$ has over the shape of $\alpha (\varphi )$ for an axis shape with $m=0$. It is clear that increasing $k$ results in a function closer to that required for omnigenity but at the cost of a sharp behaviour around $\varphi _{\mathrm {max}}$ to preserve periodicity.

Figure 2. Behaviour of the function $\alpha (\varphi )$ over one period, for the case $m=0$. The dotted line shows the fully omnigenous case, which is non-continuous at $\varphi _{\mathrm {max}}$. Solid lines correspond to $\alpha$ given by (4.6) with values for the parameter $k$ ranging from 1 to 8.

The equilibrium constructed using this $\alpha$ will be approximately omnigenous as long as the last term in (3.25) is small relative to the first term, i.e.

(4.7)\begin{equation} \left| {\rm \pi}\left( m - {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} /N \right) \left(\frac{\varphi-\varphi_{\mathrm{min} }^{i}}{{\rm \pi}/N} \right)^{2k+1} \right| \ll \left|{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \left( \varphi - \varphi_{\mathrm{min}}^{(i)} \right)\right|, \end{equation}

which can be further simplified and rearranged as

(4.8)\begin{equation} \left( \frac{\left| \varphi - \varphi_{\mathrm{min}}^{(i)}\right|}{{\rm \pi}/N} \right)^{2k} \ll \frac{1}{\left| Nm/{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} -1\right|}. \end{equation}

The left-hand side of this expression includes $| \varphi - \varphi _{\mathrm {min}}^{(i)}|$, which is always smaller than ${\rm \pi} /N$ in a well. Hence, the condition for omnigenity being achieved closely everywhere in the well, i.e. also for $\varphi -\varphi _{\mathrm {min}} \approx {\rm \pi}/N$, reduces to

(4.9)\begin{equation} \left| Nm/{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} -1 \right| \ll 1. \end{equation}

This implies $N m \sim {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$, which, as expected can only be satisfied for nearly rational ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$. Although a rational value of ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$ is inconsistent with confinement, it may be advantageous to seek nearly rational values as a strategy for finding almost omnigenous configurations. We also note that the $m=0$ case, in which (4.9) is not technically satisfied, is however interesting because the limit of small ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$ implies that the absolute size of $\alpha$ will remain small, as required for approximate omnigenity.

A yet smoother choice of the function $\alpha (\varphi )$ with continuous derivatives up to third order at $\varphi _{\mathrm {max}}$ can be achieved by adding an extra term,

(4.10)\begin{equation} \alpha_{{\mathrm{II}}}(\varphi) = {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} (\varphi-\varphi_{\mathrm{min}}^{i})+{\rm \pi} (2\,m i + \tfrac{1}{2}) + a(\varphi-\varphi_{\mathrm{min}}^{i})^{2k+1} + b(\varphi-\varphi_{\mathrm{min}}^{i})^{2p+1}. \end{equation}

For the configurations shown in this work we use $\alpha (\varphi )$ as described in (3.25), which appears sufficiently smooth for the cases we have studied, but more details about $\alpha _{\mathrm {II}}(\varphi )$ are discussed in Appendix A.

5 Axis shape

The shape of the magnetic axis is perhaps the most important input for the construction of stellarator configurations since the axis properties seem to strongly affect the success of the construction, as measured by the accuracy of the approximation at finite aspect ratio.

The need to have an axis with points of zero curvature has already been discussed, but, in addition, low-curvature axes are attractive because they improve the accuracy of the near-axis approximation (indeed, the original work of Garren and Boozer defined the expansion parameter in terms of the maximum curvature of the magnetic axis), and are also associated with a low-amplitude first-order magnetic field

(5.1)\begin{equation} B_1 = B_0 \bar{d} \kappa^{s} \cos(\theta - \alpha), \end{equation}

where the definition of $\bar {d}$ has been substituted in (3.17). Note that $\bar {d}$ cannot be made too small without causing large elongation, and therefore minimising $\kappa ^s$ is an effective strategy for improving omnigenity at finite $\epsilon$.

It is also desirable to limit axis torsion, which enters directly into (2.6) for $\sigma$, as will be clear in the following discussion. At first order, the cross-sections of the plasma boundary at constant $\varphi$ are elliptical, with elongation $e$ defined as the ratio between the semi-major and semi-minor axes. Landreman & Sengupta (Reference Landreman and Sengupta2018) derived an expression for elongation in terms of $X_{1}$ and $Y_{1}$ (B 4). Using (2.4) and (2.5), we obtain an expression for $e$ dependent on the input parameters of the construction, namely $\bar {d}$, $B_{0}$ and $\sigma$,

(5.2)\begin{equation} e(\varphi) = \frac{B_{0}\bar{d}^{2}}{4} + \frac{1}{B_{0}\bar{d}^2}(1+\sigma^2) + \sqrt{\left( \frac{B_{0}\bar{d}^{2}}{4}\right)^2 + \frac{(\sigma^2-1)}{2} + \left(\frac{1}{B_{0}\bar{d}^2}(1+\sigma^2)\right)^{2}}. \end{equation}

We can simplify the previous expression by introducing

(5.3)\begin{equation} \bar{e}(\varphi) = \frac{B_{0}\bar{d}^{2}}{2}, \end{equation}

leading to

(5.4)\begin{equation} e = \frac{\bar{e}}{2} + \frac{1}{2\bar{e}}(1+\sigma^2) + \sqrt{ \frac{\bar{e}^{2}}{4} + \frac{(\sigma^2-1)}{2} + \frac{(1+\sigma^2)}{4\bar{e}^2}}. \end{equation}

One particular interesting limit is the idealised case of constant elongation. This can be achieved by choosing $\sigma$ constant. Given that, due to stellarator symmetry, $\sigma (0)=0$ then $\sigma$ must be zero everywhere, which transforms (5.4) into

(5.5)\begin{equation} e = \frac{1}{2\bar{e}}(\bar{e}^{2} +1 + \left|\bar{e}^{2}-1\right|), \end{equation}

this will result in a plasma boundary with constant elongation as long as $\bar {e}$ is independent of the toroidal angle $\varphi$. Finding plasma equilibria with constant elongation can be achieved by choosing $d(\varphi )$ appropriately so that it cancels the toroidal dependence of $\bar {e}$, namely

(5.6)\begin{equation} \bar{d}(\varphi) \propto \sqrt{\frac{2}{B_{0}(\varphi)}}. \end{equation}

Now, we can introduce the conditions that lead to constant elongation in (2.6)

(5.7)\begin{equation} ({\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} -\alpha^{\prime})\left( 1 +\bar{e} \right)= \frac{2G_{0}\bar{e}}{B_{0}} \tau, \end{equation}

where we have specialised to the case of vanishing current density on axis, $I=0$, the standard situation for QI stellarators (Helander & Nührenberg Reference Helander and Nührenberg2009; Helander, Geiger & Maaßberg Reference Helander, Geiger and Maaßberg2011). As omnigenity requires ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} - \alpha ^\prime \approx 0$ (see (3.25)), (5.7) implies that $\tau \approx 0$ is necessary for solutions with constant elongation. This shows that axes with low torsion are compatible with simple equilibrium boundary shapes.

The possibility of low elongation, in particular, may be unexpected, given what has been found with traditional QI optimisation, and also from considerations of basic geometry on the magnetic axis, where the identity

(5.8)\begin{equation} {\boldsymbol{n}} \kappa = \boldsymbol{\nabla}_\perp \ln B \end{equation}

follows from force balance and Ampere's law, noting that $\boldsymbol {\nabla } p = 0$ at $\psi = 0$, and defining $\boldsymbol {\nabla }_\perp = \boldsymbol {\nabla } - {\boldsymbol {tt}}\boldsymbol {\cdot } \boldsymbol {\nabla }$ (Landreman & Sengupta Reference Landreman and Sengupta2018). This expression seems to imply that plasma volumes that are compressed in the direction of the curvature (normal) vector, i.e. high-elongation configurations, will exhibit only weak variation of the field strength in the perpendicular direction, and therefore in $\theta$. This logic, although rough, is reflected in the near-axis construction of QI configurations, where it is typical to find elongated solutions. In § 8, we show how this tendency can be overcome by choice of axis torsion.

5.1 Axis construction

A space curve's shape is entirely determined by its curvature $\kappa$ and torsion $\tau$. From these two quantities, it is possible to calculate the tangent, normal and binormal vectors $(\boldsymbol {t}, \boldsymbol {n}, \boldsymbol {b})$ using the Frenet–Serret formulae

(5.9)\begin{equation} \left.\begin{gathered} \frac{{\rm d}\boldsymbol{t}}{{\rm d}s} = \kappa(s)\boldsymbol{n}, \\ \frac{{\rm d}\boldsymbol{n}}{{\rm d}s} ={-}\kappa(s)\boldsymbol{t} + \tau(s)\boldsymbol{b}, \\ \frac{{\rm d}\boldsymbol{b}}{{\rm d}s} ={-}\tau(s) \boldsymbol{n} , \end{gathered}\right\} \end{equation}

where $s$ is the arc length, used for parametrising the curve. Then, numerical integration of the tangent vector, $\boldsymbol {t} = \textrm {d}\boldsymbol {r}/\textrm {d}s$, yields the curve described by $\boldsymbol {r}$. Unfortunately, prescribing periodic $\kappa$ and $\tau$ is not sufficient for finding closed curves. A parameter scan is thus needed to find curves that can be used as magnetic axes. This is the approach used to find the axes curves described in § 7.

Although a Frenet description seems optimal for controlling torsion and curvature precisely, a truncated Fourier series is advantageous for its simplicity, smoothness and the fact that such curves are automatically closed. Smoothness is especially important for obtaining solutions that remain accurate at lower aspect ratio. Note that sharp derivatives $d/d\varphi \sim 1/\epsilon$ invalidate the NAE, which is a limit that is especially felt at high field period number. A method for generating simple axis curves, represented by a relatively small number of Fourier coefficients, is briefly outlined here and described in greater detail in Appendix B.

We represent the magnetic axis in cylindrical coordinates as

(5.10)\begin{equation} \boldsymbol{x} = \hat{\boldsymbol{R}} R(\phi) + \hat{\boldsymbol{z}} z(\phi). \end{equation}

The usual Fourier representation for a stellarator-symmetric axis is

(5.11)\begin{gather} R(\phi) = \sum_{n=0}^{n_\mathrm{max}} R_c(n) \cos(n N \phi), \end{gather}

(5.12)\begin{gather}z(\phi) = \sum_{n=1}^{n_\mathrm{max}} z_s(n) \sin(n N \phi). \end{gather}

A local form is used to establish conditions on the derivatives of these functions about a point of stellarator symmetry (also coinciding with an extrema of $B_0(\phi )$), that can be then used to generate a linear system of equations for the Fourier coefficients.

Conditions on torsion and curvature, specifically zeros of different orders, are imposed locally and converted into constraints on the derivatives of the axis components $R$ and $z$, and then applied to a truncated Fourier representation. This results in a set of linear conditions on the Fourier coefficients that can be solved numerically, or by computer algebra. The orders of the zeros, together with the set of Fourier coefficients, define a space that can be used for further optimisation.

As a very simple example of curves with points of vanishing curvature at multiples of ${\rm \pi} /N$, one may consider a symmetric class of curves, as in Plunk et al. (Reference Plunk, Landreman and Helander2019)

(5.13)\begin{gather} R = 1 + R_c(2) \cos(2 N \phi), \end{gather}

(5.14)\begin{gather}z = z_s(2) \sin(2 N \phi). \end{gather}

Owing to the fact that only the even mode numbers are retained ($z_s(1) = R_c(1) = 0$), the condition for zeros of curvature at first order need only be applied at $\phi = 0$, i.e. ${\textrm {d}^2R}/{\textrm {d}\varphi }|_{0} = R(0)$, and the result is

(5.15)\begin{equation} R_c(2)={-}\frac{1}{4 N^2+1}, \end{equation}

which matches (8.3a) of Plunk et al. (Reference Plunk, Landreman and Helander2019) for the case $N = 1$. The coefficient $Z_s(2)$ is free to be adjusted to satisfy other desired requirements. Note that additional Fourier modes may be retained to define a near-axis QI optimisation space, as done by Jorge et al. (Reference Jorge, Plunk, Drevlak, Landreman, Lobsien, Camacho and Helander2022).

5.2 Controlling torsion

Noting that a curve of zero torsion lies within a plane, it would seem straightforward to realise a stellarator-symmetric axis shape of low torsion by simply reducing the magnitude of its $z$ component, thereby constraining the curve to lie close to the $x$–$y$ plane. We can do this with the single-parameter curve defined by (5.13)–(5.14), by letting the parameter $z_s(2)$ tend to zero. Unfortunately this limit is not well-behaved, as shown in figure 3. Although torsion goes to zero almost everywhere, it tends to infinity in the neighbourhood of the zeros of curvature.

Figure 3. Torsion behaviour with decreasing value of $z_{s}(2)$: (a) singularity in torsion at points of stellarator symmetry as $z_{s}(2)$ approaches zero; (b) maximum torsion at different values of $z_{s}(2)$.

Another approach to minimising torsion is, somewhat paradoxically, to take the limit of large $z_s(2)$. Such elongated axis shapes are nearly planar around the points of zero curvature, and have their torsion peaked midway between these points. The torsion is, however, small in magnitude due to the large values of curvature around such points. As figure 4 confirms, the lowering of torsion by this method is accomplished only at the cost of raising the maximum value of curvature. Furthermore, the maximum value of torsion obtained at a fixed value of maximum curvature ($\mathrm {max}(\kappa ) = 3$) grows linearly with field period number, as shown in the second panel of figure 4. Requiring the maximum value of curvature to remain bounded when increasing the number of field periods results in the maximum of torsion increasing, as shown in figure 5. Finding closed curves with torsion and curvature remaining under a certain value gets more difficult when increasing the number of field periods. This might be a reason why finding good solutions with $N>1$ is challenging for the optimisation procedure described in Jorge et al. (Reference Jorge, Plunk, Drevlak, Landreman, Lobsien, Camacho and Helander2022).

Figure 4. (a) Maximum torsion versus maximum curvature for $N=2,4,8$ (blue, orange, green). Curvature and torsion are individually maximised over $\phi$ for fixed values of $z_{s}(2)$. Values used correspond to large-$z_{s}(2)$ regime, $z_{s}(2)\gtrsim 0.4/N^2$. (b) Maximum torsion versus $N$ for $z_{s}(2)$ chosen such that $\mathrm {max}_{\varphi }(\kappa ) = 3$.

Figure 5. (a) Curvature and (b) torsion versus $\phi$ for different field period numbers, $N=2,4,6,8,10$. Here $z_{s}(2)$ is chosen such that $\mathrm {max}_{\varphi }(\kappa ) = 3$.

6 QI construction, two-field-period example

The construction of a two-field-period configuration is now described. The axis shape was chosen using (5.15) for the case $N=2$, yielding

(6.1)\begin{gather} R = 1 - \tfrac{1}{17} \cos{(4\phi)}, \end{gather}

(6.2)\begin{gather}z = 0.3921 \sin{(2\phi)} + 4.90\times10^{{-}3}\sin{(4\phi)}. \end{gather}

The $z$-coefficients were chosen to limit the curvature and torsion to tolerable levels. The normal vector $\boldsymbol {n}^{s}$ does not complete any full rotation around the axis, hence $m=0$. This curve contains points of zero curvature as shown in figure 6. The location of these points coincide with the extrema of the on-axis magnetic field strength, which is chosen as

(6.3)\begin{equation} B_{0} = 1 + 0.15\cos{(2\varphi)}. \end{equation}

Figure 6. Unsigned curvature $\kappa$ (solid), torsion $\tau$ (dashed) and on-axis magnetic field intensity $B_0$ (dotted) profiles for one field period of the axis described by (6.1) and (6.2). The curvature has zeros at points of extrema of $B_{0}$ as required by the theory.

In general, the values of the toroidal coordinates $\varphi$, the Boozer angle, and $\phi$, the cylindrical angle, are not the same. However, thanks to stellarator symmetry, they coincide at the extrema points of $B_0(\varphi )$, which thus also correspond to the zeros of the axis curvature.

The choice of $d(\varphi )$ has an important effect on the elongation of the plasma boundary as can be seen from (2.4). We observed that keeping it proportional to $\kappa ^{s}(\varphi )$ helped to reduce the elongation to manageable levels. For this example it was chosen as $d(\varphi ) = 1.03\ \kappa ^{s}$. The parameter $k$ entering (4.6) controls the deviation from omnigenity and was set to $k=2$.

To construct numerical solutions, we first find the signed Frenet–Serret frame quantities of the axis: $\kappa ^s$, $\tau$ and $(\boldsymbol {t},\boldsymbol {n}^{\boldsymbol {s}},\boldsymbol {b}^{\boldsymbol {s}})$. Then, the relation $\varphi (\phi )$, as well as $G_{0}$, need to be found along the axis. This is done by iteratively solving (8.1) of Plunk et al. (Reference Plunk, Landreman and Helander2019)

(6.4a)\begin{equation} \frac{{\rm d}\varphi}{{\rm d}\phi}= \frac{B_0}{|G_0|}\left|\frac{{\rm d} \boldsymbol{r}}{{\rm d} \phi} \right|, \quad |G_0| = \frac{1}{2{\rm \pi} N}\int_0^{2{\rm \pi}}{\rm d}\phi \,B_0(\varphi(\phi))\left| \frac{{\rm d} \boldsymbol{r}}{{\rm d} \phi} \right|. \end{equation}

We then proceed to solve (2.6), self-consistently with ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$, for one field period, i.e. in the region $\varphi \in [0,2{\rm \pi} /N]$.

A boundary is constructed with aspect ratio $A=10$, where $A$ can be expressed in terms of the distance from the axis as

(6.5)\begin{equation} A=\sqrt{\frac{\overline{B_{0}}}{2}}\frac{R_c(0)}{\epsilon}, \end{equation}

where $\overline {B_0}$ is the average value of $B_{0}$. This boundary is then used to find a fixed-boundary magnetic equilibrium with the code VMEC. The strength of the magnetic field on the boundary is shown in figure 7. The rotational transform profile obtained with VMEC is shown in figure 8 and coincides with the value calculated numerically from (2.6) ${{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} = 0.107}$. The maximum elongation of the flux surface cross-sections (taken at $\phi =\textrm {const}.$), as defined by (5.4) is $e_{\mathrm {max}}=4.4$.

Figure 7. Intensity of the magnetic field on the plasma boundary for a two field period, $A=10$ configuration: (a) side view and (b) top view.

Figure 8. (a) Rotational transform, ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$, and (b) effective ripple, $\epsilon _{\mathrm {eff}}$, for an $N=2$, $A=10$ configuration.

The effective ripple, $\epsilon _{\mathrm {eff}}$, is a simple and convenient parameter that characterises low-collisionality neoclassical transport of electrons (Beidler et al. Reference Beidler, Allmaier, Isaev, Kasilov, Kernbichler, Leitold, Maaßberg, Mikkelsen, Murakami and Schmidt2011). We calculate it using the procedure described by Drevlak et al. (Reference Drevlak, Heyn, Kalyuzhnyj, Kasilov, Kernbichler, Monticello, Nemov, Nührenberg and Reiman2003) in 16 radial points, and find an $\epsilon _{\mathrm {eff}}$ below 1 % up to mid-radius and lower than 2 % everywhere in the plasma volume, see figure 8.

The boundary shape was generated for different aspect ratios, i.e. different distances from the axis. A comparison between the contours of constant magnetic field strength obtained from the NAE (2.1) and those obtained from the VMEC calculation and transformed to boozer coordinates using BOOZ_XFORM (Sanchez et al. Reference Sanchez, Hirshman, Ware, Berry and Spong2000) is shown in figure 9. The difference between both results decreases with increasing the aspect ratio, as expected because the expansion is performed in the distance from the axis. For the largest aspect ratio, 160, the difference is almost imperceptible. The root-mean-squared difference between both magnetic fields is calculated for each aspect ratio and the results are shown in figure 10. The scaling with the aspect ratio is as expected from a first-order expansion, proportional to $1/A^2$.

Figure 9. Contours of the magnetic field intensity on the plasma boundary, $B$, are shown for increasing values of aspect ratio, starting at $A=20$ (a,b). Dotted lines are the contours obtained from the construction and solid lines from VMEC and BOOZ_XFORM. Contours on (b,d, f) correspond to the first-order correction to $B$ and (a,c,e) to the total magnetic field.

Figure 10. Root-mean-squared difference between the magnetic field from the construction and that calculated by VMEC, for different aspect ratios. The dashed line shows the expected scaling $\propto 1/A^{2}$.

7 A family of constant torsion QI solutions

As illustrated by the previous example, we are capable of directly constructing QI equilibria with low electron neoclassical transport, as measured by $\epsilon _{\mathrm {eff}}$, at aspect ratios comparable to existing devices. Performing an optimisation procedure in the space of QI solutions described by the NAE has led to the discovery of $N=1$ configurations with excellent confinement properties as shown in Jorge et al. (Reference Jorge, Plunk, Drevlak, Landreman, Lobsien, Camacho and Helander2022). However, finding configurations with more than one field period and good confinement has proved challenging, even when using optimisation procedures. In order to explore the role of the axis shape in causing this behaviour, we choose a family of axis shapes with $N=2,3,4$, constant torsion and equal per-field-period maximum curvature. Constant torsion was chosen for simplicity in order to obtain smooth solutions for $\sigma$ from (2.6).

The axis shapes chosen are closed curves with minimal bending energy and constant torsion as described in Pfefferlé et al. (Reference Pfefferlé, Gunderson, Hudson and Noakes2018). Their curvature and torsion are

(7.1)\begin{gather} \tau = 4 X(p)K(p) N/L, \end{gather}

(7.2)\begin{gather}\kappa^{2}(s)=\kappa_{0}^{2}\left[ 1 - \mathrm{sn}^{2}\left(\frac{\kappa_{0}}{2p}(s-0.5L),p\right)\right], \end{gather}

where $\mathrm {sn}(a,b)$ is the Jacobi elliptic sine function, $L$ the curve length, $N$ the number of field periods, $p$ a parameter yet to be chosen and $\kappa _{0}$ the maximum curvature. Here $K(p)$ is the complete elliptical integral of first kind and $X^{2}(p)= 2E(p)/K(p) -1$, where $E(p) = E({\rm \pi} /2,p)$ is the incomplete elliptic integral of the second kind.

The maximum curvature is chosen as $\kappa _0=1.6$ for all cases. Given a number of field periods $N$, the parameter $p$ is scanned until a closed curve is found. We find $p_{2}=0.4527$, $p_{3}=0.5927$ and $p_{4}=0.6580$, for two, three and four field periods, respectively. The curvature per field period has the same maximum and similar toroidal dependence for all three cases and has zeros at multiples of ${\rm \pi} /N$, as required by the construction. The necessary torsion for a closed curve increases with $N$ as seen in figure 11. The Frenet–Serret formulae (5.9) are then used to find the curve described by these values of $\kappa$ and $\tau$, as described in § 5. For the three curves, $m=1$.

Figure 11. (a) Curvature, $\kappa$, and torsion, $\tau$, with respect to the toroidal angle $\phi$ of the magnetic axes used in the construction of the different $N$ field-period solutions. Here $\kappa (\varphi )$ and $\tau (\varphi )$ are given by (7.2) and (7.1). The maximum curvature, $\kappa _0 = 1.6$, is the same for the three cases. (b) The $\sigma$ profiles obtained by solving (2.6) for each of the curves described by $\kappa$ and $\tau$ on the left.

Using these curves as axis shapes, three configurations were constructed following the method described in § 6. All per-period parameters and functions used in the construction are kept the same as in § 6; $d(\varphi ) = 1.03\kappa ^{s}$, the parameter $k$ entering (4.6) was set as $k=2$, and the magnetic field on axis as

(7.3)\begin{equation} B_{0} = 1 + 0.15\cos(N\varphi). \end{equation}

The main differences between these configurations is the value of the torsion $\tau$ and the number of field periods. As we are solving (2.6) for $\sigma$ in a single period, we can compare the per-field-period solutions with respect to a scaled angular variable $\varphi /N$, as seen in figure 11(b). The solution for $\sigma$ in one field period for each of these cases is very similar and only deviates in those regions where the curvature values differ, as expected. Consequently, the values found for the per-period rotational transform, ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _{N}={\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} /N$, are also very similar: ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _{\mathrm {2}}= 0.4565$, ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _{\mathrm {3}}= 0.4674$ and ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _{\mathrm {4}}= 0.4707$.

The resulting boundary shape is used to find the magnetic field strength on the boundary using the MHD equilibrium code VMEC and BOOZ_XFORM. The result is compared with the magnetic field from the construction (2.1), and shown for $A=40$ in figure 12, where it is clear that the approximation deteriorates with increasing number of field periods.

Figure 12. Contours of the magnetic field intensity on the plasma boundary, $B$, are shown for increasing values of $N$, the number of field periods, starting at (a) $N=2$, (b) $N=3$ and (c) $N=4$. Dotted lines show the contours obtained from the construction and solid lines those from the calculation with VMEC and BOOZ_XFORM. Aspect ratio was set to $A=40$ for all cases.

In the same manner as in the previous section, we calculate the root-mean-square difference between the intensity of $B(\theta,\varphi )$ in the boundary as calculated by VMEC and as expected from the construction. This calculation is done for each configuration at different aspect ratios, and is shown in figure 13. The scaling for all cases is proportional to $1/A^2$, as expected, but the magnitude of the difference increases with the number of field periods, indicating a deterioration of the approximation with increasing $N$. The same behaviour is observed in the effective ripple, getting significantly worse for the case with four field periods and being optimal, below $\epsilon _{\mathrm {eff}}= 8\times 10^{-3}$, for two field periods, as is evident from figure 14.

Figure 13. Root-mean-squared difference between the magnetic field from the construction and that calculated by VMEC for configurations with $N=2,3,4$ and different values of $A$. The dashed lines show the expected scaling $\propto 1/A^{2}$.

Figure 14. Effective ripple, $\epsilon _{\mathrm {eff}}$, profiles for configurations with $N=2,3,4$ field periods.

The only apparent significant differences between the solutions constructed in this section are the number of field periods and the value of the torsion. Stellarator designs constructed through conventional optimisation have $N = 4,5,6$ and larger, including W7-X, a QI design (very approximately) with 5 field periods. This indicates there should be no fundamental obstacle to obtaining a good approximation to QI fields at larger values of $N$. Thus, with the hope to find such fields with the near-axis framework, we are further motivated to explore magnetic axes with low torsion.

8 A three-field-period, low-torsion-axis example

Following the discussion at the beginning of § 5, and also the results of the previous section, we are motivated to consider axis curves with low torsion to both control elongation and better approximation to QI at larger $N$. Using the procedure described in § 5, we therefore investigate closed curves with $N=3$ and zeros of first order in the curvature at extrema of $B(\varphi )$, and zeros of torsion at second order at $\varphi _{\mathrm {min}}^{i}$. We choose one axis shape from this class that fulfils

(8.1)\begin{equation} \mathrm{max_{\varphi}}(\tau) < \mathrm{max_{\varphi}}(\kappa) \end{equation}

for all toroidal points. The Fourier coefficients describing this axis are

(8.2)\begin{align} R & = 1 + 9.07\times10^{{-}2} \cos{(3\phi)} - 2.06\times10^{{-}2} \cos{(6\phi)}\nonumber\\ & \quad - 1.11\times10^{{-}2} \cos{(9\phi)} - 1.64\times10^{{-}3} \cos{(12\phi)}, \end{align}

(8.3)\begin{align} z & = 0.36 \sin{(3\phi)} + 2.0\times10^{{-}2}\sin{(6\phi)}+ 1.0\times10^{{-}2}\sin{(9\phi)}. \end{align}

Its curvature and torsion are shown as functions of the toroidal angle $\phi$ in figure 15, and $m=0$.

Figure 15. Curvature $\kappa$ (solid), torsion $\tau$ (dashed) and magnetic field on-axis $B_0$ (dotted) with respect to the toroidal angle $\phi$ of the magnetic axis used in the construction of an $N=3$ equilibria. Curvature has zeros of first order at multiples of ${\rm \pi} /3$, corresponding to extrema of $B_0$. Torsion has a second-order zero at minima of $B_0$, $\varphi _{\mathrm {min}}^{i}$.

Using this axis, an $N=3$, QI boundary was constructed for $A=20$, with a magnetic field on-axis given by

(8.4)\begin{equation} B_{0} = 1 + 0.25\cos{(3\varphi)}, \end{equation}

$d(\varphi ) = 1.03\kappa ^{s}$ and the parameter $k$ entering equation (4.6) set to $2$. The resulting plasma boundary and the magnetic field strength on it as found by VMEC are shown in figure 16, where we can observe straight sections around $\varphi _{\mathrm {min}}$, as a consequence of the axis choice. Figure 17 shows the contours of $|B|$ as calculated by VMEC and from the near-axis construction in Boozer coordinates. At this aspect ratio most of the contours still close poloidally, as necessary for quasi-isodynamicity. The magnetic field contours are less straight around the point of maximum $B_{0}$, which is expected because this is the region where $\alpha (\varphi )$ deviates from a perfectly omnigenous form.

Figure 16. Intensity of the magnetic field on the plasma boundary for a three field period, $A=20$ configuration: (a) side view and (b) top view.

Figure 17. Contours of the magnetic field intensity on the boundary, for $N=3$, $A=20$, as calculated by VMEC (solid lines) and from the near-axis construction (dotted lines).

Zero torsion around the bottom of the magnetic well was chosen, together with $d(\varphi )$, as an attempt to reduce the elongation of the plasma boundary, as described in § 5. In figure 18 we see the cross-sections of the plasma boundary for different values of toroidal angle on the left and the evolution of elongation with the toroidal angle. The maximum elongation for this configuration is $e_\mathrm {max}=3.2$, and includes regions where the elongation is nearly $e=1$, the case of a circular cross-section. We also observe that the elongation remains low around the region where torsion vanishes. These facts seem to validate our approach for controlling elongation, and demonstrate that high elongation is not a necessary sacrifice for achieving good omnigenity in QI stellarators, as might be feared from previous results, for example W7-X (Grieger et al. Reference Grieger, Lotz, Merkel, Nührenberg, Sapper, Strumberger, Wobig, Burhenn, Erckmann and Gasparino1992) with a maximum elongation around $e_{\mathrm {max}}\gtrsim 4.5$, or the NAE configuration described by Jorge et al. (Reference Jorge, Plunk, Drevlak, Landreman, Lobsien, Camacho and Helander2022) and QIPC (Subbotin et al. Reference Subbotin, Mikhailov, Shafranov, Isaev, Nührenberg, Nührenberg, Zille, Nemov, Kasilov and Kalyuzhnyj2006), both with maximum elongations higher than $e_{\mathrm {max}}\gtrsim 5$.

Figure 18. (a) Cross sections of the plasma boundary for different values of the toroidal angle $\phi$, and (b) elongation with respect to the cylindrical toroidal angle $\phi$.

The rotational transform profile obtained with VMEC is shown in figure 19(a) and coincides with the value calculated numerically from (2.6), ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} = 0.375$. The effective ripple, $\epsilon _{\mathrm {eff}}$, was found to be lower than 1 % throughout the plasma volume, another indication of the closeness of the solution to omnigenity. The value at different distances from the axis is shown in figure 19(b).

Figure 19. (a) Rotational transform, ${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$, and (b) effective ripple $\epsilon _{\mathrm {eff}}$