Reduction of electrostatic turbulence in a quasi-helically symmetric stellarator via critical gradient optimization

Abstract

We present a stellarator configuration optimized for a large threshold (‘critical gradient’) for the onset of the ion temperature gradient (ITG) driven mode, which achieves the largest critical gradient we have seen in any stellarator. Above this threshold, gyrokinetic simulations show that the configuration has low turbulence levels over an experimentally relevant range of the drive strength. The applied optimization seeks to maximize the drift curvature, leading to enhanced local-shear stabilization of toroidal ITG modes, and the associated turbulence. These benefits are combined with excellent quasi-symmetry, yielding low neoclassical transport and vanishingly small alpha particle losses. Analysis of the resulting configuration suggests a trade-off between magnetohydrodynamic (MHD) and ITG stability, as the new configuration possesses a vacuum magnetic hill.

1. Introduction

The excitation of the ion temperature gradient (ITG) mode in magnetically confined fusion devices leads to turbulence that is responsible for energy losses, which reduce the plasma confinement needed for potential fusion performance. For instance, it has been argued that during operation of the Wendelstein 7-X (W7-X) stellarator in electron heating scenarios, the ITG mode leads to so-called ‘ion temperature clamping’ (Beurskens et al. 2021, 2022), thus preventing the heating of ions in the plasma core above 2 keV. While resolving the location and extent of the fine-scale (ion Larmor radius) fluctuations is difficult in experiments, implementation of diagnostics using phase contrast imaging (Bähner et al. 2021) and Doppler reflectometry (Carralero et al. 2021) has helped to characterize the dependence of W7-X experimental temperature and density profiles on such fluctuations.

To lessen the negative effects of the ITG mode, a possible strategy to follow involves the manipulation of the density and electron profiles, as implemented in W7-X using pellet injections (Bozhenkov et al. 2020; Pablant et al. 2020). A stellarator magnetic field may also be shaped to reduce microturbulence losses during steady-state operation. In the case of electron-temperature gradient turbulence, it is argued that multiple field periods are stabilizing via reduction of the parallel connection length (Plunk et al. 2019), while trapped-electron mode turbulence can be ameliorated by requiring that the parallel adiabatic invariant J achieves its maximum on the magnetic axis (Proll et al. 2012), the so-called ‘Maximum-J’ property (Helander et al. 2012; Mackenbach, Proll & Helander 2022; Sánchez et al. 2023).

While most previous works have exploited magnetic field shaping to lower the rate of the ITG transport scaling (‘stiffness’) (see, e.g. Mynick, Pomphrey & Xanthopoulos 2010; Nunami, Watanabe & Sugama 2013; Xanthopoulos et al. 2014; Hegna, Terry & Faber 2018; Jorge & Landreman 2021; Stroteich et al. 2022; Jorge et al. 2023), the goal of the present work is to increase the linear onset of ITG modes (‘critical gradient’) (Roberg-Clark, Plunk & Xanthopoulos 2021), since for configurations with stiff transport, it can determine the radial ion temperature profile (Baumgaertel, Hammett & Mikkelsen 2013). Even though the nature of the marginally unstable fluctuations may differ between low-magnetic-shear stellarators (Bhattacharjee et al. 1983; Zocco et al. 2018, 2022) (Floquet-like) and tokamaks (Terry, Anderson & Horton 1982; Biglari, Diamond & Rosenbluth 1989; Romanelli 1989; Jenko, Dorland & Hammett 2001) (Toroidal-like), it seems that the size of the so-called ‘drift curvature’, a factor appearing in the gyrokinetic equation, can be used to predict the critical gradient for most optimized stellarator configurations (Roberg-Clark, Plunk & Xanthopoulos 2022). Here we leverage this predictive capability to generate a new quasi-helically symmetric configuration with low neoclassical transport and large ITG critical gradient. It turns out that turbulent losses from ITG modes above this threshold are suppressed when compared with a well-known quasi-helically symmetric configuration, in particular near the critical gradient. We attribute this enhanced stability, at least in part, to the high onset gradient of localized, toroidal ITG modes, which are further damped by flux expansion and local shear effects.

