Stellarator equilibrium axis-expansion to all orders in distance from the axis for arbitrary plasma beta

Abstract

A systematic theory of the asymptotic expansion of the magnetohydrostatics (MHS) equilibrium in the distance from the magnetic axis is developed to include arbitrary smooth currents near the magnetic axis. Compared with the vacuum and the force-free system, an additional magnetic differential equation must be solved to obtain the pressure-driven currents. It is shown that there exist variables in which the rest of the MHS system closely mimics the vacuum system. Thus, a unified treatment of MHS fields is possible. The mathematical structure of the near-axis expansions to arbitrary order is examined carefully to show that the double-periodicity of physical quantities in a toroidal domain can be satisfied order by order. The essential role played by the leading-order Birkhoff–Gustavson normal form in solving the magnetic differential equations is highlighted. Several explicit examples of vacuum, force-free and MHS equilibrium in different geometries are presented.

1 Introduction

Understanding the magnetohydrostatics (MHS) equilibrium in magnetically confined devices is paramount. The theory of axisymmetric MHS equilibrium, applicable to tokamaks, is well-developed at this point. The axisymmetric problem can be reduced to solving a nonlinear elliptic partial differential equation (PDE) called the Grad–Shafranov equation (Grad 1971; Freidberg 1982). In contrast, the three-dimensional (3-D) theory cannot be reduced to a single elliptic PDE (Bauer, Betancourt & Garabedian 2012b; Helander 2014; Weitzner 2014). Analytical understanding is obstructed by the intrinsic nonlinearity of the MHS system, coupled with the 3-D nature of the 3-D equilibrium, which is generic for stellarators. These obstacles hinder an intuitive understanding of important physical effects, such as the deformations of the vacuum plasma surfaces due to pressure-driven currents or the dependence of magnetic shear on the geometry. Large-scale numerical computation has taught us much about stellarator equilibrium, but simulations can be expensive. Furthermore, the computational error can be significant, particularly near the magnetic axis (Panici et al. 2022). Therefore, obtaining an analytical understanding of 3-D stellarator MHS equilibrium is beneficial for computational studies and their applications to experiments.

