2.1. General problem set-up

Let us set up the equilibrium problem by introducing the governing set of equations that the magnetic field $\mathbf {B}$ must satisfy. First, for the field to represent a magnetic field, it must be solenoidal, i.e. (i) $\nabla \cdot \mathbf {B}=0$ . In addition, we shall consider the field to have nested toroidal surfaces labelled by $\psi$ (the toroidal flux enclosed by the flux surfaces over $2\pi$ ) tangent everywhere to the field; that is, (ii) $\mathbf {B}\cdot \nabla \psi =0$ . An appropriate degree of smoothness of this toroidal foliation of space is assumed (especially in the neighbourhood of the field axis, where the near-axis expansion will ensue) (Burby, Duignan & Meiss Reference Burby, Duignan and Meiss2021; Duignan & Meiss Reference Burby, Duignan and Meiss2021). Finally, we must impose the equilibrium condition, which we do in its simplest form (the static limit of MHD (Kruskal & Kulsrud Reference Kruskal and Kulsrud1958; Wesson Reference Wesson2011; Freidberg Reference Freidberg2014)); (iii) $\mathbf {j}\times \mathbf {B}=\nabla p$ , where $\mu _0\mathbf {j}=\nabla \times \mathbf {B}$ is the current density. It is possible to generalise the treatment beyond MHD with isotropic pressure (Rodríguez & Bhattacharjee, Reference Rodríguez and Bhattacharjee2021), but we do not do that here.

We understand an equilibrium solution to be the field $\mathbf {B}$ (and its associated flux surfaces $\psi$ ) as a function of $\mathbf {r}\in \mathbb {R}^3$ , the position in the laboratory frame, that satisfies this set of equations. This makes it natural to directly express $\mathbf {B}=\mathbf {B}(\mathbf {r})$ and $\psi =\psi (\mathbf {r})$ , and solve our set of equations directly in this form. This approach is known as the direct-coordinate approach, which was pioneered by Mercier (Reference Mercier1962) and Solov’ev & Shafranov (Reference Solov’ev and Shafranov1970), and subsequently developed and used by various authors (Shafranov & Yurchenko Reference Shafranov and Yurchenko1968; Lortz & Nührenberg, Reference Lortz and Nührenberg1978; Jorge et al. Reference Jorge and Landreman2020 a,b; Duignan & Meiss Reference Burby, Duignan and Meiss2021; Sengupta et al. Reference Sengupta, Rodríguez, Jorge and Bhattacharjee2024). This approach is, however, not ideally suited to dealing with optimised stellarators. The reason is that many physical properties of the field are not a direct consequence of the full field geometry, but rather only $B=|\mathbf {B}|$ (see for example the neoclassical behaviour or single-particle dynamics (Boozer Reference Boozer1983)), and thus an approach that incorporates $|\mathbf {B}|$ more organically is preferred. This brings us to the alternative option, which we favour in this paper, known as the inverse-coordinate approach (Garren & Boozer, Reference Garren and Boozer1991 Reference Garren and Boozerb ; Landreman & Sengupta Reference Landreman and Sengupta2019), in which $\mathbf {B}$ is described as a function of Boozer coordinates $\{\psi, \theta, \varphi \}$ (D’haeseleer et al. Reference D’haeseleer, Hitchon, Callen and Shohet2012; Boozer Reference Boozer1981), where $\theta$ and $\varphi$ are the poloidal and toroidal angles, respectively.

The use of Boozer coordinates is particularly convenient as it enables a succinct representation of $\mathbf {B}$ in a form that ensures the following:

(2.1a) \begin{align} \mathbf{B}=\nabla \psi \times \nabla \theta +\iota (\psi )\nabla \varphi \times \nabla \psi , \end{align}

(2.1b) \begin{align} = G(\psi )\nabla \varphi +I(\psi )\nabla \theta +B_\psi (\psi, \theta, \varphi )\nabla \psi . \\[6pt] \nonumber \end{align}

The first contravariant form introduces into the problem the rotational transform, $\iota$ , while the latter brings the Boozer currents $I$ and $G$ , as well as the covariant $B_\psi$ . The solution of our problem is now a consistent set of functions $\{\iota, G,I,B_\psi \}$ as functions of Boozer coordinates.

However, we are overlooking a key element in the construction. By adopting a representation in Boozer coordinates, the inverse-coordinate approach loses a direct link to real space, the laboratory frame. If one was given the field at some coordinate triplet $\{\psi _1,\theta _1,\varphi _1\}$ , one would actually not know where in space this point is. The solution is thus incomplete. To connect Boozer coordinates to real space we define the following auxiliary functions $\{X,Y,Z\}$ , such that

(2.2) \begin{align} \mathbf {r}=\mathbf {r}_{{axis}}+X\hat {\boldsymbol {\kappa }}+Y\hat {\boldsymbol {\tau }}+Z\hat {\textbf{t}}, \end{align}

