2. Linear gyrokinetic equation

Following Plunk et al. (Reference Plunk, Helander, Xanthopoulos and Connor2014), we use the standard gyrokinetic system of equations (Brizard & Hahm Reference Brizard and Hahm2007) to describe electrostatic fluctuations destabilized along a thin flux tube tracing a magnetic field line. The ballooning transform (Dewar & Glasser Reference Dewar and Glasser1983) and twisted slicing representation (Roberts & Taylor Reference Roberts and Taylor1965) are used to separate out the fast perpendicular (to the magnetic field) scale from the slow parallel scale. The magnetic field in field-following (Clebsch) representation reads $\boldsymbol {B}=\boldsymbol {\nabla } \psi \times \boldsymbol {\nabla } \alpha$, where $\psi$ is a toroidal flux surface label and $\alpha =\vartheta - \iota \phi$ labels the magnetic field line on the surface, with $q$ the safety factor, $\vartheta$ the poloidal angle and $\phi$ the toroidal angle. The perpendicular wave vector is then expressed as $\boldsymbol {k_{\perp }} = k_{\alpha } \boldsymbol {\nabla } \alpha + k_{\psi } \boldsymbol {\nabla } \psi$, where $k_{\alpha }$ and $k_{\psi }$ are constants, so the variation of $\boldsymbol {k_{\perp }}(l)$ stems from that of the geometric quantities $\boldsymbol {\nabla } \alpha$ and $\boldsymbol {\nabla } \psi$, with $\ell$ the field-line-following (arc length) coordinate.

We assume Boltzmann-distributed (adiabatic) electrons, thus solving for the perturbed ion distribution $g_{i}({v_{\parallel }},{v_{\perp }},\ell,t)$, defined to be the non-adiabatic part of $\delta f_{i}$ ($\delta f_{i}=f_{i}-f_{i0})$ with $f_{i}$ the ion distribution function and $f_{i0}$ a Maxwellian. The electrostatic potential is $\phi (\boldsymbol {\ell })$, and ${v_{\parallel }}$ and ${v_{\perp }}$ are the particle velocities parallel and perpendicular to the magnetic field, respectively.

The linear gyrokinetic equation for the ions is written

(2.1)\begin{equation} i{v_{{\parallel}}} \frac{\partial g}{\partial \ell} + (\omega - \tilde{\omega}_d)g = \varphi J_0(\omega - \widetilde{\omega}_*)f_0,\end{equation}

with the following definitions: $J_0 = J_0(k_{\perp }v_{\perp }/\varOmega ) = J_0(k_{\perp }\rho \sqrt {2}v_{\perp }/v_{{T}})$; the ion thermal velocity is $v_{{T}} = \sqrt {2T/m}$ and the thermal ion Larmor radius is $\rho = v_{{T}}/(\varOmega \sqrt {2})$; $n$ and $T$ are the background ion density and temperature; $q$ is the ion charge; $\varphi = q\phi /T$ is the normalized electrostatic potential; $\varOmega =q B/m$ is the cyclotron frequency, with $B=|\boldsymbol {B}|$ the magnetic field strength. Assuming Boltzmann electrons, the quasi-neutrality condition is

(2.2)\begin{equation} \int d^3{\boldsymbol v} J_0 g = n(1 + \tau) \varphi,\end{equation}

where $\tau = T/(ZT_e)$ with the charge ratio defined as $Z = q/q_e$. The equilibrium distribution is the Maxwellian

(2.3)\begin{equation} f_0 = \frac{n}{(v_{{T}}^2\pi)^{3/2}}\exp({-}v^2/v_{{T}}^2), \end{equation}

and we introduce the velocity-dependent diamagnetic frequency

(2.4)\begin{equation} \widetilde{\omega}_* = \omega_*^{\mathrm{T}} \left[\frac{v^2}{v_{{T}}^2} - \frac{3}{2}\right], \end{equation}

