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

(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 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)

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

(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. 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

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

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

(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 (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.