In this work, we shall focus on an asymptotic expansion called the near-axis expansion (NAE) developed by several authors (Mercier 1964; Solov'ev & Shafranov 1970; Bernardin & Tataronis 1985; Garren & Boozer 1991; Weitzner 2016). The basic idea is to expand the MHS equilibrium equations in a power series in the distance from the magnetic axis. The small parameter is $\kappa a$, where $\kappa$ is the curvature of the magnetic axis, and $a$ is a typical length scale of the order of the minor radius in the direction normal to the magnetic axis. Recently, there has been a surge of interest in applying the NAE to study quasisymmetry (Landreman & Sengupta 2018, 2019; Jorge, Sengupta & Landreman 2020a; Rodríguez, Sengupta & Bhattacharjee 2023), quasiisodynamic stellarators (Plunk, Landreman & Helander 2019; Mata, Plunk & Jorge 2022; Mata & Plunk 2023; Rodríguez & Plunk 2023), energetic particles (Figueiredo et al. 2023), MHS stability (Landreman & Jorge 2020; Kim, Jorge & Dorland 2021; Rodríguez 2023) and turbulence (Jorge & Landreman 2020, 2021; Rodríguez & Mackenbach 2023).

The development of the NAE theory has proceeded along two approaches: direct and indirect. In the direct approach (Mercier 1964; Solov'ev & Shafranov 1970; Bernardin & Tataronis 1985; Jorge et al. 2020a; Jorge, Sengupta & Landreman 2020b; Duignan & Meiss 2021), one uses a coordinate system pioneered by Mercier, and expands in the distance from the magnetic axis. In the inverse approach, one expands in the toroidal flux-surface label (Garren & Boozer 1991; Weitzner 2016; Landreman & Sengupta 2018, 2019; Landreman, Sengupta & Plunk 2019; Rodríguez & Bhattacharjee 2021; Rodríguez, Sengupta & Bhattacharjee 2021, 2022).

Both the direct and the indirect approaches have distinct advantages and disadvantages. The inverse coordinate approach is best suited to address concepts such as quasisymmetry and omnigeneity, which impose constraints on the magnitude of the magnetic field. However, the direct approach offers better physical insight. This is mainly because the flux-surface shaping factors, such as ellipticity and triangularity, can be chosen independently in the direct approach but not in the indirect approach. Moreover, in the indirect approach, the lowest order involves the solution of a nonlinear Riccati equation. In the direct coordinate approach, we can avoid solving a nonlinear ordinary differential equation (ODE) by using Mercier's formula (Solov'ev & Shafranov 1970; Mercier 1974), which also follows from Floquet theory (Duignan & Meiss 2021). The expressions for the rotational transform (Mercier 1964; Helander 2014) and its derivatives are also physically more intuitive in the direct approach.

An outstanding issue that arises while constructing NAE solutions order by order is the maintenance of the smoothness of the solutions of the MHS equilibrium. If we assume that the magnetic field satisfies MHS with smooth, nested toroidal pressure surfaces, it can be shown rigorously (Burby, Duignan & Meiss 2021) that the NAE can be assured to be smooth. This rigorous result has been extended recently to smooth vacuum and force-free fields assuming smooth nested flux surfaces (Perrella, Duignan & Pfefferlé 2023). These deep mathematical results hold even in the case that iota ($\iota$) is irrational. Although we are guaranteed that smooth NAE solutions will exist, several difficulties arise when constructing the NAE solution explicitly order by order.

The first problem we encounter is the classic problem of the occurrence of singular currents in stellarators on rational surfaces because of a lack of axisymmetry (Mercier 1964; Grad 1967; Helander 2014; Loizu et al. 2015; Weitzner 2016). These singularities are fundamentally tied to the magnetic differential equations (MDE) (Newcomb 1959), whose solution can be singular everywhere, $\iota$, the rotational transform, is rational. Since rationals are dense, avoiding such singularities with non-constant $\iota$ is challenging (Grad 1967; Boozer 1981b; Weitzner 2014). In the NAE expansion, thankfully, the rational singularities are governed only by $\iota _0$, the on-axis rotational transform (Mercier 1964; Solov'ev & Shafranov 1970; Duignan & Meiss 2021). Choosing it to be sufficiently irrational helps with the practical implementation of lower-order NAE (Landreman & Sengupta 2019; Jorge et al. 2020b). However, since any irrational number can be approximated closely by rationals, the problem of singular divisors appearing in higher orders can only be partially resolved mathematically in this approach. An alternative is to assume that $\iota _0$ is close to a given rational number and construct special non-singular solutions (Weitzner 2016; Jaquiery & Sengupta 2019) or use the near-resonant normal form theory (Duignan & Meiss 2021).

We remark that there could be other potential sources of non-smoothness as well. Even in the vacuum limit, the smoothness of the solution within the NAE framework is not obvious (Jorge et al. 2020b). When currents are present, the situation is complicated by another factor: the possibility of occurrence of logarithmic terms of the form $\rho ^n (\log \rho )^m, n>m$ that, if present, make the equilibrium weakly singular near the axis ($\rho =0$) (Weitzner 2016). Perhaps the easiest way to understand the physical origin of these logarithmic terms is to think of the plasma in the large aspect ratio picture as a thin current-carrying loop (Freidberg 1982; Jackson 1999). The logarithmic terms are then unavoidable in the poloidal flux function (Freidberg 1982). However, this simplified description of plasma needs to be modified near the axis to avoid the singularities on the axis. As shown in Greene, Johnson & Weimer (1971), removing such weakly singular terms is tricky, even at lower orders. These issues are expected to worsen in stellarators where the currents interact with the 3-D geometry. In this work, we develop a systematic way to avoid such singular terms within the NAE formalism, showing that the assumption of smoothness of the NAE can be satisfied order by order.

The current work generalizes the direct NAE approach of Jorge et al. (2020b) by including force-free and plasma-beta-driven currents, and by discussing both the structure of the flux surface and the field line label to all orders. The outline of the paper is as follows. In § 2, we introduce the basic formalisms developed by Mercier (Mercier 1964, 1974) and Weitzner (Weitzner 2014, 2016) and show how they can be combined into one, which we call the Mercier–Weitzner formalism. We assume the rotational transform to be sufficiently irrational. With this assumption, we can solve the MDEs order by order without encountering resonances on rational surfaces. We discuss the NAE of the vacuum fields in § 3, force-free fields in § 4 and finally, MHS fields with finite plasma beta in § 5.

Our strategy is to construct a new variable that mimics the vacuum scalar potential but fully incorporates the effects of currents. The NAE using this variable shows a structure very close to the vacuum problem. Under the assumption of smooth fields and currents near the axis, we can show that power-series solutions to ideal MHS can be constructed to all orders. We do not attempt to prove convergence of the series. We show that the NAE calculations are simplified to a great extent through the use of the leading-order Birkhoff–Gustavson Hamiltonian normal form variables (Bernardin & Tataronis 1985; Duignan & Meiss 2021). Normal forms in the present work show up as a basis that diagonalizes the magnetic differential operator (MDO), thereby simplifying the calculations of the solutions of MDEs ubiquitous in MHS.

In § 7, we present a numerical example that compares the higher-order finite-beta near-axis asymptotic expansions with a finite aspect-ratio equilibrium generated using the VMEC code (Hirshman & Whitson 1983). We present our conclusions and open questions in § 8. Finally, a variety of special cases are discussed in the appendices.

2 The Mercier–Weitzner near-axis formalism

Our goal in this work is to extend the direct approach to the near-axis formalism, as developed by Mercier (Mercier 1964) and Solov’ev & Shafranov (Solov'ev & Shafranov 1970), to include finite plasma beta. Mercier's original approach (Mercier 1974) indeed considers currents and finite beta. However, the methodology to calculate beyond the two lowest orders is not transparent. Solov'ev & Shafranov (1970) extend Mercier's work significantly but primarily for vacuum fields. Effects of plasma beta, such as on the flux surface shapes, are treated through expansion subsidiary to the NAE. Only recently, a complete description of NAE in direct coordinates that included both vacuum and force-free equilibria was obtained (Duignan & Meiss 2021).

Here we provide a unified framework to include vacuum, force-free as well as finite beta ideal MHS. We shall use Weitzner's approach (Weitzner 2014, 2016) to obtain finite-beta MHS equilibrium equations. In this approach, currents are included through the current potential obtained through the solution of an MDE. In the following, we describe the NAE in direct coordinates using the current potential approach due to Weitzner.

2.1 The Weitzner formulation of ideal MHS with scalar pressure

We begin with the discussion of the formalism developed in Weitzner (2014, 2016) 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. 2012; Helander 2014) 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. (2012)). 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 1980, 1981a,b) 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 1971; Weitzner 2016) 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 (2014), 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 1964), which are also known as the direct coordinates (Jorge et al. 2020b).

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 1970; Duignan & Meiss 2021), 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. 2019; Mata & Plunk 2023). 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 (1964) 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 1987; Garren & Boozer 1991) 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. (2020b) and Duignan & Meiss (2021), 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 (1964) and Solov'ev & Shafranov (1970), 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 (2016). 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 The NAE of vacuum fields in Mercier–Weitzner formalism

