A single-field-period quasi-isodynamic stellarator

Abstract

A single-field-period quasi-isodynamic stellarator configuration is presented. This configuration, which resembles a twisted strip, is obtained by the method of direct construction, that is, it is found via an expansion in the distance from the magnetic axis. Its discovery, however, relied on an additional step involving numerical optimization, performed within the space of near-axis configurations defined by a set of adjustable magnetic field parameters. This optimization, completed in 30 s on a single CPU core using the SIMSOPT code, yields a solution with excellent confinement, as measured by the conventional figure of merit for neoclassical transport, effective ripple, at a modest aspect ratio of eight. The optimization parameters that led to this configuration are described, its confinement properties are assessed and a set of magnetic field coils is found. The resulting transport at low collisionality is much smaller than that of W7-X, and the device needs significantly fewer coils because of the reduced number of field periods.

1. Introduction

In 1951, Lyman Spitzer proposed the stellarator concept, where magnetic field lines are twisted by deforming a torus to break axisymmetry (Spitzer 1958). Spitzer's idea involved the shape of a figure eight, where he was able to calculate its corresponding rotational transform, and was realized experimentally with the Model A, B and C stellarators (Stix 1998). During the next several decades, the stellarator concept saw important developments in the calculation of magnetohydrodynamic (MHD) equilibria with good nested flux surfaces, high-$\beta$ stability properties, coil optimization, reduced neoclassical transport and other improvements (Gates et al. 2017). More recently, the figure-eight configuration proposed by Spitzer has continued to inspire advancements in the design of stellarators and linked mirrors (Feng et al. 2021).

