Intrinsic turbulence stabilization in a stellarator

The magnetic surfaces of modern stellarators are characterized by complex, carefully optimized shaping and exhibit locally compressed regions of strong turbulence drive. Massively parallel computer simulations of plasma turbulence reveal, however, that stellarators also possess two intrinsic mechanisms to mitigate the effect of this drive. In the regime where the length scale of the turbulence is very small compared to the equilibrium scale set by the variation of the magnetic field, the strongest fluctuations form narrow band-like structures on the magnetic surfaces. Thanks to this localization, the average transport through the surface is significantly smaller than that predicted at locations of peak turbulence. This feature results in a numerically observed upshift of the onset of turbulence on the surface towards higher ion temperature gradients as compared with the prediction from the most unstable regions. In a second regime lacking scale separation, the localization is lost and the fluctuations spread out on the magnetic surface. Nonetheless, stabilization persists through the suppression of the large eddies (relative to the equilibrium scale), leading to a reduced stiffness for the heat flux dependence on the ion temperature gradient. These fundamental differences with tokamak turbulence are exemplified for the QUASAR stellarator [G. H. Neilson et al., IEEE Trans. Plasma Sci. 42, 489 (2014)].

Over the last several decades, magnetic confinement fusion research has been dominated by two concepts, the tokamak and the stellarator. Most of the effort has gone into tokamaks, which employ axisymmetric magnetic fields to confine the plasma, but it has always been recognized that stellarators enjoy certain advantages although they use more complicated fields, and therefore face a number of technical challenges [1]. Compared to tokamaks, the physics of stellarator plasmas is less understood, and the greatest uncertainty concerns the nature of the plasma turbulence. This topic is of particular interest to the Wendelstein 7-X stellarator, which has recently been commissioned and completed its first-plasma milestone [2,3]. Its performance is expected to depend crucially on the turbulence present in the plasma, which will be explored thoroughly in upcoming experimental campaigns.
The great challenge for both tokamaks and stellarators is to provide stable and robust plasma confinement with minimal energy losses. The latter are caused by plasma turbulence and by so-called neoclassical transport, which results from the random walk taken by plasma particles moving along complicated orbits in the magnetic field whilst colliding with each other [4]. Neoclassical losses are modest in tokamaks but tend to be prohibitively large in stellarators, unless the magnetic field geometry is optimized to reduce the neoclassical transport [5], in which case the turbulent loss channel becomes relatively more important. Both for W7-X and other modern stellarators, such as QUASAR (based on the NCSX quasi-axisymmetric design from Princeton Plasma Physics Laboratory, which to date has not been assembled) [6], the result of the optimization is a magnetic field with quite complex geometry [7]. The geometric properties of the field lines vary greatly over the magnetic surfaces, and the latter are pushed against each other in some locations, resulting in locally large temperature gradients and consequently strong drive for turbulence. The natural question that arises is how this geometric complexity affects the turbulent transport and how the latter compares to that in tokamaks.
Turbulence in tokamaks and stellarators is generally thought to be caused by low-frequency plasma instabilities, such as the electrostatic, collisionless ion temperature gradient (ITG) driven mode [8,9]. A theoretical method to understand and, ideally, predict the behavior of turbulence is furnished by "gyrokinetics" [10], according to which, the fully kinetic description of the plasma is reduced by one dimension, owing to the fast gyration of the charged particles around the magnetic field lines. Despite this simplification, the coupled system of nonlinear partial differential equations is five-dimensional (plus time) and can only be solved by numerical codes, such as GKV-X [11], GS2[12] and GENE [13], with the help of modern supercomputers. We employ the massively parallel GENE code, which is able to address both tokamak (covering the entire plasma radius) [14] and stellarator (covering the entire magnetic surface, while radially local) [15] configurations. For all simulations shown, the plasma ions are treated gyrokinetically while the electrons are assumed to have a Boltzmann (adiabatic) distribution (the inclusion of gyrokinetic electrons is at present intractable due to the immense computational resources required). The dimensionless parameter expressing the strength of the ITG turbulence drive is 1/L Ti = −1/T i dT i /ds (T i is the ion temperature). Here, s = Φ tor /Φ lcfs ∈ [0, 1] is the radial coordinate (Φ tor is the toroidal magnetic flux, with Φ lcfs its value at the outermost surface); for all simulations we have selected s = 0.5, which corresponds to a radius where tur-bulence is expected to be particularly active. Although several physical mechanisms for ITG turbulence stabilization have been reported in tokamak-related gyrokinetic simulations, for instance involving electromagnetic effects [16,17] or geometric effects like triangularity [18], here we address intrinsic stabilization effects that are due solely to the nonaxisymmetric magnetic geometry. Specifically, the pivotal point in the present work is that, in a stellarator, the transport averaged over a magnetic surface can be significantly lower than the transport produced locally at the most unstable regions on the surface.
The effects we consider encompass a wide range of values for the dimensionless parameter ρ * = ρ/a, where ρ is the ion gyroradius and a the minor radius of the torus. In the "small-ρ * regime" characterized here by ρ * = 1/250, the turbulence scale is much smaller than the variation of the equilibrium magnetic field. Turbulence thus appears localized on the surface, in that the strongest fluctuations accumulate along thin band-like structures, as shown in Fig. 1 for the QUASAR configuration. This localization FIG. 1. (Color) Root-mean-squared normalized density fluctuations caused by ITG turbulence (1/LT i = 2) on the magnetic surface of the QUASAR stellarator in the small-ρ * regime. The strongest fluctuations appear localized around the magnetic field line with α = π/3 (black), and similarly for the line with α = −π/3 (not shown), due to the symmetry of the configuration. The number of grid points in the five-dimensional phase space (x, y, z, v , µ) is 128 × 128 × 128 × 64 × 10. The simulation box size for the radial direction is Lx = 147ρ and for the binormal direction Ly = 204ρ. Also, Lz = 2π for the parallel direction (parameterized by the poloidal angle). The simulation box size for the parallel velocity is −3vt ≤ v ≤ 3vt and for the magnetic moment 0 ≤ µ ≤ 9Ti/B0. Here, vt = 2Ti/mi is the thermal velocity and B0 = 2Φ lcfs /a 2 .
is responsible for a reduction of the transport through the surface as compared with its local maximum, since large relatively quiescent regions bring down the average level. Figure 2 presents ion heat flux timetraces in gyroBohm units, Q gB = c s P (ρ * ) 2 (c s = T i /m i is the ion sound speed and P is the plasma pressure). The timetraces for the most unstable magnetic field line with α = π/3 and the surface average are displayed. [The coordinate α = ζ − qθ labels a magnetic field line cover-ing one poloidal turn; θ, ζ are the poloidal and toroidal angles and q is the safety factor, defining the amount of winding of the magnetic field line on the surface.] Also shown in Fig. 2 is the timetrace from a location of weak turbulence, at the magnetic line with α = 0. It is found that the resulting transport on the surface lies between the two field line levels, demonstrating the averaging effect. There is an important observation concerning the small-ρ * regime, namely that, despite the deviation of the surface-averaged transport from the strongest local one, ITG turbulence still behaves locally with respect to a field line. Indeed, if we use a slender computational domain surrounding a single field line, a so-called "flux tube" [20], then locality is to be understood in the sense that the transport stemming from a magnetic field line on the surface is equal to that produced separately, using a flux tube around this line. We stress, however, that locality in the stellarator does not imply that the transport calculated from a single flux tube can predict the surface-averaged transport. on the magnetic surface of the QUASAR stellarator in the small-ρ * regime. The timetraces correspond to the field line with α = π/3 (strong turbulence region), the field line with α = 0 (weak turbulence region) and the surface average, demonstrating the stabilization (averaging) effect. Also shown is the timetrace from a separate flux tube simulation surrounding the field line with α = π/3, suggesting the local nature of stellarator turbulence in small-ρ * conditions.
To appreciate the critical role of nonaxisymmetry for turbulence stabilization, we can compare Fig. 1 with Fig.  3, corresponding to a tokamak with similar aspect ratio (3.87) as QUASAR (4.47). Evidently, in the small-ρ * regime, no localization of the ITG turbulence exists in the tokamak, since the axisymmetric geometry induces an almost even distribution of the fluctuations on the outboard side of the surface. As a result, the surfaceaveraged heat flux is essentially the same as the heat flux measured at various field lines, implying a lack of this stabilization for the tokamak (see Fig. 4). The ben- eficial turbulence localization in QUASAR could be lost, however, in the "large-ρ * regime" characterized here by ρ * = 1/125, where the ITG density fluctuations occupy the entire outboard domain of the configuration (see Fig.  5). In this regime, the distribution of turbulent intensity resembles that of a tokamak. Be that as it may, ITG turbulence still enjoys stabilization, as another mechanism applies, once ρ * increases, via interaction of turbulent eddies with the rapidly-varying magnetic field in a 3D configuration. A theoretical explanation of this mechanism can be obtained in the fluid limit, by examining the linear drive of the electrostatic fluctuations φ(α) of the