The NAE of vacuum fields in direct coordinates has been treated extensively in Jorge et al. (2020b). In the Mercier–Weitzner formalism, $K=0$ corresponds to the vacuum limit, and the analysis is similar to Jorge et al. (2020b), which also uses direct coordinates. However, a key point of departure of our present work from all other previous developments of NAE is the treatment of the Clebsch variable, which represents the field line label, $\alpha$. The primary motivation for focusing on $\alpha$ is that the rotational transform and its derivatives (which includes the magnetic shear) are related directly to the secular part of $\alpha$ (through (2.11b)). Thus, a systematic construction of $\alpha$ is crucial.

Our goal in this section is to study the structure of the Mercier–Weitzner equations in the vacuum limit. First, we show how the condition of nested flux surfaces leads to MDEs. We then explicitly employ the NAE and elucidate the properties of the solutions obtained, paying particular attention to the secular and periodic structure of the field line label, $\alpha$.

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 1958; Goedbloed 1983; Courant & Hilbert 2008; Weitzner 2014) 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 (2014). 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 1970) 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 (2021). 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 2014).

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. (2020b). 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. (2020b), 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 1959) 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 (1970), 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 (1970) 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 (2021) and Jorge et al. (2020b). More concretely, ‘analyticity’ implies (Duignan & Meiss 2021; Solov'ev & Shafranov 1970)

(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 (2021), 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.

3.6.6 Summary of Steps I–IV

We now briefly summarize the main results obtained in § 3.6. We have individually analysed the structure of $(\varPhi,\psi,\alpha )$, proving that the necessary periodicity conditions (3.28a,b) for obtaining physically meaningful solutions can be satisfied order by order. In § 3.6.1 we have discussed the $O(n+2)$ Laplace equation and the structure of its solution $\varPhi ^{(n+2)}$. We have also given a prescription, (3.70), with which the free functions (the homogeneous solution to the Poisson equation) can be expressed in terms of the flux-surface shaping given by (3.71a,b).

In § 3.6.2, we study the general MDE and obtain some of its fundamental properties. These properties allow one to construct homogeneous solutions of the MDEs and annihilators for the MDEs. Later, we use the homogeneous solution for $\psi ^{(n+2)}$ and the annihilator to enforce the periodicity constraint (3.42). The last important property is that the poloidal harmonics of the solution of the MDE truncate if the forcing term has finite poloidal harmonics. Using these properties we found that at $O(n)$, $\varPhi ^{(n+2)}$ will have at most $(n+2)$ harmonics. The structure of $\alpha ^{(n)}$ is more complicated, and solving the MDE for $\alpha$ (3.27) is far more involved than MDE for $\psi$ given by (3.61). Fortunately, we do not have to solve the MDE for $\alpha$.

Finally, in § 3.6.5, we demonstrate how one can obtain $\alpha ^{(n)}$ once $\varPhi ^{(n+2)},\psi ^{(n+2)}$ has been constructed. This brings us to the end of the analysis of the structure of vacuum equations in the Mercier–Weitzner formalism. In § 6.1 and Appendix H, we discuss analytically tractable and physically interesting limits of 3-D vacuum fields that possess nested surfaces. We can obtain explicit solutions order by order using the NAE formalism we have developed thus far.

4 The NAE of force-free fields

For force-free systems, the current $\boldsymbol {J}$ and magnetic field $\boldsymbol {B}$ are colinear, i.e.

(4.1)\begin{equation} \boldsymbol{J}=\lambda \boldsymbol{B}, \quad \boldsymbol{J}=\boldsymbol{\nabla}\times \boldsymbol{B}. \end{equation}

Since the divergence of $\boldsymbol {J}$ is zero, we have

(4.2)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \lambda =0. \end{equation}

We further assume that the force-free fields possess nested flux surfaces. Therefore, (4.2) implies that $\lambda$ must be a function of only $\psi$ and $\alpha$. Physically, the $\alpha$ dependence of $\lambda$ is special. It is only possible if the rotational transform is a constant rational number. Any amount of magnetic shear will lead to

(4.3)\begin{equation} \lambda= \lambda(\psi), \end{equation}

which we assume from now on.

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 (2016), 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 1964)

(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 (1974), 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 The MHS equilibrium with finite $\beta \sim p'/B^2$

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 1982)

(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 1982; Helander 2014):

(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 2012a; Weitzner 2014). 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 Construction of analytical solutions in the Mercier–Weitzner formalism

We have completed the development of the Mercier–Weitzner formalism of the direct NAEs. In this section, we shall provide explicit analytical solutions to low orders for vacuum, force-free and MHS fields that can be obtained within the Mercier–Weitzner formalism.

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 2010). 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 2010)

