2. Plasma equilibria and supporting coils

This paper builds upon the construction of the plasma equilibrium as a stationary point of the MRxMHD energy functional, which we now describe.Footnote ⁵

The computational domain, $\varOmega \subset \mathbb {R}^{3}$, is considered to be a solid torus with computational boundary ${{\mathcal {D}}} \equiv \partial \varOmega$. The latter is a prescribed toroidal surface, which unless explicitly stated otherwise is held fixed throughout the calculation. The computational domain contains the plasma volume, ${{\mathcal {V}}}\subset \varOmega$, a smaller solid torus enclosed by the plasma boundary ${{\mathcal {S}}}\equiv \partial \mathcal {V}$. The plasma volume is further partitioned into $v = 1, \dots , N_V$ nested toroidal subvolumes, ${{\mathcal {V}}}_v$, which are separated by a set of ideal interfaces, $\mathcal {I}_v$ so that $\partial \mathcal {V}_v={{\mathcal {I}}}_v\cup {{\mathcal {I}}}_{v-1}$, with the outermost interface coinciding with the plasma boundary, ${{\mathcal {I}}}_{N_V} = {{\mathcal {S}}}$. The innermost subregion is a simple (solid) torus, and each other subregion is a toroidal annulus which may be thought of as a ‘hollow’ torus. In each ${{\mathcal {V}}}_v$, the magnetic field is written as the curl of the magnetic vector potential, $\boldsymbol {B}_v = \boldsymbol {\nabla } \times {\boldsymbol {A}}_v$. In the innermost toroidal region, the toroidal flux, computed as a poloidal loop integral of ${\boldsymbol {A}}_1$ on ${{\mathcal {I}}}_1$, is constrained to match a prescribed value. In each annular region, both the enclosed toroidal and poloidal fluxes, computed, respectively, as the difference between the poloidal and toroidal loop integrals of ${\boldsymbol {A}}_v$ on ${{\mathcal {I}}}_v$ and ${{\mathcal {I}}}_{v-1}$, are constrained. In each region, the magnetic helicity, defined as the volume integral of ${\boldsymbol {A}}_v \cdot \boldsymbol {B}_v$, is also constrained. The magnetic field in each region is taken to be tangential to that region's boundary, $\boldsymbol {n}\cdot {\boldsymbol {B}}_v=0$ on ${{\mathcal {I}}}_v$ where $\boldsymbol {n}$ is a unit outward normal vector, but otherwise the topology of the field lines is unconstrained and magnetic islands and irregular field lines are allowed in each ${\mathcal {V}}_v$.

For free-boundary calculations, by which we mean that the plasma boundary is allowed to move, we include the vacuum region ${{\mathcal {V}}}_+=\varOmega {\setminus} {{\mathcal {V}}}$, with the inner boundary of ${{\mathcal {V}}}_+$ coinciding with ${{\mathcal {S}}}$ and its outer boundary with ${{\mathcal {D}}}$. It is on ${{\mathcal {D}}}$ that the normal plasma field plus the normal external field is required to equal the normal total field.

By construction, there are no currents inside ${{\mathcal {V}}}_+$. We may employ a scalar potential and write the vacuum field as $\boldsymbol {B}_{+} = \boldsymbol {\nabla } \varPhi$.Footnote ⁶ To obtain a unique solution, constraints are imposed on the linking and net toroidal plasma currents, which are computed as loop integrals of $\boldsymbol {B}_+$. We take $\boldsymbol {B}_+$ to be tangential to the plasma surface, $\boldsymbol {n} \boldsymbol {\cdot } \boldsymbol {\nabla } \varPhi |_{{{\mathcal {S}}}}=0$ where $\boldsymbol {n}$ is a unit outward normal vector on ${{\mathcal {S}}}$. The normal field on the computational boundary ${{\mathcal {D}}}$, $B_{T,n} \equiv \bar {\boldsymbol {n}} \boldsymbol {\cdot } \boldsymbol {\nabla } \varPhi |_{{{\mathcal {D}}}}$ where $\bar {\boldsymbol {n}}$ is normal to ${{\mathcal {D}}}$, is the sum, $B_{P,n} + B_{E,n}$, of the plasma field and the external field, neither of which are necessarily known a priori; generally, $B_{T,n} \ne 0$.

The equilibria that we consider are stationary points of the following energy functional:

(2.1)\begin{equation} {{\mathcal{F}}} \equiv \sum_{v=1}^{N_V} {{\mathcal{F}}}_v + {{\mathcal{F}}}_+, \end{equation}

where

(2.2)\begin{gather} {{\mathcal{F}}}_v \equiv \int_{{\mathcal{V}}_v} \left( \frac{p_v}{\gamma-1} + \frac{B_v^{2}}{2} \right) \mathrm{d} v - \frac{\mu_v}{2} \left[ \int_{{\mathcal{V}}_v} {\boldsymbol{A}}_v \cdot \boldsymbol{B}_v \, \mathrm{d} v - H_{v} \right], \end{gather}

(2.3)\begin{gather}{{\mathcal{F}}}_+{\equiv} \int_{{{\mathcal{V}}}_+} \tfrac{1}{2} \boldsymbol{\nabla} \varPhi \boldsymbol{\cdot} \boldsymbol{\nabla} \varPhi \, \mathrm{d} v -\int_{{{\mathcal{D}}}} \bar{\varPhi} \boldsymbol{B}_T \cdot \mathrm{d} \bar{\boldsymbol{s}}, \end{gather}

with $p_v$ the pressure in the $v$-th region,Footnote ⁷ $\gamma$ is the specific heat ratio (adiabatic index) and $H_v$ is the prescribed value of the magnetic helicity in ${{\mathcal {V}}}_v$. The constraints on the toroidal and poloidal fluxes in the plasma volumes, ${{\mathcal {V}}}_v$, are implied; that is, the fields $\boldsymbol {B}_v$ that we consider in (2.2) are restricted to the class of divergence-free vector fields that are tangential to the interfaces and have the prescribed values of the toroidal and poloidal fluxes. The helicity constraints are enforced explicitly using the Lagrange multipliers $\mu _v$. The loop integral constraints on $\boldsymbol {\nabla }\varPhi$ are also implied: the allowed $\varPhi$ in (2.3) is restricted to the class of multivalued scalar functions with prescribed loop integrals of $\boldsymbol {\nabla }\varPhi$. We write $\bar {\varPhi }\equiv \varPhi |_{{{\mathcal {D}}}}$ to denote the scalar potential evaluated on ${{\mathcal {D}}}$ to distinguish it from $\varPhi \equiv \varPhi |_{{{\mathcal {S}}}}$ on ${{\mathcal {S}}}$. Similarly, we write $\mathrm {d}\bar {\boldsymbol {s}}$ to denote an infinitesimal surface element on ${{\mathcal {D}}}$.

We may call ${{\mathcal {F}}}$ the free-plasma-boundary, fixed-computational-boundary, generalized- boundary-condition MRxMHD energy functional. This is to accommodate the common understanding in the magnetic confinement community that a free-boundary equilibrium calculation allows the plasma boundary to move, which it does in this case, and the common terminology in the broader mathematical community that refers to the ‘boundary’ as boundary of the computational domain, which is fixed in this case. We refer to the boundary conditions as ‘generalized’ because the total normal field, $B_{T,n}$, is not constrained to be zero. For brevity, we hereafter refer to ${{\mathcal {F}}}$ simply as the energy functional.

Using the virtual-casing principle (Shafranov & Zakharov Reference Shafranov and Zakharov1972), the normal plasma field, $B_{P,n}$, produced by currents within the plasma volume is computed on the computational boundary as

(2.4)\begin{equation} B_{P,n}(\bar{\boldsymbol{x}}) \equiv \left( - \frac{1}{4{\rm \pi}} \int_{{{\mathcal{S}}}} \frac{\left( \boldsymbol{B}_T|_{{{\mathcal{S}}}_+} \times \mathrm{d} \boldsymbol{s} \right) \times \boldsymbol{r}}{r^{3}} \right) \cdot \bar{\boldsymbol{n}}(\bar{\boldsymbol{x}}), \end{equation}

where ${{\boldsymbol {r}}} \equiv \boldsymbol {x} - \bar {\boldsymbol {x}}$ for $\boldsymbol {x} \in {{\mathcal {S}}}$ and $\bar {\boldsymbol {x}}\in {{\mathcal {D}}}$, and we use $\mathrm {d} \boldsymbol {s}$ to denote a surface element of ${{\mathcal {S}}}$. Here, $\boldsymbol {B}_T|_{{{\mathcal {S}}}_+}=\boldsymbol {\nabla } \varPhi |_{{{\mathcal {S}}}_+}$ is the total magnetic field immediately outside of ${{\mathcal {S}}}$. To accommodate for the possible existence of sheet currents lying on the plasma boundary, we must use in (2.4) the total magnetic field lying immediately outside the plasma boundary; that is, we must use the magnetic field on the inner boundary of the vacuum region. When that information is not available, we can use the magnetic field immediately inside the plasma boundary. In § 4, we shall exploit the circumstance that the evaluation of the plasma field on the computational boundary is a non-local, linear operator of the scalar potential, $\varPhi$.

