Combined plasma-coil optimization algorithms

Combined plasma-coil optimization approaches for designing stellarators are discussed and a new method for calculating free-boundary equilibria is proposed. Four distinct categories of stellarator optimization, two of which are novel approaches, are the fixed-boundary optimization, the generalized fixed-boundary optimization, the quasi free-boundary optimization, and the free-boundary (coil) optimization. These are described using the multi-region relaxed magnetohydrodynmics (MRxMHD) energy functional, the Biot-Savart integral, the coil-penalty functional and the virtual casing integral, and their derivatives. The proposed free-boundary equilibrium calculation differs from existing methods in how the boundary-value problem is posed, and for the new approach it seems that there is not an associated energy minimization principle because a non-symmetric functional arises. We propose to solve the weak formulation of this problem using a spectral-Galerkin method, and this will reduce the free-boundary equilibrium calculation to something comparable to a fixed-boundary calculation. In our discussion of combined plasma-coil optimization algorithms, we emphasize the importance of the stability matrix.


Introduction
The design space for stellarators is larger than that of tokamaks because stellarators exploit three-dimensional (3D) magnetic fields, by which it is meant that there is no continuous (e.g. rotational) symmetry, whereas tokamaks are notionally axisymmetric (2D) (Helander 2014). The larger design space allows more freedom in the geometry of the plasma boundary. Geometry affects several important plasma properties such as stability and transport (Helander et al. 2012), including turbulent transport (Hegna et al. 2018). If the search includes 3D configurations, generally, one may expect that it is more likely to find a feasible fusion device.
A particularly advantageous feature of stellarators is that the rotational-transform, which is essential for confinement in toroidal configurations, can be provided externally (Mercier 1964), either by current-carrying coils or by permanent magnets. This reduces or even eliminates the need of generating toroidal plasma currents, which can lead to problematic disruptions.
Along with advantages, 3D configurations can have drawbacks. A feasible fusion device must possess good particle confinement and the plasma equilibrium must be supported by an external set of coils or magnets that do not pose unrealistic engineering challenges (Grieger et al. 1992). Stellarators are guaranteed neither. Plasmas in early stellarator designs were not well confined because of neoclassical particle losses caused by unconfined orbits (Galeev et al. 1969;Beidler et al. 2011). Geometrically distorted coils complicated the engineering of the cancelled NCSX experiment (Neilson et al. 2010) and caused delays in the construction of Wendelstein 7-X (Riße 2009).
Nonetheless, with careful exploitation of the large design space of 3D configurations confinement in stellarators can be significantly improved, e.g., by using quasi-symmetric (Nührenberg & Zille 1988;Boozer 1995;Henneberg et al. 2020) or quasi-isodynamic fields (Gori et al. 1996;Landreman & Catto 2012). Perfectly quasi-isodynamic fields have vanishing bootstrap current, which allows for simultaneous optimization of both neoclassical transport and small toroidal net current (Helander & Nührenberg 2009), which can be beneficial for stability and is necessary if an island divertor is to be employed. The 3D optimization effort is encouraged by recent Wendelstein 7-X results (Klinger et al. 2019), which was optimized for neoclassical transport (Nührenberg & Zille 1986). There are new methods in coil-design, such as using permanent magnets Zhu et al. 2020;Landreman & Zhu 2020) and designing coils with generous tolerances (Lobsien et al. 2018(Lobsien et al. , 2020.
With these developments in understanding and designing 3D configurations, the question is not whether we should search the large 3D configuration space but how. With an order of magnitude more degrees of freedom, give or take, it is significantly more difficult to search the 3D configuration space. To simultaneously achieve the desired properties of a feasible stellarator, we need efficient optimization approaches, and we need to understand the solution space.
An essential component of stellarator optimization is the evaluation of equilibrium properties. † Equilibrium calculations are conveniently divided into two types, namely fixed-boundary and free-boundary. ‡ The fixed-boundary approach requires the geometry of the plasma boundary to be provided as input information (Hirshman & Whitson 1983;Bauer et al. 1984;Hudson et al. 2012). Experience suggests that fixed-boundary calculations are faster and more robust than their free-boundary counterpart. In reality, however, the geometry of the plasma boundary is not known a priori. Freeboundary equilibrium calculations require as input the external magnetic field, and the self-consistent plasma boundary is determined as part of the equilibrium calculation (Hirshman et al. 1986;Hudson et al. 2020). Existing free-boundary equilibrium codes invariably perform additional iterations compared to their fixed-boundary analogs. The present algorithm in free-boundary SPEC , for example, performs a Picard-style iteration over the fixed-boundary code in order to determine the plasma boundary that is consistent with the external field. The speed of existing fixed-boundary codes is advantageous for stellarator optimization because the equilibrium is typically computed at each iteration.
One might ask, what is the "best" approach for stellarator optimization? We do not expect that there will be a be-all and end-all of stellarator optimization algorithms. We should embrace a variety of methods. In this paper, we seek to elucidate the mathematical structure of the various inter-related calculations that underpin the integrated stellarator optimization problem, which we hereafter call the combined plasma-coil optimization. The overall numerical efficiency of the combined plasma-coil optimization problem cannot be understood without considering how the equilibrium calculation, the coil calculation, and the optimization calculation communicate with each other. This needs to be understood from both a mathematical and an algorithmic perspective. † We consider vacuum approximation or large aspect ratio expansions to be types of equilibrium calculations.
‡ Here the word "boundary" refers to the plasma boundary; later we shall generalize this terminology to include equilibrium calculations with a fixed computational boundary but with a plasma boundary that moves during the calculation.
Existing algorithms can loosely be sorted into two categories, which may be called the "two-step" method and the "direct-coil" method, which we now describe. ¶ The two-step approach separates the equilibrium optimization calculation from the problem of finding coils (Nührenberg & Zille 1986). In the first step, the plasma boundary is varied to optimize the desired equilibrium properties, such as rotational transform profile, MHD stability, etc. (Beidler et al. 1990;Zarnstorff et al. 2001;Strickler et al. 2002a;Henneberg et al. 2019;Bader et al. 2020). In the second step, the geometry of an external set of coils that provides the required external field is determined. The coil geometry must meet certain engineering criteria; the coils cannot have too small curvature and they must not intersect each other, for example. This stellarator optimization can be performed with STELLOPT , ROSE , or SIMSOPT (Simons Hidden Symmetries Collaboration, https://github.com/hiddenSymmetries/simsopt), which are all under active development. The two-step approach was used to design Wendelstein 7-X (Beidler et al. 1990), CHS-qa (Okamura et al. 2001), NCSX (Nelson et al. 2003), HSX (Canik et al. 2007), ESTELL (Drevlak et al. 2013), and CFQS (Shimizu et al. 2018).
An advantage of the two-step approach is that it is primarily the fusion-relevant performance of the plasma that is most important; so, it is reasonable to optimize the plasma performance first without compromise; that is, without regard to the complexity of the coils. After all, a stellarator that produces fusion power with complicated coils is incomparably more commercially valuable than a stellarator with simple coils that does not. An advantage of using fixed-boundary calculations to optimize the equilibrium properties is the availability of fast and robust fixed-boundary equilibrium codes.
Using fixed-plasma-boundary calculations for the plasma equilibrium leads to an inverse problem for the coil geometry: the required magnetic field is known, and the Biot-Savart law must be inverted to obtain the coil geometry. Many different codes have been developed to solve the inverse problem for the coils, e.g., NESCOIL (Merkel 1987), ONSET (Drevlak 1998), COILOPT (Brown et al. 2015), FOCUS (Zhu et al. 2017), REG-COIL (Landreman 2017), OMIC (Singh et al. 2020), andFOCUSADD (McGreivy et al. 2020). There are codes that solve the inverse problem for a set of permanent magnets Zhu et al. 2020;Landreman & Zhu 2020).
A disadvantage of the two-step approach is that not all plasma boundaries can be produced by a discrete set of coils or magnets that are easy to build. If there is no consideration of the coil and magnet complexity until after the plasma equilibrium has been chosen, complicated coils can result. It may be the case that small changes in the plasma geometry that result in small improvements in the plasma performance may also result in large disadvantageous changes in the coil geometry. Integrated algorithms that incorporate both plasma performance and coil complexity are therefore desirable (Boozer 2015;Landreman & Boozer 2016;Drevlak et al. 2019).
Measures of coil complexity can be included in the optimization cost function using the two-step method Carlton-Jones et al. 2020). This comes at the cost of computing the coil geometry or some approximation of it at every stage of the fixed-boundary plasma optimization.
Another approach for the combined plasma-coil design is the direct coil optimization using a free-boundary equilibrium code (Strickler et al. 2002b;Hudson et al. 2002). The direct coil optimization approach avoids the need for the inverse coil code because the coil geometry is taken to be the independent degree of freedom in the optimization, ¶ Note that this is loose terminology and there are numerous variations of these methods. We shall give a more precise description of various algorithms in Sec. 5.
i.e., the coil geometry is assumed to be known. The external magnetic field that is required as input by the free-boundary equilibrium calculation is computed using the Biot-Savart formula. Metrics of coil complexity can easily be included in the optimization cost function.
There are direct coil optimization approaches that dispense with the equilibrium calculation and instead optimize the properties of the vacuum field, which may be thought of the trivial (zero plasma-current, zero pressure) plasma equilibrium state. This should be sufficient for designs that do not target high plasma pressure; however, finite-β effects are important: the pressure-induced Shafranov shift (Shafranov 1963) can affect stability and confinement. A related direct coil optimization uses a near-axis analytic approximation to the equilibrium state (Giuliani et al. 2020), which might be sufficient for large aspect-ratio configurations but perhaps not so for "compact" stellarators. Also, recent work has exploited analytical gradient information and adjoint methods for a variety of optimization problems Paul et al. 2018;Antonsen et al. 2019;Paul et al. 2019;Carlton-Jones et al. 2020).
In this paper, we review and categorize different combined plasma-coil optimization. We discuss free-boundary calculations and propose an alternative method to calculate the virtual-casing self-consistent vacuum magnetic field, which reduces the cost of a freeboundary calculation to that of a fixed-boundary calculation. In Sec. 2, the multi-region relaxed magnetohydrodynamic (MRxMHD) energy functional, the coil-penalty functional and the virtual casing integral are described. In Sec. 3, we present the relevant variational calculus for the free-boundary equilibrium calculation, the coil geometry optimization and the virtual casing integral that are employed in the combined plasma-coil optimization algorithms discussed in Sec. 5. In Sec. 4, we discuss the construction of the magnetic field using a (weak) Galerkin method. In Sec. 5, we examine various fixed-and free-boundary algorithms for the combined plasma-coil optimization that differ in which quantity, which we denote with z, is chosen to be the independent degree-of-freedom, and derive the derivatives of the plasma equilibrium and the coil geometry with respect to z in terms of derivatives of the plasma-energy functional, the coil-penalty functional, and the virtualcasing integral. † Finally, in Sec. 6, we discuss the results of this paper.

