Gyrokinetic studies of trapped electron mode turbulence in the Helically Symmetric eXperiment stellarator

The development of quasi-symmetry¹ has provided stellarators with a viable path to magnetic confinement fusion as an alternative to the tokamak. The Helically Symmetric eXperiment (HSX) is a four period modular stellarator that possesses quasi-helical symmetry and has been shown to reduce neoclassical transport levels to values similar to those observed in tokamaks.² The residual transport observed in HSX discharges is turbulent. In this work, we employ gyrokinetic simulation techniques to study the properties of trapped electron mode turbulence in HSX geometry.

The density-gradient-driven, collisionless trapped electron mode (TEM)³ has been shown to be a viable candidate to explain the electron thermal fluxes observed in many tokamak experiments with peaked density profiles^4,5 and is detrimental to achieving confinement in fusion plasmas. Fundamentally, the density-gradient-driven TEM exists for drift-wave frequencies $ω<ωb$⁠, where ω_b is the electron bounce frequency on a flux surface. It is destabilized by density perturbations of the trapped particle population in the presence of bounce-averaged bad magnetic curvature. For tokamaks, this criterion is typically satisfied in a single poloidally localized region near the outboard midplane.

Over the past decade, there has been significant work and understanding gained in the gyrokinetic simulation of stellarator plasmas as it relates to the ion temperature gradient (ITG) and electron temperature gradient (ETG) modes^6–13 and subsequent turbulent transport optimization.^14–19 However, only recently have computational capabilities reached the point where TEMs in stellarators can be simulated effectively. As such, analysis of TEM-driven turbulence in stellarators is in the early stages, with most prior work relating to the W7-X stellarator.^12,20–23

The TEM destabilization condition is easily met in HSX, as will be shown in Sec. II, due to the overlap of regions of bad curvature and magnetic wells, a property that generally occurs for quasi-symmetric stellarator configurations.^8,24 One would then expect TEMs to be linearly unstable in HSX for some density gradient, which has been shown in the linear gyrokinetic growth rate calculations performed in Ref. 25. The present work represents the first examination of the nonlinear state of TEM-driven turbulence in HSX, and together with the work published in Ref. 23 the first examination of TEM turbulence in neoclassical-transport-optimized stellarators.

Previous work on the density-gradient-driven TEM in tokamaks has shown that zonal flows are the primary nonlinear saturation mechanism.^26–28 Zonal flows can provide a saturation mechanism through flow shear^29,30 or by providing a route through which an instability can couple and transfer energy to damped modes.^30–36 The nonlinear simulations of HSX presented in this paper show the development of zonal flows consistent with the standard picture in tokamaks, and we show that the zonal flows are the nonlinear saturation mechanism.

The paper is structured as follows: In Sec. II, the details of HSX equilibria are introduced and important numerical considerations in applying gyrokinetics to HSX are discussed. In Sec. III, we present linear gyrokinetic simulation results, which show the existence of multiple, distinct ion-scale modes propagating in the electron direction. The nonlinear simulation results in Sec. IV show that the excitation of an ion-propagating direction coherent structure not predicted by linear simulations is associated with the development of zonal flows. Results are summarized and discussed in Sec. V.

We investigate trapped particle drift-waves through the gyrokinetic framework.^37,38 The Gene code³⁹ is used to solve the coupled gyrokinetic Vlasov-Maxwell system of equations in both linear and nonlinear simulations. Simulations are performed in flux tube geometry,^40,41 where, by definition, the background equilibrium pressure and its gradients are constant. The GIST code⁴² is used to calculate the metric elements of the 3-D magnetic field configuration of HSX in a flux tube, where the metric elements are now solely a function of the parallel coordinate along the magnetic field. In a periodic modular stellarator such as HSX, the flux tube is only computed for one field period and quasi-periodic boundary conditions are applied in the parallel direction.

