Skip to main content Accessibility help


  • Access
  • Cited by 2



MathJax is a JavaScript display engine for mathematics. For more information see
      • Send article to Kindle

        To send this article to your Kindle, first ensure is added to your Approved Personal Document E-mail List under your Personal Document Settings on the Manage Your Content and Devices page of your Amazon account. Then enter the ‘name’ part of your Kindle email address below. Find out more about sending to your Kindle. Find out more about sending to your Kindle.

        Note you can select to send to either the or variations. ‘’ emails are free but can only be sent to your device when it is connected to wi-fi. ‘’ emails can be delivered even when you are not connected to wi-fi, but note that service fees apply.

        Find out more about the Kindle Personal Document Service.

        Equilibrium 𝛽-limits in classical stellarators
        Available formats

        Send article to Dropbox

        To send this article to your Dropbox account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Dropbox.

        Equilibrium 𝛽-limits in classical stellarators
        Available formats

        Send article to Google Drive

        To send this article to your Google Drive account, please select one or more formats and confirm that you agree to abide by our usage policies. If this is the first time you use this feature, you will be asked to authorise Cambridge Core to connect with your <service> account. Find out more about sending content to Google Drive.

        Equilibrium 𝛽-limits in classical stellarators
        Available formats
Export citation


A numerical investigation is carried out to understand the equilibrium $\unicode[STIX]{x1D6FD}$ -limit in a classical stellarator. The stepped-pressure equilibrium code (Hudson et al., Phys. Plasmas, vol. 19 (11), 2012) is used in order to assess whether or not magnetic islands and stochastic field-lines can emerge at high $\unicode[STIX]{x1D6FD}$ . Two modes of operation are considered: a zero-net-current stellarator and a fixed-iota stellarator. Despite the fact that relaxation is allowed (Taylor, Rev. Mod. Phys., vol. 58 (3), 1986, pp. 741–763), the former is shown to maintain good flux surfaces up to the equilibrium $\unicode[STIX]{x1D6FD}$ -limit predicted by ideal-magnetohydrodynamics (MHD), above which a separatrix forms. The latter, which has no ideal equilibrium $\unicode[STIX]{x1D6FD}$ -limit, is shown to develop regions of magnetic islands and chaos at sufficiently high $\unicode[STIX]{x1D6FD}$ , thereby providing a ‘non-ideal $\unicode[STIX]{x1D6FD}$ -limit’. Perhaps surprisingly, however, the value of $\unicode[STIX]{x1D6FD}$ at which the Shafranov shift of the axis reaches a fraction of the minor radius follows in all cases the scaling laws predicted by ideal-MHD. We compare our results to the High-Beta-Stellarator theory of Freidberg (Ideal MHD, 2014, Cambridge University Press) and derive a new prediction for the non-ideal equilibrium $\unicode[STIX]{x1D6FD}$ -limit above which chaos emerges.

1 Introduction

In stellarators, the maximum achievable $\unicode[STIX]{x1D6FD}$ is most probably set by the equilibrium and not by its stability (Helander et al. 2012). In fact, magnetic surfaces are not guaranteed to exist in three-dimensional magnetohydrodynamics (MHD) equilibria without a continuous symmetry (Meiss 1992). While stellarators can be designed to possess magnetic surfaces in vacuum (Hanson & Cary 1984; Cary & Hanson 1986; Hudson & Dewar 1997; Pedersen et al. 2016), the necessary existence of plasma currents that maintain force-balance at finite plasma pressure engenders the potential destruction of magnetic surfaces at sufficiently high $\unicode[STIX]{x1D6FD}$ and can thus lead to the loss of confinement (Drevlak, Monticello & Reiman 2005; Suzuki et al. 2008).

The equilibrium $\unicode[STIX]{x1D6FD}$ -limit is not fully understood since it requires the accurate computation of three-dimensional MHD equilibria, which generally consist of an intricate combination of magnetic surfaces, magnetic islands and magnetic field-line chaos. The stepped-pressure equilibrium code (SPEC) was developed as one possible approach to fulfil this highly non-trivial task (Hudson et al. 2012), although there are a few more ongoing parallel efforts (Suzuki et al. 2006; Hirshman, Sanchez & Cook 2011). SPEC has been rigorously verified in axisymmetry (Hudson et al. 2012), in slightly perturbed configurations (Loizu et al. 2015a , b , 2016a ) and more recently in stellarator geometries (Loizu, Hudson & Nührenberg 2016b ).

With a view to progressing towards an understanding of the $\unicode[STIX]{x1D6FD}$ -limit in advanced, fusion-relevant stellarator experiments, we focus on a classical stellarator geometry with a simple pressure pedestal and perform a basic numerical study of its equilibrium $\unicode[STIX]{x1D6FD}$ -limit. The simplified geometry allows us to use the High-Beta-Stellarator model (Freidberg 2014) to guide our investigation. This paper leads to the distinction between ‘ideal’ and ‘relaxed’ equilibrium $\unicode[STIX]{x1D6FD}$ -limits, for which we derive analytical expressions that push our theoretical understanding forward and validate the numerical calculations. Although this study is not aimed at providing experimental predictions yet, the potential consequences of the results presented herein for real experimental situations is briefly discussed at the end, together with some of the open questions.

2 Model and control parameters

We consider the fixed-boundary problem of a finite $\unicode[STIX]{x1D6FD}$ equilibrium in a classical $l=2$ stellarator (Freidberg 2014). Namely, we must provide (i) the geometry of the boundary, e.g. via the Fourier coefficients of the cylindrical coordinates defining the boundary surface, $\{R_{mn},Z_{mn}\}$ ; (ii) the pressure profile as a function of the enclosed toroidal magnetic flux, $p(\unicode[STIX]{x1D6F9})$ ; and (iii) an additional profile, e.g. the rotational transform, , or the net toroidal current, $I_{\unicode[STIX]{x1D711}}(\unicode[STIX]{x1D6F9})$ . The total toroidal magnetic flux enclosed by the boundary, $\unicode[STIX]{x1D6F9}_{\text{edge}}$ , is also provided but its value is irrelevant and only acts as a global scale-factor for the magnetic field strength.

2.1 Boundary

The simplest boundary representation that can model an $l=2$ stellarator is that of a rotating ellipse with no toroidally averaged elongation. Namely,