where we neglect background density variation and define $\omega _*^{\mathrm {T}} = (Tk_{\alpha }/q)d\ln T/d\psi$. The magnetic drift frequency is $\tilde {\omega }_d = {\boldsymbol v}_d\boldsymbol {\cdot } {\boldsymbol k}_{\perp }$ and the magnetic drift velocity is ${\boldsymbol v}_d = \widehat {\boldsymbol b}\times ((v_{\perp }^2/2)\boldsymbol {\nabla } \ln B + {v_{\parallel }}^2\boldsymbol {\kappa })/\varOmega$, where $\boldsymbol {\kappa } = \widehat {\boldsymbol b}\boldsymbol {\cdot } \boldsymbol {\nabla }\widehat {\boldsymbol b}$. We take $\boldsymbol {\nabla } \ln B = \boldsymbol {\kappa }$, the zero $\beta$ approximation, for simplicity. We then let

(2.5)\begin{equation} \tilde{\omega}_d = \frac{\boldsymbol{k_{{\perp}}} \boldsymbol{\cdot} (\hat{\boldsymbol{b}} \times \boldsymbol{\kappa})v^{2}_{T}}{\varOmega} \left[\frac{{v_{{\parallel}}}^2}{v_{{T}}^2} + \frac{{v_{{\perp}}}^2}{2v_{{T}}^2}\right] = \omega_d(\ell) \left[\frac{{v_{{\parallel}}}^2}{v_{{T}}^2} + \frac{{v_{{\perp}}}^2}{2v_{{T}}^2}\right], \end{equation}

where the velocity-independent drift frequency $\omega _d(\ell )$ generally varies along the field line. Positive values of $\omega _d$ correspond to ‘bad’ curvature, i.e. are destabilizing for ITG modes, assuming $\omega _*^{\mathrm {T}}$ is negative. Assuming $k_{\psi }=0$ for simplicity, we then define the drift curvature $K_{d}$ by writing

(2.6)\begin{equation} \omega_d(\ell) \propto K_{d}(\ell) \equiv a^2{\boldsymbol{\nabla}}\alpha \boldsymbol{\cdot} \hat{\boldsymbol{b}} \times \boldsymbol{\kappa}, \end{equation}

where $K_{d}$ contains the purely geometric variation of the drift frequency and $a$ is the minor radius of the flux surface at the edge. For the purpose of analysing gyrokinetic simulation results, we define the metrics $g^{yy}=a^{2}s_{0}/(q_{0}^{2})( \boldsymbol {\nabla } \alpha )^{2}$, with $s=\psi /\psi _{edge}$ the toroidal flux normalized to its value at the last closed flux surface ($q_{0}$ is the safety factor on a particular surface $s=s_{0}$), and $g^{xx}=a^{2}/(4s_{0})(\boldsymbol {\nabla } s)^{2}$. Finally, we define the poloidal wavenumber $k_{y}=(q_{0}/\sqrt {s_{0}})k_{\alpha }$.

We neglect particle trapping so ${x_{\parallel }}$ and ${x_{\perp }}$ do not depend on $\ell$. Most evidence indicates that the ITG mode uniformly responds to changes in the parameters $\tau$ and $\mathrm {d}n/\mathrm {d}\psi$, and is stabilized by increasing either of them, except for a relatively small region of parameter space where positive density gradients can destabilize the mode. For simplicity, we thus set $\tau =1$ and $\boldsymbol {\nabla } n=0$.

2.1. Estimating the ITG mode critical gradient

As reported by Roberg-Clark et al. (Reference Roberg-Clark, Plunk and Xanthopoulos2022), the geometric dependence of the linear ITG critical gradient can be estimated in a simple way. The drift curvature profile ((2.6)) on a particular magnetic field line (at radial location $s = s_{0}$ and field line $\alpha =\alpha _{0}$) is fitted with quadratic curves in regions of bad curvature, which we refer to as ‘drift wells.’ The fitting acts as an effective coarse-graining of the geometry, producing a smoothed curvature profile as mentioned in the preceding section. Using the fitted profile, a value of the predicted critical gradient is produced, corresponding to the drift well with the smallest critical gradient. The formula reads

