Global gyrokinetic simulations of ITG turbulence in the configuration space of the Wendelstein 7-X stellarator

We study the effect of turbulent transport in different magnetic configurations of the Weldenstein 7-X stellarator. In particular, we performed direct numerical simulations with the global gyrokinetic code GENE-3D, modeling the behavior of Ion Temperature Gradient turbulence in the Standard, High-Mirror, and Low-Mirror configurations of W7-X. We found that the Low-Mirror configuration produces more transport than both the High-Mirror and the Standard configurations. By comparison with radially local simulations, we have demonstrated the importance of performing global nonlinear simulations to predict the turbulent fluxes quantitatively.


I. INTRODUCTION
Turbulent transport is one of the major uncertainties in the design of any magnetic fusion device. This is particularly true for stellarators, given the many degrees of freedoms available for the shape of the magnetic geometry [1,2]. Indeed, with the recent progress in the optimization for neoclassical transport in stellarators, as demonstrated in the planning, construction, and operation of Wendelstein 7-X (W7-X) [3,4], it has become clear that turbulent transport is the main mechanism limiting the confinement time [5]. Therefore, in order to further improve stellarators as a viable long-term alternative to tokamaks for future power plants [6][7][8], we need to understand, predict, and control turbulent transport.
In this context, the present paper aimed to understand how Ion Temperature Gradient (ITG) turbulence [9], which appears to be one of main micro-instabilities limiting the performance of W7-X in standard scenarios [10], is affected by different magnetic configurations which are experimentally available. This knowledge might provide guidance to prepare and improve future campaigns of W7-X.
We used the gyrokinetic simulation framework [11,12] from which much progress has been made towards understanding turbulence transport in stellarators [13] and performed one of the first global gyrokinetic turbulence studies of W7-X. In particular, we used the gyrokinetic code GENE-3D [14], a newly developed global version of the well established gyrokinetic code GENE [15] for stellarator geometries, to perform linear and nonlinear simulations of W7-X for different magnetic configurations.
The paper is organized as follows. In Section II, we briefly describe the gyrokinetic code GENE-3D. In Section III, the different W7-X magnetic configurations and equilibrium profiles used in this work are introduced. We then present linear simulation results in Section IV, followed by nonlinear turbulent results in Section V. A comparison with radially local simulations is performed in Section VI, demonstrating the importance of global nonlinear simulations to predict turbulent a) Electronic mail: alejandro.banon.navarro@ipp.mpg.de fluxes quantitatively. Finally, conclusions are given in Section VII.