ITG mode with frequency ω,
where T e is the temperature of the electrons, ω T * = (qT i /e) d ln T i /dΦ tor , and the magnetic drift frequency ω d = ∇α · (b × ∇ ln B)2T i /(eB) introduces nonaxisymmetry via its field line dependence (b is the magnetic field unit vector and B the magnetic field intensity). Due to 3D magnetic shaping, geometric features such as magnetic shear and curvature cause the magnetic drift to vary in α (unlike in a tokamak), so that it exhibits areas both favoring (ω d < 0) and suppressing (ω d > 0) local instability. To model this effect, we can consider an analytically tractable form of ω d (α) given by a periodic, piecewise constant function, alternating between regions of positive and negative values, with the extent of these regions representing the equilibrium scale associated with the variation of ω d . The solutions of Eq. 1 are in this case waves (∂ α → ik α ), but since ω d multiplies the derivative term, finiteness of ∂ 2 φ/∂ 2 α requires that zeros in the ITG eigenmode φ(α) occur where zeros of ω d exist. Thus, in analogy to a standing wave on a string, a minimum allowable wave number k α is imposed, here determined by the equilibrium scale of the magnetic drift. In short, the size of the mode is limited by the size of the unstable region where it resides. At small ρ * this is of no consequence, since the length scale of the turbulence is much smaller than the equilibrium scale. At sufficiently large ρ * though, this scale separation is lost, and the lower bound on k α implies that large eddies, that would carry most of the fluctuation energy, are suppressed. In a tokamak, turbulence stabilization due to ρ * is expected to be much weaker compared to the stellarator (radial effects that can additionally influence the tokamak transport are not taken into account [21]): although axisymmetric configurations formally also im-pose a geometric constraint on the scale of turbulence (k α ), this is merely set by toroidal periodicity, rather than smaller-scale 3D magnetic shaping via ω d (α). As a consequence, significant turbulence stabilization for the tokamak requires much larger gradients and ρ * values compared to the stellarator. This is evidenced by 6, which serves as a summary of the effect of the turbulence stabilization mechanisms in QUASAR and the tokamak. For the stellarator, ion heat flux scaling curves are shown as a function of the ion temperature gradient, for the two values of ρ * , together with the curve corresponding to the most unstable magnetic field line (α = π/3) on the surface. We note that stabilization in the large-ρ * regime primarily causes a favorable reduction of the transport stiffness (the rate at which the transport increases), while in the small-ρ * regime, stabilization is mainly responsible for the upshift observed by comparing the threshold of the surface-averaged transport to that of the most active field line. We emphasize that this effect is solely related to the nonaxisymmetric geometry, and note, as Fig. 4 suggests, that no such upshift should be reported in tokamak gyrokinetic simulations, as the heat flux levels for the surface and any magnetic field line are practically identical (these data points do not appear in Fig. 6 to avoid clutter).
The newly observed upshift is comparable in magnitude but otherwise unrelated to the well-known Dimits shift [22,23] in the context of ITG turbulence, which refers to the difference between the linear and nonlinear thresholds; the upshift in Fig. 6 instead involves two nonlinear scalings. In fact, a typical Dimits shift is also found in QUASAR, by comparing the local linear stabil- ity threshold (≈ 1.5 in Fig. 7) with the nonlinear threshold (≈ 1.65 in Fig. 6) measured on the corresponding field line. As the nonlinear scaling is derived from the surface simulation, we conclude that zonal flow activity should be also present on the entire magnetic surface. Interestingly, however, the linear threshold of the surface mode (≈ 1.85 in Fig. 7) appears to coincide with the small-ρ * nonlinear threshold of the surface transport. This appar-   contradiction is resolved by noting that the global linear mode stabilization should not be related to nonlinear stabilization (heat flux reduction), as the global linear mode involves communication across the entire flux surface, contrary to the local behavior of turbulence in the small-ρ * regime. In this context, we also warn against inferring a direct connection between the upshift of the nonlinear threshold and the corresponding linear one. As demonstrated in Fig. 8, even in the large-ρ * regime, the scale of the corresponding linear mode is very different from the turbulence scale. For that matter, one can observe that the heat flux reduction is disproportionately large relative to the growth rate reduction.
Summarizing, we have highlighted two different mechanisms by which a stellarator is able to stabilize ITG-driven turbulence on the surface, solely thanks to nonaxisymmetry. The first mechanism applies in the small-ρ * regime, where ITG fluctuations are localized in the vicinity of certain magnetic field lines on the surface, whereas the second is relevant in the large-ρ * regime, where the fluctuations are distributed more evenly on the surface. The first stabilization mechanism is due to the averaging between regions of strong and weaker turbulence and, relative to the threshold of the most turbulent field line, additionally involves a newly observed upshift. The transport difference stemming from the surface and the most unstable location might serve as an additional figure-of-merit for the characterization of turbulence optimization in a stellarator [24]. In addition, due to the very small electron gyroradius, it is highly likely that turbulence caused by the electron temperature gradient will also experience stabilization in the small-ρ * regime. The other type of stabilization occurs due to the suppression of large energy-containing eddies by the fast variation of the equilibrium magnetic field, and leads to a reduction of the transport stiffness. Because underlying mechanisms are quite generic, both stabilization mechanisms are expected to act also in other types of stellarators, where the "poloidal" localization of turbulence is either theoretically predicted [15] or experimentally verified [25]. Finally, we note that 3D geometrical turbulence stabilization could also be introduced in a tokamak if magnetic perturbations break the axisymmetry, thereby enhancing the local shear [26], and confining turbulence structures in a stellarator-like fashion.
The GENE simulations were performed on the Helios (Japan) and Hydra (Germany) supercomputers. We appreciate technical support by J. Geiger and Y. Turkin with the magnetic equilibria.