2 Derivation

To begin, we expand the magnetic field around a vacuum magnetic field $\boldsymbol {\nabla }\chi$, where $\nabla ^2\chi = O(\epsilon ^2)$:

(2.1)\begin{equation} \boldsymbol{B} = \left(1 + \frac{B_1}{B_v}\right)\boldsymbol{\nabla}\chi + \boldsymbol{B}_{{\perp} 1} + O(\epsilon^2).\end{equation}

We will order relative to the $\boldsymbol {\nabla }\chi$ term: $B_v = |\boldsymbol {\nabla }\chi | = O(1)$, $B_1 = O(\epsilon )$ and $\boldsymbol {B}_{\perp 1} = O(\epsilon )$, where $\epsilon = \max |\boldsymbol {B} - \boldsymbol {\nabla }\chi |/B_v \ll 1$ and $\boldsymbol {B}_{\perp 1}\perp \boldsymbol {\nabla }\chi$. In addition, the derivative along the vacuum field must be ordered as $\boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla } = O(\epsilon )$, whereas $\boldsymbol {\nabla }_\perp = \boldsymbol {\nabla } - B_v^{-2}\boldsymbol {\nabla }\chi\boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla } = O(1)$. The derivative ordering can be justified as follows. The the perpendicular length scale will be $\sim r_0$, whereas the length scale along $\boldsymbol {\nabla }\chi$ will be $\sim 2{\rm \pi} R_0$, where $r_0$ and $R_0$ are the minor and major radii, respectively. Thus, the ratio of the length scales must be less than $1/2{\rm \pi}$, and no more than 0.1 for any realistic aspect ratio.

Taking the divergence of (2.1), we have $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B} = \boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla }(B_1/B_v) + \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {B}_{\perp 1} = {\boldsymbol {\nabla }_\perp \boldsymbol {\cdot }\boldsymbol {B}_{\perp 1}} + O(\epsilon ^2)$. Since $\boldsymbol {\nabla }_\perp$ is two-dimensional, a stream function $A = O(\epsilon )$ can be introduced: $\boldsymbol {B}_{\perp 1} = \boldsymbol {\nabla } A\times \boldsymbol {\nabla }\chi$. Note that Strauss uses the symbol $\psi$ for his stream function, but we will use $A$ instead to avoid confusion with the flux surface label.

Following Strauss, we can write the current density as follows:

(2.2)\begin{equation} \boldsymbol{j} = \frac{\boldsymbol{B}\times\boldsymbol{\nabla} p}{B^2} + j_\parallel\frac{\boldsymbol{B}}{B} = \frac{\boldsymbol{\nabla}\chi\times\boldsymbol{\nabla} p}{B_v^2} - \frac{\boldsymbol{B}}{\mu_0}\varDelta^* A + O(\epsilon^2),\end{equation}

where $\varDelta ^* = B_v^{-2}\boldsymbol {\nabla }\boldsymbol {\cdot }(B_v^2\boldsymbol {\nabla }_\perp )$, as defined by Strauss (Reference Strauss1997), and the expression for $j_\parallel$ is easily obtained by dotting the curl of (2.1) with $\boldsymbol {\nabla }\chi$ and using the identity $\boldsymbol {\nabla } f\boldsymbol {\cdot }\boldsymbol {\nabla }\times \boldsymbol {U} = -\boldsymbol {\nabla }\boldsymbol {\cdot }(\boldsymbol {\nabla } f\times \boldsymbol {U})$. An alternate expression can be obtained by taking the curl of (2.1):

(2.3)\begin{equation} \boldsymbol{j} = \frac{1}{\mu_0}\boldsymbol{\nabla}\left(\frac{B_1}{B_v}\right)\times\boldsymbol{\nabla}\chi - \frac{\boldsymbol{\nabla}\chi}{\mu_0}\nabla^2 A + \frac{1}{\mu_0}(\boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{\nabla} A - \frac{1}{\mu_0}(\boldsymbol{\nabla} A\boldsymbol{\cdot}\boldsymbol{\nabla})\boldsymbol{\nabla}\chi. \end{equation}

The last two terms are both $O(\epsilon ^2)$. This becomes more obvious for the last term if we write it as $(\boldsymbol {\nabla } A\boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {\nabla }\chi = \boldsymbol {\nabla }(\boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla } A) - (\boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla })\boldsymbol {\nabla } A$. Equating the perpendicular components of the current at lowest order, we obtain

(2.4)\begin{equation} \left[\boldsymbol{\nabla}\left(\frac{B_1}{\mu_0 B_v}\right) + \frac{1}{B_v^2}\boldsymbol{\nabla} p\right] \times\boldsymbol{\nabla}\chi = 0.\end{equation}

Unless the second term is a gradient, both terms are linearly independent and must be individually zero at order $\epsilon$ (see Appendix A for details). If the second term is a gradient, we must have either $B_v = B_v(p)$ or $B_v = B_0 + O(\epsilon )$ where $B_0 = \mathrm {const}$; however, the former cannot be true as $|\boldsymbol {B}|$ can only be a flux function in axisymmetry (Schief Reference Schief2003), and since $|\boldsymbol {B}| = B_v + O(\epsilon ^2)$, $B_v$ cannot be a flux function either. Thus, to satisfy the equation, we must either ensure that both terms are individually zero, i.e. $p = O(\epsilon ^2)$ and $B_1 = 0$, or let $B_v = B_0 + O(\epsilon )$. The former corresponds to the equilibrium limit of the Strauss equations, and so we focus on the latter. The equation then yields:

(2.5)\begin{equation} \frac{B_0 B_1}{\mu_0} + p = 0.\end{equation}