The paper is structured as follows. In § 2, we define the linear gyrokinetic system used to analyse ITG mode properties. Section 3 describes the optimization procedure used to generate the new configuration, while § 4 presents results of the optimization, along with a description of the local shear effect. We conclude with an outlook on the apparent competition between ITG and MHD stability in § 5.

2. Linear gyrokinetic equation

Following Plunk et al. (2014), we use the standard gyrokinetic system of equations (Brizard & Hahm 2007) to describe electrostatic fluctuations destabilized along a thin flux tube tracing a magnetic field line. The ballooning transform (Dewar & Glasser 1983) and twisted slicing representation (Roberts & Taylor 1965) 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. (2022), 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. (2022) 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. 2021) 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 1983) 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 2022), 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. 2021). 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. 2022). 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 2022). Similarly to Landreman & Paul (2022), 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]$.

4. Results

Figure 1 shows a surface plot of the outermost flux surface of the resulting optimized configuration, which we dub ‘HSK’ (helically symmetric kompakt stellarator). The contours of $B$ are plotted in the Boozer angle plane at $s=0.25$ in figure 2(a) indicating the quasi-helical symmetry. The highly compact configuration has roughly triangular cross-sections (cuts at constant toroidal angle; figure 3) and noticeably lacks a ‘bean-shaped’ cross-section at toroidal angle $\phi =0$, which most optimized stellarator configurations possess. The spatial separation between the flux surfaces is substantial, indicating low surface compression, while the lack of an indentation on the inboard side at $\phi =0$ appears to be a result of not optimizing for a vacuum magnetic well. The configuration instead has a sizeable magnetic hill, i.e. $V^{\prime \prime }(\psi )$, the second derivative of the surface volume, is positive at all radii. Heuristically speaking, this result is expected from configurations with average bad curvature, i.e. the average value of the drift curvature over the entire surface is positive. HSK has a rotational transform varying from $\iota = 1.3$ to $1.8$ (demonstrating significant global magnetic shear; figure 2b) and values of the neoclassical transport coefficient $\epsilon _{\rm eff}$ (Nemov et al. 1999) ranging from $0.4\,\%$ on axis to $1.6\,\%$ at the last closed flux surface (figure 2c), implying excellent neoclassical confinement. Using a symplectic integrator for guiding centre trajectories (Albert, Kasilov & Kernbichler 2020) in this configuration but rescaled to the same average $B$ and minor radius (5.3 Tesla and 1.85 metres, respectively) as the ARIES-CS reactor (Mau et al. 2008), we find no collisionless alpha particle losses after 10 ms for particles launched from $s=0.25$, suggesting excellent quasi-symmetry (figure not shown). The test particle orbits were simulated in vacuum with no coils included in the final equilibrium, which would presumably affect the results.

Figure 1. Boundary surface of HSK.

Figure 2. Properties of HSK. (a) Contours of $B$ with the trajectory of a field line (black curve near the centre) overlaid in Boozer coordinates. (b) Rotational transform as a function of radius. (c) Neoclassical transport coefficient $\epsilon _{\rm eff}$ as a function of radius.

Figure 3. Surface cuts of HSK taken at constant toroidal angle, $\phi =\text {const}$.

Linear gyrokinetic simulations in local flux tube geometry using the GENE code (Jenko et al. 2000) are performed to determine the critical gradient in the same manner as by Roberg-Clark et al. (2021, 2022), i.e. by reducing the applied temperature gradient until a single, marginally unstable mode remains. An unusually large resolution in $\mu$ (i.e. ${v_{\perp }}$ space) of $n_{\mu }=32$ was required for numerical convergence of the critical gradient, likely because of the extreme values of curvature near the outboard midplane that imply significant linear phase mixing in ${v_{\perp }}$ through the $\boldsymbol {\nabla } B$ drift term in (2.5). The linear ITG critical gradient at $(s_{0}=0.5, \alpha _{0}=0)$ ($k_{y}\rho =0.3$), while not reaching $2$ as predicted by the fitting model ((2.7)), nonetheless attains the value of $a/L_{T,{\rm crit}}=1.75$, the highest we have seen in any stellarator (figure 4). Thus the simultaneous optimization for quasi-helical symmetry, aspect ratio and ITG linear critical gradient was successful. We discuss how such a large $R_{\mathrm {eff}}$ is achieved in the next section.