The difference between the total normal magnetic field and the plasma normal field on ${{\mathcal {D}}}$, which we call the required external normal field, must be provided by an external source. A plasma cannot create its own magnetic bottle (Shafranov Reference Shafranov1966). This quantity, $D_n \equiv B_{T,n}-B_{P,n}$, is very closely related to $B_{E,n}$, the provided external field, but they may differ. The difference is what we later call the ‘coil error’.

For clarity of exposition regarding the externally applied field and to facilitate discussion of the various combined plasma–coil optimization algorithms considered in § 5, we imagine a typical and easily differentiable example; namely, that the external magnetic field is provided by a set of $i = 1, \ldots , N_C$ current-carrying one-dimensional curves (filaments) with geometry $\boldsymbol {c}_i$ that produce a magnetic field given by the Biot–Savart law,

(2.5)\begin{equation} \boldsymbol{B}_E(\bar{\boldsymbol{x}}) ={-}\sum_i \frac{\mu_0 I_i}{4{\rm \pi}}\int_{\boldsymbol{c}_i} \frac{\boldsymbol{r} \times \mathrm{d} {{\boldsymbol{l}}}}{r^{3}}, \end{equation}

where $I_i$ is the current in the $i$-th coil and $\mathrm {d} {{\boldsymbol {l}}}$ is the differential line segment. For brevity, we shall hereafter ignore the dependency on the magnitude of the currents and we set $I_i = 4{\rm \pi} /\mu _0$. It is a simple matter to extend the following to allow for the variation of the coil currents. It is also a simple matter to extend the following to the case where a surface potential on a prescribed winding surface provides the external field (Merkel Reference Merkel1987; Landreman Reference Landreman2017), or to the case where a ‘finite-build’ approximation is used for the coils (Li et al. Reference Li, Liu, Xu, Shimizu, Kinoshita, Okamura, Isobe, Xiong, Luo and Cheng2020; McGreivy et al. Reference McGreivy, Hudson and Zhu2021; Singh et al. Reference Singh, Kruger, Bader, Zhu, Hudson and Anderson2020), or when permanent magnets are included (Helander et al. Reference Helander, Drevlak, Zarnstorff and Cowley2020; Landreman & Zhu Reference Landreman and Zhu2020; Zhu et al. Reference Zhu, Hammond, Brown, Gates, Zarnstorff, Corrigan, Sibilia and Feibush2020).

If the required external field, $D_{n}$, is given, the geometry of the coils may be chosen to minimize a suitable ‘coil-penalty’ functional; for example (Merkel Reference Merkel1987; Landreman Reference Landreman2017; Hudson et al. Reference Hudson, Zhu, Pfefferlé and Gunderson2018; Zhu et al. Reference Zhu, Hudson, Song and Wan2018),

(2.6)\begin{equation} {{\mathcal{E}}} \equiv \varphi_2 + \omega L, \end{equation}

where $\varphi _2 \equiv \int _{{{\mathcal {D}}}} (Q^{2} \,\mathrm {d}\bar {s})/2$ is the commonly used quadratic-flux error functional, $Q = D_n - B_{E,n}$ where $B_{E,n}[\boldsymbol {c},{{\mathcal {D}}}]$ is the normal component of $\boldsymbol {B}_E$ on ${{\mathcal {D}}}$ produced by the filamentary coils, where we use $\boldsymbol {c}$ to denote the geometry of all the coils; i.e. $\boldsymbol {c} \equiv \{ \boldsymbol {c}_i, i = 1, \dots N_C \}$. We include a regularization term, $L = L[\boldsymbol {c}]$, which we take to be the total length of the coils, and $\omega$ is a scalar penalty.

The geometry of the optimal coils is given by $\delta {{\mathcal {E}}}/\delta \boldsymbol {c} = 0$. This equation may be solved using descent algorithms (Zhu et al. Reference Zhu, Hudson, Song and Wan2017; Hudson et al. Reference Hudson, Zhu, Pfefferlé and Gunderson2018) or Newton-style methods (Zhu et al. Reference Zhu, Hudson, Song and Wan2018). Practically, it is not possible to obtain a coil set that exactly produces the required normal field and generally $\varphi _2\ne 0$. This is why the required external normal field, $D_n$, must be distinguished from the actual external normal field, $B_{E,n}$. Upon solving $\delta {{\mathcal {E}}}/ \delta \boldsymbol {c} = 0$, $\varphi _2$ measures what we call the coil error.Footnote ⁸

In the above, we have described the three basic components of the combined plasma–coil optimization algorithms that we consider in § 5; namely, the energy functional, ${{\mathcal {F}}}$, the Biot–Savart law, the coil-penalty functional, ${{\mathcal {E}}}$, and the evaluation of the normal plasma field, $B_{P,n}$, on ${{\mathcal {D}}}$ using the virtual-casing functional. In the following section, § 3, we present the first and second variations of ${{\mathcal {F}}}$, ${{\mathcal {E}}}$ and $B_{P,n}$ that will be required in § 5. In § 4, we elaborate upon the numerical construction of the magnetic fields and show that, instead of prescribing the total normal magnetic field, $B_{T,n}$, on ${{\mathcal {D}}}$, computing the equilibrium and then determining the coil geometry that provides (approximately) the required external normal field, we may directly consider the case where the external normal field, $B_{E,n}$, on ${{\mathcal {D}}}$ is given. We can solve for the virtual-casing self-consistent vacuum magnetic field using a Galerkin method. In § 5, we show that the derivatives described in § 3 and the Galerkin construction of the magnetic fields may be combined to construct a variety of combined plasma–coil stellarator optimization algorithms.

3. Variation of the energy, coil-penalty and virtual-casing functionals

In this section, we present the variations of the energy, coil-penalty and virtual casing functionals that are used in § 5. The variational calculus provides useful insight into the physics and mathematics of the different optimization problems, as we shall comment upon in § 6.

Requesting that the first variations of the energy functionals vanish results in weak formulations of the various coupled boundary value problems. When the boundary data is ‘smooth’ enough, it is expected that weak solutions coincide with (strong) solutions. In what follows, we assume the fields are sufficiently regular to apply integration by parts and derive local Euler–Lagrange equations. Conversely, we convert any linear elliptic partial differential equations subject to boundary conditions into variational problems (weak formulation) simply by integrating against arbitrary variations. In the case where the bilinear functional thus obtained is symmetric (self-adjoint), the variational problem can be associated with an energy minimization principle. However, as we will see in § 3.4, the weak formulation of certain boundary value problems do not necessarily result in self-adjoint bilinear functionals, in which case there is no immediate energy minimization principle.

The weak formulation is the basis of our numerical representation of solutions, in particular leading to a spectral-Galerkin method where our fields are approximated by a finite linear combination of Fourier modes and orthogonal polynomials. The coefficients of the linear combination become our degrees of freedom and the weak formulation boils down to a matrix inversion.

3.1. Energy functional

In the plasma volumes we write $\boldsymbol {B}_v = \boldsymbol {\nabla }\times \boldsymbol {A}_v$, and thus $\delta \boldsymbol {B}_v = \boldsymbol {\nabla }\times \delta \boldsymbol {A}_v$. We use the fact that the fields in the vicinity of the interfaces obey ideal MHD, so that the variation of the magnetic field on the interfaces can be expressed with an displacement vector $\boldsymbol {\xi }_v$, $\delta \boldsymbol {B}_v = \boldsymbol {\nabla }\times (\boldsymbol {\xi }_v\times \boldsymbol {B}_v)$, i.e. $\delta \boldsymbol {A}_v = \boldsymbol {\xi }_v\times \boldsymbol {B}_v + \nabla g_v$, where $g_v$ is a gauge potential. In MRxMHD, the mass is constrained in each volume $\mathcal {V}_v$ by the isentopic ideal-gas constraint, $p_v = a_v / V_v^{\gamma }$, where $V_v$ is the volume of the $v$-th region, $V_v = \int _{\mathcal {V}_v} \,\mathrm {d} v$ and $a_v$ is a constant. The first variation with respect to deformations of the volume boundary is given by $\delta p_v = - \gamma p_v (\int _{{{{\mathcal {I}}}}_v} \boldsymbol {\xi }_v\cdot \mathrm {d} {\boldsymbol {s}}-\int _{{{{\mathcal {I}}}}_{v-1}} \boldsymbol {\xi }_{v-1}\cdot \mathrm {d} {\boldsymbol {s}})/V_v$.

Using $\boldsymbol {B}_v \cdot {\boldsymbol {n}}_v = 0$ on the interfaces, where $\boldsymbol {n}_v$ is the normal vector of the $v$-th interface, one can derive

(3.1)\begin{equation} \delta {\mathcal{F}} = \sum_{v=1}^{N_v}\left(\int_{{\mathcal{V}}_v} \frac{\delta{\mathcal{F}}_v}{\delta \boldsymbol{A}_v} \cdot \delta \boldsymbol{A}_v \,\mathrm{d} v + \int_{{{\mathcal{I}}}_v} \frac{\delta {\mathcal{F}}_v}{\delta \boldsymbol{\xi}_{v}} \cdot \boldsymbol{\xi}_{v} \,\mathrm{d} s\right) + \delta \mathcal{F}_+, \end{equation}

where

(3.2)\begin{gather} \frac{\delta{\mathcal{F}}_v}{\delta \boldsymbol{A}_v} = \boldsymbol{\nabla}\times\boldsymbol{B}_v - \mu_v \boldsymbol{B}_v, \end{gather}