In a tokamak, one flux tube is identical to all other flux tubes on the same magnetic surface, and thus a single flux tube can be used to represent a magnetic surface. This is not the case in stellarators,⁴³ where different flux tubes can experience different curvatures and trapping regions along the field line, with resulting changes to drift-wave behavior. The two representative flux tubes used in this work at the half-radius in the toroidal flux coordinate $Ψ/Ψ0≃0.5(r/a≃0.71)$ are identified by the shape of the plasma cross-section at zero poloidal angle. The normalized normal curvature and normalized magnetic field strength of the flux tubes are given in Fig. 1, corresponding to a bean shaped cross-section, and Fig. 2, corresponding to a triangular cross-section. The flux tubes will hereafter be identified as HSX-b and HSX-t, respectively, and are the up-down symmetry planes one-half field period away from each other. As seen in Figs. 1 and 2, a feature of both flux tubes is the correlation of the regions of bad (negative) curvature and magnetic field minima (trapped particle regions), which is similar to tokamak geometry. However, unlike a tokamak, the majority of particles trapped in these wells are not toroidally, but rather helically trapped.

The geometry of HSX is such that the averaged magnetic field shear is very small across the plasma radius. In flux tube simulations, periodicity at the parallel boundary requires⁴¹ 

where the background magnetic shear is defined as $ŝ=r0ιdιdx|r0$⁠, where $d/dx$ denotes the radial derivative of the rotational transform profile ι and r₀ is the minor radius of the flux surface under investigation. L_x designates the radial box size, $N$ is a positive integer, and $kyminρs$ is the minimum binormal mode number simulated. If $ŝ$ becomes small, the box size in the radial direction can become quite large, in contrast to many tokamak applications with $ŝ∼O(1)$⁠. At the half-toroidal flux surface, the background magnetic shear is $ŝ=−0.045$ and decreases as $r/a→0$⁠, which leads to typical box sizes $Lx≳200ρs$⁠, where ρ_s is the ion sound gyroradius. As the numerical resolution in the radial direction is set such that $O(ρs)$ features are captured, the large radial box size provides a significant computational constraint for flux-tube simulations.

Two different TEM parameter sets will be investigated. As HSX operates with primarily ECRH heating, the ions have a flat temperature profile for a majority of the minor radius. As such, there will be no ion temperature gradient in the following simulations. The first parameter set will be called the “standard” TEM and was selected to capture the essential physics of the density-gradient-driven TEM when both a density and electron temperature gradient are present but the density gradient is larger (⁠$ηe=dlnTe/dlnne<1$⁠). It uses the normalized density gradient $a/Ln=2$ and normalized electron temperature gradient $a/LTe=1$⁠, where a is the averaged minor radius, $Lξ=1ξdξdx$ is the inverse length scale for any scalar ξ, and $Te0/Ti0=1$ with $mi/me=1836$⁠, where $Ts0, ms$ is the background temperature and mass of species $s$⁠. The standard TEM parameter set will be applied to simulations in both flux tubes.

The second parameter set was selected to examine the sensitivity of the TEM to only the driving density gradient. This parameter set uses four different density gradient values, $a/Ln=1,2,3,4$ and no electron temperature gradient. This parameter set will only be used in simulations done for HSX-b. The normalized electron pressure is $β=8πne0Te0/B02=0.05$%, where $ne0$ is the background electron density.

The linear growth rates, frequencies, and mode structures presented in this section have been calculated by two different methods. Initial value simulations are used to obtain only the fastest growing mode at each $kyρs$⁠. Eigenvalue simulations calculate both the dominant and subdominant modes at each $kyρs$⁠. The three-dimensional nature of the magnetic geometry makes for a rich eigenmode landscape characterized by many subdominant, but growing modes for each $kyρs$⁠.

The growth rates have been checked for convergence with resolution, by performing a series of simulations with increasingly fine resolution until the growth rates show no change. Two factors make these simulations substantially more expensive than tokamak simulations: (1) the magnetic field along the field line has more structure, and as such, stellarator simulations may require between 3 and 10 times as much resolution in the parallel direction;^9,23 and (2) the low $ŝ$ allows linear modes at low $kyρs$ to become very extended along the field line and require large numbers of radial connections to resolve, corresponding to increasingly fine radial grids in real space. For linear simulations, we used a numerical grid size of $Nz×Nv∥×Nμ=64×36×8$ with $17≤Nx≤97$⁠, where μ is the magnetic moment. The parallel hyper-diffusion coefficient, see Ref. 44, was set to $ϵz=4$⁠.