(2.1) $$\begin{eqnarray}\displaystyle \left.\begin{array}{@{}c@{}}R(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})=R_{00}+R_{10}\cos \unicode[STIX]{x1D703}+R_{11}\cos (\unicode[STIX]{x1D703}-N_{p}\unicode[STIX]{x1D711})\\[5.0pt] Z(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})=Z_{00}+Z_{10}\sin \unicode[STIX]{x1D703}+Z_{11}\sin (\unicode[STIX]{x1D703}-N_{p}\unicode[STIX]{x1D711}),\end{array}\right\} & & \displaystyle\end{eqnarray}$$

with $Z_{00}=0$ , $Z_{10}=-R_{10}$ and $Z_{11}=R_{11}$ . For our $\unicode[STIX]{x1D6FD}$ -limit study, the main parameters of interest in (2.1) are the major radius, $R_{00}$ , and the number of field periods, $N_{p}$ . In fact these can be used to vary independently the inverse aspect ratio, $\unicode[STIX]{x1D716}$ , and the vacuum rotational transform, , which are predicted to determine the ideal equilibrium $\unicode[STIX]{x1D6FD}$ -limit. We therefore choose to fix the other parameters to $R_{10}=1$ and $R_{11}=0.25$ . Two examples of such boundaries with different values of $N_{p}$ are shown in figure 1.

The inverse aspect ratio is

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D716}=\frac{r_{\text{eff}}}{R_{00}},\end{eqnarray}$$

where the effective minor radius is $r_{\text{eff}}=\sqrt{r_{\text{max}}r_{\text{min}}}$ , with $r_{\text{max}}=R_{10}+R_{11}=1.25$ and $r_{\text{min}}=R_{10}-R_{11}=0.75$ , respectively the major and minor axis of the rotating ellipse. The vacuum rotational transform can be estimated analytically (Helander 2014) as


For example, for $N_{p}=5$ we get .

Figure 1. Boundary of a classical $l=2$ stellarator with $N_{p}=5$ (a) and $N_{p}=10$ (b) field periods. The inverse aspect ratio is $\unicode[STIX]{x1D716}=0.1$ and the colour represents the amplitude of the vacuum magnetic field on the boundary as computed from SPEC.

2.2 Pressure profile

We model a pressure pedestal by assuming that all the pressure gradient is concentrated on a single flux-surface, namely $p(\unicode[STIX]{x1D6F9})=p_{0}$ for $\unicode[STIX]{x1D6F9}\leqslant \unicode[STIX]{x1D6F9}_{a}$ and $p(\unicode[STIX]{x1D6F9})=0$ for $\unicode[STIX]{x1D6F9}\geqslant \unicode[STIX]{x1D6F9}_{a}$ . This step in the pressure is naturally described by the SPEC code: two Taylor-relaxed volumes (Taylor 1974) separated by an ideal-interface supporting a pressure step $\unicode[STIX]{x27E6}p\unicode[STIX]{x27E7}=p(\unicode[STIX]{x1D6F9}_{a}^{+})-p(\unicode[STIX]{x1D6F9}_{a}^{-})=p_{0}$ , in correspondence to which a jump in $B$ must arise according to $\unicode[STIX]{x27E6}p+(B^{2}/(2\unicode[STIX]{x1D707}_{0}))\unicode[STIX]{x27E7}=0$ . This implies, by virtue of Ampère’s law, the presence of a surface current density, $\boldsymbol{j}=\boldsymbol{n}\times \unicode[STIX]{x27E6}\boldsymbol{B}\unicode[STIX]{x27E7}~\unicode[STIX]{x1D6FF}(\boldsymbol{x}-\boldsymbol{x}_{a})$ , where $\unicode[STIX]{x1D6FF}(x)$ is the Dirac delta-function, $\boldsymbol{x}=\boldsymbol{x}_{a}$ defines the points on the surface across which the jump occurs and $\boldsymbol{n}$ is the unit vector normal to that surface. This current density is simply a weak representation of the pressure-driven current density (diamagnetic and Pfirsch–Schlüter), by which we mean that it shall be interpreted only in the integral sense (see § 2.3).

For our basic $\unicode[STIX]{x1D6FD}$ -limit study, we choose to fix the value $\unicode[STIX]{x1D6F9}_{a}=0.3\unicode[STIX]{x1D6F9}_{\text{edge}}$ and use the freedom in $p_{0}$ to control the value of $\unicode[STIX]{x1D6FD}$ , which we define here as $\unicode[STIX]{x1D6FD}=2\unicode[STIX]{x1D707}_{0}p_{0}/B_{0}^{2}$ , where $B_{0}=B(\unicode[STIX]{x1D6F9}=0)$ .

We would like to remark that the SPEC code can handle pressure profiles with arbitrarily many interfaces, hence arbitrarily many pressure jumps, thereby allowing the equilibrium to approach the ideal-MHD limit (Dennis et al. 2013a ; Loizu et al. 2015b ). In this study, however, we are looking for the ‘worst case scenario’ and hence we consider a single, localized pressure pedestal. By ‘worst case scenario’ we mean that ‘maximally relaxed’ or ‘minimally constrained’ equilibrium states are sought. In fact, when the magnetic helicity is assumed to be conserved only globally in a certain volume (Taylor 1986), magnetic reconnection is ‘maximal’ within that volume: there are no more possible reconnecting events that could lower the plasma potential energy. If the pressure gradient was distributed among many surfaces, the available volume for relaxation would be smaller; and in the ideal-MHD limit with all surfaces supporting a pressure gradient, no reconnection would be possible. To conclude this discussion: whenever good flux surfaces are to be found despite the allowed relaxation, it shall be understood that no possible relaxation mechanism can destroy these surfaces. When, on the other hand, islands and chaotic field-lines are produced, it shall be inferred that this is the ‘worst case scenario’ and the potential destruction of flux-surfaces is subject to the available relaxation and healing mechanisms. A similar approach was taken (Dennis et al. 2013b ) in order to reproduce self-organized helical states in reversed field pinches.

2.3 Zero-net-current versus fixed-iota

The SPEC code calculates MHD equilibria as extrema of the Multiregion, Relaxed MHD (MRxMHD) energy functional (Hole, Hudson & Dewar 2007; Hudson, Hole & Dewar 2007). In essence, the energy functional is the same as in conventional ideal MHD equilibrium theory (Kruskal & Kulsrud 1958), but the constraints under which the function is extremized are different. While in ideal-MHD the magnetic topology is continuously constrained, in MRxMHD the topology is only discretely constrained, thus allowing for partial relaxation. More precisely, the plasma is partitioned into a finite number, $N_{V}$ , of nested volumes, $V_{v}$ , that undergo Taylor relaxation. These volumes are separated by $N_{V}-1$ interfaces that are constrained to remain magnetic surfaces during the energy minimization process. For the $\unicode[STIX]{x1D6FD}$ -limit study at hand, we have $N_{V}=2$ volumes separated by one ideal-interface. The location and shape of this interface is unkown a priori and determined self-consistently by a force-balance condition. MRxMHD equilibrium states satisfy

