2.2.1. Stellarator optimization

The basic concepts of the modern strategy for stellarator optimization were initiated in Garching (Chodura et al. Reference Chodura, Dommaschk, Herrnegger, Lotz, Nührenberg and Schlüter1981; Nührenberg & Zille Reference Nührenberg and Zille1988) in the 1980s. The state of a plasma with a definite pressure and current profile is fixed by a choice of the outermost surface (Bauer et al. Reference Bauer, Betancourt and Garabedian1984).

A standard expression for the outermost surface is the Garabedian representation (Bauer et al. Reference Bauer, Betancourt and Garabedian1984). A Garabedian surface is described in $(R_{g},{\it\varphi},Z_{g})$ cylindrical coordinates, $\boldsymbol{x}_{g}({\it\theta},{\it\varphi})=R_{g}\hat{\boldsymbol{R}}({\it\varphi})+Z_{g}\hat{\boldsymbol{Z}}$ , by expressing $R_{g}({\it\theta},{\it\varphi})$ and $Z_{g}({\it\theta},{\it\varphi})$ using complex notation,

(2.7) $$\begin{eqnarray}R_{g}+\text{i}Z_{g}=\mathop{\sum }_{mn}{\it\Delta}_{mn}\text{e}^{\text{i}(n{\it\varphi}-(m-1){\it\theta})},\end{eqnarray}$$

where $R_{g}$ , $Z_{g}$ and ${\it\Delta}_{mn}$ are all real. A given set of ${\it\Delta}_{mn}$ defines a surface.

An optimized stellarator is found by minimizing a quantity ${\it\chi}^{2}$ as a function of the shape parameters ${\it\Delta}_{mn}$ . Each physical quantity that is to be optimized has a target value, ${\rm\pi}_{{\it\alpha}}^{(target)}$ , and the value ${\rm\pi}_{{\it\alpha}}({\it\Delta}_{mn})$ achieved in a plasma state defined by a set of shape parameters ${\it\Delta}_{mn}$ . The definition of ${\it\chi}^{2}({\it\Delta}_{mn})$ is

(2.8) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\chi}^{2}({\it\Delta}_{mn})\equiv \mathop{\sum }_{{\it\alpha}}\frac{{\it\chi}_{{\it\alpha}}^{2}}{{\it\sigma}_{{\it\alpha}}^{2}}, & \displaystyle\end{eqnarray}$$

(2.9) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\chi}_{{\it\alpha}}^{2}({\it\Delta}_{mn})\equiv ({\rm\pi}_{{\it\alpha}}({\it\Delta}_{mn})-{\rm\pi}_{{\it\alpha}}^{(target)})^{2}. & \displaystyle\end{eqnarray}$$

${\it\sigma}_{{\it\alpha}}$

${\it\sigma}_{{\it\alpha}}$

The indeterminacy of the weighting ${\it\sigma}_{{\it\alpha}}$ of different properties is well known in everyday life. When searching for a new home, a safe neighbourhood, a short commute, comfortable size and low cost are all desirable features, but they cannot all be optimized. Despite agreement on desirable features, their relative weighting in determining the optimal home will vary from family to family.

Each of the ${\it\chi}_{{\it\alpha}}^{2}$ represents a simple property. Two that will be used as examples are (1) the effective ripple in the magnetic field strength on a magnetic surface (Nemov et al. Reference Nemov, Kasilov, Kernbichler and Heyn1999), which is a measure of the enhancement of neoclassical transport, and (2) the resonant current driven on a particular rational magnetic surface due to an ideal MHD perturbation. In a tokamak this current can be calculated by the IPEC code (Park et al. Reference Park, Schaffer, Menard and Boozer2007) and measures the extent to which that perturbation would drive an island if the resonant current were allowed to relax. A related calculation has been done for stellarators (Nührenberg et al. Reference Nührenberg, Boozer and Hudson2009). For both examples, the target function is zero.

A third example of a ${\it\chi}_{{\it\alpha}}^{2}$ is a method for determining the most efficient external magnetic field consistent with the plasma having a given set of physics properties. Define functions ${\rm\pi}_{{\it\alpha}}=\text{e}^{2t_{{\it\alpha}}}{\it\Phi}_{{\it\alpha}}^{2}$ . The ${\it\Phi}_{{\it\alpha}}$ are obtained from the external field required to support the equilibrium, $\boldsymbol{B}_{x}\boldsymbol{\cdot }\hat{\boldsymbol{n}}({\it\Delta}_{mn})$ , by an expansion in the $\boldsymbol{F}_{{\it\alpha}}^{(t)}$ , which are the left eigenvectors of the transfer matrix $\unicode[STIX]{x1D64F}$ . To exclude highly inefficient magnetic fields, which have an index greater than $i_{c}$ , a ${\rm\pi}_{{\it\alpha}}^{2}({\it\Delta}_{mn})$ can be defined by the square of the amplitude of the field components that are to be excluded, ${\rm\pi}_{{\it\alpha}}^{2}=\text{e}^{2t_{i_{c}}}\sum _{i>i_{c}}{\it\Phi}_{i}^{2}$ . The code STELLOPT/COILOPT combines the stellarator physics and coil optimization in a single code (Strikler et al. Reference Strikler, Hirshman, Spong, Cole, Lyon, Nelson, Williamson and Ware2004) but does so by optimizing coefficients in a Fourier representation of the coils. STELLOPT/COILOPT could be modified to optimize in the space of complete and non-redundant set external magnetic fields.

The standard method for calculating stellarator equilibria is the VMEC code (Hirshman & Whitson Reference Hirshman and Whitson1983), and that code coupled with one of many insights tied to the name of Shafranov, the virtual casing principle (Shafranov & Zakharov Reference Shafranov and Zakharov1972), gives the external magnetic field $\boldsymbol{B}_{x}\boldsymbol{\cdot }\hat{\boldsymbol{n}}({\it\Delta}_{mn})$ required to support a plasma state defined by a set of ${\it\Delta}_{mn}$ . The normal field $\boldsymbol{B}_{x}\boldsymbol{\cdot }\hat{n}$ can be expanded in the left eigenfunctions of the transfer matrix to obtain dimensionless amplitude coefficients,

