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

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.


Introduction
An essential requirement of magnetic confinement is that the plasma must be in a stable magnetohydrodynamic (MHD) equilibrium.Since the confining magnetic field in a stellarator is mainly produced by external current-carrying coils or magnets, it does not suffer from the current-driven instabilities intrinsic to tokamaks.For a stellarator, a set of nested toroidal flux surfaces with a common magnetic axis such that the magnetic field lines are tangential to the surfaces is desired.However, the lack of continuous symmetry in a stellarator can lead to a breakdown of nested flux surfaces and the formation of island chains and ergodic field-line regions.Without symmetry, the trapped particles can drift out radially and be lost.Both of these effects are severely detrimental to confinement.Furthermore, the lack of symmetry and the nonlinear nature of ideal MHD equations seriously complicate our understanding of stellarator equilibrium.Consequently, compared with tokamaks, MHD equilibrium theory for generic stellarators still has major unanswered questions.
A hidden symmetry of the magnetic field strength, called quasisymmetry (QS), can ameliorate some of the difficulties mentioned above.While the magnetic field remains fully three-dimensional (3-D), QS requires that the field strength, B, be independent of one of the angular coordinates.The continuous symmetry of B then leads to a conserved canonical angular momentum (due to Noether's theorem), ensuring neoclassical particle confinement of the same quality as that of a tokamak (Helander 2014;Landreman & Catto 2010).Moreover, exact QS1 leads to nested flux surfaces (Burby et al. 2020;Rodríguez, Helander & Bhattacharjee 2020).Indeed, one of the benefits of understanding the space of near-axisymmetric QS is that it provides a valuable perspective for 3-D error-field correction of tokamaks (Park et al. 2021), which are known to be sensitive to 3-D perturbations (Park, Boozer & Glasser 2007a;Park et al. 2009).If the 3-D perturbations of an axisymmetric tokamak restore QS, they can help mitigate harmful transport effects (Park et al. 2007b(Park et al. , 2018)).
Unfortunately, exact QS is difficult, if not impossible, to obtain.Analytical results (Garren & Boozer 1991b, a;Plunk & Helander 2018;Landreman & Sengupta 2019;Jorge, Sengupta & Landreman 2020) seem to indicate that imposing QS on ideal MHD with scalar pressure leads in general to an overdetermined problem, with no global (or volumetric) solution.In particular, formal series expansions carried out in the distance from the axis, the so-called near-axis expansion (NAE), show that, beyond the second order, there are more equations than unknowns.A large-scale numerical optimization approach has successfully provided multiple practical designs (Anderson et al. 1995;Zarnstorff et al. 2001;Najmabadi et al. 2008;Ku & Boozer 2010;Bader et al. 2019;Landreman & Paul 2022).However, the many degrees of freedom arising from the three-dimensionality of stellarators lead to computational challenges in optimizing them.A well-known problem for multi-parameter optimization is getting stuck in local minima in parameter space, and the deviations from exact QS cannot generally be made arbitrarily small.
Significant analytical and numerical progress has been made in recent times in understanding the QS constraint, which has been achieved with very high precision (Landreman & Paul 2022).Near-axis expansions have provided helpful analytical insights into new previously unknown configurations and provided excellent initial guesses for further numerical optimization.It has also been possible to map out the entire quasisymmetric phase space using second-order NAEs.The important question of how the magnetic field strength shapes magnetic flux surfaces can be addressed within the NAE framework.For a local equilibrium, one is free to prescribe the flux-surface shape and the magnetic field strength (Boozer 2002;Skovoroda 2007).However, for a global quasisymmetric equilibrium, the flux-surface shapes are highly constrained (Jorge et al. 2020;Rodriguez 2022).Second-order NAE theory allows one to understand and explore this relationship (Rodriguez 2022;Rodríguez, Sengupta & Bhattacharjee 2023).
Despite the several advantages of NAE, a significant drawback, when applied to QS, is that physical quantities can be approximated only by low-order polynomials in the radial variable.The reason lies in the overdetermined nature of QS.Therefore, only plasma profile functions, such as the pressure, magnetic shear and bootstrap current that can be sufficiently well described by low-order polynomials, can be treated within the NAE framework for QS.A more serious drawback is that the overdetermination problem (discussed above) becomes an obstacle to calculating global magnetic shear, which shows up at the order at which the NAE approach breaks down.In the absence of reliable information on the global magnetic shear, MHD stability analyses, such as those pertaining to ballooning or interchange modes, become unreliable.Therefore, an alternative to the NAE is needed for the equilibrium and stability analysis of quasisymmetric devices.
In this work, we present an alternative approach to QS, which provides analytical insight akin to NAE but allows us to compute approximate global equilibrium solutions with more general profile functions.The expansion parameter in our reduced MHD model is the ratio of the equilibrium perpendicular and parallel length scales (L ⊥ /L ) 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 ⊥ /L ).We further restrict the lowest-order magnetic field to be a purely toroidal vacuum field.With this restriction, we can only treat quasi-axisymmetric configurations.Plunk & Helander (2018) have shown that exact volumetric vacuum QS cannot be obtained close to vacuum axisymmetry under quite general conditions.Without any contradiction with Plunk & Helander (2018), we are able to realize approximate volumetric QS in our model since we satisfy the ideal MHD force balance and QS only to the lowest non-trivial orders.We are also consistent with the CDG (Constantin-Drivas-Ginzberg) theorem (Constantin, Drivas & Ginsberg 2021), which states that approximate volumetric QS can be obtained if the force-balance condition in ideal MHD is modified by the addition of a small force, allowing one to satisfy two of the weak QS conditions exactly.We have imposed the near-axisymmetry restriction to leverage the simplicity of the lowest-order axisymmetric geometry but shall relax this restriction in a future study.
The reduced MHD structure of our model coincides with the traditional large-aspectratio expansion (LAE) of ideal MHD (Hazeltine & Meiss 2003).An LAE model for stellarators, accurate to the first order in the inverse aspect ratio, with the plasma beta ordered as the inverse aspect ratio, is known in the literature (Wakatani 1998;Freidberg 2014) as the high-beta stellarator (HBS) model.(The rotational transform is finite in this model.)Unlike earlier large-aspect-ratio stellarator models (Greene & Johnson 1961;Strauss 1980;Strauss & Monticello 1981), the HBS does not require a large toroidal periodicity mode number, N. Given the high interest in quasisymmetric stellarators, we impose the QS constraint on the HBS and assess some of the important consequences for MHD equilibrium and stability.
The structure of the paper is as follows.In § 2, after discussing the implications of QS for ideal MHD equilibria, we provide a derivation of the quasisymmetric Grad-Shafranov equation (QS-GSE) in reduced MHD.We then discuss the quasisymmetric HBS (QS-HBS) model in § 3 and, in particular, we derive a special form of the QS-GSE.We present an analytic solution to the QS-GSE and its numerical verification in § 4. We discuss ballooning and interchange stability of the analytical equilibria in § 6.We conclude with a discussion of the implications of our results and possible generalizations in § 7.

