2. Vacuum fields with nested surfaces

We consider a generic toroidal coordinate system $(\psi ,\theta ,\phi )$ where $\psi$ is the toroidal flux and $\theta$ and $\phi$ are the poloidal and toroidal angles that are 2${\rm \pi}$ periodic. Denoting the magnetic scalar potential and the field line label by $\varPhi$ and $\alpha$, respectively, we represent the vacuum magnetic field with surfaces as

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

The three components of the above equation yield

(2.2a)\begin{gather} \sqrt{g}\begin{pmatrix} \varPhi_{,\theta}\\ \varPhi_{,\phi} \end{pmatrix}= \begin{pmatrix} g_{\theta\theta} & g_{\theta\phi}\\ g_{\theta\phi} & g_{\phi\phi} \end{pmatrix} \begin{pmatrix} -\alpha_{,\phi} \\ \alpha_{,\theta} \end{pmatrix}, \end{gather}

(2.2b)\begin{gather}\sqrt{g}\varPhi_{,\psi} = g_{\psi\phi} \alpha_{,\theta} - g_{\psi \theta} \alpha_{,\phi}, \end{gather}

where $\sqrt {g} = (\boldsymbol {\nabla } \psi \times \boldsymbol {\nabla } \theta \cdot \boldsymbol {\nabla } \phi )^{-1}$ is the Jacobian and $g_{ij} = \boldsymbol {r}_{,x_i} \cdot \boldsymbol {r}_{,x_j}$ is the metric tensor. The first two equations are the ‘in-surface’ derivatives of $\varPhi$, whereas the third equation gives the ‘out-of-surface’ derivative of $\varPhi$. The magnetic differential operator is given by

(2.3)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} = \frac{1}{\sqrt{g}}\left( \alpha_{,\theta}\partial_\phi - \alpha_{,\phi}\partial_\theta\right)= \frac{1}{\sqrt{g}} \{\alpha,\:\:\}_{(\theta,\phi)}, \end{equation}

where $\{\,f,g\}_{(\theta ,\phi )}$ denotes a Poisson bracket of $f$ and $g$ with respect to $(\theta ,\phi )$. The field strength is computed from

(2.4)\begin{equation} B^2 = \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} \varPhi= \frac{1}{\sqrt{g}}\left(\varPhi_{,\phi} \alpha_{,\theta} - \varPhi_{,\theta} \alpha_{,\phi}\right). \end{equation}

3. The QS constraint on vacuum fields

In the literature various definitions of QS can be found. (A clear comparison is provided in Landreman (Reference Landreman2019b).) Here we discuss the following form of QS when the rotational transform is irrational, namely

(3.1)\begin{equation} \text{QS}: \quad \boldsymbol{B}\times\boldsymbol{\nabla} \psi\cdot \boldsymbol{\nabla} B-F(\psi) \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} B =0. \end{equation}

This form of QS has been extensively discussed previously in the literature (Helander & Simakov Reference Helander and Simakov2008; Helander Reference Helander2014). This version of QS has been recently termed ‘weak QS’ (Burby et al. Reference Burby, Kallinikos and MacKay2020; Rodriguez, Helander & Bhattacharjee Reference Rodriguez, Helander and Bhattacharjee2020; Constantin et al. Reference Constantin, Drivas and Ginsberg2021). However, we note that under the assumption of MHD equilibrium, the strong and weak forms of QS coincide. As we will make the assumption of a vacuum magnetic field, a subset of MHD equilibrium models, our QS definition can be considered ‘strong QS’. We note that $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } B$ vanishes whenever $\boldsymbol {B} \times \boldsymbol {\nabla } \psi \cdot \boldsymbol {\nabla } B$ does (Landreman & Catto Reference Landreman and Catto2012), and their ratio remains finite. In the above expression, $F(\psi )$ is a non-vanishing, smooth function of flux. Under the vacuum field constraint, an explicit expression for $F(\psi )$ is provided in (3.7).

In $(\varPhi ,\psi ,\alpha )$ coordinates with $B^2=\boldsymbol {\nabla }\varPhi \cdot \boldsymbol {\nabla } \psi \times \boldsymbol {\nabla }\alpha$,

(3.2a,b)\begin{equation} \boldsymbol{B}\times\boldsymbol{\nabla}\psi\cdot \boldsymbol{\nabla}=\boldsymbol{\nabla} \varPhi\times \boldsymbol{\nabla} \psi\cdot \boldsymbol{\nabla}=B^2 \partial_\alpha,\quad\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} = \boldsymbol{\nabla} \psi\times \boldsymbol{\nabla} \alpha\cdot \boldsymbol{\nabla}=B^2 \partial_\varPhi. \end{equation}

We assume that $B^2 \ne 0$ in the neighborhood of the surface. Since $B^2\neq 0$, QS implies

(3.3)\begin{equation} B=B(\varPhi +F(\psi)\alpha,\psi). \end{equation}

We define

(3.4)\begin{equation} \eta= \varPhi + F(\psi) \alpha, \end{equation}

such that $B=B(\eta ,\psi )$. Note that a similar function can be constructed even if the rotational transform is rational by replacing $\psi$ by $\oint \textrm {d} l/B$, where the integral is carried out along a closed field line (Elbarmi, Sengupta & Weitzner Reference Elbarmi, Sengupta and Weitzner2020). The streamlines of $\eta$ must close on a constant $\psi$ surface. This implies that

(3.5)\begin{equation} \eta \propto M \theta - N \phi + \tilde{\eta}(\theta,\phi) \end{equation}

where $M$ and $N$ are integers that specify the helicity of the QS. We note that $\varPhi$ and $\alpha$ also generally have secular terms:

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

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

where tilde denotes periodic functions. This implies that $\bar {I} + F(\psi ) \propto M$ and $-\bar {G} + F(\psi )\iota \propto N$:

(3.7)\begin{equation} F(\psi) = \frac{\bar{G} + N/M \bar{I}}{\iota(\psi) -N/M}.\end{equation}

While $\bar {I}$ vanishes in a global vacuum field, the surface about which we expand may enclose currents that are not in the neighbourhood of the surface. As we are performing a local expansion, this does not preclude the existence of currents asymptotically away from the surface. Thus we retain both secular terms in the scalar potential. In terms of the variable $\eta$ one can write $B^2$ using (2.4) and (3.4) as

(3.8)\begin{equation} B^2= \frac{1}{\sqrt{g}}\{\alpha,\varPhi\}_{\theta,\phi}=\frac{1}{\sqrt{g}}\{\alpha,\eta\}_{\theta,\phi}. \end{equation}

The ‘in-surface’ derivatives of $\varPhi$, $\alpha$ and $\eta$ can then be related to each other using (2.2a) and (3.4) as