(2.10) $$\begin{eqnarray}\displaystyle a_{i}\equiv \frac{1}{B_{0}}\oint \boldsymbol{B}_{x}\boldsymbol{\cdot }\hat{n}F_{i}^{(t)}w\,\text{d}{\it\theta}\,\text{d}{\it\varphi}=\frac{{\it\Phi}_{i}}{B_{0}A}. & & \displaystyle\end{eqnarray}$$

2.2.2. Sensitivity of an optimum

Once an optimum stellarator configuration is found, the ${\it\Delta}_{mn}$ can be varied by a small amount, ${\it\delta}{\it\Delta}_{mn}$ , away from their optimal values. A small change in the shape of the plasma, a given ${\it\delta}{\it\Delta}_{mn}$ , implies a small change in the external magnetic field required to support the plasma equilibrium, which can be expressed as ${\it\delta}\boldsymbol{B}_{x}\boldsymbol{\cdot }\hat{\boldsymbol{n}}$ on a fixed toroidal surface just outside the plasma. Writing the external magnetic field as ${\it\delta}\boldsymbol{B}_{x}\boldsymbol{\cdot }\hat{\boldsymbol{n}}=w\boldsymbol{f}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D731}/A$ , (2.5), changes away from the optimum have the form

(2.11) $$\begin{eqnarray}{\it\delta}{\it\chi}^{2}=\unicode[STIX]{x1D731}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\boldsymbol{\cdot }\unicode[STIX]{x1D731}.\end{eqnarray}$$

$\unicode[STIX]{x1D64E}$ is the sensitivity matrix and is what is known as the Hessian matrix in optimization theory.

The sensitivity matrix is symmetric, its eigenvalues have real positive values, $\text{e}^{s_{i}}$ , and its eigenvectors $F_{i}^{(s)}({\it\theta},{\it\varphi})$ can be ordered by the magnitude of the associated eigenvalue – the eigenvector with the largest eigenvalue first.

2.2.3. Sensitivity of individual ${\it\chi}_{{\it\alpha}}^{2}$

Each ${\it\chi}_{{\it\alpha}}^{2}$ , which is associated with a simple property, has an eigenfunction $f_{{\it\alpha}}({\it\theta},{\it\varphi})$ , which gives the unique distribution of external magnetic field ${\it\delta}\boldsymbol{B}_{x}\boldsymbol{\cdot }\hat{\boldsymbol{n}}\propto wf_{{\it\alpha}}$ that produces a small change in that property – whether small changes are quadratic or linear in the perturbation amplitude.

The dependence of ${\it\chi}_{{\it\alpha}}^{2}$ on the perturbation amplitude can be quadratic, as is the case for the singular current on a particular rational surface in a tokamak. The singular current is proportional to the driving perturbation and has the form $I_{{\it\alpha}}=\sum _{i}\unicode[STIX]{x1D61D}_{i}{\it\Phi}_{i}$ , where ${\it\alpha}$ denotes a particular rational surface. The matrix $\unicode[STIX]{x1D61D}_{i}$ is a matrix vector, which has only one non-zero singular value and that singular value is associated with a particular eigenfunction $f_{{\it\alpha}}({\it\theta},{\it\varphi})$ .

Even near a minimum of ${\it\chi}^{2}$ , the dependence of a particular ${\it\chi}_{{\it\alpha}}^{2}$ can be linear in the perturbation amplitude, ${\it\chi}_{{\it\alpha}}^{2}=\sum _{i}\unicode[STIX]{x1D61D}_{i}{\it\Phi}_{i}$ . This occurs when a change in the ${\it\Delta}_{mn}$ about the optimum of necessity makes one property better and another worse. The matrix row vector $\unicode[STIX]{x1D61D}_{i}$ has only one non-zero singular value, which is associated with a particular eigenfunction $f_{{\it\alpha}}({\it\theta},{\it\varphi})$ .

The eigenfunctions $f_{{\it\alpha}}({\it\theta},{\it\varphi})$ from different ${\it\chi}_{{\it\alpha}}^{2}$ are not generally independent. Indeed, the symmetric matrix $\unicode[STIX]{x1D64A}$ of overlap integrals

(2.12) $$\begin{eqnarray}\unicode[STIX]{x1D616}_{{\it\alpha}{\it\beta}}\equiv \oint f_{{\it\alpha}}f_{{\it\beta}}w\,\text{d}{\it\theta}\,\text{d}{\it\varphi}\end{eqnarray}$$

must have an eigenvalue that is zero for ${\it\chi}_{{\it\alpha}}^{2}$ to have a linear dependence on the perturbation to a plasma state in which ${\it\chi}^{2}$ is at a minimum. When the $f_{{\it\alpha}}$ from different properties are orthogonal, $\unicode[STIX]{x1D616}_{{\it\alpha}{\it\beta}}={\it\delta}_{{\it\alpha}{\it\beta}}$ , and all eigenvalues are unity. The eigenvalues of $\unicode[STIX]{x1D64A}$ will be written as $\text{e}^{-o_{i}}$ and the eigenfunctions as $F_{i}^{(o)}({\it\theta},{\it\varphi})$ . The coupling between two properties, ${\it\alpha}$ and ${\it\beta}$ , is given by the overlap integral $\unicode[STIX]{x1D616}_{{\it\alpha}{\it\beta}}$ between their eigenfunctions.

The number of non-zero eigenvalues of $\unicode[STIX]{x1D64A}$ is the number of the functions $F_{i}^{(o)}$ that are independent. When an eigenvalue of $\unicode[STIX]{x1D64A}$ is small, properties coupled in that eigenfunction are tightly coupled – one cannot be changed without changing another. The weightings ${\it\sigma}_{{\it\alpha}}$ in the optimization only determine where in ${\it\Delta}_{mn}$ space the ${\it\chi}^{2}$ optimized plasma state is located. The choice of the ${\it\sigma}_{{\it\alpha}}$ has no direct effect on the coupling between properties, and in particular what types of further optimization of certain properties are possible, in the neighbourhood of that plasma state.

