Optimizing the HSX stellarator for microinstability by coil-current adjustments

The optimization of helically symmetric experiment (HSX) for reduced microinstability has been achieved by examining a large set of configurations within a neighborhood of the standard operating configuration. This entailed generating a database of more than 106 magnetic-field configurations for HSX by varying the currents in external coils. Using a set of volume-averaged metrics and gyrokinetic simulations, this database has helped to identify a set of configurations that can be used to regulate trapped-electron-mode stability in HSX. This set of configurations is also found to correlate flux-surface elongation and triangularity with an increase in magnetic-well depth, an increase in rotational transform, and low neoclassical heat-flux relative to the standard quasi-helically-symmetric configuration. These results demonstrate sensitivity of plasma behavior in response to changes in a 3D magnetic field to both neoclassical and gyrokinetic models, and the experimental potential in HSX to explore turbulence optimization. This perturbative optimization approach is not unique to HSX, and can readily be deployed on existing fusion devices to identify novel magnetic-fields to be used in turbulence-optimization experiments.


Introduction
One important aspect towards the development of a magnetically confined fusion power plant is the reduction of heat losses arising from transport perpendicular to the confining * Author to whom any correspondence should be addressed.
Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. magnetic field. In stellarators with quasi-symmetric fields, the cross-field transport induced by neoclassical effects has been reduced to levels comparable to tokamaks in existing optimized stellarator experiments [1][2][3] and predicted to be reduced in new stellarator designs [4,5]. As a consequence, transport due to turbulent fluctuations has become an important focus of study in recent years [6][7][8]. One advantage of the stellarator concept is that the 3D shaping of the magnetic field provides a vast optimization space for targeting configurations with reduced turbulent heat losses [9][10][11][12]. However, to date, turbulence-optimized configurations have been obtained primarily for idealized or synthetic geometries, and no such configurations have been validated against experiments. The goal of the present work is to identify experimental candidate configurations that have reduced microinstability growth rates by modifying currents in the existing coil set of the helically symmetric experiment (HSX) stellarator [13]. This lays the groundwork for experimental validation of turbulence optimization in stellarator design.
In the case of HSX, the standard quasi-helically-symmetric magnetic-field configuration, referred to as QHS, has been shown to reduce neoclassical transport [2,3] to levels at which turbulent fluctuations driven by trapped-electron modes (TEMs) dominate heat losses [14,15]. Moreover, a TEM optimization study has been performed in which a reduction in the TEM transport was achieved in simulations by modifying the magnetic geometry [16]. This reduction, however, came at the expense of the quasi-helical symmetry and is not experimentally realizable with the existing coil set on HSX. Regarding the ion-temperature gradient (ITG) mode, the experimentally observed gradients do not exceed the necessary critical gradient for destabilization [17,18]. Therefore ITGs are neglected in the present analysis.
In this article, a set of experimental candidate configurations based on the HSX magnetic-field coils are identified that, when compared numerically to QHS, simultaneously exhibit comparable or even lower neoclassical heat fluxes, increased magnetic-well depth, avoid large stochastic regions of the magnetic field, and reduce TEM growth rates at long wavelength. These effects are all shown to correlate with increased flux-surface elongation and triangularity, which provides an experimentally relevant set of parameters that can be used to investigate the possibility of regulating TEM transport in future HSX discharges.
To identify these correlations, an optimization approach has been implemented in which a database of >10 6 magnetic field geometries was produced by adjusting individual coil currents within 10% of the default operating current. These configurations were then filtered using a set of volume-averaged metrics intended to predict modifications in the TEM behavior while preserving the helical symmetry of the QHS configuration. While this optimization procedure is more computationally expensive than alternate approaches [4,19,20], it does present several key advantages. First, it enables exploration of trends in configuration space by uniformly sampling a domain around a fixed operating point, allowing one to identify properties that might otherwise be obscured when using a gradient-descent method. Second, by identifying these properties, reduced models can be compared and benchmarked to test their efficacy regarding magnetic-field optimization. The first advantage is demonstrated in the present work by the identification of fluxsurface elongation and triangularity as geometric parameters that can be used in HSX to regulate TEM growth rates, while the second will be the subject of future study.
The remainder of this paper is structured as follows. In section 2, a set of volume-averaged metrics are motivated and defined that aim to identify configurations with modified trapped-particle dynamics relative to QHS. Then, an existing helical symmetry metric is discussed and adapted for the database. Section 3 describes how these metrics are used to reduce the full set of configurations to a set of 100 experimental candidate configurations. These configurations are then used to identify flux-surface elongation and triangularity as the geometric parameters of interest regarding TEM stability. The neoclassical tranport of the 100 configurations is analyzed in section 4, and the heat flux is shown to be comparable to the heat flux in QHS. The magnetic-well depth and rotationaltransform profiles are considered in section 5 and used to demonstrate that an increase in elongation and triangularity has other potential benefits. In section 6, TEM growth rates are obtained from gyrokinetic simulations, with a clear reduction in growth rates observed with increasing flux-surface elongation. Finally, in section 7, conclusions and an outlook are provided.

