2.1. Naive globalization of phase I

A first naive approach to globalization is to perturb the geometry and currents of flat electromagnetic coils, as well as the parameters of the magnetic axis to generate a set of initial starting points. These starting points are used to initialize phase I, where the local gradient-based optimization algorithm BFGS is used to find a local minimum of the objective in the neighbourhood of the perturbed initial guess. The physics and engineering quantities that we target are $\iota _{\text {target}}=0.9$ using $n_{\text {fp}}, n_{\text {coils per hp}}=2$, $d_{\min }=0.1\,{\rm m}$, $\kappa _{\max }=5\,{\rm m}^{-1}$ $\kappa _{\text {msc}}=5\,{\rm m}^{-2}$, the target coil length $L_{\text {target}}$ is varied between 40 and 80 m. For each value of coil length, we solve (2.1) with distinct initial guesses 16 times. The initial guesses are obtained by perturbing the Fourier coefficients of initially flat coils with normally distributed noise with zero mean and standard deviation $\epsilon$. One has substantial freedom to decide many Fourier coefficients to perturb, as well as the standard deviation of the perturbation and how it decays. If $\epsilon$ is too large then the local optimizer might be sent into an uninteresting region of coil parameter space. If $\epsilon$ is too small, and it might not fully explore the set of feasible coil designs. For the experiment here, we perturb the currents, $\bar {\eta }$, and only the first two Fourier harmonics of the coils, and magnetic axis. Extensive (but tedious and computationally intensive) tuning revealed that the standard deviation of the noise that discovered devices with lowest quasi-symmetry error was $\epsilon = 0.01\,{\rm m}$. We find that the quality of the minima found depends strongly on the choice of $\epsilon$, the number of Fourier modes that are perturbed, and that there is no way to know a priori that this value works well for other stellarator designs.

In figure 2, we plot the trade-off between the on-axis quasisymmetry error and total coil length. We find, as expected, that lower quasisymmetry errors are attainable with longer coils. It is clear that making no attempt at globalization to solve this optimization problem might result in a sub-optimal coil set. The problem of multiple minima becomes more pronounced for devices with longer coil lengths because there are more, intricate ways for the coils to arrange themselves. Finally, at around $50\,{\rm m}$, there is a device that appears as an outlier with remarkably low near-axis quasisymmetry error. This hints that the naive approach is missing interesting solution branches.

Figure 2. (a) We plot the on-axis quasisymmetry error with respect to the total coil length used, where we observe that there are multiple minima. (b) We plot four distinct local minima when the total coil length is $49.81\,{\rm m}$, where the colours of the coils correspond to those of the squares on the left. The coil sets’ associated magnetic axes are not plotted to avoid clutter in the image.

2.2. The TuRBO globalization of phase I

An alternative to the naive approach is to use a global optimization algorithm. The coil set it finds is then given as the starting point in phase I to a local gradient-based optimization algorithm BFGS. This global-to-local approach was explored in Glas et al. (Reference Glas, Padidar, Kellison and Bindel2022), where stochastic global optimization of coil geometries was done using the DTuRBO algorithm (Padidar et al. Reference Padidar, Zhu, Huang, Gardner and Bindel2021), an extension of the gradient-free global optimization algorithm TuRBO (Eriksson et al. Reference Eriksson, Pearce, Gardner, Turner, Poloczek, Wallach, Larochelle, Beygelzimer, d'Alché-Buc, Fox and Garnett2019) that takes into account gradients. The TuRBO algorithm solves the minimization problem

(2.5)\begin{equation} \left. \begin{aligned} & \min_{\boldsymbol{x} \in \mathbb{R}^N}\ f(\boldsymbol{x})\\ & \ell_i \leq x_i \leq u_i,\quad \text{for } i = 1 \ldots N, \end{aligned} \right\} \end{equation}

where $N$ is the dimension of the optimization problem, $\ell _i, u_i$ are respectively lower and upper bounds on the components of the control variables. That is, the algorithm attempts to find the global minimum $f(\boldsymbol {x})$ on an $N$-dimensional hyper-rectangle. The TuRBO algorithm completes a preliminary exploration phase by quasi-uniformly sampling the hyper-rectangle to find interesting areas of parameter space. Then, it uses trust-region-based Bayesian optimization algorithms to further refine the solution.