An area in which improved optimization may be important is the confinement of energetic alpha particles. Power loss is not the issue: roughly half the alphas are on passing-particle trajectories, which are well confined when magnetic surfaces exist. The issue is controlling the damage to the wall surrounding the plasma from bombardment by energetic alpha particles. An extremely low rate of alpha loss may not be required. It may be possible to steer the lost alphas into a pool of liquid metal, which would eliminate the material-damage issue.

An important question is how collections of properties can be independently controlled. This question is addressed by the singular value decomposition of the matrix

(2.13) $$\begin{eqnarray}\unicode[STIX]{x1D60A}_{ij}\equiv \oint \text{e}^{-t_{i}}F_{i}^{(t)}\text{e}^{-o_{j}}F_{j}^{(o)}w\,\text{d}{\it\theta}\,\text{d}{\it\varphi}.\end{eqnarray}$$

When a singular value $\text{e}^{-c_{\ell }}$ is small compared to unity, the right eigenvector ${\mathcal{F}}_{k}^{(c)}$ , which is a linear combination of the $F_{j}^{(o)}$ can only be driven with low efficiency. When $\text{e}^{-c_{\ell }}$ is small compared to other singular values, extreme accuracy is required in the ratio of currents in different coils in order to separately drive this particular collection of properties.

2.3.1. Flexibility for many configurations

Once a desirable plasma state is found that requires only efficiently produced external magnetic fields, the optimization of the design for supporting $p_{s}\gg 1$ plasma states should be carried out. This can be done by finding the minimal set of efficiently produced distributions $n_{x}$ that are required to support a large set of plasma states. The external magnetic fields required to support each of these plasma states can be designated on an envelope toroidal surface, which is defined so it lies outside the plasmas in all of the states but close to a plasma in at least one of the states. Each plasma state $p$ has a set of distributions ${\mathcal{F}}_{i}^{p}$ of the efficiency-ordered fields on the envelope surface required to support it, with each at an amplitude $a_{i}^{p}$ , (2.10). The issue is how many coils that are independently controllable, $n_{x}$ , are required to provide most of the required magnetic field to support the $p_{s}$ plasma states. The remaining fields are either provided by trim coils or perhaps some parts of the remaining fields are sufficiently small not to require control.

The error in the required magnetic field to support the $p\text{th}$ plasma state which is defined by expansion coefficients $\boldsymbol{a}_{p}$ , (2.10), can be written as

(2.14) $$\begin{eqnarray}{\it\delta}\boldsymbol{B}^{(p)}\boldsymbol{\cdot }\hat{\boldsymbol{n}}=w(\boldsymbol{a}_{p}-\unicode[STIX]{x1D63C}\boldsymbol{\cdot }\boldsymbol{D})^{\dagger }\boldsymbol{\cdot }\boldsymbol{{\mathcal{F}}}^{(t)},\end{eqnarray}$$

where each element in the matrix column vector $\boldsymbol{D}$ represents the B-normal distribution on the envelope surface that is produced by a single coil. The matrix $\unicode[STIX]{x1D608}_{ij}$ is the amplitude of the $i\text{th}$ efficiency-ordered normal field distribution produced by the $j\text{th}$ element of $\boldsymbol{D}$ . With the basic coils, one wishes to minimize the additional fields that must be supplied by auxiliary coils. That is the coils should be chosen so $\boldsymbol{D}$ can be adjusted to minimize the root-mean-square (RMS) error in the normal magnetic field on the control surface. A weighting matrix, such as the sensitivity matrix, can be added, but for simplicity it will be assumed that one wishes to minimize the (RMS) field error ${\mathcal{E}}_{p}$ in each of the $p_{s}$ plasma states, where

(2.15) $$\begin{eqnarray}{\mathcal{E}}_{p}^{2}\equiv (\boldsymbol{a}_{p}-\unicode[STIX]{x1D63C}\boldsymbol{\cdot }\boldsymbol{D})^{\dagger }\boldsymbol{\cdot }(\boldsymbol{a}_{p}-\unicode[STIX]{x1D63C}\boldsymbol{\cdot }\boldsymbol{D}).\end{eqnarray}$$

When $\boldsymbol{D}$ is chosen to minimize ${\mathcal{E}}_{p}^{2}$ for the $p\text{th}$ plasma state,

(2.16) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{D}=\unicode[STIX]{x1D63C}_{s}^{-1}\boldsymbol{\cdot }\boldsymbol{a}_{p}, & \displaystyle\end{eqnarray}$$

(2.17) $$\begin{eqnarray}\displaystyle & \displaystyle {\mathcal{E}}_{p}^{2}=\boldsymbol{a}_{p}^{\dagger }\boldsymbol{\cdot }(\unicode[STIX]{x1D64B}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D64B})\boldsymbol{\cdot }\boldsymbol{a}_{p},\quad \text{where} & \displaystyle\end{eqnarray}$$

(2.18) $$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x1D64B}\equiv \mathbf{1}-\unicode[STIX]{x1D63C}\boldsymbol{\cdot }\unicode[STIX]{x1D63C}_{s}^{-1}, & \displaystyle\end{eqnarray}$$

$\mathbf{1}$

$\unicode[STIX]{x1D63C}_{s}^{-1}$

$\unicode[STIX]{x1D63C}$

