Asymptotic quasisymmetric high-beta three-dimensional magnetohydrodynamic equilibria near axisymmetry

Abstract

Quasisymmetry (QS), a hidden symmetry of the magnetic field strength, is known to support nested flux surfaces and provide superior particle confinement in stellarators. In this work, we study the ideal magnetohydrodynamic (MHD) equilibrium and stability of high-beta plasma in a large-aspect-ratio stellarator. In particular, we show that the lowest-order description of a near-axisymmetric equilibrium vastly simplifies the problem of three-dimensional quasisymmetric MHD equilibria, which can be reduced to a standard elliptic Grad–Shafranov equation for the flux function. We show that any large-aspect-ratio tokamak, deformed periodically in the vertical direction, is a stellarator with approximate volumetric QS. We discuss exact analytical solutions and numerical benchmarks. Finally, we discuss the ideal ballooning and interchange stability of some of our equilibrium configurations.

1 Introduction

An essential requirement of magnetic confinement is that the plasma must be in a stable magnetohydrodynamic (MHD) equilibrium. Since the confining magnetic field in a stellarator is mainly produced by external current-carrying coils or magnets, it does not suffer from the current-driven instabilities intrinsic to tokamaks. For a stellarator, a set of nested toroidal flux surfaces with a common magnetic axis such that the magnetic field lines are tangential to the surfaces is desired. However, the lack of continuous symmetry in a stellarator can lead to a breakdown of nested flux surfaces and the formation of island chains and ergodic field-line regions. Without symmetry, the trapped particles can drift out radially and be lost. Both of these effects are severely detrimental to confinement. Furthermore, the lack of symmetry and the nonlinear nature of ideal MHD equations seriously complicate our understanding of stellarator equilibrium. Consequently, compared with tokamaks, MHD equilibrium theory for generic stellarators still has major unanswered questions.