(2.4) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\times \boldsymbol{B}=\unicode[STIX]{x1D707}_{v}\boldsymbol{B}\quad \text{in the volumes} & \displaystyle\end{eqnarray}$$
(2.5) $$\begin{eqnarray}\displaystyle & \displaystyle \left[\hspace{-2.0pt}\left[p+\frac{B^{2}}{2\unicode[STIX]{x1D707}_{0}}\right]\hspace{-2.0pt}\right]=0\quad \text{on the interface} & \displaystyle\end{eqnarray}$$

for $v=1,2$ . In addition to providing the enclosed toroidal fluxes in each volume ( $\unicode[STIX]{x1D6F9}_{a}$ and $\unicode[STIX]{x1D6F9}_{\text{edge}}$ ), the solution to equation (2.4) requires one more parameter if the volume is a topological torus (the innermost volume) and two more parameters if the volume is an annulus (the outer volume). Hence we must provide a total of 3 parameters to determine the equilibrium solution at a given value of $\unicode[STIX]{x1D6FD}$ .

If we want to enforce a zero net-toroidal-current, $I_{\unicode[STIX]{x1D711}}=0$ , we can impose $\unicode[STIX]{x1D707}_{1}=\unicode[STIX]{x1D707}_{2}=0$ and then iterate on the total enclosed poloidal flux, $\unicode[STIX]{x1D713}_{p}$ , until the surface current has no net toroidal component. At each iteration step, the net toroidal surface current can be easily calculated as

(2.6) $$\begin{eqnarray}I_{\unicode[STIX]{x1D711}}^{\text{CS}}=\int _{0}^{2\unicode[STIX]{x03C0}}\unicode[STIX]{x27E6}\boldsymbol{B}\unicode[STIX]{x27E7}\boldsymbol{\cdot }\boldsymbol{e}_{\unicode[STIX]{x1D703}}\,\text{d}\unicode[STIX]{x1D703}\end{eqnarray}$$

by virtue of Ampère’s law. Here $\boldsymbol{e}_{\unicode[STIX]{x1D703}}=\unicode[STIX]{x2202}_{\unicode[STIX]{x1D703}}\boldsymbol{x}$ and $\boldsymbol{x}(\unicode[STIX]{x1D703},\unicode[STIX]{x1D711})$ parametrizes the surface. The iterative procedure can be implemented via a Newton method and brings $I_{\unicode[STIX]{x1D711}}^{\text{CS}}$ down to machine precision in a few steps. We refer to this mode of operation as zero-net-current.

If we want to constrain the rotational transform, , we can force it to remain constant on both sides of the ideal-interface,  and at the edge, . Once again, this can be achieved by iterating on the values of $\unicode[STIX]{x1D707}_{1,2}$ and $\unicode[STIX]{x1D713}_{p}$ . We refer to this mode of operation as fixed-iota.

We would like to remark that while the zero-net-current mode guarantees $I_{\unicode[STIX]{x1D711}}=0$ , it does not guarantee that the rotational transform remains constant and in particular we expect . Conversely, the fixed-iota mode guarantees that  remains constant on certain surfaces but in general we expect $I_{\unicode[STIX]{x1D711}}\neq 0$ , in particular at the location of the pressure-gradient.

3 High- $\unicode[STIX]{x1D6FD}$ equilibria and Shafranov shift

Figure 2 shows Poincaré plots of the equilibrium magnetic field at different values of $\unicode[STIX]{x1D6FD}$ for both the zero-net-current stellarator and the fixed-iota stellarator. In both cases there is a Shafranov shift that increases with $\unicode[STIX]{x1D6FD}$ . However, the Shafranov shift of the axis, $\unicode[STIX]{x1D6E5}_{\text{ax}}$ , increases with $\unicode[STIX]{x1D6FD}$ much faster in the zero-net-current stellarator. This can be explained as follows. For small $\unicode[STIX]{x1D6FD}$ , one expects  (Miyamoto 2005). As shown in § 4 decreases with $\unicode[STIX]{x1D6FD}$ in the zero-net-current stellarator, thus amplifying the effect.

Figure 2. Poincaré section ( $\unicode[STIX]{x1D711}=0^{\circ }$ ) of the equilibrium magnetic field at different values of $\unicode[STIX]{x1D6FD}$ . (a) Zero-net-current stellarator. (b) Fixed-iota stellarator. Here $N_{p}=5$ and $\unicode[STIX]{x1D716}=0.1$ . Indicated in red are the boundary surface and inner interface supporting the pressure pedestal.

It is useful to define the quantity

(3.1) $$\begin{eqnarray}\unicode[STIX]{x1D6FD}_{0.5}\equiv \unicode[STIX]{x1D6FD}\left(\unicode[STIX]{x1D6E5}_{\text{ax}}=\frac{r_{\text{eff}}}{2}\right),\end{eqnarray}$$

namely, the value of $\unicode[STIX]{x1D6FD}$ at which the Shafranov shift of the axis reaches half of the minor radius. According to ideal-MHD equilibrium theory (Miyamoto 2005), $\unicode[STIX]{x1D6FD}_{0.5}$ is predicted to scale as


for large aspect ratios, $\unicode[STIX]{x1D716}\ll 1$ , and slowly varying , which is true for $N_{p}\sim 1$ . A scan in both $R_{00}$ and $N_{p}$ has been carried out in order to assess how $\unicode[STIX]{x1D6FD}_{0.5}$ scales in the numerical MHD calculations. Figure 3 shows the result of this scan. Despite the fact that SPEC allows for plasma relaxation, the scaling law (3.2) is very well reproduced in both modes of operation; except at high values of $\unicode[STIX]{x1D716}\sim 1/R_{00}$ and $N_{p}$ , for which the scaling law ceases to be valid. Moreover, the values of $\unicode[STIX]{x1D6FD}_{0.5}$ are much higher in the fixed-iota stellarator, by a factor of about 6. This reflects, once again, the fact that the Shafranov shift increases faster in the zero-net-current stellarator.

