2.1. The equilibrium problem

We define an equilibrium solution as a magnetic field with nested flux surfaces that satisfies the MHS equilibrium equation $\boldsymbol{j}\times \boldsymbol{B}=\boldsymbol{\nabla }p$ , where $\boldsymbol{j}=\mu _0^{-1}\boldsymbol{\nabla }\times \boldsymbol{B}$ is the plasma current density, $p$ is the plasma pressure, and $\mu_0$ is the vacuum permeability. The magnetic field must, of course, be a solenoidal vector field, $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{B}=0$ , and we define a toroidal magnetic flux $2\pi \psi$ whose level sets define nested flux surfaces satisfying $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi =0$ . Both the global equilibrium solver and the near-axis description attempt the construction of equilibrium fields of this form.

Within this set of equations, we confer on the latter two a more prominent role. That is, we prioritise the magnetic field being solenoidal and having nested flux surfaces, which we impose on the field in an exact form. This is done by adopting the inverse-coordinate form of the problem (Bauer, Betancourt & Garabedian Reference Bauer, Betancourt and Garabedian1978; Hirshman & Whitson Reference Hirshman and Whitson1983). Instead of describing the magnetic field $\boldsymbol{B}$ explicitly as a function of space, we instead focus on describing the magnetic flux surfaces on which the field lives, as well as the form in which field lines wrap over them.

To construct such a description, we start by writing the magnetic field in its Clebsch form (D’haeseleer et al. Reference D’haeseleer, Hitchon, Callen and Shohet2012)

(2.1) \begin{equation} \boldsymbol{B}=\boldsymbol{\nabla }\psi \times (\boldsymbol{\nabla }\vartheta -\iota \boldsymbol{\nabla }\varphi ), \end{equation}

where $\vartheta$ and $\varphi$ are poloidal and toroidal angles (in straight field line coordinates with non-vanishing Jacobian), and $\iota$ (a function of $\psi$ ) is the rotational transform. A field written in this form satisfies both field requirements exactly. Defining a position vector $\boldsymbol{x}=\boldsymbol{x}(\psi ,\vartheta ,\varphi )$ , and taking that set of straight field line coordinates as our set of independent coordinates, the magnetic field may then be expressed as

(2.2) \begin{equation} \boldsymbol{B}=\mathcal{J}_\vartheta ^{-1}\left (\frac {\partial \boldsymbol{x}}{\partial \varphi }+\iota \frac {\partial \boldsymbol{x}}{\partial \vartheta }\right )\!, \end{equation}

where $\mathcal{J}_\vartheta =\partial _\psi \boldsymbol{x}\times \partial _\vartheta \boldsymbol{x}\boldsymbol{\cdot }\partial _\varphi \boldsymbol{x}$ is the Jacobian of the straight field line coordinate system. Full knowledge of flux surfaces in straight field line coordinates (i.e. $\boldsymbol{x}$ ), along with the rotational transform, uniquely defines the magnetic field. Alternatively, one could provide information about the average toroidal current rather than the rotational transform (see Appendix A).

The central task is then to find a function $\boldsymbol{x}(\psi ,\vartheta ,\varphi )$ , which we may represent in a multitude of ways. The NAE, in its inverse-coordinate form pioneered by Garren & Boozer (Reference Garren and Boozer1991c ), takes as a reference of this $\boldsymbol{x}$ the magnetic axis. This is a natural choice given that the NAE is concerned with providing a consistent asymptotic description of the field near the magnetic axis (i.e. $\psi =0$ ). Explicitly,

(2.3) \begin{equation} \boldsymbol{x}(\psi ,\vartheta ,\varphi )=\boldsymbol{r}_0+X(\psi ,\vartheta ,\varphi )\hat {\boldsymbol \kappa }+Y(\psi ,\vartheta ,\varphi )\hat {\boldsymbol \tau }+Z(\psi ,\vartheta ,\varphi )\hat {\boldsymbol{t}}, \end{equation}

where $\boldsymbol{r}_0$ describes the magnetic axis, and the vector triad $\{\hat {\boldsymbol \kappa },\hat {\boldsymbol \tau },\hat {\boldsymbol{t}}\}$ corresponds to the Frenet–Serret basis (Frenet Reference Frenet1852; Animov Reference Animov2001) (namely, the normal, binormal and tangent vectors defined with respect to the magnetic axis). Here, we specialise to regular magnetic axes with no, or only few, isolated flattening points (i.e. vanishing curvature points) in which such a frame (or its signed version (Carroll, Köse & Sterling Reference Carroll, Köse and Sterling2013; Plunk et al. Reference Plunk, Landreman and Helander2019; Camacho Mata & Plunk Reference Camacho Mata and Plunk2023; Rodríguez et al. Reference Rodríguez, Plunk and Jorge2025)) exists. This is enough to describe the asymptotic behaviour of all optimised configurations; namely, quasisymmetric and quasi-isodynamic ones (Gori, Lotz & Nührenberg Reference Gori, Lotz and Nührenberg1997; Plunk et al. Reference Plunk, Landreman and Helander2019). In the NAE, then, finding a magnetic field corresponds to finding an asymptotic form of the functions $X,\,Y$ and $Z$ , which describe the shape of flux surfaces, alongside the rotational transform and the axis shape. This description is performed in a special class of straight field line coordinates, Boozer coordinates $\{\psi ,\vartheta _{B},\varphi _{B}\}$ , which eases the treatment of optimised stellarators due to the simple involvement of $|\boldsymbol{B}|$ (Boozer Reference Boozer1981), which imposes additional constraints in the construction.

In the case of the global equilibrium description, the numerical solver DESC, which we shall focus on, describes flux surfaces not with respect to the axis (unlike new developments (Hindenlang et al. Reference Hindenlang, Plunk and Maj2025)), but in cylindrical coordinates so that

(2.4) \begin{equation} \boldsymbol{x}_{\mathrm{DESC}}=R \hat {\boldsymbol{R}} + Z \hat {\boldsymbol{Z}}. \end{equation}