Database metrics
The collisionless TEM is a drift-wave instability that is driven by the resonant interaction of magnetically trapped electrons and poloidally propagating drift waves [21][22][23][24][25]. A necessary condition for this resonance to occur depends on two frequencies, namely the electron diamagnetic frequency ω * e = (T e k α /e)dn e /dψ, and the magnetic electron drift frequency ω de = k ⊥ · v de . Here, v de is the electron drift velocity and k ⊥ = k ψ ∇ψ + k α ∇α is the wave vector perpendicular to a magnetic field-line, with ψ and α defined in the Clebsch-representation of the magnetic field B = ∇ψ × ∇α, so that ψ is the toroidal magnetic flux used to label nested flux surfaces and α labels field lines within a flux surface. Furthermore, T e , n e and e are the electron temperature, density, and electric charge, respectively.
Resonance requires that ω de ω * e > 0, where the overline denotes a bounce-average, so that both the bounce-averaged trapped-electron drift and the electron diamagnetic drift are propagating in the same direction [24,26]. If the electron density profile is monotonically decreasing, then ω * e < 0. Therefore, the resonance condition requires that ω de < 0. The sign of ω de in the low-β, high aspect-ratio limits is determined primarily by the sign of the curvature-drive C d at the mirror points, as described in appendix A. From the local 3D equilibrium model [27,28], C d can be defined as: Here κ n = κ ·n and κ g = κ · (b ×n), where κ =b · ∇b, b = B/B, andn = ∇ψ /|∇ψ|, so that κ n and κ g are the normal and geodesic curvatures, respectively. The term D is defined in: where s is the local magnetic shear, φ is the toroidal angle in a straight field-line coorindate system, and ι-′ = dι-/dψ, where ι-is the rotational transform. Therefore, D describes the integrated variation of the local magnetic-shear relative to the fluxsurface averaged shear along a magnetic field-line. The relationship between C d and ω de is derived for a quasiaxisymmetric equilibrium in [28], with a similar derivation for the quasi-helically symmetric case provided in appendix A. The salient feature of the relation is that, in the paradigm of a collisionless TEM instability that requires ω de < 0, C d < 0 describes a destabilizing curvature-drive in the electron diamagnetic direction while C d > 0 describes a stabilizing curvature drive in the ion diamagnetic direction.
In figure 1, the magnetic field strength and C d are plotted along a magnetic field-line in the QHS configuration. The magnetic field strength and curvature drive are shown in blue and orange, respectively, while λ is the poloidal coordinate in a straight field-line coordinate system. Included in the figure is a plot of just the κ n /|∇ψ| term from the curvature drive, shown as the dotted orange curve. Given that |∇ψ| > 0, the two orange curves show that the sign of C d is determined primarily by κ n . It can also be observed that the regions where κ n < 0 strongly overlap with the magnetic wells. Because the trapped particle density peaks near the bottom of a magnetic well, this overlap results in a considerable free energy source for TEMs, thereby making TEMs the dominate microinstability in HSX [15,16]. Therefore, one possible approach to regulate the TEM stability in HSX is to identify magnetic-field configurations that exhibit a modified distribution of κ n with respect to the magnetic wells. This can be done by considering a set of volume-averaged metrics, where the volume average is denoted as: with √ g the flux-surface Jacobian, V the total enclosed volume, and θ and ζ are the poloidal and toroidal fluxcoordinates, respectively. A volume-averaged magnetic-field curvature metric is then defined as: where B max = B max (ψ) is the maximum magnetic-field strength on each flux-surface. B κ can be interpreted as a volume average of κ n , weighted by the trapping-well depth.
Denoting the value of B κ for the QHS configuration with a * superscript, it is found that B * κ ≈ −5.11 × 10 −3 T/m, where the negative sign suggests a destabilizing curvature drive with respect to the trapping-well weighting.
In addition to the correlation between negative normalcurvature and trapping-well depth, it is also important to consider the weighting of κ n with respect to 1/|∇ψ|, since it is the product of these two terms that appear in the curvature drive. Therefore, a second metric is defined as: and describes a volume average of the normal curvature's contribution to C d . For QHS, it is found that C * n = −7.21 × 10 −3 (T · m 2 ) −1 , indicating a volume-averaged destabilizing contribution to C d .
Neither B κ nor C n account for changes in the neoclassical transport properties of a given configuration. Therefore, a third helical symmetry metric is defined as: where ψ LCFS is the toroidal magnetic-flux at the last closed flux surface (LCFS), B nm are the Boozer modes [29] of the magnetic-field strength on a given flux surface, n and m are the toroidal and poloidal mode numbers, respectively, and the summation is performed over all symmetry-breaking modes. Note that the symmetry modes in HSX are all integer multiples of (n, m) = (4, 1). Therefore, low values of Q indicate good helical symmetry, and, by extension, low neoclassical fluxes, with the quantity for the QHS configuration being Q * ≈ 8.52 × 10 −3 .

