2 Ideal magnetohydrodynamics with the quasisymmetry constraint

We begin with the well-known model for a plasma equilibrium, the ideal MHD equations

(2.1a)\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{B}=0 , \end{gather}

(2.1b)\begin{gather}\boldsymbol{J}\times \boldsymbol{B} = \boldsymbol{\nabla} p , \end{gather}

(2.1c)\begin{gather}\boldsymbol{J}=\boldsymbol{\nabla} \times \boldsymbol{B}. \end{gather}

Here, $\boldsymbol {B}$, $\boldsymbol {J}$ and $p$ denote the magnetic field vector, the current density and the plasma pressure, respectively. We have used units such that the permeability of free space, $\mu _0=1$. We will assume that pressure is constant on a set of nested toroidal flux surfaces labelled by $\psi$, i.e.

(2.2a,b)\begin{equation} p=p(\psi), \quad \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\psi=0. \end{equation}

We shall further assume that $\boldsymbol {\nabla } \psi$ is non-zero almost everywhere except on a finite set of closed field lines that include the magnetic axis and the separatrix.

The QS constraint can be imposed on the ideal MHD equilibrium in many different but equivalent ways. Here, we shall make use of the following form, which is most commonly referred to as the two-term form (Helander Reference Helander2014; Burby et al. Reference Burby, Kallinikos and MacKay2020; Rodríguez et al. Reference Rodríguez, Helander and Bhattacharjee2020):

(2.3)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla} B= \frac{1}{F(\psi)}\boldsymbol{B} \times \boldsymbol{\nabla} \psi \boldsymbol{\cdot} \boldsymbol{\nabla} B . \end{equation}

The flux function $\psi$ denotes the toroidal flux, $B$ is the magnetic field strength and $F(\psi )$ is a flux function that can be shown (Helander Reference Helander2014) to be related to the rotational transform $\iota (\psi )$, through

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

Here, $N/M$ is the helicity, and $G, I$ are poloidal and toroidal currents. For quasiaxisymmetry (QA), which is the subject of this work, $N=0$. The currents can be obtained by integrating $\boldsymbol {B}$ along the constant $\vartheta$ (poloidal) and constant $\varphi$ (toroidal) angles. These angles can be chosen as the Boozer angles for convenience. Thus,

(2.5a,b)\begin{equation} G(\psi)=\frac{1}{2{\rm \pi}}\oint_\vartheta \boldsymbol{B}\boldsymbol{\cdot} {\rm d}\boldsymbol{r}, \quad I(\psi)=\frac{1}{2{\rm \pi}}\oint_\varphi \boldsymbol{B}\boldsymbol{\cdot} {\rm d}\boldsymbol{r} . \end{equation}

Another relevant quantity of interest is the QS vector $\boldsymbol {u}$ (Burby et al. Reference Burby, Kallinikos and MacKay2020; Rodríguez et al. Reference Rodríguez, Helander and Bhattacharjee2020), which may be defined by the following two equivalent relations:

(2.6a,b)\begin{equation} \boldsymbol{u} = \frac{\boldsymbol{\nabla} \psi\times \boldsymbol{\nabla} B}{\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{\nabla} B} ,\quad \boldsymbol{u} = \frac{1}{B^2}\left( F(\psi)\boldsymbol{B} -\boldsymbol{B} \times \boldsymbol{\nabla} \psi \right) . \end{equation}

We note here that the QS vector of Burby et al. differs by a factor of iota in the quasi-axisymmetric case assumed here. Therefore, the two-term QS relation is equivalent to $\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {B} =F(\psi )$. The QS vector $\boldsymbol {u}$ defines curves that lie on constant flux surfaces along which $B$ does not change since $\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla } B=0$ and $\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla } \psi =0$. Furthermore, $\boldsymbol {u}$ lines are closed since $B$ is single valued.

Let us now demonstrate how QS helps ensure that pressure-driven singular currents do not generally form on rational surfaces. We can obtain the ideal current from (2.1b) in the form

(2.7)\begin{equation} \boldsymbol{J}=j_{||} \boldsymbol{B} + \frac{\boldsymbol{B}\times \boldsymbol{\nabla} p}{B^2}, \end{equation}

where the parallel component $j_{||}$ is not yet determined. To determine $j_{||}$, we use the fact that the current must be divergence free due to (2.1c). This leads to the following consistency condition:

(2.8)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} j_{||} + \boldsymbol{B} \times \boldsymbol{\nabla} p\boldsymbol{\cdot} \boldsymbol{\nabla} \left( \frac{1}{B^2} \right)=0. \end{equation}

Equation (2.8) is a magnetic differential equation for $j_{||}$. Single valuedness of $j_{||}$ leads to the following Newcomb constraint on every closed field line:

(2.9)\begin{equation} \oint {\rm d}\ell\,\boldsymbol{B} \times \boldsymbol{\nabla} p\boldsymbol{\cdot} \boldsymbol{\nabla} \left( \frac{1}{B^2} \right)=0. \end{equation}

Without a continuous symmetry such as axisymmetry, the Newcomb condition is generally not satisfied on all rational surfaces, permitting the formation of singular currents in 3-D ideal MHD equilibria.

The QS condition, (2.3), implies that