The pseudo-inverse $\unicode[STIX]{x1D63C}_{s}^{-1}$ is defined by SVD, which means the form $\unicode[STIX]{x1D63C}=\unicode[STIX]{x1D650}\boldsymbol{\cdot }\unicode[STIX]{x1D668}_{A}\boldsymbol{\cdot }\unicode[STIX]{x1D651}^{\dagger }$ , where $\unicode[STIX]{x1D650}^{-1}=\unicode[STIX]{x1D650}^{\dagger }$ , $\unicode[STIX]{x1D651}^{-1}=\unicode[STIX]{x1D651}^{\dagger }$ , and $\unicode[STIX]{x1D668}_{A}$ is diagonal with all elements greater than or equal to zero. An SVD exists for essentially any matrix. The pseudo-inverse $\unicode[STIX]{x1D63C}_{s}^{-1}\equiv \unicode[STIX]{x1D650}\boldsymbol{\cdot }\unicode[STIX]{x1D668}_{{\mathcal{A}}}^{-1}\boldsymbol{\cdot }\unicode[STIX]{x1D651}^{\dagger }$ , where for any diagonal elements $(s_{A})_{i}\geqslant s_{min}$ , the element $(\unicode[STIX]{x1D668}_{{\mathcal{A}}}^{-1})_{i}\equiv 1/(s_{A})_{i}$ , but if $(s_{A})_{i}<s_{min}$ , then $(\unicode[STIX]{x1D668}_{{\mathcal{A}}}^{-1})_{i}\equiv 0$ . When $s_{min}$ is made smaller, more equilibria can be supported by a given set of coils, but when $s_{min}$ is too small, the currents may become too large and accuracy requirements on the currents may be too great to be feasible.

Let $i_{x}$ be the number of efficiency-ordered distributions required to support all $p_{s}$ plasma states, then $\unicode[STIX]{x1D63C}$ is an $i_{x}\times i_{D}$ matrix, where $i_{D}$ is the number of elements in $\boldsymbol{D}$ , which is the number of independent drivable fields in the main coil set. The matrix $\unicode[STIX]{x1D63C}$ must have $i_{x}-i_{D}$ singular values that are zero, when $i_{x}>i_{D}$ . What one wishes to do is to choose the components of $\unicode[STIX]{x1D63C}$ for fixed $i_{x}$ and $i_{D}$ to make the maximum over $p$ of ${\mathcal{E}}_{p}$ as small as possible. ${\mathcal{E}}_{p}$ represents the strength of the field that must be produced by trim coils. When ${\mathcal{E}}_{p}$ is too large, either $i_{D}$ must be made larger or the range of plasma states reduced.

2.3.2. Error-field mitigation

Reducing deviations of the externally provided magnetic field from the design field by extreme precision of construction increases the cost, time, and the uncertainty of building fusion experiments. These three factors can be minimized by including field error mitigation in the design of an experiment. When the magnetic field errors are too large, too many trim coils are required for practical control and too much current can be required in the trim coils for a practical design. The basic theory of error-field mitigation given here follows § 12.3 of Boozer (Reference Boozer2015) and is related to the theory of machine flexibility, § 2.3.1. The theory of plasma control through modification of the magnetic field distributions to which plasma is sensitive is essentially identical to the theory of error-field mitigation.

The error in the magnetic field ${\it\delta}\boldsymbol{B}_{m}$ produced by the main equilibrium coils is determined in the plasma region by the normal component ${\it\delta}\boldsymbol{B}_{m}\boldsymbol{\cdot }\hat{\boldsymbol{n}}$ to a control surface just on the plasma side of the coils. That is,

(2.19) $$\begin{eqnarray}{\it\delta}\boldsymbol{B}_{m}\boldsymbol{\cdot }\hat{\boldsymbol{n}}=w\mathop{\sum }_{i}{\it\delta}{\it\Phi}_{i}^{(m)}f_{i}({\it\theta}_{c},{\it\varphi}_{c}),\end{eqnarray}$$

where the $f_{i}({\it\theta}_{c},{\it\varphi}_{c})$ are an orthonormal set of functions on the coil control surface.

The perturbation to the externally produced field on plasma boundary is then

(2.20) $$\begin{eqnarray}{\it\delta}\unicode[STIX]{x1D731}_{x}=\unicode[STIX]{x1D64F}\boldsymbol{\cdot }{\it\delta}\unicode[STIX]{x1D731}_{m}+\unicode[STIX]{x1D648}\boldsymbol{\cdot }\boldsymbol{J},\end{eqnarray}$$

where $\unicode[STIX]{x1D64F}$ is a transfer matrix, which relates the magnetic perturbations on the plasma surface to the perturbations near the coils, and $\unicode[STIX]{x1D648}$ is a mutual inductance matrix that gives the effect of currents $\boldsymbol{J}$ in the control coils.

The plasma sensitivity is measured by ${\it\delta}\unicode[STIX]{x1D731}_{x}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\boldsymbol{\cdot }{\it\delta}\unicode[STIX]{x1D731}_{x}$ . The currents $\boldsymbol{J}$ in the trim coils should be chosen to minimize ${\it\delta}\unicode[STIX]{x1D731}_{x}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\boldsymbol{\cdot }{\it\delta}\unicode[STIX]{x1D731}_{x}$ . The answer is $\unicode[STIX]{x1D648}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\boldsymbol{\cdot }\unicode[STIX]{x1D648}\boldsymbol{\cdot }\boldsymbol{J}=-\unicode[STIX]{x1D648}^{\dagger }\boldsymbol{\cdot }S\boldsymbol{\cdot }\unicode[STIX]{x1D64F}\boldsymbol{\cdot }{\it\delta}\unicode[STIX]{x1D731}_{m}$ . The optimal trim-coil currents are then

(2.21) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{J}=-\boldsymbol{{\mathcal{C}}}_{k}\boldsymbol{\cdot }{\it\delta}\unicode[STIX]{x1D731}_{m},\quad \text{ where} & \displaystyle\end{eqnarray}$$

(2.22) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\mathcal{C}}}_{k}\equiv (\unicode[STIX]{x1D648}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\boldsymbol{\cdot }\unicode[STIX]{x1D648})_{k}^{-1}\boldsymbol{\cdot }\unicode[STIX]{x1D648}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\boldsymbol{\cdot }\unicode[STIX]{x1D64F}, & \displaystyle\end{eqnarray}$$

$k$

$k$

$\unicode[STIX]{x1D648}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\boldsymbol{\cdot }\unicode[STIX]{x1D648}$