where $\mathbf {r}_{\mathrm {axis}}$ describes the magnetic axis of the equilibrium, and the triad $\{\hat {{\textit{t}}},\hat {\boldsymbol {\kappa }}, \hat {\boldsymbol {\tau }}\}$ is the Frenet–Serret (FS) triad of the axis (Frenet Reference Frenet1852; Animov Reference Animov2001). That is, the unit tangent, normal and binormal vectors, respectively. The complete description of the equilibrium requires solving for these space functions $\{X,Y,Z\}$ alongside the field; it will then be through (2.2) that the shape of flux surfaces in space will be described (think of $\mathbf {r}$ at constant $\psi$ ). The functions $\{X,Y,Z\}$ are also directly involved in (2.1). This comes apparent when using the dual relations $\partial _{q_i}\mathbf {r}=\mathcal {J}\epsilon _{ijk}\nabla q_j\times \nabla q_k$ , for ${\textit{q}}$ Boozer coordinates (Garren & Boozer, Reference Garren and Boozer1991 Reference Garren and Boozerb ; Hazeltine & Meiss Reference Hazeltine and Meiss2003; D’haeseleer et al. Reference D’haeseleer, Hitchon, Callen and Shohet2012).

With the problem set-up in this form, we may now turn to the asymptotic treatment. In the context of the near axis, the ordering ‘parameter’ of the asymptotic treatment is some measure of the distance from the magnetic axis, which we will have to construct and shall refer to as $r$ . The near-axis description consists of an expansion of the problem (i.e. the functions and governing equations) in powers of $r$ and the ordered solution of the ensuing hierarchy. Given that this distance $r$ is more a coordinate than an externally imposed constant, one may always choose a distance sufficiently close to the axis in which the asymptotic description holds. This distinguishes the near-axis from other similar asymptotic approaches such as large aspect-ratio expansions (Greene, Johnson & Weimer Reference Greene, Johnson and Weimer1971; Freidberg Reference Freidberg2014; Solov’ev & Shafranov, Reference Solov’ev and Shafranov1970; Cowley et al. Reference Cowley, Kaw, Kelly and Kulsrud1991). With the closed axis enclosing zero toroidal flux $\psi =0$ , it is natural to choose as a pseudoradial coordinate $r=\sqrt {2\psi /\bar {B}}$ , where $\bar {B}$ is a normalising reference magnetic field (assuming $\psi \geqslant 0$ ).

2.2. Shaping the magnetic axis

Before explicitly looking at the expansion in $r$ , we must start with the choice of an appropriate reference axis. Its shape is rather important, as it does not only serve as reference, but also provides the basis in (2.2) for the construction of the field. Its FS framing is chosen as a convenient basis (Serret, Reference Serret1851; Frenet Reference Frenet1852; Mathews & Walker Reference Mathews and Walker1964; Animov Reference Animov2001), but we must comment on the subtleties that arise with that choice. The 3-D non-intersecting (simple) space curves with nowhere vanishing tangent (regular) (Moffatt & Ricca Reference Moffatt and Ricca1992; Oberti & Ricca Reference Oberti and Ricca2016) nor curvature (geometric Def. 4.2.4 Fuller Jr, Reference Fuller and Edgar1999 or non-degenerate (Pohl Reference Pohl1968) curves), have a uniquely defined FS frame everywhere along the curve satisfying (Mathews & Walker Reference Mathews and Walker1964).

(2.3) \begin{align} \frac {\mathrm {d}\mathbf {r}_{{axis}}}{\mathrm {d}\ell }=\hat {{\textit{t}}},\quad \frac {\mathrm {d}\hat {{\textit{t}}}}{\mathrm {d}\ell }=\kappa \hat {\boldsymbol {\kappa }},\quad \frac {\mathrm {d}\hat {\boldsymbol {\kappa }}}{\mathrm {d}\ell }=-\kappa \hat {{\textit{t}}} + \tau \hat {\boldsymbol {\tau }},\quad \frac {\mathrm {d}\hat {\boldsymbol {\tau }}}{\mathrm {d}\ell }=-\tau \hat {\boldsymbol {\kappa }}, \end{align}

where $\kappa$ and $\tau$ are the curvature and torsion of the curve, and $\ell$ is the length along it (with all $\hat {\cdot }$ vectors being normalised). All quantities on the right-hand side are generally functions of $\ell$ .

Although this set of space curves exhausts the possibilities for QS fields (Garren & Boozer, Reference Garren and Boozer1991 Reference Garren and Boozera ; Landreman & Sengupta Reference Landreman and Sengupta2019; Rodríguez et al. Reference Rodríguez2022), the construction of equilibria with poloidally closed $|\mathbf {B}|$ contours requires that curvature vanishes somewhere along the axis. The reason is that any amount of field line bending (in this case the magnetic axis) necessarily supports a finite gradient of $|\mathbf {B}|$ normal to the axis, $\nabla _\perp \ln B=\kappa \hat {\boldsymbol {\kappa }}$ , and thus some non-zero poloidal variation (2.21.5, Wesson (Reference Wesson2011)). We must then specialise on degenerate curves with vanishing curvature points, which may be called flattening points and shall take to be smooth. Although we shall assume such points to be isolated, they can be more or less flat, depending on how many $\ell$ -derivatives of $\kappa$ vanish there. The order of the first non-vanishing derivative is called the order of the zero. We shall additionally specialise on flattening points that are also points of stellarator symmetry (Dewar & Hudson Reference Dewar and Hudson1998), which simplifies the treatment significantly.

In such curves, the FS frame is discontinuous at flattening points. If the order of the zero is even, though, one may show that the frame has a removable discontinuity there. If it is instead of odd order, the frame undergoes a $\pi$ -rotation about the axis across the point. To form a continuous frame, then, a signed frame needs to be defined (the $\beta$ -frame of Carroll, Köse & Sterling (Reference Carroll, Köse and Sterling2013) and Plunk et al. (Reference Plunk, Landreman and Helander2019)). The frame is constructed by solving the regular FS equations, (2.3), starting from some non-degenerate point on the curve and moving along the curve, flipping the FS frame every time a zero of odd order is traversed. The result is a smooth, continuous frame in $\ell \in [0,L)$ , where $L$ is the total length of the axis and the binormal and normal (as well as the curvature) are $\pm$ their FS form. Note that this signed frame still satisfy the FS equations (so (2.2) and (2.3) can be interpreted to involve the signed frame), and thus the near-axis construction in its usual form, as in LS, may be used. We will consider the signed frame throughout this paper. See appendix A for details and proofs of these statements.

