Hostname: page-component-cd9895bd7-gxg78 Total loading time: 0 Render date: 2024-12-23T06:25:07.619Z Has data issue: false hasContentIssue false
  • English
  • Français

Exploration of the parameter space of quasisymmetric stellarator vacuum fields through adjoint optimisation

Published online by Cambridge University Press:  20 December 2024

Richard Nies Open the ORCID record for Richard Nies [Opens in a new window] ,
Elizabeth J. Paul Open the ORCID record for Elizabeth J. Paul [Opens in a new window] ,
Dario Panici Open the ORCID record for Dario Panici [Opens in a new window] ,
Stuart R. Hudson Open the ORCID record for Stuart R. Hudson [Opens in a new window]  and
Amitava Bhattacharjee Open the ORCID record for Amitava Bhattacharjee [Opens in a new window]
Richard Nies*
Affiliation:
Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08543, USA Princeton Plasma Physics Laboratory, Princeton, NJ 08540, USA
Elizabeth J. Paul
Affiliation:
Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY 10027, USA
Dario Panici
Affiliation:
Department of Mechanical Engineering, Princeton University, Princeton, NJ 08543, USA
Stuart R. Hudson
Affiliation:
Princeton Plasma Physics Laboratory, Princeton, NJ 08540, USA
Amitava Bhattacharjee
Affiliation:
Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08543, USA Princeton Plasma Physics Laboratory, Princeton, NJ 08540, USA
*
†Email address for correspondence: [email protected]
Rights & Permissions [Opens in a new window]

Abstract

Optimising stellarators for quasisymmetry leads to strongly reduced collisional transport and energetic particle losses compared with unoptimised configurations. Although stellarators with precise quasisymmetry have been obtained in the past, it remains unclear how broad the parameter space is where good quasisymmetry may be achieved. We study the range of aspect ratios and rotational transform values for which stellarators with excellent quasisymmetry on the boundary can be obtained. A large number of Fourier harmonics is included in the boundary representation, which is made computationally tractable by the use of adjoint methods to enable fast gradient-based optimisation and by the direct optimisation of vacuum magnetic fields, which converge more robustly compared with solutions from magnetohydrostatics. Several novel configurations are presented, including stellarators with record levels of quasisymmetry on a surface, three field period quasiaxisymmetric stellarators with substantial magnetic shear, and compact quasisymmetric stellarators at low aspect ratios similar to tokamaks.

Keywords

