Ion temperature gradient turbulence modification in quasi-axisymmetry

The numerical optimization of stellarator configurations for improved transport properties¹ has yielded devices with tokamak levels of neoclassical transport^2–4 and suggests that devices with reduced turbulent transport may be developed as well.⁵ Stellarators have been shown to benefit from three-dimensional effects on a flux surface, causing a reduction of turbulence levels, compared to the values located around the largest normal curvature.⁶ The QUASAR stellarator provides a setting in which adiabatic Ion Temperature Gradient (ITG) turbulence may be explored.^7,8 In this work, a set of five experimentally realizable equilibria are developed which scan magnetic shear $(dι/dρ)$ at fixed edge rotational transform (ι). This produces configurations which vary from traditional low-shear stellarator configurations to the ones with a “tokamak-like” rotational transform profiles. The resulting variation in ITG turbulence properties is investigated using linear and nonlinear gyrokinetic calculations, suggesting new modifications to existing ITG turbulence proxy functions.

The QUASAR stellarator is a re-envisioning of the cancelled National Compact Stellarator Experiment (NCSX).^9–11 In QUASAR, the NCSX coil set and vacuum vessel are reused with the possibility of additional coils. For this work, the auxiliary ohmic and vertical field coil sets, as envisioned for NCSX, are utilized. It should be noted that the optimized modular stellarator coils and toroidal field coils have already been manufactured. Thus, the QUASAR facility is a quasi-axisymmetric stellarator capable of Ohmic current drive and heating. The device is envisioned to be run with stable discharge conditions for 3 to 5 s.

This study showcases the unique capabilities of QUASAR, provides an experimental path towards validating turbulence predictions in stellarators, and hints at the possibility of a fully transport optimized stellarator reactor. Unlike modern high performance stellarators (LHD and W7-X), the QUASAR stellarator has a large configuration space.¹² This is in part due to its extensive coil set which includes three optimized stellarator coils, toroidal field coils, poloidal field coils, and ohmic coil. The ability to drive Ohmic current is shared by only one other stellarator in the world, the Compact Toroidal Hybrid.¹³ The low levels of predicted neoclassical transport allow QUASAR to serve as a stellarator turbulence test-stand. Moreover, the quasi-axisymmetric nature of the device should allow comparison with existing tokamak turbulence studies.¹⁴ Finally, if predictions of a turbulence optimized stellarator can be validated in experiment, the stellarator reactor design may become preferable to tokamak configurations. In Sec. II, we discuss the tools and methods used in this work. In Sec. III, we present the results of this analysis followed by a discussion of the results in Sec. IV.

An exploration of ITG turbulence in the QUASAR device is conducted using the VMEC,¹⁵ STELLOPT,^16,17 and GENE^18,19 codes. The VMEC code is used to solve for free-boundary ideal MHD equilibria with continuously nested flux surfaces in the QUASAR stellarator. Care is taken to tailor profiles which are consistent with a β ∼ 2% ohmic discharge. The STELLOPT code is used to optimize the equilibrium inputs (coil currents, current, and enclosed flux) so that ITG turbulence can be evaluated through a shear scan in the device. Edge rotational transform and plasma volume serve as targets. The resulting equilibria are then evaluated at multiple radial locations to provide estimates of ITG driven heat flux. Such an analysis highlights the unique capabilities of the QUASAR facility to explore turbulence and its optimization.

The VMEC code is used in this work to calculate the free boundary three dimensional ideal MHD equilibrium of the QUASAR device. This is achieved by minimizing the ideal MHD energy functional

subject to the nested flux surface constraint $B→·∇ψ=0$⁠, where p is the plasma pressure, γ is the adiabatic index, $B→$ is the magnetic field, ψ is the toroidal flux, μ₀ is the permittivity of free space, and integration is over the domain of the plasma equilibrium. Free boundary simulations are made possible through consideration of the vacuum field energy on the plasma boundary.²⁰ Assuming temperature profiles of the form $Tk(s)=Tk0(1−s2)$ (where s = ψ/ψ_edge is the normalized toroidal flux and index k = i, e denotes the ion and electron species), a β = 2% equilibrium can be constructed. In this scenario, T_i0 = 3.5 keV, T_e0 = 1.5 keV, and the electron density is assumed to be constant at 5 × 10¹⁹ m⁻³. These profiles are chosen to mimic an ohmic discharge in the experiment. As a result, the current density profile has a form $I′(s)=(1−s2)$ consistent with the electron temperature. All profiles are specified using the 100 point Akima spline parameterization.²¹ The net toroidal current is left as a free parameter allowing edge iota to be held fixed, while the magnetic shear is varied.