(3.3)\begin{gather}\frac{\delta{\mathcal{F}}_v}{\delta \boldsymbol{\xi}_{v}} ={-} \left[ \! \left[ p_v + \frac{B_v^{2}}{2} \right] \! \right] \boldsymbol{n}_v. \end{gather}

In vacuum, the first variation in ${{\mathcal {F}}}_+$ with respect to variations in $\varPhi$ is

(3.4)\begin{equation} \delta {{\mathcal{F}}}_+{=} \int_{\mathcal{V}} \frac{\delta {{\mathcal{F}}}_+}{\delta \varPhi} \delta \varPhi \, \mathrm{d} v + \int_{{{\mathcal{D}}}}\delta \bar\varPhi \left( \boldsymbol{\nabla} \bar{\varPhi} - \boldsymbol{B}_T \right) \cdot \mathrm{d}\bar{{\boldsymbol{s}}} - \int_{{{\mathcal{S}}}} \delta \varPhi \,\boldsymbol{\nabla} \varPhi \cdot \mathrm{d} {\boldsymbol{s}.} \end{equation}

Setting the functional derivative,

(3.5)\begin{equation} \frac{\delta {{\mathcal{F}}}_+}{\delta \varPhi} ={-}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\nabla} \varPhi, \end{equation}

to zero leads to the Laplace equation. The other terms provide corresponding Neumann boundary conditions.

The second variation with respect to deformations of the interface geometry of ${\mathcal {F}}_v$ is given by

(3.6)\begin{align} \delta^{2} {\mathcal{F}}_v & = \int_{{\mathcal{V}}_v} \delta \boldsymbol{A}_v \boldsymbol{\cdot} ( \boldsymbol{\nabla}\times\boldsymbol{\nabla}\times\delta\,\,\tilde{\!\!\boldsymbol{A}}_v-\mu \boldsymbol{\nabla}\times \delta\,\,\tilde{\!\!\boldsymbol{A}}_v )\,\mathrm{d} v + \gamma p_v \frac{\displaystyle\int_{{{\mathcal{I}}}_v} \boldsymbol{\xi}_v\cdot \mathrm{d} \boldsymbol{s} \int_{{{\mathcal{I}}}_v} \tilde{\boldsymbol{\xi}}_v\cdot \mathrm{d} {\boldsymbol{s}}}{V_v}\nonumber\\ &\quad - \int_{{{\mathcal{I}}}_v} ( \boldsymbol{B}_v\boldsymbol{\cdot}\boldsymbol{\nabla}\times \delta\,\,\tilde{\!\!\boldsymbol{A}})\boldsymbol{\xi}_v \cdot\mathrm{d}{\boldsymbol{s}} - \int_{{{\mathcal{I}}}_v} \delta \boldsymbol{A}_v \boldsymbol{\cdot}(\mu_v \boldsymbol{B}_v - \boldsymbol{\nabla} \times \boldsymbol{B}_v)\tilde{\boldsymbol{\xi}}_v \cdot\mathrm{d}{\boldsymbol{s}} \nonumber\\ &\quad - \int_{{{\mathcal{I}}}_v} \tilde{\boldsymbol{\xi}}_v\boldsymbol{\cdot}\boldsymbol{n}_v \boldsymbol{\nabla}\boldsymbol{\cdot}\left(\left(p_v+\frac{B_v^{2}}{2}\right)\boldsymbol{\xi}_v\right) \mathrm{d} s, \end{align}

where the tilde ($\sim$) indicates that the second variation can differ from the first. Here, for simplicity, we allowed only one of the two interfaces to vary. For a complete expression (see (3.8)), one has to be careful since cross-terms, including integrals over two different interfaces, appear in the $\gamma p_v$ term.

In vacuum, the second variation with respect to the scalar potential is given by

(3.7)\begin{equation} \delta^{2}{{\mathcal{F}}}_+{=} -\int_\mathcal{V} \delta \varPhi \, \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{\nabla} \delta \tilde{\varPhi} \, \mathrm{d} v - \int_{{{\mathcal{D}}}} \delta \bar{\varPhi}\, \bar{\boldsymbol{n}}\boldsymbol{\cdot}\boldsymbol{\nabla} \delta \tilde{\bar{\varPhi}} \, \mathrm{d} \bar{s}. \end{equation}

Collecting all the different terms for the second variation, we obtain

(3.8)\begin{align} \delta^{2} {\mathcal{F}} &= \sum_v\left(\int_{{\mathcal{V}}_v} \delta \boldsymbol{A}_v \boldsymbol{\cdot} ( \boldsymbol{\nabla}\times\boldsymbol{\nabla}\times\delta\,\,\tilde{\!\!\boldsymbol{A}}_v-\mu_v \boldsymbol{\nabla}\times \delta\,\,\tilde{\!\!\boldsymbol{A}}_v )\,\mathrm{d} v \right. \nonumber\\ &\quad +\gamma\left( \frac{p_v}{V_v}+\frac{p_{v+1}}{V_{v+1}}\right) \left(\int_{{{\mathcal{I}}}_v} \boldsymbol{\xi}_v\cdot \mathrm{d} \boldsymbol{s}\right) \left(\int_{{{\mathcal{I}}}_v} \tilde{\boldsymbol{\xi}}_v\cdot \mathrm{d} \boldsymbol{s} \right) \nonumber\\ &\quad -\gamma \frac{p_v}{V_v} \left[\left(\int_{{{\mathcal{I}}}_{v-1}}\tilde{\boldsymbol{\xi}}_{v-1} \cdot \mathrm{d} \boldsymbol{s} \right) \left(\int_{{{\mathcal{I}}}_v}\boldsymbol{\xi}_v\cdot \mathrm{d} \boldsymbol{s}\right) + \left(\int_{{{\mathcal{I}}}_{v-1}}\boldsymbol{\xi}_{v-1}\cdot \mathrm{d} \boldsymbol{s} \right) \left(\int_{{{\mathcal{I}}}_v}\tilde{\boldsymbol{\xi}}_v\cdot \mathrm{d} \boldsymbol{s}\right)\right]\nonumber\\ &\quad \left.+ \int_{{{\mathcal{I}}}_v} \left[ \! \left[ \boldsymbol{B} \boldsymbol{\cdot} ((\boldsymbol{B}\times\boldsymbol{\nabla})\times\boldsymbol{n}) \right] \! \right] \tilde{\boldsymbol{\xi}}_v \boldsymbol{\cdot} \boldsymbol{n} \,\boldsymbol{\xi}_v \cdot\mathrm{d}\boldsymbol{s} \right) + \delta^{2} {{\mathcal{F}}}_+, \end{align}

where we used the equilibrium expressions; we included variations of both interfaces; we used $\boldsymbol {B} \boldsymbol {\cdot } ((\boldsymbol {B}\times \boldsymbol {\nabla })\times \boldsymbol {n})=\boldsymbol {B}^{2} \boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {n} - \boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla } \boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {B}$; and we set $\boldsymbol {\xi }_v$ and $\tilde {\boldsymbol {\xi }}_v$ proportional to the normal vector $\boldsymbol {n}$ since only normal variations of the interfaces affect the plasma energy. The second variation of the energy functional determines the MRxMHD stability of the equilibrium state (Kumar et al. Reference Kumar, Qu, Hole, Wright, Loizu, Hudson, Baillod, Dewar and Ferraro2021).Footnote ⁹

3.2. Coil-penalty functional

The coil-penalty functional, ${{\mathcal {E}}}$ defined in (2.6) is considered to be a functional of the coil geometry, $\boldsymbol {c}$, the required normal magnetic field, $D_n$, on the computational boundary, ${{\mathcal {D}}}$, and on ${{\mathcal {D}}}$ itself, which is, however, kept fixed here; i.e. ${{\mathcal {E}}}={{\mathcal {E}}}[\boldsymbol {c},D_n]$. We may formally write the first variation, $\delta {{\mathcal {E}}}$, with respect to variations $\delta \boldsymbol {c}$, and $\delta D_n$ as

(3.9)\begin{equation} \delta {{\mathcal{E}}} = \sum_{i} \int_0^{2{\rm \pi}} \frac{\delta {{\mathcal{E}}}}{\delta \boldsymbol{c}_i} \cdot \delta \boldsymbol{c}_i \, \mathrm{d} \alpha + \int_{{{\mathcal{D}}}} \frac{\delta {{\mathcal{E}}}}{\delta D_n} \cdot \delta D_n \, \mathrm{d} \bar{s}, \end{equation}

where $\alpha$ is an angle-like curve parameter, $\boldsymbol {c}_i(\alpha +2{\rm \pi} )=\boldsymbol {c}_i(\alpha )$. The functional derivatives are

(3.10)\begin{gather} \frac{\delta {{\mathcal{E}}}}{\delta \boldsymbol{c}_i} = \boldsymbol{c}_i' \times \left( \int_{{{\mathcal{D}}}} {{\boldsymbol{R}}}_{i,n} Q \,\mathrm{d} s + \omega \, \boldsymbol{\kappa}_i \right), \end{gather}

(3.11)\begin{gather}\frac{\delta {{\mathcal{E}}}}{\delta D_n} = Q, \end{gather}

