2.1 Mercier’s coordinate system

In this section, leveraging the work in Mercier (Reference Mercier1964) and Solov’ev & Shafranov (Reference Solov’ev and Shafranov1970), we construct an orthogonal coordinate system associated with a particular field line of force $\boldsymbol{r}_{0}(s)$, which is taken to be the magnetic axis curve. We let $L$ denote the total length of $\boldsymbol{r}_{0}$, such that $0\leqslant s\leqslant L$. The unit tangent vector $\boldsymbol{t}$ is defined as

(2.1)$$\begin{eqnarray}\boldsymbol{t}=\boldsymbol{r}_{0}^{\prime }(s),~0\leqslant s\leqslant L.\end{eqnarray}$$

Using the fact that $\boldsymbol{t}^{\prime }(s)$ is orthogonal to $\boldsymbol{t}$, the unit normal vector $\boldsymbol{n}$ is defined as $\boldsymbol{n}=\boldsymbol{t}^{\prime }(s)/\unicode[STIX]{x1D705}$ with $\unicode[STIX]{x1D705}=|\boldsymbol{t}^{\prime }(s)|=|(\boldsymbol{t}\boldsymbol{\cdot }\unicode[STIX]{x1D735})\boldsymbol{t}|$ the curvature, while the unit binormal vector $\boldsymbol{n}$ obeys $\boldsymbol{b}=\boldsymbol{t}\times \boldsymbol{n}$. The triad $(\boldsymbol{t},\boldsymbol{n},\boldsymbol{b})$ forms a right-handed system of orthogonal unit vectors, usually called the Frenet–Serret frame (Spivak Reference Spivak1999), which obey the following set of first-order differential equations

(2.2)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{t}^{\prime }(s)=\unicode[STIX]{x1D705}\boldsymbol{n}, & \displaystyle\end{eqnarray}$$

(2.3)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{n}^{\prime }(s)=-\unicode[STIX]{x1D705}\boldsymbol{t}+\unicode[STIX]{x1D70F}\boldsymbol{b}, & \displaystyle\end{eqnarray}$$

(2.4)$$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{b}^{\prime }(s)=-\unicode[STIX]{x1D70F}\boldsymbol{n}, & \displaystyle\end{eqnarray}$$

with $\unicode[STIX]{x1D70F}$ the torsion. Explicit expressions for the curvature and torsion when the curve $\boldsymbol{r}_{0}$ is parametrized in terms of a parameter $t$ other than the arclength $s$ (e.g. the toroidal angle $\unicode[STIX]{x1D6F7}$ in cylindrical coordinates) can be obtained using

(2.5)$$\begin{eqnarray}\unicode[STIX]{x1D705}(t)={\displaystyle \frac{|\boldsymbol{r}^{\prime }(t)\times \boldsymbol{r}^{\prime \prime }(t)|}{|\boldsymbol{r}^{\prime }(t)|^{3}}},\end{eqnarray}$$

and

(2.6)$$\begin{eqnarray}\unicode[STIX]{x1D70F}(t)={\displaystyle \frac{\left(\boldsymbol{r}^{\prime }(t)\times \boldsymbol{r}^{\prime \prime }(t)\right)\boldsymbol{\cdot }\boldsymbol{r}^{\prime \prime \prime }(t)}{|\boldsymbol{r}^{\prime }(t)\times \boldsymbol{r}^{\prime \prime }(t)|^{2}}}.\end{eqnarray}$$

It can be shown that any curve (such as the magnetic axis) can be described only with $\unicode[STIX]{x1D705}$ and $\unicode[STIX]{x1D70F}$ (see e.g. Spivak (Reference Spivak1999)). Given $\unicode[STIX]{x1D705}$ and $\unicode[STIX]{x1D70F}$, the Frenet–Serret frame can then be found using (2.2)–(2.4) and a set of initial conditions. An example of a magnetic axis curve is shown in figure 1, together with the Frenet–Serret unit vectors, namely the tangent (blue), normal (green) and binormal (red) unit vectors.

Figure 1. Example of a magnetic surface using Mercier’s coordinates, where the radial coordinate $\unicode[STIX]{x1D70C}$ is prescribed as a function of $s$ and $\unicode[STIX]{x1D714}=\unicode[STIX]{x1D703}+\unicode[STIX]{x1D6FE}$ in order to obtain a toroidal shape with elliptical cross-section. The magnetic axis curve is shown in black, together with the Frenet–Serret unit vectors, namely the tangent (cyan), normal (grey) and binormal (orange) unit vectors. The blue and yellow lines in the outermost surface denote the curves with constant $\unicode[STIX]{x1D703}$ and $s$, respectfully.

We now let $\unicode[STIX]{x1D70C}$ denote the distance between an arbitrary point $\boldsymbol{r}$ to the axis $\boldsymbol{r}_{0}$ in the plane $(\boldsymbol{n},\boldsymbol{b})$ orthogonal to the axis, and let $\unicode[STIX]{x1D703}$ be the angle measured from the normal to $\boldsymbol{r}-\boldsymbol{r}_{0}$. The radius vector $\boldsymbol{r}$ can then be written as

(2.7)$$\begin{eqnarray}\boldsymbol{r}=\boldsymbol{r}_{0}(s)+\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}\boldsymbol{n}(s)+\unicode[STIX]{x1D70C}\sin \unicode[STIX]{x1D703}\boldsymbol{b}(s).\end{eqnarray}$$

We note that the set of coordinates $(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703},s)$ can be made orthogonal by introducing the angle $\unicode[STIX]{x1D714}$, defined by

(2.8)$$\begin{eqnarray}\unicode[STIX]{x1D714}=\unicode[STIX]{x1D703}+\unicode[STIX]{x1D6FE}(s),~\unicode[STIX]{x1D6FE}(s)=\int _{0}^{s}\unicode[STIX]{x1D70F}(s^{\prime })\,\text{d}s^{\prime },\end{eqnarray}$$

with $\unicode[STIX]{x1D6FE}$ the integrated torsion. An example of a magnetic surface constructed using Mercier’s formalism is shown in figure 1, where the radial coordinate $\unicode[STIX]{x1D70C}$ was chosen to be a function of $s$ and $\unicode[STIX]{x1D714}$ in order to obtain poloidal cross-sections with an elliptical shape. The explicit expression of $\unicode[STIX]{x1D70C}(s,\unicode[STIX]{x1D714})$ for the case of an elliptical expression can be found in § 3. The determinant of the metric tensor is given by

(2.9)$$\begin{eqnarray}\sqrt{g}=\unicode[STIX]{x1D70C}h_{s}\end{eqnarray}$$

in both the $(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703},s)$ and $(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D714},s)$ coordinate systems. However, the metric tensor $\unicode[STIX]{x1D628}_{ij}$ is diagonal only when expressed in terms of the $(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D714},s)$ coordinates, with $\unicode[STIX]{x1D628}_{ij}=\text{diag}(1,\unicode[STIX]{x1D70C}^{2},h_{s}^{2})$. In (2.9), we defined $h_{s}$ as

(2.10)$$\begin{eqnarray}h_{s}=1-\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos (\unicode[STIX]{x1D714}-\unicode[STIX]{x1D6FE}).\end{eqnarray}$$

In the following, we assume that the plasma boundary is close enough to the axis such that $\unicode[STIX]{x1D70C}<1/\unicode[STIX]{x1D705}$ leading to a non-zero Jacobian across the whole plasma volume. Finally, we introduce a Cartesian coordinate system in the $(\boldsymbol{n},\boldsymbol{b})$ plane by defining $x=\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}$ and $y=\unicode[STIX]{x1D70C}\sin \unicode[STIX]{x1D703}$, yielding

(2.11)$$\begin{eqnarray}\boldsymbol{r}=\boldsymbol{r}_{0}(s)+x\boldsymbol{n}(s)+y\boldsymbol{b}(s).\end{eqnarray}$$

We note that the position vector in (2.11) coincides with the definition in the inverse coordinate method (Garren & Boozer Reference Garren and Boozer1991a) only to lowest order in $\unicode[STIX]{x1D716}$. A key difference between the direct method pursued here and the inverse method in (Garren & Boozer Reference Garren and Boozer1991a) is that the latter includes a finite contribution in the $\boldsymbol{t}$ direction in (2.11), leading to distinct vectors $\boldsymbol{r}-\boldsymbol{r}_{0}(s)$ between both the direct and the inverse method at the same point $\boldsymbol{r}$.

2.2 Power series expansion

In this section, we show how to construct the asymptotic series related to a physical quantity $f$ in terms of $\unicode[STIX]{x1D716}$ and $\unicode[STIX]{x1D70C}$. As shown below, such construction imposes conditions on the series coefficients of $f$. In § 3, we prove that the magnetic field satisfies such conditions in the vacuum case at all orders. In what follows, we normalize $\unicode[STIX]{x1D70C},x$ and $y$ to the characteristic perpendicular scale $a$ and the quantities $\unicode[STIX]{x1D705},\unicode[STIX]{x1D70F}$ and $s$ to the characteristic parallel scale $R$. The magnetic field is normalized to a constant $\overline{B}=\int _{0}^{L}B_{0}(s)\,\text{d}s/L$, with $B_{0}(s)$ the magnetic field on axis and $\unicode[STIX]{x1D713}$ is normalized to $\overline{B}R^{2}$. The asymptotic expansion of a function $f$ is then constructed by noting that any analytic function $f$ has a Taylor expansion near the origin $(x,y)=(0,0)$ of the form

(2.12)$$\begin{eqnarray}f(x,y,s)=\mathop{\sum }_{p=0}^{\infty }\mathop{\sum }_{j=0}^{\infty }f_{pj}(s)x^{p}y^{j}\unicode[STIX]{x1D716}^{p+j}.\end{eqnarray}$$

Similarly, using (2.11), we write $f$ as