In § 4, the fundamental, macroscopic differences between the two modes of operation are explained in terms of the High-Beta-Stellarator (HBS) model developed in Freidberg (2014). In § 5, we attempt to describe and predict the $\unicode[STIX]{x1D6FD}$ -induced generation of islands and chaotic field-lines in the fixed-iota stellarator (figure 2).

Figure 3. Scaling of $\unicode[STIX]{x1D6FD}_{0.5}$ with the inverse aspect ratio, $\unicode[STIX]{x1D716}\sim 1/R_{00}$ (a), and with the vacuum iota, $\,\unicode[STIX]{x1D704}\text{-}_{v}\sim N_{p}$ (b). Black circles are for the fixed-iota stellarator. Magenta stars are for the zero-net-current stellarator. The dashed lines have slope 1 (a) and 2 (b).

Figure 4. Rotational transform as a function of toroidal magnetic flux from both SPEC (black stars) and VMEC (solid blue line) at $\unicode[STIX]{x1D6FD}=0.15\,\%$ . For comparison, the vacuum transform is also shown (dashed magenta line). Here $N_{p}=5$ , $\unicode[STIX]{x1D716}=0.1$ and the vertical dashed line indicates the location of the pressure pedestal.

4 Ideal $\unicode[STIX]{x1D6FD}$ -limit and the HBS theory

The HBS model for a classical stellarator developed in Freidberg (2014) and briefly summarized in appendix A predicts that the rotational transform at the plasma edge, , evolves with $\unicode[STIX]{x1D6FD}$ and plasma current as


where  is the transform produced by the net toroidal current,




where $a$ is the effective minor radius of the plasma edge and $\unicode[STIX]{x1D716}_{a}=a/R$ . For our system, we have $a=\sqrt{\unicode[STIX]{x1D6F9}_{a}/\unicode[STIX]{x1D6F9}_{\text{edge}}}r_{\text{eff}}$ and thus $\unicode[STIX]{x1D716}_{a}=\sqrt{\unicode[STIX]{x1D6F9}_{a}/\unicode[STIX]{x1D6F9}_{\text{edge}}}\unicode[STIX]{x1D716}$ .

In the context of the HBS theory, the zero-net-current stellarator can be analyzed by taking . Equation (4.1) then implies that  decreases with increasing $\unicode[STIX]{x1D6FD}$ . This is visible in figure 4, where the profile  obtained from SPEC at finite $\unicode[STIX]{x1D6FD}$ is shown and compared to the vacuum transform. A jump in the rotational transform self-consistently develops on the ideal interface supporting the pressure gradient, namely at $\unicode[STIX]{x1D6F9}_{a}=0.3\unicode[STIX]{x1D6F9}_{\text{edge}}$ . The ideal MHD, variational moments equilibrium code, or VMEC (Hirshman & Whitson 1983), was also run for this case with a pressure pedestal of small but finite width (the calculation requires a rather high radial resolution, with about 3000 flux surfaces) and shown to produce essentially the same transform profile.

Figure 5. Rotational transform at the plasma edge, $\,\unicode[STIX]{x1D704}\text{-}_{a}^{+}=\,\unicode[STIX]{x1D704}\text{-}(\unicode[STIX]{x1D6F9}_{a}^{+})$ , as a function of $\unicode[STIX]{x1D6FD}$ , from SPEC calculations (stars), from the HBS prediction with Solov’ev pressure profile, equation (4.1) (dashed line) and from the HBS prediction adapted to a stepped-pressure profile, equation (A 10) (dash-dotted line). Two cases are shown: $N_{p}=3$ and $N_{p}=5$ at fixed $\unicode[STIX]{x1D716}=0.1$ .

In figure 5, the value of  is shown as a function of $\unicode[STIX]{x1D6FD}$ and compared to the HBS prediction, equation (4.1), showing fairly good agreement (notice that there are no free parameters). The agreement is even more remarkable if we notice that the HBS model assumes a circular cross-section and uses Solov’ev pressure profiles. In appendix A, we show that adapting the theory to the case of a stepped-pressure profile produces very similar predictions (see figure 5). This reflects the fact that the macroscopic equilibrium depends on integral quantities (e.g., the total plasma pressure) and not so much on the details of the profiles.

The point where , which from (4.1) happens when $\unicode[STIX]{x1D708}=1$ , corresponds to the emergence of a separatrix (see, e.g., figure 2) and this defines the ideal $\unicode[STIX]{x1D6FD}$ -limit, namely


For example, for the $N_{p}=5$ case depicted in figure 5 we have $\unicode[STIX]{x1D716}_{a}\approx 0.05$ and , thus (4.4) gives $\unicode[STIX]{x1D6FD}_{\text{lim}}\approx 0.4\,\%$ , in good agreement with the SPEC calculations. We note that the separatrix forming at $\unicode[STIX]{x1D6FD}>\unicode[STIX]{x1D6FD}_{\text{lim}}$ seems to possess two X-points that connect to the current sheet established on the ideal-surface (see middle-left panel in figure 2); a situation that is reminiscent of the Waelbroeck ribbon forming in the reconnecting kink-tearing mode in tokamaks (Waelbroeck 1989; Zakharov, Rogers & Migliuolo 1993).

For the fixed-iota stellarator, we can impose  in the HBS model. This leads to an expression for the value of the plasma current that is necessary to keep  constant. One obtains (Freidberg 2014)



(4.6) $$\begin{eqnarray}H=\frac{\unicode[STIX]{x1D6FD}}{\unicode[STIX]{x1D6FD}_{\text{lim}}}.\end{eqnarray}$$

Figure 6 shows the net toroidal surface-current, $I_{\unicode[STIX]{x1D711}}$ , self-consistently generated in SPEC equilibria as a function of $\unicode[STIX]{x1D6FD}$ and compares it to the HBS prediction, equation (4.5), showing good agreement (again, there are no free parameters). Since the predicted current is entirely pressure-driven, we do not expect currents to develop in the relaxed volumes and we have checked that for these SPEC calculations the values of $\unicode[STIX]{x1D707}_{v}$ remain small $({<}10^{-2})$ , such that the corresponding volume currents are more than 10 times smaller than the pressure-driven surface-current. For large $H\gg 1$ , one has $I_{\unicode[STIX]{x1D711}}\sim \sqrt{\unicode[STIX]{x1D6FD}}$ and the HBS model predicts that no $\unicode[STIX]{x1D6FD}$ -limit is reached because the plasma current keeps rising and preventing the separatrix to form. From SPEC equilibrium calculations, however, where plasma relaxation is allowed, we observe that magnetic islands and chaotic field-lines emerge at sufficiently high $\unicode[STIX]{x1D6FD}$ , thereby providing a ‘non-ideal $\unicode[STIX]{x1D6FD}$ -limit’.