A series of ideal equilibria which scan magnetic shear at fixed edge rotational transform are generated using the STELLOPT code. The modified Levenberg-Marquardt algorithm^22,23 is used to target an edge rotational transform value of $ι=0.616$ for five equilibria at different levels of auxiliary toroidal field coil current. In order to simplify the optimization, the three modular coil types of the stellarator coil-set are treated as one coil system, holding their relative strengths fixed. Table I depicts the target parameters used in the optimization and their relative weightings. In Sec. III, the resulting five equilibria are described.

The end goal of this work is to demonstrate that the QUASAR facility can be used to test adiabatic ITG turbulence. This is done through evaluation of these equilibria using the gyrokinetic code GENE. Recently, the STELLOPT code has been interfaced to the GIST code,²⁴ allowing it to directly calculate the geometric quantities necessary for linear and nonlinear runs of GENE. A comparison between the standalone GIST computation and that incorporated into STELLOPT can be seen in Fig. 1. This has allowed the development of ITG and TEM proxy functions for turbulence^5,25 and the inclusion of linear GENE calculations directly in the optimizer. This last feature has been utilized in this work to simplify calculations of the radial profiles of linear growth rates. Standalone nonlinear flux tube runs are performed for select flux tubes. Such simulations are still too computationally expensive to directly incorporate in the optimizer.

In Fig. 1, various magnetic and geometric quantities are benchmarked between the STELLOPT and GIST implementation. The contravariant metric elements are defined as

where s = Φ/Φ_edge is the normalized toroidal flux, $∇α=q′θ*∇s+q∇θ*−∇ϕ$ (prime is derivation with respect to s), a is the effective minor radius, θ* is the VMEC poloidal coordinate, ϕ is the VMEC toroidal coordinate, $B̂=|B|/|Φedge/πa2|$⁠, and $ι=1/q$ is the rotational transform. All these values are in excellent agreement between the two codes. Differences in the equations presented here are due to the utilization of PEST coordinates (as opposed to Boozer in Ref. 24) The quantities labeled L2 and L1 are curvature terms related to the normal and geodesic curvature (respectively), while dBdt refers to the mirror term.

These quantities serve as input for the GENE gyrokinetic simulation runs. All gene simulations utilized adiabatic electrons and a gradient scale length of $a/LTi=3$ unless otherwise noted. This simulation grid was chosen to 128 points the parallel direction, 128 radial Fourier modes, and 96 poloidal Fourier modes. The parallel and perpendicular velocity space grids had 32 and 16 grid points, respectively. The radial and binormal extents of the flux tube were chosen to be 125 gyro radii, and the parallel extent is equal to 2πa. In velocity space, the extent of the velocity grids is 6v_therm and $9Ti/(Φedge/πa2)$ in the parallel and perpendicular directions. The extent of the grid in the radial x direction is automatically adjusted to fulfill a quantization constraint on the grid (involving shear, x grid size, and the largest resolved gyroradius). In cases of a low magnetic shear, it is appropriate to set the shear to zero and then assume parallel boundary conditions in the direction z (along the magnetic field line). This avoids unnecessarily expensive computations with large grids in the x direction.

In this work, five ohmic discharge type scenarios, varying shear but holding edge rotational transform fixed, are developed using the STELLOPT code. These scenarios are then evaluated for their adiabatic ITG turbulence properties using proxy functions, linear growth rates, and nonlinear flux tube simulations. The variation of the turbulent heat flux is found to exceed factors of two at multiple radii in the plasma.