where $\boldsymbol {c}_i^{\prime }$ is the derivative of the $i$-th coil with respect to $\alpha$. We have written ${\textsf {R}}_i \equiv 3 \, \boldsymbol {r}_i \boldsymbol {r}_i / r_i^{5} - {{\sf I}} / r_i^{3}$ and defined ${{\boldsymbol {R}}}_{i,n} = {\textsf {R}}_i {\cdot } \boldsymbol {n}$, where $\boldsymbol {r}_i$ is the displacement between a point on the $i$-th coil and the evaluation point on ${{\mathcal {D}}}$, and ${{\sf I}}$ is the identity tensor or idemfactor. The curvature of the $i$-th coil is $\boldsymbol {\kappa }_i \equiv \boldsymbol {c}_i' \times \boldsymbol {c}_i'' / \boldsymbol {c}_i'^{3}$ and the functional derivatives are generalized expressions of the ones presented in Dewar, Hudson & Price (Reference Dewar, Hudson and Price1994) and Hudson et al. (Reference Hudson, Zhu, Pfefferlé and Gunderson2018). We used the expression for the variation of the normal component of the magnetic field with respect to the coil geometry

(3.12)\begin{equation} \delta B_{n} = \oint_i \frac{\delta B_n}{\delta \boldsymbol{c}_i} \cdot \delta\boldsymbol{c}_i \,\mathrm{d}{\alpha_i}, \end{equation}

with the functional derivatives

(3.13)\begin{equation} \frac{\delta B_n}{\delta \boldsymbol{c}_i} =\boldsymbol{c}_i' \times {\boldsymbol{R}}_{i,n}. \end{equation}

The required second variations of ${{\mathcal {E}}}$ are (Hudson et al. Reference Hudson, Zhu, Pfefferlé and Gunderson2018)

(3.14)\begin{gather} \delta^{2} {{\mathcal{E}}} = \sum_{i,j} \oint_i \oint_j \delta \boldsymbol{c}_i \cdot \frac{\delta^{2} {{\mathcal{E}}}}{\delta \boldsymbol{c}_i \delta\boldsymbol{c}_j} \cdot \delta \boldsymbol{c}_j \,\mathrm{d}{\alpha_i} \,\mathrm{d}{\alpha_j}, \end{gather}

(3.15)\begin{gather}\delta^{2} {{\mathcal{E}}} = \sum_i \oint_i \,\mathrm{d}{\alpha} \oint_{{{\mathcal{D}}}} \,\mathrm{d} \bar{s} \, \delta \boldsymbol{c}_i \cdot \frac{\delta^{2} {{\mathcal{E}}}}{\delta \boldsymbol{c}_i \delta D_n} \, \delta D_n, \end{gather}

where

(3.16)\begin{gather} \frac{\delta^{2} {{\mathcal{E}}}}{\delta \boldsymbol{c}_i\delta\boldsymbol{c}_j} = \oint_{{{\mathcal{D}}}} (\boldsymbol{c}_i' \times {\boldsymbol{R}}_{i,n}) (\boldsymbol{c}_j' \times {\boldsymbol{R}}_{j,n}) \,\mathrm{d} \bar s, \quad (\mbox{for}\ i \ne j), \end{gather}

(3.17)\begin{gather}\frac{\delta^{2} {{\mathcal{E}}}}{\delta\boldsymbol{c}_i \delta D_n } = \boldsymbol{c}_i' \times {{\boldsymbol{R}}}_{i,n}. \end{gather}

For the case where $i=j$ the variation cannot be written in such a compact form but also does not provide much insight. The other second variations of ${{\mathcal {E}}}$ are not required in § 5 and are omitted for brevity.

We only need the variation with respect to the geometry of the computational boundary, ${{\mathcal {D}}}$, when $D_n=-B_{P,n}=-\boldsymbol {B}_{P}\cdot \bar {\boldsymbol {n}}$. Since there is an explicit dependence with respect to $\bar{\boldsymbol {n}}$, we calculate the combined variation of ${{\mathcal {E}}}$ and $B_{P,n}$ with respect to variations in ${{\mathcal {D}}}$, which allows us to express the results in compact form. The first variation is

(3.18)\begin{equation} \delta {{\mathcal{E}}} = \int_{{{\mathcal{D}}}} \frac{\delta {{\mathcal{E}}}}{\delta {{\mathcal{D}}}} \cdot \delta {{\mathcal{D}}} \, \mathrm{d} \bar{s}, \end{equation}

where

(3.19)\begin{equation} \frac{\delta {{\mathcal{E}}}}{\delta {{\mathcal{D}}}} = \left( -\frac{1}{2}Q^{2} \boldsymbol{\nabla}\boldsymbol{ \cdot} \bar{\boldsymbol{n}} - \left( \boldsymbol{B}_{P,s} + \boldsymbol{B}_{E,s} \right) \boldsymbol{\cdot} \boldsymbol{\nabla} Q \right) \bar{\boldsymbol{n}}, \end{equation}

where we used (12) from Dewar et al. (Reference Dewar, Hudson and Price1994). The notation $\boldsymbol {f}_{s} = ( {{\sf I}} - \bar {\boldsymbol {n}}\bar {\boldsymbol {n}})\cdot \boldsymbol {f}$ denotes the surface projection of an arbitrary vector or tensor, $\boldsymbol {f}$, onto a surface, which in this case is the computational boundary. The second variation of ${{\mathcal {E}}}$ with respect to ${{\mathcal {D}}}$ and $\boldsymbol {c}_i$ is

(3.20)\begin{equation} \delta^{2} {{\mathcal{E}}} = \sum_i \oint_{\boldsymbol{c}_i} \delta \boldsymbol{c}_i \cdot \oint_{{{\mathcal{D}}}} \frac{\delta^{2} {{\mathcal{E}}}}{\delta \boldsymbol{c}_i \delta {{\mathcal{D}}}} \cdot \delta {{\mathcal{D}}} \, \mathrm{d} \bar{s} \, \mathrm{d}\alpha, \end{equation}

where

(3.21)\begin{equation} \frac{\delta^{2} {{\mathcal{E}}}}{\delta \boldsymbol{c}_i \delta {{\mathcal{D}}}} = \boldsymbol{c}_i^{\prime} \times \left( Q {\textsf{R}}_i \boldsymbol{\cdot} \boldsymbol{H} - {\textsf {R}}_{i,s} \boldsymbol{\cdot} \boldsymbol{\nabla} {\mathcal{H}} + (\boldsymbol{B}_{P,s}+\boldsymbol{B}_{E,s}) \boldsymbol{\cdot} \boldsymbol{\nabla} {{\boldsymbol{R}}}_{i,n} \right) \bar{\boldsymbol{n}}, \end{equation}

where the mean curvature of ${{\mathcal {D}}}$ is $\boldsymbol {H} \equiv - \bar {\boldsymbol {n}} \, (\boldsymbol {\nabla }\boldsymbol {\cdot }\bar {\boldsymbol {n}} )$. Equation (3.21) is a more general form of (12) of Hudson et al. (Reference Hudson, Zhu, Pfefferlé and Gunderson2018) where $D_n=0$.

3.3. Plasma normal field

The normal plasma field, $B_{P,n}$, on ${{\mathcal {D}}}$, as defined in (2.4), is considered to be a functional of the plasma boundary, ${\mathcal {S}}$, the total (tangential) magnetic field, $\boldsymbol {B}_T|_{{{\mathcal {S}}}}$ on ${{\mathcal {S}}}$, and on ${{\mathcal {D}}}$, which is kept fixed; i.e. $B_{P,n}=B_{P,n}[{\mathcal {S}},\boldsymbol {B}_T|_{{{\mathcal {S}}}}]$. The normal component of the magnetic field produced by plasma currents, $B_{P,n}$, given in (2.4) may be written as

(3.22)\begin{equation} B_{P,n} ={-} \frac{1}{4{\rm \pi}} \int_{{{\mathcal{S}}}} \boldsymbol{B}_{T}|_{{{\mathcal{S}}}} \cdot {\textsf {M}} \cdot \mathrm{d}\boldsymbol{s}, \end{equation}

where the antisymmetric tensor, ${\textsf {M}} \equiv ( {\boldsymbol {r}} \, \bar {\boldsymbol {n}} - \bar {\boldsymbol {n}} \, {\boldsymbol {r}} ) / r^{3}$, is defined, where ${\boldsymbol {r}} = \boldsymbol {x}-\bar {\boldsymbol {x}}$ for $\boldsymbol {x} \in {{\mathcal {S}}}$ and $\bar {\boldsymbol {x}}\in {{\mathcal {D}}}$. The first variation setting $\delta {{\mathcal {D}}} =0$ is

(3.23)\begin{equation} \delta B_{P,n} = \int_{{{\mathcal{S}}}} \frac{\delta B_{P,n}}{\delta {{\mathcal{S}}}} \cdot \delta {{\mathcal{S}}} \, \mathrm{d} s + \int_{{{\mathcal{S}}}} \frac{\delta B_{P,n}}{\delta \boldsymbol{B}_T|_{{{\mathcal{S}}}} } \cdot \delta \boldsymbol{B}_T|_{{{\mathcal{S}}}} \, \mathrm{d} s, \end{equation}

where