It is then the task of the solver to find the appropriate $R$ and $Z$ functions as a function of $\{\psi ,\,\vartheta ,\,\varphi \}$ . Given that the description is in cylindrical coordinates, the toroidal angle is naturally chosen to be the cylindrical one, $\phi$ .Footnote ¹ Because we insist on this specific form of the toroidal angle, it requires a very particular choice of the poloidal angle $\vartheta$ in order to guarantee the coordinate system to be a straight field line coordinate system. This coordinate system is known as PEST coordinates $\{\psi ,\vartheta _{\mathrm{PEST}},\phi \}$ . However, the global solver does not know at the starting point what $\vartheta _{\mathrm{PEST}}$ is, as this depends on the final magnetic field solution. Hence, in practice, the solver must introduce a stream function $\lambda$ such that $\vartheta _{\mathrm{PEST}}=\theta +\lambda$ , where $\theta$ can be any arbitrary poloidal angle used by the equilibrium solver. In that case, the magnetic field, as in (2.2) is written as

(2.5a) \begin{equation} \boldsymbol{B}=\mathcal{J}_\theta ^{-1}\left [(1+\partial _\theta \lambda )\frac {\partial \boldsymbol{x}}{\partial \phi }+(\iota -\partial _\phi \lambda )\frac {\partial \boldsymbol{x}}{\partial \theta }\right ]\!. \end{equation}

Where now the Jacobian is of the computational coordinates $\mathcal{J}_\theta =\partial _\psi \boldsymbol{x}\times \partial _\theta \boldsymbol{x}\boldsymbol{\cdot }\partial _\varphi \boldsymbol{x}\!$ . In the case of the global solver then, a magnetic field is uniquely described by finding functions $R,\,Z$ and $\lambda$ as a function of $\{\psi ,\,\theta ,\,\phi \}$ .

So far, we have only represented a magnetic field but have not yet solved the equilibrium. We must still find a consistent set of flux surfaces $\boldsymbol{x}$ (described in appropriate coordinates), such that the resulting magnetic field is in equilibrium. We must deform flux surfaces until the field comes to equilibrium; i.e. it minimises the force residual $\boldsymbol{F}=\boldsymbol{j}\times \boldsymbol{B}-\boldsymbol{\nabla }p$ for a prescribed function $p(\psi )$ . In the context of the NAE, the consistent construction of flux surfaces and Boozer coordinates can be carried out systematically by setting $\boldsymbol{F} = \boldsymbol{0}$ order by order, and detailed descriptions of how to do this may be found in Landreman & Sengupta (Reference Landreman and Sengupta2019) and Rodríguez et al. (Reference Rodríguez, Plunk and Jorge2025). We content ourselves with knowing that, at the end of the day, we will be able to have an asymptotic form of a $X,Y,Z$ rotational transform and the axis shape.

The equilibrium solver DESC seeks equilibria by numerically and iteratively minimising the force residual $\boldsymbol{F}$ directly. Note that this is not the only existing practical approach, with the important case of the VMEC code approaching the problem by minimising energy (Hirshman & Whitson Reference Hirshman and Whitson1983). To be more explicit about the force residual, consider first the covariant form of the magnetic field

(2.6) \begin{equation} \boldsymbol{B}=B_\theta \boldsymbol{\nabla }\theta +B_\phi \boldsymbol{\nabla }\phi +B_\psi \boldsymbol{\nabla }\psi . \end{equation}

There might appear to be new information in the problem, but all of the covariant components in (2.6) can be directly computed from the geometry by simply leveraging the dual relations (D’haeseleer et al. Reference D’haeseleer, Hitchon, Callen and Shohet2012, (2.3.12)–(2.3.13)) $\boldsymbol{\nabla }q_i=\epsilon _{ijk}\partial _{q_j}\boldsymbol{x}\times \partial _{q_k}\boldsymbol{x}$ and (2.5a ). With this, the force residual is (see Panici et al. (Reference Panici, Conlin, Dudt, Unalmis and Kolemen2023) as well)

