Direct construction of optimized stellarator shapes. Part 3. Omnigenity near the magnetic axis

Abstract

The condition of omnigenity is investigated, and applied to the near-axis expansion of Garren & Boozer (Phys. Fluids B, vol. 3 (10), 1991a, pp. 2805–2821). Due in part to the particular analyticity requirements of the near-axis expansion, we find that, excluding quasi-symmetric solutions, only one type of omnigenity, namely quasi-isodynamicity, can be satisfied at first order in the distance from the magnetic axis. Our construction provides a parameterization of the space of such solutions, and the cylindrical reformulation and numerical method of Landreman & Sengupta (J. Plasma Phys., vol. 84 (6), 2018, 905840616); Landreman et al. (J. Plasma Phys., vol. 85 (1), 2019, 905850103), enables their efficient numerical construction.

1 Introduction

The ‘near-axis expansion’, an expansion in the distance from the magnetic axis, has been investigated as a method to characterize the space of magnetic equilibria, and determine how well the condition of quasi-symmetry can be maintained across the plasma volume (Garren & Boozer 1991a , b ). More recently, it is being reconsidered as a method to rapidly generate such equilibria, without resorting to costly numerical optimization.

The present work, the third paper in a series (Landreman & Sengupta 2018; Landreman, Sengupta & Plunk 2019), is concerned with extending these ideas to the broader class of omnigenous fields, which confine collisionless particle orbits, but do not necessarily satisfy quasi-symmetry. By relaxing this requirement, a large amount of freedom is gained, and the solution space can be considered as an exhaustive catalogue of ‘good’ equilibria (i.e. all those that confine trapped particle orbits). We note that some analytical work of this type provided a basis for the development of stellarator configurations (Gori, Lotz & Nührenberg 1996; Nührenberg 2010), but was limited to treatment of deeply trapped particles in a particular (quasi-isodynamic) class of equilibria.

Two related properties, characterizing optimal stellarator equilibria, are quasi-symmetry and omnigenity. The relationship between them is somewhat paradoxical, as found by Cary & Shasharina (1997): quasi-symmetry is formally a subclass of omnigenity – all quasi-symmetric fields are also omnigenous. However, those omnigenous fields that break quasi-symmetry are, under reasonable conditions, also non-analytic, so that under the restriction of analyticity, it appears, disappointingly, that these classes are the same. However, by replacing a non-analytic solution with a sufficiently smooth approximation, one can indeed approach omnigenity while sharply violating quasi-symmetry, explaining why it is indeed fruitful to simply search for equilibria that confine particles, without necessarily imposing the satisfaction of quasi-symmetry.

These facts set the stage for our near-axis expansion, where analyticity is indeed assumed in the two dimensions perpendicular to the magnetic axis. It should be no surprise, therefore, that we find two variants of omnigenity (those with toroidally and helically closed contours of magnetic field strength) cannot be achieved without breaking quasi-symmetry – one, however, cannot rule out their existence entirely, as they might be found via a less restrictive construction.

We do, however, find that omnigenity can be approximately satisfied (at first order in the distance from the magnetic axis) for the class of equilibria with poloidally closed contours of magnetic field strength, so-called quasi-isodynamic equilibria. Constructing such equilibria presents some technical difficulties (as compared with quasi-symmetric equilibria), but the potential benefit in terms of freedom is quite substantial – three additional free functions of the toroidal angle, as compared with quasi-symmetric construction!

The near-axis construction of optimal stellarator equilibria can thus be summarized as follows. Quasi-axisymmetric and quasi-helically symmetric equilibria can be directly constructed, but there are no other consistent omnigenous equilibria with toroidally or helically closed contours of field strength. On the other hand, one cannot construct quasi-poloidally symmetric equilibria near the magnetic axis, but we find that one can indeed find approximately quasi-isodynamic ones, i.e. quasi-symmetry must be broken to achieve this kind of omnigenity. Therefore, it seems that quasi-helical symmetry, quasi-axial symmetry and quasi-isodynamicity are the three natural ways to achieve good stellarator equilibria with a near-axis approximation. These are the only types of equilibria with collisionless orbit confinement and therefore small neoclassical transport. It should be noted, however, that quasi-isodynamic equilibria do not share other neoclassical properties of quasi-symmetric ones. The bootstrap current vanishes identically (Helander & Nührenberg 2009; Landreman & Catto 2012) and the transport is not intrinsically ambipolar, which severely restricts plasma rotation (Helander 2007; Helander & Simakov 2008).

The paper is organized as follows. In § 2 we review the definition and geometric interpretation of omnigenity, and introduce a coordinate mapping, based on that of Cary & Shasharina (1997). We use the mapping as a tool to construct omnigenous magnetic-field-strength functions. In particular, in § 5, the mapping is expanded for nearly quasi-symmetric fields, which, due to its invertibility in this case, means that such fields can be parameterized by a free function representing the perturbation. Although this expansion may have more general application, it is particularly useful for generating equilibria using the near-axis expansion, since such equilibria must necessarily be nearly quasi-symmetric. This is discussed in § 6, where the near-axis and omnigenous forms of the magnetic field strength are checked for consistency. In this section two classes of omnigenity are shown to yield inconsistencies, and are therefore eliminated from consideration, leaving only the third, quasi-isodynamicity. After demonstrating consistency in the forms of magnetic field strength, the remaining conditions needed to construct an equilibrium solution are then examined in § 7. In § 8 it is demonstrated numerically that the equilibrium equation can indeed be solved, yielding solutions that satisfy omnigenity to the degree predicted by the theory. The results are summarized and discussed in § 9.

2 Omnigenity

The most fundamental optimization for a stellarator is the shaping of the magnetic field to confine collisionless trapped particle orbits. This is achieved when there is no ‘radial’ excursion of such particles during a single orbit, as expressed in terms of the change in the coordinate labelling the magnetic flux surfaces, here $\unicode[STIX]{x1D713}$ ( $2\unicode[STIX]{x03C0}\unicode[STIX]{x1D713}$ is the toroidal magnetic flux). Denoting the magnetic drift velocity of a particle as $\boldsymbol{v}_{d}$ and using $\text{d}\unicode[STIX]{x1D713}/\text{d}t=\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ , the condition is written as $\unicode[STIX]{x0394}\unicode[STIX]{x1D713}=0$ , where, for a given magnetic well, is found to be

(2.1) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x0394}\unicode[STIX]{x1D713} & = & \displaystyle \oint \text{d}t\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713},\end{eqnarray}$$

(2.2) $$\begin{eqnarray}\displaystyle & = & \displaystyle 2\int _{l_{1}}^{l_{2}}\frac{\text{d}l}{v_{\Vert }}\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713},\end{eqnarray}$$

(2.3) $$\begin{eqnarray}\displaystyle & = & \displaystyle \mathop{\sum }_{\unicode[STIX]{x1D6FE}}2\unicode[STIX]{x1D6FE}\int _{B_{\text{min}}}^{1/\unicode[STIX]{x1D706}}\frac{\text{d}B}{v_{\Vert }\unicode[STIX]{x2202}B/\unicode[STIX]{x2202}l}\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713},\end{eqnarray}$$

(2.4) $$\begin{eqnarray}\displaystyle & = & \displaystyle \mathop{\sum }_{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D6FE}\int _{B_{\text{min}}}^{1/\unicode[STIX]{x1D706}}\text{d}B\,h(B,\unicode[STIX]{x1D706})Y(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},B,\unicode[STIX]{x1D6FE}).\end{eqnarray}$$

In the first line, the quantity is represented as a time average. In the second line, it is rewritten as an integral over the arc length $l$ between two bounce points $l_{1}$ and $l_{2}$ . In the third line, the arc length is replaced by the magnetic field strength $B$ as a variable of integration, and the sign $\unicode[STIX]{x1D6FE}\equiv \text{sign}(\unicode[STIX]{x2202}B/\unicode[STIX]{x2202}l)$ is introduced which identifies the two sides of the magnetic well, and the upper bound of integration is written explicitly in terms of the quantity $\unicode[STIX]{x1D706}=v_{\bot }^{2}/(v^{2}B)$ . To obtain the final line, we use the expression for the drift velocity of a particle of mass $m$ and charge $q$ , $\boldsymbol{v}_{d}=m\hat{\boldsymbol{b}}/qB\times (v_{\Vert }^{2}\unicode[STIX]{x1D73F}+(v_{\bot }^{2}/2)\unicode[STIX]{x1D735}\ln B)$ , where $\unicode[STIX]{x1D73F}=\hat{\boldsymbol{b}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\hat{\boldsymbol{b}}$ is magnetic curvature, $\hat{\boldsymbol{b}}=\boldsymbol{B}/B$ , and note that $\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }(\boldsymbol{B}\times \unicode[STIX]{x1D73F})=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }(\boldsymbol{B}\times \unicode[STIX]{x1D735}\ln B)$ . Thus, the integrand may be separated into a purely geometric factor $Y$ (depending generally on the flux surface label $\unicode[STIX]{x1D713}$ , field line label $\unicode[STIX]{x1D6FC}$ and the magnetic field strength $B$ ),

(2.5) $$\begin{eqnarray}Y=\frac{\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B}{\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}B},\end{eqnarray}$$

and an additional factor $h=(2mv/qB)(1-\unicode[STIX]{x1D706}B/2)/\sqrt{1-\unicode[STIX]{x1D706}B}$ depending on the magnetic field strength and the pitch angle. The expression for $\unicode[STIX]{x0394}\unicode[STIX]{x1D713}$ given by (2.4) can be viewed as an integral transformation from $B$ space to $\unicode[STIX]{x1D706}$ space (the velocity $v$ appears only as a multiplicative factor), which can be inverted, implying that the condition $\unicode[STIX]{x0394}\unicode[STIX]{x1D713}=0$ is satisfied if and only if

(2.6) $$\begin{eqnarray}\mathop{\sum }_{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D6FE}Y=0.\end{eqnarray}$$

For a quasi-symmetric field, $Y$ is constant on flux surfaces (Helander & Simakov 2008) and this condition is trivially satisfied. Equation (2.6) is a purely geometric condition, in the sense that it is a condition on the spatial dependence of the magnetic field strength. For the purposes of this paper, we will consider the consequences of this condition on the magnetic field in Boozer poloidal and toroidal coordinates, $\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D711}$ respectively. Using (2.5), we write the magnetic field in the numerator in the covariant form $\boldsymbol{B}_{\text{cov}}=G\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}+I\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+K\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ , and the field in contravariant form in the denominator, , obtaining

