Stellarator turbulence at electron gyroradius scales

Electromagnetic gyrokinetic simulations of electron-temperature-gradient-driven modes on electron gyroradius scales are performed in the geometry of an advanced stellarator fusion experiment, Wendelstein 7-AS. Based on linear simulations, a critical electron-temperature-gradient formula is established which happens to agree quite well with a previously derived formula for tokamaks in the appropriate limit. Nonlinear simulations are used to study the turbulence and transport characteristics which are dominated by the presence of high-amplitude radially elongated vortices or `streamers'. The role of Debye shielding effects is also examined.


Introduction
It has long been known that the cross-field particle and heat transport in magnetically confined fusion plasmas generally exceeds predictions based on collision-induced processes by up to two orders of magnitude. There is strong experimental evidence that these anomalous losses can generally be attributed to various kinds of fine-scale turbulence (at perpendicular spatial scales somewhat larger than the ion gyroradius) driven by density and/or temperature gradients [1]. For several reasons it is useful to distinguish between edge and core turbulence where the edge region is mainly characterized by much shorter length scales of the density and temperature profiles, L edge ⊥ ∼ 1 cm versus L core ⊥ ∼ 1 m. First and foremost, the microinstabilities underlying the turbulent transport processes differ in both regions. Whereas edge turbulence seems to be mainly driven by electromagnetic ion-temperature-gradient (ITG) modes, drift Alfvén waves (electromagnetic electron drift waves), and ballooning modes [2,3], core turbulence is believed to be caused by electrostatic ITG modes and trapped electron modes [4]. Several experimental studies indicate, however, that this 'standard model' of turbulence in fusion plasmas is incomplete. For example, high-performance plasma discharges with internal transport barriers 35.2 exhibit fairly high electron heat transport (way beyond the neoclassical level) despite the fact that all fine-scale turbulence is suppressed by strong poloidal shear flows [5]. This discrepancy might be due to plasma turbulence driven by electron-temperature-gradient (ETG) modes on hyperfine (i.e., electron gyroradius) scales.
Assuming a perfect isomorphism between ETG and ITG turbulence, χ ETG e is smaller than χ ITG i by the square root of the ion-to-electron mass ratio if the temperatures and normalized temperature gradients of the two species are comparable. However, nonlinear gyrokinetic simulations of tokamak plasmas have shown recently that ETG turbulence can yield experimentally relevant electron heat flux levels since even in the electrostatic and adiabatic limit there are important differences in nonlinear dynamics between ETG and ITG modes which break this isomorphism [6,7]. If this surprising finding also holds for stellarators, it might have implications for the performance of next-generation experiments with minimized neoclassical transport like Wendelstein 7-X [8]. In the absence of generic model magnetohydrodynamic equilibria for stellarators, we will focus specifically on Wendelstein 7-AS [9] which can be seen as a prototype for a certain important class of magnetic fusion devices called 'helical advanced stellarators' [10]. A typical flux surface shape for Wendelstein 7-AS is shown in figure 1 together with part of the confining magnetic field coil system. The rest of this paper is organized as follows. In section 2, we derive a compact algebraic formula for the critical ETG of ETG modes which is compared to a corresponding tokamak formula. In section 3, we present nonlinear gyrokinetic simulations of ETG turbulence and transport, highlighting the role of streamers and examining Debye shielding effects. We conclude with a summary in section 4.

Nominal parameters
Our nominal physical parameters are based on the Wendelstein 7-AS discharge no 49689 at an effective minor radius of 6 cm, i.e., at about 1/3 of the effective minor radius of the last closed flux surface. The electron density and temperature are given by n e = 2.2 × 10 19 m −3 and T e = 2.3 keV, respectively; their normalized gradients are R/L ne = −(R/n e )(dn e /dr) = 0 and R/L Te = −(R/T e )(dT e /dr) = 10 where R = 2.0 m is the major plasma radius. The 35.3 magnetic field strength is B = 2.5 T, the rotational transform is -ι = 1/q ≈ 1/3, and the global magnetic shear isŝ = −(r/ι)(d-ι/dr) ≈ 0. The nominal value for τ = Z eff (T e /T i ) was chosen to be unity although the experimental value for this particular discharge was of the order of 10. The reason for this choice will become clear in section 3.2. Moreover, the ratio of the electron Debye length to the electron gyroradius is λ De /ρ e = 1.7 and the electron plasma beta is given by β e = 8πn e T e /B 2 = 3.2 × 10 −3 .

Brief characterization of ETG modes
In the electrostatic and adiabatic limit, the linear dynamics of ETG modes (see, e.g., [11]- [13]) bears strong resemblance to that of ITG modes, with the roles of electrons and ions reversed. At typical perpendicular wavenumbers k ⊥ ρ e < ∼ 1 k ⊥ ρ i and frequencies ω ∼ v te /L Te , toroidal ETG modes are linearly unstable when R/L Te exceeds a critical value. Here, ρ e and ρ i are, respectively, the gyroradii of thermal electrons and ions, and v te is the electron thermal speed. From the gyrokinetic Poisson equation it can be shown that the response of the ions is adiabatic on sub-ion-gyroradius scales, i.e.,ñ i /n i = −Z i (T e /T i )(eΦ/T e ). In the presence of additional impurity ion species, Z i can simply be replaced by Z eff -assuming that the temperatures of the main ions and impurities are comparable-and enters only via the combined quantity τ = Z eff (T e /T i ). Because of ion adiabaticity, ETG modes do not cause significant particle or ion heat transport.

Some linear simulation details
In order to study the dependence of the critical value of R/L Te on various plasma parameters we employ a linear electromagnetic gyrokinetic code, ALEGC, which is described in [14]. The system of gyrokinetic Vlasov-Maxwell equations [15,16] for electrons and adiabatic ions is solved with an explicit numerical scheme in a general three-dimensional magnetic field geometry. Since the latter is given and many plasma parameters (such as β e , the mass ratio m i /m e , and the electron-ion collision frequency ν ei ) are found to have a negligible effect on (R/L Te ) crit , the only relevant quantities left to vary are τ and R/L ne . The role of Debye shielding effects will be discussed in section 3.4. For now, we set the electron Debye length to zero.
The method used to derive critical gradients for any given set of physical parameters is the same as in [17]: we choose k x = 0 (i.e., the ballooning parameter θ 0 = 0) and a set of k y s (poloidal wavenumbers) in the region of the fastest-growing linear ETG modes (for flat density profiles; the growth rate spectra reach their maximum values around k y ρ e ∼ 0.3). Then we vary R/L Te and find the linear growth rates as a function of R/L Te for each value of k y . Linear extrapolation and subsequent minimization over all k y s yields the final result, (R/L Te ) crit . The two-dimensional parameter space was covered by a grid with τ ∈ (0.5, 1, 1.5, 2, 4) and (R/L n ) ∈ (1, 2.5, 5). For each set of parameters, we used 3-7 values of k y to find the growth rate spectrum, and evaluated this spectrum for 3-5 values of R/L Te . In total, approximately 250 linear ALEGC runs were carried out. The numerical parameters (particularly the time step, the box size, and the number of grid points in the parallel direction) were checked regularly to ensure convergence.
A typical linear mode structure (real and imaginary part of the perturbed electrostatic potential) of an ETG mode in Wendelstein 7-AS geometry (using the nominal parameters and choosing k y ρ e = 0.3) is shown in figure 2. For comparison, the result for an equivalent tokamak is also shown, with the area below the curve kept constant. The field-line-following coordinate z (basically equivalent to the toroidal angle) has its origin at the outboard side of the elliptical plane since there the magnetic curvature driving toroidal ETG modes reaches its maximum value. Accordingly, the parallel mode structure exhibits a rather pronounced ballooning structure, i.e., a localization around z = 0, with some helical perturbations.

Critical gradient formula
Compared to the axisymmetric case of a tokamak, the magnetic field structure of a stellarator is complicated by helical deformations of the flux surface geometry. Local variations of magnetic field strength, field line curvature, and shear alter the character of the geometric coefficients in the gyrokinetic equation significantly from simple toroidicity. Moreover, in contrast to the case for tokamaks, the global magnetic shear is chosen to be close to zero across the entire plasma domain in order to prevent the formation of magnetic islands at low-order rational values of the rotational transform -ι. Nevertheless, local deviations of radially neighbouring field lines occur due to the deformations of the flux surface varying in a helical period of Wendelstein 7-AS between triangular and elliptical. These specific properties do not, in general, allow one to apply results found in tokamak geometry to a stellarator. It is thus quite surprising that the analysis of our linear simulation data leads to a critical gradient formula which agrees quite well with a previously derived formula for tokamaks [17]. The linear threshold obtained in this work (neglecting plasma shaping effects which turn out to have only minor influence) is given by with In Wendelstein 7-AS, the global magnetic shearŝ is negligibly small, but the local magnetic shear varies along the field line. Making an ansatz analogous to that in the tokamak case and 35.5 * Figure 3. Linear threshold (R/L Te ) * crit of the ETG instability according to the regression formula, equations (3) and (4), as a function of the linear thresholds (R/L Te ) GK crit obtained directly from linear gyrokinetic simulations.
setting B = 0, the data from our parameter studies are used to evaluate A and C. A least-squares fit yields the following critical gradient formula for ETG modes: with A = 1.23 ± 0.1 and The quality of this regression analysis is illustrated in figure 3.

Discussion
Remarkably, our result is in good agreement with the tokamak formula, equations (1) and (2), a result which is consistent with recent work by Jost et al [18] who performed global linear gyrokinetic simulations of ITG modes with adiabatic electrons in the core of quasi-symmetric stellarators, showing that ITG modes are not significantly affected by local magnetic shear. Since in the electrostatic, adiabatic limit, linear ITG and ETG modes are analogous (apart from global geometry and finite-λ De effects), one can indeed expect the two studies to yield qualitatively similar results.
On the other hand, our finding is in contrast to the results from simulations of drift waves driven by cold-ion density gradients for plasma edge parameters, where it was found that local magnetic shear in a stellarator is an effective damping mechanism when the global shear is low [19,20]. This difference can be understood by remembering that drift waves are dissipative modes, depending crucially on the parallel electron dynamics which in turn reacts to changes in the magnetic geometry, whereas ITG and ETG modes are reactive and of basically adiabatic nature. Local shear variations should thus have less impact on core than on edge instabilities which is exactly what we observe.

35.6
A possible explanation for the (small) deviation of the Wendelstein 7-AS critical gradient formula from the one for tokamaks is rather likely to be found in the influence of (elliptic and triangular) flux surface deformation, in analogy to the plasma shaping effects observed in the tokamak formula [17]. One should also note, however, that the tokamak formula has been derived forŝ ≥ 0.2, and therefore needs to be extrapolated toŝ = 0, which may lead to additional changes.

Nonlinear flux-tube simulations in stellarator geometry
The fact that we focus on hyperfine perpendicular scales in a configuration with an edge rotational transform -ι a ≈ 1/3 and nearly zero global magnetic shear means that it is justified to employ the well-established method of flux-tube simulations, using periodic boundary conditions in the field-line-following coordinate after exactly three toroidal turns (assuming -ι = 1/q = 1/3) and neglecting perpendicular variations of the geometric coefficients. This is a serious simplification in view of the need for global simulations of ion-scale turbulence in arbitrary stellarator geometries [18].
Due to the broken toroidal symmetry in a stellarator, it is in general not sufficient to represent an entire flux surface by just one flux tube like in a tokamak. However, the fivefold discrete symmetry together with the rotational transform of one third ensure that a fair proportion of the flux surface is already covered by just one representative flux tube, intersecting the outboard side of the elliptical plane. The situation is visualized in figure 4   the three-dimensional magnetic field structure are generated by means of the Gourdon code [21] as a Fourier expansion in straight-field-line coordinates.

Some nonlinear simulation details
To study ETG turbulence in Wendelstein 7-AS geometry, we employ a massively parallel code, GENE, which solves the nonlinear electromagnetic gyrokinetic Vlasov-Maxwell equations [15,16] on a fixed grid in five-dimensional phase space. Numerical details concerning GENE can be found in [17]. The simulation is carried out in a magnetic flux tube with perpendicular dimensions of L x = L y = 128ρ e where ρ e is the gyroradius of thermal electrons. The ions are taken to respond adiabatically, an approximation which is well justified on sub-ion-gyroradius scales. Trapped electron effects are neglected for similar reasons. We use 64 grid points in both perpendicular directions as well as in the field-line-following coordinate to resolve the strong parallel variations of the geometric coefficients as illustrated in figure 2. 30 × 10 points are employed in velocity space.
It turns out that in addition to that, a long simulation time is necessary to reach a quasi-static dynamical equilibrium state of the turbulent system. The simulation shown below took 150 h on 128 processors of a Cray T3E-600. In view of this immense numerical effort we refrained from effecting systematic parameter variations around the nominal values. In accordance with the linear runs presented above and corresponding studies in tokamak geometry [22], we expect a smooth rise of χ ETG e with R/L Te increasing from its critical value, (R/L Te ) crit ≈ 2.5, to its nominal value, R/L Te = 10.
At this point, a few remarks concerning our choice of the nominal value of τ seem to be in place. The aim of the present work is not to simulate and explain previous experiments, but rather to explore scenarios for Wendelstein 7-AS in which ETG modes might be detected in future experiments. Choosing a large value of the parameter τ would, according to equations (3) and (4), raise the critical temperature gradient substantially, making it more likely that other electron thermal transport mechanisms dominate. Moreover, the electron heat transport would be reduced (holding the distance of R/L Te from its critical value fixed) as has been found in corresponding tokamak simulations. For example, for the plasma parameters given in [6], using τ = 4 instead of τ = 1 leads to a reduction of χ e by more than a factor of 4. Therefore, one expects that for discharges like No 49689 with dominant electron heating and τ ∼ 10, it would be hard if not impossible to directly observe ETG modes experimentally. Therefore we based our simulation parameters on previously measured electron temperature profiles, but used a higher ion temperature and thus a lower value of τ , namely τ = 1. Experimentally, this can be achieved by adding some direct ion heating (e.g., by neutral beam injection). In addition, this choice reflects better the typical conditions expected in a possible stellarator reactor.

Nonlinear simulation results for λ De /ρ e = 0
As has been shown recently by nonlinear gyrokinetic simulations, ETG turbulence in tokamak plasmas may exhibit high-amplitude 'streamers', i.e., radially elongated vortices, which are able to greatly enhance the turbulent radial transport [6,7]. To achieve experimentally relevant levels of the electron heat conductivity, χ e , this feature is vital. That is because the intrinsically set mixing length level for χ ETG e is smaller than that for ITG turbulence by the square root of the electron-to-ion mass ratio (assuming similar temperatures and normalized temperature gradients for both species). With the latter being typically of the order of the experimentally observed 1 m 2 s −1 , this indicates that χ ETG e without streamers is not expected to contribute significantly 35.8 to the turbulent electron heat transport under normal tokamak conditions. In the presence of streamers, however, the hyperfine-scale ETG modes might not be negligible. Here, we present the first simulations of ETG turbulence in stellarator geometry investigating this issue.
Time traces of the electrostatic and electromagnetic components of χ e as defined, e.g., in [23], are shown in figure 5. Here, we expressed both contributions in MKS units, using the nominal parameters given above. Following an initial exponential growth phase which is dominated by the fastest-growing linear ETG mode, the simulation starts to saturate around t ≈ 300L Te /v te . The initial saturation level is maintained until t ≈ 600L Te /v te . Then, the electrostatic component of χ ETG e (due to perpendicular E ×B advection) more than doubles, and the electromagnetic component (due to fast electron motion along perturbed magnetic field lines) no longer fluctuates around zero. The spatially averaged fluctuation amplitude ofÃ increases by a factor of four. Obviously, this slow relaxation to a turbulent equilibrium state involves nonlinear electromagnetic effects. (Note that linear electromagnetic effects are negligible for ETG modes as was mentioned above.) The transport is clearly dominated by its electrostatic part, with the electromagnetic component contributing only some 2% to the total χ ETG e ∼ 2.5 m 2 s −1 . This result is in accordance with [24] where it was found that, despite possible electromagnetic effects on various kinds of tokamak turbulence, the induced transport is predominantly electrostatic in nature except close to the ideal ballooning limit.
The fairly high transport level is indicative of streamer activity. Indeed, this is confirmed by a snapshot ofφ,ñ e , andT e at t = 900L Te /v te presented in figure 4. Hence we have shown that the presence of high-amplitude streamers with its associated transport enhancement is not dependent on a particular kind of magnetic geometry (i.e., that of a tokamak). In physical units, n e /n e0 ∼ 0.2% andT e /T e0 ∼ 0.25% in this case, which significantly exceed the mixing length expectations,ñ e /n e0 ∼T e /T e0 ∼ ρ e /L Te ∼ 0.023%. An additional transport boost is provided by the finite vortex aspect ratio as can be understood, e.g., in terms of a simple random walk argument: The transport level can be estimated using the formula χ ∼ ∆x 2 /∆t where ∆x and ∆t are typical radial (spatial) and temporal scales of the turbulence. From this it is clear that a radial elongation larger than unity increases the turbulent transport (holding the poloidal reference scale fixed). In figure 6, the radial autocorrelation function of φ is shown as a function of ∆x in units of 35.9  ρ e . Its FWHM value turns out to be 21ρ e which is about three times larger than the value obtained in a typical (tokamak) case without streamers which is shown for comparison. The average vortex aspect ratio was computed to be 2.0 (0.9) in the simulation with (without) streamers.
We also note in passing that the radial correlation length of 21ρ e is <1/6 of the radial width of the simulation domain. Under these circumstances there is no influence of the nature of the boundary conditions on the results. In particular, the existence and nature of the streamers is not affected. Runs with varying box sizes (up to 512ρ e ) in tokamak geometry (presented in [22]) confirm that. Although Dirichlet boundary conditions (zeroing out the fluctuations at the radial boundaries) are also an implemented option in the GENE code, they are rarely used since they lead to problems like quasi-linear flattening of the background profiles, gradually shutting off the turbulence. Periodic radial boundary conditions are therefore both justified and advantageous. In figure 7, the time-averaged k y spectra [25] of two quantities are presented:ñ 2 e and T 2 e , denoted, respectively, by n and t. Both the electron density and temperature spectra have pronounced peaks around k y ρ e ∼ 0.2 (despite the fact that the linearly most unstable modes are at k y ρ e ∼ 0.4) and fall off with a near power law for k y ρ e > ∼ 0.3, reflecting the self-similarity of the turbulent ETG system at small scales. magnetic field in teslas and n e19 is the electron density in units of 10 19 m −3 ). In general, we have λ De /ρ e ∼ 1 for typical core plasma conditions. The role of finite-λ De effects can be understood in the following way. λ De enters the gyrokinetic equations via the modified Poisson equation which for k ⊥ ρ e 1 k ⊥ ρ i can be written as Increasing k ⊥ λ De is thus expected to have a similar effect to increasing τ . In both cases, for a given density perturbation the potential drops, weakening the instability. This is also reflected by the critical gradient formula, equations (3) and (4), which states that (R/L Te ) crit increases with increasing τ . Obviously, Debye shielding acts predominantly on the scales which are smaller than λ De . This means that the high-k y part of the linear growth rate spectrum is stabilized first as can be seen in figure 8. For τ > 1, we find that finite-λ De effects are reduced, in accordance with equation (5).
As we have shown in figure 7, the turbulent k y -spectra tend to peak around k y ρ e ∼ 0.15−0.2. Therefore it is reasonable to expect Debye shielding to affect the turbulence only in as much it affects those long-wavelength ETG modes. Going from λ De /ρ e = 0 to λ De /ρ e = 1.7, the linear growth rates in that k y -range go down by about a factor of 2-3 (see figures 8 and 9). As can be inferred from figure 10, this change is mirrored by the nonlinear simulation. The transport is still dominated by electrostatic streamers, with the electromagnetic contribution being negligible. The k y -spectra are found to be very similar to the ones displayed in figure 7, except that they tend to peak at slightly lower values of k y .

Summary
In conclusion, we have examined some main characteristics of toroidal ETG modes in the geometry of the advanced stellarator fusion experiment Wendelstein 7-AS by means of gyrokinetic simulations. We established a critical ETG formula which agrees quite well with a previously derived formula for tokamaks in the appropriate limit, indicating that toroidal ETG modes are not significantly affected by local magnetic shear. The single most important plasma parameter determining (R/L Te ) crit is the modified temperature ratio τ = Z eff (T e /T i ).

35.12
Our nonlinear gyrokinetic simulation for the nominal parameters and λ De /ρ e = 0 yielded an electron heat flux level of χ ETG e ∼ 2.5 m 2 s −1 . The transport was clearly dominated by its electrostatic part (due to perpendicular E × B advection), with the electromagnetic component (due to fast electron motion along perturbed magnetic field lines) contributing only some 2% to the total χ e . This surprisingly high value of χ ETG e is caused by radially elongated turbulent eddies with large saturation amplitudes. Hence, for the first time, it is shown that the occurrence of high-amplitude streamers is not dependent on a particular kind of magnetic geometry (i.e., that of a tokamak). Fluctuation spectra and levels were also measured, together with the average streamer aspect ratio.
Debye shielding mainly affects the part of the linear growth rate spectrum which is characterized by k ⊥ λ De > ∼ 1. For our nominal parameters (λ De /ρ e = 1.7), Debye shielding is found to reduce the relevant linear growth rates as well as the turbulent transport level by about a factor of two, with the long-poloidal-wavelength streamer activity remaining dominant.