2.1 The Weitzner formulation of ideal MHS with scalar pressure

We begin with the discussion of the formalism developed in Weitzner (Reference Weitzner2014, Reference Weitzner2016) to study 3-D (non-symmetric) ideal MHS equilibria with scalar pressure. We assume throughout that the magnetic field $\boldsymbol {B}$ satisfies the ideal MHS equations

(2.1a--c)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{B} =0, \quad \boldsymbol{J}=\boldsymbol{\nabla} \times \boldsymbol{B} , \quad \boldsymbol{J}\times \boldsymbol{B} = \boldsymbol{\nabla} p, \end{equation}

where $\boldsymbol {J}$ represents the current, and $p$ is the scalar pressure. We have used units such that $\mu _0$ is set to unity. We shall assume that the magnetic field possess a set of nested toroidal flux surfaces, denoted by the flux label $\psi$, which are also isosurfaces of the pressure. Furthermore, we assume that $\boldsymbol {\nabla } \psi$ is non-zero everywhere except at the magnetic axis.

To describe the magnetic field vector, we now use the following contravariant and covariant representation:

(2.2a)\begin{gather} \boldsymbol{B} =\boldsymbol{\nabla} \psi \times \boldsymbol{\nabla} \alpha, \end{gather}

(2.2b)\begin{gather}\boldsymbol{B} = \boldsymbol{\nabla} \varPhi + K \boldsymbol{\nabla} \psi. \end{gather}

Equation (2.2a) is the standard Clebsch representation (D'haeseleer et al. Reference D'haeseleer, Hitchon, Callen and Shohet2012; Helander Reference Helander2014) of the magnetic field with $\psi$ and $\alpha$ denoting the toroidal magnetic flux and field line label, respectively. The Clebsch form manifestly satisfies $\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {B}=0$ and $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\psi =0$. The second representation of the magnetic field, (2.2b), is the so-called Boozer–Grad representation (see chapter 7.2 of D'haeseleer et al. (Reference D'haeseleer, Hitchon, Callen and Shohet2012)). We call $\varPhi$ the magnetic scalar potential and $K$ the current potential, since, in the limit $K\to 0$, or in the case of $K= K(\psi )$, the curl of $\boldsymbol {B}$ (2.2b) vanishes. The quantities $\varPhi, K,\alpha$ need not be single-valued in a torus. Note also that the separation of $\varPhi,K$ in the representation for $\boldsymbol {B}$ in (2.2b) is not unique since a shift of $\varPhi \to \varPhi + \varphi (\psi )$ results in a shift $K\to K-\varphi '(\psi )$ without changing $\boldsymbol {B}$.

For the sake of comparison, we first look at the standard representation given by

(2.3a)\begin{gather} \boldsymbol{B}=\boldsymbol{\nabla}\psi\times \boldsymbol{\nabla} \theta_B +\iota(\psi)\boldsymbol{\nabla} \phi_B \times \boldsymbol{\nabla} \psi, \end{gather}

(2.3b)\begin{gather}\boldsymbol{B}=I_B(\psi)\boldsymbol{\nabla}\theta_B + G_B(\psi)\boldsymbol{\nabla}\phi_B + B_\psi(\psi,\theta_B,\phi_B)\boldsymbol{\nabla}\psi. \end{gather}

Here, $(\psi,\theta _B,\phi _B)$ represents the standard Boozer coordinates (Boozer Reference Boozer1980, Reference Boozer1981a,Reference Boozerb) with $2{\rm \pi} \psi$ equal to the toroidal flux, $\iota (\psi )$ equal to the rotational transform and $(\theta _B,\varPhi _B)$ represent the poloidal and toroidal Boozer angles. The quantities $I_B(\psi ),G_B(\psi )$ and $B_\psi$ are all single-valued functions. The Boozer representation is an important special case of (2.2) for the choice

(2.4a--c)\begin{align} \varPhi= I_B(\psi)\theta_B + G_B(\psi) \phi_B, \quad K= B_\psi-\left( I_B'(\psi)\theta_B + G_B'(\psi) \phi_B \right), \quad \alpha=\theta_B -\iota(\psi)\theta_B. \end{align}

In this work, we use the Boozer–Grad coordinates instead of the usual Boozer coordinates.

Returning to the covariant representation of the magnetic field as given in (2.2b), we find that the chief advantage of this form is that the current, $\boldsymbol {J}=\boldsymbol {\nabla } \times \boldsymbol {B}$, has a Clebsch form

(2.5)\begin{equation} \boldsymbol{J}=\boldsymbol{\nabla} K\times \boldsymbol{\nabla} \psi, \end{equation}

which manifestly satisfies

(2.6)\begin{equation} \boldsymbol{J} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi=0. \end{equation}

Furthermore, the ideal MHS force balance

(2.7)\begin{equation} \boldsymbol{J}\times \boldsymbol{B}=\boldsymbol{\nabla} p, \end{equation}

takes a particularly simple form in terms of an MDE for $K$:

(2.8)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} K = p'(\psi). \end{equation}

The solution to (2.8) can be represented as

(2.9a--c)\begin{equation} K= \bar{K}(\psi,\alpha)+ G, \quad \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \bar{K} =0, \quad \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} G = p' , \end{equation}

where $\bar {K}(\psi,\alpha )$ is the homogeneous solution of the MDE and hence only a function of $\psi$ and possibly $\alpha$.

Let us collect the basic equations in the Weitzner formalism. These are obtained by equating the contravariant and the covariant representations of $\boldsymbol {B}$ as given in (2.2), together with (2.8),

(2.10a)\begin{gather} \boldsymbol{B} =\boldsymbol{\nabla} \psi \times \boldsymbol{\nabla} \alpha = \boldsymbol{\nabla} \varPhi + K \boldsymbol{\nabla} \psi, \end{gather}

(2.10b)\begin{gather}\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} K = p'(\psi). \end{gather}

Note that (2.10a), which we will call the basic MHS equations, is a set of three coupled nonlinear PDEs, while (2.10b) determines $K$. Together, they form a complete set of equations for the variables $(\varPhi,\psi,\alpha,K)$. We shall now discuss the boundary conditions that must be satisfied such that (2.10) has physically meaningful solutions in a torus.

In a torus, all physical quantities, such as magnetic fields and currents, must satisfy periodic boundary conditions in the two angles $(\theta,\phi )$, where $\theta$ is a poloidal angle, and $\phi$ is a toroidal angle. Unlike standard Boozer coordinates $(\theta _B,\phi _B)$, we do not require $(\theta,\phi )$ to be straight field line coordinates. This provides greater flexibility and freedom since we do not have any additional constraints on the angles. As mentioned, the functions $(\alpha,\varPhi,K)$ in a toroidal domain are multivalued. However, the multivalued part of the $\{\varPhi,\alpha,K\}$ system is not arbitrary due to the physical requirement that the magnetic field and the current must be single-valued quantities.

Following (Grad Reference Grad1971; Weitzner Reference Weitzner2016) we observe that in order for $\boldsymbol {B},\boldsymbol {J}$ (as given by (2.2a) and (2.5)) to be single-valued, only the $\boldsymbol {\nabla }\psi$ covariant components of $\boldsymbol {\nabla } \alpha$ and $\boldsymbol {\nabla } K$ can be multivalued. Furthermore, from (2.2b) we find that the secular terms of $\varPhi$ and $K$ must be related to each other such that $\boldsymbol {B}$ is single-valued.

As shown in Weitzner (Reference Weitzner2014), single valuedness of $\boldsymbol {B},\boldsymbol {J}$ in a torus implies that $\varPhi,\alpha$ and $K$ must be of the form

(2.11a)\begin{gather} \varPhi(\psi,\theta,\phi)= I(\psi)\theta + G(\psi)\phi + \tilde{\varPhi}(\psi,\theta,\phi), \end{gather}

(2.11b)\begin{gather}\alpha(\psi,\theta,\phi)= f(\psi)(\theta - \iota(\psi)\phi) + \tilde{\alpha}(\psi,\theta,\phi), \end{gather}

(2.11c)\begin{gather}K(\psi,\theta,\phi)={-}I'(\psi) \theta -G'(\psi) \phi +\tilde{K}(\psi,\theta,\phi), \end{gather}

where, $(\tilde {\varPhi },\tilde {\alpha },\tilde {K})$ denotes functions periodic in both $\theta$ and $\phi$, and $\iota (\psi )$ denotes the rotational transform of the magnetic field. The form (2.11) builds in the nestedness of flux surfaces. The single-valued flux functions $I(\psi ),G(\psi ),\iota (\psi )$ denote the same quantities as in standard Boozer coordinates (2.4a–c). We briefly discuss the physical significance of these flux functions below.

Without loss of generality, we can choose $2{\rm \pi} \psi$ to be the toroidal flux, as is common in the stellarator literature. As shown in Appendix A, with this choice, we can normalize $\alpha$ such that $f(\psi )=1$. The physical interpretation of $I(\psi ),G(\psi )$ is obtained by considering the net current, $\oint \boldsymbol {J}\boldsymbol {\cdot } \boldsymbol {dS}$, through a poloidal and a toroidal circuit. From (2.5) and (2.11c) it follows that $G'(\psi ),I'(\psi )$ are proportional to the net poloidal and toroidal current. In the vacuum limit ($K=0,p'=0$), the net toroidal current is zero, while the net poloidal current is a constant due to external currents. Consequently, for vacuum fields, we have $I(\psi )=0$ and $G(\psi )=G_0$ is a constant.

Equations (2.10), subject to the conditions (2.11), constitute the Weitzner formalism. The MHS equations (2.10a), (2.10b) are highly nonlinear, and in a generic stellarator, analytical progress is severely hindered by lack of any apparent continuous symmetry. Furthermore, the periodicity constraint (2.11) imposes a non-trivial requirement. In order to gain valuable physical insight, one, therefore, resorts to asymptotic analysis in some small parameter.

2.2 Mercier's NAE

Mercier developed one such asymptotic expansion scheme to study the behaviour near the magnetic axis by utilizing the distance from the magnetic axis as the expansion parameter. In the following, we shall briefly describe Mercier's near-axis formalism and Mercier's coordinates (Mercier Reference Mercier1964), which are also known as the direct coordinates (Jorge et al. Reference Jorge, Sengupta and Landreman2020b).

2.2.1 Mercier coordinates

The Mercier coordinates are defined with respect to the magnetic axis, which is a closed magnetic field line. Although the magnetic axis can be elliptic or hyperbolic (Solov'ev & Shafranov Reference Solov'ev and Shafranov1970; Duignan & Meiss Reference Duignan and Meiss2021), we will only treat the elliptic case here. We will assume that the magnetic axis is a 3-D smooth closed curve of total length $L$, described by the Frenet–Serret frame

(2.12a--d)\begin{equation} \dfrac{{\rm d}\boldsymbol{r}_0}{{\rm d}\ell}=\boldsymbol{t},\quad \dfrac{{\rm d}\boldsymbol{t}}{{\rm d}\ell}=\kappa\boldsymbol{n}, \quad \dfrac{{\rm d}\boldsymbol{n}}{{\rm d}\ell}=\tau\boldsymbol{b}-\kappa \boldsymbol{t}, \quad \dfrac{{\rm d}\boldsymbol{b}}{{\rm d}\ell}={-}\tau\boldsymbol{n}, \end{equation}

where $\boldsymbol {r}_0$ is the position vector of the magnetic axis, $\ell$ is the arclength along the axis, the functions $\kappa (\ell ),\tau (\ell )$ are the curvature and torsion of the axis, and the vectors $(\boldsymbol {t},\boldsymbol {n},\boldsymbol {b})$ are the standard tangent, normal and binormal vectors in the Frenet–Serret frame. As is well known, Frenet frames are singular near points of $\kappa =0$. Since magnetic axes with vanishing curvatures are crucial for QI stellarators, signed Frenet frames are typically used to describe such curves (Plunk et al. Reference Plunk, Landreman and Helander2019; Mata & Plunk Reference Mata and Plunk2023). Here, for simplicity, we do not consider such cases. However, the extension to a signed Frenet frame is straightforward.

Now, let us consider a tube of radius $\rho$ with the magnetic axis as its axis. We can define a poloidal angle $\theta$ such that a point on the tube is described by

(2.13)\begin{equation} \boldsymbol{r}=\boldsymbol{r}_0 + \boldsymbol{\rho}, \quad \boldsymbol{\rho}= x\; \boldsymbol{n}(\ell) + y \;\boldsymbol{b}(\ell), \quad x=\rho \cos \theta,y=\rho \sin{\theta}. \end{equation}

The coordinate system $(\rho,\theta,\ell )$ is not orthogonal when the torsion of the axis, $\tau (\ell )$, is non-zero. To obtain an orthogonal coordinate system, Mercier (Reference Mercier1964) replaced $\theta$ by

(2.14)\begin{equation} \omega =\theta +\int\tau \,{\rm d}\ell. \end{equation}

The line and volume elements, $({\rm d}s^2,{\rm d}V)$, of the Mercier coordinates are as follows:

(2.15a--c)\begin{equation} {\rm d}s^2\equiv |{\rm d}\boldsymbol{r}|^2=({\rm d}\rho)^2+(\rho \,{\rm d}\omega)^2+(h\,{\rm d}\ell)^2, \quad {\rm d}V=\rho h\, {\rm d}\rho\, {\rm d}\omega\, {\rm d}\ell, \quad h=1-\kappa \rho \cos{\theta}. \end{equation}

In other words, the $(\rho,\omega,\ell )$ coordinates are orthogonal with $(1,\rho,h)$ as the scale factors, and the metric tensor is given by

(2.16a,b)\begin{equation} g_{ij}=\text{diag}\left( 1,\rho^2,h^2 \right), \quad \sqrt{g}\equiv\sqrt{\text{Det}(g_{ij})}=\rho h. \end{equation}

Associated with the Mercier coordinates are the following orthonormal vectors $(\boldsymbol {\hat {\rho }},\boldsymbol {\hat {\omega }},\boldsymbol {\hat {\ell }})$ such that

(2.17a--c)\begin{align} \boldsymbol{\hat{\rho}}\equiv \boldsymbol{\nabla} \rho= \cos{\theta}\, \boldsymbol{n}+\sin{\theta}\, \boldsymbol{b}, \quad \boldsymbol{\hat{\omega}}\equiv \rho \boldsymbol{\nabla} \omega={-}\sin{\theta}\,\boldsymbol{n}+\cos{\theta}\,\boldsymbol{b}, \quad \boldsymbol{\hat{\ell}}\equiv h \boldsymbol{\nabla} \ell=\boldsymbol{t}. \end{align}

There is a certain freedom in choosing the poloidal and toroidal angles in the direct coordinates $(\rho,\omega,\ell )$. While the toroidal angle follows from the arclength $\ell$, one can generally add any arbitrary function $\delta (\ell )$, where $\delta '(\ell )$ is single-valued, to $\theta$ to obtain another poloidal angle without changing the Jacobian of transformation. Besides $(\theta,\phi =2{\rm \pi} \ell /L)$ another possible choice is $(u,\phi )$, where

(2.18a,b)\begin{equation} u=\theta+\delta(\ell)= \omega-\int \tau \,{\rm d}\ell +\delta(\ell) ,\quad \phi = 2{\rm \pi}\frac{\ell}{L}. \end{equation}

Let us now describe the Weitzner system described in § 2.1 using the Mercier coordinates. To connect to the contravariant form of $\boldsymbol {B}$ (2.2a) in the Weitzner formalism, we shall make a choice $(\theta,\phi )$ for defining angular coordinates. To connect to the covariant form of $\boldsymbol {B}$ given by (2.2b) in the Weitzner formalism, we shall utilize the orthonormal vectors $(\boldsymbol {\hat {\rho }},\boldsymbol {\hat {\omega }},\boldsymbol {\hat {\ell }})$. Using (2.15a–c), we can write the magnetic field in the following component form:

(2.19)\begin{gather} \boldsymbol{B}=B_\rho \boldsymbol{\hat{\rho}}+B_\omega \boldsymbol{\hat{\omega}}+ B_\ell \boldsymbol{\hat{\ell}}, \end{gather}

(2.20a--c)\begin{gather} B_\rho= \varPhi_{,\rho}+K \psi_{,\rho}, \quad B_\omega= \frac{1}{\rho}\left( \varPhi_{,\omega}+K \psi_{,\omega}\right), \quad B_\ell= \frac{1}{h}\left( \varPhi_{,\ell}+K \psi_{,\ell}\right). \end{gather}

Here and elsewhere, we shall use the notation $X_{,\ell }$ to denote the partial derivative of $X$ with respect to $\ell$ holding the other two Mercier coordinates $(\rho,\omega )$ fixed.

2.2.2 Near-axis expansions

The fundamental idea underlying the NAE is to solve the MHS equations using perturbation theory, where the expansion parameter scales like the distance from the magnetic axis, $\rho$. We can formally define the dimensionless expansion parameter $\epsilon \ll 1$ as

(2.21)\begin{equation} \epsilon=\kappa_{\text{max}}a, \end{equation}

where $\kappa _{\text {max}}$ is the maximum axis curvature and $0\leq \rho < a$. With this definition, we have $\kappa \rho <1$ throughout in the domain of validity of the NAE, which is necessary for the metric $\sqrt {g}=\rho h$ to be positive everywhere. Expanding the quantities $(\psi,\alpha,\varPhi,K)$ in $\epsilon$ and substituting them in the basic equilibrium equations (2.10), one can solve the MHS equations order by order in $\epsilon \ll 1$.

It is usually assumed (Kuo-Petravic & Boozer Reference Kuo-Petravic and Boozer1987; Garren & Boozer Reference Garren and Boozer1991) that the physical quantities are sufficiently regular near the magnetic axis such that one can carry out a regular power series expansion in positive powers of $\rho$. Thus, any function of the form $\mathfrak {f}(\rho,\omega,\ell )$ is expanded as follows:

(2.22)\begin{equation} \mathfrak{f}(\rho,\omega,\ell)=\sum_{n=0}^\infty (\epsilon \rho)^n\, \mathfrak{f}^{(n)}(\omega,\ell). \end{equation}

The function $\mathfrak {f}^{(n)}(\omega,\ell )$ can be determined order by order by substituting (2.22) into the MHS equations (2.10). Moreover, if $\mathfrak {f}$ has a well-defined Taylor series near the magnetic axis, we would expect $\mathfrak {f}^{(n)}$ to be a polynomial of order $n$ in $(x,y)=\rho (\cos \theta,\sin \theta )$. In Jorge et al. (Reference Jorge, Sengupta and Landreman2020b) and Duignan & Meiss (Reference Duignan and Meiss2021), it has been demonstrated that for vacuum fields with nested surfaces, the functions $(\varPhi -G_0 \phi ),\psi$ are of such form. Therefore, for such regular functions, we have

(2.23)\begin{align} \mathfrak{f}(x,y,\ell)& =\sum_{n=0}^\infty \epsilon^n \sum_{m=0}^n \mathfrak{c}^{(n)}(\ell) x^m y^{n-m}\nonumber\\ & = \sum_{n=0}^\infty (\epsilon\rho)^n\:\sum_{m=0}^n\left( \mathfrak{f}^{(n)}_{mc}(\ell) \cos{\left( m \omega +\delta_m(\ell)\right)}+\mathfrak{f}^{(n)}_{ms}(\ell) \sin{\left( m \omega +\delta_m(\ell)\right)}\right), \end{align}

where $\left ( \mathfrak {c}^{(n)}(\ell ),\mathfrak {f}^{(n)}_{mc}(\ell ),\mathfrak {f}^{(n)}_{ms}(\ell ),\delta _m(\ell ) \right )$ are functions of $\ell$ that must be determined by substituting (2.23) into (2.10). We note that even though the series (2.23) looks like a regular Taylor series for $\mathfrak {f}(x,y,\ell )$ near the axis, the series may be divergent. The NAE, like many other perturbative expansions, is written formally here, with no claims regarding the convergence of the series.

We shall substitute the formal series (2.23) into the basic MHS equations (2.10) and collect powers of $\rho$ and the poloidal harmonics $(\cos {m \ell },\sin {m \ell })$. This then reduces the nonlinear PDE system (2.10) to a nonlinear ODE system for $(\mathfrak {f}^{(n)}_{mc}(\ell ),\mathfrak {f}^{(n)}_{ms}(\ell ))$, which is a significant step forward. Moreover, (2.23) already includes periodicity in $\omega$, and only the periodicity in $\ell$ remains to be imposed on the ODEs.

Before proceeding further, we should point out some issues that can lead to a lack of regularity of the NAE. Firstly, although the periodicity in $\ell$ can always be satisfied if the on-axis rotational transform is sufficiently irrational, the coefficients $(\mathfrak {f}^{(n)}_{mc}(\ell ),\mathfrak {f}^{(n)}_{ms}(\ell ))$ can formally exhibit singular behaviour when the on-axis rotational transform is close to a rational number. As shown in Mercier (Reference Mercier1964) and Solov'ev & Shafranov (Reference Solov'ev and Shafranov1970), the form of the singularities on a $(m,n)$ rational surface is $(\iota _0-n/m)^{-1}$, where $\iota _0$ is the on-axis rotational transform. In this work, we shall assume that $\iota _0$ is sufficiently irrational to formally avoid such resonances in the vicinity of the axis. Secondly, we note that the regular power series expansion (2.22) is not always guaranteed to be correct. The possibility of having weak logarithmic singularities on the axis in the form of $\rho ^n (\log \rho )^m, n>m>0$ has been pointed out in Weitzner (Reference Weitzner2016). Furthermore, functions such as $\alpha$ can not be represented in the analytic power series given by (2.23). In this work, we demonstrate that for a large class of equilibria, one can still construct a power series in $\rho$ valid to all orders in $\rho$ without encountering logarithmically singular terms. We also pay special attention to the multivalued function $\alpha, K$ and discuss their structure.