Although continuity and differentiability of the frame are guaranteed as we move along the axis, because the curve is closed, this signed frame is not guaranteed to be periodic. In fact, if the sum of the order of all the flattening points is odd, then the frame at 0 and $L$ will have a flip respect to each other. In that case, we refer to the axis as a half-helicity curve, reminding ourselves the definition of the helicity as ‘the number of times the normal to the axis encircles the axis in a complete toroidal turn’ (see appendix A for more precise definitions and also Camacho-Mata & Plunk (Reference Camacho-Mata and Plunk2023)). In practice, for a field that has a single magnetic well per field period and is QI, the curve will have two distinct flattening points: one where $B_0=B_{{min}}$ and another where $B_0=B_{{max}}$ (see § 3 for the specific implications of QI that lead to this). The flattening point at $B_{{min}}$ must be odd by omnigeneity, which makes the order of the top decide the helicity of the axis: half-helicity if even, integer helicity if odd. Because the frame itself is non-periodic in the half-helicity case, the magnetic field described in (2.2) would be unphysical unless the functions $X$ and $Y$ are what we define as half-periodic functions (see appendix B.1). That is, functions that satisfy $f(\ell )=-f(\ell +L)$ and thus can be understood to be periodic in the extended $\ell \in [0,2L)$ domain. A consistent treatment of these half-periodic axes is possible and explicitly discussed in detail in appendix B.1, where the necessary subtle modifications to the second-order construction are detailed.

The choice of a space curve with appropriate stellarator symmetric flattening points as our magnetic axis is the starting point of the near-axis construction. Some of these curves may be constructed straightforwardly following the prescriptions in Plunk et al. (Reference Plunk, Landreman and Helander2019), Rodríguez et al. (Reference Rodríguez2022), Jorge et al. (Reference Jorge, Plunk, Drevlak, Landreman, Lobsien, Camacho Mata and Helander2022), Camacho Mata et al. (Reference Camacho Mata, Plunk and Jorge2022) and Camacho-Mata & Plunk (Reference Camacho-Mata and Plunk2023). A direct way of doing so by specifying the curvature and the torsion will be presented in a future publication.

Once we have a magnetic axis, we must specify how the magnetic field magnitude varies along it; namely, we must provide a $L$ -periodic function $B_0(\ell )$ . The same way as the field strength in a straight magnetic mirror can be externally curated, the function $B_0(\ell )$ must also be provided to the construction. We shall specialise, for simplicity, on fields with a single distinct trapping well in each field period along the magnetic axis, which implies a unique $B_{{min}}$ and $B_{{max}}$ (with the ratio $\Delta =(B_{{max}}-B_{{min}})/(B_{{max}}+B_{{min}})$ defined as the mirror ratio). As mentioned above, these extremal points must match the flattening points of the axis curvature. This is a necessary consequence of having poloidal contours (‘pseudosymmetry’) (Rodríguez & Plunk, Reference Rodríguez and Plunk2023; Skovoroda Reference Skovoroda2005; Landreman & Catto Reference Landreman and Catto2012); i.e. unless this was true, the radial drift would be non-vanishing at these points, and deeply and barely trapped particles would drift away. Due to stellarator symmetry, $B_0$ must be an even function about these points.

In order to sustain such a magnetic field line, we must thread the magnetic axis, through Ampere’s law, with a poloidal current $G_0$ ((6.6.2), D’haeseleer et al. (Reference D’haeseleer, Hitchon, Callen and Shohet2012)). Finally, and because the axis itself is a magnetic field line, there is a strict connection between the length along the field line $\ell$ and $\varphi$ (the Boozer angle), by virtue of being straight field line coordinates, $\mathrm {d}\ell /\mathrm {d}\varphi =|G_0|/B_0$ (see (A20) in LS).

2.3. First-order near-axis construction

Once we have our axis, we move on to a description of the field in its neighbourhood; that is, we must now explicitly consider the expansion in $r=\sqrt {2\psi /\bar {B}}$ , and look at the leading $O(r)$ parts of the field. Here $\bar {B}$ is a reference normalising magnetic field value, typically some average of $B_0$ . Because of its radial-like nature, to avoid coordinate singularities on axis, the expansion in $r$ requires a careful coupling between powers of $r$ and harmonics of $\theta$ (Mercier Reference Mercier1964; Kuo-Petravic & Boozer Reference Kuo-Petravic and Boozer1987; Landreman & Sengupta Reference Landreman and Sengupta2018). In particular, any function $f$ must take the form $f=\sum _{n=0}^\infty r^n f_n$ and $f_n=\sum _{m=0}^n (f_{nm}^c(\varphi )\cos m\chi +f_{nm}^s(\varphi )\sin m\chi )$ , where the latter sum is over even or odd numbers depending on the parity of $n$ , and $\chi =\theta -M\varphi$ with $2M\in \mathbb {Z}$ the helicity of the signed frame. This defines the subscript notation to be repeatedly utilised throughout the paper.Footnote ¹ Considering lower powers of $r$ then implies keeping a small number of poloidal harmonics, which is the key to the strength of the method. Resolving the resulting equations in $\chi$ -harmonics and powers of $r$ , the problem reduces to a hierarchy of equations on $\varphi$ . The detailed accounts of the formal expansion of (2.1) and equilibrium can be found in many works, in particular LS (originally in Garren & Boozer (Reference Garren and Boozer1991 Reference Garren and Boozerb )). We now present the structure of the problem to first order in $r$ , including the equations that need to be solved, their physical meaning and the inputs needed to complete the description. The solution of the problem to $O(r)$ is summarised in a schematic way in figure 1, which we shall follow closely in the description that follows.