(2.23) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\delta}\unicode[STIX]{x1D731}_{x}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\boldsymbol{\cdot }{\it\delta}\unicode[STIX]{x1D731}_{x}={\it\delta}\unicode[STIX]{x1D731}_{m}^{\dagger }\boldsymbol{\cdot }\boldsymbol{{\mathcal{S}}}_{k}\boldsymbol{\cdot }{\it\delta}\unicode[STIX]{x1D731}_{m}, & \displaystyle\end{eqnarray}$$

(2.24) $$\begin{eqnarray}\displaystyle & \displaystyle \boldsymbol{{\mathcal{S}}}_{k}\equiv (\unicode[STIX]{x1D64F}-\unicode[STIX]{x1D648}\boldsymbol{\cdot }\boldsymbol{{\mathcal{C}}}_{k})^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\boldsymbol{\cdot }(\unicode[STIX]{x1D64F}-\unicode[STIX]{x1D648}\boldsymbol{\cdot }\boldsymbol{{\mathcal{C}}}_{k}). & \displaystyle\end{eqnarray}$$

The number of singular values $k$ retained in $\boldsymbol{{\mathcal{C}}}_{k}$ can be interpreted as the effective number of trim coils that can be used for error-field mitigation. The condition number of $\boldsymbol{{\mathcal{C}}}_{k}$ , which is the ratio of the largest to the smallest retained singular values, should be as small as possible since it measures how accurately the coil currents must be chosen to independently control $k$ error fields. The required construction accuracy is generally set by having sufficiently few field errors that require mitigation. That is, ${\it\delta}\unicode[STIX]{x1D731}_{x}^{\dagger }\boldsymbol{\cdot }\unicode[STIX]{x1D64E}\boldsymbol{\cdot }{\it\delta}\unicode[STIX]{x1D731}_{x}$ must be sufficiently small for a practical value of $k$ , where $k$ is both the effective number of control coils and the number of error fields requiring correction.

The singular values of $\boldsymbol{{\mathcal{C}}}_{k}$ generally increase exponentially from one to the next, so the allowed number of errors that can require correction has a relatively sharp upper limit. It should be noted that neither $\boldsymbol{{\mathcal{S}}}_{k}$ nor $\boldsymbol{{\mathcal{C}}}_{k}$ has a dependence on the actual magnetic field errors, so an error-field mitigation system can be designed without knowing the as-built coil errors. Once an error-field mitigation system is designed, the largest tolerable as-built errors can be determined both in magnitude and form.

3.2.1. Resonant divertors

In a resonant divertor, the external magnetic field has a resonance on an outer magnetic surface that has a rotational transform ${\it\iota}_{d}$ . In a tokamak, the resonance is ${\it\iota}_{d}=0$ , and the magnetic island that forms is axisymmetric. The W7-X stellarator has a resonant divertor; the resonance is ${\it\iota}_{d}=5/5$ , which means that the external magnetic field has five periods toroidally and has a component that has five periods poloidally. The $5/5$ component splits the ${\it\iota}=5/5$ surface to form a magnetic island, which diverts the magnetic field lines (Strumberger Reference Strumberger1996; Feng et al. Reference Feng, Sardei, Grigull, McCormick, Kisslinger and Reiter2006). The same divertor plate locations can be used if the resonant surface is chosen to be at ${\it\iota}_{d}=5/6$ or ${\it\iota}_{d}=5/4$ .

The advantages of a magnetic-island or resonant divertor are that (1) the divertor region is compact staying very close to the body of the plasma and (2) the divertor has a long connection length because of the slowness with which magnetic field lines encircle a narrow island.

The resonant divertor has the following two disadvantages.

1. (1) The formation of a magnetic island requires a fixed resonant rotational transform ${\it\iota}_{d}$ be maintained at the plasma edge, which is difficult when the plasma evolution naturally has a large variation in the net toroidal current.

2. (2) The intersection points of the magnetic field lines with the surrounding structure in a resonant divertor are separated by only a small distance from the main plasma body, despite the many toroidal transits required by a magnetic field line to go from the plasma edge to the surrounding structure.

3.2.2. Non-resonant divertors

A non-resonant divertor is produced by using the external magnetic field to force a sharp edge on the plasma surface. The axisymmetric sharp edge on the surface of a tokamak, the X-line, cannot be crossed by a magnetic field line, which forces the rotational transform to be zero on the separatrix. This is despite the fact that the toroidal current enclosed by the separatrix produces a poloidal field that, away from the X-line, would be associated with a rotational transform ${\it\iota}\sim 1/3$ .

The magnetic field lines near a sharp edge, even an edge that is just a segment of a line, resemble those near at tokamak X-line. That is, magnetic field lines are carried off to the walls, and non-resonant terms in the magnetic field line Hamiltonian, when sufficiently large, can produce a divertor configuration. Non-resonant divertors have been discussed for quasi-axisymmetric stellarators (Mioduszewski et al. Reference Mioduszewski, Owen, Spong, Fenstermacher, Koniges, Rognlien, Umansky, Grossman and Kugel2007; Mau et al. Reference Mau, Kaiser, Grossman, Raffray, Want, Lyon, Maingi, Ku and Zarnstorff2008).

Few studies have been made of non-resonant stellarator divertors, but, as in a tokamak, a non-resonant divertor can be insensitive to the net plasma current, and the diverted field lines can intersect the divertor structures far from the main plasma body.

Sharp edges are a requirement for a non-resonant divertor but may or may not be present in a resonant divertor. Sharp edges are present in the W7-X resonant divertor. Nevertheless, an island for a resonant divertor can be located on what would otherwise be a smooth though rational magnetic surface. That is, the location of the place from which field lines are diverted can be determined either by the global geometry defining sharp edges, as in a non-resonant divertor, or by a particular rational value of the rotational transform ${\it\iota}_{d}$ , as in resonant divertor.

3.5.1. Position vector for resonant divertors