In the following, we discuss various limits of ideal MHS: vacuum, force-free and, finally, finite beta MHS equilibrium. We start with the vacuum limit because its mathematical structure is fundamental and the easiest to analyse. Subsequently, we will show that the force-free and MHS fields can be formulated in a form very similar to the vacuum fields.

3.1 Basic equations and the need for an alternative set of equations

In the vacuum limit, the basic equilibrium equations (2.10) reduce to

(3.1)\begin{equation} \boldsymbol{\nabla} \psi\times \boldsymbol{\nabla} \alpha = \boldsymbol{\nabla} \varPhi. \end{equation}

Projecting (3.1) onto the orthonormal vectors $(\boldsymbol {\hat {\rho }},\boldsymbol {\hat {\omega }},\boldsymbol {\hat {\ell }})$ defined in (2.17a–c), we get

(3.2a)\begin{gather} \left( \psi_{,\omega}\alpha_{,\ell}- \psi_{,\ell}\alpha_{,\omega}\right) = \rho h\,\varPhi_{,\rho}, \end{gather}

(3.2b)\begin{gather}\left( \psi_{,\ell}\alpha_{,\rho}- \psi_{,\rho}\alpha_{,\ell}\right) = \frac{h}{\rho}\varPhi_{,\omega}, \end{gather}

(3.2c)\begin{gather}\left( \psi_{,\rho}\alpha_{,\omega}- \psi_{,\omega}\alpha_{,\rho}\right) = \frac{\rho }{h} \varPhi_{,\ell}. \end{gather}

Equation (3.2) is a closed set of equations that, in principle, can be solved for the three unknowns $(\varPhi,\psi,\alpha )$. However, the coupling between the various quantities makes the system (3.2) difficult to tackle analytically in the direct coordinate approach to NAE. Therefore, instead of (3.2), we examine alternative equations that decouple $(\varPhi,\psi,\alpha )$ as much as possible.

To arrive at an alternative set of equations, we now look at the consistency conditions that directly follow from (3.2). To derive the first consistency condition, we differentiate (3.2b) and (3.2c) with respect to $\omega$ and $\ell$, respectively, and add them. Using the commutativity of the cross-derivatives for the functions $\psi$ and $\alpha$, we get a consistency condition

(3.3)\begin{equation} {-}\partial_\rho\left( \psi_{,\omega}\alpha_{,\ell}- \psi_{,\ell}\alpha_{,\omega}\right) = \partial_\omega \left( \frac{h}{\rho}\varPhi_{,\omega} \right)+ \partial_\ell \left( \frac{\rho }{h} \varPhi_{,\ell}\right) . \end{equation}

It follows from (3.2a) that (3.3) is nothing but the Laplace equation for $\varPhi$,

(3.4)\begin{equation} \Delta \varPhi = \frac{1}{\rho h}\frac{\partial}{\partial \rho }\left(\rho h \frac{\partial \varPhi}{\partial \rho}\right) +\frac{1}{\rho^2 h}\frac{\partial}{\partial \omega }\left( h \frac{\partial \varPhi}{\partial \omega}\right)+ \frac{1}{ h}\frac{\partial}{\partial \ell }\left( \frac{1}{h} \frac{\partial \varPhi}{\partial \ell}\right) =0, \end{equation}

which also follows from $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {B} = \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {\nabla } \varPhi =0$.

Two other consistency conditions directly follow from (3.1) and are given by

(3.5a)\begin{gather} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \psi =\boldsymbol{\nabla} \varPhi\boldsymbol{\cdot} \boldsymbol{\nabla} \psi=0, \end{gather}

(3.5b)\begin{gather}\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \alpha =\boldsymbol{\nabla} \varPhi\boldsymbol{\cdot} \boldsymbol{\nabla} \alpha=0 . \end{gather}

One might suppose that the complete set (3.2) could be replaced by the new set of (3.4) and (3.5), which are the MDEs for $\psi$ and $\alpha$. The Laplace equation only involves $\varPhi$ and the scale factors and is ideally suited for determining $\varPhi$. Once $\varPhi$ is known, we might determine the Clebsch variables $(\psi,\alpha )$ through (3.5a) and (3.5b). The benefit of using the characteristic equations, (3.5), is that $\psi$ and $\alpha$ can be decoupled. However, several mathematical intricacies need to be addressed before replacing the complete set (3.2), which we now highlight.

Firstly, (3.2) is a third-order system because it has three first-order derivatives in $(\rho,\omega,\ell )$, while the Laplace equation is second order and so is (3.5) because of the two first derivatives. Therefore, the replacement of (3.2) by (3.4) and (3.5) raises the order by one since we have taken an additional derivative of the equations, which might lead to the inclusion of spurious solutions. Furthermore, the conditions (3.4), (3.5a) and (3.5b) are not completely independent when the equations are solved order by order. To see this, let us assume that the equations that govern the angular derivatives of $\varPhi$, (3.2b) and (3.2c), have been solved. It is then straightforward to see that the radial derivative of $\varPhi$ given by (3.2a) and characteristic equations (3.5a) and (3.5b) are not independent but are related through

(3.6)\begin{equation} \varPhi_{,\rho}- \frac{1}{\rho h}\left( \psi_{,\omega}\alpha_{,\ell}- \psi_{,\ell}\alpha_{,\omega}\right)=\frac{1}{\psi_{,\rho}}\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \psi =\frac{1}{\alpha_{,\rho}}\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \alpha. \end{equation}

As discussed in Appendix B, in order to retain the correct vacuum and MHS characteristics (Friedrichs & Kranzer Reference Friedrichs and Kranzer1958; Goedbloed Reference Goedbloed1983; Courant & Hilbert Reference Courant and Hilbert2008; Weitzner Reference Weitzner2014) we cannot simply replace the complete set (3.2) by the consistency conditions (3.4), (3.5a) and (3.5b). Fortunately, a workaround is possible that allows a partial decoupling of the variables $\varPhi,\psi,\alpha$. In the following sections, we will provide a recipe to carry out the finite-beta NAE efficiently and self-consistently to arbitrary orders by choosing the right set of fundamental equations. We will justify all the necessary mathematical steps below.

3.2 Alternative set of basic equations

We first observe that the consistency condition (3.3) implies that if the $\omega$ and $\ell$ components of the basic vacuum equations (3.1), i.e. (3.2b) and (3.2c), are satisfied, the solution of the Laplace equation satisfies the $\rho$ component, (3.2a), up to an arbitrary function of $(\omega,\ell )$. Therefore, (3.2a) needs to be imposed to eliminate this arbitrariness. As can be seen from (3.6) (and further discussed later on), the equation for the radial derivative of $\varPhi$ (3.2a) is identical to the MDE for $\psi$ (3.5a) order by order in $\rho$. Furthermore, it will be shown later that once (3.4) and (3.5a) (equivalently (3.2a)) are satisfied, the (3.2b) and (3.2c) are compatible, in the sense that either can be used to solve for $\alpha$. We will find it more straightforward to use (3.2b) to solve for $\alpha$ up to a function of $\ell$ by integrating with respect to $\omega$. The reason behind choosing (3.2b) over (3.2c) is that the $\varPhi$ and $\psi$ can be shown to have a finite number of poloidal harmonics which depend only on the order of expansion, whereas no such restriction exists for the toroidal harmonics.

Finally, to determine $\alpha$ fully, we need only the poloidal average of either (3.2c) or (3.5b). We use the latter because of the relevance of the MDE to calculating rotational transform and its higher derivatives. Note that the only MDE we solve is the MDE for $\psi$, (3.5a). We are not solving the MDE for $\alpha$; we are simply fixing the $\ell$ dependent part of $\alpha$ that was left undetermined because of the partial integration with respect to $\omega$. The procedure will be explained in detail in §§ 3.3–3.6.

In summary, we have seen that equating the covariant and the contravariant representations of $\boldsymbol {B}$, (2.2), in the vacuum limit leads to a complete set of (3.2). However, the variables $(\varPhi,\psi,\alpha )$ are coupled nonlinearly, making (3.2), not the best set of equations to start the NAE programme in direct coordinates. The alternative is to solve the Laplace equation for $\varPhi$ (3.4) and the MDE for $\psi$ (3.5a). We use the $\varPhi _{,\omega }$ equation, (3.2b), to obtain $\alpha$ up to a function of $\ell$, which is then determined from the poloidally averaged MDE for $\alpha$ (3.5b). The Laplace equation does not introduce extraneous solutions because we are solving the MDE for $\psi$, which is equivalent to solving (3.2a), as can be seen from (3.6). Furthermore, the Laplace equation is also the necessary solvability condition for the function $\alpha$.

3.3 Derivation of the vacuum NAE equations: the lowest order

We now expand the quantities $(\varPhi,\psi,\alpha )$ in power series of $\rho$ in the general form

(3.7a)\begin{gather} \varPhi(\rho,\omega,\ell)=\varPhi^{(0)}(\omega,\ell)+\rho\, \varPhi^{(1)}(\omega,\ell)+\rho^2 \varPhi^{(2)}(\omega,\ell)+\rho^3 \varPhi^{(3)}(\omega,\ell)+\cdots, \end{gather}

(3.7b)\begin{gather}\alpha(\rho,\omega,\ell)=\alpha^{(0)}(\omega,\ell)+\rho\, \alpha^{(1)}(\omega,\ell)+\rho^2 \alpha^{(2)}(\omega,\ell)+\rho^3 \alpha^{(3)}(\omega,\ell)+\cdots, \end{gather}

(3.7c)\begin{gather}\psi(\rho,\omega,\ell)=\rho^2 \psi^{(2)}(\omega,\ell)+\rho^3 \psi^{(3)}(\omega,\ell)+\rho^4 \psi^{(4)}(\omega,\ell)\cdots. \end{gather}

The analysis of the lowest-order quantities is important for the subsequent development of the analysis. Therefore, we now discuss in some detail the lowest-order structure. Firstly, note that the expansion of $\psi$ starts at $\rho ^2$ since we need $\psi$ to be zero on the axis. We also need $\psi _{,\rho }$ to vanish on the axis to ensure that the poloidal field, $B_\omega =(1/\rho )\varPhi _\omega \approx \psi _{,\rho }\alpha _{,\ell }$, is also zero on the magnetic axis. Furthermore, the vanishing of the poloidal field on the axis also leads to the $O(1/\rho ^2)$ and $O(1/\rho )$ of the vacuum system (3.2c) to be satisfied provided

(3.8a,b)\begin{equation} \varPhi^{(0)}(\omega,\ell)=\varPhi_0(\ell),\quad\varPhi^{(1)}(\omega,\ell)=0. \end{equation}

We assume that the strength of the magnetic field on the axis is positive and single-valued, i.e.

(3.9a,b)\begin{equation} B_0(\ell)=\varPhi_0'(\ell) >0, \quad B_0(\ell +L)=B_0(\ell). \end{equation}

To make $B_0(\ell )$ single-valued, $\varPhi ^{(0)}$ must be of the form

(3.10a,b)\begin{equation} \varPhi^{(0)}(\ell)=\int B_0(\ell) \,{\rm d}\ell = \bar{B}_0 \ell + \tilde{\varPhi}_0(\ell),\quad \bar{B}_0=\frac{1}{L}\int_0^L B_0(\ell)\,{\rm d}\ell, \end{equation}

where, $\tilde {\varPhi }_0(\ell )$ is periodic in $\ell$.

Expanding the vacuum system (3.2) to $O(1)$, we obtain the following coupled equations for $(\varPhi ^{(2)},\psi ^{(2)},\alpha ^{(0)})$:

(3.11a)\begin{gather} 2\varPhi^{(2)} =\alpha_{,\ell}^{(0)} \psi^{(2)}_{,\omega}-\alpha_{,\omega}^{(0)} \psi^{(2)}_{,\ell}, \end{gather}

(3.11b)\begin{gather}2\psi^{(2)} \alpha_{,\omega}^{(0)} = \varPhi_0'(\ell) , \end{gather}

(3.11c)\begin{gather}2\psi^{(2)} \alpha_{,\ell}^{(0)} ={-}\varPhi^{(2)}_{,\omega}. \end{gather}

Multiplying (3.11a) by $2\psi ^{(2)}$ and using ((3.11b), (3.11c)) we get the following linear homogeneous equation for $\psi ^{(2)}$:

(3.12)\begin{equation} \left( \varPhi_0'(\ell)\partial_\ell + \varPhi^{(2)}_{,\omega} \partial_\omega+4 \varPhi^{(2)} \right)\psi^{(2)}=0. \end{equation}

To solve (3.12) we need $\varPhi ^{(2)}$. The compatibility condition $\partial _\ell \alpha _{,\omega }^{(0)}=\partial _\omega \alpha _{,\ell }^{(0)}$ with the use of ((3.11b), (3.11c)) leads to the following Poisson equation for $\varPhi ^{(2)}$:

(3.13)\begin{equation} (\partial^2_{\omega}+2^2)\varPhi^{(2)}+\varPhi''_0=0. \end{equation}

Since (3.13) is a simple harmonic oscillator, its solution can be easily seen to be of the form

(3.14)\begin{equation} \varPhi^{(2)}={-}\frac{1}{4}\varPhi''_0(\ell) + \varPhi_0'(\ell) \left( f^{(2)}_{c2}(\ell)\cos{(2u(\omega,\ell))} +f^{(2)}_{s2}(\ell)\sin{(2u(\omega,\ell))} \right), \end{equation}

where, $f^{(2)}_{c2}(\ell ),f^{(2)}_{s2}(\ell )$ are the ‘integration constants’ with respect to the $\partial _\omega$ operator. We have also chosen $u$, defined by (2.18a,b), as the helical angle. Note that from the form of $\varPhi ^{(0)}$ as given in (3.10a,b), $\varPhi ^{(2)}$ can be made single-valued in $\ell$ by ensuring that the integration constants are periodic. The function $\varPhi ^{(2)}_{,\omega }$ obtained from (3.14) is harmonic with zero net poloidal average.

The solution of the equation for $\psi ^{(2)}$, (3.12), is of the form

(3.15)\begin{equation} \psi^{(2)}= \varPhi_0'(\ell)(a(\ell)+b(\ell)\cos{(2 u)}), \end{equation}

which can be checked by substituting (3.15) in (3.12) and collecting the poloidal harmonics $1,\cos {2u},\sin {2u}$. Equating the coefficients of the various poloidal harmonics to zero, we get

(3.16a,b)\begin{equation} a(\ell)=c_0 \cosh{\eta(\ell)}, \quad b(\ell)=c_0 \sinh{\eta(\ell)}, \end{equation}

together with

(3.17a,b)\begin{equation} f^{(2)}_{c2}(\ell)={-}\frac{b'(\ell)}{4 a(\ell)}={-}\frac{1}{4}\eta'(\ell),\quad f^{(2)}_{s2}(\ell)= \frac{b(\ell)}{2a(\ell)}\left( \delta'(\ell) -\tau(\ell)\right). \end{equation}

Here $\eta (\ell )$ is an arbitrary function of $\ell$. We require that $\left ( f^{(2)}_{c2}(\ell ),f^{(2)}_{s2}(\ell ) \right )$ are single-valued functions by requiring $\eta (\ell ),\delta '(\ell )$ to be single-valued and periodic in $\ell$. With this choice, we make both $\varPhi ^{(2)}, \psi ^{(2)}$ single-valued functions. Later, we determine the constant $c_0$ by imposing the choice of $\psi$ as the toroidal flux. To see that this is indeed the only solution, we can divide the equations for the $\omega,\ell$ derivatives of $\alpha ^{(0)}$, (3.11b) and (3.11c), by $2\psi ^{(2)}$ and eliminate $\alpha ^{(0)}$ through cross-differentiation. The resultant consistency condition

(3.18)\begin{equation} \partial_\ell \left( \frac{\varPhi_0'}{\psi^{(2)}}\right)+ \partial_\omega \left( \frac{\varPhi^{(2)}_{,\omega}}{\psi^{(2)}}\right)=0, \end{equation}

is a linear PDE for $(1/\psi ^{(2)})$ with periodic boundary conditions in $\omega,\ell$. Since the expression for $\psi ^{(2)}$ satisfies (3.18) and the periodic boundary conditions, it must be the unique solution.

For the physical interpretation of the various functions $\eta (\ell ),\delta (\ell )$, we refer to the excellent review article, Helander (Reference Helander2014). In short, the flux-surface shape obtained from the equation for $\psi ^{(2)}$ (3.15) is a rotating ellipse. The function $\eta (\ell )$ controls the ellipticity, while $\delta (\ell )$ measures its angle of rotation respect to the Frenet frame. We note that the choices of the free-functions in $f^{(2)}_{c2},f^{(2)}_{s2}$ in (3.17a,b) were made in terms of the shaping parameters of the lowest-order flux-surface. Later, we make choices such that the interpretation of higher-order free-functions $\varPhi ^{(n+2)}_{c},\varPhi ^{(n+2)}_{s}$ can be made in terms of higher-order flux surface shapes $\psi ^{(n+2)}_{c},\psi ^{(n+2)}_{s}$.

It is possible to find special coordinates such that $\psi ^{(2)}$ is a circle in these coordinates. Since $\varPhi _0'=B_0>0$ and $a(\ell )>b(\ell )\geq 0$, it follows that $\psi ^{(2)}$ is non-vanishing. As a result, we can define ‘rotating coordinates’ (Solov'ev & Shafranov Reference Solov'ev and Shafranov1970) or the leading-order ‘normal form’ coordinates

(3.19a,b)\begin{equation} X_{\mathcal{N}}= {\rm e}^{\eta/2}\rho \cos u , \quad Y_{\mathcal{N}}= {\rm e}^{-\eta/2}\rho \sin u, \end{equation}

such that $\psi ^{(2)}$, given by

(3.20)\begin{equation} \psi^{(2)} \rho^2= \frac{1}{2}\left( X^2_{\mathcal{N}}+Y_{\mathcal{N}}^2\right) \varPhi_0', \end{equation}

is a circle in the $(x_\mathcal {N}\equiv X_\mathcal {N}/\rho,y_\mathcal {N}\equiv Y_\mathcal {N}/\rho )$ coordinates. The benefits of using the ‘normal form’ coordinates have been noted in Duignan & Meiss (Reference Duignan and Meiss2021). As shown in Appendix F, these coordinates are particularly useful in solving the coupled ODEs that result from the MDE for higher-order flux-surface shapes.

With $\varPhi ^{(2)}$ and $\psi ^{(2)}$ in hand let us focus on obtaining $\alpha ^{(0)}$. Substituting the expressions for $\psi ^{(2)}$, (3.15), together with the expressions for $a,b$, (3.16a,b), in the equation for $\alpha ^{(0)}$, (3.11b), and integrating with respect to $\omega$ we get

(3.21)\begin{equation} \alpha^{(0)}=\frac{1}{2c_0}\arctan{\left( {\rm e}^{-\eta(\ell)}\tan{u}\right)}+\mathfrak{a}^{(0)}(\ell). \end{equation}

To determine the unknown function $\mathfrak {a}^{(0)}(\ell )$ we only need the poloidally averaged (3.11c). Since $\varPhi ^{(2)}_{,\omega }$ has zero average, the poloidally averaged $\alpha ^{(0)}_{,\ell }$ equation (3.11c) reads

(3.22)\begin{equation} \oint \frac{{\rm d}u}{2{\rm \pi}}\left( 2\psi^{(2)} \alpha_{,\ell}^{(0)} \right)=0. \end{equation}

Using the expressions for $\alpha ^{(0)},\psi ^{(2)}$ as given in (3.21) and (3.15), we can show that (3.22) yields

(3.23)\begin{equation} {-}2a(\ell)\, {\mathfrak{a}^{(0)}}'(\ell)=\partial_\ell u=\delta'(\ell)-\tau(\ell). \end{equation}

Therefore,

(3.24)\begin{equation} \alpha^{(0)}=\frac{1}{2c_0}\arctan{\left( {\rm e}^{-\eta(\ell)}\tan{u}\right)}-\int \,{\rm d}\ell \frac{\delta'(\ell)-\tau(\ell)}{2 c_0 \cosh{\eta(\ell)}}. \end{equation}

Let us briefly summarize the findings of NAE up to $O(1)$. From the requirement that the poloidal magnetic field and the flux surface label vanish on the axis, we can determine $\varPhi ^{(0)},\varPhi ^{(1)}$ to be as given in (3.8a,b). The vacuum system to $O(1)$ gives three coupled equations for $\varPhi ^{(2)},\psi ^{(2)},\alpha ^{(0)}$. We have shown that these equations can be decoupled to yield a Poisson equation for $\phi ^{(2)}$, (3.13) and an MDE for $\psi ^{(2)}$, (3.12). The solutions $\varPhi ^{(2)},\psi ^{(2)}$, given in (3.14) and (3.15), can then be used to determine $\alpha ^{(0)}$, which is given in (3.24) up to an overall constant. With the help of these lowest-order expressions, we now carry out the expansion to $O(n)$ in the next section, where we demonstrate that at $O(\rho ^n)$, an analogous decoupling of the higher-order quantities $(\varPhi ^{(n+2)},\psi ^{(n+2)},\alpha ^{(n)})$ is possible.

3.4 Derivation of the vacuum NAE equations: $n{\text {th}}$ order

The system (3.2) to $O(n)$ leads to the following closed set of equations:

(3.25a)\begin{gather} \left( {\alpha_{,\omega}^{(0)}} \partial_{\ell}-{\alpha_{,\ell}^{(0)}}\partial_{\omega}\right)\psi^{(n+2)}+ \left( {\alpha^{(n)}_{,\omega}} \partial_{\ell}-{\alpha^{(n)}_{,\ell}}\partial_{\omega}\right)\psi^{(2)}={-}(n+2) \varPhi^{(n+2)}+\cdots, \end{gather}

(3.25b)\begin{gather}\left( 2\psi^{(2)} \right)^{n/2+1} \partial_{\omega}\left(\dfrac{\alpha^{(n)}}{\left( 2\psi^{(2)} \right)^{n/2}}\right) = {\alpha_{,\omega}^{(0)}}\left(- (n+2)\:\psi^{(n+2)} \right) +\cdots, \end{gather}

(3.25c)\begin{gather}-\left( 2\psi^{(2)} \right)^{n/2+1} \partial_{\ell}\left(\dfrac{\alpha^{(n)}}{\left( 2\psi^{(2)} \right)^{n/2}}\right) = {\alpha_{,\ell}^{(0)}}\left( (n+2)\:\psi^{(n+2)} \right)+\varPhi^{(n+2)}_{,\omega}+\cdots, \end{gather}

where $\cdots$ refers to terms involving the lower-order known quantities.

Although (3.25) is linear in $(\varPhi ^{(n+2)},\psi ^{(n+2)},\alpha ^{(n)})$, it still couples these functions in a complicated manner. As discussed in § 3.1, an alternate approach to solving the system (3.25) is the system comprising of the Poisson equation for $\varPhi ^{(n+2)}$, the MDE for $\psi ^{(n+2)}$, (3.25b), and finally, the averaged MDE for $\alpha ^{(n)}$ to determine $\alpha ^{(n)}$ completely. In this section, we show how these alternative equations are equivalent to the complete set (3.25).

Expanding the equation relating the $\varPhi _{,\rho }$ equation and the MDEs for $\psi,\alpha$, (3.6), to $O(n)$ we observe that the equation for $\psi ^{(n+2)}$ as given in (3.25a) is equivalent to the $O(n)$ MDE for $\psi ^{(n+2)}$. However, the solution of the MDE for $\psi ^{(n+2)}$ does not uniquely determine $\psi ^{(n+2)}$ since a homogeneous solution that satisfies $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\psi$ to $O(n+2)$ can always be added to it. We show later that imposing the choice of $2{\rm \pi} \psi$ to be the toroidal flux, order by order, determines the homogeneous piece of $\psi$.

To solve the MDE for $\psi ^{(n+2)}$ or equivalently (3.25a) we need $\varPhi ^{(n+2)}$ as evident from the right-hand side of (3.25a). We have demonstrated in § 3.1 that the Laplace equation for $\varPhi$ arises as an exact compatibility condition of (3.2b) and (3.2c), provided (3.2a) holds. The same procedure works in higher orders as well. To $O(n)$, elimination of $\alpha ^{(n)}$ from (3.25b) and (3.25c) through cross-differentiation, i.e. $\partial _\ell$ (3.25c)$+ \partial _\omega$ (3.25b), followed by elimination of the derivatives of $\psi ^{(n+2)}$ using (3.25a), leads to an equation for $\varPhi ^{(n+2)}$ of the form

(3.26)\begin{equation} \left( \partial^2_\omega +(n+2)^2 \right) \varPhi^{(n+2)}=F^{(n+2)}_\varPhi, \end{equation}

where $F^{(n+2)}_\varPhi$ denotes known inhomogeneous forcing terms from lower-order quantities. Note that (3.26) is simply the Laplace equation for $\varPhi$ to $O(n+2)$ that leads to a Poisson equation for $\varPhi ^{(n+2)}$.

Now that the equation for its angular derivatives (3.25b) and (3.25c) are compatible thanks to (3.26), we use (3.25b) to solve for $\left ( \alpha ^{(n)}\right )/\left ( 2\psi ^{(2)}\right )^{n/2}$ up to an arbitrary function of $\ell$ say $\mathfrak {\bar {a}}^{(n)}(\ell )$. To determine $\mathfrak {\bar {a}}^{(n)}(\ell )$ we need the poloidally averaged MDE for $\alpha ^{(n)}$. The MDE for $\alpha ^{(n)}$ can be obtained through an algebraic combination of (3.25b) and (3.25c) that eliminates $\psi ^{(n+2)}$. It reads

(3.27)\begin{equation} \left( 2\psi^{(2)}\right)^{n/2+1}\left( {\alpha_{,\omega}^{(0)}} \partial_{\ell}-{\alpha_{,\ell}^{(0)}}\partial_{\omega}\right) \left(\dfrac{\alpha^{(n)}}{\left( 2\psi^{(2)}\right)^{n/2}}\right) ={-}{\alpha_{,\omega}^{(0)}}\left( \varPhi^{(n+2)}_{,\omega}+\dots \right)+\cdots. \end{equation}

Thus, we have derived an alternate set of equations which decouples $\varPhi ^{(n+2)},\psi ^{(n+2)},\alpha ^{(n)}$. The next step is to impose the periodicity constraint (2.11). To understand the nature of the constraints (2.11) in a torus, it turns out that inverse coordinates, where $\psi$ is used as the radial coordinate, are ideally suited. Connecting the direct and the inverse coordinates enables us to establish a direct relation between higher derivatives of the rotational transform, $\iota$, and the secular pieces of $\alpha ^{(n)}$. Therefore, in § 3.5, we discuss the map between the direct and the inverse to understand the periodicity constraints better.

In what follows, our strategy would be to solve the following equations for $(\varPhi ^{(n+2)},\psi ^{(n+2)},\alpha ^{(n)})$:

1. (I) the Poisson equation (3.26) for $\varPhi ^{(n+2)}$;

2. (II) the MDE for $\psi ^{(n+2)}$ (3.25a);

3. (III) the $\alpha _{,\omega }$ equation (3.25b) to solve for $\alpha ^{(n)}$ up to a function $\mathfrak {\bar {a}}^{(n)}(\ell )$;

4. (IV) the poloidal $\omega$ average of the MDE for $\alpha$ (3.27) to determine $\mathfrak {\bar {a}}^{(n)}(\ell )$.

Adding currents and plasma beta does not appreciably change the overall structure, as we show later. Therefore, a thorough understanding of the mathematical structure of the vacuum system (I–IV) and its solution is of merit and will be discussed in detail in the next few sections.

We discuss boundary conditions for Steps I–IV, derived from the periodicity requirements of physical quantities in the poloidal and toroidal angles, in § 3.5. Then we discuss the solutions to the Poisson equation (3.26) in § 3.6.1, the MDE for $\psi ^{(n+2)}$ in § 3.6.3 and in § 3.6.5 we provide the details of obtaining $\alpha$ following Steps III and IV.

3.5 Derivation of the vacuum NAE equations: boundary conditions in toroidal geometry and periodicity constraints

To obtain solutions of the equations (Steps I–IV) in a torus, we need to impose double periodicity in the two angles and regularity in the radial coordinate near the axis. In § 2, we presented the secular and periodic pieces of $\varPhi,\alpha$ in (2.11) in the inverse coordinates. The structure of the functions of $\varPhi,\psi,\alpha$ follows from the physical requirement that $\boldsymbol {J},\boldsymbol {B}$ be periodic in a torus. Since the Mercier direct coordinates do not use the flux-surface label as a coordinate, separating secular and periodic parts of the multivalued functions $\varPhi,\alpha$ is challenging. In particular, the form of $\alpha$ is crucial for understanding the MDEs. Hence, before discussing how to solve the equations, we need to understand the periodicity constraints in the direct coordinates.

With the choice of $2{\rm \pi} \psi$ to be the toroidal flux, we get the following expression for vacuum $(\alpha,\varPhi )$ from (2.11):

(3.28a,b)\begin{equation} \alpha=\theta - \iota(\psi) \phi +\tilde{\alpha}(\psi,\theta,\phi), \quad \varPhi=G_0 \phi +\tilde{\varPhi}(\psi,\theta,\phi). \end{equation}

To work out the NAE in inverse Mercier coordinates $(\psi,\omega,\phi =\ell / L)$ one needs to invert the function $\psi (\rho,\omega,\ell )$ to obtain $\rho (\psi,\omega,\ell )$. Equivalently one needs to invert the NAE for $\psi$ (3.7c) to obtain an expression for $\rho$ as a power series in $\sqrt {2\psi }$. We can asymptotically invert the power series of $\psi$ in $\rho$ to obtain

(3.29)\begin{equation} \rho=\sqrt{\frac{\psi}{\psi^{(2)}}}-\frac{\psi^{(3)}}{2\psi^{(2)}}\left( \frac{\psi}{\psi^{(2)}}\right) -\frac{1}{8}\left( 4\frac{\psi^{(4)}}{\psi^{(2)}}-5\left( \frac{\psi^{(3)}}{\psi^{(2)}}\right)\right) \left(\frac{\psi}{\psi^{(2)}}\right)^{3/2}+O(\psi^2). \end{equation}

Let us first focus on the field-line label $\alpha$. The NAE of $\alpha$ in the inverse (denoted by subscripts) and the direct (denoted by superscripts) coordinates,

(3.30a)\begin{gather} \alpha(\rho,\omega,\ell)= \alpha^{(0)}(\omega,\ell) +\rho \alpha^{(1)}(\omega,\ell)+\rho^2 \alpha^{(2)}(\omega,\ell)+O(\rho^3), \end{gather}

(3.30b)\begin{gather}\alpha(\psi,\omega,\ell)=\alpha_0(\omega,\ell) +\sqrt{\psi} \alpha_1(\omega,\ell)+\psi\, \alpha_2(\omega,\ell)+O(\psi^{3/2}), \end{gather}

are, therefore, connected through the relations

(3.31)\begin{equation} \left. \begin{aligned} \alpha_0 & = \alpha^{(0)},\quad \alpha_1= \frac{\alpha^{(1)}}{\sqrt{\psi^{(2)}}}, \quad \alpha_2=\frac{\alpha^{(2)}}{\psi^{(2)}}-\frac{1}{2}\frac{\alpha^{(1)}}{\sqrt{\psi^{(2)}}}\frac{\psi^{(3)}}{\left( \psi^{(2)}\right)^{3/2}},\\ \alpha_n & = \frac{\alpha^{(n)}}{(\psi^{(2)})^{n/2}}+\sum_{m=1}^{n-1} \frac{\alpha^{(m)}}{\left( \psi^{(2)} \right)^{m/2}} \mathcal{P}_{m(n-1)}, \end{aligned} \right\} \end{equation}

where, $\mathcal {P}_{m(n-1)}$ is a polynomial of order $(n-1)$ whose arguments are $\left ( \psi ^{(j)}/\left ( \psi ^{(2)}\right )^{j/2} \right )$ with $2\leq j\leq m+1$. Expanding $\iota (\psi )$ in a Taylor series in $\psi$ near the axis,

(3.32)\begin{equation} \iota(\psi)= \iota_0 +\iota_2 \psi + \dots +\iota_{2k}\psi^{(2k)} +O(\psi^{(2k+1)}), \end{equation}

and using the form of $\alpha$ in the inverse coordinates as given in (2.11), we find that

(3.33a--d)\begin{align} \alpha_0=\theta-2{\rm \pi}\iota_0 \frac{\ell}{ L}+\tilde{\alpha}_0, \quad \alpha_1 = \tilde{\alpha}_1, \quad \alpha_2 ={-}2{\rm \pi}\iota_2 \frac{\ell}{ L} + \tilde{\alpha}_2,\quad \alpha_{2k}={-}2{\rm \pi}\iota_{2k} \frac{\ell}{L}+\tilde{\alpha}_{2k}, \end{align}

where $\tilde {\alpha }_k, k=0,1,2$ denotes functions periodic in the angles. Averaging the $\omega,\ell$ derivatives of $\alpha _{2k}$ from (3.33a–d) over $\omega$ and $\ell$, and exploiting the periodicity of $\tilde {\alpha }_k$, we find that

(3.34a,b)\begin{equation} \oint \frac{{\rm d}\omega}{2{\rm \pi}} \partial_\omega \alpha_{2k}=\delta_{2k},\quad 2{\rm \pi}\iota_{2k}={-}\oint \frac{{\rm d}\omega}{2{\rm \pi}}\int_0^L \,{\rm d}\ell\, \partial_\ell \left( \alpha_{2k}\right), \end{equation}

where, $\delta _{2k}=0$ for $k>0$ and $1$ for $k=0$.

We are now in a position to determine the constant $c_0$ in (3.16a,b). To do so, we impose the choice of $2{\rm \pi} \psi$ as the toroidal flux (equivalently $f(\psi )=1$ in (2.11b)), order by order, by insisting that the secular $\theta$ piece only appears in $\alpha ^{(0)}$ as shown in (3.33a–d). Here, we note that (3.33a–d) implies that the average of $\alpha _{,\omega }^{(0)}$ must be unity.

Using the form of $\alpha ^{(0)}$ obtained in (3.24) and averaging over $\omega$, we obtain the average $\alpha _{,\omega }^{(0)}$ in the following form:

(3.35)\begin{equation} \oint \frac{{\rm d}\omega}{2{\rm \pi}} \alpha_{,\omega}^{(0)}= \frac{1}{2a(\ell)} \oint \frac{{\rm d}\omega}{2{\rm \pi}}\frac{1}{1+\varepsilon(\ell) \cos{2u}}. \end{equation}

Here, we have used (3.11b) for $\alpha ^{(0)}_{,\omega }$, (3.15) for the functional form of $\psi ^{(2)}$, and (3.16a,b) to define

(3.36)\begin{equation} \varepsilon(\ell) \equiv \frac{b(\ell)}{a(\ell)}= \tanh{\eta(\ell)}. \end{equation}

Replacing the $\omega$ integral by a $u$ integral, using the identity

(3.37)\begin{equation} \oint \frac{{\rm d}u}{2{\rm \pi}}\frac{1}{1+\varepsilon \cos{2u}}=\frac{1}{\sqrt{1-\varepsilon^2}}=\cosh{\eta(\ell)}, \end{equation}

and finally imposing the condition that the average of $\alpha _{,\omega }^{(0)}$ must be unity, we find that

(3.38)\begin{equation} c_0=\tfrac{1}{2}, \quad a(\ell)= \cosh\eta(\ell). \end{equation}

Therefore, the expression for $\alpha ^{(0)}$, (3.24), now reads

(3.39)\begin{equation} \alpha^{(0)}=\arctan{\left( {\rm e}^{-\eta(\ell)}\tan{u}\right)}-\int \,{\rm d}\ell \frac{\delta'(\ell)-\tau(\ell)}{2a(\ell)}. \end{equation}

The on-axis rotational transform can be obtained from (3.39) by utilizing

(3.40)\begin{equation} \iota_0= \alpha(\ell,\theta)-\alpha(\ell+L,\theta)+N = \frac{v(L)-(\delta(L)-\delta(0))}{2{\rm \pi}}+N, \end{equation}

where we have defined

(3.41)\begin{equation} v(\ell)=\int_0^\ell \,{\rm d}s \frac{u'(s)}{\cosh{\eta(s)}}= \int_0^\ell \,{\rm d}s \frac{\delta'(s)-\tau(s)}{\cosh{\eta(s)}} . \end{equation}

The extra integer $N$ comes from the rotation of the Frenet frame (Helander Reference Helander2014).

We now explicitly derive the periodicity conditions in the direct coordinates. Firstly, every $\psi ^{(m+2)}, m=0,1,2,\ldots$ needs to be single-valued so that the flux-surface label $\psi$ is single-valued. Secondly, by substituting the relation between the components of $\alpha$ in the inverse representation $\alpha _n$ and the direct representation $\alpha ^{(n)}$ as given in (3.33a–d) into the identities (3.34a,b), we get

(3.42)\begin{equation} \left. \begin{gathered} \oint \frac{{\rm d}\omega}{2{\rm \pi}} \partial_\omega\left( \frac{\alpha^{(n)}}{(\psi^{(2)})^{n/2}}\right)=0 \quad (n>0), \\ \oint \frac{{\rm d}\omega}{2{\rm \pi}}\oint\frac{{\rm d}\ell}{2{\rm \pi}}\partial_\ell \left( \frac{\alpha^{(n)}}{(\psi^{(2)})^{n/2}} +\sum_{m=1}^{n-1} \frac{\alpha^{(m)}}{\left( \psi^{(2)} \right)^{m/2}} \mathcal{P}_{m(n-1)}\right) ={-}\iota_{2k} \quad (n=2k). \end{gathered} \right\} \end{equation}

Returning to the discussion of the multivalued structure of $\varPhi$, we note that the analysis is almost identical to that of $\alpha$. The analogues of the relation between inverse and direct components (3.30), and the order by order form in inverse-coordinates (3.31) for $\varPhi$ are obtained by simply replacing $\alpha$ by $\varPhi$. An important difference between $\varPhi$ and $\alpha$ occurs because in vacuum $G_0$ in (3.28a,b) is a constant. Therefore, we have

(3.43)\begin{equation} \varPhi_0=G_0 \frac{\ell}{ L}+\tilde{\varPhi}_0, \quad \varPhi_{n}=\tilde{\varPhi}_{n} \quad \text{for all }n>0, \end{equation}

which implies that except for $\varPhi ^{(0)}$, all other $\varPhi ^{(n)}$ are periodic since $\psi ^{(2)}$ is periodic. Single valuedness of $\varPhi ^{(n)}, n>0$ follows from the Laplace equation as will be shown in § 3.6.1.

To summarize, we have deduced the periodicity constraints on $(\varPhi ^{(n+2)},\psi ^{(n+2)},\alpha ^{(n)})$ in this section by connecting to the inverse coordinates. The flux label $\psi ^{(n+2)}$ must be periodic in both angles on a torus. In the inverse coordinates, the general periodicity requirement on $\varPhi,\alpha$ is given by (3.28a,b) and order by order the requirement is (3.33a–d) together with (3.43). In the direct coordinates, the periodicity constraints are given by (3.42) and single valuedness of $\varPhi ^{(n+2)}$ for all $n>0$.

3.6 Solutions of the vacuum NAE equations to $O(n)$

Our goal now is to analyse the harmonic contents of $(\varPhi ^{(n+2)},\psi ^{(n+2)},\alpha ^{(n)})$. The structure of the vacuum potential $\varPhi$ was analysed in detail in Jorge et al. (Reference Jorge, Sengupta and Landreman2020b). We briefly summarize some of the relevant results here. To analyse the harmonic content of $\psi ^{(n+2)}$ and $\alpha ^{(n)}$, we study the structure and properties of MDEs in the next sections.

3.6.1 Solution of Step I: the Poisson equation for $\varPhi ^{(n+2)}$

The Laplace equation to $O(n)$ is

(3.44)\begin{align} & \left( \partial^2_{\omega}+(n+2)^2 \right)\varPhi^{(n+2)}\nonumber\\ & \quad ={-}\sum_{m=0}^{n} \kappa^{m} \cos^{m} \theta\left[\vphantom{\frac{(m+1)(m+2)}{2}}\kappa\left(\partial_\omega \varPhi^{(n-m+1)} \sin \theta- (n-m+1)\varPhi^{(n-m+1)}\cos \theta\right)\right.\nonumber\\ & \qquad +\left.(m+1)\partial^2_{\ell}\varPhi^{(n-m)}+ \frac{(m+1)(m+2)}{2}\partial_\ell(\kappa \cos \theta) \partial_\ell\varPhi^{(n-m-1)}\right]. \end{align}

As shown in Jorge et al. (Reference Jorge, Sengupta and Landreman2020b), the right-hand side of (3.44) has at most harmonic terms of order $(n)$. If the $(n+2)$th poloidal harmonics were present in the forcing term of (3.44), the inversion of (3.44) would lead to linear (hence secular) $\omega$ terms in $\varPhi ^{(n+2)}$. Absence of $(n+2)$ forcing poloidal harmonics implies that $\varPhi ^{(n+2)}$ has poloidal up to $n+2$, with the coefficients of $(n+2){\textrm {th}}$ harmonics being free.

Therefore, we can write down the solution of $\varPhi ^{(n+2)}$ as

(3.45)\begin{equation} \varPhi^{(n+2)}= \sum_j \left( \varPhi^{(n+2)}_{cj}(\ell)\cos(j u)+ \varPhi^{(n+2)}_{sj}(\ell)\sin(j u)\right), \end{equation}

where, $u$ is defined by (2.18a,b). For even $n$, $j=2r$ and $r$ varies from $0$ to $(n/2+1)$. For odd $n$, $j=2r+1$ and $r$ varies from $0$ to $(n+1)/2$. The functions $\varPhi ^{(n+2)}_{c(n+2)}(\ell ),\varPhi ^{(n+2)}_{s(n+2)}(\ell )$ are the ‘free functions’, and the rest are completely determined by the lower-order quantities.

We now ensure that derivatives of $\varPhi ^{(n+2)}$ are single-valued since they are physically related directly to the components of the magnetic field. The right-hand side is manifestly periodic in the poloidal angle. Therefore, periodicity in the poloidal angle is maintained order by order. Let us now check for periodicity in $\ell$. The presence of $\varPhi ^{(n-m+1)}$ without derivatives in the forcing term on the right-hand side of (3.45) might suggest that $\varPhi ^{(n+2)}$ is not periodic since $\varPhi ^{(0)}$ is multivalued in $\ell$ after all. However, $\varPhi ^{(0)}$ never appears without any derivative in the forcing term. To have $\varPhi ^{(0)}$, we must have $m=n+1$, which is outside the summation range. Furthermore, $\varPhi ^{(1)}=0$ implies that for any $n$, the lowest non-trivial term of the form $\varPhi ^{(n-m+1)}$ is $\varPhi ^{(2)}$ for $m=n-1$, which is periodic. Therefore, the forcing term in (3.44) is a linear superposition of derivatives of $(\varPhi ^{(0)},\varPhi ^{(2)},\ldots,\varPhi ^{(m)})$ and $\varPhi ^{m}$ terms with $m>0$. By induction $\varPhi ^{(n+2)}$ must be single-valued for all $n>0$ as required by (3.43).

This completes our solution of Step I. The solution of Step II is more involved and will be developed over the next few sections.

3.6.2 Solution of Step II: analysis of the general MDE

To solve the MDE for $\psi$ in Step II, we need to discuss some properties of the MDE. In this section, we shall derive these properties for a general vacuum MDE.

The MDE (Newcomb Reference Newcomb1959) for a quantity $\chi$ takes the general form

(3.46)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \chi +F_\chi=0,\quad \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\equiv B_\rho \partial_\rho + B_\omega \frac{1}{\rho}\partial_\omega + B_\ell \frac{1}{h}\partial_\ell, \end{equation}

where $F_\chi$ is the forcing (or the inhomogeneous) term. Expressing the components of $\boldsymbol {B}$ in (3.46) in terms of derivatives of the scalar potential $\varPhi$, the MDO can be written as

(3.47)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}= \varPhi_{,\rho} \partial_\rho + \frac{1}{\rho^2}\varPhi_{,\omega}\partial_\omega +\frac{1}{h^2}\varPhi_{,\ell}\partial_\ell . \end{equation}

Expanding in powers of $\rho$ and using the expressions for $\varPhi ^{(0)}, \varPhi ^{(1)}$ from (3.8a,b), we find that (3.46) takes the form

(3.48)\begin{equation} (\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla})^{(n)}_0 \chi^{(n)} +\mathfrak{f}^{(n)}_\chi=0, \end{equation}

where

(3.49)\begin{equation} (\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla})^{(n)}_0\equiv \left( \varPhi_0'(\ell)\partial_\ell + \varPhi^{(2)}_{,\omega} \partial_\omega+2 n \varPhi^{(2)} \right). \end{equation}

The term $\mathfrak {f}^{(n)}_\chi$ consists of $\chi ^{(n-1)}$ and other lower-order quantities obtained from both $F_\chi$ and $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } \chi$ expanded to $O(n)$. Note also that the definition of the lowest order MDO depends on the order $n$ due to the $\partial _\rho$ operator in (3.47). Clearly, for $n=2$, $\psi ^{(2)}$ is a homogeneous solution of (3.48), as claimed earlier.

Unsurprisingly, we will find the MDO (3.49) repeatedly in this work. We shall therefore record a few essential properties of the MDO before we continue with the analysis of the MDE (3.48).

The first important property to note is that an alternate form of the MDO is given by

(3.50)\begin{equation} (\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla})^{(n)}_0\left( \boldsymbol{\cdot} \right)= 2\left( \psi^{(2)}\right)^{n/2+1}\left\{\alpha^{(0)},\frac{1}{ \left( \psi^{(2)}\right)^{n/2}}\left( \boldsymbol{\cdot} \right) \right\}_{(\omega,\ell)}, \end{equation}

where, $\{f,g\}_{(\omega,\ell )}$ denotes a standard Poisson bracket of functions $f$ and $g$ with respect to $(\omega,\ell )$. The form (3.50) follows from the definition of MDO (3.49), the equation for $\psi ^{(2)}$ (3.12) and the relations (3.11) that are satisfied by $(\varPhi ^{(2)},\alpha ^{(0)},\psi ^{(2)})$.

The Poisson-bracket form of the MDO, (3.50), allows us to readily find the homogeneous solutions of the MDE (3.48). Thus, if $\chi ^{(n)}_h$ is a homogeneous solution of (3.48), it must be of the form

(3.51)\begin{equation} \chi^{(n)}_h=\left( \psi^{(2)}\right)^{n/2}\mathcal{C}_h\left( \alpha^{(0)}\right). \end{equation}

When the rotational transform is irrational, $\mathcal {C}_h$ must be a constant if $\chi ^{(n)}_h$ represents a physical quantity.

The next important property of the MDO (3.50) is that it possesses an annihilator (an operator to whose kernel $(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla })^{(n)}_0$ belongs). To obtain the annihilator, we use the following property of the Poisson bracket