HSX coil-current database
In this section, the HSX coil-current database is introduced, and the metrics defined in the previous section will be used to identify configurations that exhibit a modified distribution of normal curvature throughout the plasma volume while also preserving good neoclassical properties. Moreover, these configurations will be shown to exhibit increased flux-surface elongation and triangularity. As depicted in figure 2, the HSX magnetic field is generated by 48 3D-shaped main coils and 48 planar auxiliary coils. Due to the four-fold symmetry in HSX, both sets of the 48 main and auxiliary coils can be subdivided into four sets of twelve, each set spanning a single field period. Moreover, accounting for stellarator symmetry, only six of the twelve main and six of the twelve auxiliary coils are unique. In the standard QHS configuration, all main coils are powered at 100% their typical operating current, while the auxiliary coils are not powered.
To assess the geometric flexibility in the HSX magnetic field for reducing TEM turbulence, a coil-current database has been produced using the 3D magnetohydrodynamic (MHD) equilibrium code VMEC [30], which was run in free-boundary mode to calculate the magnetic-field structure assuming zero equilibrium β and fixed ψ LCFS . Here β = 2µ 0 p/B 2 is the ratio of the thermal pressure to the magnetic field pressure.
Each configuration in the database is produced via coilcurrent modifications to the six unique main and six unique auxiliary coils. Main coil currents are reduced in steps of 1% of the typical operating current, with a maximal reduction of 10%. This maximal modification is based on the expected limitation of what could be experimentally realizable in HSX, while the step size results in a computationally tractable number of equilibria that can be calculated. Moreover, the reduction of the current is only considered for two coils at a time. This is because, at present, all main coils are connected in series, so reducing the current in any one of them will require a shunt resistor to by-pass that particular coil. Restricting the analysis to no more than two main coil-current modifications at a time limits the effort needed in the event such a configuration is pursued in future experimental campaigns. The choice in coil-current permutation, combined with the eleven coil-current levels, defines a total of 1561 unique configurations. In addition to the main coils, all six auxiliary coil currents are permuted across states of (−10, 0, +10)% of the corresponding default main coil current. Hence, an additional 3 6 = 729 auxiliary coil-current configurations are obtained, producing a total of 1561 × 729 = 1 137 969 unique configurations. Discretizing the coil-current space in this way allows for a sampling of the stellarator-symmetric equilibrium space that may be explored in future HSX experiments.
To sort the database, all metrics introduced in section 2 are normalized to the corresponding calculation for the standard QHS configuration. In all figures, the QHS geometry is identified by the star-shaped datum point. In figure 3, C n /|C * n | is plotted with respect to B κ /|B * κ |, and Q/Q * is shown in color. It can be observed that a region with relatively low Q/Q * bisects the distribution of C n /|C * n | and B κ /|B * κ | configurations. It is important to consider that as |C n | becomes large, the κ n curvature drive is expected to increase in the electron diamagnetic direction. Also, as |B κ | becomes small, the regions where κ n < 0 are expected to move away from regions with deep trapping wells. Therefore, following the set of configurations that span the valley of low Q/Q * from the top left of the figure to the bottom right, one expects that configurations will exhibit increased curvature drive in the electron diamagnetic direction, but that curvature drive will be shifted away from the magnetic trapping wells. What the net effect of these trends are on the TEM growth rates is unclear, however a clear region in the metric space from which to sample a subset of configurations for gyrokinetic simulations is presented.
To further investigate this region, a set of 100 configurations were selected. The selected configurations are shown in figure 3 as the square data points. From this set, it is observed that flux-surface elongation and triangularity, denoted as K and δ, respectively, are well correlated with B κ and C n . This can be seen in figure 4, where both K/K * and δ/δ * are plotted as functions of B κ /|B * κ | and C n /|C * n |, where K * = 1.50 and δ * = 0.66 are the corresponding values in the QHS configuration, and the symmetry breaking ratios are shown in color. This demonstrates that flux-surface shapes in HSX can be significantly modified while preserving relatively good helical symmetry. Note that configurations with a large elongation ratio K/K * indicate flux surfaces that are, on average, elongated in the direction perpendicular to the displacement of the helical magnetic axis relative to the center of the torus and compressed in the direction parallel to that displacement. A quantitative definition for K is provided in appendix B. Furthermore, configurations with a large triangularity ratio δ/δ * exhibit an increased absolute value of triangularity in a rotated reference frame, which is described in appendix C. Three exemplary configurations are plotted in figure 5, where poloidal cross-sections are shown of the flux surfaces in a configuration with a low elongation ratio in figure 5(a), QHS in (b), and a highly elongated configuration in (c).
The correlations in figure 4 suggest that the distribution of normal curvature throughout a plasma volume is modified in HSX by elongating the plasma boundary and increasing the triangularity. For simplicity, the 100 configurations will be plotted throughout the remainder of this paper based on their elongation ratio alone. However, because both shaping parameters are strongly correlated, the individual impact of each parameter on an HSX equilibrium cannot be discerned. The impact of the shaping modifications on TEM growth rates are presented in section 5, but first a few considerations regarding the experimental viability of these configurations are presented.