(2.13)$$\begin{eqnarray}f(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D714},s)=\mathop{\sum }_{n=0}^{\infty }f_{n}(\unicode[STIX]{x1D714},s)\unicode[STIX]{x1D716}^{n}\unicode[STIX]{x1D70C}^{n}.\end{eqnarray}$$

The equality between the two expansions in (2.12) and (2.13) yields

(2.14)$$\begin{eqnarray}f_{n}(\unicode[STIX]{x1D714},s)=\mathop{\sum }_{p=0}^{n}f_{pn-p}(s)\cos ^{p}(\unicode[STIX]{x1D703})\sin ^{(n-p)}(\unicode[STIX]{x1D703})=\mathop{\sum }_{p=0}^{n}f_{np}^{c}(s)\cos (p\unicode[STIX]{x1D703})+f_{np}^{s}(s)\sin (p\unicode[STIX]{x1D703}).\end{eqnarray}$$

Equation (2.14) shows that $f_{n}$ can be written as a Fourier series in terms of $\cos p\unicode[STIX]{x1D703}$ and $\sin p\unicode[STIX]{x1D703}$ with only even or odd values of $p$ in the range $0\leqslant p\leqslant n$ depending on whether $n$ is even or odd (Kuo-Petravic & Boozer Reference Kuo-Petravic and Boozer1987). This is proven in appendix A. In the following, when dealing with Mercier’s coordinates $(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D714},s)$, analyticity is shown in the form of (2.14) rather than (2.12).

We expand the components of the magnetic field $\boldsymbol{B}$ as $\boldsymbol{B}=B_{\unicode[STIX]{x1D70C}}\boldsymbol{e}_{\unicode[STIX]{x1D70C}}+B_{\unicode[STIX]{x1D714}}\boldsymbol{e}_{\unicode[STIX]{x1D714}}+B_{s}\boldsymbol{e}_{s}$ with $\boldsymbol{e}_{\unicode[STIX]{x1D70C}}=\cos \unicode[STIX]{x1D703}\boldsymbol{n}+\sin \unicode[STIX]{x1D703}\boldsymbol{b}$, $\boldsymbol{e}_{\unicode[STIX]{x1D714}}=-\sin \unicode[STIX]{x1D703}\boldsymbol{n}+\cos \unicode[STIX]{x1D703}\boldsymbol{b}$ and $\boldsymbol{e}_{s}=\boldsymbol{t}$ using (2.13), i.e.

(2.15)$$\begin{eqnarray}\displaystyle B_{\unicode[STIX]{x1D70C}}=\mathop{\sum }_{n=0}^{\infty }B_{\unicode[STIX]{x1D70C}n}(\unicode[STIX]{x1D714},s)\unicode[STIX]{x1D716}^{n}\unicode[STIX]{x1D70C}^{n}, & & \displaystyle\end{eqnarray}$$

and similarly for $B_{\unicode[STIX]{x1D714}},B_{s}$ and $\unicode[STIX]{x1D713}$. The fields $\boldsymbol{B},\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D6FC}$ are split into a vacuum and a finite $\unicode[STIX]{x1D6FD}$ component as

(2.16)$$\begin{eqnarray}\boldsymbol{B}=\boldsymbol{B}^{0}+\boldsymbol{B}^{1},\end{eqnarray}$$

with $\boldsymbol{B}^{0}$ the vacuum magnetic field in the absence of a plasma and $\boldsymbol{B}^{1}$ a linear perturbation in $\unicode[STIX]{x1D6FD}$. A similar split is applied to $\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D6FC}$. The linear perturbations satisfy

(2.17)$$\begin{eqnarray}\left|{\displaystyle \frac{\boldsymbol{B}^{1}}{\boldsymbol{B}^{0}}}\right|\sim {\displaystyle \frac{\unicode[STIX]{x1D713}^{1}}{\unicode[STIX]{x1D713}^{0}}}\sim {\displaystyle \frac{\unicode[STIX]{x1D6FC}^{1}}{\unicode[STIX]{x1D6FC}^{0}}}\sim \unicode[STIX]{x1D6FD}\ll 1.\end{eqnarray}$$

We first consider the vacuum case and obtain a set of equations for the coefficients $B_{\unicode[STIX]{x1D70C}n}^{0},B_{\unicode[STIX]{x1D714}n}^{0},B_{sn}^{0},\unicode[STIX]{x1D713}_{n}^{0}$ and $\unicode[STIX]{x1D6FC}_{n}^{0}$.

3.1 Vacuum magnetic scalar potential

We now expand $\unicode[STIX]{x1D719}$ using (2.13) and collect terms of the same order in $\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}$ in (3.3) in order to obtain a single equation for $\unicode[STIX]{x1D719}_{n}$. In the following, we define $\dot{\unicode[STIX]{x1D719}}=\unicode[STIX]{x2202}_{\unicode[STIX]{x1D714}}\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D719}^{\prime }=\unicode[STIX]{x2202}_{s}\unicode[STIX]{x1D719}$ for the partial derivatives with respect to $\unicode[STIX]{x1D714}$ and $s$, respectively. We first focus on the $O(\unicode[STIX]{x1D716}^{2})$ term of (3.3)

(3.6)$$\begin{eqnarray}\ddot{\unicode[STIX]{x1D719}}_{2}+4\unicode[STIX]{x1D719}_{2}=-B_{0}^{\prime }(s).\end{eqnarray}$$

The homogeneous solution of (3.6) can be written as a linear combination of $\sin 2\unicode[STIX]{x1D714}$ and $\cos 2\unicode[STIX]{x1D714}$ terms, while the particular solution is given by $-B_{0}^{\prime }/4$, i.e.

(3.7)$$\begin{eqnarray}\unicode[STIX]{x1D719}_{2}=-{\displaystyle \frac{B_{0}^{\prime }}{4}}+c_{1}\sin 2\unicode[STIX]{x1D714}+c_{2}\cos 2\unicode[STIX]{x1D714}.\end{eqnarray}$$

For convenience, and to later obtain a direct relation between the coefficients of $\sin 2\unicode[STIX]{x1D714}$ and $\cos 2\unicode[STIX]{x1D714}$ and the geometric parameters of the flux surface function, we introduce the functions $\unicode[STIX]{x1D6FF}(s),\unicode[STIX]{x1D702}(s)$ and $\unicode[STIX]{x1D707}(s)=\tanh \unicode[STIX]{x1D702}(s)$, and write $\unicode[STIX]{x1D719}_{2}$ as

(3.8)$$\begin{eqnarray}\unicode[STIX]{x1D719}_{2}={\displaystyle \frac{B_{0}}{2}}\left[(\ln B_{0}^{-1/2})^{\prime }+\unicode[STIX]{x1D707}u^{\prime }\sin 2u-{\displaystyle \frac{\unicode[STIX]{x1D702}^{\prime }}{2}}\cos 2u\right],~u=\unicode[STIX]{x1D714}-\unicode[STIX]{x1D6FE}(s)+\unicode[STIX]{x1D6FF}(s).\end{eqnarray}$$

The integration constants $\unicode[STIX]{x1D6FF}(s)$ and $\unicode[STIX]{x1D702}(s)$ that characterize the $O(\unicode[STIX]{x1D716}^{2})$ scalar potential $\unicode[STIX]{x1D719}$ are arbitrary (periodic) functions of $s$ can be used to impose additional constraints on the magnetic field $\boldsymbol{B}$, such as quasisymmetry or omnigeneity.

Focusing on the $O(\unicode[STIX]{x1D716}^{3})$ term of (3.3), we obtain the following equation for $\unicode[STIX]{x1D719}_{3}$

(3.9)$$\begin{eqnarray}9\unicode[STIX]{x1D719}_{3}+\ddot{\unicode[STIX]{x1D719}}_{3}=-B_{0}\cos \unicode[STIX]{x1D703}\unicode[STIX]{x1D705}^{\prime }+\unicode[STIX]{x1D705}[2\unicode[STIX]{x1D719}_{22}^{c}\cos (\unicode[STIX]{x1D703}+2\unicode[STIX]{x1D6FF})+2\unicode[STIX]{x1D719}_{22}^{s}\sin (\unicode[STIX]{x1D703}+2\unicode[STIX]{x1D6FF})-B_{0}\unicode[STIX]{x1D70F}\sin \unicode[STIX]{x1D703}-2B_{0}^{\prime }\cos \unicode[STIX]{x1D703}],\end{eqnarray}$$

with $\unicode[STIX]{x1D719}_{22}^{c}=-B_{0}\unicode[STIX]{x1D702}^{\prime }/4$ and $\unicode[STIX]{x1D719}_{22}^{s}=B_{0}\unicode[STIX]{x1D707}u^{\prime }/2$. Similarly to (3.8), we write $\unicode[STIX]{x1D719}_{3}$ as

(3.10)$$\begin{eqnarray}\unicode[STIX]{x1D719}_{3}=\unicode[STIX]{x1D719}_{31}^{c}\cos u+\unicode[STIX]{x1D719}_{31}^{s}\sin u+\unicode[STIX]{x1D719}_{33}^{c}\cos 3u+\unicode[STIX]{x1D719}_{33}^{s}\sin 3u.\end{eqnarray}$$

Plugging the form of (3.10) in (3.9), we obtain for the particular solution

(3.11)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{31}^{c}(s)=-{\displaystyle \frac{\unicode[STIX]{x1D705}B_{0}}{8}}\left({\displaystyle \frac{5}{2}}{\displaystyle \frac{B_{0}^{\prime }}{B_{0}}}+{\displaystyle \frac{k^{\prime }}{k}}+{\displaystyle \frac{\unicode[STIX]{x1D702}^{\prime }}{2}}\cos 2\unicode[STIX]{x1D6FF}-\unicode[STIX]{x1D707}u^{\prime }\sin 2\unicode[STIX]{x1D6FF}\right), & \displaystyle\end{eqnarray}$$