(2.10)\begin{equation} \boldsymbol{B} \times \boldsymbol{\nabla} p\boldsymbol{\cdot} \boldsymbol{\nabla} \left( \frac{1}{B^2} \right)=\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\left( F(\psi) \frac{p'(\psi)}{B^2} \right), \end{equation}

which satisfies the Newcomb constraint. This also allows us to solve for $j_{||}$ up to a flux function $H'(\psi )$

(2.11)\begin{equation} j_{||} ={-}H'(\psi) -F(\psi) \frac{p'(\psi)}{B^2}. \end{equation}

We note here that the pressure-driven term in (2.11) constitutes the Pfirsch–Schluter current, whereas the flux function, $H'(\psi )$, encompasses all other non-pressure-driven current sources. Since a $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }$ has been lifted to obtain $j_{||}$ from (2.8), singular Dirac-delta currents (Loizu et al. Reference Loizu, Hudson, Bhattacharjee, Lazerson and Helander2015; Rodríguez & Bhattacharjee Reference Rodríguez and Bhattacharjee2021; Huang et al. Reference Huang, Hudson, Loizu, Zhou and Bhattacharjee2022, Reference Huang, Zhou, Loizu, Hudson and Bhattacharjee2023) are possible as well on rational surfaces. We plan to address the singular currents in the future.

In the case of axisymmetry, $j_{||}$ can be written as an elliptic Laplacian-like operator acting on the poloidal flux, leading to a Grad–Shafranov (GS) equation. We would now like to derive a QS-GS equation in analogy with axisymmetry. The general QS-GSE (Burby et al. Reference Burby, Kallinikos and MacKay2020) (discussed in Appendix B) being too complicated, we shall use the reduced model for ideal MHD derived by Strauss (Reference Strauss1997) in the remainder of this section. Strauss’ model assumes that the magnetic field is approximately given by

(2.12a,b)\begin{equation} \boldsymbol{B} = \boldsymbol{B}_v + \boldsymbol{\nabla} A\times \boldsymbol{B}_v, \quad \boldsymbol{B}_v= \boldsymbol{\nabla} \chi, \end{equation}

where $\boldsymbol {B}_v$ is a vacuum field, $A$ is an $O( L_\perp /L_\parallel )$ function that is analogous to the poloidal flux in the axisymmetric case. The fundamental assumption is that the gradients along $\boldsymbol {B}_v$ are smaller than those perpendicular to it. The form (2.12a,b) ensures that $\boldsymbol {B}$ is divergence free to first order in $L_\perp /L_\parallel$. We shall assume that plasma beta is first order in $L_\perp /L_\parallel$.

Calculating $j_{||}$ by taking the curl of (2.12a,b) and substituting in (2.11) then leads to the following elliptic partial differential equation (PDE)  for $A$:

(2.13a,b)\begin{equation} -\!\varDelta^* A + F(\psi) \frac{p'(\psi)}{B^2} +H'(\psi)=0, \quad \varDelta^* A= \frac{1}{B_v^2}\boldsymbol{\nabla}\boldsymbol{\cdot} \left( B_v^2 \boldsymbol{\nabla} A \right). \end{equation}

Thus, (2.11) is equivalent to an elliptic QS-GSE, which reduces to the standard GS in the axisymmetric limit, as shown in Appendix B. This follows immediately from $\boldsymbol {B}_v \propto \boldsymbol {\nabla } \phi, B_v^2 \propto 1/R^2, A\propto -\psi$ in standard axisymmetric cylindrical coordinates. Thus, the $(F/B^2) p', H'$ terms reduce to the standard $R^2 p', I I'$ terms in the axisymmetric Grad–Shafranov equation.

3 The quasisymmetric high-beta stellarator model

In the previous section, we have shown that perfect QS helps to integrate out one hyperbolic characteristic of ideal MHD (Grad Reference Grad1971). Therefore, the Hamada condition (Helander Reference Helander2014) is automatically satisfied, which leads to a Pfirsch–Schluter current that is non-singular on rational flux surfaces. However, there remains another hyperbolic characteristic of ideal MHD that is related to $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\psi =0$. In the absence of a continuous symmetry, nested flux surfaces do not exist, and islands and chaotic fields are formed generically throughout the plasma volume.

In this section, we shall primarily focus on the HBS model by choosing the lowest-order vacuum field to be purely toroidal, which leads to an analytically tractable model. We shall show that the assumption of QS also helps support nested flux surfaces.

3.1 Derivation of the QS-HBS model

Henceforth, we assume that the aspect ratio is large and that the leading-order magnetic field is a purely toroidal vacuum field. It is then convenient to use the standard right-handed cylindrical coordinates $(R, Z, \phi )$, with unit vectors $(\boldsymbol {e_R},\boldsymbol {e}_{\boldsymbol {Z}},\boldsymbol {e}_{\boldsymbol {\phi }})$

(3.1a-c)\begin{equation} \boldsymbol{r}=R\boldsymbol{e_R}+ Z \boldsymbol{e}_{\boldsymbol{Z}}, \quad \boldsymbol{e}_{\boldsymbol{\phi}}=\boldsymbol{e_R}\times \boldsymbol{e}_{\boldsymbol{Z}}, \quad \boldsymbol{\nabla} \phi =\frac{1}{R}\boldsymbol{e}_{\boldsymbol{\phi}}. \end{equation}

The inverse aspect ratio $\epsilon$, defined as the ratio of the minor radius $r_0$, and the major radius $R_0$, is our expansion parameter, i.e.

(3.2)\begin{equation} \epsilon\equiv \frac{r_0}{R_0} \ll 1. \end{equation}

However, assuming a purely toroidal vacuum field to lowest order implies that the magnetic axis will remain close to a planar ring whose normal does not rotate around itself. Therefore, we can only access QA for which $N=0$ (Landreman & Sengupta Reference Landreman and Sengupta2019; Rodríguez et al. Reference Rodríguez, Sengupta and Bhattacharjee2023) and

(3.3)\begin{equation} F(\psi)= \frac{G(\psi)}{\iota(\psi)} \quad \text{(quasiaxisymmetry)}. \end{equation}

Following Freidberg (Reference Freidberg2014), we now expand the various physical quantities in formal power series of $\epsilon$,

(3.4a-d)\begin{align} \boldsymbol{B} =\boldsymbol{B}_0+\epsilon \boldsymbol{B}_1 +\dots, \quad B =B_0 +\epsilon B_1 +\dots,\quad p=\epsilon p_1 +\dots \quad \boldsymbol{J} = \epsilon \boldsymbol{J}_1+\dots, \end{align}

assuming $B_0$ to be a constant, and the plasma beta and currents are first order in $\epsilon$. All quantities such as $\boldsymbol {B}_1,B_1,p_1$ etc. are assumed to be $O(1)$. It is convenient to normalize all length scales with the minor radius $r_0$. Thus

(3.5a-c)\begin{equation} R=\frac{1}{\epsilon}+x ,\quad Z=y, \quad \boldsymbol{\nabla} \phi = \boldsymbol{e}_{\boldsymbol{\phi}}\,\epsilon(1+\epsilon x)^{{-}1}, \end{equation}

where $(x,y)$ are order-unity coordinates along $\boldsymbol {e_R}$ and $\boldsymbol {e}_{\boldsymbol {Z}}$. Since the lowest-order magnetic field is assumed to be a toroidal vacuum field, $G\sim R_0 B_0$, and

(3.6a,b)\begin{equation} \boldsymbol{B} = \frac{1}{\epsilon}G_{{-}1} \boldsymbol{\nabla} \phi + O(\epsilon)= B_0 \boldsymbol{e}_{\boldsymbol{\phi}}+ O(\epsilon), \quad G_{{-}1}=B_0. \end{equation}

The magnetic field is thus divergence free and curl free to $O(\epsilon )$. Moreover, $\boldsymbol {B}_0\boldsymbol {\cdot } \boldsymbol {\nabla }= O(\epsilon )$ from (3.5a–d) as required from $L_\perp /L_\parallel \sim \epsilon$. Thus, the (normalized) gradients naturally split into

(3.7a-c)\begin{equation} \boldsymbol{\nabla} =\boldsymbol{\nabla}_{\boldsymbol{\perp}} + \epsilon \boldsymbol{\nabla}_\phi,\quad \boldsymbol{\nabla}_{\boldsymbol{\perp}} \equiv \boldsymbol{e}_{\boldsymbol{R}}\partial_x + \boldsymbol{e}_{\boldsymbol{Z}}\partial_y, \quad \boldsymbol{\nabla}_\phi \equiv \boldsymbol{e}_{\boldsymbol{\phi}}(1+\epsilon x)^{{-}1} \partial_\phi. \end{equation}

This completes the description of the lowest-order magnetic field and its associated coordinate system. Next, we analyse the ideal MHD system (2.1) and the QS condition (2.3), order by order. Since $B=B_0+O(\epsilon )$ and $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\sim \epsilon$, the QS condition only imposes a constraint at $O(\epsilon ^2)$ and higher. Except for the treatment of QS, all other details of the expansion are given in Freidberg (Reference Freidberg2014), so we only include the essential results here. We note that one must be careful with the operator $\boldsymbol {\nabla }_\phi$ as defined in (3.7a–c) since it generates $x$ dependent terms in higher orders to account for the $1/R$ term in $\boldsymbol {\nabla }\phi$.

Proceeding to $O(\epsilon )$, we find

(3.8a)\begin{gather} \boldsymbol{\nabla}_{\boldsymbol{\perp}} \boldsymbol{\cdot} \boldsymbol{B}_1=0 , \end{gather}

(3.8b)\begin{gather}\boldsymbol{J}_1\times \boldsymbol{B}_0 = \boldsymbol{\nabla}_{\boldsymbol{\perp}} p_1, \end{gather}

(3.8c)\begin{gather}\boldsymbol{J}_1=\left( \boldsymbol{\nabla} \times \boldsymbol{B}\right)_1. \end{gather}

It can be shown Freidberg (Reference Freidberg2014) that the divergence-free condition, (3.8a), leads to

(3.9)\begin{equation} \frac{\boldsymbol{B}_1}{B_0} ={-}\left( x +b_{\phi 1}\right)\boldsymbol{e}_{\boldsymbol{\phi}} + \boldsymbol{\nabla} \left(\frac{A_1}{B_0}\right) \times \boldsymbol{e}_{\boldsymbol{\phi}}, \end{equation}

where the $A_1$ term defines a poloidal magnetic field with $A_1$ being the streamfunction. The $b_{\phi 1}$ term on the right of the above equation, which originates from the pressure-induced corrections to the $1/R$ vacuum field, can be determined in terms of plasma beta using the force-balance condition (3.8b) and the definition of current ${J}_1$ (3.8c)

(3.10)\begin{equation} b_{\phi 1}(\psi)= \frac{p_1(\psi)}{B_0^2}. \end{equation}

Taking into account the orthogonality of the two terms in (3.9), it follows that the field strength $B\approx B_0+\epsilon B_1$, where

(3.11)\begin{equation} \frac{B_1}{B_0}={-} \left( x + b_{\phi 1}(\psi) \right). \end{equation}

Finally, taking the curl of $\boldsymbol {B}$ to first order, we obtain the current

(3.12a,b)\begin{equation} \frac{\boldsymbol{J}_1}{B_0}={-}\boldsymbol{e}_{\boldsymbol{\phi}} \nabla_\perp^2 \left(\frac{A_1}{B_0}\right)-\boldsymbol{\nabla}_{\boldsymbol{\perp}} b_{\phi 1}\times \boldsymbol{e}_{\boldsymbol{\phi}}, \quad \nabla_\perp^2 \equiv \partial_x^2 +\partial_y^2. \end{equation}

As discussed in the previous section, the parallel component of $\boldsymbol {J}_1$ can be determined through a magnetic differential equation of the form (2.8). This requires us to go to second order in $\epsilon$.

From the second order, we find the two fundamental HBS equations

(3.13a)\begin{gather} \boldsymbol{\nabla}_{||} \psi=0 , \end{gather}

(3.13b)\begin{gather}\boldsymbol{\nabla}_{||} \left( \nabla_\perp^2 \left( \frac{A_1}{B_0}\right)\right)={-}2 \frac{{\rm d} b_{\phi 1}}{{\rm d}\psi}\partial_y \psi , \end{gather}

where we have used

(3.14a,b)\begin{equation} \boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}=\epsilon B_0 \boldsymbol{\nabla}_{||}, \quad \boldsymbol{\nabla}_{||}\equiv \partial_\phi -\left\{\left(\frac{A_1}{B_0} \right), \ \right\}_{(x,y)}. \end{equation}

The Poisson bracket appearing in (3.14a,b) is defined in the usual way as

(3.15)\begin{equation} \left\{\mathcal{A}_1, \ \right\}_{(x,y)}= \boldsymbol{e}_{\boldsymbol{\phi}}\times \boldsymbol{\nabla}_{\boldsymbol{\perp}} \mathcal{A}_1 \boldsymbol{\cdot} \boldsymbol{\nabla}_{\boldsymbol{\perp}}. \end{equation}

The system (3.13) is much simpler than the full ideal MHD system but is still highly nonlinear.

Now, let us turn to the QS condition (2.3). We can obtain $G(\psi )$ and $F(\psi )$ from their definitions (2.5a,b) and (3.3) as

(3.16a,b)\begin{equation} \frac{G(\psi)}{B_0}= \oint \frac{{\rm d}\phi}{2{\rm \pi}}\,R\frac{\boldsymbol{B}}{B_0} \boldsymbol{\cdot} \boldsymbol{e}_{\boldsymbol{\phi}} \approx \frac{1}{\epsilon} -b_{\phi 1}, \quad \frac{F(\psi)}{B_0}= \frac{1}{\epsilon\, \iota(\psi)} \left( 1 -\epsilon b_{\phi 1}\right) . \end{equation}

The QS condition (2.3) then implies that

(3.17)\begin{gather} \boldsymbol{e}_{\boldsymbol{\phi}}\times \boldsymbol{\nabla}_{\boldsymbol{\perp}} \psi_F \boldsymbol{\cdot} \boldsymbol{\nabla}_{\boldsymbol{\perp}} B_1 = B_0\boldsymbol{\nabla}_{||} B_1, \end{gather}

(3.18)\begin{gather}\psi_F=\int {\rm d}\psi \,\frac{B_0}{\epsilon F(\psi)}= \int {\rm d}\psi\,\iota(\psi)(1+\epsilon b_{\phi_1}). \end{gather}

It follows from the expression (3.16a,b) of $F(\psi )$ that, to the lowest order in $\epsilon$, $\psi _F$ is equal to the poloidal flux

(3.19)\begin{equation} \psi_F \approx \psi_p\equiv \int \iota(\psi) \,{\rm d}\psi. \end{equation}

Now, substituting $B_1$ from (3.11) into (3.17), and using (3.14a,b) and (3.15) we can rewrite the QS constraint as

(3.20)\begin{equation} \partial_y \left( \psi_F+A_1 \right) =0. \end{equation}

Lifting the partial $y$ derivative in (3.20), we obtain

(3.21)\begin{equation} \psi_F + A_1= B_0 \mathfrak{a}(x,\phi) . \end{equation}

Equations (3.13) and (3.21) constitute our quasisymmetric HBS system. To proceed further, it is convenient to introduce normalized variables $(\beta, \varPsi,\mathcal {A})$ such that

(3.22a-c)\begin{equation} \varPsi=\frac{\psi_F}{B_0} ,\quad \mathcal{A}=\frac{ A_1}{ B_0} , \quad \beta=\frac{p_1}{B_0^2} . \end{equation}

In terms of the normalized variables, together with the definitions

(3.23a,b)\begin{equation} \boldsymbol{\nabla}_{||} \equiv \partial_\phi -\{\mathcal{A},\ \}_{(x,y)}, \quad \nabla_\perp^2 \equiv \partial_x^2 +\partial_y^2, \end{equation}

the QS-HBS model reads

(3.24a)\begin{gather} \boldsymbol{\nabla}_{||} \varPsi =0, \quad \text{(nested flux-surface condition)} \end{gather}

(3.24b)\begin{gather}\boldsymbol{\nabla}_{||} \nabla_\perp^2 \mathcal{A} ={-}2 \frac{{\rm d}\beta}{{\rm d}\varPsi}\partial_y \varPsi, \quad \text{(parallel current condition)} \end{gather}

(3.24c)\begin{gather}\varPsi +\mathcal{A} =\mathfrak{a}(x,\phi). \quad \text{(relation of $\varPsi$ and $A$ from QS)} \end{gather}

Once these equations are solved, we can calculate $\boldsymbol {B},\boldsymbol {J}$ to $O(\epsilon )$ through

(3.25a)\begin{gather} \frac{\boldsymbol{B}}{B_0} = \boldsymbol{e}_{\boldsymbol{\phi}} +\epsilon \left( -(x+\beta)\boldsymbol{e}_{\boldsymbol{\phi}}+\boldsymbol{\nabla}_{\boldsymbol{\perp}} \mathcal{A} \times \boldsymbol{e}_{\boldsymbol{\phi}} \right), \end{gather}

(3.25b)\begin{gather}\frac{B}{B_0} =1-\epsilon \left( x + \beta \right), \end{gather}

(3.25c)\begin{gather}\frac{\boldsymbol{J}}{B_0} =\epsilon \left( -\boldsymbol{e}_{\boldsymbol{\phi}} \nabla_\perp^2 \mathcal{A} - \boldsymbol{\nabla}_{\boldsymbol{\perp}} \beta \times \boldsymbol{e}_{\boldsymbol{\phi}}\right). \end{gather}

This completes our derivation of the basic equations that govern near-axisymmetric quasisymmetric reduced MHD. These are nonlinear equations for $(\varPsi,\mathcal {A})$, with single valuedness of $(\varPsi,\mathcal {A})$ as boundary conditions in the $(x,y,\phi )$ coordinates. The pressure term $\beta =p_1/B_0^2$ is a free input function. The function $\mathfrak {a}(x,\phi )$ is also an input function, but there are some constraints on it that we will discuss in § 3.3 and Appendix A.

In the next section, we proceed with the QS-HBS system (3.24) and obtain a QS-GSE for the function $\varPsi$.

3.2 Derivation of the quasisymmetric Grad–Shafranov equation

The QS-GSE was obtained in full generality in Burby et al. (Reference Burby, Kallinikos and MacKay2020). However, it is hard to use and somewhat opaque due to the complexities arising from the 3-D nature of the metric coefficients and the additional consistency conditions. The LAE allows us to cut through these complications and reduce the QS-GSE to a simple form, thus highlighting the essential role of geometry and QS.

We begin by noting that the definition of $\boldsymbol {\nabla }_{||}$ (3.23a,b), and the relation between $\varPsi$ and $\mathcal {A}$ (3.24c) yields the identity

(3.26)\begin{equation} \boldsymbol{\nabla}_{||} x ={-}\partial_y \varPsi, \end{equation}

which is the large-aspect-ratio limit of the two-term QS formula (2.3). Identity (3.26) implies that the parallel current equation (3.24b) can be written as

(3.27)\begin{equation} -\nabla_\perp^2 \mathcal{A} + 2 x\frac{{\rm d}\beta}{{\rm d}\varPsi} + H'(\varPsi)=0. \end{equation}

To obtain a single nonlinear equation for $\varPsi$ we now eliminate $\mathcal {A}$ from (3.27) using (3.24c). We then obtain the QS-GSE

(3.28)\begin{equation} \nabla_\perp^2 \varPsi+ 2 x\frac{{\rm d}\beta(\varPsi)}{{\rm d}\varPsi} + H'(\varPsi)-\mathfrak{a}_{,xx}=0. \end{equation}

Here and elsewhere, we use the notation $f_{,x}$ to denote the $x$ derivative of $f$.

Next, we substitute $\mathcal {A}=\mathfrak {a} -\varPsi$ in the $\boldsymbol {\nabla }_{||} \varPsi =0$ equation to find

(3.29)\begin{equation} \varPsi_{,\phi} -\mathfrak{a}_{,x}\varPsi_{,y}=0. \end{equation}

Using the method of characteristics, we can obtain the following exact solution to (3.29) in the form of a travelling wave (TW)

(3.30a,b)\begin{equation} \varPsi=\varPsi(x,\xi), \quad \xi= y + \int {\rm d}\phi\,\mathfrak{a}_{,x}. \end{equation}

To understand the physical meaning of the quantity $\xi$, we now construct the QS vector $\boldsymbol {u}$ as defined in (2.6a,b). As shown in Appendix B

(3.31a,b)\begin{equation} \boldsymbol{u}= \frac{1}{\epsilon\, \iota}\left( \boldsymbol{e}_{\boldsymbol{\phi}} (1+\epsilon x)-\epsilon \mathfrak{a}_{,x}\boldsymbol{e}_{\boldsymbol{Z}} \right), \quad \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla} = \frac{1}{\iota} \left( \partial_\phi - \mathfrak{a}_{,x}\partial_y\right) . \end{equation}

