Stellarator optimization for good magnetic surfaces at the same time as quasisymmetry

The geometry of a stellarator can be optimized to improve confinement and stability. Many of the figures of merit that are commonly included in stellarator optimization, such as neoclassical transport and magnetohydrodynamic (MHD) stability, are most conveniently calculated under the assumption that nested toroidal magnetic surfaces exist. Moreover, most of the physics codes in the stellarator community use the data format of the Variational Moments Equilibrium Code (VMEC),¹ in which the existence of nested magnetic surfaces is assumed. However, magnetic surfaces are not guaranteed to exist in stellarators, as magnetic islands and chaotic field regions may be present. As confinement across islands and chaotic regions is severely diminished, avoiding or controlling these phenomena is a high priority in stellarator experiments^2–4 The extent of islands and chaos is another function of the magnetic geometry that can be optimized.^5,6 Islands and chaos can be diagnosed using other magnetic field representations or MHD equilibrium codes, but then physics objectives that presume the existence of surfaces are not straightforward to calculate. In stellarator optimization, it would therefore be ideal to have the best of both worlds, taking advantage of codes that assume the existence of surfaces and use the VMEC data format, while at the same time controlling islands. In this paper, we demonstrate one approach to achieving these goals.

The approach here involves simultaneously using two magnetic field representations, one that assumes the existence of surfaces and one that does not. In particular, we use VMEC and the Stepped Pressure Equilibrium Code (SPEC),^7,8 which both compute three-dimensional MHD equilibria. VMEC has been widely used in the stellarator community since the 1980s, so a large number of transport and stability codes are available that analyze the numerical solutions it produces. VMEC operates by minimizing the MHD energy subject to two constrained radial profiles and the constraint of nested magnetic surfaces. Hence, magnetic islands and chaos cannot be represented. SPEC is a newer code in which the toroidal domain is divided into a number of nested annular regions, with the pressure constant in each region. The magnetic field is constrained to be tangent to the toroidal boundary surfaces of each region, and pressure balance is enforced across these interfaces. Within each of the regions, there is no constraint that magnetic surfaces exist, and so islands and chaos can be represented. VMEC and SPEC are both applicable to equilibria with or without plasma pressure and current; in this work, vacuum fields will be considered to expedite the SPEC calculations. While VMEC has been used as part of stellarator optimization for decades, the present work is the first time SPEC has been used in optimization.

In the new optimization approach presented here, both VMEC and SPEC are run at each evaluation of the objective function. If any islands are present in the SPEC solution, SPEC and VMEC necessarily will not agree exactly on the magnetic field. A measure of the magnetic island width is computed from the SPEC solution and included as a penalty in the objective function. Motivated by Refs. 5 and 6, the measure we use is the residue,⁹ a real number that can be computed for any periodic field line, which is zero when the island width vanishes. Due to this island width penalty, the islands are eliminated during the optimization so that the VMEC and SPEC representations agree by the end. At the same time, the objective function also includes quantities derived from the VMEC solution. In particular, here we minimize the deviation from quasisymmetry. Quasisymmetry is a continuous symmetry in the field strength that ensures guiding-center confinement.¹⁰ The quasisymmetry objective is based on a conversion of the VMEC solution to Boozer coordinates.¹¹ At the start of the optimization, this VMEC-derived part of the objective function may be somewhat inaccurate because the islands were ignored in its computation. However, by the end of the optimization, since the islands have been eliminated by the residue penalty, the quasisymmetry term is accurate. In this way, we can take advantage of calculations that are based on the existence of magnetic surfaces and of the many available codes that postprocess VMEC equilibrium files, while fully accounting for the possibility of magnetic islands.

Compared to earlier island healing calculations done for the National Compact Stellarator Experiment (NCSX) project,^12–15 there are several differences in the approach here, such as use of different equilibrium codes and physics objectives and different measures of island width. Whereas we optimize island widths and other quantities concurrently, the NCSX approach did not use optimization; instead a nonlinear system of equations was solved in which quantities other than island width (such as MHD instability growth rates) were held fixed.¹⁵ 

Here, we will consider the first stage of the two-stage optimization procedure used for the design of experiments such as W7-X¹⁶ and HSX.¹⁷ In this first stage, the parameter space for optimization is the shape of a toroidal boundary magnetic surface. Specifically, the parameter space consists of the Fourier modes ${Rm,n, Zm,n}$ of the boundary toroidal surface

where $ϕ$ is the standard toroidal angle, θ is any poloidal angle, $nfp$ is the number of field periods, and we have assumed stellarator symmetry. We exclude the major radius $R0,0$ from the parameter space, in order to fix the spatial scale. Here, we will not consider the second stage of the two-stage approach, in which coils are optimized to produce the boundary surface resulting from the first stage. In the future, the method of this paper could also be used in a single-stage optimization, in which the parameter space consists of coil shapes, and free-boundary equilibria are used.