(2.7)

where the toroidal and poloidal currents $I(\unicode[STIX]{x1D713})$ and $G(\unicode[STIX]{x1D713})$ have been introduced. This expression can be interpreted as the ratio of two (linearly independent) components of the gradient of $B$ in the $\unicode[STIX]{x1D703}$ – $\unicode[STIX]{x1D711}$ plane. Equation (2.6) states that this ratio must be equal at the two (bounce) points of equal magnetic field strength, i.e.  $\unicode[STIX]{x1D6FE}=\pm 1$ . This fact has the simple geometric interpretation that the gradients of $B$ must be anti-parallel, i.e. the contours of constant magnetic field strength must be parallel at those two points. Equivalently, the angular separation between points of equal magnetic field strength (within a given magnetic well) must be the same for all field lines. A number of corollaries follow intuitively, for instance that the contours of magnetic field strength must all close in the same sense (either poloidally, toroidally or helically), and the contour of the global maximum must be straight (for irrational rotational transform).

3 Cary–Shasharina (CS) mapping

It was shown by Cary & Shasharina (1997) that the condition of omnigenity can be satisfied by functions $B(\unicode[STIX]{x1D709},\unicode[STIX]{x1D701})$ that satisfy a set of geometric conditions. Here, $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D701}$ are angular coordinates (defined in terms of Boozer angles) defined so that contours of constant magnetic field strength close in the $\unicode[STIX]{x1D709}$ direction. We will use the angle $\unicode[STIX]{x1D709}$ to denote either the poloidal or toroidal coordinate, $\unicode[STIX]{x1D703}$ or $\unicode[STIX]{x1D711}$ , depending on the class of symmetry. The idea with the CS mapping is to introduce a new coordinate $\unicode[STIX]{x1D702}$ , to replace the $\unicode[STIX]{x1D701}$ coordinate, so that lines of constant magnetic field strength will appear straight in the $\unicode[STIX]{x1D702}$ – $\unicode[STIX]{x1D709}$ plane.

The expression for magnetic field strength given by Cary & Shasharina (1997) is

(3.1) $$\begin{eqnarray}B=B_{0}(1+\unicode[STIX]{x1D716}_{r}\cos (\unicode[STIX]{x1D702})),\end{eqnarray}$$

where $|\unicode[STIX]{x1D716}_{r}|<1$ ensures that $B$ has no zeros; other expressions would also work, but this one has the benefit that it is analytic in $\unicode[STIX]{x1D702}$ . The coordinate mapping for $\unicode[STIX]{x1D702}$ is defined via

(3.2) $$\begin{eqnarray}\unicode[STIX]{x1D701}=\unicode[STIX]{x1D702}+F(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702}).\end{eqnarray}$$

When expressed in this manner, an explicit form of $F$ can be given, from which the condition of omnigenity is enforced

(3.3)

where we note that   is the rotational transform in the  $\unicode[STIX]{x1D709}$ – $\unicode[STIX]{x1D701}$ plane (number of $\unicode[STIX]{x1D701}$ transits per $\unicode[STIX]{x1D709}$ transit). The contour of maximum field strength is straight, i.e. $\unicode[STIX]{x1D702}=0$ at $\unicode[STIX]{x1D701}=0$ , independent of $\unicode[STIX]{x1D709}$ , so we have the additional constraint

(3.4) $$\begin{eqnarray}f(\unicode[STIX]{x1D709},0)=0.\end{eqnarray}$$

The function $\unicode[STIX]{x0394}\unicode[STIX]{x1D701}(\unicode[STIX]{x1D702})$ specifies the angular separation, as measured along a magnetic field line, between a contour of constant $\unicode[STIX]{x1D702}$ , say $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}_{0}$ , and the corresponding contour where the magnetic field has the same strength, i.e.  $\unicode[STIX]{x1D702}_{1}=2\unicode[STIX]{x03C0}-\unicode[STIX]{x1D702}_{0}$ ; this function is therefore defined for $\unicode[STIX]{x03C0}\leqslant \unicode[STIX]{x1D702}\leqslant 2\unicode[STIX]{x03C0}$ , and has the properties $\unicode[STIX]{x0394}\unicode[STIX]{x1D701}(\unicode[STIX]{x03C0})=0$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D701}(2\unicode[STIX]{x03C0})=2\unicode[STIX]{x03C0}$ . The condition of omnigenity is that this angular separation only depends on the magnetic field strength, in this case labelled by $\unicode[STIX]{x1D702}$ .

The above construction only considers the relatively simple case of a single magnetic well. Parra et al. (2015) demonstrated that more complicated fields are possible, by providing an explicit construction for a two-well case. It is not surprising that more general omnigenous magnetic fields, with arbitrarily complex arrangement of magnetic wells are also allowed, and a general construction is given in the following section.

4 General omnigenous fields

One can think of the CS mapping as representing the general deformation of a one-dimensional prototype magnetic-field-strength function, i.e. (3.1), such that quasi-symmetry is broken while preserving omnigenity. Let us introduce a general function $\bar{B}(\unicode[STIX]{x1D702})$ , and propose a procedure to continuously deform it. This function is assumed to satisfy

(4.1) $$\begin{eqnarray}\bar{B}(\unicode[STIX]{x1D702})>0,\end{eqnarray}$$

be $2\unicode[STIX]{x03C0}$ periodic and reach its global maximum at $\unicode[STIX]{x1D702}=0$ . We would also like this function to possess some degree of smoothness (so that smoothness imparted in the mapping for $\unicode[STIX]{x1D702}$ may be inherited by the magnetic field itself) but it need not be analytic; for instance piecewise polynomials may be convenient.

Figure 1. Example magnetic field with several magnetic wells. The numbered intervals can be identified as the trapping domains, i.e. from left to right we have $D_{1L}$ , $D_{2L}$ , $D_{3L}$ , $D_{3R}$ , $D_{4L}$ , $D_{4R}$ , $D_{2R}$ , $D_{5L}$ , $D_{5R}$ , $D_{1R}$ .

The function $\bar{B}(\unicode[STIX]{x1D702})$ , an example of which is given in figure 1, implies a set of boundary points in $\unicode[STIX]{x1D702}$ , which define what we will call trapping domains. Let us itemize these points: First there are the locations of local minima and maxima, starting with the global maximum at $\unicode[STIX]{x1D702}=0$ . Then, each local maximum that is less than the global maximum implies a number of additional boundary points associated with values of $\unicode[STIX]{x1D702}$ where the magnetic field reaches the level of that local maximum. We will denote these boundary points, in the order they appear along the field line, as $\unicode[STIX]{x1D702}_{0},\unicode[STIX]{x1D702}_{1},\unicode[STIX]{x1D702}_{2}\ldots \,$ , and we will number the trapping domains, defined by neighbouring points, in the order of increasing $\unicode[STIX]{x1D702}$ , but assign the same index to the two domains that have matching bounce points; see figure 1. Thus, the left-hand domain, where $\bar{B}^{\prime }<0$ will be numbered $i$ , and denoted $D_{iL}$ while the complimenting domain, with $\bar{B}^{\prime }>0$ will have the same index and be denoted $D_{iR}$ . We denote the unions of left and right trapping domains as $D_{L}$ and $D_{R}$ .

Next a function $\unicode[STIX]{x1D702}_{b}(\unicode[STIX]{x1D702},i)$ must be constructed, which returns the left bounce point for a given point in a right-hand domain. That is, for any point $\unicode[STIX]{x1D702}$ contained in a right-hand domain $D_{iR}$ , the function $\unicode[STIX]{x1D702}_{b}(\unicode[STIX]{x1D702},i)$ returns the corresponding bounce point in the left-hand part of domain $D_{iL}$ . (The quantity corresponding to $\unicode[STIX]{x1D702}_{b}$ in the original CS mapping is $2\unicode[STIX]{x03C0}-\unicode[STIX]{x1D702}$ .) Note the explicit dependence on the index $i$ , since we must distinguish between domains at locations that exactly coincide with boundary points. For example, for the case depicted in figure 1, $\unicode[STIX]{x1D702}_{b}(\unicode[STIX]{x1D702}_{6},4)=\unicode[STIX]{x1D702}_{4}$ , while $\unicode[STIX]{x1D702}_{b}(\unicode[STIX]{x1D702}_{6},2)=\unicode[STIX]{x1D702}_{2}$ . Note that it may be advantageous to use a simple method for constructing $\bar{B}(\unicode[STIX]{x1D702})$ , so that the function $\unicode[STIX]{x1D702}_{b}$ might be defined analytically.

As before, the mapping is given by (3.2), but now the function $F$ must be defined over the various trapping domains,

(4.2)

where $\unicode[STIX]{x0394}\unicode[STIX]{x1D701}(\unicode[STIX]{x1D702},i)$ prescribes the angular difference between field lines in magnetic coordinate $\unicode[STIX]{x1D701}$ . Note that this can be confirmed by evaluating (3.2) at two bounce points in a well, and taking the difference. That is, take ( $\unicode[STIX]{x1D709}_{1}$ , $\unicode[STIX]{x1D701}_{1}$ , $\unicode[STIX]{x1D702}_{1}$ ) in $D_{iL}$ and ( $\unicode[STIX]{x1D709}_{2}$ , $\unicode[STIX]{x1D701}_{2}$ , $\unicode[STIX]{x1D702}_{2}$ ) in $D_{iR}$ . Since these are bounce points, we have the identities $\unicode[STIX]{x1D702}_{b}(\unicode[STIX]{x1D702}_{2},i)=\unicode[STIX]{x1D702}_{1}$ , and . The difference of the two equations thus gives $\unicode[STIX]{x1D701}_{2}-\unicode[STIX]{x1D701}_{1}=\unicode[STIX]{x0394}\unicode[STIX]{x1D701}(\unicode[STIX]{x1D702}_{2},i)$ .

While (4.2) is a relatively simple and direct generalization of the CS mapping, it is cumbersome to enforce continuity of $\unicode[STIX]{x1D701}$ at the boundary points. In part this is due to the fact that the function $\unicode[STIX]{x0394}\unicode[STIX]{x1D701}$ must be constructed in a consistent manner, for instance so that the angular distance spanning a well is consistent with the distance prescribed in the external trapping regions, resulting in a hierarchy of additional constraints.