Neoclassical analysis
One important consequence of running VMEC in free boundary mode with a fixed ψ LCFS is that the volume enclosed by each configuration is related to the total current in all external coils. Therefore, if the total current is reduced, the equilibrium volume will increase as the magnetic field strength is reduced. For the set of 100 configurations, the volume increase is approximately linear with elongation, going from 0.32 to 0.41 m 3 , while the the volume-averaged magneticfield decreases from 1.11 to 0.86 T, possibly having a deleterious effect on plasma confinement. Moreover, from an experimental perspective, it is import to verify explicitly that the neoclassical transport in the 100 configurations is comparable to or does not substantially exceed that in QHS, otherwise any reduction in TEM transport may be obfuscated in experiment by an increase in neoclassical transport.
To address these concerns, the SFINCS code [31,32] was used to model the neoclassical heat flux across all 100 configurations. SFINCS solves the 4D linearized drift-kinetic equations, (16) and (27) in [32], with a linearized Fokker-Planck collision operator on an individual flux surface with a specified radial electric field E r . A scan over E r can then be performed to find the field strength at which the electron and ion particle fluxes are equivalent, thereby satisfying the ambipolarity condition. Typically, there are three roots associated with this condition. In HSX the ion root is the one with the weakest radial electric field and is typically observed [17]. In low collisionality, the cross-field neoclassical transport in the ion root is more sensitive to the quality of the helical symmetry than the electron root. For this reason, it is the ion-root heat flux that is considered in the present analysis.
The SFINCS code requires as inputs the magnetic field geometry, which is provided by VMEC, and density and temperature profiles for both the ion and electron species. Using experimental HSX profiles results in highly singular behavior in the particle fluxes due to the proximity of the resonant electric field with the ambipolar electric field [32], a consequence of the low T i /T e ratio. HSX, however, is currently undergoing an upgrade in which a 70 GHz gyrotron will be installed, potentially allowing for higher electron densities and hotter ions. Figure 6 shows the predicted density and temperature profiles for the upgraded HSX. Using these profiles, the electronand ion-root particle fluxes were calculated in SFINCS for the QHS configuration and compared to results from the PENTA code [33,34]. The particle fluxes calculated in the two codes were found to be in agreement within a few percent across the radial domain of the plasma, indicating the increase in the T i /T e ratio in the temperature profiles is sufficient to avoid the singularities in the SFINCS particle fluxes. Note that these simulations are merely meant to assess the quality of the neoclassical confinement of a given configuration compared to QHS and not for validation.
The flux-surface-averaged ion-root electron heat-flux ⟨Q e ⟩ was calculated for all 100 configurations and QHS. Figure 7 shows ⟨Q e ⟩ as a function of elongation at four radial positions in the plasma, with the symmetry-breaking ratio shown in color. It can be observed that the neoclassically predicted heat flux for K/K * ≲ 1.05 is comparable to, or in some cases better than, the heat flux predicted for QHS. For elongation ratios larger than 1.05, a sharp increase in ⟨Q e ⟩ is observed for r/a ⩽ 0.6, which is consistent with the corresponding increase in the symmetry-breaking ratio. This indicates that configurations  with K/K * ≲ 1.05 make good experimental candidates to test what impact plasma elongation has on turbulent transport in HSX.