Figure 1. Inverse-coordinate near-axis construction to first order. Diagram depicting the key elements in the near-axis description of an equilibrium to first order, and the sequential order (left to right) in which to proceed. The construction starts with a shape for the magnetic axis and the magnetic field strength along it, both of which are given as inputs. To construct the neighbouring flux surfaces one must then provide the leading-order variation of the normal shaping $X_1$ , which directly relates to the field $B_1$ (this can go both ways). With that, one may then solve a differential equation on an auxiliary function $\sigma$ , which is all that is needed to complete the flux surface description through $Y_1$ . The encircled functions are the functions and parameters involved in the near-axis description, with the label (A xx) denoting the equations from Landreman & Sengupta (Reference Landreman and Sengupta2019) needed to find them.

We already started (leftmost part of figure 1) the near-axis construction by the choice of an appropriate magnetic axis shape and the magnetic field strength on it. Naturally, when moving to $O(r)$ the field gains a finite flux surface build around the axis. Because of the analyticity constraint on the functions $X$ and $Y$ of (2.2), these surfaces must be elliptical in the FS frame. There is a certain amount of freedom available (although limited) in the problem to choose how to shape this elliptical surface (see figure 1). In particular, we must provide the function $X_1$ which specifies the distance along the normal from the flux surface to the axis as a function of $\chi$ and $\varphi$ . It is convenient to define it in terms of two $\phi$ -periodic functions (Garren & Boozer, Reference Garren and Boozer1991 Reference Garren and Boozerb ; Plunk et al. Reference Plunk, Landreman and Helander2019) $\bar {d}(\varphi )\gt 0$ and $\alpha _1(\varphi )$ , such thatFootnote ²

(2.4) \begin{align} X_1=\bar {d}\cos \left (\chi -\alpha _1\right ). \\[-24pt] \nonumber \end{align}

The shaping $X_1$ is of particular physical significance: bringing the surface closer (further) to (from) the axis will (with a fixed gradient of $B$ , from $\nabla _\perp (B^2/2)=B^2\boldsymbol {\kappa }$ ) lead to a larger (lower) magnetic field strength on the surface. Thus, the shaping brings in a direct control of the on-surface variation of $|\mathbf {B}|$ . Within the context of the near axis, this formally translates to $B_1=\kappa B_0X_1$ ((A22) of LS). Note that this formulation is different from the standard one, where the magnetic field magnitude is provided as a control input to the problem (Garren & Boozer, Reference Garren and Boozer1991 Reference Garren and Boozera ; Landreman & Sengupta Reference Landreman and Sengupta2019; Rodríguez et al. Reference Rodríguez and Mackenbach2023; Landreman Reference Landreman2022). We need to do so because of the presence of ‘insensitive’ straight sections in which $B_1$ has an intrinsic form (in particular, $B_1=0$ wherever $\kappa =0$ ). A similar observation will appear at higher orders.

Specifying the distance from the surface to the axis along the normal only constitutes one part in the description of the shape of flux surfaces. To complete it, we must know the behaviour along the binormal as well. Finding that $Y_1$ constitutes the last steps of the first-order near-axis construction (see figure 1). To find it, we must remember that we are describing a flux surface, and thus any cross-section of the surface must be thread by a constant toroidal flux. This means that as the field strength changes in $\varphi$ , so must the cross-sectional area, so that $X_1Y_1\sim 1/B_0$ . This justifies the condition $\bar {d}\gt 0$ , as flux surfaces would become infinitely elongated in the event of $\bar {d}$ ever vanishing. The exact relation between $X_1$ and $Y_1$ is, however, not so simple, because besides the conservation of flux, one must make sure that the resulting surface is consistent with the solenoidal magnetic field that lives on it. The consequence of this careful balance is a strong relation between the shaping of the elliptic flux surfaces (including both the elongation and the rotation), the shape of the axis and the rotational transform, $\iota _0$ .

At a formal level, this careful balancing act reduces to the solution of a first-order nonlinear Riccati ( (Polyanin & Zaitsev Reference Polyanin and Zaitsev2017), § 1.4) differential equation ((A26) in LS) on a periodic auxiliary function $\sigma (\varphi )$ (taken in the stellarator symmetric case to satisfy $\sigma (0)=0$ ) and the rotational transform on axis, $\iota _0$ ,

(2.5) \begin{align} \frac {\mathrm {d}\sigma }{\mathrm {d}\varphi }+\left (\bar {\iota }_0-\frac {\mathrm {d}\alpha _1}{\mathrm {d}\varphi }\right )\left [\left (\bar {d}^2\frac { B_0}{\bar {B}}\right )^2+1+\sigma ^2\right ]-2\bar {d}^2\frac { B_0}{\bar {B}}\frac {\mathrm {d}\ell }{\mathrm {d}\varphi }\left (\frac {I_2}{\bar {B}}-\tau \right )=0. \end{align}