First, we examine the linear behavior for both HSX-b and HSX-t. The growth rates and real frequencies for the standard TEM case in both flux tubes are presented in Fig. 3. The magnitude of the maximum growth rates are different (⁠$γ=0.643$ in HSX-b vs. $γ=0.443$ in HSX-t), and the peak occurs at slightly different $kyρs$ (⁠$kyρs=2.1$ in HSX-b vs. $kyρs=1.6$ in HSX-t). The difference in growth rate magnitudes can be explained by examining Figs. 1 and 2. With z in Figs. 1 and 2 as the parallel coordinate, one can see in Fig. 1 (HSX-b) a region centered around z = 0 corresponding to a magnetic field minimum and unfavorable curvature. The situation is reversed in HSX-t, where modes localized near z = 0 see favorable curvature. The observed faster growth at high-$kyρs$ in HSX-b is more strongly ballooning, i.e. its structure is mainly confined to k_x = 0, whereas its HSX-t counterpart peaks at a finite $|kx|=±4,6π$⁠.

The peak growth rates in HSX are at larger $kyρs$ compared with density-gradient-driven TEM in tokamaks,^26,28 where $kyρs(γmax)≲0.7$⁠. They are comparable to the $kyρs$ of peak growth rates in W7-X,²³ where $kyρs(γmax)≈1.6$⁠. Applying a simple quasilinear estimate for the heat flux, $QkyQL∝γk/(kyρs)2$ (assuming k_x = 0) suggests that scales at $kyρs<1$⁠, where both HSX-b and HSX-t have similar linear growth rates, will likely dominate transport. Thus, despite there being a significant difference in the maximum linear growth rates between the two flux tubes, the nonlinear behavior of the flux tubes may be similar.

The real frequencies show electron propagation direction modes for $kyρs≤2$⁠. This feature of the linear modes is robust against variations in β, where the low-$kyρs$ linear modes maintain electron propagation direction. No kinetic ballooning modes are expected, as $β=0.05$% is well below marginal ideal MHD ballooning stability point for HSX.⁴⁵ The real frequencies in Fig. 3 show discontinuities in $kyρs$ for $kyρs<1$⁠, which is indicative of the coexistence of different modes. These modes will be classified in Sec. III C. At the high-$kyρs$ end of the spectrum, the real frequencies for both flux tubes cross the $ωr=0$ boundary shortly after the peak growth rate, switching from electron to ion propagation direction. This type of mode has been observed in linear TEM simulations for a variety of configurations from stellarators²³ to tokamaks^26,28,46 and the reversed-field pinch,⁴⁷ and is referred to as the “ubiquitous” mode.⁴⁸ Despite the high growth rates in these modes, the high $kyρs$ modes contribute a relatively small amount to the overall transport, as shown in Sec. IV, due to the quasilinear $k⊥−2$ scaling.

To isolate the role of the density gradient, we set the electron temperature gradient to zero and only vary the density gradient. The growth rates and real frequencies for the density gradient scan in HSX-b are presented in Fig. 4. Without an electron temperature gradient, the linear growth rates and frequencies for the density gradient scan are similar to that of the standard TEM, providing evidence that the standard TEM is a density-gradient driven TEM with little sensitivity to $a/LTe$⁠, and is a general property of HSX plasmas, consistent with previous work.^8,24,25 The high $kyρs$ ubiquitous mode transitions to a different mode at even higher $kyρs$ as $a/Ln$ is increased. The behavior of the different modes at low $kyρs$ becomes clearer, as the different mode boundaries are delineated by discontinuities in the real frequencies. Examination of the potential mode structures, shown for the standard TEM parameters in Sec. III C, reveal that the dominant mode at each $kyρs$ for each $a/Ln$ are the same as the modes at each $kyρs$ for the standard TEM. At $kyρs=0.1$⁠, the growth rate for each $a/Ln$ is small, which was also observed in the standard TEM simulations, indicating it may not be necessary to resolve lower $kyρs$ in nonlinear simulations.