fusion plasmaplasma confinementplasma devices
Type
Research Article
Information
Journal of Plasma Physics , Volume 90 , Issue 6 , December 2024 , 905900620
Check for updates
Creative Commons
Creative Common License - CCCreative Common License - BY
This is an Open Access article, distributed under the terms of the Creative Commons Attribution licence (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted re-use, distribution and reproduction, provided the original article is properly cited.
Copyright
Copyright © The Author(s), 2024. Published by Cambridge University Press

1 Introduction

Magnetic confinement fusion in toroidal geometry revolves around two main concepts, the axisymmetric tokamak and the three-dimensionally shaped stellarator (Spitzer Reference Spitzer1958). The axisymmetry of the tokamak leads to greater engineering simplicity and good confinement of particles, at the cost of difficulty maintaining a stable steady-state plasma due to the need for a large plasma current. Stellarators avoid these issues by mostly relying on the external coils to provide the confining magnetic field, but they are generally plagued by poor neoclassical and energetic particle confinement at reactor-relevant low collisionality. However, the confinement may be returned to tokamak-like levels by making stellarators omnigenous (Hall & McNamara Reference Hall and McNamara1975) through careful optimisation of the three-dimensional geometry, as has been verified experimentally in the Wendelstein-7X device (Beidler et al. Reference Beidler, Smith, Alonso, Andreeva, Baldzuhn, Beurskens, Borchardt, Bozhenkov, Brunner and Damm2021).

As a subset of omnigeneity, quasisymmetry (QS) (Nührenberg & Zille Reference Nührenberg and Zille1988) recovers tokamak-like guiding-centre dynamics through a symmetry in the magnetic field strength $B=B(s, M\theta _B - N \zeta _B)$, with the flux surface label $s$, the Boozer poloidal and toroidal angles $\theta _B$ and $\zeta _B$, and the QS helicities $M$ and $N$. Depending on the helicities, QS configurations may be either quasiaxisymmetric (QA, $N=0$), quasihelical (QH, $M, N \,{\neq}\, 0$) or quasipoloidal (QP, $M=0$). Several QS stellarators have been designed, e.g. the QA NCSX (Zarnstorff et al. Reference Zarnstorff, Berry, Brooks, Fredrickson, Fu, Hirshman, Hudson, Ku, Lazarus and Mikkelsen2001) and MUSE (Qian et al. Reference Qian, Chu, Pagano, Patch, Zarnstorff, Berlinger, Bishop, Chambliss, Haque and Seidita2023), and the QH HSX (Anderson et al. Reference Anderson, Almagri, Anderson, Matthews, Talmadge and Shohet1995), the latter two of which were constructed. Recently, configurations with record levels of QS have been obtained in optimisations of vacuum fields (Landreman & Paul Reference Landreman and Paul2022) and finite-pressure equilibria (Landreman, Buller & Drevlak Reference Landreman, Buller and Drevlak2022).

Most QS optimisations to date have relied on derivative-free algorithms (Drevlak et al. Reference Drevlak, Brochard, Helander, Kisslinger, Mikhailov, Nührenberg, Nührenberg and Turkin2013; Bader et al. Reference Bader, Drevlak, Anderson, Faber, Hegna, Likin, Schmitt and Talmadge2019; Henneberg et al. Reference Henneberg, Drevlak, Nührenberg, Beidler, Turkin, Loizu and Helander2019; Henneberg, Drevlak & Helander Reference Henneberg, Drevlak and Helander2020), or on gradient-based algorithms using finite differences to evaluate the gradient (Landreman & Paul Reference Landreman and Paul2022; Landreman et al. Reference Landreman, Buller and Drevlak2022). Due to the large parameter space of the stellarator boundary representation, these schemes are computationally expensive, which has motivated the introduction of automatic differentiation techniques (Dudt et al. Reference Dudt, Conlin, Panici and Kolemen2023) and adjoint methods (Paul, Landreman & Antonsen Reference Paul, Landreman and Antonsen2021). Moreover, the accuracy of gradients evaluated through finite differences can be sensitive to the step size, requiring careful scans in the numerical parameters to ensure success in the optimisation, e.g. avoiding artificial local minima.

Although many configurations with good QS have been obtained in the past, the degree of compatibility of QS with other objectives remains unclear.Footnote 1 Furthermore, slower optimisation algorithms limit the amount of shaping that can be taken advantage of in the optimisation. To remedy these gaps, we here perform QS optimisation of stellarators using adjoint methods. The fast optimisation facilitates the exploration of a large parameter space in the stellarator boundary shape, and the compatibility of QS with aspect ratio and rotational transform targets. We optimise vacuum magnetic fields, instead of considering the vacuum limit (zero plasma current and pressure) of magnetohydrostatics (MHS) with the imposition of nested flux surfaces, as computed e.g. by the codes VMEC (Hirshman, van Rij & Merkel Reference Hirshman, van Rij and Merkel1986) and DESC (Dudt & Kolemen Reference Dudt and Kolemen2020). The vacuum magnetic field solution converges more robustly than the solution from MHS due to the linearity of Laplace's equation, allowing us to go to high resolutions reliably.

In § 2, our objective function and optimisation methods are introduced. In §§ 3 and 4, we then present the obtained quasisymmetric configurations with varying aspect ratio and rotational transform, respectively. In § 5, we show how the optimised configurations may be further refined, e.g. by optimising for QS in the full volume (§ 5.1), or by improving the nestedness of flux surfaces in a QA with substantial magnetic shear (§ 5.2). Finally, we present our conclusions in § 6.

2 Methods

The objective function $f$ considered in the optimisations presented herein contains three targets

(2.1)\begin{equation} f = f_\mathrm{QS}^\star{+} w_\iota f_\iota + w_A f_A, \end{equation}

with weights $w_\iota$ and $w_A$. The objective includes a target for QS on the boundary $f_\mathrm {QS}^\star$, defined in (A4), another for the edge rotational transform $\iota _e$

(2.2)\begin{equation} f_\iota = \frac{1}{2}\left( \iota_e - \iota_T \right)^2, \end{equation}

with prescribed $\iota _T$, and finally an aspect ratio objective,

(2.3)\begin{equation} f_A = \frac{1}{2}\left( A - A_T \right)^2, \end{equation}

with prescribed $A_T$. The aspect ratio definition (see e.g. Landreman & Sengupta Reference Landreman and Sengupta2019) follows that employed in the VMEC code (Hirshman et al. Reference Hirshman, van Rij and Merkel1986)

(2.4)\begin{equation} A = \frac{R_\mathrm{maj}}{a_\mathrm{min}}, \quad \text{with} \ R_\mathrm{maj} = \frac{\mathcal{V}}{2{\rm \pi}^2 a_\mathrm{min}^2}, \quad a_\mathrm{min} = \sqrt{\frac{\bar{S}}{\rm \pi}}, \quad \bar{S} = \frac{1}{2{\rm \pi}}\int_0^{2{\rm \pi}}\,\mathrm{d}\zeta\, S(\zeta), \end{equation}

with $S(\zeta )$ the cross-sectional area of the boundary at fixed toroidal angle $\zeta$, and $\mathcal {V}$ the volume enclosed by the boundary.

The optimisation is performed over a set of harmonics $R_{mn}$ and $Z_{mn}$ describing the boundary, here assuming stellarator symmetry

(2.5a)\begin{gather} R(\theta, \zeta) = \sum_{m, n} R_{mn} \cos\left( m\theta - n N_\mathrm{fp} \zeta \right), \end{gather}
(2.5b)\begin{gather}Z(\theta, \zeta) = \sum_{m, n} Z_{mn} \sin\left( m\theta - n N_\mathrm{fp} \zeta \right), \end{gather}

with $N_\mathrm {fp}$ the number of field periods of the device, leading to a discrete symmetry under the transformation $\zeta \rightarrow \zeta + 2{\rm \pi} /N_\mathrm {fp}$. In the optimisation, the series in (2.5) are truncated at some chosen $m_\mathrm {max}$ and $n_\mathrm {max}$. Furthermore, $R_{00}=1$ is kept constant to set the length scale of the problem, as the objective (2.1) is invariant under a uniform rescaling of all $R_{mn}$ and $Z_{mn}$ coefficients.

The derivatives of the QS and rotational transform targets with respect to the boundary coefficients is obtained using adjoint methods. The corresponding adjoint equations and shape gradient were originally derived in Nies et al. (Reference Nies, Paul, Hudson and Bhattacharjee2022), and are here slightly modified to make the QS objective dimensionless, as shown in Appendix A. The derivative of the aspect ratio target (2.3) is derived in Appendix B.

Each optimisation begins with small values of $m_\mathrm {max}$ and $n_\mathrm {max}$ before gradually increasing them in optimisation ‘stages’, letting the optimisation run in each such stage until progress has halted. The Fourier space resolution of the solutions to the Laplace equation for the vacuum magnetic field, to the straight-field-line equation, and to their adjoint equations, is always higher than the boundary resolution and is also increased at each stage. We found it generally optimal to start the optimisation with a single poloidal harmonic $m_\mathrm {max}=1$, but with a larger number of toroidal harmonics, e.g. $n_\mathrm {max} = 3$. This may be due to the connection between the axis shape and QS (Rodríguez, Sengupta & Bhattacharjee Reference Rodríguez, Sengupta and Bhattacharjee2022b). Owing to the use of adjoint methods, a large parameter space in the boundary representation (2.5) may be explored, such that our optimisations typically involve a dozen stages or more, incrementally going up to $m_\mathrm {max}\sim n_\mathrm {max} \sim 10-20$. Despite the high resolution, a typical optimisation requires only $\mathcal {O}(10^2)$ CPU hours, with $\mathcal {O}(10^3)$ iterations of the optimiser over the multiple stages. All optimisations are performed using the Broyden-Fletcher-Goldfarb-Shanno  algorithm (Fletcher Reference Fletcher1987) implemented in the scipy package (Virtanen et al. Reference Virtanen, Gommers, Oliphant, Haberland, Reddy, Cournapeau, Burovski, Peterson, Weckesser and Bright2020), with the initial boundary taken to be that of a simple rotating ellipse. The Laplace equation for the vacuum magnetic field, as well as the corresponding adjoint problem (Nies et al. Reference Nies, Paul, Hudson and Bhattacharjee2022), are solved with the SPEC code (Hudson et al. Reference Hudson, Dewar, Dennis, Hole, McGann, von Nessi and Lazerson2012; Qu et al. Reference Qu, Pfefferlé, Hudson, Baillod, Kumar, Dewar and Hole2020). The rotational transform on the boundary is obtained by solving the straight-field-line equation $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\alpha = 0$, with the field-line label $\alpha = \theta - \iota \zeta + \lambda (\theta, \zeta )$, where $\lambda$ is a single-valued function of the angular coordinates.

After an optimised vacuum field boundary has been obtained, it may then be input to VMEC (Hirshman et al. Reference Hirshman, van Rij and Merkel1986) or DESC (Dudt & Kolemen Reference Dudt and Kolemen2020; Conlin et al. Reference Conlin, Dudt, Panici and Kolemen2023; Dudt et al. Reference Dudt, Conlin, Panici and Kolemen2023; Panici et al. Reference Panici, Conlin, Dudt, Unalmis and Kolemen2023), solving MHS force balance while setting the pressure and current profiles to vanish. This procedure yields a good approximation to the vacuum field, provided the flux surfaces are nested. This is generally the case for the QH configurations, as well as the QA configurations with $N_\mathrm {fp}=2$, but not for those with $N_\mathrm {fp}=3$ due to their larger magnetic shear, see §§ 3 and 4. The DESC solution is used only for post-processing, to evaluate volume properties and the magnetic axis. The VMEC code is used in conjunction with the SIMSOPT (Landreman et al. Reference Landreman, Medasani, Wechsung, Giuliani, Jorge and Zhu2021a) optimisation suite in § 5, to further optimise two configurations.

Although $f_\mathrm {QS}^\star$ from (A4) is used in the optimisation, for post-processing the degree of QS is also evaluated using the infinity norm of the symmetry-breaking modes

(2.6)\begin{equation} {\left|B_{mn}\right|}_\infty \equiv \max_{n/m\neq N/M} B_{mn}, \end{equation}

which may be evaluated on the boundary from either the SPEC or the DESC solution. Here, $B_{mn}$ are the coefficients of the field strength when Fourier expanded in the Boozer angles $\theta _B$ and $\zeta _B$, i.e. $B = \sum B_{mn} \cos (m\theta _B - n N_\mathrm {fp} \zeta _B)$. For the SPEC vacuum field solution, the Boozer coordinate transformation required to evaluate (2.6) follows the procedure presented in Appendix C. From the DESC solution, we further define a metric of QS in the volume as

(2.7)\begin{equation} \left\langle \frac{{\left|B_{mn}\right|}_\infty}{B_{00}} \right\rangle = \int_0^1\,\mathrm{d}s\, \frac{{\left|B_{mn}\right|}_\infty}{B_{00}}, \end{equation}

where $s$ is a normalised flux coordinate ranging from $s=0$ on the magnetic axis to $s=1$ on the boundary. We note that, although perfect QS would result in both $f_\mathrm {QS}^\star =0$ and ${\left |B_{mn}\right |}_\infty = 0$, the two measures of QS employed in this study do not perfectly correlate as QS is approached (Rodríguez, Paul & Bhattacharjee Reference Rodríguez, Paul and Bhattacharjee2022a), such that larger ${\left |B_{mn}\right |}_\infty$ values than expected might be obtained in the configurations optimised for low $f_\mathrm {QS}^\star$. Moreover, we observe differences in the level of QS on the boundary given by the DESC and vacuum field solutions, see §§ 3 and 4. These may be due to the different algorithms employed, or small inaccuracies in the equilibrium solutions becoming important as very small QS deviations are reached.

3 Aspect ratio scan

We first investigate how the choice of aspect ratio $A$ affects the QS optimisation. We consider QA configurations with $N_\mathrm {fp}=2$ and $N_\mathrm {fp}=3$, with a boundary rotational transform target $\iota _T = 0.42$, as had been used for the precise QA configuration obtained by Landreman & Paul (Reference Landreman and Paul2022). Optimisations for QA at $N_\mathrm {fp}=1$ and $N_\mathrm {fp}=4$ generally performed poorly, consistent with previous studies (Landreman Reference Landreman2022), and are thus not shown here. For the QH configurations, $N_\mathrm {fp}=3$ and $N_\mathrm {fp}=4$ are considered, with rotational transform targets of $\iota _T \approx 1.3$. Although it is crucial only for the QA to prescribe a rotational transform target, as the optimisation would otherwise tend towards an axisymmetric solution with $\iota =0$, we also included an $\iota$ target for the QH to make the comparison between different aspect ratio values more meaningful. We do not include here QH configurations with either $N_\mathrm {fp}= 2$, due to their larger deviations from QS in our optimisation, or with $N_\mathrm {fp}\geq 5$, as these could only be optimised at larger aspect ratio values.

For the QA configuration with $N_\mathrm {fp}=2$ and $A=6$, the optimisation was pushed to very high resolutions to reach record levels of QS on the boundary. This optimisation was performed as a proof of principle that the QS error on the boundary could be reduced substantially, although a large number of boundary harmonics are required in practice, and an increased computational cost is incurred. We note that reducing the QS error from ${\left |B_{mn}\right |}_\infty /B_{00} \sim 10^{-5}$ to ${\left |B_{mn}\right |}_\infty /B_{00}\sim 10^{-10}$ does not substantially change the boundary shape, or the average QS error in the volume.

The contours of the magnetic field strength on the boundary of the configurations with the lowest aspect ratio obtained for each case are displayed in figure 1. Configurations with good QA (i.e. those for which the boundary magnetic field strength contours look approximately straight in Boozer coordinates) could be found down to $A=2.6$ and $A=4.5$ for $N_\mathrm {fp}=2$ and $N_\mathrm {fp}=3$, respectively. For the QH configurations, the aspect ratio could be reduced to $A=3.6$ and $A=3$ for $N_\mathrm {fp} = 3$ and $N_\mathrm {fp}=4$, respectively. As attested by the contours in Boozer coordinates (which are perfectly straight in the limit of exact QS), these configurations have excellent QS on the boundary. Only the low aspect ratio $N_\mathrm {fp}=4$ QH configuration has visible deviations of the contours from straight lines. For the QA configurations, the three-dimensional shaping is visibly strongest on the inboard side, in agreement with the analytical results of Plunk & Helander (Reference Plunk and Helander2018) for low aspect ratio QA stellarators close to axisymmetry.

Figure 1. Contours of the magnetic field strength $B$ on the boundary of quasisymmetric stellarators at small aspect ratio $A$: top view (a–d), side view (e–h) and in Boozer coordinates (i–l).

The QS error and rotational transform profiles for all configurations in the aspect ratio scan are shown in figure 2 as a function of the normalised flux $s$. Further properties of the configurations are listed in table 1. As expected, the QS is best on the boundary, down to record normalised values of ${\left |B_{mn}\right |}_\infty /B_{00} \sim 10^{-8}$ in the DESC solution, and ${\left |B_{mn}\right |}_\infty /B_{00} \sim 10^{-10}$ in the SPEC solution. The QS error also remains low throughout the volume in most cases, ${\left |B_{mn}\right |}_\infty /B_{00} \lesssim 10^{-2}$ at small aspect ratio, and ${\left |B_{mn}\right |}_\infty /B_{00} \lesssim 10^{-3}$ at large aspect ratio. The good QS levels throughout the volume suggest that, in practice, the optimisation of QS on a boundary is achieved by approximating a globally QS configuration, see § 5.1.

Figure 2. Volume properties of stellarators optimised for QS on the boundary, with varying aspect ratio $A$. The QS error (a) and rotational transform (b) profiles were obtained using the DESC code. The QS error reaches record levels on the boundary ($s=1$), and generally remains relatively low in the core. The rotational transform profiles are flat for all configurations except for $N_\mathrm {fp}=3$ QA.

Table 1. Properties of stellarators optimised for QS at varying aspect ratio $A$: QS figure of merit $f_\mathrm {QS}^\star$, maximum symmetry-breaking mode on boundary from vacuum solution (SPEC) and MHS solution (DESC), volume-averaged QS error in MHS solution, difference between $\iota$ on axis and at the edge $\Delta \iota = \iota (s=1)-\iota (s=0)$ and minimum value of axis curvature $\kappa _\mathrm {min}$.

The QA $N_\mathrm {fp}=2$ and $N_\mathrm {fp}=3,4$ QH configurations all have small magnetic shear, as attested by the flatness of the rotational transform profiles, in accordance with previous QS optimisations of vacuum fields (Landreman & Paul Reference Landreman and Paul2022). However, the $N_\mathrm {fp}=3$ QA stellarators attain substantial magnetic shear, with $\iota$ varying from $0.42$ on the boundary up to $0.75$ on axis. In those cases, the magnetic shear $\mathrm {d}\iota /\mathrm {d}\psi$ is positive (tokamak-like) and approximately constant, such that $\Delta \iota = \iota (s=1)-\iota (s=0) \sim \mathrm {d}\iota /\mathrm {d}\psi \;\Delta \psi \propto 1/A^2$, as the change in flux $\Delta \psi \sim B r^2 \sim (B R^2) / A^2$.

The substantial magnetic shear of the $N_\mathrm {fp}=3$ QA configurations means the rotational transform generally crosses low-order rational values in the core, leading to island chains and chaotic regions. In those cases, the DESC solution may not be trustworthy, as it assumes nested flux surfaces a priori. However, the integrability may be improved in an auxiliary optimisation, as shown in § 5.2. We note that large chaotic regions are found for the lowest aspect ratio $N_\mathrm {fp}=3$ QA with $A=4.5$. One may hypothesise that the aspect ratio is here limited by the chaotic region increasing in size until it encounters the boundary.

4 Rotational transform scan

We now vary the boundary rotational transform target $\iota _T$ in the QS optimisation. We again consider $N_\mathrm {fp}=2$ and $N_\mathrm {fp}=3$ QA, with a fixed aspect ratio of $A=6$. We also optimise $N_\mathrm {fp}=4$ QH stellarators at a higher aspect ratio $A=8$. The aspect ratio targets were chosen to be identical to the precise $N_\mathrm {fp}=2$ QA and $N_\mathrm {fp}=4$ QH configurations of Landreman & Paul (Reference Landreman and Paul2022). Good QS could be obtained for rotational transform values up to $\iota _e \sim 0.7$ and $\iota _e \sim 0.8$ for the $N_\mathrm {fp}=2$ and $N_\mathrm {fp}=3$ QA configurations, respectively, and in the range $\iota _e \sim 1$ to $\iota _e \sim 2$ for the $N_\mathrm {fp}=4$ QH.

The contours of the magnetic field strength for the QA configurations with highest rotational transform obtained ($\iota _e = 0.65$ and $\iota _e = 0.82$ for $N_\mathrm {fp}=2$ and $N_\mathrm {fp}=3$, respectively) are shown in figure 3, alongside the contours for the QH with the lowest and highest rotational transform values obtained ($\iota _e=1.12$ and $\iota _e=1.97$). The contours of the field strength in Boozer coordinates are straight to the naked eye, attesting to the high level of QS.

Figure 3. Contours of the magnetic field strength $B$ on the boundary of quasisymmetric stellarators at varying edge rotational transform $\iota _e$: top view (a–d), side view (e–h) and in Boozer coordinates (i–l).

The QS error in the volume and the rotational transform profile for all configurations obtained in the rotational transform scan are shown in figure 4. Further properties of the configurations are listed in table 2. The QS error is again very low at the edge in most configurations, down to ${\left |B_{mn}\right |}_\infty /B_{00} \sim 10^{-8}$. The QS error in the core also remains low, with ${\left |B_{mn}\right |}_\infty /B_{00} \lesssim 10^{-3}$ for most QA configurations and the high $\iota$ QH, while ${\left |B_{mn}\right |}_\infty /B_{00} \lesssim 10^{-2}$ for the high $\iota$ $N_\mathrm {fp}=3$ QA and lower $\iota$ QH. The magnetic shear is again found to be small for the $N_\mathrm {fp}=2$ QA, and for the $N_\mathrm {fp}=4$ QH configurations at lower $\iota _e \lesssim 1.75$. Substantial magnetic shear is found for the $N_\mathrm {fp}=3$ QA configurations, with a nearly constant difference between the transform in the core and in the edge, $\Delta \iota =\iota (s=1)-\iota (s=0)\approx 0.2$. Noticeable magnetic shear is also found for the highest $\iota _e=1.97$ QH configuration.

Figure 4. Volume properties of stellarators optimised for QS on the boundary, with varying edge rotational transform $\iota _e$. The QS error (a) and rotational transform (b) profiles were obtained using the DESC code. The QS error reaches record levels on the boundary ($s=1$), and generally remains relatively low in the core. The rotational transform profiles are flat for all configurations except for $N_\mathrm {fp}=3$ QA and the high $\iota _e=1.97$ QH.

Table 2. Properties of stellarators optimised for QS at varying edge rotational transform $\iota _e$: QS figure of merit $f_\mathrm {QS}^\star$, maximum symmetry-breaking mode on boundary from vacuum solution (SPEC) and MHS solution (DESC), volume-averaged QS error in MHS solution, difference between $\iota$ on axis and at the edge $\Delta \iota = \iota (s=1)-\iota (s=0)$ and minimum value of axis curvature $\kappa _\mathrm {min}$.

The magnetic axes of the configurations at varying rotational transform are shown in figure 5. For the QA configurations, as the rotational transform is increased, the axis develops regions of low curvature, as may also be seen from the straightening of the axis, and of high torsion. The straightening of the magnetic axis with increasing $\iota$ is in agreement with a near-axis model of QS (Rodríguez et al. Reference Rodríguez, Sengupta and Bhattacharjee2022b) that showed how higher values of $\iota$ move QA configurations towards a phase transition where the configuration would be quasi-isodynamic, which requires the axis to have regions of vanishing curvature. For the QH configurations, the phase transition is approached as $\iota$ is decreased, corresponding again to a straightening of the magnetic axis at lower $\iota$ in figure 5, although it is less pronounced than for the QA cases.

Figure 5. Magnetic axes’ top view (a), curvature (b) and torsion (c), of the stellarators optimised for QS on the boundary with varying edge rotational transform $\iota _e$. The DESC solution is used to obtain the magnetic axis. Regions of reduced curvature develop at high $\iota$ for the QA, and at low $\iota$ for the QH configurations, as predicted by Rodríguez et al. (Reference Rodríguez, Sengupta and Bhattacharjee2022b). The legend for the various colours and shadings is given in figure 4.

5 Optimisation refinement for volume properties

We here use the SIMSOPT (Landreman et al. Reference Landreman, Medasani, Wechsung, Giuliani, Jorge and Zhu2021a) optimisation suite in conjunction with the VMEC code Hirshman et al. (Reference Hirshman, van Rij and Merkel1986) to refine the optimisation of two configurations with good QS on the boundary. We first demonstrate improvement of QS in the volume in § 5.1, and then show how a $N_\mathrm {fp}=3$ QA can be made integrable in § 5.2. We note that the VMEC code struggles to obtain a converged equilibrium for the QH configurations at high resolutions, such that only QA configurations were considered.

5.1 Volume quasisymmetry

We here optimise for QS in the volume of the $N_\mathrm {fp}=2$ QA configuration with aspect ratio $A=6$ and high rotational transform $\iota = 0.65$. The optimisation follows the procedure of Landreman & Paul (Reference Landreman and Paul2022). Because of the large number of harmonics employed in the adjoint optimisation ($n_\mathrm {max} = 11$ and $m_\mathrm {max} = 7$ for the $\iota _T = 0.65$ case under consideration), only a small number of steps in the SIMSOPT optimisation are computationally feasible due to the use of finite differences to evaluate the objective function's derivative.

After only $15$ iterations of the optimiser, the QS in the volume is reduced by approximately an order of magnitude, although the QS on the boundary is degraded, as shown in figure 6. The QS error is now highest on the boundary, similar to previous optimisation results (e.g. Landreman & Paul Reference Landreman and Paul2022). The improvement of the QS in the volume does not require large changes in the boundary shape, as demonstrated by the small changes to the boundary cross-sections. This further supports the hypothesis that the stellarators optimised for boundary QS are close to solutions with good volume QS. Finally, we note that the already small magnetic shear is further reduced by the global QA optimisation, decreasing from $\Delta \iota = \iota (s=1)-\iota (s=0)=0.012$ to $\Delta \iota = 0.002$.

Figure 6. Two field-period QA with large edge rotational transform $\iota _e = 0.65$ before (Boundary QA) and after (Global QA) auxiliary optimisation for QS in the volume. The QS error (a), evaluated using the VMEC code, shows an improvement in the global QS level, at the cost of higher QS error on the boundary. The optimisation did not require a large modification to the stellarator shape, as attested by the cross-sections of the boundary (b), shown for $\zeta /(2{\rm \pi} / 5N_\mathrm {fp}) \in \{0,1,2,3,4\}$.

5.2 Flux surface nestedness

Generally, the QA configurations with $N_\mathrm {fp}=3$ are not integrable, i.e. they do not possess nested magnetic flux surfaces, as the magnetic shear is large enough to cause the $\iota$ profile to cross low-order rational values. Depending on the configuration at hand, this leads to magnetic island chains, or stochastic regions. In contrast, the QH and $N_\mathrm {fp}=2$ QA configurations all have nested flux surfaces, as they have small magnetic shear and the targeted $\iota$ was chosen so as to avoid low-order rationals.

Even for the non-integrable $N_\mathrm {fp}=3$ QA configurations, the integrability may be improved in a subsequent optimisation, as shown here for the $N_\mathrm {fp}=3$ QA configuration with high edge rotational transform $\iota _e=0.72$ and aspect ratio $A=6$. The integrability is targeted by minimising Greene's residue (Greene Reference Greene1979) on a set of rational values of $\iota$ seen to cause islands in the Poincaré plot, as in Landreman, Medasani & Zhu (Reference Landreman, Medasani and Zhu2021b). An objective targeting QS in the volume is further included.

The Greene residue is targeted for the $\iota = 3/4$ and $\iota = 6/7$ resonances. Although the boundary coefficients go up to $m_\mathrm {max}=n_\mathrm {max}=10$, the optimisation space is limited to a small number of coefficients in the boundary representation (2.5), up to $m = 2$ and $n=4$ to make the optimisation computationally tractable.

Before the integrability optimisation, the Poincaré plots clearly reveal the presence of multiple magnetic island chains, see figure 7(a). These are removed by the optimisation, as shown in figure 7(b). Furthermore, the QS error in the core could be reduced, although the QS error on the boundary is degraded, see figure 7(c). The changes in the QS and integrability do not require a substantial modification to the boundary shape, as demonstrated by the cross-sections in figure 7(d).

Figure 7. Optimisation for integrability and volume QS of $N_\mathrm {fp}=3$ QA with high edge rotational transform $\iota _e=0.72$. The integrability optimisation leads to the disappearance of the magnetic island chains, such that the final configuration has optimised configuration. Furthermore, the volume QS error (evaluated using the VMEC code) is also reduced. The required change in the stellarator shape is minimal, as attested by the cross-sections of the boundary at $\zeta /(2{\rm \pi} /5N_\mathrm {fp}) \in \{0,1,2,3,4\}$. (a) Poincaré before integrability optimisation. (b) Poincaré after integrability optimisation. (c) Quasisymmetry error. (d) Boundary cross-sections at varying $\zeta$.

6 Conclusions

We leveraged the computational efficiency of adjoint methods to optimise for QS on the boundary of stellarator vacuum fields, reaching record levels of QS on the boundary. For a $N_\mathrm {fp}=2$ quasi-axisymmetric configuration with $A=6$ and $\iota _e$ similar to (Landreman & Paul Reference Landreman and Paul2022), where a volume symmetry-breaking mode amplitude of ${\left |B_{mn}\right |}_\infty /B_{00} \approx 3 \times 10^{-5}$ was achieved, we could optimise QS to ${\left |B_{mn}\right |}_\infty /B_{00} \approx 2 \times 10^{-10}$ on the boundary. The configurations optimised for boundary QS appear close to solutions with precise global QS, such that the boundary QS optimisation may be followed by a global QS optimisation, the latter requiring only small changes to the boundary, see § 5.1.

We were also able to obtain QA configurations at $N_\mathrm {fp}=3$, which possess substantial magnetic shear compared with the QH and $N_\mathrm {fp}=2$ QA configurations (see also Landreman & Paul Reference Landreman and Paul2022). An increased magnetic shear could be beneficial, e.g. to reduce the turbulent transport. However, the larger magnetic shear generally leads to the rotational transform crossing resonances, causing magnetic island chains and stochastic regions. We showed in § 5.2 how these may be removed in an auxiliary optimisation, restoring flux surface nestedness.

We studied in § 3 the range of aspect ratio values for which good QS could be obtained, leading to designs with small aspect ratios similar to those of tokamaks, e.g. an $A=2.6$ QA with $N_\mathrm {fp}=2$, and QH configurations with $A=3.6$ for $N_\mathrm {fp}=3$, and $A=3$ for $N_\mathrm {fp}=4$. Reaching such low aspect ratios could prove crucial for a stellarator fusion power plant, as a large aspect ratio might require a prohibitively large major radius for a prescribed volume target. However, configurations with tight aspect ratios have degraded levels of QS and generally have a large amount of shaping, for which economical coils could be precluded.

We further investigated in § 4 the range of $\iota$ values compatible with good QS on the boundary. Quasiaxisymmetric configurations could be found up to $\iota \sim 0.7-0.8$, while the QH configurations at $N_\mathrm {fp}=4$ had a range of $\iota$ between $\sim 1$ and $2$. The value of $\iota$ is crucial in many ways, e.g. to avoid prompt losses of energetic particles (Paul et al. Reference Paul, Bhattacharjee, Landreman, Alex, Velasco and Nies2022), as the orbit width scales inversely with ${\left |\iota -N/M\right |}$. For the QA configurations, good QS is most easily obtained at low rotational transform values, as these are closer to the axisymmetric case. For the QH configurations, good QS is most readily obtained at intermediate values of $\iota _e \sim 1.5$, though the high $\iota _e \sim 2$ configurations might prove interesting due to their increased magnetic shear.

We note that, although the inclusion of thermal pressure and plasma current can lead to significant changes to the QS solutions (e.g. Landreman et al. Reference Landreman, Buller and Drevlak2022), vacuum configurations can provide useful guidance for finite-pressure equilibria (Boozer Reference Boozer2019). Furthermore, although the present study was mostly limited to optimisation of QS on the boundary, this might be sufficient to guarantee good confinement by effectively creating an edge transport barrier (for the guiding-centre trajectories). Moreover, reducing the confinement in the plasma core may be desirable to avoid impurity and ash accumulation. If global QS is desired, however, the configurations presented generally possess good levels of QS even in the core, which may be further improved by a subsequent global QS optimisation, see § 5.1.

In this study, we focused on the compatibility of QS with aspect ratio and rotational transform targets. Future work ought to consider the trade-offs between QS and other targets typically considered in the design of a stellarator, such as MHD stability, turbulent transport optimisation or engineering feasibility. Furthermore, adjoint-based fixed-boundary optimisation could be combined with coil optimisation for an efficient derivative-based single-stage approach simultaneously optimising the plasma and the coils, as studied by Jorge et al. (Reference Jorge, Goodman, Landreman, Rodrigues and Wechsung2023).

Acknowledgements

We thank the anonymous referees for their comments and suggestions. R.N. thanks E. Rodriguez, S. Buller and M. Landreman for helpful discussions.

Editor P. Helander thanks the referees for their advice in evaluating this article.

Funding

This work was supported by the U.S. DOE (grants DE-AC02-09CH11466, DE-SC0016072, DE-AC02-76CH03073, DE-SC0024548, DE-SC0022005, Field Work Proposal No. 19); and by the Simons Foundation/SFARI (A.B., 560651).

Declaration of interests

The authors report no conflict of interest.

Data availability statement

All configurations presented in this study may be found at DOI:10.34770/rt5a-sy08 (Nies et al. Reference Nies, Paul, Panici, Hudson and Bhattacharjee2024).

Appendix A. Normalised quasisymmetry objective

We introduce two normalisations in our QS objective, compared with Nies et al. (Reference Nies, Paul, Hudson and Bhattacharjee2022). As a reminder, our objective is based on the definition of QS as

(A1)\begin{equation} \frac{\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla}\psi\times\boldsymbol{\nabla} B}{\boldsymbol{B}\boldsymbol{\cdot}\boldsymbol{\nabla} B} ={-}\frac{MG + NI}{N-\iota M}, \end{equation}