The rotational transform, here $\bar {\iota }_0=\iota _0-M$ , is a scalar parameter that must be chosen to guarantee periodicity of $\sigma$ depending on the toroidal current $I_2$ assumed (often set to zero) (Landreman & Sengupta Reference Landreman and Sengupta2018).Footnote ³ This equation has been explored in detail by other authors (Landreman, Sengupta & Plunk Reference Landreman, Sengupta and Plunk2019; Rodríguez et al. Reference Rodríguez and Mackenbach2023) and $\sigma$ given a geometric interpretation of (roughly) signifying the rotation of the ellipses respect to the FS frame (Camacho Mata et al. Reference Camacho Mata, Plunk and Jorge2022; Rodríguez Reference Rodríguez2023). Once the $\sigma$ -equation is solved, $Y_1$ can be directly found ((A25) in LS),

(2.6) \begin{align} Y_1=\frac {\bar {B}}{\bar {d}B_0}\left [\sin\! (\chi -\alpha _1)+\sigma \cos\! (\chi -\alpha _1)\right ] \end{align},

and the first-order construction is completed.

2.4. Second-order near-axis construction

Much of the second-order construction can be understood as a natural continuation to the first order, where similar physics govern the way to proceed. There is, however, as indicated in figure 2, one significant difference: for the first time in the construction, the description of the field involves the pressure gradient, $p_2$ , directly. Sufficiently close to the axis, the magnetic field is force-free, which is why $p$ was not involved at first order. The fulfilment of pressure balance is guaranteed by accommodating the magnetic field component $B_\psi =\mathbf {B}\cdot \mathbf {e}_\psi$ (Boozer Reference Boozer1981). Formally, this involves the solution of simple first-order ordinary differential equations (ODEs) ((A51)–(A52) in LS), which may be done straightforwardly using periodic, vanishing endpoint boundary conditions. In addition, to balance a finite pressure gradient it is necessary to have a finite net current gradient. A formal statement of that is that $G_2+\iota _0 I_2$ (related to the radial variation of the net poloidal and toroidal currents), can be found algebraically at this point in the construction ((A50) in LS). Thus, there is one free constant in the second-order construction, with either $G_2$ being determined by $p_2$ or vice versa.

Figure 2. Inverse-coordinate near-axis construction at second order. Diagram representing the key elements in the typical second-order near-axis equilibrium construction (left to right) to be taken as continuation of figure 1. The construction at second order starts by introducing the pressure gradient, $p_2$ , explicitly, which $B_\psi$ and the current $G_2$ must balance. The $Z_2$ shaping is then uniquely determined to remain consistent. Providing the $\theta$ -dependent shaping of $X_2$ as input (or the, for poloidal $|\mathbf {B}|$ fields less convenient, alternative $B_2$ ), the rest of the construction is uniquely determined, much like at first order. First one solves two coupled ODEs for the rigid displacement of flux surfaces (i.e. for $X_{20}$ and $Y_{20}$ ), to finally complete the magnetic field and flux surface shaping in the binormal direction. The encircled functions are the key elements needed to describe the field, with the label (A xx) denoting the equations from Landreman & Sengupta (Reference Landreman and Sengupta2019) needed to find these quantities. The small circles on the edge of the functions denote that a derivative of lower-order quantities is required to compute the function. The star, in turn, denotes that the function is a solution of an ODE with a derivative of lower-order quantities in its inhomogeneous term.

Once pressure is included in the problem, we solve for $Z_2$ (see figure 2). Formally, finding $Z_2$ only involves operations of near-axis quantities already known, including their toroidal derivatives (see (A 27–29) in LS). This component describes the non-planar shape of cross-sections of constant $\varphi$ . In that sense, it describes how the Boozer coordinate $\varphi$ changes in the neighbourhood of the magnetic axis for it to remain consistent with the near-axis magnetic field.

With the toroidal coordinate resolved, we now meet a step analogous to that with $X_1$ shaping at first order (see figure 2). We have the freedom to impose some shaping along the normal to the axis at second order. In particular, the second harmonic components $X_{2s}$ and $X_{2c}$ , which affect directly the triangularity and up-down symmetry of the cross-sections (Rodríguez Reference Rodríguez2022, Reference Rodríguez2023). As it occurred at first order, the choice of this distance will, through the curvature, directly affect the corresponding harmonics of the second-order magnetic field strength $B_2$ ((A 35–36) in LS),

(2.7a) \begin{align} B_{2c} = \kappa B_0X_{2c}-\hat {\mathcal {T}}_c+\frac {\left(B_{1c}^{\mathrm {QI}} \right)^2}{2B_0}\cos 2\alpha _1, \\[-24pt] \nonumber \end{align}

(2.7b) \begin{align} B_{2s}=\kappa B_0X_{2s}-\hat {\mathcal {T}}_s+\frac {\left(B_{1c}^{\mathrm {QI}} \right)^2}{2B_0}\sin 2\alpha _1, \\[6pt] \nonumber \end{align}

where $\hat {\mathcal {T}}_c$ and $\hat {\mathcal {T}}_s$ are known functions of first-order quantities (see appendix D), and $B_{1c}^{\mathrm {QI}}=B_0 \bar {d}\kappa$ . As it occurred at first order, at the straight sections of the field, the second-order shaping has no effect on the behaviour of the field (this favours, once again, the choice of $X_2$ harmonics ( $X_{2c}$ and $X_{2s}$ ) as inputs, as opposed to $B_2$ ones). At those points, the field becomes solely determined by the first-order solution. Elsewhere, the first-order solution also has a direct effect (physically from the need to preserve divergencelessness and tangent fieldlines), but present a larger degree of freedom. The choice $X_{2c}=0=X_{2s}$ , which we will refer to as minimal shaping, makes that underlying structure manifest at second order.