The eigenmodes in Figs. 3 and 4 are distinguished under the assumption that the real frequencies of a particular mode branch should be a continuous function in $kyρs$⁠. One can distinguish four different dominant linear modes in HSX-b but only two different dominant linear modes in HSX-t. The potential structure in the ballooning representation for the different mode branches of HSX-b and HSX-t are shown in Figs. 5 and 6, respectively, as a function of the ballooning angle $θ=z+2πp$⁠, where p is an integer.⁴⁹ For HSX-b, there is significant variation in the structure between modes. At $kyρs=0.1$⁠, mode B1 in Fig. 5 has a strongly extended envelope in θ with maxima symmetrically located at roughly $θ∼±10π$ in the ballooning angle. On top of the extended envelope, there are spikes in the mode amplitude within each radial connection (every $2π$ in θ), which correspond to the localized bad curvature at z = 0 in Fig. 1. These spikes are well-resolved. The broad envelope and narrow spikes indicate two disparate scales in the modes at low $kyρs$⁠. The combination of low $kyρs$ and small $ŝ$ leads to a slab-like mode that extends along the field line, where the larger scale is related to the magnetic shear $ŝ$ and the small scale is related to the localization via bad curvature. This result is consistent with what occurs in a tokamak, where lower $kyρs$ modes have broader ballooning structure,⁵⁰ particularly in the $ŝ≪1$ limit.^49,51

For $0.2≤kyρs≤0.6$⁠, corresponding to mode B2 of Fig. 5, the mode still has an extended structure, but the envelope peaks more narrowly at $θ≈±20π$⁠, which corresponds to a finite-k_x mode. Mode B3, dominant for $0.7≤kyρs≤0.9$⁠, is considerably more localized but possesses tearing parity, $Φ$ has odd parity, unlike the other modes. The mode peaks in the first radial connection, $±2π$⁠, indicating it is a finite-k_x mode. For $kyρs>0.9$⁠, the ubiquitous mode is now highly localized to θ = 0. This mode is responsible for the maximum growth rate and the structure of the mode does not change as the real frequency changes from electron to ion propagation direction.

Some of the modes in HSX-t show a markedly different character than the modes at equivalent $kyρs$ for HSX-b. The mode T1a at $kyρs=0.1$ has two distinct, extended lobes that peak at $θ=±9π$⁠. Unlike the mode B1 of Fig. 5, the amplitude in the central $−π<θ<π$ range, corresponding to k_x = 0, is small. The T2 modes of Fig. 6, dominant for $kyρs≥0.8$ have the same character and parity, with peak amplitudes on either side at $θ=±5π$⁠. The T2 modes are identified with the ubiquitous mode for HSX-t and are distinct from the B3 mode due to lack of tearing parity. There does not appear to be an equivalent mode with tearing parity in HSX-t.

The sensitivity of these modes to the driving gradients was examined by independently applying a 10% variation to the density and temperature gradients of the standard TEM parameter set. For each of the different modes in Figs. 3 and 4, the growth rates have a positive correlation with increasing density and electron temperature gradient, but consistently show a stronger dependence on the density gradient. Increases in β generally result in moderately decreasing growth rates, hence partial stabilization, for all modes except B1.

With all of these considerations in mind, we have classified the B modes in HSX-b and the T modes in HSX-t according to Table I. We will call the B1 mode TEM1 and mode B2 TEM2. The B3 modes are called “tearing-parity” TEM (TTEM) and behave much like the ubiquitous mode B4. Finally, the B4 modes are ubiquitous modes (UTEM). HSX-t shows only two distinct modes, with the T1 mode being a more traditional TEM and the T2 mode UTEM for that flux tube.

To examine the behavior of the dominant eigenmodes in Figs. 3 and 4, we computed the first five most unstable linear eigenmodes for the standard TEM at every $kyρs$⁠. The results are shown in Fig. 7 with the five eigenvalues denoted as EV#1 to EV#5. An aspect of the eigenvalue decomposition that is pervasive in HSX simulations and markedly different from tokamak results is the clustering of linear eigenmodes. The first four most subdominant modes in the range $kyρs≤0.6$ all have growth rates that are only slightly smaller than that of the dominant mode and very similar, real frequencies. As such, the precise ordering of the four subdominant modes was not investigated for $kyρs≤0.6$⁠.