This is simply the lowest-order radial pressure balance condition; similar expressions appear in most reduced MHD models that assume $p = O(\epsilon )$ (Freidberg Reference Freidberg2014; Zocco, Helander & Weitzner Reference Zocco, Helander and Weitzner2020; Kruger, Hegna & Callen Reference Kruger, Hegna and Callen1998). Taking the divergence of (2.2) then gives the generalized version of Strauss’ (26) in the equilibrium limit:

(2.6)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\left(\frac{j_\parallel}{B}\right) + \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\frac{\boldsymbol{B}\times\boldsymbol{\nabla} p}{B^2}\right) = \frac{-1}{\mu_0}\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\varDelta^* A + \boldsymbol{\nabla}\left(\frac{1}{B_v^2}\right) \boldsymbol{\cdot}(\boldsymbol{\nabla}\chi\times\boldsymbol{\nabla} p) + O(\epsilon^3) = 0.\end{equation}

The only other equation of Strauss’ reduced MHD that is non-trivial in the equilibrium limit can be written simply as $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } p = 0$.

Equation (2.6) can be further simplified by exploiting the two-term quasisymmetry condition, the imposition of which is equivalent to demanding that $|\boldsymbol {B}| = B(\psi, M\theta - N\phi )$ (Helander Reference Helander2014): $(\boldsymbol {B}\times \boldsymbol {\nabla }\psi )\boldsymbol {\cdot }\boldsymbol {\nabla } B = F(\psi )\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } B$, where $\psi$ is the toroidal flux. Applying the ordering, this condition can be written in one of two equivalent ways:

(2.7a)\begin{gather} (\boldsymbol{\nabla}\chi\times\boldsymbol{\nabla}\varPsi)\boldsymbol{\cdot}\boldsymbol{\nabla} B_v = \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} B_v, \end{gather}

(2.7b)\begin{gather}(\boldsymbol{\nabla} B_v\times\boldsymbol{\nabla}\chi)\boldsymbol{\cdot}\boldsymbol{\nabla}(\varPsi + A) = \boldsymbol{\nabla}\chi\boldsymbol{\cdot}\boldsymbol{\nabla} B_v, \end{gather}

where an alternate flux surface label $\varPsi = \int d\psi /F(\psi )$ is defined to absorb the $F(\psi )$ factor. Note that $F(\psi )$ is $O(\epsilon ^{-1})$ (Helander Reference Helander2014), so $\varPsi = O(\epsilon )$. Using (2.7a), the second term in (2.6) can be rewritten as $p'(\boldsymbol {\nabla }\chi \times \boldsymbol {\nabla }\varPsi )\boldsymbol {\cdot }\boldsymbol {\nabla } B_v^{-2} = p'\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } B_v^{-2}$, which results in both terms in the equation having $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }$ acting on something. Since $p'$ commutes with $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }$, we can remove the $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }$ operator, obtaining a Grad–Shafranov-like equation:

(2.8)\begin{equation} \varDelta^* A - \mu_0\frac{{\rm d} p}{{\rm d}\varPsi}\left(\frac{1}{B_v^2} - \frac{1}{B_0^2}\right) = H(\varPsi),\end{equation}

where $H$ and $\mu _0 p'/B_0^2$ are functions of $\varPsi$ that appear due to the integration. Although $H$ is arbitrary, the $B_0$ term cannot be absorbed into it as it is necessary to cancel the $O(1)$ part of $1/B_v^2$ so that all terms are of the same order. Further, note that $H$ is a degree of freedom that allows us to specify the toroidal current profile; its axisymmetric equivalent is the $FF'$ term in the standard Grad–Shafranov equation.

To make further progress, $\chi$ has to be specified explicitly. The constraint $B_v = B_0 + O(\epsilon )$, which, as we have seen, arises from the high-$\beta$ assumption, requires that $\chi = B_0 l$. This can be seen by writing $B_v^2$ in orthogonal Mercier coordinates $(\rho,\omega,l)$ aligned to a magnetic axis with curvature $\kappa$ and torsion $\tau$ (Solov'ev & Shafranov Reference Solov'ev and Shafranov1970; Jorge et al. Reference Jorge, Sengupta and Landreman2020):

(2.9)\begin{equation} B_v^2 = \left(\frac{\partial\chi}{\partial\rho}\right)^2 + \frac{1}{\rho^2}\left(\frac{\partial\chi}{\partial\omega}\right)^2 + \frac{1}{h^2}\left(\frac{\partial\chi}{\partial l}\right)^2 = B_0^2 + O(\epsilon), \end{equation}

where $h = 1 - \kappa \rho \cos \theta$, $\theta = \omega - \int _0^l\tau dl$. On the axis, we must have $\partial \chi /\partial \rho |_{\rho =0} = \lim _{\rho \to 0}\rho ^{-1}\partial \chi /\partial \omega = 0$, since the axis is a field line. Thus, the constant $B_0^2$ can only come from the last term. It then follows that $\kappa \rho = O(\epsilon )$, so that $1/h^2 = 1 + O(\epsilon )$; otherwise, $\chi$ would have depended on $\rho$ and $\omega$ at the lowest order, which is a contradiction. With that in mind, we can integrate at the zeroth order, obtaining $\chi = B_0 l + O(\epsilon )$. We can ignore the $O(\epsilon )$ correction to $\chi$ without loss of generality; as will be shown at the end of the section, it can be absorbed into the $\boldsymbol {\nabla } A\times \boldsymbol {\nabla }\chi$ term. Note that we can now give $\epsilon$ a more intuitive meaning: given that $\rho \sim r_0$ and $\kappa,\tau \sim 1/R_0$, we should order $r_0 = O(1)$ and $R_0 = O(\epsilon )$, and hence we have a large aspect ratio expansion with $\epsilon$ acting as the inverse aspect ratio. Finally, it should be emphasized that the only arbitrary assumption made in this derivation is that $p = O(\epsilon )$. The rest of the ordering, while not rigorous, is still justified by heuristic arguments.