II. THE GLOBAL GYROKINETIC STELLARATOR CODE GENE-3D
The results presented in this paper are produced with GENE-3D, the global nonlinear version of the GENE code which supports stellarator geometries. GENE-3D has been benchmarked linearly against EUTERPE [16] in W7-X geometry and nonlinearly against the global version of GENE [17] in tokamak geometry. In addition, GENE-3D has been recently applied to study gyrokinetic turbulence in configurations derived by novel optimization techniques in W7-X [18]. In the following section, we outline the main features of the code, but for a detailed description, we refer the reader to Ref. [14].
GENE-3D solves the gyrokinetic Vlasov equation coupled self-consistently to Maxwell's equations on a fixed grid in five-dimensional phase space (plus time), consisting of two velocity coordinates v (velocity parallel to the magnetic field), and µ (magnetic moment), and three magnetic fieldaligned coordinates x, y and z, defined as The radial coordinate (x) is based on the normalized toroidal flux ρ tor = Φ tor /Φ edge , where Φ tor is the toroidal flux and Φ edge its value at the last closed flux surface. The bi-normal coordinate (y) selects a field line α = qθ − ϕ on a given flux-surface. Here, q is the safety factor, θ is the poloidal PEST angle [19] and ϕ is the geometrical toroidal angle. In the definition of y, the constant C y = x 0 /|q 0 |, with q 0 being the safety factor at a reference position x 0 , is used to have y as a length rather than an angle-like coordinate. Finally, the parallel coordinate (z) denotes a position along the magnetic field line. The sign of the poloidal magnetic field σ Bp is introduced so that the unit vector in the parallel direction is always in the direction of the magnetic field.
In these coordinates, the equilibrium magnetic field can be expressed as with C(x) = (xB 0 )/(|q(x)|C y ), where B 0 denotes the magnetic field strength on the magnetic axis. The magnetic field is a solution of a three-dimensional ideal MHD equilibrium with nested flux surfaces, computed by the Galerkin Variational Equilibrium Code GVEC. In addition to the field and the q(x) profile, GVEC supplies, through an interface, the geometrical information necessary to map the equilibrium to the field-aligned grid. It also provides the mapping into Cartesian coordinates to visualize the simulation data produced by GENE-3D.
To solve the gyrokinetic Vlasov equation numerically, GENE-3D uses the δf approach [12]. In this approach, the gyrocenter distribution function F σ of species σ is split into a background F 0σ and a first order perturbation part F 1σ : where the background distribution function is assumed to be a local Maxwellian in GENE-3D.
Keeping only first-order terms in the perturbed distribution function, the resulting electrostatic gyrokinetic Vlasov equation solved in GENE-3D reads: with The last term on the left-hand side in equation (6) couples neoclassical and turbulence transport, and it may affect the longterm evolution of the system in the presence of collisions [20]. A linearized Landau-Boltzmann collision operator is currently implemented for C[F 1σ ]. In the present paper, neoclassical contributions and collisions are neglected for simplicity. The drift velocities in equation (6) are defined as: Here, Ω σ = (q σ B 0 )(m σ c) is the gyrofrequency of a species σ with charge q σ and mass m σ ,b 0 = B 0 /B 0 and c is the speed of light. Furthermore, φ 0 is the equilibrium electrostatic potential which can be employed to consider externally imposed (long-wavelength) radial electric field effects [21,22] andφ 1 is the gyroaveraged perturbed electrostatic potential, defined asφ where X is the gyrocenter position and the gyroradius vector r(α) is orthogonal to the local magnetic field. The perturbed electrostatic potential is calculated selfconsistently from Poisson's equation, which, when written in terms of the perturbed particle density n 1σ , reads: where the left-hand side is neglected in GENE-3D (quasineutral limit). Finally, in this work, we treat the electrons as an adiabatic species. In this approximation, the adiabaticity relation on a flux-surface reads: with e being the (positive) elementary charge and · FS denoting a flux surface average [23] and V being the volume enclosed by that flux surface. Equation (6) is solved numerically by discretizing the distribution function on the aforementioned five-dimensional grid. This allows one to reduce the original hyperbolic integro-differential system of equations to a system of ordinary differential equations, which are then explicitly integrated in time. GENE-3D currently uses a 4th order explicit Runge-Kutta (RK4) integrator; the timestep is computed at the initialization and maximized during a run to ensure optimal stability [24]. Because of spatial dependencies in all three spatial coordinates, fourth-order centered finite difference schemes are used to calculate derivatives. This is in contrast to other GENE versions that can rely on spectral (Fourier) methods in certain directions. A zero Dirichlet boundary condition is used in the the radial direction while periodic boundary conditions are applied in the bi-normal direction. The conventional twist-and-shift [25] boundary condition is applied in the parallel direction after having followed a field line for one poloidal turn. The velocity space is discretized using a regular equidistant grid for v (assuming as well a zero Dirichlet boundary condition and employing fourth-order centered finite differences) while a Gauss quadrature scheme, with Gauss-Legendre weights and knots, is used in the µ direction.
The Arakawa scheme [26] is used to evaluate the nonlinear term (the sixth term on the left-hand side of equation (6)) in order to ensure free energy conservation [27]. Finally, the gyroaverage operator, which does not vary during a simulation since it depends only on the background temperature profile and magnetic geometry, is discretized using finite element in the direction perpendicular to the magnetic field using bicubic piecewise polynomials as basis. This allows us to write Poisson's equation as a system of linear equations that is solved using the PETSCs library [28][29][30].