(3.52)\begin{equation} \left\{\alpha^{(0)},\chi\right\}_{(\omega,\ell)} \equiv \alpha_{,\omega}^{(0)} \chi_{,\ell}-\alpha_{,\ell}^{(0)} \chi_{,\omega}= \partial_\ell \left( \alpha_{,\omega}^{(0)} \chi \right)-\partial_\omega \left( \alpha_{,\ell}^{(0)} \chi \right). \end{equation}

Here we have used the fact that $\alpha _{,\omega }^{(0)},\alpha _{,\ell }^{(0)}$ are single-valued in order to commute the partial derivatives $\partial _\ell,\partial _\omega$ acting on $\alpha ^{(0)}$ in (3.52). Averaging (3.52) over $(\omega,\ell )$ we get

(3.53)\begin{equation} \oint \,{\rm d}\ell\oint \frac{{\rm d}\omega}{2{\rm \pi}} \left\{\alpha^{(0)},\chi \right\}_{(\omega,\ell)} =0, \end{equation}

provided $\chi$ is single-valued. Therefore, the operator

(3.54)\begin{equation} \langle \left( \boldsymbol{\cdot} \right) \rangle_{(n)} = \oint \frac{{\rm d}\ell}{2{\rm \pi}} \oint \frac{{\rm d}\omega}{2{\rm \pi}}\frac{1}{ \left( \psi^{(2)}\right)^{n/2+1}} \left( \boldsymbol{\cdot} \right) \end{equation}

is an annihilator for the MDO (3.50). Note that the factor of $2{\rm \pi}$ below $\textrm {d}\ell$ is because (3.50) has a derivative with respect to $\ell$. When written in terms of $\phi =2{\rm \pi} \ell /L$, the factor of $2{\rm \pi}$ turns out to be the usual measure.

Finally, we must discuss the poloidal harmonics generated when the MDO acts on the $n{\textrm {th}}$ harmonics. This property is essential in order to solve the MDE order by order. In the inverse coordinate approach, one can use analyticity to show that there is an upper bound on the poloidal harmonics at each order. In the direct coordinate approach, the upper bound on the poloidal harmonics is not evident since various powers of $\left ( \psi ^{(2)}\right )$ appear in the MDO (3.50) and its annihilator (3.54).

To deduce the poloidal harmonic content of $\chi ^{(n)}$ that satisfies the MDE (3.48), we will need another important property of the MDO (3.49). Let us consider a generic periodic function $\chi ^{(n)}_j$ of the form

(3.55)\begin{equation} \chi^{(n)}_j= \chi^{(n)}_{cj}(\ell) \cos{(j u(\omega,\ell))}+ \chi^{(n)}_{sj}(\ell) \sin{(j u(\omega,\ell))}, \end{equation}

with $u(\omega,\ell )$ defined in (2.18a,b). The integer $j$ such that $2\leq j \leq n$, is even (odd) if n is even (odd). Using (3.14) and (3.15) we can show that

(3.56a)\begin{gather} (\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla})^{(n)}_0 \chi^{(n)}_j = C^{(n)}_{cj}(\ell) \cos{j u}+C^{(n)}_{sj}(\ell) \sin{j u} +(n-j) C^{(n)}_{(+)}(\omega,\ell)+(n+j) C^{(n)}_{(-)}(\omega,\ell), \end{gather}

(3.56b)\begin{gather}C^{(n)}_{cj}(\ell)= \varPhi_0' \left( \frac{{\rm d}\chi^{(n)}_{cj}}{{\rm d}\ell} +n(\delta'-\tau)\chi^{(n)}_{sj}-\frac{n}{2}\frac{\varPhi''_0}{\varPhi_0'}\chi^{(n)}_{cj}\right), \end{gather}

(3.56c)\begin{gather}C^{(n)}_{sj}(\ell)= \varPhi_0' \left( \frac{{\rm d}\chi^{(n)}_{sj}}{{\rm d}\ell} -n(\delta'-\tau)\chi^{(n)}_{cj}-\frac{n}{2}\frac{\varPhi''_0}{\varPhi_0'}\chi^{(n)}_{sj}\right), \end{gather}

(3.56d)\begin{gather}C^{(n)}_{(+)}=\varPhi_0' \left[\left( f^{(2)}_{c2}\chi^{(n)}_{cj}- f^{(2)}_{s2}\chi^{(n)}_{sj} \right)\cos{(j+2)u}+\left( f^{(2)}_{c2}\chi^{(n)}_{sj}+ f^{(2)}_{s2}\chi^{(n)}_{cj} \right)\sin{(j+2)u} \right], \end{gather}

(3.56e)\begin{gather}C^{(n)}_{(-)}=\varPhi_0' \left[\left( f^{(2)}_{c2}\chi^{(n)}_{cj}+ f^{(2)}_{s2}\chi^{(n)}_{sj} \right)\cos{(j-2)u} -\left( f^{(2)}_{s2}\chi^{(n)}_{cj}- f^{(2)}_{c2}\chi^{(n)}_{sj} \right)\sin{(j-2)u}\right]. \end{gather}

We note, in particular, that for $j=n$, there are no terms with $(n+2)$ poloidal harmonics. Thus, for a function of the form (3.55), i.e. with $j{\textrm {th}}$ poloidal harmonics in $u$, where $j$ has the same parity as $n$ and $0\leq j\leq n$, $(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla })^{(n)}_0 \chi ^{(n)}_j$ is a function with poloidal harmonics strictly between $j-2$ and $\min {(j+2,n)}$.