Therefore, we propose a simplification, without loss of generality, namely to consider the angular distance as being prescribed by $\bar{B}(\unicode[STIX]{x1D702})$ itself. That is, we take $\unicode[STIX]{x0394}\unicode[STIX]{x1D701}(\unicode[STIX]{x1D702},i)=\unicode[STIX]{x0394}\unicode[STIX]{x1D702}(\unicode[STIX]{x1D702},i)=\unicode[STIX]{x1D702}-\unicode[STIX]{x1D702}_{b}(\unicode[STIX]{x1D702},i)$ . We arrive then at the simple form

(4.3)

Note now that the function $f(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702})$ need only be defined over the left-hand domain $D_{L}$ . We emphasize that there is no loss of generality in this choice of $\unicode[STIX]{x0394}\unicode[STIX]{x1D701}$ , as angular separation between bounce points can be controlled during the construction of $\bar{B}(\unicode[STIX]{x1D702})$ . This choice also has the benefit of separating out the freedom associated with the corresponding narrower class of quasi-symmetric magnetic fields, from the freedom associated with breaking quasi-symmetry.

What remains is merely to enforce continuity of $\unicode[STIX]{x1D701}$ (thus $F$ ) at boundary points (one may also wish to enforce continuity of some number of derivatives). Some of these conditions are trivial: by construction, continuity of $F$ is already satisfied at local minima ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D702}=0$ at such points). Also, at boundary points between two left-hand domains, $F$ will be continuous as long as $f$ is continuous. The other boundary points must satisfy shifted conditions due to the translation by  in the $\unicode[STIX]{x1D709}$ coordinate, but this can clearly be satisfied since the function $f$ is completely arbitrary, aside from the requirement that $f(\unicode[STIX]{x1D709},0)=0$ and the invertibility requirement $\unicode[STIX]{x2202}F/\unicode[STIX]{x2202}\unicode[STIX]{x1D702}>-1$ (i.e. that the Jacobian of the variable transformation is non-zero).

Thus it is demonstrated that arbitrarily complex magnetic well structures are compatible with omnigenous $B$ , and the construction of such functions is only constrained by a set of matching conditions, and an invertibility requirement. In the remainder of the paper, the definition for $F$ given in (4.3) will be used.

4.1 Analyticity

Analyticity of the magnetic field $B$ in its spatial coordinates, will, under reasonable conditions, imply analyticity of the function $F$ in the coordinate $\unicode[STIX]{x1D702}$ (Krantz & Parks 2002). Using their mapping, Cary & Shasharina (1997) showed that (assuming irrational rotational transform) analyticity of $F$ in $\unicode[STIX]{x1D702}$ can only be satisfied for quasi-symmetric fields, in particular continuity of derivatives of the magnetic field at $\unicode[STIX]{x1D702}=0,2\unicode[STIX]{x03C0}$ , where the field maximum is located, require that the corresponding derivatives of $f$ are zero,

(4.4) $$\begin{eqnarray}\left.\frac{\text{d}^{n}f}{\text{d}\unicode[STIX]{x1D702}^{n}}\right|_{\unicode[STIX]{x1D702}=0}=0.\end{eqnarray}$$

Thus if $f$ is to be analytic, it must, generally speaking, be independent of $\unicode[STIX]{x1D702}$ , i.e. it must be quasi-symmetric. The interesting cases, therefore, are non-analytic. Note that other boundary points, for instance local minima of the magnetic field (say $\unicode[STIX]{x1D702}_{\text{min}}=\unicode[STIX]{x1D702}_{j}$ ) are also locations of potential discontinuity in derivatives. Therefore, it may be useful to establish the analyticity constraints there as well. For example, by imposing $(\unicode[STIX]{x2202}F/\unicode[STIX]{x2202}\unicode[STIX]{x1D702})_{\unicode[STIX]{x1D702}_{\text{min}}-}=(\unicode[STIX]{x2202}F/\unicode[STIX]{x2202}\unicode[STIX]{x1D702})_{\unicode[STIX]{x1D702}_{\text{min}}+}$ , using (4.3), we arrive at the condition

(4.5)

Thus it would seem that analyticity at $\unicode[STIX]{x1D702}_{\text{min}}$ demands a particular relationship between $\unicode[STIX]{x0394}\unicode[STIX]{x1D702}$ and $f$ at boundary points. However, given that analyticity must already be abandoned at $\unicode[STIX]{x1D702}=0$ , it is not obvious why one should demand analyticity at other locations such as $\unicode[STIX]{x1D702}_{\text{min}}$ . Nevertheless, conditions such as (4.4) and (4.5), can be used to enforce continuity of arbitrary derivatives of the magnetic field strength, if desired.

4.2 Invertibility

It would seem that the mapping (4.3) allows one to generate a generic omnigenous field-strength-function using the free functions $f$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D701}$ . However, a complication arises in that its invertibility is not guaranteed. That is, although $\unicode[STIX]{x1D701}$ is uniquely defined in the $\unicode[STIX]{x1D709}$ – $\unicode[STIX]{x1D702}$ plane, it is not necessarily true that a unique value of $\unicode[STIX]{x1D702}$ can be obtained as a function of $\unicode[STIX]{x1D709}$ and $\unicode[STIX]{x1D701}$ . Note, for instance, that two contours of two values of $\unicode[STIX]{x1D702}$ could intersect, which cannot be allowed since the magnetic field strength must be single valued.

The issue is reminiscent of the problem involving the specification of a general coordinate mapping $\boldsymbol{x}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$ , which can generate self-intersecting magnetic surfaces. Practically, this presents a challenge for the parametrization of the space of omnigenous magnetic field strengths, i.e. devising a general method to construct omnigenous functions $B(\unicode[STIX]{x1D709},\unicode[STIX]{x1D701})$ . However, since the change of coordinates is fairly simple (it essentially one-dimensional, $\unicode[STIX]{x1D701}\rightarrow \unicode[STIX]{x1D702}$ ), the Jacobian is readily computed, and the condition that it is non-zero could be asserted, i.e.

(4.6) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D701}}{\text{d}\unicode[STIX]{x1D702}}=1+\frac{\text{d}F}{\text{d}\unicode[STIX]{x1D702}}>0,\quad \text{for all }\unicode[STIX]{x1D709}\text{ and }\unicode[STIX]{x1D702},\end{eqnarray}$$

which translates, via (4.3), into a condition on $f$ ; note that requiring the derivatives of $f$ to be small would suffice. This might be suitable for numerically generating well-behaved mappings. However, for present purposes, we find a more convenient approach is to expand about known omnigenous fields, the quasi-symmetric ones, relying on the fact that sufficiently small deviations will preserve the topology of the lines of constant $\unicode[STIX]{x1D702}$ , and therefore ensure the invertibility of the mapping.

5 Omnigenity near quasi-symmetry

In this section we derive an approximate form of omnigenity, expanding about a quasi-symmetric field. Besides the near-axis expansion, this form could have other applications, such as numerical optimization or direct solutions applying omnigenity on a single surface (Plunk & Helander 2018). Let us express the mapping as being composed of a quasi-symmetric and non-quasi-symmetric parts,

(5.1) $$\begin{eqnarray}\displaystyle & \displaystyle F=F_{0}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D716}F_{1}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702}), & \displaystyle\end{eqnarray}$$

(5.2) $$\begin{eqnarray}\displaystyle & \displaystyle f=f_{0}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D716}f_{1}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702}). & \displaystyle\end{eqnarray}$$

The mapping is then defined by

(5.3) $$\begin{eqnarray}\unicode[STIX]{x1D701}=\unicode[STIX]{x1D702}+F_{0}(\unicode[STIX]{x1D702})+\unicode[STIX]{x1D716}F_{1}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D702}).\end{eqnarray}$$

We can now consider the deviation to be small, $\unicode[STIX]{x1D716}\ll 1$ , so that the coordinate $\unicode[STIX]{x1D702}$ may be expanded as

(5.4) $$\begin{eqnarray}\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}_{0}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D702}_{1}.\end{eqnarray}$$

The mapping at zeroth order is then

(5.5) $$\begin{eqnarray}\unicode[STIX]{x1D701}=\unicode[STIX]{x1D702}_{0}+F_{0}(\unicode[STIX]{x1D702}_{0}),\end{eqnarray}$$

where

(5.6) $$\begin{eqnarray}F_{0}(\unicode[STIX]{x1D702})=\left\{\begin{array}{@{}l@{}}f_{0}(\unicode[STIX]{x1D702}),\quad \text{for }\unicode[STIX]{x1D702}\in D_{iL}\quad \\ f_{0}(\unicode[STIX]{x1D702}_{b}(\unicode[STIX]{x1D702})),\quad \text{for }\unicode[STIX]{x1D702}\in D_{iR}.\quad \end{array}\right.\end{eqnarray}$$

Henceforth, and without loss of generality, we will assume $f_{0}=0$ , so that $F_{0}=0$ and $\unicode[STIX]{x1D702}_{0}=\unicode[STIX]{x1D701}$ . This amounts to absorbing all the freedom of the zeroth-order mapping into the definition of the zeroth-order function $\bar{B}_{0}(\unicode[STIX]{x1D702}_{0})=\bar{B}_{0}(\unicode[STIX]{x1D701})$ . Note that $\bar{B}_{0}$ also determines the function $\unicode[STIX]{x1D702}_{b}$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D702}$ . The first-order terms of (5.3) yield

(5.7) $$\begin{eqnarray}\unicode[STIX]{x1D702}_{1}=-F_{1}(\unicode[STIX]{x1D709},\unicode[STIX]{x1D701}),\end{eqnarray}$$

where

(5.8)

The magnetic field to first order can be expressed as

(5.9) $$\begin{eqnarray}\bar{B}\approx \bar{B}_{0}(\unicode[STIX]{x1D701})+\unicode[STIX]{x1D716}[\unicode[STIX]{x1D702}_{1}\bar{B}_{0}^{\prime }(\unicode[STIX]{x1D701})+\bar{B}_{1}(\unicode[STIX]{x1D701})],\end{eqnarray}$$

where we have introduced the small parameter $\unicode[STIX]{x1D716}$ , and introduced $\bar{B}_{1}$ , which can account for changes to the functional form of $\bar{B}$ (small modifications of well depths, heights, locations, etc.).

Thus, the problem of specifying an omnigenous field has been reduced to specifying the zeroth-order quasi-symmetric magnetic field, and adding an arbitrary, small deviation to the CS mapping via the perturbation $f_{1}$ . There is a great deal of freedom in specifying $f_{1}$ : it must only satisfy the appropriate continuity conditions, and $f_{1}|_{\unicode[STIX]{x1D702}=0}=0$ , and, optionally, continuity in derivatives of a desired order.

