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. Reference Beurskens2021, Reference Beurskens2022), 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. Reference Bähner2021) and Doppler reflectometry (Carralero et al. Reference Carralero2021) 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. Reference Bozhenkov2020; Pablant et al. Reference Pablant2020). 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. Reference Plunk2019), 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. Reference Proll, Helander, Connor and Plunk2012), the so-called ‘Maximum-J’ property (Helander et al. Reference Helander2012; Mackenbach, Proll & Helander Reference Mackenbach, Proll and Helander2022; Sánchez et al. Reference Sánchez, Velasco, Calvo and Mulas2023).
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 Reference Mynick, Pomphrey and Xanthopoulos2010; Nunami, Watanabe & Sugama Reference Nunami, Watanabe and Sugama2013; Xanthopoulos et al. Reference Xanthopoulos, Mynick, Helander, Turkin, Plunk, Jenko, Görler, Told, Bird and Proll2014; Hegna, Terry & Faber Reference Hegna, Terry and Faber2018; Jorge & Landreman Reference Jorge and Landreman2021; Stroteich et al. Reference Stroteich, Xanthopoulos, Plunk and Schneider2022; Jorge et al. Reference Jorge, Dorland, Kim, Landreman, Mandell, Merlo and Qian2023), the goal of the present work is to increase the linear onset of ITG modes (‘critical gradient’) (Roberg-Clark, Plunk & Xanthopoulos Reference Roberg-Clark, Plunk and Xanthopoulos2021), since for configurations with stiff transport, it can determine the radial ion temperature profile (Baumgaertel, Hammett & Mikkelsen Reference Baumgaertel, Hammett and Mikkelsen2013). Even though the nature of the marginally unstable fluctuations may differ between low-magnetic-shear stellarators (Bhattacharjee et al. Reference Bhattacharjee, Sedlak, Similon, Rosenbluth and Ross1983; Zocco et al. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018, Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022) (Floquet-like) and tokamaks (Terry, Anderson & Horton Reference Terry, Anderson and Horton1982; Biglari, Diamond & Rosenbluth Reference Biglari, Diamond and Rosenbluth1989; Romanelli Reference Romanelli1989; Jenko, Dorland & Hammett Reference Jenko, Dorland and Hammett2001) (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 Reference Roberg-Clark, Plunk and Xanthopoulos2022). 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. (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
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
where $\tau = T/(ZT_e)$ with the charge ratio defined as $Z = q/q_e$. The equilibrium distribution is the Maxwellian
and we introduce the velocity-dependent diamagnetic frequency
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
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
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
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
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
where
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]$.
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. Reference Nemov, Kasilov, Kernbichler and Heyn1999) 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 Reference Albert, Kasilov and Kernbichler2020) 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. Reference Mau, Kaiser, Grossman, Raffray, Wang, Lyon, Maingi, Ku and Zarnstorff2008), 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.
Linear gyrokinetic simulations in local flux tube geometry using the GENE code (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000) are performed to determine the critical gradient in the same manner as by Roberg-Clark et al. (Reference Roberg-Clark, Plunk and Xanthopoulos2021, Reference Roberg-Clark, Plunk and Xanthopoulos2022), 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.
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. Reference Talmadge, Anderson, Anderson, Deng, Guttenfelder, Likin, Lore, Schmitt and Zhai2008), another QH optimized configuration (figure 5b). Rotational transform and neoclassical losses for HSX are shown in figure 6.
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 Reference Helander2014). This effect, distinct from flux expansion, is associated with a stabilization effect (Waltz & Boozer Reference Waltz and Boozer1993) 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. Reference Angelino2009; Helander & Plunk Reference Helander and Plunk2021; Plunk & Helander Reference Plunk and Helander2022; Stroteich et al. Reference Stroteich, Xanthopoulos, Plunk and Schneider2022).
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 Reference Sugama1999; Terry, Baver & Gupta Reference Terry, Baver and Gupta2006; Faber et al. Reference Faber, Pueschel, Terry, Hegna and Roman2018), which can form the basis for nonlinear saturation of ITG turbulence (Pueschel et al. Reference Pueschel, Faber, Citrin, Hegna, Terry and Hatch2016; Hegna et al. Reference Hegna, Terry and Faber2018; Mckinney et al. Reference Mckinney, Pueschel, Faber and Hegna2019). 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 Reference Plunk, Xanthopoulos and Helander2017; Mckinney et al. Reference Mckinney, Pueschel, Faber and Hegna2019), 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.
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. Reference Beurskens2021) and in global gyrokinetic simulations of W7-X (Bañón Navarro et al. Reference Bañón Navarro, Di Siena, Velasco, Wilms, Merlo, Windisch, LoDestro, Parker and Jenko2023). 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. (Reference Hegna2022) 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. Reference Bañón Navarro, Roberg-Clark, Plunk, Fernando, Di Siena, Wilms and Jenko2023) has shown that the apparent nonlinear threshold for ITG modes is appreciably larger in HSK than HSX across the plasma volume.
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. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018) 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 Reference Jenko and Dorland2002; Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022). 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. Reference Jenko, Dorland and Hammett2001) 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. 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.
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. Reference Roberg-Clark, Plunk and Xanthopoulos2022). Rather than reducing the magnitudes of both drift curvature ((2.6)) and flux surface compression (Mynick et al. Reference Mynick, Pomphrey and Xanthopoulos2010; Xanthopoulos et al. Reference Xanthopoulos, Mynick, Helander, Turkin, Plunk, Jenko, Görler, Told, Bird and Proll2014), 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. Reference Helander, Bird, Jenko, Kleiber, Plunk, Proll, Riemann and Xanthopoulos2015; Proll et al. Reference Proll, Plunk, Faber, Görler, Helander, McKinney, Pueschel, Smith and Xanthopoulos2022), non-zero density gradients (Thienpondt et al. Reference Thienpondt2023), trapped-electron mode turbulence (Faber et al. Reference Faber, Pueschel, Proll, Xanthopoulos, Terry, Hegna, Weir, Likin and Talmadge2015; Mackenbach et al. Reference Mackenbach, Proll and Helander2022) or finite beta effects, which are generally stabilizing (Pueschel, Kammerer & Jenko Reference Pueschel, Kammerer and Jenko2008; Zocco, Helander & Connor Reference Zocco, Helander and Connor2015). 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 Reference Landreman2022; Landreman & Paul Reference Landreman and Paul2022). 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. Reference Mishchenko, Borchardt, Hatzky, Kleiber, Könies, Nührenberg, Xanthopoulos, Roberg-Clark and Plunk2023), 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. Reference Fujiwara, Kawahata, Ohyabu, Kaneko, Komori and Yamada2001) and stellarators such as TJ-II (De Aguilera et al. Reference De Aguilera2015) and W7-AS (Geiger et al. Reference Geiger, Weller, Zarnstorff, Nührenberg, Werner and Kolesnichenko2004; Weller et al. Reference Weller2006) 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.
Acknowledgements
The authors thank the anonymous reviewers for providing helpful suggestions to improve the paper. We thank C. Nührenberg for checking MHD stability of the HSK configuration, R. Jorge for providing the NEAT code to track alpha particle losses, and M. Landreman and B. Medasani for help with the use of SIMSOPT. We thank P. Helander, A. Zocco and M. Landreman for contributing useful ideas to this work. Computing resources at the Cobra cluster at IPP Garching and the Marconi Cluster were used to perform the simulations.
Editor Alex Schekochihin thanks the referees for their advice in evaluating this article.
Declaration of interests
The author reports no conflict of interest.
Funding
This research was supported by a grant from the Simons Foundation (No. 560651 to G.T.R.-C.). This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200 – EUROfusion). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them.