Now let us define a function $\chi ^{(n)}=\sum _j \chi ^{(n)}_j$ such that

(3.57)\begin{equation} \chi^{(n)}= \begin{cases} \displaystyle\sum_{r=0}^{ n/2}\chi^{(n)}_{2r} , & n \text{ is even},\\ \\ \displaystyle\sum_{r=0}^{(n-1)/2 }\chi^{(n)}_{2r+1} , & n \text{ is odd}, \end{cases} \end{equation}

where, $\chi ^{(n)}_j$'s for all $j,n$ are of the form (3.55), $(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla })^{(n)}_0\chi ^{(n)}$ will have even (odd) poloidal harmonics only up to $n$ for even (odd) values of $n$. If we are given a function $\mathfrak {f}^{(n)}_\chi$, such that

(3.58)\begin{equation} \mathfrak{f}^{(n)}_\chi= \begin{cases} \displaystyle\sum_{r=0}^{ n/2} \left( F^{(n)}_{c(2r)} \cos{2r u} + F^{(n)}_{s(2r)} \sin{2r u} \right) , & n\text{ is even},\\ \\ \displaystyle\sum_{r=0}^{(n-1)/2 }\left( F^{(n)}_{c(2r+1)} \cos{(2r+1) u} + F^{(n)}_{s(2r+1)} \sin{(2r+1) u} \right) , & n \text{ is odd}, \end{cases} \end{equation}

the solution of

(3.59)\begin{equation} (\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla})^{(n)}_0 \chi^{(n)}+\mathfrak{f}^{(n)}_\chi=0 \end{equation}

is obtained by solving the coupled ODEs obtained by equating to zero the sum of the $(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla })^{(n)}_0 \chi ^{(n)}$ terms given in (3.56) and the $\mathfrak {f}^{(n)}_\chi$ terms given in (3.58).

We want to emphasize that without the property that no $(n+2)$ poloidal harmonics are generated when the MDO acts on the $n{\textrm {th}}$ harmonics, there would be no guarantee of an upper bound on the number of poloidal harmonics in $u$. Thanks to this property, a forcing term $\mathfrak {f}^{(n)}_\chi$ that has $n$ poloidal harmonics will lead to a solution $\chi ^{(n)}$, which also has at most $n$ poloidal harmonics.

3.6.3 Solution of Step II: solution of the MDE for $\psi ^{(n+2)}$

As we discussed in § 3.4, equation (3.25a) can be rewritten as a MDE for $\psi ^{(n+2)}$ by eliminating $\alpha ^{(n)}$ algebraically using (3.25b) and (3.25c). To show that explicitly, we need the following identity obtained from the $n{\textrm {th}}$ order NAE vacuum equations (3.25b), (3.25c):

(3.60)\begin{align} \left( 2\psi^{(2)}\right)\{\alpha^{(n)},\psi^{(2)}\}_{(\omega,\ell)}& =(n+2)\psi^{(n+2)}\{\alpha^{(0)},\psi^{(2)}\}_{(\omega,\ell)}+\psi^{(2)}_{,\omega}\varPhi^{(n+2)}_{,\omega}\ldots \nonumber\\ & =2(n+2)\psi^{(n+2)}\varPhi^{(2)} +\psi^{(2)}_{,\omega}\varPhi^{(n+2)}_{,\omega}\ldots , \end{align}

where, in the last step, we have used the $O(1)$ relation (3.11) that gives the Poisson bracket of $\alpha ^{(0)},\psi ^{(2)}$ in terms of $\varPhi ^{(2)}$.

Multiplying (3.25a) by $2\psi ^{(2)}$, using (3.11) and (3.60), we now rewrite (3.25a) as

(3.61)\begin{align} \left( \varPhi_0'(\ell)\partial_\ell + \varPhi^{(2)}_{,\omega} \partial_\omega+ 2(n+2) \varPhi^{(2)} \right)\psi^{(n+2)}+ \left( 2\psi^{(2)} (n+2)\varPhi^{(n+2)}+ \psi^{(2)}_{,\omega}\varPhi^{(n+2)}_{,\omega}\right)+\cdots=0, \end{align}

which is clearly of the form

(3.62)\begin{align} (\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla})^{(n+2)}_0 \psi^{(n+2)}+ F^{(n+2)}_\psi=0, \quad F^{(n+2)}_\psi= \left( 2\psi^{(2)} (n+2)\varPhi^{(n+2)}+ \psi^{(2)}_{,\omega}\varPhi^{(n+2)}_{,\omega}\right)+\cdots, \end{align}

where, $(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla })^{(n+2)}_0$ is defined by (3.49).

Firstly, following the discussion in § 3.6.2, we note that (3.62) determines $\psi ^{(n+2)}$ only up to the homogeneous solution,

(3.63)\begin{equation} {\psi_H}^{(n+2)}= \left( {2\psi^{(2)}}\right)^{({n+2})/{2}} \mathcal{Y}_H^{(n+2)}(\alpha^{(0)}). \end{equation}

For irrational rotational transform, $\mathcal {Y}_H^{(n+2)}$ must be a constant. The value of the constant will be determined in a subsequent section by imposing the constraint that $2{\rm \pi} \psi$ is the net toroidal flux.

To determine the particular solution, $\psi _P^{(n+2)}$, of the MDE (3.62), we first need to investigate the harmonic structure of the forcing term $F^{(n+2)}_\psi$. Let us convince ourselves that $F^{(n+2)}_\psi$ does not have any harmonics higher than $(n+2)$. We will start with the highest harmonics of $\varPhi ^{(n+2)}$ beating with $\psi ^{(2)}$ according to (3.62). Straightforward algebra using (3.15) shows that

(3.64)\begin{align} & 2\psi^{(2)} (n+2)\varPhi^{(n+2)}+ \psi^{(2)}_{,\omega}\varPhi^{(n+2)}_{,\omega} \nonumber\\ & \quad = 2(n+2)\varPhi_0'\left[ a(\ell) \left( \varPhi^{(n+2)}_{c(n+2)}\cos{(n+2)u}+ \varPhi^{(n+2)}_{c(n+2)}\sin{(n+2)u} \right) \right.\nonumber\\ & \qquad \left.+b(\ell)\left( \varPhi^{(n+2)}_{c(n+2)}\cos{n u}+ \varPhi^{(n+2)}_{s(n+2)}\sin{n u} \right) \right]+\cdots . \end{align}

The rest of the terms are similar, as shown inductively in Appendix D. Given this form for the forcing term, the solution of $\psi ^{(n+2)}$ can then be deduced from the solution of $\chi ^{(n)}$ as given in (3.57) with $n$ replaced by $(n+2)$.

We shall now briefly analyse the coupled ODEs that result from equating the poloidal harmonics of the MDE (3.62). A more complete description will be presented in the next section. Our goal in the remainder of this section is to motivate a choice for the free-functions, $\varPhi ^{(n+2)}_{c(n+2)}, \varPhi ^{(n+2)}_{s(n+2)}$, that appear in the forcing term (3.64), in terms of the $O(n+2)$ flux-surface shaping functions $\psi ^{(n+2)}_{c(n+2)}, \psi ^{(n+2)}_{s(n+2)}$.

We first note that using the expressions for $\varPhi ^{(2)}$ as given in (3.14) together with the definitions (3.17a,b), the MDO $(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla })^{(n+2)}_0$ defined in (3.49) can be explicitly written out as follows:

(3.65)\begin{equation} (\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla})^{(n+2)}_0=\varPhi_0' \left( \partial_\ell -\frac{n+2}{2}\frac{\varPhi_0''}{\varPhi_0'}+\varOmega_0 \mathcal{L}_1 + \frac{\eta'}{2}\mathcal{L}_2 \right). \end{equation}

Here,

(3.66)\begin{equation} \varOmega_0 = \frac{u'}{\cosh\eta} , \quad u'=\delta-\tau, \end{equation}

and the operators $\mathcal {L}_1,\mathcal {L}_2$ are given by

(3.67a,b)\begin{align} \mathcal{L}_1=\sinh\eta \left( \cos{2u}\:\partial_\omega + (n+2)\sin{2u}\right), \quad \mathcal{L}_2= \left( \sin{2u}\:\partial_\omega - (n+2)\cos{2u}\right). \end{align}

Using (3.65) to solve the MDE (3.62), we note that the $\varPhi _0''$ term in (3.65) can be removed by defining $Y^{(n+2)}$ such that

(3.68)\begin{equation} \left. \begin{aligned} \psi^{(n+2)} & =\left( \varPhi_0'\right)^{({n+2})/{2}} Y^{(n+2)}, \\ Y^{(n+2)} & = \sum_{m=0,1}^{n}Y^{(n+2)}_{c(m+2)}\cos{(m+2)u}+ Y^{(n+2)}_{s(m+2)}\sin{(m+2)u}. \end{aligned} \right\} \end{equation}

The lower limit on the sum is zero (one) for even (odd) orders. It is straightforward to show that the operators $\mathcal {L}_1,\mathcal {L}_2$ , in general, couple the $m{\textrm {th}}$ harmonics to $m\pm 2$ harmonics. This leads to a complicated set of coupled linear ODEs for the variables $Y^{(n+2)}_{c m},Y^{(n+2)}_{s m}$. However, $\mathcal {L}_1,\mathcal {L}_2$ acting on the $(n+2)$ harmonics of $\psi ^{(n+2)}$ yields only harmonics of order $n$. From (3.45) we can deduce that the free-functions $\varPhi ^{(n+2)}_{c(n+2)}, \varPhi ^{(n+2)}_{s(n+2)}$ also appear with the $(n+2)$ harmonics. We shall now focus only on the $\varPhi ^{(n+2)}_{c(n+2)}, \varPhi ^{(n+2)}_{s(n+2)}$ equations. Substituting the forms of $(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla })^{(n+2)}_0$ and $\psi ^{(n+2)}$ , (3.65) and (3.68) we find that

(3.69)\begin{align} (\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla})^{(n+2)}_0 Y^{(n+2)} = \left( \varPhi_0' \right)^{({n+2})/{2}}\, \varPhi_0'\frac{\partial}{\partial \ell}\left( Y^{(n+2)}_{c(n+2)} \cos{(n+2)u}+Y^{(n+2)}_{s(n+2)}\sin{(n+2)u} \right)+\cdots, \end{align}

where ‘$\cdots$’ represent the rest of the terms including the $(n+2)$ harmonics that are generated by $\mathcal {L}_1,\mathcal {L}_2$ acting on $\varPhi ^{(n+2)}_{c n}, \varPhi ^{(n+2)}_{s n}$. The forcing terms with $(n+2)$ harmonics are of the form (3.64). We shall now make a choice for the free-functions $\varPhi ^{(n+2)}_{c(n+2)}, \varPhi ^{(n+2)}_{s(n+2)}$ that would allow us to solve the MDE (3.62) partially to obtain $Y^{(n+2)}_{c(n+2)}, Y^{(n+2)}_{s(n+2)}$. (The solution is partial because the coupling to $\varPhi ^{(n+2)}_{c n}, \varPhi ^{(n+2)}_{s n}$ is not included.)

Motivated by the particular choice for the free functions in Solov'ev & Shafranov (Reference Solov'ev and Shafranov1970), we now choose the free function as

(3.70)\begin{align} \varPhi^{(n+2)}_{\text{free}}& \equiv \left( \varPhi^{(n+2)}_{c(n+2)}\cos{(n+2)u}+ \varPhi^{(n+2)}_{c(n+2)}\sin{(n+2)u} \right)\nonumber\\ & \quad ={-}\frac{\left( \varPhi_0'(\ell)\right)^{(n+2)/2}}{2a(\ell)(n+2)^2}\frac{\partial}{\partial \ell}\left( Q^{(n+2)}(\ell)\cos{(n+2)u}+P^{(n+2)}(\ell)\sin{(n+2)u} \right). \end{align}

Note that our choice differs from Solov'ev & Shafranov (Reference Solov'ev and Shafranov1970) by an extra factor of ${\varPhi _0'}^{n/2}/2a(\ell )$.

The MDE for $\psi ^{(n+2)}$ (3.62) upon using the expressions (3.69) and (3.70) then leads to

(3.71a,b)\begin{equation} Y^{(n+2)}_{s(n+2)}=\frac{1}{n+2}P^{(n+2)}+\cdots,\quad Y^{(n+2)}_{c(n+2)}=\frac{1}{n+2}Q^{(n+2)}+\cdots, \end{equation}

where ‘$\cdots$’ represent the terms obtained from $Y^{(n+2)}_{c n},Y^{(n+2)}_{s n}$ that we have ignored. Thus, at each order, $(n+2)$, the choice of free function (3.70) corresponds to a choice of the $(n+2){\textrm {th}}$ poloidal of the flux-surface shape.

3.6.4 Solution of Step II: decoupling the ODEs obtained from the MDE for $\psi ^{(n+2)}$

In the last section, we showed how the MDE for $\psi ^{(n+2)}$ can be expanded in a finite number of poloidal Fourier harmonics, given in (3.68). Collecting the various harmonics, the MDE for $\psi ^{(n+2)}$, which is a PDE, can be reduced to a coupled set of ODEs for the amplitudes, $Y^{(n+2)}_{c m},Y^{(n+2)}_{s m}$, where $0\leq m \leq n+2$. The coupling between $m$ and $m\pm 2$ harmonics can be seen from (3.56) to involve the quantities $f^{(2)}_{c2},f^{(2)}_{s2}$. We now address how to decouple these ODEs. To do so, we need to discuss the leading-order normal form representation of $\psi ^{(n+2)}$ instead of the poloidal harmonic representation.

The poloidal harmonic structure of $\varPhi ^{(n+2)}$ and $\psi ^{(n+2)}$ are the same order by order, both having a maximum poloidal frequency of $(n+2)$ at order $(n+2)$. This observation is directly related to the ‘analyticity’ (in the sense of regularity near the axis) of vacuum magnetic fields assumed by Mercier and others. Rigorous proof of the analyticity condition for $\varPhi,\psi$ can be found in Duignan & Meiss (Reference Duignan and Meiss2021) and Jorge et al. (Reference Jorge, Sengupta and Landreman2020b). More concretely, ‘analyticity’ implies (Duignan & Meiss Reference Duignan and Meiss2021; Solov'ev & Shafranov Reference Solov'ev and Shafranov1970)

(3.72)\begin{align} \psi^{(n+2)}\rho^{n+2}=\mathcal{P}_{(n+2)}\left( X_\mathcal{N},Y_\mathcal{N} \right), \end{align}

where, $\mathcal {P}_{(n+2)}$ is a homogeneous polynomial of order $(n+2)$ in $(X_\mathcal {N},Y_\mathcal {N})$ with coefficients that are functions of $\ell$. Here, the leading-order normal form variables $(X_\mathcal {N},Y_\mathcal {N})$ satisfy (3.19a,b). At $O(n)$, the polynomial $\mathcal {P}_{(n+2)}$ can have at most $(n+2)$ poloidal harmonics since the leading-order normal form variables only have first harmonics.

In the following, we use the following variables obtained from the leading-order normal form variables $(X_\mathcal {N},Y_\mathcal {N})$:

(3.73a,b)\begin{equation} x_\mathcal{N}\equiv\frac{X_\mathcal{N}}{\rho}= {\rm e}^{\eta/2}\cos u, \quad y_\mathcal{N}\equiv\frac{Y_\mathcal{N}}{\rho}={\rm e}^{-\eta/2}\sin u. \end{equation}

From the expression of $\psi ^{(2)}$ in the leading-order normal form variables (3.20), we find that

(3.74a--c)\begin{equation} x_\mathcal{N}= \sqrt{\frac{2\psi^{(2)}}{\varPhi_0'}}\cos\varTheta, \quad y_\mathcal{N}= \sqrt{\frac{2\psi^{(2)}}{\varPhi_0'}}\sin\varTheta, \quad \tan\varTheta = {\rm e}^{-\eta}\tan u. \end{equation}

Following Duignan & Meiss (Reference Duignan and Meiss2021), we now employ complex variables

(3.75a,b)\begin{equation} z_\mathcal{N}=x_\mathcal{N} + {\rm i}\, y_\mathcal{N} = \sqrt{\frac{2\psi^{(2)}}{\varPhi_0'}}{\rm e}^{{\rm i}\varTheta}, \quad \bar{z}_\mathcal{N}=x_\mathcal{N}-{\rm i}\, y_\mathcal{N} = \sqrt{\frac{2\psi^{(2)}}{\varPhi_0'}}{\rm e}^{-{\rm i}\varTheta}, \end{equation}

to rewrite the expression of $\psi ^{(n+2)}$ in the complex leading-order normal form coordinates (3.72) as

(3.76)\begin{equation} \psi^{(n+2)}= \sum_{m=0}^{n+2} \varPsi^{(n+2)}_{m (n+2-m)}(\ell) \bar{z}_\mathcal{N}^m {z_\mathcal{N}}^{(n+2-m)}. \end{equation}

The connection between the amplitudes $\varPsi ^{(n+2)}_{m r}(\ell )$ and the amplitudes $\psi ^{(n+2)}_{cm},\psi ^{(n+2)}_{sm}$ follows from the definitions (3.73a,b) and (3.75a,b). Note that the reality condition of $\psi ^{(n+2)}$ implies that $\varPsi ^{(n+2)}_{m (n+2-m)}= {\bar {\varPsi }}^{(n+2)}_{(n+2-m) m}$.

We now substitute the complex leading-order normal form expression for $\psi ^{(n+2)}$ (3.76) into the MDE for $\psi ^{(n+2)}$ (3.62). For that we first need to evaluate $(\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla })^{(n+2)}_0 \psi ^{(n+2)}$.

In Appendix C, we prove the following identities:

(3.77a)\begin{gather} (\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla})^{(n+2)}_0 \zeta_{m\mathcal{N}}^{(n+2)} = \varPhi_0' \left( {\rm i}\, \varOmega_0 (n+2-2m)-\frac{n+2}{2}\frac{\varPhi_0''}{\varPhi_0'}\right) \zeta_{m\mathcal{N}}^{(n+2)},\quad \zeta_{m\mathcal{N}}^{(n+2)}=\bar{z}_\mathcal{N}^m\: {z}^{n+2-m}_\mathcal{N}, \end{gather}