(3.9)\begin{equation} \sqrt{g}\begin{pmatrix} \varPhi_{,\theta}\\ \varPhi_{,\phi} \end{pmatrix}= \sqrt{g}\begin{pmatrix} \eta_{,\theta} - F \alpha_{,\theta} \\ \eta_{,\phi} - F \alpha_{,\phi} \end{pmatrix} =\begin{pmatrix} g_{\theta\theta} & g_{\theta\phi}\\ g_{\theta\phi} & g_{\phi\phi} \end{pmatrix} \begin{pmatrix} -\alpha_{,\phi} \\ \alpha_{,\theta} \end{pmatrix}, \end{equation}

which yields

(3.10)\begin{equation} \sqrt{g}F\begin{pmatrix} \eta_{,\theta} \\ \eta_{,\phi} \end{pmatrix}=\begin{pmatrix} g_{\theta\theta} & g_{\theta\phi} + F \sqrt{g}\\ g_{\theta\phi} - F\sqrt{g} & g_{\phi\phi} \end{pmatrix} \begin{pmatrix} -F\alpha_{,\phi} \\ +F\alpha_{,\theta} \end{pmatrix}.\end{equation}

This matrix equation can be inverted in order to obtain an equation for $\eta$:

(3.11)\begin{equation} \frac{\sqrt{g}F}{\mathcal{J}^2 + F^2 \sqrt{g}^2}\begin{pmatrix} g_{\phi\phi} & -g_{\theta\phi} - F \sqrt{g}\\ -g_{\theta\phi} + F\sqrt{g} & g_{\theta\theta} \end{pmatrix} \begin{pmatrix} \eta_{,\theta} \\ \eta_{,\phi} \end{pmatrix} = \begin{pmatrix} -F\alpha_{,\phi} \\ +F\alpha_{,\theta} \end{pmatrix}, \end{equation}

where we have defined the surface Jacobian to be $\mathcal {J} = \sqrt {g_{\theta \theta } g_{\phi \phi } - g_{\theta \phi }^2}$.

4. Near-surface expansions

Let us now expand about a surface $\psi = \psi _0$ such that

(4.1a)\begin{gather} \boldsymbol{r}(\psi,\theta,\phi) = \boldsymbol{r}^{(0)}(\theta,\phi) + \delta \psi\boldsymbol{r}^{(1)}(\theta,\phi) + {O}(\delta \psi^2), \end{gather}

(4.1b)\begin{gather}\alpha(\psi,\theta,\phi) = \alpha^{(0)}(\theta,\phi) + \delta \psi \alpha^{(1)}(\theta,\phi) + {O}(\delta\psi^2), \end{gather}

(4.1c)\begin{gather}\eta(\psi,\theta,\phi) = \eta^{(0)}(\theta,\phi) + \delta \psi \eta^{(1)}(\theta,\phi) + {O}(\delta\psi^2), \end{gather}

(4.1d)\begin{gather}\varPhi(\psi,\theta,\phi) = \varPhi^{(0)}(\theta,\phi) + \delta \psi \varPhi^{(1)}(\theta,\phi) + {O}(\delta\psi^2), \end{gather}

(4.1e)\begin{gather}F(\psi) = F^{(0)} + \delta \psi F^{(1)} + {O}(\delta\psi^2), \end{gather}

where $\delta \psi = \psi - \psi _0$. Note that we are using a standard Taylor expansion in linear power of $\psi$, as opposed to the expansion in $\sqrt {\psi }$ employed in the near-axis expansions. As we are expanding about a non-degenerate surface away from the magnetic axis (i.e. the two-dimensional surface is assumed to not degenerate to a one-dimensional curve), the coordinate singularity is avoided and a regular Taylor series suffices. We also assume enough smoothness in the various quantities appearing in (4.1) such that the expansions are valid.

At $\mathcal {O}((\delta \psi )^0)$ the vacuum equations (2.2) yield

(4.2a)\begin{gather} \sqrt{g}^{(1)} \begin{pmatrix} \varPhi^{(0)}_{,\theta}\\ \varPhi^{(0)}_{,\phi} \end{pmatrix}= \begin{pmatrix} g^{(0)}_{\theta\theta} & g^{(0)}_{\theta\phi}\\ g^{(0)}_{\theta\phi} & g^{(0)}_{\phi\phi} \end{pmatrix} \begin{pmatrix} -\alpha^{(0)}_{,\phi} \\ \alpha^{(0)}_{,\theta} \end{pmatrix}, \end{gather}

(4.2b)\begin{gather}\sqrt{g}^{(1)} \varPhi^{(1)}= \boldsymbol{r}^{(1)} \cdot \left( \boldsymbol{r}^{(0)}_{,\phi} \alpha^{(0)}_{,\theta} - \boldsymbol{r}^{(0)}_{,\theta} \alpha^{(0)}_{,\phi} \right), \end{gather}

where

(4.3a)\begin{gather} \sqrt{g}^{(1)}= \boldsymbol{r}^{(1)} \cdot \boldsymbol{r}^{(0)}_{,\theta} \times \boldsymbol{r}^{(0)}_{,\phi}, \end{gather}

(4.3b)\begin{gather}g^{(0)}_{x_ix_j} = \boldsymbol{r}^{(0)}_{,x_i} \cdot \boldsymbol{r}^{(0)}_{,x_j}, \end{gather}

for $x_i, x_j \in \{\theta ,\phi \}$.

We can define the angles in the neighbourhood of the surface such that the metric tensor simplifies (Boozer Reference Boozer2002; Imbert-Gerard, Paul & Wright Reference Imbert-Gerard, Paul and Wright2019):

(4.4a)\begin{equation} \boldsymbol{r}^{(1)} \cdot \boldsymbol{r}_{,\theta}^{(0)} = \boldsymbol{r}^{(1)} \cdot \boldsymbol{r}_{,\phi}^{(0)} = 0. \end{equation}

Defining $\rho ^{(1)} = |\boldsymbol {r}^{(1)}|$ and $\mathcal {J}^{(0)} = \sqrt {g}^{(1)}/\rho ^{(1)}$, we find that with such a choice of coordinates

(4.5)\begin{equation} \mathcal{J}^{(0)} = |\boldsymbol{r}^{(0)}_{,\theta} \times \boldsymbol{r}^{(0)}_{,\phi}|. \end{equation}

This eliminates the ‘out-of-surface’ equation for $\varPhi _{,\psi }$ (2.2b) locally but not globally. So our vacuum equations simplify to

(4.6)\begin{equation} \rho^{(1)}\mathcal{J}^{(0)}\begin{pmatrix} \varPhi^{(0)}_{,\theta}\\ \varPhi^{(0)}_{,\phi} \end{pmatrix} = \begin{pmatrix} g^{(0)}_{\theta\theta} & g^{(0)}_{\theta\phi}\\ g^{(0)}_{\theta\phi} & g^{(0)}_{\phi\phi} \end{pmatrix} \begin{pmatrix} -\alpha^{(0)}_{,\phi} \\ \alpha^{(0)}_{,\theta} \end{pmatrix}. \end{equation}

As shown in Appendix A, a global toroidal coordinate system $(\psi ,\theta ,\phi )$ can also be utilized to bring the ‘in-surface’ vacuum equations to the same form. In order to carry out a fully global calculation, one must retain the equation for $\varPhi _{,\psi }$ (2.2b). In a local expansion, this equation only gives higher-order corrections to $\varPhi$ and will not be considered in this work. To enforce QS, we need to rewrite the vacuum equations in terms of the variable $\eta$. Using