It is difficult to define sensible lower and upper bounds in (2.5) on the unknowns when they are the Fourier coefficients of the coils. There should be some decay in the size of the boxes with increasing mode number, but it unclear how this decay should be chosen. In Glas et al. (Reference Glas, Padidar, Kellison and Bindel2022), manufacturing errors were projected onto a Fourier basis, from which reasonable bounds on the Fourier coefficients were inferred. In this work, we propose a different way of determining the box constraints for deterministic coil design problems, however, our technique is generic enough to be used in the stochastic context as well.

2.2.1. Bounding boxes

It is straightforward to constrain the coil currents to $-1 \leq \mu _0 I_i \leq 1$, where $\mu _0 = 4{\rm \pi} \times 10^{-7}$ is the magnetic constant. The parameter $\overline \eta$ can also be constrained by $0\leq \bar {\eta } \leq 2$.

We still use a Fourier series to represent the coils and axis, but change the degrees of freedom from Fourier coefficients to spatial coordinates, linked to each other by the discrete Fourier transform. In this way, physically meaningful box constraints can be provided to TuRBO. These positions, or anchor points, in cylindrical coordinates, $(r_i, \theta _i, z_i)$, are constrained to the box $[0, 1 + R_{\textrm {minor}}] \times [\theta _i-\Delta \theta /2, \theta _i + \Delta \theta /2] \times [-R_{\textrm {minor}}, R_{\textrm {minor}}]$, where $R_{\textrm {minor}} = L_{\textrm {target}}/2{\rm \pi}$, and $\Delta \theta = \rho (2{\rm \pi} /2n_{\textrm {fp}} n_{\textrm {coils per hp} })$, $\theta _i = (2{\rm \pi} /2n_{\textrm {fp}}n_{\textrm {coils per hp} })(i+1/2)$, $R_{\textrm {minor}}$ is the radius of the perfectly circular coil with length $L_{\textrm {target}}$, $\Delta \theta$ defines a cylindrical sector that the coil can occupy and the $\rho$ multiplier ensures that the coil bounding boxes overlap somewhat. Randomly sampling $2N_{f,c}+1$ points on the box will frequently generate coils with complex geometries (figure 3a). These geometries are not particularly useful and we would like TuRBO to avoid wasting time in these uninteresting areas of parameter space. These configurations can be avoided, while still using the same anchor points, by re-indexing them so that the points are ordered counterclockwise about their barycentre in the $RZ$ plane (figure 3b). This unravels the coils. Finally, these coil coordinates are converted to Cartesian coordinates, projected onto a Fourier basis, then used to evaluate the near-axis objective (2.4).

Figure 3. The randomly generated coil in cylindrical coordinates in (a) is too complex. By $Re$-indexing the nodes by ordering them about their barycentre in the $RZ$ plane, the coil untangles in (b). (c) A set of coils, the magnetic axis and their associated bounded boxes in physical space are illustrated. The full stellarator is not shown, but can be obtained by applying symmetries.

We also generate $(R_i, \phi _i, Z_i)$ positions, or anchor points, for the magnetic axis, where $(R_i, Z_i)$ are constrained to the box $[1-1/(1+n^2_{\textrm {fp}}), 1+ 1/(1+n^2_{\textrm {fp}})] \times [-0.2, 0.2]$, and $\phi _i$ is a uniformly spaced cylindrical angle in $2{\rm \pi} /n_{\textrm {fp}}$. The bounding boxes for $R_i$ and $Z_i$ are informed by the fact that quasi-axisymmetric magnetic axes have zero helicity. Helicity is an integer associated with the axis geometry that measures how many poloidal transits the axis normal vector makes as the axis traces one toroidal revolution. It was shown in Rodriguez, Sengupta & Bhattacharjee (Reference Rodriguez, Sengupta and Bhattacharjee2022) that magnetic axes of the form

(2.6)\begin{gather} R(\phi) = 1+a\cos(n_{\text{fp}}\phi), \end{gather}

(2.7)\begin{gather}Z(\phi) = b\sin(n_{\text{fp}}\phi), \end{gather}