Plasma equilibria and supporting coils
This paper builds upon the construction of the plasma equilibrium as a stationary point of the multi-region relaxed magnetohydrodynamic (MRxMHD) energy functional, which we now describe. ‡ The computational domain, Ω ⊂ R 3 , is considered to be a solid torus with computational boundary D ≡ ∂Ω. The latter is a prescribed toroidal surface, which unless explicity stated otherwise is held fixed throughout the calculation. The computational domain contains the plasma volume, V ⊂ Ω, a smaller solid torus enclosed by the plasma boundary S ≡ ∂V. The plasma volume is further partitioned into v = 1, . . . , N V nested toroidal sub-volumes, V v , which are separated by a set of ideal interfaces, I v so that ∂V v = I v ∪ I v−1 , with the outermost interface coinciding with the plasma boundary, I NV = S. The innermost sub-region is a simple (solid) torus, and each other sub-region † If the geometry of the coils is chosen to be the independent degree-of-freedom, then the Biot-Savart integral, Eqn. 2.5, is required instead of the coil-penalty functional.
‡ We have chosen to use the MRxMHD functional and not the ideal MHD functional because the ideal MHD functional encounters problems near rational surfaces. However it is straightforward to extend the work of this paper to use ideal MHD. In fact, we believe that most of the results of this paper can be extrapolated to ideal MHD.
is a toroidal annulus which may be thought of as a "hollow" torus. In each V v , the magnetic field is written as the curl of the magnetic vector potential, B v = ∇ × A v . In the innermost toroidal region, the toroidal flux, computed as a poloidal loop integral of A 1 on 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 A v on I v and I v−1 , are constrained. In each region, the magnetic helicity, defined as the volume integral of A v ·B v , is also constrained. The magnetic field in each region is taken to be tangential to that region's boundary, n · B v = 0 on I v where n is a unit outward normal vector, but otherwise the topology of the fieldlines is unconstrained and magnetic islands and irregular field lines are allowed in each V v .
For free-boundary calculations, by which we mean that the plasma boundary is allowed to move, we include the vacuum region V + = Ω \ V, with the inner boundary of V + coinciding with S and its outer boundary with D. It is on 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 V + . We may employ a scalar potential and write the vacuum field as B + = ∇Φ. † To obtain a unique solution, constraints are imposed on the linking and net toroidal plasma currents, which are computed as loop integrals of B + . We take B + to be tangential to the plasma surface, n · ∇Φ| S = 0 where n is a unit outward normal vector on S. The normal field on the computational boundary D, B T,n ≡n · ∇Φ| D wheren is normal to 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 = 0.
The equilibria that we consider are stationary points of the following energy functional, with p v the pressure in the v-th region ‡, γ is the specific heat ratio (adiabatic index), and H v is the prescribed value of the magnetic helicity in V v . The constraints on the toriodal and poloidal fluxes in the plasma volumes, V v , are implied; that is, the fields B v that we consider in Eqn. 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 µ v . The loop integral constraints on ∇Φ are also implied: the allowed Φ in Eqn. 2.3 is restricted to the class of multi-valued scalar functions with prescribed loop integrals of ∇Φ. We writeΦ ≡ Φ| D to denote the scalar potential evaluated on D to distinguish it from Φ ≡ Φ| S on S. Similarly, we write ds to denote an infinitesimal surface element on D.
We may call F the free-plasma-boundary fixed-computational-boundary generalized- † We may also use the magnetic vector potential to describe the magnetic field in V+, and this can have advantages ), but we will not discuss this possibility here.
‡ Here, the pressure is understood in units of µ0 × [J/m 3 ] where µ0 = 4π × 10 −7 [H/m] is the vacuum permeability 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 F simply as the energy functional.
Using the virtual-casing principle (Shafranov & Zakharov 1972), the normal plasma field, B P,n , produced by currents within the plasma volume is computed on the computational boundary as where r ≡ x −x for x ∈ S andx ∈ D, and we use ds to denote a surface element of S. Here, B T | S+ = ∇Φ| S+ is the total magnetic field immediately outside of S. To accommodate for the possible existence of sheet currents lying on the plasma boundary, we must use in Eqn. 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 Sec. 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, Φ.
The difference between the total normal magnetic field and the plasma normal field on 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 1966). This quantity, D n ≡ 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 Sec. 5, we imagine a typical and easily differentiable example; namely, that the external magnetic field is provided by a set of i = 1, . . . , N C current-carrying one-dimensional curves (filaments) with geometry c i that produce a magnetic field given by the Biot-Savart law, where I i is the current in the i-th coil and dl 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π/µ 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 1987;Landreman 2017), or to the case where a "finite-build" approximation is used for the coils (Singh et al. 2020;McGreivy et al. 2020;Li et al. 2020), or when permanent magnets are included Zhu et al. 2020;Landreman & Zhu 2020). 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 1987;Landreman 2017;Zhu et al. 2018;Hudson et al. 2018), is the normal component of B E on D produced by the filamentary coils, where we use c to denote the geometry of all the coils; i.e., c ≡ {c i , i = 1, . . . N C }. We include a regularization term, L = L[c], which we take to be the total length of the coils, and ω is a scalar penalty.
The geometry of the optimal coils is given by δE/δc = 0. This equation may be solved using descent algorithms (Zhu et al. 2017;Hudson et al. 2018) or Newton-style methods . Practically, it is not possible to obtain a coil set that exactly produces the required normal field and generally ϕ 2 = 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 δE/δc = 0, ϕ 2 measures what we call the coil error. † In the above, we have described the three basic components of the combined plasma-coil optimization algorithms that we consider in Sec. 5; namely, the energy functional, F , the Biot-Savart law, the coil-penalty functional, E, and the evaluation of the normal plasma field, B P,n , on D using the virtual-casing functional. In the following section, Sec. 3, we present the first and second variations of F , E and B P,n that will be required in Sec. 5. In Sec. 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 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 D is given. We can solve for the virtual-casing self-consistent vacuum magnetic field using a Galerkin method. In Sec. 5, we show that the derivatives described in Sec. 3 and the Galerkin construction of the magnetic fields may be combined to construct a variety of combined plasma-coil stellarator optimization algorithms.

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 Sec. 5. The variational calculus provides useful insight into the physics and mathematics of the different optimization problems, as we shall comment upon in Sec. 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 to an energy minimization principle. However, as we will see in Sec. 3.4, the weak formulation of certain boundary value problems do not necessarily result in self-ajoint 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 † Note: this is not a "numerical" error, in the sense that this error will not decrease as the numerical resolution of the discrete set of coils increases. The coil error results from the inevitable requirement that the external field is provided by a discrete set of coils (or permanent magnets Zhu et al. 2020;Landreman & Zhu 2020)) which also are constraint by other (engineering) properties. 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.