6 Compatibility of omnigenous and near-axis forms of $B$

In the near-axis expansion, the magnetic field is nearly quasi-symmetric in a trivial sense, because the zeroth-order contribution (i.e. the on-axis magnetic field) can only depend on the Boozer toroidal angle. The field is therefore simultaneously quasi-axisymmetric, quasi-helically symmetric and quasi-poloidally symmetric on the axis if the on-axis field is constant; if the field strength varies with the Boozer toroidal angle, it cannot be quasi-axisymmetric or quasi-helically symmetric, but it is still automatically quasi-poloidally symmetric. At a small distance from the axis, we are concerned with whether omnigenity can still be satisfied, when quasi-symmetry is broken. The first step is to check for consistency between two forms of the magnetic field strength, namely the omnigenous form and the near-axis form.

In § 6.1 we consider the case of constant on-axis magnetic field, and conclude that the two forms are inconsistent, implying that omnigenous non-quasi-symmetric equilibria of this type do not exist – at least not any that are consistent with the axis expansion to first order. The question of consistency for the case of varying on-axis field strength is then treated in § 6.2, where more optimistic conclusions are drawn.

Following Garren & Boozer (1991a ), we write the general form of the magnetic field to first order in the distance from the axis

(6.1) $$\begin{eqnarray}B(\unicode[STIX]{x1D716},\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})\approx B_{\text{a}}(\unicode[STIX]{x1D711})(1+\unicode[STIX]{x1D716}d(\unicode[STIX]{x1D711})\cos [\unicode[STIX]{x1D703}-\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D711})]),\end{eqnarray}$$

where $\unicode[STIX]{x1D716}\ll 1$ can be taken as a measure of the distance from the axis. Please note that $\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D711})$ is not to be confused with the field line label in the Clebsch representation of the magnetic field. We will compare this form to the general form of a weakly non-quasi-symmetric omnigenous magnetic field, given by (5.9), for the three symmetry classes.

6.1 Impossibility of omnigenity with constant on-axis magnetic field

Let us expand the magnetic field $\bar{B}(\unicode[STIX]{x1D702})$ , introduced in (4.1), about its constant on-axis value,

(6.2) $$\begin{eqnarray}\bar{B}\approx \bar{B}_{0}+\unicode[STIX]{x1D716}\bar{B}_{1}(\unicode[STIX]{x1D702}).\end{eqnarray}$$

Note that we do not employ the expansion of § 5 here, which is inapplicable as the first-order field $\bar{B}_{1}$ determines the magnetic well structure non-perturbatively when the zeroth-order field $\bar{B}_{0}$ is constant.

Within § 6.1, we will allow this to describe a possibly non-omnigenous field, that nevertheless does satisfy a more relaxed condition called ‘pseudo-symmetry’. Pseudo-symmetry means that the contours of magnetic field all have a common topology – i.e. they all close poloidally, all close toroidally, or all close helically – which is a necessary but not sufficient requirement for omnigenity (Subbotin et al. 1999; Helander & Nührenberg 2011). To assert pseudo-symmetry we assume that (3.2) still applies, but with $F$ not necessarily satisfying (4.3). Equating the two forms of the magnetic field, equations (6.1)–(6.2), we obtain at dominant order $\bar{B}_{0}=B_{\text{a}}$ , i.e. a constant, and at first order

(6.3) $$\begin{eqnarray}\frac{\bar{B}_{1}(\unicode[STIX]{x1D702})}{B_{\text{a}}}=d(\unicode[STIX]{x1D711})\cos (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D711})).\end{eqnarray}$$

6.1.1 Toroidally or helically closed lines of magnetic field strength

Let us consider the cases of toroidally or helically closed contours of constant magnetic field, taking $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D703}-N\unicode[STIX]{x1D711}$ and $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D711}$ , where $N=0$ for the toroidal case. Let us write $\unicode[STIX]{x1D6FC}=\widetilde{\unicode[STIX]{x1D6FC}}+N\unicode[STIX]{x1D711}$ , so that (6.3) becomes

(6.4) $$\begin{eqnarray}\bar{B}_{1}(\unicode[STIX]{x1D702})=d(\unicode[STIX]{x1D711})B_{\text{a}}\cos (\unicode[STIX]{x1D701}-\widetilde{\unicode[STIX]{x1D6FC}}(\unicode[STIX]{x1D711})).\end{eqnarray}$$

First, we note that pseudo-symmetric fields are consistent with (6.4), for instance $\bar{B}_{1}\propto \cos (\unicode[STIX]{x1D701}+c\sin (\unicode[STIX]{x1D711}))$ , with $|c|<1$ . This example is not, however, omnigenous since the contour of maximum field strength (i.e.  $\unicode[STIX]{x1D702}=0$ ) is not straight. In fact, omnigenous examples are not possible, which can be proved as follows. First we can determine $d(\unicode[STIX]{x1D711})$ by evaluating (6.4) at the global maximum, $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D701}=0$ , yielding

(6.5) $$\begin{eqnarray}d(\unicode[STIX]{x1D711})=\frac{\bar{B}_{1}(0)}{B_{\text{a}}\cos (\widetilde{\unicode[STIX]{x1D6FC}})}.\end{eqnarray}$$

Equation (6.4) then becomes

(6.6) $$\begin{eqnarray}\bar{B}_{1}(\unicode[STIX]{x1D702})=\bar{B}_{1}(0)\frac{\cos (\unicode[STIX]{x1D701}-\widetilde{\unicode[STIX]{x1D6FC}})}{\cos (\widetilde{\unicode[STIX]{x1D6FC}})}.\end{eqnarray}$$

We can show that the only omnigenous fields consistent with this form are quasi-symmetric. Note first that if $\bar{B}_{1}$ is omnigenous but not quasi-symmetric, then $\widetilde{\unicode[STIX]{x1D6FC}}$ must be a non-constant function of $\unicode[STIX]{x1D711}$ . However, the function $1/\cos (\widetilde{\unicode[STIX]{x1D6FC}})$ must then reach a global maximum at some value of $\unicode[STIX]{x1D711}$ . The coordinate $\unicode[STIX]{x1D701}$ can be independently optimized to maximize the function $\cos (\unicode[STIX]{x1D701}-\widetilde{\unicode[STIX]{x1D6FC}})$ , i.e.  $\unicode[STIX]{x1D701}=\widetilde{\unicode[STIX]{x1D6FC}}$ or $\unicode[STIX]{x1D701}=\widetilde{\unicode[STIX]{x1D6FC}}+\unicode[STIX]{x03C0}$ , depending on the desired sign. This means that the magnetic field reaches its maximum at a point in the $\unicode[STIX]{x1D701}$ – $\unicode[STIX]{x1D711}$ plane, which is inconsistent with the required topology of contours of constant $\unicode[STIX]{x1D702}$ . Therefore we reach a contradiction and conclude that $\widetilde{\unicode[STIX]{x1D6FC}}$ must be a constant, i.e. the only possible omnigenous fields with toroidally or helically closed contours of magnetic field strength are quasi-symmetric.

6.1.2 Poloidally closed lines of magnetic field strength

Finally, we consider the case of poloidally closed lines of constant field strength, taking $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D711}$ and $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D703}$ . We can show that pseudo-symmetric fields are altogether impossible. Fixing $\unicode[STIX]{x1D702}=\unicode[STIX]{x1D702}_{c}$ such that $\bar{B}_{1}(\unicode[STIX]{x1D702}_{c})\neq 0$ , we consider a poloidal transit $\unicode[STIX]{x1D703}\rightarrow \unicode[STIX]{x1D703}+2\unicode[STIX]{x03C0}$ . While the left-hand side of (6.3) remains constant, the argument of $\cos$ on the right-hand side must increase by $2\unicode[STIX]{x03C0}$ since the $\unicode[STIX]{x1D711}$ returns to its initial value at the end of the transit. Therefore, by the intermediate value theorem, this argument must pass through a zero of $\cos$ , making the right-hand side zero, contradicting the assumption that $\bar{B}_{1}\neq 0$ . Thus we find the only fields that satisfy (6.3) are $\bar{B}_{1}=0$ .

6.2 The case of varying on-axis magnetic field strength (quasi-isodynamic fields)

We now consider the on-axis magnetic field strength to vary with $\unicode[STIX]{x1D711}$ (i.e.  $\unicode[STIX]{x2202}B_{\text{a}}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}\neq 0$ almost everywhere). This restricts us to the class of poloidally closed magnetic fields (e.g. quasi-poloidal symmetric and quasi-isodynamic fields). We thus henceforth fix the angles of the CS mapping to be $\unicode[STIX]{x1D701}=\unicode[STIX]{x1D711}$ and $\unicode[STIX]{x1D709}=\unicode[STIX]{x1D703}$ . For this case we find it is possible to construct first-order omnigenous fields, at least in a certain approximate sense. Note that quasi-poloidal symmetry cannot be achieved at first order because it requires $d$ , a quantity proportional to the axis curvature, to be zero. A magnetic axis with zero curvature everywhere cannot be closed.

We now set the two forms of the magnetic field strength equal, i.e. by (5.9) and (6.1), and use $\unicode[STIX]{x1D702}_{0}=\unicode[STIX]{x1D711}$ to obtain at zeroth order in $\unicode[STIX]{x1D716}$

(6.7) $$\begin{eqnarray}B_{\text{a}}(\unicode[STIX]{x1D711})=\bar{B}_{0}(\unicode[STIX]{x1D711}),\end{eqnarray}$$

and, at first order,

(6.8) $$\begin{eqnarray}B_{\text{a}}(\unicode[STIX]{x1D711})d\cos (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D6FC})=\unicode[STIX]{x1D702}_{1}\bar{B}_{0}^{\prime }(\unicode[STIX]{x1D711})+\bar{B}_{1}(\unicode[STIX]{x1D711}).\end{eqnarray}$$

The first thing we notice is that although the left-hand side depends on $\unicode[STIX]{x1D703}$ , the term $\bar{B}_{1}$ on the right-hand side depends only on $\unicode[STIX]{x1D711}$ , therefore, it must balance and cancel with part of the other term on the right-hand side. Formally, we may filter out the $\unicode[STIX]{x1D703}$ -averaged part and ignore it since, by (6.8), it cannot affect the magnetic field strength, so it represents non-physical freedom in the function $\bar{B}(\unicode[STIX]{x1D702})$ . That is, we assume $\unicode[STIX]{x1D702}_{1}$ to henceforth have no $\unicode[STIX]{x1D703}$ -average, and set $\bar{B}_{1}=0$ without loss of generality. (This is entirely consistent with (5.8)–(5.7) since the filtering commutes with the $\unicode[STIX]{x1D701}$ -shift.)

