1 Introduction
Three-dimensional magnetohydrodynamic (MHD) equilibrium solutions with topologically toroidal magnetic flux are conventionally specified by the plasma current, pressure profile and the shape of a flux surface (Bauer, Betancourt & Garabedian Reference Bauer, Betancourt and Garabedian1984). In a vacuum, merely imposing the surface shape is sufficient; this corresponds to a Neumann boundary condition, since the magnetic field, which is then equal to the gradient of a scalar field, has no component perpendicular to the surface, and so the normal derivative of the scalar field is zero. It seems plausible that another kind of differential constraint might be substituted at the boundary, instead of surface shape, to constitute an alternative boundary condition (or an alternative ‘problem specification’). The condition of quasi-symmetry, a symmetry of the magnetic field strength expressed in Boozer angles (Nührenberg & Zille Reference Nührenberg and Zille1988), is one such possibility. This property ensures that the magnetic field confines collisionless particle orbits. In fact, it was argued by Garren & Boozer (Reference Garren and Boozer1991a ) that, although it is generally not possible to satisfy quasi-symmetry across the entire volume, it should be possible to do so exactly on one surface, since there is a sufficient amount of freedom in the solution. However, a proof of the existence of such solutions has not been given, nor is it known how many such solutions exist. Our understanding of such issues is even more incomplete for the general class of ‘omnigenous’ equilibria, which have good particle confinement without (necessarily) having quasi-symmetry.
A subclass of quasi-symmetry, known as quasi-axisymmetry (QAS) (Nührenberg et al. Reference Nührenberg, Sindoni, Lotz, Troyon, Gori and Vaclavik1994), is particularly interesting, as it can be thought of as a generalization of axisymmetry. It has been suggested that QAS configurations might be accessed by non-axisymmetric shaping of tokamak equilibria (Boozer Reference Boozer2008). Practically speaking, non-axisymmetric shaping is a possible avenue for mitigating major problems faced with tokamaks, related to stability and disruptions (Ku & Boozer Reference Ku and Boozer2009). The stabilizing influence of non-axisymmetric shaping has been verified experimentally in recent years (Hartwell et al. Reference Hartwell, Knowlton, Hanson, Ennis and Maurer2017), but the shaping was not designed to preserve QAS. In light of the success of optimized stellarator design, i.e. with the Helically Symmetric Experiment and Wendelstein 7-X, it seems realistic to expect that stabilization by non-axisymmetric shaping of tokamaks could be achieved while also preserving the good confinement of particle orbits. This further motivates us to reconsider the theory of QAS magnetic fields.
In the present work, we perform an expansion about axisymmetry, establishing the existence of QAS magnetic fields (i.e. satisfying QAS on a single surface), characterizing the space of these solutions and relating them to the more general class of ‘omnigenous’ solutions. A numerical method to explicitly find QAS solutions is also provided. QAS is different from other forms of quasi-symmetry, since there is a (known) case (namely true axisymmetry) where the QAS condition is exactly satisfied globally, and for this reason it has been suggested that QAS equilibria may be obtained by continuous deformation of tokamak equilibria. (Equilibria with helical or cylindrical symmetry are globally quasi-symmetric, but such equilibria do not close toroidally.) We find that QAS can indeed be satisfied at a single surface, to all orders in the expansion, lending support to the conclusion of Garren & Boozer (Reference Garren and Boozer1991a ), although not establishing it exactly (i.e. non-perturbatively). Note here that QAS is satisfied on the outer surface but the magnetic field solution is global, satisfying the magnetostatic equations all the way to the magnetic axis. This approach can be contrasted with local equilibrium theory, which solves the magnetostatic equilibrium equations only in the region close to a given flux surface (Hegna Reference Hegna2000; Boozer Reference Boozer2002).
In our expansion, the deviation from axisymmetry is controlled by the (arbitrary) small parameter
$\unicode[STIX]{x1D716}$
. We consider the vacuum case, because it is significantly simpler than the non-vacuum case; it also has the convenient feature that any obtained rotational transform can be attributed entirely to three-dimensional (3-D) shaping, since none can arise from the zeroth-order axisymmetric solution.
Our findings can be summarized as follows. Due to axisymmetry at zeroth order, the linear first-order equations have an ignorable coordinate, and thus a free parameter, the toroidal mode number
$N$
. We treat this as a parameter of first-order deformation, and note that it can then be interpreted as the number of field periods of the stellarator, since higher-order contributions only involve harmonics of this fundamental mode number. From the first-order equations, a single linear second-order partial differential equation (PDE) for a single scalar field is obtained, whose solution completely determines the magnetic field. QAS is expressed as a simple differential constraint, and it is demonstrated that it cannot be satisfied globally, except for the trivial case where axisymmetry is preserved by the perturbation. The constraint of QAS can, however, be applied on a single surface as a non-standard boundary condition for the PDE, in which case the problem takes the form of an ‘oblique derivative’ problem. The solution of this mathematical problem is generally non-unique, but for our case (a homogeneous elliptic PDE), exactly two linearly independent solutions are guaranteed to exist. We conclude that there is a certain freedom in these single-surface QAS solutions: (i) the freedom of an arbitrary shape for the zeroth-order axisymmetric flux surfaces; (ii) the toroidal mode number
$N$
of the deformation (toroidal field period); and (iii) two complex amplitudes for the independent solutions of the oblique derivative problem.
At second order, the magnetostatic equation, with the QAS condition, is again reformulated as an oblique derivative problem, but in this case the system is non-homogeneous. The problem is known to have solutions, but the number of solutions cannot generally be predicted. The externally induced rotational transform (i.e. due to the non-axisymmetric deformation) is obtained at second order (but in terms of first-order quantities), i.e. it scales as
$\unicode[STIX]{x1D716}^{2}$
– this is the amount that can be attributed to the non-axisymmetric deformation. At third and higher order, the problems takes a form similar to the second-order problem, and therefore a solution is guaranteed to exist at each order. We conclude that QAS can be satisfied at all orders in the expansion.
At the surface where QAS is satisfied, a local magnetostatic condition is derived, from which a general conclusion can be drawn about QAS magnetic fields, namely that the non-axisymmetric perturbation in the field is confined to the inboard side at low aspect ratio, and the effect is amplified at high toroidal mode number. This may be a useful consideration in the design of QAS devices. Finally, we observe that, as an alternative to QAS, the condition of omnigeneity (also referred to as ‘generalized quasi-symmetry’ Landreman & Catto Reference Landreman and Catto2012) can be imposed on the field strength
$B_{0}+\unicode[STIX]{x1D716}B_{1}$
as a function of Boozer angles on a single magnetic surface. It is noted that the QAS solution represents a homogeneous solution to this general problem, and therefore represents an arbitrary freedom in the solution.
The paper is organized as follows. In § 2, the basic equations, assumptions and notation are introduced. In § 3, the problem formulation is stated, and the expansion about axisymmetry is performed. The main theoretical results are presented here, although some details are contained in the appendices. In § 4, a numerical solution of the first-order QAS problem is demonstrated, and a comparison with the widely used variational moments equilibrium code (VMEC) (Hirshman & Whitson Reference Hirshman and Whitson1983) is made to confirm the validity of the solutions. In § 5, the extension of the method to treat generalized QAS (omnigeneity) is discussed. Section 6 contains further discussion and conclusions.
2 Preliminaries
The MHD equilibrium equations are
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\times \boldsymbol{B}=\unicode[STIX]{x1D707}_{0}\,\boldsymbol{j}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{j}\times \boldsymbol{B}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\frac{\text{d}p}{\text{d}\unicode[STIX]{x1D713}}. & \displaystyle\end{eqnarray}$$
We assume topologically toroidal flux surfaces, labelled by
$\unicode[STIX]{x1D713}$
, such that
$p=p(\unicode[STIX]{x1D713})$
. To solve these equations, we use a similar approach as previous works (Garren & Boozer Reference Garren and Boozer1991a
,Reference Garren and Boozer
b
; Hegna Reference Hegna2000; Boozer Reference Boozer2002; Weitzner Reference Weitzner2014). Boozer angles (Boozer Reference Boozer1981) are denoted
$\unicode[STIX]{x1D703}$
and
$\unicode[STIX]{x1D711}$
. The contravariant form of
$\boldsymbol{B}$
is written
where 
is the rotational transform, and
$2\unicode[STIX]{x03C0}\unicode[STIX]{x1D713}$
is the toroidal magnetic flux. This form of
$\boldsymbol{B}$
satisfies zero divergence and assumes flux-surface geometry. The covariant form is written
$$\begin{eqnarray}\boldsymbol{B}_{\text{cov}}=G(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}+I(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+K(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})\unicode[STIX]{x1D735}\unicode[STIX]{x1D713},\end{eqnarray}$$
where
$G$
and
$I$
are poloidal and toroidal currents, respectively. This form is a consequence of
$\boldsymbol{j}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=0$
(2.3), and Ampere’s law (2.1), which itself implies
$\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{j}=0$
; see e.g. Helander (Reference Helander2014).
The basic strategy to find an equilibrium is to assert
$\boldsymbol{B}_{\text{con}}=\boldsymbol{B}_{\text{cov}}$
together with force balance (2.3), relying on the fact that these forms of the magnetic field incorporate (2.1) and (2.2) as well as the assumption of magnetic flux surfaces. Either the magnetic coordinates
$\unicode[STIX]{x1D713}$
,
$\unicode[STIX]{x1D703}$
and
$\unicode[STIX]{x1D711}$
can be considered as unknown functions of spatial coordinates (‘direct formulation’), or the coordinate mapping
$\boldsymbol{x}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$
can be considered as an unknown function of magnetic coordinates (‘inverse formulation’). The direct formulation is appropriate for the axisymmetric problem, where symmetry is defined relative to the geometric toroidal angle, while the inverse formulation is appropriate for higher orders in the expansion, where axisymmetry is broken and the condition of QAS can be enforced as a symmetry in the Boozer toroidal angle
$\unicode[STIX]{x1D711}$
.
In what follows we assume vacuum conditions,
$\text{d}G/\text{d}\unicode[STIX]{x1D713}=0$
, and
$I=K=0$
. The vacuum case has the curious feature that flux surfaces are not uniquely defined for the zeroth-order axisymmetric solution, since the field has no rotational transform (actually, we will find the rotational transform in our expansion is also zero at first order). However, the non-axisymmetric perturbation introduces a non-zero rotational transform, which then makes the surfaces unique.
3 Quasi-axisymmetry near axisymmetry
In the following, to more easily apply the condition of QAS, we employ the ‘inverse’ formulation, where unknown of the theory is the coordinate mapping
$\boldsymbol{x}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$
, and the magnetostatic equations are written in terms of various derivatives, denoted as
$\boldsymbol{e}_{1}=\unicode[STIX]{x2202}\boldsymbol{x}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}$
,
$\boldsymbol{e}_{2}=\unicode[STIX]{x2202}\boldsymbol{x}/\unicode[STIX]{x2202}\unicode[STIX]{x1D703}$
and
$\boldsymbol{e}_{3}=\unicode[STIX]{x2202}\boldsymbol{x}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}$
. These equations can be translated into equations involving the metrics via the usual identities (reviewed in appendix A). The equation
$\boldsymbol{B}_{\text{con}}=\boldsymbol{B}_{\text{cov}}$
yields what we will call the magnetostatic equation
From this, a ‘local’ magnetostatic constraint (Skovoroda Reference Skovoroda2007) (involving only derivatives within the surface) can be constructed:
where we define
$\unicode[STIX]{x1D628}_{ij}\equiv \boldsymbol{e}_{i}\boldsymbol{\cdot }\boldsymbol{e}_{j}$
. Alternatively, this follows directly from
$(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \boldsymbol{B}_{\text{cov}})\boldsymbol{\cdot }\boldsymbol{B}_{\text{con}}=0$
. A noteworthy consequence of this constraint is discussed in appendix E.
3.1 The condition of quasi-axisymmetry (QAS)
Quasi-symmetry can be stated as the condition that the magnetic field strength is a function of only a single helicity of the Boozer angles (Nührenberg & Zille Reference Nührenberg and Zille1988), i.e.
$B=B(\unicode[STIX]{x1D713},M\unicode[STIX]{x1D703}-N\unicode[STIX]{x1D711})$
where
$M$
and
$N$
are integers. QAS is the
$N=0$
case, i.e. simply the condition that the magnetic field strength is independent of the Boozer angle
$\unicode[STIX]{x1D711}$
, i.e.
$\unicode[STIX]{x2202}B/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}=0$
. Defining
$E=G^{2}/B^{2}$
and using (B 3), we can express the quantity
$E$
in terms of the surface metrics:
We can then restate QAS as
$\unicode[STIX]{x2202}E/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}=0$
, i.e.
3.2 The expansion about axisymmetry
We write the coordinate mapping
$\boldsymbol{x}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$
as a series expansion in the small parameter
$\unicode[STIX]{x1D716}$
,
$$\begin{eqnarray}\boldsymbol{x}=\boldsymbol{x}_{0}+\unicode[STIX]{x1D716}\boldsymbol{x}_{1}+\unicode[STIX]{x1D716}^{2}\boldsymbol{x}_{2}+\cdots \,,\end{eqnarray}$$
where
$\boldsymbol{x}_{0}$
corresponds to the zeroth-order axisymmetric magnetic field. For simplicity, we will consider the current
$G$
as fixed to its zeroth-order values (there is no loss of generality as these functions can always be adjusted at zeroth order). We however allow the deformation to modify 
, since we would like todetermine how much ‘external’ rotational transform can be achieved by the 3-D shaping.
where we have used the fact that a vacuum axisymmetric magnetic field has no rotational transform 
. We also introduce thefollowing notation for expanding the metrics:
$$\begin{eqnarray}\unicode[STIX]{x1D628}_{ij}=\unicode[STIX]{x1D628}_{ij}^{(0)}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D628}_{ij}^{(1)}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D628}_{ij}^{(2)}+\cdots \,.\end{eqnarray}$$
3.2.1
$O(\unicode[STIX]{x1D716}^{0})$
The coordinate mapping at zeroth order is written
$$\begin{eqnarray}\boldsymbol{x}_{0}=\hat{\boldsymbol{R}}R_{0}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})+\hat{\boldsymbol{z}}Z_{0}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713}).\end{eqnarray}$$
Note that the unit vectors
$\hat{\boldsymbol{\unicode[STIX]{x1D719}}}$
,
$\hat{\boldsymbol{R}}$
and
$\hat{\boldsymbol{z}}$
are defined with respect to the Boozer angle
$\unicode[STIX]{x1D711}$
, (this choice will simplify the vector algebra at higher order)
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{\boldsymbol{R}}=\hat{\boldsymbol{x}}\cos (\unicode[STIX]{x1D711})+\hat{\boldsymbol{y}}\sin (\unicode[STIX]{x1D711}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \hat{\boldsymbol{\unicode[STIX]{x1D719}}}=\hat{\boldsymbol{y}}\cos (\unicode[STIX]{x1D711})-\hat{\boldsymbol{x}}\sin (\unicode[STIX]{x1D711}). & \displaystyle\end{eqnarray}$$
We emphasize that these are not the usual cylindrical basis vectors, as the Boozer toroidal angle only corresponds to the geometric toroidal angle at zeroth order; see discussion at the end of § C.1. Substituting (3.8) into (3.1), and projecting along the
$\hat{\boldsymbol{\unicode[STIX]{x1D719}}}$
,
$\hat{\boldsymbol{R}}$
and
$\hat{\boldsymbol{z}}$
directions, yields only a single non-trivial constraint
$$\begin{eqnarray}G\left\{Z_{0},R_{0}\right\}=R_{0},\end{eqnarray}$$
where we define
$$\begin{eqnarray}\left\{A,B\right\}\equiv \frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}-\frac{\unicode[STIX]{x2202}B}{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}}\frac{\unicode[STIX]{x2202}A}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}.\end{eqnarray}$$
Of course, QAS is satisfied at this order,
$\unicode[STIX]{x2202}E_{0}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}=0$
, as
$$\begin{eqnarray}E_{0}=R_{0}^{2}.\end{eqnarray}$$
Note that (3.11) does not constrain the choice of the zeroth-order magnetic field very much. In fact, we may freely choose not only the shape of the outer flux surface, but the shape of all the interior surfaces as well; see appendix C.2. If the reader is bothered by the fact that such surfaces are not proper flux surfaces, in the sense that, because 
, they are notunique, we remark that such non-uniqueness does not cause any mathematical inconsistencies at this order (the contravariant form of
$\boldsymbol{B}$
still correctly describes the field), and the uniqueness property will indeed by attained at higher order where the rotational transform is non-zero.
3.2.2
$O(\unicode[STIX]{x1D716}^{1})$
So as to preserve the unit vectors, as mentioned above, we do not perturb the geometric angle
$\unicode[STIX]{x1D719}$
. Instead, we include a (third) component of
$\boldsymbol{x}_{1}$
in the
$\hat{\boldsymbol{\unicode[STIX]{x1D719}}}$
direction, i.e. we write
$$\begin{eqnarray}\boldsymbol{x}_{1}=\hat{\boldsymbol{R}}R_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713},\unicode[STIX]{x1D711})+\hat{\boldsymbol{z}}Z_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713},\unicode[STIX]{x1D711})+\hat{\boldsymbol{\unicode[STIX]{x1D719}}}\unicode[STIX]{x1D6F7}_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713},\unicode[STIX]{x1D711}).\end{eqnarray}$$
As
$\unicode[STIX]{x1D711}$
is an ignorable coordinate in the first-order magnetostatic equations, we will assume
$$\begin{eqnarray}\displaystyle & \displaystyle R_{1}=\bar{R}_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})+2\text{Re}[\hat{R}_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})\exp (\text{i}N\unicode[STIX]{x1D711})], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Z_{1}=\bar{Z}_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})+2\text{Re}[\hat{Z}_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})\exp (\text{i}N\unicode[STIX]{x1D711})], & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D6F7}_{1}=\bar{\unicode[STIX]{x1D6F7}}_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})+2\text{Re}[\hat{\unicode[STIX]{x1D6F7}}_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})\exp (\text{i}N\unicode[STIX]{x1D711})], & \displaystyle\end{eqnarray}$$
where
$N\neq 0$
is an integer, which will be interpreted as the field period number, since, as we will see, higher-order corrections will only involve harmonics of this mode number. Note that this is not the only way we can construct an
$N$
-period stellarator, as additional harmonics (
$2N$
,
$3N$
, etc.) could be included already at this order, but we make the above choice for simplicity. The
$\unicode[STIX]{x1D711}$
-averaged part of the local magnetostatic constraint ((3.2), 
) yields
Periodicity of
$\bar{\unicode[STIX]{x1D6F7}}_{1}/R_{0}$
yields the solubility constraint
implying 
. Then, from (3.18), we find
$\bar{\unicode[STIX]{x1D6F7}}_{1}=C_{1}R_{0}$
, where
$C_{1}(\unicode[STIX]{x1D713})$
is an arbitrary constant. The only remaining useful information that can be obtained from the
$\unicode[STIX]{x1D711}$
-average of (3.1) comes from its
$\hat{\boldsymbol{\unicode[STIX]{x1D719}}}$
component, which gives a relationship between
$\bar{R}_{1}$
and
$\bar{Z}_{1}$
:
$$\begin{eqnarray}\bar{R}_{1}=G\left(\left\{\bar{Z}_{1},R_{0}\right\}+\left\{Z_{0},\bar{R}_{1}\right\}\right).\end{eqnarray}$$
We are thus free to set
$\bar{\unicode[STIX]{x1D6F7}}_{1}=\bar{R}_{1}=\bar{Z}_{1}=0$
at this point; note the axisymmetric part of the first-order solution could simply be absorbed into the zeroth-order solution, and therefore does not represent interesting freedom in the QAS solution. Furthermore, these components do not affect what follows, so there is no loss of generality. For the
$\unicode[STIX]{x1D711}$
-dependent part of the first-order solution, we take components of the magnetostatic equation (3.1) in the
$\hat{\boldsymbol{\unicode[STIX]{x1D719}}}$
,
$\hat{\boldsymbol{R}}$
and
$\hat{\boldsymbol{z}}$
directions to obtain
$$\begin{eqnarray}\displaystyle & \displaystyle \text{i}N\hat{\unicode[STIX]{x1D6F7}}_{1}+\hat{R}_{1}=G\left(\left\{\hat{Z}_{1},R_{0}\right\}+\left\{Z_{0},\hat{R}_{1}\right\}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{i}N\hat{R}_{1}-\hat{\unicode[STIX]{x1D6F7}}_{1}=G\left\{\hat{\unicode[STIX]{x1D6F7}}_{1},Z_{0}\right\}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{i}N\hat{Z}_{1}=G\left\{R_{0},\hat{\unicode[STIX]{x1D6F7}}_{1}\right\}. & \displaystyle\end{eqnarray}$$
A single second-order equation for
$\hat{\unicode[STIX]{x1D6F7}}_{1}$
can be obtained from this system by substituting (3.22) and (3.23) into (3.21).
$$\begin{eqnarray}(N^{2}-1)\hat{\unicode[STIX]{x1D6F7}}_{1}+2G\left\{Z_{0},\hat{\unicode[STIX]{x1D6F7}}_{1}\right\}=G^{2}\left(\left\{R_{0},\left\{R_{0},\hat{\unicode[STIX]{x1D6F7}}_{1}\right\}\right\}+\left\{Z_{0},\left\{Z_{0},\hat{\unicode[STIX]{x1D6F7}}_{1}\right\}\right\}\right).\end{eqnarray}$$
We are further able to simplify the resulting equation significantly by changing coordinates from (
$\unicode[STIX]{x1D713}$
,
$\unicode[STIX]{x1D703}$
) to (
$R_{0}$
,
$Z_{0}$
). Defining
$$\begin{eqnarray}P_{1}(R_{0},Z_{0})=\hat{\unicode[STIX]{x1D6F7}}_{1}(\unicode[STIX]{x1D713}(R_{0},Z_{0}),\unicode[STIX]{x1D703}(R_{0},Z_{0})),\end{eqnarray}$$
and using (3.11) (see also appendix D), we obtain finally the simple equation
$$\begin{eqnarray}(N^{2}-1)P_{1}=R_{0}^{2}\unicode[STIX]{x1D6E5}_{0}^{\ast }P_{1},\end{eqnarray}$$
where we encounter the familiar Grad–Shafranov operator
$$\begin{eqnarray}\unicode[STIX]{x1D6E5}_{0}^{\ast }=R_{0}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}R_{0}}\left(\frac{1}{R_{0}}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}R_{0}}\right)+\frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}Z_{0}^{2}}.\end{eqnarray}$$
The QAS condition is enforced by setting
$\unicode[STIX]{x2202}E_{1}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}=0$
. Using
$E_{1}=\unicode[STIX]{x1D628}_{33}^{(1)}$
we obtain
$$\begin{eqnarray}\text{i}NR_{0}(\text{i}N\hat{\unicode[STIX]{x1D6F7}}_{1}+\hat{R}_{1})=0,\end{eqnarray}$$
i.e.
$\text{i}N\hat{\unicode[STIX]{x1D6F7}}_{1}+\hat{R}_{1}=0$
. Using (3.22) we obtain from this a condition on
$P_{1}$
,
$$\begin{eqnarray}(N^{2}-1)P_{1}+R_{0}\frac{\unicode[STIX]{x2202}P_{1}}{\unicode[STIX]{x2202}R_{0}}=0.\end{eqnarray}$$
It can be easily verified that (except for a trivial
$N=1$
case corresponding to rotation; see appendix F) the only way to satisfy this condition across the entire plasma volume is to have
$P_{1}=0$
. Global QAS is thereby proved impossible for vacuum fields, supporting the conclusions of Garren & Boozer (Reference Garren and Boozer1991a
). If, however, we require (3.29) to be satisfied on a single surface, denoted
$\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}_{bc}$
, then (3.29) constitutes an ‘oblique’ (or more specifically, ‘tangential–oblique’) boundary condition for the elliptic homogeneous second-order equation (3.26). This problem was first studied by Poincaré (Poincaré Reference Poincaré1910), and is therefore sometimes called ‘Poincaré’s problem’. The chief difficulty in the problem is due to the fact that the direction of the derivative is tangential to the boundary surface
$\unicode[STIX]{x1D713}_{bc}$
at a discrete set of points. In the present case, because (3.26) is homogeneous, and the direction of the derivative does not rotate as the boundary curve is traced, it is known that exactly two linearly independent solutions exist (Vekua Reference Vekua1953; see also ch. III, § 23 of Miranda Reference Miranda1970), aside from the trivial null solution
$P_{1}=0$
.
For completeness we state the equations for
$\hat{R}_{1}$
and
$\hat{Z}_{1}$
in terms of
$P_{1}$
. From (3.22) to (3.23) we have
$$\begin{eqnarray}\displaystyle & \displaystyle \text{i}N\hat{R}_{1}=P_{1}-R_{0}\frac{\unicode[STIX]{x2202}P_{1}}{\unicode[STIX]{x2202}R_{0}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{i}N\hat{Z}_{1}=-R_{0}\frac{\unicode[STIX]{x2202}P_{1}}{\unicode[STIX]{x2202}Z_{0}}. & \displaystyle\end{eqnarray}$$
3.2.3
$O(\unicode[STIX]{x1D716}^{2})$
We note that, at second order, the finite-toroidal-number components of the magnetostatic equation will generally be non-zero for
$n=\pm 2N$
due to source terms that are quadratic in first-order quantities. We are free to choose the other components to be zero, and write the coordinate mapping in the same form as at first order, i.e.
$\boldsymbol{x}_{2}=\hat{\boldsymbol{R}}R_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713},\unicode[STIX]{x1D711})+\hat{\boldsymbol{z}}Z_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713},\unicode[STIX]{x1D711})+\hat{\boldsymbol{\unicode[STIX]{x1D719}}}\unicode[STIX]{x1D6F7}_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713},\unicode[STIX]{x1D711})$
, with
$$\begin{eqnarray}\displaystyle R_{2}=\bar{R}_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})+2\text{Re}[\hat{R}_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})\exp (2\text{i}N\unicode[STIX]{x1D711})], & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle Z_{2}=\bar{Z}_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})+2\text{Re}[\hat{Z}_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})\exp (2\text{i}N\unicode[STIX]{x1D711})], & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D6F7}_{2}=\bar{\unicode[STIX]{x1D6F7}}_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})+2\text{Re}[\hat{\unicode[STIX]{x1D6F7}}_{2}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})\exp (2\text{i}N\unicode[STIX]{x1D711})]. & & \displaystyle\end{eqnarray}$$
To find 
we can use the local magnetostatic constraint (3.2) at second order, i.e. 
. Taking the
$\unicode[STIX]{x1D711}$
-average, we obtain
Integrating over
$\unicode[STIX]{x1D703}$
(at constant
$\unicode[STIX]{x1D713}$
) we obtain 
as a solubility constraint for
$\bar{\unicode[STIX]{x1D6F7}}_{2}$
:
Following essentially the same procedure used at first order, we can obtain a single elliptic partial differential equation for
$P_{2}(R_{0},Z_{0})=\hat{\unicode[STIX]{x1D6F7}}_{2}(\unicode[STIX]{x1D713}(R_{0},Z_{0}),\unicode[STIX]{x1D703}(R_{0},Z_{0}))$
,
$$\begin{eqnarray}\displaystyle & & \displaystyle (4N^{2}-1)P_{2}-R_{0}^{2}\unicode[STIX]{x1D6E5}_{0}^{\ast }P_{2}\nonumber\\ \displaystyle & & \displaystyle \quad =2\text{i}NR_{0}\left\{\hat{Z}_{1},\hat{R}_{1}\right\}_{R_{0}Z_{0}}+R_{0}^{2}\left(\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}Z_{0}}\left\{\hat{R}_{1},\hat{\unicode[STIX]{x1D6F7}}_{1}\right\}_{R_{0}Z_{0}}+\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}R_{0}}\left\{\hat{\unicode[STIX]{x1D6F7}}_{1},\hat{Z}_{1}\right\}_{R_{0}Z_{0}}\right),\end{eqnarray}$$
where
$\left\{A,B\right\}_{R_{0}Z_{0}}=\unicode[STIX]{x2202}_{R_{0}}A\unicode[STIX]{x2202}_{Z_{0}}B-\unicode[STIX]{x2202}_{R_{0}}B\unicode[STIX]{x2202}_{Z_{0}}A$
. The condition of QAS is stated as
$$\begin{eqnarray}(2N^{2}-1)P_{2}+R_{0}\frac{\unicode[STIX]{x2202}P_{2}}{\unicode[STIX]{x2202}R_{0}}=-R_{0}\left\{\hat{\unicode[STIX]{x1D6F7}}_{1},\hat{Z}_{1}\right\}_{R_{0}Z_{0}}+\frac{\text{i}N}{R_{0}}\left[(\text{i}N\hat{R}_{1}-\hat{\unicode[STIX]{x1D6F7}}_{1})^{2}+(\text{i}N\hat{Z}_{1})^{2}\right].\end{eqnarray}$$
which is again an oblique boundary condition. Equation (3.37) with boundary condition (3.38) is solvable (Vekua Reference Vekua1953; Miranda Reference Miranda1970).
3.2.4
$O(\unicode[STIX]{x1D716}^{n})$
,
$n>2$
At higher order, there will be more equations to solve due to the nonlinear interaction of lower-order contributions. For instance, at third order, there will be non-zero
$n=-3N,-N,0,N,3N$
components due to the source terms arising from the product of first- and second-order quantities. However, these problems will be of the type found at second order, and therefore must have solutions. We conclude that QAS (at a single surface) can be satisfied to all orders in the expansion.
4 Numerical solution of quasi-axisymmetric magnetic fields at first order
We turn now to the numerical solution of (3.26), imposing the oblique boundary condition, equation (3.29), corresponding to the condition of QAS. Note that, upon solving for
$P_{1}$
, the functions
$\hat{R}_{1}$
and
$\hat{Z}_{1}$
can then be immediately calculated from (3.30)–(3.31) and, fixing
$\unicode[STIX]{x1D716}$
, the total mapping can be constructed as
$$\begin{eqnarray}\boldsymbol{x}\approx \hat{\boldsymbol{R}}(R_{0}+2\unicode[STIX]{x1D716}\text{Re}[\hat{R}_{1}\exp (\text{i}N\unicode[STIX]{x1D711})])+\hat{\boldsymbol{z}}(Z_{0}+2\unicode[STIX]{x1D716}\text{Re}[\hat{Z}_{1}\exp (\text{i}N\unicode[STIX]{x1D711})])+2\hat{\boldsymbol{\unicode[STIX]{x1D719}}}\unicode[STIX]{x1D716}\text{Re}[\hat{\unicode[STIX]{x1D6F7}}_{1}\exp (\text{i}N\unicode[STIX]{x1D711})].\end{eqnarray}$$
Due to the change of variables to (
$R_{0}$
,
$Z_{0}$
) space, the only knowledge of the zeroth-order solution needed is the shape of the outer flux surface, which defines the computational domain. If the first-order deformation is specified, along with the amplitude
$\unicode[STIX]{x1D716}$
, then
$\boldsymbol{x}$
is determined by (4.1), in the coordinates
$R_{0}$
,
$Z_{0}$
and
$\unicode[STIX]{x1D711}$
. We will henceforth drop the subscripts and work directly with the following equations for
$P$
in the
$R$
–
$Z$
plane:
$$\begin{eqnarray}\displaystyle & \displaystyle (N^{2}-1)P+3R\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}R}=(\hat{\boldsymbol{R}}\unicode[STIX]{x2202}_{R}+\hat{\boldsymbol{z}}\unicode[STIX]{x2202}_{Z})\boldsymbol{\cdot }\left(\hat{\boldsymbol{R}}R^{2}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}R}+\hat{\boldsymbol{z}}R^{2}\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}Z}\right), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle (N^{2}-1)P+R\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}R}=0. & \displaystyle\end{eqnarray}$$
Note that (4.2) is simply (3.29), with the diffusive term rewritten in conservative form.
The numerical method used to solve (4.2) with boundary condition (4.3) is based on a standard finite element method approach to a Neumann boundary value problem. As usual, the function
$P$
is expressed as a sum of basis functions that have finite support over mesh elements. Multiplying (4.2) by a basis function and integrating over the domain, one obtains a matrix equation containing the so-called stiffness matrix. Part of the stiffness matrix is due to the flux at the boundary (arising from the conservative term on the right-hand side of (4.2)), and involves the normal derivative of
$P$
. This flux term is evaluated using our boundary condition, equation (4.3). That is, the derivative term in (4.3) is rewritten in terms of components that are tangential and normal to the boundary:
$$\begin{eqnarray}(N^{2}-1)P+\unicode[STIX]{x1D6FE}(l)\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}n}+\unicode[STIX]{x1D6FD}(l)\frac{\unicode[STIX]{x2202}P}{\unicode[STIX]{x2202}l}=0.\end{eqnarray}$$
The functions
$\unicode[STIX]{x1D6FE}$
and
$\unicode[STIX]{x1D6FD}$
depend on the shape of the boundary. Equation (4.4) is divided by
$\unicode[STIX]{x1D6FE}$
to obtain an expression for the normal derivative, and the resulting equation is discretized over a one-dimensional mesh whose nodes match the boundary nodes of the finite element mesh. The result is used to evaluate the relevant contribution the stiffness matrix obtained from the discretization of (4.2).
The only difficulty in this procedure arises from the fact that the coefficient of the normal derivative,
$\unicode[STIX]{x1D6FE}(l)$
, passes through zero at a discrete set of boundary points (those where the normal vector points in the
$\hat{\boldsymbol{z}}$
direction) and so the flux is not determined by (4.4) at those points. To overcome this, we only need to avoid these points when choosing the nodes of the mesh that lie on the boundary.
The final result of the above procedure is a homogeneous matrix equation. The two eigenfunctions with eigenvalues that tend toward zero (with increasing resolution) are identified as the solutions. The solutions can (due to their linear independence) be chosen such that one is even and the other is odd in a suitably defined geometric angle
$\unicode[STIX]{x1D703}$
. A sum of the two solutions can be taken, and the relative phase and amplitude can be adjusted to maximize the rotational transform.
Let us illustrate the above procedure with the simple example of a circular boundary centred at
$R=R_{c}$
with radius
$1$
. This case yields
$\unicode[STIX]{x1D6FE}(l)=[R_{c}+\cos (l)]\cos (l)$
, and
$\unicode[STIX]{x1D6FD}(l)=-[R_{c}+\cos (l)]\sin (l)$
. The resulting mesh and two independent solutions are depicted in figure 1.
Figure 1. Solution for circular boundary with
$N=2$
chosen. (a) A visualization of a circular second-order mesh, with boundary nodes indicated with blue dots, and the two points of tangency indicated with red dots. The numerical solution, including mesh generation, is done with the recently introduced finite element method capabilities of Mathematica 10. (b,c) Two independent solutions found for
$P(R_{0},Z_{0})$
. Note that due to the up–down mirror symmetry of the solution, the two solutions can be chosen as even (b) and odd (c).
The even and odd solutions are superposed to create a total solution
$P=AP_{e}+BP_{o}$
, where
$A$
and
$B$
are arbitrary complex constants. To obtain finite rotational transform for the present example, it is necessary to assume a non-zero phase shift between these amplitudes. Taking
$A=1$
and
$B=\exp (-\text{i}\unicode[STIX]{x03C0}/2)$
and fixing
$\unicode[STIX]{x1D716}$
the surface shape, according to (4.1), is plotted in figure 2.
Figure 2. Visualization of boundary surface shape up to first order, assuming circular zeroth-order shape (same case as figure 1). Here
$N=2$
and
$\unicode[STIX]{x1D716}=2.0$
; note that the eigenfunctions are normalized to a small numerical absolute value, as shown in figure 1, so the total deformation is small despite the fact that
$\unicode[STIX]{x1D716}$
is not. The toroidal cuts and lines of constant poloidal and toroidal angle are included for purely stylistic reasons.
Figure 3. Demonstration that QAS is satisfied to first order in size of non-axisymmetric perturbation; the measure of deviation
$Q$
, defined by (4.5), is computed only at the surface where the QAS condition is applied. Here the field period number is
$N=2$
and the zeroth-order flux-surface shape is circular with an aspect ratio of
$4$
. Solution of (4.2) with boundary condition (4.3) leads to the expected scaling with
$\unicode[STIX]{x1D716}^{2}$
, as shown in (a). In (b), a ‘control’ (non-QAS) deformation is used for comparison (
$\unicode[STIX]{x1D6F7}_{1}=0$
,
$R_{1}=0.4\cos (\unicode[STIX]{x1D717}-2\unicode[STIX]{x1D711})$
and
$Z_{1}=0.4\sin (\unicode[STIX]{x1D717}-2\unicode[STIX]{x1D711})$
, where
$\unicode[STIX]{x1D717}$
is a geometric poloidal angle), and it is found that QAS is broken at first order in
$\unicode[STIX]{x1D716}$
, as expected. Visualization of surface shape (for largest value of
$\unicode[STIX]{x1D716}$
) is included within each plot. (Note that
$\unicode[STIX]{x1D716}$
is somewhat arbitrary as the solutions are linear and subject to arbitrary normalization.)
To independently confirm the result, the surface shape, determined by (4.1), can be used as input for the VMEC code. This shape specification has a free parameter
$\unicode[STIX]{x1D716}$
, which can be varied to test the degree to which the result satisfies the QAS condition. Since the solution is only correct to first order, the ‘exact’ solution obtained from VMEC will be expected to deviate from quasi-axisymmetry at second order, on the appropriate flux surface. Therefore, we conduct a series of runs with increasing strength of perturbation,
$\unicode[STIX]{x1D716}$
. The resulting solution is then passed to a separate code ‘BOOZ_XFORM’ (Sanchez et al.
Reference Sanchez, Hirshman, Ware, Berry and Spong2000), which computes the components of
$|\boldsymbol{B}|$
in Boozer angles separately on each flux surface. Examining these components for the outer surface, we find that the ratio of the non-QAS part to the full field, does indeed scale as
$\unicode[STIX]{x1D716}^{2}$
. For comparison, the scaling of a control (non-QAS, non-axisymmetric) deformation is also computed, and found to be linear in
$\unicode[STIX]{x1D716}$
as expected; see figure 3. In this figure the quantity
$Q$
measures the deviation from QAS, and is defined in terms of the Fourier component of the magnetic field magnitude expressed in Boozer coordinates
$\hat{B}_{mn}$
, where
$m$
and
$n$
are the poloidal and toroidal mode numbers respectively:
$$\begin{eqnarray}Q=\frac{\displaystyle \left(\mathop{\sum }_{m,n\neq 0}|\hat{B}_{mn}|^{2}\right)^{1/2}}{\left(\displaystyle \mathop{\sum }_{m,n}|\hat{B}_{mn}|^{2}\right)^{1/2}}.\end{eqnarray}$$
5 Omnigeneity
Omnigeneity (Hall & McNamara Reference Hall and McNamara1975) is the more general condition of having zero average radial particle drift. This can be shown to be equivalent to the geometric condition that points of equal magnetic field strength have constant separation in Boozer angles, i.e. the separation (
$\unicode[STIX]{x0394}\unicode[STIX]{x1D703}$
,
$\unicode[STIX]{x0394}\unicode[STIX]{x1D711}$
) is independent of the field line chosen (see also Landreman & Catto (Reference Landreman and Catto2012) for further discussion of this). It has been argued that exactly omnigenous configurations, must be quasi-symmetric if they are to be analytic (Cary & Shasharina Reference Cary and Shasharina1997). The same authors point out, however, that analytic fields can be specified that are arbitrarily close to being omnigenous, while also being far from quasi-symmetric. One can conclude that omnigeneity is a useful target for optimization. Indeed, it is fair to say, roughly speaking, that nearly quasi-symmetric configurations represent a small subspace of nearly omnigenous configurations. Therefore, it may be useful to be able to impose omnigeneity instead of quasi-symmetry as a boundary condition in our expansion near axisymmetry. As will be shown, we can impose an arbitrary modification to the magnetic field strength,
$E_{1}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$
. This function can, in particular, be chosen such that
$E_{0}+\unicode[STIX]{x1D716}E_{1}$
satisfies the condition of omnigeneity. We write
$E_{1}$
as a Fourier series,
$$\begin{eqnarray}E_{1}=\bar{E}_{1}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703})+\mathop{\sum }_{n\neq 0}\hat{E}_{1}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},n)\exp (\text{i}n\unicode[STIX]{x1D711}),\end{eqnarray}$$
and first consider the
$\unicode[STIX]{x1D711}$
-independent part,
$\bar{E}_{1}$
. Note that the magnetostatic equation, (3.1), was already solved for the QAS problem, and yields
$\bar{\unicode[STIX]{x1D6F7}}_{1}=C_{1}R_{0}$
and 
. From (3.3), we then find that
$\bar{E}_{1}=2R_{0}\bar{R}_{1}$
, i.e.
$$\begin{eqnarray}\bar{R}_{1}=\frac{\bar{E}_{1}}{2R_{0}}.\end{eqnarray}$$
We now can rewrite (3.20), which relates
$\bar{Z}_{1}$
to
$\bar{R}_{1}$
, using (3.11), in terms of independent variables
$Z_{0}$
and
$R_{0}$
as follows
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\bar{Z}_{1}}{\unicode[STIX]{x2202}Z_{0}}=\bar{R}_{1}-\frac{\unicode[STIX]{x2202}\bar{R}_{1}}{\unicode[STIX]{x2202}R_{0}}.\end{eqnarray}$$
Thus
$\bar{Z}_{1}$
is determined by
$\bar{R}_{1}$
via (5.3), completing the
$\unicode[STIX]{x1D711}$
-independent part of the solution. For the
$\unicode[STIX]{x1D711}$
-dependent part of the solution, we essentially only need to add a non-zero contribution to the right-hand side of (3.29), noting, however, that the resulting equation must be satisfied for all toroidal mode numbers. We therefore modify the notation
$N\rightarrow n$
, and obtain the following PDE for
$P_{1}(R_{0},Z_{0},n)$
, with
$n\neq 0$
:
$$\begin{eqnarray}(n^{2}-1)P_{1}=R_{0}^{2}\unicode[STIX]{x1D6E5}_{0}^{\ast }P_{1}.\end{eqnarray}$$
The boundary condition filling the role of (3.29), now with an inhomogeneity proportional to
$\hat{E}_{1}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},n)$
, becomes (for
$n\neq 0$
):
$$\begin{eqnarray}R_{0}(n^{2}-1)P_{1}+R_{0}^{2}\frac{\unicode[STIX]{x2202}P_{1}}{\unicode[STIX]{x2202}R_{0}}=-\text{i}n\hat{E}_{1}.\end{eqnarray}$$
We note first, that the existence of at least one solution is guaranteed for this non-homogeneous oblique derivative problem (Vekua Reference Vekua1953; see also ch. III, § 23 of Miranda Reference Miranda1970). This problem is more difficult to solve than the QAS problem, since the equations for all necessary Fourier components of
$E_{1}$
must be solved. Furthermore, it seems necessary to solve for the
$\unicode[STIX]{x1D703}(R_{0},Z_{0})$
to be able to construct a function
$E_{1}$
such that the total field is omnigenous; this was not necessary for the QAS problem since only the shape of the zeroth-order boundary flux surface was needed to solve the first-order system.
We may draw an interesting conclusion at this point, namely that the QAS problems, which can be indexed by
$n\neq 0$
, are valid homogeneous solutions of the omnigenous problem, i.e. the corresponding deformations do not modify the magnitude of the magnetic field (by construction), so they can be added to a solution of the general first-order problem (5.4) and (5.5), and the result will still satisfy (5.4) across the volume, and (5.5) on the boundary. The QAS deformations therefore represent an arbitrary freedom for omnigenous solutions.
6 Conclusions
One of the perhaps surprising outcomes of this work is the abundance of QAS solutions. Instead of being isolated, they form a continuum in the vicinity of axisymmetric solutions. The parameter space for this continuum consists of the (i) the space of axisymmetric flux-surface shape functions, (ii) the discrete parameter
$N$
, defining the stellarator field period number and (iii) the two complex amplitudes of the solutions of the oblique derivative problem. From one perspective, the size of this space might be expected: the axisymmetric part of our solution satisfies QAS in a trivial sense, and globally. So much of the freedom in the solution corresponds to the part of the solution that is known to be exceptional. If, as is true for the near-axis expansion of (Garren & Boozer Reference Garren and Boozer1991a
), the results found in the vacuum case are also valid for the non-vacuum case, a striking and simple conclusion could be drawn – namely, that a large space of QAS stellarator configurations are accessible via suitably small deformations of any tokamak.
We reach another interesting conclusion, following the discussion of § 5, namely that an alternative problem specification is possible, by fixing the perturbation in magnetic field strength at the boundary flux surface, as a function of Boozer angles. A solution is guaranteed to exist for this problem, and we note that an arbitrary QAS deformation may be added to this solution, such that the total solution still satisfies the field strength boundary condition. Furthermore, the field strength may be chosen such that a generalized form of QAS (omnigeneity) is satisfied, significantly broadening the space of ‘good’ configurations that can be found in the neighbourhood of axisymmetry. Although a deformation can only satisfy exact QAS at one surface, it remains an open question whether omnigeneity can be satisfied throughout the volume.
As a final note, it is worth remarking on how our perturbed solutions, which can be thought of as nearly two-dimensional, relate to fully 3-D magnetic fields, which are known to suffer from singularities associated with rational rotational transform that can break up flux surfaces. Since our treatment begins with the assumption of flux surfaces, it is fair to ask to what extent this assumption could be broken in actual configurations unconstrained by this assumption. However, we can argue that this will not be a problem, for at least sufficiently small deviations from axisymmetry. This is because the obtained rotational transform is necessarily small, scaling as
$\unicode[STIX]{x1D716}^{2}$
, and so the rational surfaces must be very high order, i.e. 
with
$m\gg 1$
. Such surfaces, as argued by Rosenbluth et al. (Reference Rosenbluth, Sagdeev, Taylor and Zaslavski1966), give rise to island chains that have exponentially small width.
Acknowledgements
Thanks to S. Henneberg for helpful conversations, and S. Lazerson and J. Geiger for assistance with the VMEC code. This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 under grant agreement no. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.
Appendix A. Magnetic geometry
Given arbitrary toroidal coordinates (
$\unicode[STIX]{x1D713}$
,
$\unicode[STIX]{x1D703}$
,
$\unicode[STIX]{x1D711}$
), whose handedness is
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703})>0$
, it is the convention to work with local basis vectors, (
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$
,
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}$
,
$\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}$
) and (
$\boldsymbol{e}_{1}=\unicode[STIX]{x2202}\boldsymbol{x}/\unicode[STIX]{x2202}\unicode[STIX]{x1D713}$
,
$\boldsymbol{e}_{2}=\unicode[STIX]{x2202}\boldsymbol{x}/\unicode[STIX]{x2202}\unicode[STIX]{x1D703}$
,
$\boldsymbol{e}_{3}=\unicode[STIX]{x2202}\boldsymbol{x}/\unicode[STIX]{x2202}\unicode[STIX]{x1D711}$
). The metric components are defined in the usual way
$$\begin{eqnarray}\unicode[STIX]{x1D628}_{ij}\equiv \boldsymbol{e}_{i}\boldsymbol{\cdot }\boldsymbol{e}_{j},\end{eqnarray}$$
and the Jacobian for these coordinates is
$$\begin{eqnarray}J=\frac{1}{\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D711})}=\boldsymbol{e}_{1}\boldsymbol{\cdot }(\boldsymbol{e}_{2}\times \boldsymbol{e}_{3}).\end{eqnarray}$$
Additionally, assigning
$(u_{1},u_{2},u_{3})=(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$
, we see that
$\boldsymbol{e}_{i}\boldsymbol{\cdot }\unicode[STIX]{x1D735}u_{j}=\unicode[STIX]{x1D6FF}_{ij}$
, and the following identities are easily verified
$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{e}_{1}=J(\unicode[STIX]{x1D735}u_{2}\times \unicode[STIX]{x1D735}u_{3}),\quad \boldsymbol{e}_{2}=J(\unicode[STIX]{x1D735}u_{3}\times \unicode[STIX]{x1D735}u_{1}),\quad \boldsymbol{e}_{3}=J(\unicode[STIX]{x1D735}u_{1}\times \unicode[STIX]{x1D735}u_{2}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}u_{1}=\frac{\boldsymbol{e}_{2}\times \boldsymbol{e}_{3}}{J},\quad \unicode[STIX]{x1D735}u_{2}=\frac{\boldsymbol{e}_{3}\times \boldsymbol{e}_{1}}{J},\quad \unicode[STIX]{x1D735}u_{3}=\frac{\boldsymbol{e}_{1}\times \boldsymbol{e}_{2}}{J}. & \displaystyle\end{eqnarray}$$
Appendix B. Useful forms of
$J$
and
$B$
Taking
$B^{2}=\boldsymbol{B}_{\text{cov}}\boldsymbol{\cdot }\boldsymbol{B}_{\text{con}}$
gives
Taking
$B^{2}=\boldsymbol{B}_{\text{con}}\boldsymbol{\cdot }\boldsymbol{B}_{\text{con}}$
gives
Using (B 1)–(B 2) we can express the magnetic field strength ‘locally’ (in terms of only surface metrics),
Using (B 1)–(B 2) we can also express the Jacobian locally,
Appendix C. General axisymmetric fields in Boozer coordinates
The condition of axisymmetry can be stated as the condition that the
$\hat{\boldsymbol{R}}$
,
$\hat{\unicode[STIX]{x1D753}}$
and
$\hat{\boldsymbol{z}}$
components of
$\boldsymbol{B}$
are independent of
$\unicode[STIX]{x1D719}$
. This implies that the flux surfaces must be themselves axisymmetric, so
$\hat{\unicode[STIX]{x1D753}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}=0$
, i.e.
$\unicode[STIX]{x2202}_{\unicode[STIX]{x1D719}}\unicode[STIX]{x1D713}=0$
. Then from
$\unicode[STIX]{x2202}_{\unicode[STIX]{x1D719}}(\hat{\unicode[STIX]{x1D753}}\boldsymbol{\cdot }\boldsymbol{B}_{\text{cov}})=0$
we obtain
$$\begin{eqnarray}G\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D711}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{2}}+I\frac{\unicode[STIX]{x2202}^{2}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}^{2}}=0.\end{eqnarray}$$
Integrating, using
$\unicode[STIX]{x1D711}=\unicode[STIX]{x1D719}+\unicode[STIX]{x1D708}$
, and periodicity in
$\unicode[STIX]{x1D719}$
we obtain
$$\begin{eqnarray}G\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}+I\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}=0.\end{eqnarray}$$
Now, taking the
$\unicode[STIX]{x2202}_{\unicode[STIX]{x1D719}}(\hat{\boldsymbol{R}}\boldsymbol{\cdot }B_{\text{con}})=0$
, we likewise obtain
Note that (C 2)–(C 3) are linearly independent (as guaranteed by (B 1)), so we can conclude that
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}=\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D708}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}=0.\end{eqnarray}$$
C.1 Axisymmetric fields in the direct formulation
Switching to the direct formulation, let us consider magnetic coordinates
$\unicode[STIX]{x1D713}$
,
$\unicode[STIX]{x1D703}$
,
$\unicode[STIX]{x1D711}$
as functions of cylindrical coordinates
$R$
,
$Z$
and
$\unicode[STIX]{x1D719}$
. Here
$Z$
is the distance along the axis of symmetry (of the zeroth-order solution),
$R$
is the distance from the origin in the plane perpendicular to
$Z$
and
$\unicode[STIX]{x1D719}$
is the geometric toroidal angle measuring the rotation about the
$Z$
axis. Writing
$\boldsymbol{B}_{\text{con}}=\boldsymbol{B}_{\text{cov}}$
we have
In a vacuum axisymmetric field we have 
and
$\text{d}G/\text{d}\unicode[STIX]{x1D713}=0$
, i.e.
$$\begin{eqnarray}G\unicode[STIX]{x1D735}\unicode[STIX]{x1D711}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}.\end{eqnarray}$$
By axisymmetry,
$\unicode[STIX]{x1D713}$
and
$\unicode[STIX]{x1D703}$
are independent of
$\unicode[STIX]{x1D719}$
(see appendix C), whereas
$\unicode[STIX]{x1D711}$
can generally be expressed as
$$\begin{eqnarray}\unicode[STIX]{x1D711}=\unicode[STIX]{x1D719}+\unicode[STIX]{x1D708}(R,Z).\end{eqnarray}$$
The right-hand side of (C 6) has only a toroidal component, so the
$\hat{\boldsymbol{R}}$
and
$\hat{\boldsymbol{z}}$
components of the equation require that
$\unicode[STIX]{x1D708}$
is a constant, which we set to zero. The
$\hat{\unicode[STIX]{x1D753}}$
component of (C 6) yields the only non-trivial constraint,
$$\begin{eqnarray}\frac{G}{R}=\left\{\unicode[STIX]{x1D703},\unicode[STIX]{x1D713}\right\}_{RZ},\end{eqnarray}$$
where
$\left\{A,B\right\}_{RZ}=\unicode[STIX]{x2202}_{R}A\unicode[STIX]{x2202}_{Z}B-\unicode[STIX]{x2202}_{R}B\unicode[STIX]{x2202}_{Z}A$
. To summarize,
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}(R,Z), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}(R,Z), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D711}=\unicode[STIX]{x1D719}. & \displaystyle\end{eqnarray}$$
We note that (C 11) can be interpreted (within the context of the expansion about axisymmetry) to mean that the zeroth-order geometric toroidal angle is equal to the Boozer toroidal angle, i.e.
$\unicode[STIX]{x1D719}_{0}=\unicode[STIX]{x1D711}$
; the geometric angle
$\unicode[STIX]{x1D719}$
can be defined exactly in terms of the coordinate mapping via the relation
$\tan (\unicode[STIX]{x1D719})=(\hat{\boldsymbol{y}}\boldsymbol{\cdot }\boldsymbol{x})/(\hat{\boldsymbol{x}}\boldsymbol{\cdot }\boldsymbol{x})$
. One can observe that this geometric angle will acquire corrections at higher order,
$\unicode[STIX]{x1D719}_{1}$
,
$\unicode[STIX]{x1D719}_{2}$
, etc., due to the modifications to the coordinate mapping. However, we can avoid solving for these corrections by defining the unit vectors
$\hat{\boldsymbol{R}}$
and
$\hat{\boldsymbol{\unicode[STIX]{x1D719}}}$
in terms of the Boozer angle
$\unicode[STIX]{x1D711}$
(3.9)–(3.10) and allowing the perturbation of the coordinate mapping to have a component in the
$\hat{\boldsymbol{\unicode[STIX]{x1D719}}}$
direction; see (3.14). We hope these conventions do not cause confusion, but they significantly simplify the vector algebra that follows, since it avoids unnecessary complications associated with expanding the conventionally defined unit vector,
$\hat{\boldsymbol{R}}(\unicode[STIX]{x1D719})=\hat{\boldsymbol{x}}\cos (\unicode[STIX]{x1D719})+\hat{\boldsymbol{y}}\sin (\unicode[STIX]{x1D719})$
.
C.2 Freedom at zeroth order according to (3.11)
Consider a set of nested arbitrarily shaped axisymmetric surfaces, labelled by the variable
$\unicode[STIX]{x1D70C}$
; let these correspond to constant-
$\unicode[STIX]{x1D713}$
surfaces. The function
$\unicode[STIX]{x1D713}(\unicode[STIX]{x1D70C})$
is to be determined, as follows. Let the variable
$l$
be the arc length parameterizing a contour of constant
$\unicode[STIX]{x1D70C}$
, such that
$(\hat{\unicode[STIX]{x1D753}}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D70C})\boldsymbol{\cdot }\unicode[STIX]{x1D735}l>0$
. Then (3.11) becomes
$$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D713}}{\text{d}\unicode[STIX]{x1D70C}}\frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}{\unicode[STIX]{x2202}l}=\frac{G}{R_{0}|\unicode[STIX]{x1D735}\unicode[STIX]{x1D70C}|}\equiv F(\unicode[STIX]{x1D70C},l)\end{eqnarray}$$
using
$\oint \unicode[STIX]{x2202}\unicode[STIX]{x1D703}/\unicode[STIX]{x2202}l\,\text{d}l=-2\unicode[STIX]{x03C0}$
we find
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D713}=\frac{1}{2\unicode[STIX]{x03C0}}\int _{0}^{\unicode[STIX]{x1D70C}}\text{d}\unicode[STIX]{x1D70C}^{\prime }\oint \text{d}l\,F(\unicode[STIX]{x1D70C}^{\prime },l), & & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle \unicode[STIX]{x1D703}=\unicode[STIX]{x1D703}_{0}(\unicode[STIX]{x1D70C})-2\unicode[STIX]{x03C0}\int _{0}^{l}\text{d}l^{\prime }\frac{F(\unicode[STIX]{x1D70C},\ell ^{\prime })}{\oint \text{d}l^{\prime \prime }F(\unicode[STIX]{x1D70C},l^{\prime \prime })}, & & \displaystyle\end{eqnarray}$$
where we use that the toroidal flux on axis
$\unicode[STIX]{x1D713}(\unicode[STIX]{x1D70C}=0)$
is zero. It is thus shown that the choice of an arbitrary set of nested axisymmetric surfaces is consistent with (3.11), and this choice determines the spatial dependence of the Boozer coordinates
$\unicode[STIX]{x1D713}$
and
$\unicode[STIX]{x1D703}$
up to an arbitrary angular shift
$\unicode[STIX]{x1D703}_{0}(\unicode[STIX]{x1D70C})$
, which can be inverted to find
$R_{0}$
and
$Z_{0}$
as functions of
$\unicode[STIX]{x1D713}$
and
$\unicode[STIX]{x1D703}$
.
Appendix D. Using
$R_{0}$
and
$Z_{0}$
as coordinates
The differential operators of (3.24) take the form of Poisson brackets involving the functions
$R_{0}$
and
$Z_{0}$
. Such operators can be interpreted as 2-D advection along contours of constant
$R_{0}$
and
$Z_{0}$
, which suggests that it could be fruitful to reformulate the equations, using
$R_{0}$
and
$Z_{0}$
as the independent variables, especially given that these variables are well-behaved coordinates. Furthermore, it will be useful, for the purposes of numerically solving the first-order system, to be able to reuse the domain defined already at zeroth order.
The following relations, found using (3.11), may make the derivation of (3.26) more transparent. For any function
$f$
, we have
$$\begin{eqnarray}\displaystyle & \displaystyle \left\{Z_{0},f\right\}=\frac{R_{0}}{G}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}R_{0}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \left\{R_{0},f\right\}=-\frac{R_{0}}{G}\frac{\unicode[STIX]{x2202}f}{\unicode[STIX]{x2202}Z_{0}}. & \displaystyle\end{eqnarray}$$
Appendix E. Consequence of the local magnetostatic constraint on the shape of QAS configurations
A differential constraint involving only in-surface derivatives can be derived directly from the condition
$\boldsymbol{B}_{\text{con}}\boldsymbol{\cdot }(\boldsymbol{B}_{\text{cov}}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713})=0$
. This, which we called the ‘local magnetostatic constraint’, is given in (3.2) and repeated here:
At first order in the expansion about axisymmetry, (E 1) yields for
$N\neq 0$
$$\begin{eqnarray}R_{0}\frac{\unicode[STIX]{x2202}\hat{\unicode[STIX]{x1D6F7}}_{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}+\frac{\unicode[STIX]{x2202}R_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}(\text{i}N\hat{R}_{1}-\hat{\unicode[STIX]{x1D6F7}}_{1})+\text{i}N\frac{\unicode[STIX]{x2202}Z_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\hat{Z}_{1}=0.\end{eqnarray}$$
This equation, combined with the condition of QAS, (3.28), gives
$$\begin{eqnarray}\frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}(\hat{\unicode[STIX]{x1D6F7}}_{1}R_{0}^{N^{2}-1})=-\text{i}N\frac{\unicode[STIX]{x2202}Z_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}R_{0}^{N^{2}-1}\hat{Z}_{1},\end{eqnarray}$$
which can be integrated over
$\unicode[STIX]{x1D703}$
to give the solubility constraint
$$\begin{eqnarray}\oint \text{d}\unicode[STIX]{x1D703}\frac{\unicode[STIX]{x2202}Z_{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D703}}\mathop{R}_{0}^{N^{2}-1}\hat{Z}_{1}=0.\end{eqnarray}$$
This constitutes a rather strong constraint on the behaviour of the function
$\hat{Z}_{1}(\unicode[STIX]{x1D703})$
at low aspect ratio and large field period number
$N$
since in that limit the factor
$R_{0}^{N^{2}-1}$
varies strongly in magnitude. This leads to the qualitative behaviour around the surface where QAS is satisfied (indeed, observed in our numerical solutions) of the perturbation favouring the ‘inboard’ side of the device (locations of small
$R$
), i.e. having greater amplitude there; the surface shape on the outboard side, in such cases, remains mostly unperturbed from its axisymmetric shape. At small
$N$
(i.e.
$N=2$
) and large aspect ratio (i.e. when the variation in the overall magnitude of
$R_{0}$
is small) this effect is less noticeable.
Appendix F. Non-existence of global QAS magnetic fields
We show here that the QAS condition can be satisfied globally only for the special cases of
$N=1$
, and, in that case, the deformation merely corresponds to trivial transformations that preserve axisymmetry. We demonstrate this at first order in the expansion. It is easy to integrate the QAS condition (it can be interpreted as a first-order ordinary differential equation since it is differential in only the
$R_{0}$
variable) to obtain a global solution:
$$\begin{eqnarray}P_{1}=\bar{P}_{1}(Z_{0})R_{0}^{1-N^{2}}.\end{eqnarray}$$
Substituting this solution into (3.26), we obtain
$$\begin{eqnarray}\bar{P}_{1}^{\prime \prime }/\bar{P}_{1}=N^{2}(1-N^{2})R_{0}^{-2}.\end{eqnarray}$$
The left-hand side of this equation does not depend on
$R$
whereas the right-hand side clearly does, unless it is zero, i.e.
$N=0,1$
(the
$N=0$
case is uninteresting). Taking
$\bar{P}_{1}^{\prime \prime }=0$
, we obtain
$$\begin{eqnarray}\bar{P}_{1}=A_{0}+A_{1}Z_{0}.\end{eqnarray}$$
The first term corresponds to a translation in the
$x$
–
$y$
plane, and the second term corresponds to a tilt, i.e. the solution corresponds to a axisymmetry preserving deformation. The solutions for
$\hat{R}_{1}$
and
$\hat{Z}_{1}$
, obtained from (3.22)–(3.23), are also consistent with this conclusion. The numerical solution also confirms this analytic solution. Therefore, the two solutions to the boundary value problem in fact satisfy the differential constraint globally, not just at the boundary, but correspond to trivial and uninteresting deformations.