with the magnetic field $\boldsymbol {B}$, the toroidal flux $\psi$ and the enclosed poloidal and toroidal currents, $G$ and $I$, respectively. We further consider a vacuum magnetic field

(A2)\begin{equation} \boldsymbol{B}=G\boldsymbol{\nabla}(\zeta+\omega) \equiv G \boldsymbol{\breve{B}}, \end{equation}

with $I=0$, $\omega$ a single-valued scalar function and $\mathcal {S}$ the stellarator boundary. The QS objective in Nies et al. (Reference Nies, Paul, Hudson and Bhattacharjee2022) was defined as

(A3)\begin{equation} f_\mathrm{QS} = \int_{\mathcal{S}} \,\mathrm{d}S\; \frac{1}{2}v_\mathrm{QS}^2, \quad \text{with}\ v_\mathrm{QS} = \boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla} \breve{B} - \boldsymbol{\breve{B}}\times\boldsymbol{\breve{g}}_\psi\boldsymbol{\cdot}\boldsymbol{\nabla} \breve{B} \left( \iota - N/M\right), \end{equation}

where $\boldsymbol {\breve {B}}=\boldsymbol {B}/G$, and $\boldsymbol {\breve {g}}_\psi$ is a vector defined on the boundary that reduces to $\boldsymbol {\nabla }\psi /G$ in the limit of an integrable magnetic field, see Nies et al. (Reference Nies, Paul, Hudson and Bhattacharjee2022). In this study, we consider the objective

