On the effect of flux-surface shaping on trapped-electron modes in quasi-helically symmetric stellarators

Using a novel optimization procedure it has been shown that the Helically Symmetric eXperiment (HSX) stellarator can be optimized for reduced trapped-electron-mode (TEM) instability [M.J.~Gerard et al., \textit{Nucl.~Fusion} \textbf{63} (2023) 056004]. Presently, with a set of 563 experimental candidate configurations, gyrokinetic simulations are performed to investigate the efficacy of available energy $E_\mathrm{A}$, quasi-helical symmetry, and flux-surface shaping parameters as metrics for TEM stabilization. It is found that lower values of $E_\mathrm{A}$ correlate with reduced growth rates, but only when separate flux-surface shaping regimes are considered. Moreover, configurations with improved quasi-helical symmetry demonstrate a similar reduction in growth rates and less scatter compared to $E_\mathrm{A}$. Regarding flux-surface shaping, a set of helical shaping parameters is introduced that show increased elongation is strongly correlated with reduced TEM growth rates, however, only when the quasi-helical symmetry is preserved. Using a newly derived velocity-space-averaged TEM resonance operator, these trends are analyzed to provide insights into the physical mechanism of the observed stabilization. For elongation, stabilization is attributed to geometric effects that reduce the destabilizing particle drifts across the magnetic field. Regarding quasi-helical symmetry, the TEM resonance in the maximally resonant trapping well is shown to increase as the quasi-helical symmetry is broken, and breaking quasi-helical symmetry increases the prevalence of highly resonant trapping wells. While these results demonstrate the limitations of using any single metric as a linear TEM proxy, it is shown that quasi-helical symmetry and plasma elongation are highly effective metrics for reducing TEM growth rates in helical equilibria.