For simplicity, we will consider a configuration with no plasma current or pressure. This choice minimizes the computational cost because a single radial domain can be used in SPEC. In the future, the procedure here could be applied to configurations with nonzero plasma current and pressure, using multiple radial domains in SPEC.

The numerical example is carried out using SIMSOPT, a new software framework for stellarator optimization.^18,19 The optimization is driven in python, using the default algorithm (trust region reflective) for nonlinear least squares optimization from the scipy package.²⁰ Gradients for the minimization are calculated with forward finite differences, using MPI for concurrent function evaluations.

The initial state is an axisymmetric circular-cross section torus, and the number of field periods $nfp$ is set at two. We first carry out a preliminary optimization without SPEC or residues. The reason for this preliminary optimization is that the ι profile evolves significantly at the beginning, causing resonances to enter and leave the domain. In this case, it is awkward to include residues in the objective function. If a particular rational ι is outside the computational domain, the corresponding residues cannot be calculated. If the residues in the objective are set to zero in this case, the objective would have a discontinuous change whenever the rational surface entered or left the computational domain during the optimization, which will likely cause poor performance for many optimization algorithms. Therefore, no residues are included in the objective function for the preliminary optimization, which is

Here, A is the effective aspect ratio as defined in VMEC, ι₀ and ι_a are the rotational transform at the magnetic axis and edge, $Bm,n$ is the amplitude of the $cos (mϑ−nφ)$ Fourier mode of the field strength in Boozer coordinates $(ϑ,φ)$ for the flux surface with normalized toroidal flux s = 0.5, and only $n≠0$ modes are included in the sum. All terms in the objective are computed from the VMEC solution, with the $Bm,n$ values computed by postprocessing of the VMEC solution with the BOOZ_XFORM code.²¹ The aspect ratio is included in the objective because if not, the quasisymmetry term can be reduced to zero by increasing the aspect ratio to infinity. The rotational transform terms are included in the objective because if there are no constraints on ι, true axisymmetry is an optimum. The factor of 2 in the quasisymmetry term of (2) is chosen based on experience to give the best optimum. Using a script similar to the one described shortly for the combined VMEC-SPEC optimization, the optimization is performed in a series of three steps, as the size of the parameter space and resolution parameters are increased.

For the combined VMEC-SPEC optimization, the objective function is

where R₁ and R₂ are the residues⁹ for two periodic field lines on the rational surface. These two curves are located using Newton's method by searching for field lines with the desired periodicity to the right and left of the magnetic axis along the line Z = 0 starting at the $ϕ=0$ plane. When an island chain is present, R₁ may correspond to the O-point and R₂ to the X-point, or vice versa. A value of the residue in the interval $(0, 1)$ indicates that the periodic field line is an O-point, whereas if a residue is negative or > 1, nearby field lines diverge from it, as in an X-point. Field lines on a good flux surface with rational ι have a residue of zero. Since residues are dimensionless and of order unity, a residue can be considered small if its absolute value is $≪1$⁠.

The residues in Eq. (3) are computed from the SPEC solution by integrating along the periodic field lines to compute the tangent map. The weight factor of 2 in the residue terms of (3) is chosen by experimentation to yield a good optimum.

We expect that the addition of the residue terms in Eq. (3) should have relatively little impact on minimization of the other objective terms. This expectation is based in part on Refs. 12–15, where it was found that the condition of zero island width could be imposed while leaving a high dimensional null space that could be exploited for imposing other objectives or constraints. Moreover, in other numerical experiments in which only squared residues are minimized, we find the optimum is not unique: many different optima without islands can be obtained, depending on which boundary Fourier modes are allowed to vary in the optimization. Therefore, the condition that the residues vanish appears to be less constraining than the condition of quasisymmetry.

The SIMSOPT python driver script to define and solve the minimization problem is shown in Fig. 1, as several features are noteworthy. On line 14, the Spec object is configured to use the same boundary Surface object as the Vmec instance. Therefore, when the shape of this single surface is modified during the optimization, the outputs of both VMEC and SPEC change accordingly. The objective function (3) is specified in lines 16–34. Also, since the optimization problem is defined with a script, any other desired scripting elements can be included. Here, this capability is used to define a series of three optimization stages, in which the size of the parameter space (the maximum m and n values of the ${Rm,n, Zm,n}$ to vary) is increased at each step, along with the numerical resolution parameters of the codes. The former is valuable to avoid getting stuck in a poor local minimum, and the latter improves computational efficiency. In lines 36–41, it can be seen that for the first step, the boundary amplitudes ${Rm,n, Zm,n}$ are varied for $m=0…3$ and $n=−3⋯3$⁠. In the second step, the maximum m and $|n|$ are increased to 4, and in the third step, the maximum m and $|n|$ are increased to 5. (For the preliminary optimization, the corresponding maximum mode numbers were 1, 2, and 3.)