have zero helicity when $a < 1/(1+n^2_{\textrm {fp}})$. The situation is more delicate if additional Fourier harmonics are used to represent the axis (as we do here), but in the case of a single dominant harmonic, this bound on $a$ is a well-informed estimate. The box associated with $Z_i$ was chosen to prevent large excursions from the $Z=0$ plane, although the size of this box was somewhat arbitrary and other values might have been used.

There is one final detail that ensures that the resulting axis curve is stellarator symmetric. After sampling $N_{f,a}+1$ points from the bounding box associated with the radial coordinates $R_i$, they are arranged in an $2N_{f,a}+1$-sized array as follows:

(2.8)\begin{equation} [R_0-s, R_1-s, \ldots, R_{N_{f,a}}-s, R_{N_{f,a}}-s, \ldots, R_1-s], \end{equation}

where $s$ is a shift to ensure the array has mean 1. After the shift, the radial positions of the magnetic axis may lie slightly outside the original bounding box. Similarly, we sample $N_{f,a}$ points from the bounding box associated with the vertical coordinates $Z_i$, and arrange them in a $2N_{f,a}+1$-sized array as follows:

(2.9)\begin{equation} [0, Z_1, \ldots, Z_{N_{f,a}}, -Z_{N_{f,a}}, \ldots, -Z_1]. \end{equation}

The stellarator symmetric $R$ and $Z$ harmonics of the magnetic axis are obtained by discrete Fourier transform of these arrays and provided to the near-axis objective (2.4).

The bounding boxes of different coils can overlap each other and all the bounding boxes of the coils necessarily overlap the one of the magnetic axis, since we do not search for windowpane coils. Thus, it is possible for the generated coils to be linked, or for the coils to not be linked with the magnetic axis. In addition, it is still possible that the magnetic axis sampled from the bounding boxes to have non-zero helicity. It is for this reason that we add penalty terms to the objective $f_{\textrm {axis}}$ that avoid these undesirable configurations

(2.10)\begin{align} f(\boldsymbol{c} , \boldsymbol{I}, \boldsymbol{a}, \bar{\eta}) = f_{\text{axis}}(\boldsymbol{c} , \boldsymbol{I}, \boldsymbol{a} , \bar{\eta}) + w\left[\sum_{i}\sum_{j>i} l_{i,j}(\boldsymbol{c}) + \sum_{i}|l_{i,\text{axis}}(\boldsymbol{c}, \boldsymbol{a})-1| + h(\boldsymbol{a})\right], \end{align}

where we abuse notation here and use $\boldsymbol {c}, \boldsymbol {a}$ to represent the vector of anchor points on the coils and axis, respectively. These positions can be converted to the Fourier coefficients by using a discrete Fourier transform. The value $l_{i,j}$ is the absolute value of the linking number between coils $i$ and $j$, and $l_{i, \textrm {axis}}$ is the absolute value of the linking number between coil $i$ and the magnetic axis. For the linking number calculation, we use the implementation in Bertolazzi, Ghiloni & Specogna (Reference Bertolazzi, Ghiloni and Specogna2019) since it returns a linking number with certified accuracy by tracking rounding errors in the calculation. Finally, $h(\boldsymbol {a})$ is the absolute value of the axis helicity. The weight $w=10^7$ is chosen to make the terms that it multiplies greater than $f_{\textrm {axis}}$. The penalty terms multiplied by $w$ are zero when the design variables describe a stellarator with unlinked coils that each wrap once poloidally around a magnetic axis with zero helicity. Since the objective is no longer differentiable due to the linking number and helicity calculation, we use TuRBO rather than DTuRBO.

We also designed stellarators using an axis bounding box independent of $n_{\textrm {fp}}$, $(R_i, Z_i) \in [0.85, 1.15] \times [-0.15, 0.15]$, and neglecting to include the helicity penalty. In a postprocessing step, devices that had an axis with non-zero helicity were discarded. This approach risks discarding devices and wasting computational resources but it is capable of producing comparable results.

2.3. Illustration of TuRBO globalization followed by phase I

