Stellarator optimization for nested magnetic surfaces at finite β and toroidal current

In toroidal geometries, three-dimensional (3D) magnetohydrodynamic (MHD) equilibria are, in general, a mix of nested magnetic surfaces, magnetic islands, and magnetic field line chaos.^1–3 In the plasma core, the latter two topologies are usually detrimental to confinement, i.e., the radial transport of particles and energy is generally greater than in regions of nested magnetic surfaces. In addition to other desirable properties, a common target of stellarator optimization is to increase the volume occupied by magnetic surfaces.⁴ 

The equilibrium approach for optimizing the volume occupied by nested flux surfaces requires three tools: (i) a fast 3D equilibrium code that does not assume nested flux surfaces, (ii) a numerical diagnostic that provides a measure of integrability of the magnetic field,^5–8 and (iii) an optimization algorithm. Many 3D MHD equilibrium codes have been developed, for example, VMEC,^9,10 SIESTA,^11,12 HINT,^13,14 PIES,^15,16 DESC,¹⁷ or SPEC.¹⁸ Among these, DESC and VMEC assume nested magnetic surfaces everywhere in the plasma and are not suitable for describing equilibria with magnetic islands and field line chaos when used as stand-alone codes. The SIESTA, HINT, and PIES codes allow the generation of these field line topologies, but are comparatively slow and are as of yet not suitable as a main MHD solver in a stellarator optimization loop. On the other hand, SPEC can compute 3D equilibria with nested magnetic surfaces, magnetic islands, and field line chaos. It is based on an energy principle and has, thus, the potential to be as fast as VMEC, i.e., sufficiently fast to be considered as the main MHD solver of an optimization loop. SPEC, however, is constrained to solve equilibria with stepped-pressure profiles. In addition, SPEC has been recently extended to allow free-boundary calculations¹⁹ and to allow the prescription of net toroidal current profiles.²⁰ Numerical work has also improved its robustness and speed.²¹ 

Most physics codes developed by the stellarator community are based on the assumption of nested flux surfaces and thus require, for example, a VMEC equilibrium as input. In a recently published paper, Landreman et al.²² showed that by optimizing the plasma boundary geometry, SPEC can be used in combination with VMEC to obtain self-consistent vacuum configurations where both codes are in agreement, ensuring good magnetic surfaces in the region of interest. This allows us then to trust in any auxiliary codes that assume nested flux surfaces in this configuration and to safely optimize for other metrics.

In this paper, we extend the work by Landreman et al.²² by showing that finite-β SPEC equilibria with non-zero net toroidal currents can also be optimized to reduce the volume occupied by magnetic islands and field line chaos in a reasonable amount of time. Additionally, we explore the use of parameter spaces other than the plasma boundary geometry that could be of interest. Indeed, we leverage new capabilities of SPEC to show that the volume of magnetic surfaces in a stellarator can be maximized by optimizing the injected toroidal current profile or the coil configuration—two experimentally relevant “knobs.”²³ For the optimization, we follow Landreman et al.²² and use the SIMSOPT framework,²⁴ which, in particular, can construct an objective function based on Greene's residues⁷ of some selected rational surfaces.

Since optimizations of both fixed- and free-boundary equilibria are considered in this paper, we first describe how these are computed by SPEC. More details about the SPEC numerical algorithm can be found in Hudson et al.^18,19 and references therein.

SPEC computes equilibria where the pressure profile is stepped, i.e., a finite number of nested surfaces support pressure discontinuities. Between these interfaces, there is a finite volume where the pressure is constant, and the force-free magnetic field is described by a Taylor state.^25,26 These interfaces thus describe a set of N_vol nested annular regions, called SPEC volumes or simply volumes in what follows.

Using the standard cylindrical coordinate system $(R,ϕ,Z)$⁠, a plasma boundary surface $ΓPB$ is parametrized by $R(θ,ϕ), Z(θ,ϕ)$⁠, where θ is an as of yet undetermined poloidal angle. A fixed-boundary SPEC equilibrium is then determined by $ΓPB$⁠, the number of volumes N_vol and the toroidal flux they enclose ${ψl}l={1,…,Nvol}$⁠, the pressure in each volume ${pl}l={1,…,Nvol}$⁠, the net toroidal current flowing at the volumes' interfaces ${Iϕ,ls}l={1,…,Nvol−1}$⁠, and the net toroidal current flowing in each volume ${Iϕ,lv}l={1,…,Nvol}$⁠. Interface currents represent all equilibrium pressure-driven currents, such as diamagnetic, Pfirsch–Schlüter, and bootstrap currents, as well as δ-function current densities arising when an ideal interface is positioned on a resonance,²⁷ while volume currents represent externally driven currents, such as electron cyclotron current drive (ECCD), neutral beam current drive (NBCD), or Ohmic current.²⁰ Note that SPEC can be run with different constraints, e.g., the rotational transform—this will however not be covered in the present paper. As an output of a fixed-boundary calculation, SPEC gives the magnetic field in each plasma volume as well as the geometry of each interface.