From (3.31a,b) we find that $\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla } \xi =\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla } x=0$. Thus, the symmetry vector $\boldsymbol {u}$ lies on surfaces of constant $\xi$ and $x$, thereby ensuring that $\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {\nabla } \varPsi (x,\xi )=0$. Therefore, $\xi$ denotes one of the Clebsch variables for $\boldsymbol {u}$, the other being $x$. Moreover, since the symmetry lines must close on themselves on a torus, $\xi$ must be periodic in $\phi$.

We now look into the relation of the QS vector, $\boldsymbol {u}$, with the killing vectors of 3-D Euclidean space that generate rotations and translations (Burby et al. Reference Burby, Kallinikos and MacKay2020). Following Burby et al. (Reference Burby, Kallinikos and MacKay2020), we define

(3.32)\begin{equation} \boldsymbol{u}_{\text{HBS}}\equiv \iota \boldsymbol{u}, \end{equation}

which is equivalent to choosing the poloidal flux instead of toroidal flux in the expression for $\boldsymbol {u}$ (2.6a,b). Comparing the symmetry vectors corresponding to axisymmetry, helical symmetry and the QA symmetry in the HBS

(3.33a-c)\begin{equation} \boldsymbol{u}_{\text{AS}}= R\boldsymbol{e}_{\boldsymbol{\phi}}, \quad \boldsymbol{u}_{\text{HS}}= R\boldsymbol{e}_{\boldsymbol{\phi}}-l\,\boldsymbol{e}_{\boldsymbol{Z}}, \quad \boldsymbol{u}_{\text{HBS}}\approx R\boldsymbol{e}_{\boldsymbol{\phi}} - \mathfrak{a}_{,x} \boldsymbol{e}_{\boldsymbol{Z}}, \end{equation}