III. W7-X MAGNETIC CONFIGURATION SPACE AND EQUILIBRIUM PROFILES
Wendelstein 7-X is the first large, superconducting machine of the HELIAS (HELIcal axis Advanced Stellarator) [31] type, optimized for strongly reduced neoclassical transport and low bootstrap current. W7-X has a major radius of R 5.5 m, an aspect ratio of A = R/a 10 and 5 field periods. It has a very flexible coil current system composed of 5 modular and 2 planar (tilted) superconducting magnetic coils in each half period. This coil system allows a large variety of magnetic configurations [32,33].
In this work, we investigate the influence of different magnetic field geometries on Ion Temperature Gradient (ITG) turbulence. For this purpose, we selected three different magnetic configurations: the Standard configuration (EIM) and two configurations with different magnetic mirror ratios with respect to the Standard configuration; the High-Mirror (KJM), and the Low-Mirror (AIM) configurations. The magnetic equilibria were calculated with the VMEC code [34,35]. The change in the magnetic mirror ratio slightly changes the rotational transform (see figure 1).
To isolate the effect of the magnetic configuration on ITG turbulence, we used the same equilibrium profiles in all previously described magnetic geometries. We also assumed equal ion and electron temperatures. The profiles and their respective normalized logarithmic gradients, with the latter defined as a/L X = −a d ln X where X is the ion temperature (T i ) or the electron density (n e ), are shown in figure 2. For simplic- Low-Mirror (AIM). The High-Mirror configuration has a magnetic mirror ratio of 6%, while is reduced to 3.7% for the lower-mirror configuration. ity, the profiles were not taken from an actual W7-X discharge, but were produced with the DKES code [36,37] assuming an idealized W7-X scenario characterized by a broad range of radial positions in which the ion logarithmic temperature gradient is large and the ratio η i = L ni /L Ti is greater than one, as shown in figure 2. In such a scenario, we expect ITG modes to be strongly driven [9], and trapped electron physics, which is neglected in the adiabatic electron response employed in the simulations, to be less important in capturing the dependence of ITG turbulence on the different magnetic configurations correctly [38].

IV. LINEAR GLOBAL ITG SIMULATIONS
Before performing turbulence simulations, we first want to characterize the linear modes present in the different W7-X magnetic configurations. Therefore, in this section, we compare the growth rate, the toroidal mode numbers, and the mode structure for the different cases considered in this paper. The numerical setup for linear simulations includes the following: neglect of the nonlinear term so that only ITG modes are captured; a radial domain from ρ tor ∈ [0.1, 0.9] with zero Dirichlet boundary condition; a buffer zone of 5% of the radial domain with a Krook damping operator with a coefficient of 1.0 v i /a significantly larger than the expected growth rate of the instability at each side of the radial domain to avoid nu- We display in figure 3 the growth rates (γ) and peak toroidal mode numbers (n is the mode number corresponding to the peak in the Fourier spectrum of φ 1 ) for the different magnetic configurations. We observe that the growth rates and the peak toroidal mode numbers are very similar for all the different W7-X configurations. In particular, the ITG mode is characterized by a maximum growth rate between 0.14 − 0.16 v i /a and a peak toroidal mode number around n ≈ 180.
We also found that the mode structure is very similar for all configurations. We show this in figure 4 by comparing the square amplitude of the electrostatic potential φ 2 1 vs. the radial coordinate (top left) and vs. the PEST poloidal angle (top right) -averaged over the remaining coordinates. The mode is localized around the outboard mid-plane (θ ≈ −0.5) and peaks at ρ tor ≈ 0.75. This radial location is at a position where a/L Ti is large, and η i = L ni /L Ti is at a local maximum (see figure 2), thus, favoring strongly unstable ITG modes. Finally, the three-dimensional mode structure in this flux-surface is also shown in figure 4 (bottom) but only for the Standard EIM configuration to avoid redundancy. We observe that the mode is highly localized in the flux-surface, which agrees with previously published linear results of W7-X by EUTERPE [16,39,40]. To summarize this section, we found that ITG modes modeled with an adiabatic electron response do not vary significantly for the different W7-X configurations studied in the present paper. They are characterized by large toroidal mode numbers (n ≈ 180), and their mode structure is highly localized both radially and in the flux-surface.