Figure 6. Net toroidal plasma current as a function of $\unicode[STIX]{x1D6FD}$ , from SPEC calculations (stars) and from (4.5) for both Solov’ev profiles (dashed line) and stepped-pressure profile (dash-dotted line). Two cases are shown: $\unicode[STIX]{x1D716}=1/10$ and $\unicode[STIX]{x1D716}=1/20$ at fixed $N_{p}=5$ .

5 Non-ideal $\unicode[STIX]{x1D6FD}$ -limit and emergence of chaos

We can quantify the emergence of chaos by calculating the fractal dimension of the field-lines on the Poincaré section as a function of $\unicode[STIX]{x1D6FD}$ (Meiss 1992). More precisely, we can evaluate the so-called box-counting dimension, or Hausdorff dimension,

(5.1) $$\begin{eqnarray}D=\lim _{L\rightarrow 0}\left|\frac{\log (N)}{\log (L)}\right|,\end{eqnarray}$$

where $L$ is the size of the boxes and $N$ is the number of boxes containing at least one point of the magnetic field-line on the Poincaré section. If the field-line traces a magnetic surface, or even a magnetic island, one expects $D=1$ . If the magnetic field-line trajectory is chaotic, however, it fills up a certain ‘area’ in the Poincaré section and $D>1$ is expected. We remark that the accurate evaluation of $D$ requires a large number of toroidal transits, $N_{\text{trans}}$ , when generating the Poincaré section via field-line-tracing. Satisfactory convergence was found at values of about $N_{\text{trans}}>2\times 10^{4}$ . The value of $D$ for each field-line $i=1,\ldots ,n_{\text{lines}}$ can be plotted against an effective toroidal flux, $\unicode[STIX]{x1D6F9}_{i}\sim \unicode[STIX]{x1D70C}_{i}^{2}$ , which is proportional to the square of the radial coordinate, $\unicode[STIX]{x1D70C}_{i}$ , obtained from the interpolation of the inner and outer ideal surfaces.

Figure 7 shows the calculated fractal dimension as a function of the toroidal flux in equilibria of increasing $\unicode[STIX]{x1D6FD}$ . First, we observe that for sufficiently low $\unicode[STIX]{x1D6FD}$ we obtain $D(\unicode[STIX]{x1D6F9})=1$ , as expected, because magnetic surfaces are preserved in the entire volume. Second, we notice that for sufficiently high $\unicode[STIX]{x1D6FD}$ there are regions in which $D(\unicode[STIX]{x1D6F9})>1$ for $\unicode[STIX]{x1D6F9}>\unicode[STIX]{x1D6F9}_{a}$ . Third, the value of $D$ seems to be almost-binary, taking values at either $D\approx 1$ or $D\approx 1.6$ . Fourth, the regions with $D\approx 1.6$ correspond to what appears to be stochastic regions in the corresponding Poincaré section. Finally, the volume occupied by these regions increases with $\unicode[STIX]{x1D6FD}$ . These observations suggest that $D$ is a good proxy for the emergence of chaos, which greatly simplifies the task of probing the ‘non-ideal $\unicode[STIX]{x1D6FD}$ -limit’. In fact, we can now define the volume of chaos, $V_{\text{chaos}}$ , in the system as

(5.2) $$\begin{eqnarray}V_{\text{chaos}}=V_{\text{tot}}\mathop{\sum }_{i=1}^{n_{\text{lines}}}\frac{(\unicode[STIX]{x1D6F9}_{i}-\unicode[STIX]{x1D6F9}_{i-1})}{\unicode[STIX]{x1D6F9}_{\text{edge}}}{\mathcal{H}}(D(\unicode[STIX]{x1D6F9}_{i})-D_{\text{crit}}),\end{eqnarray}$$

where $V_{\text{tot}}$ is the total volume defined by the fixed-boundary, $n_{\text{lines}}$ is the number of traced field-lines, $\unicode[STIX]{x1D6F9}_{i}-\unicode[STIX]{x1D6F9}_{i-1}$ measures the enclosed toroidal flux between two neighbouring field lines, and ${\mathcal{H}}$ is the Heaviside function, with ${\mathcal{H}}=0$ for $D<D_{\text{crit}}=1.5$ and ${\mathcal{H}}=1$ otherwise. Figure 8 shows the profile of $V_{\text{chaos}}(\unicode[STIX]{x1D6FD})$ calculated for three different values of $N_{p}$ . Clearly, the emergence of chaos occurs at some critical value of $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{\text{chaos}}$ , which we define as the non-ideal equilibrium $\unicode[STIX]{x1D6FD}$ -limit. The question remains: can we theoretically predict the value of $\unicode[STIX]{x1D6FD}_{\text{chaos}}$ ?

Figure 7. Fractal dimension of the magnetic field lines in a Poincaré section as a function of toroidal flux. The pressure pedestal is at $\unicode[STIX]{x1D6F9}/\unicode[STIX]{x1D6F9}_{\text{edge}}=0.3$ . Different curves are for different values of $\unicode[STIX]{x1D6FD}$ . All equilibria have $N_{p}=5$ and $\unicode[STIX]{x1D716}=0.1$ .

Figure 8. Volume of chaos as a function of $\unicode[STIX]{x1D6FD}$ for $N_{p}=3$ , $N_{p}=5$ and $N_{p}=10$ , at fixed $\unicode[STIX]{x1D716}=0.1$ . Vertical dashed lines indicate the predicted transition to chaos at $\unicode[STIX]{x1D6FD}=\unicode[STIX]{x1D6FD}_{\text{chaos}}$ given by (5.4).