I. INTRODUCTION
Reaching burning plasma conditions in magnetic confinement fusion (MCF) devices requires high densities, high temperatures, and long energy-confinement times.However, most fusion devices show a strong degradation of the confinement time during high-temperature and high-density operation.This reduction can be attributed primarily to the cross-field transport arising from turbulence at the gyroradius scale, which is driven by microinstabilities occurring at comparable scales [1].
One technique for reducing turbulent transport is to modify the magnetic field geometry [2][3][4][5][6][7][8][9].In the stellarator approach to fusion, this can be achieved with the coils used to generate the magnetic field.The first stellarator in the world to demonstrate the efficacy of magnetic-field shaping for reducing neoclassical transport is the Helically Symmetric eXperiment (HSX) [10,11].HSX is a four-fold symmetric device that can generate a plasma with approximately 1.2 m major radius, 0.12 m minor radius, and plasma volume around 0.44 m 3 .It was designed to be able to generate different magnetic topologies by varying the current in its external coils.Notably, the device can generate configurations with both quasihelical symmetry and broken quasi-helical symmetry, though the standard configuration, referred to as the QHS configuration, is quasi-helically symmetric [12,13].Due to the successful neoclassical optimization of HSX, the dominant source of transport is attributed to turbulence.With densities around 10 18 m −3 and electron-cyclotron heating, the electrons in HSX are much hotter than the ions.Moreover, peaked density profiles are typically observed [14].Therefore, the ion-temperature-gradient (ITG) mode tends not to contribute significantly to the cross-field transport, though ITG modes in HSX geometries have been studied computationally [15,16].Alternatively, simulation work has shown that the trappedelectron mode (TEM) is destabilized in HSX, and likely drives the anomalous transport in the device [17,18].
Recently, it has been demonstrated computationally that the collisionless TEM growth rates can be reduced in HSX by elongating the plasma while preserving the quasi-helical symmetry [19].Those results are based on a set of 100 configurations selected from a database of over 1 million unique magnetohydrodynamic (MHD) HSX-based equilibria.This database was produced by numerically modifying coil currents in the 48 auxiliary and 48 main coils of the experiment.This flexibility in HSX provides a useful domain for testing the efficacy of different TEM optimization metrics while simultaneously investigating the relationship between the magnetic field geometry and the governing TEM resonance.Therefore, in the present paper, the number of configurations investigated in linear gyrokinetic simulation is increased to include configurations with broken quasi-helical symmetry.Moreover, the simulation results are compared against several TEM metrics, and the physical mechanism governing the observed trends is discussed.This process, while consistent with results from Ref. [19], reveals new insights into the optimization of quasi-helically symmetric equilibria for reduced TEM growth rates.
The collisionless TEM is destabilized by a resonance between the precessional drift of the trapped electrons and the electrostatic drift wave.This has been reported extensively for tokamak geometries [20][21][22] and more recently has been shown to apply to stellarators [23][24][25].The stability condition for TEMs is well characterized by the equation for the energy transfer rate between the trapped electrons and the electrostatic drift wave [4,[23][24][25].Near marginal stability (γ → 0 + ), this equation can be expressed as which was first derived in Ref. [23].Throughout this paper, an overbar is used to denote a bounce average, which is defined in Eq. (A3).Then, e is the elementary charge, T e0 the background electron temperature, ℓ the magnetic-field-line arc length, B the magneticfield strength, ϕ the bounce-averaged electrostatic potential, f e0 the zeroth-order electron distribution-function, v the particle velocity, and δ the Dirac-delta distribution-function.
There are also three frequencies in the formula, of which ω is the electrostatic driftwave frequency, ω de is the electron bounce-averaged-drift frequency, and ω T * e the electrondiamagnetic-drift frequency.Note that the e subscript will be dropped for the remainder of this paper.
The diamagnetic-drift frequency is defined as where n 0 denotes the background electron density, E = mv 2 /2 is the kinetic energy with electron mass m and velocity magnitude v, ψ is the toroidal magnetic flux, and α is the Clebsh angle, which locally defines the magnetic field as B = ∇ψ × ∇α.Then, k α is the wavenumber in the ∇α direction.The bounce-averaged-drift frequency is defined as where v d is the electron-drift velocity and k ψ is the wavenumber in the ∇ψ direction.The two components of the bounce-averaged-drift frequency are further identified as the bounceaveraged precessional-drift frequency ω α = k α v d • ∇α and the bounce-averaged radial-drift frequency ω ψ = k ψ v d • ∇ψ.In terms of a particle's second adiabatic invariant J = mv ∥ dℓ, with v ∥ its parallel velocity, the precessional-and radial-drift frequencies can be expressed as where the subscripts denote variables being held constant, with µ = mv 2 ⊥ /2B the magnetic moment and v ⊥ the magnitude of the particle's perpendicular velocity, relative to the magnetic field.
When P < 0, Eq. ( 1) describes a destabilizing transfer of energy from the trapped electrons into the drift wave via resonant interaction.Therefore, the destabilization of a TEM requires that ω(ω − ω T * ), or equivalently ω d (ω d − ω T * ), be negative in regions along a magnetic field line where the trapped-particle fraction is high and the bounce-averaged electrostatic eigenfuction is peaked.The other terms in the expression then provide a fluxsurface and velocity-space weighting of this condition.Note that in so-called omnigenous magnetic fields, where the bounce-averaged radial drifts of trapped particles goes to zero (ω ψ = 0), then ω d = ω α .Therefore, to produce an unstable TEM, an omnigenous magnetic field requires that the diamagnetic drift and the bounce-averaged precessional drift propagate in the same direction.To be sure, no magnetic field is perfectly omnigenous; however, quasiomnigenous configurations with ω ψ /ω α ≪ 1 do exist, in which case the destabilization of TEMs is due to unfavorable precessional drifts in the electron-diamagnetic direction.Since the QHS configuration is quasi-helically symmetric, and therefore quasi-omnigenous, this is also the case for the collisionless TEM in the QHS configuration.
Calculation of the energy transfer described in Eq. ( 1) requires knowledge of the driftwave frequency and the bounce-averaged electrostatic potential, both of which typically require explicit calculation as done in the gyrokinetic code Gene [26,27].For stellarator optimization, the more computationally efficient calculation of the available energy E A could be used as a proxy for TEM growth rates.E A is a TEM metric that defines an upper bound on the amount of energy that is available in the trapped-electron population to drive fluctuations [28][29][30][31].As described in Ref. [30], E A can be calculated in flux-tube geometry as where ∆ψ and ∆α are the length scales over which energy is available and Note that the flux-tube E A can be normalized as where ρ * = ρ s /a with ρ s the ion gyroradius at the ion sound speed, a is the effective-minorradius at the last closed flux surface (LCFS), and is the thermal energy of a plasma in a flux tube to leading order in an expansion in the directions perpendicular to the magnetic field [31].
It is important to note the significance of ∆ψ and ∆α in the calculation of E A .In a gyrokinetic simulation, one may resolve fluctuation dynamics in a volumetrically minimized domain by making a flux-tube approximation [32].Importantly, in this local approximation, increasing the simulation domain beyond this minimally resolved domain does not increase fluctuation energies.Therefore, for E A to relate to these local fluctuations, ∆ψ and ∆α must correspond to the length scale of those fluctuations.Following arguments presented in Ref. [31], these free parameters in the available energy are defined as ∆ψ = ρ 2 * ψ edge and ∆α = ρ * / √ s, where s = ψ/ψ edge and ψ edge is the toroidal magnetic-flux at the LCFS, so that s labels the flux surface being considered.
The concept of available energy has been used to argue for a possible connection between reduced TEM turbulence and reduced neoclassical radial particle transport in the collisionless regime.This conjecture is motivated by the observation that a Maxwellian distribution can only satisfy the minimal-E A state if the magnetic field is omnigenous [29,30].This is observed in Eq. ( 6) as a monotonically increasing dependence of E A on ω 2 ψ .Since quasiomnigeneity is achieved in HSX through quasi-helical symmetry, TEM growth rates will be compared against the E A and the quasi-helical symmetry of various HSX equilibria.
Quasi-helical symmetry can be quantified as where N and M are the number of toroidal and poloidal modes in a straight field-line coordinate system and B n,m are the mode amplitudes of the magnetic field strength on a flux surface [33].Due to the n = 4m symmetry in the QHS magnetic field, Q quantifies the extent to which the symmetry is broken in the magnetic field's so-called Boozer spectrum.
Recent analytic work has also shown the resiliency of so-called maximum-J configurations towards TEM stability [23,24], where maximum-J means that ∂J /∂ψ < 0 throughout an equilibrium.Notably, this property is locally satisfied for deeply trapped particles in quasiisodynamic configurations like Wendelstein 7-X (W7-X) and is achieved over a broader set of trapped-electron pitch-angles in new configurations presented in Refs.[34,35].The effect of this maximum-J condition is to reduce the TEM resonance by confining the electron distribution function to a domain in which trapped electrons dominantly precess in the iondiamagnetic direction [23].Despite the absence of TEMs, gyrokinetic simulations in W7-X geometries still show finite growth rates for modes propagating the electron-diamagnetic direction [36].These modes have been identified as the universal instability (UI), which can be destabilized by the passing-or trapped-electron population, and these modes dominantly appear in situations when the TEM has been sufficiently stabilized [37].
In this paper, a set of 563 unique HSX equilibria that span a broad range of coil-current perturbations is analyzed for their TEM growth rates.This allows for the testing of existing TEM metrics and analytic theories.Section II provides details on the set of equilibria that have been selected, with a description of how E A and quasi-helical symmetry can be used in conjunction with a set of helically defined flux-surface shaping parameters.In Sec.III, results from linear gyrokinetic simulations are shown to be consistent with findings in Ref. [19], where the lowest growth rates were observed in highly elongated configurations.However, it is observed here that elongation must not come at the expense of quasi-helical symmetry, which, when broken, leads to an increase in growth rates.Moreover, no universal correlation between TEM growth rates and E A is observed.Instead, both E A and quasi-helical symmetry are shown to scale with growth rates when considered in particular shaping regimes.A velocity-space-averaged TEM resonance operator is defined in Sec.IV and compared against the growth rates.From this analysis, it is found that the reduction in growth rates is consistent with the energy transfer rate defined in Eq. ( 1), but that the stabilization does not result in a transition from a dominant TEM to UI.By analyzing trapped-electron drifts in helically linked magnetic trapping wells, the growth-rate reduction with an increase in elongation is shown to be the result of an increase in the rotational transform and a modification in the magnetic-field-strength spectrum that preserves quasi-helical symmetry.Moreover, the dependence of growth rates on quasi-helical symmetry is shown to be caused by an increase in the prevalence of highly resonant trapping wells.These results show that elongation can be used along with either E A or a quasi-helical-symmetry metric as a promising combination of metrics for the optimization of helical equilibria for reduced TEM growth rates.Moreover, this further demonstrates the efficacy of the optimization technique introduced in Ref. [19] of exploring the configuration space surrounding the operating point of an existing experiment.

II. EQUILIBRIUM SELECTION A. Selection Metrics
In this section, a set of new shaping parameters is defined for helical equilibria, which will be used to characterize flux-surface geometries throughout the HSX coil-current database.
The shaping parameters are based on those defined in Ref. [38] for an axisymmetric system.
Therein, definitions are provided for flux-surface elongation κ, triangularity δ, and squareness ζ, all of which are based on the extrema along a flux surface at a single toroidal angle in a cylindrical coordinate system.The primary adaptation made to extend these definitions to a non-axisymmetric system is to define a rotated reference frame from which the flux-surface shapes can be calculated.The rotated reference frame is defined as where {R, Z} and {R ′ , Z ′ } are the radial and vertical coordinates in the cylindrical and rotated reference frames, respectively, and N fp and φ are the number of field periods and toroidal angle, respectively.This rotated reference frame is selected so that, relative to the magnetic axis, the radial-like direction êR ′ = ∇R ′ /|∇R ′ | always points towards the low-field side of a toroidal cross-section, making the shaping parameters analogous to those of the axisymmetric system.Figure 1 demonstrates this rotated reference frame in three toroidal cross-sections of the QHS equilibrium, with the magnetic field strength shown in color and the rotated basis vectors shown after translation to the magnetic axis.
It should be noted that in Ref. [38], a distinction is made between upper and lower values of elongation and triangularity.For squareness, a further distinction between upper/lower and inner/outer is made.For our purposes, all such distinctions are calculated and then averaged together at each toroidal angle.Moreover, due to the lack of axisymmetry, each shaping parameter is a function of toroidal angle.Thus, an integrated toroidal average for each shaping parameter is used as the final parameter under consideration.Though helpful in characterizing flux-surface shapes, it should be realized that because of this toroidal averaging, these shaping parameters do not uniquely describe a flux-surface geometry.

B. Selection Process
An initial set of 600 configurations, of which 563 will be retained for gyrokinetic analysis, is selected from the database of > 10 6 MHD equilibria, with the database described in Ref. [19].These configurations are selected using the helical shaping and symmetry metrics.To ensure all metrics are calculated along flux surfaces with the same real-space position, different fluxsurface labels are selected for consideration in each configuration.With flux-surface label s = ψ/ψ edge = (r/a) 2 , where r is the effective-minor-radius of the flux surface, the standard flux-surface is taken to be the QHS surface at s q = 0.5.For each i th configuration it is then required that r i = r q .This leads to the salient relation s i = (a q /a i ) 2 s q , where s i labels each flux surface to be considered.
The three shaping parameters, the quasi-helical symmetry metric, and normalized fluxtube E A are all calculated along the s i flux surfaces and normalized to their corresponding quantity on the s q = 0.5 flux surface, with the latter denoted by a star superscript.Since both triangularity and squareness can take on either positive or negative values, the absolute value of the corresponding metric in the QHS configuration is used in the normalization.
In this way, the sign of the normalized metrics is determined by the configuration being considered.For the E A calculation, the α = 0 field line is selected and followed for four poloidal turns with 512 grid points along the direction of the magnetic field vector.The GIST code [39] is used to calculate these flux-tube geometries.Due to the inclusion of ω T * in Eq. ( 6), E A depends on the electron temperature and density gradients across the flux tube.These gradients are taken to be the typical experimentally observed quantities in QHS [10,14], defined as a/L T e = a/L ne = 2 for the electron temperature and density gradients, respectively.Note, L χ is the characteristic gradient scale length defined as L χ = −χ/(dχ/dr).
The procedure for down-sampling the configurations then goes as follows.Using the three helical-shaping parameters, along with the quasi-helical symmetry metric, the database of configurations is partitioned into a set of selection bins that cover all metric values calculated from the equilibria.From each bin, a random set of configurations is selected in proportion to the number of configurations in that bin.This prevents over-or under-sampling any particular region of the configuration space.Importantly, configurations with rotational transforms that cross the 4/4 or 8/7 resonance are excluded, as such configurations can form large islands inside the plasma confinement region which contravene the VMEC solutions (see Sec. 5 in Ref. [19]).

C. Equilibrium Reconstruction
To calculate an ideal MHD equilibrium, VMEC requires an initial guess for the LCFS geometry and the enclosed toroidal magnetic flux ψ edge .When constructing the original database, ψ edge was held fixed, and the LCFS of the QHS configuration was used as the initial guess.Then VMEC was run in free-boundary mode, meaning the initial guess for the LCFS was allowed to change based on how the coil-current information modified the vacuum magnetic field, with the LCFS expanding or contracting in volume to achieve the specified ψ edge .Moreover, solutions were provided using eight toroidal and eight poloidal mode numbers, allowing for computational expediency when generating the database.
To ensure a robust geometric comparison in the proceeding gyrokinetic analysis, higher spectral resolutions are desired.However, VMEC often struggles to find free-boundary solutions in HSX geometries with poloidal or toroidal mode numbers greater than twelve.
Alternatively, running VMEC in fixed-boundary mode allows for much higher mode resolution, but requires that a reliable LCFS be supplied as input since this boundary will be fixed as the solution is identified.Therefore, the LCFS of each down-sampled configuration is identified using a magnetic field-line-following (FLF) code in the corresponding vacuum-field configuration.This LCFS is then used to re-generate the VMEC equilibrium in fixed-boundary mode with 24 and 18 toroidal and poloidal mode numbers, respectively.The use of vacuum FLF is justified as the plasma pressure p in HSX is low enough to warrant the β = 0 assumption, where β = 2µ 0 p/B 2 with µ 0 the permeability of free space, making β the ratio of thermal pressure to the magnetic-field pressure.
Initially, the same ψ edge is used in all high-resolution VMEC calculations, and the fluxsurface-averaged magnetic field strength ⟨B⟩ i is calculated on each s i surface.Then, VMEC is run a second time in fixed-boundary mode with the same mode resolution but with ψ edge,i = (⟨B⟩ q /⟨B⟩ i )ψ edge,q .This results in a set of s i surfaces with the same flux-surfaceaveraged field strength.As a final check, the high-resolution s i surfaces are compared against the corresponding vacuum-field surface produced by FLF.This is to ensure the VMEC surfaces are not distorted by the presence of magnetic islands near those surfaces.
In total, 563 configurations have been retained of the 600 configurations initially sampled.Regarding the shaping parameters, these results demonstrate the intuitive notion that elongation is determined almost exclusively by low-order modes, while triangularity and squareness are related to progressively higher-order modes.Less intuitive, however, is that quasi-helical symmetry, like elongation, is determined by low-order modes, while E A depends strongly on moderate-to high-order modes.This latter result is an important observation for the use of E A as an optimization metric, showing that it requires high-resolution equilibria for an accurate calculation.This does not, however, cast doubt on the process for downsampling the equilibria, since the high-resolution quantities are observed to span a range of metric values comparable to that found in the low-resolution equilibria.Therefore, the approach of investigating the equilibrium space around an existing configuration using the rapid generation of magnetic equilibria remains a valid optimization procedure.Lastly, the dependence of quasi-helical symmetry and E A on low-and high-order modes, respectively, may provide some benefit regarding TEM optimization.For example, one could target quasi-helical symmetry by modifying the low-order modes and then perform a fine-tuning of the resultant configuration by targeting higher-order modes to enhance TEM stability.
Additional support for this notion is presented towards the end of Sec.IV.

III. LINEAR GYROKINETIC ANALYSIS
A set of linear gyrokinetic simulations across all 563 configurations and QHS are performed using the Gene code [26,27].The binormal wavenumbers investigated include k y ρ s = 0.1, 0.4, 0.7, and 1.This choice of wavenumbers is motivated by nonlinear simulations in QHS, where it has been observed that the heat flux spectrum peaks at k y < 1 [18].
The flux-tube geometries used in these simulations are selected from the α = 0 field line and generated using the GIST code [39].Convergence testing was done for 10 configurations  Relative to the largest observed growth rates, the reduction factor is 29 and 3.1 for k y ρ s = 0.1 and 0.4, respectively.These results are consistent with the gyrokinetic results from Ref. [19] regarding the stabilizing effect of plasma elongation.However, these new results further demonstrate that the reduction in γ with increasing κ is limited by the ability to increase κ while preserving the quasi-helical symmetry.This is exemplified by the increase in growth rates at high elongation and increasingly broken quasi-helical symmetry.
In Fig. 4, a set of outlier configurations is present at all wavenumbers and is distinguished from the bulk configurations with square markers.Notably, in Fig. 4 E A can be observed, exemplified by the wide range of growth rates at any given value of E A .
However, if elongation is considered, one observes two separate branches of configurations, with each branch delineated by elongations higher or lower than that observed in the QHS flux surface.These separate branches exhibit a less scattered trend, showing reduced growth rates for lower E A .This suggests that, as a linear metric, E A may be limited in its ability to identify a globally optimized configuration, but that within a particular shaping regime, it may be used as a local optimization metric.For example, if one were to begin optimization with an initial configuration appearing in the top-right quadrant of Fig. 6(b), and then smoothly deform the flux surface to reduce the E A , one would likely arrive at a configuration that has a lower growth rate.However, there could still be a considerable disparity between that growth rate and another growth rate from a configuration with a very similar E A .
Furthermore, the set of outlier configurations with large squareness is identified by their exceptionally low E A and their uncharacteristically large growth rates at some k y .The reason for the anomalously large growth rates in relation to the E A is unknown.However, it is possible that these configurations constitute a third shaping regime, distinct from the high-and low-elongation regimes, in which separate scaling relations could exist between γ and E A .
To relate the E A of these configurations to their quasi-helical symmetry, Fig. 7 shows the growth rates as a function of the symmetry-breaking ratio, with elongation shown in color.and 7, similar trends can be identified; however, the scaling between γ and Q in Fig. 7 appears less scattered.A certain degree of consistency between E A and quasi-helical symmetry is to be expected, given that minimizing ω ψ is one method by which E A can be reduced, as seen in Eq. ( 6).What is more surprising is that in the comparison of Figs. 6   and 7, the quasi-helical symmetry appears to be a better predictor of the relative growth rate than E A .Of course, the set of outlier configurations is still present, appearing in Fig. 7 as the set of configurations with the lowest elongations and Q/Q * ≤ 5.

IV. TEM RESONANCE ANALYSIS A. Deriving the velocity-space-averaged TEM resonance operator
To better understand these results, and to compare them to existing analytic theory, a velocity-space-averaged resonance operator is defined and calculated from the linear gyrokinetic data.This resonance operator is motivated by Eq. ( 1), restated here for convenience Note that ω(ω − ω T * ) is the salient quantity regarding the influence of trapped-electron drifts on the TEM stability of a particular flux tube.Therefore, the velocity-space-averaged resonance operator is defined.By evaluating this quantity, one may characterize the contribution of the trappedelectron drifts to the destabilizing transfer of energy described by P < 0. The evaluation will be carried out by first making a series of substitutions to express the integral in terms of particle energy and pitch angle, and then normalizing the equation variables to be consistent with Gene.The evaluation of the Dirac-delta distribution function will then be made possible with the drift-wave frequency ω calculated in Gene, which will result in the reduction of Eq. ( 12) to a field-line integral.
The velocity differential can be transformed to a particle energy and pitch angle differential as described in Ref. [41] as where ε = E/T 0 is the normalized particle energy, θ b is the poloidal angle where a trapped particle is bounced -therefore relating it to the pitch angle -and the phase-space-density factor is with B = B min /B(θ b ), B min the minimum magnetic-field strength over the poloidal domain, and B ′ = dB/dθ.Importantly, p(θ b ) is proportional to the trapped-particle density, and peaks when B = B min .Then, a Maxwellian is used for the zeroth-order distribution function, To express Eq. ( 12) in a manner consistent with Gene, first consider that the drift-wave frequency and wavenumber normalizations are where the 0 subscript denotes the variable is evaluated along the field line at the center of the flux tube, q 0 is the safety factor, and ω, kx , and ky are the corresponding dimensionless values in Gene [26,32].Then, the diamagnetic-drift frequency is with ωT * = ky [(a/L ne ) + (a/L T e )(z − 3/2)], and the bounce-averaged-drift frequencies are 18) Note that ωy = ε ky vdy and ωx = ε kx vdx , with normalized drift velocities vdy and vdx defined in Eq. ( A21)-(A22).A general discussion on how these quantities are calculated from a GIST output is presented in Appendix A. Finally, the full bounce-average-drift frequency is Thus defined, the root of the Dirac-delta distribution in Eq. ( 12) can be identified as with ξ = T 0 /(2eq 0 ψ edge ).Then, substituting Eq. ( 13)-( 15) into Eq.( 12) and applying the Gene normalizations, Eq. ( 12) is reduced to the field-line integral with H(ε ′ ) the Heaviside function and Importantly, radial wavenumber resonances in Gene are related to the ballooning representation of the electrostatic eigenfunction [32].This relation is shown graphically in Fig. 8, where the square of the bounce-averaged-eigenfunction amplitude in the QHS configuration at ky = 0.1 is shown as a function of the poloidal-bounce angle in ballooning space.This choice of wavenumber is selected because in the derivation of Eq. ( 1) it is assumed that γ → 0 + , which is best fulfilled at low k y .Therefore, in the subsequent analysis, the focus will be on the ky = 0. From Fig. 8, a significant eigenfunction amplitude at finite-radial wavenumber can be observed, with the peak found at |j| = 2.This means the radial drifts cannot be neglected in Eq. ( 22) a priori, especially when considering configurations with broken quasi-helical symmetry.However, as will be shown, the TEM resonance is predominantly determined by the precessional, rather than radial, drifts.Therefore, to understand how the precessional drifts influence the TEM resonance, it is helpful to define the velocity-space-averaged bounceaverage precessional drift as Here, Eqs. ( 13), (15), and ( 18) are used to reduce the velocity-space average to a field-line integral covering the poloidal domain.