To illustrate the global-to-local procedure, we examine the TuRBO globalization procedure followed by the phase I optimization for a single stellarator. We use a batch size of 100 and restrict TuRBO to a budget of 15 000 function evaluations, of which 1000 are used in the initial exploration phase; TuRBO's internal Gaussian process model is only trained using the most recent 1000 function evaluations for speed. Finally, the underlying Fourier discretization uses $N_{f,c}=N_{f,a}=2$ for the coils and axis during this global exploration phase. The coil set and magnetic axis from the global search are used as the initial guess to the BFGS algorithm to polish the solution. During this phase, the number of Fourier modes used is increased to $N_{f,c}=N_{f,a}=6$. Before BFGS obtains a good enough approximation to the inverse Hessian, especially during the first few iterations, the line search might accept a large enough step such that the coils become interlocked or unlinked with the magnetic axis. To avoid this, we use a linking number calculation to trigger the line search.

In figure 4 we provide the objective function evaluations over the course the TuRBO globalization then phase I of figure 1. The TuRBO phase of the run has large variations in the objective function, varying over several orders of magnitude. When the objective is larger than $10^7$, this likely corresponds to a device that has an axis with non-zero helicity, linked coils or coils unlinked with the axis; TuRBO quickly finds an interesting region of parameter space and it could have used a smaller budget. The notable jump in the objective function value around 11 000 function evaluations occurs when TuRBO no longer can make progress in the trust region, so it is discarded and another one is initialized. The best configuration found by TuRBO is provided and labelled Design 4.A. The coil set at the end of phase I (Design 4.B) has changed substantially from those obtained after the TuRBO phase.

Figure 4. The value of the objective function when applying TuRBO then BFGS. For the TuRBO phase, the objective values are evaluations of (2.10), which includes the non-differentiable penalties and the ‘iteration #’ is actually number of function evaluations. For the BFGS phase, the objective values are evaluations of (2.2), and ‘iteration #’ refers to the objective function values after an accepted line search. The red line is the best function evaluation so far. (b) We plot the best stellarator found by TuRBO, given by Design 4.A, provided as the initial guess to BFGS in phase I. Then, the stellarator found by BFGS after increasing the number Fourier modes is given by Design 4.B.

2.4. Comparison of devices from TuRBO and naive globalization after phase I

In this section, we compare the performance of TuRBO with the more naive approach by solving (2.1) for many target coil length values when $\iota =0.9$, $n_{\textrm {fp}}, n_{\textrm {coils per hp}}=2$. To this end, we execute the TuRBO globalization followed by phase I, as illustrated in § 2.3 and figure 4. In both naive and TuRBO approaches, the workflow is executed 16 times per target coil length. The quasisymmetry error of the resulting configurations are provided in figure 5. TuRBO finds devices that outperform those found by the naive approach, and we plot two devices from each algorithm, called Designs 5.A and 5.B. We find that Design 5.B is from a genuinely different solution branch. To drive this home, we use it as an initial guess, we solve (2.1) for multiple different target coil lengths. This solution branch persists on both sides, and was evidently even missed by TuRBO. This highlights the difficulty of finding global minima.

Figure 5. (a) We provide the trade-off between axis quasisymmetry error and total coil length used for local minima found using the naive approach and TuRBO with $\iota =0.9$, and $n_{\textrm {fp}}, n_{\textrm {coils per hp}}=2$. (b) We have plotted two coil sets, called Designs 5.A and 5.B, corresponding respectively to the red, and green square markers. Using Design 5.B as an initial guess, we solve (2.1) for different target coil lengths, resulting in the devices corresponding to the blue crosses. The axis QA error refers to the terms in (2.2).

With very modest input from the practitioner, TuRBO does a good job of localizing performant minima. We find that the choice of the overlap factor $\rho$ only modestly affected the quality of the minimizers found. In contrast, one has to decide on how many Fourier modes to perturb and the size of each perturbation to use with the naive approach. This can have a large effect on the quality of the minimizers found.

4.1. Illustration of phases II and III

In this section, we give more details on the final two phases of the coil design algorithm to illustrate how it works. The configurations obtained from the near-axis formulation (phase I) only attempt to find coils that target quasisymmetry on the magnetic axis. The formulation does not control anything about the magnetic field away from the axis and, as a result, magnetic surfaces may not exist on a large volume (Lee et al. Reference Lee, Paul, Stadler and Landreman2022). The BoozerLS surfaces (phase II) are robust and capable of healing generalized chaos (Giuliani et al. Reference Giuliani, Wechsung, Cerfon, Landreman and Stadler2023), although it can be computationally expensive. Finally, the BoozerExact (phase III) surfaces polish the coil set for precise QA and are computationally much cheaper.