Energy functional
In the plasma volumes, we write B v = ∇ × A v , and thus δB v = ∇ × δA v . We use that 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 In MRxMHD, the mass is constrained in each volume V v by the isentopic ideal-gas is a constant. The first variation with respect to deformations of the volume boundary is In vacuum, the first variation in F + with respect to variations in Φ is Setting the functional derivative, 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 F v is given by where the tilde (∼) 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 Eqn. 3.8), one has to be careful since cross-terms, including integrals over two different interfaces, appear in the γp v term. In vacuum, the second variation with respect to the scalar potential is given by: Collecting all the different terms for the second variation, we obtain where we used the equilibrium expressions; we included variations of both interfaces; we used B · ((B × ∇) × n) = B 2 ∇ · n − B · ∇n · B; and we set ξ v andξ v proportional to the normal vector 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 . †

Coil-penalty functional
The coil-penalty functional, E defined in Eqn. 2.6 is considered to be a functional of the coil geometry, c, the required normal magnetic field, D n , on the computational boundary, D, and on D itself, which is however kept fixed here; i.e., E = E[c, D n ]. We may formally write the first variation, δE, with respect to variations δc, and δD n as where α is an angle-like curve parameter, c i (α + 2π) = c i (α). The functional derivatives are where c ′ i is the derivative of the i-th coil with respect to α. We have written R i ≡ 3 r i r i /r 5 i − I/r 3 i and defined R i,n = R i · n, where r i is the displacement between a point on the i-th coil and the evaluation point on D, and I is the identity tensor or idemfactor. The curvature of the i-th coil is κ i ≡ c ′ i × c ′′ i /c ′3 i and the functional derivatives are generalized expressions of the ones presented in (Dewar et al. 1994;Hudson et al. 2018). We used the expression for the variation of the normal component of the magnetic field † For ideal MHD the first and second variations of the energy functional are well known (Freidberg 2014). Codes like TERPSICHORE and CAS3D (Anderson et al. 1990;Schwab 1993) can compute the second variation of the ideal MHD functional.
with respect to the coil geometry (3.12) with the functional derivatives (3.13) The required second variations of E are ) 14) ( 3.17) 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 E are not required in Sec. 5 and are omitted for brevity. We only need the variation with respect to the geometry of the computational boundary, D, when D n = −B P,n = −B P ·n. Since there is an explicit dependence with respect ton, we calculate the combined variation of E and B P,n with respect to variations in D, which allows us to express the results in compact form. The first variation is where we used Eqn. 12 from (Dewar et al. 1994). The notation f s = (I −nn) · f denotes the surface projection of an arbitrary vector or tensor, f , onto a surface, which in this case is the computational boundary. The second variation of E with respect to D and c i is where the mean curvature of D is H ≡ −n (∇ ·n). Eqn. 3.21 is a more general form of Eqn. 12 of  where D n = 0.
3.3. Plasma normal field The normal plasma field, B P,n , on D, as defined in Eqn. 2.4, is considered to be a functional of the plasma boundary, S, the total (tangential) magnetic field, B T | S on S, and on D, which is kept fixed; i.e., B P,n = B P,n [S, B T | S ]. The normal component of the magnetic field produced by plasma currents, B P,n , given in Eqn. 2.4 may be written as where the anti-symmetric tensor, M ≡ (rn −n r)/r 3 , is defined, where r = x −x for x ∈ S andx ∈ D. The first variation setting δD = 0 is where we used again Eq. 12 from (Dewar et al. 1994) for the first expression. Note that only M and R depend onx. The relevant variation with respect to the computational boundary D is presented in the previous section. The second variations are not required in the following.

