Hostname: page-component-cd9895bd7-gbm5v Total loading time: 0 Render date: 2024-12-25T18:19:15.721Z Has data issue: false hasContentIssue false

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

Published online by Cambridge University Press:  05 April 2024

Wrick Sengupta*
Affiliation:
Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08543, USA Princeton Plasma Physics Laboratory, Princeton, NJ 08540, USA
Nikita Nikulsin
Affiliation:
Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08543, USA Princeton Plasma Physics Laboratory, Princeton, NJ 08540, USA
Rahul Gaur
Affiliation:
Department of Mechanical and Aerospace Engineering, Princeton University, Princeton, NJ 08544, USA
Amitava Bhattacharjee
Affiliation:
Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08543, USA Princeton Plasma Physics Laboratory, Princeton, NJ 08540, USA
*
Email address for correspondence: wsengupta@princeton.edu

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.

Type
Research Article
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
Copyright © The Author(s), 2024. Published by Cambridge University Press

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 Reference Helander2014; Landreman & Catto Reference Landreman and Catto2010). Moreover, exact QSFootnote 1 leads to nested flux surfaces (Burby et al. Reference Burby, Kallinikos and MacKay2020; Rodríguez, Helander & Bhattacharjee Reference Rodríguez, Helander and Bhattacharjee2020). 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. Reference Park, Yang, Logan, Hu, Zhu, Zarnstorff, Nazikian, Paz-Soldan, Jeon and Ko2021), which are known to be sensitive to 3-D perturbations (Park, Boozer & Glasser Reference Park, Boozer and Glasser2007a; Park et al. Reference Park, Boozer, Menard, Garofalo, Schaffer, Hawryluk, Kaye, Gerhardt and Sabbagh2009). If the 3-D perturbations of an axisymmetric tokamak restore QS, they can help mitigate harmful transport effects (Park et al. Reference Park, Schaffer, Menard and Boozer2007b, Reference Park, Jeon, In, Ahn, Nazikian, Park, Kim, Lee, Ko and Kim2018).

Unfortunately, exact QS is difficult, if not impossible, to obtain. Analytical results (Garren & Boozer Reference Garren and Boozer1991b, Reference Garren and Boozera; Plunk & Helander Reference Plunk and Helander2018; Landreman & Sengupta Reference Landreman and Sengupta2019; Jorge, Sengupta & Landreman Reference Jorge, Sengupta and Landreman2020) 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. Reference Anderson, Almagri, Anderson, Matthews, Talmadge and Shohet1995; Zarnstorff et al. Reference Zarnstorff, Berry, Brooks, Fredrickson, Fu, Hirshman, Hudson, Ku, Lazarus and Mikkelsen2001; Najmabadi et al. Reference Najmabadi, Raffray, Abdel-Khalik, Bromberg, Crosatti, El-Guebaly, Garabedian, Grossman, Henderson and Ibrahim2008; Ku & Boozer Reference Ku and Boozer2010; Bader et al. Reference Bader, Drevlak, Anderson, Faber, Hegna, Likin, Schmitt and Talmadge2019; Landreman & Paul Reference Landreman and Paul2022). 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 Reference Landreman and Paul2022). 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 Reference Boozer2002; Skovoroda Reference Skovoroda2007). However, for a global quasisymmetric equilibrium, the flux-surface shapes are highly constrained (Jorge et al. Reference Jorge, Sengupta and Landreman2020; Rodriguez Reference Rodriguez2022). Second-order NAE theory allows one to understand and explore this relationship (Rodriguez Reference Rodriguez2022; Rodríguez, Sengupta & Bhattacharjee Reference Rodríguez, Sengupta and Bhattacharjee2023).

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 (Reference Plunk and Helander2018) have shown that exact volumetric vacuum QS cannot be obtained close to vacuum axisymmetry under quite general conditions. Without any contradiction with Plunk & Helander (Reference Plunk and Helander2018), 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 Reference Constantin, Drivas and Ginsberg2021), 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 Reference Hazeltine and Meiss2003). 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 Reference Wakatani1998; Freidberg Reference Freidberg2014) as the high-beta stellarator (HBS) model. (The rotational transform is finite in this model.) Unlike earlier large-aspect-ratio stellarator models (Greene & Johnson Reference Greene and Johnson1961; Strauss Reference Strauss1980; Strauss & Monticello Reference Strauss and Monticello1981), 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 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.5ad) 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.7ac) 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.33ac), 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.37ac) 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.

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 Reference Freidberg2014). 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 (Reference Freidberg2014) 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 (Reference Freidberg2014), the nonlinearity and three-dimensionality present a severe hindrance to constructing analytical equilibria. An exception is the Greene–Johnson (Greene & Johnson Reference Greene and Johnson1961) limit of the HBS model, $N\sim \epsilon ^{-1/2}\gg 1$, for a classical stellarator (Wakatani Reference Wakatani1998; Freidberg Reference Freidberg2014; Loizu et al. Reference Loizu, Hudson, Nührenberg, Geiger and Helander2017; Baillod et al. Reference Baillod, Loizu, Qu, Arbez and Graves2023), 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. (Reference Burby, Kallinikos and MacKay2020). 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 (Reference Plunk, Landreman and Helander2019), Plunk (Reference Plunk2020), and the CDG theorem (Constantin et al. Reference Constantin, Drivas and Ginsberg2021).

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. Reference Miller, Chu, Greene, Lin-Liu and Waltz1998). 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 Reference Spong2005) and experimentally demonstrated (Gerhardt et al. Reference Gerhardt, Talmadge, Canik and Anderson2005), 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 Reference Freidberg2014; Loizu et al. Reference Loizu, Hudson, Nührenberg, Geiger and Helander2017; Baillod et al. Reference Baillod, Loizu, Qu, Arbez and Graves2023). 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. Reference Park, Yang, Logan, Hu, Zhu, Zarnstorff, Nazikian, Paz-Soldan, Jeon and Ko2021).