Solving (6.8) for $\unicode[STIX]{x1D702}_{1}$ and using (5.7), we obtain

(6.9) $$\begin{eqnarray}F_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})=-\frac{\bar{B}_{0}(\unicode[STIX]{x1D711})d(\unicode[STIX]{x1D711})}{\bar{B}_{0}^{\prime }(\unicode[STIX]{x1D711})}\cos (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D711})).\end{eqnarray}$$

Now we observe that the omnigenous form of $F_{1}$ , equation (5.8), expresses a sort of symmetry that can be written as

(6.10)

Applying this condition to the expression for $F_{1}$ given in (6.9), we obtain an equation of the form $a\cos (\unicode[STIX]{x1D703}+t_{1})=b\cos (\unicode[STIX]{x1D703}+t_{2})$ from which it follows that $a=\pm b$ and $t_{1}=t_{2}+2\unicode[STIX]{x03C0}n+(\unicode[STIX]{x03C0}\mp \unicode[STIX]{x03C0})/2$ ; we choose the upper sign convention and take $n=0$ , obtaining

(6.11) $$\begin{eqnarray}d(\unicode[STIX]{x1D711})=\unicode[STIX]{x1D702}_{b}^{\prime }(\unicode[STIX]{x1D711})d(\unicode[STIX]{x1D702}_{b}(\unicode[STIX]{x1D711})),\quad \text{for }\unicode[STIX]{x1D711}\in D_{iR},\end{eqnarray}$$

and

(6.12)

To obtain (6.11), note that $\bar{B}_{0}(\unicode[STIX]{x1D702}_{b}(\unicode[STIX]{x1D711}))=\bar{B}_{0}(\unicode[STIX]{x1D711})$ implies $\bar{B}_{0}^{\prime }(\unicode[STIX]{x1D711})=\unicode[STIX]{x1D702}_{b}^{\prime }(\unicode[STIX]{x1D711})\bar{B}_{0}^{\prime }(\unicode[STIX]{x1D702}_{b}(\unicode[STIX]{x1D711}))$ . Equation (6.12) can be thought of as providing a way to construct the function $\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D711})$ for $\unicode[STIX]{x1D711}\in D_{iR}$ , given its dependence within $D_{iL}$ . Equation (6.11) provides a similar prescription for $d(\unicode[STIX]{x1D711})$ . Note that, assuming irrational , the function $\unicode[STIX]{x1D6FC}$ cannot be consistently defined to force periodicity. That is, evaluating (6.12) at the global maximum, $\unicode[STIX]{x1D711}=2\unicode[STIX]{x03C0}$ , we obtain , which can only be a multiple of $2\unicode[STIX]{x03C0}$ if  is an integer. This does not actually affect the periodicity of the first-order magnetic field, since, from (6.8) and $\bar{B}_{0}^{\prime }(0)=0$ , we must have $d(0)=0$ , implying that it is zero at $\unicode[STIX]{x1D711}=0$ (and therefore is periodic). In fact this must be true at all extrema of $\bar{B}_{0}$ ,

(6.13) $$\begin{eqnarray}d=0,\quad \text{at all local extrema,}\end{eqnarray}$$

which also follows from (6.11) if $d$ is to be continuous at these locations. Even if the field strength is periodic, this does not automatically imply periodicity in derivatives of the magnetic field strength, but this would be too much to hope for, given that omnigenous solutions are generally non-analytic (Cary & Shasharina 1997). As we will see, periodicity is also not be so easily enforced for the full solution, but in the next section we will propose a way to repair these discontinuities.

7 Constructing quasi-isodynamic equilibria

In the previous section, we have demonstrated that the first order magnetic field strength of the near-axis expansion can be made consistent with the condition of omnigenity. This alone is not sufficient to solve the equilibrium problem. We must also demonstrate that there is a consistent solution of the equilibrium equation introduced by Garren & Boozer (1991a ).

To apply the results of Garren & Boozer (1991a ) to the present problem, however, we must account for the fact that the curvature of the magnetic axis will be zero at some toroidal locations (this follows from (6.13) as shown below). The conventional Frenet frame is discontinuous at points of zero curvature, when the normal and binormal unit vectors change their signs, which can be considered the typical case since it corresponds to first-order zeros of the curvature. Fortunately, we can employ a modified Frenet frame that accounts for these sign flips; a formal development of such frames is given by Carroll, Kose & Sterling (2013). Under reasonable conditions, this frame will be continuous everywhere, and its derivatives satisfy precisely the same equations as the traditional Frenet frame, so that the derivations of Garren & Boozer (1991a ) proceed identically, with only the trivial modification of replacing the curvature with the ‘signed curvature’  $\unicode[STIX]{x1D705}^{s}(\unicode[STIX]{x1D711})=s(\unicode[STIX]{x1D711})\unicode[STIX]{x1D705}(\unicode[STIX]{x1D711})$ , and reinterpreting the components of the coordinate mapping using the modified frame. Note that the sign $s(\unicode[STIX]{x1D711})$ takes values $+1$ or $-1$ , and switches at locations of zero curvature, and the normal and binormal unit vectors of the modified Frenet frame are given by multiplying those of the traditional Frenet frame by $s(\unicode[STIX]{x1D711})$ .

The equation that must be solved to find an equilibrium solution is a first-order nonlinear ordinary differential equation (ODE) for the quantity $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D711})$ , which can be interpreted as follows with the modified Frenet frame: the coordinate mapping is expressed to first order as $\boldsymbol{x}\approx \boldsymbol{r}_{0}+X_{1}\boldsymbol{n}^{s}+Y_{1}\boldsymbol{t}^{s}$ , where $\boldsymbol{r}_{0}(\unicode[STIX]{x1D711})$ specifies the magnetic axis, and, via its derivatives, also the signed normal and binormal unit vectors, $\boldsymbol{n}^{s}(\unicode[STIX]{x1D711})$ and $\boldsymbol{t}^{s}(\unicode[STIX]{x1D711})$ . The iterative solution of the equilibrium equations leads to expressions for the components $X_{1}$ and $Y_{1}$ , so that the quantity $\unicode[STIX]{x1D70E}$ is seen to be related to the amplitude of the in-phase part of the binormal component of the coordinate mapping (see Garren & Boozer (1991a ) and Landreman & Sengupta (2018) for more details),

(7.1) $$\begin{eqnarray}\displaystyle & \displaystyle X_{1}=\bar{d}(\unicode[STIX]{x1D711})\cos [\unicode[STIX]{x1D703}-\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D711})], & \displaystyle\end{eqnarray}$$

(7.2) $$\begin{eqnarray}\displaystyle & \displaystyle Y_{1}=\frac{2}{\bar{d}(\unicode[STIX]{x1D711})}\{\sin [\unicode[STIX]{x1D703}-\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D711})]+\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D711})\cos [\unicode[STIX]{x1D703}-\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D711})]\}, & \displaystyle\end{eqnarray}$$

where $\bar{d}$ is related to $d$ as

(7.3) $$\begin{eqnarray}d(\unicode[STIX]{x1D711})=\bar{d}(\unicode[STIX]{x1D711})\unicode[STIX]{x1D705}^{s}(\unicode[STIX]{x1D711}).\end{eqnarray}$$

Note that from this and (6.13) we see that the curvature must indeed be zero at extrema of $B_{a}$ . Although $d$ generally changes sign at points of zero curvature, $\bar{d}$ can be continuous and non-zero at these places. Note that $\bar{d}$ cannot have any zeros, for $Y_{1}$ to remain finite. Because of this, $\unicode[STIX]{x1D705}^{s}$ and $d$ must have zeros of the same order. Furthermore, the periodicity of $X_{1}$ cannot be achieved by setting the coefficient, $\bar{d}$ , of (7.1) to zero. Instead, periodicity of $X_{1}$ and $Y_{1}$ implies the constraints

(7.4) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{d}(2\unicode[STIX]{x03C0})=\bar{d}(0), & \displaystyle\end{eqnarray}$$

(7.5) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6FC}(2\unicode[STIX]{x03C0})-\unicode[STIX]{x1D6FC}(0)=2\unicode[STIX]{x03C0}N,\quad \text{for some integer }N, & \displaystyle\end{eqnarray}$$

(7.6) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D70E}(2\unicode[STIX]{x03C0})=\unicode[STIX]{x1D70E}(0), & \displaystyle\end{eqnarray}$$

where the first should already be satisfied if $d$ , satisfies (6.11). Given irrational , the second of these conditions is not consistent with omnigenity (6.12), which implies . We will propose a possible way to overcome this in the next section.

Finally, the equation that must be solved for $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D711})$ is

(7.7)

where $P=1+\bar{d}^{4}/4$ , $Q=-G_{0}\bar{d}^{2}(\unicode[STIX]{x1D70F}+I_{2}/2)$ , with $G_{0}$ , $I_{2}$ and $\unicode[STIX]{x1D70F}(\unicode[STIX]{x1D711})$ being related to the poloidal current, toroidal current and torsion of the magnetic axis, respectively.

Let us now summarize the problem of finding an omnigenous magnetic field at first order: first one can freely specify the on-axis magnetic field strength $\bar{B}_{0}(\unicode[STIX]{x1D711})$ , an axis shape, (having zero curvature at the extrema of $\bar{B}_{0}$ ) and a function $d(\unicode[STIX]{x1D711})$ satisfying (6.11). A function $\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D711})$ can also be arbitrarily specified on $D_{L}$ but its dependence on $D_{R}$ is constrained by (6.12) which depends on . Therefore the solution of (7.7) for  $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D711})$ and , must be done consistently; see § 8 for the demonstration of this. In principle, this would complete the solution. However, as we have discussed, there is a conflict between (6.12) and the periodicity constraint on $\unicode[STIX]{x1D6FC}$ , equation (7.5), and in the following section we propose a practical way to address this.

7.1 Controlled approximation of omnigenous fields

Because an omnigenous magnetic field (that is not quasi-symmetric) is necessarily non-analytic, such a field can only ultimately be physically realized by a smooth approximation. One approach, proposed by Cary & Shasharina (1997), is to truncate a series representation of an exactly omnigenous field, resulting in a smooth but slightly non-omnigenous one. Another idea is to specify a sufficiently smooth magnetic field that violates omnigenity in a controlled way – this, as we show below, can also simultaneously help resolve the conflict between omnigenity and periodicity, which was encountered in the previous section.