(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$).

7 Numerical implementation

In this section, we present some numerical examples of the second-order near-axis equilibria constructed numerically using the formalism developed in the previous sections. We verify aspects of this implementation by generating a plasma boundary at a finite aspect ratio, using it to construct a global equilibrium using the VMEC code (Hirshman & Whitson 1983), and comparing some of its properties with the near-axis equilibrium. We detail this procedure in what follows.

The examples presented here are constructed with an axis with constant torsion and a nearly circular cross-section, which we chose for simplicity. The details on the axis are provided in Appendix H.4. The choice of nearly circular cross-sections, discussed in detail in Appendix H.4, leads to the toroidal flux function,

(7.1)\begin{equation} \psi = \rho^2\left(\frac{1}{2}+\rho \psi_3\right), \end{equation}

where

(7.2)\begin{equation} \psi_3 = \frac{3 \kappa}{8}\cos u + \sigma^{(3)}_{c1}\cos u + \sigma^{(3)}_{s1}\sin u, \end{equation}

with $\sigma ^{(3)}_{c1}$ and $\sigma ^{(3)}_{s1}$ given by (H42).

In addition, we choose $\varPhi _0=l$ so that the magnetic field is constant on the axis, $B_0=1$, a feature easily checked in the finite aspect ratio equilibrium constructed. Given the emphasis in this paper on the near-axis construction for finite current and pressure, we choose finite values for these: $\lambda _0=1.2$ and $p_2=-0.1$.