From general MHD equilibria, we know that physical effects such as the Shafranov shift should come about here, influenced by the pressure gradient but also, and to a large extent, by the self-consistent shaping of the field ((Wesson Reference Wesson2011), § 3.7). This naturally brings us to the next step in the second-order near-axis construction (see figure 2), which is to self-consistently solve for $X_{20}$ and $Y_{20}$ . These two elements are measures of Shafranov shift in the normal and binormal direction (i.e. the relative rigid displacement of flux surfaces with radius) (Landreman Reference Landreman2021; Rodríguez Reference Rodríguez2022, Reference Rodríguez2023). To find these functions consistently with the rest of pieces of the construction requires the solution of two linear coupled first-order ODEs on $\varphi$ . These may be solved numerically in a straightforward way ((A41)–(A48) in LS), the considerations for half-helicity fields requiring us to work with half-periodic functions. With that solution and the second-order shaping inputs, the construction of the surface at second order is completed algebraically finding $Y_{2s}$ and $Y_{2c}$ to preserve, as we learnt at first order, the constant flux assumption as well as the appropriate shaping of the field lines over them ((A32)–(A33) in LS).

Finally, we must compute the total magnetic field strength on the surface, of which only the second harmonics had been directly found after providing the normal shaping of the flux surfaces. The missing component is $B_{20}$ , i.e. the average $\psi$ -derivative of $|\mathbf {B}|$ on the surface. Of course, we could not compute this before knowing the complete shaping of the surfaces, and that is why its evaluation comes at last. Its construction is algebraic ((A34) in LS), similar to (2.7a)–(2.7 b ). Despite coming last in our solution method, the element $B_{20}$ plays a key role in physics, for instance MHD stability (Landreman & Jorge Reference Landreman and Jorge2020; Rodríguez Reference Rodríguez2023) and particle precession (Rodríguez & Mackenbach, Reference Rodríguez and Mackenbach2023; Rodríguez et al. Reference Rodríguez, Helander and Goodman2024).

This concludes the essential elements of the near-axis equilibrium construction at second order, for a field with poloidally closed $|\mathbf {B}|$ contours. The above holds regardless of whether the axis is half or integer helicity, though the former case does have subtleties related to periodicity and continuity, as discussed in appendix B.1.

4.1. Minimally shaped vacuum configuration

Let us start our benchmark by looking at a simple case which we take from the recent publication of Camacho Mata et al. (Reference Camacho Mata, Plunk and Jorge2022). We consider the $N=2$ field in that paper, which is approximately QI at first order with a standard buffer region with $k=2$ . The axis is defined to have an integer helicity with first-order zeroes of curvature, parametrised by

(4.2a) \begin{align} R\,[\text {m}] &= 1-\frac {1}{17}\cos 2\phi, \\[-6pt] \nonumber \end{align}

(4.2b) \begin{align} Z\,[\text {m}] &= \frac {0.8}{2.04}\sin 2\phi +\frac {0.01}{2.04}\sin 4\phi, \\[9pt] \nonumber \end{align}

where the form of $R$ guarantees the correct behaviour of flattening points (Rodríguez et al. Reference Rodríguez, Sengupta and Bhattacharjee2022; Camacho Mata et al. Reference Camacho Mata, Plunk and Jorge2022). The magnetic field on axis is simply described by $B_0\,[\text {T}]=1+0.15\cos 2\varphi$ , and $\bar {d}=0.73$ is assumed constant. We extend this field into second order following the minimal shaping assumption; that is, $X_{2c}=0=X_{2s}$ . For this choice, the second-order field is a mere extension of the first order, and investigating the properties of this second-order completion may be used in the future as a tool to investigate the properties and suitability of the first-order near-axis construction, all within the near axis framework.

Figure 3. Global equilibrium construction from the second-order near-axis field in § 4.1. The figure presents the second-order near-axis construction based on the $N=2$ field of Camacho Mata et al. (Reference Camacho Mata, Plunk and Jorge2022) presented in § 4.1. (a) Cross-sections at constant cylindrical angle in a half-period and 3-D rendering at $A=10$ . The dotted line traces the magnetic axis, with crosses representing the intersection with the cross-sections. (b) Root-mean-square difference in the second-order near-axis magnetic field between the ideal near-axis description and the finite aspect ratio global equilibrium construction with VMEC. A reference quadratic scaling $\propto 1/A^2$ is given, from which the convergence of the error deviated at larger aspect ratios. (c) Explicit comparison between the different poloidal components of the second-order magnetic field magnitude from the global VMEC equilibrium at different aspect ratios, and the ideal near-axis value (black dotted line).