To ensure the appropriate behaviour of $\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D711})$ , equation (7.5), we introduce small matching regions near $\unicode[STIX]{x1D711}=0,2\unicode[STIX]{x03C0}$ , where $\unicode[STIX]{x1D6FC}$ will be defined such that condition (7.5) is satisfied, while abandoning the condition of omnigenity there,

(7.8) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D711})=\left\{\begin{array}{@{}ll@{}}\unicode[STIX]{x1D6FC}_{\text{I}}(\unicode[STIX]{x1D711}),\quad & \text{for }0\leqslant \unicode[STIX]{x1D711}\leqslant \unicode[STIX]{x1D6FF},\\ \unicode[STIX]{x1D6FC}_{\text{II}}(\unicode[STIX]{x1D711}),\quad & \text{for }\unicode[STIX]{x1D6FF}\leqslant \unicode[STIX]{x1D711}\leqslant 2\unicode[STIX]{x03C0}-\unicode[STIX]{x1D6FF},\\ \unicode[STIX]{x1D6FC}_{\text{III}}(\unicode[STIX]{x1D711}),\quad & \text{for }2\unicode[STIX]{x03C0}-\unicode[STIX]{x1D6FF}\leqslant \unicode[STIX]{x1D711}\leqslant 2\unicode[STIX]{x03C0},\end{array}\right. & & \displaystyle\end{eqnarray}$$

In region II where omnigenity is satisfied, equation (6.12) implies that $\unicode[STIX]{x1D6FC}_{\text{II}}$ may be written in terms of the function $a(\unicode[STIX]{x1D711})$ , defined on $D_{L}$ ,

(7.9)

Introducing , extending the definition of $\unicode[STIX]{x1D702}_{b}(\unicode[STIX]{x1D711})$ to $D_{L}$ such that it is its own inverse, $\unicode[STIX]{x1D702}_{b}=\unicode[STIX]{x1D702}_{b}^{-1}$ , we can write $\unicode[STIX]{x1D6FC}$ in a more symmetric form,

(7.10)

where we have used that $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}(\unicode[STIX]{x1D711})=\unicode[STIX]{x1D711}_{b}(\unicode[STIX]{x1D711})-\unicode[STIX]{x1D711}$ , so $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}(\unicode[STIX]{x1D711}_{b}(\unicode[STIX]{x1D711}))=-\unicode[STIX]{x0394}\unicode[STIX]{x1D711}(\unicode[STIX]{x1D711})$ .

To make further progress, let us make some simplifying assumptions. We assume symmetric, single-well $\bar{B}_{0}(\unicode[STIX]{x1D711})$ , so that $\unicode[STIX]{x1D711}_{b}(\unicode[STIX]{x1D711})=2\unicode[STIX]{x03C0}-\unicode[STIX]{x1D711}$ , and $\unicode[STIX]{x0394}\unicode[STIX]{x1D711}(\unicode[STIX]{x1D711})=2\unicode[STIX]{x1D711}-2\unicode[STIX]{x03C0}$ . Equation (7.7) can be written as

(7.11) $$\begin{eqnarray}\unicode[STIX]{x1D70E}^{\prime }+\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D711})(\unicode[STIX]{x1D70E}^{2}+P)+Q=0,\end{eqnarray}$$

where . The functions  $\unicode[STIX]{x1D6FC}_{\text{I}}$ and $\unicode[STIX]{x1D6FC}_{\text{III}}$ must be chosen such that $\unicode[STIX]{x1D6FC}(2\unicode[STIX]{x03C0})-\unicode[STIX]{x1D6FC}(0)=2\unicode[STIX]{x03C0}N$ , but are otherwise free. However, a particularly simple prescription is found as follows. We note that $\unicode[STIX]{x1D6FC}_{\text{II}}$ , as expressed in (7.10), is composed of a secular piece (), and an even period piece, which need not be corrected to satisfy periodicity. Thus, we extend $\unicode[STIX]{x1D6FC}_{\text{II}}$ into the buffer regions, and add a correction to the secular part. That is, we write $\unicode[STIX]{x1D6FC}_{\text{I}}=\unicode[STIX]{x1D6FC}_{\text{II}}+\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}_{\text{I}}$ and $\unicode[STIX]{x1D6FC}_{\text{III}}=\unicode[STIX]{x1D6FC}_{\text{II}}+\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}_{\text{III}}$ , where continuity of $\unicode[STIX]{x1D6FC}$ requires that the correction goes to zero at the interface with the buffer ( $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}_{\text{I}}(\unicode[STIX]{x1D6FF})=0$ , and $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}_{\text{III}}(2\unicode[STIX]{x03C0}-\unicode[STIX]{x1D6FF})=0$ ), and the periodicity constraint is satisfied via $\unicode[STIX]{x1D6FC}_{\text{I}}(0)=-N\unicode[STIX]{x03C0}$ and $\unicode[STIX]{x1D6FC}_{\text{III}}(2\unicode[STIX]{x03C0})=N\unicode[STIX]{x03C0}$ . Taking $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}_{\text{I}}(\unicode[STIX]{x1D711})$ and $\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}_{\text{III}}(\unicode[STIX]{x1D711})$ to be linear functions of $\unicode[STIX]{x1D711}$ then uniquely determines them, resulting in

(7.12)

The function $\unicode[STIX]{x1D6FE}$ is then simply

(7.13)

where $\unicode[STIX]{x1D6FE}_{\text{o}}(\unicode[STIX]{x1D711})$ is the odd function

(7.14) $$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6FE}_{\text{o}}(\unicode[STIX]{x1D711})=\left\{\begin{array}{@{}ll@{}}-c^{\prime }(\unicode[STIX]{x1D711}),\quad & \text{for }\unicode[STIX]{x1D702}\in [0,\unicode[STIX]{x03C0}]\\ c^{\prime }(2\unicode[STIX]{x03C0}-\unicode[STIX]{x1D711}),\quad & \text{for }\unicode[STIX]{x1D702}\in [\unicode[STIX]{x03C0},2\unicode[STIX]{x03C0}].\end{array}\right. & & \displaystyle\end{eqnarray}$$

We note that the integer $N$ , introduced in (7.5), must be consistently set, taking into account the number of times the coordinate mapping rotates around the magnetic axis during a toroidal transit, so that $\unicode[STIX]{x1D703}$ behaves as a proper poloidal angle – i.e. it increases by $2\unicode[STIX]{x03C0}$ with a poloidal transit, but is periodic in the toroidal direction.

Note that  only appears in the definition of  $\unicode[STIX]{x1D6FE}(\unicode[STIX]{x1D711})$ in the matching regions, so we may consider solving (7.11) independently in the omnigenous region. Then,  can be tuned in the matching regions to achieve $2\unicode[STIX]{x03C0}$  periodicity of $\unicode[STIX]{x1D70E}$ . Assuming that, for small $\unicode[STIX]{x1D6FF}$ , the functions $P$ , $Q$ and $\unicode[STIX]{x1D6FE}$ can be considered constant in the matching regions, the proof of the existence and uniqueness of this value of  in fact follows easily from the results of Landreman et al. (2019). A general existence and uniqueness theorem, i.e. allowing general $\unicode[STIX]{x1D6FC}_{\text{I}}$ , $\unicode[STIX]{x1D6FC}_{\text{III}}$ and $\bar{B}_{0}$ , appears more difficult to prove.

8 Numerical solutions

8.1 Algorithm

We now describe how the equations of § 6.2 can be solved in practice numerically. The inputs to the calculation are the shape of the magnetic axis, the functions $B_{0}=\bar{B}_{0}(\unicode[STIX]{x1D711})$ and $d(\unicode[STIX]{x1D711})$ , a chosen width for the buffer regions $\unicode[STIX]{x1D6FC}_{\text{I}}$ and $\unicode[STIX]{x1D6FC}_{\text{III}}$ and a finite value of aspect ratio. The outputs of the calculation are ,  $B_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})=B_{0}(\unicode[STIX]{x1D711})d\cos (\unicode[STIX]{x1D703}-\unicode[STIX]{x1D6FC})$ and the shape of the elliptical flux surfaces.

The first step in the numerical solution is to compute $G_{0}$ and $\unicode[STIX]{x1D711}(\unicode[STIX]{x1D719})$ along the axis, where $G_{0}$ is the on-axis value of the coefficient $G$ in $\boldsymbol{B}_{\text{cov}}$ , and $\unicode[STIX]{x1D719}$ is the standard toroidal angle (i.e. the azimuthal angle for cylindrical coordinates). These calculations can be done by iteratively solving (2.20)–(2.22) in Landreman & Sengupta (2018),

(8.1) $$\begin{eqnarray}\displaystyle \frac{\text{d}\unicode[STIX]{x1D711}}{\text{d}\unicode[STIX]{x1D719}}=\frac{\ell ^{\prime }B_{0}}{|G_{0}|},\quad |G_{0}|=\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{2\unicode[STIX]{x03C0}}\text{d}\unicode[STIX]{x1D719}B_{0}(\unicode[STIX]{x1D711}(\unicode[STIX]{x1D719}))\ell ^{\prime }, & & \displaystyle\end{eqnarray}$$

where $\ell ^{\prime }=|\text{d}\boldsymbol{r}/\text{d}\unicode[STIX]{x1D719}|$ is the arclength increment. Beginning with the guess $\unicode[STIX]{x1D711}\approx \unicode[STIX]{x1D719}$ , fixed-point (Picard) iteration converges quickly. The sign of $G_{0}$ (i.e. whether $\boldsymbol{B}$ points towards increasing or decreasing $\unicode[STIX]{x1D711}$ ) is a free input.

Next we solve the equation for $\unicode[STIX]{x1D70E}$ , equation (7.7), using Newton’s method. This solution procedure is similar to the one for quasi-symmetry described in § 3 of (Landreman et al. 2019), but with a few modifications. The discrete unknowns include values of $\unicode[STIX]{x1D70E}$ on a uniform grid of $N_{\unicode[STIX]{x1D719}}$ points in the domain $\unicode[STIX]{x1D719}\in [0,2\unicode[STIX]{x03C0}/n_{fp}]$ , where $n_{fp}$ is the number of identical field periods. The $\text{d}\unicode[STIX]{x1D70E}/\text{d}\unicode[STIX]{x1D711}$ term in (7.7) is discretized using the pseudospectral differentiation matrix. For computing omnigenous solutions, it is convenient if none of the grid points lie exactly where the curvature vanishes; otherwise (7.7) would require evaluating $d/\unicode[STIX]{x1D705}=0/0$ . To avoid points of vanishing curvature at both $\unicode[STIX]{x1D719}=0$ and $\unicode[STIX]{x1D719}=\unicode[STIX]{x03C0}/n_{fp}$ (half-period), the grid points are shifted by one third of the grid spacing: $\unicode[STIX]{x1D719}_{j}=(j-2/3)2\unicode[STIX]{x03C0}/n_{fp}$ for $j=1\ldots N_{\unicode[STIX]{x1D719}}$ .