Ideal magnetohydrodynamics with the quasisymmetry constraint
We begin with the well-known model for a plasma equilibrium, the ideal MHD equations Here, B, 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, μ 0 = 1.We will assume that pressure is constant on a set of nested toroidal flux surfaces labelled by ψ, i.e.
We shall further assume that ∇ψ is non-zero almost everywhere except on a finite set of closed field lines that include the magnetic axis and the separatrix.
The QS constraint can be imposed on the ideal MHD equilibrium in many different but equivalent ways.Here, we shall make use of the following form, which is most commonly referred to as the two-term form (Helander 2014;Burby et al. 2020;Rodríguez et al. 2020): 3) The flux function ψ denotes the toroidal flux, B is the magnetic field strength and F(ψ) is a flux function that can be shown (Helander 2014) to be related to the rotational transform . (2.4) 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 B along the constant ϑ (poloidal) and constant ϕ (toroidal) angles.These angles can be chosen as the Boozer angles for convenience.Thus, Another relevant quantity of interest is the QS vector u (Burby et al. 2020;Rodríguez et al. 2020), which may be defined by the following two equivalent relations: 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 u • B = F(ψ).The QS vector u defines curves that lie on constant flux surfaces along which B does not change since u • ∇B = 0 and u • ∇ψ = 0. Furthermore, 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 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) 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) 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 which satisfies the Newcomb constraint.This also allows us to solve for j || up to a flux function H (ψ) (2.11) We note here that the pressure-driven term in (2.11) constitutes the Pfirsch-Schluter current, whereas the flux function, H (ψ), encompasses all other non-pressure-driven current sources.Since a B • ∇ has been lifted to obtain j || from (2.8), singular Dirac-delta currents (Loizu et al. 2015;Rodríguez & Bhattacharjee 2021;Huang et al. 2022Huang et al. , 2023) ) are possible as well on rational surfaces.We plan to address the singular currents in the future.
In the case of axisymmetry, j || can be written as an elliptic Laplacian-like operator acting on the poloidal flux, leading to a Grad-Shafranov (GS) equation.We would now like to derive a QS-GS equation in analogy with axisymmetry.The general QS-GSE (Burby et al. 2020) (discussed in Appendix B) being too complicated, we shall use the reduced model for ideal MHD derived by Strauss (1997) in the remainder of this section.Strauss' model assumes that the magnetic field is approximately given by where B v is a vacuum field, A is an O(L ⊥ /L ) function that is analogous to the poloidal flux in the axisymmetric case.The fundamental assumption is that the gradients along B v are smaller than those perpendicular to it.The form (2.12a,b) ensures that B is divergence free to first order in L ⊥ /L .We shall assume that plasma beta is first order in L ⊥ /L .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: 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 B v ∝ ∇φ, B 2 v ∝ 1/R 2 , A ∝ −ψ in standard axisymmetric cylindrical coordinates.Thus, the (F/B 2 )p , H terms reduce to the standard R 2 p , II terms in the axisymmetric Grad-Shafranov equation.