Acknowledgements

The authors would like to thank F. P. Diaz, V. Duarte, H. Weitzner, E. J. Paul, M. Zarnstroff, J. Loizu, A. Baillod, H. Zhu, R. Jorge, E. Rodríguez and M. Landreman for stimulating discussions and helpful suggestions.

Editor Per Helander thanks the referees for their advice in evaluating this article.

Funding

This research was supported by a grant from the Simons Foundation/SFARI (560651, AB), and the Department of Energy Award No. DE-SC0024548.

Declaration of interests

The authors report no conflict of interest.

Data availability

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Appendix A. Consistency of QS-GSE and ${\nabla}$$_{||}\varPsi =0$

We shall now investigate the consistency of the following set of equations:

(A1a)\begin{gather} \nabla_\perp^2 \varPsi+ 2 x\frac{{\rm d}\beta(\varPsi)}{{\rm d}\varPsi} + H'(\varPsi)-\mathfrak{a}_{,xx}=0, \quad \nabla_\perp^2 = \partial_x^2+\partial_y^2 , \end{gather}
(A1b)\begin{gather}\mathcal{L} \varPsi=0, \quad \mathcal{L}= \left( \partial_\phi - \mathfrak{a}_{,x}\partial_y\right) , \end{gather}

obtained from (3.28) and (3.29). The consistency condition follows from the fact that $\mathcal {L}$ on (A1a) should vanish identically. We shall use the following identities, which follow from straightforward, direct calculations using (A1b):

(A2a)\begin{gather} \mathcal{L} \nabla_\perp^2 \varPsi= \left( \partial_x^3 \mathfrak{a}\right) \partial_y \varPsi +2 \left( \partial_x^2 \mathfrak{a} \right)\partial^2_{xy}\varPsi, \end{gather}
(A2b)\begin{gather}\mathcal{L} \left( 2 x\frac{{\rm d}\beta(\varPsi)}{{\rm d}\varPsi} + H'(\varPsi) \right) =0, \end{gather}
(A2c)\begin{gather}\mathcal{L} \mathfrak{a}_{,xx}= \partial_\phi (-\mathfrak{a}_{,xx}). \end{gather}

Using the above identities, we can rewrite the consistency condition as an equation for $\mathfrak {a}$

(A3a,b)\begin{equation} \left( \varPsi_{,y} \partial_x +2 \varPsi_{,xy}-\partial_\phi\right) \zeta=0, \quad \zeta=\mathfrak{a}_{,xx}. \end{equation}

Using the fact that $\mathfrak {a}$ is $y$-independent, we can rewrite (A3a,b) as

(A4)\begin{equation} \partial^2_{x y} \left( \sqrt{\zeta} \varPsi\right) - \partial_\phi \sqrt{\zeta}=0. \end{equation}

Averaging over $\phi$ and using periodicity of $\mathfrak {a}, \zeta$ in $\phi$ we obtain

(A5)\begin{equation} \partial^2_{xy} \oint {\rm d}\phi \,\sqrt{\zeta}\;\varPsi =0 \quad \Leftrightarrow \quad \oint {\rm d}\phi\, \sqrt{\zeta}\,\varPsi = \text{const}. \end{equation}

We note that $\varPsi$ depends on all three variables whereas $\mathfrak {a}, \zeta$ depends only on $(x,\phi )$. Thus, satisfying (A3a,b) or the constraint (A5), with a non-zero $\zeta$, is in general not possible. Thus, a possible solution is $\mathfrak {a}_{,xx}=0$, which implies that $\mathfrak {a}$ is linear in $x$. Thus, we obtain (3.34).

Appendix B. Comparison with Burby–Kallinikos–MacKay

Burby et al. (Reference Burby, Kallinikos and MacKay2020) obtained the following generalized QS-GSE:

(B1)\begin{equation} \nabla^2 \psi - \frac{\boldsymbol{u} \times \boldsymbol{v}}{u^2}\boldsymbol{\cdot} \boldsymbol{\nabla} \psi + \frac{\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{v}}{u^2}\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{B} -\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{J} =0, \end{equation}

where $\boldsymbol {u}$ is the QS vector, $\boldsymbol {v}=\boldsymbol {\nabla } \times \boldsymbol {u}$. Note that in (B1), $\psi$ denotes $\boldsymbol {A}\boldsymbol {\cdot } \boldsymbol {u}$ with $\boldsymbol {A}$ denoting the vector potential. For MHD equilibrium, $\boldsymbol {u}\boldsymbol {\cdot } \boldsymbol {B} =F(\psi )$ and (B1) takes the form

(B2)\begin{equation} \boldsymbol{\nabla}^2 \psi - \frac{\boldsymbol{u} \times \boldsymbol{v}}{u^2}\boldsymbol{\cdot} \boldsymbol{\nabla} \psi + \frac{\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{v}}{u^2}F +FF'-u^2 p'(\psi)=0. \end{equation}

For axisymmetry and helical symmetry, we have, respectively

(B3)\begin{gather} \left.\begin{gathered} \boldsymbol{u}=R^2\boldsymbol{\nabla}\phi,\quad u=R,\quad \boldsymbol{v}={-}2 \boldsymbol{e}_{\boldsymbol{Z}},\quad \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{v}=0,\quad \boldsymbol{u}\times \boldsymbol{v}=2 R \boldsymbol{e}_{\boldsymbol{R}}, \\ \boldsymbol{u}=R^2\boldsymbol{\nabla}\phi-l \boldsymbol{\nabla} Z,\quad u=\sqrt{R^2+l^2},\quad \boldsymbol{v}={-}2 \boldsymbol{e}_{\boldsymbol{Z}}, \quad \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{v} = 2 l,\quad \boldsymbol{u}\times \boldsymbol{v}=2 R \boldsymbol{e}_{\boldsymbol{R}}. \end{gathered}\right\} \end{gather}