We expect that chaos emerges due to the overlapping of magnetic islands according to the Chirikov criterion (Chirikov 1979). That requires, however, predicting which resonances appear first and how does the width of the corresponding islands grow with $\unicode[STIX]{x1D6FD}$ . While we intend to investigate this question in more detail in the future, we derive here a criterion based on the following general idea. The expected width of an island generated by a resonance at  is expected to scale as  (Boozer 2004). Thus, new resonances can appear at finite $\unicode[STIX]{x1D6FD}$ due to (i) a change in the rotational transform, or (ii) a change in the harmonic content of $B$ inevitably caused by the Shafranov shift. As $\unicode[STIX]{x1D6FD}$ increases, the net toroidal current increases and modifies the rotational transform by increasing its value in the region $\unicode[STIX]{x1D6F9}_{a}<\unicode[STIX]{x1D6F9}<\unicode[STIX]{x1D6F9}_{\text{edge}}$ (data not shown). The rising of  allows new resonances to appear (and with lower values of $m$ ). At this point we make the following hypothesis: the emergence of chaos may occur when the perturbations in the poloidal field due to finite toroidal current are comparable to the vacuum poloidal field. Namely, chaos may emerge when . From the HBS theory we know that  increases with $\unicode[STIX]{x1D6FD}$ according to (4.5), hence applying our constraint we have


and hence


which gives $\unicode[STIX]{x1D6FD}_{\text{chaos}}\approx 0.5\,\%$ for $N_{p}=3$ , $\unicode[STIX]{x1D6FD}_{\text{chaos}}\approx 1.4\,\%$ for $N_{p}=5$ , and $\unicode[STIX]{x1D6FD}_{\text{chaos}}\approx 6.5\,\%$ for $N_{p}=10$ , thus in excellent agreement with the observed transition to chaos (see figure 8). More importantly, equation (5.4) predicts that the non-ideal equilibrium $\unicode[STIX]{x1D6FD}$ -limit scales exactly as the ideal equilibrium $\unicode[STIX]{x1D6FD}$ -limit but with a larger factor in front, of the order of $\sqrt{12}\approx 3.5$ .

6 Discussion

The equilibrium $\unicode[STIX]{x1D6FD}$ -limit in a classical stellarator has been thoroughly investigated via numerical calculations that have guided our analytical understanding. A classical stellarator with zero net-toroidal-current possesses an equilibrium $\unicode[STIX]{x1D6FD}$ -limit as predicted by ideal MHD, , above which a separatrix forms due to the vanishing of the rotational transform at the plasma edge, . A classical stellarator with constant , however, has a higher equilibrium $\unicode[STIX]{x1D6FD}$ -limit, which is of non-ideal nature. In fact,  can be maintained at any value of $\unicode[STIX]{x1D6FD}$ as long as a net-toroidal-current flows in the vicinity of the pressure pedestal; when such current produces a change in transform that is comparable to the vacuum transform, , magnetic field-line chaos emerges in maximally relaxed equilibria and this occurs at . For $\unicode[STIX]{x1D6FD}>\unicode[STIX]{x1D6FD}_{\text{chaos}}$ , the volume of destroyed magnetic surfaces increases monotonically with $\unicode[STIX]{x1D6FD}$ and radially outward from the location of the pressure pedestal. We remark that this non-ideal $\unicode[STIX]{x1D6FD}$ -limit does not consider the possibility of island-healing mechanisms; on the contrary, it considers the ‘worst case scenario’ of complete relaxation. Therefore $\unicode[STIX]{x1D6FD}_{\text{chaos}}$ should be interpreted as a lower bound for the $\unicode[STIX]{x1D6FD}$ -limit of a classical stellarator where a net-toroidal-current clamps the value of . Furthermore, we would like to emphasize that a relatively small toroidal current is enough to maintain  constant and therefore to raise the $\unicode[STIX]{x1D6FD}$ -limit. For example, for a classical stellarator with $N_{p}=5$ , $\unicode[STIX]{x1D716}=0.1$ , $B_{0}\sim 1$ T and $R_{0}=10$ m, we have $I_{\unicode[STIX]{x1D711}}\approx 40$ kA at about $\unicode[STIX]{x1D6FD}\approx 2\,\%$ (figure 6). This current is easily overwhelmed by the bootstrap current.

Two questions remain to be investigated: (1) can this predictive theory be extended to more complex stellarator geometries and non-trivial pressure profiles? (2) How to incorporate the possibility of pressure-induced island healing (Bhattacharjee et al. 1995; Narushima et al. 2008; Hegna 2012) in the derivation of the equilibrium $\unicode[STIX]{x1D6FD}$ -limit? Some new ideas are needed here.


The authors would like to thank S. Lazerson, A. Alonso and J. Nührenberg for useful discussions. This work has been carried out in the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014–2018 under grant agreement no. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

Appendix A

The main model assumptions and equations in the HBS theory (Freidberg 2014) when applied to the high- $\unicode[STIX]{x1D6FD}$ equilibrium of an $l=2$ stellarator are as follows. The starting equilibrium equations are the usual ideal-MHD equations (  $\boldsymbol{j}\times \boldsymbol{B}=\unicode[STIX]{x1D735}p$ , $\unicode[STIX]{x1D735}\times \boldsymbol{B}=\unicode[STIX]{x1D707}_{0}\,\boldsymbol{j}$ , $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{B}=0$ ) and a reduction is carried out by expansion in the small parameter $\unicode[STIX]{x1D716}=a/R_{0}$ , namely the inverse aspect ratio. An ordering is assumed for the key normalized plasma parameters: the helical field amplitude, $\unicode[STIX]{x1D6FF}=|\boldsymbol{B}_{p}|/B_{\unicode[STIX]{x1D711}}\sim \unicode[STIX]{x1D716}$ , the plasma pressure, $\unicode[STIX]{x1D6FD}=2\unicode[STIX]{x1D707}_{0}p/B_{\unicode[STIX]{x1D711}}^{2}\sim \unicode[STIX]{x1D716}$ , the toroidal periodicity mode number $N\sim 1$ and the poloidal periodicity mode number, $l\sim 1$ . The resulting reduced equations are then compared to the Greene–Johnson stellarator model (Greene & Johnson 1961), which assumes a different ordering, namely $\unicode[STIX]{x1D6FF}\sim \unicode[STIX]{x1D716}^{1/2}$ , $\unicode[STIX]{x1D6FD}\sim \unicode[STIX]{x1D716}$ , $N\sim 1/\unicode[STIX]{x1D716}$ and $l\sim 1$ . A useful region of overlap is identified, with an intermediate ordering $\unicode[STIX]{x1D6FF}\sim \unicode[STIX]{x1D716}^{3/4}$ , $\unicode[STIX]{x1D6FD}\sim \unicode[STIX]{x1D716}$ , $N\sim 1/\unicode[STIX]{x1D716}^{1/2}$ and $l\sim 1$ . A substantial amount of analysis is required to reduce the HBS model to the Greene–Johnson model in the overlap region (Freidberg 2014). The final result is a Grad–Shafranov-like partial differential equation, also called the overlap equation,