This information is sufficient to describe the near-axis equilibrium and to construct a flux surface at a finite aspect ratio. Once such a surface has been constructed, we find a global equilibrium solution with VMEC. The position vector $\boldsymbol {r}$ for a surface at constant $\psi$ is given by

(7.3)\begin{equation} \boldsymbol{r} = \boldsymbol{r}_0 + \rho(\cos \theta \boldsymbol{n} + \sin \theta \boldsymbol{b}), \end{equation}

where $\rho$ is found by inverting (7.1), namely $\rho = \sqrt {2 \psi }- 2 \psi _3 \psi$, at given values of $\psi$. Surfaces of constant $\psi$ at $\phi =l=0$ are shown in figure 2 from the magnetic axis to $\psi =0.001$ (figure 2a, high-aspect ratio $\epsilon =1/50$) and from the magnetic axis to $\psi =0.02$ (figure 2b, moderate aspect ratio $\epsilon =1/11$). Note the noticeable difference in the Shafranov shift of these configurations. This difference is consistent with their effective plasma $\beta$ on axis: for the high-aspect-ratio $\beta =0.12\,\%$, and the moderate aspect ratio case, $\beta =2.69\,\%$.

Figure 2. Cross-sections of constant $\psi$ for near-axis construction. Surfaces of constant $\psi$ at $\phi =l=0$ following (7.3) in the range from $\psi =0$, the magnetic axis, to (a) $\psi =0.001$ and (b) $\psi =0.02$.

To proceed further with the comparison, we need to compute the Fourier coefficients of the finite aspect ratio flux surface. That is, the Fourier coefficients of $R(\theta, \phi _c)$ and $Z(\theta, \phi _c)$ in cylindrical coordinates $(R, Z, \phi _c)$ for the plasma surface. Note that this cylindrical coordinate system is not the Mercier coordinate system, the natural coordinate system for near-axis construction. In particular, the toroidal angle in cylindrical coordinates is not the same as the toroidal angle associated with the length $\ell$ along the magnetic axis. Thus, if we need to use the cylindrical angle off-axis $\phi _c(\psi, \theta, \phi )$ as our coordinate for the surface, we must invert its definition,