A free-boundary equilibrium is determined by a computational boundary surface $ΓCB$ surrounding the plasma and enclosed by the coils, and the Fourier harmonics V_mn of the component of the vacuum magnetic field normal to $ΓCB$⁠, i.e.,

where $n=eθ×eϕ$ is a vector normal to $ΓCB$⁠, $Bc$ is the magnetic field produced by the coils, M_pol and N_tor are, respectively, the number of poloidal and toroidal Fourier modes used in the calculation, N_fp is the number of field periods, and stellarator symmetry has been assumed. The Fourier harmonics V_mn can be obtained from the coil geometry and coil currents; changing the coil shapes, coil currents, or the computational boundary $ΓCB$ will modify the V_mn harmonics non-linearly. In addition, the calculation of a free-boundary equilibrium requires specifying the total coil current flowing through the torus hole, and, as in fixed-boundary calculations, the profiles ${ψl}l={1,…,Nvol−1},{pl}l={1,…,Nvol−1}, {Iϕ,ls}l={1,…,Nvol−1}$⁠, and ${Iϕ,lv}l={1,…,Nvol−1}$⁠. The pressure and currents in the last region (⁠$l=Nvol$⁠) are set to zero, thereby emulating a vacuum region. As an output of a free-boundary calculation, SPEC gives the plasma boundary geometry $ΓPB$ (i.e., the geometry of the outer interface of volume $Nvol−1$⁠), the magnetic field in each plasma volume as well as in the vacuum region between $ΓPB$ and $ΓCB$⁠, the total toroidal flux enclosed by $ΓCB$⁠, denoted hereafter by $Ψ̂≡∑l=1Nvolψl$⁠, and the geometry of the interfaces.

We now consider the optimization of a finite-β, finite current equilibrium in a classical stellarator geometry (rotating ellipse) calculated with SPEC, which presents regions of magnetic islands and magnetic field line chaos. This model of a classical stellarator was chosen for simplicity, as few Fourier modes are required to described the equilibrium. An experimental instance of a classical stellarator was W7-A.²⁸ However, the optimization procedure presented here does not depend on the specific choice of geometry. A strongly shaped stellarator, such as W7-X, could also in principle be optimized with the presented approach.

We construct a free-boundary equilibrium, which will be the initial state for the optimizations presented in this paper. We choose the computational boundary $ΓCB$ to be a rotating ellipse,

where $R00CB=10$ m, $R10CB=−Z10=1$ m, $R11CB=Z11=0.25$ m, and with N_fp = 5. We notice that this choice of Fourier series to represent an ellipse determines uniquely the θ angle on the computational boundary. The Fourier resolution is set to $Mpol=Ntor=6$⁠. We impose a pressure profile linear in toroidal flux ψ that we approximate with seven steps, namely, N_vol = 8, with the last volume being a vacuum. We scale the pressure profile such that the plasma averaged $⟨β⟩$⁠, defined as the value of β averaged over the volume contained by the plasma boundary $ΓPB$⁠, is $⟨β⟩=1.5%$⁠. We assume no externally induced currents by setting the total toroidal current flowing in each volume to zero, ${Iϕ,lv}l={1,…,Nvol}=0$⁠, and we assume that the plasma generates a bootstrap-like toroidal current proportional to the pressure jump at the interface $Iϕ,ls∝(pl+1−pl)$ (see Fig. 3), which sums up to a net toroidal current enclosed by the plasma of 29.5 kA. Finally, we suppose the existence of coils, with a total current of 17.1MA, such that $B·n=0$ on $ΓCB$⁠, with $B=Bc+Bp$ and $Bp$ is the magnetic field produced by the plasma currents. The condition $B·n=0$ is only valid for the specified $ΓCB$⁠, $⟨β⟩$ and profiles. Note that the existence of a feasible coil system that generates such a field is not covered in this paper—only their normal projection on the computational boundary is specified here. Other codes, such as FOCUS,²⁹ can be used to obtain such coils.