B. TEM and UI cross-phase and resonance analysis
It is additionally illuminating to consider the cross-phase between the electrostatic potential and density fluctuations, with the trapped-and passing-electron populations considered separately.The reason for this separation is that TEM and UI cross phases are typically in proximity to π/2 in the trapped-and passing-electron populations, respectively.Therefore, by comparing the cross-phase data between modes, one may identify mode branches that have a TEM, UI, or mixed cross-phase.Cross phases are commonly calculated for each population of electrons as a histogram throughout all of phase space.Therefore, a cross-phase distribution must be considered.The process by which these histograms are reduced to a scalar is described in Appendix B; however, the salient quantity presented here is M ϕ×n , which is shifted towards −1 when the cross phase has a strong UI signature, +1 when indicating a TEM, and ≈ 0 when the mode is a mixed UI/TEM.
Figure 9 shows the maximal growth rates at ky = 0.1 for all 564 configurations as a function of ⟨ω(ω − ω T * )⟩, with M ϕ×n shown in color.In panel (a), all radial wavenumbers are included in the evaluation of Eq. ( 22), while in panel (b), only the kx = 0 mode is considered, thereby excluding information regarding the radial drifts.It can be observed that the reduction in growth rates, in both panels, is correlated with a shift of ⟨ω(ω − ω T * )⟩ towards zero, demonstrating that the stabilization is consistent with the analytic theory of collisionless TEMs [23,24].
A set of outliers are presented with triangle markers in Fig. 9(b), and identified by their large negative values of ⟨ω(ω − ω T * )⟩ and relatively low-growth rates.Due to the strong UI signature of these modes, with M ϕ×n < −0.86, one does not expect the velocity-spaceaveraged TEM resonance operator to accurately capture their stabilization.However, a similar trend is observed with these outliers as is observed in the bulk configurations, namely that the growth rates are reduced as the resonance shifts towards zero.Interestingly, the outlying nature of these configurations is minimized when all finite radial wavenumbers are included, as seen in Fig. 9(a).Moreover, the cross-phase metric for these configurations is shifted towards −1 as the TEM resonance is reduced, with the lowest outlier growth rate exhibiting M ϕ×n = −0.94.This shows these modes are increasingly UI dominant as the TEM resonance is reduced.Alternatively, no clear transition from a TEM to a UI cross-phase is observed in the bulk configurations.To be sure, if one disregards the outlier configurations, the strongest UI signatures are found near these zero crossings, but so too are the strongest TEM signatures.A likely explanation for this is that the UI is stronger in the outlier configurations, making the transition to UI dominance more accessible with a reduction in the TEM resonance.Interestingly, these outliers do not exhibit any unique properties regarding their helical shaping parameters, meaning they are decidedly not the same set of outliers identified in Sec.III.This means the resonance operator accurately captures the resonance properties of the large-negative-squareness configurations, which are identified by their square markers.This further highlights the inability of E A to characterize respectively.The square markers denote the set of large-negative-squareness configurations and the triangle markers are a set of configurations with a strong UI cross phase.It is observed that the growth-rate reduction is consistent with analytic theory and that the resonance is determined primarily by the precessional, rather than radial, drifts.
TEM growth rates in various shaping regimes.
In Ref. [36], the configuration that exhibited the clearest transition from a TEM to UI cross phase was the W7-X configuration most congruous with the maximum-J condition.
Therefore, to reconcile the observed reduction in the TEM resonance operator without this corresponding transition, it is informative to consider how well these configurations satisfy ∂J /∂ψ < 0. From Eq. ( 4) it is known that the sign of ω α is determined exclusively from ∂J /∂ψ.Therefore, Eq. ( 24) can be used as a proxy for maximum-J -ness.
Figure 10 shows the ky = 0.1 growth rates as a function of {ω α }, with the cross-phase signature M ϕ×n shown in color.It is of significant note that no clear trend appears in the comparison between the growth rates and velocity-space-averaged precessional drift or the cross-phase signature.Moreover, the sign of {ω α } indicates that ∂J /∂ψ > 0 is likely for all configurations, meaning the maximum-J condition is never observed.The implications here are twofold.First, the reduction in the TEM resonance results from the detailed interaction between the trapped-particle fraction, electrostatic eigenfunction, and the bounce-averaged precessional drifts rather than a general trend towards drift reversal along the entire flux tube.Second, the drifts that persist in the destabilizing electron-diamagnetic direction are sufficient to prevent the dominance of the UI over the TEM.

