Representing the boundary of stellarator plasmas

In stellarator optimization studies, the boundary of the plasma is usually described by Fourier series that are not unique: several sets of Fourier coefficients describe approximately the same boundary shape. A simple method for eliminating this arbitrariness is proposed and shown to work well in practice.


Introduction
In optimized stellarators, the magnetic field lines usually trace out simply nested flux surfaces. Large magnetic islands or regions with chaotic field lines are avoided, at least in the plasma core, in the interest of good confinement. Kruskal and Kulsrud have shown that magnetostatic equilibria with this property (insofar as they exist) are uniquely determined by the shape of the toroidal boundary and by the plasma current and pressure profiles (Kruskal & Kulsrud 1958;Helander 2014). Instead of the current profile, that of the rotational transform can also be prescribed.
This fundamental result provides the theoretical basis for fixed-boundary magnetohydrodynamic (MHD) equilibrium calculations, which are commonly used in stellarator optimization studies. The shape of the plasma boundary is prescribed, usually as a Fourier series in poloidal and toroidal angles (Hirshman & Whitson 1983;Nührenberg & Zille 1988) R(θ, ϕ) = to a search in a space of 2(M + 1)(2N + 1) dimensions. † However, as has sometimes been remarked (Hirshman & Breslau 1998;Lee et al. 1988), this representation is not unique but contains "tangential degrees of freedom" in the limit M → ∞, N → ∞.
If a large but finite number of terms are included in the Fourier series, several very different choices of the coefficients {R mn , Z mn } correspond to approximately the same surface shape. Unless this problem is addressed, the search is therefore performed in a space of unnecessarily large dimensionality. Hirshman and co-workers devised a method called "spectral condensation" to deal with this problem, which is used internally in the VMEC and SPEC equilibrium codes to minimize the number of coefficients in the Fourier representation of all magnetic surfaces, including interior ones (Hirshman & Whitson 1983;Hirshman & Van Rij 1986;Hudson et al. 2012). However, spectral condensation is rarely used for the plasma boundary in optimization studies. The present article suggests another method of dealing with the problem of non-uniqueness of the boundary representation. This method is simpler but mathematically less sophisticated than spectral condensation. Unlike the latter, it does not correspond to a representation that is optimally economical, but it is simpler to implement numerically, requires less computation, and appears to work quite well in practice. The remainder of the present paper first describes the non-uniqueness of the representation (1.1) and how it can be eliminated, followed by examples showing how this technique simplifies the problem of optimization by eliminating a plethora of spurious and approximate minima of the target function in configuration space.
The fact that the addition of the function (θ, ϕ) to the poloidal angle θ does not change S indicates great freedom in the parameterization of the surface. If M = N = ∞ in the sum (1.1), then infinitely many choices of coefficients {R mn , Z mn } generate the same surface. Note that very different sets of coefficients can describe the same surface. If, on the other hand, M and N are finite in Eq. (1.1), so that the Fourier series of the functions R(θ, ϕ) and Z(θ, ϕ) terminate after a finite number of terms, then the corresponding series forR(θ, ϕ) andZ(θ, ϕ) will in general not terminate. † † Usually, the number is in fact slightly smaller, since negative values of n are not included in the terms with m = 0.
† For simplicity, we take to satisfy (−θ, −ϕ) = − (θ, ϕ) in order to preserve stellarator symmetry in the series (1.1). This observation has implications for the nature of the representation (1.1) when M and N are fixed, finite numbers: (i) If M and N are not too large, so that only a few terms are kept in the series (1.1), then the surface it represents will usually correspond to a unique set of coefficients {R m,n , Z m,n } or a small number of such sets. Of course, with only a few harmonics, not every surface can be represented by Eq. (1.1).
(ii) On the other hand, if many terms are included in the sum, M ∼ N 1, so that almost any stellarator-symmetric surface can be described the series Eq. (1.1), then many widely different choices of coefficients {R m,n , Z m,n } can correspond to almost the same surface. The more harmonics that are allowed in the sum, the less unique the representation becomes. These observations are confirmed by practical experience. At the beginning of a stellarator optimization run, it is usually futile to include many Fourier harmonics in the representation of the plasma boundary; the optimziation then "gets stuck" and does not proceed far from the initial state. Instead, it often proves useful to begin with only a few harmonics and gradually add more terms as the optimization proceeds, in order to allow for greater freedom in the shape of the plasma.

