2.1. Near-axis geometry

We define a magnetic axis as a $C^\infty$ closed curve $\boldsymbol{r}_0 : \mathbb{T} \to \mathbb{R}^3$ with $\boldsymbol{r}_0' \neq 0$ and a non-zero tangent magnetic field (see § 2.3 for details about the field). We define a near-axis domain about the axis with radius $\sigma$ as

(2.1) \begin{equation} \Omega _\sigma = \left \{\boldsymbol{r} \in \mathbb{R}^3 \mid \forall s \in \mathbb{T}, \ \left \lVert \boldsymbol{r} - \boldsymbol{r}_0(s) \right \rVert \lt \sigma \right \}. \end{equation}

We note that the assumption that the axis is infinitely differentiable is necessary for the near-axis expansion to formally go to arbitrary order, and we will eventually reduce the required regularity for the inputs of the regularised expansion, summarised in box (3.30).

Figure 1. Schematic of the direct near-axis Frenet–Serret coordinate frame.

We define a direct coordinate system $\boldsymbol{r} : \Omega _\sigma ^0 \to \Omega _\sigma$ where $(\rho ,\theta ,s) \in \Omega _\sigma ^0 = [0,\sigma ) \times \mathbb{T}^2$ is the solid torus as (see figure 1)

(2.2) \begin{equation} \boldsymbol{r}(\rho , \theta , s) = \boldsymbol{r}_0(s) + Q(s) \begin{pmatrix} x \\ y \\ 0 \end{pmatrix}, \qquad (x, \ y) = (\rho \cos \theta , \ \rho \sin \theta ), \end{equation}

where $Q$ is an orthonormal basis for the local coordinates at the axis with the tangent vector in the third column, i.e. $\boldsymbol{t} = \boldsymbol{r}_0' / \ell ' = Q\boldsymbol{e}_3$ , where we assume $\ell ' = \left \lVert \boldsymbol r_0' \right \rVert \gt 0$ . Both the $\boldsymbol{x} = (x,y,s)$ and the $\boldsymbol{q} = (\rho , \theta ,s)$ coordinate frames are useful, as $\boldsymbol{x}$ is non-singular with respect to the coordinate transformation, while $\boldsymbol{q}$ diagonalises the near-axis expansion operator. To perform calculus in the $\boldsymbol{q}$ basis, we require the induced metric from $\mathbb{R}^3$ . For this, we first compute the coordinate derivative

(2.3) \begin{equation} F = \frac {\mathrm{d} \boldsymbol{r}}{\mathrm{d} \boldsymbol{q}} = Q \left(\begin{array}{c@{\quad}c@{\quad}c} \cos \theta & - \rho \sin \theta & 0\\ \sin \theta & \rho \cos \theta & 0\\ 0 & 0 & \ell ' \end{array}\right) + Q K \boldsymbol{x} \boldsymbol{e}_3^T, \end{equation}

where $K$ is an antisymmetric matrix determining the derivative of $Q$ in the near-axis basis $Q'(s) = Q K$ . Using this matrix, the metric is computed as

(2.4) \begin{equation} g = F^T F. \end{equation}

For the numerical examples in this paper, we specifically consider the Frenet–Serret frame, meaning the local basis $Q$ and its derivative are defined by

(2.5) \begin{equation} Q = \left(\begin{array}{c@{\quad}c@{\quad}c} | & | & | \\ \boldsymbol{n} & \boldsymbol{b} & \boldsymbol{t} \\ | & | & | \end{array}\right), \qquad K = \ell ' \left(\begin{array}{c@{\quad}c@{\quad}c} 0 & -\tau & \kappa \\ \tau & 0 & 0 \\ -\kappa & 0 & 0 \end{array}\right)\!, \end{equation}

where $\kappa = \left \lVert \boldsymbol{t}' \right \rVert /\ell '$ is the axis curvature, $\boldsymbol{n} = \left \lVert \boldsymbol{t}' \right \rVert /(\kappa \ell ')$ is the normal vector, $\boldsymbol{b} = \boldsymbol{t} \times \boldsymbol{n}$ is the binormal vector and $\tau = -\left \lVert \boldsymbol{b}' \right \rVert /\ell '$ is the axis torsion. Alternative forms of the curvature and torsion are