(2.7)\begin{equation} F_{\rm crit}=\frac{a}{L_{T,{\rm crit}}}=2.66\left(\frac{a}{R_{\mathrm{eff}}}\right), \end{equation}

with $a/R_{\mathrm {eff}}$ the peak value of the drift curvature within the drift well. Here we have ignored the parallel stabilizing term included by Roberg-Clark et al. (Reference Roberg-Clark, Plunk and Xanthopoulos2022) as we intend to produce a configuration with a large critical gradient solely by large ‘drift curvature’; see the later discussion of geometric interpretation in § 4.1. For a finite parallel stabilization effect to take hold for the critical gradient, an unusually large amount of global shear and a relatively small drift curvature are required, thus interfering with the main effect explored here.

An alternative measure would be to find the threshold gradient relative to the major radius $R$. By focusing on $R/R_{\mathrm {eff}}$ and holding $a$ fixed, we would explore a different problem, which is how much the shaping at fixed aspect ratio can increase the strength of the curvature. In the interest of increasing core ion temperatures over the scale of the minor radius, however, $a/L_{T}$ is a useful measure. $R/R_{\mathrm {eff}}$ is a parameter which is representative of the curvature, but the ‘additional’ factor of $a/R$ in our definition is also motivated on physical grounds: lower aspect ratios improve the confinement time via the gyro-Bohm scaling of heat fluxes $\tau _{1} / \tau _{2} \sim (A_{1}/A_{2})(Q_{GB,1}/Q_{GB,2})$. Note that this ratio is arrived at assuming the surface area of each configuration goes like $aR$, with the $1/a^2$ factor coming from gyro-Bohm scaling, such that different surface areas are accounted for in this estimate. We are thus compelled to find compact configurations, as described in the following section.

3. Optimization method

We use the SIMSOPT software framework (Landreman et al. Reference Landreman, Medasani, Wechsung, Giuliani, Jorge and Zhu2021) to generate a vacuum stellarator configuration with quasi-helical (QH) symmetry and a significant linear ITG critical gradient as a result of drift curvature. The stellarator magnetic field is described by a boundary surface given in the Fourier representation

(3.1)\begin{gather} R(\vartheta,\phi)=\sum_{m,n} R_{m,n} \cos(m\vartheta-4 n\phi) , \end{gather}

(3.2)\begin{gather}Z(\vartheta,\phi)=\sum_{m,n} Z_{m,n} \sin(m\vartheta-4 n\phi), \end{gather}

where we have assumed stellarator symmetry and have set the number of field periods to $n_{fp}=4$. Global vacuum solutions are constructed at each iteration by running the VMEC (Hirshman & Whitson Reference Hirshman and Whitson1983) code, which solves the MHD equations using an energy-minimizing principle.

Optimization proceeds by treating the Fourier coefficients in (3.1) as parameters and varying them to find a least-squares minimization (using a trust region) of the specified objective function, which reads

(3.3)\begin{equation} f = f_{QS} + (A-4.10)^{2} + (F_{\rm crit} - 2.00 )^2, \end{equation}

where

(3.4)\begin{equation} f_{QS} = \sum_{s_{j}}\left\langle \left( \frac{1}{B^{3}}[({-}1-\iota)\boldsymbol{B}\times \boldsymbol{\nabla} B \boldsymbol{\cdot} \boldsymbol{\nabla} \psi - G\boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} B] \right)^{2} \right\rangle \end{equation}