The construction of the near-axis field and the comparison with the global equilibrium is shown in figure 3. The top left plots show cross-sections and a flux surface of the configuration evaluated at $r\,[\text {m}]=0.1$ , which roughly corresponds to an aspect ratio of $A=1/r=10$ . The bottom plot shows a direct comparison between the different components of $B_2$ as a function of $\varphi$ for a number of aspect ratios, showing that as the aspect ratio $A$ is increased, the agreement becomes better. Even at as low an aspect ratio as $A=5$ the second-order near-axis construction reproduces the behaviour of the equilibrium at least qualitatively. This points towards the significance and physical nature of second-order properties. The top right plot shows the level of agreement as a function of aspect ratio in a more quantitative fashion. It shows the root-mean-square difference between the VMEC and near-axis functions defined in (4.1). Although there is an initial improvement with increasing aspect ratio at small $A$ , this seems to saturate quite rapidly (at approximately $A\sim 20$ ). The reason behind this discrepancy can be traced to the behaviour of $B_{2s}$ in figure 3(c). Near the endpoints of the domain, one can see that the function $B_{2s}$ from the near-axis construction is actually discontinuous. And of course, the global equilibrium solve cannot reproduce such unphysical features. Hence the saturation of the error at a level dictated by the size of the discontinuity.

The origin of said discontinuity had already been discussed in the main text preceding this section: it is the lack of smoothness in the standard buffer. To show that this is responsible for the observed disagreement, we repeat the same numerical exercise, but this time using a smooth buffer with $k=5$ . The results are shown in figure 4, where the previous saturation of the error disappears.

Figure 4. Comparison of second-order field of the equilibrium in figure 3 with smooth buffer regions. The figure presents the comparison of the second-order near-axis field for the same configuration as in figure 3 but for the difference in the analyticity of the construction near the field edges. (a) Root-mean-square difference in the second-order near-axis magnetic field between the ideal near-axis description and the finite aspect ratio global equilibrium construction with VMEC. A reference quadratic scaling $\propto 1/A^2$ is given, from which the convergence of the error deviated at larger aspect ratios. (c) Explicit comparison between the different poloidal components of the second-order magnetic field magnitude from the global VMEC equilibrium at different aspect ratios, and the ideal near-axis value (black dotted line). The agreement near the edge is better than that of figure 3.

The disagreement between the near-axis and global equilibrium scales like $\sim 1/A^2$ (this is the difference in the second-order field $B_2$ , and not $r^2B_2$ ). This quadratic reduction of the error is expected from a construction that is accurate to second order, as the finite build of the near axis is (Landreman & Sengupta Reference Landreman and Sengupta2019). This requires the error to be at least $1/A$ , but because of analyticity of the field near the magnetic axis, the $m=0,2$ harmonics only exhibit even powers and thus the scaling with $1/A^2$ . This behaviour is observed over a wide range of aspect ratios. At the highest values of aspect ratio (in figure 4, $A=316$ ) there appears to be some level of discrepancy as VMEC struggles to converge (note, however, the high demands we are placing on the code, as we require the field to be correct to approximately 1 part in $10^8$ or more).

With this first example we show that: (i) the near-axis construction to second order for QI fields works; (ii) that the non-smooth constructions cannot be faithfully reproduced by global physical equilibria; but (iii) that their qualitative aspects can.

4.2. Shaped, finite plasma- $\beta$ configuration

We now consider a more involved example that includes finite plasma $\beta$ and non-zero explicit second-order shaping. In this case, we base the construction off the $N=3$ in Camacho Mata et al. (Reference Camacho Mata, Plunk and Jorge2022), with a smooth buffer region. The axis is defined approximately by

(4.3a) \begin{align} R [\text {m}] & = 1+ 0.09075485\cos 3\phi -0.02058279 \cos 6\phi -0.01106766\cos 9\phi \nonumber\\ & \qquad -0.001644390e \cos 12\phi, \end{align}

(4.3b) \begin{align} Z [\text {m}] & = 0.36 \sin 3\phi + 0.02\sin 6\phi + 0.01\sin 9\phi, \end{align}

which when implemented numerically one should carefully check to satisfy the vanishing curvature condition at the stellarator symmetric points. The magnetic strength on axis is $B_0\,[\text {T}]=1+0.25\cos 3\varphi$ , and $\bar {d}=0.73$ once again. In this case we choose some arbitrary non-zero values for $X_{2c}[\text {m}^{-1}]=0.1\sin 3\varphi -0.6\sin 6\varphi +0.1\sin 9\varphi$ and $X_{2s}[\text {m}^{-1}]=0.1-0.1\cos 6\varphi +0.6\cos 9\varphi$ . Finally, we put a finite toroidal current $I_2\,[\text {T/m}] = -0.9$ and pressure gradient $p_2\,[\text {Pa/m}^2]=-6\times 10^5$ . The resulting field and comparison with VMEC is presented in figure 5. Once again, the quadratic scaling of the error shows agreement between the near-axis construction and the global equilibrium solve. Note that this time the equilibrium is significantly more shaped than in the case of § 4.1, in part as a consequence of the larger number of field periods (and thus larger toroidal derivatives) but also the explicit shaping introduced at second order. This increased complexity is apparent in the comparison of the different components of $B_2$ , but also on the scaling of the error $\Delta B_{{rms}}$ itself, as the deviations from the ideal scaling occur at lower aspect ratios. This case is an example of one of the struggles in the design of near-axis QI fields, which is to avoid this extreme shaping.

Figure 5. Global equilibrium construction from the second-order near-axis field in § 4.2. The figure presents the second-order near-axis construction based on Camacho Mata et al. (Reference Camacho Mata, Plunk and Jorge2022) presented in § 4.1 including finite $\beta$ , toroidal current and additional second-order shaping. (a) Cross-sections at constant cylindrical angle in a half-period and 3-D rendering at $A=10$ . The dotted line traces the magnetic axis, with crosses representing the intersection with the cross-sections. (b) Root-mean-square difference in the second-order near-axis magnetic field between the ideal near-axis description and the finite aspect ratio global equilibrium construction with VMEC. A reference quadratic scaling $\propto 1/A^2$ is given, from which the convergence of the error deviated at larger aspect ratios. (c) Explicit comparison between the different poloidal components of the second-order magnetic field magnitude from the global VMEC equilibrium at different aspect ratios, and the ideal near-axis value (black dotted line).