Fig. 7 shows that the TTEM is the dominant mode for $0.7≤kyρs≤0.9$ (the red square) but becomes subdominant to the ubiquitous mode for $kyρs≥1$⁠. At no point in the range $0.7≤kyρs≤2.1$ does the TTEM cease to exist and all the subdominant modes at high $kyρs$ exhibit the same real frequency behavior as the ubiquitous mode. One feature of Fig. 7 is the pairing of two subdominant, unstable modes which are almost identical in both growth rate and frequency, but not mode structure. For example, in Fig. 7, for $1≤kyρs$⁠, EV#2 has tearing parity (it is the TTEM) and EV#3 has ballooning parity.

This analysis shows that HSX is unstable to TEM turbulence and there is a variety of linear modes that exist in the transport-relevant regime in HSX. HSX-b exhibits multiple TEMs in the range $kyρs≤1$⁠, while HSX-t exhibits only two different modes. The linear growth rates calculated in each flux tube are similar in this range despite differences in mode structure. A quasilinear estimate predicts that despite the large growth rates of the ubiquitous mode, it will not contribute significantly to transport. The existence of a significant number of multiple, subdominant, unstable eigenmodes at each $kyρs$ adds complexity to subsequent nonlinear analysis as multiple modes at the same $kyρs$ can contribute to the turbulence. As a consequence, energy input to the nonlinear state may be significantly greater than what is inferred from the dominant mode alone.

Using the insight gained from linear simulations, nonlinear simulations were performed on the standard TEM and density gradient scan parameter sets. The default nonlinear resolution settings use N_ky = 48 modes and $kyminρs=0.1$⁠, which sets the radial box size according to Eq. (1), where $N=1$ was chosen to yield the smallest radial box size possible. The value of $kyminρs$ was checked for convergence by halving $kyminρs$ and ensuring that the saturated flux does not change. With $kyminρs=0.1$⁠, the radial box size is $Lx=222ρs$ and the radial resolution is set at N_x = 192, but as we will show a posteriori, the nonlinear structures are considerably smaller than the radial box size. Doubling the N_x resolution does not change saturated flux values. The linear resolutions $Nz×Nv∥×Nμ=64×36×8$ also result in nonlinearly converged fluxes.

Quantities of interest, such as the heat flux, are determined, unless otherwise noted, by time averaging over the quasi-stationary state. All heat fluxes reported here have gyro-Bohm normalization, $QgB=csne0Te0(ρs∗)2$⁠, where $ρs∗=ρs/a$ with a is the average minor radius, $ne0$ is the background electron density, and $Te0$ is the background electron temperature. The fluctuating electrostatic potential $Φ$ and electron density n_e are reported as normalized values: $Φ$ in terms of $(Te0/e)ρs∗$ and n in terms of $ne0ρs∗$⁠.

The time traces of the electron heat flux for both HSX-b and HSX-t are shown in Fig. 8. From the time traces, it is apparent that the quasi-stationary electrostatic electron heat fluxes in the two flux tubes are very similar, with $⟨Qes⟩=1.75$ for HSX-b and $⟨Qes⟩=1.59$ for HSX-t. This result is consistent with the arguments made in Sec. III A, where the quasilinear estimate predicts that the $kyρs<1$ modes, where both flux tubes have similar growth rates, dominate the transport.

The electron heat fluxes in Fig. 8 have been decomposed into spectral components in $kyρs$ in Fig. 9, highlighting some moderate differences between the flux tubes while confirming the relevance of the $kyρs≤1$ region. The heat flux spectrum for HSX-t has only one peak at $kyρs=$ 0.5–0.6 and exhibits a quick decay toward high $kyρs$⁠. The spectrum in HSX-b is similar, but the peak heat flux occurs at slightly higher $kyρs=$ 0.7–0.8. The important difference between the flux tubes, however, is the presence of significant flux at $kyρs=0.1$ for HSX-b, which is not predicted by the linear growth rate results of Sec. III A.

In a numerically resolved simulation, it is desired to have only a small amount of flux in the $kyminρs$ wavenumber of the simulation. Fig. 8 shows that HSX-b exhibits no change in the integrated flux when $kyminρs$ is lowered to $kyρs=0.05$⁠. Fig. 9 demonstrates there is little flux in the $kyρs=0.05$ wavenumber but a peak in the heat flux spectrum remains at $kyρs=0.1$ (and a second peak at $kyρs=0.2$⁠). This indicates that the overall flux in this low-$kyρs$ feature is captured in either case. As such, every result hereafter that we will show uses $kyminρs=0.1$⁠.