C. Impact of elongation on TEM resonance
Motivated by these observations, the impact of elongation on the velocity-space-averaged precessional drifts will be investigated in individual trapping wells instead of along the entire magnetic field line.This is done in Fig. 11, where three individual trapping wells are identified for each flux tube, with panel (a) showing a characteristic partitioning of these wells in the QHS configuration, labeling trapping wells as i = 0, 1, or 2. These three wells are selected because they span one complete poloidal rotation.In panels (b), (c), and (d), the velocity-space-averaged precessional drift for the electrons populating each trapping well is shown as a function of elongation, with the symmetry-breaking ratios for each configuration shown in color.Recall that {ω α } < 0 is destabilizing for TEMs while {ω α } > 0 is stabilizing.
From Fig. 11(b), (c), and (d), the first thing to be observed is that the dependence of {ω α } on elongation is qualitatively different in each trapping well, an observation that is consistent with the lack of any trend in growth rates with respect to {ω α } when the latter is calculated along the entire flux-tube.The starkest comparison is between {ω α } 0 and {ω α } 1 , shown in panels (b) and (c), respectively, with the subscript denoting the i th trapping well.It can be observed that these wells exhibit opposite trends, where {ω α } 0 and {ω α } 1 are shifted towards stabilizing and destabilizing values, respectively, as elongation increases.
However, if one considers the symmetry-breaking ratio Q/Q * , then {ω α } 0 can be seen to shift toward stabilizing values at high elongation while preserving good quasi-helical symmetry in some configurations.Alternatively, the most stabilizing values of {ω α } 1 only appear in low-elongation configurations with severely broken quasi-helical symmetry.Similarly, in well 2, where the dependence on elongation is more complicated, the fact remains that more configurations with good quasi-helical symmetry, and {ω α } 2 shifted towards zero, are found at high elongation.Therefore, in HSX, if one is restricted to considering only configurations with relatively good quasi-helical symmetry, then increasing elongation will reduce the destabilizing drifts in some, but not all, magnetic trapping wells, thereby reducing linear growth rates without eliminating the TEMs.Similarly, the increase in the destabilizing drifts in well 1 with increasing elongation may contribute to the increase in growth rates at high elongation observed in Fig. 4, though this effect is difficult to distinguish from the effect of breaking the quasi-helical symmetry, which occurs simultaneously.Therefore, it remains unclear whether growth rates would continue to decrease with increasing elongation if the quasi-helical symmetry could be preserved to arbitrarily high elongation.
To investigate the drift reversal in well 0 with increasing elongation, Fig. 12  From Fig. 12 it can be observed that the trapping well in (b) and (d) is more extended in the poloidal domain.This can be attributed to an increase in the configuration's rotational transform, which results from an increase in its helical elongation (see Fig. 9 in Ref. [19]).
One significant effect of this extension can be observed in panels (c) and (d), where the maximal pitch angle for which a particle's bounce point (v ∥ = 0) occurs in a region with stabilizing curvature (L 2 > 0) is higher in the high elongation configuration.This means the high elongation configuration has fewer trapped particles that spend their entire trajectory in regions of destabilizing curvature.
Another effect of elongation is that it excites higher-order Boozer harmonics.This is observed in panels (a) and (b) of Fig. 12, where the more elongated configuration has a larger magnetic-field ripple.Importantly, the symmetry-breaking ratio in these two configurations is nearly identical, with Q/Q * = 1.772 and 1.784 for the left and right configurations, respectively.This is possible because a significant fraction of the harmonics excited in the more elongated configuration are attributed to higher-order symmetry modes such as (n, m) = (8, 2) and (12,3).Analogous to a square well, the inclusion of these higherorder harmonics results in a flattening of the magnetic field strength towards the bottom of the trapping well.This flattening then allows for higher parallel velocities through the region in which L 2 is negative, meaning the destabilizing curvature's contribution to the bounce-averaged precessional drift is minimized.
To demonstrate this effect of elongation throughout the down-sampled database, a new metric intended to capture the presence of high-order symmetry modes is introduced.Similar to Eq. ( 10), which quantifies the symmetry breaking in the magnetic-field-strength spectrum, this new quantity is defined as , where J is the number of symmetry-preserving modes.To exclude the dominant symmetry mode (n, m) = (4, 1), the summation begins with j = 2.For the QHS configuration transform.More interesting, however, is that if one considers configurations with good quasi-helical symmetry, indicated by low Q/Q * , then increasing elongation results in an increase in the amplitudes of the higher-order symmetry modes, with amplitudes increasing by 50% in the most elongated quasi-helically symmetric configurations.
These shaping effects are subtle, but their influence on the particle drifts is profound, as exemplified by the near-reversal in {ω α } shown at the top of Fig. 12(a) and (b) and more broadly in Fig. 11(b).This then supports the conjecture made in Sec.II, where it was stated that TEM stability optimization could target the higher-order modes in the magnetic spectrum while preserving the overall quasi-helical symmetry.This is an important observation that can be leveraged in future optimization.It is therefore concluded that a quasi-helical symmetry configuration will be less unstable relative to a broken-symmetry configuration because of the minimal variation in the underlying resonance across the individual trapping wells.Alternatively, in a broken-symmetry configuration, the inclusion of non-symmetric trapping wells will increase the likelihood that more strongly resonant wells will be accessible, leading to greater instability.This is in agreement with the observations made previously based on E A , namely that a reduction in E A is observed with improved quasi-helical symmetry.Since the observed dependence of E A on quasi-helical symmetry is best explained by the dependence of E A on the bounce-averaged radial particle drifts ω ψ , these results provide some evidence for the more general claim that quasi-omnigeneity might limit TEM growth rates by reducing available energy in the system.