is the quasi-symmetry residual taken from (Landreman & Paul Reference Landreman and Paul2022), with $A=R/a$ the aspect ratio output by VMEC, and $a/L_{T,{\rm crit}}$ as defined in (2.7) with the field line $(s_{0}=0.5$, $\alpha _{0}=0)$ chosen. Here $G(\psi )$ is $\mu _{0}/(2\pi )$ times the poloidal current outside the surface and we have set the toroidal current to zero, while $\langle \boldsymbol {\cdot } \rangle$ denotes a flux surface average over the individual surfaces $s=s_{j}$. We choose $B=B(\theta +4\phi _{B})$, with $\theta$ and $\phi _{B}$ the Boozer poloidal and toroidal angles, corresponding to quasi-helical symmetry with a negative axis helicity. The input equilibrium for the optimization is a ‘warm start’ example file included in SIMSOPT with poor quasi-helical symmetry, $n_{fp}=4$ and aspect ratio $A=R/a=7$, with $R$ the major radius and $a$ the minor radius, and an initial $a/L_{T,{\rm crit}}=0.50$. The target $A=4.10$ was chosen to enhance the drift curvature $K_{d}$, and to take into account that pushing the aspect ratio below the number of field periods is difficult to achieve while maintaining finite rotational transform in a helical stellarator (Roberg-Clark et al. Reference Roberg-Clark, Plunk and Xanthopoulos2021). We speculate $a/L_{T,{\rm crit}}=2.00$ to be roughly the maximum achievable linear critical gradient in the absence of very large global shear stabilization (Roberg-Clark et al. Reference Roberg-Clark, Plunk and Xanthopoulos2022). We calculate the residual $f_{QS}$ on the surfaces $s=[0,0.1,0.2,0.3,0.4,0.5]$, thus leaving the region $0.5< s<1.0$ free. We make this choice knowing that requiring a ‘precise’ degree of quasi-symmetry in the entire volume often has the effect of reducing the global magnetic shear in the configuration to vanishingly small values (Landreman & Paul Reference Landreman and Paul2022). Similarly to Landreman & Paul (Reference Landreman and Paul2022), the number of Fourier coefficient parameters are increased in a series of four steps, with the toroidal and poloidal mode numbers in the VMEC calculation chosen to be $m_{\rm pol}=n_{\rm tor}=[3,5,6]$.

5. Balancing MHD and ITG stability

HSK has a vacuum magnetic hill so the comparison to a reactor with finite plasma pressure is theoretical. However, it is revealing to see the extent to which ITG turbulent losses can be quenched by abandoning MHD stability in the optimization. This trade-off was hinted at in previous reactor studies carried out for QH stellarators (Bader et al. Reference Bader, Drevlak, Anderson, Faber, Hegna, Likin, Schmitt and Talmadge2019, Reference Bader2020). We can gain some intuition for the lack of MHD stability in HSK by comparing with the standard ‘bean-shaped’ cross-section in W7-X, which has an indentation on the inboard, and vertical, compressed surfaces on the outboard. The indentation aids the formation of a vacuum magnetic well by strongly reducing the volumetric expansion of the surfaces such that $V^{\prime \prime }(\psi )<0$ (Cooper Reference Cooper1992). The outboard of the bean section has a small drift curvature, a result of large values of $g^{\psi \psi }$ as well as a small geometric curvature (see discussion of the QIPC stellarator (Subbotin et al. Reference Subbotin2006; Beidler et al. Reference Beidler2011), which has a similar bean-shaped section to that of W7-X). Quantitatively, we also know from near-axis theory that the magnitude of the drift curvature contributes to the formation of an unstable magnetic hill (Landreman & Jorge Reference Landreman and Jorge2020), and that the expression for Mercier stability (governing large-$n$ ballooning modes near rational surfaces) in general contains destabilizing terms proportional to $1/|\boldsymbol {\nabla } \psi |$ (Landreman & Jorge Reference Landreman and Jorge2020). Thus the bean cross-section, possessing large values of $|\boldsymbol {\nabla } \psi |$, a small geometric curvature and $V^{\prime \prime }(\psi )<0$, is generally stable to MHD modes.

In contrast, HSK lacks an indentation and has expanded surfaces on the outboard side, leading to a magnetic hill and significant instability with respect to the Mercier criterion. However, it is plausible that stability of a configuration like HSK could be enhanced through the addition of an indentation and a weakening of the drift curvature on the outboard side. Optimization studies exploring this compromise between MHD and ITG stability are currently underway.