(2.6) \begin{equation} \kappa = \frac {\left \lVert \boldsymbol{r}_0' \times \boldsymbol{r}_0'' \right \rVert }{(\ell ')^3}, \qquad \tau = \frac {(\boldsymbol{r}_0' \times \boldsymbol{r}_0'') \cdot \boldsymbol{r}_0'''}{\left \lVert \boldsymbol{r}_0' \times \boldsymbol{r}_0'' \right \rVert ^2}. \end{equation}

For the Frenet–Serret coordinate system to be well-defined and non-singular on $\Omega _\sigma ^0$ , we require $\sigma ^{-1} \gt \kappa \gt 0$ . In particular, no straight segments are allowed in the Frenet–Serret frame, disallowing quasi-isodynamic (QI) stellarators (Mata et al. Reference Mata, Plunk and Jorge2022). An alternative choice for axis coordinates that allow for straight segments is Bishop’s coordinates (Bishop Reference Bishop1975; Duignan & Meiss Reference Duignan and Meiss2021).

Replacing the Frenet–Serret basis into (2.3), we obtain

(2.7) \begin{equation} F = Q \left(\begin{array}{c@{\quad}c@{\quad}c} \cos \theta & - \rho \sin \theta & -\ell ' \tau \rho \sin \theta \\ \sin \theta & \rho \cos \theta & \ell ' \tau \rho \cos \theta \\ 0 & 0 & h_s \end{array}\right),\qquad g = \left(\begin{array}{c@{\quad}c@{\quad}c} 1 & 0 & 0 \\ 0 & \rho ^2 & \ell ' \tau \rho ^2 \\ 0 & \ell ' \tau \rho ^2 & (\ell ')^2\tau ^2 \rho ^2 + h_s^2 \end{array}\right)\!, \end{equation}

where

(2.8) \begin{equation} h_s(\rho ,\theta ,s) = \ell ' (1 - \kappa \rho \cos \theta ). \end{equation}

The local volume ratio is given by

(2.9) \begin{equation} \sqrt g = \det F = \rho h_s \end{equation}

and the inverse metric is

(2.10) \begin{equation} g^{-1} = \left(\begin{array}{c@{\quad}c@{\quad}c} 1 & 0 & 0 \\ 0 & (\ell ')^2 \tau ^2 h_s^{-2} + \rho ^{-2} & -\ell ' \tau h_s^{-2} \\ 0 & -\ell ' \tau h_s^{-2} & h_s^{-2} \end{array}\right)\!. \end{equation}

To find metrics associated with the more general coordinate system (2.2), we consider transformations of the form $(\rho ,\theta ,s) \mapsto (\rho ,\omega (\theta ,s),s)$ , where $\omega (\theta ,s) =\theta + \lambda (s)$ . This transformation rotates the orthonormal frame, changing the metric to

(2.11) \begin{equation} g = \left(\begin{array}{c@{\quad}c@{\quad}c} 1 & 0 & 0 \\ 0 & \rho ^2 & \ell ' T \rho ^2 \\ 0 & \ell ' T \rho ^2 & (\ell ')^2T^2 \rho ^2 + h_s^2(\rho ,\omega -\lambda ,s), \end{array}\right)\!, \qquad T = \tau - \frac {\lambda '}{\ell '}. \end{equation}

In this way, the general set of near-axis frames can be represented by a simple replacement of $\tau$ with $T$ . A special case of this transformation is when $T=0$ , yielding (Mercier Reference Mercier1964)

(2.12) \begin{equation} \lambda = \int _{0}^s \tau (s) \ell '(s) \hspace {1pt} \mathrm{d} s. \end{equation}

In this coordinate frame, the metric becomes diagonal: $g = \textrm{Diag}(1, \rho ^2, h_s^2(\rho ,\omega -\lambda ,s))$ . The fact that $g$ is diagonalised is convenient for theoretical manipulations, but variables expressed in terms of $\omega$ are multivalued for curves with non-integer total torsion. This unfortunate consequence is important for numerical methods, as Fourier series cannot be used in $s$ and additional consistency requirements are necessary. This is part of the motivation for using the Frenet–Serret frame for the numerical examples herein.

2.2. Near-axis expansion

Now, we consider expansions of functions $A \in C^\infty (\Omega _\sigma ^0)$ about the magnetic axis. We formally expand in small distances from the axis $\rho \ll \min \kappa ^{-1}$ as

(2.13) \begin{gather} A(\rho , \theta , s) = \sum _{m=0}^{\infty }A_m(\theta ,s) \rho ^m, \qquad A_m(\theta , s) = \sum _{n=0}^m A_{mn}(s) e^{(2n-m){\mathrm{i}}\theta }, \nonumber \\ A_{mn}(s) = \sum _{k\in \mathbb{Z}} A_{mnk}e^{{\mathrm{i}} k s}. \end{gather}

This expansion is not guaranteed to converge anywhere for $A\in C^\infty$ , but it is asymptotic to $A$ near the axis, i.e.

(2.14) \begin{equation} \left | A - A_{\lt m} \right | = \mathcal{O}(\rho ^{m}), \end{equation}

where we define $A_{\lt m}$ as the partial sum

(2.15) \begin{equation} A_{\lt m} = \sum _{n=0}^{m-1} A_n \rho ^n. \end{equation}

In defining magnetic coordinates, we also find it convenient to expand $A$ in $\boldsymbol x$ as

(2.16) \begin{equation} A(x,y,s) = \sum _{\mu =0}^{\infty } \sum _{\nu = 0}^\mu A_{\mu \nu }(s) x^{\mu -\nu } y^\nu , \end{equation}

where we use Latin indices for the $\boldsymbol{q}$ frame and Greek indices for the $\boldsymbol{x}$ frame. If we require $A$ to be real-analytic on $\Omega _\sigma ^0$ , there additionally exists a $\sigma ' \lt \sigma$ so that the the asymptotic series converges uniformly on $\Omega _{\sigma '}^0$ (this does not necessarily extend to all of $\Omega _\sigma ^0$ ).

Throughout this paper, we attempt to minimise the number of complicated summation formulae resulting from the near-axis expansion (NAE). For instance, if we have two functions $A, B\in C^{\infty }$ and we want the $m$ th component of the series of $C=AB$ , we will write $C_{m} = (AB)_{m}$ , rather than

(2.17) \begin{equation} C_{m} = \sum _{\ell = 0}^m A_{m-\ell } B_\ell . \end{equation}

As expressions become increasingly complicated, this notation provides a concise description of the mathematics involved. In Appendix B, we define a number of relevant operations on series that the interested reader can use to expand the expressions within this paper. In § 4, we discuss how this is similarly convenient for the purposes of programming NAE operations. Rather than implementing residuals via complicated summation formulae, the operations in Appendix B are called, allowing for a simple framework for developing new code.

An important exception to the general rule of condensing notation is in defining any linear operators that must be inverted through the course of an asymptotic expansion. Detailed understanding of such operators are necessary for both numerical implementation and analysis on the series.

2.3. Vacuum fields

In the steady state, the vacuum magnetic field $\boldsymbol{B}\in C^\infty (\Omega _\sigma ^0)$ satisfies

(2.18) \begin{equation} \nabla \cdot \boldsymbol{B} = 0, \qquad \boldsymbol{J} = \nabla \times \boldsymbol{B} = 0, \end{equation}

where $\boldsymbol{J}$ is the plasma current density. The fundamental near-axis assumption is that $\boldsymbol{B}$ evaluated on the axis is non-zero and parallel to the magnetic axis, i.e. for some $B_0 \in C^{\infty }(\mathbb{T})$ ,

(2.19) \begin{equation} \boldsymbol{B}(0,\theta ,s) = B_0(s) \boldsymbol{t}(s). \end{equation}

Off the axis, we express the magnetic field on $\Omega _\sigma ^0$ as

(2.20) \begin{equation} \boldsymbol{B} = \nabla \phi + \tilde {\boldsymbol{B}}, \qquad \tilde {\boldsymbol{B}} = \nabla \left [ \int _{0}^s B_0(s') \ell '(s') \hspace {1pt} \mathrm{d} s'\right ] = B_0(s)\ell '(s) \nabla s, \end{equation}

where $\phi \in C^{\infty }(\Omega _\sigma ^0)$ is the magnetic potential satisfying $\left .\nabla \phi \right |_{\rho =0} = 0$ . Then, taking the divergence of (2.20), we find the magnetic scalar potential satisfies Poisson’s equation,

(2.21) \begin{equation} \Delta \phi = - \nabla \cdot \tilde {\boldsymbol{B}}. \end{equation}

By construction, the field $\boldsymbol{B}$ in (2.20) is locally the gradient of some function, so it is curl-free. This means the contributions from $\tilde {\boldsymbol{B}}$ in (2.20) can locally be absorbed to recover Laplace’s equation for the potential. However, because $\Omega _\sigma$ is not simply connected and $\boldsymbol B \cdot \boldsymbol t$ is single-valued on axis, the closed-loop axis integral $\oint \boldsymbol{B} \cdot \mathrm{d} \boldsymbol \ell$ demonstrates that it is not possible for $\boldsymbol{B}$ to be globally the gradient of a single-valued function $\phi$ . In contrast, the Poisson’s equation formulation contains only single-valued functions, which is convenient both numerically and analytically.

In coordinates, we can write the magnetic field in (2.20) as

(2.22) \begin{equation} B^i = g^{ij} \frac {\partial \phi }{\partial q^j} + \tilde {B}^i, \qquad \tilde {B}^i = B_0 \ell ' g^{ij}\frac {\partial s}{\partial q^j}, \end{equation}

where we assume summation over repeated indices and $g_{ij}$ and $g^{ij}$ are the components of the metric and inverse metric, respectively. Then, Poisson’s equation in coordinates becomes

(2.23) \begin{equation} \sqrt {g}^{-1}\frac {\partial }{\partial q^i}\left ( \sqrt {g} g^{ij} \frac {\partial \phi }{\partial q^j} \right ) = - \sqrt {g}^{-1} \frac {\partial }{\partial q^i}\left (\sqrt {g} \tilde {B}^i \right ). \end{equation}

Multiplying this by $h_s$ and expanding this in the $\boldsymbol q$ coordinate system, we have

(2.24) \begin{align} &\frac {1}{\rho ^2}\left [ \rho \frac {\partial }{\partial \rho }\left (\mathsf{A} \rho \frac {\partial \phi }{\partial \rho } \right )\right . + \left .\frac {\partial }{\partial \theta }\left (\mathsf{B} \frac {\partial \phi }{\partial \theta }\right )\right ] = \nonumber \\ &- \frac {\partial }{\partial \theta }\left ( \mathsf{C} \left ( \frac {\partial \phi }{\partial s} + B_0 \ell '\right ) \right ) - \frac {\partial }{\partial s}\left ( \mathsf{C} \frac {\partial \phi }{\partial \theta } \right ) - \frac {\partial }{\partial s} \left (\mathsf{D} \left (\frac {\partial \phi }{\partial s} + B_0 \ell ' \right ) \right ), \end{align}

where

(2.25) \begin{align} \mathsf{A} &= h_s, & \mathsf{B} &= h_s + \rho ^2 h_s^{-1} (\ell ')^2 \tau ^2 , & \mathsf{C} &= -h_s^{-1} \ell ' \tau , & \mathsf{D} &= h_s^{-1}. \end{align}

From here, we can substitute the asymptotic expansions of $\phi$ and the coefficients $\mathsf{A}$ , $\mathsf{B}$ , $\mathsf{C}$ and $\mathsf{D}$ into (2.24). At each order in $\rho$ , Poisson’s equation becomes

(2.26) \begin{align} \ell ' \left (m^2 + \frac {\partial ^2}{\partial \theta ^2}\right ) \phi _{m} =& - (\nabla \cdot \tilde {B})_{m-2} -(\Delta \phi _{\lt m})_{m-2}\nonumber \\ =&- m(m-1) \mathsf{A}_1 \phi _{m-1} - \left [\frac {\partial }{\partial \theta }\left (\mathsf{B} \frac {\partial \phi _{\lt m}}{\partial \theta } \right )\right ]_m \nonumber \\ &- \left [ \frac {\partial }{\partial \theta }\left ( \mathsf{C} \left (\frac {\partial \phi _{\lt m}}{\partial s} + B_0 \ell ' \right ) \right ) + \frac {\partial }{\partial s}\left ( \mathsf{C} \frac {\partial \phi _{\lt m}}{\partial \theta } \right )\right ]_{m-2} \\ \nonumber &- \left [\frac {\partial }{\partial s} \left (\mathsf{D} \left ( \frac {\partial \phi _{\lt m}}{\partial s} + B_0 \ell '\right ) \right )\right ]_{m-2}, \end{align}

where $\mathsf{A}_1 = (h_s)_1 = -\ell ' \kappa \cos \theta$ . The right-hand side of (2.26) does not depend on orders of $\phi$ higher than $\phi _{m-1}$ , so inverting the left-hand-side operator gives an iteration for obtaining $\phi _m$ at each order.

However, the operator $\ell ' (m^2 + ({\partial ^2}/{\partial \theta ^2)})$ is singular, so we must confirm there are no secular terms. Specifically, we always have an unknown homogeneous solution at $\mathcal{O}(m)$ of the form $\phi _{m0} e^{-i m \theta } + \phi _{mm} e^{i m \theta }$ . For this, we use Fredholm’s alternative, which states that (2.26) is solvable if the right-hand side is orthogonal to the null space of the adjoint of the operator on the left-hand side. Because the operator is self-adjoint, the right-hand side must be orthogonal to $e^{\pm i m \theta }$ , i.e.

(2.27) \begin{align} &\left \langle m(m-1) \mathsf{A}_1 \phi _{m-1} + \left [\frac {\partial }{\partial \theta }\left (\mathsf{B}_{\gt 0} * \frac {\partial \phi _{\lt m}}{\partial \theta } \right )\right ]_m \right . + \left [ \frac {\partial }{\partial \theta }\left ( \mathsf{C} \left ( \frac {\partial \phi _{\lt m}}{\partial s} + B_0 \ell '\right ) \right )\right ]_{m-2}\nonumber \\ &\qquad\qquad\qquad\quad\quad\,\,\,+ \left .\left [\frac {\partial }{\partial s}\left ( \mathsf{C} \frac {\partial \phi }{\partial \theta } \right ) + \frac {\partial }{\partial s} \left (\mathsf{D} \left ( \frac {\partial \phi _{\lt m}}{\partial s} + B_0 \ell '\right ) \right )\right ]_{m-2}, \ e^{\pm i m \theta }\right \rangle =0, \end{align}

where the inner product is defined by

(2.28) \begin{equation} \left \langle f , \ g \right \rangle = \int _{0}^{2\pi } f(\theta )\overline {g}(\theta ) \hspace {1pt} \mathrm{d} \theta . \end{equation}

The $m-2$ coefficient of any analytic function is orthogonal to $e^{i m \theta }$ , so we can remove the $\mathsf{C}$ and $\mathsf{D}$ terms. The same argument allows us to remove the torsion terms in $\mathsf{B}$ . This means only contributions from $\mathsf{A}_1 = \mathsf{B}_1 = - \ell ' \kappa \cos \theta$ and $\phi _{m-1}$ survive, so we only need to verify

(2.29) \begin{align} \left \langle m(m-1) \phi _{m-1} \cos \theta + \left [\frac {\partial }{\partial \theta }\left (\frac {\partial \phi _{m-1}}{\partial \theta } \cos \theta \right )\right ]_m , \ e^{\pm i m \theta } \right \rangle = 0. \end{align}

Using the identity $\cos \theta = (e^{i \theta } + e^{-i \theta })/2$ , a quick calculation confirms the above identity holds.

Now, we consider the problem where $\phi$ is unknown. The fact that the Fredholm condition is automatically satisfied at each order implies that $\phi _{m0}$ and $\phi _{mm}$ are free parameters at each order in the near-axis expansion. So, these coefficients are an infinite-dimensional set of initial conditions for the near-axis expansion of Poisson’s equation. Intuitively, one can think of the coefficients $\phi _{m0}$ and $\phi _{mm}$ as specifying the Fourier coefficients of an infinitely thin tube about the magnetic axis. In this way, the imposition of conditions at each order compensates for the fact that partial differential equations (PDEs) typically satisfy conditions on co-dimension-1 surfaces, whereas the near-axis expansion is specified on a co-dimension-2 curve.

In addition to $\phi _{m0}$ and $\phi _{mm}$ as free parameters, we also treat $\boldsymbol{r}_0$ and $B_0$ as inputs to the near-axis problem. The requirement that the magnetic field be tangent to the axis with magnitude $B_0$ results in a constraint that $\phi _{00} = \phi _{10} = \phi _{11} = 0$ . In total, the direct vacuum near-axis problem can be written as

(2.30)

2.4. Straight field-line coordinates: leading order

Given a solution magnetic field from box (2.30), we consider the problem of finding straight field-line magnetic coordinates. We assume that the magnetic field is locally elliptic about the axis and the rotation number is irrational, so that the leading-order behaviour is rotation about the magnetic axis. This means that both hyperbolic orbits (x-points) and on-axis resonant perturbations are excluded from this work. In the language of Hamiltonian normal forms, the leading order field-line dynamics is conjugate to a non-resonant harmonic oscillator; see Burby et al. Reference Burby, Duignan and Meiss(2021) and Duignan & Meiss (Reference Duignan and Meiss2021) for a more rigorous derivation of magnetic coordinates in the near-axis expansion. We note that our process of finding coordinates is formal: we make no claims that this problem converges in the limit. However, in § 5, we find that this procedure appears to converge well numerically.

To find magnetic coordinates, we attempt to build a conjugacy between magnetic field-line dynamics $\dot {\boldsymbol{r}} = (B^s)^{-1} \boldsymbol{B}(\boldsymbol{r})$ and straight field-line dynamics $\dot {\boldsymbol{\xi }} = (- \iota (\psi ) \eta , \, \iota (\psi ) \xi , \, 1)$ , where $\boldsymbol{\xi } = (\xi , \, \eta , \, s) \in \mathbb{R}^2 \times \mathbb{T}$ are Cartesian coordinates, $\psi =\xi ^2 + \eta ^2$ is a flux-like coordinate and $\iota$ is the rotational transform. To make the connection with straight field-line coordinates precise, consider the transformation to polar coordinates $(\xi , \eta ) \to \sqrt {\psi }(\cos \gamma , \sin \gamma )$ . Then, the field-line is traced by $(\dot {\psi },\, \dot {\gamma },\, \dot {s}) = (0, \, \iota (\psi ), \, 1)$ , i.e. magnetic field lines are straight with slope $\iota$ . However, we use the Cartesian version of magnetic coordinates because it removes the coordinate singularity associated with polar coordinates, simplifying the following steps.

Figure 2. A schematic of the process of finding straight field-line coordinates. On the left, we plot the surfaces of the magnetic field $(h^x,h^y)$ on a cross-section for fixed $s$ . Moving one plot to the right, the leading correction transforms to a coordinate frame where the main elliptic component is eliminated. Going one further, the next correction accounts for the most prominent triangularity. This process continues until, in $(\xi ,\eta )$ coordinates, the magnetic surfaces are nested circles.

There are two main steps to our process of finding magnetic coordinates: the leading-order problem and the higher-order problems (see figure 2 for a sketch of the process). If we use the notation $\tilde {\boldsymbol{x}} = (x,\, y)$ and $\tilde {\boldsymbol{\xi }} = (\xi , \, \eta )$ for the out-of-plane coordinates, the leading order transformation takes the form

(2.31) \begin{equation} \tilde {\boldsymbol{\xi }}_1 = \frac {\partial \tilde {\boldsymbol{\xi }}_1}{\partial \tilde {\boldsymbol{x}}}\tilde {\boldsymbol{x}} = G_0 \tilde {\boldsymbol{x}} = \begin{pmatrix} \xi _{10} x + \xi _{11} y \\ \eta _{10} x + \eta _{11} y \end{pmatrix}. \end{equation}

We will find that the problem for $G_0$ is an eigenvalue problem for the on-axis rotation number, as is typical for the linearised dynamics about a fixed point. In the following section, we will discuss the inductive step to higher orders.

To begin, consider the contravariant form of the Cartesian near-axis magnetic field

(2.32) \begin{equation} \frac {1}{B^s} \boldsymbol B = \frac {1}{B^s} \frac {\mathrm{d} \boldsymbol{r}}{\mathrm{d} \boldsymbol{x}} \begin{pmatrix} B^x \\ B^y \\ B^s \end{pmatrix} = \frac {B^x}{B^s} \frac {\partial \boldsymbol{r}}{\partial x} + \frac {B^y}{B^s} \frac {\partial \boldsymbol{r}}{\partial y} + \frac {\partial \boldsymbol{r}}{\partial s}. \end{equation}

We would like to equate this to the straight field-line dynamics as

(2.33) \begin{equation} \frac {1}{B^s} \boldsymbol{B} = \frac {\mathrm{d} \boldsymbol{r}}{\mathrm{d} \boldsymbol{\xi }} \begin{pmatrix} - \iota (\psi ) \eta \\ \iota (\psi ) \xi \\ 1 \end{pmatrix}, \end{equation}

where $\iota$ depends smoothly upon the radial label as

(2.34) \begin{equation} \iota = \iota _0 + \iota _2 \psi + \iota _4\psi ^2 + \ldots , \end{equation}

where we emphasise $\iota _\mu = 0$ for odd $\mu$ . Multiplying both sides by $\mathrm{d} \boldsymbol{\xi } / \mathrm{d} \boldsymbol{r}$ , we find the Floquet conjugacy problem

(2.35) \begin{equation} \frac {\partial \tilde {\boldsymbol{\xi }}}{\partial s} = -G(\boldsymbol{\xi }) \tilde {\boldsymbol{h}}^{\boldsymbol{x}} +\iota (\psi (\boldsymbol{\xi })) J \tilde {\boldsymbol{\xi }}, \end{equation}

where

(2.36) \begin{equation} G = \frac {\partial \tilde {\boldsymbol{\xi }}}{\partial \tilde {\boldsymbol{x}}}, \qquad \tilde {\boldsymbol{h}}^{\boldsymbol{x}} = \begin{pmatrix} h^x \\ h^y \end{pmatrix} = \frac {1}{B^s} \begin{pmatrix} B^x \\ B^y \end{pmatrix}, \qquad J = \begin{pmatrix} 0 & -1 \\ 1 & 0 \end{pmatrix}, \end{equation}

where $J$ is known as the symplectic matrix. Our problem is to solve (2.35) for $\xi (x,y)$ , $\eta (x,y)$ and $\iota (\psi )$ .

To find the leading-order problem for (2.35), we note the magnetic field $(h^x,h^y)$ is linear at leading order:

(2.37) \begin{equation} \tilde {\boldsymbol{h}}^{\boldsymbol{x}} = \begin{pmatrix} h^x_{10} x + h^x_{11} y \\ h^y_{10} x + h^y_{11} y \end{pmatrix} + \mathcal{O}(\psi ) = H_0(s) \tilde {\boldsymbol{x}} + \mathcal{O}(\psi ). \end{equation}

Substituting this, (2.31) and $\iota = \iota _0 + \mathcal{O}(\psi )$ into (2.35), we have

(2.38) \begin{equation} \frac {\partial G_0}{\partial s} + G_0 H_0 = \iota _0 J G_0. \end{equation}

The leading order problem (2.38) is a Floquet eigenvalue problem for the linearised field-line dynamics about the magnetic axis. Assuming that the near-axis expansion is elliptic at leading order, the value of $\iota _0$ is real. Otherwise, $\iota _0$ is not real, meaning ellipticity can be numerically verified for a given input.

There are many equivalent solutions to (2.38), owing to the symmetries that if $(G_0, \iota _0)$ satisfies (2.38), then the following are also solutions:

(2.39) \begin{align} (R(ns) G_0, n + \iota _0), \qquad (J G_0, \iota _0), \qquad \left (\begin{array}{l@{\quad}c} 0 & 1 \\ 1 & 0 \end{array}\right) G_0, -\iota _0, \end{align}

where $R$ is a rotation matrix

(2.40) \begin{equation} R(\theta ) = \left(\begin{array}{l@{\quad}c} \cos \theta & - \sin \theta \\ \sin \theta & \cos \theta \end{array}\right)\!. \end{equation}

The question of which solution to choose is then a question of practicalities. For instance, one could choose the rotational transform corresponding to a non-twisting right-handed coordinate frame $\boldsymbol{\xi }_{1}$ . Here, ‘non-twisting’ means that closed coordinate lines near the axis, implicitly defined as curves $\boldsymbol r(x,y,s)$ where $\tilde {\boldsymbol{\xi }}(x,y,s) \neq 0$ is held constant, can be continuously deformed to a point on $\mathbb{R}^3 \backslash \boldsymbol r_0$ . In other words, the coordinates lines do not link with the axis. In this frame, $\iota _0$ agrees with the intuitive definition of the rotational transform as the limiting ratio of poloidal turns divided by toroidal turns of fieldlines about the axis. Other choices may have other benefits, e.g. there may be an eigenfunction that $G_0$ behaves best numerically. In this paper, we opt for an option that is easy to implement: we take the real solution where $\iota _0$ has the smallest magnitude and $\det G_0$ is positive. From here, other equivalent coordinates can easily be found by applying the transformations in (2.39).

For the scaling of $G_0$ , we choose $\psi$ to be the actual magnetic flux at leading order. The formula for the flux is

(2.41) \begin{equation} \psi = \int _{\boldsymbol{r}(D_\psi , 0)} \boldsymbol{B} \cdot \boldsymbol{t} \hspace {1pt} \mathrm{d} A = \int _{\boldsymbol{r}(D_\psi , 0)} B_s \hspace {1pt} \mathrm{d} A(\boldsymbol{r}), \end{equation}

where $D_\psi = \{(\xi ,\eta ) \mid \xi ^2 + \eta ^2 \lt \psi \}$ and

(2.42) \begin{equation} B_s = g_{sj} B^j = \frac {\partial \phi }{\partial s} + B_0 \ell '. \end{equation}

Pulling this back to the $\boldsymbol{\xi }$ frame, we have

(2.43) \begin{equation} \psi = \int _{D_\psi }B_s(\boldsymbol{x}(\boldsymbol{\xi }, s), 0) (\det G)^{-1} \hspace {1pt} \mathrm{d} A(\boldsymbol{\xi }). \end{equation}

Both $B_s$ and $G_0$ are constant in $\boldsymbol{\xi }$ to leading order, so (2.43) at leading order becomes

(2.44) \begin{equation} \int _{D_\psi }B_s(\boldsymbol{x}(\boldsymbol{\xi }, s), 0) (\det G)^{-1} \hspace {1pt} \mathrm{d} A(\boldsymbol{\xi }) = \frac {(B_{s})_0}{\det G_0} \pi \psi + \mathcal{O}(\psi ^2). \end{equation}

Setting this equal to $\psi$ , we find

(2.45) \begin{equation} \det G_0 = \pi (B_s)_0, \end{equation}

where we note that this is only possible when $(B_s)_0$ is chosen to be positive.

2.5. Straight field-line coordinates: higher order

Now that we have the leading-order behaviour, we iterate to go to higher order. To do so, first define the near-axis expansion of $\tilde {\boldsymbol{\xi }}$ near the axis as

(2.46) \begin{equation} \tilde {\boldsymbol{\xi }} = \sum _{\mu =1}^\infty \tilde {\boldsymbol{\xi }}_\mu (\tilde {\boldsymbol{x}}, s), \qquad \tilde {\boldsymbol{\xi }}_\mu = \sum _{\nu = 0}^\mu \tilde {\boldsymbol{\xi }}_{\mu \nu }(s) x^{\mu -\nu } y^\nu , \end{equation}

and define the partial sums as

(2.47) \begin{equation} \tilde {\boldsymbol{\xi }}_{\leqslant \mu } = \sum _{\mu ' = 1}^{\mu }\tilde {\boldsymbol{\xi }}_{\mu '}, \end{equation}

where at leading order, $\tilde {\boldsymbol{\xi }}_{\leqslant 1} = \tilde {\boldsymbol{\xi }}_1 = G_0 \tilde {\boldsymbol{x}}.$

At each order in the iteration, we consider the update to be a function of the previous coordinates, i.e.

(2.48) \begin{align} \tilde {\boldsymbol{\xi }}_{\leqslant \mu +1}\Big(\xi _{\leqslant \mu }, \eta _{\leqslant \mu }, s_{\leqslant \mu }\Big) &= \tilde {\boldsymbol{\xi }}_{\leqslant \mu } + \sum _{\nu =0}^\mu \tilde {\boldsymbol{\xi }}_{\mu \nu }\Big(s_{\leqslant \mu }\Big) \xi _{\leqslant \mu }^{\mu -\nu } \eta _{\leqslant \mu }^\nu \nonumber,\\ &= \tilde {\boldsymbol{\xi }}_{\leqslant \mu } + \sum _{\nu =0}^\mu \tilde {\boldsymbol{\xi }}_{\mu \nu } \xi _{1}^{\mu -\nu } \eta _{1}^\nu + \mathcal{O}\Big(\psi ^{(\mu +2)/2}\Big). \end{align}

We explicitly write the transformed toroidal coordinate $s_{\leqslant \mu }=s$ so that it is clear that $({\partial }/{\partial s})$ and $({\partial }/{\partial s_{\leqslant \mu })}$ are different operators. The purpose of performing the update in this way is primarily to make the update step as clear as possible.

To wit, the magnetic field in the new frame satisfies

(2.49) \begin{equation} \frac {1}{B^s} \boldsymbol{B} = \frac {\mathrm{d} \boldsymbol{r}}{\mathrm{d} \boldsymbol{\xi }_{\leqslant \mu }} \begin{pmatrix} h^{\xi _{\leqslant \mu }} \\ h^{\eta _{\leqslant \mu }} \\ 1 \end{pmatrix} = \frac {\mathrm{d} \boldsymbol{r}}{\mathrm{d} \boldsymbol{\xi }} \begin{pmatrix} -\iota \eta \\ \iota \xi \\ 1 \end{pmatrix}\!. \end{equation}

We can use the first equality to write

(2.50) \begin{align} \tilde {\boldsymbol{h}}^{\boldsymbol{\xi }_{\leqslant \mu }}(\boldsymbol{\xi }_{\leqslant \mu }) &= \begin{pmatrix} h^{\xi _{\leqslant \mu }} \big(\boldsymbol{\xi }_{\leqslant \mu }\big) \\[2pt] h^{\eta _{\leqslant \mu }} \big(\boldsymbol{\xi }_{\leqslant \mu }\big) \end{pmatrix} \nonumber \\ &= \frac {\partial \tilde {\boldsymbol{\xi }}_{\leqslant \mu }}{\partial \tilde {\boldsymbol{x}}} \tilde {\boldsymbol{h}}^{\boldsymbol{x}}\left(\boldsymbol{x}\big(\boldsymbol{\xi }_{\leqslant \mu }\big)\right) + \frac {\partial \tilde {\boldsymbol{\xi }}_{\leqslant \mu }}{\partial s}, \end{align}

(2.51) \begin{align} &= [\iota (\psi ) J \tilde {\boldsymbol{\xi }}]_{\leqslant \mu } + \left [ \frac {\partial \tilde {\boldsymbol{\xi }}_{\leqslant \mu }}{\partial \tilde {\boldsymbol{x}}}\tilde {\boldsymbol{h}}^{\boldsymbol{x}}(\boldsymbol{x}(\boldsymbol{\xi }_{\leqslant \mu }))\right ]_{\gt \mu }\!,\\[6pt] \nonumber\end{align}

where we use the notation $f_{\gt \mu } = \sum _{\nu \gt \mu } f_{\nu }$ . To get from the second to the third line, we have used the inductive assumption that $\boldsymbol{\xi }_{\leqslant \mu }$ matches $\boldsymbol{\xi }$ up to order $\mu$ , meaning that $\boldsymbol{\xi }_{\leqslant \mu }$ is a straight field-line coordinate system up to order $\mu$ . In this way, we will find that the update residual depends neatly upon the transformed magnetic field.

It is worth noting that (2.51) has two new operations that have not been introduced so far. The first is that we are computing $\boldsymbol{x}(\boldsymbol{\xi }_{\leqslant \mu })$ , i.e. we are inverting the coordinate transformation. The second is that we are composing functions with this inversion as $\tilde {\boldsymbol{h}}^{\boldsymbol{x}}(\boldsymbol{x}(\boldsymbol{\xi }_{\leqslant \mu }))$ . So, for any function $f(\boldsymbol{x})$ and coordinate transformation $\boldsymbol{\xi }(\boldsymbol{x})$ , we can compute the equivalent function in indirect coordinates as $f(\boldsymbol{\xi })$ using these two steps. Moreover, the inverse transformation $\boldsymbol{x}(\boldsymbol{\xi })$ can be precomputed for all transformations one wishes to perform of this type. This gives a framework to move back and forth between direct and indirect near-axis formalisms to high order. For more details on how these transformations are computed, see Appendices (B.6) and (B.7).

Given the residual field (2.51), we can use the second equality in (2.49) to find the updated equation

(2.52) \begin{equation} \frac {\partial \tilde {\boldsymbol{\xi }}}{\partial s_{\leqslant \mu }} = - \frac {\partial \tilde {\boldsymbol{\xi }}}{\partial \tilde {\boldsymbol{\xi }}_{\leqslant \mu }} \tilde {\boldsymbol{h}}^{\boldsymbol{\xi }_{\leqslant \mu }}\big(\boldsymbol{\xi }_{\leqslant \mu }\big) + \iota (\psi ) J \tilde {\boldsymbol{\xi }}. \end{equation}

Note that this has the exact same form as (2.35), except we have shifted the underlying coordinates. Substituting $\tilde {\boldsymbol{\xi }}_{\leqslant \mu +1} = \tilde {\boldsymbol{\xi }}_{\leqslant \mu } + \tilde {\boldsymbol{\xi }}_\mu$ into (2.52), we obtain

(2.53) \begin{align} \frac {\partial \tilde {\boldsymbol{\xi }}_{\mu +1}}{\partial s_{\leqslant \mu }} + \iota _0 \frac {\partial \tilde {\boldsymbol{\xi }}_{\mu +1}}{\partial \tilde {\boldsymbol{\xi }}_{\leqslant \mu }} J \tilde {\boldsymbol{\xi }}_{\leqslant \mu } - \iota _0 J \tilde {\boldsymbol{\xi }}_{\mu +1} &= - \tilde {\boldsymbol{h}}^{\boldsymbol{\xi }_{\leqslant \mu }}_{\mu +1} + \iota _\mu \psi ^{\mu /2}J \tilde {\boldsymbol{\xi }}_{\leqslant \mu } \nonumber \\ &= \boldsymbol{F}_{\mu +1}\!. \end{align}

Because the leading order problem (2.38) is an eigenvalue problem, (2.53) has the form of a higher-order correction to the eigenvalue and eigenfunction. To see what we might expect, consider the eigenvalue problem $K\boldsymbol y=\lambda M \boldsymbol y$ , where each term is expanded in a small parameter, e.g. $K=K_0 + K_1 \epsilon + K_2 \epsilon ^2 + \ldots$ for small $\epsilon$ . The analogous update equation would be

(2.54) \begin{equation} (K_0 - \lambda _0 B_0)\boldsymbol{y}_{\mu } = \lambda _\mu M_0 \boldsymbol{y}_0 + \boldsymbol{R}_\mu , \end{equation}

where $\boldsymbol{R}_\mu$ contains all of the residual terms. Assume for simplicity that $\lambda _0$ is an isolated eigenvalue. Then, there is a single secular term, which can be identified by taking the inner product of the above expression with the leading left eigenvector $\boldsymbol{z}_0$ to give $\lambda _\mu = -(\boldsymbol{z}_0^T \boldsymbol{R}_\mu )/(\boldsymbol{z}_0^T M_0 \boldsymbol{y}_0)$ . Once this is satisfied, the equation can be solved for $\boldsymbol{y}_\mu$ , where we typically choose the free component in $\boldsymbol{y}_0$ so that the norm is constant.

To perform the same steps on (2.53), we first diagonalise the left-hand-side operator by converting to polar coordinates $(\xi _{\leqslant \mu }, \eta _{\leqslant \mu }) = R_{\leqslant \mu } (\cos \Theta _{\leqslant \mu }, \sin \Theta _{\leqslant \mu })$ as

(2.55) \begin{equation} \frac {\partial \tilde {\boldsymbol{\xi }}_{\mu +1}}{\partial s_{\leqslant \mu }} + \iota _0 \frac {\partial \tilde {\boldsymbol{\xi }}_{\mu +1}}{\partial \Theta _{\leqslant \mu }} - \iota _0 J \tilde {\boldsymbol{\xi }}_{\mu +1} = \boldsymbol{F}_{\mu +1}, \end{equation}

where

(2.56) \begin{equation} \tilde {\boldsymbol{\xi }}_{\mu +1} = (R_{\leqslant \mu })^{\mu +1} \sum _{n=0}^{\mu } \sum _{\ell \in \mathbb{Z}} \tilde {\boldsymbol{\xi }}_{\mu +1,n\ell } e^{(2n-\mu ){\mathrm{i}} \theta + {\mathrm{i}} \ell s_{\leqslant \mu }}. \end{equation}

After substitution, we find that

(2.57) \begin{equation} \left ((\iota _0 (2n - \mu - 1) + \ell ) {\mathrm{i}} I - \iota _0 J \right )\tilde {\boldsymbol{\xi }}_{\mu +1,n\ell } = \boldsymbol{F}_{\mu +1,n\ell }, \end{equation}

where the $2\times 2$ left-hand-side matrix has the eigenvalues $\lambda _{\mu +1,n\ell 0} = \iota _0 {\mathrm{i}} (2n - \mu ) + {\mathrm{i}} \ell$ and $\lambda _{\mu +1,n \ell 1} = \iota _0 {\mathrm{i}} (2n - \mu - 2) + {\mathrm{i}} \ell$ with corresponding right eigenvectors $(\boldsymbol{v}_0, \boldsymbol{v}_1) = ([-{\mathrm{i}}/\sqrt {2}, 1/\sqrt {2}], [{\mathrm{i}}/\sqrt {2},1/ \sqrt {2}])$ . So, the resulting updates in the coefficients are

(2.58) \begin{equation} \tilde {\boldsymbol{\xi }}_{\mu +1,n\ell } = \sum _{k\in \{0,1\}} \frac {1}{{\mathrm{i}}(\iota _0 (2n - \mu - 2k) + \ell )} \left \langle \overline {\boldsymbol{v}}_{k} , \ \boldsymbol{F}_{\mu +1,n\ell } \right \rangle \boldsymbol{v}_{k}, \end{equation}

where $\overline {\boldsymbol{v}}_k$ are the corresponding left eigenvectors due to $(\iota _0(2n-\mu -1)+\ell ){\mathrm{i}} I - \iota _0 J$ being skew-adjoint.

There are two cases where the update (2.58) fails. The first case is when $2n - \mu - 2k = 0$ and $\ell = 0$ , occurring only when $\mu$ is even. This is the standard secularity that indicates that $\iota _\mu$ must be updated, giving the condition that $\left \langle \overline {\boldsymbol{v}}_{0} , \ \boldsymbol{F}_{\mu +1,\mu /2,0} \right \rangle = \left \langle \overline {\boldsymbol{v}}_{1} , \ \boldsymbol{F}_{\mu +1,\mu /2+1,0} \right \rangle = 0$ for single-valued solutions, where we note that these formulae are equivalent for real magnetic fields. The resulting formula is

(2.59) \begin{equation} \iota _\mu = h^{\eta _{\leqslant \mu }}_{\mu +1,\mu /2,0} - {\mathrm{i}} h^{\xi _{\leqslant \mu }}_{\mu +1,\mu /2,0}. \end{equation}

The second case of failure is when $\iota _0$ is rational, as then there are other values of $(\mu ,n,k,\ell )$ such that (2.58) is singular. This is attributed to the expanding number of $\theta$ modes at each order, where higher-order poloidal perturbations resonate with the axis. To avoid this, extra resonant terms in the higher-order magnetic field must be introduced to avoid secularity (Duignan & Meiss Reference Duignan and Meiss2021) (equivalently, this requires adding terms to the Hamiltonian normal form). Here, we assume that $\iota _0$ is irrational so that the iteration is well defined.

We note that there is still one undetermined part of the problem: what value to choose for $\langle {\overline {\boldsymbol{v}}_{0}},{\tilde {\boldsymbol{\xi }}_{\mu +1,\mu /2,0}}\rangle$ and its complex conjugate. Because this is arbitrary, we currently set this coefficient to $0$ . However, other options could be to choose this value to improve the radius of convergence or to match the flux (2.43). In summary, the straight field-line coordinate process is

(2.60)

3.1. Ill-posedness

We define a problem as ill-posed if it is not well-posed, where the standard definition of a well-posed problem is that:

1. (i) the solution exists;

2. (ii) the solution is unique;

3. (iii) the solution is continuous in the initial data.

We note that the interpretations of these statements depend on what space we require the solution to belong to and over which topology continuity is described in. For instance, it is straightforward to show existence and uniqueness in the sense of a formal power series.

Note that this proposition says nothing about the convergence of the power series to a solution off-axis; it only shows that we can find the coefficients of the power series. So, formal existence does not necessarily imply good or consistent computational results.

Another straightforward existence result for harmonic inputs is the following.

Proposition (3.2) is a useful result because it says that vacuum fields can, in principle, be written as solutions to the infinite near-axis problem. However, this is a difficult theorem to use in practice, as it is difficult to verify a priori whether the input data to the near-axis expansion agree with a solution of Poisson’s equation.

So, for a more computational approach, we must define normed spaces of inputs and outputs that agree with notions of convergence. To intuit what the correct space may be, we observe that the radial direction behaves as a ‘time-like’ variable, whereas the $\theta$ and $s$ behave more like spatial variables. That is, the near-axis PDE can be thought of as propagating surface information off the axis. This motivates a decision to separate our treatment of these coordinates. Moreover, we desire convergence in a power series in the radial variable, so we choose to treat it in an analytic manner. In contrast, the coefficients $\phi _m(\theta ,s)$ are obtained by solving a linear PDE at each order, so we treat them in a Sobolev sense as a function of $\theta$ and $s$ .

To this end, let $\alpha = (\alpha _1,\alpha _2)$ be a multi-index of degree $2$ and $\left | \alpha \right | = \sum _j \alpha _j$ . We define the $H^q$ Sobolev norm of functions $f : \mathbb{T}^2 \to \mathbb{R}$ as

(3.1) \begin{equation} \left \lVert f \right \rVert ^2_{H^q} = \sum _{\left | \alpha \right | \leqslant q} \left \lVert D_\alpha f \right \rVert ^2_{L^2}, \end{equation}

where

(3.2) \begin{equation} D_\alpha f = \frac {\partial ^{\left | \alpha \right |}f}{\partial \theta ^{\alpha _1} \partial s^{\alpha _2}}, \qquad \left \lVert f \right \rVert ^2_{L^2} = \int _{\mathbb{T}^2}\left | f \right |^2 \hspace {1pt} \mathrm{d} \mu , \end{equation}

and $\mu$ is the Lebesgue measure on $\mathbb{T}^2$ . We additionally define the $C^q$ norm of a $q$ -times differentiable function $f : \mathbb{T}^2 \to \mathbb{R}$ as

(3.3) \begin{equation} \left \lVert f \right \rVert _{C^q} = \sum _{\left | \alpha \right |\leqslant q} \sup _{(\theta ,s)\in \mathbb{T}^2} \left | D_\alpha f(\theta ,s) \right |. \end{equation}

Then, we define the following convenient near-axis function space.

Definition 3.3. Let $\sigma \gt 0$ and $W$ be a Banach space on functions on $\mathbb{T}^2$ . We define a function $f : \Omega _\sigma ^0 \to \mathbb{R}^d$ as $(\sigma ,W)$ -analytic if $f$ has a convergent near-axis expansion of the form

(3.4) \begin{equation} f(\rho ,\theta ,s) = \sum _{m=0}^\infty f_m(\theta , s)\rho ^m, \qquad f_m(\theta ,s) = \sum _{n=0}^m f_{mn}(s) e^{(2n-m){\mathrm{i}} \theta }, \end{equation}

where the norm

(3.5) \begin{equation} \left \lVert f \right \rVert _{\sigma ,W} = \sup _n \left (\sigma ^n \left \lVert f_n \right \rVert _{W}\right ) \end{equation}

is bounded. Here, convergence of the near-axis expansion is pointwise in $\rho$ and in norm in $(\theta ,s)$ , i.e. for all $\rho \lt \sigma$ ,

(3.6) \begin{equation} \lim _{M\to \infty }\left \lVert f(\rho ,\cdot ,\cdot ) - \sum _{m=0}^M f_m(\cdot ,\cdot ) \rho ^m \right \rVert _{W} \to 0. \end{equation}

Paralleling standard linear regularity theory (Evans Reference Evans2010), we will consider near-axis solutions $\phi$ to belong to a Sobolev $(\sigma ,H^q)$ -analytic space, while near-axis PDE coefficients belong to differentiable $(\sigma ,C^q)$ -analytic spaces. Because the coefficients of Poisson’s equation (2.21) are functions of the metric, control on the $(\sigma ,C^q)$ -analytic norm of the coefficients in the Frenet–Serret coordinate system will be entirely determined by the differentiability class of the axis $\boldsymbol r_0 \in C^q(\mathbb{T})$ .

To build some intuition for $(\sigma ,W)$ -analytic functions, we turn to some straightforward facts about their convergence. Let $f$ be $(\sigma ,W)$ -analytic. Then, for any $F\geqslant \left \lVert f \right \rVert _{\sigma ,W}$ , the definition of the norm $\left \lVert \cdot \right \rVert _{\sigma ,W}$ tells us that the coefficients are bounded as

(3.7) \begin{equation} \left \lVert f_m \right \rVert _{W} \leqslant F \sigma ^{-m}. \end{equation}

This means that surfaces of $f$ converge geometrically for $\rho ^* \lt \sigma$ in $W$ , i.e.

(3.8) \begin{equation} \left \lVert \left . f \right |_{\rho =\rho ^*} \right \rVert _{W} \leqslant F\sum _{m=0}^{\infty } \left (\frac {\rho ^*}{\sigma }\right )^{m} = \frac {F}{1-\rho ^*/\sigma }. \end{equation}

For a more practical statement of pointwise convergence, we have the following proposition.

A direct consequence of the preceding statements is that both $(\sigma ,C^0)$ -analytic and $(\sigma ,H^2)$ -analytic functions are continuous in 3-D, but $(\sigma ,H^1)$ -analytic functions need not be.

Now, let us return to the question of ill-posedness. To define the norm of the input, let

(3.9) \begin{equation} \phi ^{(\text{IC})} = \sum _{m=0}^\infty \rho ^m \left (\phi _{m0} e^{-{\mathrm{i}} m \theta } + \phi _{mm} e^{{\mathrm{i}} m \theta } \right ). \end{equation}

This allows us to naturally define the norm of the input functions $\phi _{m0}$ and $\phi _{mm}$ via a single $(\sigma ,{H^q})$ -analytic norm. Using this, we prove that the problem is ill-posed in the following sense.

In other words, Theorem (3.6) tells us that smooth bounded perturbations in the input lead to unbounded deviations in the solution, even if the problem is initially prepared to be convergent. The characteristic form of the unbounded perturbations are high frequency in $s$ , which grow exponentially off-axis according to their wavenumber. When the near-axis expansion is discretised, this appears to be a poor condition number for the truncated problem. This motivates us to introduce a term in the near-axis expansion that damps the behaviour of the high-frequency modes.

Before we continue, we note there is a strong connection between the ill-posedness of the near-axis expansion and the ill-posedness of coil design. It is typically the case that magnetic fields with large gradients are difficult to approximate using plasma coils far from the boundary (Kappel et al. Reference Kappel, Landreman and Malhotra2024). This is because high frequencies in coil design decay quickly towards the surface of the plasma – the opposite view of the problem of high frequencies growing outward from the axis. The effect is that it is difficult to match high frequencies on the plasma boundary, and coil design codes also require some form of regularisation (Landreman Reference Landreman2017).

3.2. Regularisation

We have just seen in Theorem (3.6) that the near-axis expansion described in box (2.30) is ill-posed. Ill-posedness can potentially cause significant problems for numerical simulations, with the primary one being that discretisation refinement – both in toroidal resolution and in the order of expansion – will not lead to convergence to a solution. Instead, the deviation from the correct answer typically increases with refinement.

In the proof of Theorem (3.6) (§ A.3), the problematic perturbations to the input have small amplitude and high poloidal wavenumber, leading to exponential growth off the axis with a rate proportional to the wavenumber. Here, we would like to damp the problematic high-wavenumber modes while maintaining the fidelity of the low-wavenumber modes. We do so by adding a high-order differential operator in the toroidal and poloidal directions, as high-order derivatives affect high wavenumbers significantly more than low wavenumbers. Specifically, we propose the following version of Poisson’s equation (2.21) with regularisation:

(3.10) \begin{equation} \Delta \phi + \frac {\rho }{\sqrt {g}}\Delta _{\perp }P \phi = - \nabla \cdot \tilde {\boldsymbol{B}}. \end{equation}

The new term contains the regularising differential operator $\Delta _\perp P$ that is at least fourth order. The operator is composed of a product of the perpendicular Laplacian defined by

(3.11) \begin{equation} \Delta _{\perp } \phi = \frac {1}{\rho } \frac {\partial }{\partial \rho }\left ( \rho \frac {\partial \phi }{\partial \rho }\right ) + \frac {1}{\rho ^2} \frac {\partial ^2\phi }{\partial \theta ^2}, \end{equation}

and the regularising differential operator $P$ satisfying the following hypotheses.

We note that in hypothesis (i), the coefficients of $P$ do not depend on $\theta$ . This means poloidal Fourier modes ‘block diagonalise’ $P$ , i.e. for $f(s) \in H^{q + 2D}$ , we have $P(f(s) e^{{\mathrm{i}} m \theta }) = g(s) e^{{\mathrm{i}} m \theta }$ for some $g(s) \in H^{q}$ . This is sufficient for $P$ to map $(\sigma ,H^{q+2D})$ -analytic functions to $(\sigma ,H^{q})$ -analytic functions, while $\theta$ dependence would make the operator non-analytic. It also means that $P$ commutes with $\Delta _{\perp }$ , i.e. $P\Delta _\perp = \Delta _\perp P$ on sufficiently regular functions. Moreover, the hypotheses (ii) and (iii) are sufficient for $1+P$ to be invertible, which is necessary at each step of the iteration.

In practice, we use the operator

(3.15) \begin{equation} P = \left (-\frac {1}{K^2}\left (\frac {L}{2\pi \ell '} \frac {\partial }{\partial s}\right )^2\right )^{D} + \left (- \frac {1}{K^2} \frac {\partial ^2}{\partial \theta ^2} \right )^{D}, \end{equation}

where $K \gt 0$ is the characteristic wavenumber of regularisation, $2D\geqslant 2$ is the algebraic order of the frequency damping and $L = \oint \ell ' \hspace {1pt} \mathrm{d} s$ is the length of the magnetic axis. We have chosen the specific form of $P$ so that it is easy to implement numerically when the curve is in arclength coordinates (and therefore $\ell '$ is constant, satisfying the smoothness requirement). Indeed, any polynomial in the arclength derivative $(\ell ')^{-1}({\partial }/{\partial s})$ and the poloidal derivative $ ({\partial }/{\partial \theta })$ is simply inverted in Fourier space. Additionally, the parameter $K$ can be tuned from large (weak damping) to small (strong damping) to adjust the regularisation strength.

With regularisation, the new iterative form of the near-axis expansion (3.10) is

(3.16) \begin{equation} \left (m^2 + \frac {\partial ^2}{\partial \theta ^2}\right )(1 + P) \phi _m = -(\nabla \cdot \tilde {\boldsymbol{B}} + \Delta \phi _{\lt m})_{m-2}, \end{equation}

where we discuss the relevant regularity in Corollary (3.9). If we use the regularisation (3.15) and $s$ to be a scaled arclength coordinate so that $\ell ' = L/2\pi$ is constant, the componentwise version of the iteration is

(3.17) \begin{align} (m^2 - (2n-m)^2)\left (1 + \left (\frac { k}{K}\right )^{2D} + \left (\frac {2n-m}{K}\right )^{2D}\right )\phi _{mnk} = \nonumber \\ -(\nabla \cdot \tilde {\boldsymbol{B}} + \Delta \phi _{\lt m})_{m-2,n-1,k} \end{align}

for $0\lt n\lt m$ . We see that in the near-axis iteration, the regularisation damps the high-order modes by dividing by high-order polynomials in the poloidal and toroidal wavenumbers.

By construction, the regularisation has another benefit: the full problem can be represented as the divergence of a perturbed magnetic field $\boldsymbol{B}_K$ . If we let $G = Diag(1, \rho ^{-2}, 0)$ , we have

(3.18) \begin{equation} B_K^i = g^{ij} \frac {\partial \phi }{\partial q^j} + \tilde {B}^i + \frac {\rho }{\sqrt {g}} G^{ij} \frac {\partial (P \phi )}{\partial q^j} . \end{equation}

The fact that there is still an underlying divergence-free field is important for the relation between the near-axis expansion and Hamiltonian mechanics. Therefore, the regularised near-axis expansion of the field $\boldsymbol{B}_K$ can be expressed by the problem

(3.19) \begin{equation} \nabla \times \boldsymbol{B}_K = \boldsymbol J_K, \qquad \nabla \cdot \boldsymbol{B}_K = 0, \end{equation}

where $\boldsymbol J_K$ is a fictitious regularising current. In contrast, the regularisation means $\boldsymbol{B} = \nabla \phi + \tilde {\boldsymbol{B}}$ is not divergence-free. To find the flux surfaces, we note that the procedure in box (2.60) in no way depends on the magnetic field being curl-free. As such, we perform the same steps for finding flux coordinates for $\boldsymbol{B}_K$ as with $\boldsymbol{B}$ .

3.3. Convergence of the regularised expansion

Consider a PDE of the form

(3.20) \begin{equation} \Delta _\perp (1+P)\phi = f + L^{(a)}\phi , \end{equation}

where $P$ satisfies Hypotheses (3.7), $f$ is $(\sigma _0,H^{q})$ -analytic for $q\geqslant 0$ and $\sigma _0\gt 0$ , and

(3.21) \begin{align} L^{(a)} &= a^{(1)}\phi + \frac {\partial a^{(2)}}{\partial \rho } \frac {\partial \phi }{\partial \rho } + \frac {1}{\rho ^2} \frac {\partial a^{(2)}}{\partial \theta }\frac {\partial \phi }{\partial \theta } + a^{(3)} \frac {\partial \phi }{\partial \theta } \nonumber \\ &\quad+ a^{(4)} \frac {\partial \phi }{\partial s} + a^{(5)} \frac {\partial ^2\phi }{\partial \theta ^2} + a^{(6)} \frac {\partial ^2\phi }{\partial \theta \partial s} + a^{(7)} \frac {\partial ^2 \phi }{\partial s^2}, \end{align}

and each $a^{(j)}$ is a $(\Sigma , C^{q})$ -analytic for $\Sigma \gt \sigma _0$ . By a straightforward application of chain rule, we see that (3.10) satisfies this form for $\boldsymbol r_0 \in H^{4+q}$ and $B_0 \in H^{1+q}$ , where $\Sigma \leqslant \min \kappa ^{-1}$ (see proof of Corollary 3.9 in Appendix A.5). Equation (3.20) also encompasses other coordinate and regularisation assumptions of the near-axis expansion, including those not defined by Frenet–Serret coordinates.

The implicit near-axis iteration for (3.20) can be written as

(3.22) \begin{equation} (\Delta _\perp (1+ P) \phi )_m = \left [f + L^{(a)}\right ]_{m}\!. \end{equation}

The form of the iteration automatically respects the Fredholm condition, so it is solvable at each order (see the proof of Theorem (3.8)). So, to solve at each order of the iteration, we can invert to find that

(3.23) \begin{equation} \phi _m = \phi ^{(\text{IC})}_m + \left( (1+P)^{-1} \Delta _\perp ^{+} \left [f + L^{(a)}\phi \right ]\right )_{m}\!, \end{equation}

where $\phi ^{(\text{IC})}$ is $(\sigma _0,H^{q+2})$ -analytic of the form (3.9) and we define the pseudoinverse $\Delta ^+_\perp$ as

(3.24) \begin{equation} \Delta _\perp ^{+} (\rho ^m e^{{\mathrm{i}} (2n-m) \theta }) = \frac {1}{(m+2)^2 - (2n-m)^2} \rho ^{m+2} e^{{\mathrm{i}} (2n-m) \theta } \qquad \text{ for }0\leqslant n \leqslant m. \end{equation}

Given this iteration, we find the following theorem for its convergence.

Theorem 3.8. Let $0\lt \sigma _0\lt \Sigma$ , $q \geqslant 1$ and $D\geqslant 1$ . Then, let $f$ be $(\sigma _0,H^q)$ -analytic, $a^{(j)}$ for $1\leqslant j \leqslant 7$ be $(\Sigma , C^q)$ -analytic, $\phi ^{(\text{IC})}$ as defined by (3.9) be $(\sigma _0,H^{q+2D})$ -analytic and $P$ satisfy Hypotheses (3.7). Then, there is a unique $(\sigma ,H^{q+2D})$ -analytic near-axis solution $\phi$ of (3.20) with initial data $\phi ^{(\text{IC})}$ . Moreover, this solution operator is continuous in $a^{(j)}$ and satisfies the bound

(3.25) \begin{equation} \left \lVert \phi \right \rVert _{\sigma ,H^{q+2D}} \leqslant C \left \lVert f \right \rVert _{\sigma _0,H^q} + C \left \lVert \phi ^{(\text{IC})} \right \rVert _{\sigma _0,H^{q+2D}}\!. \end{equation}

This theorem can be translated to the Frenet–Serret problem.

Corollary 3.9. Let $\sigma \gt 0$ , $q \geqslant 1$ and $D\geqslant 1$ . Then, let $\boldsymbol r_0 \in C^{4+q}$ in scaled arclength coordinates with $\ell ' = L/2\pi \gt 0$ and $\kappa \gt 0$ , $B_0 \in H^{1+q}$ , $\phi ^{(\text{IC})}$ be $(\sigma _0,H^{q+2D})$ -analytic, and $P$ be defined as in ( 3.15 ). Then, for some $\sigma \gt 0$ satisfying $\sigma \leqslant \sigma _0$ and $\sigma \lt \min \kappa ^{-1}$ , the near-axis solution $\phi$ of ( 3.10 ) is $(\sigma ,H^{q+2D})$ -analytic, continuous in $\boldsymbol r_0$ and satisfies

(3.26) \begin{equation} \left\lVert{\phi}\right\rVert_{\sigma,H^{q+2D}} \leq C \left\lVert{B_0}\right\rVert_{H^{1+q}} + C \left\lVert{{\phi^{(\text{IC})}}}\right\rVert_{\sigma_0,H^{q+2D}}. \end{equation}

Furthermore, the associated divergence-free field $\boldsymbol{B}_K$ is well-defined and is $(\sigma ',H^q)$ -analytic for all $\sigma '\lt \sigma$ .

There are two primary reasons why Theorem (3.8) is useful for computation. First, it gives a guarantee that the norm of the inputs controls the size of the output. For instance, in the context of stellarator optimisation, if an appropriate Sobolev penalty is put on $\boldsymbol r_0$ , $B_0$ and $\phi ^{(\text{IC})}$ , then one can expect the output to be appropriately bounded. Second, because the output is $(\sigma ,H^{q+2D})$ -analytic, it tells us that truncations of the near-axis expansion are approximations of the true solution. So, we can be more confident that finite asymptotic series are approximately correct, at least for the regularised problem.

It is worth noting that while Theorem (3.8) tells us there is a solution to the regularised problem, it does not tell us that the solution solves the original problem. To address this, we can develop an a posteriori handle on the error.

Proposition 3.10. Consider the hypotheses of Theorem ( 3.8 ) with $q\geqslant 2$ . Additionally, suppose $0\lt \sigma '\lt \sigma$ , $L=\Delta _\perp - L^{(a)}$ is a second-order negative-definite operator on $H^{2}_{0}(\Omega _{\sigma '}^0)$ , and $\phi ^{\sigma '} = \phi (\sigma ',\theta ,s)$ . Then, the solution $\phi$ of the near-axis expansion is the unique solution in $H^{q+2D}(\Omega _{\sigma '}^0)$ of the boundary value problem

(3.27) \begin{equation} P \Delta _\perp \phi + L \phi = f, \qquad \left . \phi \right |_{\rho = \sigma '} = \phi ^{\sigma '}\!. \end{equation}

Moreover, let $\tilde \phi$ be the solution to

(3.28) \begin{equation} L \tilde \phi = f, \qquad \left . \tilde \phi \right |_{\rho =\sigma '} = \phi ^{\sigma '}\!. \end{equation}

Then,

(3.29) \begin{equation} \left \lVert \phi - \tilde \phi \right \rVert _{H^{q}\big(\Omega _{\sigma '}^0\big)} \leqslant C \left \lVert P \Delta _\perp \phi \right \rVert _{H^{q-2}\big(\Omega _{\sigma '}^0\big)}\!. \end{equation}

In other words, this tells us that given a solution to the regularised problem, we can bound the distance to a non-regularised boundary value problem via the norm of $P \Delta _\perp \phi$ . This applies directly to the Frenet–Serret case because the Laplacian is negative-definite. So, given a solution, this gives us an estimate of the error. As before, we can summarise the new problem (cf. box (2.30)):

(3.30)

5.1. Equilibrium initialisation

Figure 3. Coil sets for the rotating ellipse and Landreman–Paul examples. The colour indicates the normalised $\boldsymbol{B} \cdot \boldsymbol N$ error on the outer closed flux surface.

A major task in computing high-order near-axis expansions is choosing the input coefficients $\phi ^{(\text{IC})}$ in box (3.30). While low-order expansions can often be expressed in physically intuitive variables parametrising the rotation and stretching of elliptical magnetic surfaces, it is not as intuitive how to determine the high-order coefficients of $\phi ^{(\text{IC})}$ . In practice, the best option for finding equilibria would likely be via optimisation. However, for the purposes of demonstration, we initialise our inputs by a more direct method: via magnetic coils (see figure 3).

The primary advantage of using coils for equilibrium initialisation is accuracy. The accuracy comes from the fact that the coil field can be expanded analytically about the axis, giving a direct input to the near-axis expansion. This circumvents the potentially error-prone problem of interpolating stellarator equilibria. Coils also provide an accurate ground truth to which to compare our equilibrium. In the following subsections, we use this to assess the accuracy of the near-axis expansion, both close to the axis and farther away.

The coil optimisation method employed here follows the approach described by Wechsung et al. (Reference Wechsung, Landreman, Giuliani, Cerfon and Stadler2022) and Jorge et al. (Reference Jorge, Giuliani and Loizu2024) using the code SIMSOPT (Landreman et al. Reference Landreman, Medasani, Wechsung, Giuliani, Jorge and Zhu2021). Coils are modelled as single closed 3-D filaments of current $\boldsymbol \Gamma ^{(i)} : \mathbb{T} \to \mathbb{R}^3$ . Each coil $i$ is modelled as a periodic function in Cartesian coordinates where

(5.1) \begin{equation} \boldsymbol \Gamma ^{(j)}(s)=\boldsymbol c_{0}^{(j)}+\sum _{\ell =1}^{N_F}\left [\boldsymbol c_{\ell }^{(j)}\cos (\ell \theta )+\boldsymbol s_{\ell }^{(j)}\sin (\ell \theta )\right ], \end{equation}

where each $\boldsymbol c_{\ell }^{(j)}, \boldsymbol s_{\ell }^{(j)} \in \mathbb{R}^3$ , yielding a total of $3\times (2N_F +1)$ degrees of freedom per coil. In this work, we used $N_F=12$ , with four coils per half-field period for Landreman–Paul case and eight coils per half-field period for the ellipse. The degrees of freedom for the coil shapes are then

(5.2) \begin{equation} \boldsymbol x_{\text{coils}}=[\boldsymbol c_{l}^{(j)},\boldsymbol s_{l}^{(j)}, I_j], \end{equation}

with $I_j$ the current that goes through each coil. We take advantage of stellarator and rotational symmetries to only optimise a set of $N_c$ coils per half-field period. This leads to a total of $2 \times n_{\text{fp}} \times N_c$ modular coils, where $n_{\text{fp}}$ is the number of toroidal field periods with $n_{\text{fp}} = 2$ for Landreman–Paul and $n_{\text{fp}} = 5$ for the rotating ellipse. The remaining coils are determined by symmetry. The magnetic field $ \textbf{B}_{\text{ext}}$ of each coil is evaluated using the Biot–Savart law

(5.3) \begin{equation} \boldsymbol{B}_{\mathrm{coil}}(\boldsymbol{r}) = \frac {\mu _0}{4\pi } \sum _{j=1}^{2 n_{\mathrm{fp}} N_c} \int _{0}^{2\pi } \frac {I_j \frac {\partial \boldsymbol \Gamma _j}{\partial s'} \times (\boldsymbol \Gamma _j(s') - \boldsymbol{r})}{\left | \boldsymbol \Gamma _j(s') - \boldsymbol{r} \right |^3} \hspace {1pt} \mathrm{d} s', \end{equation}

where $s'$ parametrises the coil curve. Each coil is divided into 150 quadrature points, and the cost functions used to regularise the optimisation problem use the minimum distance between two coils, the length of each coil, their curvature and mean-squared curvature (see Wechsung et al. Reference Wechsung, Landreman, Giuliani, Cerfon and Stadler2022).

Using the coil magnetic field (5.3), we find the magnetic axis $\boldsymbol{r}_0$ via a shooting method. Then, using the near-axis coordinate representation of $\boldsymbol{r}$ in (2.2), we expand the quadrature rule of (5.3) using the operations in Appendix B to find a near-axis expansion for the magnetic field $\boldsymbol{B}(\rho ,\theta ,s)$ . Given the near-axis field, it is straightforward to compute $B_0$ by

(5.4) \begin{equation} B_0(s) = \boldsymbol{t}(s) \cdot \boldsymbol{B}(0,0,s), \end{equation}

and $\phi$ is found by a near-axis expansion of the path integral (note that $\phi =0$ on the axis)

(5.5) \begin{equation} \phi (\rho ,\theta ,s) = \int _{0}^\rho (\boldsymbol{B}(\rho ',\theta ,s) - \tilde {\boldsymbol{B}}(\rho ',\theta ,s)) \cdot (\cos (\theta ) \boldsymbol{n}(s) + \sin (\theta ) \boldsymbol{b}(s)) \hspace {1pt} \mathrm{d} \rho '. \end{equation}

Finally, the coefficients $\phi _{m0}$ and $\phi _{mm}$ of this are used as input for $\phi ^{(\text{IC})}$ . The input is computed to the $N_\rho =9$ orders in $\rho$ with $N_{s,\mathrm{RE}}=100$ and $N_{,s\mathrm{LP}}=50$ Fourier modes in $s$ for the rotating ellipse and Landreman–Paul, respectively (see (4.1)), where we use the subscript ‘RE’ for the rotating ellipse and ‘LP’ for Landreman–Paul wherever necessary.

To begin our analysis of the examples, we consider the inputs to the near-axis expansion. Corollary (3.9) suggests that there are two length scales dictated by the inputs. First is the radius of curvature of the axis. Letting $\Sigma = \min _s \kappa ^{-1}(s)$ be the radius of curvature, we find

(5.6) \begin{equation} \Sigma_\textrm{RE} = 0.987, \qquad \Sigma_\textrm{LP} = 0.681. \end{equation}

The other length scale of interest is the radius of convergence $\sigma _0$ of $\phi ^{(\text{IC})}$ in the $(\sigma _0,H^{q+2D})$ -analytic norm. However, to the orders to which we compute, this radius appears to depend on the exponent $q+2D$ . To determine the most informative exponent, we consider the work by Kappel et al. (Reference Kappel, Landreman and Malhotra2024), where it was shown that the normalised gradient of the magnetic field is a strong predictor for plasma–coil separation. Because the $H^2(\mathbb{T}^2)$ norm measures the size of the second derivative of $\phi ^{(\text{IC})}_m$ (and therefore the gradient of the input magnetic field), we conjecture this corresponds to the most practical exponent.

To verify this, we first compute the minimum axis-to-coil distance $\sigma _{\mathrm{coil}}$ for both configurations to be

(5.7) \begin{equation} \sigma _{\mathrm{coil,RE}}=0.241, \qquad \sigma _{\mathrm{coil,LP}} = 0.350. \end{equation}

We note that $\sigma _{\mathrm{coil}}\lt \Sigma$ for both configurations, indicating that the distance-to-coil is the limiting factor for convergence (cf. Corollary (3.9)). In figure 4, we plot the $H^2$ norms of $\phi ^{(\text{IC})}_m$ versus $A \sigma _{\mathrm{coil}}^{-m}$ , where the coefficient $A$ is found via a best fit for both configurations. In both cases, we find there is remarkable agreement, indicating that the $H^2$ radius of convergence of $\phi ^{(\text{IC})}$ could be used as a proxy for distance-to-coils.

5.2. Magnetic field convergence

Figure 4. Plot of the coefficient norm $\lVert \phi ^{(\text{IC})}_m\rVert _{H^2}$ versus the order $m$ (markers) and best-fit lines $A \sigma _{\mathrm{coil}}^{-m}$ (lines), where $\sigma _{\mathrm{coil}}$ is the axis-to-coil distance.

Next, we compute the near-axis expansion via the procedure in box (3.30). We perform the expansion both without regularisation ( $P=0$ ) and with the regularisation operator in (3.15). For the regularised runs, we use $D=2$ throughout and vary $K$ between $10$ and $200$ to assess how the equilibrium changes between strong and weak regularisation, respectively.

As a first test of the output convergence, we compare the coil magnetic field against the unregularised expansion. Because the input is harmonic, we expect that the near-axis expansion will converge from Proposition (3.2) (unless floating-point errors overwhelm the solution, which is not observed to this order). We verify this by computing the $L^2$ magnetic field error on surfaces about the axis

(5.8) \begin{equation} \left \lVert \boldsymbol{B} - \boldsymbol{B}_{\mathrm{coil}} \right \rVert _{L^2(\partial \Omega _\sigma ^0)} = \frac {1}{4\pi ^2} \int _0^{2\pi } \int _{0}^{2\pi } \left | \boldsymbol{B}(\sigma ,\theta ,s) - \boldsymbol{B}_{\mathrm{coil}}(\sigma ,\theta ,s) \right |^2 \hspace {1pt} \mathrm{d} \theta \mathrm{d} s, \end{equation}

where $\boldsymbol{B}$ is computed from the near-axis expansion and $\boldsymbol{B}_{\mathrm{coil}}$ is computed directly from Biot–Savart.

In figure 5( a,b ), we plot the error (5.8) versus the normalised distance-from-axis $\sigma /\sigma _{\mathrm{coil}}$ . This error is computed for the approximation

(5.9) \begin{equation} \phi _{\leqslant N_\rho } \approx \sum _{m=2}^{N_\rho } \phi _m \rho ^m, \qquad \boldsymbol{B}_{\leqslant N_\rho -1} = \nabla \phi _{\leqslant N_\rho } + \tilde {\boldsymbol{B}}_{\leqslant N_\rho -1} \end{equation}

where $N_\rho$ is varied from $1$ to $8$ . For both configurations, we find that the error of the magnetic field obeys the expected power law

(5.10) \begin{equation} \left \lVert \boldsymbol{B}_{\leqslant N_\rho -1} - \boldsymbol{B}_{\mathrm{coil}} \right \rVert _{L^2(\partial \Omega _\sigma ^0)} = \mathcal{O}\big(\sigma ^{-N_\rho }\big). \end{equation}

The error curves for varying $N_\rho$ meet at $\sigma = \sigma _{\mathrm{coil}}$ , indicating that the output radius of convergence is limited by the coils. This tells us that the limit of convergence $\sigma = \sigma _0$ is achievable in Corollary (3.9).

Figure 5. (a,b) The error (5.8) as a function of the normalised distance from axis $\sigma /\sigma _{\mathrm{coil}}$ for varying orders of approximation $N_\rho$ . (c,d) The error (5.8) as a function of $\sigma /\sigma _{\mathrm{coil}}$ for varying values of the regularisation parameter $K$ ( $K=\infty$ is unregularised).

Turning to the effects of regularisation, we fix $N_\rho =9$ and plot the error (5.8) versus $\sigma /\sigma _{\mathrm{coil}}$ for varying $K$ between $10$ and $200$ in figure 5( c,d ). We also include the unregularised solution, labelled with $K=\infty$ . We find that as the regularisation becomes stronger ( $K$ decreases), the magnetic field loses fidelity near the core. We attribute to the increasing loss of accuracy of the high-wavenumber $s$ modes, while the low-wavenumber modes maintain accuracy. Then, far from the axis, the regularised error inflects to begin to agree with the rate of convergence of the unregularised solution. So, while the solution loses a high-wavenumber fidelity, the low wavenumbers maintain a similar level of accuracy. Comparing figures 5( a,b ) with figure 5( c,d ), we see that a regularised high-order expansion can achieve an equivalent error to an unregularised lower-order expansion near the axis while maintaining that fidelity far from the axis.

To address the role of regularisation more fully, however, we need to consider the fidelity of the expansion on less tuned inputs. To do this, we perturb $\boldsymbol{r}_0$ , $B_0$ and $\phi ^{(\text{IC})}$ by the random functions as

(5.11) \begin{align} \delta r_{0,i} &= \epsilon \sum _{\ell = -N_s}^{N_s} \frac {X_{i\ell }}{1 + (\ell /K_{\epsilon })^{4+q}} , \nonumber \\ \delta B_0 &= \epsilon \sum _{\ell = -N_s}^{N_s}\frac {Y_{\ell }}{1 + (\ell /K_{\epsilon })^{1+q}} , \nonumber \\ \delta \phi ^{(\text{IC})}_{mn} &= \epsilon \sum _{\ell = -N_s}^{N_s} \frac {Z_{mn\ell }}{1 + (m/K_{\epsilon })^{q+2D} + (\ell /K_{\epsilon })^{q+2D}} , && \quad n\in \{0,m\}, \end{align}

where $X_{i\ell }$ , $Y_\ell$ and $Z_{mn\ell }$ are independent and identically distributed (i.i.d.) unit normal random variables, $q = 2$ , and $\epsilon = 10^{-6}$ for the rotating ellipse and $\epsilon =10^{-4}$ for Landreman–Paul. We have chosen the regularity of the perturbation to align with the inputs in box (3.30).

Figure 6. Finite difference residual (5.12) as a function of the normalised distance from axis $\sigma /\sigma _{\mathrm{coil}}$ for the perturbed rotating ellipse and Landreman–Paul inputs (see (5.11)). For both plots, three lines are coloured and labelled, while the grey lines represent other values of $K$ interpolating between $K=10$ and $K=200$ .

In figure 6, we consider the accuracy of the solution to the perturbed problem for $K$ varying between $10$ and $200$ for both examples with $N_\rho = 9$ fixed. To measure the accuracy, we no longer have a coil set with which to compare the solution directly. So, we instead measure the residual of Poisson’s equation

(5.12) \begin{equation} \left \lVert \nabla \cdot \boldsymbol{B} \right \rVert _{L^2(\partial \Omega _\sigma ^0)} = \frac {1}{4\pi ^2} \int _0^{2\pi } \int _{0}^{2\pi } \left | \nabla \cdot (\nabla \phi + \tilde {\boldsymbol{B}}) \right |^2 \hspace {1pt} \mathrm{d} \theta \mathrm{d} s, \end{equation}

where we evaluate every derivative (including in the metric) via finite differences. For both perturbed examples, the best solution near the axis is the lightly regularised $K=200$ solution. However, beyond a certain radius between $0.3\sigma _{\mathrm{coil}}$ and $0.4 \sigma _{\mathrm{coil}}$ , more regularised solutions improve upon the less regularised ones in the finite difference metric. For our examples, we find that $K=26$ for the rotating ellipse and $K=33$ for Landreman–Paul are perhaps the best choices in practice. This figure potentially indicates a more general principle: the further from the axis one wants accuracy of the expansion, the more regularised the expansion likely has to be.

5.3. Magnetic coordinate convergence

Figure 7. Near-axis approximations of flux surfaces for varying orders of approximation $N_\rho$ (black); a Poincaré plot of the true coil magnetic field lines (colour). In the final $N_\rho = 9$ panel, we plot a circle with radius $\sigma _{\mathrm{coil,LP}}$ in red.

To compute straight field-line coordinates, we return to the unregularised Landreman–Paul configuration. Then, using the straight field-line magnetic coordinate equation from box (2.60), we compute the approximate coordinates $(\xi ,\eta )$ for $N_\rho$ varying between $2$ and $9$ , where we note the $N_\rho =2$ approximation of $\phi$ provides the leading-order field-line behaviour. To find flux surfaces, we then invert $(\xi ,\eta )$ to find the distance-to-axis coordinates $(x(\xi ,\eta ,s),y(\xi ,\eta ,s))$ , where magnetic surfaces are parametrised by $\xi ^2 + \eta ^2 = \psi$ .

In figure 7, we plot in black the computed surfaces on the $s=0$ Poincaré section for varying values of $N_\rho$ . For comparison, we plot the intersections of coil magnetic field lines in the background. At leading order, we see the surfaces are elliptical, while higher orders account for more shaping in the $\theta$ direction. Then, as the order increases beyond $5$ , the surfaces away from the core start diverging.

To investigate this divergence, we plot a red circle of constant radius $\sigma _{\mathrm{coil,LP}}$ in the $N_\rho = 9$ panel. We see that the circle appears to separate the divergent surfaces from the convergent ones. We believe this is the likely reason for the divergence; however, there are still other possibilities.

Figure 8. (a) The rotational transform and (b) the parametrisation error $R_{\mathrm{param}}$ defined in (5.13) as a function of the inboard $x$ distance (cf. figure 7) for varying orders of approximation $N_\rho$ .

To assess the errors of surfaces closer to the axis, we turn to a more quantitative measure. To do this, we first use the method from Ruth & Bindel (Reference Ruth and Bindel2024) in the SymplecticMapTools.jl package (Ruth Reference Ruth2024b ) to compute invariant circles $(x_{\mathrm{coil}}(\theta ), y_{\mathrm{coil}}(\theta ))$ and the rotational transform $\iota _{\mathrm{coil}}$ on the cross-section from the Poincaré plot trajectories. Then, as a function of the inboard distance $x$ from the axis (see figure 7), we compute rotational transform and parametrisation errors as

(5.13) \begin{align} R_{\iota } &= \left | \iota - \iota _{\mathrm{coil}} \right |, \nonumber \\ R_{\mathrm{param}} &= \left [\int _{0}^{2\pi } (x(\psi (x,0,0),\theta ,0)-x_{\mathrm{coil}}(\theta ))^2 + (y(\psi (x,0,0),\theta ,0)-y_{\mathrm{coil}}(\theta ))^2 \hspace {1pt} \mathrm{d} \theta \right ]^{1/2}. \end{align}

In figure 8, we plot both errors with varying $N_\rho$ . In both cases, the rotational transform and parametrisation converge to high accuracy near the core. However, they begin to diverge before the the outermost surface, agreeing with the visual divergence in figure 7.