(A 1) $$\begin{eqnarray}\unicode[STIX]{x1D735}_{\bot }^{2}\unicode[STIX]{x1D713}=-J-\frac{\text{d}\unicode[STIX]{x1D6FD}}{\text{d}\unicode[STIX]{x1D713}}(x-\langle x\rangle )+\unicode[STIX]{x1D735}_{\bot }^{2}\left(\mathop{\sum }_{1}^{\infty }\frac{i}{n}\boldsymbol{e}_{\unicode[STIX]{x1D711}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\!\bot }A_{n}^{\ast }\times \unicode[STIX]{x1D735}_{\!\bot }A_{n}\right)\!,\end{eqnarray}$$

where $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})$ is the normalized poloidal flux, $J(\unicode[STIX]{x1D713})$ is the normalized toroidal current density, $\unicode[STIX]{x1D6FD}(\unicode[STIX]{x1D713})$ is the normalized plasma pressure and $A_{n}(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})$ are the vector potentials associated with the vacuum helical fields. The operator $\unicode[STIX]{x1D735}_{\!\bot }$ is the gradient in the polar $(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})$ coordinates, $x=\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}$ , and $\langle ~\rangle$ denotes an average over the poloidal angle as $\langle Q\rangle =(\oint Q(r,\unicode[STIX]{x1D703})\,\text{d}l_{p}/|\unicode[STIX]{x1D735}_{\!\bot }\unicode[STIX]{x1D713}|)/(\oint \text{d}l_{p}/|\unicode[STIX]{x1D735}_{\!\bot }\unicode[STIX]{x1D713}|)$ , with $l_{p}(\unicode[STIX]{x1D703})$ the poloidal arc length.

Equation (A 1) can be solved by choosing a single $l=2$ vacuum field helicity,

(A 2) $$\begin{eqnarray}A_{n}(\unicode[STIX]{x1D70C},\unicode[STIX]{x1D703})=\frac{\unicode[STIX]{x1D6E5}_{2}}{4}\unicode[STIX]{x1D70C}^{2}\text{e}^{\text{i}2\unicode[STIX]{x1D703}},\end{eqnarray}$$

and Solov’ev profiles for the free-functions,

(A 3) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D6FD}}{\text{d}\unicode[STIX]{x1D713}}=C & \displaystyle\end{eqnarray}$$
(A 4) $$\begin{eqnarray}\displaystyle & \displaystyle J=A+C\langle x\rangle , & \displaystyle\end{eqnarray}$$

where $\unicode[STIX]{x1D6E5}_{2}$ , $C$ and $A$ are constants. This reduces the overlap equation to

(A 5) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D713}=-A-2\unicode[STIX]{x1D6E5}_{2}^{2}-C\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703}.\end{eqnarray}$$

Equation (A 5) with boundary condition $\unicode[STIX]{x1D713}(\unicode[STIX]{x1D70C}_{a},\unicode[STIX]{x1D703})=0$ (thus assuming a fixed, circular boundary) is easily solved analytically. The solution can be written as

(A 6)

where  and $\unicode[STIX]{x1D708}$ are defined in § 4. Finally, one can evaluate the rotational transform at the plasma edge, , with

(A 7)

which is exactly (4.1).

If a stepped-pressure profile is considered instead of a Solov’ev profile, for example by replacing (A 3) with

(A 8) $$\begin{eqnarray}\frac{\text{d}\unicode[STIX]{x1D6FD}}{\text{d}\unicode[STIX]{x1D713}}=C\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D713}-\unicode[STIX]{x1D713}_{a}),\end{eqnarray}$$

and assuming that the plasma-vacuum surface, $\unicode[STIX]{x1D713}=\unicode[STIX]{x1D713}_{a}$ , remains approximately circular, then (A 1) reduces to

(A 9) $$\begin{eqnarray}\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D713}=-A-2\unicode[STIX]{x1D6E5}_{2}^{2}-\frac{C}{\unicode[STIX]{x1D713}_{0}^{\prime }(\unicode[STIX]{x1D70C}_{a})}\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{a})\unicode[STIX]{x1D70C}\cos \unicode[STIX]{x1D703},\end{eqnarray}$$

where we have written $\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D713}-\unicode[STIX]{x1D713}_{a})=\unicode[STIX]{x1D6FF}(\unicode[STIX]{x1D70C}-\unicode[STIX]{x1D70C}_{a})/\unicode[STIX]{x1D713}_{0}^{\prime }(\unicode[STIX]{x1D70C}_{a})$ with $\unicode[STIX]{x1D713}_{0}(\unicode[STIX]{x1D70C})$ the solution to equation (A 1) in the $\unicode[STIX]{x1D6FD}=0$ limit. Equation (A 9) can be easily solved for $\unicode[STIX]{x1D70C}<\unicode[STIX]{x1D70C}_{a}$ and for $\unicode[STIX]{x1D70C}>\unicode[STIX]{x1D70C}_{a}$ , and an appropriate matching determines the global solution. The resulting final expression for the rotational transform at the plasma edge, , is almost the same as (A 7),

(A 10)


(A 11) $$\begin{eqnarray}\unicode[STIX]{x1D706}=\left(\frac{1+\unicode[STIX]{x1D70C}_{a}^{2}}{2\unicode[STIX]{x1D70C}_{a}}\right)\unicode[STIX]{x1D708}.\end{eqnarray}$$

In particular, the prediction for the $\unicode[STIX]{x1D6FD}$ -limit in the zero-net-current case is therefore

(A 12)

For the SPEC calculations presented in this paper, we have $\unicode[STIX]{x1D70C}_{a}^{2}=\unicode[STIX]{x1D6F9}_{a}/\unicode[STIX]{x1D6F9}_{\text{edge}}=0.3$ , and hence $\unicode[STIX]{x1D6FD}_{\text{lim}}$ is slightly lower (by a factor of 0.84) but very similar to the one obtained for Solov’ev profiles (see figures 5 and 6 for a comparison).