4.3. Half-helicity case

The same systematic construction can be followed for the special QI class of fields with half-helicity axes. We construct an $N=3$ half-helicity field with a smooth $k=5$ buffer. When dealing with half-helicity axes, and especially when high control is required on the curvature and torsion, as is indeed needed in order to control the extreme shaping that is the tendency of QI fields, it is convenient to define the axis by directly parametrising the curvature and the torsion. In that case, one needs to reconstruct the axis by integrating FS equations, (2.3). The well-known drawback of doing so is that curvature and torsion must be chosen carefully to guarantee closure of the axis, which must be achieved numerically to a sufficiently high degree of accuracy. Failing to do so will lead to discontinuities (or lack of smoothness in the best of cases) in the near-axis construction of flux surfaces, and thus will lead to an unphysical field.Footnote ¹⁰ Taking this into account, we consider in this case,

(4.4a) \begin{gather} \kappa [\text {m}^{-1}]=-\frac {8.0725268}{2}\left [1+\cos \left (6\pi \frac {\ell }{L}\right )\right ]\sin \left (3\pi \frac {\ell }{L}\right )\sin \left (6\pi \frac {\ell }{L}\right ), \\[-24pt] \nonumber \end{gather}

(4.4b) \begin{gather} \tau [\text {m}^{-1}]= 1.34174997-0.133333\cos \left (6\pi \frac {\ell }{L}\right ), \\[6pt] \nonumber \end{gather}

and $L [\text {m}]= 2\pi$ , which gives a third-order zero at the bottom of the well and second order at the top. The magnetic field $B_0$ and the shaping $\bar {d}$ are not presented explicitly in the text, but can be found in the repository associated with this paper, as they have been chosen carefully to limit the amount of elongation and shaping of the construction. The steps informing the choice of inputs and the use of the FS approach will be detailed in a future publication. We complete the inputs to the near-axis construction by choosing $X_{2c}=0=X_{2s}$ .

Figure 6. Global equilibrium construction from the second-order near-axis field in § 4.3. The figure presents a second-order near-axis construction for a half-helicity, three period field. (a) Cross-sections at constant cylindrical angle in a half-period and 3-D rendering at $A=10$ . The dotted line traces the magnetic axis, with crosses representing the intersection with the cross-sections. (b) Root-mean-square difference in the second-order near-axis magnetic field between the ideal near-axis description and the finite aspect ratio global equilibrium construction with VMEC. A reference quadratic scaling $\propto 1/A^2$ is given, from which the convergence of the error deviated at larger aspect ratios. The VMEC solver struggles at larger aspect ratios. (c) Explicit comparison between the different poloidal components of the second-order magnetic field magnitude from the global VMEC equilibrium at different aspect ratios, and the ideal near-axis value (black dotted line).

The results of the comparison with VMEC are shown in figure 6, where we get the expected agreement. This shows that the treatment of the half-helicity scenario (as discussed in § 2 and appendix B.1) and is indeed the correct one.

4.4. The QI optimised minimally shaped field

Finally, we present a field which has been optimised to be omnigeneous at second order near the bottom of the well following § 3. In this proof-of-principle exercise, we allowed the $Z$ harmonics as well as $\bar {d}$ to vary as our degrees of freedom, and penalise deviations from omnigeneity at second order in the central 20 % of the domain. The starting point of the construction was the same as that of § 4.2, namely the $N=3$ from Camacho Mata et al. (Reference Camacho Mata, Plunk and Jorge2022). The resulting $Z$ component of the axis,

(4.5) \begin{align} Z [\text {m}] \approx 0.397890 \sin 3\phi + 0.0244030\sin 6\phi + 0.00292413\sin 9\phi, \end{align}

and

(4.6) \begin{align} \bar {d} \approx 0.642233+0.000346819\cos 6\varphi +0.00114030\cos 9\varphi &+0.00366299\cos 12\varphi \nonumber \\ &-0.00230278\cos 15\varphi . \end{align}

Figure 7. Global equilibrium construction from the second-order near-axis field in § 4.4. The figure presents the second-order near-axis construction based on Camacho Mata et al. (Reference Camacho Mata, Plunk and Jorge2022) presented in § 4.4 that has been reshaped to satisfy the QI criterion at second order near $B_{{min}}$ . (a) Cross-sections at constant cylindrical angle in a half-period and 3-D rendering at $A=10$ . The dotted line traces the magnetic axis, with crosses representing the intersection with the cross-sections. (b) Root-mean-square difference in the second-order near-axis magnetic field between the ideal near-axis description and the finite aspect ratio global equilibrium construction with VMEC. A reference quadratic scaling $\propto 1/A^2$ is given, from which the convergence of the error deviated at larger aspect ratios. (c) Explicit comparison between the different poloidal components of the second-order magnetic field magnitude from the global VMEC equilibrium at different aspect ratios, and the ideal near-axis value (black dotted line). The shaded region indicates the $\varphi$ range optimised for QI at second order.

Once again, the comparison with VMEC shows agreement (see figure 7) over a range of aspect ratios, struggling to converge at very large aspect ratios, likely driven by the high degree of toroidal variation.