(3.12)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{31}^{s}(s)={\displaystyle \frac{\unicode[STIX]{x1D705}B_{0}}{8}}\left(-\unicode[STIX]{x1D70F}+{\displaystyle \frac{\unicode[STIX]{x1D702}^{\prime }}{2}}\sin 2\unicode[STIX]{x1D6FF}+\unicode[STIX]{x1D707}u^{\prime }\cos 2\unicode[STIX]{x1D6FF}\right), & \displaystyle\end{eqnarray}$$

with both $\unicode[STIX]{x1D719}_{33}^{c}(s)$ and $\unicode[STIX]{x1D719}_{33}^{s}(s)$ integration constants from the homogeneous solution.

To obtain an expression for the solution of $\unicode[STIX]{x1D719}_{n}$ at arbitrary order in $\unicode[STIX]{x1D716}$, we first define the Laplacian operator $D$ in polar $(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D714})$ coordinates multiplied by $\unicode[STIX]{x1D70C}^{2}$ as

(3.13)$$\begin{eqnarray}D\unicode[STIX]{x1D719}=\left[\unicode[STIX]{x1D70C}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}\left(\unicode[STIX]{x1D70C}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}\right)+{\displaystyle \frac{\unicode[STIX]{x2202}^{2}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}^{2}}}\right]\unicode[STIX]{x1D719}=\mathop{\sum }_{n=2}^{\infty }\left(n^{2}\unicode[STIX]{x1D719}_{n}+\ddot{\unicode[STIX]{x1D719}}_{n}\right)\unicode[STIX]{x1D716}^{n}\unicode[STIX]{x1D70C}^{n}.\end{eqnarray}$$

Laplace’s equation, equation (3.3), can then be written as

(3.14)$$\begin{eqnarray}D\unicode[STIX]{x1D719}=-\mathop{\sum }_{n=2}^{\infty }\unicode[STIX]{x1D716}^{n}\unicode[STIX]{x1D70C}^{n}\left[{\displaystyle \frac{\unicode[STIX]{x1D705}}{h_{s}}}\left(\dot{\unicode[STIX]{x1D719}}_{n-1}\sin \unicode[STIX]{x1D703}-(n-1)\unicode[STIX]{x1D719}_{n-1}\cos \unicode[STIX]{x1D703}\right)+{\displaystyle \frac{\unicode[STIX]{x1D719}_{n-2}^{\prime \prime }}{h_{s}^{2}}}+{\displaystyle \frac{(\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703})^{\prime }}{h_{s}^{3}}}\unicode[STIX]{x1D719}_{n-3}^{\prime }\right],\end{eqnarray}$$

where we defined $\unicode[STIX]{x1D719}_{-1}=0$. We then expand the inverse powers of $h_{s}=1-\unicode[STIX]{x1D716}\unicode[STIX]{x1D705}\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}$ in a power series in $\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}$ and collect terms of the same power in $\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}$, yielding for $n\geqslant 3$

(3.15)$$\begin{eqnarray}\displaystyle n^{2}\unicode[STIX]{x1D719}_{n}+\ddot{\unicode[STIX]{x1D719}}_{n} & = & \displaystyle -\mathop{\sum }_{m=2}^{n}\unicode[STIX]{x1D705}^{n-m}\cos \unicode[STIX]{x1D703}^{n-m}\left[\vphantom{\int }\unicode[STIX]{x1D705}\left(\dot{\unicode[STIX]{x1D719}}_{m-1}\sin \unicode[STIX]{x1D703}-(m-1)\unicode[STIX]{x1D719}_{m-1}\cos \unicode[STIX]{x1D703}\right)\right.\nonumber\\ \displaystyle & & \displaystyle \left.+(n-m+1)\unicode[STIX]{x1D719}_{m-2}^{\prime \prime }+{\displaystyle \frac{(n-m+1)(n-m+2)}{2}}(\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703})^{\prime }\unicode[STIX]{x1D719}_{m-3}^{\prime }\right]\!.\qquad\end{eqnarray}$$

We note that (3.15) has the form of a periodically driven harmonic oscillator with natural frequency $n$. Similarly to $\unicode[STIX]{x1D719}_{2}$ and $\unicode[STIX]{x1D719}_{3}$, we decompose $\unicode[STIX]{x1D719}_{n}$ into its Fourier harmonics as

(3.16)$$\begin{eqnarray}\unicode[STIX]{x1D719}_{n}=\mathop{\sum }_{p=0}^{n}\unicode[STIX]{x1D719}_{np}^{c}(s)\cos pu+\unicode[STIX]{x1D719}_{np}^{s}(s)\sin pu.\end{eqnarray}$$

At each order $n$, the expansion of (3.16) can be plugged in (3.15), yielding a set of one-dimensional differential equations for the coefficients $\unicode[STIX]{x1D719}_{np}^{c}(s)$ and $\unicode[STIX]{x1D719}_{np}^{s}(s)$. As the right-hand side of (3.15) only contains frequencies in $\unicode[STIX]{x1D714}$ up to $(n-2)$ (shown in appendix A), the solutions of $\unicode[STIX]{x1D719}_{np}^{c}$ and $\unicode[STIX]{x1D719}_{np}^{s}$ in (3.16) for $0\leqslant p\leqslant n-2$ are determined by the lower-order $O(\unicode[STIX]{x1D716}^{n-1}\unicode[STIX]{x1D70C}^{n-1})$ particular solutions. The two remaining functions $\unicode[STIX]{x1D719}_{nn}^{c}(s)$ and $\unicode[STIX]{x1D719}_{nn}^{s}(s)$ are then integration parameters from the solution of the homogeneous equation in (3.16). Finally, we remark that the analyticity condition of (2.14) for $\unicode[STIX]{x1D719}$ can be derived from the solution of Laplace’s equation, equation (3.3), as shown in appendix A.

3.2 Vacuum magnetic toroidal flux surface function

We now determine the expression for the normalized vacuum toroidal flux $\unicode[STIX]{x1D713}^{0}(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D714},s)$. Using (1.4) and (3.1), we determine $\unicode[STIX]{x1D713}^{0}$ via

(3.17)$$\begin{eqnarray}\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}^{0}=0.\end{eqnarray}$$

Expanding the gradient operator using (3.2), we rewrite (3.17) as

(3.18)$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}+{\displaystyle \frac{1}{\unicode[STIX]{x1D70C}^{2}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}+{\displaystyle \frac{\unicode[STIX]{x1D716}^{2}}{h_{s}^{2}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}s}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{0}}{\unicode[STIX]{x2202}s}}=0.\end{eqnarray}$$

An expression for the power series coefficients $\unicode[STIX]{x1D713}_{n}^{0}$ of $\unicode[STIX]{x1D713}^{0}$ can then be obtained by expanding both $\unicode[STIX]{x1D719}$ and $\unicode[STIX]{x1D713}^{0}$ in (3.18) in powers of $\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}$ using (2.13) and

(3.19)$$\begin{eqnarray}\unicode[STIX]{x1D713}^{0}=\mathop{\sum }_{n=0}^{\infty }\unicode[STIX]{x1D713}_{n}^{0}(\unicode[STIX]{x1D714},s)\unicode[STIX]{x1D716}^{n}\unicode[STIX]{x1D70C}^{n},\end{eqnarray}$$

yielding

(3.20)$$\begin{eqnarray}2n\unicode[STIX]{x1D719}_{2}\unicode[STIX]{x1D713}_{n}^{0}+\dot{\unicode[STIX]{x1D719}}_{2}\dot{\unicode[STIX]{x1D713}}_{n}^{0}+B_{0}\unicode[STIX]{x1D713}_{n}^{0^{\prime }}=F_{n}^{0},\end{eqnarray}$$

where $F_{0}^{0}=0$ and

(3.21)$$\begin{eqnarray}\displaystyle F_{n}^{0} & = & \displaystyle -\mathop{\sum }_{m=0}^{n-1}\left[\vphantom{\mathop{\sum }_{f=m}^{n}}(n+2-m)m\unicode[STIX]{x1D719}_{n+2-m}\unicode[STIX]{x1D713}_{m}^{0}+\dot{\unicode[STIX]{x1D719}}_{n+2-m}\dot{\unicode[STIX]{x1D713}}_{m}^{0}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\unicode[STIX]{x1D713}_{m}^{0^{\prime }}\mathop{\sum }_{f=m}^{n}(n-f+1)\unicode[STIX]{x1D705}^{n-f}\cos \unicode[STIX]{x1D703}^{n-f}\unicode[STIX]{x1D719}_{f-m}^{\prime }\right],\end{eqnarray}$$

for $n>0$. Although a formal solution of (3.20) can be obtained using the method of characteristics (as shown in appendix B), here, we focus on deriving a one-dimensional system of differential equations for $\unicode[STIX]{x1D713}_{n}^{0}$ with the coefficients $\unicode[STIX]{x1D719}_{np}^{c}$ and $\unicode[STIX]{x1D719}_{np}^{s}$ as sources. As an aside, we note that both $\unicode[STIX]{x1D713}^{0}$ and $\unicode[STIX]{x1D6FC}^{0}$ obey the constraint in (3.17). The distinction between the two is made by requiring $\unicode[STIX]{x1D713}^{0}$ to obey the analyticity condition, equation (2.14).

Using (3.20) and the analyticity condition, equation (2.14), the lowest-order solutions for $\unicode[STIX]{x1D713}^{0}$ are then given by $\unicode[STIX]{x1D713}_{0}^{0^{\prime }}=\unicode[STIX]{x1D713}_{1}^{0}=0$. The constant $\unicode[STIX]{x1D713}_{0}^{0}$ is set to zero by requiring that $\unicode[STIX]{x1D713}^{0}$ vanishes on the magnetic axis. Focusing on the $n=2$ case in (3.20), the equation for $\unicode[STIX]{x1D713}_{2}^{0}=\unicode[STIX]{x1D713}_{20}^{0}+\unicode[STIX]{x1D713}_{22}^{0s}\sin 2u+\unicode[STIX]{x1D713}_{22}^{0c}\cos 2u$ can be written as