(4.7a)\begin{gather} \varPhi^{(0)} = \eta^{(0)} - F^{(0)} \alpha^{(0)}, \end{gather}

(4.7b)\begin{gather}\rho^{(1)}\mathcal{J}^{(0)} B^2(\eta^{(0)}) = \{\alpha^{(0)},\varPhi^{(0)}\}_{\theta,\phi}= \{\alpha^{(0)},\eta^{(0)} \}_{\theta,\phi} \end{gather}

to eliminate $\varPhi$, we obtain the following set of equations:

(4.8a)\begin{align} \begin{pmatrix} -F^{(0)}\alpha_{,\phi}^{(0)} \\ +F^{(0)}\alpha_{,\theta}^{(0)} \end{pmatrix} &= \frac{\rho^{(1)}F^{(0)}}{\mathcal{J}^{(0)} \left(1 + \left(F^{(0)}\rho^{(1)}\right)^2 \right)} \nonumber\\ &\quad\times\begin{pmatrix} g_{\phi\phi}^{(0)} & -g_{\theta\phi}^{(0)} - F^{(0)} \rho^{(1)} \mathcal{J}^{(0)}\\ -g_{\theta\phi}^{(0)} + F^{(0)}\rho^{(1)} \mathcal{J}^{(0)} & g_{\theta\theta}^{(0)} \end{pmatrix} \begin{pmatrix} \eta_{,\theta}^{(0)} \\ \eta_{,\phi}^{(0)} \end{pmatrix}, \end{align}

(4.8b)\begin{gather} \rho^{(1)}\mathcal{J}^{(0)} B^2(\eta^{(0)}) = \eta^{(0)}_{,\phi}\alpha^{(0)}_{,\theta}-\eta^{(0)}_{,\theta}\alpha^{(0)}_{,\phi}, \end{gather}

using (3.11). For the remainder of the text, the superscripts will be dropped for brevity. If the surface is axisymmetric, we can further choose the toroidal angle $\phi$ to coincide with the geometric toroidal angle. Therefore,

(4.9a,b)\begin{equation} g_{\phi\theta}=0,\quad \mathcal{J}=\sqrt{g_{\theta\theta}g_{\phi\phi}}. \end{equation}

We now make a coordinate transformation to an angle $\varTheta$ such that

(4.10)\begin{equation} \varTheta = \int \textrm{d} \theta \frac{\sqrt{g_{\theta \theta}}}{\sqrt{g_{\phi\phi}}}. \end{equation}

To see that $\varTheta$ is indeed an angle, we first note that the integrand $\sqrt {g_{\theta \theta }/g_{\phi \phi }}$ is a periodic function of only $\theta$ due to axisymmetry. We can therefore Fourier expand the integrand in harmonics of $\theta$ and integrate over $\theta$ to obtain $\varTheta$. The function $\varTheta (\theta )$ is clearly of the form $\varTheta =l_0\theta + \tilde {\varTheta }(\theta )$, where $l_0$ is the average of $\sqrt {g_{\theta \theta }/g_{\phi \phi }}$ and $\tilde {\varTheta }$ is periodic in $\theta$. Whereas physical quantities are required to have $2{\rm \pi}$ periodicity in the angle $\theta$, periodicity in $\varTheta$ is not necessarily $2{\rm \pi}$ and depends on the lowest-order surface. For numerical analysis we will use to $\theta$ and demand $2{\rm \pi}$ periodicity.

We rewrite the vacuum equations (4.6) and the QS condition (4.7b):

(4.11a)\begin{gather} \rho \begin{pmatrix} \varPhi_{,\varTheta}\\ \varPhi_{,\phi} \end{pmatrix}= \begin{pmatrix} -\alpha_{,\phi} \\ \alpha_{,\varTheta} \end{pmatrix}, \end{gather}

(4.11b)\begin{gather}\rho g_{\phi\phi} B^2 =\{\alpha,\varPhi\}_{\varTheta,\phi}. \end{gather}

The first two equations have the structure of the generalized Cauchy–Riemann condition and reduce to standard Cauchy–Riemann if $\rho =1$. The reduction to the above Cauchy–Riemann condition is the main advantage of using $\varTheta$ over $\theta$ as a coordinate. For the QS problem, $\rho$ must be self-consistently determined along with $\varPhi$ and $\alpha$. We note that $\varPhi$ and $\alpha$ are therefore pseudo-analytic functions (Bers Reference Bers1953).

We similarly write the equations in terms of $\eta$ and $\alpha$ (4.8) with this change of variables as

(4.12a)\begin{gather} \begin{pmatrix} -F\alpha_{,\phi} \\ +F\alpha_{,\varTheta} \end{pmatrix} = \frac{\rho F}{1 + \left(F\rho\right)^2 } \begin{pmatrix} 1 & -\rho F\\ \rho F & 1 \end{pmatrix} \begin{pmatrix} \eta_{,\varTheta}\\ \eta_{,\phi} \end{pmatrix}, \end{gather}

(4.12b)\begin{gather}\rho B^2(\eta) G(\varTheta) = \eta_{,\phi}\alpha_{,\varTheta}-\eta_{,\varTheta}\alpha_{,\phi}, \end{gather}

where we have defined $G(\varTheta ) = g_{\phi \phi }$.

The system (4.12) is a coupled PDE for the three unknowns ($\eta ,\alpha ,\rho$). The inputs are the functional form of the magnetic field strength, $B(\eta )$, and the geometric factor, $G(\varTheta )$, which encodes the surface shape. The helicity of the QS can also be chosen, which fixes the helicity of the streamlines of $\eta$. The outputs from (4.12) determine the field lines (through $\alpha (\varTheta ,\phi )$), the field strength (through $\eta (\varTheta ,\phi )$) and the Jacobian factor ($\rho (\varTheta ,\phi )$).

5. A parabolic PDE to enforce exact QS on a surface

We now simplify (4.12) to get a single PDE for $\eta$. Solving for $\rho F$ in terms of $\eta$ and $B(\eta )$ we obtain

(5.1)\begin{equation} 1+(\rho F)^2= \frac{\eta_{,\varTheta}^2+\eta_{,\phi}^2}{G(\varTheta)B(\eta)^2}. \end{equation}

Eliminating $\alpha$ by equating mixed partial derivatives, we now get a single nonlinear PDE for $\eta$:

(5.2)\begin{equation} (\partial^2_{\varTheta}+\partial^2_{\phi})\eta + \frac{2}{\tan{2\varOmega}}(\eta_{,\varTheta}\partial_\varTheta+ \eta_{,\phi}\partial_\phi)\varOmega+2(\eta_{,\varTheta}\partial_\phi-\eta_{,\phi}\partial_\varTheta)\varOmega =0, \end{equation}

which we refer to as the parabolic QS equation. Here,