The nonlinear real frequencies vs. $kyρs$ are shown in Fig. 10 for HSX-b. The linear frequencies have been overlaid and the magnitude of the nonlinear frequencies at higher $kyρs$ is roughly a factor of 3 larger than the linear frequencies. The modes corresponding to dominant transport peak (⁠$kyρs≥0.5$⁠) have negative frequency, which indicates they are TEM. The difference in magnitude between the linear and nonlinear frequencies could be the result of either subdominant modes with higher frequencies or of three wave coupling between different unstable modes. Without any ion temperature gradient, no ITG modes are present in these simulations, however, the modes at low $kyρs$⁠, near to the low-$kyρs$ transport peak, all propagate in the ion direction. It is important to reiterate that no linear subdominant mode with positive frequency was found at these $kyρs$⁠.

The contours of the electrostatic potential and electron density for HSX-b and HSX-t are shown in Figs. 11(a) and 11(b), respectively. These contours are a representative snapshot taken at the final time step of the simulation and are not averaged over the saturated phase. There is zonal flow activity in both flux tubes, as evidenced by the well-defined, alternating vertical bands of positive and negative potential. Zonal density structures are relatively weak in the present case, indicating that they do not play an important role in nonlinear saturation, consistent with corresponding density-gradient-driven tokamak results.^27,28 The density fluctuations do show that the turbulence has a somewhat different character between the two flux tubes. The density fluctuations in HSX-b (Fig. 11(a)) away from x = 0 are anisotropic in the x-y plane with radial elongation not exceeding the zonal flow width, which is consistent with one requirement for eddy shearing by the zonal flows.^29,30 The HSX-triangle density fluctuations (Fig. 11(b)) are more isotropic but are still appear sheared along the zonal flow boundaries.

The appearance of the zonal flow corresponds with the turnover of the heat flux at $t≈15a/cs$ in Fig. 8 and the development of the quasi-stationary state. This indicates that the zonal flows may play a role as the nonlinear saturation mechanism in these simulations; however, the time-averaged, fluctuation-driven $E×B$ shearing rate $ωE=⟨d2Φky=0/dx2⟩≈0.35$ in both flux tubes and is on the order of the respective linear growth rates. This shearing rate has not been corrected for the finite-frequency effects of Ref. 52. Typically, uncorrected values of ω_E as defined are consistent with the rule of thumb established in Refs. 10 and 53, where $ωE≃10γmax$ should be satisfied if the zonal flows are to be an important nonlinear saturation mechanism. A different definition of the shearing rate was used in Refs. 28 and 54, where time scales shorter than an eddy lifetime and spatial scales smaller than the radial correlation length were filtered from fluctuating zonal potential, as a result of which their subsequent shearing rate has to be compared directly to the linear growth rate, without the application of finite-frequency corrections or a factor of $≃10$⁠. To confirm equivalence of both approaches, our present formalism was applied to the data of Ref. 54, and our definition was found to yield a shearing rate value about eight times as large as the other definition, consistent with the aforementioned rule of thumb.

For the present case, the values of ω_E reported here indicate that the zonal flows may not play as big a role as a saturation mechanism. A simulation performed with the zonal flows artificially removed showed that the TEM turbulence still saturates in the absence of the zonal flows but with a quasi-stationary flux level 2.7 times larger. This result confirms that zonal flows are important to the saturation of the TEM in HSX; this is consistent with previous results of density-gradient driven TEMs in tokamaks,^27,53 where the absence of zonal flows leads to larger transport. As a consequence, it is concluded that the use of the shearing rate to determine saturation via zonal flows does not apply in the present case, perhaps due to low magnetic shear or the specifics of the saturation.

In the contours of HSX-b in Fig. 11, there is a large scale coherent structure at x = 0 in the fluctuating potential and density that is not present in HSX-t. This coherent structure is not static and drifts in the -y (ion) direction. Filtering out the $0.1≤kyρs≤0.3$ components in Fig. 11 removes this structure and reveals an underlying zonal flow smaller in amplitude than the coherent structure and similar in amplitude to the zonal flows at other radial positions in the box. The large flux seen in the flux spectrum of Fig. 9 at $kyρs≈0.1$ is a product of this coherent structure.