Considerations on HSX experimental viability
In this section, preliminary steps are taken to examine the experimental viability of the 100 configurations previously identified. These considerations include calculations of the magnetic-well depth and the formation of magnetic islands around low-order rational surfaces.
Though HSX is typically operated at low β, changes in the magnetic-field geometry have been shown to produce significant changes to the ideal-MHD β stability properties as predicted by both ballooning mode theory and the Mercier criterion [35]. Moreover, in [35], it has been shown that increasing the magnetic-well depth leads to an increase in the onset β for ideal-MHD instability. Therefore, the normalized  Figure 8 shows the normalized well depth as a function of flux-surface elongation at various r/a. It can be observed that configurations with increased elongation exhibit deeper magnetic wells, and that a magnetic hill forms around r/a = 0.8 in configurations with K/K * ≲ 0.97. Therefore, at high elongation, one expects more robust MHD stability properties.
Because the elongated axis rotates poloidally as one moves around the torus, increasing the elongation will also increase the rotational transform ι-= ι/2π. This effect can be seen clearly in figure 9, which shows ι-at various r/a as a function of elongation. When ι-crosses a low-order rational value, magnetic islands can form, which may significantly erode plasma confinement [37]. Therefore, it is important to avoid such  rational values. Figure 9 indicates the 4/4 and 8/7 rationals with two horizontal lines, from which ι-can be observed to cross each low-order rational at low and high elongation ratios, respectively.
In HSX, the 4/4 resonance produces the largest magnetic islands of all low-order rationals. This is demonstrated in figure 10(a), where a Poincaré plot is shown for a configuration with K/K * ≈ 0.96 at a toroidal angle of π/4, with the LCFS of the VMEC equilibrium shown as the thick red line, and the interior wall of the vacuum vessel shown as the black line. A large 4/4 island chain is observed just outside the LCFS of the VMEC equilibrium, which significantly distorts the flux surfaces traced out in the Poincaré plot. Because VMEC assumes nested flux surfaces, this distortion is not captured in the VMEC equilibrium. The VMEC rotational transform is shown as a function of the normalized radius in figure 10(b), where ι-is observed to approach the 4/4 resonance at the edge but not cross it. Therefore, when considering configurations with the lowest available elongation, the proximity of ι-to the 4/4 resonance at the plasma edge means many such configurations are experimentally untenable, as they are likely to exhibit large magnetic islands just outside the LCFS that significantly distort the flux surfaces in the bulk.
The 8/7 resonance can be found in some configurations with 1.04 ≲ K/K * ≲ 1.05, as can be observed in figure 9, but this resonance has considerably less impact on flux-surface quality. This can be seen in figure 11(a), where a Poincaré plot is shown for a configuration with K/K * ≈ 1.05, and with the rotational transform provided in figure 11(b). It can be observed that ι-is sufficiently above the 8/7 resonance in the plasma core to avoid the formation of any large-scale magnetic islands. The 12/10 and 16/13 resonances do produce stochastic field lines near the plasma edge, but this is tolerable because the layer of stochasticity appears sufficiently narrow to avoid affecting transport in the core. This demonstrates that configurations with K/K * ≳ 1.05 can be found with good flux surfaces in the plasma core by allowing for a narrow region of stochasticity at the plasma edge.
In summary, there exist multiple experimental candidate configurations with high elongation ratios that exhibit increased magnetic-well depth and good flux surfaces in the plasma core.