(5.3)\begin{equation} \tan\varOmega =\rho F= \frac{\sqrt{\eta_{,\varTheta}^2+\eta_{,\phi}^2-G(\varTheta) B(\eta)^2}}{G(\varTheta)B(\eta)^2}. \end{equation}

The above equation when further simplified takes the form

(5.4)\begin{align} & (\rho F \eta_{,\varTheta}+\eta_{,\phi})^2\partial^2_\phi \eta+ 2(\rho F \eta_{,\varTheta}+\eta_{,\phi})(\eta_{,\varTheta}-\rho F \eta_{,\phi})\partial_\varTheta\partial_\phi \eta + (\eta_{,\varTheta}-\rho F \eta_{,\phi})^2\partial^2_\varTheta \eta \nonumber\\ &\quad =(1+(\rho F)^2)B^3\left( B'(\eta)(1-(\rho F)^2) + \frac{G'(\varTheta)}{2G(\varTheta)}\left(\eta_{,\varTheta}(1-(\rho F)^2)-2\rho F \eta_{,\phi}\right)\right). \end{align}

The right-hand side of the $\eta$ equation shows that the ‘source’ terms are the gradients of the magnetic field strength $B(\eta )$ and the geometry factor $G(\varTheta )$. We shall assume that these two given source terms are well defined and smooth enough for the solution $\eta (\varTheta ,\phi )$ to be meaningful. Given these two source terms, $B(\eta )$ and $G(\varTheta )$, we can in principle solve the nonlinear $\eta$ equation with the requirement that

(5.5)\begin{equation} \eta = a_1 \varTheta -a_2 \phi + \tilde{\eta}(\varTheta,\phi), \end{equation}

where $\tilde {\eta }$ is doubly periodic. We once again note that $\tilde {\eta }$ is periodic with respect to $\varTheta$, but with a periodicity that is generally not $2{\rm \pi}$. The transformed coordinate, $\varTheta$, could alternatively be rescaled to be $2{\rm \pi}$-periodic without impacting our results. The change of coordinates from $\theta$ to $\varTheta$ has no effect on the decomposition of $\eta$ into a secular and a periodic part in (5.5).

We note that the PDE that $\eta$ satisfies is of parabolic type in general geometry (see Appendix B) in contrast to the generalized Grad–Shafranov equation (Burby et al. Reference Burby, Kallinikos and MacKay2020) which is elliptic. As has been noted by Burby et al. (Reference Burby, Kallinikos and MacKay2020), the Grad–Shafranov equation is only part of the complete set of QS constraints: there are several additional non-trivial constraints ((94)–(97) in Burby et al. (Reference Burby, Kallinikos and MacKay2020)) that need to be satisfied besides the elliptic Grad–Shafranov equations. Addressing the compatibility of the local parabolic QS equation approach and the Grad–Shafranov approach with the additional constraints is left as a future work.

The repeated characteristic of the parabolic equation is given by

(5.6)\begin{equation} \frac{\textrm{d}\varTheta}{\textrm{d}\phi}= \frac{\eta_{,\varTheta}-\rho F \eta_{,\phi}}{\rho F \eta_{,\varTheta}+\eta_{,\phi}} ={-}\frac{F\alpha_{,\phi}}{F\alpha_{,\varTheta}}. \end{equation}

The last expression follows from using the vacuum equations.

Therefore, constant $\alpha$ lines are the repeated characteristics for the $\eta$ equation.

We shall now make a few changes of variables that considerably simplify the nonlinear parabolic equation. Firstly, to bring the parabolic equation to the standard form we shall transform to $(\varTheta ,F\alpha )$ coordinates. To further simplify the equations, in particular the source terms, we change variables from ($\eta ,\varTheta$) to ($\mathcal {U},\vartheta$), where

(5.7a,b)\begin{equation} \mathcal{U}= \int\frac{\textrm{d}\eta}{B(\eta)},\quad \vartheta= \int \textrm{d}\varTheta G(\varTheta). \end{equation}

Since $G$ is a physical function periodic in the angle $\theta$ and $\varTheta = l_0 \theta + \tilde {\varTheta }(\theta )$, Fourier expanding and integrating in $\theta$, we conclude that $\vartheta$ is indeed an angle of the form $\vartheta = l_1 \theta +\tilde {\vartheta }(\theta )$, where $l_1$ is a constant and $\tilde {\vartheta }$ is a periodic function of $\theta$. For arbitary $l_1$, $\vartheta$ does not have a $2{\rm \pi}$ period. The parabolic QS equation (5.4) in these variables reads

(5.8)\begin{equation} \mathcal{U}_{,\vartheta\vartheta}+\mathcal{U}_{,\vartheta}\left( \mathcal{U}_{,\vartheta}^2-1\right)\left( \frac{G'(\vartheta)}{2G(\vartheta)}+\frac{B'(\mathcal{U})}{B(\mathcal{U})}\frac{\mathcal{U}_{,\vartheta}}{1-B \mathcal{U}_{,F\alpha}}\right)=0.\end{equation}

The quantity $\rho$ is given by

(5.9)\begin{equation} \rho F ={\pm}\frac{\sqrt{\mathcal{U}_{,\vartheta}^2-1}}{1-B \mathcal{U}_{,F\alpha}}. \end{equation}

The positive and negative signs of $\rho$ represent outward and inward surface expansions corresponding to the assumption that the flux increases or decreases away from the axis. Equations (5.8) and (5.9) constitute the principal result of this work. The solution of (5.8) yields $\eta (\theta ,\phi )$ and thus the quasisymmetric field strength through $B(\eta )$. The necessary deformation of the axisymmetric geometry to support quasisymmetric fields is then provided by (5.9). Owing to $\vartheta$ having a period different from $2{\rm \pi}$, the domain of $\vartheta$ is $[0,l_1(2{\rm \pi} ))$, where $l_1$ depends on the surface about which we are expanding.

We now make a few observations from the structure of (5.8) and (5.9). Firstly, for $\rho$ to be real, we must have

(5.10)\begin{equation} |\mathcal{U}_{,\vartheta}|> 1. \end{equation}

The case $\mathcal {U}_{,\vartheta }=1$ is not allowed since $\rho =0$ implies that the flux surfaces overlap. Physically, $\rho =0$ would correspond to the formation of islands. Secondly, if the gradients of $G$ and $B$ are zero, $\eta =B \mathcal {U}$ and the general solution of $\mathcal {U}$ is linear in the angle $\vartheta$:

(5.11)\begin{equation} \mathcal{U} = u_1(F\alpha) \vartheta + u_2(F\alpha). \end{equation}

Since $\eta$ must satisfy the boundary condition (5.5), $u_1(F\alpha )$ must be a constant and $u_2(F\alpha )$ must be a periodic function of $F\alpha$. Furthermore, $u_1\neq 1$ to avoid $\rho =1$. We note that the case $B'=0$ corresponds to local isodynamic fields (Palumbo Reference Palumbo1968).