In addition to the unknowns $\unicode[STIX]{x1D70E}(\unicode[STIX]{x1D719}_{j})$ there is one more unknown: . Corresponding to this additional unknown is one additional equation, reflecting the freedom to specify $\unicode[STIX]{x1D70E}(0)$ (as discussed in the Appendix of Landreman et al. (2019)). In the common case of stellarator symmetry, this extra condition is $\unicode[STIX]{x1D70E}(0)=0$ . While $\unicode[STIX]{x1D711}=0$ is not one of the grid points, this condition can nonetheless be imposed by interpolating $\unicode[STIX]{x1D70E}$ from the grid points $\unicode[STIX]{x1D719}_{j}$ (using pseudospectral interpolation).

To apply Newton’s method, the Jacobian is needed. One block of the Jacobian corresponds to the derivative of (7.7) with respect to $\unicode[STIX]{x1D70E}$ , which is straightforward. Another block corresponds to the derivative of (7.7) with respect to , and in contrast to the quasi-symmetric case, here we must account for the fact that $\unicode[STIX]{x1D6FC}$ depends on . Fortunately, in both the central region and buffer regions, $\unicode[STIX]{x1D6FC}$ depends on  as 

(8.2)

for some functions  and  $\unicode[STIX]{x1D6FC}_{0}(\unicode[STIX]{x1D711})$ . Therefore the Jacobian block corresponding to the derivative of (7.7) with respect to  is given by .

Finally, once $\unicode[STIX]{x1D70E}$ with self-consistent  has been found, the result can be converted to cylindrical coordinates and a finite-aspect-ratio configuration can be generated using any of the methods described in § 4 of Landreman et al. (2019). For results shown below, we use the method of § 4.1. The resulting toroidal boundary shape in cylindrical coordinates can then be supplied to an equilibrium code such as VMEC (Hirshman & Whitson 1983).

Instead of solving equation (7.7) for $\unicode[STIX]{x1D70E}$ , it is also possible to directly solve the equations for the near-axis flux surface shape in cylindrical coordinates derived in Landreman & Sengupta (2018). This approach is conceptually valuable since it makes no use of the Frenet frame, so one need not worry about discontinuities of the Frenet representation. The numerical solution can be obtained using Newton’s method, as detailed in § 4.3 of Landreman et al. (2019). In contrast to the quasi-symmetric case, for the omnigenous case $B_{1}$ depends on  through $\unicode[STIX]{x1D6FC}$ . The block of the Jacobian corresponding to the derivative of the residual with respect to  can again be evaluated by noting (8.2). We have implemented both this cylindrical coordinate approach and the Frenet ( $\unicode[STIX]{x1D70E}$ ) approach, and verified the results from the two methods converge towards each other as $N_{\unicode[STIX]{x1D719}}\rightarrow \infty$ . For instance, in the example of the next section with $N_{\unicode[STIX]{x1D719}}=101$ grid points, the two approaches yield the same value of   through 9 digits.

8.2 Example

An example of an omnigenous non-quasi-symmetric configuration constructed using this procedure is shown in figure 2–5. We begin by choosing the axis shape

(8.3a,b ) $$\begin{eqnarray}\displaystyle R_{0}(\unicode[STIX]{x1D719})=1-0.2\cos (2\unicode[STIX]{x1D719}),\quad z_{0}(\unicode[STIX]{x1D719})=0.35\sin (2\unicode[STIX]{x1D719}), & & \displaystyle\end{eqnarray}$$

shown in figure 2(a), which has vanishing curvature at two points: $\unicode[STIX]{x1D719}=0$ and $\unicode[STIX]{x03C0}$ . The Frenet normal and binormal vectors both display a discontinuous reversal in direction at these points. When viewed from above, e.g. from a viewpoint with $R=0$ and $z>0$ , the axis has the shape of a racetrack oval, with the points of vanishing curvature at the middle of each straightaway. The oval is twisted out of the $z=0$ plane so it has net torsion. We choose the on-axis field strength to be $B_{0}(\unicode[STIX]{x1D711})=1+0.1\cos \unicode[STIX]{x1D711}$ . Note that while the axis shape is symmetric under rotation by $\unicode[STIX]{x03C0}$ , as in a device with two field periods, the field strength does not have this symmetry, so the number of field periods is one.

Figure 2. Elements of the numerical construction include (a) the shape of the magnetic axis, (b) the function $d(\unicode[STIX]{x1D711})$ , (c) the phase $\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D711})$ and (d) the first-order (in $\unicode[STIX]{x1D716}$ ) variation of the field strength, $B_{1}$ .

Figure 3. The shape of the constructed omnigenous configuration, for aspect ratio 20. (a) Boundary shape in three dimensions, with colour indicating $B$ computed by VMEC. A bird’s eye view and side view are shown. (b) Cross-sections of the configuration at 8 values of toroidal angle $\unicode[STIX]{x1D719}$ .

The elongation of the final configuration is sensitive to the choice of $d(\unicode[STIX]{x1D711})$ . It can be challenging to find a function $d(\unicode[STIX]{x1D711})$ for which the elongation does not grow to impractically large values, $\gg 10$ . For the configuration in the figures we use

(8.4) $$\begin{eqnarray}\displaystyle d(\unicode[STIX]{x1D711})=1.08\sin (\unicode[STIX]{x1D711})+0.26\sin (2\unicode[STIX]{x1D711})+0.46\sin (3\unicode[STIX]{x1D711}), & & \displaystyle\end{eqnarray}$$

where the coefficients were chosen to keep the elongation down at tolerable levels, ${\sim}5$ . Further tailoring of the $d(\unicode[STIX]{x1D711})$ function might yield lower values of elongation.

We choose $\unicode[STIX]{x1D6FC}=-3\unicode[STIX]{x03C0}/2$ at $\unicode[STIX]{x1D711}=0$ and $\unicode[STIX]{x1D6FC}=\unicode[STIX]{x03C0}/2$ at $\unicode[STIX]{x1D711}=2\unicode[STIX]{x03C0}$ . In region II, . We allocate the first 10 % and last 10 % of the toroidal domain to the buffer regions: $\unicode[STIX]{x1D6FC}_{\text{I}}$ occupies $\unicode[STIX]{x1D711}\in [0,\unicode[STIX]{x03C0}/5]$ and $\unicode[STIX]{x1D6FC}_{\text{III}}$ occupies $\unicode[STIX]{x1D711}\in [9\unicode[STIX]{x03C0}/5,2\unicode[STIX]{x03C0}]$ . In the buffer regions we choose $\unicode[STIX]{x1D6FC}$ to be a fourth-order polynomial in $\unicode[STIX]{x1D711}$ , enforcing continuity of $\unicode[STIX]{x1D6FC}$ and its first two derivatives at the domain boundaries $\unicode[STIX]{x1D711}=0$ , $\unicode[STIX]{x03C0}/5$ , and $9\unicode[STIX]{x03C0}/5$ . (When $\unicode[STIX]{x1D6FC}$ does not have a continuous derivative at these boundaries, the VMEC solution exhibits numerical oscillations near these points.)

The shape of the resulting configuration is shown in figure 3, for aspect ratio 20. It can be seen that the shape resembles a Möbius strip. There is a pronounced twist in the shape near the region of maximum field strength. The toroidal extent of this twist corresponds directly to the width of the buffer regions. The constructed boundary shape is provided as input to the VMEC equilibrium code (Hirshman & Whitson 1983). According to the near-axis construction, the rotational transform is . The rotational transform computed by VMEC is similar (0.70 on axis, 0.72 at the edge), and it converges to the value predicted by the construction as the aspect ratio is increased.

Given a VMEC solution, we then evaluate $B$ as a function of Boozer coordinates using the code BOOZ_XFORM (Sanchez et al. 2000), to compare the achieved $B$ to the $B$ predicted by the near-axis construction. We find that $B$ for the numerical equilibrium inside the constructed boundary (according to BOOZ_XFORM) converges to the desired function $B(r,\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$ as the aspect ratio $A$ increases. One aspect of this convergence can be seen in figure 4, which shows the convergence of the on-axis $B$ to the target function $1+0.1\cos \unicode[STIX]{x1D711}$ . For $A\geqslant 80$ , differences from the target function are barely discernible on the scale of the figure. Furthermore, the convergence of $B$ at the boundary can be seen in figure 5. The total $B$ is shown in the left column, with the target $B_{0}$ subtracted in the right column, and the three rows show three increasing values of aspect ratio. As $A$ increases, the total $B$ converges to the desired function, and $B-B_{0}$ converges to the desired $B_{1}$ . Finally, figure 6 shows how the difference between the target and achieved magnetic field strength scales with $A$ . The difference between the target and achieved $B$ is measured by the root-mean-square (r.m.s.) difference $[\int \!\text{d}\unicode[STIX]{x1D703}\!\int \!\text{d}\unicode[STIX]{x1D711}(B_{\text{VMEC}}-B_{\text{construction}})^{2}]^{1/2}$ . In this formula, $B$ at the boundary is used. This r.m.s. difference is found to scale as $1/A^{2}$ , as expected since the construction is done here through $O(\unicode[STIX]{x1D716})$ .

Figure 4. As the aspect ratio $A$ used for the construction increases, $B$ on the axis of the numerical VMEC equilibrium inside the constructed boundary (solid coloured curves) converges to the desired target function (dotted grey curve).

Figure 5. As the aspect ratio $A$ used for the construction increases, $B$ for the numerical VMEC equilibrium inside the constructed boundary (solid curves) converges to the desired target function (dotted curves). Panels (a,c,e) shows the total $B$ at the boundary; (b,d,f) shows $B-(1+0.1\cos \unicode[STIX]{x1D711})$ at the boundary.

Figure 6. The difference between $B$ predicted by the construction and $B$ computed by VMEC $+$ BOOZ_XFORM inside the constructed boundary, as measured by the root-mean-square $[\int \text{d}\unicode[STIX]{x1D703}\int \text{d}\unicode[STIX]{x1D711}(B_{\text{VMEC}}-B_{\text{construction}})^{2}]^{1/2}$ , scales as $A^{-2}$ . This scaling is consistent with the fact that the construction here is carried out through $O(\unicode[STIX]{x1D716})$ .

8.3 Measuring deviation from omnigenity via $1/\unicode[STIX]{x1D708}$ transport

The collisionless confinement of trapped particle orbits also has a direct effect on neoclassical transport, which can be observed in the so-called $1/\unicode[STIX]{x1D708}$ transport, determined by the effective helical ripple, $\unicode[STIX]{x1D716}_{\text{eff}}$ (Nemov et al. 1999),

(8.5) $$\begin{eqnarray}\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}=\left(\frac{\unicode[STIX]{x03C0}R^{2}}{2^{7/2}}\right)\lim _{L_{s}\rightarrow \infty }\left(\int _{0}^{L_{s}}\frac{\text{d}l}{B}\right)\left(\int _{0}^{L_{s}}\frac{\text{d}l|\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}|}{B}\right)^{-2}\int _{B_{\text{min}}/B_{\text{r}}}^{B_{\text{max}}/B_{\text{r}}}\text{d}x\mathop{\sum }_{j}\frac{H_{j}^{2}(x)}{I_{j}(x)},\end{eqnarray}$$