(3.24)\begin{gather} \frac{\delta B_{P,n}}{\delta {{\mathcal{S}}}} ={-} \frac{1}{4{\rm \pi}} \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{B}_T|_{{{\mathcal{S}}}} \cdot {\textsf {M}} ), \end{gather}

(3.25)\begin{gather}\frac{\delta B_{P,n}}{\delta \boldsymbol{B}_T|_{{{\mathcal{S}}}}} ={-} \frac{1}{4{\rm \pi}} {\textsf {M}} \cdot \boldsymbol{n}, \end{gather}

where we used again (12) from Dewar et al. (Reference Dewar, Hudson and Price1994) for the first expression. Note that only ${\textsf {M}}$ and ${\textsf {R}}$ depend on $\bar {\boldsymbol {x}}$. The relevant variation with respect to the computational boundary ${{\mathcal {D}}}$ is presented in the previous section. The second variations are not required in the following.

3.4. Supplied external field

To enable efficient free-boundary optimization algorithms, which we describe in § 5, we revisit the energy functional; in particular, we pay attention to the boundary condition for the normal field on ${{\mathcal {D}}}$.

It is perfectly legitimate to prescribe the total normal field, $B_{T,n}$, on ${{\mathcal {D}}}$. Indeed, every fixed-boundary equilibrium calculation does as much. As described in § 2, the total normal field is the sum of the plasma normal field and the external normal field, $B_{T,n} = B_{P,n} + B_{E,n}$, neither of which are known a priori in fixed-boundary calculation.

It seems unlikely that there are any circumstances for which we will a priori know $B_{P,n}$, except for the trivial equilibrium, the vacuum, for which $B_{P,n}=0$. After all, this is why the equilibrium calculation is required, to determine the plasma currents and the magnetic field that they produce.

In contrast, a particularly important calculation arises when $B_{E,n}$ is known. In this case, if we were to proceed with the approach of specifying $B_{T,n}$ on ${{\mathcal {D}}}$, then some type of ‘self-consistent’ iteration, for example, must be implemented to determine the $B_{T,n}$ that satisfies the matching condition, namely that $B_{T,n} - B_{P,n}[\boldsymbol {B}_T|_\mathcal {S}] = B_{E,n}$, where $B_{P,n}[\boldsymbol {B}_T|_{{{\mathcal {S}}}}]$ may be considered to be a linear, non-local operator acting on the tangential total field on the plasma boundary $\boldsymbol {B}_T|_{{{\mathcal {S}}}}$, which is only known after the equilibrium has been computed. In the current implementation of free-boundary SPEC (Hudson et al. Reference Hudson, Loizu, Zhu, Qu, Nührenberg, Lazerson, Smiet and Hole2020), for example, a Picard iteration over $B_{T,n}$ is required.

Here, we present a more direct strategy. The boundary value problem in the vacuum region is to find $\boldsymbol {B}_+ = \boldsymbol {\nabla }\varPhi$ such that $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {B}_+ = \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {\nabla } \varPhi = 0$ on ${{\mathcal {V}}}_+$ satisfying boundary conditions $\boldsymbol {n}\boldsymbol {\cdot }\boldsymbol {\nabla }\varPhi =0$ on ${{\mathcal {S}}}$ and $\bar{\boldsymbol {n}}\boldsymbol {\cdot } \boldsymbol {\nabla }\bar {\varPhi } - B_{P,n}[\boldsymbol {\nabla }\varPhi |_{{{\mathcal {S}}}}]= B_{E,n}$ on ${{\mathcal {D}}}$. This translates into the weak problem of finding $\varPhi$ such that

(3.26)\begin{equation} \int_{{{\mathcal{V}}}_+} \boldsymbol{\nabla}\psi \boldsymbol{\cdot} \boldsymbol{\nabla}\varPhi \,\mathrm{d} v + \int_{{{\mathcal{D}}}} \int_{{{\mathcal{S}}}} \frac{(\mathrm{d} \bar{\boldsymbol{s}}\times \boldsymbol{\nabla}\bar{\psi}) \boldsymbol{\cdot}(\mathrm{d} \boldsymbol{s}\times \boldsymbol{\nabla}\varPhi)}{|\bar{\boldsymbol{x}}-\boldsymbol{x}|} = \int_{{{\mathcal{D}}}} \bar{\psi} \, B_{E,n} \,\mathrm{d} \bar{s}, \quad \forall \psi, \end{equation}

where $\psi$ is an arbitrary (test) function. An important solvability condition is that the normal field of the coils is consistent with a divergence-free field, $\int _{{{\mathcal {D}}}} B_{E,n} \,\mathrm {d} \bar {s} =0$.

As it stands, it is not possible to cast this weak problem into an energy minimization problem because the second term on the left-hand side of (3.26) spoils the required symmetry of the functional. The first variation of a quadratic energy functional necessarily results in symmetric bilinear forms.