(3.22)$$\begin{eqnarray}\unicode[STIX]{x1D6F9}_{2}^{0^{\prime }}=A_{2}^{0}\unicode[STIX]{x1D6F9}_{2}^{0},\end{eqnarray}$$

with $\unicode[STIX]{x1D6F9}_{2}^{0}=B_{0}(s)^{-1}[\unicode[STIX]{x1D713}_{20}^{0},\unicode[STIX]{x1D713}_{22}^{0s},\unicode[STIX]{x1D713}_{22}^{0c}]^{\text{T}}$ and $A_{2}^{0}$ given by

(3.23)$$\begin{eqnarray}A_{2}^{0}=\left(\begin{array}{@{}ccc@{}}0 & -2\unicode[STIX]{x1D707}u^{\prime } & \unicode[STIX]{x1D702}^{\prime }\\ -2\unicode[STIX]{x1D707}u^{\prime } & 0 & 2u^{\prime }\\ \unicode[STIX]{x1D702}^{\prime } & -2u^{\prime } & 0\\ \end{array}\right).\end{eqnarray}$$

The system of equations in (3.22) can be simplified by introducing the transformation

(3.24)$$\begin{eqnarray}\unicode[STIX]{x1D6F9}_{2}^{0}=T_{2}\unicode[STIX]{x1D70E}_{2}^{0},\end{eqnarray}$$

with $\unicode[STIX]{x1D70E}_{2}^{0}=[\unicode[STIX]{x1D70E}_{21}^{0},\unicode[STIX]{x1D70E}_{22}^{0},\unicode[STIX]{x1D70E}_{23}^{0}]^{\text{T}}$ and $\unicode[STIX]{x1D64F}_{2}$ the matrix

(3.25)$$\begin{eqnarray}\unicode[STIX]{x1D64F}_{2}=2\left(\begin{array}{@{}ccc@{}}-\text{sinh}\,\unicode[STIX]{x1D702} & \text{i} & \cosh \unicode[STIX]{x1D702}\\ -\text{sinh}\,\unicode[STIX]{x1D702} & -\text{i} & \cosh \unicode[STIX]{x1D702}\\ \cosh \unicode[STIX]{x1D702} & 0 & -\text{sinh}\,\unicode[STIX]{x1D702}\\ \end{array}\right).\end{eqnarray}$$

The quantities $\unicode[STIX]{x1D70E}_{2}^{0}$ then satisfy the following decoupled system of equations:

(3.26)$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{2}^{0^{\prime }}=\left(\begin{array}{@{}ccc@{}}2{\displaystyle \frac{\text{i}u^{\prime }}{\cosh \unicode[STIX]{x1D702}}} & 0 & 0\\ 0 & -2{\displaystyle \frac{\text{i}u^{\prime }}{\cosh \unicode[STIX]{x1D702}}} & 0\\ 0 & 0 & 0\end{array}\right)\unicode[STIX]{x1D70E}_{2}^{0}.\end{eqnarray}$$

The solution of (3.26) can then be given in terms of the integral

(3.27)$$\begin{eqnarray}v(s)=\int _{0}^{s}{\displaystyle \frac{u^{\prime }(x)}{\cosh \unicode[STIX]{x1D702}(x)}}\,\text{d}x=\int _{0}^{s}\sqrt{1-\unicode[STIX]{x1D707}(x)^{2}}[\unicode[STIX]{x1D6FF}^{\prime }(x)-\unicode[STIX]{x1D70F}(x)]\,\text{d}x,\end{eqnarray}$$

as

(3.28)$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{2}^{0}(s)=[\unicode[STIX]{x1D70E}_{21}^{00}\text{e}^{2\text{i}v(s)},\unicode[STIX]{x1D70E}_{22}^{00}\text{e}^{-2\text{i}v(s)},\unicode[STIX]{x1D70E}_{23}^{00}]^{\text{T}},\end{eqnarray}$$

with $\unicode[STIX]{x1D70E}_{21}^{00},\unicode[STIX]{x1D70E}_{22}^{00}$ and $\unicode[STIX]{x1D70E}_{23}^{00}$ constants. The flux surface function $\unicode[STIX]{x1D713}_{2}^{0}$ can then be found using (3.24).

A more streamlined method to obtain the lowest-order flux surface function $\unicode[STIX]{x1D713}_{2}^{0}$ can be found by noting that the free parameter $\unicode[STIX]{x1D6FF}$ in the analyticity condition in (2.14) can be used to set $\unicode[STIX]{x1D713}_{22}^{0s}=0$, i.e. the expansion coefficients of the field $\unicode[STIX]{x1D719}$ are chosen in such a way that the $\sin 2u$ terms in $\unicode[STIX]{x1D713}_{2}^{0}$ vanish. From (3.22), the remaining $\unicode[STIX]{x1D713}_{20}^{0}$ and $\unicode[STIX]{x1D713}_{22}^{0c}$ terms are then given by

(3.29)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D702}^{\prime }={\displaystyle \frac{\unicode[STIX]{x1D713}_{022}^{c^{\prime }}}{\unicode[STIX]{x1D713}_{20}^{0}}}={\displaystyle \frac{\unicode[STIX]{x1D713}_{20}^{0^{\prime }}}{\unicode[STIX]{x1D713}_{22}^{0c}}}, & \displaystyle\end{eqnarray}$$

(3.30)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D707}=\tanh \unicode[STIX]{x1D702}={\displaystyle \frac{\unicode[STIX]{x1D713}_{22}^{0c}}{\unicode[STIX]{x1D713}_{20}^{0}}}. & \displaystyle\end{eqnarray}$$

Defining $\unicode[STIX]{x1D713}_{22}^{0c}=\sinh \unicode[STIX]{x1D702}$ and $\unicode[STIX]{x1D713}_{20}^{0}=\cosh \unicode[STIX]{x1D702}$, the vacuum magnetic flux surface function $\unicode[STIX]{x1D713}^{0}=\unicode[STIX]{x1D713}_{2}^{0}\unicode[STIX]{x1D70C}^{2}\unicode[STIX]{x1D716}^{2}+O(\unicode[STIX]{x1D716}^{3})$ can then be written as

(3.31)$$\begin{eqnarray}\unicode[STIX]{x1D713}_{2}^{0}={\displaystyle \frac{B_{0}\unicode[STIX]{x03C0}}{\sqrt{1-\unicode[STIX]{x1D707}^{2}}}}\left(1+\unicode[STIX]{x1D707}\cos 2u\right)=B_{0}\unicode[STIX]{x03C0}\left(e^{\unicode[STIX]{x1D702}}\cos ^{2}u+\text{e}^{-\unicode[STIX]{x1D702}}\sin ^{2}u\right).\end{eqnarray}$$

The multiplicative constant in (3.31) is chosen such that $\unicode[STIX]{x1D713}$ equals the toroidal magnetic flux, i.e.

(3.32)$$\begin{eqnarray}\unicode[STIX]{x1D713}={\displaystyle \frac{1}{L}}\int (\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s)\,\text{d}V,\end{eqnarray}$$

with $\text{d}V$ the volume element, which in the $(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},s)$ coordinate system reads $\text{d}V=\text{d}\unicode[STIX]{x1D713}\,\text{d}\unicode[STIX]{x1D703}\,\text{d}s/(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}\boldsymbol{\cdot }\unicode[STIX]{x1D735}s)$. We note that for a circular cross-section with $\unicode[STIX]{x1D707}=0$, $\unicode[STIX]{x1D713}^{0}=\unicode[STIX]{x03C0}B_{0}\unicode[STIX]{x1D70C}^{2}$.

The analysis that led to the solution in (3.28) can be extended to higher orders. Explicit expressions for the third- and fourth-order solutions can be found in appendix C. Furthermore, as shown in appendix D, the system of equations for $\unicode[STIX]{x1D6F9}_{n}^{0}$, with $\unicode[STIX]{x1D6F9}_{n}^{0}$ the column vector with the coefficients of the Fourier expansion of $B_{0}^{-n/2}\unicode[STIX]{x1D713}_{n}^{0}$ in $u$ for arbitrary $n$ with $\unicode[STIX]{x1D713}_{n}^{0}$ expanded as

(3.33)$$\begin{eqnarray}\unicode[STIX]{x1D713}_{n}^{0}=\mathop{\sum }_{p=0}^{n}\unicode[STIX]{x1D713}_{np}^{0c}\cos pu+\unicode[STIX]{x1D713}_{np}^{0s}\sin pu,\end{eqnarray}$$

can be cast into the following form

(3.34)$$\begin{eqnarray}\unicode[STIX]{x1D6F9}_{n}^{0^{\prime }}=\unicode[STIX]{x1D63C}_{n}\unicode[STIX]{x1D6F9}_{n}^{0}+\unicode[STIX]{x1D63D}_{n}^{0},\end{eqnarray}$$

with $\unicode[STIX]{x1D63C}_{n}$ and $\unicode[STIX]{x1D63D}_{n}^{0}$ square matrices with periodic coefficients. An analysis of the properties of $\unicode[STIX]{x1D63C}_{n}$ and $\unicode[STIX]{x1D63D}_{n}^{0}$ using the methods of Floquet theory is left for a future study.

The form of (3.26), (C 6) and (C 10) suggests the existence of a matrix $\unicode[STIX]{x1D64F}_{n}$ at arbitrary order in $n$ such that the vectors $\unicode[STIX]{x1D748}_{n}^{0}$ defined by

(3.35)$$\begin{eqnarray}\unicode[STIX]{x1D713}_{n}^{0}=\unicode[STIX]{x1D64F}_{n}\unicode[STIX]{x1D748}_{n}^{0},\end{eqnarray}$$