The quasisymmetric high-beta stellarator model
In the previous section, we have shown that perfect QS helps to integrate out one hyperbolic characteristic of ideal MHD (Grad 1971).Therefore, the Hamada condition (Helander 2014) is automatically satisfied, which leads to a Pfirsch-Schluter current that is non-singular on rational flux surfaces.However, there remains another hyperbolic characteristic of ideal MHD that is related to B • ∇ψ = 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.

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, φ), with unit vectors (e R , e Z , e φ ) The inverse aspect ratio , defined as the ratio of the minor radius r 0 , and the major radius R 0 , is our expansion parameter, i.e.
However, assuming a purely toroidal vacuum field to lowest order implies that the magnetic axis will remain close to a planar ring whose normal does not rotate around itself.Therefore, we can only access QA for which N = 0 (Landreman & Sengupta 2019;Rodríguez et al. 2023) and Following Freidberg (2014), we now expand the various physical quantities in formal power series of , assuming B 0 to be a constant, and the plasma beta and currents are first order in .All quantities such as 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 where (x, y) are order-unity coordinates along e R and e Z .Since the lowest-order magnetic field is assumed to be a toroidal vacuum field, G ∼ R 0 B 0 , and The magnetic field is thus divergence free and curl free to O( ).Moreover, 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( ) and B • ∇ ∼ , the QS condition only imposes a constraint at O( 2 ) and higher.Except for the treatment of QS, all other details of the expansion are given in Freidberg ( 2014), so we only include the essential results here.We note that one must be careful with the operator ∇ φ as defined in (3.7a-c) since it generates x dependent terms in higher orders to account for the 1/R term in ∇φ.
Proceeding to O( ), we find It can be shown Freidberg (2014) that the divergence-free condition, (3.8a), leads to where the A 1 term defines a poloidal magnetic field with A 1 being the streamfunction.The b φ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) Taking into account the orthogonality of the two terms in (3.9), it follows that the field strength B ≈ B 0 + B 1 , where Finally, taking the curl of B to first order, we obtain the current As discussed in the previous section, the parallel component of J 1 can be determined through a magnetic differential equation of the form (2.8).This requires us to go to second order in .
From the second order, we find the two fundamental HBS equations where we have used The Poisson bracket appearing in (3.14a,b) is defined in the usual way as (3.15) 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(ψ) and F(ψ) from their definitions (2.5a,b) and (3.3) as The QS condition (2.3) then implies that It follows from the expression (3.16a,b) of F(ψ) that, to the lowest order in , ψ F is equal to the poloidal flux 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 Lifting the partial y derivative in (3.20), we obtain Equations (3.13) and (3.21) constitute our quasisymmetric HBS system.To proceed further, it is convenient to introduce normalized variables (β, Ψ, A) such that In terms of the normalized variables, together with the definitions the QS-HBS model reads Once these equations are solved, we can calculate B, J to O( ) through This completes our derivation of the basic equations that govern near-axisymmetric quasisymmetric reduced MHD.These are nonlinear equations for (Ψ, A), with single valuedness of (Ψ, A) as boundary conditions in the (x, y, φ) coordinates.The pressure term β = p 1 /B 2 0 is a free input function.The function a(x, φ) 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 Ψ .