Equation (B2) reduces to the standard axisymmetric and helically symmetric Grad–Shafranov equations, respectively.

We now derive the expression for $\boldsymbol {u}$. We use the following identities, obtained in §§ 3.1 and 3.2:

(B4a-d)\begin{equation} \frac{{\rm d}\varPsi}{{\rm d}\psi}=\frac{1}{\epsilon F}, \quad \frac{F}{B_0}=\frac{1}{\epsilon \iota}(1-\epsilon \beta), \quad \frac{B}{B_0}=1-\epsilon(x+\beta), \quad \boldsymbol{\nabla}_{||} x ={-}\partial_\xi\varPsi(x,\xi). \end{equation}

Starting with the definition for $\boldsymbol {u}$ (2.6a,b), we find that

(B5)\begin{equation} \boldsymbol{u} = \frac{\boldsymbol{\nabla} \varPsi \times \boldsymbol{\nabla} B}{B \left( \epsilon \boldsymbol{\nabla}_{||} B\right)}\left( \frac{{\rm d}\varPsi}{{\rm d}\psi} \right)^{{-}1}= \frac{F/B_0}{B/B_0} \frac{\boldsymbol{\nabla} \varPsi \times \boldsymbol{\nabla} x}{\boldsymbol{\nabla}_{||} x}. \end{equation}

Using the identities (B4ad) together with the definition of $\xi$ (3.30a,b), $\boldsymbol {u}$ takes the form

(B6)\begin{equation} \boldsymbol{u}=\frac{1}{\epsilon \iota}(1-\epsilon \beta)(1+\epsilon (x+\beta))\frac{\boldsymbol{\nabla} \xi\times \boldsymbol{\nabla} x}{\boldsymbol{\nabla}_{||} x}\varPsi_{,\xi} = \frac{1}{\epsilon \iota}\left( \boldsymbol{e}_{\boldsymbol{\phi}}(1+\epsilon x)-\epsilon \mathfrak{a}_{,x}\boldsymbol{e}_{\boldsymbol{Z}}\right). \end{equation}

The expression for $\boldsymbol {u}_{\textrm {HBS}}=\iota \boldsymbol {u}$ is

(B7)\begin{equation} \boldsymbol{u}_{\text{HBS}}= R^2 \boldsymbol{\nabla} \phi -\mathfrak{a}_{,x}\boldsymbol{\nabla} Z = \boldsymbol{e}_{\boldsymbol{\phi}}\left( \frac{1}{\epsilon}+x\right) -\mathfrak{a}_{,x}\boldsymbol{e}_{\boldsymbol{Z}}, \end{equation}

which follows from (B6) and the definition $R=\epsilon ^{-1}+x$.

We now show that (B2) for HBS takes the same form as the QS-GSE (3.36). Substituting the following identities into (B2):

(B8a)\begin{gather} |\boldsymbol{u}_{\text{HBS}}|=\frac{1}{\epsilon}+x,\quad \boldsymbol{v}_{\text{HBS}}={-}2 \boldsymbol{e}_{\boldsymbol{Z}}-\epsilon \mathfrak{a}_{x\phi}\boldsymbol{e}_{\boldsymbol{R}},\quad FF'={-}\frac{1}{\epsilon}\beta' +O(1) \end{gather}
(B8b)\begin{gather}\frac{\boldsymbol{u}_{\text{HBS}}\boldsymbol{\cdot} \boldsymbol{v}_{\text{HBS}}}{|\boldsymbol{u}_{\text{HBS}}|^2} ={-} 2\mathfrak{a}_{,x} \epsilon^2, \quad \frac{\boldsymbol{u}_{\text{HBS}}\times \boldsymbol{v}_{\text{HBS}}}{|\boldsymbol{u}_{\text{HBS}}|^2} = 2 \epsilon \boldsymbol{e}_{\boldsymbol{R}}, \end{gather}

we find that both the $\boldsymbol {u}\times \boldsymbol {v}$ and the $\boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {v}$ terms in (B2) are negligible to the lowest order in $\epsilon$. Also, the $FF'$ and the pressure term cancel to the lowest order, leaving the $2x\beta '$ term in (3.36). As a result, we get (3.36) with the identification that $H'(\varPsi )$ is the $O(1)$ term of the expression $FF'$.

Finally, let us show that $\boldsymbol {u}_{\textrm {HBS}}$ is not a killing vector. As discussed in Burby et al. (Reference Burby, Kallinikos and MacKay2020), a vector $\boldsymbol {u}$ is killing if, and only if, the vector $\boldsymbol {w}$ is identically zero, where