V. NONLINEAR GLOBAL ITG SIMULATIONS
Having characterized linear ITG modes for all configurations, we proceed to the study of nonlinear simulations of ITG turbulence. There are two ways of performing global nonlinear simulations and reaching a quasi-steady state: the flux-driven approach, in which fixed sources (e.g., of energy) are used and cause the profiles to slowly evolve toward a quasi-steady state; and the gradient-driven approach, which employs sources and sinks to keep the plasma profiles close to the initial ones, so that the turbulence drive, characterized by the pressure gradients, is maintained constant during the simulation. The flux-driven setup can be considered closer to the experimental situation, but it requires longer simulations. The gradient-driven approach needs shorter runs to reach a quasi-steady state and is, therefore, less computational demanding. In this work, we adopted the gradient-driven approach. The simulation setup used is as follows: the resolution used is the same as in the linear case, {x, y, z, v , µ} = {192, 256, 128, 48, 12} -we have performed convergence studies to ensure the validity of this choice; a Krook-type heat source with a relaxation rate of 0.02 v i /a was used to main- where F pc 1i is the perturbed part of the ion particle distribution function, and A is the flux-surface area. In this figure, the shaded region corresponds to an estimate of the standard deviation of the set of means of consecutive temporal subdomains of the saturated state; and the horizontal dashed lines show their radial average over the whole region. In particular, we observe that the Low-Mirror configuration produces approximately 1.7 times more radially averaged heat flow than the High-Mirror configuration and 1.4 times higher than the Standard configuration. The difference in the heat fluxes is mainly localized in the region around ρ tor = 0.65 − 0.75. In this region, the heat flow can reach values around 30 MW for the Low-Mirror configuration and around 20 MW for the other two configurations. These values are larger than the currently available heating power of W7-X. Although a quantitative estimate of the flux amplitudes would require a more comprehensive plasma description, e.g. kinetic electron physics, and/or electromagnetic fluctuations, these results may suggest that the profiles used in this work could not be achieved experimentally in W7-X.
Note that for the Low-Mirror configuration to produce more heat flux, while maintaining the profiles, more heat must be injected by the Krook-type heat source term. This is indeed the case as shown in figure 6, where we plot the profiles of the logarithmic ion temperature gradient (top) and the heat source term (bottom). We observe that the ion gradients are well maintained during the simulation for all cases. However, the net heat source input, given by the sum of the positive contributions of the heat source term, is indeed higher for the Low-Mirror configuration. Although the amplitudes are different between the configurations, the spatial structure of the heat flux is very similar between them. As we show in figure 7, the spatial structure of the heat flux at ρ tor = 0.7 is maximal around θ ≈ 0 and at α ≈ 0 (and 2π/5 due to periodicity) for all cases. We then compare this structure with that of the bi-normal curvature, defined as K y = − (B 0 × ∇B 0 ) ·ŷ/B 2 0 . As expected by the fact that ITG modes are driven unstable by negative values of bi-normal curvature, we observe that the heat flux is at a maximum in the regions of bad curvature, i.e., where K y is at a minimum. Taking into account that at this position the bi-normal curvature has a very similar structure (see figure 7(bottom)) for all the magnetic configurations, the differences in transport must come from a different geometrical term.
Furthermore, in order to compare with the linear case, we display in figure 8 the three-dimensional representation of the time average heat flux, but only for the Standard configura- tion. We observe that most transport is localized in a narrow region at the outboard mid-plane. This localization is, however, broader than the one presented by the linear mode structure (see figure 4). This difference can be explained by analyzing the ion heat flux spectra at ρ tor = 0.7. As shown in figure 9, the heat flux spectra are dominated by a broad range of toroidal modes with a maximum around n ≈ 70. This maximum is therefore a larger scale mode than the most dominant mode in linear simulations (n ≈ 180). In fact, the most unstable mode in linear simulations does not contribute significantly to transport. This indicates that it is crucial to perform nonlinear simulations to characterize the system correctly.
Finally, we compare the radial profile of E ×B shear flows, defined as: since they are one of the important phenomena responsible for the saturation of turbulence [41,42]. We observe that the radial profiles of the shear flows (figure 10) are practically identical (and smaller than the linear growth rate) for all the configurations, despite the higher heat flux in the Low-Mirror configuration. Therefore, this result suggests that shear flows To summarize the section, we showed that ITG turbulent transport depends on the magnetic configuration used in W7-X. This is contrary to what we observed in linear simulations. In particular, we found that the Low-Mirror configuration produces more heat flux than both the Standard and the High-Mirror configurations. This happens despite similar levels of E × B shear flows, implying that the shear flows are less important for the saturation of turbulence in the Low-Mirror configuration. Nevertheless, we found that the spatial structure of the heat flux is very similar for all configurations, being localized in the regions of bad curvature, although it is broader than in the linear case. The discrepancy between these results may be explained by the fact that the heat flux spectra are dominated by larger scale modes than in the corresponding linear simulations.