Derivation of the quasisymmetric Grad-Shafranov equation
The QS-GSE was obtained in full generality in Burby et al. (2020).However, it is hard to use and somewhat opaque due to the complexities arising from the 3-D nature of the metric coefficients and the additional consistency conditions.The LAE allows us to cut through these complications and reduce the QS-GSE to a simple form, thus highlighting the essential role of geometry and QS.
We begin by noting that the definition of ∇ || (3.23a,b), and the relation between Ψ and A (3.24c) yields the identity 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 To obtain a single nonlinear equation for Ψ we now eliminate A from (3.27) using (3.24c).We then obtain the QS-GSE Here and elsewhere, we use the notation f ,x to denote the x derivative of f .
Next, we substitute A = a − Ψ in the ∇ || Ψ = 0 equation to find Using the method of characteristics, we can obtain the following exact solution to (3.29) in the form of a travelling wave (TW) To understand the physical meaning of the quantity ξ , we now construct the QS vector u as defined in (2.6a,b).As shown in Appendix B From (3.31a,b) we find that u • ∇ξ = u • ∇x = 0. Thus, the symmetry vector u lies on surfaces of constant ξ and x, thereby ensuring that u • ∇Ψ (x, ξ) = 0. Therefore, ξ denotes one of the Clebsch variables for u, the other being x.Moreover, since the symmetry lines must close on themselves on a torus, ξ must be periodic in φ.
We now look into the relation of the QS vector, u, with the killing vectors of 3-D Euclidean space that generate rotations and translations (Burby et al. 2020).Following Burby et al. (2020), we define which is equivalent to choosing the poloidal flux instead of toroidal flux in the expression for u (2.6a,b).Comparing the symmetry vectors corresponding to axisymmetry, helical symmetry and the QA symmetry in the HBS we find that they have the form of a linear combination of a rotation in φ and a translation in Z.We note that, in a straight cylinder with helical symmetry, it is customary to use the poloidal angle ϑ to denote the angle of rotation, and φ is along Z.We have chosen not to do so in order to express the close relationship between these vectors.From (3.33a-c), we observe that the translation in Z is zero for axisymmetry, constant l for helical symmetry, and periodic (with zero average) in HBS.As shown in Appendix B, although the axisymmetry and helical symmetry vectors are killing, the HBS symmetry vector is not killing for a generic φ dependent ∂ x a(x, φ).Thus, the QA symmetry of the HBS model is indeed a hidden symmetry and not an isometry of the 3-D Euclidean space.

Quasisymmetric Grad-Shafranov equation and consistency conditions
In the case of axisymmetry, a = 0, which implies Ψ ,φ = 0.The QS-GSE (3.28) then reduces to the axisymmetric GSE as expected.However, in the non-symmetric case with a = 0, it is not obvious that the QS-GSE equation (3.28), and the ∇ || Ψ = 0 condition (3.29), are consistent.The conflict lies in the fact that the non-symmetric terms can enter Ψ only through 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 a(x, φ).