As a sanity check, we will demonstrate that $\chi = B_0 l$ only breaks the divergence-free condition at second order, which is acceptable, since this is a first-order model. Using the Laplace operator in orthogonal Mercier coordinates (see (3.4) of Sengupta et al. (Reference Sengupta, Rodriguez, Jorge, Landreman and Bhattacharjee2024b)), we have

(2.10)\begin{align} \nabla^2\chi = \frac{1}{\rho h}\frac{\partial}{\partial\rho} \left(\rho h\frac{\partial\chi}{\partial\rho}\right) + \frac{1}{\rho^2 h} \frac{\partial}{\partial\omega}\left(h\frac{\partial\chi}{\partial\omega}\right) + \frac{1}{h}\frac{\partial}{\partial l}\left(\frac{1}{h}\frac{\partial\chi}{\partial l}\right) ={-}\frac{B_0}{h^3}\frac{\partial h}{\partial l} = O(\epsilon^2), \end{align}

since $h = 1 + O(\epsilon )$ and $\partial /\partial l = B_0^{-1}h^2\boldsymbol {\nabla }\chi \boldsymbol {\cdot }\boldsymbol {\nabla } = O(\epsilon )$.

Throughout the rest of the paper, we will use non-orthogonal Mercier coordinates $(x,y,l)$, where $x = \rho \cos \theta$ and $y = \rho \sin \theta$. The metric tensor is as follows (Solov'ev & Shafranov Reference Solov'ev and Shafranov1970):

(2.11)\begin{equation} g^{ik} = \frac{1}{h^2}\begin{pmatrix} h^2 + \tau^2 y^2 & -\tau^2 xy & \tau y \\ -\tau^2 xy & h^2 + \tau^2 x^2 & -\tau x \\ \tau y & -\tau x & 1 \end{pmatrix},\quad h = 1 - \kappa x. \end{equation}

Using the metric tensor and $\chi = B_0 l$, which implies $B_v = B_0/h$, (2.7b) can be written as

(2.12)\begin{equation} B_0^2\frac{\partial}{\partial y}(\varPsi + A) ={-}B_0^2\left(\tau y + \frac{1}{\kappa} \frac{{\rm d}\kappa}{{\rm d} l}x\right). \end{equation}

The non-orthogonal Mercier coordinates have simplified the equation such that it can now be integrated:

(2.13)\begin{equation} \varPsi + A ={-}\frac{\tau y^2}{2} - \frac{1}{\kappa}\frac{{\rm d}\kappa}{{\rm d} l}xy + a(x,l),\end{equation}

where $a$ is an arbitrary function. This allows us to write (2.8) as an equation for $\varPsi$:

(2.14)\begin{equation} \varDelta_\perp\varPsi + \tau - \frac{\partial^2 a}{\partial x^2} = \frac{2\mu_0}{B_0^2}\frac{{\rm d} p}{{\rm d}\varPsi}\kappa x - H(\varPsi),\end{equation}

where we have used $B_v = B_0|\boldsymbol {\nabla } l| = B_0/(1 - \kappa x)$ and $\varDelta ^* = \varDelta _\perp + O(\epsilon )$, and dropped terms of order $\epsilon ^2$ and higher.

Although we now have an equation for $\varPsi$, solving it is non-trivial, as $\varPsi$ is overconstrained, due to also having to satisfy the equation $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\varPsi = 0$, which, after inserting $A$ from (2.13), can be written in Mercier coordinates as

(2.15)\begin{equation} \frac{\partial\varPsi}{\partial l} - \frac{1}{\kappa}\frac{{\rm d}\kappa}{{\rm d} l}x \frac{\partial\varPsi}{\partial x} + \left(\frac{1}{\kappa}\frac{{\rm d}\kappa}{{\rm d} l}y- \tau x - \frac{\partial a}{\partial x}\right)\frac{\partial\varPsi}{\partial y}= 0,\end{equation}

where we have inserted the expression for $A$ from (2.13). The characteristics of this equation are given by

(2.16a)\begin{equation} \kappa x = C_1;\quad \frac{y}{\kappa} ={-} C_1\int\frac{\tau}{\kappa^2}{\rm d} l - \int\left.\frac{1}{\kappa}\frac{\partial a}{\partial x}\right|_{x=C_1/\kappa}{\rm d} l + C_2.\end{equation}

For $\varPsi$ to satisfy (2.15), it must be a function of only $C_1$ and $C_2$, as defined in (2.16a,b). One must then ensure that the function also satisfies (2.14) at all values of $l$.

Having derived the model, we now show that adding an $O(\epsilon )$ correction to $\chi$ will not change it. Thus, if the assumption that the $B_0 l$ is the dominant term in $\chi$ holds, the model is fully general. Indeed, suppose that we add a correction $\chi _1 = O(\epsilon )$ to $\chi$, then at order $\epsilon$, this term can be absorbed into the $\boldsymbol {\nabla } A\times \boldsymbol {\nabla }\chi$ term by replacing $A\mapsto A + \tilde {A}$, where $\boldsymbol {\nabla }\chi _1 = B_0\boldsymbol {\nabla }\tilde {A}\times \boldsymbol {\nabla } l$. The fact that such an $\tilde {A}$ exists can be seen by writing out the contravariant $x$ and $y$ components of this relation:

(2.17a)\begin{equation} \frac{\partial\chi_1}{\partial x} = B_0\frac{\partial\tilde{A}}{\partial y} + O(\epsilon^2),\quad \frac{\partial\chi_1}{\partial y} ={-}B_0\frac{\partial\tilde{A}}{\partial x} + O(\epsilon^2). \end{equation}

These are just the Cauchy–Riemann equations; since $\chi _1$ must satisfy the Laplace equation at order $\epsilon$, a corresponding $\tilde {A}$ is guaranteed to exist.

3 Consistency with the first-order near-axis model

The simplest ansatz for $\varPsi$ that will work in this model is a quadratic polynomial in $C_1$ and $C_2$, which are defined in (2.16a,b):

(3.1)\begin{align} \varPsi & = s_1 C_1^2 + s_2 C_1 C_2 + s_3 C_2^2 = s_1 (\kappa x)^2 + s_2 (\kappa x)\left(\frac{y}{\kappa} + \kappa x\varLambda\right) + s_3 \left(\frac{y}{\kappa} + \kappa x\varLambda\right)^2 \nonumber\\ & = (s_1 + s_2\varLambda + s_3\varLambda^2)(\kappa x)^2 + (s_2 + 2s_3\varLambda)xy + s_3\left(\frac{y}{\kappa}\right)^2, \end{align}

where $\varLambda = \int \kappa ^{-2}(\tau + a_2)\,{\rm d} l$; and $a$ was chosen as $a = a_2(l)x^2/2$. This represents a plasma with an elliptic cross-section that rotates about the magnetic axis as $l$ varies. Inserting this expression into the Grad–Shafranov equation (2.14), we obtain an equation expressing $a_2$ in terms of $\varLambda$. Note that, when inserting (3.1) into the Grad–Shafranov equation, all terms on the left-hand side of (2.14) will depend only on $l$, whereas terms on the right-hand side can depend on all three variables. Thus, we further assume $p = O(\epsilon ^2)$ and $H(\varPsi ) = H_0 = \mathrm {const}$; this is consistent with the near-axis model, where pressure enters only at the second order. Proceeding, we have

(3.2)\begin{equation} 2(s_1 + s_2\varLambda + s_3\varLambda^2)\kappa^2 + 2\frac{s_3}{\kappa^2} + \tau - a_2 + H_0 = 0.\end{equation}

Finally, we can combine ${\rm d}\varLambda /{\rm d} l = (\tau + a_2)/\kappa ^2$, which follows from the definition of $\varLambda$, with the above equation, resulting in a Riccati equation for $\varLambda$:

(3.3)\begin{equation} \frac{{\rm d}\varLambda}{{\rm d} l} - 2(s_1 + s_2\varLambda + s_3\varLambda^2) - \frac{2s_3}{\kappa^4} - \frac{2\tau}{\kappa^2} - \frac{H_0}{\kappa^2} = 0.\end{equation}

The last step is to enforce periodicity on $\varLambda$. Averaging the above equation over $l$ removes the derivative term, resulting in the following constraint:

(3.4)\begin{equation} -2(s_1 + s_2\langle\varLambda\rangle + s_3\langle\varLambda^2\rangle) - 2s_3 \left\langle\frac{1}{\kappa^4}\right\rangle - 2\left\langle\frac{\tau}{\kappa^2}\right\rangle - H_0\left\langle\frac{1}{\kappa^2}\right\rangle = 0,\end{equation}

where $\langle f(l)\rangle = L^{-1}\int _0^L f(l)\mathrm {d}l$, with $L$ being the axis length. Thus, only two of the $s_i$ constants are free, with the remaining one determined by the above constraint. In practice, since $\varLambda$ is not known a priori, (3.3) and (3.4) must be solved iteratively: first, one makes an initial guess for the unknown $s_i$, then (3.3) is solved for $\varLambda$, which allows one to calculate the corrected value of the unknown $s_i$ from (3.4); this cycle is repeated until convergence is achieved. This same method is used in the near-axis model to solve the $\sigma$-equation while simultaneously finding the correct rotational transform that is compatible with periodicity (Landreman, Sengupta & Plunk Reference Landreman, Sengupta and Plunk2019).

Now consider the lowest order near-axis expression for $\varPsi$ (Jorge et al. Reference Jorge, Sengupta and Landreman2020):

(3.5)\begin{align} \varPsi & \simeq \frac{\psi}{F_0} \simeq \frac{{\rm \pi} B_0}{F_0}\rho^2 [\cosh\eta + \sinh\eta\cos 2(\theta + \delta)] \nonumber\\ & = \frac{{\rm \pi} B_0}{F_0}[(\cosh\eta + \sinh\eta\cos 2\delta)x^2 - 2\sinh\eta\sin 2\delta~xy +(\cosh\eta - \sinh\eta\cos 2\delta)y^2], \end{align}

where $x = \rho \cos \theta$, $y = \rho \sin \theta$, $F_0$ is the value of $F(\psi )$ on axis and $\eta$ and $\delta$ are functions of $l$. Note that $F_0$, as defined here, is the inverse of the $F_0$ of Jorge et al. (Reference Jorge, Sengupta and Landreman2020). In quasisymmetry, we have $\cosh \eta - \sinh \eta \cos 2\delta = \bar {\eta }^2/\kappa ^2$ (Jorge et al. Reference Jorge, Sengupta and Landreman2020). Since $\bar {\eta } = \mathrm {const}$Footnote ¹ , the $y^2$ term of (3.1) agrees with that of (3.5) and $s_3 = {\rm \pi}B_0\bar {\eta }^2/F_0$. Further, taking $4(\text {coef. of }x^2)(\text {coef. of }y^2) - (\text {coef. of } xy)^2$ from (3.5), we see that the $l$-dependent terms cancel and we are left with just a constant, $4{\rm \pi} ^2 B_0^2/F_0^2$. Likewise, when using the corresponding coefficients from (3.1), the terms with $\varLambda$ cancel and we get $4s_1 s_3 - s_2^2 = 4{\rm \pi} ^2 B_0^2/F_0^2$.

