Direct construction of stellarator-symmetric quasi-isodynamic magnetic configurations

Abstract

We develop the formalism of the first-order near-axis expansion of the magnetohydrodynamic equilibrium equations described by Garren & Boozer (Phys. Fluids B, vol. 3, issue 10, 1991, pp. 2805–2821) and Plunk et al. (J. Plasma Phys., vol. 85, issue 6, 2019; J. Plasma Phys., vol. 87, issue 6, 2021) for the case of a quasi-isodynamic, $N$-field-period, stellarator-symmetric, single-well magnetic field equilibrium. The importance of the magnetic axis shape is investigated, and we conclude that control of the curvature and torsion is crucial to obtain omnigenous configurations with finite aspect ratio and low effective ripple, especially for a higher number of field periods. For this reason a method is derived to construct classes of axis shapes with favourable curvature and torsion. Solutions are presented, including a three-field-period configuration constructed at an aspect ratio of $A=20$, with a maximum elongation of $e=3.2$ and an effective ripple under $1\,\%$, which demonstrates that high elongation is not a necessary feature of quasi-isodynamic stellarators.

1 Introduction

In recent years, the W7-X experiment has demonstrated that intricately optimised stellarators can be successfully built and operated (Pedersen et al. 2018; Beidler et al. 2021). Stellarators have attractive qualities, such as little net toroidal current and capacity for steady-state operation, that make them suitable for reactors. However, unlike in tokamaks, confinement is not inherently good due to the three-dimensional geometry of the magnetic field. Stellarator magnetic fields need to be designed carefully to ensure neoclassical transport is sufficiently low. A sufficient condition for good orbit confinement in a stellarator is that of omnigenity,

(1.1)\begin{equation} \int (\boldsymbol{v}_{d}\boldsymbol{\cdot} \boldsymbol{\nabla} \psi) \,{\rm d}t = 0, \end{equation}

where $\boldsymbol {v}_{d}$ is the total drift, and the integration is done over the bounce time of a trapped particle. Magnetohydrodynamic (MHD) equilibria for which trapped particle motion fulfils (1.1) have collisionless orbits that are radially confined. A subset of omnigenous magnetic fields are those that satisfy quasi-symmetry, in which the strength of the magnetic field is symmetric (axially, helically or poloidally) in magnetic coordinates. An example of omnigenous fields that are not quasi-symmetric are those called quasi-isodynamic (QI), which have poloidally closed contours of the magnetic field strength, as well as vanishing bootstrap currents at low collisionality (Helander & Nührenberg 2009; Helander 2014). The present work concerns this class of stellarators.

Identification of good configurations, e.g. those with low neoclassical transport and easily buildable coils, is traditionally done through a two-step optimisation procedure. In the first step, a plasma boundary is deformed and the physical properties of the resulting equilibrium are assessed using various codes until the desired plasma properties have been found. Then, in a second optimisation process, coils that reproduce the desired plasma boundary are sought.¹ The equilibrium optimisation procedure requires the use of computationally expensive codes in each step and an initial plasma boundary to be prescribed (Henneberg et al. 2021b). Such optimisation procedures generally identify local minima (Henneberg, Helander & Drevlak 2021a), and it is therefore possible that other, lower, minima may be found elsewhere in the parameter space. In addition, QI optimisation often results in highly shaped plasma boundaries. In order to generate such plasma shapes, complicated coils may be necessary, which are difficult and expensive to build. It is not clear whether this complexity is inherent to particular QI equilibria or if it may be overcome by a more exhaustive or differently initialised search of parameter space. The near-axis approach, pursued here, has the potential to do just this.

This method, first introduced by Garren & Boozer (1991) and revitalised by Landreman & Sengupta (2018) and Landreman, Sengupta & Plunk (2019), allows the construction of omnigenous solutions of the MHD equations at first order in the distance from the axis, using Boozer coordinates. Starting with an axis shape and a set of functions of the toroidal angle, this near-axis expansion (NAE) allows the systematic and efficient construction of plasma boundaries corresponding to omnigenous equilibria. Different procedures apply for the case of QI and for quasi-symmetry, with the latter described in Landreman et al. (2019) and used successfully combined with an optimisation procedure in Giuliani et al. (2022). The case of quasi-isodynamicity was first discussed by Plunk, Landreman & Helander (2019) and is described in greater detail in § 2.

In § 3 we derive all the constraints on the input functions for the case of a stellarator-symmetric, single-magnetic-well² solution with multiple field periods. We then proceed to find expressions for the geometric functions required as input for the near-axis construction and show that, for this particular case, the number of free parameters is reduced.

Exactly omnigenous fields are necessarily non-analytic as shown by Cary & Shasharina (1997). Hence, omnigenity must necessarily be broken to achieve analytical solutions. In the NAE we do this in a controlled way, through the careful definition of one of the geometrical input functions, $\alpha (\varphi )$. A new way to express this function, more smoothly than done previously by Plunk et al. (2019), is proposed in § 4.

Very recently, single-field-period QI configurations with excellent confinement properties have been found by optimising within the space of these NAE solutions in Jorge et al. (2022). Configurations with more field periods have proved challenging to find, even when using the optimiser.

To investigate the causes of this difficulty, we analyse possible choices of the magnetic axis shape. In contrast to traditional optimisation, where a plasma boundary is the starting point and a magnetic axis has to be found as part of the equilibrium calculation, the NAE requires the magnetic axis shape as input. In § 5 we discuss the freedom in the shape of the axis and its mathematical description. We shed some light on the difficulties associated with finding closed curves compatible with the NAE, specifically the need to find closed curves with points of zero curvature at prescribed locations, as well as the role curvature and torsion play on the quality of the approximation. We then describe a method for generating stellarator-symmetric axis shapes with points of zero curvature and torsion at prescribed orders.