(B9a,b)\begin{equation} \boldsymbol{w}={-}\boldsymbol{u} \times \boldsymbol{v} + \boldsymbol{\nabla} u^2, \quad \boldsymbol{v} =\boldsymbol{\nabla} \times \boldsymbol{u}. \end{equation}

Evaluating $\boldsymbol {w}_{\textrm {HBS}}$ we find it to be non-zero since

(B10)\begin{equation} \boldsymbol{w}_{\text{HBS}}={-}\mathfrak{a}_{,x\phi} \boldsymbol{e}_{\boldsymbol{Z}}= Y''(\phi) \boldsymbol{e}_{\boldsymbol{Z}}\neq 0. \end{equation}

Appendix C. Various coordinate transformations

Navigating between different coordinate systems is essential for simplifying PDEs and obtaining solutions. In this work, we find the following coordinate systems, besides the HBS coordinates $(x,y,\phi )$, to be advantageous:

(C1a-d)\begin{equation} (x,\xi,\phi), \quad (\varPsi,\xi,\phi), \quad (\varPsi,\alpha,\phi), \quad \text{and}\quad (\varPsi,\alpha,\ell). \end{equation}

We begin with the transformation from $(x,y,\phi )$ to $(x,\xi = y +Y(\phi ),\phi )$,

(C2a-c)\begin{equation} \left.\partial_x\right|_{(\kern0.7pt y,\phi)}=\left.\partial_x\right|_{(\xi,\phi)}, \quad {\partial_y}|_{(x,y)}={\partial_\xi}|_{(x,\phi)}, \quad {\partial_\phi}|_{(x,y)}= {\partial_\phi}|_{(x,\xi)}+Y'(\phi)\partial_\xi. \end{equation}

The $\boldsymbol {\nabla }_{||}$ operator in $(x,y,\phi )$ is

(C3)\begin{equation} \boldsymbol{\nabla}_{||} = \partial_\phi -Y'(\phi)\partial_y +\{\varPsi,\ \}_{(x,y)}, \end{equation}

which follows from $\mathfrak {a}_{,x}=Y'(\phi )$, and the definition of $\boldsymbol {\nabla }_{||}$ given in (3.23a,b). Using (C2ac), the expression for $\boldsymbol {\nabla }_{||}$ can be rewritten as

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

Thus, $\boldsymbol {\nabla }_{||}\varPsi = \partial _\phi \varPsi =0$ in the $(x,\xi,\phi )$ coordinates, pointing to the 2-D nature of $\varPsi$.

We now exploit the two-dimensionality of $\varPsi$ to transform to the $(\varPsi,\xi,\phi )$ coordinates using

(C5a-c)\begin{equation} \left.\partial_x\right|_{(\kern0.7pt y,\phi)}=\varPsi_{,x}\left.\partial_\varPsi\right|_{(\xi,\phi)}, \quad {\partial_\xi}|_{(x,y)}={\partial_\xi}|_{(x,y)}+\varPsi_{,\xi}\left.\partial_\varPsi\right|_{(\xi,\phi)}, \quad {\partial_\phi}|_{(x,\xi)}= {\partial_\phi}|_{(\varPsi,\xi)}. \end{equation}

The parallel derivative $\boldsymbol {\nabla }_{||}$ (C4), now takes the form

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

It follows from the fact that $\varPsi _{,x}$ is only a function of $\varPsi,\xi$, $\boldsymbol {\nabla }_{||} \alpha =0$ and $\epsilon \boldsymbol {\nabla }_{||}\ell =1$ that

(C7a,b)\begin{equation} \alpha = \phi - \int_\varPsi \frac{{\rm d}\xi}{\varPsi_{,x}},\quad \epsilon \ell = \phi +L(\varPsi). \end{equation}

To check the validity of (C7a,b), we now check that we can recover $\boldsymbol {B}$ as given in (3.25a) from $\boldsymbol {\nabla }\alpha \times \boldsymbol {\nabla }\varPsi$. From (C7a,b), we have

(C8)\begin{equation} \boldsymbol{\nabla} \alpha \times \boldsymbol{\nabla} \varPsi ={-}\boldsymbol{\nabla} \varPsi \times \boldsymbol{\nabla} \phi +\boldsymbol{\nabla} \varPsi \times \boldsymbol{\nabla} \xi \frac{1}{\varPsi_{,x}}. \end{equation}

We can simplify the last term in (C8) in $(x,y,\phi )$ coordinates using $\varPsi =\varPsi (x,\xi ), \xi =y+Y(\phi ), \mathfrak {a}_{,x}=Y'(\phi )$. We obtain

(C9)\begin{align} \boldsymbol{\nabla} \varPsi \times \boldsymbol{\nabla} \xi \frac{1}{\varPsi_{,x}}& = \varPsi_{,x}\boldsymbol{\nabla} x \times \boldsymbol{\nabla} (\kern0.7pt y+Y(\phi)) \frac{1}{\varPsi_{,x}},\nonumber\\ & = \boldsymbol{e}_{\boldsymbol{\phi}}+Y'(\phi)\boldsymbol{\nabla} x \times \boldsymbol{\nabla} \phi,\nonumber\\ & = \boldsymbol{e}_{\boldsymbol{\phi}}+\boldsymbol{\nabla} \mathfrak{a} \times \boldsymbol{\nabla} \phi. \end{align}

Therefore

(C10)\begin{equation} \boldsymbol{\nabla} \alpha \times \boldsymbol{\nabla} \varPsi =\boldsymbol{e}_{\boldsymbol{\phi}}+\boldsymbol{\nabla} (\mathfrak{a}-\varPsi) \times \boldsymbol{\nabla} \phi = \boldsymbol{e}_{\boldsymbol{\phi}}+\boldsymbol{\nabla}_{\boldsymbol{\perp}} \mathcal{A} \times \boldsymbol{\nabla} \phi, \end{equation}