V. CONCLUSIONS
Using an HSX coil-current database [19], a set of 563 HSX equilibria has been compared against the QHS configuration to identify trends between macroscopic shaping parameters and TEM stability.It has been found that the lowest TEM growth rates occur in helically elongated geometries that preserve the quasi-helical symmetry of the QHS configuration.
Moreover, the reduction in growth rates with improved quasi-helical symmetry is shown to be consistent with predictions from available-energy calculations.However, both E A and quasi-helical symmetry fail to predict the reduction with respect to elongation.
Alternatively, the velocity-space-averaged resonance operator has been shown to accurately capture the reduction in TEM growth rates.However, this resonance calculation requires knowledge of the drift-wave frequency and eigenmode structure, precluding this quantity as an easy-to-evaluate TEM metric.The quantity does demonstrate utility in revealing the physical mechanisms that contribute to the reduction in growth rates.For example, it has been found that the reduction in growth rates with increasing elongation occurs by reducing the TEM resonance without supporting a maximum-J flux surface.Moreover, breaking the quasi-helical symmetry is shown to broaden the distribution of trapping-well resonances to include a greater number of highly destabilizing trapping wells.
Regarding general trends towards optimization, it has been found that the elongated configurations of HSX produce lower TEM growth rates because of the extension of the trapping wells in the poloidal domain as well as the inclusion of high-order symmetry harmonics in their Boozer spectrum.The latter suggests that if a configuration has been optimized for quasi-helical symmetry, then further optimization to increase the amplitude of high-order symmetry harmonics may be beneficial for reducing TEM growth rates.However, it is not clear, at present, how such an optimization would be performed since the sign of a particular symmetry mode will determine whether that mode steepens or flattens a magnetic well, and Eq. ( 25) does not preserve such information.Therefore, it may be easier to target the excitation of these modes by simply targeting increased plasma elongation, with a strong penalty for configurations that break the quasi-helical symmetry.Notably, consistency between the quasi-helical symmetry metric and E A , as seen when comparing Figs. 6 and 7, suggests that a similar approach could be to target both increased elongation and reduced The next step is to perform nonlinear simulations to verify key findings, which will be reported in a future publication.Regarding the experimental validation of these observations, plans are underway to measure density and temperature fluctuations and to calculate heat fluxes from power balance analysis in HSX geometries identified from the coil-current database.
The equilibria considered throughout this paper are all vacuum field configurations.As such, p′ = 0 and L 2 is identified, from Eq. (A21), as the curvature drive for the precessional drift.Note, −L 1 and L 2 are provided as a function of z in the GIST output.

