## 1. Introduction

In 1951, Lyman Spitzer proposed the stellarator concept, where magnetic field lines are twisted by deforming a torus to break axisymmetry (Spitzer Reference Spitzer1958). 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 Reference Stix1998). 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.* Reference Gates, Boozer, Brown, Breslau, Curreli, Landreman, Lazerson, Lore, Mynick and Neilson2017). 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.* Reference Feng, Yu, Jiang and Fu2021).

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 Reference Gori, Lotz and Nuhrenberg1994), (2) quasi-axisymmetric (QA) and (3) quasi-helically symmetric (QH) ones (Boozer Reference Boozer2015). 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 Reference Helander2014). 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.* Reference D'haeseleer, Hitchon, Callen and Shohet1991), 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 Reference Garren and Boozer1991; Landreman, Sengupta & Plunk Reference Landreman, Sengupta and Plunk2019; Plunk, Landreman & Helander Reference Plunk, Landreman and Helander2019), while only quasi-symmetric stellarators have been obtained using an approach based on Mercier coordinates (Jorge, Sengupta & Landreman Reference Jorge, Sengupta and Landreman2020*b*, Reference Jorge, Sengupta and Landreman*a*). The largest stellarator in operation, Wendelstein 7-X (Grieger *et al.* Reference Grieger, Beidler, Harmeyer, Lotz, KißLinger, Merkel, Nührenberg, Rau, Strumberger and Wobig1992), 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.* Reference Gori, Lotz and Nuhrenberg1994; Nührenberg *et al.* Reference Nührenberg, Lotz, Merkel, Nührenberg, Schwenn, Strumberger and Hayashi1995; Cary & Shasharina Reference Cary and Shasharina1997; Nührenberg Reference Nührenberg2010) and the bootstrap current vanishes identically at low collisionality (Helander & Nuhrenberg Reference Helander and Nuhrenberg2009; Helander, Geiger & Maaßberg Reference Helander, Geiger and Maaßberg2011; Landreman & Catto Reference Landreman and Catto2012). 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 Reference Landreman and Sengupta2019; Jorge *et al.* Reference Jorge, Sengupta and Landreman2020*a*; Landreman Reference Landreman2021; Landreman & Jorge Reference Landreman and Jorge2020), 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.* (Reference Geiger, Beidler, Feng, Maassberg, Marushchenko and Turkin2015) 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.* Reference Landreman, Medasani, Wechsung, Giuliani, Jorge and Zhu2021*a*; Landreman, Medasani & Zhu Reference Landreman, Medasani and Zhu2021*b*), 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 Reference Hirshman, van RIJ and Merkel1986), the BOOZ_XFORM code (Sanchez *et al.* Reference Sanchez, Hirshman, Ware, Berry and Spong2000), the NEO code (Nemov *et al.* Reference Nemov, Kasilov, Kernbichler and Heyn1999), the Stepped Pressure Equilibrium Code (SPEC) (Qu *et al.* Reference Qu, Pfefferlé, Hudson, Baillod, Kumar, Dewar and Hole2020; Hudson *et al.* Reference Hudson, Dewar, Dennis, Hole, McGann, von Nessi and Lazerson2012), the coil optimization suite ONSET (Drevlak Reference Drevlak1998) and the SIMPLE code (Albert, Kasilov & Kernbichler Reference Albert, Kasilov and Kernbichler2020). 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.* (Reference Plunk, Landreman and Helander2019), 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 Reference Camacho Mata, Plunk and Jorge2022).

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$:

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 Reference Boozer1981), 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

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

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 Reference Carroll, Kose and Sterling2013) 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 (Reference Landreman2019).

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

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 Reference Landreman and Sengupta2019)

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

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.* Reference Plunk, Landreman and Helander2019).

As shown by Cary & Shasharina (Reference Cary and Shasharina1997), 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.* Reference Plunk, Landreman and Helander2019). 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

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.* (Reference Camacho Mata, Plunk and Jorge2022).

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

and parametrize the axis curve $\boldsymbol r_0$ as

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.* (Reference Plunk, Landreman and Helander2019), with an added term proportional to the curvature of the magnetic axis:

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.* (Reference Landreman, Sengupta and Plunk2019) 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.* (Reference Plunk, Landreman and Helander2019). 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.* Reference Plunk, Landreman and Helander2019). 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:

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 (Reference Landreman2021)), 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.* (Reference Camacho Mata, Plunk and Jorge2022). 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:

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 (Reference Jorge2022).

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.* (Reference Plunk, Landreman and Helander2019), 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.

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 (Reference Rodríguez, Sengupta and Bhattacharjee2022).

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.* Reference Beidler, Grieger, Herrnegger, Harmeyer, Kisslinger, Lotz, Maassberg, Merkel, Nuehrenberg and Rau1990; Nemov *et al.* Reference Nemov, Kasilov, Kernbichler and Heyn1999), 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.

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.* (Reference Reiman, Hirshman, Hudson, Monticello, Rutherford, Boozer, Brooks, Hatcher, Ku and Lazarus2007), 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.* Reference Yamazaki, Yanagi, Ji, Kaneko, Ohyabu, Satow, Morimoto, Yamamoto and Motojima1993; Neilson *et al.* Reference Neilson, Gruber, Harris, Rej, Simmons and Strykowsky2010). In contrast to VMEC, the SPEC is able to compute magnetic islands (Hudson *et al.* Reference Hudson, Dewar, Dennis, Hole, McGann, von Nessi and Lazerson2012). 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.

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.* (Reference Geiger, Beidler, Feng, Maassberg, Marushchenko and Turkin2015) 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.

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 Reference Merkel1987), 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 Reference Drevlak1998) 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 (Reference Lobsien, Drevlak and Sunn2018) 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.* (Reference Lobsien, Drevlak, Kruger, Lazerson, Zhu and Pedersen2020), 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

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

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

## 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.* (Reference Plunk, Landreman and Helander2019) and Camacho Mata *et al.* (Reference Camacho Mata, Plunk and Jorge2022), 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 Reference Helander, Proll and Plunk2013). 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.

## Acknowledgements

*Editor Peter Catto thanks the referees for their advice in evaluating this article.*

## Funding

R.J. was supported by a grant by Alexander-von-Humboldt-Stiftung, Bonn, Germany, through a research fellowship. K.C.M. and J.-F.L. were supported by a grant from the Simons Foundation (560651).

## Declaration of interests

The authors report no conflict of interest.

## Data availability statement

The data that support the findings of this study are openly available in Zenodo at https://doi.org/10.5281/zenodo.7040688, reference number 7040688.