Next, we show that (3.3) is equivalent to the $\sigma$-equation of Jorge et al. (Reference Jorge, Sengupta and Landreman2020). A relation between $\varLambda$ and $\sigma$ can be obtained by comparing the coefficients of $xy$ in (3.1) and (3.5). Since $\sigma = \sinh \eta \sin 2\delta$ (Jorge et al. Reference Jorge, Sengupta and Landreman2020), we have

(3.6)\begin{equation} \varLambda ={-}\frac{{\rm \pi} B_0}{F_0 s_3}\sigma - \frac{s_2}{2s_3}. \end{equation}

Inserting this into (3.3) and replacing the $s_i$ variables with near-axis variables, we get

(3.7)\begin{equation} \frac{{\rm d}\sigma}{{\rm d} l} + \frac{2{\rm \pi} B_0}{F_0} \left(1 + \sigma^2 + \frac{\bar{\eta}^4}{\kappa^4}\right) + 2\tau\frac{\bar{\eta}^2}{\kappa^2} + H_0\frac{\bar{\eta}^2}{\kappa^2} = 0.\end{equation}

In the case where $H_0 = 0$, this equation reduces to (29) of Jorge et al. (Reference Jorge, Sengupta and Landreman2020).

To conclude this section, we show that when $s_2 = 0$, there is a simple relationship between the characteristics and Boozer angles. Choosing a value for $s_2$ is similar to gauge fixing, since a non-zero $s_2$ simply represents a deviation of the ellipse from the upright position in the $(C_1,C_2)$ plane. The orientation of the ellipse in real space is determined by $\sigma$, which is unaffected by changes in $s_2$, as long as the initial condition for $\varLambda$ is changed accordingly. The coordinates $x$ and $y$ can be written in terms of $C_1,C_2$ as

(3.8a)\begin{equation} x = \frac{C_1}{\kappa},\quad y = \kappa C_2 - \kappa\varLambda C_1. \end{equation}

Meanwhile, (1) and (54) from Jorge et al. (Reference Jorge, Sengupta and Landreman2020) give

(3.9)\begin{equation} \boldsymbol{r} = \boldsymbol{r}_0 + x\boldsymbol{n} + y\boldsymbol{b} = \boldsymbol{r}_0 + \sqrt{\frac{\psi}{{\rm \pi} B_0}}\left[\frac{\bar{\eta}}{\kappa} \cos\vartheta\boldsymbol{n} + \frac{\kappa}{\bar{\eta}}(\sin\vartheta + \sigma\cos\vartheta)\boldsymbol{b}\right]. \end{equation}

Equating the two and using the relations between $\varLambda$ and $\sigma$ as well as $s_3$ and $\bar {\eta }$, we arrive at the following:

(3.10a)\begin{equation} C_1 = \frac{F_0}{{\rm \pi} B_0}\sqrt{s_3\varPsi}\cos\vartheta,\quad C_2 = \sqrt{\frac{\varPsi}{s_3}}\sin\vartheta, \end{equation}

where $\vartheta = \theta _B - N\phi _B$, with $\theta _B$ and $\phi _B$ being the Boozer angles.

4 New classes of solutions

We now present three classes of exact and approximate solutions to (2.14) and (2.15) that fall outside the scope of the near-axis model. We will begin with a cubic polynomial solution with constant $H$ and then proceed to a non-polynomial solution and a solution with variable $H$.

4.1 Cubic solution with finite pressure

Consider the following cubic polynomial in $C_1$ and $C_2$:

(4.1)\begin{align} \varPsi & = s_1 C_1^2 + s_2 C_1 C_2 + s_3 C_2^2 + r C_1^3 \nonumber\\ & = (s_1 + s_2\varLambda + s_3\varLambda^2)(\kappa x)^2 + (s_2 + 2s_3\varLambda)xy + s_3\left(\frac{y}{\kappa}\right)^2 + r(\kappa x)^3. \end{align}

Inserting this ansatz into (2.14), the following is obtained:

(4.2)\begin{equation} 2(s_1 + s_2\varLambda + s_3\varLambda^2)\kappa^2 + \frac{2s_3}{\kappa^2} + 6r\kappa^3 x + \tau - a_2 = \frac{2\mu_0 p_1}{B_0^2}\kappa x - H_0,\end{equation}

where we have additionally assumed that $dp/d\varPsi = p_1 = \mathrm {const}$, $H(\varPsi ) = H_0 = \mathrm {const}$ and $a = a_2(l)x^2/2$. Only two terms in the above equation depend on $x$; equating those two terms, we get

(4.3)\begin{equation} r = \frac{\mu_0 p_1}{3B_0^2\kappa^2},\end{equation}

while the remaining terms match (3.2). Thus, a solution is obtained by first finding a quadratic solution, as discussed in § 3, and then adding a cubic term $r(\kappa x)^3$ with $r$ given by (4.3). Equation (4.3) also imposes the constraint that $\kappa = \mathrm {const}$, but, since $\tau$ is unconstrained, we can still get non-planar closed curves. At this point, we can see why $rC_1^3$ is the only cubic term that we can include in (4.1). If we were to add cubic terms that involve $C_2$, then (4.2) would also have terms that depend on $y$ and we would end up with an additional equation that would constrain $\tau$, meaning that we would not be able to ensure that the axis is closed.