Supplied external field
To enable efficient free-boundary optimization algorithms, which we describe in Sec. 5, we revisit the energy functional; in particular, we pay attention to the boundary condition for the normal field on D.
It is perfectly legitimate to prescribe the total normal field, B T,n , on D. Indeed, every fixed-boundary equilibrium calculation does as much. As described in Sec. 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 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 [B T | S ] = B E,n , where B P,n [B T | S ] may be considered to be a linear, nonlocal operator acting on the tangential total field on the plasma boundary B T | S , which is only known after the equilibrium has been computed. In the current implementation of free-boundary SPEC , 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 B + = ∇Φ such that ∇ · B + = ∇ · ∇Φ = 0 on V + satisfying boundary conditions n · ∇Φ = 0 on S andn · ∇Φ − B P,n [∇Φ| S ] = B E,n on D. This translates into the weak problem of finding Φ such that where ψ is an arbitrary (test) function. An important solvability condition is that the normal field of the coils is consistent with a divergence-free field, D B E,n ds = 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 Eqn. 3.26 spoils the required symmetry of the functional. The first variation of a quadratic energy functional necessarily results in symmetric bi-linear forms.
Energy principles do exist for exterior Neumann problems in the form of boundary integral methods (Giroire & Nédélec 1978). Their formulation however requires the plasma boundary and the computational boundary to coincide. The NESTOR code is a good example (Merkel 1986), as well as BIEST (Malhotra et al. 2019). There are several disadvantages and challenges in implementing these energy principles. First, as seen in Eqn. 3.26 in the limit D → 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 Eqn. 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 †.
In the following section, Sec. 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 selfconsistent vacuum field is obtained as the solution to a set of linear equations, given below in Eqn. 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.

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. 2012Qu et al. 2020). In this paper, we restrict our attention to the problem of constructing B + .
A standard numerical approach for finding stationary points of F + is to represent the scalar potential using a mixed Chebyshev-Fourier representation; e.g., Φ(s, θ, ζ) = I θ + G ζ + l,m,n Φ l,m,n T l (s) exp(imθ − inζ), where (s, θ, ζ) 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 integrals/2π. When this representation is inserted into Eqn. 2.2, and assuming that S and D do not change, the problem of calculating the vacuum field amounts to solving where φ represents the vector of independent degrees-of-freedom, namely the Φ l,m,n and the Lagrange multipliers that may be used to enforce the constraints, and A is a † Readers familiar with VMEC will recognize this as the task of computing the magnetic field on a Cartesian grid from coils via Biot-Savart law, known as mgrid file.
symmetric matrix whose elements are the second derivatives of F + with respect to these degrees-of-freedom and involve volume integrals of coordinate metric elements and the Chebyshev-Fourier functions. The matrix A depends on the geometry of S and D; i.e., A = A[S, D]. The right-hand-side vector, B T , contains the Fourier harmonics of B T,n . The system of linear equations given in Eqn. 4.1 is essentially a discrete realization of Laplace's equation, ∇·∇Φ = 0, with the supplied boundary conditions and loop integrals, written in matrix form. We can determine how the vacuum magnetic field varies with S, D and B T,n using matrix perturbation methods, where δA = ∂A/∂S · δS + ∂A/∂D · δD. From Eqn. 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 B + | S+ . The immediate outside of S is the inner boundary of the vacuum region, V + , and the magnetic field in V + , namely B + , is the gradient of the scalar potential. The Fourier harmonics, B P,n , of the plasma field on D are linear in the degrees-of-freedom of the scalar potential; that is, we may write where B P represents the vector of Fourier harmonics of B P,n , and B is a matrix derived from inserting the mixed Chebyshev-Fourier representation for Φ into Eqn. 2.4. We note that B depends only on the geometry of S and D; i.e., B = B [S, D].
Provided that the external field B E is "tame" (in particular D B E,n ds = 0), we put forward the following propositions: (i) for toroidal annular domains of practical interest, with the constraint that B T · n = 0 on the inner boundary, S, there exists a normal magnetic field on the outer boundary, 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 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 Eqn. 4.3 with Eqn. 4.1; i.e.
where we defined the Laplace-virtual-casing matrix, L vc ≡ A − B, and B E is that part of the right-hand-side vector given in Eqn. 4.1 that does not depend on the plasma currents. † 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) where δB = ∂B/∂S · δS + ∂B/∂D · δD.

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 Sec. 5 below, we simplify as follows. We re-define the energy functionals, F v , in the plasma volumes to (4.6) where here and hereafter the magnetic field, 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 ∇ × B v = µ v B v . This equation is known as the Beltrami equation and is obtained, see Eqn. 3.2, as the Euler Lagrange equation δF V /δA v = 0. We note that B v depends only on the geometry of the adjacent interfaces; i.e., The energy functional in the vacuum region is redefined simply as where B + = ∇Φ is restricted to be the unique magnetic field with the prescribed loop integrals (for the enclosed currents) and is tangential to S, and that obeys ∇ · ∇Φ = 0 in V + . This field is obtained as the solution to either Eqn. 4.1 if the total normal field is given or Eqn. 4.4 if the external normal field is given. With the enclosed currents and D held constant, we note that B + depends only on S and the boundary condition; i.e., With these conventions hereafter implied, we omit the explicit dependence of F on the vector and scalar potentials and the Lagrange multipliers, and write To distinguish the case when B T,n is provided and Eqn. 4.1 is used to compute the vacuum field, we use F t to denote the energy functional in this case, whereas we instead denote this quantity by F e when B E,n is provided and Eqn. 4.4 is used. We use x to denote the geometry of the v = 1, . . . N V interfaces, which includes the plasma boundary; i.e., x NV = S.
The interface geometry, and by implication the equilibrium magnetic field, is defined by δF /δx = 0, which from Eqn. 3.3 is the equilibrium equation, F ≡ [[p + B 2 /2]] = 0 across the ideal interfaces. Minima of F with respect to x can be found using descent-style algorithms, ∂x/∂τ = −δF /δx, where τ 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., x = x[B T,n , D] or x = x[B E,n , D]. We shall write the total (tangential) magnetic field, B T | S , on S, which is required for the virtual casing calculation of B P,n , as a function of x; i.e., The second derivatives of the restricted energy functions, Eqn. 4.6 and Eqn. 4.7, with respect to variations in the interface geometry, 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 Eqn. 4.2 and Eqn. 4.5 with the expressions presented in Sec. 3.1. ‡

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. † The energy functional also depends on the mass and entropy profiles, the flux profiles, and the helicity profiles; but, in this article we assume that these are given. In practice, these must be chosen to be consistent with an assumed model of transport; or, so that the parallel current or rotational transform profiles match experimentally measured profiles, for example, or so that pressure gradients do not coincide with rational surfaces.
‡ All the required derivatives have been implemented in SPEC (Hudson et al. 2012, 2020).
When we are required to construct the geometry of the coils as extrema of the coilpenalty functional, E, we assume that the coil geometry satisfies δE/δ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 [c, D].
When the normal plasma field, B P,n , on 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, S − , or immediately outside, S + , the plasma boundary, depending on whether B T is determined using a fixed-boundary calculation, for which the total magnetic field outside 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, B T | S on S, is obtained as an output of the equilibrium calculation, δF /δx = 0.
In the following, we consider choosing the independent degrees-of-freedom, z, in the combined plasma-coil optimization to be (i) the plasma boundary, z ≡ S; (ii) the total normal field on D, z ≡ B T,n ; (iii) the required normal field on D, z ≡ D n ; and (iv) the coil geometry, z ≡ c. For each of these choices, the descent direction depends on the derivatives of the plasma-energy functional, 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, P[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 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, C[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 nonplanar-ness of the coils, as measured by C = i τ 2 /2 dl where τ 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., It is straightforward to extend the following treatment to include such properties.
For the combined plasma-coil optimization, we imagine a differentiable cost-function, 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 T , P and C until one has stated what these quantities are, this is left to a future article. In the following, we present expressions for ∂x/∂z and ∂c/∂z for the different choices of z mentioned above.

Fixed-boundary optimization
We can consider the independent degree-of-freedom in the optimization to be the plasma boundary, z ≡ S. The free-plasma-boundary equilibrium calculation described in Sec. 2 reduces to a fixed-plasma-boundary calculation by eliminating the vacuum region; that is, we choose D = S and B T,n = 0. In this subsection only, because we are choosing D to be coincident with S to facilitate a fixed-boundary equilibrium calculation, D will move during the optimization.
The equilibrium, x, satisfies ∂F t (x, 0, S)/∂x = 0. We compute B P,n (S, B T | S− , S) from the virtual-casing integral. † With B T,n = 0, the required normal field is D n = −B P,n (S, B T | S− , S). The coil geometry satisfies ∂E(c, D n , S)/∂c = 0. Expanding and collecting terms, we obtain and where The analogous functional derivatives for ∂ 2 E/∂c∂c and ∂ 2 E/∂c∂D are given in Eqn. 3.16 and Eqn. 3.21, those for ∂ 2 F t /∂x∂x, ∂ 2 F t /∂x∂D are given in Eqn. 3.6 and Eqn. 3.7, and those for ∂B P,n /∂S and ∂B P,n /∂B T | S can be found in Eqn. 3.24, and Eqn. 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 ϕ 2 = ϕ 2 (D n , c, D) The required derivative is The analogous functional derivative, ∂ϕ 2 /∂D, is given in Eqn. 3.19.

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, z ≡ B T,n on D. D is assumed to lie outside the plasma boundary and to remain fixed during the calculation. The equilibrium satisfies ∂F t (x, B T,n , D)/∂x = 0. The required normal field is D n = B T,n − B P,n (S, B T | S+ , D). The coil geometry satisfies ∂E(c, D n , D)/∂c = 0. Expanding and collecting terms, we † If D = S, there is a singularity in the integrand of the virtual casing equation. This can be overcome by subtracting an integral which posses the same type of singularity and adding the analytic expression of that subtracted integral (Merkel 1986). (5.9) and where (5.10) The analogous functional derivatives for ∂E/∂c∂D n , ∂B P,n /∂S, and ∂B P,n /∂B T | S are given in Eqn. 3.17,Eqn. 3.24,and Eqn. 3.25. Note the use of B T | S− in Eqn. 5.5 and B T | S+ in Eqn. 5.10. The derivative of the quadratic-flux error is The analogous functional derivative for ∂φ 2 /∂D n is given in Eqn. 3.11.

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, z ≡ D n on 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 ∂F e (x, D n , D)/∂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, Eqn. 4.4. The coil geometry satisfies ∂E(c, D n , D)/∂c = 0.
Expanding and collecting terms, we obtain and for the quadratic-flux error (5.14)

Free-boundary (coil) optimization
In this final case, the independent degree-of-freedom in the optimization is the coil geometry, z ≡ c. The equilibrium satisfies ∂F e (x, B E,n , D)/∂x = 0, where B E,n = B E,n [c, D] is computed using the Biot-Savart integral.
We obtain The analogous functional derivative for ∂B E,n /∂c is given in Eqn. 3.13. The coil geometry is given directly, c = z.
Since the input to the equilibrium calculation, B E,n , is computed directly from the actual coils there is no coil error, ϕ 2 = 0.

Discussion
In this paper, we have categorized and investigated four different combined plasmacoil optimization approaches and proposed a novel method for improving free-boundary equilibrium calculation. We began by summarizing all the functional derivatives of the MRxMHD energy, † the coil-penalty, and the virtual-casing integral needed for a combined plasma-coil optimization. We emphasized that absence of an energy-principle formulation with an explicit boundary condition on the normal component of the external magnetic field on the computational boundary if it does not coincide with the plasma boundary. This construction of the field in the vacuum region is advantageous for the free-boundary optimization approach. For this special case, we proposed solving for the required field using a weak formulation. Finally, we explicitly stated which derivatives are necessary for the four distinct combined plasma-coil optimization algorithms: fixedboundary optimization, generalized fixed-boundary optimization, quasi-free-boundary optimization, and free-boundary (coil) optimization. To the best of our knowledge, all existing stellarator optimization algorithms can be grouped in either the fixed-boundary optimization or the free-boundary optimization approach. The collection of these four distinct optimization algorithms helps to clarify certain intrinsic problems of combined plasma-coil optimization, as we will discuss later in this section.
The novel proposed approach for calculating the virtual-casing self-consistent vacuum field will reduce the cost of the free-boundary calculation to something comparable to that of a fixed-boundary calculation. A full-Newton method, with analytic derivatives, can be used for both. In the proposed approach, in Eqn. 4.4 appears the Laplace-virtualcasing matrix L vc = A − B, which is an important matrix that deserves some discussion. This matrix, in a vague sense, "connects" the plasma to the external magnetic field. Since it is not symmetric, it cannot directly follow from an energy principle, although there are round-about ways to force this to happen (e.g. by introducing adjoint degrees of freedom). Investigation of the well-posedness of the Laplace-virtual-casing problem is the subject of future publication.
On the basis of the summary of the four distinct optimization algorithms presented in Sec. 5, we now discuss the differences and potential advantages and disadvantages of the various optimization algorithms described above. The independent degree-of-freedom is (i) a two-dimensional surface, S(θ, ζ), where θ and ζ are poloidal and toroidal angles, embedded in three-dimensional space for the fixed-boundary approach; (ii) a scalar function, B T,n (θ, ζ), with the constraint that the net-flux is zero for the generalized fixed-boundary approach; (iii) a similar function, D n (θ, ζ), for the quasi free-boundary approach; and (iv) a discrete set of one-dimensional curves embedded in three-dimensional space, c i (θ) for i = 1, . . . , N C , for the free-boundary (coil) approach. Note that a surface is really just one scalar function of two angles: even though some codes represent a surface using two functions, namely R and Z, this introduces a tangential null space, which is removed by exploiting spectral condensation (Hirshman & Meier 1985;Lee et al. 1988; † We could have used the ideal MHD functional instead of the MRxMHD functional. In fact, we believe that this would have simplified some of the algebra; however, we chose MRxMHD because it can treat rational surfaces better than ideal MHD. Hirshman & Breslau 1998). Also, in (iv) we have restricted attention to when the external field is provided by a finite number of closed current filaments; other representations can be included.
For the quasi free-boundary approach, the theory of "efficient fields" can be used when D n is the independent degree-of-freedom (Landreman & Boozer 2016) to determine which Fourier spectrum of the normal component of the magnetic field corresponds to distant coils. Using only efficient fields, the optimization space can be effectively reduced. For the free-boundary (coil) algorithm, we expect that a concise parameterization of the coils can be introduced. For the other approaches, we suggest to use some measure of the coil error in total cost function: the quadratic-flux error would be a good target, and targets weighting resonant normal field errors are probably better.
In each presented case, the ∂ 2 F /∂x∂x matrix needs to be inverted to calculate ∂x/∂z. This is the equilibrium stability matrix, and near marginal stability its eigenvalues vanish. This suggests that the optimization can encounter singular points. Physically, this has to do with the fact that if a plasma is only marginally stable the equilibrium can be change substantially by even a small perturbation. This presumably would induce a large change in the coil geometry, whose optimal shape is therefore very uncertain. Marginal stability is a particularly relevant and important case. Toroidal confinement of fusion plasmas typically improves as the plasma pressure and or currents increase, but pressure and currents ultimately are responsible for instabilities. An economically viable fusion reactor would presumably operate at the highest fusion-performance parameters but also sufficiently far from bifurcation boundaries in parameter space so that small variations could be controlled.
Similarly, in all but the last case, the ∂ 2 E/∂c∂c matrix needs to be inverted to calculate ∂c/∂z. Small changes in the equilibrium may result in large changes to the coils if this matrix has a small eigenvalue.
Algorithms for combined plasma-coil optimization might become problematic near marginal stability boundaries, particularly if the response in either the plasma or the coils to small variations is not well understood. Depending on the plasma and coil properties chosen, there might be bifurcation boundaries in the optimization parameter space, which determine which initial conditions will converge to which local minima.
Marginal plasma stability is a lower-dimensional subspace of the full space of configurations, and it can be defined when there is a direction in which there is no variation, i.e. there are eigenvalues that are zeros. A similarly defined marginal stability subspace exists in the coil space. The stability boundary of the combined plasma-coil optimization, in the parameter space denoted by z, will be a combination of both. Depending on the plasma, P(x), and coil, C(c), properties chosen, there may exist subspaces for which ∂P/∂x = 0 and ∂C/∂c = 0, and there may be additional bifurcation boundaries in T (z), depending on what T (z) is chosen to be. For robust multi-objective functional optimization, it would be advantageous to understand separately the stability boundaries in the plasma, coil and optimization spaces for a variety of relevant plasma and coil properties.
In this paper, we focused on the ∂x/∂z and ∂c/∂z part of Eqn. 5.1, because these derivatives appear independently of which properties one wants to tailor in the optimization. We neglected the actual cost-function with the equilibrium, P, and coil, C, properties and their derivatives because these quantities are problem specific; for example, a reactor-grade plasma has different design criteria than a research experiment. Recently, much work has been devoted to the development of adjoint methods and automatic differentiation to calculate derivatives of the plasma and coil properties, e.g. for the rotational transform, neoclassical transport and for the volume averaged β For both the equilibrium calculation and the coil calculation, both descent-and Newton-style methods can be used. Descent methods can be employed for combined plasma-coil optimization. An interesting question is whether or not all these descent calculations can be performed simultaneously; for example, where τ is an arbitrary integration parameter, and α and β can be chosen for numerical stability so that, for example, the "inner" equilibrium and coil geometry calculations are sufficiently converged so that the "outer" optimization calculation is provided with sufficiently accurate information. The equations in Eqn. 6.1 can be thought of as a dynamical system; and, like most dynamical systems, the location and stability of the fixed points of Eqn. 6.1 gives important information about the "dynamics" (Strogatz 2018), which in this case is the convergence and stability properties of the combined plasma-coil optimization.

Acknowledgements
The lead author would like to thank J. Loizu, E. Paul, and B. Shanahan for helpful conversations. This work has been carried out in the framework of the EUROfusion Consortium and has received funding from the Euratom research and training programme 2014-2018 and 2019-2020 under grant agreement no. 633053. It was also supported by a grant from the Simons Foundation (560651, PH). The views and opinions expressed herein do not necessarily reflect those of the European Commission.