Bhattacharjee, A., Hayashi, T., Hegna, C. C., Nakajima, N. & Sato, T. 1995 Theory of pressure-induced islands and self-healing in three-dimensional toroidal magnetohydrodynamic equilibria. Phys. Plasmas 2 (3), 883.
Boozer, A. H. 2004 Physics of magnetically confined plasmas. Rev. Mod. Phys. 76 (4), 10711141.
Cary, J. R. & Hanson, J. D. 1986 Stochasticity reduction. Phys. Fluids 29, 2464.
Chirikov, B. V. 1979 A universal instability of many-dimensional oscillator systems. Phys. Rep. 52 (5), 263379.
Dennis, G. R., Hudson, S. R., Dewar, R. L. & Hole, M. J. 2013a The infinite interface limit of multiple-region relaxed magnetohydrodynamics. Phys. Plasmas 20 (3), 032509.
Dennis, G. R., Hudson, S. R., Terranova, D., Franz, P., Dewar, R. L. & Hole, M. J. 2013b Minimally constrained model of self-organized helical states in reversed-field pinches. Phys. Rev. Lett. 111 (5), 15.
Drevlak, M., Monticello, D. & Reiman, A. 2005 PIES free boundary stellarator equilibria with improved initial boundary conditions. Nucl. Fusion 45, 731740.
Freidberg, J. P. 2014 Ideal MHD. Cambridge University Press.
Greene, J. M. & Johnson, J. L. 1961 Determination of hydromagnetic equilibria. Phys. Fluids 4, 875.
Hanson, J. D. & Cary, J. R. 1984 Elimination of stochasticity in stellarators. Phys. Fluids 27, 767.
Hegna, C. C. 2012 Plasma flow healing of magnetic islands in stellarators. Phys. Plasmas 19 (5), 056101.
Helander, P. 2014 Theory of plasma confinement in non-axisymmetric magnetic fields. Rep. Prog. Phys. 77 (8), 087001.
Helander, P., Beidler, C. D., Bird, T. M., Drevlak, M., Feng, Y., Hatzky, R., Jenko, F., Kleiber, R., Proll, J. H. E., Turkin, Y. et al. 2012 Stellarator and tokamak plasmas: a comparison. Plasma Phys. Control. Fusion 54 (12), 124009.
Hirshman, S. P., Sanchez, R. & Cook, C. R. 2011 SIESTA: a scalable iterative equilibrium solver for toroidal applications. Phys. Plasmas 18 (6), 062504.
Hirshman, S. P. & Whitson, J. C. 1983 Steepest-descent moment method for three-dimensional magnetohydrodynamic equilibria. Phys. Fluids 26 (1), 35533568.
Hole, M. J., Hudson, S. R. & Dewar, R. L. 2007 Equilibria and stability in partially relaxed plasma–vacuum systems. Nucl. Fusion 47 (8), 746753.
Hudson, S. R. & Dewar, R. L. 1997 Manipulation of islands in a heliac vacuum field. Phys. Lett. A 226, 8592.
Hudson, S. R., Dewar, R. L., Dennis, G., Hole, M. J., McGann, M., Von Nessi, G. & Lazerson, S. 2012 Computation of multi-region relaxed magnetohydrodynamic equilibria. Phys. Plasmas 19 (11), 112502.
Hudson, S. R., Hole, M. J. & Dewar, R. L. 2007 Eigenvalue problems for Beltrami fields arising in a three-dimensional toroidal magnetohydrodynamic equilibrium problem. Phys. Plasmas 14 (5), 052505.
Kruskal, M. D. & Kulsrud, R. M. 1958 Equilibrium of a magnetically confined plasma in a toroid. Phys. Fluids 1 (4), 265.
Loizu, J., Hudson, S., Bhattacharjee, A. & Helander, P. 2015a Magnetic islands and singular currents at rational surfaces in three-dimensional magnetohydrodynamic equilibria. Phys. Plasmas 22 (2), 022501.
Loizu, J., Hudson, S. R., Bhattacharjee, A., Lazerson, S. & Helander, P. 2015b Existence of three-dimensional ideal-magnetohydrodynamic equilibria with current sheets. Phys. Plasmas 22 (9), 090704.
Loizu, J., Hudson, S. R., Helander, P., Lazerson, S. A. & Bhattacharjee, A. 2016a Pressure-driven amplification and penetration of resonant magnetic perturbations. Phys. Plasmas 23 (5), 055703.
Loizu, J., Hudson, S. R. & Nührenberg, C. 2016b Verification of the SPEC code in stellarator geometries. Phys. Plasmas 23 (11), 112505.
Meiss, J. D. 1992 Symplectic maps, variational principles and transport. Rev. Mod. Phys. 64 (3), 795848.
Miyamoto, K. 2005 Plasma Physics and Controlled Nuclear Fusion. Springer.
Narushima, Y., Watanabe, K. Y., Sakakibara, S., Narihara, K., Yamada, I., Suzuki, Y., Ohdachi, S., Ohyabu, N., Yamada, H. & Nakamura, Y. 2008 Dependence of spontaneous growth and suppression of the magnetic island on beta and collisionality in the LHD. Nucl. Fusion 48 (7), 075010.
Pedersen, T. S., Otte, M., Lazerson, S., Helander, P., Bozhenkov, S., Biedermann, C., Klinger, T., Wolf, R. C., Bosch, H. S.& The Wendelstein 7-X Team 2016 Confirmation of the topology of the Wendelstein 7-X magnetic field to better than 1:100 000. Nat. Commun. 7, 13493.
Suzuki, Y., Nakajima, N., Watanabe, K., Nakamura, Y. & Hayashi, T. 2006 Development and application of HINT2 to helical system plasmas. Nucl. Fusion 46 (11), L19L24.
Suzuki, Y., Watanabe, K. Y., Sakakibara, S., Nakajima, N. & Ohyabu, N. 2008 Theoretical studies of equilibrium beta limit in heliotron plasmas. In Proceedings of the 22nd IAEA Fusion Energy Conference, Geneva, pp. 919. Nuclear Fusion.
Taylor, J. B. 1974 Relaxation of toroidal plasma and generation of reverse magnetic fields. Phys. Rev. Lett. 33 (19), 11391141.
Taylor, J. B. 1986 Relaxation and magnetic reconnection in plasmas. Rev. Mod. Phys. 58 (3), 741763.
Waelbroeck, F. L. 1989 Current sheets and nonlinear growth of the $m=1$ kink-tearing mode. Phys. Fluids B 1 (12), 23722380.
Zakharov, L., Rogers, B. & Migliuolo, S. 1993 The theory of the early nonlinear stage of $m=1$ reconnection in tokamaks. Phys. Fluids B 5 (7), 24982505.