The configurations before and after the final optimization are shown in Figs. 2–4. In the figures, the configuration resulting from the preliminary optimization used as input for the combined VMEC-SPEC optimization is labeled “Before optimization.” Figure 3 shows that this initial configuration has a significant island chain at the $ι=2/5$ resonance. Therefore, VMEC and SPEC do not agree on the internal flux surface shapes near the islands for this configuration. It can also be seen in Fig. 3 that the optimization has successfully eliminated the islands. Indeed, the magnitudes of the residues have been reduced from $2×10−3$ to $2×10−6$⁠. The two codes agree very well on the internal surface shapes by the end of the optimization. Therefore, calculations for the final configuration based on the VMEC solution, such as the Boozer coordinate transformation, can be trusted.

Figure 5 displays the rotational transform profiles at the beginning and end of the optimization. For this figure, ι was computed by following field lines in the SPEC solution, starting from an array of points on the inboard midplane between the magnetic axis and computational boundary. A flat region in this ι profile for the initial configuration reflects the presence of a magnetic island. The ι profile ranges from 0.39 to 0.42 for both the initial and optimized configurations, so the islands were not eliminated by shifting the resonance out of the domain, but rather by tuning of the resonant field.

The final configuration also has extremely good quasiaxisymmetry, especially on the s = 0.5 surface where symmetry was optimized. This can be seen in the straight horizontal contours of $|B|$ in Fig. 6. Any deviation from symmetry is not perceptible in the figure. By contrast, analogous figures of $|B|$ on a surface for previously published quasisymmetric configurations have almost always shown visible ripples in the $|B|$ contours. Examples include Figs. 5 and 6 of Ref. 23, Fig. 2 of Ref. 24, and Fig. 4 in Ref. 25. The only previously published configurations we are aware of without clear curvature of the $|B|$ contours are those of Fig. 21 in Ref. 26, which are a much higher aspect ratio (⁠$≥78$⁠).

The quality of quasisymmetry in the optimized configuration can also be seen in Fig. 7, which shows the dependence of the $|Bm,n|$ amplitudes on minor radius. At s = 0.5 where symmetry was optimized, there is a striking notch in the $n≠0$ modes where their amplitude becomes extremely small, $<10−5$ of the mean field. This finding supports the conjecture that it may be possible to obtain quasisymmetry exactly on an isolated magnetic surface.²⁷ At the plasma edge, where symmetry-breaking $Bm,n$ modes are largest, the largest nonsymmetric mode remains smaller than all of the symmetric modes for $m=0–3$⁠.

In this work, we have demonstrated a method for optimizing a stellarator's geometry to eliminate magnetic islands while simultaneously optimizing other objectives that assume the existence of magnetic surfaces. This makes it possible for optimizations to include physics codes that use equilibria from the VMEC code or other equilibrium codes that assume the existence of nested magnetic surfaces, while also ensuring good surface quality.

Although islands are eliminated in the method here, rational surfaces remain in the plasma. Therefore, error fields have the potential to reintroduce any island chains that were eliminated. For an experiment, one should therefore consider the sensitivity of the configuration to coil position errors or other error fields. The most meaningful calculation of this sensitivity would involve computing coil shapes for the final configuration, and then determining the island size that results from various coil shape errors or the addition of other error fields.

The method of this paper can be applied in the future to configurations with nonuniform pressure by using multiple radial domains in SPEC. While quasi-axisymmetry was the main objective in the example here, the method is equally applicable to other objectives.

Another extension of this work could be to use a measure of flux surface quality or island width other than residues. In principle, any such method could be used in the procedure of this paper. Several such alternative measures include Mather's^28,29 $ΔW$⁠, the estimate by Cary and Hanson,³⁰ and converse KAM.^29,31

We gratefully acknowledge discussions with and assistance from Aaron Bader, David Bindel, Benjamin Faber, Andrew Giuliani, Stuart Hudson, Rogerio Jorge, Thomas Kruger, Zhisong Qu, Jonathan Schilling, Florian Wechsung, and Mike Zarnstorff. This work was supported by a grant from the Simons Foundation (560651, M.L.). B.M. and C.Z. are supported by the U.S. Department of Energy under Contract No. DE-AC02-09CH11466 through the Princeton Plasma Physics Laboratory.

The data that support the findings of this study are openly available in Zenodo at http://doi.org/10.5281/zenodo.5035515, Ref. 22.