Gyrokinetic analysis
Having identified viable candidates to assess anomalous transport experimentally, the next step is to predict which configurations are likely to exhibit strong or weak microinstability and turbulence. To study the effect of plasma shaping on linear TEM growth rates, gyrokinetic simulations have been performed across all 100 configurations with the GENE code [38,39], which solves the coupled gyrokinetic Vlasov-Maxwell system of equations. For computational efficiency, flux-tube simulations are performed here with the necessary geometry information obtained from VMEC and the GIST code [40]. The field lines along which the flux tubes were generated lie on the ψ/ψ LCFS = 0.5(r/a = 0.71) flux surface, and each flux-tube covers four poloidal turns in the parallel domain; the present study focuses on the α = 0 field line, which corresponds to a flux tube centered at the outboard midplane at toroidal ζ = 0. This flux surface was selected because it lies within the region where the TEM transport is thought to dominate over the neoclassical transport in HSX, and at four poloidal turns, gyrokinetic simulations have been shown to be insensitive to the choice in α in HSX [41,42].  Linear initial-value simulations were performed, which provide the real frequency ω and growth rate γ for the most unstable mode at a specified binormal wavenumber k y . Both ω and γ are normalized by the ratio of the ion sound speed c s to the effective minor radius a, while k y is normalized to the inverse of the ion sound gyroradius ρ s . Convergence testing was done for 15 configurations sampled uniformly across elongation ratios at k y = 0.1, 0.4, 0.7, and 1, with all grid resolutions numerically converged to within 5% deviation in γ. The highest required resolutions from the set of 15 configurations were used as the resolutions for all configurations. The grid parameters are N x × N z × N v∥ × N µ = 15 × 512 × 64 × 8, where N x and N z are the number of radial and parallel grid points, respectively, and N v∥ and N µ are the number of parallel and perpendicular velocity-space grid points, respectively. The velocity-space box sizes used were L v∥ = 3 and L µ = 9. Physical input parameters were specified as β = 0, ν ei = 0, T i /T e = 0.2, a/L ne = a/L ni = a/L Te = 2 and a/L Ti = 0, where ν ei is the collision frequency and L χ = −(χ/a)/(dχ/dr), with χ the zeroth-order density or temperature of either electron or ion species. These parameters were selected to best approximate the plasma parameters in HSX experimental discharges at this radius [3,17].
All simulations were performed with kinetic ions and electrons, and the most unstable mode is identified as a collisionless TEM [15]. Figure 12(a) shows the growth rate γ of the fastest-growing mode as a function of elongation at four different k y values. Note that a significant reduction in the growth rate can be observed when comparing the highly elongated configurations to those with low elongation. Moreover, a reduction is also observed relative to QHS. Throughout the elongation range, this stabilization can be observed to scale strongly with increasing elongation, except for moderately elongated configurations at k y ⩾ 0.7. Considering that previous studies of nonlinear heat-flux spectra in HSX typically show that the bulk of the transport occurs at long wavelength (k y < 1) [15], this result is expected to lead to a measurable impact on the microinstability induced drive for turbulent transport, as growth rates at relevant scales can be reduced by elongating the plasma. Moreover, the region of K/K * ≈ 1.05, which features low neoclassical heat fluxes and good flux surfaces in the plasma core, may be advantageous from the point of reduced turbulent transport since a clear reduction in the linear growth rates is predicted for all k y < 1. Figure 12(b) shows ω as a function of elongation across four values of k y , demonstrating that all modes propagate in the electron diamagnetic direction (ω < 0). The discontinuities in ω across elongation ratios is a consequence of regime transitions resulting from a substantial number of subdominant modes. Therefore, subtle changes in the magnetic field geometry will change which individual mode dominates, compare [15,41].
Comparing figure 12 with figure 4, the reduction in growth rate with increasing elongation and triangularity means the growth rate is also reduced with decreasing |B κ | and increasing |C n |. This suggests that the shift of negative κ n away from the magnetic trapping-wells, resulting in lower values of |B κ |, has a more significant stabilizing effect than the volumeaveraged increase in the κ n curvature drive in the electron diamagnetic direction described by the increase in |C n |. Therefore, in HSX, TEM stability appears to depend more on the the drifts of the deeply trapped particles than the volume averaged drifts. While this result justifies the use of these metrics in the perturbative optimization approach, it is unclear how they would perform in a gradient-descent based optimization. Furthermore, in a general stellarator equilibrium, these metrics should be modified to include the additional terms in the curvature drive equation (1), since C d will not, in general, be dominated by the κ n /|∇ψ| component.