There are currently three types of optimized stellarator configurations that appear to have the potential to become future fusion reactors: (1) quasi-isodynamic (QI) (Gori, Lotz & Nuhrenberg 1994), (2) quasi-axisymmetric (QA) and (3) quasi-helically symmetric (QH) ones (Boozer 2015). In these configurations, termed omnigenous, the second (or longitudinal) adiabatic invariant, $J$, is a function of the toroidal magnetic flux only, $J = J(\psi )$, which can be shown to lead to confined guiding-centre orbits (Helander 2014). The difference between the three cases lies mainly in the way that the contours of constant magnetic field strength $|\boldsymbol B|$ close when plotted on magnetic flux surfaces: in QI they close poloidally, in QA toroidally and in QH helically. Additionally, while QA and QH belong to the class of quasi-symmetric stellarators, where all contours of $|\boldsymbol B|$ are straight in Boozer and Hamada coordinates (D'haeseleer et al. 1991), in QI only the contours of the maxima of $|\boldsymbol B|$ are required to be straight. All three types of configurations can be ‘directly constructed’ using the near-axis expansion, based on an approach using Boozer coordinates (Garren & Boozer 1991; Landreman, Sengupta & Plunk 2019; Plunk, Landreman & Helander 2019), while only quasi-symmetric stellarators have been obtained using an approach based on Mercier coordinates (Jorge, Sengupta & Landreman 2020b, a). The largest stellarator in operation, Wendelstein 7-X (Grieger et al. 1992), is of the QI type, though only very approximately given the fact that the magnetic field strength at local maxima and minima varies substantially over the flux surface, there is substantial neoclassical transport and losses of fast ions and the bootstrap current are non-negligible. This type of configuration is relatively insensitive to the plasma pressure since the Shafranov shift is small (Gori et al. 1994; Nührenberg et al. 1995; Cary & Shasharina 1997; Nührenberg 2010) and the bootstrap current vanishes identically at low collisionality (Helander & Nuhrenberg 2009; Helander, Geiger & Maaßberg 2011; Landreman & Catto 2012). Unlike quasi-symmetric configurations, QI ones may thus be operated with essentially no net toroidal current experimentally, even with finite plasma pressure.

In this work, we construct a QI configuration using a near-axis expansion framework based on Boozer coordinates, which reduces the computational effort considerably compared with other approaches (the advantages and limitations of the near-axis expansion framework are detailed in the next section). Unlike nearly all previous stellarator designs, our configuration only has a single field period and a racetrack shape that resembles Spitzer's original idea. In contrast to the latter, however, it has carefully shaped flux surfaces in order to satisfy the stringent requirements of quasi-isodynamicity. As shown recently (Landreman & Sengupta 2019; Jorge et al. 2020a; Landreman 2021; Landreman & Jorge 2020), the near-axis expansion not only can be used as a tool to find configurations with enhanced particle confinement, but it also has the potential to find configurations that are Mercier stable, if the expansion is carried to second order. Our approach, based on a first-order near-axis expansion, relies on a careful choice of these parameters such that good flux surfaces are found even for low to medium aspect ratios, i.e. for $R/a$ between 6 and 10, which are within the realm of the near-axis expansion framework. In the past, QI stellarators have usually had aspect ratios of at least 10 (as for example the W7-X A configuration of Geiger et al. (2015) at $\beta =0$).

The QI configuration found here is obtained by varying the input parameters for the near-axis expansion, namely the magnetic axis and the zeroth- and first-order magnetic fields, and optimizing the resulting output such that a given objective function is minimized. This optimization is performed using the SIMSOPT code (Landreman et al. 2021a; Landreman, Medasani & Zhu 2021b), which is able to perform gradient-based optimizations with MPI parallelization of a finite-difference method and can easily be generalized to a Variational Moments Equilibrium Code (VMEC)-based optimization approach that uses near-axis configurations as initial conditions. After a suitable set of parameters is found, a finite-aspect-ratio configuration is constructed and its properties are assessed using several numerical tools commonly used in stellarator optimization studies: the VMEC (Hirshman, van RIJ & Merkel 1986), the BOOZ_XFORM code (Sanchez et al. 2000), the NEO code (Nemov et al. 1999), the Stepped Pressure Equilibrium Code (SPEC) (Qu et al. 2020; Hudson et al. 2012), the coil optimization suite ONSET (Drevlak 1998) and the SIMPLE code (Albert, Kasilov & Kernbichler 2020). Their main functions and results are described in the following sections. The method used to construct near-axis QI fields is that of Plunk et al. (2019), although we choose different forms of some of the free functions that parametrize the solution space. The most notable example is a function that controls the deviation from omnigenity near maxima of the magnetic field strength (see (2.8)); this change, which improves the smoothness of solutions, will be explained in more detail in an upcoming publication (Camacho Mata, Plunk & Jorge 2022).

This paper is organized as follows. In § 2 the near-axis expansion formalism is outlined, in particular its application to QI fields and the corresponding numerical implementation. In § 3 the physics-based figures of merit are described and their analytical expressions in the near-axis expansion formalism are shown. The resulting optimized configuration is shown in § 4, and our conclusions follow in the final section.

2. The near-axis expansion

The near-axis expansion solves the equilibrium MHD equations by performing an expansion in the inverse aspect ratio $\epsilon$:

(2.1)\begin{equation} \epsilon=\frac{a}{R}\ll 1, \end{equation}

where $a$ and $R$ are measures of the minor and major radius of the device, respectively. We note that, although the construction is based on an expansion on $\epsilon$, it is able to describe the core region of any configuration, including those with low aspect ratio. We employ the near-axis expansion using Boozer coordinates $(\psi,\theta,\varphi )$ (Boozer 1981), with $\psi$ the toroidal magnetic flux divided by $2{\rm \pi}$, $\theta$ the poloidal angle and $\varphi$ the toroidal angle, writing the magnetic field vector as

(2.2)\begin{align} \boldsymbol B & = \boldsymbol{\nabla} \psi \times \boldsymbol{\nabla} \theta + \iota(\psi) \boldsymbol{\nabla} \varphi \times \boldsymbol{\nabla} \psi \end{align}

(2.3)\begin{align} & = \beta(\psi,\theta,\varphi) \boldsymbol{\nabla}\psi + I(\psi) \boldsymbol{\nabla} \theta + G(\psi) \boldsymbol{\nabla} \varphi, \end{align}

where $\iota =\iota (\psi )$ is the rotational transform, $I(\psi )$ is $\mu _0/(2{\rm \pi} )$ times the toroidal current enclosed by the flux surface, $G(\psi )$ is $\mu _0/(2{\rm \pi} )$ times the poloidal current outside the flux surface and $\beta$ is related to the Pfirsch–Schlüter current.

In the near-axis expansion method for QI configurations, the position vector $\boldsymbol {r}$ is written as

(2.4)\begin{equation} \boldsymbol{r}(r,\theta,\varphi) = \boldsymbol{r}_0(\varphi) +X(r,\theta,\varphi) \boldsymbol{n}^s(\varphi) +Y(r,\theta,\varphi) \boldsymbol{b}^s(\varphi) +Z(r,\theta,\varphi) \boldsymbol{t}^s(\varphi), \end{equation}

where $r=\sqrt {2 \psi /\bar {B}}$ is a radial-like variable, $\bar {B}$ is a constant reference field strength and $\boldsymbol {r}_0(\varphi )$ is the magnetic axis curve parametrized using $\varphi$. We employ a modified Frenet–Serret frame (Carroll, Kose & Sterling 2013) where the curvature is replaced by the signed curvature $\kappa ^s=s \kappa$, where $s(\varphi )$ takes values of $\pm 1$ and switches at locations of zero axis curvature (a characteristic of QI configurations). The normal and binormal vectors are also multiplied by $s$. Note that the Frenet–Serret formulas are invariant under such substitution. We therefore write the signed Frenet–Serret frame as $(\boldsymbol {t},\boldsymbol {n}^s=s \boldsymbol n,\boldsymbol {b}^s=s \boldsymbol b)$, where $(\boldsymbol {t},\boldsymbol n, \boldsymbol b)$ are the tangent, normal and binormal unit vectors of the Frenet–Serret frame of the magnetic axis. In the following, the arc length along the axis is denoted by $\ell$ with $0\le \ell < L$, the axis curvature by $\kappa (\varphi )$ and the axis torsion by $\tau (\varphi )$ with the sign convention of Landreman (2019).

We expand the magnetic field and the position vector only up to first order in $r$ as (Garren & Boozer 1991)

(2.5)\begin{equation} B=B_0[1+r d \cos(\theta-\alpha)], \end{equation}

where $B_0=B_0(\varphi )$ is the magnetic field on axis, $d=d(\varphi )$ is a free function describing the first-order magnetic field strength, while $\alpha =\alpha (\varphi )$ is an angle-like variable also describing the first-order magnetic field strength. This yields (Landreman & Sengupta 2019)

(2.6)\begin{equation} X = \frac{r d}{\kappa^s} \cos(\theta-\alpha),\ Y = \frac{r \kappa^s }{d}\frac{\bar{B}}{B_0}\left[\sin(\theta-\alpha)+\sigma \cos(\theta-\alpha)\right],\ Z = 0, \end{equation}

where $G= G_0 = L/\int _0^{2{\rm \pi} }(d\varphi /B_0)$ and $\sigma$ is a solution of

(2.7)\begin{equation} \frac{{\rm d} \sigma}{{\rm d}\varphi}+\gamma\left(1+\sigma^2+\frac{B_0^2}{\bar{B}^2}\frac{{\rm d}^4}{(\kappa^s)^4}\right)-2\left(\frac{I_2}{\bar{B}}-\tau\right)\frac{G_0 d^2}{\bar{B} (\kappa^s)^2}=0, \end{equation}

with $I=r^2 I_2$ and $\gamma =\iota -\alpha '(\varphi )$. As the plasma pressure only appears at second order in the expansion, the configurations considered here are effectively force-free configurations. Incidentally, the fact that the axis curvature should vanish at points where $B_0'(\varphi )=0$ stems from the fact that QI fields need to have $d=0$ at all local extrema combined with $Y$ being proportional to $r \kappa /d$ (Plunk et al. 2019).

As shown by Cary & Shasharina (1997), QI fields are necessarily non-analytic. Furthermore, for the near-axis expansion case, the omnigenity condition leads to the relation $\alpha (2{\rm \pi} )-\alpha (0)=2 {\rm \pi}\iota$, while periodicity of the magnetic field in (2.5) requires $\alpha (2{\rm \pi} )-\alpha (0)=2 {\rm \pi}N$ (Plunk et al. 2019). To alleviate this conflict between omnigenity and periodicity, we choose the function $\alpha$ such that omnigenity is violated in a controlled way, by writing it as

(2.8)\begin{equation} \alpha = \iota\left(\varphi-\frac{\rm \pi}{N_{fp}}\right)\left[1+\frac{\rm \pi}{\iota}\left(\frac{N-\iota/N_{fp}}{({\rm \pi}/N_{fp})^{2k+1}}\right)\left(\varphi-\frac{\rm \pi}{N_{fp}}\right)^{2k}\right]+{\rm \pi}\left(2N+\frac{1}{2}\right), \end{equation}

where $N_{fp}$ is the number of field periods of the device. The integer $k$ in (2.8) effectively controls the spatial distribution along $\varphi$ of the deviation from omnigenity. For a more detailed discussion and a comparison with other forms of $\alpha$ including the exact omnigenous form, see Camacho Mata et al. (2022).

We consider magnetic fields with a single minimum along the magnetic axis of the form

(2.9)\begin{equation} B_0 = B_{00}+B_{01} \cos(N_{fp}\varphi), \end{equation}

and parametrize the axis curve $\boldsymbol r_0$ as

(2.10)\begin{equation} \boldsymbol r_0(\phi) = \sum_l R_l \cos(N_{fp} l \phi) \boldsymbol e_R + \sum_l Z_l \sin(N_{fp} l \phi) \boldsymbol e_Z, \end{equation}

where $\phi$ is the standard cylindrical toroidal angle, $\boldsymbol e_R$ and $\boldsymbol e_Z$ are the standard cylindrical coordinate unit vectors and stellarator symmetry is assumed. The form for $d$ chosen here is similar to the approach taken by Plunk et al. (2019), with an added term proportional to the curvature of the magnetic axis:

(2.11)\begin{align} d = d_\kappa \kappa + \sum_l \bar{d}_l \sin(N_{fp} l \varphi), \end{align}

where $d_\kappa$ and $\bar {d}_i$ are constants.

The equation for $\sigma$, (2.7), is solved using Newton iteration with a pseudo-spectral collocation discretization together with the constraint $\sigma (\phi =0)=\sigma _0$. In the following, the parameter $\sigma _0$ is set to $0$ in order to enforce stellarator symmetry. If $\phi =0$ is not one of the grid points, this condition is imposed by interpolating $\sigma$ using pseudo-spectral interpolation. A uniform grid of $N_\phi$ points is used, with $\phi _j=(j+j_{{\rm shift}}-1)2{\rm \pi} /(N_\phi N_{fp})$ for $j=1 \ldots N_\phi$ and $j_{{\rm shift}}=0$ and $1/3$ in the quasi-symmetric and QI cases, respectively. The value of $j_{{\rm shift}}=1/3$ is used to avoid points along the axis where the curvature reaches zero but different values of $j_{{\rm shift}}$ could be used to achieve the same effect. The discrete unknowns include the values of $\sigma$ on the $\phi$ grid and $\iota _0$. As a boundary condition, we impose periodicity in $\sigma$, with $\sigma (0)=\sigma (2{\rm \pi} )=\sigma _0$, yielding a single value of $\iota _0$ and the function $\sigma (\phi )$ as solution. Finally, a conversion to cylindrical coordinates is performed in order to create VMEC and SPEC input files and for visualization purposes using the nonlinear method described in Landreman et al. (2019) and by choosing a particular value for the radius $r$.

3. Optimization method

The optimization is performed with the SIMSOPT code, using the trust region reflective algorithm for nonlinear least squares problems from the scipy package using the Python programming language. The input parameters for the optimization are the axis shape coefficients $(R_i,Z_i)$ except $R_0=1$, which is fixed, the magnetic field on axis, where we fix $B_{00}=1$ and vary $B_{01}$, and the scalars $d_\kappa$ and $\bar {d}_i$ present in the first-order magnetic field function $d$ in (2.11). The additional parameters $N_{fp}$, the number of field periods, and the exponent $k$ in (2.8) were varied manually. As it was seen that values of $k=3$ and $N_{fp}=1$ yielded configurations with consistently lower elongation and neoclassical transport (as measured by $\epsilon _{\text {eff}}$; see below), these parameters were then held fixed during the construction of the configuration presented here.

The optimization is performed in a series of steps. As an initial condition, we employ $R_0=1$, $R_2=-0.2$, $Z_2=0.35$, $(R_1,Z_1)=0.0$, $B_{01}=0.16$ (a value used later as reference), $d_\kappa =0.5$, $\bar {d}_0=0$ and $\bar {d}_1=0.01$. The values of $R_0, R_2$ and $Z_2$ are those of the configuration in Plunk et al. (2019). Then, SIMSOPT is called inside a loop over the number of axis coefficients, starting at only two coefficients up to 12 sequentially, that is, there are a total of 12 free axis coefficients, $R_2, R_4, R_6, R_8, R_{10}, R_{12}, Z_2, Z_4, Z_6, Z_8, Z_{10}, Z_{12}$, and therefore there are 6 steps, where first $R_2, Z_2$ is allowed to vary, then $R_4, Z_4$, then $R_6, Z_6$ and so forth. We note that, as more axis coefficients are introduced, the previous axis coefficients are still varied within the optimization. By choosing even mode numbers only, we obtain a two-field-period axis shape with two points of zero curvature at the points of the extrema of $|\boldsymbol B|$ ($\varphi =0$ and $\varphi ={\rm \pi}$), as required by the omnigenity condition (Plunk et al. 2019). The grid resolution is $N_\phi =131$ and is increased by 20 each time the number of axis coefficients increases. Each iteration inside the loop is run until the change of the cost function is smaller than $10^{-4}$, which usually takes between 10 and 40 steps to be achieved. The optimization process takes less than 30 s using a single CPU core.

The objective function has the following form:

(3.1)\begin{align} f_{QI}& =\frac{w_{L_{\boldsymbol{\nabla} B}}}{L_{\boldsymbol{\nabla} \boldsymbol B}^2}+w_{\mathrm{E}} |\mathrm{E}|^2+w_{R_0}[\text{min}(R_{{\rm axis}})-R_{\text{min}}]^2+w_{Z_0}[\text{max}(Z_{{\rm axis}})-Z_{\text{max}}]^2\nonumber\\ & \quad +w_d |d|^2+w_{\bar{d}_s}\sum_i |\bar{d}_i|^2+w_{\bar{B}_{01}}|\bar{B}_{01}-\bar{B}^{\text{ref}}_{01}|^2+w_{d'(0)}|d'(0)|^2+w_{\alpha}|\alpha-\alpha_0|^2. \end{align}

In (3.1), $L_{\boldsymbol {\nabla } \boldsymbol B}$ is the scale length associated with the Frobenius norm $|\boldsymbol {\nabla } \boldsymbol B|$ of the $\boldsymbol {\nabla } \boldsymbol B$ tensor (equation (3.11) in Landreman (2021)), given by $L_{\boldsymbol {\nabla } \boldsymbol B} = B_0 \sqrt {2 / |\boldsymbol {\nabla } \boldsymbol B|}$, $\min (R_{{\rm axis}})$ and $\max (Z_{{\rm axis}})$ are the value of the minimum cylindrical radial and maximum vertical coordinates, respectively, of the axis curve $\boldsymbol r_0$, $|d|$ and $|\mathrm {E}|$ are the $L2$ norm of the discretized functions $d$ and the elongation $\mathrm {E}=a/b$ associated with the first-order magnetic field in (2.5) where $a$ and $b$ are the semi-major and semi-minor axes of the elliptical flux surface cross-section, respectively, and $\alpha _0=\iota (\varphi -{\rm \pi} )$ is an exactly omnigenous version of the function $\alpha$ in (2.8). The terms in (3.1) have three main goals: (1) select configurations with small deviations from QI (min of $\alpha -\alpha _0$, $\bar {d}$, $d$ and $d'(0)$); (2) select axis shapes with small elongation and aspect ratio (reduction of $1/L_{\boldsymbol {\nabla } \boldsymbol B}$, $\min (R_{{\rm axis}})$, $\max (Z_{{\rm axis}})$ and $E$); and (3) penalize high mirror ratios (reduction of $\bar {B}_{01}-\bar {B}_{01}^{{\rm ref}}$). For a more in-depth assessment of the relation between axis shapes and the resulting elongation and the difficulty of obtaining solutions with low elongation in the near-axis expansion framework, we refer the reader to Camacho Mata et al. (2022). As parameters in (3.1) we choose $R_{\text {min}}=Z_{\text {max}}=0.4$ and $\bar {B}_{01}^{\text {ref}} = 0.16$. The weights used in the optimization process are the following: $w_{L_{\boldsymbol {\nabla } B}}=0.03,w_e=0.4/N_{\phi },w_{R_0}=w_{Z_0}=30,w_d=20/N_{\phi },w_{\bar {d}_s}=100,w_{\bar {B}_{01}}=200,w_{d'(0)}=2,w_{\alpha }=60$. Such values were found by first scaling the weights such that every term in the objective function has an order of magnitude of unity when a reasonable equilibrium is found and then they are fine tuned to ensure the three main goals described before.

4. Results

The optimization procedure resulted in the following parameters:

(4.1)\begin{align} R_0 (\phi)\ [\text{m}] & = 1 - 0.40566229 \cos(2\phi)+0.07747378 \cos (4 \phi)-0.00780386 \cos(6 \phi),\end{align}

(4.2)\begin{align} Z_0 (\phi)\ [\text{m}] & ={-}0.24769666 \sin(2\phi)+0.06767352 \sin(4 \phi)-0.00698062\sin(6\phi)\nonumber\\ & \quad -6.8162\times 10^{{-}4}\sin(8\phi)-1.451\times 10^{{-}5}\sin(10\phi)-2.83\times 10^{{-}6}\sin(12\phi), \end{align}

(4.3)\begin{align} B_0(\phi)\ [\text{T}] & = 1+0.16915531 \cos(\varphi),\end{align}

(4.4)\begin{align} d (\phi) & = 0.00356311 \sin(2\varphi)+0.00020159\sin(4\varphi)-0.00121786\sin(6\varphi)\nonumber\\ & \quad -0.00011629\sin(8\varphi)-8.3\times 10^{{-}7}\sin(10\varphi)+3.20\times 10^{{-}6}\sin(12\varphi)\nonumber\\ & \quad +0.51837838\kappa. \end{align}

This configuration has a rotational transform on axis of $\iota _0 = 0.680$, a maximum elongation of $\mathrm {E}=5.100$, a derivative of $d$ on axis of $d'(0)=0.094$ and a total value of the objective function of $f_{QI}=32.618$. The coefficients $R_8$, $R_{10}$ and $R_{12}$ are smaller than $10^{-8}$ and were therefore set to zero as they have a negligible contribution to the properties of the equilibrium.

Data for the magnetic configurations are available at Jorge (2022).

To create a plasma boundary, we chose the radial near-axis coordinate as $r=1/9$, which leads to a calculated minor and major radius of Aminor_p = 0.125 (corresponding VMEC output parameter) and Rmajor_p = 0.993 (corresponding VMEC output parameter), respectively, leading to an aspect ratio of $\epsilon$ = 7.944. The resulting boundary can be seen in figure 1, while the magnetic field and the magnetic axis can be seen in figure 2. While the axis shape is similar that found in Plunk et al. (2019), namely a racetrack oval with the points of vanishing curvature at the middle of each straight section and the surfaces resembling a twisted strip, there is no sudden twist near the region of maximum field strength. This not only results in better convergence of the VMEC at lower aspect ratios, but has a practical consequence of simplifying the coil shapes needed to produce such a configuration. As is evident from figure 2, not all contours of $|\boldsymbol B|$ close poloidally, but as we shall see below, the neoclassical transport is nevertheless very small.

Figure 1. Shape of the idealized constructed configuration from the near-axis expansion. (a) Bird's eye and side views of the boundary shape in three dimensions with an aspect ratio of 8. (b) Cross-sections of the configuration at eight values of toroidal angle $\phi$.

Figure 2. Elements of the near-axis construction. (a) Magnetic field on the boundary using (2.5). (b) Shape of the magnetic axis with white spheres in the locations of zero curvature. In red we depict the signed normal vector, in blue the signed binormal vector and in green the tangent vector.

We then use this boundary shape as input to VMEC, a magnetic field equilibrium code. The VMEC uses an inverse moment method where the cylindrical coordinates $R$ (radial) and $Z$ (vertical) are expanded in a double Fourier series involving a poloidal angle and the cylindrical toroidal angle. The resulting rotational transform profile is shown in figure 3(a). The rotational transform predicted by the near-axis expansion is $\iota =0.680$ while VMEC finds a relatively linear rotational transform with an on-axis value of $\iota =0.671$ and a value of $\iota =0.685$ at the edge. We note that, while the rotational transform from the near-axis expansion is fixed to the on-axis value in figure 3, the inclusion of magnetic shear is possible by performing an expansion to third order as shown in Rodríguez, Sengupta & Bhattacharjee (2022).

Figure 3. (a) Profile of the rotational transform $\iota$ from VMEC (blue) and $\iota$ on axis from the near-axis expansion (red). Contours of constant magnetic field strength in Boozer coordinates at $s=0.17$ (b), $s=0.5367$ (c) and $s=0.9033$ (d).

Based on the VMEC result, we use the BOOZ_XFORM code to find the magnetic field strength $|B|$ as a function of Boozer coordinates. This allows us to draw contours of constant $|B|$ and compare with the predicted contours in figure 2. The resulting properties of $|B|$ computed with BOOZ_XFORM are shown in figure 3. In figure 3(b–d), contours of constant $|B|$ are shown at $s=0.17$, $0.5367$ and $0.9033$, where $s=\psi /\psi _b$ with $\psi _b$ the toroidal magnetic flux at the boundary.

Next, we calculate the effective helical ripple $\epsilon _{\text {eff}}$ (Beidler et al. 1990; Nemov et al. 1999), which quantifies the direct effect of the radial magnetic drift of trapped-particle orbits on neoclassical transport in the so-called 1/$\nu$ transport regime. This parameter vanishes for perfectly omnigeneous configurations and is a function of the radial position $r^2$. In figure 4, the profile of $\epsilon _{\text {eff}}$ is shown for the present configuration, calculated using the NEO code. The $\epsilon _{\text {eff}}$ found here is substantially lower than that for W7-X, even for the lower aspect ratio employed here. For such levels of $\epsilon _{\text {eff}}$, the collisional transport is expected to almost certainly be weaker than the turbulent transport for this configuration.

Figure 4. The magnitude $\epsilon _{\text {eff}}$ of the $1/\nu$ transport for the constructed configuration (labelled as QI ${\rm NFP} = 1$) and for the standard configuration of W7-X in fixed boundary mode.

While VMEC can be used to analyse the approximate magnetic field inside a given boundary, it cannot be used to directly examine magnetic islands. As discussed by Reiman et al. (2007), low-aspect-ratio configurations are particularly prone to have magnetic islands, which are caused by undesirable radial field components at surfaces with rational rotational transform. As the presence of island chains and chaotic field regions in the core significantly degrades confinement, practically useful MHD equilibria should avoid or limit this phenomenon (Yamazaki et al. 1993; Neilson et al. 2010). In contrast to VMEC, the SPEC is able to compute magnetic islands (Hudson et al. 2012). It does so using the formalism of multi-region relaxed MHD to divide the computational domain into a number of nested annular regions. While the magnetic field is required to be tangential to the boundary of each region, there is no requirement that magnetic surfaces exist within them. The SPEC can thus be used to assess the presence of magnetic islands in the interior of such regions. For the optimized configurations found here, this procedure is particularly simple since the pressure vanishes and only a single region needs to be considered, which corresponds to the vacuum problem exactly (in the numerical sense). The Poincaré plot of such a calculation using SPEC is shown in figure 5, where it is seen that no major resonances are encountered and that flux surfaces are smoothly nested from the magnetic axis to the computational boundary.

Figure 5. Poincaré plot computed from the SPEC solution at $\phi =0$.

Next, we assess the confinement of 3.5 MeV alpha particles, which would be generated in a fusion reactor. The fraction of lost particles is evaluated using the drift-orbit code SIMPLE. For this study, we scale our magnetic configuration to a reactor size with minor radius of 1.7 m and an average on-axis magnetic field strength of $B_{0,0}=5.7$ T. Figure 6 shows the loss fraction of fast particles following collisionless guiding-centre drift orbits. A total of 5000 test particles, equally distributed on each flux surface, were launched with uniformly distributed pitch angles and traced for 0.2 s, typical for the collisional slowing-down time in a fusion reactor, or until they cross the $s = 1$ boundary surface and are considered lost. As shown in figure 6, the loss fraction is approximately 3.8 % for the particles starting at a flux surface with $s = 0.06$ and for the ones starting at a flux surface with $s = 0.25$ (approximately half the radius) it is approximately $7.2\,\%$. For comparison, figure 6 also shows the loss fraction in a W7-X configuration scaled to the same minor radius and magnetic field. For this study, the vacuum W7-X standard configuration is used, which corresponds to the A configuration of Geiger et al. (2015) at $\beta =0$. Its loss fraction is higher than the single-field-period QI for all particles at $t=0.2$ s except those that are started at $s=0.9$. We note, however, that at such long time scales, collisional effects might start to play a role, making collisionless simulations less reliable. Indeed, for times $t$ smaller than 0.01 s, the single-field-period configuration provides a higher loss fraction mainly due to prompt losses. The source mechanism for such prompt losses will be the subject of future studies.

Figure 6. The fast particle loss fraction at flux surfaces between $s = 0.06$ and $s = 0.9$ for a scaled configuration to a minor radius of 1.7 m and a magnetic field $B_{0,0}=5.7\,{\rm T}$. (a) Present one-field-period QI configuration. (b) W7-X.

Finally, a set of 30 coils that approximately reproduces the toroidal magnetic surface in figure 1 was obtained. The goal here is to show that simple coils can be found for the proposed one-field-period configuration. This was done by first scanning the space of current potentials that approximate the target magnetic field. Their distance to the plasma boundary as well as the number of Fourier modes describing the analytic current can be varied. Then, the contours of the current potential, found using the NESCOIL code (Merkel 1987), with four toroidal and three poloidal modes that lie on a current-carrying surface 30 cm away from the plasma were transformed into 15 modular coils per half-period and optimized using the ONSET code (Drevlak 1998) into a constructable shape shown in figure 7. The coils are parametrized with three-dimensional splines independent of any constraining surface. The coil optimization technique used is similar to that of Lobsien, Drevlak & Sunn (2018) with the difference that a starting point was chosen that originates from a solution with higher Fourier modes. This makes the approximation of the target magnetic field better but the complexity of the starting coil configuration is prohibitive. Therefore, a first design step of the nonlinear coil optimization focused on reducing the coil complexity as well as lowering the field error was taken. This procedure is also described in Lobsien et al. (2020), which uses the same notation for the penalty values as described below. In the optimization, the average and maximal curvatures were reduced and the distance between adjacent coils was increased as well as the distance between coils and the plasma boundary. The minimum distance between coil centrelines is 0.338 m. The final optimization step focused on properties of the vacuum magnetic field, ensuring that no low-order rational values of $\iota$ are present inside the plasma boundary and that $\iota$ increases with minor radius. The shapes of the magnetic surfaces defined by following field lines are highly sensitive to the iota profile due to the low shear nature of this configuration. We show in figure 8 a Poincaré plot resulting from the coil configuration in figure 7, which produces an island chain surrounding the last closed flux surface and could potentially be used for an island divertor. We note that the appearance of an island chain outside the plasma boundary was not targeted directly. As seen in figure 9, the approximation of the target magnetic field defined by the VMEC solution converges to a normalized maximum field error of 2 % and a normalized average field error below 1 % after multiple optimization steps. The normalized field error is defined locally as

(4.5)\begin{equation} q_{le}(x) = \frac{|\boldsymbol{B}\boldsymbol{\cdot} \boldsymbol{n}|}{|\boldsymbol{B}|}, \end{equation}

where $\boldsymbol n$ is the normal unit vector perpendicular to the plasma boundary and the maximum local field error is $\max q_{le}$. The average field error is defined globally as

(4.6)\begin{equation} q_{ae}(x) = \frac{\int_A q_{le}\,{\rm d}A}{A}, \end{equation}

where $A$ is the area of the plasma boundary.

Figure 7. Coil shapes for the magnetic configuration obtained here.

Figure 8. Poincaré plot obtained with the coils shown in figure 7 at the location $\phi =0$.

Figure 9. Relative field error $q_{ae}$ between the magnetic field produced by the coil shapes and the VMEC boundary shape (vertical axis) as a function of the optimization step (horizontal axis).

5. Conclusion

In summary, a new QI configuration that exhibits a number of favourable properties has been found using a first-order near-axis expansion approach. It has relatively weak shaping (implying relatively simple geometry) other than strong elongation, one field period, has low neoclassical transport and can be realized with relatively simple coils. No attempt was made to ensure favourable MHD properties, such as stability and small Shafranov shift at finite plasma pressure. In contrast to standard optimization procedures where a plasma boundary is varied, this design was found by varying the degrees of freedom of the near-axis expansion, namely the magnetic axis and the lowest-order magnetic field strength. This procedure builds on the work of Plunk et al. (2019) and Camacho Mata et al. (2022), using a newly developed approach to perform a controlled approximation to omnigenity.

In future work, we intend to reduce the neoclassical transport further and improve its fast-particle confinement properties by extending the optimization procedure to second order in the inverse aspect ratio $\epsilon$. Carrying the expansion to this order will allow us to study configurations with finite plasma pressure and to optimize for other relevant quantities such as magnetic well and the maximum-$J$ property (Helander, Proll & Plunk 2013). Indeed, this configuration is characterized by having a magnetic hill instead of a magnetic well as such measure is not available as an explicit target for the optimization at first order in the expansion. The particular optimization procedure outlined here can also be applied to other types of omnigenous configurations such as QA and QH ones, and will be the subject of future work.