As shown in Appendix
where ā(φ), Y(φ) are periodic functions of φ.Note that, in axisymmetry (∂ φ = 0), a must be of the form 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, a can be set to zero in the axisymmetric limit with no loss of generality.
The rather strong restriction on a(x, φ), given by (3.34), implies that we need to solve the following equations for Ψ , subject to single valuedness of Ψ and its derivatives: Thus, we have a tokamak-like GSE subject to a TW deformation.The profile functions β(Ψ ), H(Ψ ) 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 → y + Y(φ) = ξ .Alternatively, one can shift from (x, y, φ) coordinates to (x, ξ, φ) coordinates, such that The details of the transformation are provided in Appendix C. Since the QS-GSE equation (3.36a)only contains (x, y) derivatives and a ,x is only a function of φ, the transformed QS-GSE reads Let us now compare our system with similar systems discussed in Landreman & Sengupta (2018), Rodríguez et al. (2023), Plunk & Helander (2018), Plunk (2020) and Burby et al. (2020).An important consequence of the TW nature (3.36b) of Ψ is that the quasisymmetric deformation to the axisymmetric equilibrium is a periodic displacement purely in the vertical e Z direction.Thus, the system (3.36)does not have a rotating ellipse-like solution near the axis (Landreman & Sengupta 2018).It can be shown self-consistently (Rodríguez et al. 2023) that, as one approaches axisymmetry, the function σ (φ), which controls the rotation of the ellipse, goes to zero in agreement with our results.The non-rotating aspect is also markedly different from Plunk & Helander (2018) and Plunk (2020), where the quasisymmetric perturbations have a harmonic exp (iNφ) phase factor (with N being an integer) in both e R and e Z directions.
A detailed comparison of (3.36) with the general QS-GSE system of Burby et al. (2020) is carried out in Appendix B. To summarize, we find that in the large-aspect-ratio limit, we recover the generalized Grad-Shafranov equations for QS derived by Burby et al. (2020).In general geometry, it is not enough to solve the GSE, and Burby et al. (2020) had to impose several extra compatibility constraints.On the other hand, the only constraint we had to impose to ensure compatibility is ∂ 2 x a(x, φ) = 0. (The extra constraints presumably appear at higher order, which is beyond the scope of the current work.)3.4.Quasisymmetric high-beta stellarator as an integrable system A magnetic field that satisfies MHD equilibrium conditions and possesses nested pressure and magnetic flux surfaces is, in principle, integrable (in the sense described precisely in Burby, Duignan & Meiss 2021).Implicitly, Hamada's condition is assumed to be satisfied by any MHD equilibrium with nested surfaces (Helander 2014).However, the integrability of a model obtained through a formal asymptotic expansion of 3-D ideal MHD equations is not always guaranteed.In particular, Freidberg's HBS model, which consists of two nonlinear coupled PDEs for A, Ψ , 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, ξ, φ) coordinates to construct action-angle coordinates.
First, we note that the definition of ∇ || in (3.37a-c) implies that the pair (x, ξ) satisfies Moreover, the 2-D nature of Ψ is explicit in the (x, ξ, φ) coordinates since We can interpret the ∇ || operator as a total derivative with respect to φ along the magnetic field line.Therefore, we can cast (3.39a,b) as Hamilton's equations of motion where ξ is the coordinate, x is the conjugate momentum and Ψ (x, ξ) is the Hamiltonian.The condition (3.40) implies that the Hamiltonian is independent of time.
Next, we perform a canonical transform from (x, ξ, φ) to (I, ϑ, φ) such that Requiring that Ψ be independent of ϑ, we arrive at the action-angle coordinates, where It then follows that, to this order of accuracy, the action I is nothing but the toroidal flux ψ, ϑ is the canonical straight-field-line poloidal angle and ι is the rotational transform The following more explicit form of ι can be derived from (3.43a,b) using Ψ ,x x ,Ψ = 1: (3.44) 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 ξ .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.

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, (Ψ, α, ) where α is the field-line label and is the arclength along B.
The expression for the field-line label α can be derived from the condition B • ∇α = 0, which implies To simplify (3.45), we further change coordinates to (Ψ, ξ, φ).We have provided the necessary details of the transformation to the straight-field-line coordinates in Appendix C.
In (Ψ, ξ, φ) coordinates, (3.45) takes the simplified form whose solution can be immediately obtained as (3.47) Here, we have used the fact that Ψ is independent of φ at a fixed ξ .We can easily check that α = φ − ι −1 ϑ from the Hamilton's equations for the field lines We note that there is a sign difference between (3.49) and what is typically used in the stellarator literature (D'haeseleer et al. 2012;Helander 2014) owing to the definition of α in (3.47).However, the benefit of this form of α and B is that it smoothly goes over to the tokamak expressions when the deviation from axisymmetry is negligible.To that end, we can recast α as It is important to clarify that (B φ /B 0 − 1) is only a function of x(Ψ, ξ ) and β(Ψ ), and independent of φ.Through the canonical transform (3.43a,b), we then obtain b φ = b φ (Ψ, ϑ).Furthermore, since b φ = 1 + O( ), α = φ − q(Ψ )ϑ + O( ).Thus, ϑ is indeed the straight-field-line poloidal angle as discussed in § 3.4.
Similarly, the expression for can be derived from the condition B • ∇ = B, or equivalently in (Ψ, ξ, φ) coordinates (3.53) The −1 factor in the expression for 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 (3.53) also follows from the fact that in (Ψ, α, φ) coordinates ∇ || = ∂ φ .The function L(Ψ ) is a homogeneous solution of the operator ∇ || .