The behaviour of the boundary between confined and non-confined magnetic field lines in a resonant divertor is essentially described by the magnetic field line Hamiltonian, (3.1), and only a moderate number of the constants ${\it\alpha}_{mn}$ need be non-zero. The position vector can be chosen to be an extrapolated Garabedian surface, (2.7), $\boldsymbol{x}_{p}({\it\psi}_{t},{\it\theta},{\it\varphi})=R\hat{\boldsymbol{R}}+Z\hat{\boldsymbol{Z}}$ , where

(3.4) $$\begin{eqnarray}R+\text{i}Z=\mathop{\sum }_{m,n}{\it\Delta}_{mn}\left(\frac{{\it\psi}_{t}}{{\it\psi}_{c}}\right)^{m/2}\text{e}^{\text{i}(n{\it\varphi}-(m-1){\it\theta})}.\end{eqnarray}$$

3.5.2. Position vector for non-resonant divertors

The position vector and the Garabedian surface, (2.7), have fundamentally different roles in the modelling of non-resonant versus resonant divertors. For non-resonant divertors, the position vector can be chosen to have a simple form. The Garabedian surface determines the constants that appear in the magnetic field line Hamiltonian. A field line of the Hamiltonian should coincide with the Garabedian surface. That is, $\boldsymbol{x}_{p}({\it\psi}_{t}={\it\psi}_{g},{\it\theta},{\it\varphi})$ must equal $\boldsymbol{x}_{g}({\it\theta},{\it\varphi})$ for some choice of the function ${\it\psi}_{g}({\it\theta},{\it\varphi})$ .

A consistent position vector can be defined using ordinary $(R,{\it\varphi},Z)$ cylindrical coordinates, $\boldsymbol{x}_{p}({\it\psi}_{t},{\it\theta},{\it\varphi})=R({\it\psi}_{t},{\it\theta},{\it\varphi})\hat{\boldsymbol{R}}({\it\varphi})+Z({\it\psi}_{t},{\it\theta},{\it\varphi})\hat{\boldsymbol{Z}}$ . The $R$ and $Z$ positions are real but can be defined by the complex form

(3.5) $$\begin{eqnarray}R+\text{i}Z=R_{a}({\it\varphi})+\text{i}Z_{a}({\it\varphi})+{\it\Delta}_{0,0}\sqrt{\frac{{\it\psi}_{t}}{{\it\psi}_{c}}}\text{e}^{\text{i}({\it\theta}+{\it\delta}{\it\theta})},\end{eqnarray}$$

where the constant ${\it\psi}_{c}$ is a characteristic value of the toroidal flux and ${\it\Delta}_{0,0}$ , a constant from the Garabedian representation, is the minor radius. The function ${\it\delta}{\it\theta}({\it\psi}_{t},{\it\theta},{\it\varphi})$ is required for the equation $\boldsymbol{x}_{p}({\it\psi}_{g},{\it\theta},{\it\varphi})=\boldsymbol{x}_{g}({\it\theta},{\it\varphi})$ to have a solution for the toroidal magnetic flux ${\it\psi}_{t}={\it\psi}_{g}({\it\theta},{\it\varphi})$ on the Garabedian surface, (2.7). The equation for the axis of the poloidal angle is

(3.6) $$\begin{eqnarray}R_{a}+\text{i}Z_{a}=\mathop{\sum }_{n=0}^{\infty }{\it\Delta}_{1,n}\text{e}^{\text{i}n{\it\varphi}}.\end{eqnarray}$$

3.6.1. Expressions for ${\it\psi}_{g}$ and ${\it\delta}{\it\theta}$

The position vector defined in (3.5) is consistent with the Garabedian position vector, (2.7), only when the functions ${\it\psi}_{g}({\it\theta},{\it\varphi})$ and ${\it\delta}{\it\theta}({\it\theta},{\it\varphi})$ have a certain form.

The equation $\boldsymbol{x}_{p}({\it\psi}_{g},{\it\theta},{\it\varphi})=\boldsymbol{x}_{g}({\it\theta},{\it\varphi})$ , where ${\it\psi}_{t}={\it\psi}_{g}({\it\theta},{\it\varphi})$ on the Garabedian surface, implies

(3.7) $$\begin{eqnarray}\sqrt{\frac{{\it\psi}_{g}}{{\it\psi}_{c}}}\text{e}^{\text{i}{\it\delta}{\it\theta}}=\mathop{\sum }_{m\neq 0}\frac{{\it\Delta}_{m,n}}{{\it\Delta}_{0,0}}\text{e}^{\text{i}(n{\it\varphi}-m{\it\theta})}.\end{eqnarray}$$

This equation for ${\it\psi}_{g}({\it\theta},{\it\varphi})$ and ${\it\delta}{\it\theta}({\it\theta},{\it\varphi})$ can be solved by letting

(3.8) $$\begin{eqnarray}\displaystyle & \displaystyle C({\it\theta},{\it\varphi})\equiv 1+{\it\Sigma}^{\prime }\frac{{\it\Delta}_{m,n}}{{\it\Delta}_{0,0}}\cos (n{\it\varphi}-m{\it\theta}), & \displaystyle\end{eqnarray}$$

(3.9) $$\begin{eqnarray}\displaystyle & \displaystyle S({\it\theta},{\it\varphi})\equiv {\it\Sigma}^{\prime }\frac{{\it\Delta}_{m,n}}{{\it\Delta}_{0,0}}\sin (n{\it\varphi}-m{\it\theta}), & \displaystyle\end{eqnarray}$$

${\it\Sigma}^{\prime }$

$(m,n)$

$(0,0)$

$(1,n)$

(3.10) $$\begin{eqnarray}\displaystyle & \displaystyle \cos {\it\delta}{\it\theta}=\frac{C}{\sqrt{C^{2}+S^{2}}}, & \displaystyle\end{eqnarray}$$

(3.11) $$\begin{eqnarray}\displaystyle & \displaystyle \sin {\it\delta}{\it\theta}=\frac{S}{\sqrt{C^{2}+S^{2}}}, & \displaystyle\end{eqnarray}$$