The coherent structure observed in Fig. 11(a) shows a significant dependence on several parameters: the presence of zonal flows, β, and the background magnetic shear. In the zonal-flow-free simulation, the coherent structure was also absent, suggesting that the evolution of the coherent structure may be moderated by the zonal flow. Fig. 12 shows that the low-$kyρs$ peak in transport is eliminated when β is changed from $β=5×10−4$ to $β=5×10−5$⁠. The time trace (not shown) of heat flux has slightly higher quasi-stationary flux levels for the lower-β case, consistent with the linear result where growth rates throughout most of the $kyρs$ range are slightly stabilized with increasing β. Simulations performed with an artificially increased background magnetic shear (not shown) have reduced the low-$kyρs$ flux and lack a nonlinear coherent structure at x = 0, while the integrated flux levels are not affected as much.

The cross phases of $Φ$ with n_e, $Te∥$⁠, and $Te⊥$⁠, shown in Fig. 13, show a clear difference between linear and nonlinear simulations at low $kyρs$ in HSX-b. In all three phase relations, the potential is out of phase with the other quantities at $kyρs≤0.2$⁠, corresponding to the coherent structure. As $kyρs$ increases, the linear and nonlinear phases agree reasonably well for $kyρs>0.3$⁠. The $Φ$ vs. n_e and $Φ$ vs. $Te⊥$ cross phases, Figs. 13(a) and 13(c), follow the same trend as the trapped-particle population in Ref. 55 for temperature-gradient-driven turbulence, where the phase angle increases with increasing $kyρs$⁠. This indicates that despite the linear-nonlinear frequency mismatch, the turbulence is of TEM type. Next, this interpretation shall be confirmed by comparing dependencies in the driving density gradient.

The nonlinear electrostatic electron heat fluxes for the density gradient scan parameters, where $a/LTe=0$⁠, are presented in Fig. 14. These values are on the order of the observed experimental fluxes obtained from calculations of the power deposition profile,⁵⁶ where $3≲a/Ln≲4$ and $2≲⟨Qes exp ⟩/QgB≲7$ in the region $0.4≤r/a≤0.6$⁠. This suggests that TEM turbulence is responsible for the observed heat fluxes in HSX.

The scaling of the nonlinear saturated heat flux with driving gradient is shown in Fig. 15 and compared with the linear growth rates at $kyρs=0.7$⁠, approximately where the nonlinear heat flux spectrum peaks in Fig. 16. A linear fit has been applied to the nonlinear flux data for $a/Ln≤3$⁠. An upshift in critical density gradient is observed for the nonlinear saturated fluxes, consistent with the TEM nonlinear upshift discovered in tokamak simulations,²⁶ with a critical density gradient at $a/Ln≈0.8$ as compared with a critical density gradient at $a/Ln≈0.2$ predicted by the linear growth rate scaling. The fact that there is an upshift again provides evidence that zonal flows play a role in moderating the turbulence.²⁸ 

The electron heat flux spectrum presented in Fig. 16 shows that the fraction of the flux contributed by the low-$kyρs$ wavenumbers increases with driving density gradient. At $a/Ln=4$⁠, the fraction of flux for $kyρs≤0.3$ is approximately 22% of total integrated flux. Again, the nonlinear frequencies for the low-$kyρs$ modes are ion frequencies, consistent with Fig. 10.

The increase in $a/Ln$ leads to an increase both in the variability of the heat flux and in the transport in the lowest-$kyρs$ wavenumbers, the cause of which can be determined by examining the contours of electrostatic potential, shown for all four $a/Ln$ in Fig. 17. The potential structures in Figs. 17(a) and 17(b) are similar to both Fig. 11(a) and Fig. 11(b) for standard TEM parameters, showing clear zonal flows. However, there is no coherent structure observed for $a/Ln=1$ and subsequently no low-$kyρs$ transport peak in Fig. 16. For $a/Ln=2$ (Fig. 17(b)), the coherent structure is located near the radial boundary of the box and is responsible for the low-$kyρs$ transport peak in Fig. 16. For $a/Ln=3$ and $a/Ln=4$⁠, Figs. 17(c) and 17(d) show a large amplitude potential structure centered at x = 0 that drifts in the ion direction and locally dominates the zonal flows.