(2.7) \begin{align} \boldsymbol{F}&=\left [\mathcal{J}_\theta ^{-1}\mu _0^{-1}(\partial _\phi B_\psi -\partial _\psi B_\phi )(1+\partial _\theta \lambda )-\mathcal{J}_\theta ^{-1}\mu _0^{-1}(\iota -\partial _\phi \lambda )(\partial _\psi B_\theta -\partial _\theta B_\psi )-p'\right ]\boldsymbol{\nabla }\psi \nonumber \\& \quad -\mathcal{J}_\theta ^{-1}\mu _0^{-1}(\partial _\theta B_\phi -\partial _\phi B_\theta )((1+\partial _\theta \lambda )\boldsymbol{\nabla }\theta -(\iota -\partial _\phi \lambda )\boldsymbol{\nabla }\phi ). \end{align}

The $\boldsymbol{\nabla }\psi$ component can be interpreted as the radial force balance equation, with the helical component representing $\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi$ . Thus, the global solver is set to find $\{R,\,Z,\,\lambda \}$ as a function of $\{\psi ,\theta ,\phi \}$ that make (2.7) vanish.

The near-axis and global treatments of the problem thus have a key difference not only in the coordinate representation used, but also in their treatment of boundary conditions. One must provide certain inputs to the equilibrium problem to describe different equilibria. In the NAE case, this prescription occurs outwards; that is, one must provide the shape of the axis, and then, order by order, some features of the flux surfaces away from it. Which features must be provided depends on the type of stellarator under consideration, and the order of the expansion of interest (see Landreman & Sengupta Reference Landreman and Sengupta2019; Rodríguez et al. Reference Rodríguez, Plunk and Jorge2025). The global solution fundamentally attempts the opposite: one specifies the shape of the boundary, and the solution is constructed inwards (what is known as a fixed-boundary solution). This is most straightforwardly understood by imagining a vacuum field in the form of Laplace’s equation. From this difference it should be clear that there is not a one-to-one correspondence between the two approaches. A finite truncated near-axis construction does not, in principle, uniquely describe a global field solution, but rather simply constrains its core behaviour. We are now in a position to connect both.

2.2. Near-axis expansion in DESC

In the NAE we learned that the description of flux surfaces $\boldsymbol{x}(\psi ,\vartheta _B,\varphi _B)$ is given by $\{X,Y,Z\}$ , (2.3). The forms of $\{X,Y,Z\}$ are particularly simple, as may be directly appreciated by considering their Taylor–Fourier expansion (Garren & Boozer Reference Garren and Boozer1991c ; Landreman & Sengupta Reference Landreman and Sengupta2019)

(2.8) \begin{equation} f(\psi ,\theta _{B},\varphi _{B})=\sum _{l=0}^\infty r^l\sum _{m=0}^l{}^{\prime }\left (f_{lm}^C(\varphi _{B})\cos m\vartheta _{B}+f_{lm}^S(\varphi _B)\sin m\vartheta _{B}\right )\!, \end{equation}

where $r=\sqrt {2\psi /\bar {B}}$ is a pseudo-radial coordinate, $\bar {B}$ is a representative normalisation value of the magnetic field and the sum $\sum '$ denotes sum over even or odd orders depending on the parity of $l$ . Because the poloidal dependence grows in a controlled manner linked to the powers of $r$ , the flux surfaces gain complexity order by order; they start by being elliptical, then acquire triangularity, and so on.

To impose this asymptotic form of the field on the global equilibrium solver, we must be able to relate this asymptotic description to the representation in DESC, and constrain the latter according to the former. The structure in (2.8) is a consequence of maintaining the field description singularity free close to the magnetic axis ( $\psi =0$ ) (Kuo-Petravic & Boozer Reference Kuo-Petravic and Boozer1987; Garren & Boozer Reference Garren and Boozer1991b ; Landreman & Sengupta Reference Landreman and Sengupta2018). This requirement was indeed observed in the implementation of DESC, where the global solver employs a Zernike basis to represent $\{R,Z,\lambda \}$ , crucially coupling the powers of DESC’s chosen pseudo-radial coordinate $\rho =\sqrt {\psi /\psi _{\mathrm{edge}}}$ to the poloidal modes $\theta$ . This makes the treatment in DESC close to the form used in the near-axis description.

To make that connection, however, we must consider the differences carefully. One is the use of different re-scaled pseudo-radial coordinates; we may relate $\rho$ to the NAE $r$ by $\rho =r\sqrt {\bar {B}/2\psi _{\mathrm{edge}}}$ . The second point of difference is the use of a different poloidal angle. While the near-axis description uses the poloidal Boozer angle, $\vartheta _{B}$ , the global equilibrium is described in a general poloidal angle $\theta$ . The simplest comparison between the two descriptions is then made when $\theta$ is forced to match the Boozer angle, at least to leading order near the axis. When imposing the near-axis constraints on $\theta$ as if it were Boozer coordinates, the solver is forced to use a very particular form of $\lambda$ , and thus can be overly constraining. Defining $\nu =\varphi _{B}-\phi$ as the function that maps the Boozer toroidal angle to the cylindrical one, our choice of $\theta$ dictates that

(2.9) \begin{gather} \lambda =-\iota \nu . \end{gather}

This simplest choice of poloidal angle could be modified to accommodate a more general and efficient $\theta$ . Details about how this would be incorporated are shown in Appendix B, but we shall nevertheless consider the simplest scenario in the main body of the paper, which should suffice as a proof of principle.

With the coordinates in place, what remains from the comparison is the different radial–poloidal structure of the NAE expansion and the Zernike basis. In the equilibrium solver a generic function is written as (ignoring the toroidal coordinate for now)

(2.10) \begin{equation} f(\rho ,\theta )=\sum _{l=0}^\infty \sum _{m=-l}^l{}^{\prime }f_{lm}\mathcal{Z}_l^m(\rho ,\theta ), \end{equation}

where the $\sum ^{\prime }$ again denotes a sum over the same parity as $l$ , and the Zernike polynomials are

(2.11) \begin{equation} \mathcal{Z}_l^m(\rho ,\theta )=\mathcal{R}_l^{|m|}\begin{cases} \cos (m\theta ) & m\geqslant 0, \\ \sin (|m|\theta ) & m\lt 0, \end{cases} \end{equation}

where the radial parts are the shifted Jacobi polynomials, given as

(2.12) \begin{equation} \mathcal{R}_l^{|m|}=\sum _{k=0}^{(l-|m|)/2}\frac {(-1)^k(l-k)!}{k!\left ((({l+m})/{2})-k\right )!\left ((({l-m})/{2})-k\right )!}\rho ^{l-2k}. \end{equation}

The coefficients $\{f_{lm}\}$ for $\{R,Z,\lambda \}$ are the set which DESC solves for in order to construct the equilibrium field. To make a one-to-one connection to the NAE, (2.8), we need to unravel this sum and pick the different powers of $\rho$ . Generally, a function $f$ defined by $f_{lm}$ coefficients can be written in terms of the equivalent near-axis coefficients as

(2.13) \begin{equation} f_{l|m|}^{C/S}=\left (\frac {\rho }{r}\right )^l\times \begin{cases} \begin{aligned} &(-1)^{l/2}\sum _{k=l/2}^\infty \Bigg [ (-1)^k\binom {k+l/2}{l}\\&\hspace {1in} \times \binom {l}{(l-|m|)/2}f_{2k,m}\Bigg ]\quad |m| \mathrm{\,is\,even} ,\\ &(-1)^{(l+1)/2}\sum _{k=(l+1)/2}^\infty \Bigg [(-1)^k\binom {k+(l-1)/2}{l}\\&\hspace {1in}\times \binom {l}{(l-|m|)/2}f_{2k-1,m}\Bigg ] \quad |m| \mathrm{\,is\,odd}, \\ \end{aligned} \end{cases} \end{equation}

and the sign of $m$ has been taken to represent, as is the case in the DESC notation, the sine or cosine terms $S/C$ (negative and positive respectively).Footnote ² The above only holds for $l\geqslant |m|$ , and $m$ sharing the same parity as $l$ . Otherwise, the contribution will be zero. Note how, for large $k$ at constant $l$ and $m$ , the factor in front of the Zernike components goes like $\sim k^l$ . This shows the exponential contribution from the higher-order Zernike modes, which must therefore be conveniently truncated. See details on the derivation in Appendix C. Thus, imposing the near-axis behaviour constitutes imposing constraints on linear combinations of the DESC degrees of freedom.

As a final remark, we note that what have been referred to as coefficients $f_{lm}$ so far are generally functions of $\phi$ . In the pseudo-spectral representation of DESC, we may then treat the toroidal dependence with a Fourier series in $\phi$ , with the $n$ Fourier components corresponding to $f_{lmn}$ , with $n\geqslant 0$ representing cosines, and $n\lt 0$ sines.

3.1. Zeroth order: the magnetic axis

To leading order in the near-axis description, flux surfaces reduce to a single closed curve at the origin ( $\psi =0$ ):Footnote ³ the magnetic axis. In practice, such a closed curve is described as a set of Fourier harmonics (in $\phi$ ) in cylindrical coordinates $\{R,\,Z\}$ . That is,

(3.1a) \begin{align} R(\phi ) & =R_0+\sum _{n=1}^N \left(R_{n}^C\cos n\phi +R_{n}^S\sin n\phi \right)\!,\\[-12pt]\nonumber \end{align}

(3.1b) \begin{align} Z(\phi ) & =\sum _{n=1}^N \left(Z_{n}^C\cos n\phi +Z_n^S\sin n\phi \right)\!,\end{align}

where the Fourier coefficients are constant. This form of describing the axis is natural to NAE codes (such as pyQSC), and equilibrium codes (such as DESC). Following the discussion on the Zernike basis, using (C2), the relevant constraints for the Zernike modes are then

(3.2a) \begin{align} R_{n}^{C/S}=\sum _{k=0}^\infty (-1)^k R_{2k,0,\pm |n|},\\[-10pt]\nonumber\end{align}

(3.2b) \begin{align} Z_{n}^{C/S}=\sum _{k=0}^\infty (-1)^k Z_{2k,0,\pm |n|}.\end{align}

We shall thus have $2(2N+1)$ constraints to impose onto the Zernike harmonics in order to constrain the magnetic axis to match that of the NAE. No further geometric consideration is needed here.

3.2. First order: elliptical surfaces

At first order things start to become more interesting, and we need to explicitly transform the constraints into cylindrical coordinates (the language of DESC) using the information in our magnetic-axis frame. This includes transforming the Boozer toroidal angle to the cylindrical angle $\phi$ . This transformation is routinely performed numerically by codes like pyQSC to represent near-axis magnetic fields in the laboratory frame, and was previously explored in Landreman & Sengupta (Reference Landreman and Sengupta2018). Here, we present an alternative simple geometric treatment of the problem.

Start by defining the flux surface as $\boldsymbol{x}=\boldsymbol{r}_0+\boldsymbol{x}_1$ , where $\boldsymbol{x}_1\propto \rho$ and thus is in the asymptotic sense small. That is, we construct flux surfaces near the axis through a small displacement $\boldsymbol{x}_1$ from it, which is a function of $\rho ,\,\theta$ and $\phi$ . Our task is to find $R$ and $Z$ corresponding to these points. Start with $R$ , the major radius, and consider the projection of the problem to a constant $Z$ plane, which we define as the $\pi$ plane (see the diagram in figure 1). Define a point on the magnetic axis by $R_0$ (the radial distance to a point along the axis), and $\phi _0$ (the cylindrical angle of said point). The displacement $\boldsymbol{x}_1$ is then a function of that point (see figure 1), whose projection normal to $z$ we define to be $\boldsymbol{x}_1^\pi$ . Define $\phi _1$ to be the angle that $\boldsymbol{x}_1^\pi$ makes with $\hat {\boldsymbol{R}}_0$ , which satisfies $\cos \phi _1=\hat {\boldsymbol{R}}_0\times\boldsymbol{x}_1^\pi /|\boldsymbol{x}_1^\pi |$ . Here, $\hat {\boldsymbol{R}}_0$ corresponds to the radial vector direction evaluated at the axis; that is, $\hat {\boldsymbol{R}}_0=\hat {\boldsymbol{R}}(\phi =\phi _0$ ). We then define $\delta \phi$ to be the angle which, when added to $\phi _0$ , results in the cylindrical toroidal angle of the point on the flux surface $\boldsymbol{x}_0+\boldsymbol{x}_1$ .

Figure 1.

Diagram illustrating the key element for the geometric transformation at first order. Schematic diagram showing the position of a point on the surface (at radial distance $R$ and angle $\phi _0+\delta \phi$ ), in reference to other quantities. These include the position along the magnetic axis (radial position $R_0$ and angle $\phi _0$ ), and $\rho \boldsymbol{x}_1$ from the axis to the point, projected onto the $R,\,\phi _c$ plane (superindex $\pi$ ).

Application of the sine rule gives

(3.3) \begin{equation} \delta \phi \approx \frac {x_1^\phi }{R_0}, \end{equation}

where $x_1^\phi =|\boldsymbol{x}_1^\pi |\sin \phi _1$ and only the leading piece in $\rho$ is being kept (so, $R \approx R_0$ ) and $\delta \phi \ll 1$ has been assumed. We define $x_1^j=\boldsymbol{x}_1\boldsymbol{\cdot }\hat {j}$ for $j=R,\,\phi ,\,z$ as the projections of $\boldsymbol{x}_1$ along the cylindrical basis vectors, which when expressed in terms of the Frenet–Serret frame of the axis will involve various projections of the triad $\{\hat {\boldsymbol \kappa },\hat {\boldsymbol \tau },\hat {\boldsymbol{t}}\}$ .

Applying the cosine rule, and defining $R-R_0=\delta R$ , we find that to leading order in $\rho$

(3.4) \begin{equation} \delta R\approx x_1^R. \end{equation}

Therefore, the radial position as a function of the geometric angle at the point on the surface,Footnote ⁴ is given as

(3.5) \begin{equation} R(\phi )\approx R_0(\phi -\delta \phi )+\delta R(\phi -\delta \phi )\approx R_0(\phi )+\left [\delta R(\phi )-\delta \phi \partial _\phi R_0\right ]+O(\rho ^2), \end{equation}

where all the quantities in brackets are evaluated at $\phi$ , and $R_0$ is a known function (from the shape of the magnetic axis). The quantity in square brackets is what will become the constraint on the Zernike basis, and in the asymptotic sense is correct to $O(\rho )$ .

Then, plugging in our expressions for $\delta \phi$ and $\delta R$ , the first-order $\rho$ terms in $R$ are

(3.6) \begin{equation} R_1=x_1^R-x_1^\phi \frac {\partial _\phi R_0}{R_0}. \end{equation}

This form preserves the simple $\theta$ -harmonic content of $\boldsymbol{x}_1$ , which from the near axis has the simple form $\boldsymbol{x}_1=\cos \theta (X_1^c\hat {\boldsymbol \kappa }+Y_1^c\hat {\boldsymbol \tau })+\sin \theta (X_1^s \hat {\boldsymbol \kappa } + Y_1^s\hat {\boldsymbol \tau } )$ . Thus, (3.6) can be written as $R_1=\mathcal{R}_{1,1}(\phi )\cos \theta +\mathcal{R}_{1,-1}(\phi )\sin \theta$ , where the form of the coefficients follows directly from (3.6) and the form of $\boldsymbol{x}_1$ (these are derived explicitly in Appendix D). The subscripts in $\mathcal{R}_{nm}$ refer to the radial and poloidal orders of the term, respectively. To evaluate these we need the projections of the axis normal and binormal onto the cylindrical basis, information which is known in the context of the NAE (numerically calculated by codes like pyQSC).

For $Z$ , analogously, we get

(3.7) \begin{equation} Z_1 = x_1^z-x_1^\phi \frac {\partial _\phi Z_0}{R_0}, \end{equation}

where $Z_0(\phi )$ is the $Z$ of the magnetic axis that is known fully as a function of $\phi$ . Once again, one may collect terms to construct $Z_1=\mathcal{Z}_{1,1}(\phi )\cos \theta +\mathcal{Z}_{1,-1}(\phi )\sin \theta$ .

Using (C4), we may (using DESC notation) write

(3.8a) \begin{align} \frac {r}{\rho }\mathcal{R}_{1,1,n} &= -\sum _{k=1}^M(-1)^kkR_{2k-1,1,n},\\[-10pt]\nonumber\end{align}

(3.8b) \begin{align} \frac {r}{\rho }\mathcal{R}_{1,-1,n} & = -\sum _{k=1}^M(-1)^kkR_{2k-1,-1,n},\end{align}

for all $n\in [-N,N]$ representing Fourier components in $\phi$ , and the negative sign corresponds to the sine $\phi$ components. The left-hand sides represent the NAE components, which may be obtained through the NAE, and thus serve as constraints on the DESC Fourier–Zernike mode amplitudes, which appear on the right-hand side. The constraints for $Z$ have exactly the same form but with $\mathcal{Z}$ .

3.2.1. Straight field line angle

Recall that, in order to complete the description of a field in DESC, it is important to find the streamfunction $\lambda$ that transforms the angle $\theta$ to the PEST poloidal coordinate. In the approach here, and given that the near-axis expressions in Boozer coordinates are being used as constraints, we have (at least asymptotically) $\theta =\vartheta _B$ , and thus $\lambda =-\iota \nu$ , (2.9). Hence, by constraining $R$ and $Z$ as above, and in order for DESC to find an equilibrium, $\lambda$ is in practice constrained.

To find what this $\lambda$ is expected to be, we need to find $\phi$ , the cylindrical angle, as a function of $\phi _{B}$ and $\vartheta _{B}$ within the NAE. We may write $\phi =\phi ^{(0)}(\varphi _B)+r\phi ^{(1)}(\vartheta _B,\varphi _B)+\cdots$ , in which the leading-order expression is (Landreman & Sengupta Reference Landreman and Sengupta2019 (A20))

(3.9) \begin{equation} \lambda ^{(0)}=-\iota _0\int _0^\phi \left (\frac {B_0}{G_0}\frac {\mathrm{d}\ell }{\mathrm{d}\phi }-1\right )\mathrm{d}\phi ', \end{equation}

where $\ell$ is the length along the magnetic axis, $B_0$ the magnetic field magnitude, $G_0$ the poloidal Boozer current and $\iota _0$ the rotational transform, all the zero subscripts indicating evaluation on axis.

The first correction $\phi ^{(1)}$ is directly related to $\delta \phi$ in (3.3), so that $\nu ^{(1)}=-\delta \phi _1=-x_1^\phi /R_0$ , parametrised with the cylindrical angle on axis. Parametrised in terms of the standard cylindrical angle it would read

(3.10) \begin{equation} \lambda ^{(1)}=\iota _0\frac {x_1^\phi }{R_0}\left (1-\frac {\lambda ^{(0)}{}'}{\iota _0}\right )\!, \end{equation}

where the prime denotes a toroidal derivative.

3.2.2. Interpretation of the first-order transformation

The above might appear to be a purely formal artefact resulting from the change of frame, but it has a straightforward geometric interpretation. The expressions in (3.6) and (3.7) show the difference in describing the shape of cross-sections in the near-axis frame (cutting flux surfaces normal to the magnetic axis) and the ‘laboratory frame’ (where the cuts are made at constant cylindrical toroidal angle). The elliptical cross-sections that live in the plane normal to the axis will generally have some inclination with respect to the cylindrical basis, by virtue of the axis being non-planar. Hence, a cut in the laboratory frame introduces an additional projection.

Figure 2.

Definition of the slant angle $\alpha$ . Diagram showing the definition of the angle $\alpha$ measuring the inclination of the magnetic axis at the origin ( $\phi =0$ ) with the ‘laboratory’ cylindrical coordinate system. The symbols have their usual meaning.

Take as an illustrating example the cross-section at a stellarator-symmetric point (Dewar & Hudson Reference Dewar and Hudson1998), which corresponds to an up–down-symmetric ellipse in the plane normal to the magnetic axis. Because under the map $(\phi ,\theta )\rightarrow (-\phi ,-\theta )$ it must be the case that $R,Z\rightarrow R,-Z$ , it follows that $R$ and $Z$ must be even and odd functions, respectively. In particular, this means that $\partial _\phi R_0=0=\partial _\phi ^2 Z_0$ at $\phi =0$ . Define a slant angle $\alpha$ as the angle that the plane normal to the magnetic axis makes with respect to the cylindrical toroidal angle; i.e. the magnetic axis is rotated about $\hat {R}$ by an angle $\alpha$ (see figure 2). Then we may write the components of $\hat {\boldsymbol \tau }$ as $\tau ^z=\cos \alpha$ and $\tau ^\phi =\sin \alpha$ , and because of the symmetry, $\hat {\kappa }=-\hat {R}$ . The transformation in (3.6) and (3.7) then reduces to a scaling of the ellipse in the binormal ( $Y$ ) direction

(3.11) \begin{align} \bar {Y}_1=Y_{1}\bigg (\tau ^z-\underbrace {\frac {\partial _\phi Z_0}{R_0}}_{-\tau ^\phi /\tau ^z}\tau ^\phi \bigg )=\frac {Y_{1}}{\cos \alpha }. \end{align}

Due to the inclination of the axis, the ellipticity of the cross-section changes (see an example in figure 3), becoming stretched in the binormal direction. Away from this point of stellarator symmetry, the axis also presents an inclination with respect to the radial coordinate, leading to an additional reshaping of the elliptical shape.

Figure 3.

Geometric deformation of cross-sections from the near-axis to the laboratory frame. Example of the change in the cross-sections due to the geometric effects of going from the frame of the axis (broken lines) to the laboratory frame (solid line). These correspond to the cross-sections of the ‘precise QH’ stellarator configuration from Landreman & Paul (Reference Landreman and Paul2022) at one of its stellarator-symmetric points. The left corresponds to the cross-section evaluated at $r=0.01$ , while the shape to the right is at $r=0.05$ . The shape to the left shows the enhancement of elongation in the vertical direction due to the inclination of the axis, and the right the change of triangularity but immutability of the centre of the cross-section.

3.3. Second order: triangular surfaces

The same approach detailed in the previous section for the first order can be extended to higher orders following the same geometric construction as that in figure 1. The details of how to proceed to second order are presented in Appendix E. Second order is perhaps the highest order of interest, given that near-axis constructions are hardly performed beyond it (Garren & Boozer Reference Garren and Boozer1991a ; Landreman & Sengupta Reference Landreman and Sengupta2019; Rodríguez et al. Reference Rodríguez, Plunk and Jorge2025).

The derivation at second order proceeds by writing $\boldsymbol{x}=\boldsymbol{r}_0+\rho \boldsymbol{x}_1+\rho ^2\boldsymbol{x}_2$ , where $\boldsymbol{x}_2=X_2\hat {\boldsymbol \kappa }+Y_2\hat {\boldsymbol \tau }+Z_2\hat {\boldsymbol{t}}$ . Using the same considerations as in figure 1

(3.12a) \begin{align} R_2 & =-\frac {1}{2}\partial _\phi ^2R_0\left (\frac {x_1^\phi }{R_0}\right )^2-\left (\frac {x_2^\phi }{R_0}-\frac {x_1^Rx_1^\phi }{R_0^2}\right )\partial _\phi R_0\nonumber \\[5pt]& \quad -\frac {x_1^\phi }{R_0}\partial _\phi \left (x_1^R-\frac {x_1^\phi }{R_0}\partial _\phi R_0\right ) +\left (x_2^R+\frac {(x_1^\phi )^2}{2R_0}\right )\!,\\[-5pt]\nonumber \end{align}

(3.12b) \begin{align} Z_2 &=-\frac {1}{2}\partial _\phi ^2 Z_0\left (\frac {x_1^\phi }{R_0}\right )^2-\left (\frac {x_2^\phi }{R_0}-\frac {x_1^Rx_1^\phi }{R_0^2}\right )\partial _\phi Z_0-\frac {x_1^\phi }{R_0}\partial _\phi \left ( x_1^z-\frac {x_1^{\phi }}{R_0}\partial _\phi Z_0\right )+x_2^z,\end{align}

where we may express their poloidal harmonic content $R_2=R_{2,0}+R_2^c\cos 2\theta +R_2^s\sin 2\theta$ (and similarly for $Z$ ) by appropriate use of multiple angle formulae. Note that the complexity has significantly increased in going to higher order, which will tend to make higher harmonic toroidal content proliferate.

Once we have these geometric transformations, we may then impose them as constraints on the Zernike basis, by simply considering the general expression in (2.13) to write

(3.13a) \begin{align} f_{2,0}^{\mathrm{NAE}}=-\left (\frac {\rho }{r}\right )^2\sum _{k=1}^\infty (-1)^kk(k+1)f_{2k,0}, \end{align}

(3.13b) \begin{align} f_{2,\pm 2}^{\mathrm{NAE}}=-\left (\frac {\rho }{r}\right )^2\frac {1}{2}\sum _{k=1}^\infty (-1)^kk(k+1)f_{2k,\pm 2}, \end{align}

where the latter corresponds to the cosine piece for the + sign and the sine for −.

3.3.1. Interpretation of the second-order transformation

To illustrate the geometric meaning of this transformation, let us consider again the shape of the cross-section at the stellarator-symmetric point. At second order, the cross-section gains shaping beyond ellipticity (Landreman Reference Landreman2021; Rodríguez Reference Rodríguez2023). First, it acquires some level of triangularity, $\delta$ , a measure of left–right asymmetry of the cross-section, defined as the distance between the vertical turning points of the cross-section with respect to the centre along the symmetry line of the cross-section. In the near-axis description

(3.14) \begin{equation} \delta =2\left (\frac {Y_{2s}}{Y_{1s}}-\frac {X_{2c}}{X_{1c}}\right )\!, \end{equation}

where $Y$ and $X$ are along the normal and binormal, but an equivalent definition could be adapted to the shape with respect to the laboratory frame, with $Z$ instead of $Y$ and $R$ instead of $X$ .

The other important ingredient is the Shafranov shift, the relative displacement of the centres of the cross-sections from one surface to the next. Following Rodríguez (Reference Rodríguez2023), in the normal direction (in this case the radial one as well)

(3.15) \begin{equation} \Delta _X=X_{20}+X_{2c}. \end{equation}

The transformations in (3.12a ) and (3.12b ) for the stellarator-symmetric cross-section read

(3.16a) \begin{align} R_2=(X_{20}+\varLambda )+(X_{2c}-\varLambda )\cos 2\theta , \end{align}

(3.16b) \begin{align} Z_2 = Y_{20}\tilde {\tau }+ \left (Y_{2s}\tilde {\tau }-\tau ^\phi Y_{1s}\left [\frac {Y_{1c}'}{2R_0}\tilde {\tau }+\frac {X_{1c}}{2R_0}\left (\kappa _Z'-(\kappa _R+\kappa _\phi ')\frac {Z_0'}{R_0}\right )\right ]\right )\sin 2\theta , \end{align}

where $\tilde {\tau }=1/\tau ^z$ and $\varLambda =Y_{1s}^2\tau ^\phi [R_0(\tau ^\phi -2\tau _R')+\tau ^zR_0'']/4R_0^2$ . From the top equation, we may read the Shafranov shift $\Delta _R$ , which the transformation leaves unchanged, i.e. $\Delta _X=\Delta _R$ . The triangularity, however, does change (see Rodríguez Reference Rodríguez2024) provided $\alpha \neq 0$ and there is a rotation of the elliptical cross-section about the axis and across the stellarator-symmetric point. An example of this is illustrated in figure 3. There is as a result an offset in the value of triangularity when going from the axis to the laboratory frame. However, changing triangularity in the Frenet frame at second order directly affects triangularity in the laboratory frame linearly.

4.1. Verification tests between global solutions and NAE

Before proceeding to test the above, we introduce a number of measures to assess the correct behaviour of the constrained equilibrium. Because of the final step in the construction above, we shall achieve equilibrium to the desired level of accuracy, but the relaxation of the near-axis constraints may lead to deviations from its expected asymptotic behaviour if something went wrong. That is why it is fitting to check certain aspects of the field at the final step and compare those with the near-axis expectation. Depending on the order in the expansion in $\rho$ considered, different quantities may be considered.

To make the comparison as impartial as possible, we do not compare field quantities that correspond to the shape of flux surfaces which we are directly enforcing by the constraint, but rather, derived ones. We define the average discrepancies between on-axis quantities as

(4.2) \begin{equation} \unicode{x0394} f_0 = \left \langle \frac {\left |f_0^{\mathrm{DESC}} - f_0^{\mathrm{NAE}}\right |}{f^{\mathrm{NAE}}_0}\right \rangle _{\phi }, \end{equation}

where $\langle \ldots \rangle _\phi =\int _0^{2\pi }\ldots \mathrm{d}\phi /(2\pi )$ . The $\mathtt {DESC}$ quantities are defined as $f_0^{\mathrm{DESC}}=f^{\mathrm{DESC}}(\rho =0)$ .

For assessment of the first-order constrained fields, we consider $f = \iota ,\,B_0$ and $\lambda$ . In a correctly constrained equilibrium we expect $\unicode{x0394} \iota _0\rightarrow 0$ , as the rotational transform on axis is fully determined by the elliptical cross-sections about the axis (Mercier Reference Mercier1964; Landreman & Sengupta Reference Landreman and Sengupta2018). However, we must be careful here, because this is only correct when $\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi =0$ exactly. Because DESC only achieves this approximately, as part of the optimisation program, the rotational transform will have some dependence on higher-order shaping of $R_2$ and $Z_2$ (see Appendix A). Hence, $\unicode{x0394} \iota _0$ shall provide both a measure of the correct implementation of the constraints, as well as the compatibility with a global solution with $\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi =0$ . This potential dependence on higher orders does not apply to the magnetic field on axis, which the modulation of the ellipse areas at first-order controls. The third measure is a check on the behaviour of the $\theta$ coordinate in DESC, which by the constraints on the shaping, and in order to satisfy force balance, should tend towards becoming a Boozer coordinate. As shown in § 3.2.1, this asymptotic behaviour should manifest in the streamfunction, and in particular on the on-axis value, (3.9). We expect these normalised quantities then to be small in a successful implementation. To compare the verification errors with something, we shall present the NAE-constrained equilibrium errors alongside the errors found when constructing the global equilibrium not through the near-axis constraint, but when proceeding in the standard finite radius near-axis fixed-boundary equilibrium fashion (Landreman & Sengupta Reference Landreman and Sengupta2019).

At second order, we add the comparison of $f=V''$ , the so-called vacuum well of the configuration (Greene Reference Greene1998; Landreman & Jorge Reference Landreman and Jorge2020). This quantity reflects an important aspect of the MHD stability of the magnetic field, and it directly involves the second-order NAE. If a NAE was designed to be stable, we would like the global construction to preserve this curated feature, and thus it is natural to look at $\unicode{x0394} V_0''$ .

In addition to these scalar measures, we shall also consider an additional quantitative diagnostic away from the magnetic axis: the scaling of the magnetic field magnitude deviations from the ideal near-axis behaviour with $\rho$ . In particular, we shall look at the radial dependence of the quasisymmetry error (as measured by the normalised Boozer $f_B$ metric (Rodríguez et al. Reference Rodríguez, Paul and Bhattacharjee2022a ; Dudt et al. Reference Dudt, Conlin, Panici and Kolemen2023))Footnote ⁶ for quasisymmetric configurations, which we expect to scale (ideally) as $\mathcal{O}(\rho ^2)$ for the first-order constrained problem, and $O(\rho ^3)$ when constrained to second order. Note that this scaling will only be observed if the field is sufficiently quasisymmetric at second order, as one must generally allow for some deviations from exact quasisymmetric at second order (Garren & Boozer Reference Garren and Boozer1991a ; Landreman & Sengupta Reference Landreman and Sengupta2019).

5.1. First-order constraint

We first consider the equilibrium solutions imposing the near-axis constraints of both zeroth and first order; namely, (3.2) and (3.8). We present three equilibria representing different classes of optimised (omnigenous) stellarators for the benchmark: a quasi-axisymmetric (QA) and a quasi-helically symmetric (QH) configuration (near-axis constructions of the configurations presented in Landreman & Paul (Reference Landreman and Paul2022), which we refer to as ‘precise QA’ and ‘precise QH’, respectively), and a quasi-isodynamic (QI) configuration (the single field-period vacuum field described in § 8.2 of Plunk et al. Reference Plunk, Landreman and Helander2019). In all cases, we are considering stellarator-symmetric vacuum fields, although note that the constraints derived are equally valid for non-symmetric equilibria.

We start by summarising some of the key verification quantities in table 1 and figure 4.Footnote ⁷ From the comparison of cross-sections, it is clear that the NAE-constrained solutions match the NAE axis and elliptical shapes to leading order, which is not true of the fixed-boundary cases, especially for the QH and QI equilibria. These configurations are more strongly shaped than their QA counterpart, which results in a more limited radius of validity of the NAE. As a result, using the finite radius built surface as a boundary condition leads to an equilibrium that significantly departs from the original near-axis behaviour. Note that the fields are being constructed at a rather low aspect ratio, which enhances this effect. This discrepancy is also apparent in the rotational transform and the behaviour of $|\boldsymbol{B}|$ . The constrained equilibria correctly capture the on-axis rotational transform orders of magnitude better, and in the quasisymmetric fields the symmetry breaking error $f_{B}$ (Garren & Boozer Reference Garren and Boozer1991a ; Rodríguez et al. Reference Rodríguez, Paul and Bhattacharjee2022a ; Dudt et al. Reference Dudt, Conlin, Panici and Kolemen2023) shows the correct scaling $\sim \rho ^2$ . In the case of the QI configuration, the sign of a correct behaviour of $|\boldsymbol{B}|$ near the axis is shown by looking at its value on axis. Note that away from the axis, the global solver completes the equilibrium solution, e.g. in the case of rotational transform, introducing a small global shear.

Although the accuracy achieved in a quantity such as $\iota _0$ is orders of magnitude better than the conventional fixed-boundary case, which would need an extremely high ( $R_0/a\gt 160$ ) aspect ratio to match it (as shown in Landreman et al. (Reference Landreman, Sengupta and Plunk2019), see also Appendix G), $\unicode{x0394} \iota _0$ is not down to machine precision. The imperfect force balance achieved in the equilibrium solution is likely the culprit. As indicated above, unless $\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{\nabla }\psi =0$ exactly, then the rotational transform will depend on second-order shaping, which not being constrained, will affect the precision of the agreement with the NAE. The larger and more complex the higher-order shaping tends to be, the easier it is for this discrepancy to be larger. This explains the relatively worse performance of the QI equilibrium, whose increased shaping can be seen in the heightened difficulty in representing the QI NAE behaviour with the cylindrical toroidal angle, as explained in Appendix F. Despite these limitations, the accuracy is practically sufficient to show that the global equilibrium behaves like the NAE near axis.

Table 1.

Quantitative verification of equilibrium field. Table including the verification measures introduced in § 4.1 comparing the near-axis-constrained (NAE) and the fixed-boundary (Surf) equilibria. The first column shows the aspect ratio of the equilibria considered. As expected, the near-axis-constrained solution performs significantly better than the fixed-boundary one.

Figure 4.

Verification of first-order NAE-constrained equilibria. The figure shows a comparison of DESC first-order NAE-constrained equilibrium solutions (in red) against DESC near-axis fixed-boundary solutions (in green) and the ideal near-axis field (broken line). From left to right the configurations correspond to the precise QA, precise QH and QI introduced in the text. (Top) Cross-sections at $\phi =0$ . (Middle) Rotational transform profile showing the correct matching of the NAE-constrained field. (Bottom) The left plots show the Boozer quasisymmetry metric for the QS solutions (measure of $|\boldsymbol{B}|$ error). The expected quadratic scaling is observed in the NAE-constrained equilibria. The right plot shows the on-axis magnetic field strength for the QI solution (for which $f_B$ is not a meaningful measure).The QA and QI global equilibria were solved in DESC at radial ( $L$ ), poloidal ( $M$ ) and toroidal ( $N$ ) spectral resolutions of $L=9,\,M=9,\,N=24$ , and the near axis using pyQSC and pyQIC. The QH equilibrium was solved at a higher resolution of $L=10,\,M=10,\,N=24$ , which was necessary to achieve good force balance for the fixed-surface solution.

5.2. Second-order constraint

Let us now repeat the construction, but this time using near-axis constructions to second order. We use for that purpose the near-axis versions of the four highly quasisymmetric equilibria in Landreman & Paul (Reference Landreman and Paul2022). Besides precise QA and QH, which we shall now construct to second order, we also consider NAE models of the other two, namely ‘precise QA + well’ and ‘precise QH + well’ configurations (Landreman & Paul Reference Landreman and Paul2022). We now construct NAE-constrained equilibria in DESC, constraining the $\mathcal{O}(\rho ^2)$ $R$ and $Z$ behaviour of the NAE solution. We shall also fix in this case the on-axis $\lambda$ , which was found to help the equilibrium solve while not being overconstraining. We shall be particularly fixated with the DESC solution retaining some of the higher-order features of the equilibrium. Specifically, we place renewed focus on the scaling of the quasisymmetric error and the magnetic well $V''$ (Landreman & Jorge Reference Landreman and Jorge2020). The latter is particularly important, as MHD stability is often found to be at odds with other physics and engineering objectives in stellarator optimisation (Rodríguez Reference Rodríguez2023; Conlin et al. Reference Conlin, Kim, Dudt, Panici and Kolemen2024), but with this constraint one may enforce some measure of stability on the MHS equilibrium solution.

Figure 5.

Verification of second-order NAE-constrained equilibria. The figure shows a comparison of DESC second-order NAE-constrained equilibrium solutions (in red) against DESC fixed-boundary solutions (in green) and the ideal near-axis field (broken line). The plots correspond to the equilibria introduced in the text, the ‘precise QA’ and ‘precise QA + well’ on the left, the ‘precise QH’ and ‘precise QH + well’ on the right. (Top) Cross-sections at $\phi =0$ . (Middle upper) Rotational transform profile showing the correct matching of the NAE-constrained field. (Middle lower) Boozer quasisymmetry metric for the QS solutions (measure of $|\boldsymbol{B}|$ error). The NAE-constrained equilibria generally adhere closely to the expected cubic scaling. (Bottom) Magnetic well parameter, where the NAE-constrained solutions are closer to the NAE value on axis, which in the ‘QA + well’ case corresponds to stability. Each global equilibrium was solved in DESC at radial ( $L$ ), poloidal ( $M$ ) and toroidal ( $N$ ) spectral resolutions of $L=15,\,M=10,\,N=25$ , and the near-axis solutions using pyQSC.

The results for the second order are given in figure 5 and summarised quantitatively in table 2. Once again, the flux surfaces and rotational transform near the axis are better matched by the NAE-constrained equilibrium than by the fixed surface. The Boozer quasisymmetry error ( $f_{B}$ ) (Garren & Boozer Reference Garren and Boozer1991a ; Dudt et al. Reference Dudt, Conlin, Panici and Kolemen2023; Rodríguez et al. Reference Rodríguez, Paul and Bhattacharjee2022a ) scales like $\mathcal{O}(\rho ^3)$ , as expected for a second-order NAE field, highly quasisymmetric at second order (Landreman & Sengupta Reference Landreman and Sengupta2019). Finally, and most importantly, the magnetic well parameter is matched to the NAE solution to a good degree. See Appendix G for another example which obtains stability.

Table 2.

Table with quantitative verification of equilibrium field. Table including the verification measures introduced in § 4.1 comparing the near-axis-constrained (NAE) and the fixed-boundary (Surf) equilibria. The first column shows the aspect ratio of the equilibria considered. As expected, the near-axis-constrained solution performs significantly better than the fixed-boundary one.