3.1. Relation to QP symmetry

We start by examining LQP surfaces. By definition, field lines are geodesics of the surface if the geodesic curvature $K_2$ vanishes. Physically, this is because the ‘acceleration’ term $\partial \boldsymbol{\hat {t}}/\partial l$ has no component in the tangent plane of the surface, spanned by $(\boldsymbol{\hat {t}},\boldsymbol{\hat {b}})$ , as can be seen from (2.3a ). By (2.5a ), this implies that $\eta _l=0$ , and we can therefore change coordinates from $(l,\alpha )$ to $(l',\alpha )$ by shifting the origin of the field lines,

(3.2) \begin{align} l' := l - \int ^\alpha \text{d}\alpha ' \, \eta (\alpha ',\psi ). \end{align}

It is noted that in a toroidal geometry, one has to exercise care in extending the $(l,\alpha )$ coordinates globally, but this transformation introduces no novel issues, since it is merely a statement of where we choose to define $l=0$ (Imbert-Gerard et al. Reference Imbert-Gerard, Paul and Wright2024). For instance, if $2\pi \psi$ is the enclosed toroidal flux then we can restrict $\alpha \in [0,2\pi )$ and pick the lower limit of the integral in (3.2) to be $0$ . Equation (3.2) then implies that $\eta '=0$ , where $\eta '$ is defined as in (2.2) but now with respect to $(l',\alpha ,\psi )$ coordinates. This choice of coordinates is a particular example of geodesic coordinates (see, e.g. Eisenhart Reference Eisenhart2013). The coordinates $l'$ and $\alpha$ are orthogonal on the surface, since $(\partial _\alpha \boldsymbol{r})_{l'} = \rho ' \boldsymbol{\hat {b}}$ is perpendicular to $(\partial _{l'}\boldsymbol{r})_\alpha = \boldsymbol{\hat {t}}$ . Henceforth, we work in such coordinates, so we drop the prime on $l',\eta '$ , etc. Thus, the geodesic condition on LQP surfaces simply becomes $\eta =0$ .

Additionally, by definition, the current normal to the LQP surfaces vanishes. Thus, by (2.8), we conclude that

(3.3) \begin{align} \left (\mu \rho \right )_\alpha = 0 \end{align}

or equivalently that $\partial B/\partial \alpha = 0$ . This symmetry of the field strength is simply quasisymmetry. Indeed, we define the quasisymmetry vector as the covariant basis vector $\boldsymbol{u}=\partial _\alpha \boldsymbol{r} = \rho \boldsymbol{\hat {b}}$ . Then, using the quasisymmetry definition of Rodríguez et al. (Reference Rodríguez, Helander and Bhattacharjee2020), it follows that locally

(3.4a) \begin{align} \boldsymbol{u}\boldsymbol{\cdot } \boldsymbol{\nabla }B &= 0,\\[-6pt]\nonumber \end{align}

(3.4b) \begin{align} \boldsymbol{u}\times \boldsymbol{B} &= \boldsymbol{\nabla }\psi ,\\[-6pt]\nonumber \end{align}

(3.4c) \begin{align} \boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} &= 0. \end{align}

Equation (3.4a ) follows from $B=1/\mu \rho$ and (3.3). Equation (3.4b ) follows from $\boldsymbol{u}=\rho \boldsymbol{\hat {b}}$ , $\boldsymbol{B}=\boldsymbol{\hat {t}}/\mu \rho$ , $\boldsymbol{\hat {b}}\times \boldsymbol{\hat {t}}=\boldsymbol{\hat {n}}$ and $\boldsymbol{\nabla }\psi = \boldsymbol{\hat {n}}/\mu$ . Note that (3.4b ) usually has an additional minus sign, but we define $\boldsymbol{u}$ such that $\rho \gt 0$ and $(\boldsymbol{\hat {t}},\boldsymbol{\hat {n}},\boldsymbol{\hat {b}})$ forms a right-handed coordinate system. Equation (3.4c ) follows from $\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = \partial _\alpha (\mu \rho )/ \mu \rho =0$ . It is notable that we did not have to assume quasisymmetry; it emerged naturally from assuming geodesic field lines and zero current normal to the surfaces and revealed itself in geodesic coordinates.

We now demonstrate why QP surfaces are natural choices for LQP surfaces. We prove that, given a toroidal surface with zero enclosed toroidal current and pointwise zero normal current, the field lines are geodesics on that surface if and only if the magnetic field strength on the surface has QP symmetry.

Firstly, we prove the forward direction: geodesic field lines imply QP symmetry. We have already shown that LQP surfaces are quasisymmetric, so they possess a symmetry vector $\boldsymbol{u}$ and a corresponding symmetry coordinate $\alpha$ . On a toroidal surface where $B$ is not constant on a flux surface, the integral curves of $\boldsymbol{u}$ must close, according to (3.4a ). The integral of $\boldsymbol{B}$ around a closed $\boldsymbol{u}$ integral curve (denoted $\partial S$ ) must vanish because $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{B}=0$ by the choice of geodesic coordinates, so

(3.5) \begin{align} 0=\oint _{\partial \mathcal{S}} \text{d}\alpha \, \frac {\partial \boldsymbol{r}}{\partial \alpha }\boldsymbol{\cdot }\boldsymbol{B}= \iint _{\mathcal{S}} \text{d}\boldsymbol{S} \boldsymbol{\cdot }\boldsymbol{J}, \end{align}

where the second equality follows from Stokes’ theorem. If $\partial \mathcal{S}$ closed toroidally then $\iint _S \text{d}\boldsymbol{S}\boldsymbol{\cdot }\boldsymbol{J}$ is the external poloidal current, which does not vanish. If $\partial \mathcal{S}$ closed poloidally then $\iint _S \text{d}\boldsymbol{S}\boldsymbol{\cdot }\boldsymbol{J}$ is the enclosed toroidal current, which does vanish. So, $\partial S$ as a poloidal curve is a possible solution. As a final possibility, if $\partial \mathcal{S}$ closed helically then $\iint _S \text{d}\boldsymbol{S}\boldsymbol{\cdot }\boldsymbol{J}$ is a linear combination between external poloidal and enclosed toroidal currents, which also does not vanish. So, $\partial \mathcal{S}$ must close poloidally. This yields QP symmetry.

The poloidal current exterior to the flux surface does not vanish, whereas the toroidal current enclosed by the flux surface does vanish.

We now show the backward direction: QP symmetry implies geodesic field lines, given zero enclosed toroidal current and $\boldsymbol{J}\boldsymbol{\cdot }\boldsymbol{\hat {n}}=0$ . By quasisymmetry and the vanishing of the normal current, we have $(\eta /\mu \rho )_l=0$ . For an irrational surface, this means that $\eta /\mu \rho$ is constant on the surface, since field lines get arbitrarily close to every point on the surface. For rational surfaces, assuming that the rotational transform is not constant and we have QP in an infinitesimal volume around the surface, we also have the fact that $\eta /\mu \rho$ is constant by continuity. Then, we integrate the magnetic field along the poloidal symmetry curve $\partial \mathcal{S}$ . Using the fact that the enclosed toroidal current vanishes and Stokes’ theorem, this integral vanishes:

(3.6) \begin{align} 0=\oint _{\partial \mathcal{S}} \text{d}\alpha \, \frac {\partial \boldsymbol{r}}{\partial \alpha }\boldsymbol{\cdot }\boldsymbol{B}= \oint _{\partial \mathcal{S}} \text{d}\alpha \, \frac {\eta }{\mu \rho }. \end{align}

As $\eta /\mu \rho$ is constant, this implies that $\eta =0$ , and consequently, that we have geodesic field lines.

In Appendix C we make the connection with Boozer (Reference Boozer1981) coordinates, which can be defined since $\boldsymbol{J}\boldsymbol{\cdot }\boldsymbol{\hat {n}}=0$ (Rodríguez et al. Reference Rodríguez, Sengupta and Bhattacharjee2021). We also show that, under the further assumption of MHS force balance, the parallel current on a QP surface vanishes whenever the enclosed toroidal current vanishes. Thus, we conclude from this subsection that QP symmetry is the natural global analogue for LQP surfaces. However, LQP surfaces, defined locally, constitute a larger class of possible surfaces, toroidal or not (e.g. magnetic mirrors).

3.2. The surface equations

In this subsection we present a reduced, decoupled model for single LQP flux surfaces. We first discuss what it means to have defined a single LQP surface in the language of the moving frame.

A magnetic field is exactly described by $\eta ,\rho ,\lambda ,\mu ,\nu ,\boldsymbol{K},\boldsymbol{\chi }$ and $\boldsymbol{\sigma }$ as functions of $(l,\alpha ,\psi )$ , provided that the geometric constraints of § 2 are obeyed. In contrast, a single flux surface, taken as an embedded surface in real space, is precisely described by $\eta ,\rho ,\boldsymbol{K}$ and $\boldsymbol{\chi }$ as functions of $(l,\alpha )$ , provided that (2.4a ) and (2.5) are satisfied. This result of classical differential geometry is referred to as Bonnet’s theorem (see, e.g. § 3.6 in Struik Reference Struik1950). A magnetic field requires additional information because it is a 3D bundle of field lines, whereas a single flux surface is a 2D manifold.

Even if we consider only a single LQP flux surface, finding an LQP surface is in principle a 3D problem. This is because defining the field strength and current on a flux surface requires information about field lines away from that surface. For the field strength, given $B= 1/(\mu \rho ) = 1/(\boldsymbol{r}_l\boldsymbol{\cdot }\boldsymbol{r}_\alpha \times \boldsymbol{r}_\psi )$ , we need to know field lines in the infinitesimal neighbourhood of the surface, i.e. we need to know $\boldsymbol{r}(l,\alpha ,\psi _0+\delta \psi )$ to first order in $\delta \psi$ , where $\psi _0$ labels the relevant flux surface.

In general, describing currents on a flux surface requires information about field lines further from the surface, i.e. $\boldsymbol{r}(l,\alpha ,\psi _0+\delta \psi )$ to two orders in $\delta \psi$ . This is manifest in the expression (2.8) for the current as terms like $\mu _\psi$ depend on $\boldsymbol{r}_{\psi \psi }$ (see (2.2b )). However, when the field lines are geodesics on the surface, one requires information only about field lines in the immediate infinitesimal neighbourhood, i.e. to first order in $\delta \psi$ , to define the currents normal to the surface and parallel to the field lines. This is because when $\eta =0$ on the flux surface, which is the geodesic condition in geodesic coordinates as described in § 3.1, there are no terms involving $\boldsymbol{r}_{\psi \psi }$ contributing to the normal and parallel current, as can be seen from (2.8).

Therefore, to define a single LQP surface, we need to define a particular $\boldsymbol{r}(l,\alpha ,\psi _0+\delta \psi )$ to first order in $\delta \psi$ . In other words, a single LQP surface is defined by field lines on and directly adjacent to that surface. Geometrically, defining such field lines is equivalent to providing $\lambda ,\mu ,\nu $ and $\boldsymbol{\sigma }$ on the surface as functions of $(l,\alpha )$ at fixed $\psi =\psi _0$ and $\boldsymbol{K},\boldsymbol{\chi },\rho $ and $\eta$ as functions of $(l,\alpha ,\psi )$ for $\psi$ in the immediate neighbourhood of $\psi _0$ , provided that (2.4), (2.6) and (2.7) hold at $\psi =\psi _0$ and (2.5) holds for $\psi =\psi _0+\delta \psi$ for infinitesimal, arbitrary $\delta \psi$ . Solving these equations for those variables merely provides the requisite structure; an LQP surface must also satisfy the physical constraints of having (i) geodesic field lines, (ii) zero normal current, and (iii) zero parallel current:

(3.7a) \begin{align} \eta &= 0 \quad \text{for } \psi =\psi _0+\delta \psi ,\\[-10pt]\nonumber \end{align}

(3.7b) \begin{align} (\mu \rho )_\alpha &= 0\quad \text{for } \psi =\psi _0,\\[-10pt]\nonumber \end{align}

(3.7c) \begin{align} \lambda _\alpha &= 0 \quad \text{for } \psi =\psi _0. \end{align}

The first equation corresponds to having geodesic field lines and choosing $(l,\alpha )$ to be geodesic coordinates (as discussed in § 3.1). We do not just impose $\eta =0$ on the relevant flux surface, but also in the neighbourhood of the flux surface. This reflects the requirement that a finite volume of surfaces possess the physical benefits of geodesic field lines. The second equation corresponds to zero normal current and is equivalent to quasisymmetry (also discussed in § 3.1). Given the previous two equations, the third equation corresponds to zero parallel current on the relevant flux surface.

We note one crucial simplification arising from field lines being geodesics. By (2.5a ), we deduce that $K_2=0$ . Thus, the Darboux frame on the surface reduces to the Frenet–Serret frame for the field lines, since (2.3a ) becomes the Frenet–Serret formulae (Serret Reference Serret1851; Frenet Reference Frenet1852). To comply with the usual notation, we replace ‘ $\boldsymbol{K}$ ’ with the curvature $\kappa := -K_3$ (equal to the normal curvature) and torsion $\tau := -K_1$ . We note that the curvature defined in this way can be negative, precisely when the principal normal of the field lines is anti-aligned with the surface normal.

We have determined the relevant variables and equations for finding a single LQP surface using the moving frame formalism. In Appendix D we manipulate and reduce these equations to obtain the main result of this paper: the geometry of a single LQP flux surface is determined by the solution to the surface equations

(3.8a) \begin{align} \kappa _\alpha &= -2\rho _l \tau - \rho \tau _l,\\[-10pt]\nonumber \end{align}

(3.8b) \begin{align} \tau _\alpha &= \rho _l \kappa + \left ( \frac {\rho _{ll}-\rho \tau ^2}{\kappa } \right )_l,\\[-14pt]\nonumber \end{align}

(3.8c) \begin{align} \left ( \frac {\kappa }{\rho }\right )_\alpha &= 0. \end{align}

For the reader who skipped to these equations, we reiterate that $\kappa$ is the signed curvature of the field lines, $\tau$ is the torsion of the field lines, $(l,\alpha )$ are a particular choice of field-line-following coordinates where $\partial _\alpha \boldsymbol{r}=\boldsymbol{u}$ is the quasisymmetry vector, and $\rho =|\boldsymbol{u}|$ is the ‘orthogonal distance’ between the field lines.

The surface equations have three remarkable properties. Firstly, the surface equations depend only on the geometry of the flux surface and the field lines on that flux surface (through the variables $(\rho ,\kappa ,\tau )$ ); they do not depend on field lines away from the flux surface. How can this be if the field strength and currents on the flux surface depend on field lines away from the flux surface? This is because of a second remarkable property: given a solution to the surface equations (representative of an embedded surface foliated by field lines), there exist (at least locally) some neighbouring field lines that would confer upon the surface the relevant LQP properties of zero normal and parallel currents. We show this mathematically by providing a set of equations that (along with the surface equations) are equivalent to the original problem for finding an LQP surface, and then showing that these equations admit a solution. These equations are derived alongside the surface equations in Appendix D and are presented here:

(3.9a) \begin{align} \eta &= 0,\\[-10pt]\nonumber \end{align}

(3.9b) \begin{align} K_2 &=0,\\[-10pt]\nonumber \end{align}

(3.9c) \begin{align} \chi _1 &= -\frac {\rho _{ll}-\rho \tau ^2 }{\kappa },\\[-10pt]\nonumber \end{align}

(3.9d) \begin{align} \chi _2 &= \rho _l,\\[-10pt]\nonumber \end{align}

(3.9e) \begin{align} \chi _3 &= \rho \tau ,\\[-10pt]\nonumber \end{align}

(3.9f) \begin{align} \lambda _\alpha &= 0,\\[-10pt]\nonumber \end{align}

(3.9g) \begin{align} \lambda _l &= \kappa \mu ,\\[-10pt]\nonumber \end{align}

(3.9h) \begin{align} \nu _l &= \frac {\rho _l}{\rho }\nu - 2 \tau \mu ,\\[-10pt]\nonumber \end{align}

(3.9i) \begin{align} \mu _\alpha &= - \frac {\rho _\alpha }{\rho }\mu ,\\[-10pt]\nonumber \end{align}

(3.9j) \begin{align} \sigma _1 &= - \frac {\rho _\alpha }{\rho ^2} \mu - \tau \lambda - \frac {\nu }{\rho }\frac {\rho _{ll}-\rho \tau ^2}{\kappa },\\[-10pt]\nonumber \end{align}

(3.9k) \begin{align} \sigma _2 &= \frac {\rho _l}{\rho }\nu -\tau \mu ,\\[-10pt]\nonumber \end{align}

(3.9l) \begin{align} \sigma _3 &= -\mu _l -\kappa \lambda + \tau \nu ,\\[-10pt]\nonumber \end{align}

(3.9m) \begin{align} \rho _\psi &= \nu _\alpha + \rho _{l} \lambda + \frac {\rho _{ll}-\rho \tau ^2}{\kappa } \mu ,\\[-10pt]\nonumber \end{align}

(3.9n) \begin{align} \eta _\psi &= 0,\\[-10pt]\nonumber \end{align}

(3.9o) \begin{align} \kappa _\psi &= \mu _{ll} + (\kappa ^2 + 3 \tau ^2)\mu + \kappa _l \lambda + \frac {\kappa _\alpha }{\rho }\nu ,\\[-10pt]\nonumber \end{align}

(3.9p) \begin{align} \tau _\psi &= \frac {\tau _\alpha }{\rho }\nu + \tau _l \lambda + \frac {\rho _\alpha }{\rho ^2}\mu _l + \left ( \left ( \frac {\rho _\alpha }{\rho ^2}\right )_l - 2 \frac {\tau }{\rho } \frac {\rho _{ll}-\rho \tau ^2}{\kappa } \right )\mu ,\\[-10pt]\nonumber \end{align}

(3.9q) \begin{align} K_{2\psi } &= 0,\\[-10pt]\nonumber \end{align}

(3.9r) \begin{align} \chi _{1\psi } &= -\left (\frac {\rho _{ll}-\rho \tau ^2}{\kappa } \right )_\psi ,\\[-10pt]\nonumber \end{align}

(3.9s) \begin{align} \chi _{2\psi } &= \rho _{l \psi },\\[-10pt]\nonumber \end{align}

(3.9t) \begin{align} \chi _{3 \psi } &= (\rho \tau )_\psi . \end{align}

These equations apply at fixed $\psi =\psi _0$ . Taken together with the surface equations, they enforce the relevant physical and geometric conditions. We now justify why these equations are solvable, given a solution $(\rho ,\kappa ,\tau )$ to the surface equations. First, we consider the variables $(\lambda ,\mu ,\nu )$ . They are determined by the homogeneous, linear differential equations (3.9f )–(3.9i ). By the Frobenius (Reference Frobenius1877) theorem, there exists a solution $(\lambda ,\mu ,\nu )$ to these differential equations. This is because equations (3.9f ) and (3.9g ) are compatible: both equations agree that $\lambda _{\alpha l}=0$ , since $(\kappa \mu )_\alpha =0$ by the surface equation (3.8c ) and the quasisymmetry condition, presented here as (3.9i ). To uniquely specify a solution to these equations, certain boundary conditions must be specified: $\mu$ as a function of $l$ at fixed $\alpha$ (which corresponds to choosing the field strength; see § 3.4), $\nu$ as a function of $\alpha$ at fixed $l$ and $\lambda$ at a particular point. There is less freedom to specify boundary conditions in $\lambda$ because it is controlled by two equations, (3.9f ) and (3.9g ). The remaining functions $(\eta ,K_2,\boldsymbol{\chi },\boldsymbol{\sigma },\rho _\psi ,\eta _\psi ,\tau _\psi ,\boldsymbol{K}_\psi ,\boldsymbol{\chi }_\psi )$ are algebraically and compatibly determined in terms of $\rho ,\kappa ,\tau ,\lambda ,\mu$ and $\nu$ by (3.9a )–(3.9e ) and (3.9j )–(3.9t ).

A solution to these equations, alongside a solution to the surface equations, represents a single LQP flux surface. This is because, as discussed, a single LQP surface is defined by field lines in an infinitesimal volume around that surface, or equivalently, by the functions $\lambda ,\mu ,\nu ,\sigma ,\boldsymbol{K}, \boldsymbol{\chi },\rho ,\eta ,\boldsymbol{K}_\psi , \boldsymbol{\chi }_\psi ,\rho _\psi $ and $\eta _\psi$ at fixed $\psi =\psi _0$ as functions of $(l,\alpha )$ . Thus, the problem of finding an LQP flux surface has been reduced to solving the surface equations, since the remaining part of the problem is solvable given a solution to the surface equations. It is surprising that we do not have an overdetermined system, despite the numerous equations involved. We discuss the physical reasons for this in § 3.4.

A third and final remarkable property of the reduction to the surface equations is that $(\rho ,\kappa ,\tau )$ completely determine the surface and the field lines on it. Intuitively, this is clear: $(\kappa ,\tau )$ determine the field lines individually and $\rho$ tells us the spacing between the field lines. Formally, we can recall Bonnet’s theorem (stated at the start of this subsection) and note that $\boldsymbol{\chi }$ is algebraically determined by $(\rho ,\kappa ,\tau )$ through (3.9c )–(3.9e ). The relevant compatibility conditions (2.4a ) and (2.5) are satisfied identically when $\boldsymbol{\chi }$ is expressed in this way, with the exception of two components of (2.4a ), which are equivalent to the first two surface equations. Thus, the first two surface equations are purely geometric conditions, as viewed by geodesics. Curiously, the physics of vanishing normal and parallel current reduces to the simple equation (3.8c ).

To summarise this subsection, the problem of finding an LQP surface has been reduced from a 3D problem to a 2D problem, since we need only solve the surface equations for the geometry of an embedded surface. An LQP surface is defined by its properties of field strength and currents, which depend on field lines away from the surface. However, once we have a solution to the surface equations, we can locally guarantee that there exist neighbouring field lines that would give that surface the relevant properties of an LQP flux surface, viz. zero parallel and normal currents. We can make this guarantee because the variables that define the neighbouring field lines can be found by solving a well-posed set of equations. In this sense, the surface equations are decoupled and complete. Appendix E makes a connection with the language of classical differential geometry for the reader who is fluent in such a language.

Finally, we comment on MHS force balance (3.1), which depends on the current in the binormal direction (perpendicular to the field lines and tangential to the flux surface). This component of current has not yet been considered and need not vanish, unlike the parallel and normal components of current. By substituting the expression (2.8) for the current into MHS force balance (3.1), we find that

(3.10) \begin{align} \mu _\psi = -\frac {\rho _\psi }{\rho }\mu + p'(\psi ) \rho ^2 \mu ^3 - \rho \mu ^2 \left ( \frac {\lambda }{\mu \rho }\right )_l. \end{align}

Because we look for a solution on a single surface, (3.10) does not lead to overdetermination, since $\mu _\psi$ was unconstrained by the system of equations in (3.9). Physically, this is because the current in the binormal direction is affected by field lines even further from the surface. Formally, the presence of $\mu _\psi$ in the expression for this component of the current implies that one requires $\boldsymbol{r}(l,\alpha ,\psi _0+\delta \psi )$ to two orders in $\delta \psi$ to define $\boldsymbol{J}\boldsymbol{\cdot }\boldsymbol{\hat {b}}$ on the flux surface $\psi =\psi _0$ . For a single LQP flux surface, these ‘far away’ field lines do not affect the parallel or normal currents, so we can choose them to satisfy MHS force balance on the relevant surface: hence, there is no overdetermination by imposing force balance when considering a single LQP flux surface. This is a fortunate and unique scenario: typically in studies of MHS equilibria, imposing MHS force balance on a quasisymmetric equilibrium creates overdetermination (Garren & Boozer Reference Garren and Boozer1991; Sengupta et al. Reference Sengupta, Nikulsin, Buller, Madan, Paul, Nies, Kaptanoglu, Hudson and Bhattacharjee2024). Moreover, (3.10) implies that $(\mu \rho )_{\psi \alpha }=0$ , which is just a statement of quasisymmetry in the neighbourhood of the surface. Therefore, all considerations thus far, including the surface equations, are consistent with a finite-beta equilibrium. However, it must be cautioned that (3.10) does determine the allowable size of $\delta \psi$ for which the present description is accurate.

3.3. A sketch for solving the surface equations through optimisation

In this subsection we sketch how the surface equations could be solved as an optimisation problem for generating LQP flux surfaces. This could be used as part of a larger optimisation for generating a practical QPS field or for quickly and efficiently generating flux surfaces to study trends in QP equilibria. This approach would be intrinsically faster than usual MHD equilibrium solvers (e.g. VMEC, Hirshman & Whitson Reference Hirshman and Whitson1983 or DESC, Dudt & Kolemen Reference Dudt and Kolemen2020) because we do not need to directly solve for MHD equilibrium within a 3D volume. For concreteness, we take the example of a toroidal QP flux surface, but a similar process could be designed for the generation of asymmetric magnetic mirror flux surfaces with the desirable properties of an LQP surface.

Computationally, we work in the coordinate system $(\vartheta ,\varphi )$ , where $\varphi$ is some particular toroidal angle coordinate and $\vartheta$ is some particular poloidal angle coordinate. The flux surface is parameterised by the single-valued functions $R(\vartheta ,\varphi )$ and $Z(\vartheta ,\varphi )$ , where $R$ and $Z$ correspond to the usual cylindrical coordinates with respect to an axis passing through the centre of the torus. In practice, one could work in discrete Fourier modes of $(\vartheta ,\varphi )$ , e.g. $R_{mn}$ , up to some chosen level of precision.

On the toroidal surface, we also parameterise a foliation of closed poloidal loops. In an exact QP solution, these would correspond to lines of constant $l$ (in the geodesic coordinates used throughout this paper). Equivalently, in an exact QP solution, they are the field lines for the field $\boldsymbol{u}=\partial _\alpha \boldsymbol{r}$ or the loops along which the magnetic field strength would be constant (though we do not solve for the magnetic field strength in this 2D problem). We take the poloidal loops to be level sets of some function $\varLambda (\vartheta ,\varphi )$ , which is necessarily periodic in $\vartheta$ . The field lines are taken as perpendicular to these poloidal loops. For specificity, we could take $\varLambda$ to be the distance along a chosen field line from some reference poloidal loop. In the solved problem where the field lines are geodesics, this $\varLambda$ becomes independent of the chosen field line and is equivalent to our choice of arc length coordinate $l$ .

We generate LQP flux surfaces by solving the surface equations as an optimisation problem. The free parameters are the surface, through $R(\vartheta ,\varphi )$ and $Z(\vartheta ,\varphi )$ , and the poloidal loops, through the function $\varLambda (\vartheta ,\varphi )$ . The function $\varLambda (\vartheta ,\varphi )$ is constrained to have level sets that close poloidally. The cost function, $f[\varLambda ,R,Z]$ , has two parts, $f_1[\varLambda ,R,Z]$ and $f_2[\varLambda ,R,Z]$ , where $f_1[\varLambda ,R,Z]=0$ ensures that the field lines are geodesics and if $f_2[\varLambda ,R,Z]=0$ also, the third surface, (3.8c ), is solved. We now consider these two parts in more detail.

Firstly, we enforce that the field lines are geodesics. The tangent vector to the field lines $\boldsymbol{\hat {t}}$ is obtained by calculating the orthogonal to the poloidally closed level sets of $\varLambda$ , a calculation that depends on $R$ and $Z$ . Then, using the tangent vector field $\boldsymbol{\hat {t}}$ , the geodesic curvature can be calculated, in a way that likewise depends on $R$ and $Z$ and is also calculable via well-established methods. The field lines are geodesics if and only if the geodesic curvature is zero everywhere, so we take this part of the cost function to be

(3.11) \begin{align} f_1[\varLambda ,R,Z] = \int _0^{2 \pi } \text{d}\vartheta \int _0^{2 \pi } \text{d}\varphi \,\sqrt {g}_{\text{surf}} \, |\kappa _g|, \end{align}

where $\kappa _g \equiv K_2$ is the geodesic curvature and the surface Jacobian $\sqrt {g}_{\text{surf}}$ ensures the integral is a surface average.

The second part of the cost function is connected with the third surface equation (3.8c ) and is somewhat more subtle. We define the field line label $\tilde {\alpha }$ to be constant along the magnetic field lines and go from $0$ to $2 \pi$ (for specificity) along the closed poloidal loops. We then take $s$ to be the arc length along the closed poloidal loop and define $\tilde {\rho } := \partial s/\partial \tilde {\alpha }$ . For the third surface equation, we must also calculate the normal curvature $\kappa _n \equiv -K_3$ from $\boldsymbol{\hat {t}}$ . When the geodesic curvature is zero, the normal curvature is the same as the curvature $\kappa$ .

We take the second part of the cost function to be

(3.12) \begin{align} f_2[\varLambda ,R,Z] = \int _0^{2 \pi } \text{d}\vartheta \int _0^{2 \pi } \text{d}\varphi \,\sqrt {g}_{\text{surf}}\, \left |\left ( \ln \frac {\kappa _n}{\tilde {\rho }}\right )_{l \tilde {\alpha }}\right |. \end{align}

To justify this expression, we show that precisely when

(3.13) \begin{align} \left ( \ln \frac {\kappa }{\tilde {\rho }}\right )_{l \tilde {\alpha }}=0, \end{align}

there exists a transformation $\alpha =f(\tilde {\alpha })$ such that the third surface equation (3.8c ) is satisfied. To show this, we suppose that (3.13) is true. Then,

(3.14) \begin{align} \left ( \ln \frac {\kappa }{\tilde {\rho }}\right )_{ \tilde {\alpha }}=g(\tilde {\alpha }) \end{align}

for some particular $g$ . We define $\alpha = f(\tilde {\alpha })$ . It follows that $\rho = \tilde {\rho } ({\text{d}f}/{\text{d} \tilde {\alpha }})$ . We choose $f$ to be

(3.15) \begin{align} f(\tilde {\alpha })= \int _0^{\tilde {\alpha }} \text{d}a \, \exp \left ( \int _0^a \text{d}b \, g(b)\right )\!, \end{align}

where $f$ has all the necessary features of a transform between periodic coordinates. First, it is monotonic; second, for $\tilde {\alpha }$ increasing by a period, which is $2\pi$ here for specificity, $f$ increases by a fixed amount. This relies on $\int _0^{2\pi } \text{d} \tilde {\alpha } \; g(\tilde {\alpha })=0$ , which follows from (3.14) and the fact that $\kappa$ and $\tilde {\rho }$ are single valued. Since $f$ solves the differential equation

(3.16) \begin{align} \frac {\text{d}^2 f}{\text{d} \tilde {\alpha }^2} = \frac {\text{d}f}{\text{d} \tilde {\alpha }} g(\tilde {\alpha }), \end{align}

the third surface equation (3.8c ) is satisfied. In the 2D picture, there is a ‘gauge’ freedom $\alpha \rightarrow f(\alpha )$ (corresponding to $\rho \rightarrow \rho /f'(\alpha )$ ) in the orthogonal coordinate $\alpha$ , but the third surface equation (3.8c ) still imposes a non-trivial constraint, through (3.13).

We do not need to verify the first two surface equations, (3.8a ) and (3.8b ). These are identically satisfied when the field lines are geodesics, with $\kappa =\kappa _n$ and $\rho \rightarrow \tilde {\rho }$ . Indeed, these two surface equations are purely geometric constraints and are satisfied by any geodesic foliation of an embedded surface. In this set-up, the geometry is encoded through $R(\vartheta ,\varphi )$ and $Z(\vartheta ,\varphi )$ .

We also note that in practice, solving for $f_2[\varLambda ,R,Z]=0$ makes sense only if the field lines are sufficiently close to geodesics (i.e. $\kappa _n \gg \kappa _g$ ). In a computational optimisation problem, this could be achieved by including the $f_2$ part only when $f_1$ is sufficiently small. Moreover, there could be numerical issues with the logarithm in (3.12). If geodesic curvature is close to zero everywhere (i.e. $f_1$ is sufficiently small), then normal curvature is likely to only be zero on isolated points or curves, lest the surface have flat regions. In any case, the integrand can be transformed to regularise these possible issues, e.g. by using $|( \tanh (\ln ( {\kappa _n}/{\tilde {\rho }})))_{l \tilde {\alpha }}|$ instead.

In practice, once a toroidal surface has been selected by this method, it can be used as the boundary for a complete equilibrium solver. In vacuum, the boundary shape completely determines the magnetic field inside the boundary. However, because of the locality of the assumptions behind LQP surfaces, there is no guarantee that there are nested QP flux surfaces within the boundary of a surface selected by this method. In this case, a second optimisation, using the equilibrium solver, could be used to fine tune the boundary surface to give the desired global properties. This approach has reasons to be better than existing methods for selecting initial conditions for QP symmetry. For instance, a heuristic, near-axis construction (such as that devised by Goodman et al. Reference Goodman, Camacho Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach, Beidler and Helander2023 and applied to QP surfaces in Haque et al. Reference Haque, Paul, Giuliani, Campagna, Madan and Sengupta2026) necessarily contains deviations from exact QP symmetry, since QP symmetry is not possible through the near-axis expansion. The current framework is not limited to the near axis. Nonetheless, it remains as future work to verify whether this approach can indeed lead to improved QP quality.

3.4. What is so special about LQP surfaces?

There are two unique features of LQP surfaces that lead to the drastic simplification of the surface equations. Indeed, such a reduction cannot be found generally. If we were to impose only quasisymmetry as a physical constraint (i.e. only (3.7b )), the moving frame formalism would lead to equations that are highly coupled and overdetermined. In such cases, describing the magnetic field through a moving frame is not particularly useful.

The first special feature is that, for LQP surfaces, the surface geometry is decoupled from the magnetic field strength. In the case of general quasisymmetry with zero normal current, $\boldsymbol{u}\boldsymbol{\cdot }\boldsymbol{B}=\eta /\mu \rho$ is constant on a flux surface (Rodríguez et al. Reference Rodríguez, Helander and Bhattacharjee2020). This couples a quantity that influences the geometry of the flux surface, $\eta$ , with a quantity that is determined by field lines in the neighbourhood of the surface, namely the field strength $1/\mu \rho$ . Local QP surfaces avoid this coupling by having geodesic field lines, which implies that $\eta =0$ everywhere. Notably, given a single LQP surface, it is always possible to construct another LQP surface with the same flux surface shape but a different magnetic field strength profile (so long as $\partial B/\partial \alpha =0)$ . This is because we can take a complete solution to the LQP equations (3.8) and (3.9), and construct another solution through $\mu \rightarrow g(l) \mu$ for arbitrary function $g$ .

There is a second way in which LQP surfaces are special. Without assuming geodesic field lines, the general expression for the parallel current in the moving frame representation is

(3.17) \begin{align} J_\parallel = \frac {1}{\mu \rho } \left ( \left ( \frac {\eta }{\mu \rho } \right )_\psi -\left ( \frac {\lambda }{\mu \rho } \right )_\alpha \right ) + \frac {\eta }{\mu \rho } \left ( \left ( \frac {\lambda }{\mu \rho } \right )_l -\left ( \frac {1}{\mu \rho } \right )_\psi \right ), \end{align}

where we have continued to assume that the normal current vanishes. As an aside, the second set of terms is the Pfirsch–Schlüter current, which vanishes for geodesic field lines ( $\eta =0$ ). For $\eta =0$ , the expression for the parallel current ceases to depend on $\mu _\psi$ . This means that, for the particular case of geodesic field lines, the parallel current is determined only by field lines in the immediate neighbourhood of the surface, not by field lines at second order from the surface. Despite being more restrictive – there are fewer field lines that can contribute to having zero parallel current – this fact leads to mathematical simplification, since the equations are decoupled from field lines that are not in the immediate neighbourhood of the surface. In other words, in the general non-geodesic case, we have to consider a ‘thicker’ volume around the surface in order to evaluate the parallel current, which means there are more variables coupled into the equations (though optimistically, this also represents increased freedom). It is also this fact that means our single LQP surface can always be force balanced without additional overdetermination, since the binormal current $\boldsymbol{J}\boldsymbol{\cdot }\boldsymbol{\hat {b}}$ can be manipulated without affecting the parallel current, by altering field lines further from the surface through $\mu _\psi$ , in accordance with (3.10).

4.1. The magnetic field strength on toroidal Hasimoto surfaces

We have thus far ignored the magnetic field strength profile. This was because, as discussed in § 3.4, the question of field strength decouples from the question of surface geometry. However, with some additional weak assumptions, we can make further statements about the magnetic field strength on Hasimoto surfaces. We first establish a connection between Hasimoto surfaces and an attractive property for practical stellarators, termed the flat mirror property. We then derive an exact governing differential equation that determines the magnetic field strength profile in terms of the surface geometry quantities, $(\rho ,\kappa ,\tau )$ .

First, we show that a Hasimoto surface in vacuum is a flat mirror, as defined and studied in Velasco et al. (Reference Velasco, Calvo, Sánchez and Parra2023). Velasco et al. show that such magnetic fields exhibit small radial transport of energy and good confinement of bulk and fast ions. A flat mirror is defined as a QI magnetic field satisfying

(4.6) \begin{align} \partial _\psi \left ( B_{\text{max}} -B_{00}\right )=0 , \end{align}

where $B_{00}=\iint \text{d}\theta \text{d}\phi \, B/4\pi ^2$ , with $(\theta , \phi$ ) as Boozer angles, defined in Appendix C, and $B_{\text{max}}$ is the maximum of the magnetic field strength on the surface.

In vacuum, we have $(\lambda /\mu \rho )_l = (1/\mu \rho )_\psi$ , as can be seen from the expression for the current (2.8) or equivalently by setting $p'(\psi )=0$ in the expression for MHS force balance (3.10). Using the fact that $B=1/\mu \rho$ is independent of $\alpha$ , we integrate this to find that

(4.7) \begin{align} \lambda = \frac {1}{B(l,\psi )} \int ^l \text{d}l' \partial _\psi B(l',\psi ) + \frac {f(\psi )}{B(l,\psi )}, \end{align}

where $f$ is the integration constant and is independent of $\alpha$ because we have $\lambda _\alpha = 0$ , as in (3.9f ). We now consider the transformation from $(l,\alpha ,\psi )$ coordinates to $(\phi ,\alpha ,\psi )$ , where $\phi$ is the toroidal Boozer angle. We thus have

(4.8) \begin{align} \left ( \frac {\partial }{\partial \psi }\right )_l &= \left ( \frac {\partial }{\partial \psi }\right )_\phi + \left (\frac {\partial \phi }{\partial \psi }\right )_l \left ( \frac {\partial }{\partial \phi }\right )_\psi\\[-10pt]\nonumber \end{align}

(4.9) \begin{align} &= \left ( \frac {\partial }{\partial \psi }\right )_\phi + \left ( \frac {1}{G} \int \text{d}l' \partial _\psi B(l',\psi ) \right ) \left ( \frac {\partial }{\partial \phi }\right )_\psi . \end{align}

The second equality follows from noting that in vacuum, the magnetic field can be written as $\boldsymbol{B}=G\boldsymbol{\nabla }\phi$ , where $G$ is constant. In other words, $G\phi$ is the magnetic potential for a vacuum magnetic field. Rearranging this expression and using $\partial _\phi = (G/B) \partial _l$ gives

(4.10) \begin{align} \left ( \frac {\partial }{\partial \psi }\right )_\phi = \left ( \frac {\partial }{\partial \psi }\right )_l - \lambda \left ( \frac {\partial }{\partial l}\right )_\psi , \end{align}

where we have absorbed the constant $f(\psi )$ into the choice of origin between $l$ and $\phi$ .

Now, we return to the fact that the current in the $\partial _\alpha \boldsymbol{r}$ direction vanishes: $(\lambda /\mu \rho )_l = (1/\mu \rho )_\psi$ . We combine this with (3.9g ) to find that

(4.11) \begin{align} \frac {\kappa }{\rho } &= (\partial _\psi - \lambda \partial _l ) B\\[-10pt]\nonumber \end{align}

(4.12) \begin{align} &= \left ( \frac {\partial B}{\partial \psi }\right )_\phi , \end{align}

where we have used (4.10) in the second equality. For a foliation of Hasimoto surfaces, $\kappa /\rho =\gamma (\psi )$ , so we find that

(4.13) \begin{align} B(\phi ,\psi ) = B_0(\psi ) + B_1(\phi ), \end{align}

where we use the fact that $B$ depends on only $\phi$ and $\psi$ . This arises from (3.7b ) and is equivalent to QP symmetry in the vacuum case considered here. We find that the magnetic field strength is separable into a radially varying part $B_0(\psi )$ and a surface-varying part $B_1(\phi )$ , which is constant across flux surfaces. As an even stronger constraint on the field strength, this satisfies the flat mirror condition (4.6). Technically, for the benefits of a flat mirror to be gained, it is required that $\partial _\psi B_{00}\gt 0$ (Velasco et al. Reference Velasco, Calvo, Sánchez and Parra2023). However, this can often be achieved by adding a finite plasma pressure, due to the plasma behaving diamagnetically. We have indeed concluded that vacuum foliations of Hasimoto surfaces, which are a subset of QP magnetic fields, are flat mirrors. It is important to note that this property does not follow from QP fields generally, but just this particular class of QP solutions. The optimisation problem sketched in § 3.3 can be easily modified to look for Hasimoto surfaces, which then provides a 2D method for investigating flat mirror QP configurations.

In addition to the flat mirror result, we now derive a differential equation for the magnetic field strength in terms of the surface geometry quantities. In principle, this allows us to express quantities of practical interest, such as the mirror ratio, in terms of the surface geometry. We assume that we have a toroidal Hasimoto surface and that $(\rho /\kappa )_{l \psi }=0$ , which is consistent with assuming that $\rho \propto \kappa$ in the neighbourhood of the Hasimoto surface. We do not assume here that we have vacuum. We first note that (3.9h ) can be written as

(4.14) \begin{align} \left (\frac {\nu }{\rho }\right )_l = -2 \frac {\mu \tau }{\rho }. \end{align}

Also, from (3.8c ), (3.9m ), (3.9o ), we have

(4.15) \begin{align} \left ( \frac {\nu }{\rho }\right )_\alpha = ( \partial _\psi - \lambda \partial _l) \ln \frac {\rho }{\kappa } + \left ( \frac {4 \tau ^2}{\kappa } + \kappa - \frac {\rho _{ll}}{\rho \kappa }\right ) \frac {1}{\rho B} + \frac {1}{\kappa } \left ( \frac {1}{\rho B}\right )_{ll}, \end{align}

where we have used $B=1/\mu \rho$ , which follows from the definition of the field-line-following coordinate system described in § 2. We now equate the $\alpha$ derivative of (4.14) with the $l$ derivative of (4.15) – another compatibility condition. This gives

(4.16) \begin{align} - \left ( \frac {2 \tau }{B \rho ^2}\right )_\alpha = \left [\left ( \frac {4 \tau ^2}{\kappa } + \kappa - \frac {\rho _{ll}}{\rho \kappa }\right ) \frac {1}{\rho B} + \frac {1}{\kappa } \left ( \frac {1}{\rho B}\right )_{ll} \right ]_l, \end{align}

where we have used the fact that $(\rho /\kappa )_l=0$ and $(\rho /\kappa )_{l \psi }=0$ , consistent with field lines obeying the vortex filament equation on the surface and in the neighbourhood of the surface. We now integrate this equation around an integral curve of $\boldsymbol{u}=\partial _\alpha \boldsymbol{r}$ , which is closed in toroidal geometry. Thus, as $2 \tau /B\rho ^2$ is single valued, the left-hand side vanishes and we have

(4.17) \begin{align} \left (\frac {1}{B}\right )_{ll}\oint \text{d}\alpha \frac {1}{\rho \kappa } -2 \left (\frac {1}{B}\right )_l \oint \text{d}\alpha \frac {\rho _l}{\kappa \rho ^2} + \frac {1}{B} \oint \text{d}\alpha \frac {1}{\rho \kappa }\left ( 4\tau ^2 + \kappa ^2 - \frac {2\rho _{ll}}{\rho } + \frac {2 \rho _l^2}{\rho ^2} \right ) = C, \end{align}

where $C$ is some constant and we have used $\partial _\alpha B = 0$ . The notable feature of this equation is that it establishes a direct relationship between the magnitude of the magnetic field and the geometry of the field lines, even though the exact nature of this relationship is obscure.

In this subsection we have used $\kappa /\rho = \gamma (\psi )$ not only on the surface but also in its neighbourhood. In doing so, extra consistency conditions arise by combining equations (3.9m ) and (3.9o ). We have not considered the consequences of enforcing this consistency.

5.1. Quality of QP symmetry

We define the QP error on a surface of a configuration as

(5.1) \begin{align} \text{QP error}=\frac {\int _0^{2 \pi } \text{d}\vartheta \int _0^{2 \pi } \text{d}\varphi \sqrt {g}_{\text{surf}} ({1}/{B^3}) \left | \boldsymbol{B}\times \boldsymbol{\nabla }B \boldsymbol{\cdot } \boldsymbol{\nabla }\psi \right |}{\int _0^{2 \pi } \text{d}\vartheta \int _0^{2 \pi } \text{d}\varphi \sqrt {g}_{\text{surf}}}, \end{align}

where $2\pi \psi$ is the toroidal magnetic flux and $\int _0^{2 \pi } \text{d}\vartheta \int _0^{2 \pi } \text{d}\varphi \sqrt {g}_{\text{surf}} (\ldots )$ is an integral over the surface, given that $(\vartheta ,\varphi )$ are the computational poloidal and toroidal angle coordinates (not necessarily straight field-line coordinates) and $\sqrt {g}_{\text{surf}} = | \partial \boldsymbol{r}/\partial \vartheta \times \partial \boldsymbol{r}/\partial \varphi |$ is the surface Jacobian. This QP error is based on the commonly used ‘two-term’ formulation of quasisymmetry (see, e.g. Rodríguez et al. Reference Rodríguez, Paul and Bhattacharjee2022), but there is only one term because we aim for QP configurations that have zero enclosed toroidal current. The QP objective function used during the optimisations was slightly different to this QP error definition and is spelled out in Appendix G.

For the five configurations, we present the quality of QP symmetry as a function of the surface in figure 5. We also make the comparison with the magnetic field design for the proposed QPS experiment (Nelson et al. Reference Nelson2003; Spong et al. Reference Spong, Hirshman, Lyon, Berry and Strickler2005) and the vacuum, precise QA configuration of Landreman & Paul (Reference Landreman and Paul2022). For the latter, the error is appropriately adjusted to be

(5.2) \begin{align} \text{QA error}=\frac {\int _0^{2 \pi } \text{d}\vartheta \int _0^{2 \pi } \text{d}\varphi \sqrt {g}_{\text{surf}} ({1}/{B^3}) \left | \boldsymbol{B}\times \boldsymbol{\nabla }B \boldsymbol{\cdot }\boldsymbol{\nabla }\psi + ({G}/{\iota })\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{\nabla }B\right |}{\int _0^{2 \pi } \text{d}\vartheta \int _0^{2 \pi } \text{d}\varphi \sqrt {g}_{\text{surf}}}, \end{align}

where $2 \pi G$ is the external poloidal current.

In figure 5 we see that the level of QP error is significantly higher than that of the precise QA configuration, despite the DESCconfigurations representing an improvement over QPS. We do not intend for this work to represent the forefront of possible QP solutions (see Haque et al. (Reference Haque, Paul, Giuliani, Campagna, Madan and Sengupta2026) instead), as these solutions were generated for the purpose of corroborating the surface equations. Nonetheless, it seems that optimising for a QP configuration using standard methods is more difficult than for a QA or quasi-helical configuration. For solutions outside the region of validity of the near axis, this difficulty has not been explained. It might be due to a limitation of current optimisation methods, in which case the scheme of § 3.3 could lead to improved QP solutions. It may also be due to fundamental difficulties with QP configurations. Indeed, as we detail in § 5.3, there is a fundamental tension between the third surface equation and toroidicity that can lead to cusps in the surface.

Figure 5.

The QP error (defined in (5.1)) for six different magnetic field configurations as a function of the flux surface and the QA error (defined in (5.2)) for one additional configuration. Configurations 1–3 were produced using DESC and configurations 4 and 5 were produced using SPEC (QP error plotted only on the boundary where a flux surface necessarily exists). Here `QPS’ corresponds to the magnetic field design for the QPS (Nelson et al. Reference Nelson2003; Spong et al. Reference Spong, Hirshman, Lyon, Berry and Strickler2005) and ‘precise QA’ corresponds to the configuration found in Landreman & Paul (Reference Landreman and Paul2022). Configurations 1–3 show an improvement over QPS. Despite QP error being only targeted on the boundary, there is a lower level of QP error towards the axis.

5.2. How the surface equations are tested

The first two surface equations, (3.8a ) and (3.8b ), are purely geometric and are satisfied by any embedded surface foliated by geodesic field lines. Thus, to verify the first two surface equations, it suffices to ensure that the geodesic curvature is negligible on the boundary surface. As it turns out, this is a direct target for the optimiser as the geodesic curvature is proportional to the LQP error on a surface:

(5.3) \begin{align} \frac {1}{B^3} \boldsymbol{B}\times \boldsymbol{\nabla }B \boldsymbol{\cdot }\boldsymbol{\nabla }\psi = \rho \kappa _g \end{align}

(cf. (5.1) and the optimiser objective functions presented in Appendix G). The fact that $\kappa _g \rightarrow 0$ leads to QP was explained in § 3.1: a surface with geodesic field lines, zero net enclosed toroidal current and zero normal current is necessarily QP symmetric. Despite being a direct target of the optimiser, deviations from zero geodesic curvature do exist and are significant, as demonstrated in § 5.3.

The third surface equation, (3.8c ), is assessed directly. For this, we need to interpret $\kappa$ , $\rho$ and $\partial _\alpha$ for our numerical equilibria. We take $\kappa$ to be the normal curvature, $\kappa _n$ . For calculating $\rho$ , we use the transformation to Boozer coordinates

(5.4a) \begin{align} \theta &= \vartheta + \tilde {\lambda } + \iota \omega ,\\[-6pt]\nonumber \end{align}

(5.4b) \begin{align} \phi &= \varphi + \omega , \end{align}

where $\tilde {\lambda }(\vartheta ,\varphi )$ gives the transformation to straight field-line coordinates and $\omega (\vartheta ,\varphi )$ gives the transformation from those coordinates to Boozer coordinates. Calculating $\tilde {\lambda }$ and $\omega$ from numerical equilibria is standard practice when provided with flux surfaces (since ${\boldsymbol{\nabla }}\psi$ can be straightforwardly defined). For the SPEC equilibria, there is no obvious globally defined $ \boldsymbol{\nabla }\psi$ since there are no nested flux surfaces; see Nies et al. (Reference Nies, Paul, Hudson and Bhattacharjee2022) for how $\tilde {\lambda }$ and $\omega$ are calculated in this case. In vacuum, the magnetic field can be written as $\boldsymbol{B}=G\boldsymbol{\nabla }(\omega + \varphi )$ , where $G$ is constant. Given the surface Jacobian $\sqrt {g}_{\text{surf}} = | \partial \boldsymbol{r}/\partial \vartheta \times \partial \boldsymbol{r}/\partial \varphi |$ , we then take $\rho$ to be

(5.5) \begin{align} \rho = \frac {\sqrt {g}_{\text{surf}} B}{G | (1+\omega _\varphi ) (1+ \tilde {\lambda }_\vartheta ) - \omega _\vartheta (\tilde {\lambda }_\varphi - \iota ) |}, \end{align}

where $\iota$ is the rotational transform. This agrees with the theoretical definition of $\rho$ when $\kappa _g=0$ , as shown in Appendix H.

Finally, for exact QP configurations, we have

(5.6) \begin{align} \left ( \frac {\partial }{\partial \theta } \right )_\phi = 0 \quad \iff \quad \left (\frac {\partial }{\partial \alpha }\right )_l=0, \end{align}

so the third surface equation is verified if and only if $\kappa /\rho$ depends only on $\phi$ . In figure 6 we illustrate what the plots of $\kappa _g$ and $\kappa /\rho$ on the boundary should look like for the surface equations to be verified.

Figure 6.

A schematic example of what it looks like for the surface equations to be verified or falsified on a flux surface. For $\phi \lt \pi /2$ , there is agreement with the surface equations, since the geodesic curvature is zero and $\kappa /\rho$ only depends on $\phi$ . For $\phi \gt \pi /2$ , there is disagreement as the geodesic curvature is non-zero and $\kappa /\rho$ does not have vertical contours.

5.3. The appearance of cusps

Figure 7.

We plot quantities for the first DESC configuration: (a) $\kappa /\rho$ and (b) the geodesic curvature $\kappa _g$ , both on the boundary surface as a function of Boozer angles, and (c) toroidal cross-sections of the boundary flux surface. Lengths are measured with respect to the average major radius. The black line in (a) and (b) is a field line. We see agreement with the third surface equation, except the regions with large geodesic curvature (i.e. where QP breaking is significant).

Figure 7. Long description

Panel A: A contour plot shows the geodesic curvature on the boundary surface as a function of Boozer angles. The x-axis represents the Boozer angle phi, and the y-axis represents the Boozer angle theta. The color scale on the right indicates the geodesic curvature values, ranging from -24 to 30. A black line represents a field line. Panel B: A heatmap displays the geodesic curvature on the boundary surface as a function of Boozer angles. The x-axis represents the Boozer angle phi, and the y-axis represents the Boozer angle theta. The color scale on the right indicates the geodesic curvature values, ranging from -2.16 to 2.16. A black line represents a field line. Panel C: A line plot shows toroidal cross-sections of the boundary flux surface. The x-axis represents the major radius R, and the y-axis represents the minor radius z. Different colors represent different values of the Boozer angle phi, as indicated by the legend on the right.

In this subsection we consider the first DESC magnetic field configuration. For this configuration, the vacuum target was essentially achieved. In figure 7(a) we see that $\kappa /\rho$ varies primarily horizontally in agreement with the third surface equation, except for regions of sharp vertical change around $\theta \approx 1$ and $\theta \approx 4.5$ for $\phi \in (0.75,2.5)$ . This does not falsify the third surface equation because these are precisely the regions where the geodesic curvature deviates significantly from zero, as can be seen from figure 7(b). After all, zero geodesic curvature was an assumption behind the derivation of the surface equations.

In fact, we can understand these deviations from the surface equations through the surface equations themselves. The curvature $\kappa$ is negative whenever the field line principal normal is opposed to the surface normal. Based on intuition from tokamaks, the inboard and outboard sides of the torus typically have opposite signs of $\kappa$ . However, the third surface equation tells us that $\kappa /\rho$ should be constant on poloidal loops, so $\kappa$ should not be able to change sign, since $\rho \gt 0$ . The optimiser rectifies this difficulty by isolating the regions where $\kappa$ changes sign to thin cusps, with $\kappa /\rho$ independent of $\theta$ only away from these cusps. The cusps are visible in figure 7(c) at $\phi \approx 0.5 \pi$ as the tips on the top and bottom of the cross-section. These cusps are also visible in figure 4 where the surface ‘curves around’.

We note that the cusps also appear to follow the field lines. This can be understood through MHS force balance. As the field is vacuum, it is necessarily force free, so the tension force from the field lines can only be balanced by magnetic pressure. The magnetic pressure can only grow so large given the limited field strength variation, so the magnetic tension force has to be minimised by not having field lines go across the cusp.

5.4. The effect of parallel current

We now consider the second DESC magnetic field configuration. In figure 8(a) we observe noticeable departure from the third surface equation, since the contours are largely not close to vertical. There are two related reasons for this.

Figure 8.

For the second DESC configuration: (a) $\kappa /\rho$ , (b) the geodesic curvature $\kappa _g$ , and (c) the parallel current $J_\parallel$ . All quantities are plotted on the boundary surface as a function of Boozer angles. Lengths are measured with respect to the average major radius and the current is normalised by the mean magnetic field on the boundary. There are two related reasons for deviations from the third surface equation: areas of high geodesic curvature and large parallel current.

Firstly, consider the region $\phi \in (0.4,0.8)$ . This region of the surface is similar to the configuration considered in the last subsection: there are approximately straight contours separated by ‘cusps’, which is where the geodesic curvature is largest (see figure 8 b). The geodesic curvature deviates somewhat from zero even away from these cusps, which explains why the vertical contours are slightly curved and angled. As in the last subsection, we can explain the cusps by the fact that $\kappa$ typically has to change sign from the inboard to the outboard side, in contradiction with the third surface equation.

Now, consider the remaining portion of the surface: $\phi \in (0,0.4)$ and $\phi \in (0.8,2\pi /5)$ . The contours are not vertical here at all, despite the geodesic curvature being small. Here, it is another assumption of the surface equations that breaks down: the parallel current is non-zero. Zero current density was a target for the optimiser but was not sufficiently attained everywhere in this configuration.

Even without considering the vacuum target, having non-zero parallel current is somewhat unexpected. This is because of the theoretical result presented at the end of Appendix C: assuming QP symmetry, force-balance and zero net toroidal current, the parallel current should vanish pointwise. The magnetic field configuration presented in figure 8 does indeed have zero net toroidal current and is force balanced (as current perpendicular to the field lines is indeed negligible). Thus, it is deviations from the QP configuration (i.e. deviations from zero geodesic curvature) the leads to non-zero parallel currents. Curiously, the region where geodesic curvature is significant ( $\phi \in (0.4,0.8)$ ) is different to the region where parallel currents are significant: hence, this is a non-local effect.

5.5. When are field lines sufficiently geodesic?

In this subsection we use the third DESC magnetic field configuration to point out that the third surface equation is satisfied (locally) only if the field lines are sufficiently geodesic, in the sense that the geodesic curvature should be much smaller than the normal curvature.

Figure 9.

For the boundary surface of the third DESC configuration: (a) $\kappa /\rho$ , (b) the geodesic curvature $\kappa _g$ , and (c) the logarithm of the geodesic curvature divided by normal curvature $\kappa _n$ . For the third surface equation to be true, the geodesic curvature should be small compared with the normal curvature, which is not the case within the overlaid black/white rectangles.

Figure 9. Long description

Panel A: A contour plot displays the ratio of geodesic curvature to normal curvature. The x-axis represents the variable phi ranging from 0 to 3.5, and the y-axis represents the variable theta ranging from 0 to 6. The color scale on the right indicates values from -28.8 to 36.0. Several black rectangles are overlaid on the plot, highlighting specific regions. Panel B: A heatmap shows the geodesic curvature. The x-axis represents the variable phi ranging from 0 to 3.5, and the y-axis represents the variable theta ranging from 0 to 6. The color scale on the right indicates values from -2.0 to 2.0. Several black rectangles are overlaid on the heatmap, highlighting specific regions. Panel C: A contour plot displays the logarithm of the geodesic curvature divided by normal curvature. The x-axis represents the variable phi ranging from 0 to 3.5, and the y-axis represents the variable theta ranging from 0 to 6. The color scale on the right indicates values from -7.2 to 9.0. Several black rectangles are overlaid on the plot, highlighting specific regions.

We consider two different parts of the boundary surface of this configuration. The region $\phi \in (1.3,\pi )$ is similar to the central region of the first configuration considered in § 5.3: the third surface equation is violated on field-line-following cusps (see figure 9 a), which is where the geodesic curvature (and thus QP error) is large, as per figure 9(b). Away from these cusps, the contours of $\kappa /\rho$ are nearly vertical.

We now consider the region $\phi \in (0,1.3)$ , where the magnetic field strength is largest. Apart from distinct blobs at $(\phi ,\theta )\approx (0.7,3)$ and at $(\phi ,\theta )\approx (0.7,5.5)$ , the geodesic curvature is relatively small in an absolute sense, as per figure 9(b). These blobs can also be thought of as cusps, and also represent areas where the third surface equation does not hold (contours of $\kappa /\rho$ are not vertical). However, even away from these blobs, the contours are not close to vertical, particularly in the regions outlined by the black rectangles in figure 9(a). As seen in figure 9(b), the geodesic curvature is small in these regions, in an absolute sense. However, as seen from figure 9(c), the geodesic curvature is much larger than the normal curvature in these regions, so the field lines cannot be considered geodesics in these regions. Given that the geodesicity of field lines was an assumption of the surface equations, deviation from the third surface equation is to be expected in these regions.

Given that the normal curvature does not directly enter into the definition of QP error (cf. (5.1) and (5.3)), this observation suggests that it is theoretically possible to have a surface of high (but imperfect) apparent QP quality without having that surface conform to the surface equations, as long as the normal curvature is sufficiently small. This would be most readily achieved for a high aspect ratio stellarator, where it is conceivable for normal curvature to remain small everywhere. However, this would constitute a Pyrrhic victory. As discussed in § 1, many practical advantages stem from the radial drift of the guiding centres being much smaller than the poloidal drift. If geodesic curvature is larger than normal curvature, the radial drift is larger than the poloidal drift, so one loses these practical advantages.

Figure 10.

For the first SPEC configuration, we plot (a) the boundary surface in real space, (b) $\kappa _g$ on the boundary, and (c) $\kappa /\rho$ on the boundary. Lengths are measured in terms of the average major radius and the black line in panels (b,c) is a field line. We observe a very large mirror ratio, narrow pinch points and double wells.

5.6. High mirror ratios and narrow pinch points

In this subsection we focus on the two SPEC equilibria. By construction, these are vacuum equilibria. Thus, we need not worry about non-zero currents undermining the surface equations, as was the case in § 5.4. This comes at the cost of not necessarily having nested flux surfaces, which permits pathologically narrow pinch points in the boundary surface, seen in figures 10(a) and 11(a).

In figures 10 and 11 we plot $\kappa /\rho$ and $\kappa _g$ on the boundary surface. Due to the regions of large geodesic curvature, there are deviations from the third surface equation as the contours of $\kappa /\rho$ are not vertical. There are vertical contours in select regions, particularly in figure 10(c).

For completeness, we note an additional potential reason for deviation from the third surface equation. In deriving the surface equations, we assumed not only that the field lines are geodesics on the flux surface, but also that field lines are geodesics in the neighbourhood of the surface (see (3.7a ) and surrounding discussion). In other words, we assumed QP symmetry in the neighbourhood of the surface. Thus, in these SPEC configurations, non-geodesic field lines adjacent to the flux surface can contribute to deviations of the third surface equation. This differs from the DESC configurations, where imposing nested flux surfaces naturally leads to a high level of QP symmetry on neighbouring surfaces.

One can observe in figures 10(a) and 11(a) that the SPEC configurations resemble a system of ‘linked magnetic mirrors’. This observation is consistent with the surface equations, since magnetic mirrors are solutions to the surface equations with poloidally closing orthogonals to the field lines.

We also observe narrow pinch points and extremely high mirror ratios in the SPEC configurations, as seen in figures 10(a) and 11(a). Even with a prescribed number of field periods, the optimiser tended to ‘double’ the field periods in these configurations. In other words, the optimisation would occasionally produce a ‘double well’, with two minima of the magnetic field strength within a single field period.

Based on the similarity with magnetic mirrors, we can explain the appearance of narrow pinch points, high mirror ratios and double wells through the surface equations. As it turns out, magnetic mirrors with narrow pinch points can more readily bend around (as is necessary to form a toroidal surface) and generate a non-zero rotational transform. This claim is justified in Appendix I, where we perturb surfaces of revolution through the surface equations, form an eigenvalue problem (given periodic boundary conditions) and show that the density of eigenvalues increases given a narrow pinch point.

Figure 11.

For the second SPEC configuration, we plot (a) the boundary surface in real space, (b) $\kappa _g$ on the boundary, and (c) $\kappa /\rho$ on the boundary. Lengths are measured in terms of the average major radius and the black line in panels (b, c) is a field line. We observe a very large mirror ratio, narrow pinch points and double wells.

Figure 11. Long description

Panel A: A 3D rendering of a boundary surface in real space. The surface is toroidal with a complex, twisted shape. The color gradient represents the magnetic field strength, with a scale bar on the right indicating values ranging from 10^1 to 10^2. Panel B: A contour plot of the parameter kappa_g on the boundary surface. The x-axis is labeled phi and ranges from 0.0 to 1.0. The y-axis is labeled theta and ranges from 0 to 6. The color scale on the right ranges from -4.32 to 4.32. A black line represents a field line. Panel C: A contour plot of the parameter kappa/rho on the boundary surface. The x-axis is labeled phi and ranges from 0.0 to 1.0. The y-axis is labeled theta and ranges from 0 to 6. The color scale on the right ranges from -204.0 to 25.5. A black line represents a field line.

In addition to this mathematical argument, we can understand physically why narrow pinch points assist in having a non-zero rotational transform. The field lines are geodesics, which are paths of inertial motion on the surface. Akin to a ballerina pulling their arms in to spin faster, a narrower pinch point will assist the geodesic field lines ‘spinning around’ the surface, generating a rotational transform. To consider an explicit example, consider a surface of revolution $\boldsymbol{r}(\phi ,x)=(x,y(x)\text{cos}(\phi ),y(x)\text{sin}(\phi ))$ . The geodesics are given by the solutions to

(5.7) \begin{align} 1=(1+y'(x))^2\dot {x}^2+y(x)^2\dot {\phi }^2, \quad J=y^2\dot {\phi }, \end{align}

for some constant $J$ , which is the angular momentum here. One can see that for smaller $y$ at constant $J$ , $\dot {\phi }$ has to be larger. This physical interpretation does not extend to explaining why narrower pinch points also assist in the surface being able to ‘bend around’.

Finally, there are other possible explanations for the appearance of large mirror ratios and narrow pinch points in optimised QPS. In the near-axis approximation, the curvature of the magnetic axis is proportional to the level of QP breaking at first order from the axis, as per (A4). The optimiser, in its pursuit to satisfy the QP metric, may swamp this necessary poloidal variation of the magnetic field with a large toroidal variation, leading to high mirror ratios. In addition, the near-axis expansion also indicates a preference for highly elongated surfaces, which also corresponds to high mirror ratios and narrow pinch points. High mirror ratios are physically undesirable: for instance, Haque et al. (Reference Haque, Paul, Giuliani, Campagna, Madan and Sengupta2026) find that higher mirror ratios lead to increased fast ion losses, even with an apparently improved QP quality. Furthermore, Dommaschk (Reference Dommaschk1994) provides a theoretical lower bound for the mirror ratio in QPS. Thus, it is practically important to assess whether it is possible to achieve a high quality of QP symmetry without high mirror ratios. Appendix I illustrates that it is consistent with the surface equations to have solutions with high mirror ratios when QP surfaces are close to surfaces of revolution, but it does not rule out the possibility of low-mirror-ratio QP solutions, which may be possible when the near-axis expansion breaks down.

To summarise this section, we have found that the surface equations are approximately obeyed where the LQP error, local parallel currents and geodesic curvature are sufficiently small. It is important to note that these numerical solutions deviate significantly from ideal vacuum QP surfaces. These deviations often occur in local patches. These local patches serve as a reference case for where the surface equations need not be true and indeed are violated. The LQP error is proportional to the geodesic curvature, and we found that for the surface equations to be approximately satisfied, the geodesic curvature must be much smaller than the normal curvature. This is a physically motivated condition, since the ratio between geodesic curvature and normal curvature corresponds to the ratio between radial curvature drift and poloidal curvature drift. We observe a tendency for cusps to form on the surface, corresponding to areas of high geodesic curvature and deviations from the third surface equation. These cusps appear because of the optimiser’s predicament: the third surface equation implies that $\kappa$ has a constant sign along the poloidal loops, but typical toroidal surfaces have a $\kappa$ switching sign between the inboard and outboard sides. This represents a fundamental difficulty in achieving exact QP configurations that apply outside of the regime of validity of the near-axis framework. Finally, for QP surfaces that resemble linked magnetic mirrors, we found that narrow pinch points and high mirror ratios were consistent with the surface equations. It remains as important future work to verify whether high mirror ratios and narrow pinch points are essential to QP surfaces.