The algorithm for phase II is given in Algorithm 1. The first step of the algorithm is to use the magnetic axis obtained from phase I to compute an approximate magnetic surface from (4.10). Then, the optimization problem (4.5) is solved, but we do not attempt to fully converge the coil sets and only complete a total of 300 iterations of BFGS. Each time BFGS is restarted, an attempt is made to increase the number of Fourier modes used to represent the surface. Starting from $m_\textrm {pol},n_\textrm {tor}=2$, we try to compute a BoozerLS surface. If the solve succeeds, we use the surface as an initial guess and increase $m_\textrm {pol}, n_\textrm {tor}$ by 1. If the surface solve fails after increasing the surface degree, then the degree is reverted to the previous one. If the solve succeeds, the continuation in degree continues until ${m_\textrm {pol},n_\textrm {tor}=4}$. Progressively increasing the number of Fourier modes makes the procedure more robust, as we find self-intersecting surfaces are less likely to occur.

The weight $w_r$ is chosen such that the Boozer residual term is an order of magnitude larger than the non-quasisymmetry penalty. This term favours nested flux surfaces and improves the robustness of the optimization algorithm. However, it also competes with the QA error and so the quality of quasisymmetry obtained at this stage can be limited. Since the stellarators obtained here are not fully converged and strongly depend on the value of $w$, we do not analyse the physics properties of the stellarators here yet. In figure 7, we provide Poincaré plots of a device just after the near-axis optimization, and then as BoozerLS surfaces are added to optimize for nested flux surfaces. The first panel shows that the near-axis optimization algorithm does not guarantee nested surfaces far away from the magnetic axis. The final two panels illustrate that as the surfaces are added, the volume with nested flux surfaces increases.

Figure 7. Poincaré plots at $\phi =0$ after phases I and II for aspect ratios $AR=20, 10$. The point where the magnetic axis intersects the $\phi =0$ plane corresponds to the red dot in (a). Cross-sections of the surfaces on which QA is optimized are red in (b,c).

Algorithm 1 Phase II: BoozerLS

Algorithm 2 Phase III: BoozerExact