where the summation of the magnetic well index $j$ includes all wells in the interval $[0,L_{s}]$ containing points where $B=xB_{\text{r}}$ is satisfied, $R$ is a reference major radius, $B_{\text{r}}$ is a reference magnetic field strength, and

(8.6) $$\begin{eqnarray}\displaystyle & \displaystyle I_{j}=\mathop{\sum }_{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D6FE}\int _{B_{\text{min}}^{j}}^{xB_{\text{r}}}\frac{\text{d}B}{B\unicode[STIX]{x2202}B/\unicode[STIX]{x2202}s}\sqrt{1-\frac{B}{B_{\text{r}}x}}, & \displaystyle\end{eqnarray}$$

(8.7) $$\begin{eqnarray}\displaystyle & \displaystyle H_{j}=\frac{1}{x}\int _{B_{\text{min}}^{j}}^{xB_{\text{r}}}\frac{\text{d}B}{B^{2}}\sqrt{x-\frac{B}{B_{\text{r}}}}\left(\frac{4B_{\text{r}}}{B}-\frac{1}{x}\right)\mathop{\sum }_{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D6FE}Y. & \displaystyle\end{eqnarray}$$

Here we again encounter the factor $Y$ of (2.5) that arises in the calculation of the bounce-averaged radial excursion $\unicode[STIX]{x0394}\unicode[STIX]{x1D713}$ . This expression is complicated, but we note that it is essentially the square of the average of $\sum _{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D6FE}Y$ . For a perfectly omnigenous field, we therefore expect $\unicode[STIX]{x1D716}_{\text{eff}}=0$ . To test the quality of one of our constructed solutions, we generate an ‘exact’ numerical solution by using the magnetic surface shape at a finite radial position as an input for the VMEC code. Because we expect omnigenity to be satisfied only at first order (and only in the limit that the buffer regions are small), we expect $\sum _{\unicode[STIX]{x1D6FE}}\unicode[STIX]{x1D6FE}Y\sim O(\unicode[STIX]{x1D716}^{2})$ . Noting the factor of $|\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}|^{-2}\sim O(\unicode[STIX]{x1D716}^{-2})$ in (8.5), we find $\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}\sim O(\unicode[STIX]{x1D716}^{2})$ .

This predicted scaling is borne out numerically, as shown in figure 7. There, the numerical construction is carried out for several values of the buffer region size $\unicode[STIX]{x1D6FF}$ and the boundary aspect ratio $A$ . For each case, the equilibrium is computed using the VMEC code, and then the radial profile of $\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}$ is evaluated using the NEO code (Nemov et al. 1999). The numerical results show a scaling $\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}\propto \unicode[STIX]{x1D713}\propto \unicode[STIX]{x1D716}^{2}$ as expected.

Figure 7. For sufficiently large aspect ratio $A$ and small buffer region width $\unicode[STIX]{x1D6FF}$ , the $1/\unicode[STIX]{x1D708}$ transport magnitude $\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}$ for constructed configurations is found numerically to be proportional to toroidal flux, as expected from the analytic calculation in § 8.3.

Figure 8 demonstrates that the construction for omnigenity results in reduced $\unicode[STIX]{x1D716}_{\text{eff}}$ compared to non-optimized configurations of otherwise similar geometry. In particular, we compare the configuration of § 8.2 to configurations with the same magnetic axis shape but circular cross-section in the plane perpendicular to the axis. We consider two types of these latter configurations. In the first, shown in green in figure 8, the radius of the circular cross-sections is independent of toroidal angle, leading to (nearly) constant $B$ along the axis. In the second type of non-optimized configuration, shown in blue in figure 8, the radius of the circular cross-sections is made to vary with toroidal angle as $\propto 1/\sqrt{1+0.1\cos \unicode[STIX]{x1D711}}$ , so $B_{0}(\unicode[STIX]{x1D711})$ is matched to that of the constructed omnigenous configurations. Results are shown for two values of aspect ratio $A$ , 10 (solid curves) and 80 (dashed). For each configuration, a numerical equilibrium is computed with the VMEC code and $\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}$ is then computed using the NEO code (Nemov et al. 1999). At each value of aspect ratio, the constructed omnigenous configuration has smaller $\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}$ than either of the non-optimized configurations. Thus, the procedure of § 8 does appear to be a practical way to generate finite-aspect-ratio configurations with reduced $1/\unicode[STIX]{x1D708}$ transport.

Figure 8. For given aspect ratio $A$ , the omnigenous construction (red) results in lower $1/\unicode[STIX]{x1D708}$ transport magnitude $\unicode[STIX]{x1D716}_{\text{eff}}^{3/2}$ compared to configurations with the same magnetic axis shape but circular cross-section (green and blue).

9 Conclusion

We have demonstrated that it is possible to directly construct approximately quasi-isodynamic magnetic equilibria near the magnetic axis, with low computational cost, as compared to conventional optimization. These solutions are valid to first order in the distance from the magnetic axis, and also satisfy omnigenity (zero bounce-averaged radial drift) at that order for all particle orbits except a small fraction that have bounce points in the neighbourhood of the point of maximum magnetic field strength. This unconfined fraction, $f_{\text{no}}$ , can, in principle, be made arbitrarily small, at the cost of introducing sharp behaviour in the magnetic field. It is however unnecessary to make $f_{\text{no}}$ much smaller than the square root of the plasma collisionality, since the effective scattering frequency into and out of the unconfined region in velocity space is proportional to $f_{\text{no}}^{2}$ .

Our findings imply that quasi-isodynamic fields have the only possible magnetic-field-strength topology that can be achieved near the magnetic axis that satisfies omnigenity while breaking quasi-symmetry. The present work therefore naturally complements the quasi-axial and quasi-helical cases explored in the two previous papers of this series (Landreman & Sengupta 2018; Landreman et al. 2019), giving a comprehensive set of tools for constructing stellarator equilibria optimized for collisionless particle confinement near the magnetic axis – recall that quasi-poloidal symmetry is not possible. A noteworthy advantage of quasi-isodynamic construction is the additional freedom, affording a much bigger space of possible configurations than can be achieved with quasi-symmetric construction. With the latter case, the magnetic field strength must be constant on the magnetic axis, and the main freedom is in the shape of the magnetic axis. With a quasi-isodynamic construction, the on-axis field strength is a free function, and there are two additional free functions of the toroidal angle, corresponding (roughly speaking) to the angle of poloidal rotation of the elliptical cross section ( $\unicode[STIX]{x1D6FC}(\unicode[STIX]{x1D711})$ ) and elongation of the ellipse ( $d(\unicode[STIX]{x1D711})$ ), which must only satisfy certain symmetry requirements.

We have found that the case of constant magnetic field strength is theoretically forbidden for quasi-isodynamic fields. Correspondingly, in our numerical constructions, we encountered a limit to how small the magnetic mirror amplitude could be made. This demonstrates a certain separation between optimization lines, i.e. local optimization would not likely stumble upon a quasi-axisymmetric configuration, while searching in the neighbourhood of a quasi-isodynamic one.

We studied the $1/\unicode[STIX]{x1D708}$ neoclassical transport associated with some examples of equilibria obtained numerically, and confirmed the expected theoretical scaling of that transport in the distance from the magnetic axis. We noted that quasi-symmetry here seems to have an advantage due to the trapped particle fraction tending to zero on the magnetic axis – the magnetic field cannot vary in $\unicode[STIX]{x1D711}$ there – so that $1/\unicode[STIX]{x1D708}$ transport (determined solely by trapped particles) tends more sharply to zero toward the magnetic axis. Both quasi-symmetric and quasi-isodynamic solutions, however, can be constructed with sufficiently small transport to give good (fusion relevant) stellarator design candidates.

To satisfy a certain periodicity condition in our solutions, it was necessary to introduce locations where the geometric condition of omnigenity is broken. This was done at locations of maximum magnetic field strength, to ensure that all but the marginally trapped particles remain well confined. We note however that the perturbation in the magnetic field strength goes to zero at these locations, mitigating the effect of the non-omnigenous component of the magnetic field on marginally trapped particles. Note that, generally, there exist very weakly trapped particles that drift a macroscopic distance during a bounce time, spending disproportionate time in the neighbourhood of the maximum magnetic field, but such particles cannot benefit from the cancellation of radial drift at the two ends of the magnetic well, afforded by the condition omnigenity as applied here. Perfectly confining them therefore requires a different strategy, namely satisfying zero radial drift locally, i.e. at points close to the location of maximum B (note that $\boldsymbol{v}_{d}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=0$ exactly at the maximum). However, as discussed previously, it is unnecessary to perfectly confine these particles.

We note that it appears possible to achieve smallness of the regions of unconfined orbits without incurring a cost in the smoothness of the magnetic field, if the rotational transform is close to an integer. It is a curious coincidence that Wendelstein 7-X stellarator does indeed have a near-unity rotational transform, though it was motivated by a different design principle (the island divertor concept).

Possible future continuation of work includes a stellarator design study, where traditional optimization methods are applied with initial states given by our directly constructed solutions; such an approach has the advantage of drastically reducing the parameter space that needs to be searched. It may also be possible to extend the expansion to higher order in the distance from the magnetic axis, though this would also require extension of the expansion of the CS map, as performed in § 5.