FIG. 1 .
FIG.1.Three toroidal cross-sections are shown for the QHS equilibrium, with the helically rotated basis vectors shown after translation to the magnetic axis.It is within this reference frame that the helical shaping parameters are calculated.Importantly, êR ′ is constructed to point towards the low-field side of a given cross-section.

Figure 2
Figure 2 shows the quasi-helical symmetry metric in (a) and the flux-tube E A in (b), each as a function of elongation, with triangularity shown in color.The set of 600 down-sampled configurations is shown as filled squares in the foreground, while all other configurations in the database are shown as hollow squares in the background.Otherwise, the layering of configurations in each set is random.The black cross-hairs denote the metric values of the QHS flux surface in each figure.Note that thousands of configurations exist with comparable, or better, quasi-helical symmetry than QHS, and available energies are observed that exhibit a 10% reduction from QHS.In Fig. 2(a) the down-sampled configurations are seen to probe the full extent of the metric values, except for gaps found in regions with more extreme values of elongation as Q is increased.These gaps result from the exclusion of configurations with rotational transforms near the aforementioned resonances, meaning these excluded configurations are unlikely to make good experimental candidates.Moreover, in Fig. 2(b) the down-sampled configurations are observed to span the range of ÊA / Ê * A from 0.89 to 1.45 in both κ and δ shaping regimes.This demonstrates that the set of 600 configurations provides a representative sampling of the metric values produced in the For reference, the high-resolution QHS metrics evaluate to κ * = 1.471, δ * = 0.0523, ζ * = −0.0901,Q * = 0.0089, and Ê * A = 0.0270.With the metric quantities recalculated for the new equilibrium solutions, it is informative to compare them to their original values.Such a comparison is shown in Fig. 3, where metric quantities are compared in the fixed-and free-boundary mode solutions.Note that the metric quantities calculated from the high-resolution fixed-boundary solutions are plotted along the vertical axis while the metrics calculated from the low-resolution free-boundary solutions are plotted along the horizontal axis.The dashed black line indicates perfect agreement between metric quantities calculated from both solutions.The metrics being considered are the quasi-helical symmetry metric (a), flux-tube E A (b), elongation (c), and triangularity (d).It can be observed that both the quasi-helical symmetry and elongation values are nearly identical when calculated from either equilibrium solution.Alternatively, both the E A and triangularity show considerably more scatter between their high-and low-resolution values.Note that squareness, not shown in the figure, is found to have no discernible trend between the high-and low-resolution values since the low-resolution equilibria fail to capture large values of squareness.