We will refer to this equilibrium as the initial free-boundary equilibrium; its associated magnetic topology is shown via its Poincare section and rotational transform, plotted on the top left panel of Fig. 1 and on Fig. 2, respectively. Note that the toroidal flux $Ψ$ in volume l is defined via interpolation between $ψl−1$ and ψ_l. The discontinuities observed in the rotational transform profile are due to SPEC's stepped-pressure equilibrium model—since the magnetic field is generally discontinuous across the interfaces, and so is the rotational transform.

Free-boundary equilibria can be emulated by a fixed-boundary calculation if the toroidal current and the pressure in the last volume, $Iϕ,Nvolv$ and $pNvol$⁠, respectively, are set to zero. This corresponds to an equilibrium where a perfectly conducting wall $Γw$⁠, parametrized by $Rw(θ,ϕ)$ and $Zw(θ,ϕ)$⁠, is located at the last interface, $l=Nvol$⁠, and the plasma boundary $ΓPB$ is the interface $l=Nvol−1$⁠, which is allowed to move. The exact same equilibrium as the initial free-boundary equilibrium can thus be generated with fixed-boundary conditions if we set $Γw=ΓCB$⁠, given by Eqs. (2) and (3). The difference is, however, that for any value of $⟨β⟩$ and choice of profiles ${Iϕ,lv}, {Iϕ,ls}$⁠, the condition $B·n=0$ remains valid (perfectly conducting wall). In this case, the total toroidal flux enclosed by $Γw$ is imposed to be $Ψ̂$⁠. We will refer to this equilibrium as the perfectly conducting wall equilibrium. Note that the free-boundary equilibrium has effectively no conducting wall (no wall limit).

We now consider different degrees of freedom depending on the type of initial equilibrium. For the free-boundary equilibrium, the parameter space is a selected set of the Fourier harmonics V_mn, which emulate an optimization of the coil geometry and coil currents, as would be done in a “single-step stellarator optimization.”^4,30 Indeed, assuming that a coil system exists that generates any considered V_mn, optimizing the coils would be equivalent to optimizing the V_mn harmonics. For the perfectly conducting wall equilibrium, we consider two different parameter spaces. The first is the current flowing in the plasma volumes, ${Iϕ,lv}l={1,…,Nvol}$⁠, which is the externally induced net toroidal plasma current. This parameter space is relevant, for example, for scenario design, where one could desire to heal magnetic islands and chaos for a given plasma geometry and coil system. The second parameter space is the geometry of the wall expressed as Fourier series,

The degrees of freedom are then a selected set of Fourier harmonics ${Rmnw,Zmnw}$⁠.

The objective functions for each optimization are based on Greene's residues.⁷ The Greene's residue R is a quantity that can be computed for any periodic orbit, with R = 0 when the island width is zero, $0<R<1$ for an O-point, and R > 1 or R < 0 for an X-point. The objective function is

where R_i is the residue for a field line on the $ith$ targeted island chain, and x is the considered degrees of freedom. The downside of this approach is that each resonant field line has to be selected by hand before starting the optimization. If a new resonance becomes excited during the optimization, it will not be included in the objective function. From the initial equilibrium, Fig. 1 (top left), we can identify a set of resonant rational surfaces and their associated residue, listed in Table I. We use the SIMSOPT framework to drive the optimization, which is based on the default scipy.optimize³¹ python algorithm for non-linear least squares optimization.

We start by optimizing the initial free-boundary equilibrium. Residues with indices 1–3 in Table I are used to build the objective function according to Eq. (6). These residues correspond to the islands and subsequent chaos present in the vacuum region, right outside the plasma edge. This choice is somewhat arbitrary and mainly due to the fact that this particular equilibrium presents most of the chaos in that region. Indeed, the core confinement of particles and heat may not be directly affected by the topology of field lines in the vacuum region, but it is important to control the presence of magnetic islands or field line chaos in this region for the stellarator heat exhaust system. This is particularly important in stellarators with an island divertor, such as W7-X.²³ 

The parameter space is a selected set of Fourier modes ${Vmn}, (m,n)={(6,n)}n={−3,…,3}$⁠. In principle, more residues listed in Table I could be targeted if more ${Vmn}$ modes were used as degrees of freedom. This is, however, computationally expensive and was not done for this proof-of-principle calculation.