upon using $\varPsi +\mathcal {A}=\mathfrak {a}$.

Comparing with the expression of $\boldsymbol {B}$ from (3.25a), we find that the $\boldsymbol {e}_{\boldsymbol {\phi }}(x+\beta )$ terms are mixing. Including the $O(\epsilon )$ correction to $\alpha$ such that

(C11)\begin{equation} \alpha = \phi - \int_\varPsi {\rm d}\xi \,\frac{B_\phi/B_0}{\varPsi_{,x}}, \end{equation}

remits this problem. We need not just $\alpha$ but also $\boldsymbol {\nabla } \alpha$ to be accurate to $O(\epsilon )$. The first term in $\alpha$, i.e. $\phi$, has a gradient of $O(\epsilon )$. Thus, we also need the $O(\epsilon )$ piece from $B_\phi /B_0$ to contribute to $\boldsymbol {\nabla } \alpha$ to $O(\epsilon )$.

Appendix D. On-axis rotational transform from the integrated torsion

The axis of the QS-HBS system is a circle deformed periodically in the $Z$ direction. Let us now obtain an expression for total torsion ($\oint \tau \,\textrm {d}\ell$), representing the torsion contribution to the on-axis rotational transform formula due to Mercier.

Using cylindrical coordinates $(R, Z, \phi )$, we can denote the axis by the following space curve:

(D1)\begin{equation} \boldsymbol{r}=R_0 (\boldsymbol{e}_{\boldsymbol{R}}+\mathfrak{f}(\phi)\boldsymbol{e}_{\boldsymbol{Z}}). \end{equation}

Here, the first term represents a circle of constant radius $R_0$, and $\mathfrak {f}$ is a periodic deformation of the same. Using standard relations

(D2a,b)\begin{equation} \partial_\phi \boldsymbol{e}_{\boldsymbol{R}}=\boldsymbol{e}_{\boldsymbol{\phi}} , \quad \partial_\phi \boldsymbol{e}_{\boldsymbol{\phi}}={-}\boldsymbol{e}_{\boldsymbol{R}}, \end{equation}

we find the arclength $\ell$, and the unit tangent vector to be