satisfy a decoupled system of first-order differential equations. Assuming the existence of $\unicode[STIX]{x1D64F}_{n}$ for $n>4$, the decoupled system of equations for $\unicode[STIX]{x1D748}_{n}^{0}=\unicode[STIX]{x1D64F}_{n}^{-1}\unicode[STIX]{x1D713}_{n}^{0}$, a column vector with entries $\unicode[STIX]{x1D70E}_{nm}^{0}$ can be written for arbitrary $n$ in the following simplified form:

(3.36)$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{nm}^{0^{\prime }}(s)-\text{i}mv^{\prime }(s)\unicode[STIX]{x1D70E}_{nm}^{0}=F_{nm}^{0},\end{eqnarray}$$

with $m$ an odd (even) integer if $n$ is odd (even) and $-n\leqslant m\leqslant n$ and $v(s)$ given by (3.27). The problem of determining the parameters associated with the flux surface function is then reduced to the solution of (3.36). A general solution of (3.36) is given by

(3.37)$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{nm}^{0}(s)=\text{e}^{\text{i}mv(s)}\int _{s_{0}}^{s}F_{nm}^{0}(x)\text{e}^{-\text{i}mv(x)}\,\text{d}x.\end{eqnarray}$$

We now require $\unicode[STIX]{x1D713}^{0}$ (hence $\unicode[STIX]{x1D70E}_{nm}^{0}$) to be periodic on $s$ with period $L$, i.e. we impose the periodicity condition $\unicode[STIX]{x1D70E}_{nm}^{0}(s+L)=\unicode[STIX]{x1D70E}_{nm}^{0}(s)$, which sets the constant of integration $s_{0}$. The periodic solution of (3.36) is then given by (Mercier Reference Mercier1964; Solov’ev & Shafranov Reference Solov’ev and Shafranov1970)

(3.38)$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{nm}^{0}(s)={\displaystyle \frac{\text{e}^{\text{i}mv(s)}}{\text{e}^{-\text{i}mv(L)}-1}}\int _{s}^{s+L}\text{e}^{-\text{i}mv(x)}F_{nm}^{0}(x)\,\text{d}x.\end{eqnarray}$$

Analysis of (3.38) shows that a periodic solution of $\unicode[STIX]{x1D70E}_{nm}^{0}$ (hence $\unicode[STIX]{x1D713}_{n}^{0}$) yields a resonant denominator when

(3.39)$$\begin{eqnarray}v(L)=2\unicode[STIX]{x03C0}{\displaystyle \frac{l}{m}},\end{eqnarray}$$

with $l$ an integer. We therefore conclude that, for a solution of $\unicode[STIX]{x1D70E}_{nm}^{0}$ to exist, either the rotational transform on axis is an irrational number, or the numerators in (3.38) vanish for a rational on-axis rotational transform. In § 5, we relate the parameter $v(L)$ to the rotational transform $\unicode[STIX]{x1D704}$ and show that (3.39) is satisfied when the magnetic field lines in the vicinity of the magnetic axis close on themselves after one or more circuits along $s$, i.e. (3.39) is the condition for the existence of rational surfaces.

3.3 Vacuum field line label

The vacuum field line label $\unicode[STIX]{x1D6FC}^{0}$ is found by equating (1.4) and (3.1), yielding the following set of three coupled equations:

(3.40)$$\begin{eqnarray}\displaystyle & \displaystyle h_{s}\unicode[STIX]{x1D70C}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}^{0}}{\unicode[STIX]{x2202}s}}-{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{0}}{\unicode[STIX]{x2202}s}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}^{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}, & \displaystyle\end{eqnarray}$$

(3.41)$$\begin{eqnarray}\displaystyle & \displaystyle h_{s}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}=\unicode[STIX]{x1D70C}\left({\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{0}}{\unicode[STIX]{x2202}s}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}^{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}-{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}^{0}}{\unicode[STIX]{x2202}s}}\right), & \displaystyle\end{eqnarray}$$

(3.42)$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D70C}^{2}}{h_{s}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D719}}{\unicode[STIX]{x2202}s}}=\unicode[STIX]{x1D70C}\left({\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}^{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}-{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}^{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}\right). & \displaystyle\end{eqnarray}$$

Expanding $\unicode[STIX]{x1D6FC}^{0}$ in $\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}$ as

(3.43)$$\begin{eqnarray}\unicode[STIX]{x1D6FC}^{0}=\mathop{\sum }_{n}\unicode[STIX]{x1D6FC}_{n}^{0}(\unicode[STIX]{x1D714},s)\unicode[STIX]{x1D716}^{n}\unicode[STIX]{x1D70C}^{n},\end{eqnarray}$$

we find that, to lowest order in $\unicode[STIX]{x1D716}$, equations (3.40)–(3.42) reduce to

(3.44)$$\begin{eqnarray}\displaystyle & \displaystyle 2\unicode[STIX]{x1D719}_{2}=\dot{\unicode[STIX]{x1D713}}_{2}^{0}\unicode[STIX]{x1D6FC}_{0}^{0^{\prime }}-\unicode[STIX]{x1D713}_{2}^{0^{\prime }}\dot{\unicode[STIX]{x1D6FC}}^{0}, & \displaystyle\end{eqnarray}$$

(3.45)$$\begin{eqnarray}\displaystyle & \displaystyle \dot{\unicode[STIX]{x1D719}}_{2}=-2\unicode[STIX]{x1D713}_{2}^{0}\unicode[STIX]{x1D6FC}_{0}^{0^{\prime }}, & \displaystyle\end{eqnarray}$$

(3.46)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D719}_{0}^{\prime }=2\unicode[STIX]{x1D713}_{2}^{0}\dot{\unicode[STIX]{x1D6FC}}_{0}^{0}. & \displaystyle\end{eqnarray}$$

We note that, by eliminating $\unicode[STIX]{x1D6FC}_{0}^{0^{\prime }}$ and $\dot{\unicode[STIX]{x1D6FC}}_{0}^{0}$ in (3.44), using (3.45) and (3.46) we obtain (3.20) with $n=2$. Solving (3.46) for $\unicode[STIX]{x1D6FC}_{0}^{0}$ and plugging the result in (3.45), we find that

(3.47)$$\begin{eqnarray}\unicode[STIX]{x1D6FC}^{0}={\displaystyle \frac{1}{2\unicode[STIX]{x03C0}}}\left[\arctan \left(\text{e}^{-\unicode[STIX]{x1D702}}\tan u\right)-v(s)\right]+O(\unicode[STIX]{x1D716}).\end{eqnarray}$$

As expected, in contrast with $\unicode[STIX]{x1D713}^{0}$, $\unicode[STIX]{x1D6FC}^{0}$ does not obey the analyticity condition in (2.14).

In order to obtain the vacuum field line label $\unicode[STIX]{x1D6FC}^{0}$ to arbitrary order, we expand $\unicode[STIX]{x1D6FC}^{0}$ in powers of $\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}$ and obtain the following formulas for (3.41) and (3.42) and $n>0$

(3.48)$$\begin{eqnarray}\displaystyle \hspace{-24.0pt}-2(\unicode[STIX]{x1D713}_{2}^{0})^{n/2+1}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}}\left(\unicode[STIX]{x1D6FC}_{n}^{0}(\unicode[STIX]{x1D713}_{2}^{0})^{-n/2}\right) & = & \displaystyle \dot{\unicode[STIX]{x1D719}}_{n+2}-\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703}\dot{\unicode[STIX]{x1D719}}_{n+1}\nonumber\\ \displaystyle \hspace{-24.0pt} & & \displaystyle -\,\mathop{\sum }_{m=0}^{n-1}[m\unicode[STIX]{x1D6FC}_{m}^{0}\unicode[STIX]{x1D713}_{n+2-m}^{^{\prime }0}-(n+2-m)\unicode[STIX]{x1D6FC}_{m}^{^{\prime }0}\unicode[STIX]{x1D713}_{n+2-m}^{0}],\end{eqnarray}$$

(3.49)$$\begin{eqnarray}\displaystyle \hspace{-24.0pt}2(\unicode[STIX]{x1D713}_{2}^{0})^{n/2+1}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}\left(\unicode[STIX]{x1D6FC}_{n}^{0}(\unicode[STIX]{x1D713}_{2}^{0})^{-n/2}\right) & = & \displaystyle \mathop{\sum }_{m=0}^{n}\unicode[STIX]{x1D719}_{m}^{\prime }(\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703})^{n-m}\nonumber\\ \displaystyle \hspace{-24.0pt} & & \displaystyle +\,\mathop{\sum }_{m=0}^{n-1}[m\unicode[STIX]{x1D6FC}_{m}^{0}\dot{\unicode[STIX]{x1D713}}_{n+2-m}-(n+2-m)\dot{\unicode[STIX]{x1D6FC}}_{m}^{0}\unicode[STIX]{x1D713}_{n+2-m}^{0}].\end{eqnarray}$$

As done for $\unicode[STIX]{x1D6FC}_{0}^{0}$, at each order, equation (3.49) can be used to find an analytical expression for $\unicode[STIX]{x1D6FC}_{n}^{0}$ up to an additive function of $s$, which is set by (3.48).

4.1 MHD current density vector

We start by deriving the forms of the perpendicular and parallel current densities, $\boldsymbol{J}_{\bot }$ and $J_{\Vert }$ respectively. In the following, to simplify the notation, we define the coefficients of the inverse expansion for $|\boldsymbol{B}^{0}|^{2}$ as

(4.7)$$\begin{eqnarray}{\displaystyle \frac{1}{|\boldsymbol{B}^{0}|^{2}}}={\displaystyle \frac{1}{|\unicode[STIX]{x1D735}\unicode[STIX]{x1D719}|^{2}}}=\mathop{\sum }_{n=0}^{\infty }{\displaystyle \frac{{\mathcal{B}}_{n}}{B_{0}(s)^{2}}}\unicode[STIX]{x1D716}^{n}\unicode[STIX]{x1D70C}^{n},\end{eqnarray}$$