(A4)\begin{equation} f_\mathrm{QS}^\star{=} \sqrt{\int_{\mathcal{S}} \,\mathrm{d}S\; w_\mathrm{QS}^2} \quad \text{with} \ w_\mathrm{QS} = v_\mathrm{QS} / \breve{B}^2. \end{equation}

The normalisation to $w_\mathrm {QS}$ makes the objective function dimensionless (as $w_\mathrm {QS}\sim 1/L$), and should lend itself better to QP optimisation, as $f_\mathrm {QS}$ is dominated by contributions from high-field regions and QP configurations typically exhibit significant variation of $B$ on the surface. Finally, the square root is motivated by the fact that the Fourier expansion of $w_\mathrm {QS}$ in Boozer coordinates (Rodríguez et al. Reference Rodríguez, Paul and Bhattacharjee2022a) is

(A5)\begin{equation} w_\mathrm{QS} ={-}{\rm i} \frac{1}{G} \sum_{n,m} \left(n-m N/M\right) B_{mn} \exp({{\rm i}(m\theta_B-n\zeta_B)}), \end{equation}

so that we can expect our objective to scale linearly in the size of the symmetry-breaking Fourier components.

We now compute the shape derivative of the new QS objective (A4). First, note that

(A6)\begin{equation} \delta f_\mathrm{QS}^\star{=} \frac{1}{f_\mathrm{QS}^\star} \delta\left( \int_{\mathcal{S}} \,\mathrm{d}S\; \frac{1}{2} w_\mathrm{QS}^2 \right), \end{equation}