Spectral condensation
Spectral condensation exploits the non-uniqueness of the poloidal angle by minimizing the "spectral width", which measures the spectral extent of R mn and Z mn , under the constraint of not changing the geometry of the surface S. The spectral width is defined by Hirshman & Meier (1985) and Hirshman & Breslau (1998) as where p 0 and q > 0 are constants. † To first order in the perturbation introduced in the previous section,R(θ, ϕ) = R(θ, ϕ) + δR(θ, ϕ),Z(θ, ϕ) = Z(θ, ϕ) + δZ(θ, ϕ), and the Fourier coefficients of δR and δZ are given by The first-order variation in the spectral width thus becomes m,n m p (R 2 mn + Z 2 mn ), X(θ, ϕ) = m,n m p (m q − M )R m,n cos(mθ − nϕ), † In the SPEC code (Hudson et al. 2012) it is defined without the normalization, i.e.
Note that the m = 0 terms do not contribute to the Fourier sums, and that the spectral width assumes its minimum when I(θ, ϕ) = 0. This constraint is imposed to the requisite accuracy by Fourier expanding I(θ, ϕ) and requiring a number m * of the Fourier coefficients I mn to vanish, thus effectively removing m * degrees of freedom from the representation (Hirshman & Meier 1985).

An explicit boundary representation
The method of spectral condensation is optimal in the sense that it minimizes the spectral width of the representation, but it adds conceptual and computational complexity. The number of coefficients in the representation (1.1) remains high, although the constraints I mn = 0 effectively restrict the search to a submanifold of lower dimensionality, and the system of equations corresponding to these constraints must in general be solved numerically. In this and the next section we explore a simpler and more explicit construction as a possible alternative.
There are, of course, infinitely many ways of making the choice of poloidal angle unique, some of which have been proposed before, see e.g. Hirshman & Breslau (1998) and Carlton-Jones et al. (2020). A particularly simple choice could be to express the vertical coordinate as If the functions a and b are Fourier decomposed, one finds that this representation is of the same form as Eq. (1.1) but with only two poloidal harmonics, which are equal to In each poloidal cross section of the plasma surface, the poloidal angle θ thus defined is the polar angle of the horizontal projection on a circle with a diameter equal to the vertical extent of the surface, see Fig. 1. This choice of representation, which removes the superfluous degrees of freedom, can produce all surface shapes without multiple minima and maxima in the vertical coordinate Z in each poloidal cross section. The vast majority of all stellarators considered to date possess this property. However, Eq. (4.2) suffers from another and more serious shortcoming: it cannot economically represent a classical stellarator. Such devices have an elliptical poloidal cross section that rotates co-or counter-clockwise with increasing toroidal angle ϕ. This can be seen by introducing a rotating coordinate system, where α is a constant determining the rate of rotation. For a classical stellarator, it is equal to half the number of toroidal periods of the device, α = N/2, so that the cross section rotates by 180 degrees in one period. In these coordinates, a surface with rotating elliptical boundary is represented by where A and B denote the semi-axes. In our original coordinates, we obtain Hence it is clear that Z 1,−1 = Z 1,1 in contradiction to Eq. (4.3), which can only mean that the poloidal coordinate θ used in Eq. (4.2) cannot coincide with the corresponding one in Eq. (4.4). Although the former representation can describe any stellarator-symmetric surface, it needs many harmonics R m,n for a surface with rotating elliptical cross section. Close to the magnetic axis, most stellarators have this property, making this shortcoming serious indeed. Fortunately, it is easily overcome by applying a representation similar to Eq. (4.1) in the rotating coordinate system. This leads to the prescription In our experience, and as we shall see in the next section, this simple prescription works very well in practice.

Numerical examples
In this section, we explore a few examples of increasing complexity and realism, comparing our recipe with the conventional approach.