A hidden symmetry of the magnetic field strength, called quasisymmetry (QS), can ameliorate some of the difficulties mentioned above. While the magnetic field remains fully three-dimensional (3-D), QS requires that the field strength, $B$, be independent of one of the angular coordinates. The continuous symmetry of $B$ then leads to a conserved canonical angular momentum (due to Noether's theorem), ensuring neoclassical particle confinement of the same quality as that of a tokamak (Helander 2014; Landreman & Catto 2010). Moreover, exact QS¹ leads to nested flux surfaces (Burby et al. 2020; Rodríguez, Helander & Bhattacharjee 2020). Indeed, one of the benefits of understanding the space of near-axisymmetric QS is that it provides a valuable perspective for 3-D error-field correction of tokamaks (Park et al. 2021), which are known to be sensitive to 3-D perturbations (Park, Boozer & Glasser 2007a; Park et al. 2009). If the 3-D perturbations of an axisymmetric tokamak restore QS, they can help mitigate harmful transport effects (Park et al. 2007b, 2018).

Unfortunately, exact QS is difficult, if not impossible, to obtain. Analytical results (Garren & Boozer 1991b, a; Plunk & Helander 2018; Landreman & Sengupta 2019; Jorge, Sengupta & Landreman 2020) seem to indicate that imposing QS on ideal MHD with scalar pressure leads in general to an overdetermined problem, with no global (or volumetric) solution. In particular, formal series expansions carried out in the distance from the axis, the so-called near-axis expansion (NAE), show that, beyond the second order, there are more equations than unknowns. A large-scale numerical optimization approach has successfully provided multiple practical designs (Anderson et al. 1995; Zarnstorff et al. 2001; Najmabadi et al. 2008; Ku & Boozer 2010; Bader et al. 2019; Landreman & Paul 2022). However, the many degrees of freedom arising from the three-dimensionality of stellarators lead to computational challenges in optimizing them. A well-known problem for multi-parameter optimization is getting stuck in local minima in parameter space, and the deviations from exact QS cannot generally be made arbitrarily small.

Significant analytical and numerical progress has been made in recent times in understanding the QS constraint, which has been achieved with very high precision (Landreman & Paul 2022). Near-axis expansions have provided helpful analytical insights into new previously unknown configurations and provided excellent initial guesses for further numerical optimization. It has also been possible to map out the entire quasisymmetric phase space using second-order NAEs. The important question of how the magnetic field strength shapes magnetic flux surfaces can be addressed within the NAE framework. For a local equilibrium, one is free to prescribe the flux-surface shape and the magnetic field strength (Boozer 2002; Skovoroda 2007). However, for a global quasisymmetric equilibrium, the flux-surface shapes are highly constrained (Jorge et al. 2020; Rodriguez 2022). Second-order NAE theory allows one to understand and explore this relationship (Rodriguez 2022; Rodríguez, Sengupta & Bhattacharjee 2023).

Despite the several advantages of NAE, a significant drawback, when applied to QS, is that physical quantities can be approximated only by low-order polynomials in the radial variable. The reason lies in the overdetermined nature of QS. Therefore, only plasma profile functions, such as the pressure, magnetic shear and bootstrap current that can be sufficiently well described by low-order polynomials, can be treated within the NAE framework for QS. A more serious drawback is that the overdetermination problem (discussed above) becomes an obstacle to calculating global magnetic shear, which shows up at the order at which the NAE approach breaks down. In the absence of reliable information on the global magnetic shear, MHD stability analyses, such as those pertaining to ballooning or interchange modes, become unreliable. Therefore, an alternative to the NAE is needed for the equilibrium and stability analysis of quasisymmetric devices.

In this work, we present an alternative approach to QS, which provides analytical insight akin to NAE but allows us to compute approximate global equilibrium solutions with more general profile functions. The expansion parameter in our reduced MHD model is the ratio of the equilibrium perpendicular and parallel length scales ($L_\perp /L_\parallel$) to the lowest-order magnetic field. We order plasma beta, defined by the ratio of the plasma and the magnetic pressure, to be of $O(L_\perp /L_\parallel )$. We further restrict the lowest-order magnetic field to be a purely toroidal vacuum field. With this restriction, we can only treat quasi-axisymmetric configurations. Plunk & Helander (2018) have shown that exact volumetric vacuum QS cannot be obtained close to vacuum axisymmetry under quite general conditions. Without any contradiction with Plunk & Helander (2018), we are able to realize approximate volumetric QS in our model since we satisfy the ideal MHD force balance and QS only to the lowest non-trivial orders. We are also consistent with the CDG (Constantin–Drivas–Ginzberg) theorem (Constantin, Drivas & Ginsberg 2021), which states that approximate volumetric QS can be obtained if the force-balance condition in ideal MHD is modified by the addition of a small force, allowing one to satisfy two of the weak QS conditions exactly. We have imposed the near-axisymmetry restriction to leverage the simplicity of the lowest-order axisymmetric geometry but shall relax this restriction in a future study.

The reduced MHD structure of our model coincides with the traditional large-aspect-ratio expansion (LAE) of ideal MHD (Hazeltine & Meiss 2003). An LAE model for stellarators, accurate to the first order in the inverse aspect ratio, with the plasma beta ordered as the inverse aspect ratio, is known in the literature (Wakatani 1998; Freidberg 2014) as the high-beta stellarator (HBS) model. (The rotational transform is finite in this model.) Unlike earlier large-aspect-ratio stellarator models (Greene & Johnson 1961; Strauss 1980; Strauss & Monticello 1981), the HBS does not require a large toroidal periodicity mode number, $N$. Given the high interest in quasisymmetric stellarators, we impose the QS constraint on the HBS and assess some of the important consequences for MHD equilibrium and stability.

The structure of the paper is as follows. In § 2, after discussing the implications of QS for ideal MHD equilibria, we provide a derivation of the quasisymmetric Grad–Shafranov equation (QS-GSE) in reduced MHD. We then discuss the quasisymmetric HBS (QS-HBS) model in § 3 and, in particular, we derive a special form of the QS-GSE. We present an analytic solution to the QS-GSE and its numerical verification in § 4. We discuss ballooning and interchange stability of the analytical equilibria in § 6. We conclude with a discussion of the implications of our results and possible generalizations in § 7.

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 2014; Burby et al. 2020; Rodríguez et al. 2020):

(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 2014) 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. 2020; Rodríguez et al. 2020), 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. 2015; Rodríguez & Bhattacharjee 2021; Huang et al. 2022, 2023) 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. 2020) (discussed in Appendix B) being too complicated, we shall use the reduced model for ideal MHD derived by Strauss (1997) 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 1971). Therefore, the Hamada condition (Helander 2014) 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 2019; Rodríguez et al. 2023) and

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

Following Freidberg (2014), 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 (2014), 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 (2014) 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. (2020). 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. 2020). Following Burby et al. (2020), 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 (2018), Rodríguez et al. (2023), Plunk & Helander (2018), Plunk (2020) and Burby et al. (2020). 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 2018). It can be shown self-consistently (Rodríguez et al. 2023) 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 (2018) and Plunk (2020), 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. (2020) 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. (2020). In general geometry, it is not enough to solve the GSE, and Burby et al. (2020) 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 2021). Implicitly, Hamada's condition is assumed to be satisfied by any MHD equilibrium with nested surfaces (Helander 2014). 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. 2012; Helander 2014) 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 1983), 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 1997). 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. 2023) the ballooning eigenmode equation (Connor, Hastie & Taylor 1978; Dewar & Glasser 1983)

(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 1958). Multiple studies by Redi et al. (2001, 2002) 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.

7 Discussion and conclusion