Decreasing $kyminρs$ to 0.05 (equivalent to doubling the box size in the y direction) for the $a/Ln=4$ simulation shows the coherent structure at x = 0 retains a similar scale as in Fig. 17(d) and that the quasi-stationary flux level is similar to the $kyminρs=0.1$ simulation, indicating that the $a/Ln=4$ simulations presented here are sufficiently resolved to capture the turbulent dynamics at low $kyρs$⁠. The bursts and fluctuations seen in Fig. 14 for $a/Ln=4$ are also present in the higher resolution simulation due to the large size and amplitude of the coherent structure at x = 0.

The same observations about the zonal flow as the saturation mechanism as for the standard TEM parameters apply here. For $a/Ln≤3$⁠, the shearing rate ω_E is always smaller than the maximum linear growth rate and only slightly larger than the maximum linear growth rate for $a/Ln=4$⁠. Therefore, zonal-flow-based turbulent saturation may rely on energy transfer to damped modes rather than shearing. The associated density fluctuations to Fig. 17 are comparable to the potential fluctuations but zonal density as a saturation mechanism is again precluded by the absence of zonal structures in the density.

We have performed an in-depth gyrokinetic study of density-gradient-driven TEM turbulence in the HSX stellarator. Linear simulations show that HSX profiles are unstable to TEMs across a large $kyρs$ range, with the HSX flux tube corresponding to the “bean” shape having larger growth rates than the HSX flux tube corresponding to the triangular cross section. In the investigated range $0.1≤kyρs≤6$⁠, we have found four distinct dominant TEM-type modes in HSX-b and two in HSX-t. The TEMs at low $kyρs$ have a two-scale ballooning structure, with an extended envelope whose width is influenced by the low shear and localized spikes in amplitude due to localized regions of overlapping bad curvature and particle trapping. At higher $kyρs$⁠, the modes become substantially more localized, and an analysis of the five most unstable modes for each $kyρs$ in HSX-b shows that there is a distinct tearing parity TEM that transitions from dominant to subdominant as $kyρs$ exceeds 0.9. HSX-t does not exhibit a dominant tearing parity TEM at any point.

Nonlinear TEM simulations show that the two flux tubes have similar quasi-stationary flux levels, in agreement with quasilinear estimates based on the similarity of growths for $kyρs<1$⁠. The $kyρs$ flux spectra between the two flux tubes only differ in that there is a peak in the flux at $kyρs=0.1$ in HSX-b, which does not occur for HSX-t and is not predicted by a quasilinear estimate. This low-$kyρs$ transport peak is the consequence of a robust coherent structure in HSX-b, which is further destabilized by $a/Ln$ and β and can be eliminated by increasing the background magnetic shear. The characteristic frequencies and phase relations associated with the coherent structure have no equivalents in any of the performed linear calculations. Examination of the potential contours shows that the nonlinear coherent structures are radially localized to the resonant magnetic surface at the center of the simulation domain.

Zonal flows were present in all nonlinear simulations, consistent with the tokamak-based expectations for density-gradient-driven TEM. Artificially removing the zonal flows confirms that the TEM is saturated by zonal flow interactions; however, the $E×B$ shearing rates are much smaller than what one would expect in turbulence that saturates through zonal flow shearing. Instead, energy transfer between different modes and moderated by the zonal flows may be responsible for nonlinear saturation. This, however, has not been confirmed and will be left for future investigation.

Simulated turbulent heat fluxes are comparable with experimental observations. While additional study along these lines will be necessary—e.g., investigating the impact of non-unity temperature ratio—it is concluded that density-gradient-driven TEM microturbulence is an excellent candidate in explaining heat fluxes in the HSX stellarator.

The authors would like to acknowledge H. E. Mynick and D. R. Ernst for useful discussions. This work was supported by U.S. Department of Energy Grant Nos. DE-FG02-93ER54222 and DE-SC0006103. This research used resources of the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231. The authors are extremely grateful to D. R. Mikkelsen and F. Jenko for generously providing computing resources.