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.
Reference Helander, Beidler, Bird, Drevlak, Feng, Hatzky, Jenko, Kleiber, Proll and Turkin2012). In fact, magnetic surfaces are not guaranteed to exist in three-dimensional magnetohydrodynamics (MHD) equilibria without a continuous symmetry (Meiss Reference Meiss1992). While stellarators can be designed to possess magnetic surfaces in vacuum (Hanson & Cary Reference Hanson and Cary1984; Cary & Hanson Reference Cary and Hanson1986; Hudson & Dewar Reference Hudson and Dewar1997; Pedersen et al.
Reference Pedersen, Otte, Lazerson, Helander, Bozhenkov, Biedermann, Klinger, Wolf and Bosch2016), 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 Reference Drevlak, Monticello and Reiman2005; Suzuki et al.
Reference Suzuki, Watanabe, Sakakibara, Nakajima and Ohyabu2008).
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.
Reference Hudson, Dewar, Dennis, Hole, McGann, Von Nessi and Lazerson2012), although there are a few more ongoing parallel efforts (Suzuki et al.
Reference Suzuki, Nakajima, Watanabe, Nakamura and Hayashi2006; Hirshman, Sanchez & Cook Reference Hirshman, Sanchez and Cook2011). SPEC has been rigorously verified in axisymmetry (Hudson et al.
Reference Hudson, Dewar, Dennis, Hole, McGann, Von Nessi and Lazerson2012), in slightly perturbed configurations (Loizu et al.
Reference Loizu, Hudson, Bhattacharjee and Helander2015a
,Reference Loizu, Hudson, Bhattacharjee, Lazerson and Helander
b
, Reference Loizu, Hudson, Helander, Lazerson and Bhattacharjee2016a
) and more recently in stellarator geometries (Loizu, Hudson & Nührenberg Reference Loizu, Hudson and Nührenberg2016b
).
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 Reference Freidberg2014) 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 Reference Freidberg2014). 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,
$$\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
$$\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 Reference Helander2014) as
For example, for
$N_{p}=5$
we get 
.
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 Reference Taylor1974) 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. Reference Dennis, Hudson, Dewar and Hole2013a ; Loizu et al. Reference Loizu, Hudson, Bhattacharjee, Lazerson and Helander2015b ). 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 Reference Taylor1986), 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. Reference Dennis, Hudson, Terranova, Franz, Dewar and Hole2013b ) 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 Reference Hole, Hudson and Dewar2007; Hudson, Hole & Dewar Reference Hudson, Hole and Dewar2007). In essence, the energy functional is the same as in conventional ideal MHD equilibrium theory (Kruskal & Kulsrud Reference Kruskal and Kulsrud1958), 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
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D735}\times \boldsymbol{B}=\unicode[STIX]{x1D707}_{v}\boldsymbol{B}\quad \text{in the volumes} & \displaystyle\end{eqnarray}$$
$$\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
$$\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 Reference Miyamoto2005). As shown in § 4, 
decreases with
$\unicode[STIX]{x1D6FD}$
in the zero-net-current stellarator, thus amplifying the effect.
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
$$\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 Reference Miyamoto2005),
$\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 (Reference Freidberg2014). 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).
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).
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 (Reference Freidberg2014) 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,
and
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 Reference Hirshman and Whitson1983), 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.
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 Reference Waelbroeck1989; Zakharov, Rogers & Migliuolo Reference Zakharov, Rogers and Migliuolo1993).
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 Reference Freidberg2014)
where
$$\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’.
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 Reference Meiss1992). More precisely, we can evaluate the so-called box-counting dimension, or Hausdorff dimension,
$$\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
$$\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}}$
?
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$
.
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 Reference Chirikov1979). 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 Reference Boozer2004). 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.
Reference Bhattacharjee, Hayashi, Hegna, Nakajima and Sato1995; Narushima et al.
Reference Narushima, Watanabe, Sakakibara, Narihara, Yamada, Suzuki, Ohdachi, Ohyabu, Yamada and Nakamura2008; Hegna Reference Hegna2012) in the derivation of the equilibrium
$\unicode[STIX]{x1D6FD}$
-limit? Some new ideas are needed here.
Acknowledgements
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 Reference Freidberg2014) 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 Reference Greene and Johnson1961), 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 Reference Freidberg2014). The final result is a Grad–Shafranov-like partial differential equation, also called the overlap equation,
$$\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,
$$\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,
$$\begin{eqnarray}\displaystyle & \displaystyle \frac{\text{d}\unicode[STIX]{x1D6FD}}{\text{d}\unicode[STIX]{x1D713}}=C & \displaystyle\end{eqnarray}$$
$$\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
$$\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
where 
and
$\unicode[STIX]{x1D708}$
are defined in § 4. Finally, one can evaluate the rotational transform at the plasma edge, 
, with
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
$$\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
$$\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),
where
$$\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
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).