with $B_{0}(s)=|\boldsymbol{B}^{0}|_{\unicode[STIX]{x1D70C}=0}$ the vacuum magnetic field modulus on axis and ${\mathcal{B}}_{0}$=1.

Using (4.5), the three components $(J_{\bot s},J_{\bot \unicode[STIX]{x1D70C}},J_{\bot \unicode[STIX]{x1D714}})$ of the perpendicular current can be written as

(4.8)$$\begin{eqnarray}\displaystyle \hspace{-26.39996pt} & \displaystyle B_{0}(s)^{2}J_{\bot sl+2}=\mathop{\sum }_{g=0}^{l}\mathop{\sum }_{p=0}^{g}{\mathcal{B}}_{l-g}\left[(p+2)\unicode[STIX]{x1D719}_{p+2}\dot{\unicode[STIX]{x1D713}}_{0g-p+2}-(g-p+2)\unicode[STIX]{x1D713}_{0g-p+2}\dot{\unicode[STIX]{x1D719}}_{p+2}\right]\!, & \displaystyle\end{eqnarray}$$

(4.9)$$\begin{eqnarray}\displaystyle \hspace{-26.39996pt} & \displaystyle B_{0}(s)^{2}J_{\bot \unicode[STIX]{x1D70C}l+1}=\mathop{\sum }_{f=0}^{l}\mathop{\sum }_{g=0}^{f}\mathop{\sum }_{p=0}^{g}{\mathcal{B}}_{l-g}\left[\dot{\unicode[STIX]{x1D719}}_{p}\unicode[STIX]{x1D713}_{0g-p+2}^{\prime }-\unicode[STIX]{x1D719}_{p}^{\prime }\dot{\unicode[STIX]{x1D713}}_{0g-p+2}\right](\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703})^{l-f}\!, & \displaystyle\end{eqnarray}$$

(4.10)$$\begin{eqnarray}\displaystyle \hspace{-26.39996pt} & \displaystyle B_{0}(s)^{2}J_{\bot \unicode[STIX]{x1D714}l+1}=\mathop{\sum }_{f=0}^{l}\mathop{\sum }_{g=0}^{f}\mathop{\sum }_{p=0}^{g}{\mathcal{B}}_{l-g}\!\left[(g-p+2)\unicode[STIX]{x1D713}_{0g-p+2}\unicode[STIX]{x1D719}_{p}^{\prime }-p\unicode[STIX]{x1D713}_{0g-p+2}^{\prime }\unicode[STIX]{x1D719}_{p}\right]\!(\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703})^{l-f}\!. & \displaystyle\end{eqnarray}$$

We remark that, as expected, the perpendicular current vanishes on axis, i.e. $J_{\bot s0}=J_{\bot \unicode[STIX]{x1D70C}0}=J_{\bot \unicode[STIX]{x1D714}0}=0$.

The power series expansion coefficients $J_{\Vert n}$ of $J_{\Vert }=\sum _{n}J_{\Vert n}\unicode[STIX]{x1D716}^{n}\unicode[STIX]{x1D70C}^{n}$ are obtained by expanding the terms included in (4.6) in power of $\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}$, yielding

(4.11)$$\begin{eqnarray}\displaystyle 2n\unicode[STIX]{x1D719}_{2}J_{\Vert n}+\dot{\unicode[STIX]{x1D719}}_{2}\dot{J}_{\Vert n}+B_{0}J_{\Vert n}^{\prime }={\displaystyle \frac{G_{n}}{B_{0}(s)^{2}}}-H_{n}, & & \displaystyle\end{eqnarray}$$

with $H_{0}=G_{0}=0$ and the coefficients $H_{n}$ and $G_{n}$ given by

(4.12)$$\begin{eqnarray}H_{n}=\mathop{\sum }_{g=0}^{n-1}\left[(n-g+2)\unicode[STIX]{x1D719}_{n-g+2}gJ_{\Vert g}+\dot{\unicode[STIX]{x1D719}}_{n-g+2}\dot{J}_{\Vert g}+\unicode[STIX]{x1D719}_{n-g}^{\prime }J_{\Vert g}^{\prime }+(\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703})^{n-g}\mathop{\sum }_{m=0}^{g}\unicode[STIX]{x1D719}_{g-m}^{\prime }J_{\Vert m}^{\prime }\right]\!,\end{eqnarray}$$

and

(4.13)$$\begin{eqnarray}\displaystyle G_{n} & = & \displaystyle \mathop{\sum }_{g=1}^{n}\left\{\mathop{\sum }_{p=0}^{g-1}\mathop{\sum }_{r=0}^{n-g}\unicode[STIX]{x1D719}_{r}^{\prime }(\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703})^{n-g-r}\left[(g-p)\dot{\unicode[STIX]{x1D713}}_{0p+2}{\mathcal{B}}_{g-p}-(p+2)\dot{\unicode[STIX]{x1D713}}_{0p+2}\dot{{\mathcal{B}}}_{g-p}\right]\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\mathop{\sum }_{l=1}^{g}\mathop{\sum }_{p=0}^{l-1}(\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703})^{g-l}\left[\left({\mathcal{B}}_{l-p-1}^{\prime }-2{\displaystyle \frac{B_{0}^{\prime }}{B_{0}}}{\mathcal{B}}_{l-p-1}\right)(p+2)\unicode[STIX]{x1D713}_{0p+2}\dot{\unicode[STIX]{x1D719}}_{n-g+1}\right.\right.\nonumber\\ \displaystyle & & \displaystyle \left.\left.+\,{\mathcal{B}}_{l-p}\left(\unicode[STIX]{x1D713}_{0p+1}^{\prime }\unicode[STIX]{x1D719}_{n-g+1}-(l-p)\unicode[STIX]{x1D713}_{0p+1}\dot{\unicode[STIX]{x1D719}}_{n-g+1}\right)\vphantom{\mathop{\sum }_{g=1}^{n}}\right]\right\},\end{eqnarray}$$

respectively, for $n>0$. We note that, as both $J_{\Vert }$ and $\unicode[STIX]{x1D713}^{0}$ are obtained by solving a magnetic differential equation, $J_{\Vert n}$ obeys an advection equation similar to the one of $\unicode[STIX]{x1D713}_{2}^{0}$, equation (B 1). Therefore, similarly to (3.34), the components of the parallel current column vector ${\mathcal{J}}_{\Vert n}$ with ${\mathcal{J}}_{\Vert n}=[J_{\Vert n0}J_{\Vert n2}^{c}J_{\Vert n2}^{s}~\cdots ]^{\text{T}}$ for even $n$ and ${\mathcal{J}}_{\Vert n}=[J_{\Vert n1}^{c}J_{\Vert n1}^{s}J_{\Vert n3}^{c}J_{\Vert n3}^{s}~\cdots \,]^{\text{T}}$ for odd $n$, can be shown to obey

(4.14)$$\begin{eqnarray}{\mathcal{J}}_{\Vert n}^{\prime }=A_{n}^{0}{\mathcal{J}}_{\Vert n}+C_{n}^{0},\end{eqnarray}$$

with $C_{n}^{0}$ a source term dependent on the components $G_{n}$ and $H_{n}$. Using the same transformation matrix $\unicode[STIX]{x1D64F}_{n}$ as in (3.35), the components $\unicode[STIX]{x1D70E}_{nm}^{1}$ of $\unicode[STIX]{x1D70E}_{n}^{1}$, where $\unicode[STIX]{x1D70E}_{n}^{1}$ is given by

(4.15)$$\begin{eqnarray}{\mathcal{J}}_{\Vert n}=\unicode[STIX]{x1D64F}_{n}\unicode[STIX]{x1D70E}_{n}^{1},\end{eqnarray}$$

can be shown to satisfy (3.36) with a different source term, namely

(4.16)$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{nm}^{1^{\prime }}(s)-\text{i}mv^{\prime }\unicode[STIX]{x1D70E}_{nm}^{1}=D_{nm}^{1},\end{eqnarray}$$

with $D_{nm}^{1}=\unicode[STIX]{x1D64F}_{n}^{-1}C_{n}^{0}$. The solution of $\unicode[STIX]{x1D70E}_{nm}^{1}$ is then of the form of (3.38).

For $f=0$, equation (4.11) determines the current on axis $J_{\Vert 0}$ to be a constant. For $f=1$, the equation for the coefficients $\unicode[STIX]{x1D706}$ and $\unicode[STIX]{x1D707}$ of the first-order current $J_{\Vert 1}=\unicode[STIX]{x1D706}\sin u+\unicode[STIX]{x1D707}\cos u$ can be written as

(4.17)$$\begin{eqnarray}\unicode[STIX]{x1D70E}_{1}^{\prime }-\text{i}v^{\prime }(s)\unicode[STIX]{x1D70E}_{1}=-B_{0}^{3/4}(s)4\unicode[STIX]{x1D705}\left[\text{e}^{\unicode[STIX]{x1D702}/2}(1-\unicode[STIX]{x1D707}\cos \unicode[STIX]{x1D6FF})+\text{i}\text{e}^{-\unicode[STIX]{x1D702}/2}(1+\unicode[STIX]{x1D707})\sin \unicode[STIX]{x1D6FF}\right]\!,\end{eqnarray}$$

with $\unicode[STIX]{x1D70E}_{1}=B_{0}(s)^{-1/2}(\unicode[STIX]{x1D706}\text{e}^{\unicode[STIX]{x1D702}/2}-\text{i}\unicode[STIX]{x1D707}\text{e}^{-\unicode[STIX]{x1D702}/2})$.

4.2 MHD magnetic field vector