Conclusions
Optimizing magnetic fields for reduced turbulent transport is at the forefront of modern stellarator design. In this study, a novel perturbative optimization approach around a fixed operating point has been applied to HSX. This has led to the identification of a set of magnetic-field configurations with good neoclassical properties that can be used to experimentally test the optimization of HSX with respect to reduced TEM transport. To that end, a database of more than 10 6 HSX stellarator magnetic equilibria has been generated by modifying individual main and auxiliary coil currents. A set of volumeaveraged metrics was introduced and used to identify configurations in the database that exhibit substantial changes in curvature drive of their precessional drifts. Using these metrics, along with a helical symmetry metric, a set of 100 configurations was identified that spans a wide range of elongation and triangularity values while keeping the deviations from helical symmetry low. The neoclassical heat flux was then calculated for all 100 configurations with SFINCS, which demonstrated that magnetic configurations with increased elongation and triangularity exist that do not significantly degrade the neoclassical transport properties relative to the standard QHS configuration. Therefore, these shaping parameters have been identified as macroscopic features in HSX magnetic-field geometries that can be used to modify the ω de curvature drive while preserving good neoclassical confinement.
The 100 configurations were also shown to have greater well-depth with increasing shape parameters, suggesting that the critical β below which a configuration retains ideal-MHD stability is increased at higher shaping values. Additionally, it was demonstrated that the rotational transform in configurations with low shaping values is close enough to the 4/4 resonance to produce a large region with a stochastic magnetic field, while configurations with high shaping can be found that avoid the 8/7 resonance and produce good flux surfaces in the plasma core. This demonstrates that configurations with high elongation and triangularity ratios make good experimental candidates for future HSX discharges.
Finally, using the GENE code, linear gyrokinetic simulations were performed to calculate the growth rate of the most unstable eigenmode across the set of 100 configurations. Results show a significant reduction in growth rate at wavenumbers k y < 1 with increasing elongation and triangularity. This reduction is hypothesized to result from the shift of regions with negative curvature drive away from the magnetic trapping-wells, thereby reducing the destabilizing contribution of the trapped-electron drifts to the drift wave. Note that in nonlinear simulations, the heat-flux spectrum tends to peak at k y < 1, so the reduction in growth rates of k y ⩽ 1 studied here may indeed translate to a corresponding reduction in nonlinear transport. At present, nonlinear simulations of select configurations are ongoing, and will be reported in a separate publication, along with additional analysis regarding the observed linear stabilization.
This work demonstrates the efficacy of a coil-currentmodification approach towards the optimization of a 3D magnetic field around a fixed operating point, such as existing experiments. The plasma behavior has been shown to be sensitive to subtle changes in the 3D geometry in both neoclassical and gyrokinetic models. During the process of down-selection based on various metrics, many configurations remain as good experimental candidates to validate the approach taken here on future HSX discharges. Such configurations can be used to investigate how trapped-electron drifts affect TEM-driven turbulent heat losses, paving the way for a more robust optimization of future experiments with respect to TEM turbulence. Importantly, this perturbative optimization approach is not unique to HSX, and could be deployed on any existing device, providing a wealth of new experiments in turbulence optimization.
for √ g = (∇ψ · ∇λ × ∇φ) −1 . Substituting equation (A3) into equation (A2) results in: Defining v 2 ∥ = v 2 − 2µB/m e , where µ is the magnetic moment, then differentiating equation (A4) with respect to ψ and K results in: and where ι-′ = ∂ι-/∂ψ. Following [28], one derives: Here, µ 0 is the permeability of free space, p is the plasma pressure, B λ and B λ are the contra-and covariant poloidal components of the magnetic field, and h is a flux-surface geometry term defined in equation (20) of [27]. In a large aspect ratio limit, terms that scale with h and B λ B λ /B 2 are negligible, while in the low-β limit, ∂p/∂ψ can be neglected. Dropping those terms and substituting equation (A7) into equation (A5) it is found that: Equations (A8) and (A6) provide a set of integrable equations to calculate ω de in accordance with equation (A1). Importantly, when integrating along a field line dη = (N − ι-)dφ and dλ = ι-dφ, while the limits of integration are the φ (or λ) value at which a particle is mirrored, making ω de a function of pitch angle. Moreover, because v ∥ = 0 at a particle's mirror point, the integration in equation (A8) is weighted heavily towards the value of C d at those points, where C d is defined in equation (1). Therefore, it is primarily the sign of C d near the bounce-points that determines the sign of ω de .

