1. Introduction
Omnigeneity (Hall & McNamara Reference Hall and McNamara1975; Nührenberg & Zille Reference Nührenberg and Zille1988; Cary & Shasharina Reference Cary and Shasharina1997), zero-average radial drift of the guiding centre orbits, is essential for particle confinement in three-dimensional (3-D) magnetic fields (Helander Reference Helander2014). One route to omnigeneity is quasisymmetry (QS, Boozer Reference Boozer1983), in which contours of the magnetic field strength are flat in a special coordinate system known as Boozer angles. QS with magnetohydrostatic (MHS) force balance, assuming isotropic pressure, is known to be an overdetermined system (Garren & Boozer Reference Garren and Boozer1991b ). However, recent strides have shown that extremely precise QS can be obtained by numerical optimisation (Landreman & Paul Reference Landreman and Paul2022; Wechsung et al. Reference Wechsung, Landreman, Giuliani, Cerfon and Stadler2022), reinforcing the confidence that QS is a practical approach to plasma confinement. Because quasisymmetric equilibria have been found only numerically, our understanding of the qualitative behaviour of quasisymmetric equilibria is still inchoate. Considerable effort has been applied to the near-axis theory of QS (Garren & Boozer Reference Garren and Boozer1991 Reference Garren and Boozer a ,Reference Garren and Boozer b ; Landreman Reference Landreman2019; Rodríguez & Bhattacharjee, Reference Rodríguez and Bhattacharjee2021b ,Reference Rodríguez and Bhattacharjee c ) and to weakly non-axisymmetric equilibria at finite aspect ratio (Plunk & Helander Reference Plunk and Helander2018; Henneberg et al. Reference Henneberg, Drevlak, NÃijhrenberg, Beidler, Turkin, Loizu and Helander2019; Plunk Reference Plunk2020). Relatively little work has been done to understand QS analytically far from the axis (Sato Reference Sato2022). One reason for this gap is the lack of a standard Grad–Shafranov equation (of the elliptic type), well known for axisymmetric equilibria. In 3-D, even assuming the existence of nested surfaces, the Grad–Shafranov equation is non-standard and has mixed characteristics (Burby, Kallinikos & MacKay Reference Burby, Kallinikos and MacKay2020).
In this work, we consider the limit of large aspect ratio and weak perturbations away from axisymmetry to extract a tractable, semi-analytical model that reproduces many of the qualitative features of numerically optimised QA configurations. The regime explored here is inaccessible by near-axis theory and difficult to analyse using the Variational Moments Equilibrium Code (VMEC, Hirshman & Whitson Reference Hirshman and Whitson1983) or similar tools because of sharp features in the plasma boundary which appear in many models, but in general are difficult to resolve with a Fourier representation (Schuett & Henneberg Reference Schuett and Henneberg2024; Gaur et al. Reference Gaur, Panici, Elder, Landreman, Unalmis, Elmacioglu, Dudt, Conlin and Kolemen2025). We therefore stress the complementary role that analytical theory must play in enhancing, explaining and providing new directions for purely numerical optimisation.
Our approach is stimulated by a recent 3-D quasisymmetric, high-beta stellarator (QS-HBS) model developed by Sengupta et al. (Reference Sengupta, Nikulsin, Gaur and Bhattacharjee2024b
) and extended by Nikulsin et al. (Reference Nikulsin, Sengupta, Jorge and Bhattacharjee2024), based on the HBS model of Freidberg (Freidberg Reference Freidberg2014). Freidberg’s model assumes large aspect ratio (so inverse aspect ratio
$\epsilon \ll 1$
) and ‘high’ plasma beta
$\beta \equiv 2 \mu _0 B^2/p \sim \epsilon$
(with scalar thermal pressure
$p$
and magnetic field
$B$
). The magnetic field, normalised by its value
$B_0$
at the magnetic axis, is toroidal at
$O(1)$
, with an
$O(\epsilon )$
poloidal part generated by the toroidal current. The QS-HBS model of Sengupta et al. (Reference Sengupta, Nikulsin, Gaur and Bhattacharjee2024b
) imposed QA on Freidberg’s HBS model and further required the magnetic axis to be a circle of radius
$R_0 \equiv 1/\epsilon$
plus an
$O(1)$
deviation from axisymmetry (where all distances are in units of the minor radius), and showed that the QA error is
$O(\epsilon ^2)$
if the non-axisymmetric perturbation to the axis is purely vertical. The generalisation by Nikulsin et al. (Reference Nikulsin, Sengupta, Jorge and Bhattacharjee2024) showed that the axis may also have some radial perturbation so long as the axis torsion is
$O(\epsilon ^2)$
, and the QA error then remains
$O(\epsilon ^2)$
with zero flux surface shaping. With all these conditions, the QS-HBS model relaxes the QS overdetermination problem and yields a Grad–Shafranov equation of elliptic type.
In § 2, we use this model to explore an interesting class of global solutions using a modified form of a method pioneered by Palumbo (Reference Palumbo1968) and extended by Bishop & Taylor (Reference Bishop and Taylor1986) and Hernandes & Clemente (Reference Hernandes and Clemente2009) to construct isodynamic axisymmetric solutions, i.e. flux surfaces on which
$|\boldsymbol{B}|$
is constant. In § 3, we show that the resulting quasiaxisymmetric configurations exhibit many of the qualitative features observed in finite-aspect-ratio systems such as those found numerically by Landreman and Paul (Henneberg et al. Reference Henneberg, Drevlak, NÃijhrenberg, Beidler, Turkin, Loizu and Helander2019; Landreman & Paul Reference Landreman and Paul2022). These include very low magnetic shear, D-shaped equilibria for small
$\iota$
and inboard cusp formation as the maximum plasma volume is approached. We stress that these features of our model cannot be understood solely from the existing near-axis framework. Cusp formation in our model appears as the plasma reaches the maximum allowable volume, which differs from surface self-intersection in the near-axis theory, and which simply indicates the breakdown of the near-axis approximation (Landreman Reference Landreman2021). Similarly, the magnetic shear can be calculated from the near-axis theory (Rodríguez et al. Reference Rodríguez, Sengupta and Bhattacharjee2023), but such calculations are only strictly valid in the vicinity of the axis. In § 4, we exploit the decoupled axis and plasma shaping to show that non-negligible rotational transform may be generated from weak axis perturbations, and we use VMEC to verify the asymptotic behaviour of the QA error. We conclude by considering some possible applications of this work to quasiaxisymmetric stellarator design.
2. Analytic solutions of the quasisymmetric Grad–Shafranov equation: an extended Palumbo approach
Throughout this work, we use standard cylindrical coordinates
$(R, \phi , Z)$
and local Cartesian coordinates
$(x, \xi )$
centred on the magnetic axis
$(R_0(\phi ), Z_0(\phi ))$
such that
$x=R - R_0(\phi )$
and
$\xi = Z- Z_0(\phi )$
, and define the perpendicular gradient
$\boldsymbol{\nabla} _{\perp }\equiv (\partial _x \text{, }\partial _{\xi })^T$
.The quasisymmetric Grad–Shafranov equation (QS-GSE) of the QS-HBS model takes the form (in normalised units)
\begin{equation} {\nabla} _{\perp }^2 \varPsi + 2x \frac {{\rm d} \beta }{{\rm d} \varPsi } + \frac {{\rm d} H}{{\rm d} \varPsi } = 0{} \end{equation}
and the toroidal current density
$J_{\phi }$
, in terms of the profiles
$\beta (\varPsi )$
and
$H(\varPsi )$
, is
\begin{equation} \frac {{J_{\phi }}}{B_0} = -\epsilon \left (2x \frac {{\rm d} \beta }{{\rm d} \varPsi } + \frac {{\rm d} H}{{\rm d} \varPsi }\right )\! . \end{equation}
We shall now employ an alternative approach based on an extension (Bishop & Taylor Reference Bishop and Taylor1986; Hernandes & Clemente Reference Hernandes and Clemente2009) of Palumbo’s work on isodynamic equilibrium (Palumbo Reference Palumbo1968) to obtain exact solutions to the QS-GSE (2.1). Our approach uses
$(x,\varPsi )$
as independent variables and solves for
$\xi (x,\varPsi )$
. Therefore, it is an inverse coordinate approach to solve the QS-GSE. The primary advantage of the inverse-coordinate approach stems from the fact that
$\xi$
enters the QS-GSE only through the
${\nabla} _\perp ^2$
term since the pressure term and the current terms are only functions of
$(x,\varPsi )$
. Moreover, this approach shows that the MHD profile functions cannot be arbitrary if we impose QS.
2.1. An inverse-coordinate approach to solve the QS-GSE
Using an ingenious method, Palumbo (Reference Palumbo1968) obtained MHD equilibria with isodynamic magnetic fields, where the magnetic field strength is a constant on a flux surface, i.e.
$|\boldsymbol{B}|=B(\psi )$
. Extensions of Palumbo’s approach were carried out by Bishop & Taylor (Reference Bishop and Taylor1986), Hernandes & Clemente (Reference Hernandes and Clemente2009) and Hernandes & Nogueira (Reference Hernandes and Nogueira2013).
Following the approach outlined by Hernandes & Clemente (Reference Hernandes and Clemente2009), we now carry out the inverse-coordinate analysis starting with the variable transform from
$(x,\xi )$
to
$(x,\varPsi )$
using the identities
\begin{align} {\left .{\partial _x}\right |_{\xi }}= \left .{\partial _x}\right |_{\varPsi } + u\; {\left .{\partial _\varPsi }\right |_{x}} , \quad \left .{\partial _\xi }\right |_{x}= v\; \partial _\varPsi , \quad u\equiv \varPsi _{,x}, \quad v\equiv \varPsi _{,\xi }. \end{align}
Since
$(u,v)$
are components of
$\boldsymbol{\nabla _\perp } \varPsi$
, commutation of partial derivatives
$\partial _\xi u = \partial _x v$
, which can be written in the following equivalent forms (Hernandes & Clemente Reference Hernandes and Clemente2009):
\begin{align} \left ( \partial _x + u\;\partial _\varPsi \right ) v^2 =2 v^2 \partial _\varPsi u \quad \Leftrightarrow \quad \partial _x \left ( \frac {1}{v}\right ) + \partial _\varPsi \left ( \frac {u}{v}\right )=0. \end{align}
We proceed to solve the QS-GSE (2.1), treating
$u$
and
$v$
as the unknowns. It now reads
\begin{align} \partial _x u + \partial _\varPsi \varLambda =0, \end{align}
where we have defined
\begin{align} \varLambda \equiv \frac {1}{2}\eta +2 x\beta (\varPsi ) +H(\varPsi ), \quad \eta (x,\varPsi )\equiv u^2+v^2. \end{align}
Once the QS-GSE (2.5) is solved subject to the consistency condition (2.4), we can obtain
$\xi$
(and thus the flux surface geometry) by numerically integrating
\begin{align} \xi = - \int \frac {u}{v}\; {\rm d}x \quad \text{or equivalently }\quad \xi = \int \frac {1}{v}\; {\rm d}\varPsi . \end{align}
The two integrals are equal to each other, as can be seen from the second form of the consistency condition (2.4).
The system (2.5)–(2.6) is a highly nonlinear system of partial differential equations for
$(u,v)$
. Fortunately, analytical progress can be made based on the following observation. Integrating (2.5) in
$x$
, we get
\begin{align} u(x,\varPsi )=u_0(\varPsi )-x H'(\varPsi )-x^2 \beta '(\varPsi ) -\frac {1}{2}\int {\rm d}x\; \partial _\varPsi \eta (x,\varPsi ) . \end{align}
The form of
$u$
, as given in (2.8), suggests that one could seek a solution where
$|\boldsymbol{\nabla _\perp }\varPsi |^2=\eta (x,\varPsi )$
is a polynomial in
$x$
with coefficients that are functions of
$\varPsi$
, as done by Hernandes & Clemente (Reference Hernandes and Clemente2009). For the cases that Palumbo (Reference Palumbo1968) and Bishop & Taylor (Reference Bishop and Taylor1986) considered, the polynomial nature of
$\eta$
is a direct consequence of the property of
$|\boldsymbol{B}|$
, but in the model considered here, it is an assumption made for tractability.
We shall now seek a solution where
$|\boldsymbol{\nabla _\perp }\varPsi |^2=\eta$
is a quadratic polynomial in
$x$
,
\begin{align} \eta (x,\varPsi ) = \eta _0(\varPsi )+x\; \eta _1(\varPsi ) +x^2 \eta _2{(\varPsi )}. \end{align}
Here,
$\eta _i, i=0,1,2$
are functions of
$\varPsi$
to be determined by solving the nonlinear system (2.5)–(2.4). The quadratic nature of
$\eta$
is suggested by the fact that in (2.8), the non-trivial terms that involve
$\beta ', H'$
are, at most, quadratic in
$x$
. Substituting (2.9) into (2.8), we get
\begin{align} u(x,\varPsi )&= u_0(\varPsi ) + x\: u_1(\varPsi ) +x^2 u_2(\varPsi ),\ \end{align}
where
\begin{align} u_1 &= -\left ( \frac {\eta '_0}{2}+H'\right ), \quad u_2 = -\left ( \frac {\eta '_1}{4}+\beta '\right ) \! . \end{align}
We note that to maintain the quadratic nature of
$u$
, we have chosen
$\eta _2$
to be a constant.
With
$\eta$
and
$u$
available, we can now calculate
$\varLambda , v^2$
. From the definition of
$\varLambda$
(2.6), we find that
$\varLambda$
is a quadratic in
$x$
of the form
\begin{align} \varLambda = \left (\frac {\eta _0}{2} +H(\varPsi )\right ) +2 x\left ( \frac {\eta _1}{4}+\beta (\varPsi ) \right ) +{\frac {\eta _2}{2}}\: x^2. \end{align}
Next, from the forms of
$\eta , u$
as given by (2.9), (2.11), we have
\begin{align} v^2=\eta -u^2, \end{align}
which can be expressed in terms of the following quartic in
$x$
:
\begin{align} v^2&= \sum _{n=0}^{n=4}V_n(\varPsi ) x^n, \end{align}
\begin{align} V_0&=\eta _0-u_0^2, \quad V_1=\eta _1-2u_0 u_1,\quad V_2= {\frac {\eta _2}{2}}-2 u_0 u_2 -u_1^2, \end{align}
\begin{align} V_3 &=-2 u_1 u_2 ,\quad V_4 =-u_2^2. \end{align}
We are now in a position to evaluate the consistency condition (2.4) that would help us determine the functions of
$\varPsi$
such as
$u_0,u_1,\eta _0,\eta _1$
etc. Using
$v^2=\eta -u^2$
, the consistency condition can be written in the following alternative form:
\begin{align} ( \partial _x + u \partial _\varPsi -2u_{,\varPsi }) \eta = \partial _x u^2, \end{align}
which leads to a fourth-order polynomial in
$x$
. Equating each of the coefficients of the polynomials to zero, we get five nonlinear coupled ordinary differential equations (ODEs).
Equating the coefficient of
$x^4$
to zero, we get the condition
\begin{align} u_2=\text{constant}. \end{align}
The nonlinear ODEs can be simplified by using the following variables inspired by Bishop & Taylor (Reference Bishop and Taylor1986):
\begin{align} P=\frac {\eta _0}{u_2^2}, \quad Q=\frac {\eta _1}{u_2^2}, \quad X=\frac {u_1}{u_2}, \quad \varGamma =\frac {u_0}{u_2}, \quad \gamma ={\frac {\eta _2}{2 u_2^2}} ,\quad \beta _p=\frac {\beta }{u_2^2},\quad H_p=\frac {H}{u_2^2}. \nonumber \\[2pt] \end{align}
Furthermore, the equations take a simpler form if we rescale the flux
$\varPsi$
, so we define a new independent variable
$\lambda \equiv \varPsi /u_2$
, so
\begin{align} \frac {\rm d}{{\rm d}\lambda }= u_2 \frac {\rm d}{{\rm d}\varPsi } \equiv \dot {( \;\; )}. \end{align}
In these variables, equating to zero the coefficients of
$(x^0,x,x^2,x^3)$
of the consistency condition (2.15) leads to the following coupled nonlinear ODEs for the system
$(P,Q,X,\varGamma )$
:
\begin{align} Q-2 X \varGamma + \varGamma \dot {P} -2 P \dot {\varGamma }&=0,\\[-10pt]\nonumber \end{align}
\begin{align} 2\gamma -2 X^2 + \varGamma (\dot {Q}-4) -2 Q\dot {\varGamma } + X \dot {P}-2 P \dot {X} &=0,\\[-10pt]\nonumber \end{align}
\begin{align} X(\dot {Q}-6)+ \dot {P}-2(Q \dot {X} + \gamma \dot {\varGamma }) &=0,\\[-10pt]\nonumber \end{align}
\begin{align} \dot {Q}-4-2\gamma \dot {X}&=0 .\end{align}
Admissible initial conditions for the system (2.19) are determined by enforcing regularity of
$\varPsi$
at the magnetic axis and requiring that
$v^2$
have a double root at the axis. The details of the near-axis asymptotics may be found in Appendix B. The crucial result is that all of the functions
$(P, Q, X, \varGamma )$
are completely fixed by choice of two parameters:
$X_0$
and
$\gamma$
. Once the system of ODEs (2.19) is solved numerically, using the MATLAB routine ode23tb, we can find
$H,\beta$
from (2.11) by integrating
\begin{align} \dot {H_p}= -\left ( \frac {1}{2}\dot {P} +X \right )\! , \quad \dot {\beta _p} = - \left ( 1+ \frac {1}{4}\dot Q \right )\! . \end{align}
Though
$u_2$
appears to be a free parameter, if the leading-order magnetic field
$B_0$
is normalised to unity, the poloidal field must be
$O(\epsilon )$
, which gives some restriction on
$u_2$
, viz.
\begin{align} |\boldsymbol{\nabla }\varPsi | &= \sqrt {\eta } = \sqrt {\eta _2 x^2 + \eta _2 x+ \eta _0},\\[-10pt]\nonumber \end{align}
\begin{align} &= |u_2| \sqrt {{2 \gamma } x^2 + Q x + P}, \end{align}
so if the radicand is
$O(1)$
, we are forced to choose
$u_2$
to be at most
$O(\epsilon )$
.
Next, we discuss how we can obtain
$\xi$
using (2.7). The denominator of the integrand is
$v$
, which is a square root of a fourth-order polynomial in
$x$
. Therefore, the integration with
$x$
can be carried out exactly in terms of elliptic integrals. It follows that
\begin{align} \xi = \pm X\int _{s_1(\varPsi )}^s {\rm d}s \frac {\varGamma /X^2+s(s-1)}{\sqrt {\left ( P/X^4-s Q/X^3+s^2 \gamma /X^2 \right )-(s(s-1)+\varGamma /X^2)^2 }}, \quad s\equiv -\frac {x}{X}, \end{align}
where
$s_1(\varPsi )$
is a value of
$s$
such that the denominator vanishes.
It is possible to apply the same method assuming a polynomial ansatz of greater degree for
$u$
. It is shown in Appendix A that this procedure may be carried out without issue for a cubic polynomial, but for polynomials of degree
$4$
or greater, the resulting system of ODEs is overdetermined.
The choice of a quadratic ansatz for
$u$
is post hoc justified by the observation that, to lowest order in
$\epsilon$
,
\begin{equation} {\boldsymbol{b}}\boldsymbol{\cdot }\boldsymbol{\nabla }B = -\epsilon ^2 \frac {\partial _{\xi } \varPsi B_0}{(1+\epsilon x)^2} + O(\epsilon ^3) \end{equation}
and the field strength satisfies
\begin{equation} B = B_0(1-\epsilon x) + O(\epsilon ^2), \end{equation}
so if
$v^2$
is quartic in
$x$
,
$({{\rm d} B}/{{\rm d}\ell })^2$
will be a quartic in
$B$
(where
$\ell$
is the length along the field line). Sengupta et al. (Reference Sengupta, Nikulsin, Buller, Madan, Paul, Nies, Kaptanoglu, Hudson and Bhattacharjee2024a
) found that thousands of QS configurations from the QUASR database (Giuliani Reference Giuliani2024) share the property that, in a quasiaxisymmetric field,
$( {{\rm d}B}/{{\rm d}\ell })^2$
must be either a cubic or a quartic polynomial in
$B$
. The weakly perturbed tokamaks of the QS-HBS model (Nikulsin et al. Reference Nikulsin, Sengupta, Jorge and Bhattacharjee2024; Sengupta et al. Reference Sengupta, Nikulsin, Gaur and Bhattacharjee2024b
) certainly allow deviation from this form for
$ {{\rm d}B}/{{\rm d}\ell }$
, but the Palumbo-like configurations of Bishop & Taylor (Reference Bishop and Taylor1986) and Hernandes & Clemente (Reference Hernandes and Clemente2009) do take the form
\begin{equation} \left ({\textbf{\textit{B}}}\boldsymbol{\cdot }\boldsymbol{\nabla }B^2\right )^2 = \mathbb{P}_{3,4}(B^2){,} \end{equation}
where
$\mathbb{P}_k(t)$
indicates an unspecified polynomial of degree
$k$
in the variable
$t$
.
Assuming a quartic form for
${{\rm d}B}/{{\rm d}\ell }$
, we must have
\begin{equation} v^2 = \eta - u^2 = \mathbb{P}_4(x), \end{equation}
making the simple choice of quadratic polynomials for
$\eta$
and
$u$
especially appealing.
This completes the inverse-coordinate approach to the analytic solution of the QS-GSE. The ODE system (2.19) needs to be solved to determine the profiles
$H$
and
$\beta$
that appear in the QS-GSE using (2.20). Thus, we can clearly see that the two MHD profile functions, and hence the rotational transform profile, cannot be arbitrary if we impose QS and MHS force-balance just like in the isodynamic case (Bishop & Taylor Reference Bishop and Taylor1986; Schief Reference Schief2003; Hernandes & Clemente Reference Hernandes and Clemente2009). Finally, we need to address the issue of the closure of the flux surfaces. Nestedness of flux surfaces is guaranteed by the consistency condition (2.4), but to ensure that flux surfaces are closed, we need to impose additional conditions, which we discuss in Appendix B.
3. Extended Palumbo profiles
We now explore the space of Palumbo-like profiles using the quadratic ansatz for
$\partial _x \varPsi$
described already. We calculate the plasma profiles in terms of the solution of (2.19), using (2.11) and (2.2) to find the pressure and toroidal current (figures 1 and 3, respectively), and integrating
\begin{equation} \iota ^{-1} = -\frac {1}{2\pi }\oint \frac {{\rm d}x}{v} \end{equation}
to find the rotational transform
$\iota$
. The integral over a flux surface in (3.1) is between two vertical asymptotes of the surface, i.e. two roots of the polynomial
$v^2$
. Such real roots are guaranteed to exist near the magnetic axis because we have enforced regularity of
$\varPsi$
in our choice of initial conditions; see Appendix B for more details. Thus, the integral can be efficiently approximated using Chebyshev–Gauss quadrature (Abramowitz & Stegun Reference Abramowitz and Stegun1964) or written explicitly in terms of elliptic integrals (Gradshteyn & Ryzhik Reference Gradshteyn and Ryzhik2007).
The pressure profile is nearly linear in the flux label, independent of the choice of parameters. Here,
$X_0 = 2.5$
,
$\gamma = 0.2$
and
$\epsilon = 1/6$
.
In the limit of small
$\lambda$
(using the near-axis asymptotic behaviour of (2.19) described in Appendix B),
$v$
becomes
\begin{equation} v(s) \simeq X_0^2 \sqrt {-s^2(s-1)^2 + \gamma s^2/X_0^2}, \end{equation}
with
$s = -x/X$
. The radicand has a double root at
$s = 0$
and two single roots
\begin{equation} s_{\pm } = 1\pm \sqrt {\frac {\gamma }{X_0^2}}. \end{equation}
At positive
$\lambda$
, the degeneracy is split, so the double root becomes a pair
$\{\lambda s_1, \lambda s_2\}$
, and the other two roots gain an
$\mathcal{O} (\lambda )$
correction. Substituting the asymptotic relations from Appendix B,
$s_1$
and
$s_2$
are real if
$X_0^2 \gt \gamma$
. We will see that this is equivalent to the requirement that the on-axis rotational transform be real. In the limit of small
$\lambda$
, as we integrate between the roots near zero, these corrections can be neglected:
\begin{align} & \iota ^{-1}|_{\lambda \to 0}\nonumber\\ &= \lim _{\lambda \to 0}\frac {1}{2\pi X}\oint \frac {{ d}s}{\sqrt {-(s-\lambda s_1)(s-\lambda s_2)\left(s-\left(1+\sqrt {\left({\gamma }/{X_0^2}\right)}\right)\right)\left(s-\left(1-\sqrt {\left({\gamma }/{X_0^2}\right)}\right)\right)}} \end{align}
\begin{align} &= \frac {1}{2\pi X_0}\frac {1}{\sqrt {1 - \big({\gamma }/{X_0^2}\big)}}\lim _{\lambda \to 0}\oint \frac {{\rm d}s}{\sqrt {-(s-\lambda s_1)(s - \lambda s_2)}} \end{align}
\begin{align} & = \frac {1}{X_0 \sqrt {1 - \big({\gamma }/{X_0^2}\big)}}, \end{align}
where the integral can be evaluated for any non-zero
$\lambda$
and gives
$\pm 2\pi$
. Thus, the on-axis rotational transform is simply
\begin{equation} \iota = \sqrt {X_0^2 - \gamma }. \end{equation}
It is possible to evaluate the integral (3.1) away from the magnetic axis in terms of elliptic functions, though this representation still depends on the roots of a quartic polynomial whose coefficients are the solutions to a nonlinear set of ODEs. Section 3.145 of Gradshteyn & Ryzhik (Reference Gradshteyn and Ryzhik2007) gives
\begin{equation} \iota ^{-1} = \frac {2 u_2}{\pi X \sqrt (p q)}F\left (\pi , \frac {1}{2} \sqrt {\frac {(p+q)^2+(\alpha - \beta )^2}{p q}} \right )\!, \end{equation}
where, using their notation,
\begin{align} p^2 &= (m-\alpha )^2 + n^2, \nonumber\\[5pt] q^2 &= (m-\beta )^2 + n^2, \end{align}
$\alpha$
and
$\beta$
are two real roots of
$v^2$
,
$m\pm i n$
are the complex conjugate pair of roots of
$v^2$
, and
$F$
is the incomplete elliptic function of the first kind. A similar expression holds when all four roots of
$v^2$
are real.
Further analysis of (3.8) is difficult and we have been unable to prove that the magnetic shear vanishes, even with the given formula in terms of elliptic integrals, so we resort again to numerical methods. Numerical integration using Gauss–Chebyshev quadrature, as in figure 2, shows that the
$\iota$
profile is completely flat. There is zero shear to order
$\epsilon ^2$
in these configurations. Our model is, at least approximately, in the class of ‘isochronous’ dynamical systems, though not obviously of a form exhibited in the existing literature on isochronous systems (Calogero Reference Calogero2008). We emphasise that this characteristic is robust to changes in parameters
$\gamma$
and
$X_0$
, and completely absent from the models of Bishop & Taylor (Reference Bishop and Taylor1986) or Hernandes & Clemente (Reference Hernandes and Clemente2009), both of which show the safety factor diverging at the separatrix. However, it is consistent with the observed low shear seen in optimised QA and the analysis of weakly non-asixymmetric QA of Plunk (Reference Plunk2020), which found that force balance makes the magnetic shear at most second order in the perturbation magnitude.
Integrating (3.1) with Gauss–Chebyshev quadrature shows a completely flat rotational transform throughout the plasma volume.
3.1. Cusps and singularities
We now consider the extreme limits of the equilibria produced by the system (2.19). The singularities of this nonlinear set of ODEs correspond to a singular current profile and a corresponding flux surface cusp.
First, the limit of zero rotational transform requires
$\gamma \to X_0^2$
, which is a singular limit of the ODE initial conditions (see Appendix B). It is apparent from figure 5 that this limit also creates extreme surface elongation near the axis, and D-shaped, positive-triangularity surfaces away from the axis.
At finite rotational transform, the ODEs (2.19) are well behaved near the axis. The profiles eventually develop a singularity at some finite value of
$\lambda$
, just as two roots of
$u$
and
$v$
coincide. The shape of the flux surface is locally determined by the ODE
\begin{equation} \frac {{\rm d} \xi }{{\rm d}x} = \pm \frac {u}{v}, \end{equation}
so for simple roots of
$u$
and
$v^2$
,
$u\sim v^2 \sim (x-x_0)$
implies
\begin{equation} \xi \sim (x-x_0)^{3/2}. \end{equation}
The collision of roots of
$u$
and
$v^2$
also occurs at the magnetic axis, but there we can enforce the condition
$\partial _x v^2 = 0$
to ensure local smoothness of
$\varPsi$
. For a simple root of
$v^2$
, however, the flux surfaces lose analyticity and a current singularity forms, as seen in figure 3.
A current singularity forms at the cusp, which occurs at approximately
$\lambda \simeq 4.77$
in the case
$X_0 = 2.5$
,
$\gamma = 0.2$
and
$\epsilon = 1/6$
. The eventual formation of a current singularity is insensitive to the choice of these parameters.
For most initial conditions, flux surfaces eventually form a cusp. Here,
$X_0 = 2.5$
and
$\gamma = 0.2$
.
figure 4 shows the
$\xi ^2 \sim x^3$
cusp that forms instead of an X-point, in keeping with the proven non-existence of X-points for perfect QS with non-zero rotational transform (Rodríguez & Bhattacharjee Reference Rodríguez and Bhattacharjee2021a
). Because the surface shape is independent of
$\phi$
, the cusp extends to an inboard ridge that moves in tandem with the magnetic axis. Such ridges are common features in optimised QS (Jorge, Giuliani & Loizu Reference Jorge, Giuliani and Loizu2024), but their behaviour is poorly understood. Optimisation efforts aimed at removing sharp features seem to merely push them out further from the magnetic axis, so they appear to be resilient features of QS (Paul et al. Reference Paul, Hyder, Rodríguez, Jorge and Knyazev2025).
The formation of an inboard poloidal field null at large pressure is consistent with experimental results at the poloidal beta limit, as observed in TFTR (Mauel et al. Reference Mauel1992). Conversely, in the limit of zero rotational transform, X-points are allowed (Rodríguez & Bhattacharjee Reference Rodríguez and Bhattacharjee2021a
), which we do observe as flux surfaces push towards a D shape. Though it might appear to be an artefact of our particular choice of ansatz, the index
$3/2$
in (3.11) is much more general. Locally, this index is a hallmark of a cusp catastrophe or a saddle-node or saddle-centre bifurcation (Arnold & Thomas Reference Arnold and Thomas1984; Meiss, Reference Meiss2007), though we emphasise that no bifurcation per se occurs in our model, as there is not sufficient freedom in parameter choice to permit the cusp to bifurcate into multiple critical points. The
$3/2$
cusp is one of only a handful of normal forms for catastrophes. See Appendix B for a comparison of our model to another well-known Hamiltonian system that exhibits cusp formation for some choice of parameters.
A D-shaped flux surface forms when the on-axis rotational transform is small enough. Here,
$X_0 = 1.5$
and
$\gamma = 2.17$
.
Roots of
$u$
and
$v$
as functions of the flux label. At the removable singularity,
$v_3=v_4=u_1$
, but this
$x$
-value is outside the flux surface, which is delimited by
$v_1$
and
$v_2$
. The cusp forms as
$v_2$
approaches
$u_1$
. These two root collisions are indicated by vertical dashed lines.
Maximum QS-breaking Boozer amplitude and root-mean-square norm of QS-breaking modes, fitted against
$\epsilon ^2$
for comparison.
In figure 6, we see that the flux label corresponding to the current singularity and cusp formation coincides with the collision of roots of
$u$
and
$v$
, which occurs as the denominator of (2.19) vanishes. There is an earlier, removable singularity (here, at
$\lambda \simeq 0.22$
) as the numerator also vanishes. With sufficient resolution, the integration can be continued past the first singularity, whereas the latter singularity represents a breakdown of flux surface analyticity and an infinite (albeit integrable) current density.
As observed already, in the large-aspect-ratio limit, QA configurations must satisfy
\begin{equation} \left (\frac {{\rm d}B}{{\rm d}l}\right )^2\sim \left (\frac {\partial \varPsi }{\partial y}\right )^2 \sim \text{quartic or cubic in }x, \end{equation}
so
$v \sim (x-x_0)^{1/2}$
near any simple root
$x_0$
of
$v$
. Therefore, the local form
$y^2 \sim x^3$
is the generic behaviour for a cusp in small-
$\epsilon$
, QA equilbria. See Appendix C for a more detailed discussion of the geometry of flux surfaces near the cusp.
Because we are free to choose any non-zero value of
$\iota$
throughout the plasma volume and thus in particular at the cusp, we emphasise that
$\iota$
is not necessarily rational at the cusp, further distinguishing the cusp from a standard X-point. This family of solutions may perhaps serve as a prototypical class for the study of nonresonant divertors, which have recently been of considerable interest in the stellarator community (Punjabi & Boozer Reference Punjabi and Boozer2020; Garcia et al. Reference Garcia, Bader, Frerichs, Hartwell, Schmitt, Allen and Schmitz2023; Davies et al. Reference Davies, Smiet, Punjabi, Boozer and Henneberg2025).
4. Three-dimensional equilibria
In this section, we explore permissible axis shaping and the associated tradeoffs between torsional rotational transform and QS error. The generalised QS-HBS model decouples flux surface shape, which is independent of the toroidal angle, from the magnetic axis. Using the SIMSOPT package BOOZ_XFORM (Landreman et al. Reference Landreman, Medasani, Wechsung, Giuliani, Jorge and Zhu2021) to calculate the QS error, we see that a purely vertical,
$n=2$
perturbation gives excellent quasisymmetry proportional to
$\epsilon ^2$
, as shown in figure 7 and required by the asymptotic expansion underlying the QS-HBS model. figure 8 illustrates one such vertical axis perturbation.
Plasma boundary with
$m=2$
purely vertical axis perturbation. Colour indicates
$y$
-value for visual clarity. Here, global Cartesian coordinates are used, unrelated to the local coordinates in the rest of this work.
Using the generalised HBS model recently described by Nikulsin et al. (Reference Nikulsin, Sengupta, Jorge and Bhattacharjee2024), it can be shown that the axis perturbations may include some radial component, so long as the perturbation amplitude is
$O(\epsilon )$
and the axis torsion is everywhere
$O(\epsilon ^2)$
. Here, we show that these conditions may be satisfied, generalising a similar calculation of the axis torsion by Sengupta et al. (Reference Sengupta, Nikulsin, Gaur and Bhattacharjee2024b
). The torsional contribution to the total rotational transform (Mercier et al. Reference Mercier and Luc1974) then arises at
$O(\epsilon )$
, the same order as the current-driven part.
Let the axis be defined by
\begin{align} \boldsymbol{r}=R_0 (\boldsymbol{e_R} (1 + \mathfrak{g}(\phi ))+\mathfrak{f}(\phi )\boldsymbol{e_Z}). \end{align}
Then, the curvature is
\begin{equation} \kappa = \left \lvert \frac {{\rm d} \boldsymbol{r}}{{\rm d}\ell } \times \frac {\rm {d}^2 \boldsymbol{r}}{{\rm d}\ell ^2} \right \rvert \end{equation}
and the torsion is
\begin{equation} \tau = \frac {1}{\kappa ^2} \frac {{\rm d} \boldsymbol{r}}{{\rm d}\ell } \times \frac {\rm {d}^2 \boldsymbol{r}}{{\rm d}\ell ^2} \boldsymbol{\cdot }\frac {\rm {d}^3 \boldsymbol{r}}{{\rm d}\ell ^3}. \end{equation}
It is convenient to define the arclength
\begin{align} R_0 s(\phi ) = \frac {{\rm d}\ell }{{\rm d}\phi } = \sqrt {\frac {{\rm d}\boldsymbol{r}}{{\rm d}\phi } \boldsymbol{\cdot }\frac {{\rm d} \boldsymbol{r}}{{\rm d} \phi }} = R_0 \sqrt {(1 + \mathfrak{g})^2 + \mathfrak{g}{\prime }^2 + \mathfrak{f}{\prime }^2}, \end{align}
so that the torsion becomes
\begin{equation} \tau = \frac {1}{\kappa ^2}\frac {1}{R_0^6 s^6}\frac {\partial \boldsymbol{r}}{\partial \phi } \times \frac {\partial ^2 \boldsymbol{r}}{\partial \phi ^2} \boldsymbol{\cdot }\frac {\partial ^3 \boldsymbol{r}}{\partial \phi ^3}, \end{equation}
in which we have excluded those parts of the derivatives of
$\boldsymbol{r}$
that must vanish upon taking the triple product.
Then, the torsional part of the rotational transform is
\begin{equation} \iota _{\tau } = \frac {1}{2\pi }\oint \mathrm{d}\phi s(\phi )\frac {\mathcal{A}}{\mathcal{B}}, \end{equation}
where the integrand is split for legibility into
\begin{equation} \begin{aligned} \mathcal{A} = \mathfrak{f}'''&\left ((1+\mathfrak{g})^2 - \mathfrak{g}''(1+\mathfrak{g}) + 2 \mathfrak{g}'^2\right ) - (1+\mathfrak{g})\left (\mathfrak{f}' \mathfrak{g}'' - \mathfrak{g}' \mathfrak{f}'' - \mathfrak{f}' (1+\mathfrak{g})\right ) \\[5pt] +& \left (\mathfrak{f}'' (1+\mathfrak{g}) - 2 \mathfrak{f}' \mathfrak{g}'\right )(\mathfrak{g}''' - 3 \mathfrak{g}') \end{aligned} \end{equation}
and
\begin{equation} \begin{aligned} \mathcal{B} = & \left ((1+\mathfrak{g})^2 - \mathfrak{g}'' (1+\mathfrak{g}) + 2 \mathfrak{g}'^2\right )^2 + \left (\mathfrak{f}'' (1+\mathfrak{g})-2 \mathfrak{f}' \mathfrak{g}'\right )^2 \\[5pt] &+ \left (\mathfrak{f}' \mathfrak{g}'' - \mathfrak{g}' \mathfrak{f}'' - \mathfrak{f}' (1+\mathfrak{g}) \right )^2. \end{aligned} \end{equation}
In the work that introduced the HBS model, perturbations were restricted to purely vertical displacements, i.e.
$\mathfrak{g} \equiv 0$
, and the integrated axis torsion was
$O(\epsilon ^3)$
(Sengupta et al. Reference Sengupta, Nikulsin, Gaur and Bhattacharjee2024b
). Our more complicated expression (4.6) does reduce to the torsion found previously, viz.,
\begin{equation} \tau _{\mathfrak{g} \to 0} = \frac {1}{R_0}\frac {\mathfrak{f}''' + \mathfrak{f}'}{1 + \mathfrak{f}''^2 + \mathfrak{f}'^2}. \end{equation}
Expanding in powers of
$\epsilon$
and setting
$\mathfrak{f} = \epsilon \mathfrak{f}_1$
, this expression gives
\begin{equation} \iota _{\tau , \mathfrak{g} \to 0} = -\frac {1}{4 \pi }\oint \mathrm{d}\phi \mathfrak{f}_1'^3, \end{equation}
different from the
$\mathfrak{f}_1' \mathfrak{f}_1''^2$
integrand reported in the previous work.
Allowing for non-zero
$\mathfrak{g}$
makes
$\iota _{\tau }$
formally
$O(\epsilon ^2)$
, assuming
$\mathfrak{f}$
and
$\mathfrak{g}$
are both
$O(\epsilon )$
quantities. For judiciously chosen
$\mathfrak{f}$
and
$\mathfrak{g}$
,
$\iota _{\tau }$
can be a significant contribution to the total rotational transform. For instance, taking
$\mathfrak{f}$
and
$\mathfrak{g}$
to be sums of two harmonics, i.e.
\begin{align} \mathfrak{f}(\phi )& = f_{c, m} \cos (m \phi ) + f_{c,n} \cos (n \phi )+f_{s, m} \sin (m \phi ) + f_{s,n} \sin (n \phi ),\nonumber\\[5pt] \mathfrak{g}(\phi ) &= g_{s,m} \sin (m \phi ) + g_{s,n} \sin (n \phi )+g_{c,m} \cos (m \phi ) + g_{c,n} \cos (n \phi ), \end{align}
and choosing a range of aspect ratios
$\epsilon ^{-1} \in [5, 20]$
, we numerically maximise
$\iota _{\tau }$
subject to the (somewhat arbitrary) constraints
$|a_{n,m}, b_{n,m}|\leqslant 3\epsilon$
and
$\max \{\tau (\phi )\}\leqslant 3 \epsilon ^2$
. Keeping only
$n=1$
and
$m=2$
, figure 9 shows that the axis contribution to the rotational transform can be significant, with the caveat that the axis torsion, although formally
$O(\epsilon ^2)$
, will increase the quasisymmetry error.
Numerically optimised rotational transform for two-mode perturbations (
$n=1$
and
$m=2$
) across a range of aspect ratios, subject to the constraint
$\tau _{max}\leqslant 3\epsilon ^2$
. The maximum attainable
$\iota _{\tau }$
is strongly correlated with the maximum allowable torsion.
5. Conclusion
We have constructed a novel class of solutions to the QS-GSE that arises in the limit of large-aspect-ratio, nearly axisymmetric quasiaxisymmetry. These solutions exhibit surprising qualitative similarities to QA equilibria at finite aspect ratio, including negligible shear and the formation of surface cusps as the plasma volume approaches its maximum. The free parameters
$X_0$
and
$\gamma$
in this model primarily allow one to modulate the on-axis rotational transform, which is always constant throughout the plasma volume, and the total volume of the plasma, which is determined by singularity formation in the governing ODEs.
We also note that the condition of weak non-axisymmetry means the Palumbo-QS-HBS model is a perturbed tokamak, so many features found in QS-HBS have well-known analogues in the tokamak literature. The QS-HBS ordering has
$\beta = O(\epsilon )$
and
$B_{\theta }/B_{\phi } = O(\epsilon )$
, so
$\beta _{pol} = O(1/\epsilon )$
and
$\epsilon \beta _{pol}\sim 1$
, a parameter regime that was explored in TFTR, which found that a poloidal field null forms on the inboard side at the poloidal beta limit (Mauel et al. Reference Mauel1992). Furthermore, the slender D-shaped flux surfaces in QS-HBS at
$\beta \sim \epsilon$
are reminiscent of Cowley’s high-
$\beta$
equilibria, which exhibit D-shaped surfaces and extreme Shafranov shift at
$\epsilon \beta _{pol}\gg 1$
(Cowley et al. Reference Cowley, Kaw, Kelly and Kulsrud1991). However, vanishing shear and the local breakdown of analyticity of
$\varPsi$
as a cusp (rather than an X-point) forms are not found in the aforementioned tokamak equilibria, providing suggestive, albeit inconclusive, evidence that these features are inherent to quasiaxisymmetry. In addition to these physical insights into QA, we hope that our Palumbo-like equilibria and other semi-analytical solutions to the QS-GSE will prove to be practically useful as hot starts in QA optimisation.
Acknowledgements
The authors would like to thank S. Buller, E.J. Paul, G. Plunk and M. Zarnstorff for their generous insights and the anonymous referees for their constructive criticism. Editor Per Helander thanks the referees for their advice in evaluating this article
Editor Per Helander thanks the referees for their advice in evaluating this article.
Funding
This work is supported by the DoE grant HiFiStell (until March 31, 2025) and the Simons Foundation Collaboration on Hidden Symmetries and Fusion Energy.
Declaration of interest
The authors report no conflicts of interest.
Appendix A. Higher-degree flux Ansätze
Here, we show that the method of § 2 can be applied to a cubic polynomial ansatz for
$\eta$
, but for any polynomial of degree
$4$
or greater, the resulting system of ODEs is overdetermined.
First, assuming a cubic ansatz
\begin{equation} \eta (x, \varPsi ) = \eta _3 x^3 + \eta _2(\varPsi ) x^2 + \eta _1(\varPsi )x + \eta _0(\varPsi ), \end{equation}
where again the coefficient of the highest-degree term is taken to be constant so that
$u$
is also cubic, we find the following set of ODEs:
\begin{align} 2 u_0 u_1 - \eta _1 - u_0 \eta _0' + 2 u_0' \eta _0 &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} 2 u_1^2 + 4 u_0 u_2 - 2\eta _2 - u_0 \eta _1' - u_1 \eta _0' + 2 u_0' \eta _1 + 2 u_1' \eta _0 &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} 6 u_1 u_2 - \frac {2}{3} \eta _2' u_0 - 3 \eta _3 - \eta _2' u_0 + 2\eta _2 u_0' - \eta _1' u_1 + 2\eta _1 u_1' - \eta _0' u_2 + 2 \eta _0 u_2' &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} 4 u_2^2 - \eta _2' u_1 + \frac {1}{6} \eta _0' \eta _2' + 2 \eta _1 u_2' + 2 \eta _3 u_0' - \eta _2' u_1 + 2\eta _2 u_1' - \eta _1' u_2 &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {1}{6} \eta _1' \eta _2' - \frac {4}{3} \eta _2' u_2 + 2 \eta _3 u_1' - \eta _2' u_2 + 2\eta _2 u_2' &= 0,\\[-10pt]\nonumber \end{align}
\begin{align} \frac {5}{18} (\eta _2')^2 + 2 \eta _3 u_2' &= 0. \end{align}
As
$\eta _3$
and
$u_3$
are constants, we have six equations for six unknowns. However, if
$\eta$
is quartic and we take
$\eta _4$
to be a non-zero constant, the equality
\begin{equation} u = u_0 - x H' - x^2 \beta ' - \frac {1}{2}\int {\rm d}x \left [\eta _0' + x \eta _1' + x^2 \eta _2 ' + x^3 \eta _3'\right ] \end{equation}
gives trivial relations between
$u_4$
and
$\eta _3$
. Once again applying the consistency condition (2.15), we match terms from
$O(x^0)$
up to
$O(x^8)$
, yielding nine equations for the seven unknowns
$u_0$
up to
$u_2$
and
$\eta _0$
up to
$\eta _3$
. By a similar argument, polynomial ansatzes of degree four or greater always generate overdetermined systems.
Appendix B. Details of the extended Palumbo solution
The ODE system (2.19) is highly nonlinear and can have solutions that are not regular or divergent near the axis. Thus, a proper regularity boundary condition must be imposed near the axis. Moreover, as discussed in Appendix A of Bishop & Taylor (Reference Bishop and Taylor1986), the closure condition for the flux surfaces holds in the entire volume if it holds near the axis. Therefore, we must carefully look at the near-axis behaviour of the system (2.19).
The expression for
$\xi$
given in (2.23) suggests the following new variables:
\begin{align} \overline {P}=\frac {P}{X^4}, \quad \overline {Q}=\frac {Q}{X^3}, \quad \overline {\varGamma }=\frac {\varGamma }{X^2}, \quad \overline {\gamma }=\frac {\gamma }{X^2}, \quad \overline {\xi }=\frac {\xi }{X}, \nonumber \\[5pt] \text{such that} \quad \quad \overline {\xi } = \int {\rm d}s \frac {\overline {\varGamma }+s(s-1)}{\sqrt {\left ( \overline {P}-s \overline {Q}+s^2 \overline {\gamma } \right )-\left ( s(s-1)+\overline {\varGamma }\right )^2 }}. \end{align}
Near the magnetic axis
$(x=0,\varPsi =0)$
,
$\eta =|\boldsymbol{\nabla _\perp } \varPsi |$
must vanish. Since
$\eta =u^2+v^2$
, we must have both
$u=v=0$
at
$x=0,\varPsi =0$
. Furthermore, flux surfaces of elliptic or circular cross-sections can exist near the axis provided that
$v^2$
is approximately a quadratic with two real roots, the roots coinciding on the magnetic axis. Therefore,
$\partial _x v^2$
must vanish on the axis. Thus, we get three conditions
\begin{align} u=0, \quad v=0, \quad \partial _x v^2 =0 \quad \text{at} \quad (x=0, \varPsi =0). \end{align}
From (2.9), (2.10) and (2.17), it follows that (B2) implies
\begin{align} P(0)=0, \quad \varGamma (0) =0 ,\quad Q(0)=0. \end{align}
The solution to (2.19) is very well approximated by its near-axis approximations of the form
\begin{align} X &=X_0+\lambda X_1 + \lambda ^2 X_2+O(\lambda ^3), \nonumber \\[5pt] P &= \lambda P_1 + \lambda ^2 P_2+ O(\lambda ^3), \quad \overline {P}= \lambda \overline {P}_1 + \lambda ^2 \overline {P}_2+ O(\lambda ^3), \nonumber \\[5pt] \varGamma &= \lambda \varGamma _1 + \lambda ^2 \varGamma _2+ O(\lambda ^3),\quad \overline {\varGamma }= \lambda \overline {\varGamma }_1 + \lambda ^2 \overline {\varGamma }_2+ O(\lambda ^3), \nonumber \\[5pt] Q &= \lambda Q_1 + \lambda ^2 Q_2+ O(\lambda ^3), \quad \overline {Q}= \lambda \overline {Q}_1 + \lambda ^2 \overline {Q}_2+ O(\lambda ^3), \end{align}
where
$X_i, P_i\ldots$
etc. are coefficients given by
\begin{align} X_1&=\frac {1}{2}\left ( \frac {1}{X_0^2}+\frac {3}{X_0^2-\gamma }\right ) \! , \quad X_2=-\frac {3 \gamma \big(8 X_0^4 - 7 X_0^2 \gamma + 2 \gamma ^2\big)}{16 X_0^5 (X_0^2 - \gamma )^3}, \nonumber \\[5pt] P_1&=2 X_0-\frac {2 \gamma }{X_0}, \quad P_2= \frac {\left (4 X_0^2-3 \gamma \right ) \left (\gamma ^2+10 X_0^4-2 \gamma X_0^2\right )}{4 X_0^4 \left (X_0^2-\gamma \right )^2}, \nonumber \\[5pt] \varGamma _1&=\frac { \left (\gamma +2 X_0^2\right )}{2X_0( X_0^2-\gamma )},\quad \varGamma _2=\frac {2 \gamma ^3-8 X_0^6-3 \gamma ^2 X_0^2}{16 X_0^4 \left (X_0^2-\gamma \right )^3}, \nonumber \\[5pt] Q_1&= \left (\frac {\gamma }{X_0^2}+\frac {3 \gamma }{X_0^2-\gamma }+4\right ) \! ,\quad Q_2=\frac {3 \left (\gamma ^2 \left (2 \gamma ^2+8 X_0^4-7 \gamma X_0^2\right )\right )}{8 \big(X_0^5 \big(X_0^2-\gamma \big)^3\big)}. \end{align}
The convergence of the series in
$\lambda$
clearly requires
$X_0\gt 1, \gamma \ll X_0^2$
. The limit
$\gamma \to X_0^2$
, in which the rotational transform vanishes, is singular. The quantities
$\overline {P},\overline {Q},\overline {\gamma },\overline {\varGamma }$
can be similarly expanded as
\begin{align} \overline {P} &=\frac {2 \lambda \left (X_0^2-\gamma \right )}{X_0^5}+\frac {\lambda ^2 \left (13 \gamma ^3-24 X_0^6+106 \gamma X_0^4-86 \gamma ^2 X_0^2\right )}{4 X_0^8 \left (X_0^2-\gamma \right )^2}+O(\lambda ^3), \nonumber \\[5pt] \overline {Q} &=\frac {\lambda \left( \big((3\gamma)/\big(X_0^2-\gamma\big)\big) + \big(\gamma/X_0^2\big) +4 \right)}{X_0^3}\nonumber \\[5pt]&\quad -\frac {3 \lambda ^2 \left (-2 \gamma ^4 + 64 X_0^8-80 \gamma X_0^6 + 8 \gamma ^2 X_0^4+13 \gamma ^3 X_0^2\right )}{8 \big(X_0^8 \big(X_0^2-\gamma \big)^3\big)} +O(\lambda ^3),\nonumber\end{align}
\begin{align} \overline {\varGamma } &=\frac {\lambda \left (\gamma +2 X_0^2\right )}{X_0^2 \left (2 X_0^3-2 \gamma X_0\right )}-\frac {3 \lambda ^2 \left (2 \gamma ^3+24 X_0^6-16 \gamma X_0^4-7 \gamma ^2 X_0^2\right )}{16 \big(X_0^6 \big(X_0^2-\gamma \big)^3\big)}+O(\lambda ^3), \nonumber \\[5pt] \overline {\gamma } &= \frac {\gamma }{X_0^2}-\frac {\gamma \lambda \left (3/\big(X_0^2-\gamma\big)+ \big(1/X_0^2\big)\right)}{X_0^3}+\frac {3 \gamma \lambda ^2 \!\left (11 \gamma ^2+32 X_0^4-40 \gamma X_0^2\right )}{8 X_0^6 \left (X_0^2-\gamma \right )^3}+O(\lambda ^3). \end{align}
Appendix C. Geometry of the extended Palumbo solution
In this section, we calculate in detail some geometric properties of axisymmetric extended Palumbo equilibria. For clarity, we explicitly include a factor of
$\epsilon$
in the poloidal field component, rather than subsuming this factor into the definition of
$u_2$
.
The magnetic field is, to first order in
$\epsilon$
,
\begin{equation} \boldsymbol{B} = B_0 \boldsymbol{\nabla }\phi + \epsilon \hat {\phi }\times \boldsymbol{\nabla} _\perp \varPsi , \end{equation}
so the unit vector tangent to magnetic field lines is, in cylindrical coordinates,
\begin{equation} \hat {\boldsymbol{b}} = \frac {B_0}{B(1 + \epsilon x)}\hat {\mathbf{\phi }}+\frac {\epsilon }{B} (v \hat {\boldsymbol{R}}- u \hat {\boldsymbol{Z}}) \end{equation}
and the surface normal is
\begin{equation} \hat {\boldsymbol{n}} = \frac {u \hat {\boldsymbol{R}}+v \hat {\boldsymbol{Z}}}{\sqrt {u^2 + v^2}}. \end{equation}
The geodesic curvature may be written as the triple product
\begin{equation} \kappa _g = \frac {{\rm d} \hat {\boldsymbol{b}}}{{\rm d}l}\boldsymbol{\cdot }\big (\hat {\boldsymbol{n}}\times \hat {\boldsymbol{b}}\big ). \end{equation}
The field line curvature
\begin{equation} \frac {{\rm d}\hat {\boldsymbol{b}}}{{\rm d}l} = -\epsilon \hat {\boldsymbol{R}} + O(\epsilon ^2), \end{equation}
so inserting into (C4),
\begin{equation} \kappa _g = -\epsilon \frac {v}{|\boldsymbol{\nabla }\varPsi |}. \end{equation}
A flux surface cusp forms when
$u$
and
$v^2$
share a common, simple root at
$x_0$
. Inserting local behaviour
$u\sim x-x_0$
and
$v\sim (x-x_0)^{1/2}$
, we find
\begin{equation} \begin{aligned} \kappa _g &\sim \frac {(x-x_0)^{1/2}}{\sqrt {(x-x_0)^2 + (x-x_0)}} \sim \frac {(x-x_0)^{1/2}}{(x-x_0)^{1/2}}, \end{aligned} \end{equation}
so the geodesic curvature stays finite even as the flux surfaces form a cusp. The normal curvature is
\begin{equation} \kappa _n = \frac {\boldsymbol{\kappa }\boldsymbol{\cdot }\boldsymbol{\nabla }\varPsi }{|\boldsymbol{\nabla }\varPsi |} = -\epsilon \frac {u}{|\boldsymbol{\nabla }\varPsi |} + O(\epsilon ^2), \end{equation}
which varies as
$\kappa _n \sim (x-x_0)^{1/2}$
near the cusp, vanishing while the geodesic curvature stays finite.
To find the principal curvatures, we calculate the first and second fundamental forms of the outermost flux surface, parametrised in Cartesian coordinates as
\begin{equation} \boldsymbol{r} = (R \cos \phi , R \sin \phi , Z), \end{equation}
where
$Z = Z (R, \varPsi )$
is independent of
$\phi$
in the axisymmetric case. Then, the normal vector
\begin{align} \hat {\boldsymbol{n}}&= \frac {\boldsymbol{r}_{R}\times \boldsymbol{r}_{\phi }}{|\boldsymbol{r}_{R}\times \boldsymbol{r}_{\phi }|}\nonumber\\[5pt] &= \frac {1}{\sqrt {1 + \left (\partial _R Z\right )^2}}\left (-\partial _R Z\cos \phi , -\partial _R Z\sin \phi ,1\right ), \end{align}
the first fundamental form is
\begin{equation} \mathrm{I}_p = E {\rm d}R^2 + G \,{\rm d}\phi ^2, \end{equation}
where
\begin{align} E &= \partial _R \boldsymbol{r}\boldsymbol{\cdot }\partial _R \boldsymbol{r}= \left (1 + \left (\partial _RZ\right )^2\right ),\nonumber\\[5pt] G &= \partial _{\phi } \boldsymbol{r}\boldsymbol{\cdot }\partial _{\phi } \boldsymbol{r}= R^2, \end{align}
and the second fundamental form is
\begin{equation} \mathrm{I\!I}_p = L {\rm d}R^2 + N {\rm d}\phi ^2 ,\end{equation}
where
\begin{align} L&=\boldsymbol{r}_{RR}\boldsymbol{\cdot }\hat {\boldsymbol{n}} =\frac {\partial _R^2Z}{\sqrt {1 + \left (\partial _RZ\right )^2}},\nonumber\\[8pt] N &= \boldsymbol{r}_{\phi \phi }\boldsymbol{\cdot }\hat {\boldsymbol{n}} = \frac {R \partial _R Z}{\sqrt {1 + \left (\partial _RZ\right )^2}}. \end{align}
Since the first and second fundamental forms are both diagonal, we can read off the principal curvatures
\begin{align} \kappa _1 &= \frac {L}{\sqrt {EG}} =\frac {\partial _R^2 Z}{R\left (1 + \left (\partial _RZ\right )^2\right )}\nonumber\\[5pt] &\simeq \epsilon \frac {v ({{\rm d}u}/{{\rm d}x}) - u({{\rm d}v}/{{\rm d}x}) }{u^2 + v^2} \end{align}
and
\begin{align} \kappa _2 &= \frac {N}{\sqrt {EG}} =\frac {R \partial _R Z}{R\left (1 +\left (\partial _RZ\right )^2\right )}\nonumber\\[5pt] &=\frac {uv}{u^2 + v^2} \end{align}
using
$\partial _R Z\to ({{\rm d}\xi }/{{\rm d}x})=\pm \frac {u}{v}$
.
Away from the cusp, a root of
$v^2$
that does not coincide with a root of
$u$
gives the principal curvatures
$\kappa _1\sim -\epsilon (x-x_0)^{-1/2}$
and
$\kappa _2\sim (x-x_0)^{1/2}$
, so the Gaussian curvature
$\mathcal{K}=\kappa _1\kappa _2$
is negative and
$O(\epsilon )$
, but
$\mathcal{K}$
changes sign as the cusp forms since
$\kappa _1\sim \epsilon (x-x_0)^{-1/2}$
when
$u\sim (x-x_0)$
. By continuity, there is thus some value of
$\varPsi$
at which
$\mathcal{K}\to 0$
near the cusp, right as the QS-HBS model breaks down. At the same time, the dominant toroidal field decays as
$1/R$
and the field line is forced to lie on a ridge that has the same
$O(\epsilon )$
deviation from circularity as the magnetic axis, so the magnetic field strength along the ridge only varies at
$O(\epsilon ^2)$
. Therefore, the field becomes locally isodynamic at the ridge. The connection between ridge formation in the last closed flux surface and zeros of the Gaussian curvature is the subject of ongoing research, to be presented in an upcoming publication. That QS should be unattainable when the Gaussian curvature vanishes is reminiscent of the fact, familiar from near-axis theory, that the axis curvature may not vanish in QA.
Appendix D. Analogy to the heavy symmetric top
Here, we quote results from Goldstein’s exposition (Goldstein, Poole & Safko Reference Goldstein, Poole and Safko2002) (§ 5.7) of the heavy symmetric top with one point fixed, and elucidate parallels between that system and the modified Palumbo equilibria constructed in this work. Following Goldstein’s notation, the time elapsed by an orbit
$\mathcal{U}(t)$
takes the form of an elliptic integral
\begin{equation} t = \int _{\mathcal{U}(0)}^{\mathcal{U}(t)}\frac {{\rm d}\mathcal{U}}{\sqrt {(1-\mathcal{U}^2)(\alpha - \beta \mathcal{U})-(b - a\mathcal{U})^2}} \end{equation}
and the integrand is the square root of a cubic polynomial
\begin{equation} f(\mathcal{U}) = \beta \mathcal{U}^3 - (\alpha + a^2) \mathcal{U}^2 + (2 a b - \beta )\mathcal{U} + (\alpha -b^2) \end{equation}
in the variable
$\mathcal{U} = \cos \theta$
, where
$\theta$
is the angle between the vertical and the symmetry axis of the top. The comparison to (3.1) is immediate; in our case,
$v$
is the square root of a quartic in
$x$
instead of a cubic in
$\mathcal{U}$
, but a formal similarity is apparent.
The behaviour of the top may be qualitatively understood from the root of the monomial
\begin{equation} g(\mathcal{U}) \equiv b - a\mathcal{U} \end{equation}
in relation to the roots
$\mathcal{U}_{1,2,3}$
of
$f(\mathcal{U})$
, which are all real. If
$\mathcal{U}_1 \lt \mathcal{U}'\lt \mathcal{U}_2$
, the tip of the top traces loops on the unit sphere; if
$\mathcal{U}' \gt \mathcal{U}_2$
, the top smoothly nutates without looping, as shown in figure 10. Most interestingly, if
$\mathcal{U}' = \mathcal{U}_1$
or
$\mathcal{U}' = \mathcal{U}_2$
, the loops contract and the locus develops cusps. In other words, the top locus develops cusps when
$g(\mathcal{U})$
, a polynomial of degree
$1$
, has a common root with
$f(\mathcal{U})$
, a polynomial of degree
$3$
. In the modified Palumbo case,
$\partial _x \varPsi$
is a polynomial of degree
$2$
and
$(\partial _{\xi }\varPsi )^2$
is a polynomial of degree
$4$
, but just as for the top, a cusp forms when these two polynomials share a root. We also note that the sum
\begin{equation} f(\mathcal{U})+g(\mathcal{U})^2 = (1-\mathcal{U}^2)(\alpha - \beta \mathcal{U}) \end{equation}
is analogous to
$|\boldsymbol{\nabla} _{\perp }\varPsi |^2$
, which is a quadratic in
$x$
in the modified Palumbo case. In both models, when
$g(\mathcal{U})$
or
$u$
go to zero simultaneously with
$f(\mathcal{U})$
or
$v^2$
, there is a qualitative change in the local behaviour of the integrand in the
$t$
(or
$\iota ^{-1}$
) integrals.
A crucial difference is that the higher degree in the Palumbo case allows for a pair of complex roots, which may collide outside the domain of the flux surface, creating a removable singularity instead of a physically realised cusp.
Top: paths traced by a precessing, nutating heavy symmetric top in the cases of (a) smooth (
$u'\gt u_2$
), (b) cusped (
$u'=u_2$
) and (c) looped (
$u_1\lt u'\lt u_2$
) nutation. Bottom: analogous magnetic smooth, cusped and looped magnetic surfaces as the non-axisymmetric perturbation is increased. Reproduced from Plunk (Reference Plunk2020) with permission.
X0=2.5
γ=0.2
ϵ=1/6
λ≃4.77
X0=2.5
γ=0.2
ϵ=1/6
X0=2.5
γ=0.2
X0=1.5
γ=2.17
u
v
v3=v4=u1
x
v1
v2
v2
u1
ϵ2
m=2
y
n=1
m=2
τmax⩽3ϵ2
ιτ
u′>u2
u′=u2
u1
Open access