(D3a,b)\begin{equation} \frac{{\rm d}\ell}{{\rm d}\phi}= R_0\sqrt{1+\mathfrak{f}'(\phi)^2}, \quad \boldsymbol{t}\equiv \frac{{\rm d}\boldsymbol{r}}{{\rm d}\ell} = \frac{(\boldsymbol{e}_{\boldsymbol{\phi}}+\mathfrak{f}' \boldsymbol{e}_{\boldsymbol{Z}})}{\sqrt{1+\mathfrak{f}'^2}}. \end{equation}

We shall use the following identities for the curvature and torsion, $\kappa,\tau$:

(D4a,b)\begin{equation} \kappa= \left|\frac{{\rm d}\boldsymbol{r}}{{\rm d}\ell}\times \frac{{\rm d}^2\boldsymbol{r}}{{\rm d}\ell^2} \right|, \quad \tau = \frac{1}{\kappa^2}\frac{{\rm d}\boldsymbol{r}}{{\rm d}\ell}\times \frac{{\rm d}^2\boldsymbol{r}}{{\rm d}\ell^2}\boldsymbol{\cdot} \frac{{\rm d}^3\boldsymbol{r}}{{\rm d}\ell^3}. \end{equation}

Straightforward calculation of the higher derivatives of $\boldsymbol {r}$ leads to the following expressions for $\kappa,\tau$:

(D5a,b)\begin{equation} \kappa = \frac{1}{R_0} \sqrt{\frac{1+\mathfrak{f}'^2 +\mathfrak{f}''^2}{(1+\mathfrak{f}'^2)^3}}, \quad \tau = \frac{1}{R_0} \frac{\mathfrak{f}' +\mathfrak{f}'''}{1+\mathfrak{f}'^2 +\mathfrak{f}''^2}. \end{equation}

Therefore

(D6)\begin{equation} \oint \tau \,{\rm d}\ell = \oint {\rm d}\phi \,\sqrt{1+\mathfrak{f}'^2}\frac{\mathfrak{f}' +\mathfrak{f}'''}{1+\mathfrak{f}'^2 +\mathfrak{f}''^2}. \end{equation}

In general, this is a non-vanishing quantity. However, since the deformation $\mathfrak {f}$ can be only as large as the minor radius, i.e. $\mathfrak {f}= \epsilon \mathfrak {f}_1$, the integrated torsion is not generally an $O(1)$ quantity as we now show.

We find that

(D7)\begin{equation} \oint \tau \,{\rm d}\ell = \epsilon \oint {\rm d}\phi (\mathfrak{f}_1' +\mathfrak{f}_1''') + \epsilon^3 \oint {\rm d}\phi (\mathfrak{f}_1' +\mathfrak{f}_1''')(\mathfrak{f}''_1)^2 +O(\epsilon^5). \end{equation}

The first term in the expression for torsion averages out identically. The next term is $O(\epsilon ^3)$ if the toroidal harmonic, $n\sim 1$. However, for sizable $n$, we can get an amplification factor coming from the derivatives of $\mathfrak {f}$. To estimate this factor, we note that the total torsion can be further simplified to

(D8)\begin{equation} \oint \tau \,{\rm d}\ell = \epsilon^3 \oint {\rm d}\phi\,\mathfrak{f}'_1 \mathfrak{f}_1''^2, \end{equation}

since the $\mathfrak {f}_1'''$ term is a total derivative.

We observe several facts from (D8). First, if $\mathfrak {f}_1(\phi )$ is even or simple harmonic, i.e. $\exp (\textrm {i} n \phi )$, the total torsion vanishes. The reason why a simple harmonic $\mathfrak {f}_1$ leads to a zero total torsion is that, in this case, we can always add a term $n^2 \mathfrak {f}_1''' f_1''^2$, which vanishes by itself. However, this term cancels the $\mathfrak {f}_1'f_1''^2$ term due to the simple-harmonic nature of $\mathfrak {f}_1$. So, consider a two-harmonic stellarator-symmetric deformation: $\mathfrak {f}_1 = a\sin {n\phi } + b\sin {m\phi }$. After some straightforward algebra, one obtains

(D9)\begin{align} \oint\mathfrak{f}_1'\mathfrak{f}_1''\,{\rm d}\phi & = \int_0^{2{\rm \pi}} \left[a^2 bmn\left(mn^2\sin{m\phi}\sin{2n\phi} - \frac{n^3}{2}\cos{m\phi}\cos{2n\phi}\right)\right.\nonumber\\ & \quad \left. +\, ab^2 mn\left(m^2 n\sin{n\phi}\sin{2m\phi} - \frac{m^3}{2}\cos{n\phi}\cos{2m\phi}\right)\right]{\rm d}\phi. \end{align}

Clearly, in order for the integral to be non-zero, it must be that either $m=2n$ or $n=2m$. Without loss of generality, take $m=2n$. Then

(D10)\begin{equation} \iota_\tau = \frac{1}{2{\rm \pi}}\oint\tau \,{\rm d}\ell = \frac{\epsilon^3}{2{\rm \pi}}\oint\mathfrak{f}_1'\mathfrak{f}_1''\,{\rm d}\phi = \frac{3}{2}\epsilon^3 a^2 bn^5. \end{equation}

We note the strong amplification factor in the above expression for the rotational transform due to the $n^5$ scaling. However, note that, if the ordering $\partial _\phi \sim \epsilon \nabla _\perp$ is to be valid, one must have $a\sim 1/n$ and $b\sim 1/m$; and thus $\iota _\tau = O(\epsilon ^3 n^2)$. This implies that low-amplitude large-$n$ modes could have a sizable torsional contribution to $\iota$. However, very large $n$ would lead to higher derivatives growing faster, leading to the breakdown of the LAE. Hence, $n\leqslant 3$ is perhaps the best choice for $n$, which limits the torsional contribution to $\iota _\tau \sim 0.05$ for an aspect ratio of ${\sim }5$.

Footnotes

1 In the literature, a distinction is made between the so-called strong and weak forms of QS. However, the two forms are equivalent up to scaling by a flux function for the ideal MHD force balance considered here (stated without proof in Constantin, Drivas & Ginsberg (Reference Constantin, Drivas and Ginsberg2021), but deducible from Theorems IX.11 or 13 of Burby, Kallinikos & MacKay (Reference Burby, Kallinikos and MacKay2020) modulo appropriate changes taking into account the scale factor), under some non-degeneracy conditions.

References

Anderson, F.S.B., Almagri, A.F., Anderson, D.T., Matthews, P.G., Talmadge, J.N. & Shohet, J.L. 1995 Helically symmetric experiment, (HSX) goals, design and status. Fusion Technol. 27 (3T), 273277. doi: 10.13182/FST95-A11947086.27CrossRefGoogle Scholar
Anderson, P.W. 1958 Absence of diffusion in certain random lattices. Phys. Rev. 109 (5), 1492.CrossRefGoogle Scholar
Bader, A., Drevlak, M., Anderson, D.T., Faber, B.J., Hegna, C.C., Likin, K.M., Schmitt, J.C. & Talmadge, J.N. 2019 Stellarator equilibria with reactor relevant energetic particle losses. J. Plasma Phys. 85 (5), 905850508.CrossRefGoogle Scholar
Baillod, A., Loizu, J., Qu, Z.S., Arbez, H.P. & Graves, J.P. 2023 Equilibrium $\beta$-limits dependence on bootstrap current in classical stellarators. J. Plasma Phys. 89 (5), 905890508.CrossRefGoogle Scholar
Boozer, A.H. 2002 Local equilibrium of nonrotating plasmas. Phys. Plasmas 9 (9), 37623766.CrossRefGoogle Scholar
Burby, J.W., Duignan, N. & Meiss, J.D. 2021 Integrability, normal forms, and magnetic axis coordinates. J. Maths Phys. 62, 122901.CrossRefGoogle Scholar
Burby, J.W., Kallinikos, N. & MacKay, R.S. 2020 Some mathematics for quasi-symmetry. J. Maths Phys. 61 (9), 093503.CrossRefGoogle Scholar
Connor, J.W., Hastie, R.J. & Taylor, J.B. 1978 Shear, periodicity, and plasma ballooning modes. Phys. Rev. Lett. 40 (6), 396.CrossRefGoogle Scholar
Constantin, P., Drivas, T.D. & Ginsberg, D. 2021 On quasisymmetric plasma equilibria sustained by small force. J. Plasma Phys. 87 (1), 905870111.CrossRefGoogle Scholar
Dewar, R.L. & Glasser, A.H. 1983 Ballooning mode spectrum in general toroidal systems. Phys. Fluids 26 (10), 30383052.CrossRefGoogle Scholar
D'haeseleer, W.D, Hitchon, W.N.G., Callen, J.D. & Shohet, J.L. 2012 Flux Coordinates and Magnetic Field Structure: A Guide to a Fundamental Tool of Plasma Theory. Springer Science & Business Media.Google Scholar
Freidberg, J.P. 2014 Ideal MHD. Cambridge University Press.CrossRefGoogle Scholar
Garren, D.A. & Boozer, A.H. 1991 a Existence of quasihelically symmetric stellarators. Phys. Fluids B 3 (10), 28222834.CrossRefGoogle Scholar
Garren, D.A. & Boozer, A.H. 1991 b Magnetic field strength of toroidal plasma equilibria. Phys. Fluids B 3 (10), 28052821.CrossRefGoogle Scholar
Gaur, R., Buller, S., Ruth, M.E., Landreman, M., Abel, I.G. & Dorland, W.D. 2023 An adjoint-based method for optimising mhd equilibria against the infinite-n, ideal ballooning mode. J. Plasma Phys. 89 (5), 905890518.CrossRefGoogle Scholar
Gerhardt, S.P., Talmadge, J.N., Canik, J.M. & Anderson, D.T. 2005 Experimental evidence of reduced plasma flow damping with quasisymmetry. Phys. Rev. Lett. 94 (1), 015002.CrossRefGoogle ScholarPubMed
Grad, H. 1971 Plasma containment in closed line systems. In Plasma Physics and Controlled Nuclear Fusion Research 1971. Vol. III. Proceedings of the Fourth International Conference on Plasma Physics and Controlled Nuclear Fusion Research.Google Scholar
Greene, J.M. 1997 A brief review of magnetic wells. Comments Plasma Phys. Control. Fusion 17, 389402.Google Scholar
Greene, J.M. & Johnson, J.L. 1961 Determination of hydromagnetic equilibria. Phys. Fluids 4 (7), 875890.CrossRefGoogle Scholar
Hazeltine, R.D. & Meiss, J.D. 2003 Plasma Confinement. Courier Corporation.Google Scholar
Helander, P. 2014 Theory of plasma confinement in non-axisymmetric magnetic fields. Rep. Prog. Phys. 77 (8), 087001.CrossRefGoogle ScholarPubMed
Hirshman, S.P. & Whitson, J.C. 1983 Steepest-descent moment method for three-dimensional magnetohydrodynamic equilibria. Phys. Fluids 26 (12), 35533568. https://doi.org/10.1063/1.864116CrossRefGoogle Scholar
Huang, Y.-M., Hudson, S.R., Loizu, J., Zhou, Y. & Bhattacharjee, A. 2022 Numerical study of δ-function current sheets arising from resonant magnetic perturbations. Phys. Plasmas 29 (3), 032513.CrossRefGoogle Scholar
Huang, Y.-M., Zhou, Y., Loizu, J., Hudson, S. & Bhattacharjee, A. 2023 Structure of pressure-gradient-driven current singularity in ideal magnetohydrodynamic equilibrium. Plasma Phys. Control. Fusion 65 (3), 034008.CrossRefGoogle Scholar
Jorge, R., Sengupta, W. & Landreman, M. 2020 Construction of quasisymmetric stellarators using a direct coordinate approach. Nucl. Fusion 60 (7), 076021.CrossRefGoogle Scholar
Ku, L.P. & Boozer, A.H. 2010 New classes of quasi-helically symmetric stellarators. Nucl. Fusion 51 (1), 013004.CrossRefGoogle Scholar
Landreman, M. & Catto, P.J. 2010 Effects of the radial electric field in a quasisymmetric stellarator. Plasma Phys. Control. Fusion 53 (1), 015004.CrossRefGoogle Scholar
Landreman, M. & Paul, E. 2022 Magnetic fields with precise quasisymmetry for plasma confinement. Phys. Rev. Lett. 128 (3), 035001.CrossRefGoogle ScholarPubMed
Landreman, M. & Sengupta, W. 2018 Direct construction of optimized stellarator shapes. Part 1. Theory in cylindrical coordinates. J. Plasma Phys. 84 (6), 905840616.CrossRefGoogle Scholar
Landreman, M. & Sengupta, W. 2019 Constructing stellarators with quasisymmetry to high order. J. Plasma Phys. 85 (6), 815850601.CrossRefGoogle Scholar
Loizu, J., Hudson, S.R., Bhattacharjee, A., Lazerson, S. & Helander, P. 2015 Existence of three-dimensional ideal-magnetohydrodynamic equilibria with current sheets. Phys. Plasmas 22 (9), 090704.CrossRefGoogle Scholar
Loizu, J., Hudson, S.R., Nührenberg, C., Geiger, J. & Helander, P. 2017 Equilibrium $\beta$-limits in classical stellarators. J. Plasma Phys. 83 (6), 715830601.CrossRefGoogle Scholar
Miller, R.L., Chu, M.S., Greene, J.M., Lin-Liu, Y.R. & Waltz, R.E. 1998 Noncircular, finite aspect ratio, local equilibrium model. Phys. Plasmas 5 (4), 973978.CrossRefGoogle Scholar
Najmabadi, F., Raffray, A.R., Abdel-Khalik, S.I., Bromberg, L., Crosatti, L., El-Guebaly, L., Garabedian, P.R., Grossman, A.A., Henderson, D., Ibrahim, A., et al. 2008 The ARIES-CS compact stellarator fusion power plant. Fusion Sci. Technol. 54 (3), 655672.CrossRefGoogle Scholar
Park, J.-k., Boozer, A.H. & Glasser, A.H. 2007 a Computation of three-dimensional tokamak and spherical torus equilibria. Phys. Plasmas 14 (5), 052110.CrossRefGoogle Scholar
Park, J.-k., Boozer, A.H., Menard, J.E., Garofalo, A.M., Schaffer, M.J., Hawryluk, R.J., Kaye, S.M., Gerhardt, S.P., Sabbagh, S.A., NSTX Team, et al. 2009 Importance of plasma response to nonaxisymmetric perturbations in tokamaks. Phys. Plasmas 16 (5), 056115.CrossRefGoogle Scholar
Park, J.-K., Jeon, Y.M., In, Y., Ahn, J.-W., Nazikian, R., Park, G., Kim, J., Lee, H.H., Ko, W.H., Kim, H.-S., et al. 2018 3D field phase-space control in tokamak plasmas. Nat. Phys. 14 (12), 12231228.CrossRefGoogle Scholar
Park, J.-k., Schaffer, M.J., Menard, J.E. & Boozer, A.H. 2007 b Control of asymmetric magnetic perturbations in tokamaks. Phys. Rev. Lett. 99 (19), 195003.CrossRefGoogle ScholarPubMed
Park, J.-K., Yang, S.M., Logan, N.C., Hu, Q., Zhu, C., Zarnstorff, M.C., Nazikian, R., Paz-Soldan, C., Jeon, Y.M. & Ko, W.H. 2021 Quasisymmetric optimization of nonaxisymmetry in tokamaks. Phys. Rev. Lett. 126, 125001.CrossRefGoogle ScholarPubMed
Plunk, G.G. 2020 Perturbing an axisymmetric magnetic equilibrium to obtain a quasi-axisymmetric stellarator. J. Plasma Phys. 86 (4), 905860409.CrossRefGoogle Scholar
Plunk, G.G. & Helander, P. 2018 Quasi-axisymmetric magnetic fields: weakly non-axisymmetric case in a vacuum. J. Plasma Phys. 84 (2), 905840205.CrossRefGoogle Scholar
Plunk, G.G., Landreman, M. & Helander, P. 2019 Direct construction of optimized stellarator shapes. Part 3. Omnigenity near the magnetic axis. J. Plasma Phys. 85 (6), 905850602.CrossRefGoogle Scholar
Redi, M.H., Canik, J., Dewar, R.L., Fredrickson, E.D., Cooper, W.A., Johnson, J.L. & Klasky, S. 2001 Localized ballooning modes in compact quasiaxially symmetric stellarators. Tech. Rep PPPL-3578. Princeton Plasma Physics Lab. (PPPL), Princeton, NJ (United States).Google Scholar
Redi, M.H., Johnson, J.L., Klasky, S., Canik, J., Dewar, R.L. & Cooper, W.A. 2002 Anderson localization of ballooning modes, quantum chaos and the stability of compact quasiaxially symmetric stellarators. Phys. Plasmas 9 (5), 19901996.CrossRefGoogle Scholar
Rodriguez, E. 2022 Quasisymmetry. PhD thesis. Princeton University.Google Scholar
Rodríguez, E. & Bhattacharjee, A. 2021 Islands and current singularities in quasisymmetric toroidal plasmas. Phys. Plasmas 28 (9), 092506.CrossRefGoogle Scholar
Rodríguez, E., Helander, P. & Bhattacharjee, A. 2020 Necessary and sufficient conditions for quasisymmetry. Phys. Plasmas 27 (6), 062501.CrossRefGoogle Scholar
Rodríguez, E., Sengupta, W. & Bhattacharjee, A. 2023 Constructing the space of quasisymmetric stellarators through near-axis expansion. Plasma Phys. Control. Fusion 65 (9), 095004.CrossRefGoogle Scholar
Skovoroda, A.A. 2007 Relationship between the shape of equilibrium magnetic surfaces and the magnetic field strength. Plasma Phys. Rep. 33, 523534.CrossRefGoogle Scholar
Spong, D.A. 2005 Generation and damping of neoclassical plasma flows in stellarators. Phys. Plasmas 12 (5), 056114.CrossRefGoogle Scholar
Strauss, H.R. 1980 Stellarator equations of motion. Plasma Phys. 22 (7), 733.CrossRefGoogle Scholar
Strauss, H.R. 1997 Reduced MHD in nearly potential magnetic fields. J. Plasma Phys. 57 (1), 8387.CrossRefGoogle Scholar
Strauss, H.R. & Monticello, D.A. 1981 Limiting beta of stellarators with no net current. Phys. Fluids 24 (6), 11481155.CrossRefGoogle Scholar
Wakatani, M. 1998 Stellarator and Heliotron Devices, vol. 95. Oxford University Press.CrossRefGoogle Scholar
Zarnstorff, M.C., Berry, L.A., Brooks, A., Fredrickson, E., Fu, G.-Y., Hirshman, S., Hudson, S., Ku, L.-P., Lazarus, E., Mikkelsen, D., et al. Physics of the compact advanced stellarator NCSX. Plasma Phys. Control. Fusion 43 (12A), A237A249.CrossRefGoogle Scholar
Figure 0

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

Figure 1

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 2

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

Figure 3

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

Figure 4

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 5

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

Figure 6

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.

Figure 7

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.

Figure 8

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 9

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.

Figure 10

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.