(7.4)\begin{equation} \tan \phi_c = \frac{\boldsymbol{r}_y(\psi,\theta,\phi)}{\boldsymbol{r}_x(\psi,\theta,\phi)}, \end{equation}

to compute the toroidal angle on-axis $\phi$. This can be achieved using a numerical root solver based on Newton's method. Once this inversion is achieved, we calculate the Fourier components $(RBC,ZBS)$ of $R(\theta,\phi _c)=\sum _{m,n}RBC_{m,n}\cos (m \theta - n\phi _c)$ and $Z(\theta,\phi _c)=\sum _{m,n}ZBS{m,n}\sin (m \theta -n \phi _c)$ for $0\le m \le 10$ and $0-10\le N \le 10$. The surface defined this way is now in a form to be used as an input to the equilibrium solve in VMEC. Comparing figure 2 to figure 3, we find that although the flux surfaces are nearly circular in the Frenet–Serret frame of the magnetic axis, they have significant shaping in the laboratory frame, most notably in elongation. This difference arises because of the non-circular geometry of the axis. In figure 4, we have shown the three-dimensional shapes of the outermost flux surfaces.

Figure 3. Cross-sections of the global equilibrium for $\psi =0.02$. Example of cross-sections for the global VMEC equilibrium solution for $\psi =0.02$. Cross-sections at $\phi _c=0,{\rm \pi} /4$ and ${\rm \pi} /2$ are shown at values of normalized toroidal flux from $0.1- 1.0$. The cross-sections show the non-trivial shaping of the example field used in this section.

Figure 4. Three-dimensional representation of the constructed equilibria. Plots of the boundary of the equilibria constructed from the near-axis solution for (a) $\psi =0.001$ and (b) $\psi =0.02$. Some cross-sections for the latter are shown in figure 3. The colour map represents the magnitude of the magnetic field on the surface.

The surfaces from the near-axis construction and the finite aspect ratio equilibrium are compared for several aspect ratios in figure 5. There are clear differences in the shapes of both equilibria, but these tend to decrease in the limit of increasing aspect ratio. For a more accurate matching, as discussed in Landreman & Sengupta (2019), the error incurred in setting the minor radius $\sqrt {\psi }$ to a finite value, say $a$ when constructing an equilibrium makes a one-to-one comparison between the solutions complicated. A double expansion in both small parameters, $\sqrt {\psi },a$ must be carried out to account for finite aspect ratio. Or construct the finite aspect ratio equilibrium in an alternative way (Panici et al. 2022, 2023). This is particularly important for comparing higher-order NAE features, such as the magnitude of the Shafranov shift. A more accurate and complete numerical description is left for future work.

Figure 5. Comparison of cross-sections from NAE and VMEC calculations. The plots show a comparison between the cross-sections at $\phi _c=0$ between the global equilibria computed with VMEC for $\psi =0.02,0.01,0.005,0.003$ (a–d) and the near axis solution (broken lines). The comparison of the NAE solution with the finite aspect ratio VMEC gets better with increasing aspect ratio. The flux surfaces shown correspond to $\psi =0.02,0.01,0.005,0.003,0.001$ and the magnetic axis.

The rotational transform provides another important feature of the equilibria to study. Using Mercier's formula for the rotational transform on-axis $\iota _0$ (Mercier 1964), the near-axis construction expects a value

(7.5)\begin{equation} \iota_0=\left(\frac{\lambda}{2}-\tau\right)\frac{L}{2{\rm \pi}}-N=\left(\frac{1.2}{2}+1.643\right)\frac{6{\rm \pi}}{2{\rm \pi}}-6=0.729. \end{equation}

Figure 6 shows the difference between the NAE rotational transform and the one obtained from the finite aspect ratio equilibria. The difference tends to zero quadratically in the aspect ratio, suggesting that the construction of the equilibrium is correct to first order. The differences in rotational transform may be compared in magnitude with the total global shear of the equilibria.