Figure 4. Simulation of critical gradients versus model predictions on a log-log scale for flux tubes taken from several stellarator geometries, adapted from figure 1 of Roberg-Clark et al. (2022), with the data point for HSK added. The field line $(s_{0}=0.5,\alpha _{0}=0)$ was chosen for all configurations. (a) Log-log scale and (b) linear scale.

One might reasonably expect this particular ITG optimization strategy (increasing bad curvature) to heavily exacerbate linear growth rates, and thus nonlinear transport, of ITG modes above the linear critical gradient, implying a trade-off between linear and nonlinear stability. Further results in the next section, however, oppose this intuition, revealing significant ITG stability above the critical gradient.

4.1. Mechanisms of ITG turbulence suppression

One means by which a stellarator configuration can attain large values of the drift curvature $K_{d}$ ((2.6)) is by compression of the $\alpha$ coordinate, i.e. large $|\boldsymbol {\nabla } \alpha |$. This is indeed observed for HSK, in the metrics following the magnetic field line at $(s_{0}=0.5,\alpha _{0}=0)$, where we plot $g^{yy} \propto g^{\alpha \alpha }=a^{2}|\boldsymbol {\nabla } \alpha |^{2}$ as a function of arc length $\ell$, which attains a value on the outboard midplane ($\ell =0$) larger than $2$ (figure 5a). Note that this location at the outboard midplane corresponds to $\phi =0$ in figure 1(b). For comparison, we also plot the metrics along the same field line in the HSX stellarator (Talmadge et al. 2008), another QH optimized configuration (figure 5b). Rotational transform and neoclassical losses for HSX are shown in figure 6.

Figure 5. Metrics (defined near the end of § 2) as a function of the field-line-following coordinate $\ell$ at the flux tube location $(s_{0}=0.5,\alpha _{0}=0)$, where HSK was optimized for a high critical gradient. The flux tubes were constructed with two poloidal turns. The vertical axis is in dimensionless units, such that the metrics for each configuration can be compared directly. (a) HSK; (b) HSX.

Figure 6. Rotational transform profile and neoclassical losses (effective ripple) for the HSX configuration used.

Typically in a stellarator, the vector $\boldsymbol {\nabla } \alpha$ also develops a substantial component in the direction parallel to $\boldsymbol {\nabla }\psi$, attributed to so-called ‘local shear’ of the magnetic field lines (Helander 2014). This effect, distinct from flux expansion, is associated with a stabilization effect (Waltz & Boozer 1993) due to $k_{\perp } \rho _i$ increasing, i.e. due a finite Larmor radius (FLR). Global shear may further amplify this damping through secular increase of $g^{yy}$ along the field line, although we suspect this is not the dominant effect in the case of HSK. We see that both of these effects are at play with HSK, but note for the case of HSX that the pattern near $\ell =0$ in the metrics is effectively inverted, i.e. $|K_{da}| \ll 1$, $g^{yy}<1$ and $g^{xx}>1$, which is typical of most optimized stellarator configurations.

In a simple limit, when $\boldsymbol {\nabla }\alpha \boldsymbol {\cdot } \boldsymbol {\nabla } \psi \approx 0$, and assuming weak variation of the overall magnetic field strength $B = |\boldsymbol {\nabla } \alpha \times \boldsymbol {\nabla } \psi |$, the $\alpha$-compression is directly linked to reduction of $g^{xx}\propto g^{\psi \psi }$ (i.e. $g^{xx} \sim 1/g^{yy}$). In such cases, reduced instability growths and turbulent transport can be intuitively explained via expansion of flux surfaces, resulting in an effective suppression of the applied temperature gradient, since the physical temperature gradient $|\boldsymbol {\nabla } T|={\rm d}T/{\rm d}\psi |\boldsymbol {\nabla } \psi |$ scales with $|\boldsymbol {\nabla } \psi |$ (Angelino et al. 2009; Helander & Plunk 2021; Plunk & Helander 2022; Stroteich et al. 2022).