Appendix B. Flux-surface elongation
Flux-surface elongation is defined in relation to the local magnetic axis. The magnetic axis is described in a cylindrical coordinate system whose origin is located at the center of the torus. In this reference frame, r ma = r marma , where r ma is the distance of the magnetic axis from the origin andr ma is the unit vector pointing towards the magnetic axis from the origin, both of which are functions of the azimuthal angle ϕ. Furthermore,ê ϕ is defined as the azimuthal unit vector, and z ma =r ma ×ê ϕ . Flux-surface elongation is then defined as an expansion of a flux surface in theẑ ma direction and a compression of a flux surface in ther ma direction. This results in a definition of elongation that is a function of both ϕ and ψ. This quantity is then averaged across all ϕ and ψ as: Here, the flux-surface shape is given by r = r (ψ, ϕ, θ). Following this definition, large values of K indicate configurations that are, on average, elongated in the direction perpendicular to the displacement of the helical axis relative to the center of the torus and compressed in the direction parallel to the magnetic axis.

Appendix C. Helical triangularity
In [43], upper and lower triangularity are defined for an axisymmetric equilibrium in a cylindrical coordinate system {R, ϕ, Z}. Labeling the maximum and minimum cylindrical coordinates along a flux surface with the corresponding subscripts, upper triangularity is defined as: where R geo = (R max + R min )/2.
In HSX, the direction of triangularity rotates in phase with the field period N. Therefore, defining a rotated coordinate system {R ′ , ϕ, Z ′ }, where: triangularity can be defined as in equation (C1), but in the rotated reference frame. Importantly, the triangularity switches in HSX from positive to negative δ as the flux surface goes from ϕ = 0 to ϕ = π/4. Because of this azimuthal oscillation of δ, the average δ along ϕ does not change significantly as elongation increases. Therefore, it is the azimuthal average of the absolute value of δ along the LCFS in the rotated reference frame that is presented in figures 4(c) and (d) of section 3.