Finally, note that while this cubic solution seems to be superficially similar to a second-order near-axis solution, there is a fundamental difference. Namely, the near-axis model assumes that the cubic term is a correction that is much smaller than the quadratic terms, whereas the present solution allows the cubic term to be of the same order as the quadratic ones. Also, unlike the second-order near-axis model, where the quasisymmetry error is $O(\epsilon ^3)$, the present solution is still first order, so the quasisymmetry error will be $O(\epsilon ^2)$.

4.2 Non-polynomial approximate solutions

In this subsection, we will show an approximate solution that consists of a rotating ellipse, which, as we have seen in the previous section, can be represented as a quadratic polynomial, and a non-polynomial perturbation. As we will perform a subsidiary expansion, it is convenient to rescale all quantities to be zeroth order in $\epsilon$. Thus, after performing the following replacement: $\{\kappa,\bar {\eta },\tau,a_2,H_0,\varPsi,{\rm d}/{\rm d} l\}\mapsto \epsilon \{\kappa,\bar {\eta },\tau,a_2,H_0,\varPsi,{\rm d}/{\rm d} l\}$ and $F_0\mapsto F_0/\epsilon$, we see that the $\epsilon$ parameter is cancelled in (2.14) and the rotating ellipse solution (3.1).

Now consider the case when $s_2 = 0$; then, using the expressions for the $s_i$ variables and $\varLambda$ that we found in § 3, we have $s_3 = {\rm \pi}B_0\bar {\eta }^2/F_0$, $s_1 = {\rm \pi}B_0/(F_0\bar {\eta }^2)$ and $\varLambda = -\sigma /\bar {\eta }^2$. The rotating ellipse solution can be represented as

(4.4)\begin{equation} \varPsi_{re} = \frac{{\rm \pi} B_0}{F_0}\left(\frac{1}{\bar{\eta}^2}C_1^2 + \bar{\eta}^2 C_2^2\right) = \frac{{\rm \pi} B_0}{F_0}(c_1^2 + c_2^2),\end{equation}

where the integration constants of the characteristics (2.16a,b) have been rescaled with respect to $\bar {\eta }$:

(4.5a)\begin{equation} c_1 = \frac{C_1}{\bar{\eta}} = \frac{\kappa x}{\bar{\eta}},\quad c_2 = \bar{\eta} C_2 = \frac{\bar{\eta} y}{\kappa} + \bar{\eta}\kappa x\varLambda = \frac{\bar{\eta} y}{\kappa} - \frac{\kappa x\sigma}{\bar{\eta}} \end{equation}

yielding the normal form representation. We add to (4.4) a non-polynomial perturbation $f$, and use the relation $F_0^{-1} = (\iota _0 - N)/(2{\rm \pi} B_0 r_0)$, which was found by Jorge et al. (Reference Jorge, Sengupta and Landreman2020), with $\iota _0$ being the rotational transform on axis and $N$ the helicity of the quasisymmetry. Here, the axis length $L$ was replaced with the characteristic value of the minor radius $r_0 = \epsilon L/(2{\rm \pi} )$, due to $F_0$ having been rescaled with respect to $\epsilon$. The resulting ansatz is as follows:

(4.6)\begin{equation} \varPsi = \frac{\iota_0 - N}{2r_0}\left[(c_1^2 + c_2^2) + \left(\frac{\bar{\eta}}{\kappa_0}\right)^2 f(c_2 + \sigma_0 c_1)\right],\end{equation}

where $\kappa _0$ is the minimum value of $\kappa$. Inserting this ansatz into the Grad–Shafranov equation (2.14), multiplying everything by $(\bar {\eta }/\kappa )^2$ and using $a_2 = \kappa ^2 \,{\rm d}\varLambda /{\rm d} l - \tau = -(\kappa /\bar {\eta })^2 \,{\rm d}\sigma /{\rm d} l - \tau$, the following is obtained:

(4.7)\begin{align} & \frac{\iota_0 - N}{2r_0}\left[2(1 + \sigma^2) + 2\left(\frac{\bar{\eta}}{\kappa}\right)^4 + \left(\frac{\bar{\eta}}{\kappa_0}\right)^2\left((\sigma_0 - \sigma)^2 + \left(\frac{\bar{\eta}}{\kappa}\right)^4\right)f''\right] \nonumber\\ & \quad + \frac{{\rm d}\sigma}{{\rm d} l} + \left(\frac{\bar{\eta}}{\kappa}\right)^2(2\tau + H_0) = 0. \end{align}

We can now carry out a subsidiary expansion in the limit of small $\bar {\eta }^2$, i.e. $(\bar {\eta }/\kappa _0)^2 \ll 1$, but $(\bar {\eta }/\kappa _0)^4 > \epsilon$ since we keep terms at that order in (4.6). This corresponds to the limit of a high aspect ratio stellarator with a highly elongated cross-section. Given that $\sigma \sim \bar {\eta }^2$ and $\iota _0 - N \sim \bar {\eta }^2$ (Rodriguez Reference Rodriguez2022), the two lowest orders of the above equation are as follows:

(4.8)\begin{equation} \left.\begin{gathered} O\left(\left(\frac{\bar{\eta}}{\kappa_0}\right)^2\right): \quad \frac{{\rm d}\sigma_2}{{\rm d} l} + \frac{\iota_{0,2}}{r_0} ={-}\left(\frac{\bar{\eta}}{\kappa}\right)^2(2\tau + H_0), \\ O\left(\left(\frac{\bar{\eta}}{\kappa_0}\right)^6\right): \quad \frac{{\rm d}\sigma_6}{{\rm d} l} + \frac{\iota_{0,6}}{r_0} ={-}\frac{\iota_{0,2}}{r_0}\left[\sigma_2^2 + \left(\frac{\bar{\eta}}{\kappa}\right)^4\right], \end{gathered}\right\} \end{equation}