Another factor which controls $R_{\mathrm {eff}}$ is the curvature of the magnetic field line itself. We find, for HSX and HSK at $s=0.5, \alpha =0$, that the size of the normal curvature $\kappa _{\text {norm}} = \kappa \boldsymbol {\cdot } \boldsymbol {\nabla } \psi / |\boldsymbol {\nabla } \psi |$, the dominant component at the outboard midplane, is roughly the inverse major radius $1/R$. Thus the amplification of $R_{\mathrm {eff}}$ relative to $a$ in HSK comes predominantly from the aspect ratio $a/R$ and from $|\boldsymbol {\nabla } \alpha |$. Multiplying each respective critical gradient by their aspect ratio to renormalize $R_{\mathrm {eff}}$ to the major radius, we find $R/L_{T,\text {crit}}$ is a factor of two larger in HSK compared with HSX.

We now carry out linear gyrokinetic simulations comparing HSK with HSX using the GENE code, but this time above marginal stability. We choose the flux tube $(s_{0}=0.25,\alpha _{0}=0)$ to check that the benefits of HSK are not confined to the field line optimized for ITG stability at $(s_{0}=0.5,\alpha _{0}=0)$. Figure 7(a) shows a strikingly narrow range of unstable wavenumbers for HSK in $k_{y}\rho$ at the gradient $a/L_{T}=2$. Much of the analysis of electrostatic modes in QH stellarators has focused on damped or subdominant eigenmodes (Sugama 1999; Terry, Baver & Gupta 2006; Faber et al. 2018), which can form the basis for nonlinear saturation of ITG turbulence (Pueschel et al. 2016; Hegna et al. 2018; Mckinney et al. 2019). We interpret the narrow linear growth rate spectrum of HSK at $a/L_{T}=2$ to mean that most of the linear eigenmode spectrum is stable at this gradient. It is worth mentioning that HSX sometimes has a favourable scaling of heat flux at these temperature gradients compared with other configurations (Plunk, Xanthopoulos & Helander 2017; Mckinney et al. 2019), despite larger linear growth rates. HSX was also not optimized for reduced turbulent transport, so one should not expect it to outperform turbulence-optimized configurations.

Figure 7. ITG mode linear growth rates as a function of $k_{y} \rho _{i}$ for (a) HSK and (b) HSX for the temperature gradients $a/L_{T}=2,3$ and $4$, at the flux tube location $(s_{0}=0.25,\alpha _{0}=0)$. Note the much smaller range of values for HSK.

A nonlinear gyrokinetic analysis with GENE using the same flux tubes shows that the narrow and small linear growth rates in HSK are complemented by a stark reduction in heat fluxes, a truer figure of merit for confinement. In comparison with HSX, we find heat fluxes are smaller by roughly a factor of $4$ at $a/L_{T}=2$ (figure 8b). While the reduction is not as strong at higher gradients, the curve for HSK remains below that of HSX for the full range of gradients studied. The effect of the optimization is most potent near the critical gradient, $a/L_{T} \simeq 2 - 3$. This range of gradients is likely the most relevant to experiments, as seen in the data for (Beurskens et al. 2021) and in global gyrokinetic simulations of W7-X (Bañón Navarro et al. 2023). For the same heating power, we can thus expect HSK to have a significantly larger temperature gradient compared with HSX, leading to improved confinement and higher core ion temperatures. Furthermore, HSK has the advantage of a reduced aspect ratio, which implies a favourable confinement time compared with larger aspect ratio devices with the same transport coefficients. Note that the turbulence-optimized WISTELL-C QH configuration presented by Hegna et al. (2022) is likely compared with HSX at the radial location $s=0.5$, achieving a factor of 2 reduction in heat fluxes at $a/L_{T}=4$. The same reduction factor for $a/L_{T}=4$ at $s_{0}=0.25$ is found when comparing HSK and HSX (figure 8b). While we only study the location $s = 0.25$ here, a global gyrokinetic transport analysis of HSK (Bañón Navarro et al. 2023) has shown that the apparent nonlinear threshold for ITG modes is appreciably larger in HSK than HSX across the plasma volume.

Figure 8. (a) Nonlinear heat fluxes computed by GENE at the flux tube $(s_{0}=0.25,\alpha _{0}=0)$ for HSK and HSX, varying the applied temperature gradient. (b) Same plot as in panel (a) but with the vertical axis on a log scale.