(3.12) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\psi}_{g}({\it\theta},{\it\varphi})=(C^{2}+S^{2}){\it\psi}_{c}. & \displaystyle\end{eqnarray}$$

The Fourier expansions of ${\it\psi}_{g}/{\it\psi}_{c}$ and ${\it\delta}{\it\theta}$ have the general forms

(3.13) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{{\it\psi}_{g}}{{\it\psi}_{c}}=1+2{\it\delta}({\it\theta},{\it\varphi}), & \displaystyle\end{eqnarray}$$

(3.14) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\delta}\equiv {\it\Sigma}^{\prime }{\it\delta}_{m,n}\cos (n{\it\varphi}-m{\it\theta}), & \displaystyle\end{eqnarray}$$

(3.15) $$\begin{eqnarray}\displaystyle & \displaystyle {\it\delta}{\it\theta}={\it\Sigma}^{\prime }{\it\theta}_{mn}\sin (n{\it\varphi}-m{\it\theta}). & \displaystyle\end{eqnarray}$$

$|{\it\Delta}_{mn}/{\it\Delta}_{00}|$

${\it\Sigma}^{\prime }$

${\it\delta}_{mn}$

${\it\theta}_{mn}$

${\it\Delta}_{mn}/{\it\Delta}_{00}$

3.6.2. Consistency of the Hamiltonian with the Garabedian surface

On the Garabedian surface, (2.7), the toroidal flux coordinate that appears in the standard position vector is ${\it\psi}_{t}={\it\psi}_{g}({\it\theta},{\it\varphi})$ . The constraint that the Garabedian surface be a magnetic surface implies

(3.16) $$\begin{eqnarray}\frac{\text{d}{\it\psi}_{t}}{\text{d}{\it\varphi}}=\frac{\partial {\it\psi}_{g}}{\partial {\it\theta}}\frac{\text{d}{\it\theta}}{\text{d}{\it\varphi}}+\frac{\partial {\it\psi}_{g}}{\partial {\it\varphi}}\end{eqnarray}$$

when evaluated on the Garabedian surface ${\it\psi}_{t}={\it\psi}_{g}$ . This equation and Hamilton’s equations imply

(3.17) $$\begin{eqnarray}\frac{\partial {\it\psi}_{g}}{\partial {\it\theta}}\left(\frac{\partial {\it\psi}_{p}}{\partial {\it\psi}_{t}}\right)_{g}+\left(\frac{\partial {\it\psi}_{p}}{\partial {\it\theta}}\right)_{g}+\frac{\partial {\it\psi}_{g}}{\partial {\it\varphi}}=0.\end{eqnarray}$$

The required derivatives of the magnetic field line Hamiltonian are

(3.18) $$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\partial {\it\psi}_{p}}{\partial {\it\psi}_{t}}\right)_{g}=\bar{{\it\iota}}_{g}+\mathop{\sum }_{m,n}\frac{m}{2}{\it\alpha}_{mn}(1+2{\it\delta})^{(m-2)/2}\cos (\cdots \,), & \displaystyle\end{eqnarray}$$

(3.19) $$\begin{eqnarray}\displaystyle & \displaystyle \bar{{\it\iota}}_{g}=\bar{{\it\iota}}\big({\it\psi}_{a}(1+2{\it\delta})\big), & \displaystyle\end{eqnarray}$$

(3.20) $$\begin{eqnarray}\displaystyle & \displaystyle \left(\frac{\partial {\it\psi}_{p}}{\partial {\it\theta}}\right)_{g}={\it\psi}_{a}\mathop{\sum }_{m,n}m{\it\alpha}_{mn}(1+2{\it\delta})^{m/2}\sin (\cdots \,),\hspace{14.45377pt} & \displaystyle\end{eqnarray}$$

(3.21) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\psi}_{t}^{g}}{\partial {\it\theta}}=2{\it\psi}_{a}{\it\Sigma}^{\prime }m{\it\delta}_{mn}\sin (n{\it\varphi}-m{\it\theta}), & \displaystyle\end{eqnarray}$$

(3.22) $$\begin{eqnarray}\displaystyle & \displaystyle \frac{\partial {\it\psi}_{t}^{g}}{\partial {\it\varphi}}=-2{\it\psi}_{a}{\it\Sigma}^{\prime }n{\it\delta}_{mn}\sin (n{\it\varphi}-m{\it\theta}), & \displaystyle\end{eqnarray}$$

$\cos (\cdots \,)=\cos (n{\it\varphi}-m{\it\theta})$

The solution to (3.17) for a set of ${\it\alpha}_{mn}$ can be found by an optimization technique. Optimization algorithms are made more efficient when given an approximate optimum. An approximate optimum can be found by assuming that all of the Garabedian shape coefficients $|{\it\Delta}_{mn}/{\it\Delta}_{0,0}|$ that appear in the sum ${\it\Sigma}^{\prime }$ are small, then

(3.23) $$\begin{eqnarray}{\it\alpha}_{mn}\approx 2\frac{n-\bar{{\it\iota}}_{g}m}{m}\frac{{\it\Delta}_{m,n}}{{\it\Delta}_{00}},\end{eqnarray}$$

where $\bar{{\it\iota}}_{g}=\bar{{\it\iota}}({\it\psi}_{g})\approx \bar{{\it\iota}}({\it\psi}_{c})$ . Note that the toroidal flux dependence of the model Hamiltonian must be ${\it\psi}_{t}^{m/2}$ not ${\it\psi}_{t}^{|m/2|}$ to obtain a general solution for ${\it\alpha}_{mn}$ .

3.6.3. Removing islands

As discussed in § 3.1, it is important to be able to investigate the behaviour of non-resonant divertors in the limit as the last confining magnetic surface closer approximates a magnetic surface with sharp edges. To carry out this study, it is necessary to have a way to remove magnetic islands by adjusting the ${\it\alpha}_{mn}$ . A procedure for doing this can be based on a method developed by Greene (Reference Greene1979).