(3.77b)\begin{gather}(\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla})^{(n+2)}_0 \left( \mathcal{Z}_m(\ell)\zeta_{m\mathcal{N}}^{(n+2)} \right)= \varPhi_0' \left( \mathcal{Z}_m'+ \left( {\rm i}\, \varOmega_0 (n+2-2m)-\frac{n+2}{2}\frac{\varPhi_0''}{\varPhi_0'}\right) \mathcal{Z}_m\right)\zeta_{m\mathcal{N}}^{(n+2)} , \end{gather}

which show that the complex leading-order normal form basis diagonalizes the MDO. As before, we get rid of the $\varPhi _0''$ term by defining $Y^{(n+2)}={\varPhi _0'}^{-(n+2)/2}\psi ^{(n+2)}$, equivalently by $\mathcal {Z}_m\to \mathcal {Z}_m {\varPhi _0'}^{(n+2)/2}$ in (3.77).

Therefore, the substitution of the following form of $\psi ^{(n+2)}$:

(3.78)\begin{equation} \psi^{(n+2)}=\left(\varPhi_0' \right)^{({n+2})/{2}} \sum_{m=0}^{n+2} \mathcal{Z}^{(n+2)}_{m (n+2-m)}(\ell) \bar{z}_\mathcal{N}^m {z_\mathcal{N}}^{(n+2-m)}, \quad \mathcal{Z}^{(n+2)}_{m (n+2-m)}={\bar{\mathcal{Z}}}^{(n+2)}_{(n+2-m) m}, \end{equation}

into the MDE (3.62) leads to a system of $(n+1)$ decoupled complex ODEs of the form

(3.79)\begin{equation} \dfrac{{\rm d}}{{\rm d}\ell} {\mathcal{Z}_{m (n+2-m)}^{(n+2)}} + {\rm i}\, \varOmega_0 (n+2-2m) \mathcal{Z}^{(n+2)}_{m (n+2-m)}+ \mathcal{F}^{(n+2)}_{\mathcal{Z}_m}=0, \end{equation}

where $\mathcal {F}^{(n+2)}_{\mathcal {Z}_m}$ is an appropriate forcing term. We refer to (3.79) as a reduced MDE. The above results also holds for $m=0,(n+2)=1$ as shown in Appendix C.

The complex leading-order normal form (3.72) for $\psi ^{(n+2)}$ is, therefore, of great help in systematically solving the MDE for $\psi$ by reducing it to a system of decoupled ODEs (3.79). In Appendix F we have shown explicitly how the MDE for $\psi ^{(3)}$ can be solved by leveraging the leading-order normal form structure.

3.6.5 Solutions of Steps III and IV: determination of $\alpha$

Our main goal here is to calculate $\alpha ^{(n)}$ given the solutions of $\varPhi ^{(n+2)},\psi ^{(n+2)}$, combining Steps III and IV, outlined in § 3.4. Step III consists of integrating the equation for $\alpha _{,\omega }$ with respect to $\omega$ to determine $\alpha ^{(n)}$ up to a function $\mathfrak {a}(\ell )$. Step IV involves poloidally averaging the MDE for $\alpha ^{(n)}$ to get an equation for $\mathfrak {a}(\ell )$.

We start with Step III. A power series expansion of (3.2c), rewritten as

(3.2c)\begin{align} \psi_{,\rho}\alpha_{,\omega}-\psi_{,\omega}\alpha_{,\rho}-\frac{\rho}{h}\varPhi_{,\ell}=0, \end{align}

at each order, $n$, leads to

(3.80)\begin{equation} \sum_{m=0}^n \left( (n-m) \psi^{(m+2)}_{,\omega}\alpha^{(n-m)}- (m+2) \psi^{(m+2)}\alpha^{(n-m)}_{,\omega} +\left( \kappa \cos \theta \right)^{(n-m)} \varPhi^{(m)}_{,\ell} \right) =0. \end{equation}

The function $1/h=(1-\kappa \cos \theta )^{-1}$ has been expanded in powers of $\rho$ as well.

We shall now split $\psi ^{(n+2)}$ into the homogeneous $\psi ^{(n+2)}_H$ solution and the particular solution $\psi ^{(n+2)}_P$,

(3.81a,b)\begin{equation} \psi^{(n+2)}= \psi^{(n+2)}_P+{\psi_H}^{(n+2)}, \quad {\psi_H}^{(n+2)}= \left( 2\psi^{(2)}\right)^{n/2+1}\mathcal{Y}^{(n+2)}_H. \end{equation}

Here, we assume that the on-axis rotational transform is irrational, so $\mathcal {Y}^{(n+2)}_H$ is a constant. Separating the $m=0$ and $m=n$ terms in (3.80) from the rest and using the lowest-order relation (3.11) we obtain

(3.82)\begin{equation} \partial_\omega\left(\frac{\alpha^{(n)}}{\left( 2\psi^{(2)}\right)^{n/2}} \right)= \frac{\varPhi^{(n)}_{,\ell} +\varSigma^{(n)}_1}{\left( 2\psi^{(2)}\right)^{n/2+1}}-(n+2)\left( \frac{\varPhi_0'\psi^{(n+2)}_P}{\left( 2\psi^{(2)}\right)^{n/2+2}} +\alpha_{,\omega}^{(0)} \mathcal{Y}^{(n+2)}_H\right), \end{equation}

where $\left ( \varSigma ^{(n)}_1-(\kappa \cos \theta )^n\varPhi _0'\right )$ is the sum of terms in (3.80) from $m=1$ to $(n-1)$. Imposing the poloidal integral constraint on $\alpha ^{(n)}$ as given in (3.42) and equating the average of $\alpha _{,\omega }^{(0)}$ to unity, as shown in (3.35), we find that

(3.83)\begin{equation} \mathcal{Y}^{(n+2)}_H= \oint\frac{{\rm d}\omega}{\left( 2\psi^{(2)}\right)^{n/2+1}}\frac{\left( \varPhi^{(n)}_{,\ell} +\varSigma^{(n)}_1\right)}{(n+2)} -\varPhi_0'\oint\frac{d\omega}{\left( 2\psi^{(2)}\right)^{n/2+2}}\psi^{(n+2)}_P. \end{equation}

We can now integrate (3.82) with respect to $\omega$ to obtain

(3.84)\begin{equation} \left(\frac{\alpha^{(n)}}{\left( 2\psi^{(2)}\right)^{n/2}} \right)=\bar{\mathfrak{a}}^{(n)}(\ell)+\tilde{\mathfrak{a}}^{(n)}(\omega,\ell), \end{equation}

where,

(3.85)\begin{align} \tilde{\mathfrak{a}}^{(n)}(\omega,\ell)= \int_0^\omega d\omega' \left( \frac{\varPhi^{(n)}_{,\ell} +\varSigma^{(n)}_1}{\left( 2\psi^{(2)}\right)^{n/2+1}}-(n+2)\varPhi_0'\frac{\psi^{(n+2)}_P+\psi^{(n+2)}_H}{\left( 2\psi^{(2)}\right)^{n/2+2}}\right), \quad \oint \frac{d\omega'}{2{\rm \pi}} \tilde{\mathfrak{a}}(\omega',\ell)=0. \end{align}

Note that the average of $\tilde {\mathfrak {a}}$ is zero because of the condition (3.83) that determines the homogeneous solution $\psi ^{(n+2)}_H$. This completes Step III.

We now proceed with Step IV. To determine $\bar {\mathfrak {a}}^{(n)}(\ell )$, we use the Poisson bracket form of the MDO (3.50)

(3.86)\begin{equation} \left( \alpha_{,\omega}^{(0)}\partial_\ell -\alpha_{,\ell}^{(0)}\partial_\omega\right) \left( \frac{\alpha^{(n)}}{\left( 2\psi^{(2)} \right)^{n/2}}\right) +\frac{F^{(n+2)}_\alpha}{\left( 2\psi^{(2)}\right)^{n/2+1}}=0, \end{equation}

which was derived earlier in (3.27).

Writing (3.86) in the form (3.52) and integrating with respect to $\omega$ between $(0,2{\rm \pi} )$ we obtain the poloidally averaged MDE

(3.87)\begin{equation} \partial_\ell \left( \oint \frac{{\rm d}\omega}{2{\rm \pi}} \alpha_{,\omega}^{(0)} \left( \frac{\alpha^{(n)}}{\left( 2\psi^{(2)} \right)^{n/2}}\right) \right) + \oint \frac{{\rm d}\omega}{2{\rm \pi}} \left( \frac{F^{(n+2)}_\alpha}{\left( 2\psi^{(2)}\right)^{n/2+1}} \right)=0. \end{equation}

We note here that the derivative in $\omega$ vanished because both $\alpha _{,\omega }^{(0)}$ and $\left ( \alpha ^{(n)}/\left ( 2\psi ^{(2)}\right )^{n/2}\right )$ are single-valued in $\omega$.

From (3.84) and (3.35) we then obtain the following equation for $\bar {\mathfrak {a}}^{(n)}(\ell )$:

(3.88)\begin{equation} {\bar{\mathfrak{a}}^{(n)}}{}'(\ell) + \partial_\ell \oint \frac{{\rm d}\omega}{2{\rm \pi}} \frac{\varPhi_0'}{2\psi^{(2)}}\tilde{\mathfrak{a}}^{(n)}(\omega,\ell) + \oint \frac{{\rm d}\omega}{2{\rm \pi}} \left( \frac{F^{(n+2)}_\alpha}{\left( 2\psi^{(2)}\right)^{n/2+1}} \right)=0. \end{equation}

The solution of (3.88) can be obtained through direct integration with respect to $\ell$. Note that ${\bar {\mathfrak {a}}^{(n)}}{}'(\ell )$ is not periodic in general since its average is related to the higher-order magnetic-shear as given in (3.42).

We reiterate that the form of $\alpha$ in inverse coordinates, as given in (3.28a,b), is well known. In particular, the secular terms $\theta -\iota (\psi )\phi$ remain separate from the periodic components of $\alpha$. However, in direct coordinates, the secular and the periodic terms are mixed because the $\iota (\psi )\phi$ term is now expanded in a power series of $\rho$. This explains the need to take special care in order to obtain $\alpha$ in direct coordinates.

4.1 Basic formalism for the force-free fields

In the Mercier–Weitzner formalism, force-free fields that possess nested flux surfaces can be seen to be of the form

(4.4)\begin{equation} \boldsymbol{B} =\boldsymbol{\nabla} \varPhi +\bar{K}(\psi,\alpha) \boldsymbol{\nabla} \psi, \quad \boldsymbol{B} =\boldsymbol{\nabla}\psi\times \boldsymbol{\nabla} \alpha, \end{equation}

since the current takes the form

(4.5)\begin{equation} \boldsymbol{J} = \boldsymbol{\nabla} \bar{K}\times \boldsymbol{\nabla} \psi ={-}\bar{K}_{,\alpha}\boldsymbol{B}. \end{equation}

Comparing (4.1) with (4.5) we find that

(4.6)\begin{equation} \bar{K} ={-}\lambda(\psi) \alpha. \end{equation}

Note that a homogeneous solution $\bar {K}'_0(\psi )$ could be added to (4.6). However, from (4.4) we find that the $\bar {K}'_0(\psi )$ term can be absorbed by redefining $\varPhi \to \varPhi -\bar {K}_0(\psi )$. Hence, we can set $\bar {K}_0$ to zero without loss of generality.

To carry out the NAE, we assume that $\lambda (\psi )$ is sufficiently smooth that we can expand $\lambda$ as

(4.7)\begin{equation} \lambda(\psi)=\lambda_0 + \psi \lambda_2 + \psi^2 \lambda_4+\cdots, \end{equation}

where, $\lambda _i, i=0,2,4,\ldots$ are constants. However, the usual NAE in direct coordinates (3.7) needs to be checked carefully. From the condition $\boldsymbol {\nabla }\boldsymbol {\cdot } \boldsymbol {B}=0$ we have

(4.8)\begin{equation} \Delta \varPhi + \boldsymbol{\nabla}\boldsymbol{\cdot} (\bar{K} \boldsymbol{\nabla} \psi)=0. \end{equation}

We need to show that the solution of (4.8) can be expressed as a regular power series near the axis. In the vacuum limit, we obtained (3.44), where the right-hand side did not have any harmonics that would beat with the operator on the left-hand side. However, it is not guaranteed that at each $O(\rho ^{n})$, the term $\boldsymbol {\nabla }\boldsymbol {\cdot } (\bar {K} \boldsymbol {\nabla } \psi )$ will not generate such resonant harmonics.

To further illustrate this point, let us consider the case when $\lambda =$ constant, i.e. $K=-\lambda \alpha$. The $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B}=0$ condition now reads

(4.9)\begin{equation} \Delta \varPhi -\lambda \boldsymbol{\nabla}\boldsymbol{\cdot}( \alpha \boldsymbol{\nabla} \psi)=0. \end{equation}

Expanding as in (3.7), we get

(4.10)\begin{align} \left. \begin{aligned} \Delta \varPhi & = \sum_n \rho^{n}\left( \partial^2_{\omega}+(n+2)^2\right)\varPhi^{(n+2)}+ \cdots ,\\ \boldsymbol{\nabla}\boldsymbol{\cdot}( \alpha \boldsymbol{\nabla} \psi) & = \sum_n \rho^{n}\sum_{i=0}^{n}\left( \alpha^{(i)}\psi^{(n+2-i)}(n+2)(n+2-i)+\partial_\omega\left( \alpha^{(i)}\psi_{,\omega}^{(n+2-i)}\right)\right) +\cdots . \end{aligned} \right\} \end{align}

Using the Fourier representations, $\psi ^{(n+2-i)}\sim \sum \psi ^{(n+2-i)}_m \textrm {e}^{\textrm {i} m u},\alpha ^{(i)}\sim \sum \alpha ^{(i)}_p \textrm {e}^{\textrm {i} p u}$,

(4.11)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}( \alpha \boldsymbol{\nabla} \psi) \sim \sum_n\rho^{n}\sum_{i=0}^{n+2}\sum_{m=0}^{n+2-i}\sum_{p=0}^{i}\left(\psi^{(n+2-i)}_m\alpha^{(i)}_p {\rm e}^{{\rm i} (m+p)u}c_{ni}+\cdots\right), \end{equation}

where, $c_{ni}(\ell )$ is some overall constant factor. Clearly, resonance can occur when $m+p=n+2$. If such a $(n+2)$ harmonic is generated then we cannot invert

(4.12)\begin{equation} (\partial^2_{\omega}+(n+2)^2)\varPhi^{(n+2)}\sim {\rm e}^{ {\rm i} (n+2)u}, \end{equation}

to solve for a periodic $\varPhi ^{(n+2)}$. As shown in Weitzner (Reference Weitzner2016), if such a resonance occurs, then the regular power series will fail, and weakly singular terms like $\rho ^n(\log {\rho })^m$ will appear in the Mercier expansion.

In the following section, we prove that such logarithmically singular terms can indeed be avoided under reasonable assumptions. However, the boundary cannot be arbitrarily specified and must be self-consistent with the NAE.

4.2 Definition of new variables

Equating $\boldsymbol {B} = \boldsymbol {\nabla }\varPhi -\lambda \alpha \boldsymbol {\nabla } \psi = \boldsymbol {\nabla } \psi \times \boldsymbol {\nabla } \alpha$, we get the following system of equations for the force-free magnetic fields with flux surfaces in the Mercier–Weitzner formalism:

(4.13a)\begin{gather} \varPhi_{,\rho}-\lambda \alpha \psi_{,\rho}= \frac{1}{\rho h}\left( \psi_{,\omega}\alpha_{,\ell}- \psi_{,\ell}\alpha_{,\omega}\right), \end{gather}

(4.13b)\begin{gather}\varPhi_{,\omega}-\lambda \alpha \psi_{,\omega}= \frac{\rho}{h}\left( \psi_{,\ell}\alpha_{,\rho}- \psi_{,\rho}\alpha_{,\ell}\right), \end{gather}

(4.13c)\begin{gather}\varPhi_{,\ell}-\lambda \alpha \psi_{,\ell}= \frac{h}{\rho }\left( \psi_{,\rho}\alpha_{,\omega}- \psi_{,\omega}\alpha_{,\rho}\right). \end{gather}

As discussed in the last section, unlike the vacuum limit, it is not immediately evident that resonances can be avoided in the system (4.13).

Our goal is, therefore, to bring the force-free system of (4.13) as close as possible to the vacuum system given by (3.2). We now try to define a new variable $\xi$ that will play the role analogous to $\varPhi$ in the vacuum limit. It is useful to define a variable $\varLambda$ such that

(4.14)\begin{equation} \lambda(\psi)= \varLambda'(\psi). \end{equation}

Now, we define $\xi$ such that

(4.15)\begin{equation} \xi_{,\rho}\equiv \varPhi_{,\rho}-\lambda \alpha \psi_{,\rho}= \varPhi_{,\rho}-\alpha \varLambda_{,\rho}. \end{equation}

Integrating with respect to $\rho$ we get

(4.16)\begin{equation} \xi = \varPhi + \mathcal{I},\quad \mathcal{I}\equiv{-}\int_0^\rho d\varrho\: \varLambda_{,\varrho} \alpha. \end{equation}

Note that $\xi,\varPhi$ and $\alpha$ have similar secular terms linear in the angles. Therefore, imposing the form (2.11) order by order on $\alpha$ is enough to impose the right secular structure of $\xi$.

Since, ultimately, we hope to be able to expand in power series of $\rho$, the integral $\mathcal {I}$ in (4.16), amounts to simply reshuffling the terms and multiplying the coefficients by fractions. With this definition in mind we now rewrite (4.13) in terms of $\xi$ and its derivatives as

(4.17a)\begin{gather} \xi_{,\rho}= \frac{1}{\rho h}\left( \psi_{,\omega}\alpha_{,\ell}- \psi_{,\ell}\alpha_{,\omega}\right), \end{gather}

(4.17b)\begin{gather}\xi_{,\omega}-\varLambda_{,\omega} \alpha-\mathcal{I}_{,\omega} = \frac{\rho}{h}\left( \psi_{,\ell}\alpha_{,\rho}- \psi_{,\rho}\alpha_{,\ell}\right), \end{gather}

(4.17c)\begin{gather}\xi_{,\ell}-\varLambda_{,\ell} \alpha-\mathcal{I}_{,\ell}= \frac{h}{\rho }\left( \psi_{,\rho}\alpha_{,\omega}- \psi_{,\omega}\alpha_{,\rho}\right). \end{gather}

The $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B}=0$ condition in terms of $\xi$ reads