VI. COMPARISON WITH RADIALLY LOCAL SIMULATIONS
We have performed global simulations to study the effect of ITG turbulence on different magnetic geometries of W7-X and found that the Low-Mirror configuration produces more heat flux than both the Standard and the High-Mirror configurations. In this section, we want to investigate if a similar conclusion can be reached with a reduced -and less computationally demanding -model, such as a radially local one [43,44].
For this purpose, we performed radially local simulations with GENE-3D. These simulations are conducted by selecting a flux-surface, employing constant geometry coefficients and profiles along the radial direction, but still keeping finite gradients and magnetic shear (as is usually done in local simulations in flux-tube codes), and using radially periodic boundary conditions. We selected the flux-surface where the heat flux was is at its maximum, i.e., ρ tor = 0.7, and performed two simulations: one for the High-Mirror and another for the Low-Mirror configuration. The resolution used is the same as in the global case, but the simulations are nevertheless cheaper (around 0.1 million CPU-hours) since they need less time to reach a quasi-stationary state, and are also simpler to run because they do not require sources to maintain the equilibrium profiles. We note that in this position, the normalized system size is ρ = ρ i /a ≈ 1/250, with ρ i being the ion Larmor

radius.
We observe that the trend of the Low-Mirror configuration producing more heat flux than the High-Mirror is still present in the radially local simulations ( figure 11 (top)). In addition, we also observe a close agreement between global and radially local simulations regarding the shape of the heat flux spectra ( figure 11 (bottom)). However, in both cases, the amplitudes of the fluxes are largely overestimated (up to a factor of 2.5) in local simulations.
Therefore, these results indicate that for a qualitative study of the effect of the geometry on ITG turbulence, radially local simulations are sufficient. However, for a quantitative study, such as the validation of gyrokinetic simulations against experimental results, global simulations are generally expected to be necessary. Furthermore, global nonlinear simulations will still be required in cases where effects not captured in radially local simulations, such as the shear of the radial electric field [22], are expected to play a key role in regulating turbulence.

VII. CONCLUSIONS
In the present paper, we performed linear and nonlinear global simulations of W7-X with GENE-3D to understand how ITG turbulence depends on a particular configuration used in W7-X.
We found that the Low-Mirror configuration produces more turbulent transport than both the High-Mirror and the Standard configurations. We observed that nonlinear simulations were necessary to correctly characterize how ITG turbulence is affected by the different magnetic geometries. The reason is that the most unstable mode in linear simulations is a large toroidal mode number, which does not contribute significantly to transport, and therefore it does not provide relevant information of the nonlinear system. Finally, we found that a radially local model can capture the trends observed in the global model qualitatively. There were, however, quantitative differences in the amplitudes, showing that nonlinear global simulations are generally expected to be necessary to validate the gyrokinetic results against experimental observations in W7-X.
The results presented in this paper thus represent an ongoing effort in the gyrokinetic stellarator community in moving forward to perform nonlinear global simulations [45][46][47]. In future work, we would like to perform further studies, including a more comprehensive plasma description in the simulations. In particular, we would like to treat the electrons as a kinetic species, add a radial electric field and electromagnetic effects. The goal is to perform a gyrokinetic validation study for W7-X, which will then give us confidence in our model and could help the operation and design of new machines.