and

(A7)\begin{align} & \delta\left( \int \,\mathrm{d}S\; \frac{1}{2} w_\mathrm{QS}^2 \right) = \int_{\mathcal{S}} \,\mathrm{d}S\; \left[ w_\mathrm{QS} \; \delta w_\mathrm{QS} + \frac{1}{2}\left({\delta\boldsymbol{{x}}\boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}}}\right) ({\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla} + h)w_\mathrm{QS}^2 \right] \end{align}
(A8)\begin{align} & \quad = \int_{\mathcal{S}} \,\mathrm{d}S\; \left\{ \frac{v_\mathrm{QS}}{\breve{B}^4}\delta v_\mathrm{QS} + \delta\omega \boldsymbol{\nabla}_\varGamma\boldsymbol{\cdot} \left( 2 \frac{v_\mathrm{QS}^2}{\breve{B}^6} \boldsymbol{\breve{B}} \right) + ({\delta\boldsymbol{{x}}\boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}}}) \left[ \frac{v_\mathrm{QS}}{\breve{B}^4} {\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla} v_\mathrm{QS} - 2 \frac{v_\mathrm{QS}^2}{\breve{B}^5}{\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\breve{B} + \frac{h}{2}\frac{v_\mathrm{QS}^2}{\breve{B}^4} \right] \right\}, \end{align}

where we used $\delta \breve {B} = \boldsymbol {\breve {B}}\boldsymbol {\cdot }\boldsymbol {\nabla }(\delta \omega )/\breve {B}$ and partially integrated the second term. Here, $h$ is the summed curvature, and ${\boldsymbol {\hat {{n}}}}$ is a unit vector normal to the surface. We can now proceed entirely analogously to the derivation of $\delta f_\mathrm {QS}$ in (B24) of Nies et al. (Reference Nies, Paul, Hudson and Bhattacharjee2022), normalising by $\breve {B}^{-4}$ where $v_\mathrm {QS}$ appears, leading to