After phase II terminates, we proceed to phase III, where we optimize for precise quasisymmetry. It can be shown that the island width scales like $\sqrt {|B_{n,m}|/m\iota '}$ (Hudson, Monticello & Reiman Reference Hudson, Monticello and Reiman2001), where $\iota '$ is the magnetic shear. Thus, for large enough $m$, the island width should be small. Informed by this estimate, we retain the BoozerLS formulation when the rotational transform on the surface is within 1 % of a rotational transform of the form $\iota = n_{\textrm {fp}}/m, 2n_{\textrm {fp}}/m$ for $m=1, \ldots, 15$. During this final phase, the weight $w$ is set to zero on BoozerExact surfaces, but it maintained for BoozerLS surfaces in the neighbourhood of low-order rationals. At the start of every BFGS run, we attempt to increase the degree of the surface representation as we did during phase II, up to $m_{\textrm {pol}},n_{\textrm {tor}}=10$ for BoozerExact surfaces, and 4 for the remaining BoozerLS surfaces. When all surfaces are of type BoozerExact, we cap the total number of BFGS iterations to 20 000. When there is a surface of type BoozerLS, we cap the total number of BFGS iterations 4000. The algorithm for phase III is given in Algorithm 2.

During these last two phases of coil optimization, we also increase the number of Fourier harmonics in the coils to $N_{f,c}=16$, because using too few Fourier modes in the coils can limit the attainable quality of quasisymmetry.

4.2. Comparison of devices from TuRBO and naive globalization after phase III

After the full workflow is completed (globalization, then phases I, II, then III), we compare the devices found depending on whether the TuRBO or naive globalization approach was applied just before phase I. In figure 8, we plot volume quasisymmetry error for configurations found for both globalization techniques. It appears that the highly favourable devices discovered by TuRBO do not necessarily translate to devices that have comparatively favourable volume quasisymmetry. We also find the TuRBO devices resulted in many more configurations for shorter coils than the ones found by the naive approach, although there is a larger spread of device performance. The Design 5.B does not persist upon optimizing for volume quasisymmetry, which we suspect is because the coils are too close to the magnetic axis ($16.5\,\textrm {cm}$). As a result, when the surface generated by the near-axis formula (4.10) is used at the start of phase II, the BoozerLS surface solve in (4.5b) does not converge and the optimization cannot proceed. The devices discovered by TuRBO also do not appear to outperform those found by the naive approach. Nevertheless, we emphasize that this may only be because the globalization is performed on the near-axis problem just before phase I and not directly to the volume QA optimization, just before phase II. If globalization is applied directly to the BoozerLS phase of the workflow, it is possible that new and interesting devices could still be discovered. Finally, the magnitude and form of the perturbation $\epsilon$ used in the naive approach was determined after somewhat lengthy tuning (§ 2.2). There are much fewer ad hoc choices in the bounding box approach used by TuRBO.

Figure 8. Taking the stellarators from figure 5, and using them as initial guesses for the volume QA optimization phase (§ 4), we obtain the stellarators plotted above. So that the data sets can be more easily distinguished, we use smaller marker sizes for devices corresponding to the naive algorithm.

4.3. Computational details

The run time of each phase can vary depending on the design targets, as the difficulty of the optimization problems strongly depends on the total coil length used by the device. In previous work, we used various levels of parallelism simultaneously such as MPI, OpenMP and SIMD. In this work, we take the approach of only allocating one core per stellarator and only use SIMD parallelism locally to a core when possible. The full workflow for a single device takes of the order of a day or two on a single core. However, this duration can vary substantially depending on how accurately each individual optimization problem is solved or if more cores are used per problem.

5.1. Comparison with previous devices

As a first examination, we select the devices in our database that have the closest design targets to previous configurations computed in Giuliani et al. (Reference Giuliani, Wechsung, Stadler, Cerfon and Landreman2022b), where devices have aspect ratio 6, $\iota =0.42$, $n_{\textrm {coils per hp}}=4$, total coil length $72\,\textrm {m}, 80\,\textrm {m}, 88\,\textrm {m}, 96\,\textrm {m}$ and ${n_{\textrm {fp}}=2}$. The devices in QUASR that are closest to these devices have the same number of coils, aspect ratio $6.66$ or $5$ and target mean rotational transform $0.4$. In figure 9, these devices are compared, where we observe that the QUASR devices perform comparably, although notably the precise QA device with length $96\,\textrm {m}$ is better. This might be because we only allow a total of 20 000 iterations of BFGS during phase III, while in generating the precise QA coil sets, more than double the number of iterations were allowed.

Figure 9. Comparison of devices from QUASR with aspect ratio $6.66$ and $5$, target mean rotational transform $\iota =0.4$ and the precise QA coil sets of Giuliani et al. (Reference Giuliani, Wechsung, Stadler, Cerfon and Landreman2022b) with aspect ratio $6$ and mean rotational transform $0.42$. All devices have $n_{\textrm {fp}}=2$ and $n_{\textrm {coils per hp}}=4$ here. The grey region corresponds to the Earth's background magnetic field of $50\,\mathrm {\mu }\textrm {T}$, and volume QA error refers to the mean non-QA ratio in (4.1).

5.2. A stellarator with a rotational transform profile passing through a low-order rational

Consider a device with two field periods, $n_{\textrm {coils per hp}}=3, \iota _{\textrm {target}}=0.5$ and aspect ratio 4. Since the target mean rotational transform is a low-order rational, we would expect there to be a strong likelihood of an island chain somewhere in the toroidal volume. However, the island chains are small, as shown in figure 10.

Figure 10. Poincaré plots and rotational transform profiles of a device with aspect ratio 4, mean rotational transform 0.5 and $n_{\textrm {fp}}=2$. We visually could not identify trajectories associated with islands. The horizontal blue lines on the rotational transform profiles indicate the low-order rational $\iota =1/2$. The vertical red lines on the rotational transform profiles indicate the normalized toroidal flux label on which quasisymmetry was optimized.

5.3. Impact of design targets on attainable quasisymmetry

In figure 11, we plot the volume QA error with respect to coil length for various values of $n_{\textrm {fp}}=2, 3, 4$, when optimizing for QA on a single surface of aspect ratio 20. Since this is such a high aspect ratio surface, this set-up is a close approximation to the near-axis problem in § 2, although the objective is asking for more out of the magnetic field. A notable difference with the near-axis results is that increasing the number of field periods does not appear to improve the attainable on-axis quasisymmetry, as is observed in figure 6. One possible reason for this is that the near-axis design problem (2.1) only optimizes for quasisymmetry to first order and introducing second-order near-axis penalties might explain this discrepancy. There do appear to be some highly performant devices that are outliers, so it might also be that our algorithm is just not discovering those solution branches. There is also a clear preference for $n_{\textrm {fp}}=2$ stellarators as the algorithm has difficulty finding devices with precise quasisymmetry for $n_{\textrm {fp}}=3,4$. Devices with $n_{\textrm {fp}}=1,5$ are not shown in the figure, but in both cases, the algorithm had varying difficulty finding devices with precise QA.

Figure 11. Devices on which quasisymmetry is optimized on a surface of aspect ratio 20, with a target rotational transform $\iota =0.5$. The grey region corresponds to the Earth's background magnetic field of $50\,\mathrm {\mu }\textrm {T}$, and volume QA error refers to the mean non-QA ratio in (4.1).

We also examine the trade-off between the quality of quasisymmetry, total coil length and device aspect ratio in figure 12. Increasing the target mean rotational transform does not appear to greatly affect the quality of attainable quasisymmetry, while as expected decreasing the device aspect ratio appears to have a much larger negative impact.

Figure 12. Trade-off between quality of quasisymmetry and total coil length for various device aspect ratios, and target mean rotational transform when $n_{\textrm {fp}}=2$. The grey region corresponds to the Earth's background magnetic field of $50\,\mathrm {\mu }\textrm {T}$, and volume QA error refers to the mean non-QA ratio in (4.1).

5.4. Maximum elongation

Finally, we study the values of elongation appearing in the database when $\iota _{\textrm {target}}=0.6$ and an aspect ratio of 10. High elongation increases the ratio of surface area to volume of the magnetic surfaces, moreover, it can increase bunching of the flux surfaces and this can have negative implications for stability (Goodman et al. Reference Goodman, Camacho Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach, Beidler and Helander2023). For each stellarator in our database, we compute the maximum elongation in the device by computing a surface with aspect ratio 80 in the neighbourhood of the magnetic axis. Then, we compute $N=10$ cross-sections of that surface at cylindrical angles $\phi ={\rm \pi} (i+1/2)/n_{\textrm {fp}}/N \textrm { for } i= 0, \ldots, N-1$. For each cross-section, we fit an ellipse in a least squares sense (Halır & Flusser Reference Halır and Flusser1998) and compute the ratio of the ellipse's major to minor axes. The value reported here is the maximum ratio observed at these cross-sections.

We did not target any particular value of elongation but nevertheless there do appear to be favoured values (figure 13). For $n_{\textrm {fp}}=1$ and $2$, the lowest values of quasisymmetry error occur for an elongation around 7, while for $n_{\textrm {fp}}=3$ an elongation of 4 is favoured. It is unclear whether this picture would change if a stage one optimization for QA were done instead, without coils. In Goodman et al. (Reference Goodman, Camacho Mata, Henneberg, Jorge, Landreman, Plunk, Smith, Mackenbach, Beidler and Helander2023), various QI stellarator designs were proposed, where a maximum elongation below approximately 6 was targeted. The stellarators in QUASR here illustrate that a good approximation of QA can also be found when requiring a maximum elongation below 6 too. We also observe that as expected, the lower quasisymmetry errors occur for longer coil lengths.

Figure 13. The volume quasisymmetry error with respect to the maximum elongation on a very high aspect ratio surface in the neighbourhood of the magnetic axis, for various values of $n_{\textrm {fp}}$. For all configurations plotted, $\iota _{\textrm {target}}=0.6$, where the aspect ratio of the outermost surface is 10. The grey region corresponds to the Earth's background magnetic field of $50\,\mathrm {\mu }\textrm {T}$, and volume QA error refers to the mean non-QA ratio in (4.1).