where $\sigma = \sigma _2 + \sigma _6$ with $\sigma _n\sim \bar {\eta }^n$, and $\iota _0 = N + \iota _{0,2} + \iota _{0,6}$ with $\iota _{0,n}\sim \bar {\eta }^n$. The values of $\iota _{0,n}$ are determined by enforcing periodicity on $\sigma _n$: (4.8) is averaged over $l$, which removes $d\sigma _n/dl$; the $\iota _{0,n}$ term must then be equal to the average of the right-hand side. Since we only solve (4.7) up to $O((\bar {\eta }/\kappa _0)^6)$ and $f$ does not appear until $O((\bar {\eta }/\kappa _0)^8)$, we can treat $f$ as a free function. If we were to attempt to solve this equation at $O((\bar {\eta }/\kappa _0)^8)$, then only the trivial solution $f'' = \mathrm {const}$ (i.e. $f$ is a quadratic polynomial) would be permitted.

4.3 Approximate solution with variable $H$

We conclude this section by presenting a rotating ellipse approximate solution where $H(\varPsi ) = H_0 + H_1\varPsi$. We also add a quartic term to $a$: $a = a_2(l)x^2/2 + a_4(l)x^4/24$, so the rescaled characteristics become

(4.9)\begin{equation} \left.\begin{gathered} c_1 = \frac{\kappa x}{\bar{\eta}}, \quad c_2 = \frac{\bar{\eta} y}{\kappa} + \bar{\eta}\kappa x\varLambda + \frac{\bar{\eta}\kappa^3 x^3}{6}V, \\ \varLambda = \int\frac{\tau + a_2}{\kappa^2}{\rm d} l, \quad V = \int\frac{a_4}{\kappa^4}{\rm d} l. \end{gathered}\right\} \end{equation}

Just like in the previous subsection, we rescale all quantities to be zeroth order in $\epsilon$ and let $s_2 = 0$. In addition, we order $H_1\sim \bar {\eta }^2$ and $a_4\sim \bar {\eta }^2$. It can now be seen that in the rotating ellipse ansatz (4.4), terms with $x$ or $y$ raised to a power higher than two only appear at order $(\bar {\eta }/\kappa _0)^6$ and higher. Thus, up to $O((\bar {\eta }/\kappa _0)^4)$, it is still purely a rotating ellipse.

Inserting the rotating ellipse solution (4.4) with the new characteristics (4.9) into the Grad–Shafranov equation (2.14) and multiplying everything by $(\bar {\eta }/\kappa )^2$, the following equation is obtained:

(4.10)\begin{align} & \frac{\iota_0 - N}{r_0}\left[(1 + \sigma^2) + \left(\frac{\bar{\eta}}{\kappa}\right)^4\right] + \frac{{\rm d}\sigma}{{\rm d} l} + \left(\frac{\bar{\eta}}{\kappa}\right)^2\left(2\tau - a_4\frac{x^2}{2} + H_0\right) \nonumber\\ & \quad + \frac{\iota_0 - N}{2r_0}H_1 x^2 + O\left(\left(\frac{\bar{\eta}}{\kappa_0}\right)^8\right) = 0, \end{align}

where we have again used $a_2 = -(\kappa /\bar {\eta })^2 \,\textrm {d}\sigma /\textrm {d} l - \tau$. Just as in the previous subsection, the above equation can be solved order by order. The three lowest orders are as follows:

(4.11)\begin{equation} \left.\begin{gathered} O\left(\left(\frac{\bar{\eta}}{\kappa_0}\right)^2\right): \quad \frac{{\rm d}\sigma_2}{{\rm d} l} + \frac{\iota_{0,2}}{r_0} ={-}\left(\frac{\bar{\eta}}{\kappa}\right)^2(2\tau + H_0), \\ O\left(\left(\frac{\bar{\eta}}{\kappa_0}\right)^4\right): \quad \left(\frac{\bar{\eta}}{\kappa}\right)^2 a_4 = \frac{\iota_{0,2}}{r_0}H_1, \\ O\left(\left(\frac{\bar{\eta}}{\kappa_0}\right)^6\right): \quad \frac{{\rm d}\sigma_6}{{\rm d} l} + \frac{\iota_{0,6}}{r_0} ={-}\frac{\iota_{0,2}}{r_0}\left[\sigma_2^2 + \left(\frac{\bar{\eta}}{\kappa}\right)^4\right]. \end{gathered}\right\} \end{equation}

Thus, $\sigma$ is still determined by the same equations as (4.8); in addition to that, we also have an equation to determine $a_4$.

5 Numerical example

While the model presented in this paper was derived with high-$\beta$ equilibria in mind, the only high-$\beta$ equilibrium that we have been able to find analytically is the constant curvature one presented in § 4.1. In this section, we will instead illustrate the approximate solution discussed in § 4.2 with a numerical example of a quasiaxisymmetric ($N=0$) device. This approximate solution has another important property that we wanted our model to include: non-zero shear. We were unable to get a numerical equilibrium with high $\beta$ because VMEC failed to find a solution for all of the non-planar constant curvature axes that we could find, due to their high shaping. We leave the numerical verification of the analytical solution in § 4.1, as well as the numerical search for other high-$\beta$ solutions for future work.

The example is constructed by numerically solving (4.8), and is then compared with both VMEC (Hirshman & Whitson Reference Hirshman and Whitson1983) results and the closest near-axis approximation, which is identical to taking $f=0$. The pyQSC code (Reference Landreman, Jorge, Rodriguez and DudtLandreman et al. 2020–2023) was used to solve the $\sigma$ equations and generate VMEC input files, based on the method described by Landreman & Sengupta (Reference Landreman and Sengupta2019).