we find that they have the form of a linear combination of a rotation in $\phi$ and a translation in $Z$. We note that, in a straight cylinder with helical symmetry, it is customary to use the poloidal angle $\vartheta$ to denote the angle of rotation, and $\phi$ is along $Z$. We have chosen not to do so in order to express the close relationship between these vectors.

From (3.33a–c), we observe that the translation in $Z$ is zero for axisymmetry, constant $l$ for helical symmetry, and periodic (with zero average) in HBS. As shown in Appendix B, although the axisymmetry and helical symmetry vectors are killing, the HBS symmetry vector is not killing for a generic $\phi$ dependent $\partial _x\mathfrak {a}(x,\phi )$. Thus, the QA symmetry of the HBS model is indeed a hidden symmetry and not an isometry of the 3-D Euclidean space.

3.3 Quasisymmetric Grad–Shafranov equation and consistency conditions

In the case of axisymmetry, $\mathfrak {a}=0$, which implies $\varPsi _{,\phi }=0$. The QS-GSE (3.28) then reduces to the axisymmetric GSE as expected. However, in the non-symmetric case with $\mathfrak {a}\neq 0$, it is not obvious that the QS-GSE equation (3.28), and the $\boldsymbol {\nabla }_{||}\varPsi =0$ condition (3.29), are consistent. The conflict lies in the fact that the non-symmetric terms can enter $\varPsi$ only through $\mathfrak {a}_{,x}$ in the form of a TW given in (3.30a,b). It is not immediately obvious that a TW solution will satisfy the QS-GSE (3.28) for a general $\mathfrak {a}(x,\phi )$.

As shown in Appendix A, a self-consistent solution is obtained if $\mathfrak {a}(x,\phi )$ is of the form

(3.34)\begin{equation} \mathfrak{a}_{,xx}=0 \quad \Rightarrow \quad \mathfrak{a}(x,\phi)= \bar{\mathfrak{a}}(\phi)+ x Y'(\phi), \end{equation}

where $\bar {\mathfrak {a}}(\phi ),Y(\phi )$ are periodic functions of $\phi$. Note that, in axisymmetry $(\partial _\phi =0)$, $\mathfrak {a}$ must be of the form

(3.35)\begin{equation} \mathfrak{a}= a_0 + a_1 x, \end{equation}

where $a_0,a_1$ are constants. The constants can then be absorbed through a simple redefinition of the current and pressure profiles in (3.28). Thus, $\mathfrak {a}$ can be set to zero in the axisymmetric limit with no loss of generality.