(4.18)\begin{equation} \rho^2 \Delta \xi -\frac{1}{h}\partial_\omega \left( h\left(\varLambda_{,\omega} \alpha+\mathcal{I}_{,\omega} \right) \right) -\frac{1}{h}\partial_\ell \left( \frac{\rho^2}{h}\left( \varLambda_{,\ell} \alpha+\mathcal{I}_{,\ell}\right)\right)=0. \end{equation}

Comparing (4.17), (4.18) with the vacuum limit given by (3.2) and (3.4), we find that the two differ because of the $(\varLambda _{,\omega } \alpha +\mathcal {I}_{,\omega })$ term, which is an $O(1)$ term in the $\rho$ expansion. All other deviations from the vacuum limit are higher order in $\rho$. Fortunately, a workaround is possible by noting that

(4.19)\begin{equation} \partial_\rho\left( \varLambda \alpha_{,\omega}\right)- \partial_\omega\left( \varLambda \alpha_{,\rho}\right) = \lambda \left( \psi_{,\rho}\alpha_{,\omega}-\alpha_{,\rho}\psi_{,\omega} \right). \end{equation}

Integrating once again in $\rho$, using (4.17c) and simplifying, we get

(4.20)\begin{equation} {-}\varLambda_{,\omega} \alpha-\mathcal{I}_{,\omega} =\mathcal{W},\quad \mathcal{W}\equiv\int_0^\rho \,{\rm d}\varrho \frac{\varrho}{h}\lambda \left( \xi_{,\ell}+\varXi\right), \quad \varXi\equiv{-}\left( \varLambda_{,\ell}\alpha + \mathcal{I}_{,\ell}\right). \end{equation}

With these new definitions, we can write down the force-free equations as

(4.21a)\begin{gather} \xi_{,\rho}= \frac{1}{\rho h}\left( \psi_{,\omega}\alpha_{,\ell}- \psi_{,\ell}\alpha_{,\omega}\right)=B_\rho, \end{gather}

(4.21b)\begin{gather}\xi_{,\omega}+\mathcal{W} = \frac{\rho}{h}\left( \psi_{,\ell}\alpha_{,\rho}- \psi_{,\rho}\alpha_{,\ell}\right)= \rho B_\omega , \end{gather}

(4.21c)\begin{gather}\xi_{,\ell}+\varXi= \frac{h}{\rho }\left( \psi_{,\rho}\alpha_{,\omega}- \psi_{,\omega}\alpha_{,\rho}\right) , = h B_\ell \end{gather}

where

(4.22)\begin{align} \lambda=\varLambda'(\psi),\quad \mathcal{I}\equiv{-} \int_0^\rho \,{\rm d}\varrho \varLambda_{,\varrho} \alpha,\quad \varXi\equiv{-}(\varLambda_{,\ell}\alpha + \mathcal{I}_{,\ell}), \quad \mathcal{W}\equiv \int_0^\rho \,{\rm d}\varrho \frac{\varrho}{h}\lambda \left( \xi_{,\ell}+\varXi\right). \end{align}

Finally, the Poisson equation for $\xi$ reads

(4.23)\begin{equation} \rho^2\Delta \xi +\frac{1}{h}\partial_\omega\left( h \mathcal{W} \right) +\frac{1}{h}\partial_\ell\left( \frac{\rho^2}{h}\varXi\right)=0. \end{equation}

In Appendix E, we show in detail that just like in the vacuum limit, (4.23) can be solved order by order without encountering any resonances. In short, the basic idea of the proof is to note that while the first term in (4.23) is $O(1)$ and of the form $(\partial ^2_{\omega }+(n+2)^2)\xi ^{(n+2)}$, the second and the third terms both are $O(\rho ^2)$ smaller than the first term. Hence, the second and third terms have at most $n\textrm {th}$ harmonics, thereby avoiding any resonance with the first term. Thus, the logarithmic terms are avoided. Once $\xi$ has been obtained to $O(\rho ^n)$, we can use (4.16) to determine $\varPhi$ to the same order of accuracy.

4.3 The NAE equations for force-free equilibrium

We shall now discuss the lowest-order NAE equations for the force-free equilibrium using the variables defined in the last section. Finally, we shall discuss the lowest two orders before discussing the general $O(\rho ^n)$ system of equations.

Substituting the near-axis power series expansion into the definitions of $\mathcal {I},\varXi$ given in (4.22), we find that

(4.24)\begin{align} \mathcal{I} ={-}\lambda_0 \rho^2 \alpha^{(0)} \psi^{(2)} +O(\rho^3), \quad \varXi=\rho^2 \lambda_0 \psi^{(2)} \alpha_{,\ell}^{(0)} +O(\rho^3), \quad \mathcal{W}=\frac{\lambda_0}{2}\xi^{(0)}_{,\ell}\rho^2 +O(\rho^3). \end{align}

We are now in a position to investigate the Poisson equation for $\xi$, (4.23). First, we note from (4.24) that the contribution of $\varXi$ to (4.21b) starts at $O(\rho ^4)$. Next, we find from (4.24) that the Poisson equation (4.23) matches the vacuum case exactly to $O(\rho ^{-2}), O(\rho ^{-1})$. Thus, $\xi$ and $\varPhi$ differ only from the second-order onward, i.e.

(4.25)\begin{equation} \xi^{(0)}=\varPhi_0(\ell), \quad \xi^{(1)}=0, \quad B_0=\varPhi_0'(\ell). \end{equation}

Consequently, we note that even though $\mathcal {W}\sim O(\rho ^2)$ its contribution to (4.23) through the $\partial \omega (h \mathcal {W})$ term is $O(\rho ^3)$ since it is a function of $\ell$ to lowest order. Therefore, even the $O(1)$ Poisson equation assumes the same form as in the vacuum case, namely,

(4.26)\begin{equation} \left( \partial^2_{\omega} +4 \right)\xi^{(2)} +\varPhi_0''(\ell)=0, \end{equation}

whose solution is given by

(4.27)\begin{equation} \left. \begin{aligned} \xi^{(2)} & ={-}\dfrac{1}{4}\varPhi_0''(\ell)+\varPhi_0'(z)\left(f^{(2)}_{c2}(\ell)\cos{2u}+f^{(2)}_{s2}(\ell)\sin{2u}\right),\\ u & =\theta +\delta(\ell)=\omega-\int\tau d \ell +\delta(\ell). \end{aligned} \right\} \end{equation}

The components of the magnetic field are given by

(4.28)\begin{align} B_\rho=2\rho \xi^{(2)} + O(\rho^3), \quad B_\omega = \rho \left( \xi^{(2)}_{,\omega} + \frac{\lambda_0}{2} \varPhi_0' \right), \quad B_\ell = (1+\kappa \rho \cos\theta)\varPhi_0' + O(\rho^2), \end{align}

which implies that

(4.29)\begin{equation} (\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla})_0^{(n)}= \varPhi_0' \partial_\ell + \left( \xi^{(2)}_{,\omega} + \frac{\lambda_0}{2} \varPhi_0' \right)\partial_\omega +2 n \xi^{(2)}. \end{equation}

Another implication of (4.28) is that the $O(\rho ^2)$ NAE equations are

(4.30a)\begin{gather} 2\psi^{(2)} \alpha_{,\omega}^{(0)} =\varPhi_0', \quad 2\psi^{(2)} \alpha_{,\ell}^{(0)} ={-}\left( \xi^{(2)}_{,\omega} +\frac{\lambda_0}{2} \varPhi_0'\right) , \end{gather}

(4.30b)\begin{gather}2\xi^{(2)} =\alpha_{,\ell}^{(0)} \psi^{(2)}_{,\omega}-\alpha_{,\omega}^{(0)} \psi^{(2)}_{,\ell}. \end{gather}

The lowest-order flux surface and field line label can be shown to be the one obtained by Mercier (Mercier Reference Mercier1964)

(4.31a)\begin{gather} \psi^{(2)}=\varPhi_0'(\ell)(a(\ell)+b(\ell) \cos{2u}), \end{gather}

(4.31b)\begin{gather}a(\ell)=\dfrac{1}{2} \cosh(\eta(\ell)),\quad b(\ell)=\dfrac{1}{2} \sinh(\eta(\ell)), \end{gather}