Simple axisymmetric (2D) case
Our first aim is to gain insight into the optimization space using the original arbitraryangle representation (1.1). Since only the poloidal angle θ is arbitrary, this issue can be explored in a simpler, two-dimensional setting. We choose to target an axisymmetric torus with a unit circle as the poloidal cross-section. A simple penalty function Q that is minimized by this surface is where ds = R (R 2 θ + Z 2 θ ) dθ dϕ, Here subscripts indicate partial derivatives, R θ = ∂R/∂θ, the surface area is A := ds, and the major radius is arbitrarily chosen to be R 0 = 1.5. We restrict R and Z to be axisymmetric: R mn = Z mn = 0 for all n = 0 and write R = 1.5 + m=1 R m cos(mθ), and Z = m=1 Z m sin(mθ).
If Q is not normalized to the area A, an artificial minimum exists when the area becomes small. With the normalization, the penalty function Q approaches unity in the limit of vanishing R − R 0 and Z (and thus vanishing surface area). The penalty function attains its sole minimum (Q = 0) if (R − R 0 ) 2 + Z 2 = 1. This equation is satisfied by many different choices of R m and Z m . † If all but the m = 1 Fourier coefficients vanish, i.e. if Q = Q(R 1 , Z 1 , R i = 0, Z i = 0) for i > 0, the penalty function attains the global minimum for R 1 = Z 1 = 1, see Fig. 3, and this is in general the case when only one pair of coefficients (R m , Z m ) is allowed to be non-zero, see Fig. 4. More interesting and complex behaviour is observed if Fourier harmonics with several † Although this is the only global minimum, Q becomes arbitrarily small for bounded surfaces having very large area, e.g., for highly "wrinkled" surfaces.  values of m are admitted, but one then faces the problem of graphically displaying the function Q of more than two variables. For instance, it is not easy to visualize how the cost function Q(R 1 , R 2 , Z 1 , Z 2 ) depends on all four arguments. However, some insight can be gained by plotting the mininum of Q with respect to two of the arguments as a function of the two other ones, e.g. by considering the functioñ Considering four Fourier harmonics in this way reveals the existence of three local minima, see Fig. 5. Two of these correspond to the global minimum, R 1 = Z 1 = 1, R 2 = Z 2 = 0 and R 1 = Z 1 = 0, R 2 = Z 2 = 1, and both correspond to the target surface, an axisymmetric torus with a unit circle cross section, parameterized in two different ways. The third minimum is located at R 1 = 0.0, R 2 ≈ 1.05, Z 1 = 1.15, and Z 2 = 0.0 with Q ≈ 0.133. It corresponds to a surface with zero volume but finite area, see Fig. 6. Increasing the number of poloidal harmonics m to three makes it more difficult to  : The three minima of Q(R 1 , R 2 , Z 1 , Z 2 ). In blue (global minima): axisymmetric torus with unit circle cross section described by R 1 = Z 1 = 1, R 2 = Z 2 = 0 and R 1 = Z 1 = 0, R 2 = Z 2 = 1 (R i = 0 for all i > 2). In gray (local minimum): R 1 = 0.0, R 2 ≈ 1.05, Z 1 = 1.15, and Z 2 = 0.0. locate the local minima of the cost function. Without a global optimizer, one encounters many local minima depending on the initial values chosen for the remaining Fourier coefficients. Using differential evolution, a global optimization routine, to obtain min R2,Z2,R3,Z3 Q(R 1 , R 2 , R 3 , Z 1 , Z 2 , Z 3 ) one finds a landscape broadly similar to the one found for the four-Fourier-coefficient case, but with many additional small local maxima and minima, see Fig. 7. This type of scan has to be considered with caution. Most global optimization routines do not, in practice, guarantee a global minimum but sometimes end up in local ones. In local optimization routines, this problem is of course still more acute, since the outcome generally depends on the initialization. To visualize the difficulty of finding global minima, it is useful to fix two coefficients, in the following R 1 (= 0.18) and Z 1 (= 0.4), and study how the landscape depends on the remaining ones. We note that the function min R3,Z3 Q(R 1 = 0.18, R 2 , R 3 , Z 1 = 0.4, Z 2 , Z 3 ) possesses five local minima with R, 2 and Z 2 in the range 0.0 − 1.0, see Fig. 8 and line discontinuity, as can be seen in Fig. 8. To understand the discontinuity min R3,Z3 Q(R 1 = 0.18, R 2 , R 3 , Z 1 = 0.4, Z 2 , Z 3 ), we plot Q(R 1 = 0.18, R 2 = x, R 3 , Z 1 = 0.4, Z 2 = y, Z 3 ) the function R 3 and Z 3 for selected values for x and y, Fig. 9. The number of local minima varies with x and y. For (x, y) = (0.24, 0.39) there are four local minima in the figure, for (x, y) = (0.5, 0.5) there are two of them, and for (x, y) = (1, 1) there is only one minimum. It is thus clear that a local optimizer that seeks local minima of the function Q(R 1 = 0.18, R 2 , R 3 , Z 1 = 0.4, Z 2 , Z 3 ) will find different ones depending on the starting point for R 3 and Z 3 . Abrupt changes (discontinuity) in min R3,Z3 Q(R 1 = 0.18, R 2 , R 3 , Z 1 = 0.4, Z 2 , Z 3 ) appear when a local minimum disappears and the optimizer finds a different one. The representation proposed in Sect. 4 leads to much more benign results when applied to the model problem (5.1). Restricting the optimization space to axisymmetric designs leads to R = 1.5 + m=1 R m cos(mθ), The minimized cost function min R2,Z2,R3,Z3 Q(R 1 , R 2 , R 3 , Z 1 , Z 2 , Z 3 ) with respect to R 1 and Z 1 . Differential Evolution, a global optimization routine, was used to find the minima. This time, we find that the landscape of the minimum penalty function min Q does not change much when the number of Fourier harmonics of R is increased. In the case of four Fourier harmonics, Q(R 1 , R 2 , R 3 , Z 1 ), there is one global minimum at R 1 = 1, R 2 = 0, R 3 = 0, Z 1 = 1 and a second shallow local minimum near R 1 = 0 and Z 1 = 1.0, see Fig. 10. This local minimum disappears when more Fourier harmonics are added. Importantly, the outcome is similar whether a local and global optimization algorithm is employed, see Fig. 11a and Fig. 11b, making it much easier to find the minima numerically.  : Q(R 1 = 0.18, R 2 = x, R 3 , Z 1 = 0.4, Z 2 = y, Z 3 ) with respect to R 3 and Z 3 . Figure 10: The minimized cost function min R2,R3 Q(R 1 , R 2 , R 3 , Z 1 ) with respect to R 1 and Z 1 .
(a) using a non-global optimization algorithm. (b) using Differential Evolution -a global optimization routine. Figure 11: min R2,R3,R4,R5 Q(R 1 , R 2 , R 3 , R 4 , R 5 , Z 1 with respect to R 1 and Z 1 , (a) D shape (b) bean shape Figure 12: D shape and bean shape reproduced based on Hirshman & Meier (1985) where the solution of our boundary representation overlaps with the original boundary.

Fourier representation of stellarators
We now turn to examples of explicit choices of the coefficients R 0n , Z 0n , ρ m,n and b n , corresponding to stellarator plasma boundaries that have been explored in this context in the past. We begin with examples from Hirshman & Meier (1985), who analysed shapes using spectral condensation, thus providing a convenient point of comparison with this technique. We start with a planar D-shaped boundary given by R = −0.23+0.989 cos θ+0.137 cos 2θ and Z = 1.41 sin θ − 0.109 sin 2θ after spectral condensation (Hirshman & Meier 1985). We use Fourier decomposition to obtain the coefficients in our unique boundary representation that reproduce this boundary, restricting the number of modes m to be such that all the coefficients exceed 0.01. The result is R 00 = −0.306, b 0 = 1.426 and ρ = 0.957 cos θ + 0.207 cos 2θ + 0.032 cos 3θ. The error compared to the original boundary is ≈ 0.4% although the same number of Fourier harmonics are used as in the spectral condensation technique. Hirshman and Meier also considered a bean-shaped surface given by R = −0.320 + 1.115 cos θ+0.383 cos 2θ−0.0912 cos 3θ+0.0358 cos 4θ−0.0164 cos 5θ and Z = 1.408 sin θ+ 0.154 sin 2θ−0.0264 sin 3θ after spectral condensation. Applying our representation to this case, we obtain R 00 = −0.184, b 0 = 1.419 and ρ = 1.184 cos θ+0.208 cos 2θ−0.143 cos 3θ+ 0.064 cos 4θ−0.032 cos 5θ+0.011 cos 6θ with an error of ≈ 0.15%. Our representation thus needs even fewer Fourier harmonics than this spectral condensation. † Finally, we consider a representative example from Wendelstein 7-X, where the magnetic field in the so-called standard configuration was calculated using an equilibrium solver in † Hirsman and Meier also consider a third case, a so-called belt pinch boundary, which cannot be reproduced by our boundary representation since it has multiple minima and maxima in the vertical coordinate. free-boundary mode. As shown in Fig. 13, the resulting plasma boundary can be faithfully reproduced with mode numbers m 5 and |n| 3. Thus, the representation proposed in Sec. 4 can accurately and economically reproduce relevant plasma boundary shapes, including Wendelstein 7-X and other cases studied earlier in the literature. It does not always need as few Fourier harmonics as spectral condensation, but for "reasonable" shapes it appears comparable in efficiency and avoids the need for computational optimization, which is an integral part of the spectral condensation technique.

Application to 3D stellarator optimization
Finally, we put our boundary representation to the test in a real stellarator optimization problem, where the plasma boundary serves as input for a fixed-boundary equilibrium calculation and is adjusted iteratively until a target function reflecting plasma performance has been minimized. As described in the introduction, such optimization calculations have in the past usually been performed with the ambiguous boundary representation (1.1). We start the optimization with a rotating elliptical boundary, Fig. 14. and use the optimization code ROSE (Drevlak et al. 2019) with the equilibrium code VMEC and a non-gradient, non-global optimization algorithm (Brent). The target rotational transform is chosen to be 0.25 on axis and 0.35 at the plasma boundary, and in addition we require the magnetic well to exceed a certain threshold (0.1) and the toroidal projection of the plasma boundary to be convex in every point. As usual in this type of optimization, the target function f is a weighted sum of squares, where F i is the value for criterion i,F i the corresponding target value, and w i the i'th weight, which can be adjusted to obtain various different optimal (Pareto) points. Of course, the performance of the optimization depends of on the exact choice of w i and F i as well as the initial condition, but we find that the results turn out much better, and more quickly, with the novel representation than with the standard one used in VMEC.
With the same weights chosen for both cases, we typically obtain a penalty value f that is two orders of magnitude smaller when the new boundary representation is employed. The resulting configuration is thus significantly different, and much better, than that obtained with the conventional method. An example of how the plasma boundaries differ is shown in Fig. 15. In this example, all the aims of the optimization were attained when the novel scheme was used, whereas the usual one did not succeed in achieving the prescribed rotational transform and a non-concave plasma boundary. Similar results have also been found with other, more complicated, optimization targets.

Conclusions
In summary, the usual Fourier-series representation of the plasma boundary used in stellarator optimization contains much redundancy due to the arbitrariness of the poloidal angle. This redundancy grows exponentially with the number of terms in the series and unnecessarily increases the dimensionality of the search. It causes a plethora of local  minima to appear in the optimization landscape, as can be illustrated with simple 2D examples. The situation can be remedied by making the poloidal angle unique, but some care must be taken to ensure that simple stellarator shapes can still be represented in an economical way. When this is done, local minima are eliminated and the optimization proceeds more rapidly than with the usual representation. The outcome also tends to be better, especially if a non-global optimization algorithm is used. Our specific boundary parametrization (4.5)-(4.6) is simple and intuitive, and requires less computation than the spectral condensation method, but it cannot describe stellarator boundaries with multiple maxima in the ζ-direction. Such boundaries are however highly unusual, and the representation can easily be generalized to include such shapes.

Acknowledgement
The primary author would like to thank G. Plunk 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.