While bad curvature has markedly increased in HSK relative to HSX, stiffness of ITG transport has not worsened for the range $2 \leq a/L_{T} \leq 4$. We attribute this surprising result to the fact that the onset of localized, toroidal ITG modes (Zocco et al. 2018) appears to have been increased by the combined effect of bad curvature and parallel stabilization by FLR damping, pushing it beyond the critical gradient $a/L_{T,{\rm crit}}=1.75$. It has long been suspected that these modes, rather than extended Floquet or slab-like modes, cause the worst ITG transport (Jenko & Dorland 2002; Zocco et al. 2022). The difference in these thresholds likely explains the apparent discrepancy between the linear (figure 4) and nonlinear (figure 8) results. The difference in absolute threshold between HSX and HSK is large, as seen in the linear results, but the difference in nonlinear thresholds is less stark, probably as a result of closer toroidal ITG mode thresholds. Further nonlinear stabilization of HSX, as mentioned above, likely plays a role. The toroidal threshold could be modelled by an additional stabilization term in the expression (2.7) (Jenko et al. 2001) and will be explored in future work.

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. 2019, 2020). 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 1992). 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. 2006; Beidler et al. 2011), 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 2020), 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 2020). 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.

5.1. Discussion

Our optimization strategy for stellarators seems to have simple geometric consequences. The near-ubiquitous ‘bean-shape’ cross-section, which is thought to impart significant MHD stability and is often the site of the most detrimental ITG turbulence, can be modified to acquire a more triangular shape with a point at the outboard midplane via increasing the gradient of the binormal coordinate, $|\boldsymbol {\nabla } \alpha |$. The resulting increase in ‘bad’ curvature leads to improved linear ITG mode critical gradients, while heat transport at gradients near this threshold is significantly reduced. The feared trade-off between large critical gradients and stiffness of the transport above those thresholds appears not to be a true impediment, as also hinted at by the compact W7-K configuration (Roberg-Clark et al. 2022). Rather than reducing the magnitudes of both drift curvature ((2.6)) and flux surface compression (Mynick et al. 2010; Xanthopoulos et al. 2014), one can increase drift curvature and local shear while improving both the ITG critical gradient and the stiffness of near-marginal transport. The reason appears to be that the critical gradient of localized, toroidal ITG modes is also increased, shielding the configuration from the most detrimental transport losses.

We note that the current work has not taken micro-turbulence effects into account such as kinetic electron physics of ITG modes (Helander et al. 2015; Proll et al. 2022), non-zero density gradients (Thienpondt et al. 2023), trapped-electron mode turbulence (Faber et al. 2015; Mackenbach et al. 2022) or finite beta effects, which are generally stabilizing (Pueschel, Kammerer & Jenko 2008; Zocco, Helander & Connor 2015). However, our goal here is to study a worst-case scenario for ITG modes and how it can be improved from purely geometric considerations, while also ensuring good quasi-symmetry. Our optimization strategy and gyrokinetic analysis are local, specific to certain flux tubes, and future work will address ITG stability on the entire surface, taking more locations into account to ensure the optimization succeeds globally.

The most salient compromise to emerge from this work is that between ITG and MHD stability, and future designs will likely have to prioritize which type of stability is most important. We note that a similar optimization trade-off appears to exist between quasi-symmetry and MHD stability measures such as the vacuum magnetic well (Landreman 2022; Landreman & Paul 2022). There are promising signs that Mercier stability does not have to be rigidly adhered to, e.g. in global gyrokinetic simulations of kinetic ballooning modes (Mishchenko et al. 2023), for heat fluxes at finite beta to match those of traditional MHD-optimized stellarators such as W7-X. We plan to study the stability of HSK to kinetic ballooning modes using this approach to see how detrimental electromagnetic turbulence can be when Mercier stability is strongly violated. Furthermore, there is experimental evidence from the LHD heliotron (Fujiwara et al. 2001) and stellarators such as TJ-II (De Aguilera et al. 2015) and W7-AS (Geiger et al. 2004; Weller et al. 2006) that Mercier-unstable and/or magnetic hill configurations can operate at relatively large $\beta$ values with no serious loss of confinement, with the caveat that the LHD results may be restricted to low density operational regimes. A wide range of stellarator configurations, straddling the line between rigorous MHD stability and strongly suppressed ITG turbulence, thus appears realizable.