Figure 6. Rotational transform profile of global equilibria and difference to the NAE. (a) Rotational transform $\iota$ as a function of the normalized toroidal flux $s=\psi /\psi _b$ computed by VMEC for the equilibria with $\psi =0.001$ and $0.02$. The lower aspect ratio case shows higher magnetic shear. (b) Difference between the on-axis rotational transform between the finite aspect ratio VMEC equilibria and the near axis value, as a function of the aspect ratio. The dashed line shows a scaling $\epsilon ^2$.

Finally, let us recall that the near-axis construction was made for a constant magnetic field on the axis. Thus, studying the spectrum of Fourier modes of the magnetic field strength $B$ should also be informative. Figure 7 shows the spectrum computed using the BOOZ_XFORM code (Sanchez et al. 2000). We find that, indeed, on the axis, the magnetic field is approximately constant. A mirror mode growing with decreasing aspect ratio is apparent, in line with previous findings (Landreman & Sengupta 2019; Jorge et al. 2020a), and a more sophisticated way of generating finite-aspect-ratio equilibria would be necessary (Landreman & Sengupta 2019) to suppress them.

Figure 7. Boozer spectrum of $|\boldsymbol {B}|$ for global equilibria. The plots show the Boozer spectra of $|\boldsymbol {B}|$ as a function of the normalized toroidal flux $s=\psi /\psi _b$, for the configurations constructed at $\psi =0.001$ (a) and $\psi =0.02$ (b). The spectra were computed using the BOOZXFORM code.

8 Discussion

In this work, we have presented an asymptotic theory of MHS equilibrium in a stellarator using the NAE scheme. Our approach generalizes previous work in this area by allowing MHS fields of arbitrary plasma beta. We have shown that the analyticity property of the vacuum fields, i.e. the existence of regular power-series expansions in the vicinity of the magnetic axis, can be retained if one allows for current profiles with regular power-series expansions near the axis. We have studied the vacuum limit carefully as a model for the force-free and MHS cases. We have shown that near the axis, we can transform to variables that allow us to describe the force-free and MHS systems analogous to the vacuum system.

We have demonstrated that one can carry out the NAE formally to all orders by following the same steps. The zeroth step consists of solving a MDE for the pressure-driven part of the current. We have shown how using the leading-order Birkhoff normal form considerably simplifies this calculation. The first step is to solve a Poisson equation for a scalar representing the magnetic scalar potential in the vacuum limit. The second step is to solve another MDE to obtain the flux function up to a homogeneous solution of the MDE. We fix this homogeneous solution by requiring the periodicity of physical quantities. Finally, with the help of the flux function and the solution to the Poisson equation, we can obtain the field line label through a simple integral over the poloidal angle up to a function of only the toroidal variable. We use the poloidally averaged MDE for the field line label to fix this unknown. The procedure described above is repeated in each order. We have given the complete details of the first two orders of the NAE for vacuum fields in Appendix F followed by force-free and MHS fields in Appendix G.

We have illustrated the finite beta near-axis procedure obtained through the Mercier-Weitzer formalism for various problems. We have treated the case of vacuum fields with a nearly circular cross-section in detail in § 6.1. The force-free and finite beta analogues of this example are given in Appendix H.2. Furthermore, the special cases of axisymmetry with Solov'ev profiles and magnetic axis with constant torsion are discussed in § 6.2. The analytical details of these examples are to be found in Appendices H.3 and H.4.

The discussion of physical quantities, such as the magnetic well and the magnetic shear, which play a crucial role in MHS equilibrium and stability, is out of the scope of the present work. Analytical expression of the near-axis magnetic shear has already been obtained in inverse coordinates (Rodríguez et al. 2022), while the near-axis magnetic well has been studied in Landreman & Jorge (2020) and Kim et al. (2021). In subsequent work, we shall calculate the magnetic well and the shear in direct coordinates using the formalism developed in this work.

The calculations presented here, while lengthy, move forward the NAE stellarator optimization programme significantly by presenting reduced analytical models that can be automated in fast numerical optimization suites for stellarator design and for testing numerical codes.