Greene’s method for eliminating islands in a canonical-coordinate region has three elements: (1) find the closed field line trajectories, called fixed points, of the Hamiltonian ${\it\psi}_{p}({\it\psi}_{t},{\it\theta},{\it\varphi})$ ; (2) construct the mapping matrix $\unicode[STIX]{x1D648}$ for trajectories near a fixed point; (3) adjust the appropriate coefficient ${\it\alpha}_{mn}$ in the Hamiltonian until the trace of $\unicode[STIX]{x1D648}$ equals 2.

(i) Find the fixed points. The trajectories of the model Hamiltonian, (3.1), are top–bottom symmetric in the primary cross-sections, ${\it\varphi}=0,{\rm\pi}$ . Letting $x\equiv \sqrt{{\it\psi}_{t}}\cos {\it\theta}$ be the horizontal and $y\equiv \sqrt{{\it\psi}_{t}}\sin {\it\theta}$ the vertical coordinate in a primary cross-section, the fixed points are located along the $x$ axis. To find fixed points, integrate from the starting point $(x=x_{0},y=0)$ at ${\it\varphi}=0$ to ${\it\varphi}=2{\rm\pi}M$ , where the number of toroidal transits $M$ is an integer, to obtain $x_{M}(x_{0})$ and $y_{M}(x_{0})$ . A necessary condition for an $M\text{th}$ -order fixed point is that $y_{M}(x_{0})$ change sign as $x_{0}$ is varied, with the fixed point located at the place where $y_{M}(x_{0})=0$ .

To obtain closure at an $M\text{th}$ -order fixed point, the field line trajectory advances from ${\it\theta}=0$ to ${\it\theta}=2{\rm\pi}N$ with $N$ an integer. The rotational transform of the magnetic field line that passes through the fixed point is ${\it\iota}_{f}=N/M$ .

(ii) Construct the near fixed-point map. Close to a fixed point of order $M$ , the trajectory map for going from ${\it\varphi}=0$ to ${\it\varphi}=2{\rm\pi}M$ has the form

(3.24) $$\begin{eqnarray}\left(\begin{array}{@{}c@{}}{\it\delta}x_{k+1}\\ {\it\delta}y_{k+1}\end{array}\right)=\unicode[STIX]{x1D648}\boldsymbol{\cdot }\left(\begin{array}{@{}c@{}}{\it\delta}x_{k}\\ {\it\delta}y_{k}\end{array}\right),\end{eqnarray}$$

where ${\it\delta}x_{k}$ is the $k\text{th}$ iterate of the horizontal distance, $x-x_{f}$ , of the field line from the fixed point, and ${\it\delta}y_{k}$ the vertical distance from the fixed point, and the integer $k$ is given by $kM$ being the number of toroidal transits.

The matrix is constructed by following two trajectories neighbouring the fixed point from ${\it\varphi}=0$ to ${\it\varphi}=2{\rm\pi}M$ . The first has the starting point $({\it\delta}x_{1},0)$ and ends at $({\it\delta}x_{11},{\it\delta}y_{12})$ , the second has the starting point $(0,{\it\delta}y_{2})$ and ends at $({\it\delta}x_{21},{\it\delta}y_{22})$ . Neighbouring means $|{\it\delta}x_{1}|\rightarrow 0$ and $|{\it\delta}y_{2}|\rightarrow 0$ . The mapping matrix is

(3.25) $$\begin{eqnarray}\unicode[STIX]{x1D648}=\left(\begin{array}{@{}c@{}}\displaystyle \frac{{\it\delta}x_{11}}{{\it\delta}x_{1}}\displaystyle \frac{{\it\delta}y_{12}}{{\it\delta}x_{1}}\\ \displaystyle \frac{{\it\delta}x_{21}}{{\it\delta}y_{2}}\displaystyle \frac{{\it\delta}y_{22}}{{\it\delta}y_{2}}\end{array}\right).\end{eqnarray}$$

The matrix $\unicode[STIX]{x1D648}$ should have a unit determinant, $D\equiv M_{11}M_{22}-M_{12}M_{21}=1$ , since the magnetic field line Hamiltonian is area preserving. The trace, $T\equiv M_{11}+M_{22}$ , will be the important quantity and will have an error ${\sim}|\sqrt{D}-1|$ .

(iii) Greene’s trace condition. To understand Greene’s trace condition that the trace must equal 2 to eliminate islands, first note that the trace $T$ and determinant $D$ give the eigenvalues of a $2\times 2$ matrix. That is, the matrix $\unicode[STIX]{x1D648}$ can be written in the form $\unicode[STIX]{x1D648}=\unicode[STIX]{x1D64B}\boldsymbol{\cdot }{\bf\mu}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}^{-1}$ , where ${\bf\mu}$ is diagonal with its two diagonal elements, ${\it\mu}_{\pm }=T/2\pm \sqrt{(T/2)^{2}-D}$ , the eigenvalues. For the mapping matrix, $D=1$ . As $T\rightarrow 2$ , the two eigenvalues are ${\it\mu}_{\pm }=1\pm \sqrt{T-2}$ .

When the near fixed point mapping is iterated $k$ times, the mapping matrix $\unicode[STIX]{x1D648}_{k}$ is a $k$ -fold product of the matrix $\unicode[STIX]{x1D648}$ , so $\unicode[STIX]{x1D648}_{k}=\unicode[STIX]{x1D64B}\boldsymbol{\cdot }{\bf\mu}_{k}\boldsymbol{\cdot }\unicode[STIX]{x1D64B}^{-1}$ where the diagonal elements of ${\bf\mu}_{k}$ are ${\it\mu}_{k\pm }=({\it\mu}_{\pm })^{k}=\text{e}^{\pm k\sqrt{T-2}}$ as $T\rightarrow 2$ . When $T>2$ , the trajectories exponentiate away from the fixed point, so the fixed point is the X-point of an island. When $T<2$ , the trajectories orbit the fixed point, so the fixed point is the O-point of an island. Only when $T=2$ can the fixed point lie on a magnetic surface.