We now proceed by calculating the first-order magnetic field by expanding the components of $\boldsymbol{B}^{1}$ in powers of $\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}$ as

(4.18)$$\begin{eqnarray}\boldsymbol{B}^{1}=\mathop{\sum }_{n}\boldsymbol{B}_{n}^{1}(s\unicode[STIX]{x1D714})\unicode[STIX]{x1D716}^{n}\unicode[STIX]{x1D70C}^{n}\end{eqnarray}$$

and use (4.2) to solve for $\boldsymbol{B}_{\unicode[STIX]{x1D70C}n}^{1},\boldsymbol{B}_{\unicode[STIX]{x1D714}n}^{1}$ and $\boldsymbol{B}_{sn}^{1}$. From the $\unicode[STIX]{x1D70C}$ component of (4.2) we find at lowest order the constraint constraint ${\dot{B}}_{s0}^{1}=0$. For higher order, we obtain for the $\unicode[STIX]{x1D70C}$ component

(4.19)$$\begin{eqnarray}\mathop{\sum }_{n=1}^{\infty }\unicode[STIX]{x1D716}^{n}\unicode[STIX]{x1D70C}^{n}\left[{\displaystyle \frac{\unicode[STIX]{x2202}B_{sn}^{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}-\unicode[STIX]{x1D705}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}(\cos \unicode[STIX]{x1D703}B_{sn-1}^{1})-{\displaystyle \frac{\unicode[STIX]{x2202}B_{\unicode[STIX]{x1D714}n-1}^{1}}{\unicode[STIX]{x2202}s}}\right]=A_{\unicode[STIX]{x1D70C}}^{J},\end{eqnarray}$$

with $A_{\unicode[STIX]{x1D70C}}^{j}=\unicode[STIX]{x1D6FD}h_{s}p^{\prime }(\unicode[STIX]{x1D713}^{0})(\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}J_{\bot \unicode[STIX]{x1D70C}}+J_{\Vert }\unicode[STIX]{x1D70C}\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}\unicode[STIX]{x1D70C})$, for the $\unicode[STIX]{x1D714}$ component

(4.20)$$\begin{eqnarray}\mathop{\sum }_{n=1}^{\infty }\unicode[STIX]{x1D716}^{n}\unicode[STIX]{x1D70C}^{n}\left[nB_{sn}^{1}-(n+1)\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703}B_{sn-1}^{1}-{\displaystyle \frac{\unicode[STIX]{x2202}B_{\unicode[STIX]{x1D70C}n-1}^{1}}{\unicode[STIX]{x2202}s}}\right]=A_{\unicode[STIX]{x1D714}}^{J},\end{eqnarray}$$

with $A_{\unicode[STIX]{x1D714}}^{J}=-\unicode[STIX]{x1D6FD}h_{s}p^{\prime }(\unicode[STIX]{x1D713}^{0})(\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}J_{\bot \unicode[STIX]{x1D70C}}+J_{\Vert }\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}\unicode[STIX]{x1D714})$ and for the $s$ component

(4.21)$$\begin{eqnarray}\mathop{\sum }_{n=0}^{\infty }\unicode[STIX]{x1D716}^{n}\unicode[STIX]{x1D70C}^{n}\left[(n+1)B_{\unicode[STIX]{x1D714}n}^{1}-{\displaystyle \frac{\unicode[STIX]{x2202}B_{\unicode[STIX]{x1D70C}n}^{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}\right]=A_{s}^{J},\end{eqnarray}$$

with $A_{s}^{J}=\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}\unicode[STIX]{x1D6FD}p^{\prime }(\unicode[STIX]{x1D713}^{0})(J_{\bot s}+h_{s}^{-1}J_{\Vert }\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}s)$. In order to simplify the calculations and eliminate one of the three components of (4.2), we replace (4.19) by the condition, $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}^{1}=0$, which can be written as

(4.22)$$\begin{eqnarray}{\displaystyle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}B_{\unicode[STIX]{x1D70C}}^{1})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}+{\displaystyle \frac{\unicode[STIX]{x2202}B_{\unicode[STIX]{x1D714}}^{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}=-\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}\left[{\displaystyle \frac{\unicode[STIX]{x2202}B_{s}^{1}}{\unicode[STIX]{x2202}s}}+\unicode[STIX]{x1D705}{\displaystyle \frac{\unicode[STIX]{x2202}(\cos \unicode[STIX]{x1D703}B_{\unicode[STIX]{x1D714}}^{1})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}+{\displaystyle \frac{\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703}}{\unicode[STIX]{x1D70C}}}{\displaystyle \frac{\unicode[STIX]{x2202}(\unicode[STIX]{x1D70C}^{2}B_{s}^{1})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}\right].\end{eqnarray}$$

Expanding the functions $A_{\unicode[STIX]{x1D714}}^{J}$ and $A_{s}^{J}$ in terms of powers of $\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}$, the following expressions for the perturbed magnetic field are found

(4.23)$$\begin{eqnarray}\displaystyle (n+1)^{2}B_{\unicode[STIX]{x1D70C}n}^{1}+{\displaystyle \frac{\unicode[STIX]{x2202}^{2}B_{\unicode[STIX]{x1D70C}n}^{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}^{2}}} & = & \displaystyle -(n+1)\left[{\displaystyle \frac{1}{n+1}}{\displaystyle \frac{\unicode[STIX]{x2202}A_{sn}^{J}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}+{\displaystyle \frac{\unicode[STIX]{x2202}B_{sn-1}^{1}}{\unicode[STIX]{x2202}s}}\right.\nonumber\\ \displaystyle & & \displaystyle \left.+\,\unicode[STIX]{x1D705}{\displaystyle \frac{\unicode[STIX]{x2202}(\cos \unicode[STIX]{x1D703}B_{\unicode[STIX]{x1D714}n-1}^{1})}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}+(n+2)\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703}B_{sn-1}^{1}\right],\end{eqnarray}$$

(4.24)$$\begin{eqnarray}\displaystyle B_{\unicode[STIX]{x1D714}n}^{1} & = & \displaystyle {\displaystyle \frac{1}{n+1}}\left(A_{sn}^{J}+{\displaystyle \frac{\unicode[STIX]{x2202}B_{1\unicode[STIX]{x1D70C}n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}\right),\end{eqnarray}$$

(4.25)$$\begin{eqnarray}\displaystyle B_{sn}^{1} & = & \displaystyle {\displaystyle \frac{1}{n}}\left[{\displaystyle \frac{\unicode[STIX]{x2202}B_{\unicode[STIX]{x1D70C}n-1}^{1}}{\unicode[STIX]{x2202}s}}+(n+1)\unicode[STIX]{x1D705}\cos \unicode[STIX]{x1D703}B_{sn-1}^{1}+A_{\unicode[STIX]{x1D714}n}^{J}\right].\end{eqnarray}$$

The lowest-order field $\boldsymbol{B}_{0}^{1}$ is found to satisfy $B_{\unicode[STIX]{x1D70C}0}^{1}=B_{\unicode[STIX]{x1D70C}0}^{1c}(s)\cos u+B_{\unicode[STIX]{x1D70C}0}^{1s}(s)\sin u$, $B_{\unicode[STIX]{x1D714}0}^{1}=\unicode[STIX]{x2202}B_{\unicode[STIX]{x1D70C}0}^{1}/\unicode[STIX]{x2202}\unicode[STIX]{x1D714}$ and $B_{s0}^{1}=B_{s0}^{1}(s)$, in agreement with Solov’ev & Shafranov (Reference Solov’ev and Shafranov1970). While analyticity could be proven rigorously for the vacuum case, we note that the right-hand side of the forced linear harmonic oscillator equation for $B_{\unicode[STIX]{x1D70C}n}^{1}$ in (4.23) might contain $n+1$ resonating frequencies resulting from the product $J_{\Vert }\unicode[STIX]{x2202}\unicode[STIX]{x1D719}/\unicode[STIX]{x2202}s$ in the $A_{s}^{J}$ term. These resonances can lead to the appearance of non-analytic and weakly singular terms of the form $\unicode[STIX]{x1D70C}^{n}(\log \unicode[STIX]{x1D70C})^{m}$ as discussed in Weitzner (Reference Weitzner2016).

4.3 MHD flux surface function

We now obtain $\unicode[STIX]{x1D713}^{1}$ using the linearized flux surface condition, equation (4.3). We take advantage of the fact that (4.3) is similar to (3.17), although with a non-zero source term, i.e.

(4.26)$$\begin{eqnarray}2n\unicode[STIX]{x1D719}_{2}\unicode[STIX]{x1D713}_{n}^{1}+\dot{\unicode[STIX]{x1D719}}_{2}\dot{\unicode[STIX]{x1D713}}_{n}^{1}+B_{0}\unicode[STIX]{x1D713}_{n}^{1}=F_{n}^{1},\end{eqnarray}$$

with

(4.27)$$\begin{eqnarray}F_{n}^{1}=F_{n}^{0}-{\displaystyle \frac{\unicode[STIX]{x1D716}}{n!}}\left[{\displaystyle \frac{1}{\unicode[STIX]{x1D716}^{2}}}{\displaystyle \frac{\unicode[STIX]{x2202}^{n}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}^{n}}}\left(B_{\unicode[STIX]{x1D70C}}^{1}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}+{\displaystyle \frac{B_{\unicode[STIX]{x1D714}}^{1}}{\unicode[STIX]{x1D70C}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{0}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}+\unicode[STIX]{x1D716}{\displaystyle \frac{B_{s}^{1}}{h_{s}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{0}}{\unicode[STIX]{x2202}s}}\right)\right]_{\unicode[STIX]{x1D70C}=0}.\end{eqnarray}$$

By plugging the analyticity condition for $\unicode[STIX]{x1D713}^{1}$, equation (2.14), in (4.26) and defining the column vector $\unicode[STIX]{x1D6F9}_{n}^{1}$ in an analogous manner with $\unicode[STIX]{x1D6F9}_{n}^{0}$, the following set of equations for $\unicode[STIX]{x1D713}_{n}^{1}$ is found

(4.28)$$\begin{eqnarray}\unicode[STIX]{x1D6F9}_{n}^{1^{\prime }}=\unicode[STIX]{x1D63C}_{n}\unicode[STIX]{x1D6F9}_{n}^{1}+\unicode[STIX]{x1D63D}_{n}^{1},\end{eqnarray}$$

with $\unicode[STIX]{x1D63D}_{n}^{1}$ the matrix of Fourier coefficients of $F_{n}^{1}$ at each order $n$. The components of $\unicode[STIX]{x1D63C}_{n}$ are given in appendix D.

The lowest-order solution for $\unicode[STIX]{x1D713}^{1}$ is given by

(4.29)$$\begin{eqnarray}\unicode[STIX]{x1D713}^{1}=\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}[\unicode[STIX]{x1D713}_{11}^{c}(s)\cos u+\unicode[STIX]{x1D713}_{11}^{c}(s)\sin u]+O(\unicode[STIX]{x1D716}^{2}),\end{eqnarray}$$

with $\unicode[STIX]{x1D713}_{11}^{c}$ and $\unicode[STIX]{x1D713}_{11}^{s}$ obeying the set of equations

(4.30)$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{B_{0}(s)}{2}}\unicode[STIX]{x1D713}_{11}^{c^{\prime }}=\unicode[STIX]{x1D719}_{22}^{s}\unicode[STIX]{x1D713}_{11}^{s}+(\unicode[STIX]{x1D719}_{20}+\unicode[STIX]{x1D719}_{22}^{c})\unicode[STIX]{x1D713}_{11}^{c}-B_{1\unicode[STIX]{x1D70C}}^{c}(\unicode[STIX]{x1D713}_{20}^{0}+\unicode[STIX]{x1D713}_{22}^{0c})-B_{1\unicode[STIX]{x1D70C}}^{s}\unicode[STIX]{x1D713}_{22}^{0s}, & \displaystyle\end{eqnarray}$$