The rather strong restriction on $\mathfrak {a}(x,\phi )$, given by (3.34), implies that we need to solve the following equations for $\varPsi$, subject to single valuedness of $\varPsi$ and its derivatives:

(3.36a)\begin{gather} \nabla_\perp^2 \varPsi+ 2 x\frac{{\rm d}\beta(\varPsi)}{{\rm d}\varPsi} + H'(\varPsi)=0 , \end{gather}

(3.36b)\begin{gather}\varPsi=\varPsi(x,\xi), \quad \xi= y + Y(\phi). \end{gather}

Thus, we have a tokamak-like GSE subject to a TW deformation. The profile functions $\beta (\varPsi ), H(\varPsi )$ are related to plasma pressure and currents. The solution strategy is simply solving the GSE equation in $(x,y)$ space and then performing the shift $y\to y+Y(\phi )=\xi$.

Alternatively, one can shift from $(x,y,\phi )$ coordinates to $(x,\xi,\phi )$ coordinates, such that

(3.37a-c)\begin{equation} \boldsymbol{\nabla}_{||} = \left.\partial_\phi\right|_{(x,\xi)} +\{\varPsi,\ \}_{(x,\xi)}, \quad \partial_y \to \partial_\xi, \quad \nabla_\perp^2\to \bar{\nabla}^2_\perp{\equiv} \partial_x^2+\partial_\xi^2. \end{equation}

The details of the transformation are provided in Appendix C. Since the QS-GSE equation (3.36a) only contains $(x,y)$ derivatives and $\mathfrak {a}_{,x}$ is only a function of $\phi$, the transformed QS-GSE reads

(3.38)\begin{equation} \bar{\nabla}^2_\perp \varPsi+ 2 x\frac{{\rm d}\beta(\varPsi)}{{\rm d}\varPsi} + H'(\varPsi)=0. \end{equation}

Let us now compare our system with similar systems discussed in Landreman & Sengupta (Reference Landreman and Sengupta2018), Rodríguez et al. (Reference Rodríguez, Sengupta and Bhattacharjee2023), Plunk & Helander (Reference Plunk and Helander2018), Plunk (Reference Plunk2020) and Burby et al. (Reference Burby, Kallinikos and MacKay2020). An important consequence of the TW nature (3.36b) of $\varPsi$ is that the quasisymmetric deformation to the axisymmetric equilibrium is a periodic displacement purely in the vertical $\boldsymbol {e}_{\boldsymbol {Z}}$ direction. Thus, the system (3.36) does not have a rotating ellipse-like solution near the axis (Landreman & Sengupta Reference Landreman and Sengupta2018). It can be shown self-consistently (Rodríguez et al. Reference Rodríguez, Sengupta and Bhattacharjee2023) that, as one approaches axisymmetry, the function $\sigma (\phi )$, which controls the rotation of the ellipse, goes to zero in agreement with our results. The non-rotating aspect is also markedly different from Plunk & Helander (Reference Plunk and Helander2018) and Plunk (Reference Plunk2020), where the quasisymmetric perturbations have a harmonic $\exp {({\rm i}N\phi )}$ phase factor (with $N$ being an integer) in both $\boldsymbol {e}_{\boldsymbol {R}}$ and $\boldsymbol {e}_{\boldsymbol {Z}}$ directions.

A detailed comparison of (3.36) with the general QS-GSE system of Burby et al. (Reference Burby, Kallinikos and MacKay2020) is carried out in Appendix B. To summarize, we find that in the large-aspect-ratio limit, we recover the generalized Grad–Shafranov equations for QS derived by Burby et al. (Reference Burby, Kallinikos and MacKay2020). In general geometry, it is not enough to solve the GSE, and Burby et al. (Reference Burby, Kallinikos and MacKay2020) had to impose several extra compatibility constraints. On the other hand, the only constraint we had to impose to ensure compatibility is $\partial ^2_x\mathfrak {a}(x,\phi )=0$. (The extra constraints presumably appear at higher order, which is beyond the scope of the current work.)

3.4 Quasisymmetric high-beta stellarator as an integrable system

A magnetic field that satisfies MHD equilibrium conditions and possesses nested pressure and magnetic flux surfaces is, in principle, integrable (in the sense described precisely in Burby, Duignan & Meiss Reference Burby, Duignan and Meiss2021). Implicitly, Hamada's condition is assumed to be satisfied by any MHD equilibrium with nested surfaces (Helander Reference Helander2014). However, the integrability of a model obtained through a formal asymptotic expansion of 3-D ideal MHD equations is not always guaranteed. In particular, Freidberg's HBS model, which consists of two nonlinear coupled PDEs for $\mathcal {A},\varPsi$, does not guarantee that Hamada's condition will be satisfied. Therefore, nestedness of flux surfaces can be assumed but not self-consistently demonstrated. Here, we show that adding the QS constraint on the HBS model reinforces integrability and allows us to construct explicit action-angle coordinates. We shall now use the $(x,\xi,\phi )$ coordinates to construct action-angle coordinates.

First, we note that the definition of $\boldsymbol {\nabla }_{||}$ in (3.37a–c) implies that the pair $(x,\xi )$ satisfies

(3.39a,b)\begin{equation} \boldsymbol{\nabla}_{||} \xi = \varPsi_{,x}, \quad \boldsymbol{\nabla}_{||} x ={-}\varPsi_{,\xi}. \end{equation}

Moreover, the 2-D nature of $\varPsi$ is explicit in the $(x,\xi,\phi )$ coordinates since

(3.40)\begin{equation} \boldsymbol{\nabla}_{||} \varPsi= {\partial_\phi}|_{(x,\xi)}\varPsi =0. \end{equation}

We can interpret the $\boldsymbol {\nabla }_{||}$ operator as a total derivative with respect to $\phi$ along the magnetic field line. Therefore, we can cast (3.39a,b) as Hamilton's equations of motion

(3.41a,b)\begin{equation} \frac{{\rm d}\xi}{{\rm d}\phi} = \varPsi_{,x}, \quad \frac{{\rm d}\kern0.7pt x}{{\rm d}\phi} ={-}\varPsi_{,\xi}, \end{equation}

where $\xi$ is the coordinate, $x$ is the conjugate momentum and $\varPsi (x,\xi )$ is the Hamiltonian. The condition (3.40) implies that the Hamiltonian is independent of time.

Next, we perform a canonical transform from $(x,\xi,\phi )$ to $(\mathcal {I},\vartheta,\phi )$ such that

(3.42a,b)\begin{equation} \frac{{\rm d}\vartheta}{{\rm d}\phi}= \varPsi_{,\mathcal{I}}, \quad \frac{{\rm d}\mathcal{I}}{{\rm d}\phi}={-} \varPsi_{\vartheta}. \end{equation}

Requiring that $\varPsi$ be independent of $\vartheta$, we arrive at the action-angle coordinates, where

(3.43a,b)\begin{equation} \mathcal{I} = \frac{1}{2{\rm \pi}}\oint_\varPsi x \,{\rm d}\xi, \quad \frac{{\rm d}\vartheta}{{\rm d}\phi}= \frac{{\rm d}\varPsi}{ {\rm d}\mathcal{I}}= \iota . \end{equation}

It then follows that, to this order of accuracy, the action $\mathcal {I}$ is nothing but the toroidal flux $\psi$, $\vartheta$ is the canonical straight-field-line poloidal angle and $\iota$ is the rotational transform The following more explicit form of $\iota$ can be derived from (3.43a,b) using $\varPsi _{,x}x_{,\varPsi }=1$:

(3.44)\begin{equation} \iota^{{-}1}= \frac{1}{2{\rm \pi}}\oint_\varPsi \frac{{\rm d}\xi}{\varPsi_{,x}}. \end{equation}

In summary, we have shown that nested flux surfaces from the axisymmetric configuration are preserved since they are merely subjected to a non-symmetric but periodic deformation implicitly through $\xi$. Moreover, the magnetic field-line dynamics of the QS-HBS system is that of an integrable Hamiltonian system. The action-angle coordinates of the QS-HBS corresponds to the usual straight-field-line coordinates. In other words, QS not only enables us to avoid singular currents near rational surfaces (as discussed in § 2) but also preserves the nested flux-surface structure.

3.5 Clebsch variables for the QS-HBS system

In the last section, we have discussed the straight-field-line coordinates, which are the action-angle coordinates for the QS-HBS system. Often, in the local analysis of plasma equilibrium, such as MHD and kinetic stability, it is useful to use Clebsch variables. In this section, we derive the expressions for the Clebsch variables, $(\varPsi,\alpha,\ell )$ where $\alpha$ is the field-line label and $\ell$ is the arclength along $\boldsymbol {B}$.

The expression for the field-line label $\alpha$ can be derived from the condition $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } \alpha =0$, which implies

(3.45)\begin{equation} \boldsymbol{\nabla}_{||} \alpha = \alpha_{,\phi} + \{\varPsi,\alpha\}_{(x,\xi)}=0 . \end{equation}

To simplify (3.45), we further change coordinates to $(\varPsi,\xi,\phi )$. We have provided the necessary details of the transformation to the straight-field-line coordinates in Appendix C. In $(\varPsi,\xi,\phi )$ coordinates, (3.45) takes the simplified form

(3.46)\begin{equation} \left( \partial_\phi + \varPsi_{,x}\partial_\xi\right) \alpha =0, \end{equation}

whose solution can be immediately obtained as

(3.47)\begin{equation} \alpha = \phi - \int_\varPsi\frac{ {\rm d}\xi}{\varPsi_{,x}} +O(\epsilon). \end{equation}

Here, we have used the fact that $\varPsi$ is independent of $\phi$ at a fixed $\xi$. We can easily check that $\alpha = \phi - \iota ^{-1}\vartheta$ from the Hamilton's equations for the field lines

(3.48)\begin{equation} \frac{{\rm d}\xi}{\varPsi_{,x}}=\frac{{\rm d}\kern0.7pt x}{-\varPsi_{,\xi}}= {\rm d}\phi = \iota^{{-}1}\,{\rm d}\vartheta. \end{equation}

However, as shown in Appendix C, $\boldsymbol {B}$ as given in (3.25a) is equivalent to

(3.49)\begin{equation} \boldsymbol{B}/B_0= \boldsymbol{\nabla} \alpha \times \boldsymbol{\nabla} \varPsi, \end{equation}

only if we keep the $O(\epsilon )$ correction to (3.47). The correction is needed since the $B_1/B_0$ terms do not contribute to the $\boldsymbol {\nabla }_{||}$ operator but to $\boldsymbol {B}_1/B_0$. Thus, we choose the following form of $\alpha$:

(3.50)\begin{equation} \alpha = \phi - \int_\varPsi {\rm d}\xi\,\frac{B_\phi/B_0}{\varPsi_{,x}} = \phi - \int_\varPsi {\rm d}\xi\, \frac{1-\epsilon(x+\beta)}{\varPsi_{,x}}. \end{equation}

We note that there is a sign difference between (3.49) and what is typically used in the stellarator literature (D'haeseleer et al. Reference D'haeseleer, Hitchon, Callen and Shohet2012; Helander Reference Helander2014) owing to the definition of $\alpha$ in (3.47). However, the benefit of this form of $\alpha$ and $\boldsymbol {B}$ is that it smoothly goes over to the tokamak expressions when the deviation from axisymmetry is negligible. To that end, we can recast $\alpha$ as

(3.51a-c)\begin{equation} \alpha = \phi - q(\varPsi)\int_\varPsi {\rm d}\vartheta\, b_\phi(\varPsi,\vartheta), \quad q(\varPsi)\equiv \iota^{{-}1}, \quad b_\phi = B_\phi/B_0. \end{equation}

It is important to clarify that $(B_\phi /B_0-1)$ is only a function of $x(\varPsi,\xi )$ and $\beta (\varPsi )$, and independent of $\phi$. Through the canonical transform (3.43a,b), we then obtain $b_\phi =b_\phi (\varPsi,\vartheta )$. Furthermore, since $b_\phi =1+O(\epsilon )$, $\alpha =\phi -q(\varPsi )\vartheta +O(\epsilon )$. Thus, $\vartheta$ is indeed the straight-field-line poloidal angle as discussed in § 3.4.

Similarly, the expression for $\ell$ can be derived from the condition $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\ell =B$, or equivalently in $(\varPsi,\xi,\phi )$ coordinates

(3.52)\begin{equation} \boldsymbol{\nabla}_{||} \ell =\left( \partial_\phi +\varPsi_{,x}\partial_\xi \right) \ell =\frac{1}{\epsilon}, \end{equation}

which implies that

(3.53)\begin{equation} \ell = \frac{1}{\epsilon}\left( \phi + L(\varPsi)+O(\epsilon)\right). \end{equation}

The $\epsilon ^{-1}$ factor in the expression for $\ell$ follows from the fact that we have normalized with respect to the minor radius, whereas the arclength scales with the major radius. The form of $\ell$ (3.53) also follows from the fact that in $(\varPsi,\alpha,\phi )$ coordinates $\boldsymbol {\nabla }_{||} =\partial _\phi$. The function $L(\varPsi )$ is a homogeneous solution of the operator $\boldsymbol {\nabla }_{||}$.

4 Analytic solutions of the QS-GSE: extended Soloviev profiles

We now have in our possession the QS-HBS model, which ensures both QS and force balance to first order in the LAE. As discussed in the previous section, we can start with any LAE tokamak equilibrium and deform it to obtain a quasiaxisymmetric stellarator. In this section, we show a simple analytic exact solution to the QS-HBS model that goes beyond the scope of the NAE. Using the VMEC code (Hirshman and Whitson Reference Hirshman and Whitson1983), we can verify that this solution approximates 3-D MHD equilibrium with good volumetric QA. In particular, we demonstrate that as the aspect ratio becomes large, the QS error and the differences between the numerical 3-D equilibrium from VMEC and the asymptotic analytical model tend to zero. We shall now look for a class of LAE MHD equilibria with the following profile functions $(\beta (\varPsi ),H(\varPsi ))$:

(4.1a,b)\begin{equation} \beta= \beta_0 + \beta_1 \varPsi, \quad H= H_0 + H_1 \varPsi + \frac{H_2}{2} \varPsi^2, \end{equation}

such that the QS-GSE takes the following form of a linear equation in $\varPsi$:

(4.2)\begin{equation} \bar{\nabla}^2_\perp \varPsi + 2x\beta_1 +\left( H_1 + H_2 \varPsi\right) =0. \end{equation}

Here, $\beta _i,H_j; i=0,1; j=0,1,2$ are constants that we can freely choose. If the quadratic term $H_2$ is identically zero, the profiles reduce to the so-called Soloviev profiles. It is straightforward to find a solution of (4.2) for a circular cross-section device

(4.3)\begin{equation} \varPsi= \frac{H_1}{H_2 {\rm J}_0(\sqrt{H_2})}{\rm J}_0(\sqrt{H_2}r) - \frac{H_1}{H_2} + \frac{2\beta_1}{H_2 {\rm J}_1(\sqrt{H_2})}{\rm J}_1(\sqrt{H_2}r)\cos\theta - \frac{2\beta_1 r}{H_2}\cos\theta. \end{equation}