FIG. 3 .
FIG. 3. Metric quantities are compared in the fixed-and free-boundary solutions for the ideal MHD equilibria, where (N = 24, M = 18) and (N = 8, M = 8) resolutions are used, respectively.The metric quantities under consideration are the quasi-helical symmetry (a), flux-tube E A (b), elongation (c), and triangularity (d).Note that the dashed black line indicates the hypothetical agreement between metric quantities calculated from both equilibrium solutions.

Figure 4
Figure4shows the TEM growth rates γ of the most unstable mode as a function of plasma elongation, with γ normalized to the ion sound speed c s divided by effective minor radius, and quasi-helical-symmetry ratios shown in color.Panels (a), (b), (c), and (d) correspond to wavenumbers k y ρ s = 0.1, 0.4, 0.7, and 1, respectively.The location of the QHS configuration is identified in each panel by the black dashed cross-hairs.It can be observed that the lowest growth rates, across all wavenumbers, tend to occur in configurations with high elongation and low symmetry-breaking amplitudes.Relative to QHS, the minimal growth rates are reduced by a factor of 4.8 and 1.7 at high elongation for k y ρ s = 0.1 and 0.4, respectively.

FIG. 4 .FIG. 5 .
FIG. 4. Growth rates are shown as a function of plasma elongation at binormal wavenumbers 0.1, 0.4, 0.7, and 1 in panels (a), (b), (c), and (d), respectively, with quasi-helical symmetry ratios shown in color.The black cross-hairs indicate the QHS configuration in each panel and the square markers are a set of large-negative-squareness configurations.The lowest growth rates are typically found in configurations with high elongation and good quasi-helical symmetry.