(A9)\begin{align} \delta f_\mathrm{QS}^\star & = \frac{1}{f_\mathrm{QS}^\star}\int_{\mathcal{S}} \,\mathrm{d} S\; \left\{ \delta\omega \;\boldsymbol{\nabla}_\varGamma\boldsymbol{\cdot} \left[ -\frac{v_\mathrm{QS}}{\breve{B}^4} \boldsymbol{\nabla}_\varGamma \breve{B} + \frac{\boldsymbol{\breve{B}}}{\breve{B}} \boldsymbol{\nabla}_\varGamma\boldsymbol{\cdot}\left(\frac{v_\mathrm{QS}}{\breve{B}^4} \boldsymbol{\breve{B}}\right) \right.\right.\nonumber\\ & \left.\quad + \left(\iota-N/M\right) \left( \frac{v_\mathrm{QS}}{\breve{B}^4} \; \boldsymbol{\breve{g}}_\psi\times\boldsymbol{\nabla}_\varGamma\breve{B} - \boldsymbol{\breve{B}}\; \boldsymbol{\nabla}_\varGamma\boldsymbol{\cdot}\left(\frac{1}{\breve{B}} \frac{v_\mathrm{QS}}{\breve{B}^4} \boldsymbol{\breve{B}}\times\boldsymbol{\breve{g}}_\psi\right) \right) + 2\frac{v_\mathrm{QS}^2}{\breve{B}^6} \boldsymbol{\breve{B}} \right]\nonumber\\ & \quad - \delta\iota \, \frac{v_\mathrm{QS}}{\breve{B}^4} \boldsymbol{\breve{B}}\times\boldsymbol{\breve{g}}_\psi\boldsymbol{\cdot}\boldsymbol{\nabla}_\varGamma\breve{B} \left[ \left(\iota-N/M\right) \frac{\boldsymbol{\nabla}_\varGamma\alpha\boldsymbol{\cdot}\boldsymbol{\nabla}_\varGamma\zeta}{{\left|\boldsymbol{\nabla}_\varGamma\alpha\right|}^2} + 1\right]\nonumber\\ & \quad - \delta\lambda \left(\iota-N/M\right)\boldsymbol{\nabla}_\varGamma\boldsymbol{\cdot}\left[ \frac{v_\mathrm{QS}}{\breve{B}^4} \frac{\boldsymbol{\nabla}_\varGamma\alpha}{{\left|\boldsymbol{\nabla}_\varGamma\alpha\right|}^2} \, \boldsymbol{\breve{B}}\times\boldsymbol{\breve{g}}_\psi\boldsymbol{\cdot}\boldsymbol{\nabla}_\varGamma\breve{B} \right]\nonumber\\ & \quad + ({\delta\boldsymbol{{x}}\boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}}}) \left[ \left(\iota-N/M\right) {\left|\boldsymbol{\breve{g}}_\psi\right|} \, \boldsymbol{\breve{B}}\times\boldsymbol{\nabla}\breve{B}\boldsymbol{\cdot} \left( \boldsymbol{\nabla}_\varGamma \left( \frac{v_\mathrm{QS}}{\breve{B}^4} \right) - \frac{v_\mathrm{QS}}{\breve{B}^4} \frac{\boldsymbol{\nabla}_\varGamma {\left|\boldsymbol{\nabla}_\varGamma\alpha\right|}}{{\left|\boldsymbol{\nabla}_\varGamma\alpha\right|}} \right)\right.\nonumber\\ & \quad + \frac{h}{2} \frac{v_\mathrm{QS}^2}{\breve{B}^4} - ({\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\breve{B}) \boldsymbol{\nabla}_\varGamma\boldsymbol{\cdot}\left( \frac{v_\mathrm{QS}}{\breve{B}^4} \boldsymbol{\breve{B}} \right) - \frac{v_\mathrm{QS}}{\breve{B}^4}\;\left(\boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\hat{{n}}}} - {\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\breve{B}} \right)\boldsymbol{\cdot}\boldsymbol{\nabla}_\varGamma\breve{B}\nonumber\\ & \quad \left.\left. - \frac{v_\mathrm{QS}}{\breve{B}^4} \, \boldsymbol{\breve{B}}\times\boldsymbol{\breve{g}}_\psi\boldsymbol{\cdot}\boldsymbol{\nabla} \breve{B} \left( \iota - N/M\right) \left( \frac{\boldsymbol{\nabla}_\varGamma\alpha\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_\varGamma\alpha}{{\left|\boldsymbol{\nabla}_\varGamma \alpha\right|}^2} -h \right) - 2 \frac{v_\mathrm{QS}^2}{\breve{B}^5}{\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\breve{B} \right] \right\}, \end{align}

where $\boldsymbol {\nabla }_\varGamma$ is the tangential derivative. These modifications then carry through to the adjoint equations and shape gradient, reproduced here for convenience.

First, the adjoint $q_\alpha$ to the straight-field-line equation $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\alpha =0$, with $\alpha = \theta - \iota \zeta + \lambda$ and single-valued $\lambda$, is given by

(A10a)\begin{gather} \boldsymbol{\nabla}_\varGamma \boldsymbol{\cdot} \left( q_\alpha \boldsymbol{\breve{B}} \right) ={-}\frac{1}{f_\mathrm{QS}^\star} \boldsymbol{\nabla}_\varGamma\boldsymbol{\cdot}\left[ \boldsymbol{\nabla}_\varGamma\alpha \left( \frac{v_\mathrm{QS}}{\breve{B}^4}\, \boldsymbol{\breve{B}}\times\boldsymbol{\breve{g}}_\psi\boldsymbol{\cdot}\boldsymbol{\nabla} \breve{B}\, \frac{\iota - N/M}{{\left|\boldsymbol{\nabla}_\varGamma \alpha\right|}^2}\right) \right] , \end{gather}
(A10b)\begin{gather}0=\int_{\mathcal{S}} \,\mathrm{d} S \left\{ q_\alpha \boldsymbol{\breve{B}} \boldsymbol{\cdot} \boldsymbol{\nabla}\zeta + \frac{v_\mathrm{QS}}{f_\mathrm{QS}^\star\breve{B}^4}\, \boldsymbol{\breve{B}}\times\boldsymbol{\breve{g}}_\psi\boldsymbol{\cdot}\boldsymbol{\nabla} \breve{B} \left[ \frac{\boldsymbol{\nabla}_\varGamma\alpha\boldsymbol{\cdot}\boldsymbol{\nabla}_\varGamma \zeta}{{\left|\boldsymbol{\nabla}_\varGamma \alpha\right|}^2} (\iota-N/M) + 1 \right] \right\} . \end{gather}

Second, the adjoint $q_\omega$ to the Laplace equation $\Delta \omega = 0$ may be written as

(A11a)\begin{gather} \Delta q_\omega = 0,\quad \mathrm{in}\ \mathcal{V}, \end{gather}
(A11b)\begin{gather}\boldsymbol{\nabla} q_\omega \boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}} ={-}\boldsymbol{\nabla}_\varGamma \boldsymbol{\cdot} \left\{ q_\alpha \boldsymbol{\nabla}_\varGamma \alpha + \frac{v_\mathrm{QS}}{f_\mathrm{QS}^\star\breve{B}^4} \;\boldsymbol{\nabla}_\varGamma \breve{B} - \frac{\boldsymbol{\breve{B}}}{f_\mathrm{QS}^\star\breve{B}} \boldsymbol{\nabla}_\varGamma \boldsymbol{\cdot} \left(\frac{v_\mathrm{QS}}{\breve{B}^4} \; \boldsymbol{\breve{B}}\right) - 2\frac{v_\mathrm{QS}^2}{f_\mathrm{QS}^\star\breve{B}^6} \boldsymbol{\breve{B}}\right.\nonumber\\ \left. + \frac{1}{f_\mathrm{QS}^\star}\left(\iota-N/M\right) \left[ \frac{v_\mathrm{QS}}{\breve{B}^4} \, \boldsymbol{\breve{g}}_\psi\times\boldsymbol{\nabla}_\varGamma\breve{B} - \boldsymbol{\breve{B}}\; \boldsymbol{\nabla}_\varGamma\boldsymbol{\cdot}\left(\frac{1}{\breve{B}}\frac{v_\mathrm{QS}}{\breve{B}^4} \boldsymbol{\breve{B}}\times\boldsymbol{\breve{g}}_\psi\right) \right] \right\}, \quad \text{on } \mathcal{S}, \end{gather}

with $\mathcal {V}$ the volume enclosed by $\mathcal {S}$.

Finally, the shape gradient follows as

