1 Introduction
Quasi-symmetry is a property of magnetic fields that ensures the confinement of collisionless particle orbits. Axisymmetric magnetic equilibria possess this property in a trivial sense, whereas the related class of stellarators, called quasi-axisymmetric (QAS), satisfies the symmetry in a way that is hidden to the naked eye (Nührenberg et al. Reference Nührenberg, Sindoni, Lotz, Troyon, Gori and Vaclavik1994).
The close relationship between axisymmetry and QAS suggests that the second class may be continuously connected to the first, and in particular that QAS stellarators may be obtained by deformation of axisymmetric equilibria (Boozer Reference Boozer2008). It has also been suggested that modifying tokamak equilibria by non-axisymmetric shaping might help overcome the stability issues that plague them, and a previous study, using conventional numerical optimization, has demonstrated that suitable QAS may indeed be found as deformed tokamak equilibria (Ku & Boozer Reference Ku and Boozer2009). The idea of passively stabilizing a tokamak by non-axisymmetric perturbations is also supported by a number of experimental results, when the perturbation generates a sufficient ‘external’ boost in rotational transform (W VIIA Team 1980; Pandya et al. Reference Pandya, ArchMiller, Cianciosa, Ennis, Hanson, Hartwell, Hebert, Herfindal, Knowlton and Ma2015). Solving the magnetohydrodynamic (MHD) equilibrium problem for optimal stellarator equilibria, without the use of numerical optimization algorithms (i.e. ‘direct construction’ of optimal solutions), is potentially beneficial due to the speedup offered (Landreman, Sengupta & Plunk Reference Landreman, Sengupta and Plunk2019). So far, the only ways to do this have involved approximations to the problem such as a small distance from the magnetic axis (Garren & Boozer Reference Garren and Boozer1991a,Reference Garren and Boozerb; Landreman & Sengupta Reference Landreman and Sengupta2018; Landreman et al. Reference Landreman, Sengupta and Plunk2019; Plunk, Landreman & Helander Reference Plunk, Landreman and Helander2019), or a small deviation from axisymmetry (Plunk & Helander Reference Plunk and Helander2018). However, solving these approximate problems can also lead to fundamental insights into the properties of the solutions, and the size of the solution space.
There are possible practical advantages of directly constructing QAS stellarator equilibria via the perturbation of axisymmetric equilibria, as compared to conventional optimization. For instance, the general perturbation can be constructed as a sum of independent QAS modes with different toroidal mode numbers. After pre-computation of these modes, the corresponding space of QAS equilibria can be easily scanned, without further computational cost, whereas each step of a conventional optimizer involves solving the equilibrium problem anew. Also, there is no fundamental constraint on axisymmetric equilibrium measures, such as the aspect ratio, so these may be set arbitrarily to explore new areas of stellarator design space, which may have been inaccessible with conventional optimization.
In a previous paper (Plunk & Helander Reference Plunk and Helander2018), it was proved that nearly axisymmetric magnetic fields can be constructed to satisfy the condition of quasi-axisymmetry on a single magnetic surface. These solutions, however, apply only to vacuum conditions, where the plasma itself does not contribute significantly to the magnetic field. The present work considers the more general case of finite-pressure equilibria. Formidable challenges are present in this general problem, starting with an increased complexity arising from the nonlinear coupling of multiple fields. The presence of singularities in the force balance equation makes the general problem of obtaining equilibria ill posed, even without the requirement of satisfying a special symmetry. As we will show, the issue of force balance singularities may be overcome, at least at first order in the expansion, by suitable choice of the zeroth-order rotational transform profile. The complexity of the system, however, makes it more difficult to establish the existence of solutions by the same methods employed by Plunk & Helander (Reference Plunk and Helander2018). We therefore turn to devising a method to numerically solve the system. This, as we find, gives evidence that the same problem as solved in the vacuum limit by Plunk & Helander (Reference Plunk and Helander2018), namely the problem of finding a perturbation of specified toroidal mode number 
$N$ that satisfies the condition of QAS on a single magnetic surface, is indeed well posed, at least in some practical sense.
The contents of the paper are as follows. In § 2 the basic equations and notation are established, and the ‘inverse’ MHD equilibrium problem formulation is described. In § 3, the expansion about axisymmetry is performed, and the equations are given to find perturbations satisfying QAS on a specified magnetic surface. The issue of force balance singularities is discussed, and a strategy to overcome them is described. In § 4, a numerical method is described to solve the first-order system, and a set of solutions are given, based on a zeroth-order ITER-like equilibrium. The VMEC (Hirshman & Whitson Reference Hirshman and Whitson1983) and BOOZ_XFORM (Sanchez et al. Reference Sanchez, Hirshman, Ware, Berry and Spong2000) codes are used to demonstrate that the solutions can satisfy the appropriate level of QAS as predicted by the theory.
2 Preliminaries
The MHD equilibrium equations are
\begin{gather} \boldsymbol{\nabla}\times{{\boldsymbol B}} = \mu_0\,{{\boldsymbol j}}, \end{gather}
\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}{{\boldsymbol B}} = 0. \end{gather}
\begin{gather}{{\boldsymbol j}}\times{{\boldsymbol B}} = \boldsymbol{\nabla} \psi \frac{\textrm{d} p}{\textrm{d}\psi}. \end{gather}
We assume topologically toroidal magnetic surfaces, here labelled by the flux function 
$\psi$. To solve these equations, we use a similar approach as previous works Garren & Boozer (Reference Garren and Boozer1991a,Reference Garren and Boozerb); Hegna (Reference Hegna2000); Boozer (Reference Boozer2002); Weitzner (Reference Weitzner2014). Boozer angles are denoted 
$\theta$ and 
$\varphi$. The contravariant form of 
${{\boldsymbol B}}$ is written
\begin{equation} {{\boldsymbol B}}_{\mathrm{con}} = \boldsymbol{\nabla}\psi \times \boldsymbol{\nabla}\theta - {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}(\psi) \boldsymbol{\nabla} \psi \times \boldsymbol{\nabla} \varphi,\end{equation}
where 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$ is the rotational transform, and 
$2{\rm \pi} \psi$ is the toroidal flux. This form of 
${{\boldsymbol B}}$ satisfies zero divergence, assuming a flux-surface geometry. The covariant form is written
\begin{equation} {{\boldsymbol B}}_{\mathrm{cov}} = G(\psi) \boldsymbol{\nabla}\varphi + I(\psi) \boldsymbol{\nabla}\theta + K(\psi, \theta, \varphi) \boldsymbol{\nabla} \psi.\end{equation}
This form is a consequence of 
${{\boldsymbol j}} \boldsymbol {\cdot } \boldsymbol {\nabla } \psi = 0$ (2.3), and Ampere's law (2.1); see e.g. Helander (Reference Helander2014).
The basic strategy to find an equilibrium is to assert 
${{\boldsymbol B}}_{\mathrm {con}} = {{\boldsymbol B}}_{\mathrm {cov}}$ together with force balance (2.3), relying on the fact that these forms of the magnetic field incorporate equations (2.1) and (2.2) as well as the assumption of topologically toroidal magnetic surfaces. Either the magnetic coordinates 
$\psi , \theta$ and 
$\varphi$ can be considered as the unknown functions of spatial coordinates (‘direct formulation’), or the coordinate mapping 
${\boldsymbol x} (\psi , \theta , \varphi )$ can be considered as an unknown function of the magnetic coordinates (‘inverse formulation’). Both formulations are used here.
It is convenient at zeroth order to solve the Grad–Shafranov equation (e.g. using the direct formulation). This means that we are able to specify the axisymmetric shape of the outer magnetic surface. We will also specify the current and pressure profiles at this stage, and consider them as fixed for the remainder of the calculation. We will use the indirect formulation for the problem at next order, i.e. the problem of QAS-preserving perturbations, as it casts the problem as a fixed boundary problem with QAS as the boundary condition.
2.1 Problem formulation
With the inverse formulation, the independent variables of the problem are the magnetic coordinates, and QAS is expressed as a simple constraint, 
$\partial B/\partial \varphi = 0$. Instead of using magnetic flux as a coordinate, we will use a coordinate system based on a dimensionless radial coordinate 
$\rho = \sqrt {\psi /\psi _b}$, where 
$\psi _b$ denotes the value of 
$\psi$ on the boundary surface. Note that most physical quantities are not analytic in 
$\psi$ at the magnetic axis (
$\psi = 0$), but can be expanded in 
$\rho$ (Garren & Boozer Reference Garren and Boozer1991a). This idea is motivated by considering 
$\rho$ and 
$\theta$ as polar coordinates and then assuming that a Taylor expansion can be made in the pseudo-Cartesian coordinates 
$\bar {x} = \rho \cos (\theta )$ and 
$\bar {y} = \rho \sin (\theta )$.
With the inverse formulation, the unknown of the theory is the coordinate mapping 
${\boldsymbol x}(\rho , \theta , \varphi )$, and the equilibrium equations are written in terms of various derivatives 
$\partial {\boldsymbol x}/\partial \rho$, and so forth. These equations can be translated into equations involving the metrics via the usual identities (reviewed in appendix A). Then equation 
${{\boldsymbol B}}_{\mathrm {con}} = {{\boldsymbol B}}_{\mathrm {cov}}$ becomes
\begin{equation} \rho\left(\frac{\partial {\boldsymbol x}}{\partial \varphi} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \frac{\partial {\boldsymbol x}}{\partial \theta}\right) = \bar{G} \frac{\partial {\boldsymbol x}}{\partial \rho}\times\frac{\partial {\boldsymbol x}}{\partial \theta} + \bar{I} \frac{\partial {\boldsymbol x}}{\partial \varphi}\times\frac{\partial {\boldsymbol x}}{\partial \rho} + \bar{K} \frac{\partial {\boldsymbol x}}{\partial \theta}\times\frac{\partial {\boldsymbol x}}{\partial \varphi}.\end{equation}
Force balance can be expressed as a scalar equation, since it only has a component in the 
$\boldsymbol {\nabla } \rho$ direction. We use 
${{\boldsymbol j}} = \mu _0^{-1}\boldsymbol {\nabla }\times {{\boldsymbol B}}_{\mathrm {cov}}$ and take the scaler product of (2.3) with 
$\boldsymbol {\nabla } \theta \times \boldsymbol {\nabla } \varphi$
\begin{equation} \left(\frac{\textrm{d} \bar{G}}{\textrm{d} \rho} - \frac{\partial \bar{K}}{\partial \varphi}\right) + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}\left(\frac{\textrm{d} \bar{I}}{\textrm{d} \rho} - \frac{\partial \bar{K}}{\partial \theta}\right) + \mu_0 \bar{J} \frac{\textrm{d}\bar{p}}{\textrm{d}\rho} = 0.\end{equation} We introduce the following normalized quantities: 
$\bar {G} = G/(2\psi _b), \bar {I} = I/(2\psi _b), \bar {K} = \rho K$ and 
$\bar {p} = p/(2\psi _b)^2$. Note that we define 
$\bar {J} = J/\rho$ in the limiting sense so that, although 
$J = (\boldsymbol {\nabla } \rho \boldsymbol {\cdot }\boldsymbol {\nabla }\theta \times \boldsymbol {\nabla }\varphi )^{-1}$ tends to zero with 
$\rho , \bar {J}$ does not. Finally, we will need the following expression for the regularized Jacobian:
\begin{equation} \bar{J} = \frac{\left|\dfrac{\partial {\boldsymbol x}}{\partial \varphi} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \dfrac{\partial {\boldsymbol x}}{\partial \theta}\right|^2}{\bar{G} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \bar{I}}.\end{equation}
Defining also 
$\bar {\boldsymbol {B}} = {\boldsymbol B}/(2\psi _b)$ we have the useful relation 
$\bar {J} = (\bar {G} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \bar {I})/\bar {B}^2$ ((B 1)) so that QAS can be expressed most conveniently here as
\begin{equation} \frac{\partial \bar{J}}{\partial \varphi} = 0. \end{equation}
In the present work, we look for solutions that satisfy this condition on a single magnetic surface; we will not consider here the question of whether this condition might, under special circumstances, be satisfied globally, i.e. uniformly in 
$\rho$.
3 The expansion about axisymmetry
We write the coordinate mapping 
${\boldsymbol x}(\rho , \theta , \varphi )$ as a series expansion in the small parameter 
$\epsilon$,
\begin{equation} {\boldsymbol x} = {\boldsymbol x}_0 + \epsilon {\boldsymbol x}_1 + \epsilon^2 {\boldsymbol x}_2 + \cdots, \end{equation}
where 
${\boldsymbol x}_0$ corresponds to the zeroth-order axisymmetric equilibrium. We will consider the pressure 
$\bar {p}$ and currents 
$\bar {G}$ and 
$\bar {I}$ as fixed to their zeroth-order values (there is no loss of generality as any higher-order variation in these functions can be absorbed into the zeroth-order forms). This confines our attention to axisymmetry breaking perturbations. We will, however, allow the deformation to modify 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$ and 
$\bar {K}$.
\begin{gather} \bar{K}(\rho, \varphi, \theta) = \bar{K}_0 + \epsilon \bar{K}_1 + \epsilon \bar{K}_2 + \cdots, \end{gather}
\begin{gather}{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} = {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_0 + \epsilon {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_1 + \epsilon^2 {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_2 + \cdots. \end{gather} For a nearly axisymmetric equilibrium, it is sensible to take the components of (3.11) along the cylindrical unit vectors 
$\hat {\boldsymbol R}, \hat {\boldsymbol {\phi }}, \hat {\boldsymbol z}$ (such that 
$\hat {\boldsymbol R} \times \hat {\boldsymbol {\phi }} = \hat {\boldsymbol z}$).
3.1 Order 
$\epsilon ^0$
The zeroth-order coordinate mapping is (see also appendix D)
\begin{equation} {\boldsymbol x}_0 = \hat{\boldsymbol R} R_0(\rho, \theta) + \hat{\boldsymbol z} Z_0(\rho, \theta), \end{equation}
where the cylindrical unit vectors are functions of the geometric toroidal coordinate, related to the boozer angle by 
$\varphi = \phi + \nu$, and are expanded as
\begin{equation} \phi = \varphi - \nu_0 - \nu _1 - \cdots,\end{equation}
so that 
$\phi _0 = \varphi - \nu _0$ and 
$\phi _1 = -\nu _1$, etc., and, for simplicity, the unit vectors will be defined according to the zeroth-order expression of the geometric toroidal angle,
\begin{gather} \hat{\boldsymbol R} = \hat{\boldsymbol R}(\phi_0) = \hat{\boldsymbol R}(\varphi - \nu_0), \end{gather}
\begin{gather}\hat{\boldsymbol{\phi}} = \hat{\boldsymbol{\phi}}(\phi_0) = \hat{\boldsymbol{\phi}}(\varphi - \nu_0). \end{gather}With these definitions, derivatives of the zeroth-order coordinate mapping are evaluated as
\begin{gather} \frac{\partial {\boldsymbol x}_0}{\partial \rho} = \hat{\boldsymbol R} \frac{\partial R_0}{\partial \rho} + \hat{\boldsymbol z} \frac{\partial Z_0}{\partial \rho} - \hat{\boldsymbol{\phi}} R_0 \frac{\partial \nu_0}{\partial \rho}, \end{gather}
\begin{gather}\frac{\partial {\boldsymbol x}_0}{\partial \theta} = \hat{\boldsymbol R} \frac{\partial R_0}{\partial \theta} + \hat{\boldsymbol z} \frac{\partial Z_0}{\partial \theta} - \hat{\boldsymbol{\phi}} R_0 \frac{\partial \nu_0}{\partial \theta}, \end{gather}
\begin{gather}\frac{\partial {\boldsymbol x}_0}{\partial \varphi} = \hat{\boldsymbol{\phi}} R_0. \end{gather}The zeroth-order MHD constraint is
\begin{equation} \rho\left(\frac{\partial {\boldsymbol x}_0}{\partial \varphi} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_0 \frac{\partial {\boldsymbol x}_0}{\partial \theta}\right) = \bar{G} \frac{\partial {\boldsymbol x}_0}{\partial \rho}\times\frac{\partial {\boldsymbol x}_0}{\partial \theta} + \bar{I} \frac{\partial {\boldsymbol x}_0}{\partial \varphi}\times\frac{\partial {\boldsymbol x}_0}{\partial \rho} + \bar{K}_0 \frac{\partial {\boldsymbol x}_0}{\partial \theta}\times\frac{\partial {\boldsymbol x}_0}{\partial \varphi},\end{equation}
where (3.8)–(3.10) can be substituted and the equation projected along the unit vectors 
$\hat {\boldsymbol R}, \hat {\boldsymbol {\phi }}$ and 
$\hat {\boldsymbol z}$ to obtain three coupled equations. Note that we avoid explicitly writing the lengthy equations that result, and will do likewise with others that follow, especially when they do not give any useful insight. Force balance is
\begin{equation} \frac{\textrm{d} \bar{G}}{\textrm{d} \rho} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_0 \left(\frac{\textrm{d} \bar{I}}{\textrm{d} \rho} - \frac{\partial \bar{K}_0}{\partial \theta}\right) + \mu_0 \bar{J}_0 \frac{\textrm{d}\bar{p}}{\textrm{d}\rho} = 0.\end{equation}3.1.1 Inverting the Grad–Shafranov solution
It is convenient to use Grad–Shafranov (GS) theory to obtain the zeroth-order equilibrium. This approach gives control of the axisymmetric plasma shape, and also benefits from the existing understanding of the equation and its numerical solution. A solution of the GS equation is the poloidal flux function 
$\varPsi (R, z)$ is obtained from a given pressure function 
$p$, and the poloidal flux function 
$G$. From these quantities, the corresponding profiles 
$I$ and 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0$, the current potential 
$K_0$ and coordinate mapping components 
$R_0, Z_0$ and 
$\nu _0$ can be calculated. To perform the coordinate inversion, derivatives of the GS solution are computed, so a high degree of accuracy is needed. A method is described in appendix D.
3.2 Order $\epsilon ^1$
As in Plunk & Helander (Reference Plunk and Helander2018), we do not modify the toroidal angle beyond zeroth order in the expansion (
$\nu _1 = 0$, etc., in (3.5)), but instead consider the corrections to the coordinate mapping to have a component in the 
$\hat {\boldsymbol {\phi }}$ direction, i.e.
\begin{equation} {\boldsymbol x}_1 = \hat{\boldsymbol R} R_1(\theta, \psi, \varphi) + \hat{\boldsymbol z} Z_1(\theta, \psi, \varphi) + \hat{\boldsymbol{\phi}} \varPhi_1(\theta, \psi, \varphi),\end{equation}from which it follows that
\begin{gather} \frac{\partial {\boldsymbol x}_1}{\partial \rho} = \hat{\boldsymbol R} \left(\frac{\partial R_1}{\partial \rho} + \varPhi_1 \frac{\partial \nu_0}{\partial \rho}\right) + \hat{\boldsymbol z} \frac{\partial Z_1}{\partial \rho} + \hat{\boldsymbol{\phi}} \left(\frac{\partial \varPhi_1}{\partial \rho} - R_1 \frac{\partial \nu_0}{\partial \rho}\right) \end{gather}
\begin{gather}\frac{\partial {\boldsymbol x}_1}{\partial \theta} = \hat{\boldsymbol R} \left(\frac{\partial R_1}{\partial \theta} + \varPhi_1 \frac{\partial \nu_0}{\partial \theta}\right) + \hat{\boldsymbol z} \frac{\partial Z_1}{\partial \theta} + \hat{\boldsymbol{\phi}} \left(\frac{\partial \varPhi_1}{\partial \theta} - R_1 \frac{\partial \nu_0}{\partial \theta}\right) \end{gather}
\begin{gather}\frac{\partial {\boldsymbol x}_1}{\partial \varphi} = \hat{\boldsymbol R}\left( \frac{\partial R_1}{\partial \varphi} - \varPhi_1 \right) + \hat{\boldsymbol z} \frac{\partial Z_1}{\partial \varphi} + \hat{\boldsymbol{\phi}} \left(R_1 + \frac{\partial \varPhi_1}{\partial \varphi} \right). \end{gather}
As 
$\varphi$ is an ignorable coordinate in the properly formulated first-order equilibrium equations, we will assume
\begin{gather} R_1 = \hat{R}_1(\theta, \psi) \exp(\textrm{i} N \varphi) + \mathrm{c.c.}, \end{gather}
\begin{gather}Z_1 = \hat{Z}_1(\theta, \psi) \exp(\textrm{i} N \varphi) + \mathrm{c.c.}, \end{gather}
\begin{gather}\varPhi_1 = \hat{\varPhi}_1(\theta, \psi) \exp(\textrm{i} N \varphi) + \mathrm{c.c.}, \end{gather}
\begin{gather}\bar{K}_1 = \hat{\bar{K}}_1(\theta, \psi) \exp(\textrm{i} N \varphi) + \mathrm{c.c.}, \end{gather}
with 
$N \neq 0$ an integer. The deformation is thus non-axisymmetric, and the axisymmetric (
$\varphi$-averaged) part of the local MHD constraint (C 4) is 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _1( \bar {G} g_{22}^{(0)} - \bar {I} g_{23}^{(0)} ) = 0$, from which we conclude that
\begin{equation} {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_1 = 0. \end{equation}The first-order MHD constraint is then
\begin{align} \rho\left(\frac{\partial {\boldsymbol x}_1}{\partial \varphi} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_0 \frac{\partial {\boldsymbol x}_1}{\partial \theta}\right) &= \bar{G} \left(\frac{\partial {\boldsymbol x}_0}{\partial \rho}\times\frac{\partial {\boldsymbol x}_1}{\partial \theta} + \frac{\partial {\boldsymbol x}_1}{\partial \rho}\times\frac{\partial {\boldsymbol x}_0}{\partial \theta}\right) + \bar{I} \left(\frac{\partial {\boldsymbol x}_0}{\partial \varphi}\times\frac{\partial {\boldsymbol x}_1}{\partial \rho} + \frac{\partial {\boldsymbol x}_1}{\partial \varphi}\times\frac{\partial {\boldsymbol x}_0}{\partial \rho}\right) \nonumber\\ &\quad+ \bar{K}_0 \left(\frac{\partial {\boldsymbol x}_0}{\partial \theta}\times\frac{\partial {\boldsymbol x}_1}{\partial \varphi} + \frac{\partial {\boldsymbol x}_1}{\partial \theta}\times\frac{\partial {\boldsymbol x}_0}{\partial \varphi}\right) + \bar{K}_1 \frac{\partial {\boldsymbol x}_0}{\partial \theta}\times\frac{\partial {\boldsymbol x}_0}{\partial \varphi}. \end{align}
What is needed is the 
$\exp (\textrm {i} N \varphi )$ component of this equation, obtained by substituting equations (3.17)–(3.20) into 
${\boldsymbol x}_1$, (3.13), and its derivatives, (3.14)–(3.16). The further substitution of zeroth-order expressions, (3.8)–(3.10), and projection along the cylindrical unit vectors, then yield a set of three equations for the unknowns 
$\hat {R}_1, \hat {Z}_1, \hat {\varPhi }_1, \hat {\bar {K}}_1$ in terms of the known zeroth-order solutions 
$R_0, Z_0, \nu _0$ and 
$\bar {K}_0$. The system is completed with the force balance equation,
\begin{equation} - \textrm{i} N \hat{\bar{K}}_1 - {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_0\frac{\partial \hat{\bar{K}}_1}{\partial \theta} + \hat{\bar{J}}_1 \mu_0 \frac{\textrm{d}\bar{p}}{\textrm{d}\rho} = 0.\end{equation}
The 
$\exp (\textrm {i} N \varphi )$ component of the first-order Jacobian, 
$\hat {\bar {J}}_1$, is obtained from (2.8) by substituting equations (3.9) and (3.10) and (3.15) and (3.16), into the following expression:
\begin{equation} \bar{J}_1 = 2\frac{\left(\dfrac{\partial {\boldsymbol x}_0}{\partial \varphi} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_0 \dfrac{\partial {\boldsymbol x}_0}{\partial \theta}\right)\cdot\left(\dfrac{\partial {\boldsymbol x}_1}{\partial \varphi} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_0 \dfrac{\partial {\boldsymbol x}_1}{\partial \theta}\right)}{\bar{G} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_0 \bar{I}}. \end{equation} We note that QAS implies that 
$\hat {\bar {J}}_1 = 0$, so force balance on any QAS surface reduces to 
$\textrm {i} N \hat {\bar {K}}_1 + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0 \partial \hat {\bar {K}}_1/\partial \theta = 0$, which, assuming irrational 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0$, implies
\begin{equation} \hat{\bar{K}}_1 = 0.\end{equation}
On magnetic surfaces where QAS is not satisfied, the possibility of resonances in (3.23) must be considered. It is easy to see that the equation can be uniquely solved for 
$\hat {\bar {K}}_1$, periodic in 
$\theta$, if 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0$ is not equal to a rational number. Actually, some rational numbers are resonant, and some are not, in particular there are resonances at any magnetic surface where 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0$ satisfies
\begin{equation} {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_0 = \frac{N}{m}, \end{equation}
for arbitrary integer 
$m$. One strategy to avoid resonances is to constrain 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0$ to lie between two neighbouring singular values. In that case, force balance can be considered ‘soluble’ throughout the plasma volume. Note that, assuming 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0 \sim 1$, such ‘safe’ ranges become increasingly narrow at large 
$N$, although resonances may be considered ‘high order’ in this limit, and therefore less likely to pollute the solution.
Note that (3.25) demonstrates that force balance is non-resonant on a QAS magnetic surface. Actually, this reflects a general non-perturbative property, which follows directly from exact force balance and the relationship between 
$B$ and 
$\bar {J}$ ((2.7) and (B 1)), namely that quasi-symmetry drastically simplifies the source term in force balance (Fourier series of 
$\bar {J}$ in 
$\theta$ and 
$\varphi$ has only terms of a single helicity), so that an MHD equilibrium that is quasi-symmetric globally (on all magnetic surfaces) should be free from non-trivial resonances; see also Burby, Kallinikos & MacKay (Reference Burby, Kallinikos and MacKay2019) and Rodríguez, Helander & Bhattacharjee (Reference Rodríguez, Helander and Bhattacharjee2020).
To summarize, at first order the equations to be solved are the three components of (3.22), coupled with force balance, (3.23), for the four unknown functions 
$\hat {R}_1(\rho , \theta ), \hat {Z}_1(\rho , \theta ), \hat {\varPhi }_1(\rho , \theta ), \hat {\bar {K}}_1(\rho , \theta )$. The domain is the unit disk, 
$\rho \in [0, 1]$ and the boundary condition is QAS, which translates to 
$\hat {\bar {K}}(\rho , \theta )|_{\rho = 1} = 0$. No rotational transform is obtained at this order, but the first-order solution does generally induce transform at 
${ {\textit{O}}}(\epsilon ^2)$, i.e. enters the computation of 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _2$.
3.3 Order $\epsilon ^2$
Proceeding to the next order, the equations are quite similar to before, but now include terms that are quadratic in first-order quantities. The second-order coordinate mapping is
\begin{equation} {\boldsymbol x}_2 = \hat{\boldsymbol R} R_2(\theta, \psi, \varphi) + \hat{\boldsymbol z} Z_2(\theta, \psi, \varphi) + \hat{\boldsymbol{\phi}} \varPhi_2(\theta, \psi, \varphi),\end{equation}and its derivatives are
\begin{gather} \frac{\partial {\boldsymbol x}_2}{\partial \rho} = \hat{\boldsymbol R} \left(\frac{\partial R_2}{\partial \rho} + \varPhi_2 \frac{\partial \nu_0}{\partial \rho}\right) + \hat{\boldsymbol z} \frac{\partial Z_2}{\partial \rho} + \hat{\boldsymbol{\phi}} \left(\frac{\partial \varPhi_2}{\partial \rho} - R_2 \frac{\partial \nu_0}{\partial \rho}\right) \end{gather}
\begin{gather}\frac{\partial {\boldsymbol x}_2}{\partial \theta} = \hat{\boldsymbol R} \left(\frac{\partial R_2}{\partial \theta} + \varPhi_2 \frac{\partial \nu_0}{\partial \theta}\right) + \hat{\boldsymbol z} \frac{\partial Z_2}{\partial \theta} + \hat{\boldsymbol{\phi}} \left(\frac{\partial \varPhi_2}{\partial \theta} - R_2 \frac{\partial \nu_0}{\partial \theta}\right) \end{gather}
\begin{gather}\frac{\partial {\boldsymbol x}_2}{\partial \varphi} = \hat{\boldsymbol R}\left( \frac{\partial R_2}{\partial \varphi} - \varPhi_2 \right) + \hat{\boldsymbol z} \frac{\partial Z_2}{\partial \varphi} + \hat{\boldsymbol{\phi}} \left(R_2 + \frac{\partial \varPhi_2}{\partial \varphi} \right). \end{gather}
The appearance of the nonlinear terms occurs (in the MHD constraint and force balance equation) at toroidal mode numbers 
$\pm 2N$, and also in the axisymmetric component, which now must be solved to obtain 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _2$. We note that the appearance of toroidal mode numbers 
$\pm 2N$ at second order implies a denser set of possible force balance resonances at higher orders in the expansion, i.e. 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0 = 2N/m$, which may justify further restriction on the chosen profile for 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0$. Even if the problem will only be solved at first order, higher-order resonances may occur in the exact force balance equation that must be satisfied by the full equilibrium.
In the vacuum case (Plunk & Helander Reference Plunk and Helander2018), 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _2$ was obtained as a solubility constraint of the axisymmetric component (
$\varphi$-average) of the local MHD constraint, (C 4). This result does not appear to generalize in a simple way, implying that the full system (MHD constraint plus force balance) must be solved to obtain 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _2$.\vspace*{-3pt}
4 Numerical solution
The task now is to solve the system composed of (3.22) and (3.23) for the unknowns 
$\hat {R}_1, \hat {Z}_1, \hat {\varPhi }_1, \hat {\bar {K}}_1$, subject to QAS (
$\hat {\bar {K}} = 0$) on a specified magnetic surface, typically the outermost magnetic surface (
$\rho = 1$). Note that this boundary surface need not necessarily be taken to be located at the plasma edge. It has been suggested (Henneberg et al. Reference Henneberg, Drevlak, Nührenberg, Beidler, Turkin, Loizu and Helander2019) that it may be optimal to satisfy QAS at some intermediate magnetic surface, which can be implemented here by redefinition of the coordinate 
$\rho$, or simply choosing a boundary 
$\rho < 1$. Henceforth, we assume 
$\rho = 1$ for simplicity.
The ‘pseudo-Cartesian’ coordinates 
$\bar {x} = \rho \cos (\theta )$ and 
$\bar {y} = \rho \sin (\theta )$ are used for numerical purposes, instead of the polar coordinates 
$\rho$ and 
$\theta$. These have the advantage that they do not possess the singularity of the polar coordinates (
$\rho , \theta$) as 
$\rho \rightarrow 0$, and they do not require periodicity to be enforced in 
$\theta$, or any analyticity at 
$\rho = 0$. The only advantage found in using 
$\rho$-
$\theta$ coordinates is to explicitly observe the development of non-analyticity in the solutions on resonant surfaces (satisfying 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0 = N/m$).
The finite element method is used, reformulating the problem as an equivalent ‘least squares’ problem. The least squares finite element method offers better convergence and stability properties for systems of first-order partial differential equations (Jiang Reference Jiang1998). To show how the problem is reformulated, we introduce the vector field 
${\boldsymbol u}(\bar {x}, \bar {y}) = [\hat {\bar {K}}_1, \hat {R}_1, \hat {Z}_1, \hat {\varPhi }_1]^\textrm {T}$, together with the inner product
\begin{equation} \langle {\boldsymbol v} | {\boldsymbol u} \rangle = \int \textrm{d}\bar{x} \,\textrm{d}\bar{y} {\boldsymbol v}^*\boldsymbol{\cdot}{\boldsymbol u}, \end{equation}
where 
${\boldsymbol v}^*$ denotes the complex conjugate of 
${\boldsymbol v}$, and the integral is performed over the computational domain, the unit disk 
$\varOmega$. The original first-order system of equations (i.e. the 
$\hat {\boldsymbol R}, \hat {\boldsymbol {\phi }}$ and 
$\hat {\boldsymbol z}$ components of (3.22), coupled with (3.23)) can be written as 
${\mathcal {L}} {\boldsymbol u} = 0$, where
\begin{equation} [{\mathcal{L}} {\boldsymbol u}]_i = a_{ij} u_j + \alpha_{ijk} \partial_j u_k, \end{equation}
where 
$u_j$ denotes the 
$j$th component of 
${\boldsymbol u}, \partial _1$ and 
$\partial _2$ denote 
$\partial /\partial \bar {x}$ and 
$\partial /\partial \bar {y}$, respectively, and the tensors 
$a_{ij}$ and 
$\alpha _{ijk}$ encode the coefficients of the system of equations. The adjoint of 
${\mathcal {L}}$ is denoted as 
${\mathcal {L}}^\dagger$, and is given by
\begin{equation} [{\mathcal{L}}^\dagger {\boldsymbol u}]_i = a_{ji}^* u_j - \partial_j(\alpha_{kji}^* u_k). \end{equation}With these definitions, our problem is transformed into solving the following eigenvalue problem
\begin{equation} {\mathcal{L}}^\dagger {\mathcal{L}} {\boldsymbol u} = \lambda {\boldsymbol u}, \end{equation}
subject to the QAS boundary condition 
$u_1 = 0$ on the boundary 
$\partial \varOmega$, for solutions with eigenvalue 
$\lambda \rightarrow 0$. This system is generated by computer algebra, and not explicitly written down, due to its complexity.
4.1 Examples
To demonstrate that the above numerical method works, in practice, and give a flavour of possible solutions, we consider perturbations of model tokamak equilibria, based on ITER. The ITER-like equilibria have their outer surface shape defined by the Solev'ev equilibrium given in Pataki et al. (Reference Pataki, Cerfon, Freidberg, Greengard and O'Neil2013), but scaled up so that the magnetic axis has a radial position of 
$6.68$ m and a total toroidal flux over 
$2{\rm \pi}$ of 
$15.7$ Weber. The model pressure profiles that are linear in the poloidal flux function 
$\varPsi$, as shown in figure 1, and three different constant rotational transform profiles are considered, 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0 = 0.202, 0.47$ and 
$0.98$. These values are chosen to avoid resonances for 
$N = 2, 4$ and 
$8$; see § 3.2.
Figure 1. ITER-like equilibria with constant rotational transform profiles. The pressure profile (a) is plotted versus normalized poloidal flux 
$s_p = \varPsi /|\varPsi _{{{\rm axis}}}|$ where 
$\varPsi _{{{{{\rm axis}}}}}$ is the value on axis, and 
$\varPsi$ is taken to be zero at the outermost surface. Note that the mesh on the magnetic surface is made of lines of constant geometric angle 
$\phi$ and constant Boozer angle 
$\theta$; the end cap shows lines of constant 
$\theta$ and 
$\varPsi$. (a) Pressure profile and (b) outer magnetic surface shape.
To independently evaluate the first-order QAS numerical solutions, the outer surface shape can be generated and provided to the VMEC code (Hirshman & Whitson Reference Hirshman and Whitson1983) as input for a fully nonlinear calculation, as was done in Plunk & Helander (Reference Plunk and Helander2018), and the result then passed to the BOOZ_XFORM code (Sanchez et al. Reference Sanchez, Hirshman, Ware, Berry and Spong2000) to check the level of QAS as predicted by the theory. To produce the surface shape, the perturbation amplitude is controlled via the arbitrary small parameter 
$\epsilon$ in 
${\boldsymbol x} \approx {\boldsymbol x}_0 + \epsilon {\boldsymbol x}_1$. Three such surfaces are shown in figure 2.
Figure 2. Outer magnetic surface shapes for ITER-like QAS equilibria with near-unity rotational transform, 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0 = 0.98$. Two viewing angles are shown, from the top and from the side, with the side view showing half a toroidal turn; a sample of field line segments is plotted in red. The mesh on the magnetic surface correspond to lines of constant Boozer angle, 
$\theta$ and 
$\varphi$. (a) 
$N = 2$, (b) 
$N = 4$ and (c) 
$N = 8$.
The solutions are valid to first order in the expansion, and it can therefore be expected that the error in QAS should scale as 
$\epsilon ^2$, the confirmation of which is shown for one case in figure 3, where the error is measured as follows:
\begin{equation} Q = \frac{\left(\displaystyle\sum_{m, n\neq 0}|\hat{B}_{mn}|^2\right)^{1/2}}{\left(\displaystyle\sum_{m, n}|\hat{B}_{mn}|^2\right)^{1/2}},\end{equation}
where 
$\hat {B}_{mn}$ is the Fourier coefficient of 
$|{\boldsymbol B}|$ calculated in Boozer angles by the BOOZ_XFORM. It should be noted that not all of the cases reported here match so closely with the theoretical scaling. Some have only a limited range at larger values of 
$\epsilon$ where the quadratic scaling is observed, and exhibit a weaker 
$\epsilon ^1$ scaling for smaller values of 
$\epsilon$, associated with non-QAS perturbations. Although it is expected, for instance, that numerical error in the first-order solution can introduce 
$\epsilon ^1$ scaling which must dominate at sufficiently small 
$\epsilon$, it does not appear that the linear scaling observed here is related to numerical error in the three codes being used, of the type introduced by finite resolution. It is therefore suspected that a more fundamental issue is at fault, for instance (i) the presence of force balance singularities in the fully nonlinear calculation performed by VMEC (which are formally absent from our first-order calculation of 
${\boldsymbol x}_1$), or (ii) the possibility that the problem we are solving (QAS on a single surface of non-axisymmetric perturbations) is sometimes (or always) not well posed; this issue will be investigated in future work. Nevertheless, the low observed QAS error in the solutions obtained so far indicate that the numerical method developed here should be practically useful.
Figure 3. Example of perturbed tokamak equilibrium (
$N = 2$) with ITER-like shaping. Unperturbed rotational transform is 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0 = 0.202$ at all radial locations. (a) theoretical scaling of 
$\epsilon ^2$ is well satisfied for QAS error, defined in (4.5). (b) outer surface shape visualized for case of strongest shaping (largest perturbation), with 
$1$ toroidal field period plotted, and a sample of field line segments.
One issue encountered with using the inverse representation of a magnetic equilibrium is that the coordinate mapping is not generally invertible for the magnetic coordinates. Invertibility breaks down when distinct points in the magnetic coordinate space, say (
$\rho _1, \theta _1, \varphi _1$) and (
$\rho _2, \theta _2, \varphi _2$), yield the same point in physical space, e.g. 
${\boldsymbol x}(\rho _1, \theta _1, \varphi _1) = {\boldsymbol x}(\rho _2, \theta _2, \varphi _2)$. The QAS solutions here, being based on known axisymmetric equilibria, will not suffer from this problem if the perturbation amplitude is chosen to be sufficiently small. However, the problem can be reliably encountered at large values of 
$\epsilon$, as demonstrated in figure 4. What is remarkable about this phenomenon, which is associated with the perturbation overwhelming the zeroth-order mapping, is that the theoretical QAS scaling tends to hold even as the singular point is approached, as demonstrated by figure 3. Therefore, the VMEC solutions shown here are generally chosen to correspond to a value of 
$\epsilon$ close to the singular point, but not so large as to create sharp features in the outer magnetic surface that require more than 
$10-20$ Fourier harmonics to properly resolve in VMEC.
Figure 4. As the perturbation amplitude 
$\epsilon$ is increased, the magnetic surfaces ‘reconnect’, invalidating the solution. Here is an example, with 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0 = 0.98$ and 
$N=4$, showing a single field period. For comparison, the corresponding case in table 1, and plotted in figure 2(b) corresponds to 
$\epsilon = 0.75$. (a) 
$\epsilon = 0.5$, (b) 
$\epsilon = 1.8$ and (c) 
$\epsilon = 2.5$.
Using the procedure described above, the QAS solutions, although formally only perturbative, can yield strongly shaped plasma equilibria with reasonable level of QAS, and finite ‘external’ rotational transform, as measured by the difference between the total rotational transform and that of the original axisymmetric equilibrium. This is shown by table 1, where a total of nine cases are described, corresponding to three toroidal mode numbers applied to the equilibria of three different values of constant rotational transform. Each row of the table corresponds to a single value of 
$\epsilon$ (although a sequence of values was generally calculated to investigate scaling). The fraction of external rotational transform generated by the perturbation is given in the third column
\begin{equation} \left. f^\mathrm{ext}_{{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}} = \frac{{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} - {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_0}{{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}} \right\rvert_{\rho = 1}. \end{equation}
Next is the root-mean-squared value of the modulus of the perturbed magnetic field, divided by zeroth-order field strength, with the average performed over 
$\theta$ and 
$\varphi$, denoted 
$E$
\begin{equation} E = \left\langle\frac{|\delta {\boldsymbol B}|}{|{\boldsymbol B}_0|}\right\rangle_\mathrm{rms} = \left( \frac{\epsilon^2}{4{\rm \pi} ^2} \int \, \textrm{d} \varphi \, \textrm{d} \theta \frac{\left| \dfrac{\partial {\boldsymbol x}_1}{\partial \varphi} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_0 \dfrac{\partial {\boldsymbol x}_1}{\partial\theta} \right|^2}{\left| \dfrac{\partial {\boldsymbol x}_0}{\partial \varphi} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}_0 \dfrac{\partial {\boldsymbol x}_0}{\partial\theta} \right|^2} \right)^{1/2}, \end{equation}
where we note that the above expression assumes the first-order Jacobian to be zero (QAS). This measure gives some sense of how strong the perturbation is, and may be used to estimate the size of external current distributions needed to achieve the total field. The next column provides the QAS error, 
$Q$, defined in (4.5). The chosen values of 
$\epsilon$ are arbitrary, so it is useful to calculate normalized values to compare the various solutions. For that reason we give inferred values of the magnetic perturbation measures 
$E_{10}$ and 
$E_{15}$, that would be obtained for external rotational transforms of 
$10\,\%$ and 
$15\,\%$, respectively. Analogous quantities for QAS error are denoted 
$Q_{10}$ and 
$Q_{15}$. These are calculated for each case by using the fact that 
$\delta {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} - {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0$ scales theoretically as 
$\epsilon ^2$ (confirmed for all cases), as does 
$Q$, whereas 
$E$ scales as 
$\epsilon ^1$. We stress that these values, being obtained from first-order solutions, and not benefitting from any further optimization, should only be taken as an indicator of what can be achieved by perturbing an axisymmetric equilibrium. However, what seems clear is that, despite the fact that the perturbation of the field 
$E$ scales as 
$\epsilon ^1$ whereas the external transform scales as 
$\epsilon ^2$, significant values of the latter can still be achieved at modest values of the former, as for instance shown by the 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0 = 0.43, N=4$ case where the root-mean-square field strength fraction is not much larger that the external rotational transform fraction.
Table 1. Summary of results for ITER-like QAS equilibria.
An interesting qualitative feature of the QAS perturbations is the tendency to ‘localize’ to the inboard side (e.g. at lower values of the radial coordinate 
$R$), in the sense that the amplitude of the perturbed magnetic field is larger there than on the outboard. This is quantified in the final column of the table 1, labelled ‘Loc’, where we calculate the ratio of the root-mean-squared value of the perturbed magnetic field 
$\delta {\boldsymbol B}$ on the outboard (defined such that 
$\theta = 0$) to inboard (
$\theta = {\rm \pi}$), where the average is done only over the toroidal angle 
$\varphi$. This feature is more pronounced at larger 
$N$ and lower aspect ratio, as was also observed for the vacuum case (see the appendix of Plunk & Helander Reference Plunk and Helander2018; the explanation here may be similar). We note that the high-
$N$ perturbations also only weakly penetrate radially into the plasma, as the rotational transform 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$ can be observed to fall rapidly, from the edge, to the unperturbed value 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0$.
5 Conclusion
This work gives the first set of results showing the direct construction of QAS perturbations of non-trivial axisymmetric plasma equilibria. It has been demonstrated that, despite the perturbative nature of the calculation, relatively strongly shaped stellarator equilibria may be obtained with significant external rotational transform (10 %–15 %), at a similar level of average perturbed magnetic field. This finding agrees with a previous study (Ku & Boozer Reference Ku and Boozer2009) using conventional numerical optimization. However, the method of the present work allows for a more extensive exploration of the design space, as the general QAS perturbation corresponds to a sum of modes with suitably non-resonant toroidal mode numbers. This potentially opens new avenues for exploring the concept of a stellarator–tokamak hybrid device.
One interesting initial finding is that relatively high-
$N$ perturbations (here as high as 
$N=8$) seem to efficiently produce a finite external rotational transform (e.g. 
$10\,\%$), while diminishing strongly in amplitude both radially and poloidally, and showing very good satisfaction of QAS, much less than 
$1\,\%$ error. Such perturbations might be generated by a more ‘modest’ distribution of coils localized to the inboard side of the plasma. Additionally, with the perturbation localized to the high-field side of the device, it should only significantly affect the radial drift of barely trapped particles, rendering the overall neoclassical transport especially small.
One benefit of perturbative studies like the present is the ability to characterize the size of the solution space of optimal stellarators. Similar to what was found by Plunk & Helander (Reference Plunk and Helander2018), we conclude that the freedom in QAS designs comes from (i) the zeroth-order equilibrium, which is in this case includes plasma profiles in addition to the two-dimensional shaping (e.g. triangularity, elongation, aspect ratio, etc.), and (ii) the solution space of the QAS perturbation. For the latter, there are also some differences: first, it appears that, for fixed toroidal mode number 
$N$, the solution is unique, whereas Plunk & Helander (Reference Plunk and Helander2018) found that solutions came in pairs – the latter situation may stem from the symmetry of the 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0 = 0$ scenario; we note that this limit cannot be approached within the current framework, as the resonant values of 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$ become dense near 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} = 0$. Second, the choice of toroidal mode number 
$N$ is constrained, at least formally, by the profile 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0(\rho )$, in the sense that resonances (
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0 = N/m$) must be avoided to guarantee smooth solutions at first order, with further resonances might be considered at higher order. Therefore, the realizable solution space for QAS perturbations of a given axisymmetric equilibrium may be limited to a small number of toroidal mode numbers. Such a small space might be rapidly explored to identify QAS equilibria that satisfy additional requirements.
The success demonstrated here in directly constructing QAS solutions with an inverse method, using Boozer coordinates, gives some hope that the fully nonlinear problem may be formulated and solved in a similar fashion, i.e. with a code similar to VMEC that would obtain quasi-symmetric stellarator equilibria directly, and without approximations. To accomplish this, it is necessary to identify the appropriate amount of boundary information to yield a well-posed problem; the findings of this paper should provide a useful guide in this endeavour.
Acknowledgements
Thanks to P. Helander and S. Henneberg for useful discussions, and to S. Lazerson, E. Strumberger and J. Geiger for help with VMEC. Thanks also to Y. Turkin and P. Aleynikov for help with tokamak equilibrium files.
Editor Peter Catto thanks the referees for their advice in evaluating this article.
Declaration of interests
The author reports no conflict of interest.
Appendix A. Magnetic geometry
From the coordinates (
$\rho , \theta , \varphi$) we define local basis vectors (
$\boldsymbol {\nabla }\rho , \boldsymbol {\nabla }\theta , \boldsymbol {\nabla }\varphi$) and (
$\boldsymbol e_1 = \partial {\boldsymbol x}/\partial \rho , \boldsymbol e_2 = \partial {\boldsymbol x}/\partial \theta , \boldsymbol e_3 = \partial {\boldsymbol x}/\partial \varphi$). The metric components are defined in the usual way
\begin{equation} g_{ij} \equiv \boldsymbol e_i\boldsymbol{\cdot}\boldsymbol e_j, \end{equation}and the Jacobian for these coordinates is
\begin{equation} J = \frac{1}{\boldsymbol{\nabla}\rho\boldsymbol{\cdot}(\boldsymbol{\nabla}\theta\times\boldsymbol{\nabla}\varphi)} = \boldsymbol e_1\boldsymbol{\cdot} (\boldsymbol e_2\times\boldsymbol e_3). \end{equation} Additionally, assigning 
$(u_1, u_2, u_3) = (\rho , \theta , \varphi )$, we see that 
$\boldsymbol e_i\boldsymbol {\cdot }\boldsymbol {\nabla } u_j = \delta _{ij}$, and the following identities are easily verified
\begin{gather} \boldsymbol e_1 = J(\boldsymbol{\nabla} u_2\times\boldsymbol{\nabla} u_3),\quad \boldsymbol e_2 = J(\boldsymbol{\nabla} u_3\times\boldsymbol{\nabla} u_1),\quad \boldsymbol e_3 = J(\boldsymbol{\nabla} u_1\times\boldsymbol{\nabla} u_2), \end{gather}
\begin{gather}\boldsymbol{\nabla} u_1 = \frac{\boldsymbol e_2\times\boldsymbol e_3}{J},\quad \boldsymbol{\nabla} u_2 = \frac{\boldsymbol e_3\times\boldsymbol e_1}{J},\quad \boldsymbol{\nabla} u_3 = \frac{\boldsymbol e_1\times\boldsymbol e_2}{J}. \end{gather}Appendix B. Useful forms of $J$ and $B$
Taking 
$B^2 = {{\boldsymbol B}}_{\mathrm {cov}} \boldsymbol {\cdot } {{\boldsymbol B}}_{\mathrm {con}}$ gives
\begin{equation} \bar{J} \bar{B}^2 = \bar{G} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \bar{I},\end{equation}
where we recall the definition 
$\bar {B} = B/(2\psi _b)$. Taking 
$B^2 = {{\boldsymbol B}}_{\mathrm {con}} \boldsymbol {\cdot } {{\boldsymbol B}}_{\mathrm {con}}$ gives
\begin{equation} \bar{J}^2 \bar{B}^2 = |\boldsymbol e_3 + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \boldsymbol e_2|^2 = g_{33} + 2{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} g_{23} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}^2 g_{22}.\end{equation}Using (B 1) and (B 2) we can express the magnetic field strength ‘locally’ (in terms of only surface metrics)
\begin{equation} \frac{(\bar{G}+{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \bar{I})^2}{\bar{B}^2} = g_{33} + 2 {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} g_{23} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}^2 g_{22}.\end{equation}Using (B 1) and (B 2) we can also express the Jacobian locally
\begin{equation} \bar{J} (\bar{G} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \bar{I}) = g_{33} + 2 {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} g_{23} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}^2 g_{22}.\end{equation}Appendix C. The MHD constraint
One can write the constraint equation (3.11), in different ways. Taking the 
$\boldsymbol e_1, \boldsymbol e_2$, and 
$\boldsymbol e_3$ components of this equation gives
\begin{gather} g_{13} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} g_{12} = \bar{K} \bar{J}, \end{gather}
\begin{gather}g_{23} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} g_{22} = \bar{I} \bar{J}, \end{gather}
\begin{gather}g_{33} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} g_{23} = \bar{G} \bar{J}. \end{gather}Note that (C 2) and (C 3) only involve surface metrics, and may be combined, eliminating the Jacobian, to obtain the local MHD constraint
\begin{equation} \bar{I} (g_{33} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} g_{23}) = \bar{G} (g_{23} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} g_{22}).\end{equation}
Combining the 
$\boldsymbol e_i$ components we can derive three metric constraints not explicitly involving the Jacobian 
$J$, the one above and the following two
\begin{gather} \bar{I} (g_{13} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} g_{12}) = \bar{K} (g_{23} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} g_{22}), \end{gather}
\begin{gather}\bar{K} (g_{33} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} g_{23}) = \bar{G} (g_{13} + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} g_{12}). \end{gather}Note that this system is incomplete for a vacuum field because then 
$K = 0$ and 
$I = 0$ and (C 5) provides no information. It can then be completed by including (C 3).
Appendix D. Zeroth-order axisymmetric equilibrium by direct formulation
Here we consider magnetic coordinates 
$\psi , \theta , \varphi$ as functions of cylindrical coordinates 
$R, Z$ and 
$\phi$. The condition of axisymmetry can be stated as the condition that the 
$\hat {\boldsymbol R}(\phi ), \hat {\boldsymbol {\phi }}(\phi )$ and 
$\hat {\boldsymbol z}$ components of 
${{\boldsymbol B}}$ are independent of 
$\phi$. This implies that the magnetic surfaces must be themselves axisymmetric, so 
$\hat {\boldsymbol {\phi }}\boldsymbol {\cdot } \boldsymbol {\nabla } \psi = 0$, i.e. 
$\partial _\phi \psi = 0$. Then from 
$\partial _\phi (\hat {\boldsymbol {\phi }}\boldsymbol {\cdot }{{\boldsymbol B}}_{\mathrm {cov}} )= 0$ we obtain
\begin{equation} G \frac{\partial ^2 \varphi}{\partial \phi^2} + I \frac{\partial^2 \theta}{\partial \phi^2} = 0. \end{equation}
Integrating, using 
$\varphi = \phi + \nu$, and periodicity in 
$\phi$ we obtain
\begin{equation} G\frac{\partial \nu}{\partial \phi} + I \frac{\partial \theta}{\partial \phi} = 0.\end{equation}
Now, taking the 
$\partial _\phi (\hat {\boldsymbol R}\boldsymbol {\cdot } {\boldsymbol B}_{\mathrm {con}} )= 0$, we likewise obtain
\begin{equation} \frac{\partial \theta}{\partial \phi} - {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \frac{\partial \nu}{\partial \phi} = 0.\end{equation}
Equations (D 2) and (D 3) are linearly independent (i.e. 
$G + {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} I \neq 0$ by (B 1)), so we have
\begin{equation} \frac{\partial \theta}{\partial \phi} = \frac{\partial \nu}{\partial \phi} = 0.\end{equation} For obtaining the Grad–Shafranov equation, it is convenient to use a mixed form of 
${{\boldsymbol B}}$, using the toroidal part of the covariant field and the non-toroidal part of the contravariant field
\begin{equation} {{\boldsymbol B}}_{\mathrm{mix}} = G \boldsymbol{\nabla} \phi + \boldsymbol{\nabla} \phi \times \boldsymbol{\nabla} \varPsi, \end{equation}
where 
$\varPsi (\psi )$ is the poloidal magnetic flux, and 
$\textrm {d} \varPsi / \textrm {d} \psi = {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$. Using this form, force balance and Ampere's law imply 
$\mu _0\boldsymbol {\nabla } \varPsi \, \textrm {d} p/ \textrm {d} \varPsi = (\boldsymbol {\nabla } \times {{\boldsymbol B}}_{\mathrm {mix}}) \times {{\boldsymbol B}}_{\mathrm {mix}}$, which immediately yields the Grad–Shafranov equation,
\begin{equation} \mu_0\frac{\textrm{d} p}{\textrm{d} \varPsi} = -\frac{1}{R^2}\left[G\frac{\textrm{d} G}{\textrm{d} \varPsi} + {\rm \Delta}^* \varPsi\right] \end{equation}where
\begin{equation} {\rm \Delta}^* \varPsi = R \frac{\partial }{\partial R}\left(\frac{1}{R}\frac{\partial \varPsi}{\partial R}\right) + \frac{\partial^2 \varPsi}{\partial Z^2}. \end{equation}
Although this is a complete specification of the axisymmetric field, we require the full set of coordinates to solve the perturbative problem, so we must develop the more general representations, i.e. 
${{\boldsymbol B}}_{\mathrm {cov}}$ and 
${{\boldsymbol B}}_{\mathrm {con}}$. The functions 
$\theta (R, Z)$ and 
$\nu (R,Z)$ (and 
$I$ and 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$ as functions of 
$\varPsi$) can be obtained from additional equations derived from components of the MHD constraint, 
${{\boldsymbol B}}_{\mathrm {con}} = {{\boldsymbol B}}_{\mathrm {cov}}$
\begin{equation} G(\psi) \boldsymbol{\nabla}\varphi + I(\psi) \boldsymbol{\nabla}\theta + K(\psi, \theta, \varphi) \boldsymbol{\nabla} \psi= \boldsymbol{\nabla}\psi \times \boldsymbol{\nabla}\theta - {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \boldsymbol{\nabla} \psi \times \boldsymbol{\nabla} \varphi.\end{equation} The 
$\hat {\boldsymbol {\phi }}$ component gives
\begin{equation} \frac{{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} G}{R} = \left\{ \theta - {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \nu, \varPsi \right\},\end{equation}
where 
$\left \{ A, B \right \} = \partial _R A \partial _Z B - \partial _R B \partial _Z A$, and the 
$\hat {\boldsymbol R}$ and 
$\hat {\boldsymbol Z}$ components can be combined to eliminate 
$K$, yielding
\begin{equation} \frac{1}{R} |\boldsymbol{\nabla}\varPsi|^2 = \left\{ I \theta + G \nu, \varPsi \right\}.\end{equation}
Then 
$\psi$ can be obtained from
\begin{equation} \frac{\textrm{d}\psi}{\textrm{d}\varPsi} = \frac{1}{{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}}.\end{equation}
Finally, 
$K$ can be obtained from the 
$\boldsymbol {\nabla } \psi$ component of (D 8) (i.e. the condition that 
${{\boldsymbol B}}_{\mathrm {cov}}$ has no component pointing out of the magnetic surface),
\begin{equation} K(R,Z) = \frac{1}{|\boldsymbol{\nabla}\psi|^2} ( G\boldsymbol{\nabla}\psi\boldsymbol{\cdot}\boldsymbol{\nabla}\nu + I \boldsymbol{\nabla}\psi\boldsymbol{\cdot}\boldsymbol{\nabla}\theta). \end{equation}D.1 Solving equations ( D 9) and ( D 10)
Now we would like to solve these equations for the unknowns 
$\nu$ and 
$\theta$ and 
$K$. The Boozer angle 
$\theta$ can be expressed as
\begin{equation} \theta = \vartheta + \lambda(R, Z), \end{equation}
in terms of the geometric poloidal angle 
$\vartheta$, defined in terms of the quadrant-specific 
$\arctan$ function as 
$\vartheta = \arctan (R-R_a, z - z_a)$, with 
$R_a$ and 
$z_a$ the coordinates of the magnetic axis. Defining the potentials 
$P = G \nu + I \lambda$ and 
$A = \lambda - {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} \nu$, we can obtain from (D 9) the equation
\begin{gather} \frac{\partial A}{\partial \vartheta} = -1 + \frac{{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} G}{R {\mathcal{J}}(\vartheta)}, \end{gather}
\begin{gather}\frac{\partial P}{\partial \vartheta} = -I + \frac{|\boldsymbol{\nabla} \varPsi|^2}{R {\mathcal{J}}(\vartheta)}, \end{gather}
where 
${\mathcal {J}}(\vartheta ) = -\hat {\boldsymbol {\phi }} \boldsymbol {\cdot } (\boldsymbol {\nabla } \vartheta \times \boldsymbol {\nabla }\varPsi ) = \left \{ \vartheta , \varPsi \right \}$. Note that, formally, we are changing coordinates to 
$\vartheta$ and 
$\varPsi$. The functions 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota}$ and 
$I$ may be found as solubility constraints of these two equations
\begin{gather} {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} = 2{\rm \pi} \left( \int_0^{2{\rm \pi} } \, \textrm{d} \vartheta \frac{G}{R {\mathcal{J}}(\vartheta)} \right)^{-1}, \end{gather}
\begin{gather}I = \frac{1}{2{\rm \pi} }\int_0^{2{\rm \pi} } \, \textrm{d} \vartheta \frac{|\boldsymbol{\nabla} \varPsi|^2}{R {\mathcal{J}}(\vartheta)}, \end{gather}
and the solutions for 
$A$ and 
$P$ given as
\begin{gather} A = \int_0^{\vartheta} \, \textrm{d} \vartheta^\prime \left( \frac{{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} G}{R {\mathcal{J}}(\vartheta^\prime)} - 1 \right), \end{gather}
\begin{gather}P = \int_0^{\vartheta} \, \textrm{d} \vartheta^\prime \left( \frac{|\boldsymbol{\nabla} \varPsi|^2}{R {\mathcal{J}}(\vartheta^\prime)} - I\right), \end{gather}
where note that 
$R$ and 
$|\boldsymbol {\nabla }\varPsi |^2$ are also evaluated at 
$\vartheta ^\prime$, and the choice 
$P=A=0$ at 
$\vartheta =0$ has been made. From these solutions, the desired potentials 
$\lambda$ and 
$\nu$ may be obtained as
\begin{gather} \lambda = \frac{{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} P + G A}{{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} I + G}, \end{gather}
\begin{gather}\nu = \frac{P - I A}{{\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} I + G}, \end{gather}
where we again note that 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} I + G \neq 0$ by (B 1). Finally, to obtain 
${\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} , I$ and 
$G$ as functions of 
$\rho$, we solve (D 11), and the functions 
$R(\rho , \theta )$ and 
$z(\rho , \theta )$, are obtained by inverting 
$\theta (R, z)$ and 
$\varPsi (R,z)$.
Restoring subscripts and normalizations, we are thus furnished with the functions 
$R_0(\rho , \theta ), z_0(\rho , \theta ), \nu _0(\rho , \theta ), \bar {K}_0(\rho , \theta ), {\style{display: inline-block; transform: rotate(31deg)}{\raise1.5pt{\tiny{/}}}\kern-1.4pt\iota} _0(\rho ), \bar {I}(\rho )$ and 
$\bar {G}(\rho )$ to substitute into the expressions provided in § 3.1, and proceed to the calculations, at next order, in § 3.2.
Open access