Here, ${\rm J}_n$ are order-$n$ Bessel functions of the first kind. We combine this solution with a deformation $Y(\phi ) = -0.5\sin {2\phi }$, where each poloidal plane is rigidly displaced in the vertical direction by $-Y(\phi )$.

We then use the above $\varPsi$ to find the rotational transform numerically by taking the derivative of poloidal flux with respect to toroidal flux and passing the profile to VMEC. This was repeated for four different aspect ratios, namely 5, 20/3, 10 and 20, and $H_1$ was scaled by the inverse aspect ratio as $H_1 = 0.04/\epsilon$. The values of the other parameters were held constant at $H_2 = 4.8$ and $\beta _1 = 16{\rm \pi} \times 10^{-2}$. Note that $\beta (\varPsi )$ here is normalized to be $O(1)$; the true ratio of hydrodynamic to magnetic pressure is $\epsilon \beta (\varPsi )$. To convert normalized quantities to SI units, we used the normalization factors $B_0 = 1\ \mathrm {T}$ and $r_0 = 1\ \mathrm {m}$. This choice of parameters corresponds to holding the parameters $r_0$, $B_0$, $p_1$, $H_{1,\mathrm {SI}}$ and $H_{2,\mathrm {SI}}$ constant when the QS-GSE is written in the SI system

(4.4)\begin{equation} \nabla_\perp^2\varPsi_\mathrm{SI} + \frac{2\mu_0 p_1}{R_0 B_0^2}x + \frac{H_{1,\mathrm{SI}}}{B_0} + \frac{H_{2,\mathrm{SI}}}{B_0}\varPsi_\mathrm{SI} = 0, \end{equation}

where $p_\mathrm {SI}(\varPsi _\mathrm {SI}) = p_0 + p_1\varPsi _\mathrm {SI}$ and $H_\mathrm {SI}'(\varPsi _\mathrm {SI}) = H_{1,\mathrm {SI}} + H_{2,\mathrm {SI}}\varPsi _\mathrm {SI} = H(\varPsi )$.

The outer flux surface of the analytic equilibrium is shown in figure 1, where the colors depict the field strength. The flux surfaces obtained from the analytical expression and VMEC are compared in figure 2. Note that the main difference between the VMEC flux surfaces and the analytical ones is that VMEC shows a rotating ellipse effect, which cannot be captured by the analytical model, as mentioned before. Finally, the maximum QS error in the equilibrium is defined as

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

and is plotted in figure 3 as a function of inverse aspect ratio $\epsilon$. Here, $\hat {B}_{n,m}(\psi )$ is a Fourier mode of $B = |\boldsymbol {B}|$ on flux surface $\psi$. As can be seen, it scales as $\epsilon ^2$, which is expected since the QS-HBS model is derived at order $\epsilon$.

Figure 1.

The deformed equilibrium with a circular cross-section and aspect ratio 5.

Figure 2.

Approximate analytical flux surfaces (solid blue) vs numerical flux surfaces from VMEC (dashed black) at the $\phi = 0$ poloidal plane for aspect ratios 5 (a) and 10 (b).

Figure 3.

The maximum QS error (black dots) scales as $\epsilon ^2$ (dashed blue line).

5 Numerical solutions: more complex geometry

Since (3.36a) matches the large-aspect-ratio limit of the standard axisymmetric GSE, one can numerically solve the QS-HBS model by taking any numerical tokamak equilibrium with a large aspect ratio and applying a periodic vertical deformation as given by (3.36b).

With this consideration, we use VMEC to calculate an axisymmetric equilibrium with an ITER-like cross-section at aspect ratios 4.84, 6.44, 9.65 and 28.87. The toroidal current was preserved across all aspect ratios, with the $\iota$ profile varying significantly. The aspect ratios were chosen to be close to those in the previous section except for the last one, which had to be significantly larger as the $\iota$ profile would cross 2 at aspect ratios around 20 with the chosen current profile. The presence of a low-order rational surface leads to large numerical errors in VMEC, which results in a large QS error in the deformed equilibria. In the next step, a deformation $Y(\phi ) = -0.5b\sin {2\phi }$ is applied, where $b$ is the distance from the axis to the upper tip of the outermost flux surface, and the equilibria are recomputed with VMEC. The lowest-aspect-ratio deformed equilibrium is shown in figure 4. Figure 5 shows the flux surfaces from the axisymmetric and deformed equilibria plotted on top of each other. As before, the deformed equilibria show a rotating ellipse effect, which cannot be accounted for by simply displacing the axisymmetric poloidal planes vertically by $-Y(\phi )$. The maximum QS error is plotted in figure 6 and again scales as $\epsilon ^2$, but is larger by a factor of ${\sim }2.3$ than in the circular cross-section case.

Figure 4.

The deformed equilibrium with an ITER-like cross-section and aspect ratio 4.84.

Figure 5.

Flux surfaces from the axisymmetric equilibrium (solid blue) vs from the deformed equilibrium (dashed black) at the $\phi = 0$ poloidal plane for aspect ratios 5 (a) and 10 (b).

Figure 6.

The maximum QS error (black dots) again scales as $\epsilon ^2$ (dashed blue line is $2.3\epsilon ^2$).

6 Ballooning and interchange stability

In the previous sections, we developed a model for a large-aspect-ratio stellarator with approximate volumetric quasi-axisymmetry close to axisymmetry. We have also demonstrated good agreement of our model with actual 3-D equilibria computed using the VMEC code for an aspect ratio of five and higher. These equilibria can handle plasma $\beta$ of $O(\epsilon )$, consistent with the high-$\beta$ stellarator model. Since high-$\beta$ equilibria are particularly sensitive to ideal interchange and ballooning instabilities, an MHD stability analysis is crucial. We note here that our stability analysis relies on standard tools such as VMEC. However, it is well known that VMEC can not accurately resolve current singularities. How such singularities can affect ideal MHD stability in the present model is a question we defer to the future.

In this section, we perform a stability analysis of the QS-HBS equilibria obtained here with respect to interchange and ideal ballooning modes. As is well known, interchange and ideal ballooning modes arise due to the existence of regions of unfavourable curvature, i.e. $(\boldsymbol {B}\boldsymbol {\cdot } \boldsymbol {\nabla }\boldsymbol {B})\boldsymbol {\cdot }\boldsymbol {\nabla }p > 0$, which tends to destabilize the plasma.

First, we calculate the Mercier (or interchange) stability quantified by the quantity $D_\mathrm {Merc}$ from VMEC. The Mercier stability criterion is

(6.1)\begin{equation} D_\mathrm{Merc} ={-}\frac{p^{\prime}}{\iota^{\prime}}V^{{{\dagger}}{{\dagger}}} - \frac{1}{4} > 0, \end{equation}

where the quantity $V^{{\dagger} {\dagger} }$ is related to the magnetic well (Greene Reference Greene1997). The expression for $V^{{\dagger} {\dagger} }$ contains the contributions from different terms that arise from the magnetic well, geodesic curvature, plasma current and magnetic shear. A positive $D_\mathrm {Merc}$ corresponds to stability against interchange modes, whereas a negative $D_\mathrm {Merc}$ implies instability.

The radial variation of $D_\mathrm {Merc}$ for the circular and ITER-like equilibria is plotted in figure 7. Interchange modes arise in regions of low magnetic shear. For these equilibria, the magnetic shear is the lowest near the magnetic axis, which causes the circular equilibria to become Mercier unstable. Also, we note that for the ITER-like equilibria, the Mercier stability plots do not change significantly after we impose the toroidal-symmetry-breaking perturbation.