Analytic solutions of the QS-GSE: extended Soloviev profiles
We now have in our possession the QS-HBS model, which ensures both QS and force balance to first order in the LAE.As discussed in the previous section, we can start with any LAE tokamak equilibrium and deform it to obtain a quasiaxisymmetric stellarator.In this section, we show a simple analytic exact solution to the QS-HBS model that goes beyond the scope of the NAE.Using the VMEC code (Hirshman and Whitson 1983), we can verify that this solution approximates 3-D MHD equilibrium with good volumetric QA.In particular, we demonstrate that as the aspect ratio becomes large, the QS error and the differences between the numerical 3-D equilibrium from VMEC and the asymptotic analytical model tend to zero.We shall now look for a class of LAE MHD equilibria with the following profile functions (β(Ψ ), H(Ψ )): such that the QS-GSE takes the following form of a linear equation in Ψ : Here, β 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 Here, J n are order-n Bessel functions of the first kind.We combine this solution with a deformation Y(φ) = −0.5 sin 2φ, where each poloidal plane is rigidly displaced in the vertical direction by −Y(φ).
We then use the above Ψ 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/ .The values of the other parameters were held constant at H 2 = 4.8 and β 1 = 16π × 10 −2 .Note that β(Ψ ) here is normalized to be O(1); the true ratio of hydrodynamic to magnetic pressure is β(Ψ ).To convert normalized quantities to SI units, we used the normalization factors B 0 = 1 T and FIGURE 1.The deformed equilibrium with a circular cross-section and aspect ratio 5. r 0 = 1 m.This choice of parameters corresponds to holding the parameters r 0 , B 0 , p 1 , H 1,SI and H 2,SI constant when the QS-GSE is written in the SI system 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 max Bn,m (ψ) 2 , (4.5) and is plotted in figure 3 as a function of inverse aspect ratio .Here, Bn,m (ψ) is a Fourier mode of B = |B| on flux surface ψ.As can be seen, it scales as 2 , which is expected since the QS-HBS model is derived at order .

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 ι 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 ι 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(φ) = −0.5bsin 2φ 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(φ).The maximum QS error is plotted in figure 6 and again scales as 2 , but is larger by a factor of ∼2.3 than in the circular cross-section case.

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 β of O( ), consistent with the high-β stellarator model.Since high-β 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.FIGURE 5. Flux surfaces from the axisymmetric equilibrium (solid blue) vs from the deformed equilibrium (dashed black) at the φ = 0 poloidal plane for aspect ratios 5 (a) and 10 (b).FIGURE 6.The maximum QS error (black dots) again scales as 2 (dashed blue line is 2.3 2 ).
(b) (a) 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.
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. (B • ∇B) • ∇p > 0, which tends to destabilize the plasma.
First, we calculate the Mercier (or interchange) stability quantified by the quantity D Merc from VMEC.The Mercier stability criterion is where the quantity V † † is related to the magnetic well (Greene 1997).The expression for V † † contains the contributions from different terms that arise from the magnetic well, geodesic curvature, plasma current and magnetic shear.A positive D Merc corresponds to stability against interchange modes, whereas a negative D Merc implies instability.
The radial variation of D 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.
As the equilibria analysed are stable against interchange modes for most of the volume, we then analyse their stability against the infinite-n, ideal ballooning mode.The ideal ballooning mode is another curvature-driven instability that can appear in equilibria with finite magnetic shear.To calculate the stability, we numerically solve (Gaur et al. 2023) the ballooning eigenmode equation (Connor, Hastie & Taylor 1978;Dewar & Glasser 1983) . Ballooning eigenvalue (λ) contour plots at ρ = 0.2 and ρ = 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 α = 96, n θ 0 = 96 points to accurately calculate the maximum growth rate.This process is repeated for all flux surfaces.
where λ is the ballooning eigenvalue, and g, c, f are the following normalized geometric coefficients: where α = φ − 1/ι(θ − θ 0 ) is the field-line label, θ is the PEST straight-field-line angle, φ is the cylindrical toroidal angle and θ 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ψ/a 2 N , where ψ is the toroidal flux enclosed by the boundary (without the factor of 2π).We solve (6.2) for the eigenvalue λ = γ 2 a 2 N /B 2 N , where γ a N /B N is the normalized growth rate, and the ballooning eigenfunction X.If λ > 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 θ 0 = 96 uniformly spaced values of θ 0 ∈ [−π/2, π/2] and n α = 1 on the field line α = 0.For perturbed tokamaks, we repeat the same process on n α = 96 field lines with uniformly spaced values of α ∈ [−π, π].A typical plot of the growth rate for the perturbed ITER-like case is shown in figure 8.
In figure 9, we present the ideal ballooning stability eigenvalue λ 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 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.
(b) (a) 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.the circular tokamak case.Thus, strongly shaped QS-HBS equilibria may retain favourable stability even after the toroidal symmetry is broken.
To gain insight into the ballooning stability of the ITER-like equilibria, we calculate the effective ballooning by applying the transformation X = √ gX to (6.2).The transformed ballooning equation becomes where 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 X.We can see that the stellarator equilibrium potentials have higher peaks compared with tokamaks, which cause the ballooning eigenfunction to decay rapidly, similar to the decay of a wavefunction in a potential well in Schrödinger's equation.Towards the edge, we also observe more oscillations in stellarator potentials as seen in figure 11(a) due to their 3-D shape.This decay of the eigenfunction is similar to the phenomenon of Anderson localization (Anderson 1958).Multiple studies by Redi et al. (2001Redi et al. ( , 2002) ) have demonstrated the importance of Anderson localization of ballooning modes due to perturbations that break axisymmetry of the background equilibrium.Further detailed analysis of the localization of the 3-D modes will be published elsewhere.

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 1 and high plasma beta (β ∼ O( )), are consistent with the HBS model (Freidberg 2014).With a purely toroidal vacuum field and a constant field strength as our lowest-order magnetic field, we show that the quasiaxisymmetric MHD equilibrium can be described by a Grad-Shafranov equation, much like axisymmetry.In particular, we show that approximate quasi-axisymmetric equilibria can be obtained from large-aspect-ratio axisymmetric Grad-Shafranov solutions by breaking axisymmetry through a periodic up-down deformation of the surfaces.Our study of Mercier and ballooning stability of some of these equilibria shows that the non-axisymmetric deformations do not significantly degrade the MHD stability properties.
The original HBS model Freidberg (2014) describes MHD equilibrium in large-aspect-ratio stellarators.The HBS consists of two coupled nonlinear PDEs: magnetic differential equations for the pressure surface and the parallel current.Being fully three-dimensional and nonlinear, they inherit the same problems as 3-D MHD: the non-existence of nested flux surfaces, the possibility of current singularities and the formation of magnetic islands.As discussed in Freidberg (2014), the nonlinearity and three-dimensionality present a severe hindrance to constructing analytical equilibria.An exception is the Greene-Johnson (Greene & Johnson 1961) limit of the HBS model, N ∼ −1/2 1, for a classical stellarator (Wakatani 1998;Freidberg 2014;Loizu et al. 2017;Baillod et al. 2023), where a GS equation can be obtained by an averaging over the short helical wavelength associated with large N.
In this work, we have imposed volumetric QS on the HBS model, allowing us to integrate the magnetic differential equation for the parallel current and obtain a 3-D Grad-Shafranov equation for the flux-surface label without imposing large N.The analysis of the fully 3-Dequilibria parallels that of an axisymmetric tokamak, thanks to the constraints from QS, consistent with the large-aspect-ratio limit of Burby et al. (2020).Since volumetric QS and MHD are not satisfied exactly but only to the lowest non-trivial order in , our results are also consistent with Plunk, Landreman & Helander (2019), Plunk (2020), and the CDG theorem (Constantin et al. 2021).
Our work was motivated by the success of the second-order NAEs in providing an analytic description of quasisymmetric stellarators.While the near-axis approach works for any stellarator, provided we zoom in sufficiently close to the magnetic axis, the present approach requires the stellarator to have a large aspect ratio, which is quite reasonable for most stellarators.A significant advantage of our approach over the NAE is a global description of MHD equilibrium free from any polynomial expansion in the radial variable.Therefore, magnetic shear, currents and pressure profiles can be quite general.Furthermore, we can allow any flux-surface shaping.In contrast, the second-order NAE only allows triangularity in shaping.
The QS-GSE admits several classes of exact analytical solutions.In fact, any analytic solutions of the large-aspect-ratio axisymmetric Grad-Shafranov equation can generate a quasisymmetric solution through a straightforward coordinate transform that induces a periodic up-down deformation.Utilizing this, we have generated numerical solutions starting from an axisymmetric tokamak equilibrium and deforming it accordingly.We have then verified the validity of this methodology using the VMEC code.As the aspect ratio gets larger, the deviations from the periodically up-down shifted tokamak equilibrium, and the actual stellarator equilibrium generated with the same boundary using the VMEC code, become smaller.
There are several directions in which the current work can be extended.First, our current model can only describe quasiaxisymmetric configurations close to axisymmetry due to the choice of the lowest-order vacuum field being a purely toroidal field.A more general choice for the lowest-order vacuum field is needed to describe QA and quasihelical symmetry.Second, the isomorphism between QS-HBS and axisymmetric tokamak equilibrium suggests that we can carry out a local geometry analysis for QS-HBS in close analogy to the axisymmetric Miller-geometry case (Miller et al. 1998).Near-axisymmetric local equilibria could then be used to study kinetic ballooning modes, ion temparature gradient (ITG) and electron temperature gradient (ETG) thresholds.Third, one can study subsonic flows in quasisymmetric stellarators, which have been predicted (Spong 2005) and experimentally demonstrated (Gerhardt et al. 2005), using the QS-HBS model.The QS-HBS model should offer interesting analytical insights into the structure of subsonic flows and their damping.Fourth, the HBS model in the Greene-Johnson limit has been used to study the plasma β-limit in a stellarator, which is of practical importance (Freidberg 2014;Loizu et al. 2017;Baillod et al. 2023).An analogous study can be conducted for the QS-HBS model to understand how QS impacts the β limit.Finally, we hope that the near-axisymmetric aspect of the QS-HBS model can provide an avenue to study non-axisymmetric deformations of tokamaks that preserved QS and hence can be successfully used in controlling tokamaks using 3-D perturbations (Park et al. 2021). (3.36).As a result, we get (3.36) with the identification that H (Ψ ) is the O(1) term of the expression FF .
Finally, let us show that u HBS is not a killing vector.As discussed in Burby et al. (2020), a vector u is killing if, and only if, the vector w is identically zero, where 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, φ), to be advantageous: (Ψ, ξ, φ), (Ψ, α, φ), and (Ψ, α, ). (C1a-d) We begin with the transformation from (x, y, φ) to remits this problem.We need not just α but also ∇α to be accurate to O( ).The first term in α, i.e. φ, has a gradient of O( ).Thus, we also need the O( ) piece from B φ /B 0 to contribute to ∇α to O( ).
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 ( τ d ), representing the torsion contribution to the on-axis rotational transform formula due to Mercier.
Using cylindrical coordinates (R, Z, φ), we can denote the axis by the following space curve: r = R 0 (e R + f(φ)e Z ).Straightforward calculation of the higher derivatives of r leads to the following expressions for κ, τ : (D5a,b) Therefore In general, this is a non-vanishing quantity.However, since the deformation f can be only as large as the minor radius, i.e. f = f 1 , the integrated torsion is not generally an O(1) quantity as we now show.We find that τ d = dφ(f 1 + f 1 ) + 3 dφ(f 1 + f 1 )(f 1 ) 2 + O( 5). (D7) The first term in the expression for torsion averages out identically.The next term is O( 3 ) if the toroidal harmonic, n ∼ 1.However, for sizable n, we can get an amplification factor coming from the derivatives of f.To estimate this factor, we note that the total torsion can be further simplified to since the f 1 term is a total derivative.We observe several facts from (D8).First, if f 1 (φ) is even or simple harmonic, i.e. exp(inφ), the total torsion vanishes.The reason why a simple harmonic f 1 leads to a zero total torsion is that, in this case, we can always add a term n 2 f 1 f 2 1 , which vanishes by itself.However, this term cancels the f 1 f 2 1 term due to the simple-harmonic nature of f 1 .So, consider a two-harmonic stellarator-symmetric deformation: f 1 = a sin nφ + b sin mφ.After some straightforward algebra, one obtains 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 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 ∂ φ ∼ ∇ ⊥ is to be valid, one must have a ∼ 1/n and b ∼ 1/m; and thus ι τ = O( 3 n 2 ).This implies that low-amplitude large-n modes could have a sizable torsional contribution to ι.However, very large n would lead to higher derivatives growing faster, leading to the breakdown of the LAE.Hence, n 3 is perhaps the best choice for n, which limits the torsional contribution to ι τ ∼ 0.05 for an aspect ratio of ∼5.

FIGURE 4 .
FIGURE 4. The deformed equilibrium with an ITER-like cross-section and aspect ratio 4.84.
FIGURE 11.Eigenfunction 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.