While the solution in § 4.2 allows for arbitrary $f$, here we will consider a sixth-order polynomial, $f(c) = (0.3\,\mathrm {m}^{-2})c^4 + (0.3\,\mathrm {m}^{-4})c^6$, where $c = c_2 + \sigma _0 c_1$; this is higher order than the cubic polynomials of the second-order near-axis model. We consider a four-field-period solution with $\bar {\eta } = 2.01735426\times 10^{-2}\,\mathrm {m}^{-1}$, $B_0 = 1$ T, $H_0 = 0$. Equations (4.8) result in rotational transforms $\iota _{0,2} = 0.188923174$ and $\iota _{0,6} = -0.031115615$. The axis shape given by

(5.1)\begin{align} \left.\begin{gathered} R(\phi) = 29.7794783 - 3.63597602\times 10^{{-}1}\cos 4\phi + 1.47477208\times 10^{{-}1}\cos 8\phi \\ \quad + 1.35576435\times 10^{{-}2}\cos 12\phi, \\ z(\phi) = 1.93173817\sin 4\phi + 2.38762327\times 10^{{-}2}\sin 8\phi - 7.72243217 \times 10^{{-}3}\sin 12\phi, \end{gathered}\right\} \end{align}

where $R$ and $z$ are the cylindrical coordinates in metres. The value of $(\bar {\eta }/\kappa _0)^2$ for this axis is approximately 0.3. The VMEC equilibrium, obtained by constructing the $F_0\varPsi /{\rm \pi} = 4\,\mathrm {T}\,\textrm {m}^2$ surface using (4.6) and passing it to VMEC as the boundary, is shown in figure 1. The aspect ratio of the resulting stellarator is 16.5. The assumption $\textrm {d} p/\textrm {d}\varPsi = 0$, made when constructing the base solution in § 3, and the choice of $H_0$ correspond to pressure and toroidal current profiles of zero; this was specified in the VMEC input. The total toroidal flux in the VMEC input was determined by multiplying $B_0$ by the plasma cross-section area in the plane perpendicular to the axis.

Figure 1. The outermost flux surface of the VMEC equilibrium based on the solution with a sixth-order polynomial perturbation. Colour represents $|\boldsymbol {B}|$.

Figure 2(a,b) compares the flux surfaces of the present Grad–Shafranov model, obtained from (4.6), to both the flux surfaces calculated by VMEC and the flux surfaces of the closest near-axis approximation, as given by (4.4). All flux surfaces are shown on the $\phi = 0$ poloidal plane. As expected, near the axis, the near-axis flux surfaces closely match those of the Grad–Shafranov model, but further away from the axis, the Grad–Shafranov surfaces become non-elliptic as the contribution from $f$ becomes non-negligible. Figure 2(c) compares the $\iota$ profiles in the near-axis and Grad–Shafranov models. These are computed numerically by first calculating the toroidal flux $\psi$ at each value of $\varPsi$, then finding $F = \textrm {d}\varPsi /\textrm {d}\psi$ and using the formula

(5.2)\begin{equation} F(\psi) = \frac{G(\psi) + NI(\psi)}{\iota(\psi) - N}, \end{equation}

which is given by Helander (Reference Helander2014) as an unnumbered equation. Again, these agree well with each other and the $\iota _0$ constant near the axis, but the Grad–Shafranov solution has a noticeable shear. The $\iota$ profiles calculated by VMEC (not shown here) do not match the Grad–Shafranov and near-axis profiles shown in figure 2(c), with the $\iota$ on axis in the corresponding VMEC solutions being greater by approximately 0.06 and 0.024, respectively. A similar mismatch with VMEC has been observed by Sengupta et al. (Reference Sengupta, Rodriguez, Jorge, Landreman and Bhattacharjee2024b); it likely appears because the boundary (and not the axis) is fixed in VMEC, so for a finite aspect ratio, the VMEC axis does not match the axis given by (5.1). There are also additional complications arising from VMEC having a coordinate singularity on the axis (Panici et al. Reference Panici, Conlin, Dudt, Unalmis and Kolemen2023a). Just as was reported by Sengupta et al. (Reference Sengupta, Rodriguez, Jorge, Landreman and Bhattacharjee2024b), the agreement between the $\iota$ on axis in VMEC and that predicted by the Grad–Shafranov model improves as the aspect ratio is increased. Finally, figure 3 shows that the maximum quasisymmetry error, defined as

(5.3)\begin{equation} \max_\psi \sqrt{\sum_{n\neq 0}\hat{B}_{n,m}(\psi)^2/\sum_{n,m}\hat{B}_{n,m}(\psi)^2}, \end{equation}

and calculated in VMEC equilibria with varying aspect ratios, scales as $\epsilon ^2$. Here, $\hat {B}_{n,m}(\psi )$ are the Fourier modes of $|\boldsymbol {B}|$ on flux surface $\psi$. This behaviour is to be expected, since we have only ensured quasisymmetry at order $\epsilon$.

Figure 2. (a) A comparison of the flux surfaces computed from (4.6) (lines) to the closest near-axis approximation (dots). (b) Flux surfaces computed from (4.6) (lines) compared with the flux surfaces computed by VMEC (dots). The flux surfaces are shown for $F_0\varPsi /{\rm \pi} = 0.1, 0.3, 0.5, 1, 2, 4\,\mathrm {T}\,\textrm {m}^2$. (c) $\iota$ profiles computed in the Grad–Shafranov model, the near-axis model and the constant $\iota _0$.

Figure 3. Maximum quasisymmetry error (black dots) scales as $\epsilon ^2$ (dashed blue line is $10\epsilon ^2$).