The STELLOPT optimizer was used to develop the five scenarios based on a variation of the axisymmetric toroidal field generated by the auxiliary toroidal field coils of the QUASAR device. Table II lists the free parameters for the optimization and the optimized values. All these values appear reasonable for the QUASAR experiment. Figure 2 depicts the evolution of χ² during the optimization for the 100 kA case. The Levenberg method clearly searches along two paths during the first iteration (initial equilibrium denoted as (1). The two paths equate to Gauss-Newton trajectory and the gradient descent. The descent is dominated by the edge rotational transform constraint, although all targeted quantities are minimized. This is because the contributions to χ² are fairly orthogonal with respect to their dependence on the varied parameters.

This rather simple optimization task resulted in five unique equilibria for the QUASAR device. Figure 3 depicts rotational transform profiles for these five equilibria. The magnetic shear has clearly been varied while holding the edge rotational transform fixed. While one may consider the case with zero auxiliary toroidal field (0 kA) to be the “pure stellarator case,” the rotational transform suggests a different labeling of the equilibria. Specifically, the configuration where the applied toroidal field opposes that generated by the modular coils (−100 kA) has a more stellarator-like profile. Thus, it will be referred to as the “stellarator-like” case for this work. The configuration in which the applied toroidal field is in the direction of that generated by the modular coil has a more tokamak-like profile. Thus it will be referred to as the “tokamak-like” case for this work (100 kA).

Examination of the five configurations in more detail highlights some interesting features. The neoclassical transport (Fig. 4), as characterized by the helical ripple parameter (⁠$ϵeff3/2$⁠), is low in all five configurations.²⁶ This suggests that transport for these configurations will be dominated by turbulence. Infinite-n ballooning stability analysis suggests that the configurations are mostly stable to ballooning modes (with isolated radial regions being marginally unstable, ρ > 0.8). The 50 and 100 kA cases have a small radial region near the edge where they become ballooning unstable.²⁷ Overall, these equilibria all appear to be good candidates for real experimental configurations which can test the response of turbulence to magnetic shear changes.

Before moving on to the studies of turbulent transport, a digression regarding our neglect of a bootstrap current should be presented. The neglect of a bootstrap contribution to the current profile was motivated out of convenience. Namely our equilibrium code requires that we specify a net toroidal current and a profile for dI/ds in terms of normalized toroidal flux. For an ohmic discharge, assuming a current profile in the shape of the electron temperature is a reasonable simplification. This of course does not take the place of a full transport simulation, but such tools are still in their infancy for stellarators.²⁸ In addition, using configurations with relatively low plasma β, the effects of bootstrap currents are subsequently minimized. Figure 5 depicts the bootstrap current as calculated in the collisionless limit for the 0 kA applied toroidal field case.²⁹ The plotted profile has been smoothed to remove resonant features. Inclusion of a bootstrap contribution to the current profile modifies the edge current profile, which mostly changes the edge iota (a minor correctible change). If a collisional correction is taken into account, the effect of the edge bootstrap would be further reduced. This is because collisionality (ν*) is inversely proportional to temperature, thus reducing the edge bootstrap contribution. Neglecting the bootstrap contribution is a justified simplifying assumption in this work.

Using the previously mentioned coupling between the optimizer and gyrokinetic tools, profiles of proxy and linear growth rates were evaluated. Here, we have used the $Qpr1d=〈γkg11〉$ (prox-1d) proxy,⁵ where γ_k is the growth rate. Brackets $〈 〉$ indicate averaging over the flux surface, and sums are taken over wavenumber k. Figure 6 depicts the turbulent transport proxy used to generate optimized stellarator configurations based on the QUASAR fixed boundary equilibrium.⁵ The more “tokamak-like” cases show an increase in the proxy around s = 0.75. This is also the region where these equilibria become infinite-n ballooning unstable (as calculated by the COBRA code²⁷) The presence of this feature is not surprising as the proxy has a ballooning-like character. The proxy results appear inconclusive in the core region (s < 0.25). At larger radii, the proxy function seems to indicate that as configurations become more “tokamak-like” turbulent heat flux increases. However, care should be taken as only general trend information can be extracted from this proxy and the relative amplitude comparison made here may be very approximate.

The GENE gyrokinetic code was employed to calculate various ITG growth rates. In these simulations, the electrons were treated as adiabatic. Flux tube simulations were performed to calculate the linear grow rate on the VMEC mesh for each equilibrium. Using these simulations, a profile of linear growth rates was constructed to provide a better comparison of the equilibria (Fig. 7). Here, we see that the growth rates for modes in the core (s < 0.25) are reduced as the configuration becomes more “tokamak-like.” At larger radii, this trend appears to reverse with the “stellarator-like” configuration showing the smallest growth rate. It should be noted that the ballooning like feature is no longer present for the “tokamak-like” configurations and that here we are only reporting the growth rates of the most unstable modes. Recent work suggests that proper estimation of the saturated heat flux requires an ensemble of modes to be considered.³⁰ 

Nonlinear flux tube simulations were performed for a selected set of flux tubes (with adiabatic electrons). Figure 8 (left) compares the evolution of the nonlinear heat flux for three configurations in the core (s ∼ 0.1). The more “tokamak-like” configuration shows a significantly reduced heat flux and the standard deviation in the heat is much smaller as well. Figure 8 (right) shows a similar comparison between heat fluxes but in the outer plasma region (s ∼ 0.5). As in the core, the heat flux is reduced for the more “tokamak-like” configuration. This suggests that both the proxy and most unstable linear mode fail to recover the actual behavior of the turbulence, probably due to the disparate scales involved.

The trend where more “tokamak-like” configurations have reduced turbulence levels is confirmed through further simulations. Figure 9 shows that at all radii the outboard mid plane flux tube has reduced turbulence levels for “tokamak-like” configurations. This is the flux tube which passes through the region of maximum curvature in the equilibria. In these plots, the data points are extracted as averages of the turbulent fluctuation levels after saturation. The error bars depict the standard deviation of the turbulent heat flux after saturation. Scans of the gradient length scale suggest that the relative amplitudes of the fluxes hold. Gradient length scales $(a/LTi)$ of 1, 2, and 3 were compared for a core flux tube. In each case, at least a factor of two difference in heat flux was observed. Given the ohmic nature of these scenarios, the effect seen here will mostly likely be robust to profile variation. Full surface nonlinear runs have also been conducted and corroborate these results, suggesting that the “tokamak-like” configuration has lower levels of turbulent heat-flux. Further analysis would require global simulations which are beyond the scope of this work.

In this work, a series of magnetic configurations are developed which show the possibility of testing adiabatic ITG turbulence in the QUASAR device. This series of five equilibria explore the effect of magnetic shear on turbulence where the edge iota is held fixed. Profile assumptions are made which assume ohmic heating and current drive, a unique feature of the QUASAR experiment. The auxiliary toroidal field coils of the experiment play an important role in these configurations. Nonlinear gyrokinetic simulations of adiabatic ITG turbulence confirm a consistent reduction in turbulent transport for configurations deemed to have a “tokamak-like” rotational transform profile. This suggests that, should the device be built, QUASAR could play a major role in the validation of stellarator gyrokinetic theory.

In this work, an ohmic β ∼ 2% discharge is assumed for simplicity. A full transport simulation would provide a much more convincing basis for the study. However, in this work nearly an order of magnitude changes in turbulence levels are present at all radii. This suggests that while these profile assumptions may be approximate, the response in terms of transport most likely swamps any profile discrepancies. Moreover, evaluation of the bootstrap current suggests that its contribution would be small compared to the ohmic heating profile. Additionally, we have neglected the role non-adiabatic electrons may play in the mode development via additional trapped electron drive. Recent works on Wendelstein 7-X suggest that trapped electron modes are subdominant to other modes of transport, owing to the periodic field nature of stellarators.³¹ Their role in ITG turbulence is still an open question, and we leave investigation of this effect to future works.

Differences between the proxy function, linear growth rates, and turbulent heat flux should also be addressed. The proxy functions do not claim to predict turbulent heat flux levels, but rather mimic the dependence of turbulence on magnetic field properties. These functions become essentially normalized, and the comparison presented here is illustrative of this fact. Specifically they may provide gradient information regarding what optimization parameters reduce turbulent transport, but do not provide absolute values which allow comparison of different configurations. One may ask if the optimizer had been started from the “stellarator-like” configuration and optimized using the proxy, would it create a configuration similar to the “tokamak-like” case? Such work is left to future analysis. Discrepancies between the linear growth rates are indicative of missing physics in the proxy. Recent work has suggested that the saturated heat flux may depend on a combination of linear modes.³⁰ Our analysis focused on just the most unstable mode.

Overall, this work motivates the necessity of a device like QUASAR to better explore ITG turbulence through the variation of magnetic configuration. This work also suggests that stellarators with more “tokamak-like” rotational transform profiles will have improved transport properties. This may suggest a new field of optimization for quasi-axisymmetric stellarators. Can the stellarator bootstrap current be optimized to provide similar rotational transforms? The answers to such questions are left to future work.

This manuscript is based upon the work supported by the U.S. Department of Energy, Office of Science, Office of Fusion Energy Sciences, and was authored by Princeton University under Contract No. DE-AC02-09CH11466 with the U.S. Department of Energy. The publisher by accepting the article for publication acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript or allow others to do so, for United States Government purposes.