In the next section, we shall look for exact and approximate analytical solutions of the QS equation in slab, cylindrical and toroidal geometries. We show that there are critical differences in the nature of the QS solutions for slab and cylinder geometry for which $G$ is a constant and a real axisymmetric torus for which $G=G(\vartheta )$.

6. Solutions of the QS equation in different geometries

6.1. Slab and cylindrical geometries

In slab and cylindrical geometries the geometry factor $G$ is unity (Appendix A). Hence, the only source term in the QS equation is the gradient of $B$. Exact travelling wave solutions to the $\eta$ equation (5.4) can be constructed such that $\eta =\eta (\varTheta -c \phi )$ for constant $c$. The resulting ordinary differential equation (ODE) in the variable $(\varTheta -c \phi )$ is

(6.1)\begin{equation} \eta'' -\frac{B'(\eta)}{B(\eta)}\eta'^2\left( 2-(1+c^2) \left(\frac{\eta'}{B}\right)^2\right)=0, \end{equation}

which has an exact solution given by

(6.2)\begin{equation} \int \frac{\textrm{d}\eta}{k B(\eta)^2}\sqrt{1+(1+c^2)(kB(\eta))^2}={\pm} (\varTheta-c \phi), \end{equation}

where $k$ is the integration constant. Since $B(\eta )$ is a periodic function of $\eta$, the integral on the left-hand side defines a function which has periodic terms in $\eta$ and at most a linear secular term in $\eta$ arising from the average of the integrand. This can be readily seen by Fourier expanding the integrand in $\eta$. Therefore, when the function of $\eta$ is inverted, we get a solution which satisfies the boundary condition (5.5).

Similarly, there exist exact solutions of the form $\mathcal {U}=\mathcal {U}(\vartheta -c F\alpha )$ to the QS equation (5.8) which reduces to the following ODE in the variable $(\vartheta -c F\alpha )$:

(6.3)\begin{equation} \mathcal{U}''+\mathcal{U}'\left( {\mathcal{U}'}^2-1\right)\left( \frac{B'(\mathcal{U})}{B(\mathcal{U})}\frac{\mathcal{U}'}{1+c B \mathcal{U}'}\right)=0. \end{equation}

The ODE has the exact solution

(6.4)\begin{equation} -\int \textrm{d}\mathcal{U} \frac{ \pm\sqrt{k(1-(c B(\mathcal{U}))^2)+1}+(c B(\mathcal{U}))^2}{c B(\mathcal{U}) \left({\pm}\sqrt{k (1-(c B(\mathcal{U}))^2)+1}+1\right)}= \vartheta-c F\alpha, \end{equation}

with $k$ as an integration constant. We observe that there are two possible travelling wave solutions in each of the above cases due to the quadratic nature of the equations. This can be interpreted as solutions for which the toroidal flux increases in the direction towards or away from the axis. The existence of two linearly independent solutions is also obtained in the system of equations that arises in an asymptotic expansion near axisymmetry (Plunk & Helander Reference Plunk and Helander2018; Plunk Reference Plunk2020).

6.2. Toroidal geometry with arbitrary aspect ratio and arbitrary shaping

If $G'(\vartheta )\neq 0$, exact travelling wave solutions do not exist since $G(\vartheta )$ breaks the travelling wave symmetry. We shall treat the QS equations in this case numerically in § 7. However, analytical progress can still be made in several interesting cases. In the following we give the analytical details of these cases to provide physical insight into the nature of the QS solutions and also to serve as benchmarks for the numerical solution.

6.2.1. Complete axisymmetry

The complete axisymmetric case for arbitrary aspect ratio and shaping can be tackled analytically. In this limit we have $G=G(\vartheta )$, $\mathcal {U}= \mathcal {U}( \vartheta )$, $B=B( \mathcal {U})$ and $\rho =\rho (\vartheta )$. The QS equation reduces to

(6.5)\begin{equation} \mathcal{U}''+ \mathcal{U}'( \mathcal{U}'^2-1)\left( \frac{G'}{2G}+\frac{B'}{B} \mathcal{U}' \right)=0. \end{equation}

The above equation can be integrated once and we obtain

(6.6)\begin{equation} G B^2 \left( 1- \frac{1}{ \mathcal{U}'^2} \right) = c_0, \end{equation}

where $c_0$ is the integration constant. We can rewrite the above equation in the form of an energy equation with a potential $V(\mathcal {U},\vartheta )$:

(6.7a,b)\begin{equation} \frac{1}{2} \mathcal{U}'^2+V( \mathcal{U}, \vartheta)=0, \quad V( \mathcal{U}, \vartheta)= \frac{1}{2}\frac{G B^2}{c_0-G B^2}. \end{equation}

From (5.9) it follows that

(6.8)\begin{equation} \rho F ={\pm} \frac{\sqrt{c_0}}{\sqrt{G B^2-c_0}}. \end{equation}

6.2.2. Isodynamic fields ($B=B_0,B'=0$)

The QS equation reduces to the following ODE in the variable $\vartheta$ when $B=B_0$ is constant on the flux surface:

(6.9)\begin{equation} \mathcal{U}''+\mathcal{U}'({\mathcal{U}'}^2-1)\frac{G'}{2G}=0, \end{equation}

which has an exact solution

(6.10)\begin{equation} \mathcal{U} ={\pm} \int \textrm{d}\vartheta \sqrt{\frac{G(\vartheta) B^2_0}{G(\vartheta) B^2_0-c_0}} + \varGamma F\alpha+ \bar{\mathcal{U}}(F\alpha). \end{equation}

The function $\bar {\mathcal {U}}(F\alpha )$ is an arbitrary periodic function of $F\alpha$ and $\varGamma$ is a constant. The corresponding value of $\rho F$ is

(6.11)\begin{equation} \rho F ={\pm} \sqrt{\frac{c_0}{G B^2_0-c_0}}( 1-B_0(\varGamma+\bar{\mathcal{U}}'(F\alpha)) )^{{-}1}. \end{equation}

Palumbo's original isodynamic solution (Palumbo Reference Palumbo1968) was axisymmetric, but locally we can always find non-axisymmetric isodynamic fields provided $\varGamma \neq 0$ and $\bar {\mathcal {U}}(F\alpha )$ is not a constant.

6.2.3. Localization near $B'=0$

We shall now look at the behaviour of a non-trivial non-axisymmetric solution in the neighbourhood of $B'=0$ which we assume happens along $\mathcal {U}=\mathcal {U}_m$. We shall restrict ourselves to the region $|\delta \mathcal {U}|=|\mathcal {U}-\mathcal {U}_m|\ll 1$. Even though $\delta U\ll 1$, we do not expect $\delta U'$ to be small. Therefore, we can have sharp variations in the solution near $B'=0$.

In the region near $B'=0$, we expect the solution to be close to the isodynamic case given by (6.10). We therefore make the ansatz that $\mathcal {U}$ is of the form

(6.12)\begin{equation} \mathcal{U}=\mathcal{U}_m + \varGamma F\alpha + h(\vartheta), \end{equation}

where $\varGamma$ is a constant. Taylor expanding $B(\mathcal {U})$ in powers of $\delta \mathcal {U}$ and keeping only up to the quadratic power in $\delta \mathcal {U}$ such that $B=B_0 + B_2 (\delta \mathcal {U})^2+\cdots$, and keeping all nonlinearities in $\delta \mathcal {U}'$, we get

(6.13)\begin{equation} \delta\mathcal{U}''+ \delta\mathcal{U}'( (\delta\mathcal{U}')^2-1)\left( \frac{G'}{2G}+2\frac{B_2 \delta\mathcal{U}}{B_0}\frac{ \delta\mathcal{U}'}{1-B_0\varGamma} \right)=0, \end{equation}

where prime indicates differentiation with respect to $\vartheta$. The above equation can be integrated once such that

(6.14)\begin{equation} G B_0^2\left( 1- \frac{1}{ \delta \mathcal{U}'^2} \right) \exp{\left(\frac{2B_2}{B_0}\frac{ \delta \mathcal{U}^2}{1-B_0 \varGamma}\right)} = c_0,\end{equation}

We can cast the above expression in a energy equation as follows:

(6.15a,b)\begin{align} \frac{1}{2} \delta \mathcal{U}'^2+V( \delta\mathcal{U}, \vartheta)=0, \quad V( \delta\mathcal{U}, \vartheta)= \frac{1}{2}\left(\frac{c_0}{B_0^2 G}\exp\left({-\left(\frac{2B_2}{B_0} \frac{ \delta \mathcal{U}^2}{1-B_0 \varGamma}\right)}\right)-1\right)^{{-}1}. \end{align}

We shall show later that the solution has sharp peaks near the minimum of the magnetic field strength, and we compare with the numerically obtained solution in § 7.

6.2.4. Near-isodynamic limit

We now briefly discuss a bifurcation that can happen when the background geometry $G(\vartheta )$ and magnetic field strength $B_0(\vartheta )$ are axisymmetric. If we now look for a small deformation such that magnetic field is quasisymmetric, i.e. $B=B(\eta )$, we can have the following situations:

1. (i) Both the field strength and $\rho$ are axisymmetric.

2. (ii) Both the field strength and $\rho$ are non-symmetric.

Case (i) implies that $\eta =\eta (\vartheta )$ and $\rho =\rho (\vartheta )$. We have analysed this case completely in § 6.2.1. Case (ii) implies that $\eta$ must have dependence on both $\vartheta$ and $F\alpha$, so that $\rho F$, given by (5.9), is non-symmetric. We now look for solutions of the QS equation (5.8) in the form

(6.16a)\begin{gather} \mathcal{U}=\mathcal{U}_0(\vartheta,F\alpha)+\epsilon\mathcal{U}_1(\vartheta,F\alpha)+{O}(\epsilon^2), \end{gather}

(6.16b)\begin{gather}B=B( \mathcal{U})=B_0+\epsilon B_1(\mathcal{U})+{O}(\epsilon^2), \end{gather}

(6.16c)\begin{gather}\rho=\rho_0+\epsilon \rho_1+{O}(\epsilon^2), \end{gather}

where $B_0$ is a constant background field. The expansion parameter $\epsilon$ is similar to the large aspect ratio parameter; however, we do not expand the geometry factor. We further require that to lowest order, $\mathcal {U}$ depends on $F\alpha$ only through a secular term $\varGamma F\alpha$:

(6.17)\begin{equation} \mathcal{U}_0(\vartheta,F\alpha)= \varGamma F\alpha +h_0( \vartheta). \end{equation}

For quasi-axisymmetry we must have $\varGamma =0$.

The lowest-order calculations are analogous to the case described in § 6.2.3, since $B=B_0$ is constant to this order. Denoting $\vartheta$ derivatives by primes, (5.8) to lowest order implies

(6.18)\begin{equation} h''_0+h'_0({h'_0}^2-1) \frac{G'}{2G}=0, \end{equation}

which leads to the energy equation

(6.19a,b)\begin{equation} \frac{1}{2}{h'_0}^2+V(h_0, \vartheta)=0, \quad V(h_0, \vartheta)= \frac{1}{2}\frac{G B_0^2}{c_0-G B_0^2}, \end{equation}

where the integration constant $c_0$ is a periodic function of $F\alpha$. It is convenient to separate $c_0$ into a constant part $C_0$ and a periodic part $\tilde {c}_0$. The quantity $\rho _0 F$ is given by

(6.20)\begin{equation} \rho_0 F = \sqrt{\frac{c_0}{G B_0^2-c_0}}\frac{1}{1-B_0\varGamma}. \end{equation}

To $O(\epsilon )$, we obtain the following linear equation for $\mathcal {U}_1$:

(6.21)\begin{equation} \left( \mathcal{U}'_1\sqrt{\frac{(G B^2_0-c_0)^3}{G B^2_0}}\right)'= \frac{c_0}{1-B_0\varGamma}\frac{B'_1(\mathcal{U}_0)}{B_0}\sqrt{\frac{GB_0^2}{GB^2_0-c_0}}. \end{equation}

We note that since $\mathcal {U}'_1$ is periodic and $(B_0,\varGamma ,c_0)$ are all independent of $\vartheta$, the right-hand side of (6.21) must satisfy the condition

(6.22)\begin{equation} \oint \textrm{d}\vartheta \frac{B'_1(\mathcal{U}_0)}{B_0}\sqrt{\frac{GB_0^2}{GB^2_0-c_0}}=0. \end{equation}

The above constraint is satisfied if both $B_1$ and $G$ are stellarator symmetric, for example. Solving for $\mathcal {U}_1$ we get

(6.23)\begin{equation} \mathcal{U}_1=c_1(F\alpha)\sqrt{\frac{G B^2_0}{(G B^2_0-c_0)^3}}+ \frac{c_0}{1-B_0\varGamma}\int \textrm{d}\vartheta \frac{B'_1(\mathcal{U}_0)}{B_0}\sqrt{\frac{GB_0^2}{GB^2_0-c_0}}, \end{equation}

where $c_1$ is a function of $F\alpha$ to be determined by going to the next order.

The second-order equation is once again of the form (6.21) and the periodicity requirement on $\mathcal {U}_2$ leads to a Riccati equation for $c_1$. The calculation is straightforward and we shall skip the details. The Riccati equation for $c_1$ is of the form

(6.24)\begin{equation} \frac{\textrm{d} c_1}{\textrm{d}(F\alpha)}+2 a_1 c_1 + a_2 c_1^2+a_3=0, \end{equation}

where $a_1$, $a_2$ and $a_3$ are complicated functions of $c_0(F\alpha )$. The function $c_1$ has to be a periodic function of $F\alpha$. Therefore, we must have

(6.25)\begin{equation} \oint \textrm{d} F\,\alpha \left( 2 a_1 c_1 + a_2 c_1^2+a_3 \right)=0.\end{equation}

Assuming that $c_0= C_0 + \tilde {c}_0$, we can determine $C_0$ using (6.25) for a given periodic function $\tilde {c}_0$. We note that if $\tilde {c}_0 =0$, i.e. $c_0$ is a constant, the Riccati equation (6.24) has constant coefficients and can be solved exactly as

(6.26a,b)\begin{equation} c_1=\frac{a_1-\varDelta_0 \tan \left(\varDelta_0(F\alpha-F\alpha_0)\right)}{a_2}, \quad \varDelta_0=\sqrt{a_2 a_3-a_1^2}. \end{equation}

Thus, we see that unless $\varDelta _0=0$, $c_1$ is not periodic when $c_0$ is a constant. When $\varDelta _0=0$, we obtain an axisymmetric solution.

6.2.5. A model for QS fields near an axisymmetric surface

The previous analysis suggests a strong correlation between $\rho F$ and the product of the geometry factor, $G$, and the magnetic field strength, $B^2$. We construct a simple model of QS that captures most of the interesting properties for a large class of travelling-wave-like solutions for arbitrary geometry and field strength factors:

(6.27a)\begin{gather} \mathcal{U} \approx{\pm} \int \textrm{d}\vartheta \sqrt{\frac{G B^2}{G B^2-c_0}} + \varGamma F\alpha, \end{gather}

(6.27b)\begin{gather}\rho F \approx{\pm} \sqrt{\frac{c_0}{G B^2-c_0}}\left( 1-B \varGamma \right)^{{-}1}. \end{gather}

This model is motivated by the solutions obtained in the near-isodynamic limit (6.17)–(6.20). In § 7, we verify this model numerically.

7. Numerical solutions of the QS equations

Although full toroidal geometry can be treated analytically in a few limiting cases, we obtain numerical solutions for the general case. To do so, we construct a spectral solution for $\alpha$, $\eta$ and $\rho$:

(7.1a)\begin{gather} \alpha = \theta - \iota \phi + \sum_{i=1}^{N_{\text{modes}}} \alpha_{m_i,n_i} \sin(m_i \theta- n_i \phi), \end{gather}

(7.1b)\begin{gather}\eta = \theta - N/M \phi + \sum_{i=1}^{N_{\text{modes}}} \eta_{m_i,n_i} \sin(m_i \theta - n_i \phi), \end{gather}

(7.1c)\begin{gather}\rho = \sum_{i=1}^{N_{\text{modes}}+1}\rho_{m_i,n_i}\cos(m_i\theta - n_i\phi), \end{gather}

where we have assumed stellarator symmetry and chosen a normalization such that the secular term in $\eta$ is unity. We note that $N/M$ must be a rational number, which determines the helicity of the QS. A nonlinear set of equations is obtained by inserting (7.1) into (4.8) and integrating against a set of $N_{\text {modes}}+1$ basis functions:

(7.2)\begin{equation} \{x_{m_j,n_j}\} = \{\cos(m_j \theta - n_j \phi)\}. \end{equation}

The set of $3N_{\text {modes}}+3$ unknowns is taken to be $\{\iota ,\alpha _{m,n},F,\eta _{m,n},\rho _{m,n}\}$ while the helicity $N/M$ is prescribed in addition to the functional form of $B(\eta )$. The resulting nonlinear set of equations is solved with a trust region method with a prescribed analytic Jacobian matrix. In the following sections, we verify the numerical solutions in several limits.

Once the solution is obtained, the Boozer angles can be computed as

(7.3a)\begin{gather} \vartheta_B = \theta + \frac{-N/M \tilde{\alpha} + \iota \tilde{\eta} }{\iota - N/M}, \end{gather}

(7.3b)\begin{gather}\varphi_B = \phi + \frac{\tilde{\eta} - \tilde{\alpha}}{\iota - N/M}, \end{gather}

such that only secular terms remain in the potentials:

(7.4a)\begin{gather} \varPhi(\vartheta_B,\varphi_B) = \bar{G} \varphi_B + \bar{I} \vartheta_B, \end{gather}

(7.4b)\begin{gather}\alpha(\vartheta_B,\varphi_B) = \vartheta_B - \iota \varphi_B. \end{gather}

Note that we have not assumed $\bar {I} = 0$ as is the case for a global vacuum field. Thus we allow for the existence of currents asymptotically away from the surface.

7.1. Complete axisymmetry

In the case of pure axisymmetry, we can confirm that solutions of (4.8) agree with numerical solutions of the Grad–Shafranov equation on a surface. To do so, we compute an axisymmetric equilibrium with the VMEC code (Hirshman & Whitson Reference Hirshman and Whitson1983). We consider a boundary given by

(7.5a)\begin{gather} R(\theta,\phi) = 3 + \cos(\theta), \end{gather}

(7.5b)\begin{gather}Z(\theta,\phi) = \sin(\theta), \end{gather}

and impose toroidal current and pressure profiles which are peaked near the axis,

(7.6a)\begin{gather} I_T'(\psi/\psi_0) = I_0 \exp({-(\psi/\psi_0)^2/0.05^2}), \end{gather}

(7.6b)\begin{gather}p(\psi/\psi_0) = p_0 \exp({-(\psi/\psi_0)^2/0.05^2}), \end{gather}

to provide a non-zero rotational transform. The magnetic field on the surface at $\psi /\psi _0 = 0.995$ is used to construct $\rho$, $\varPhi$, $\alpha$ and $\eta$ using (3.7) with $N = 0$ and $M = 1$. A Fourier transform is then performed to obtain the functional form of $B(\eta )$. This solution is provided as an initial guess to the numerical solver of the discretized equations (4.8), including modes with $m \le 20$, which converges in three iterations to a function tolerance of $10^{-8}$. Given the numerical solution, we compute the integration constant $c_0$ from (6.6). The relative variation of $c_0$ with respect to $\theta$ is $3.5\times 10^{-9}$, and its average over $\theta$ is used to evaluate (6.8). The value of $\rho F$ computed in this way is displayed as the ‘analytic’ result in figure 1. This is compared with the value computed from the VMEC equilibrium and that obtained from the nonlinear solution. As can be seen, the three solutions are in good agreement. The averaged relative error,

(7.7)\begin{equation} \text{error} = \frac{\displaystyle\int \textrm{d} \theta | (\rho F)^2 - ((\rho F)^{\text{model}})^2 | }{\displaystyle\int \textrm{d} \theta \, (\rho F)^2}, \end{equation}

between the VMEC and nonlinear solution is $2.49\times 10^{-5}$, and the averaged residual between the analytic and nonlinear solution is $1.27\times 10^{-8}$.

Figure 1.

Benchmark of pure axisymmetric solution (6.8).

7.2. Bifurcated solutions

We now consider solutions near an axisymmetric solution for which the field strength and $\rho$ are non-axisymmetric. In order to obtain such numerical solutions, we provide an initial guess, which is a travelling wave solution, $\eta (c\varTheta -N/M \phi )$. Here the constants $F$ and $c$ defining the initial condition are chosen such that the secular terms in $\theta$ are unity as in (7.1) and $N/M$ is taken to match the desired helicity.

With this method, we are able to obtain a class of helically symmetric solutions near an axisymmetric surface. We prescribe the field strength to be of the form

(7.8)\begin{equation} B^2(\eta) = 1 + \epsilon \cos(\eta) \end{equation}

and consider the boundary given by (7.5). We present a solution with $\epsilon = 0.3$ in figures 2 and 3. The nearby surfaces are plotted with a finite value of $\varDelta \psi = 0.05$:

(7.9a)\begin{gather} R^{{\pm}}(\theta,\phi) = R(\theta,\phi) \pm \varDelta \psi (\hat{\boldsymbol{n}} \cdot \hat{\boldsymbol{R}} ) \rho(\theta,\phi), \end{gather}

(7.9b)\begin{gather}Z^{{\pm}}(\theta,\phi) = Z(\theta,\phi) \pm \varDelta \psi (\hat{\boldsymbol{n}} \cdot \hat{\boldsymbol{z}}) \rho(\theta,\phi). \end{gather}

We note that the surface perturbations are localized on the inboard size of the torus, with a deformation that results in convexity and concavity in the neighbouring surfaces. We discuss the nature of this trend using the models in § 6.2 shortly. Although the background geometry is axisymmetric, we obtain a helical pattern of field strength. As expected, both the field lines and the contours of the field strength become straight when plotted in the Boozer coordinates $(\vartheta _B,\varphi _B)$. In contrast, $\rho$ does not exhibit a clear symmetry direction when plotted in the Boozer coordinates. However, the contours appear to be straightened due to the dependence of $\rho$ on the quantity $G B^2$, which will be discussed below.

Figure 2.

Numerical solution for (a) $\rho (\theta ,\phi )$, (b) $B^2(\theta ,\phi )$, (c) $\rho (\vartheta _B,\varphi _B)$ and (d) $B^2(\vartheta _B,\varphi _B)$ computed imposing the functional form of the field strength given by (7.8) with $\epsilon = 0.3$ near the axisymmetric surface given by (7.5). The red lines in (c,d) indicate field lines.

Figure 3.

Using the solution presented in figure 2, we plot the the surface of expansion given by (7.5) (yellow) and the nearby surfaces (red and blue) given by (7.9).

We now compare the models in § 6.2 with this numerical solution. To compare with the model near $B'(\eta ) = 0$ in § 6.2.3, we expand near the minimum of $B(\eta )$ at $\eta = \eta _0 = {\rm \pi}$ and a chosen value of $\alpha = \alpha _0$. The constant $\varGamma$ is obtained using the value of $\mathcal {U}_{,F \alpha }$ at $(\alpha _0,\eta _0)$. We then note that

(7.10)\begin{equation} B^2(\delta \mathcal{U}^2) \approx B_0^2 + 2 B_2 B_0 \delta \mathcal{U}^2, \end{equation}

so that the exponential in (6.14) can be expressed in terms of $B_0$ and $B^2$ rather than $B_2$. The constant $c_0$ is then evaluated using (6.14) with the solution for $\mathcal {U}_{,\vartheta }$ and $B^2$ at $(\alpha _0,\eta _0)$. The same expression is then used to model $\mathcal {U}_{,\vartheta }$ in the neighbourhood of $(\alpha _0,\eta _0)$, which is denoted by a black star in figure 4. This local model as well as the numerical solution are shown. For the range of $\alpha$ and $\eta$ displayed ($\alpha \in [\alpha _0-{\rm \pi} /2,\alpha _0+{\rm \pi} /2]$ and $\eta \in [\eta _0 - {\rm \pi} /20,\eta _0+{\rm \pi} /20]$), the averaged relative error,

(7.11)\begin{equation} \text{error} = \frac{\displaystyle\int \textrm{d} \alpha \int \textrm{d} \eta \left| \mathcal{U}_{,\vartheta}^2 - \left(\mathcal{U}_{,\vartheta}^{\text{model}}\right)^2 \right| }{\displaystyle\int \textrm{d} \alpha \int \textrm{d} \eta \, \mathcal{U}_{,\vartheta}^2}, \end{equation}

is $1.8\times 10^{-3}$.

Figure 4.

We compare the solution presented in figure 2 with the near-minimum (§ 6.2.3) and near-isodynamic (§ 6.2.4) models. The model is constructed by expanding about a point indicated by the black star.

We also compare with the near-isodynamic limit (6.27). We expand about the same point $(\alpha _0,\eta _0)$ and again fix $\varGamma$ to be the value of $\mathcal {U}_{,F\alpha }$ at this point. The constant $c_0$ is then fixed using the value of $\mathcal {U}_{,\vartheta }$ at this point, and (6.27) is used to model $\mathcal {U}_{,\vartheta }$ in the neighbourhood of $(\alpha _0,\eta _0)$. The result is shown in figure 4. The averaged error between the numerical solution and the model is $2.0\times 10^{-3}$.

While the near-minimum and near-isodynamic models are only valid in the neighbourhood of the minimum of the field strength, the models allow us to glean insight into the global nature of the solution. In particular, we note that the spatial dependence of $\mathcal {U}_{,\vartheta }$ and $\rho$ is largely through the combination $G B^2$. In figure 5 we present the solution for $\rho (\theta ,\phi )$ for several values of $\epsilon$. We note that as $\epsilon$ is increased, the perturbation of the nearby surfaces becomes highly localized, introducing sharp helical structures in $\rho$. For smaller values of $\epsilon$, the maximum of $\rho$ becomes localized on the inboard side where $G$ is minimized. In the limit that $\epsilon \rightarrow 0$, $\rho (\theta ,\phi )$ becomes axisymmetric. The black contours overlaid on the colour-scale plots in figure 5 are contours of $G B^2$, with the minimum of $GB^2$ coinciding with the maximum of $\rho$. While the contours of $GB^2$ do not exactly coincide with the contours of $\rho$ for $\epsilon >0$, we note that much of the structure of the numerical solution is captured. The results we present here are consistent with features of near-axisymmetric quasi-axisymmetric solutions (Plunk & Helander Reference Plunk and Helander2018; Plunk Reference Plunk2020), for which the non-axisymmetric perturbations to the field were localized to the inboard side. In figure 5(f), the major radius of the boundary given by (7.5) is increased to $R_0 = 100$. In this case, $G(\vartheta )$ approaches a constant such that the perturbation to the surfaces is no longer localized to the inboard side. The large aspect ratio toroidal solution then approaches the travelling wave (helical) solutions of the slab and cylindrical geometries.

Figure 5.

We present numerical solutions using the functional form of the field strength given by (7.8) with the displayed value of $\epsilon$ near the axisymmetric surface given by (7.5): (a) $\epsilon = 0$; (b) $\epsilon = 0.1$; (c) $\epsilon = 0.2$; (d) $\epsilon = 0.3$; (e) $\epsilon = 0.4$; (f) $\epsilon = 0.3$, $R_0 = 100$. The colour scale indicates the magnitude of $\rho (\theta ,\phi )$, with the overlaid black lines indicating the contours of $GB^2$. In (f), the boundary given by (7.5) is modified such that the major radius is $R_0 = 100$.