(4.31)$$\begin{eqnarray}\displaystyle & \displaystyle {\displaystyle \frac{B_{0}(s)}{2}}\unicode[STIX]{x1D713}_{11}^{s^{\prime }}=\unicode[STIX]{x1D719}_{22}^{s}\unicode[STIX]{x1D713}_{11}^{c}+(\unicode[STIX]{x1D719}_{20}-\unicode[STIX]{x1D719}_{22}^{c})\unicode[STIX]{x1D713}_{11}^{s}-B_{1\unicode[STIX]{x1D70C}}^{s}(\unicode[STIX]{x1D713}_{20}^{0}-\unicode[STIX]{x1D713}_{22}^{0c})-B_{1\unicode[STIX]{x1D70C}}^{c}\unicode[STIX]{x1D713}_{22}^{0s}. & \displaystyle\end{eqnarray}$$

4.4 Shafranov shift

Here, we show how to obtain the position of the magnetic axis $(x_{M},y_{M})$ once finite $\unicode[STIX]{x1D6FD}$ effects are taken into account. Although the procedure is valid for arbitrary order in the expansion parameter $\unicode[STIX]{x1D716}$, we calculate $(x_{M},y_{M})$ explicitly to lowest order in $\unicode[STIX]{x1D716}$ only. We first rewrite $\unicode[STIX]{x1D713}^{1}$ in Cartesian $(x,y)$ coordinates as

(4.32)$$\begin{eqnarray}\unicode[STIX]{x1D713}^{1}=\unicode[STIX]{x1D716}(ax+by),\end{eqnarray}$$

with the functions $a$ and $b$ given by $a=\unicode[STIX]{x1D713}_{11}^{c}\cos \unicode[STIX]{x1D6FF}+\unicode[STIX]{x1D713}_{11}^{s}\sin \unicode[STIX]{x1D6FF}$ and $b=-\unicode[STIX]{x1D713}_{11}^{c}\sin \unicode[STIX]{x1D6FF}+\unicode[STIX]{x1D713}_{11}^{s}\cos \unicode[STIX]{x1D6FF}$. The lowest-order magnetic flux surface function in $\unicode[STIX]{x1D716}$ and $\unicode[STIX]{x1D6FD}$, in $x$ and $y$ coordinates, is then given by

(4.33)$$\begin{eqnarray}\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}^{0}+\unicode[STIX]{x1D713}^{1}=ax+by+Ax^{2}+By^{2}+Cxy,\end{eqnarray}$$

with $A=(1+\unicode[STIX]{x1D707}\cos 2\unicode[STIX]{x1D6FF})/\sqrt{1-\unicode[STIX]{x1D707}^{2}}$, $B=(1-\unicode[STIX]{x1D707}\cos 2\unicode[STIX]{x1D6FF})/\sqrt{1-\unicode[STIX]{x1D707}^{2}}$ and $C=-2\unicode[STIX]{x1D707}\sin 2\unicode[STIX]{x1D6FF}/\sqrt{1-\unicode[STIX]{x1D707}^{2}}$. Setting the derivatives of $\unicode[STIX]{x1D713}$ with respect to $x$ and $y$ equal to zero, the position of the magnetic axis $(x_{M},y_{M})$ is found to be

(4.34)$$\begin{eqnarray}x_{M}={\displaystyle \frac{2aB-bC}{C^{2}-4AB}},\end{eqnarray}$$

and

(4.35)$$\begin{eqnarray}y_{M}={\displaystyle \frac{2Ab-aC}{C^{2}+4AB}}.\end{eqnarray}$$

The condition that the distortion of the magnetic surfaces be small, i.e. $x_{M}\sim y_{M}\ll a$ leads to a limit on the maximum allowed $\unicode[STIX]{x1D6FD}$. For the derivation of the plasma beta-limit $\unicode[STIX]{x1D6FD}\ll 2\unicode[STIX]{x1D716}\unicode[STIX]{x1D704}_{0}^{2}$ for the particular case of a circular magnetic axis, $\unicode[STIX]{x1D6FF}=n\unicode[STIX]{x03C0}$ and neglecting curvature effects, see Solov’ev & Shafranov (Reference Solov’ev and Shafranov1970).

4.5 MHD field line label

Finally, for the field line label $\unicode[STIX]{x1D6FC}^{1}$, using $\boldsymbol{B}^{1}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}^{1}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC}^{1}$, we derive the following set of three coupled equations:

(4.36)$$\begin{eqnarray}\displaystyle & \displaystyle h_{s}\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}B_{1\unicode[STIX]{x1D70C}}={\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}^{1}}{\unicode[STIX]{x2202}s}}-{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{1}}{\unicode[STIX]{x2202}s}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}^{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}, & \displaystyle\end{eqnarray}$$

(4.37)$$\begin{eqnarray}\displaystyle & \displaystyle h_{s}\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}B_{\unicode[STIX]{x1D714}}^{1}=\unicode[STIX]{x1D70C}\left({\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{1}}{\unicode[STIX]{x2202}s}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}^{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}-{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}^{1}}{\unicode[STIX]{x2202}s}}\right), & \displaystyle\end{eqnarray}$$

(4.38)$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D716}^{2}\unicode[STIX]{x1D70C}^{2}B_{s}^{1}=\unicode[STIX]{x1D70C}\left({\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}^{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}-{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D713}^{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}{\displaystyle \frac{\unicode[STIX]{x2202}\unicode[STIX]{x1D6FC}^{1}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D70C}}}\right). & \displaystyle\end{eqnarray}$$

Expanding $\unicode[STIX]{x1D6FC}^{1}$ in a series of powers of $\unicode[STIX]{x1D716}\unicode[STIX]{x1D70C}$ and following the approach in § 3.3, we find

(4.39)$$\begin{eqnarray}\displaystyle -(\unicode[STIX]{x1D713}_{1}^{1})^{n+1}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}s}}\left[\unicode[STIX]{x1D6FC}_{n}^{1}(\unicode[STIX]{x1D713}_{1}^{1})^{-n}\right] & = & \displaystyle B_{\unicode[STIX]{x1D714}n}^{1}-k\cos \unicode[STIX]{x1D703}B_{\unicode[STIX]{x1D714}n-1}^{1}\nonumber\\ \displaystyle & & \displaystyle -\,\mathop{\sum }_{m=0}^{n-1}\left[\unicode[STIX]{x1D713}_{n-m+1}^{1^{\prime }}m\unicode[STIX]{x1D6FC}_{m}^{1}-\unicode[STIX]{x1D713}_{n-m+1}^{1}(n-m+1)\unicode[STIX]{x1D6FC}_{m}^{1^{\prime }}\right]\end{eqnarray}$$

(4.40)$$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D713}_{1}^{1})^{n+1}{\displaystyle \frac{\unicode[STIX]{x2202}}{\unicode[STIX]{x2202}\unicode[STIX]{x1D714}}}\left[\unicode[STIX]{x1D6FC}_{n}^{1}(\unicode[STIX]{x1D713}_{1}^{1})^{-n}\right] & = & \displaystyle B_{1sn-1}\nonumber\\ \displaystyle & & \displaystyle +\,\mathop{\sum }_{m=0}^{n-1}\left[\dot{\unicode[STIX]{x1D713}}_{n-m+1}^{1}m\unicode[STIX]{x1D6FC}_{m}^{1}-\unicode[STIX]{x1D713}_{n-m+1}^{1}(n-m+1)\dot{\unicode[STIX]{x1D6FC}}_{m}^{1}\right]\!.\quad\end{eqnarray}$$

Analogously to (3.48) and (3.49), at each order in $\unicode[STIX]{x1D716}^{n}$, equation (4.40) can be used to find an analytical expression for $\unicode[STIX]{x1D6FC}_{n}^{1}$ up to an additive function of $s$, which is set by (3.49).