The Poincare section of the optimized equilibrium is plotted on the top right of Fig. 1. As expected, $ΓCB$ is no longer a flux surface since $Bc$ has been changed relative to the un-optimized equilibrium shown on the top left panel of Fig. 1. The difference between the magnetic field generated by the coils pre- and post-optimization is of the order of 1% of the total magnetic field, i.e., $δB/B∼1%$⁠. We observe that the targeted island chains are no longer visible. This optimization demonstrates that the parameter space ${Vmn}$ is suitable for stellarator optimization.

New resonances appeared close to the computational boundary $ΓCB$⁠, in particular one with mode number $(m,n)=(10,5)$⁠. The residue associated with this resonance is not included in the objective function, which explains why the optimizer was able to converge to this state. One approach to reduce the size of the $(m,n)=(10,5)$ island chain is to add the square of its associated residue to the objective function and continue the optimization from the previous optimized state. The optimizer is then able to reduce the $(m,n)=(10,5)$ island size (data not shown). This “stepped” optimization process, in which the objective function has to be modified multiple times, is avoidable if a global diagnostic to measure the extent of magnetic islands and chaos is used instead of a local diagnostic, such as the Greene's residues. This will be the topic of a future publication.

The rotational transform profile of the optimized equilibrium is plotted on Fig. 2. We observe that the rotational transform profile after optimization is of the same order as before the optimization—it still crosses the same rationals, but these rationals do not generate noticeable islands anymore!

We now look at two different optimizations of the perfectly conducting wall equilibrium. In the first one, we only target the residues in the vacuum region, i.e., residues with indices 1–3 of Table I, and the parameter space is the profile of current flowing in the plasma volumes ${Iϕ,lv}l={1,…,7}$⁠. When more residues listed in Table I are included in the objective function, the optimizer is not able to lower the objective function sufficiently to observe an effect on the island width. This might be because not enough degrees of freedom were provided. Increasing the number of volumes N_vol in SPEC, consequently increasing the number of degrees of freedom ${Iϕ,lv}l={1,…,Nvol}$⁠, is however not a guaranteed solution, since it adds additional topological constraints on the magnetic field, and some island chains might remain undetected. The optimization of the injected current profile ${Iϕ,lv}l={1,…,7}$ could, however, be combined with other parameters, for example, the coil optimization with degrees of freedom ${Vmn}$⁠, to target more residues. Note that this optimization could also be achieved in the initial free-boundary equilibrium, but the perfectly conducting wall equilibrium was considered here for simplicity.

In the second optimization, all residues listed in Table I are included in the objective function, and the geometry of the perfectly conducting wall $Γw$ is optimized. The selected degrees of freedom are the modes $Rmnw$ and $Zmnw$ with $(m,n)=(1,1), (1,2), (2,1)$⁠. Low-order Fourier modes are chosen for their capacity to penetrate deeper in the plasma.¹ In principle, higher-order Fourier modes could be used; however, the required deviation from the initial state might need to be large to minimize the objective function and thus the optimum state is difficult to find for the optimizer. Note that the poloidal angle defined by the Fourier representation of the perfectly conducting wall $Γw$ is generally not a straight-field line angle. Which Fourier mode, $Rmnw$ and $Zmnw$⁠, affects the targeted rationals is thus a non-trivial question—in general, the residue associated with the n/m rational might not be affected by the modes ${Rmnw,Zmnw}$!

Figure 1 shows the result of the optimization of ${Iϕ,lv}l={1,…,7}$ (bottom left) and the optimization of $Γw$ (bottom right). Comparing both Poincare plots with the initial equilibrium (top left), we observe that the targeted residues have indeed been minimized—the islands are now much smaller and are no longer visible in some cases.

The rotational transform profiles are plotted on Fig. 2 and the total enclosed toroidal current on Fig. 3. Again, the rotational transforms are of the same order as for the un-optimized case. Regarding the optimized current profile, the total injected current is $ΔIϕ=−5.2$ kA, less than 20% of the initial total enclosed toroidal current. Interestingly, the optimizer found an equilibrium close to the initial state and did not converge to the obvious global minimum, where a large toroidal current would be injected in the volumes to lift the rotational transform profile in order to push all the targeted rationals out of the plasma. Similarly, for the optimization of $Γw$⁠, a trivial solution would be the axisymmetric solution, with $Rmnw=Zmnw=0$ for $n≠0$⁠. Again, since a local optimizer is being used, the optimum found is a local optimum close to the initial state, and not the global optimum.