(A12)\begin{align} \mathcal{G}_\mathrm{QS}^\star & ={-}\frac{1}{f_\mathrm{QS}^\star} ({\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\breve{B}) \boldsymbol{\nabla}_\varGamma\boldsymbol{\cdot}\left( \frac{v_\mathrm{QS}}{\breve{B}^4} \boldsymbol{\breve{B}} \right) -\frac{1}{f_\mathrm{QS}^\star} \frac{v_\mathrm{QS}}{\breve{B}^4}\left(\boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\hat{{n}}}} - {\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\breve{B}} \right)\boldsymbol{\cdot}\boldsymbol{\nabla}_\varGamma\breve{B} \nonumber\\ & \quad + \frac{\iota-N/M}{f_\mathrm{QS}^\star} {\left|\boldsymbol{\breve{g}}_\psi\right|} \, \boldsymbol{\breve{B}}\times\boldsymbol{\nabla}\breve{B}\boldsymbol{\cdot} \left[ {\left|\boldsymbol{\nabla}_\varGamma\alpha\right|} \boldsymbol{\nabla}_\varGamma\left( \frac{v_\mathrm{QS}}{\breve{B}^4{\left|\boldsymbol{\nabla}_\varGamma\alpha\right|}} \right) + {\boldsymbol{\hat{{n}}}}\, \frac{v_\mathrm{QS}}{\breve{B}^4} \left( \frac{\boldsymbol{\nabla}_\varGamma\alpha\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}_\varGamma\alpha}{{\left|\boldsymbol{\nabla}_\varGamma \alpha\right|}^2} -h \right) \right] \nonumber\\ & \quad + \boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla} q_\omega + q_\alpha \left( {\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{\breve{B}} - \boldsymbol{\breve{B}}\boldsymbol{\cdot}\boldsymbol{\nabla}{\boldsymbol{\hat{{n}}}} \right)\boldsymbol{\cdot}\boldsymbol{\nabla}_\varGamma\alpha + \frac{1}{f_\mathrm{QS}^\star}\frac{h}{2} \frac{v_\mathrm{QS}^2}{\breve{B}^4} - 2 \frac{1}{f_\mathrm{QS}^\star}\frac{v_\mathrm{QS}^2}{\breve{B}^5}{\boldsymbol{\hat{{n}}}}\boldsymbol{\cdot}\boldsymbol{\nabla}\breve{B}. \end{align}

Appendix B. Shape gradient for aspect ratio objective

The definition of the aspect ratio $A$ is given in (2.4). Note that the mean cross-sectional area may be written as

(B1)\begin{equation} \bar{S} = \frac{1}{2{\rm \pi}} \int\,\mathrm{d}\zeta \int_{\mathcal{S}_\zeta} \,\mathrm{d}S = \frac{1}{2{\rm \pi}}\int_\mathcal{V}\,\mathrm{d}V\, \sqrt{\frac{g_{ss}g_{\theta\theta}-g_{s\theta}^2}{g}}, \end{equation}

with the metric elements $g_{ij}$ and its determinant $g$. The shape derivative follows as (see e.g. Walker Reference Walker2015)

(B2)\begin{equation} \delta \bar{S} = \frac{1}{2{\rm \pi}}\int_\mathcal{S}\,\mathrm{d}S\, ({\delta\boldsymbol{{x}}\boldsymbol{\cdot} {\boldsymbol{\hat{{n}}}}}) \sqrt{\frac{g_{ss}g_{\theta\theta}-g_{s\theta}^2}{g}}. \end{equation}

One may then obtain the shape derivative of the aspect ratio (2.4) as

(B3)\begin{equation} \delta A = A \left( \frac{\delta \mathcal{V}}{\mathcal{V}} - \frac{3}{2}\frac{\delta \bar{S}}{\bar{S}} \right), \end{equation}

with the enclosed volume's shape derivative (see e.g. Walker Reference Walker2015) $\delta \mathcal {V} = \int _\mathcal {S}\,\mathrm {d}S\, ({\delta \boldsymbol {{x}}\boldsymbol {\cdot } {\boldsymbol {\hat {{n}}}}})$.

Appendix C. Boozer coordinate transformation for a vacuum magnetic field

The covariant coordinates of the magnetic field in Boozer coordinates $(\psi, \theta _B, \zeta _B)$ are given by

(C1)\begin{equation} \boldsymbol{B} = G(\psi) \boldsymbol{\nabla} \zeta_B + I(\psi)\boldsymbol{\nabla} \theta_B + K(\boldsymbol{r})\boldsymbol{\nabla}\psi, \end{equation}

which reduces to $\boldsymbol {B} = G \boldsymbol {\nabla }\zeta _B$ for a vacuum magnetic field. Starting from a general set of coordinates $(s, \theta, \zeta )$, one may also write the vacuum field as

(C2)\begin{equation} \boldsymbol{B} = G\boldsymbol{\nabla}(\zeta + \omega), \end{equation}

with $\omega$ a single-valued function. The toroidal Boozer angle may thus be simply identified as

(C3)\begin{equation} \zeta_B = \zeta + \omega. \end{equation}

The poloidal Boozer angle can be obtained after solving the magnetic differential equation $\boldsymbol {B}\boldsymbol {\cdot }\boldsymbol {\nabla }\alpha = 0$ for the field-line label $\alpha = \theta - \iota \zeta + \lambda (\theta, \zeta )$, with the single-valued function $\lambda$. Then, employing the fact that Boozer coordinates are straight-field-line coordinates ($\lambda =0$), the Boozer poloidal angle may be evaluated as

(C4)\begin{equation} \theta_B = \alpha + \iota \zeta_B. \end{equation}

Footnotes

1 A notable exception is the recent study by Buller et al. (Reference Buller, Landreman, Kappel and Gaur2024) exploring the compatibility of QA with varying rotational transform values for $N_\mathrm {fp}=2$ (two field-period) stellarators.

References

Anderson, F.S.B., Almagri, A.F., Anderson, D.T., Matthews, P.G., Talmadge, J.N. & Shohet, J.L. 1995 The helically symmetric experiment, (HSX) goals, design and status. Fusion Technol. 27 (3 T), 273–277.CrossRefGoogle Scholar
Bader, A., Drevlak, M., Anderson, D.T., Faber, B.J., Hegna, C.C., Likin, K.M., Schmitt, J.C. & Talmadge, J.N. 2019 Stellarator equilibria with reactor relevant energetic particle losses. J. Plasma Phys. 85 (5), 905850508.CrossRefGoogle Scholar
Beidler, C.D., Smith, H.M., Alonso, A., Andreeva, T., Baldzuhn, J., Beurskens, M.N.A., Borchardt, M., Bozhenkov, S.A., Brunner, K.J., Damm, H., et al. 2021 Demonstration of reduced neoclassical energy transport in Wendelstein 7-X. Nature 596 (7871), 221–226.CrossRefGoogle ScholarPubMed
Boozer, A.H. 2019 Curl-free magnetic fields for stellarator optimization. Phys. Plasmas 26 (10), 102504.CrossRefGoogle Scholar
Buller, S., Landreman, M., Kappel, J. & Gaur, R. 2024 A family of quasi-axisymmetric stellarators with varied rotational transform, arXiv:2401.09021.Google Scholar
Conlin, R., Dudt, D.W., Panici, D. & Kolemen, E. 2023 The DESC stellarator code suite. Part 2. Perturbation and continuation methods. J. Plasma Phys. 89 (3), 955890305.CrossRefGoogle Scholar
Drevlak, M., Brochard, F., Helander, P., Kisslinger, J., Mikhailov, M., Nührenberg, C., Nührenberg, J. & Turkin, Y. 2013 ESTELL: a quasi-toroidally symmetric stellarator. Contrib. Plasma Phys. 53 (6), 459–468.CrossRefGoogle Scholar
Dudt, D.W., Conlin, R., Panici, D. & Kolemen, E. 2023 The DESC stellarator code suite. Part 3: quasi-symmetry optimization. J. Plasma Phys. 89 (2), 955890201.CrossRefGoogle Scholar
Dudt, D.W. & Kolemen, E. 2020 DESC: a stellarator equilibrium solver. Phys. Plasmas 27 (10), 102513.CrossRefGoogle Scholar
Fletcher, R. 1987 Practical Methods of Optimization, 2nd edn. Wiley.Google Scholar
Greene, J.M. 1979 A method for determining a stochastic transition. J. Math. Phys. 20 (6), 1183–1201.CrossRefGoogle Scholar
Hall, L.S. & McNamara, B. 1975 Three-dimensional equilibrium of the anisotropic, finite-pressure guiding-center plasma: theory of the magnetic plasma. Phys. Fluids 18 (5), 552.CrossRefGoogle Scholar
Henneberg, S.A., Drevlak, M. & Helander, P. 2020 Improving fast-particle confinement in quasi-axisymmetric stellarator optimization. Plasma Phys. Control. Fusion 62 (1), 014023.CrossRefGoogle Scholar
Henneberg, S.A., Drevlak, M., Nührenberg, C., Beidler, C.D., Turkin, Y., Loizu, J. & Helander, P. 2019 Properties of a new quasi-axisymmetric configuration. Nucl. Fusion 59 (2), 026014.CrossRefGoogle Scholar
Hirshman, S.P., van Rij, W.I. & Merkel, P. 1986 Three-dimensional free boundary calculations using a spectral Green's function method. Comput. Phys. Commun. 43 (1), 143–155.CrossRefGoogle Scholar
Hudson, S.R., Dewar, R.L., Dennis, G., Hole, M.J., McGann, M., von Nessi, G. & Lazerson, S. 2012 Computation of multi-region relaxed magnetohydrodynamic equilibria. Phys. Plasmas 19 (11), 112502.CrossRefGoogle Scholar
Jorge, R., Goodman, A., Landreman, M., Rodrigues, J. & Wechsung, F. 2023 Single-stage stellarator optimization: combining coils with fixed boundary equilibria. Plasma Phys. Control. Fusion 65 (7), 074003.CrossRefGoogle Scholar
Landreman, M. 2022 Mapping the space of quasisymmetric stellarators using optimized near-axis expansion. J. Plasma Phys. 88 (6), 905880616.CrossRefGoogle Scholar
Landreman, M., Buller, S. & Drevlak, M. 2022 Optimization of quasi-symmetric stellarators with self-consistent bootstrap current and energetic particle confinement. Phys. Plasmas 29 (8), 082501.CrossRefGoogle Scholar
Landreman, M., Medasani, B., Wechsung, F., Giuliani, A., Jorge, R. & Zhu, C. 2021 a SIMSOPT: a flexible framework for stellarator optimization. J. Open Source Softw. 6 (65), 3525.CrossRefGoogle Scholar
Landreman, M., Medasani, B. & Zhu, C. 2021 b Stellarator optimization for good magnetic surfaces at the same time as quasisymmetry. Phys. Plasmas 28 (9), 092505.CrossRefGoogle Scholar
Landreman, M. & Paul, E.J. 2022 Magnetic fields with precise quasisymmetry for plasma confinement. Phys. Rev. Lett. 128 (3), 035001.CrossRefGoogle ScholarPubMed
Landreman, M. & Sengupta, W. 2019 Constructing stellarators with quasisymmetry to high order. J. Plasma Phys. 85 (6), 815850601.CrossRefGoogle Scholar
Nies, R., Paul, E.J., Hudson, S.R. & Bhattacharjee, A. 2022 Adjoint methods for quasi-symmetry of vacuum fields on a surface. J. Plasma Phys. 88 (1), 905880106.CrossRefGoogle Scholar
Nies, R., Paul, E.J., Panici, D., Hudson, S.R. & Bhattacharjee, A. 2024 Data for “Exploration of the parameter space of quasisymmetric stellarator vacuum fields through adjoint optimisation”. Princeton Plasma Physics Laboratory, Princeton University.Google Scholar
Nührenberg, J. & Zille, R. 1988 Quasi-helically symmetric toroidal stellarators. Phys. Lett. A 129 (2), 113–117.CrossRefGoogle Scholar
Panici, D., Conlin, R., Dudt, D.W., Unalmis, K. & Kolemen, E. 2023 The DESC stellarator code suite. Part 1. Quick and accurate equilibria computations. J. Plasma Phys. 89 (3), 955890303.CrossRefGoogle Scholar
Paul, E.J., Bhattacharjee, A., Landreman, M., Alex, D., Velasco, J.L. & Nies, R. 2022 Energetic particle loss mechanisms in reactor-scale equilibria close to quasisymmetry. Nucl. Fusion 62 (12), 126054.CrossRefGoogle Scholar
Paul, E.J., Landreman, M. & Antonsen, T. 2021 Gradient-based optimization of 3D MHD equilibria. J. Plasma Phys. 87 (2).CrossRefGoogle Scholar
Plunk, G.G. & Helander, P. 2018 Quasi-axisymmetric magnetic fields: weakly non-axisymmetric case in a vacuum. J. Plasma Phys. 84 (2), 905840205.CrossRefGoogle Scholar
Qian, T.M., Chu, X., Pagano, C., Patch, D., Zarnstorff, M.C., Berlinger, B., Bishop, D., Chambliss, A., Haque, M., Seidita, D., et al. 2023 Design and construction of the MUSE permanent magnet stellarator. J. Plasma Phys. 89 (5), 955890502.CrossRefGoogle Scholar
Qu, Z.S., Pfefferlé, D., Hudson, S.R., Baillod, A., Kumar, A., Dewar, R.L. & Hole, M.J. 2020 Coordinate parameterisation and spectral method optimisation for Beltrami field solver in stellarator geometry. Plasma Phys. Control. Fusion 62 (12), 124004.CrossRefGoogle Scholar
Rodríguez, E., Paul, E.J. & Bhattacharjee, A. 2022 a Measures of quasisymmetry for stellarators. J. Plasma Phys. 88 (1), 905880109.CrossRefGoogle Scholar
Rodríguez, E., Sengupta, W. & Bhattacharjee, A. 2022 b Phases and phase-transitions in quasisymmetric configuration space. Plasma Phys. Control. Fusion 64 (10), 105006.CrossRefGoogle Scholar
Spitzer, L. 1958 The stellarator concept. Phys. Fluids 1 (4), 253.CrossRefGoogle Scholar
Virtanen, P., Gommers, R., Oliphant, T.E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., et al. 2020 SciPy 1.0: fundamental algorithms for scientific computing in Python. Nat. Meth. 17 (3), 261–272.CrossRefGoogle ScholarPubMed
Walker, S.W. 2015 The Shapes of Things: A Practical Guide to Differential Geometry and the Shape Derivative. Advances in Design and Control. Society for Industrial and Applied Mathematics.CrossRefGoogle Scholar
Zarnstorff, M.C., Berry, L.A., Brooks, A., Fredrickson, E., Fu, G.Y., Hirshman, S., Hudson, S., Ku, L.P., Lazarus, E., Mikkelsen, D., et al. 2001 Physics of the compact advanced stellarator NCSX. Plasma Phys. Control. Fusion 43 (12A), A237–A249.CrossRefGoogle Scholar
Figure 0

Figure 1. Contours of the magnetic field strength $B$ on the boundary of quasisymmetric stellarators at small aspect ratio $A$: top view (a–d), side view (e–h) and in Boozer coordinates (i–l).

Figure 1

Figure 2. Volume properties of stellarators optimised for QS on the boundary, with varying aspect ratio $A$. The QS error (a) and rotational transform (b) profiles were obtained using the DESC code. The QS error reaches record levels on the boundary ($s=1$), and generally remains relatively low in the core. The rotational transform profiles are flat for all configurations except for $N_\mathrm {fp}=3$ QA.

Figure 2

Table 1. Properties of stellarators optimised for QS at varying aspect ratio $A$: QS figure of merit $f_\mathrm {QS}^\star$, maximum symmetry-breaking mode on boundary from vacuum solution (SPEC) and MHS solution (DESC), volume-averaged QS error in MHS solution, difference between $\iota$ on axis and at the edge $\Delta \iota = \iota (s=1)-\iota (s=0)$ and minimum value of axis curvature $\kappa _\mathrm {min}$.

Figure 3

Figure 3. Contours of the magnetic field strength $B$ on the boundary of quasisymmetric stellarators at varying edge rotational transform $\iota _e$: top view (a–d), side view (e–h) and in Boozer coordinates (i–l).

Figure 4

Figure 4. Volume properties of stellarators optimised for QS on the boundary, with varying edge rotational transform $\iota _e$. The QS error (a) and rotational transform (b) profiles were obtained using the DESC code. The QS error reaches record levels on the boundary ($s=1$), and generally remains relatively low in the core. The rotational transform profiles are flat for all configurations except for $N_\mathrm {fp}=3$ QA and the high $\iota _e=1.97$ QH.

Figure 5

Table 2. Properties of stellarators optimised for QS at varying edge rotational transform $\iota _e$: QS figure of merit $f_\mathrm {QS}^\star$, maximum symmetry-breaking mode on boundary from vacuum solution (SPEC) and MHS solution (DESC), volume-averaged QS error in MHS solution, difference between $\iota$ on axis and at the edge $\Delta \iota = \iota (s=1)-\iota (s=0)$ and minimum value of axis curvature $\kappa _\mathrm {min}$.

Figure 6

Figure 5. Magnetic axes’ top view (a), curvature (b) and torsion (c), of the stellarators optimised for QS on the boundary with varying edge rotational transform $\iota _e$. The DESC solution is used to obtain the magnetic axis. Regions of reduced curvature develop at high $\iota$ for the QA, and at low $\iota$ for the QH configurations, as predicted by Rodríguez et al. (2022b). The legend for the various colours and shadings is given in figure 4.

Figure 7

Figure 6. Two field-period QA with large edge rotational transform $\iota _e = 0.65$ before (Boundary QA) and after (Global QA) auxiliary optimisation for QS in the volume. The QS error (a), evaluated using the VMEC code, shows an improvement in the global QS level, at the cost of higher QS error on the boundary. The optimisation did not require a large modification to the stellarator shape, as attested by the cross-sections of the boundary (b), shown for $\zeta /(2{\rm \pi} / 5N_\mathrm {fp}) \in \{0,1,2,3,4\}$.

Figure 8

Figure 7. Optimisation for integrability and volume QS of $N_\mathrm {fp}=3$ QA with high edge rotational transform $\iota _e=0.72$. The integrability optimisation leads to the disappearance of the magnetic island chains, such that the final configuration has optimised configuration. Furthermore, the volume QS error (evaluated using the VMEC code) is also reduced. The required change in the stellarator shape is minimal, as attested by the cross-sections of the boundary at $\zeta /(2{\rm \pi} /5N_\mathrm {fp}) \in \{0,1,2,3,4\}$. (a) Poincaré before integrability optimisation. (b) Poincaré after integrability optimisation. (c) Quasisymmetry error. (d) Boundary cross-sections at varying $\zeta$.