In this work, we have shown how the overdetermined problem of volumetric quasisymmetric MHD equilibrium can be approximately solved by utilizing the inverse aspect ratio of a stellarator as an expansion parameter. Our assumptions, namely finite rotational transform, large aspect ratio $\epsilon \ll 1$ and high plasma beta ($\beta \sim O(\epsilon )$), are consistent with the HBS model (Freidberg 2014). With a purely toroidal vacuum field and a constant field strength as our lowest-order magnetic field, we show that the quasiaxisymmetric MHD equilibrium can be described by a Grad–Shafranov equation, much like axisymmetry. In particular, we show that approximate quasi-axisymmetric equilibria can be obtained from large-aspect-ratio axisymmetric Grad–Shafranov solutions by breaking axisymmetry through a periodic up–down deformation of the surfaces. Our study of Mercier and ballooning stability of some of these equilibria shows that the non-axisymmetric deformations do not significantly degrade the MHD stability properties.

The original HBS model Freidberg (2014) describes MHD equilibrium in large-aspect-ratio stellarators. The HBS consists of two coupled nonlinear PDEs: magnetic differential equations for the pressure surface and the parallel current. Being fully three-dimensional and nonlinear, they inherit the same problems as 3-D MHD: the non-existence of nested flux surfaces, the possibility of current singularities and the formation of magnetic islands. As discussed in Freidberg (2014), the nonlinearity and three-dimensionality present a severe hindrance to constructing analytical equilibria. An exception is the Greene–Johnson (Greene & Johnson 1961) limit of the HBS model, $N\sim \epsilon ^{-1/2}\gg 1$, for a classical stellarator (Wakatani 1998; Freidberg 2014; Loizu et al. 2017; Baillod et al. 2023), where a GS equation can be obtained by an averaging over the short helical wavelength associated with large $N$.

In this work, we have imposed volumetric QS on the HBS model, allowing us to integrate the magnetic differential equation for the parallel current and obtain a 3-D Grad–Shafranov equation for the flux-surface label without imposing large $N$. The analysis of the fully 3-Dequilibria parallels that of an axisymmetric tokamak, thanks to the constraints from QS, consistent with the large-aspect-ratio limit of Burby et al. (2020). Since volumetric QS and MHD are not satisfied exactly but only to the lowest non-trivial order in $\epsilon$, our results are also consistent with Plunk, Landreman & Helander (2019), Plunk (2020), and the CDG theorem (Constantin et al. 2021).

Our work was motivated by the success of the second-order NAEs in providing an analytic description of quasisymmetric stellarators. While the near-axis approach works for any stellarator, provided we zoom in sufficiently close to the magnetic axis, the present approach requires the stellarator to have a large aspect ratio, which is quite reasonable for most stellarators. A significant advantage of our approach over the NAE is a global description of MHD equilibrium free from any polynomial expansion in the radial variable. Therefore, magnetic shear, currents and pressure profiles can be quite general. Furthermore, we can allow any flux-surface shaping. In contrast, the second-order NAE only allows triangularity in shaping.

The QS-GSE admits several classes of exact analytical solutions. In fact, any analytic solutions of the large-aspect-ratio axisymmetric Grad–Shafranov equation can generate a quasisymmetric solution through a straightforward coordinate transform that induces a periodic up–down deformation. Utilizing this, we have generated numerical solutions starting from an axisymmetric tokamak equilibrium and deforming it accordingly. We have then verified the validity of this methodology using the VMEC code. As the aspect ratio gets larger, the deviations from the periodically up–down shifted tokamak equilibrium, and the actual stellarator equilibrium generated with the same boundary using the VMEC code, become smaller.

There are several directions in which the current work can be extended. First, our current model can only describe quasiaxisymmetric configurations close to axisymmetry due to the choice of the lowest-order vacuum field being a purely toroidal field. A more general choice for the lowest-order vacuum field is needed to describe QA and quasihelical symmetry. Second, the isomorphism between QS-HBS and axisymmetric tokamak equilibrium suggests that we can carry out a local geometry analysis for QS-HBS in close analogy to the axisymmetric Miller-geometry case (Miller et al. 1998). Near-axisymmetric local equilibria could then be used to study kinetic ballooning modes, ion temparature gradient (ITG) and electron temperature gradient (ETG) thresholds. Third, one can study subsonic flows in quasisymmetric stellarators, which have been predicted (Spong 2005) and experimentally demonstrated (Gerhardt et al. 2005), using the QS-HBS model. The QS-HBS model should offer interesting analytical insights into the structure of subsonic flows and their damping. Fourth, the HBS model in the Greene–Johnson limit has been used to study the plasma $\beta$-limit in a stellarator, which is of practical importance (Freidberg 2014; Loizu et al. 2017; Baillod et al. 2023). An analogous study can be conducted for the QS-HBS model to understand how QS impacts the $\beta$ limit. Finally, we hope that the near-axisymmetric aspect of the QS-HBS model can provide an avenue to study non-axisymmetric deformations of tokamaks that preserved QS and hence can be successfully used in controlling tokamaks using 3-D perturbations (Park et al. 2021).