The normalized value of the objective function as a function of the number of iterations is plotted for each optimization in Fig. 4. The optimization of $Γw$ (orange crosses) cannot be compared with the other two, since it uses a different objective function. We see that the optimization of ${Iϕ,lv}l={1,…,7}$ (blue circles) saturates at a larger value than the optimization of ${Vmn}$ (green squares), despite optimizing the same objective function (the same residues were selected). Intuitively, this can be understood by noticing that the number Fourier harmonics of the magnetic field that are affected by the plasma current is smaller than those that can be affected by the ${Vmn}$⁠. In that sense, the space of possible perturbations of the magnetic field for influencing the islands is more constrained in the ${Iϕ,lv}l={1,…,7}$ optimization.

The optimization was run in parallel on $2n+1$ cores of Intel Broadwell processors at 2.6 GHz, where n is the number of degrees of freedom of the optimization. Each core computed a different SPEC equilibrium when evaluating the finite difference estimate of the objective function gradient. The ${Vmn}$ optimization required 165 equilibrium calculations, the ${Iϕ,lv}l={1,…,7}$ optimization 187, and the $Γw$ optimization 309. The total CPU time for execution was, respectively, $∼41$ days, $∼10$ days, and $∼21$ days for a total wall-clock time of, respectively, $∼65$⁠, $∼17$⁠, and $∼39$ h. As expected, the optimizations of the perfectly conducting wall equilibrium were faster, as fixed-boundary calculations are faster than free-boundary ones.

Note that the presented optimizations did not take advantage of the full parallelization of SPEC—only a single core computed each SPEC equilibrium, while SIMSOPT allows the user to use $M≤Nvol$ CPUs on $N≤2n+1$ SPEC instances, which would speedup the computation greatly. Nevertheless, our optimizations show that the total time required to perform a SPEC optimization is small enough to be considered in more advanced stellarator optimizations.

In this paper, we showed the first fixed- and free-boundary, multi-volume, finite β SPEC equilibrium optimizations of a classical stellarator using the SIMSOPT framework. The objective function was constructed from the Greene's residues of selected rational surfaces. Different parameter spaces were considered: either the boundary of a perfectly conducting wall surrounding the plasma, or the enclosed toroidal current profile, or the vacuum field produced by the coils were optimized.

In all three optimizations, it was possible to reduce the objective function significantly, which in turn translated to a reduction of the targeted magnetic island width. In the case of the current profile optimization, the volume occupied by magnetic islands and chaotic field lines was considerably reduced in the vacuum region surrounding the plasma. While not important for plasma confinement, the control of the magnetic field topology in the vacuum region is of paramount importance for divertor design and operation. In the case of the perfectly conducting wall optimization, the volume occupied by magnetic surfaces was increased both in the vacuum region and in the plasma volume. Finally, in the case of the optimization of the vacuum field produced by the coils, the islands present in the initial un-optimized equilibrium were reduced in size, but additional rationals and island chains emerged, since their associated residues were not included in the objective function.

Different measures of the magnetic field integrability are currently being considered to overcome the shortcomings of Greene's residues. Indeed, as observed in the coil optimization, Greene's residue is a local measure, which requires input from the user—any rational surface emerging during the optimization remains undetected. Ideally, a global measure is thus required. The volume occupied by chaotic field lines, measured from the fractal dimension of their Poincare map,⁸ could be a possible solution.

In the work by Landreman et al.,²² it has been shown that SPEC could be coupled to VMEC in order to achieve an optimization in a vacuum, where both quasi-symmetry and nested flux surfaces could be obtained. The obvious next step is then to perform a combined SPEC-VMEC finite-β optimization for good magnetic surfaces as well as other metrics, such as quasi-symmetry. This will require to deviate significantly from a classical stellarator geometry, as opposed to what has been presented in this paper, and the computational cost of the optimization will be greater due to the increased complexity of the problem. While this task does not require extensive development of the optimization tool SIMSOPT, a number of convergence issues in SPEC arising in strongly shaped geometries have to be solved. Alternative Fourier representations of the volumes interfaces³² are currently being implemented to attempt to solve this issue. This will be the topic of future investigations.

The authors would like to acknowledge the support from the SIMSOPT development team and thank Z. Qu, B. Medasani, and C. Zhu for useful discussions. This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No 101052200—EUROfusion). Views and opinions expressed are, however, those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. This work was supported by a grant from the Simons Foundation (No. 560651, ML).

The authors have no conflicts to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.