Energy principles do exist for exterior Neumann problems in the form of boundary integral methods (Giroire & Nédélec Reference Giroire and Nédélec1978). Their formulation, however, requires the plasma boundary and the computational boundary to coincide. The NESTOR code is a good example (Merkel Reference Merkel1986), as well as BIEST (Malhotra et al. Reference Malhotra, Cerfon, Imbert-Gérard and O'Neil2019). There are several disadvantages and challenges in implementing these energy principles. First, as seen in (3.26) in the limit ${{\mathcal {D}}}\to {{\mathcal {S}}}$, the integral operators tend to have singular kernels, which requires delicate numerical treatment. Second, in free-boundary calculations, the plasma boundary is varied to achieve force balance. This requires updating the numerical representation of the integral operators as well as the source term on the right-hand side of (3.26) every time the boundary changes. For this, the external field (from coils) is effectively needed in an entire volume not just on a single surface, which comes with a sizeable computational cost.Footnote ¹⁰

In the following section, § 4, we present a free-boundary equilibrium approach for a given external field that is based on a weak (Galerkin) method for constructing the required magnetic fields while keeping the computational boundary fixed and separate from the plasma boundary so that only the normal component of the external field on the computational boundary need be evaluated and only once. The virtual-casing self-consistent vacuum field is obtained as the solution to a set of linear equations, given below in (4.4). The matching condition is enforced directly in the computation of the vacuum field, and this removes the self-consistent iteration mentioned above that is otherwise required to match the equilibrium to the provided external field.

4. Galerkin method for constructing the vacuum field

Galerkin methods for constructing weak solutions to partial differential equations are well known. In the context of MRxMHD, a Galerkin method for constructing the magnetic fields for fixed- and free-boundary MRxMHD equilibria is described in Hudson et al. (Reference Hudson, Dewar, Dennis, Hole, McGann, von Nessi and Lazerson2012, Reference Hudson, Loizu, Zhu, Qu, Nührenberg, Lazerson, Smiet and Hole2020) and Qu et al. (Reference Qu, Pfefferlé, Hudson, Baillod, Kumar, Dewar and Hole2020). In this paper, we restrict our attention to the problem of constructing $\boldsymbol {B}_+$.

A standard numerical approach for finding stationary points of ${{\mathcal {F}}}_+$ is to represent the scalar potential using a mixed Chebyshev–Fourier representation; e.g. $\varPhi (s,\theta ,\zeta ) = I \theta + G \zeta + \sum _{l,m,n} \varPhi _{l,m,n} T_l(s)\exp (\textrm {i}m\theta -\textrm {i}n\zeta )$, where $(s,\theta ,\zeta )$ is a suitable toroidal coordinate system and $T_l(s)$ is the $l$-th Chebyshev polynomial, and $I$ and $G$ are given by the prescribed value of the loop $\text {integrals}/2{\rm \pi}$. When this representation is inserted into (2.2), and assuming that ${{\mathcal {S}}}$ and ${{\mathcal {D}}}$ do not change, the problem of calculating the vacuum field amounts to solving

(4.1)\begin{equation} {{\mathcal{A}}} \cdot \boldsymbol{\phi} = \boldsymbol{B}_T, \end{equation}

where $\boldsymbol {\phi }$ represents the vector of independent degrees-of-freedom, namely the $\varPhi _{l,m,n}$ and the Lagrange multipliers that may be used to enforce the constraints, and ${{\mathcal {A}}}$ is a symmetric matrix whose elements are the second derivatives of ${{\mathcal {F}}}_{+}$ with respect to these degrees-of-freedom and involve volume integrals of coordinate metric elements and the Chebyshev–Fourier functions. The matrix ${{\mathcal {A}}}$ depends on the geometry of ${{\mathcal {S}}}$ and ${{\mathcal {D}}}$; i.e. ${{\mathcal {A}}}={{\mathcal {A}}}[{{\mathcal {S}}},{{\mathcal {D}}}]$. The right-hand side vector, $\boldsymbol {B}_T$, contains the Fourier harmonics of $B_{T,n}$. The system of linear equations given in (4.1) is essentially a discrete realization of Laplace's equation, $\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {\nabla } \varPhi = 0$, with the supplied boundary conditions and loop integrals, written in matrix form. We can determine how the vacuum magnetic field varies with ${{\mathcal {S}}}$, ${\mathcal {D}}$ and $B_{T,n}$ using matrix perturbation methods,

(4.2)\begin{equation} {{\mathcal{A}}} \cdot \delta \boldsymbol{\phi} ={-} \delta {{\mathcal{A}}} \cdot \boldsymbol{\phi} + \delta \boldsymbol{B}_T, \end{equation}

where $\delta {\mathcal {A}} = \partial {{\mathcal {A}}} / \partial {{\mathcal {S}}} \cdot \delta {{\mathcal {S}}} + \partial {\mathcal {A}}/ \partial {{\mathcal {D}}} \cdot \delta {{\mathcal {D}}}$.

From (2.4), we recognize that the magnetic field produced by the plasma currents is a linear functional of the total (tangential) magnetic field on the immediate outside of the plasma boundary, namely $\boldsymbol {B}_+|_{{{\mathcal {S}}}_+}$. The immediate outside of ${{\mathcal {S}}}$ is the inner boundary of the vacuum region, ${{\mathcal {V}}}_+$, and the magnetic field in ${{\mathcal {V}}}_+$, namely $\boldsymbol {B}_+$, is the gradient of the scalar potential. The Fourier harmonics, $B_{P,n}$, of the plasma field on ${{\mathcal {D}}}$ are linear in the degrees-of-freedom of the scalar potential; that is, we may write

(4.3)\begin{equation} \boldsymbol{B}_P = {\mathcal{B}} \cdot \boldsymbol{\phi}, \end{equation}

where $\boldsymbol {B}_P$ represents the vector of Fourier harmonics of $B_{P,n}$, and ${\mathcal {B}}$ is a matrix derived from inserting the mixed Chebyshev–Fourier representation for $\varPhi$ into (2.4). We note that ${\mathcal {B}}$ depends only on the geometry of ${{\mathcal {S}}}$ and ${{\mathcal {D}}}$; i.e. ${\mathcal {B}}={\mathcal {B}}[{{\mathcal {S}}},{{\mathcal {D}}}]$.

Provided that the external field $\boldsymbol {B}_E$ is ‘tame’ (in particular $\int _{{{\mathcal {D}}}} B_{E,n}\,\mathrm {d} \bar {s}=0$), we put forward the following propositions: (i) for toroidal annular domains of practical interest, with the constraint that $\boldsymbol {B}_T \cdot \boldsymbol {n} = 0$ on the inner boundary, ${{\mathcal {S}}}$, there exists a normal magnetic field on the outer boundary, ${{\mathcal {D}}}$, that is the sum of an a priori known externally applied field plus the a priori unknown magnetic field that is produced by plasma currents inside ${{\mathcal {S}}}$, as given by the virtual-casing integral given; and (ii) that this ‘virtual-casing self-consistent’ vacuum field may be obtained as the solution to the system of linear equations resulting from combining (4.3) with (4.1); i.e.

(4.4)\begin{equation} {\mathcal{L}}_{vc}\cdot \boldsymbol{\phi} = \boldsymbol{B}_E, \end{equation}

where we defined the Laplace-virtual-casing matrix, ${\mathcal {L}}_{vc}\equiv {\mathcal {A}}-{\mathcal {B}}$, and $\boldsymbol {B}_E$ is that part of the right-hand side vector given in (4.1) that does not depend on the plasma currents.Footnote ¹¹ We shall elaborate upon these propositions in a future article.

The virtual-casing self-consistent vacuum field will change with variations in the plasma boundary according to

(4.5)\begin{equation} {\mathcal{L}}_{vc} \cdot \delta \boldsymbol{\phi} ={-} \left( \delta {{\mathcal{A}}} - \delta {\mathcal{B}}\right)\cdot \boldsymbol{\phi} + \delta \boldsymbol{B}_E, \end{equation}

where $\delta \mathcal {B} = \partial \mathcal {B} / \partial {{\mathcal {S}}} \cdot \delta {{\mathcal {S}}} + \partial \mathcal {B} / \partial {{\mathcal {D}}} \cdot \delta {{\mathcal {D}}}$.

4.1. Restricted energy functionals

To incorporate the Galerkin method for computing the magnetic fields with the energy functional and with the various combined plasma–coil optimization algorithms discussed in § 5 below, we simplify as follows. We redefine the energy functionals, ${{\mathcal {F}}}_v$, in the plasma volumes to

(4.6)\begin{equation} {{\mathcal{F}}}_v \equiv \int_{{\mathcal{V}}_v} \left( \frac{p_v}{\gamma-1} + \frac{B_v^{2}}{2} \right) \mathrm{d} v, \end{equation}

where here and hereafter the magnetic field, $\boldsymbol {B}_v$, in each region is restricted to be the unique magnetic field with the prescribed helicity and fluxes that is tangential to the boundary and obeys $\boldsymbol {\nabla } \times \boldsymbol {B}_v = \mu _v \boldsymbol {B}_v$. This equation is known as the Beltrami equation and is obtained, see (3.2), as the Euler–Lagrange equation $\delta {{\mathcal {F}}}_V / \delta \boldsymbol {A}_v = 0$. We note that $\boldsymbol {B}_v$ depends only on the geometry of the adjacent interfaces; i.e. $\boldsymbol {B}_v = \boldsymbol {B}_v[\boldsymbol {x}_{v-1},\boldsymbol {x}_v]$.Footnote ¹² The energy functional in the vacuum region is redefined simply as

(4.7)\begin{equation} {{\mathcal{F}}}_+{=} \int_{{{\mathcal{V}}}_+} \frac{B_+^{2}}{2}\, \mathrm{d} v, \end{equation}

where $\boldsymbol {B}_+ = \nabla \varPhi$ is restricted to be the unique magnetic field with the prescribed loop integrals (for the enclosed currents) and is tangential to ${{\mathcal {S}}}$, and that obeys $\boldsymbol {\nabla }\boldsymbol {\cdot }\boldsymbol {\nabla }\varPhi =0$ in ${{\mathcal {V}}}_+$. This field is obtained as the solution to either (4.1) if the total normal field is given or (4.4) if the external normal field is given. With the enclosed currents and ${{\mathcal {D}}}$ held constant, we note that $\boldsymbol {B}_+$ depends only on ${\mathcal {S}}$ and the boundary condition; i.e. $\boldsymbol {B}_+ = \boldsymbol {B}_+[{\mathcal {S}},B_{T,n}]$ or $\boldsymbol {B}_+ = \boldsymbol {B}_+[{\mathcal {S}},B_{E,n}]$.

With these conventions hereafter implied, we omit the explicit dependence of ${{\mathcal {F}}}$ on the vector and scalar potentials and the Lagrange multipliers, and write ${{\mathcal {F}}}={{\mathcal {F}}}_t[\boldsymbol {x}, B_{T,n},{{\mathcal {D}}}]$ or ${{\mathcal {F}}}={{\mathcal {F}}}_e[\boldsymbol {x}, B_{E,n},{{\mathcal {D}}}]$. To distinguish the case when ${B_{T,n}}$ is provided and (4.1) is used to compute the vacuum field, we use ${{\mathcal {F}}}_t$ to denote the energy functional in this case, whereas we instead denote this quantity by ${{\mathcal {F}}}_e$ when ${B_{E,n}}$ is provided and (4.4) is used. We use $\boldsymbol {x}$ to denote the geometry of the $v=1,\ldots N_V$ interfaces, which includes the plasma boundary, i.e. $\boldsymbol {x}_{N_V} = {{\mathcal {S}}}$.

The interface geometry, and by implication the equilibrium magnetic field, is defined by $\delta {{\mathcal {F}}} / \delta \boldsymbol {x} = 0$, which from (3.3) is the equilibrium equation, ${\boldsymbol {F}} \equiv [[p+B^{2}/2]]=0$ across the ideal interfaces. Minima of ${{\mathcal {F}}}$ with respect to ${\boldsymbol {x}}$ can be found using descent-style algorithms, $\partial \boldsymbol {x} / \partial \tau = - \delta {{\mathcal {F}}}/\delta \boldsymbol {x}$, where $\tau$ is an integration parameter, or by Newton-style methods. The geometry of the equilibrium interfaces is considered to be a function of the boundary conditions; i.e. $\boldsymbol {x} = \boldsymbol {x}[B_{T,n},{{\mathcal {D}}}]$ or $\boldsymbol {x} = \boldsymbol {x}[B_{E,n},{{\mathcal {D}}}]$. We shall write the total (tangential) magnetic field, $\boldsymbol {B}_T|_{{{\mathcal {S}}}}$, on ${{\mathcal {S}}}$, which is required for the virtual casing calculation of $B_{P,n}$, as a function of $\boldsymbol {x}$, i.e. $\boldsymbol {B}_T|_{{{\mathcal {S}}}}=\boldsymbol {B}_T|_{{{\mathcal {S}}}}[\boldsymbol {x}]$.

The second derivatives of the restricted energy functions, (4.6) and (4.7), with respect to variations in the interface geometry, $\boldsymbol {x}$, and the boundary conditions, either $B_{T,n}$ of $B_{E,n}$, which are required in the following section, may be constructed from incorporating the derivatives of the magnetic fields as calculated using (4.2) and (4.5) with the expressions presented in § 3.1.Footnote ¹³

5. Combined plasma–coil optimization algorithms

In this section, we consider the construction of the plasma equilibrium and the coil geometry simultaneously with the optimization of the plasma and coil properties.

When we are required to construct the geometry of the coils as extrema of the coil-penalty functional, ${{\mathcal {E}}}$, we assume that the coil geometry satisfies $\delta {{\mathcal {E}}} / \delta \boldsymbol {c} = 0$. When the coil geometry is taken as the independent degree-of-freedom in the optimization, the coil-penalty functional is not required and the Biot–Savart integral is used to compute ${B_{E,n}} = {B_{E,n}}[\boldsymbol {c},{{\mathcal {D}}}]$.

When the normal plasma field, $B_{P,n}$, on ${{\mathcal {D}}}$ is required and the equilibrium magnetic field is known, $B_{P,n}$ is computed using either the total (tangential) magnetic field either immediately inside, ${{\mathcal {S}}}_-$, or immediately outside, ${{\mathcal {S}}}_+$, the plasma boundary, depending on whether $\boldsymbol {B}_T$ is determined using a fixed-boundary calculation, for which the total magnetic field outside ${{\mathcal {S}}}$ may not be known, or a free-boundary calculation, for which it is. If there is a sheet current on the plasma boundary, these will differ. The total (tangential) magnetic field, $\boldsymbol {B}_T|_{{{\mathcal {S}}}}$ on ${{\mathcal {S}}}$, is obtained as an output of the equilibrium calculation, $\delta {{\mathcal {F}}} / \delta \boldsymbol {x} = 0$.

In the following, we consider choosing the independent degrees-of-freedom, $\boldsymbol {z}$, in the combined plasma–coil optimization to be (i) the plasma boundary, $\boldsymbol {z} \equiv {{\mathcal {S}}}$; (ii) the total normal field on ${{\mathcal {D}}}$, $\boldsymbol {z} \equiv B_{T,n}$; (iii) the required normal field on ${{\mathcal {D}}}$, $\boldsymbol {z} \equiv D_n$; and (iv) the coil geometry, $\boldsymbol {z} \equiv \boldsymbol {c}$. For each of these choices, the descent direction depends on the derivatives of the plasma-energy functional, ${{\mathcal {F}}} = {{\mathcal {F}}}_t[\boldsymbol {x}_v, B_{T,n},{{\mathcal {D}}}]$ or ${{\mathcal {F}}} = {{\mathcal {F}}}_e[\boldsymbol {x}_v, B_{E,n},{{\mathcal {D}}}]$ depending on whether the total or external normal field is given; the coil-penalty functional, ${{\mathcal {E}}} = {{\mathcal {E}}}[\boldsymbol {c},D_n,{{\mathcal {D}}}]$; and the virtual-casing calculation of the plasma normal field, $B_{P,n} = B_{P,n}[{{\mathcal {S}}},\boldsymbol {B}_T|_{{{\mathcal {S}}}},{{\mathcal {D}}}]$.

It is helpful to have a tangible example of what plasma and coil properties we wish to optimize. For the plasma optimization, we imagine a plasma property, ${\mathcal {P}}[\boldsymbol {x}]$, that is explicitly a scalar function of the geometry of the flux surfaces and is to be minimized. Generally, most plasma properties depend on the magnetic field, but here we are assuming that the equilibrium magnetic field depends on the geometry of the interfaces and so there is no loss of generality in treating ${\mathcal {P}}$ as an explicit function of only the geometry of the interfaces. Plasma properties of interest include the integrability of the magnetic field, quasi-symmetry properties, the rotational-transform profile, stability properties,etc.

For the coil geometry optimization, we imagine a measure of the coil complexity, ${\mathcal {C}}[\boldsymbol {c}]$, that is a scalar function of the coil geometry and is to be minimized. We might consider the total length of the coils, as shorter coils tend to be cheaper, or the non-planarness of the coils, as measured by ${\mathcal {C}} = \sum _i \oint \tau ^{2}/2 \,\mathrm {d} l$ where $\tau$ is the torsion. Other properties of interest might include a measure of the coil–plasma separation, which explicitly depends on both the coil geometry and the plasma boundary, e.g. ${\mathcal {Q}} = {\mathcal {Q}}[\boldsymbol {x},\boldsymbol {c}]$. It is straightforward to extend the following treatment to include such properties.

For the combined optimization, we imagine a differentiable cost function, $\mathcal {T}=\mathcal {T}[\mathcal {P},\mathcal {C}]$. If ${\mathcal {P}}[\boldsymbol {x}]$ and ${\mathcal {C}}[\boldsymbol {c}]$ are differentiable functions of, respectively, $\boldsymbol {x}$ and $\boldsymbol {c}$,Footnote ¹⁴ then descent algorithms can be implemented if we know the derivatives of $\boldsymbol {x}$ and $\boldsymbol {c}$ with respect to $\boldsymbol {z}$. The derivative of ${\mathcal {T}}$ is

(5.1)\begin{equation} \frac{\partial \mathcal{T}}{\partial \boldsymbol{z}} = \frac{\partial \mathcal{T}}{\partial {\mathcal{P}}} \frac{\partial {\mathcal{P}}}{\partial \boldsymbol{x}} \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{z}} + \frac{\partial \mathcal{T}}{\partial {\mathcal{C}}} \frac{\partial \mathcal{C}}{\partial \boldsymbol{c}} \frac{\partial\boldsymbol{c}}{\partial \boldsymbol{z}}. \end{equation}