Figure 7.

Mercier stability plots for the circular equilibria (a) and the ITER-like equilibria (b). The dashed curves are for the perturbed equilibria, whereas solid curves correspond to unperturbed tokamak equilibria. The red and green curves correspond to equilibria with an aspect ratio of five, and the blue and black curves correspond to an aspect ratio of ten.

As the equilibria analysed are stable against interchange modes for most of the volume, we then analyse their stability against the infinite-$n$, ideal ballooning mode. The ideal ballooning mode is another curvature-driven instability that can appear in equilibria with finite magnetic shear. To calculate the stability, we numerically solve (Gaur et al. Reference Gaur, Buller, Ruth, Landreman, Abel and Dorland2023) the ballooning eigenmode equation (Connor, Hastie & Taylor Reference Connor, Hastie and Taylor1978; Dewar & Glasser Reference Dewar and Glasser1983)

(6.2)\begin{equation} \frac{{\rm d}}{{\rm d}\theta} g \frac{{\rm d}X}{{\rm d}\theta} + c X = \lambda f X, \end{equation}

where $\lambda$ is the ballooning eigenvalue, and $g,c,f$ are the following normalized geometric coefficients:

(6.3)\begin{equation} \left.\begin{gathered} g = (\boldsymbol{b}\boldsymbol{\cdot} \boldsymbol{\nabla}\theta) \frac{\lvert \boldsymbol{\nabla}\alpha\rvert^2}{B}\\ c = \frac{1}{B^2}\frac{{\rm d}(\mu_0 p)}{{\rm d}\psi} \frac{2}{(\boldsymbol{b} \boldsymbol{\cdot}\boldsymbol{\nabla}\theta)} (\boldsymbol{b}\times (\boldsymbol{b} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{b})) \times \boldsymbol{\nabla}\alpha \\ f = \frac{1}{(\boldsymbol{b}\boldsymbol{\cdot}\boldsymbol{\nabla}\theta)} \frac{\lvert \boldsymbol{\nabla}\alpha\rvert^2}{B^3}, \end{gathered}\right\} \end{equation}

where $\alpha = \phi - 1/\iota (\theta - \theta _0)$ is the field-line label, $\theta$ is the PEST straight-field-line angle, $\phi$ is the cylindrical toroidal angle and $\theta _0$ is the ballooning parameter. All lengths in the ballooning equation are normalized using $a_{{N}}$, the effective minor radius (named Aminor_p) in the VMEC output, and magnetic field and plasma pressure are normalized using $B_{{N}} = 2 \psi /a_{{N}}^2$, where $\psi$ is the toroidal flux enclosed by the boundary (without the factor of $2 {\rm \pi}$). We solve (6.2) for the eigenvalue $\lambda = \gamma ^2 a_{{N}}^2/B_{{N}}^2$, where $\gamma a_{{N}}/B_{{N}}$ is the normalized growth rate, and the ballooning eigenfunction $X$. If $\lambda > 0$, an equilibrium is unstable against the ballooning mode; otherwise, it is stable.

For tokamak geometry, we solve (6.2) on multiple flux surfaces. For each flux surface, we scan $n_{\theta _{0}} = 96$ uniformly spaced values of $\theta _0 \in [-{\rm \pi} /2, {\rm \pi}/2]$ and $n_{\alpha } = 1$ on the field line $\alpha = 0$. For perturbed tokamaks, we repeat the same process on $n_{\alpha } = 96$ field lines with uniformly spaced values of $\alpha \in [-{\rm \pi}, {\rm \pi}]$. A typical plot of the growth rate for the perturbed ITER-like case is shown in figure 8.

Figure 8.

Ballooning eigenvalue ($\lambda$) contour plots at $\rho = 0.2$ and $\rho = 0.93$ for the perturbed aspect ratio five ITER equilibrium. The difference in the size of the ballooning unstable regions can be attributed to a large local magnetic shear in the outer region of the stellarator case. Hence, we use $n_{\alpha } = 96, n_{\theta _0} = 96$ points to accurately calculate the maximum growth rate. This process is repeated for all flux surfaces.

In figure 9, we present the ideal ballooning stability eigenvalue $\lambda$ as a function of the radial coordinate $s$. For the large aspect ratio tokamak, we observe the emergence of ballooning instability after perturbing it into a stellarator. This is not surprising since a lack of shaping will, in general, destabilize the ballooning mode. We repeat this process for the perturbed ITER-like equilibrium and plot the ballooning eigenvalue in figure 10. The tokamak equilibria are ballooning stable, whereas the deformed tokamak becomes ballooning unstable near the axis for both aspect ratios and near the edge for the aspect ratio five. However, because of the strong shaping, the ITER-like equilibria have significantly better ballooning stability properties after being perturbed compared with the circular tokamak case. Thus, strongly shaped QS-HBS equilibria may retain favourable stability even after the toroidal symmetry is broken.

Figure 9.

Ideal ballooning stability plots for the circular equilibria with an aspect ratio of ten. The dashed curves are for the perturbed equilibria, whereas solid curves correspond to unperturbed tokamak equilibria. The pressure gradient for the perturbed equilibrium is finite at the edge, which results in a finite edge growth rate.

Figure 10.

Ideal ballooning stability plots for the ITER-like equilibria with an aspect ratio of five (a) and ten (b). The solid curves correspond to the perturbed equilibria, whereas dashed curves correspond to unperturbed tokamak equilibria.

To gain insight into the ballooning stability of the ITER-like equilibria, we calculate the effective ballooning by applying the transformation $\hat {X} = \sqrt {g} X$ to (6.2). The transformed ballooning equation becomes

(6.4a,b)\begin{equation} \frac{{\rm d}^2 \hat{X}}{{\rm d}\theta^2} + \left( E-V_{\mathrm{ball}}\right) \hat{X} = 0, \quad E ={-}\frac{\lambda f}{g}, \end{equation}

where

(6.5a,b)\begin{equation} -V_{\mathrm{ball}} = \frac{c}{g} - \frac{{\rm d}}{{\rm d}\theta}\left(\frac{a}{2}\right) - \frac{a^2}{4}, \quad a = \frac{{\rm d}\log(g)}{{\rm d}\theta} \end{equation}

is called the effective ballooning potential. We plot the effective ballooning potential at the maximum growth rate values for the four growth rate plots from figure 10 in figure 11 together with the normalized eigenfunction $\hat {X}$. We can see that the stellarator equilibrium potentials have higher peaks compared with tokamaks, which cause the ballooning eigenfunction to decay rapidly, similar to the decay of a wavefunction in a potential well in Schrödinger's equation. Towards the edge, we also observe more oscillations in stellarator potentials as seen in figure 11(a) due to their 3-D shape. This decay of the eigenfunction is similar to the phenomenon of Anderson localization (Anderson Reference Anderson1958). Multiple studies by Redi et al. (Reference Redi, Canik, Dewar, Fredrickson, Cooper, Johnson and Klasky2001, Reference Redi, Johnson, Klasky, Canik, Dewar and Cooper2002) have demonstrated the importance of Anderson localization of ballooning modes due to perturbations that break axisymmetry of the background equilibrium. Further detailed analysis of the localization of the 3-D modes will be published elsewhere.

Figure 11.

Eigenfunction $\hat {X}$ and effective ballooning potential at the maximum growth rate for the original and perturbed ITER equilibria. In each panel, the eigenfunctions have been scaled by a factor of the maximum effective ballooning potential. Panels (a,c) correspond to the aspect ratio of five stellarator and ITER cases, respectively. Similarly, (b,d) correspond to the aspect ratio of ten stellarator and ITER cases, respectively.