FIG. 6 .Figs. 6
FIG. 6.Growth rates are shown as a function of the flux-tube E A at binormal wavenumbers 0.1, 0.4, 0.7, and 1 in panels (a), (b), (c), and (d), respectively, with elongation ratios shown in color.The black cross-hairs indicate the QHS configuration in each panel and the square markers are a set of large-negative-squareness configurations.A robust scaling between growth rates and E A is only observed when the different elongation branches are taken into consideration separately.

FIG. 7 .
FIG. 7. Growth rates are shown as a function of symmetry breaking ratios at binormal wavenumbers 0.1, 0.4, 0.7, and 1 in panels (a), (b), (c), and (d), respectively, with elongation ratios shown in color.The black cross-hairs indicate the QHS configuration in each panel and the square markers are a set of large-negative-squareness configurations.Similar trends as observed in Fig. 6 are observed here, though the trends with respect to their quasi-helical symmetry are less scattered.

4 FIG. 8 .
FIG.8.Bounce-averaged-electrostatic eigenfunction of the QHS configuration at ky = 0.1 as a function of poloidal-bounce angle in ballooning space.Each 32π poloidal extension of the magnetic field line about the central j = 0 section is delineated with a vertical dashed line.Sections with a finite radial wavenumber are identified by j ̸ = 0.

FIG. 9 .
FIG. 9. k y ρ s = 0.1 growth rates are plotted as a function of the velocity-space-averaged TEM resonance operator, with the resonance calculated for all radial wavenumbers in panel (a) and only for k x = 0 in panel (b).The color scale distinguishes TEMs from UI going from blue to red,

FIG. 10 .
FIG.10.k y ρ s = 0.1 growth rates are plotted as a function of {ω α }, while the color scale distinguishes TEMs from UI going from blue to red, respectively.This shows that no configuration is likely to satisfy ∂J /∂ψ < 0, which suggests the configurations may not be sufficiently TEMstabilized to observe the transition from a dominant TEM to UI.
shows the magnetic-field strength, field-line curvature L 2 , and normalized-parallel-particle velocity v∥ = v ∥ /v as a function of poloidal angle in a set of characteristic trapping wells for two different configurations.In panels (a) and (b) the magnetic field strength is shown in blue and the curvature in orange, the latter of which determines the precessional drifts as described in Eq. (A21).Then, in panels (c) and (d), normalized-parallel-velocity contours are shown, with their dependence on pitch angle Λ indicated by the various gray-scale bands.In both the top and bottom panels, the poloidal domain of the destabilizing curvature is highlighted in orange.Note that the configuration shown in panels (a) and (c) has κ/κ * = 0.989 and that shown in (b) and (d) has κ/κ * = 1.023.

FIG. 11 .
FIG. 11.In panel (a) the magnetic field strength in the QHS configuration is plotted as a function of straight poloidal and toroidal field-line coordinates, and the characteristic α = 0 field line is shown in red.The field line is partitioned into three separate trapping wells indexed by i.In panels (b)-(d) the velocity-space-averaged bounce-averaged precessional drift frequency is shown for trapping wells 0, 1, and 2 as a function of elongation, with the symmetry-breaking ratio shown in color.The black cross-hairs denote the QHS configuration.This shows how the particle drifts are affected by elongation across these three trapping wells.

FIG. 12 .
FIG. 12.The magnetic trapping well crossing the outboard mid-plane is shown for two different configurations.In (a) and (b) the magnetic field strength and curvature drive are shown in blue and orange, respectively, each plotted as a function poloidal angle.In (c) and (d), the normalizedparallel-velocity contours are shown for the corresponding trapping wells above.The parallel velocity color bands indicate the pitching angle of the trapped electrons.In both top/bottom sets the region of destabilizing curvature is highlighted in orange.The trapping well on the right has higher elongation, demonstrating the flattening of the well that leads to the increase in {ω α }, as shown atop panels (a) and (b).

FIG. 13 .
FIG.13.The higher-order symmetry metric is shown as a function of elongation, with the symmetry-breaking ratio shown in color.This shows that increasing elongation while preserving quasi-helical symmetry leads to an excitation of higher-order symmetry harmonics in the magneticfield-strength spectrum.The black cross-hairs indicate the QHS configuration and the horizontal gaps in the data are explained in the text.

FIG. 14 .
FIG.14.The maximum and minimum values of the velocity-space-averaged resonance operator, determined by calculating the resonance within all partitioned trapping wells along each flux tube, are shown as a function of the symmetry-breaking ratio.This shows that the resonance within the minimally resonant well does not change as the quasi-helical symmetry is broken, while the resonance of the maximally resonant well is increased.Therefore, a broken-symmetry configuration has an increased likelihood of producing a highly resonant trapping well, leading to greater instability over the entire flux tube.

FIG. 15 .
FIG. 15.Four histograms are shown for the TEM resonance operator across all individual magnetic trapping wells.The data compiled in each histogram stems from configurations that span four different ranges of quasi-helical-symmetry ratios.This shows that breaking the quasi-helical symmetry broadens the distribution of trapping-well resonances to include a greater number of highly destabilizing trapping wells.