(4.31c)\begin{gather}f^{(2)}_{s2}(\ell)= \dfrac{b}{2a}\left( \delta'-\tau+\frac{\lambda_0}{2}\right)\quad ,\quad f^{(2)}_{c2}(\ell)={-}\dfrac{a'}{4b}={-}\frac{\eta'}{4}, \end{gather}

(4.31d)\begin{gather}\alpha^{(0)}= \tan^{{-}1}\left( {\rm e}^{-\eta}\tan u\right)-\int \,{\rm d}\ell\left( \dfrac{\left( \delta'-\tau+\lambda_0/2\right)}{\cosh{\eta}}\right). \end{gather}

As in the vacuum case, we choose $2{\rm \pi} \psi$ as the toroidal flux. A comparison with the vacuum solutions given in § 3.3 shows minimal changes from vacuum solutions. The main difference arises because of the extra $\lambda _0/2$ factor in the expression for $\alpha _{,\ell }^{(0)}$ in (4.30a). While the structure of the equations stays the same, the forcing functions, $F_\phi,F_\psi, F_\alpha$, get modified to include the effects of force-free currents.

The on-axis rotational transform for the force-free case is given by

(4.32)\begin{equation} \left. \begin{aligned} \iota_0 & = \alpha(\ell,\theta)-\alpha(\ell+L,\theta)+N = \dfrac{v(L)-(\delta(L)-\delta(0))}{2{\rm \pi}}+N,\\ v(\ell) & = \int_0^\ell \,{\rm d}s \dfrac{\delta'(s)-\tau(s)+\dfrac{\lambda_0}{2}}{\cosh{\eta(s)}}. \end{aligned} \right\} \end{equation}

As shown by Mercier (Reference Mercier1974), pressure does not affect the on-axis rotational transform. Therefore, the on-axis rotational transform given in (4.32) also applies to the general MHS case.

We now proceed to describe the structure of the $O(\rho ^{(n+2)})$ equations. We note that just like in the vacuum case, (4.21a) gives us an MDE for $\psi ^{(n+2)}$. To obtain the MDE for $\alpha ^{(n)}$ we proceed exactly as in the vacuum limit by eliminating $\psi ^{(n+2)}$ from (4.21b) and (4.21c). Therefore, following the vacuum case, we can write down the steps involved in solving the force-free equations $(\xi ^{(n+2)},\psi ^{(n+2)},\alpha ^{(n)})$. The steps are precisely the ones we discussed in the vacuum limit:

1. (I) the Poisson equation (4.23) for $\xi ^{(n+2)}$;

2. (II) the MDE for $\psi ^{(n+2)}$ from (4.21a);

3. (III) the $\alpha _{,\omega }$ equation from (4.21c) to solve for $\alpha ^{(n)}$ up to a function $\mathfrak {\bar {a}}^{(n)}(\ell )$;

4. (IV) the poloidal $\omega$ average of the MDE for $\alpha$ to determine $\mathfrak {\bar {a}}^{(n)}(\ell )$.

We present a detailed demonstration of the above steps for $n=1$ in Appendix G.

5.1 Basic formalism for the MHS fields: current potentials

In the Mercier–Weitzner formalism, as mentioned in § 2, we have an inhomogeneous MDE of the form

(2.8)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} K =p'(\psi)\end{equation}

that must be solved to obtain the current potential $K$. In the force-free limit, we obtained the homogeneous solution $K=\bar {K}=-\lambda (\psi )\alpha$. In the presence of finite pressure gradients, in addition to $\bar {K}$, we also need to solve for the inhomogeneous solution, $G$, such that

(5.1)\begin{equation} K=\bar{K}(\psi,\alpha) +G(\rho,\omega,\ell), \quad \bar{K}(\psi,\alpha)={-}\lambda \alpha, \quad \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} G=p'(\psi). \end{equation}

Equation (5.1) implies

(5.2)\begin{equation} \boldsymbol{J}= \boldsymbol{\nabla} K\times \boldsymbol{\nabla} \psi = \lambda(\psi) \boldsymbol{B} +\boldsymbol{\nabla} G\times \boldsymbol{\nabla} \psi . \end{equation}

Thus, the quantity $G$ is the pressure-driven part of the current potential. Note that, in general, the function $\lambda$ is not the parallel current ($j_{\|}$) since

(5.3)\begin{equation} j_{\|} \equiv \frac{\boldsymbol{J} \boldsymbol{\cdot} \boldsymbol{B}}{B^2}= \lambda(\psi) -\frac{\boldsymbol{B}\times \boldsymbol{\nabla}\psi\boldsymbol{\cdot} \boldsymbol{\nabla} G}{B^2}. \end{equation}

Only in the force-free limit ($p'=0, G=0$) i $j_{\|}$ equal to $\lambda$. Also, note that $G$ is not single-valued. However, the secular parts of $G$ need to be of the same form as $K$ and $\alpha$ in (2.11).

Let us first look at the axisymmetric limit to understand better the relation between $\lambda,G$ and the parallel current. The axisymmetric magnetic field and the current in a tokamak are given by (Freidberg Reference Freidberg1982)

(5.4)\begin{equation} \boldsymbol{B} =F(\varPsi)\boldsymbol{\nabla} \zeta + \boldsymbol{\nabla} \zeta \times \boldsymbol{\nabla} \varPsi, \quad \boldsymbol{J} ={-}F'(\varPsi)\boldsymbol{B} -R^2p'(\varPsi)\boldsymbol{\nabla} \zeta, \end{equation}

where, $(R,\zeta,z)$ denotes the standard cylindrical coordinate system, $F(\varPsi )$ is the poloidal current and $\varPsi$ is the poloidal flux function that satisfies $\textrm {d}\varPsi /\textrm {d}\psi =\iota (\psi )$ together the Grad–Shafranov equation

(5.5)\begin{equation} \Delta_* \varPsi + F F' + R^2 p'(\varPsi)=0, \quad \Delta_* = R^2\boldsymbol{\nabla}\boldsymbol{\cdot} \left( \frac{1}{R^2}\boldsymbol{\nabla} \right) . \end{equation}

Furthermore, in axisymmetry, the toroidal angle dependence can only be through secular terms. Thus,

(5.6)\begin{equation} G=\bar{G}(\varPsi(R,z))\zeta + \tilde{G}(R,z), \quad \bar{G} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \zeta+\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \tilde{G}=\iota p'(\varPsi) . \end{equation}

Using the following axisymmetric relations (Freidberg Reference Freidberg1982; Helander Reference Helander2014):

(5.7)\begin{equation} j_{\|}={-}F'(\varPsi)-\frac{p'(\varPsi)F}{B^2}, \quad \boldsymbol{B}\times \boldsymbol{\nabla} \psi \boldsymbol{\cdot} \boldsymbol{\nabla} = F \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} -B^2 \partial_\zeta, \end{equation}

and (5.6), the pressure driven terms in the equation for $\lambda$, (5.3), cancel resulting in

(5.8)\begin{equation} {-}\lambda(\psi) = F'(\varPsi)+\frac{1}{\iota(\psi)}\bar{G}(\varPsi). \end{equation}

Equation (5.8) shows that in the force-free limit $-\lambda (\psi )= F'(\varPsi )$, but for finite beta, the pressure-driven current also contributes to $\lambda$ through the secular term.

5.2 The NAE of MHS fields

We now discuss the NAE of MHS fields. We postulate that $\bar {K},G$ have the following expansions:

(5.9a)\begin{gather} \bar{K}(\psi,\alpha)={-}\alpha (\lambda_0 +\lambda_2 \psi +\lambda_4 \psi^2+\cdots), \end{gather}

(5.9b)\begin{gather}G(\rho,\omega,\ell)=G^{(0)}(\ell)+G^{(1)}(\omega,\ell)\rho +G^{(2)}(\omega,\ell) \rho^2+\cdots, \end{gather}

where each quantity is finite and non-singular. By allowing power-series expansions of $\bar {K},G$, we are explicitly excluding current profiles that are singular near the axis. It is beyond the scope of this work to account for the effects of near-axis current singularities and will be left for future publications.

Carrying out the NAE assuming power series expansion in $\rho$, we get the following coupled PDE system at each order, $n$:

(5.10a)\begin{gather} (\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} G)^{(n)}=p^{(n+2)}+\cdots, \end{gather}

(5.10b)\begin{gather}(\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \psi)^{(n)}=\ldots, \end{gather}

(5.10c)\begin{gather}(\Delta \varPhi + \boldsymbol{\nabla}\boldsymbol{\cdot} (K \boldsymbol{\nabla} \psi))^{(n)}=0. \end{gather}

We have precisely two MDE for $G$ and $\psi$ and one elliptic equation for $\varPhi$. The function $\alpha$ is obtained from $\varPhi,\psi$ without solving an MDE for $\alpha$. Therefore, the PDE is of mixed elliptic–hyperbolic type as is well known (Bauer, Betancourt & Garabedian Reference Bauer, Betancourt and Garabedian2012a; Weitzner Reference Weitzner2014). However, as discussed in § 4, it is not immediately obvious that (5.10c) can be solved for $\varPhi$ without encountering logarithmic singularities. In principle, the logarithmic singularities can still arise even if the current profiles are chosen to be smooth near the axis by imposing (5.9).

5.3 Definitions of new variables

The change of variables is the same as in the force-free case, with additional terms from $G$. We shall have first the covariant and the contravariant forms of $\boldsymbol {B}$:

(5.11)\begin{equation} \boldsymbol{\nabla} \varPhi +(-\lambda \alpha +G )\boldsymbol{\nabla} \psi =\boldsymbol{\nabla} \psi\times \boldsymbol{\nabla} \alpha. \end{equation}

Then we define the new variable $\xi$, such that

(5.12)\begin{equation} \xi_{,\rho}=\varPhi_{,\rho}+(-\lambda \alpha +G)\psi_{,\rho}. \end{equation}

Integrating (5.13) with respect to $\rho$ we get

(5.13)\begin{equation} \xi = \varPhi +\mathcal{I} +\mathcal{G}, \quad \mathcal{I} \equiv{-}\int_0^\rho \,{\rm d}\varrho\;\varLambda_{,\varrho} \alpha, \quad \mathcal{G}=\int_0^\rho \,{\rm d}\varrho\;G \psi_{,\varrho}. \end{equation}

The $(\rho,\omega,\ell )$ components of (5.11) then yield

(5.14a)\begin{gather} \xi_{,\rho}= \frac{1}{\rho h}\left( \psi_{,\omega}\alpha_{,\ell}- \psi_{,\ell}\alpha_{,\omega}\right)=B_\rho , \end{gather}

(5.14b)\begin{gather}\xi_{,\omega}+\mathcal{Y} -\left( \mathcal{I}_{,\omega}+\alpha \varLambda_{,\omega}\right) = \frac{\rho}{h}\left( \psi_{,\ell}\alpha_{,\rho}- \psi_{,\rho}\alpha_{,\ell}\right)= \rho B_\omega , \end{gather}

(5.14c)\begin{gather}\xi_{,\ell}+\varXi+\mathcal{X}= \frac{h}{\rho }\left( \psi_{,\rho}\alpha_{,\omega}- \psi_{,\omega}\alpha_{,\rho}\right) = h B_\ell. \end{gather}

While the $B_\rho$ and $B_\ell$ equations above show minimal change when we go to the vacuum limit, the $B_\omega$ equation needs modifications similar to the force-free case. Utilizing the same trick, (4.19), that we used before in the force-free case, we find that

(5.15)\begin{equation} {-}\left( \mathcal{I}_{,\omega}+\alpha \varLambda_{,\omega}\right)= \int_0^\rho \,{\rm d}\varrho \frac{\varrho}{h}\lambda \left( \xi_{,\ell}+\varXi+\mathcal{X}\right). \end{equation}

Therefore, the final set of MHS equations has the form

(5.16a)\begin{gather} \xi_{,\rho}= \frac{1}{\rho h}\left( \psi_{,\omega}\alpha_{,\ell}- \psi_{,\ell}\alpha_{,\omega}\right)=B_\rho , \end{gather}

(5.16b)\begin{gather}\xi_{,\omega}+\mathcal{Y}+ \mathcal{W} = \frac{\rho}{h}\left( \psi_{,\ell}\alpha_{,\rho}- \psi_{,\rho}\alpha_{,\ell}\right)= \rho B_\omega , \end{gather}

(5.16c)\begin{gather}\xi_{,\ell}+\varXi+\mathcal{X}= \frac{h}{\rho }\left( \psi_{,\rho}\alpha_{,\omega}- \psi_{,\omega}\alpha_{,\rho}\right) = h B_\ell, \end{gather}

where

(5.16d)\begin{align} \left. \begin{aligned} \lambda & =\varLambda'(\psi),\quad \mathcal{I}\equiv{-}\int_0^\rho \,{\rm d}\varrho\; \varLambda_{,\varrho} \alpha ,\quad \varXi\equiv{-}\left( \varLambda_{,\ell}\alpha + \mathcal{I}_{,\ell}\right),\\ \mathcal{G} & =\int_0^\rho \,{\rm d}\varrho\;G \psi_{,\varrho}, \quad \mathcal{Y}\equiv G \psi_{,\omega} - \mathcal{G}_{,\omega},\quad \mathcal{X}\equiv G \psi_{,\ell} - \mathcal{G}_{,\ell},\quad \mathcal{W}=\int_0^\rho \,{\rm d}\varrho \frac{\varrho}{h}\lambda \left( \xi_{,\ell}+\varXi+\mathcal{X}\right) . \end{aligned} \right\} \end{align}

The Poisson equation for $\xi$ is now given by

(5.17)\begin{equation} \rho^2\Delta \xi +\frac{1}{h}\partial_\omega\left( h \left( \mathcal{Y}+ \mathcal{W}\right) \right) +\frac{1}{h}\partial_\ell\left( \frac{\rho^2}{h}(\varXi+\mathcal{X})\right)=0. \end{equation}

As in the force-free case, $\xi$ plays the analogous role of $\varPhi$ in the vacuum limit, and it can be shown (see Appendix E) that the Poisson equation for $\xi$ is free of any resonances.

5.4 The NAE of MHS equations

The expansions of the quantities $\mathcal {I}, \varXi$ that appear in (5.16) have already been discussed in (4.24). From the NAE, (5.9) and (2.22) we find that the quantities $\mathcal {G},\mathcal {Y},\mathcal {X},\mathcal {I},\varXi$, as defined in (5.16d), have the following forms:

(5.18)\begin{align} \left. \begin{aligned} \mathcal{G} & = \rho^2 G^{(0)}(\ell) \psi^{(2)} +O(\rho^3), \quad \mathcal{Y} =O(\rho^3),\quad \mathcal{X} ={-}\rho^2 {G^{(0)}}'(\ell) \psi^{(2)} +O(\rho^3),\\ \mathcal{I} & ={-}\lambda_0 \rho^2 \alpha^{(0)} \psi^{(2)} +O(\rho^3), \quad \varXi=\rho^2 \lambda_0 \psi^{(2)} \alpha_{,\ell}^{(0)} +O(\rho^3),\quad \mathcal{W}=\frac{\lambda_0}{2}\xi^{(0)}_{,\ell}\rho^2 +O(\rho^3). \end{aligned} \right\} \end{align}

Therefore, the deviation of the MHS system (5.16) from the vacuum system (3.2) is zero to $O(\rho ^{-2}),O(\rho ^{-1})$, just like the force-free equations (4.21). Moreover, (5.18) also implies that the expressions for the lowest-order components of $\boldsymbol {B}$ are the same as in the force-free case (4.28). As a result, MHS and force-free both have the same form of MDO given by (4.28). Consequently, the basic framework is analogous to the vacuum system and the effects of currents and pressure show up mostly through the forcing terms $F_\phi,F_\psi,F_\alpha$.

The only additional step in MHS is Step 0, which involves solving the MDE for $G$ order by order. Step 0 follows the same procedure that we have discussed for the solution of the generic MDE in § 3.6.2. To $O(\rho ^0)$, the MDE for $G$ is given by

(5.19)\begin{equation} G^{(0)}_{,\ell}\varPhi_0'(\ell)=p^{(2)}, \end{equation}

which can be easily solved to obtain the following:

(5.20)\begin{equation} G^{(0)}=p^{(2)}\int \dfrac{{\rm d}\ell}{\varPhi_0'(\ell)}, \quad K^{(0)}= G^{(0)}-\lambda_0\alpha^{(0)}. \end{equation}

The details of the general first-order solutions are given in Appendix G. Due to the similarities between the force-free and the vacuum problem, we do not explicitly show Steps I–IV.

6.1 Examples of vacuum fields

We now discuss the next to leading order, $n=1$, solutions to illustrate some of the key features of the NAE for the vacuum case. Our goal is also to clarify some of the steps in the Mercier–Weitzner formalism discussed earlier. We begin with the straightforward case where the lowest order flux-surface $\psi ^{(2)}$ has a circular cross-section. The general case (with rotating ellipse, triangularity, etc.) is described in Appendix F.

We consider the case where $a(\ell )=1/2, b(\ell )=0$ in (3.16a,b), which corresponds to a circular cross-section to the lowest order. For this case, the NAE equations to $O(1)$ (3.11) imply that

(6.1)\begin{equation} \psi^{(2)}=\frac{1}{2}\varPhi_0', \quad \alpha^{(0)}= \omega =\theta-\int \tau \,{\rm d}\ell , \quad \varPhi^{(2)}={-}\frac{\varPhi_0''}{4}. \end{equation}

To $O(\rho )$ the relevant equations for the variables $(\varPhi ^{(3)},\alpha ^{(1)},\psi ^{(3)})$ takes the form

(6.2a)\begin{gather} -(\partial^2_\omega+3^2)\varPhi^{(3)} = \cos\theta \left( \kappa' \varPhi_0'+\frac{5}{2}\kappa \varPhi_0''\right) + \sin\theta \kappa \tau \varPhi_0' , \end{gather}

(6.2b)\begin{gather}\left( \varPhi_0'\partial_\ell -\frac{3}{2}\varPhi_0''\right) \psi^{(3)} + \varPhi_0'\left( 3 \varPhi^{(3)}+\varPhi_0''\kappa \cos \theta \right) =0, \end{gather}

(6.2c)\begin{gather}- \varPhi_0'\alpha^{(1)}_{,\omega}=\left( 3\psi^{(3)}-\varPhi_0' \kappa \cos \theta\right) . \end{gather}

Equation (6.2a) is the Poisson equation for $\varPhi ^{(3)}$ whose solution can be readily seen to be

(6.3)\begin{align} \varPhi^{(3)}={-}\frac{1}{8}\cos u \left( \kappa'\varPhi_0'+\frac{5}{2}\kappa \varPhi_0'' \right) \,{-}\,\frac{\varPhi_0'}{8}\kappa \tau \sin u \,{-}\,\frac{\left( \varPhi_0'\right)^{3/2}}{9} \frac{\partial}{\partial \ell}\left( Q^{(3)} \cos{3u}+P^{(3)} \sin{3u}\right). \end{align}

Here, $u=\theta$, and we have used (3.70), the Soloviev form of the free function, in representing $\varPhi ^{(3)}$.

Next, we solve (6.2b), which is the MDE for $\psi ^{(3)}$. The first step is to eliminate the linear term in $\psi ^{(3)}$ through a change of variables,

(6.4)\begin{equation} \psi^{(3)}= \left( \varPhi_0' \right)^{3/2}Y^{(3)}, \quad Y^{(3)}_{,\ell}+\frac{1}{\left( \varPhi_0' \right)^{3/2}}\left( 3\varPhi^{(3)}+\varPhi_0'' \kappa \cos\theta\right)=0. \end{equation}

Substituting the expression for $\varPhi ^{(3)}$ (6.3) into (6.4) we find

(6.5)\begin{equation} \frac{\partial}{\partial \ell}\left( Y^{(3)} -\frac{3}{8\sqrt{\varPhi_0'}}\kappa \cos{u} -\frac{1}{3}\left( Q^{(3)} \cos{3u}+P^{(3)} \sin{3u}\right)\right)+ F^{(3)}_{\psi c1}\cos u =0, \end{equation}

with

(6.6)\begin{equation} F^{(3)}_{\psi c1}={-}\frac{1}{{\left( \varPhi_0' \right)^{3/2}}} \frac{\kappa}{8}\varPhi_0''. \end{equation}

The solution of (6.2b) is, therefore, of the form

(6.7)\begin{equation} \psi^{(3)} =\frac{3}{8}\varPhi_0'\kappa \cos u+\left(\varPhi_0'\right)^{3/2}\left( \frac{1}{3}\left( Q^{(3)} \cos{3u}+P^{(3)} \sin{3u}\right)+\left( \sigma^{(3)}_{c1}\cos u + \sigma^{(3)}_{s1}\sin u \right)\right), \end{equation}

where, $\sigma ^{(3)}_{c1}, \sigma ^{(3)}_{s1}$ satisfy

(6.8)\begin{equation} {\sigma^{(3)}_{c1}}'-\tau \sigma^{(3)}_{s1}+F^{(3)}_{\psi c1}=0, \quad {\sigma^{(3)}_{s1}}'+\tau \sigma^{(3)}_{c1}=0. \end{equation}

Note that due to our choice of the free function in (6.3), the third harmonics of $\psi ^{(3)}$ are obtained explicitly in (6.7). We can now combine (6.8) into a single complex equation

(6.9)\begin{equation} \left. \begin{aligned} & { \mathcal{Z}^{(3)}_\psi}'-{\rm i}\, \varOmega_0 \tau { \mathcal{Z}^{(3)}_\psi} + \mathcal{F}^{(3)}_\psi=0,\\ & \mathcal{Z}^{(3)}=\sigma^{(3)}_{c1} + {\rm i}\, \sigma^{(3)}_{s1}, \quad \mathcal{F}_\psi^{(3)}=F^{(3)}_{\psi c1} + {\rm i}\, F^{(3)}_{\psi s1},\quad \varOmega_0={-}\tau. \end{aligned} \right\} \end{equation}

Note that we can add a homogeneous solution to $\psi ^{(3)}$ of the form

(6.10)\begin{equation} \psi^{(3)}_H=\left( \varPhi_0' \right)^{3/2}Y^{(3)}_H. \end{equation}

We now solve (6.2c) to obtain $\alpha ^{(1)}$ up to a function $\mathfrak {a}^{(1)}(\ell )$. Substituting the solution of $\psi ^{(3)}$ and integrating (6.2c) with respect to $\omega$, we obtain

(6.11)\begin{align} \alpha^{(1)} & ={-}\frac{\kappa}{8} \sin u -\frac{\left( \varPhi_0'\right)^{1/2}}{3}\left( Q^{(3)} \sin{3u}-P^{(3)} \cos{3u}\right)\nonumber\\ & \quad +\mathfrak{a}^{(1)}(\ell)-3\left( \varPhi_0'\right)^{1/2}\left( \sigma^{(3)}_{c1}\sin u - \sigma^{(3)}_{s1}\cos u \right). \end{align}

To determine $\mathfrak {a}^{(1)}(\ell )$ we need the MDE for $\alpha ^{(1)}$,

(6.12)\begin{equation} {-}\frac{1}{2}\varPhi_0'' \alpha^{(1)}+ \varPhi_0'\alpha^{(1)}_{,\ell} +\varPhi^{(3)}_{,\omega}=0. \end{equation}

Since $\alpha ^{(1)}$ and $\varPhi ^{(3)}$ only have odd harmonics, the poloidal average of (6.12) is zero. Hence, both $Y^{(3)}_H$ and $\mathfrak {a}^{(1)}(\ell )$ are zero.

When the on-axis magnetic field, $\varPhi _0'=B_0$, is a constant, the NAE results simplify even further. For $\varPhi _0''=0$, we have $\sigma ^{(3)}_{c1}=\sigma ^{(3)}_{s1}=0$, and we obtain the following explicit forms for the lower-order quantities:

(6.13a)\begin{gather} \varPhi^{(2)}=0,\quad \psi^{(2)}=\frac{1}{2}B_0, \quad \alpha^{(0)}= \omega =\theta-\int \tau \,{\rm d}\ell , \quad u=\theta , \end{gather}

(6.13b)\begin{gather}\varPhi^{(3)}={-}B_0\frac{\partial}{\partial \ell}\left( \frac{\kappa}{8}\cos{u}+\frac{\sqrt{B_0}}{9} \left( Q^{(3)} \cos{3u}+P^{(3)} \sin{3u}\right) \right), \end{gather}

(6.13c)\begin{gather}\psi^{(3)}=B_0\left(\frac{3}{8}\kappa \cos u +\frac{\sqrt{B_0}}{3} \left( Q^{(3)} \cos{3u}+P^{(3)} \sin{3u}\right)\right), \end{gather}

(6.13d)\begin{gather}\alpha^{(1)}={-}\frac{\kappa}{8}\sin{u} -\frac{\sqrt{B_0}}{3}\left( Q^{(3)} \sin{3u}-P^{(3)} \cos{3u}\right). \end{gather}

6.2 Examples of force-free and MHS fields

We shall now compare the NAE of force-free and MHS equilibrium with the ‘one-size fits all’ Soloviev profiles (Cerfon & Freidberg Reference Cerfon and Freidberg2010). The Soloviev profiles make the Grad–Shafranov equation (5.5) linear and, therefore, analytically tractable. In the following, we shall restrict ourselves to first-order NAE, i.e. the accuracy of $\psi$ to $O(\rho ^4)$.

To the required order of accuracy, the one-size model is given by (Cerfon & Freidberg Reference Cerfon and Freidberg2010)

(6.14)\begin{equation} \varPsi(x,y)=\frac{x^4}{8}+A\left( \frac{x^2}{2}\log{x} -\frac{x^4}{8}\right) +c_1 +c_2\left( x^2\right) +c_3 \left( y^2-x^2\log{x}\right), \end{equation}

where, $A$ is a constant and $x,y$ represent the axisymmetric cylindrical $R,z$ coordinates, normalized by the major radius $R_0$. The normalized pressure and current profiles satisfy

(6.15)\begin{equation} {-}\frac{1}{B_0^2}\dfrac{{\rm d}p}{{\rm d}\varPsi}=1-A, \quad -\frac{1}{B_0^2}F(\varPsi)F'(\varPsi)=A. \end{equation}

Equivalently, the Soloviev profiles are defined by

(6.16)\begin{equation} \frac{p(\varPsi)}{B_0^2}=\frac{p_0}{B_0^2}-(1-A)\varPsi, \quad \frac{F(\varPsi)}{B_0}=\sqrt{\left( \frac{F_0}{B_0}\right)^2-2 A \varPsi}, \end{equation}

where, $p_0,F_0$ are constants. The $A=1$ limit corresponds to force-free equilibrium. The poloidal flux function $\varPsi$, given by (6.14), satisfies the following Grad–Shafranov equation:

(6.17)\begin{equation} x \partial_x (x^{{-}1}\varPsi_{,x})+\varPsi_{,yy}=(1-A)x^2+A. \end{equation}

We choose a circular cross-section to the lowest order for simplicity while allowing for arbitrary shaping at higher orders. The details of the NAE are provided in Appendix H along with other analytical examples. To $O(\rho ^4)$, the NAE matches the exact one-size results for both force-free and MHS cases. Figure 1 compares the exact one-size solution with the NAE in the case of force-free ($\beta =0\,\%$) and MHS at $\beta =5\,\%$. We see that the deviations of NAE from the exact are more on the inboard side, where the surfaces have a more significant shaping.

Figure 1. Comparison of exact axisymmetric force-free (a) and MHS (b) equilibrium at 5 % plasma $\beta$ obtained for Soloviev profiles to the NAE. The effect of $\beta$ is barely visible. The NAE matches exactly to the one-size model up to $O(\rho ^4)$. Note that the NAE deviates significantly from the exact solution only when the aspect ratio is sizable ($\approx 1/2$).