A common problem encountered when constructing NAE solutions is the tendency that the boundary cross-section becomes highly elongated. Indications that increasing elongation of the plasma boundary leads to increasingly complex coils have been found by Hudson et al. (2018). Another argument for reducing elongation is the fact that very elongated cross-sections result in small plasma volume, as seen from the ellipse area dependence on the elongation $e$, when fixing the major axis $a$, $A={\rm \pi} a b= {\rm \pi}a^2/e$. A reduction of plasma volume results in increasing transport losses and construction costs per unit plasma volume. The relation between elongation and the relevant NAE parameters is discussed in § 5, in particular the effect of torsion on shaping.

In § 6 we show how the NAE method can be used to construct a two-field-period solution that closely approximates omnigenity, and we discuss its geometric properties and neoclassical transport as measured by the effective ripple $\epsilon _{\mathrm {eff}}$. The particular case of a family of solutions with axes of constant torsion and increasing number of field periods is analysed in § 7 to show the effect of torsion on the quality of the near-axis approximation. Finally, in § 8, a three-field-period equilibrium is constructed around an axis shape specially designed to have low torsion in order to illustrate the importance of the axis shape in controlling transport and elongation.

2 Near-axis expansion

The NAE, as described by Garren & Boozer (1991), solves the equilibrium MHD equations by performing a first-order Taylor expansion about the magnetic axis in pseudo-Cartesian coordinates (Kuo-Petravic & Boozer 1987). 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 (2020).

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. (2019) and Plunk et al. (2021), 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. (2019), 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 (1991), 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 (1997) was used in Plunk et al. (2019). 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. (2019) 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. (2019) 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 (1998), 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.³ 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. (2019), 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 1983).

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 (2018) 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 2009; Helander, Geiger & Maaßberg 2011). 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 2018). 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. (2019)

(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. (2019) 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. (2022).

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. (2022).

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. (2019)

(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. 2011). We calculate it using the procedure described by Drevlak et al. (2003) 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. 2000) 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. (2022). 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. (2018). 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. 1992) with a maximum elongation around $e_{\mathrm {max}}\gtrsim 4.5$, or the NAE configuration described by Jorge et al. (2022) and QIPC (Subbotin et al. 2006), 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}}$

9 Conclusion

We have described in detail the NAE method to construct QI, stellarator symmetric, single-magnetic-well equilibria with $N\geqslant 1$ field periods, valid at first order in the distance from the magnetic axis. A new way of achieving better continuity and smoothness of these configurations, as compared with the previous method of Plunk et al. (2019), has been introduced and used to construct equilibria with $N>1$.

The problem of finding axis shapes compatible with the NAE is discussed, in particular the order of zeros in axis curvature, and the naturally arising increase of torsion for increasing number of field periods, which we argue is an underlying reason for the deterioration of the approximation at finite aspect ratio. A method to systematically describe and construct closed curves with zeros of curvature and torsion, at different orders, at toroidal locations of extrema of the magnetic field strength has also been presented.

We have demonstrated the validity of the NAE method, with a two-field-period example, by showing that the difference between the magnetic field as calculated by the NAE and that obtained using the equilibrium code VMEC falls with increasing aspect ratio, and scales as $1/A^2$, as expected from a first-order expansion. We describe the construction of this two-field-period configuration, and find that it has good confinement, as shown by $\epsilon _{\mathrm {eff}}<2\,\%$ throughout the plasma volume at aspect ratio $A=10$.

We have also constructed a family of solutions, for $N=2,3,4$, having axes with constant torsion, and very similar per-field-period initial parameters. The approximate omnigenity of these solutions deteriorates if the number of field periods is increased. This example shows the importance of reducing the maximum torsion of the axis to achieve equilibria close to omnigenity at low aspect ratio.

In the last section we have demonstrated how the choice of an axis with zero torsion around the point of minimum magnetic field strength and constrained maximum torsion enables us to find a three-field-period configuration with low elongation and small neoclassical transport. The effective ripple remains at under $1\,\%$ for an aspect ratio $A=20$, and the maximum elongation $e_{\mathrm {max}}=3.2$, demonstrating that low elongation is achievable in QI stellarators.

These configurations demonstrate that the NAE method can be used to construct magnetic equilibria with multiple field periods that maintain good confinement properties at low aspect ratios. We emphasise that these examples were all obtained without the need for numerically costly optimisation procedures.

Given the importance of the axis shape in the quality of the equilibrium, a natural next step is to reintroduce an element of optimisation, as done for the one-field-period example of Jorge et al. (2022), but appropriately restricting the search space. Specifically, we can define an optimisation space according to classes of axis curves satisfying prescribed conditions on the torsion and curvature, and search this space for configurations with attractive properties. Another interesting analysis, which can be performed thanks to the speedy calculation of solutions enabled by the NAE, is to systematically and exhaustively map the space of QI solutions and its dependence on the input functions used for the construction. Such exploration might allow physical insight to be gained into the structure of the solution space and help explain certain difficulties associated with traditional optimisation techniques.

In this work, as well as in Plunk et al. (2019), the focus has been on introducing deviations to quasi-isodynamicity near the field maxima to target deeply trapped particles. Current work is being done on breaking omnigenity at other toroidal regions in order to improve confinement of different classes of trapped particles.

The only shaping of the plasma boundary that enters at first order in the NAE is the elongation of the elliptical cross-sections. Using solutions generated at higher order, together with traditional optimisation procedures also seems a promising way for obtaining configurations with better confinement properties and stronger shaping.