Here and hereafter, for notational clarity, we assume that all quantities are discretized, so that functionals of lines, surfaces and volumes become functions of a finite set of parameters that describe those objects, and the infinite-dimensional Frechét derivatives become finite-dimensional derivatives. One cannot be more explicit regarding constructing the derivatives of ${\mathcal {T}}$, ${\mathcal {P}}$ and ${\mathcal {C}}$ until one has stated what these quantities are, this is left to a future article. In the following, we present expressions for $\partial \boldsymbol {x} / \partial \boldsymbol {z}$ and $\partial \boldsymbol {c} / \partial \boldsymbol {z}$ for the different choices of $\boldsymbol {z}$ mentioned above.

5.1. Fixed-boundary optimization

We can consider the independent degree-of-freedom in the optimization to be the plasma boundary, $\boldsymbol {z} \equiv {{\mathcal {S}}}$. The free-plasma-boundary equilibrium calculation described in § 2 reduces to a fixed-plasma-boundary calculation by eliminating the vacuum region; that is, we choose ${{\mathcal {D}}}={{\mathcal {S}}}$ and $B_{T,n}=0$. In this subsection only, because we are choosing ${{\mathcal {D}}}$ to be coincident with ${{\mathcal {S}}}$ to facilitate a fixed-boundary equilibrium calculation, ${{\mathcal {D}}}$ will move during the optimization.

The equilibrium, $\boldsymbol {x}$, satisfies $\partial {{\mathcal {F}}}_t(\boldsymbol {x},0,{{\mathcal {S}}}) / \partial \boldsymbol {x}=0$. We compute $B_{P,n}({{\mathcal {S}}},\boldsymbol {B}_T|_{{{\mathcal {S}}}_-},{{\mathcal {S}}})$ from the virtual-casing integral.Footnote ¹⁵ With $B_{T,n}=0$, the required normal field is $D_n = - B_{P,n}({{\mathcal {S}}},\boldsymbol {B}_T|_{{{\mathcal {S}}}_-},{{\mathcal {S}}})$. The coil geometry satisfies $\partial {{\mathcal {E}}}(\boldsymbol {c}, D_n,{{\mathcal {S}}})/\partial \boldsymbol {c}=0$. Expanding and collecting terms, we obtain

(5.2)\begin{gather} \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{z}} ={-} \left( \frac{\partial^{2} {{\mathcal{F}}}_t}{\partial \boldsymbol{x} \partial \boldsymbol{x}} \right)^{{-}1} \cdot \left( \frac{\partial^{2} {{\mathcal{F}}}_t}{\partial \boldsymbol{x} \partial {{\mathcal{D}}}} \right), \end{gather}

(5.3)\begin{gather}\frac{\partial \boldsymbol{c}}{\partial \boldsymbol{z}} ={-} \left( \frac{\partial^{2} {{\mathcal{E}}}}{\partial \boldsymbol{c} \partial \boldsymbol{c}} \right)^{{-}1} \cdot \left( \frac{\partial^{2} {{\mathcal{E}}}}{\partial \boldsymbol{c} \,\partial D_n} \cdot \frac{\partial D_n}{\partial \boldsymbol{z}} + \frac{\partial^{2} {{\mathcal{E}}}}{\partial \boldsymbol{c} \,\partial {{\mathcal{D}}}} \right), \end{gather}

where

(5.4)\begin{equation} \frac{\partial D_n}{\partial \boldsymbol{z}} ={-} \frac{\partial {B_{P,n}}}{\partial \boldsymbol{x}} \cdot \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{z}}, \end{equation}

and where

(5.5)\begin{equation} \frac{\partial {B_{P,n}}}{\partial \boldsymbol{x}} = \frac{\partial {B_{P,n}}}{\partial {{\mathcal{S}}}} \cdot \frac{\partial {{\mathcal{S}}}}{\partial \boldsymbol{x}} + \frac{\partial {B_{P,n}}}{\partial \boldsymbol{B}_T|_{{{\mathcal{S}}}}} \cdot \frac{\partial \boldsymbol{B}_T|_{{{\mathcal{S}}}_-}}{\partial \boldsymbol{x}}. \end{equation}

The analogous functional derivatives for $\partial ^{2} {{\mathcal {E}}} / \partial \boldsymbol {c} \partial \boldsymbol {c}$ and $\partial ^{2} {{\mathcal {E}}} / \partial \boldsymbol {c} \partial {{\mathcal {D}}}$ are given in (3.16) and (3.21), those for $\partial ^{2} {{\mathcal {F}}}_t/\partial \boldsymbol {x} \partial \boldsymbol {x}$, $\partial ^{2} {{\mathcal {F}}}_t/\partial \boldsymbol {x} \partial {{\mathcal {D}}}$ are given in (3.6) and (3.7), and those for $\partial B_{P,n} / \partial {{\mathcal {S}}}$ and $\partial B_{P,n} / \partial \boldsymbol {B}_T|_{{{\mathcal {S}}}}$ can be found in (3.24) and (3.25).

The coil geometry is constructed by minimizing the coil-penalty functional. With a finite number of discrete coils, the optimal coils will not exactly produce the required magnetic field. It may be beneficial for a combined plasma–coil optimization to minimize the quadratic-flux error $\varphi _2=\varphi _2(D_n,\boldsymbol {c},{{\mathcal {D}}})$ The required derivative is

(5.6)\begin{equation} \frac{\partial \varphi_2}{\partial \boldsymbol{z}} = \frac{\partial \varphi_2}{\partial D_n} \cdot \frac{\partial D_n}{\partial \boldsymbol{z}} + \frac{\partial \varphi_2}{\partial \boldsymbol{c}} \cdot \frac{\partial \boldsymbol{c}}{\partial \boldsymbol{z}} + \frac{\partial \varphi_2}{\partial {{\mathcal{D}}}}. \end{equation}

The analogous functional derivative, $\partial \varphi _2/\partial {{\mathcal {D}}}$, is given in (3.19).

5.2. Generalized fixed-boundary optimization

We can consider the independent degree-of-freedom in the optimization to be the total normal field on the computational boundary, ${\boldsymbol {z}} \equiv B_{T,n}$ on ${{\mathcal {D}}}$. Here ${{\mathcal {D}}}$ is assumed to lie outside the plasma boundary and to remain fixed during the calculation. The equilibrium satisfies $\partial {{\mathcal {F}}}_t(\boldsymbol {x},B_{T,n},{{\mathcal {D}}})/\partial \boldsymbol {x}=0$. The required normal field is $D_n = B_{T,n}-B_{P,n}({{\mathcal {S}}},\boldsymbol {B}_T|_{{{\mathcal {S}}}_+},{{\mathcal {D}}})$. The coil geometry satisfies $\partial {{\mathcal {E}}}(\boldsymbol {c},D_n,{{\mathcal {D}}})/\partial \boldsymbol {c}=0$. Expanding and collecting terms, we obtain

(5.7)\begin{gather} \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{z}} ={-} \left( \frac{\partial^{2} {{\mathcal{F}}}_t}{\partial \boldsymbol{x} \partial \boldsymbol{x}} \right)^{{-}1} \cdot \frac{\partial^{2} {{\mathcal{F}}}_t}{\partial \boldsymbol{x} \partial B_{T,n}}, \end{gather}

(5.8)\begin{gather}\frac{\partial \boldsymbol{c}}{\partial \boldsymbol{z}} ={-} \left( \frac{\partial^{2} {{\mathcal{E}}}}{\partial \boldsymbol{c} \partial \boldsymbol{c}} \right)^{{-}1} \cdot \frac{\partial {{\mathcal{E}}}}{ \partial \boldsymbol{c} \partial D_n} \cdot \frac{\partial D_n}{\partial \boldsymbol{z}}, \end{gather}

where

(5.9)\begin{equation} \frac{\partial D_n}{\partial \boldsymbol{z}} = 1 - \frac{\partial {B_{P,n}}}{\partial \boldsymbol{x}} \cdot \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{z}} \end{equation}

and where

(5.10)\begin{equation} \frac{\partial {B_{P,n}}}{\partial \boldsymbol{x}} = \frac{\partial {B_{P,n}}}{\partial {{\mathcal{S}}}} \cdot \frac{\partial {{\mathcal{S}}}}{\partial \boldsymbol{x}} + \frac{\partial {B_{P,n}}}{\partial \boldsymbol{B}_T|_{{{\mathcal{S}}}}} \cdot \frac{\partial \boldsymbol{B}_T|_{{{\mathcal{S}}}_+}}{\partial \boldsymbol{x}}. \end{equation}

The analogous functional derivatives for $\partial {{\mathcal {E}}}/\partial \boldsymbol {c} \partial D_n$, $\partial B_{P,n} / \partial {{\mathcal {S}}}$ and $\partial B_{P,n} / \partial \boldsymbol {B}_T|_{{{\mathcal {S}}}}$ are given in (3.17), (3.24) and (3.25). Note the use of $\boldsymbol {B}_T|_{{{\mathcal {S}}}_-}$ in (5.5) and $\boldsymbol {B}_T|_{{{\mathcal {S}}}_+}$ in (5.10). The derivative of the quadratic-flux error is

(5.11)\begin{equation} \frac{\partial \varphi_2}{\partial \boldsymbol{z}} = \frac{\partial \varphi_2}{\partial D_n} \cdot \left( 1 - \frac{\partial B_{p,n}}{\partial \boldsymbol{x}} \cdot \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{z}} \right) + \frac{\partial \varphi_2}{\partial \boldsymbol{c}} \cdot \frac{\partial \boldsymbol{c}}{\partial \boldsymbol{z}}. \end{equation}

The analogous functional derivative for $\partial \phi _2 /\partial D_n$ is given in (3.11).

5.3. Quasi-free-boundary optimization

In this case, we consider the independent degree-of-freedom in the optimization to be the required external normal field on the computational boundary, $\boldsymbol {z} \equiv D_n$ on ${{\mathcal {D}}}$. For the equilibrium calculation, we assume that the ‘actual’ external normal field is equal to the required external normal field, ${B_{E,n}} = D_n$; i.e. we assume that the equilibrium satisfies $\partial {{\mathcal {F}}}_e(\boldsymbol {x},D_n,{{\mathcal {D}}})/\partial \boldsymbol {x}=0$. The virtual casing integral is not explicitly required as it is implicitly included in the construction of the virtual-casing self-consistent vacuum field, (4.4). The coil geometry satisfies $\partial {{\mathcal {E}}}(\boldsymbol {c},D_n,{{\mathcal {D}}}) / \partial \boldsymbol {c} = 0$.

Expanding and collecting terms, we obtain

(5.12)\begin{gather} \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{z}} ={-} \left( \frac{\partial^{2} {{\mathcal{F}}}_e}{\partial \boldsymbol{x} \partial \boldsymbol{x}} \right)^{{-}1} \cdot \left( \frac{\partial^{2} {{\mathcal{F}}}_e}{\partial \boldsymbol{x} \partial B_{E,n}} \right), \end{gather}

(5.13)\begin{gather}\frac{\partial \boldsymbol{c}}{\partial \boldsymbol{z}} ={-} \left( \frac{\partial^{2} {{\mathcal{E}}}}{\partial \boldsymbol{c} \partial \boldsymbol{c}} \right)^{{-}1} \cdot \left( \frac{\partial^{2} {{\mathcal{E}}}}{\partial \boldsymbol{c} \partial D_n} \right) \end{gather}

and for the quadratic-flux error

(5.14)\begin{equation} \frac{\partial \varphi_2}{\partial \boldsymbol{z}} = \frac{\partial \varphi_2}{\partial D_n} + \frac{\partial \varphi_2}{\partial \boldsymbol{c}} \cdot \frac{\partial \boldsymbol{c}}{\partial \boldsymbol{z}}. \end{equation}

5.4. Free-boundary (coil) optimization

In this final case, the independent degree-of-freedom in the optimization is the coil geometry, $\boldsymbol {z} \equiv \boldsymbol {c}$. The equilibrium satisfies $\partial {{\mathcal {F}}}_e(\boldsymbol {x},B_{E,n},{{\mathcal {D}}})/\partial \boldsymbol {x}=0$, where $B_{E,n}=B_{E,n}[\boldsymbol {c},{{\mathcal {D}}}]$ is computed using the Biot–Savart integral.

We obtain

(5.15)\begin{equation} \frac{\partial \boldsymbol{x}}{\partial \boldsymbol{z}} ={-} \left( \frac{\partial^{2} {{\mathcal{F}}}_e}{\partial \boldsymbol{x} \partial \boldsymbol{x}} \right)^{{-}1} \cdot \frac{\partial^{2} {{\mathcal{F}}}_e}{\partial \boldsymbol{x} \partial B_{E,n}} \cdot \frac{\partial B_{E,n}}{\partial \boldsymbol{c}}. \end{equation}

The analogous functional derivative for $\partial B_{E,n}/\partial \boldsymbol {c}$ is given in (3.13). The coil geometry is given directly, $\boldsymbol {c} = \boldsymbol {z}$.

Since the input to the equilibrium calculation, ${B_{E,n}}$, is computed directly from the actual coils there is no coil error, $\varphi _2 = 0$.

A summary of the 4 different combined plasma-coil optimization algorithms can be found in Table 1.

Table 1. Overview table of the different optimization approaches.