2. Physical and numerical setting

Plasma microinstabilities are described mathematically by the gyrokinetic equation, which is an evolution equation for the gyrocentre distribution function $\delta G_s$, where $\delta G_s$ is the fluctuating part of the first-order correction to the equilibrium distribution function $f_s=F_{0s}+\delta f_s\equiv F_{0s}(1-Z_se\varphi (\boldsymbol {r},t)/T_{0,s})+\delta G_s(\boldsymbol {R}_s,\mu,\mathcal {E},t)+\mathcal {O}{\epsilon ^2}$. $\varphi$ is the electrostatic potential, $Z_s e$ the electric charge and $T_{0s}$ the equilibrium temperature for the species $s$. For an electrostatic plasma the gyrokinetic equation reads:

(2.1)\begin{align} & \left(\frac{\partial}{\partial t} +v_{{\parallel} s}\boldsymbol{\nabla}_{{\parallel}}+\boldsymbol{v}_{ds}\boldsymbol{\cdot} \boldsymbol{\nabla}\right)\delta G_s=\frac{Z_s eF_{0s}}{T_{0s}} \frac{\partial}{\partial t} \langle\varphi\rangle_{\boldsymbol{R}_s}\nonumber\\ & \quad-\frac{c}{B}\hat{\boldsymbol{b}}\boldsymbol{\cdot} \boldsymbol{\nabla}\langle\varphi\rangle_{\boldsymbol{R}_s} \times\boldsymbol{\nabla}{F_{0s}}-\frac{c}{B}\hat{\boldsymbol{b}}\boldsymbol{\cdot} \boldsymbol{\nabla}{\langle\varphi\rangle_{\boldsymbol{R}_s}}\times\boldsymbol{\nabla}{\delta G_s}. \end{align}

Collisions are neglected. The equilibrium distribution function $F_{0s}$ is taken to be Maxwellian and $\delta f_s/F_{0s}\sim k_{\parallel }/k_{\perp }\sim \rho ^*\equiv \epsilon \ll 1$, where $k_{\parallel }$ and $k_{\perp }$ are wave numbers along and across the equilibrium magnetic field $\boldsymbol {B}$. Here $\rho ^*\equiv \rho _s/L$ is the expansion parameter with $\rho _s=v_{{\rm th}s}/\varOmega _s$ the Larmor radius, $v_{{\rm th}s}=\sqrt {2T_s/m_s}$ the thermal velocity and $L$ the equilibrium space scale. The gyrocentre distribution function $\delta G_s$ depends on the gyrocentre position $\boldsymbol {R}_s=\boldsymbol {r}+\boldsymbol {v}_{\perp }\times \hat {\boldsymbol {b}}/\varOmega _s$, with $\boldsymbol {r}$ being the particle position, $\varOmega _s=Z_s e B/(m_sc)$ the gyrofrequency and $\hat {\boldsymbol {b}}=\boldsymbol {B}/B$. The velocity-space coordinates are $\mu _s=m_sv_{\perp }^2/(2B)$ and $\mathcal {E}=m_s v^2/2$. The parallel gyrocentre velocity can thus be written as $v_{\parallel s}=\pm \sqrt {2(\mathcal {E}-B\mu _s)/m_s}$, whereas $\boldsymbol {v}_{ds}$ is the gyrocentre drift velocity $\boldsymbol {v}_{ds}=(v_{\parallel }/\varOmega _s) \boldsymbol {\nabla } \times (v_{\parallel }\hat {\boldsymbol {b}})$, where $\boldsymbol {\nabla }=\partial /\partial \boldsymbol {R}_s$. The gyroaveraged electrostatic potential $\langle \varphi \rangle _{\boldsymbol {R}_s}$ is Fourier transformed with respect to $\boldsymbol {R}_s$, $\langle \varphi \rangle _{\boldsymbol {R}_s}= \sum _{\boldsymbol {k}}\langle \varphi \rangle _{\boldsymbol {R}_s, \boldsymbol {k}}\exp ({\rm i}\boldsymbol {k}\boldsymbol {\cdot } \boldsymbol {R}_s)$, where the coefficients are related to the coefficients of the transformation with respect to $\boldsymbol {r}$ through the Bessel function $\langle \varphi \rangle _{\boldsymbol {R}_s, \boldsymbol {k}}=J_0(a_s)\varphi _{\boldsymbol {k}}$, $a^2_s=2B\mu k_{\perp }^2/\varOmega _s^2$. The electrostatic potential $\varphi$ is determined by the quasineutrality condition

(2.2)\begin{equation} \sum_s Z_se\int{\rm d}^3{\boldsymbol{v}}F_{0s}\frac{Z_se\varphi}{T_{0s}}=\sum_s Z_se\int{\rm d}^3{\boldsymbol{v}}\sum_{\boldsymbol{k}} \text{e}^{{\rm i} \boldsymbol{k}\boldsymbol{\cdot} \boldsymbol{r}}J_0(a_s)\delta G_s. \end{equation}

We perform numerical studies treating both ions and electrons kinetically, but neglecting electromagnetic effects. The code used is particularly suited to this task, as it is based on a semi-implicit numerical scheme. Our principal goal is to study the behaviour of gyrokinetic turbulence close to marginality. Simulating these regimes is demanding. For weakly driven instabilities, a saturated turbulent state is difficult to reach, which translates into longer and computationally more expensive simulations than for strongly driven turbulence. For these reasons, we choose the flux tube version of the $\delta f$ gyrokinetic code stella (Barnes et al. Reference Barnes, Parra and Landreman2019). This code solves the Fourier-decomposed nonlinear gyrokinetic equation (2.1), multiplied by the normalising factor $(a^2/\rho _{\rm ref} v_{\rm th,ref})\exp (-v^2/v_{{\rm th}s}^2)/F_{0s}$, where $v_{\rm th,ref}$ denotes the thermal velocity of the user-specified reference particle species, having mass $m_{\rm ref}$ and temperature $T_{\rm ref}$. Furthermore, $v_{{\rm th}s}$ is the thermal velocity of the particle species $s$, $\rho _{\rm ref}=v_{\rm th,ref}/\varOmega _{\rm ref}$ is the Larmor radius of the reference particle species, where $\varOmega _{\rm ref}=Z_{\rm ref} e B_{\rm ref}/(m_{\rm ref} c)$ is the gyrofrequency, and $a$ and $B_{\rm ref}$ are, respectively, the reference length and reference magnetic field, which depend on the choice of the magnetic geometry.

The coordinates used by stella are $(x,y,z,v_{\parallel },\mu )$, where $(x,y)$ are, respectively, the radial and binormal coordinates in the plane perpendicular to $\hat {\boldsymbol {b}}$. The magnetic field vector is expressed as $\boldsymbol {B}=\boldsymbol {\nabla }{\alpha }\times \boldsymbol {\nabla }{\psi }$, where $\alpha = \theta - \iota \zeta$ is a coordinate that labels field lines and $\psi$ flux surfaces. Here, $\iota$ is the rotational transform, and $\theta$ and $\zeta$ denote poloidal and toroidal PEST flux coordinates (Grimm, Dewar & Manickam Reference Grimm, Dewar and Manickam1983), respectively. In terms of these coordinates, $x$ and $y$ are defined as

(2.3)\begin{equation} \begin{array}{c@{}} {x = r-r_0},\\ {y = r_0(\alpha-\alpha_0)}, \end{array}\end{equation}

where $r_0$ and $\alpha _0$ are input parameters that specify the radial position and the field line at the centre of the simulation domain. In particular, $r$ is defined as $r\equiv a\sqrt {\psi _\textrm {t}/\psi _{\textrm {t},\textrm {LCFS}}}$, where $\psi _\textrm {t}$ is the toroidal magnetic flux and $\psi _{\textrm {t},\textrm {LCFS}}$ its value at the last closed flux surface. The perpendicular wave number $\boldsymbol {k}_{\perp }$ can thus be decomposed as $\boldsymbol {k}_{\perp }=k_x\boldsymbol {\nabla }{x} +k_y\boldsymbol {\nabla }{y}$, where $k_x$ is the radial wave number and $k_y$ the binormal wave number. Furthermore, $z$ is a coordinate measuring the location along the magnetic field and, for stellarator simulations, is equal to the toroidal angle $\zeta$ measured in radians.

All the simulations presented in this paper are run in stellarator geometry. The equilibrium is passed to stella as an input file obtained through the equilibrium solver code VMEC (Hirshman & Whitson Reference Hirshman and Whitson1983). The simulations are performed either in the W7-X high-mirror (KJM) configuration (Geiger et al. Reference Geiger, Beidler, Feng, Maaßberg, Marushchenko and Turkin2014) or in the low-mirror (AIM) configuration. The main difference between these two configurations is their mirror ratio, i.e. the ratio between the maximum and minimum magnetic field strength on the magnetic axis. The high-mirror configuration is of particular interest for experiments as it is characterised by a small bootstrap current.

3. Stability diagrams

We begin by constructing a stability diagram for the various electrostatic microinstabilities that could harm plasma confinement, in the spirit of Alcusón et al. (Reference Alcusón, Xanthopoulos, Plunk, Helander, Wilms, Turkin, von Stechow and Grulke2020) but carefully monitoring the wave number dependence of the various modes. Selecting the most unstable eigenvalue of microinstability spectra is not always the most conclusive procedure, since the relevance of an instability for transport depends on the scale at which the instability is active. We perform a scan over the normalised ion temperature gradient $a/L_\textrm {Ti}$ and normalised density gradient $a/L_\textrm {n}$, where $1/L_\textrm {Ti}=-T^{-1}_{0\textrm {i}}\, \textrm {d}T_{0\textrm {i}}/{\textrm {d} x}$ and $1/L_\textrm {n}=-n^{-1}_{0}\, \textrm {d} n_{0}/{\textrm {d} x}$ are the characteristic gradients length scales and $a$ is the average minor radius of the plasma. The scan is performed in the range $[0.0, 3.0]$ for both gradients, with a step size of 0.5. Since we are simulating two kinetic species, we expect to observe ITG modes, density-driven trapped electron modes (TEMs) (Proll et al. Reference Proll, Helander, Connor and Plunk2012; Helander et al. Reference Helander, Bird, Jenko, Kleiber, Plunk, Proll, Riemann and Xanthopoulos2015), ion-driven trapped electron modes (iTEMs) (Proll et al. Reference Proll, Mynick, Xanthopoulos, Lazerson and Faber2015; Plunk, Connor & Helander Reference Plunk, Connor and Helander2017a), and universal instabilities (Coppi & Pegoraro Reference Coppi and Pegoraro1977; Smolyakov et al. Reference Smolyakov, Yagi and Kishimoto2002; Helander & Plunk Reference Helander and Plunk2015; Landreman et al. Reference Landreman, Antonsen and Dorland2015). We consider equal temperatures $\tau \equiv T_\textrm {i}/T_\textrm {e}=1$ and density gradients. We set the electron temperature gradient to zero. We acknowledge the importance that electron temperature gradients can have at the ion scale in predominantly electron-heated plasmas, such as in W7-X (Wilms, Bañón Navarro & Jenko Reference Wilms, Bañón Navarro and Jenko2023). However, in this work we decide to limit ourselves to at most density-gradient-driven electron modes and their effect on ion physics. The role of electron temperature gradients in W7-X will be addressed in dedicated studies (Zocco et al. Reference Zocco, Podavini, Wilms, Bañón Navarro and Jenko2024).

The radial location chosen is $r_0/a=0.7$, where electrostatic fluctuations have been detected in W7-X (Bähner et al. Reference Bähner, Alcusón, Hansen, von Stechow, Grulke, Windisch, Smith, Huang, Edlund and Porkolab2021). The chosen equilibrium is a high-mirror (KJM) W7-X magnetic configuration. At this radial location, for the chosen equilibrium, the local value of the global magnetic shear is $\hat {s}=-(r/\iota ) \,\textrm {d} \iota /\textrm {d} r=-0.1249$ and the safety factor $q=1/\iota =1.103$. The resolution of the simulations is set by the input parameters: $N_z\times N_{k_y}\times N_{k_x}\times N_{v_{\parallel }}\times N_{\mu }=384\times 70\times 1\times 64\times 24$. The maximum magnetic moment $\hat {\mu }_\textrm {max}$ is set by the maximum perpendicular velocity $\hat {v}_{\perp, \textrm {max}}=3$. The maximum parallel velocity is $|\hat {v}_{\parallel, \textrm {max}}|=3$, where $\hat {\mu }_s=\mu _s B_\textrm {ref}/(2T_s)$, $\hat {v}_{\perp s}=v_{\perp s}/v_{\textrm {th}s}$ and $\hat {v}_{\parallel s}=v_{\parallel s}/v_{\textrm {th}s}$. The number of binormal wave numbers corresponds to a scan in the range $k_y\rho _\textrm {i}\in [0.05,3.0]$, while the radial wave number $k_x\rho _\textrm {i}$ is set to zero. We choose to run the simulations in the flux tube denoted by $\alpha _0=0$, centred around $(\theta, \zeta )=(0,0)$, at the bean-shaped cross section, which has been characterised extensively (González-Jerez et al. Reference González-Jerez, Xanthopoulos, García-Regaña, Calvo, Alcusón, Navarro, Barnes, Parra and Geiger2022) and where turbulence is found to be significant (Helander et al. Reference Helander, Beidler, Bird, Drevlak, Feng, Hatzky, Jenko, Kleiber, Proll and Turkin2012). The flux tube length covers three full toroidal turns of the device for reasons to be discussed later.

We construct three types of stability diagrams. This is because care must be taken to avoid overlooking instabilities with growth rates smaller than, but comparable to, the fastest growing one, especially at long wavelengths. Different instabilities are recognised by the sign of their normalised frequency, which may exhibit sudden jumps as the binormal wave number $k_y\rho _\textrm {i}$ is varied, since the code always selects the fastest growing instability. The three stability diagrams reflect the three main $k_y\rho _\textrm {i}$ ranges where instabilities are found, as displayed in figure 1. At the longest wavelengths, $k_y\rho _\textrm {i}<1.0$, we find instabilities with negative frequency, which corresponds to rotation in the electron diamagnetic direction. At shorter wavelengths, we find two distinct instabilities having positive frequency, thus rotating in the ion diamagnetic direction. They are unstable for different $k_y\rho _\textrm {i}: k_y\rho _\textrm {i}\lesssim 1.5$ and $1.5\lesssim k_y\rho _\textrm {i}\lesssim 2.5$, respectively. We thus separate them into different diagrams, shown in figure 2(a–c), which are devised to show the normalised growth rate for the fastest growing mode in a selected $k_y\rho _\textrm {i}$ range. Each pixel corresponds to a single simulation having a specific value of the parameter $\{a/L_\textrm {Ti}, a/L_\textrm {n}\}$.

Figure 1. Normalised growth rate $\gamma a/v_\textrm {thi}$ and frequency $\omega a/v_\textrm {thi}$ spectra as a function of the normalised binormal wave number $k_y\rho _\textrm {i}$ for simulations with density gradient $a/L_\textrm {n}=3.0$ but different ion temperatuere gradients $a/L_\textrm {Ti}$.

Figure 2. Stability diagrams showing the normalised growth rate $\gamma a/v_\textrm {thi}$ as a function of the ion temperature gradient $a/L_\textrm {Ti}$ and density gradient $a/L_\textrm {n}$. Each diagram shows instabilities having different normalised frequency $\omega a/v_\textrm {thi}$ signs and active at different $k_y\rho _\textrm {i}$: (a) $\omega a/v_\textrm {thi}<0$, $k_y\rho _\textrm {i}<1.0$; (b) $\omega a/v_\textrm {thi}>0$, $k_y\rho _\textrm {i} \lesssim 1.5$; and (c) $\omega a/v_\textrm {thi}>0$, $1.5 \lesssim k_y\rho _\textrm {i} \lesssim 2.5$.

The first diagram (figure 2a) is populated by instabilities found at large scales and rotating in the electron diamagnetic direction. A critical value of the density gradient $a/L_\textrm {n}\approx 1.5$ is necessary to destabilise them, whereas the ion temperature gradient plays a stabilising role. In order to identify these instabilities, we study the parallel structure of the modulus of the normalised electrostatic potential $|\hat {\varphi }_{\boldsymbol {k}}|$, where $\hat {\varphi }_{\boldsymbol {k}}=(a/\rho _\textrm {i})(e/T_\textrm {i})\varphi _{\boldsymbol {k}}$. In particular, we compare it with two quantities useful for the instabilities identification: the value of the normalised magnetic field $B/B_\textrm {ref}$ and the normalised $\boldsymbol {\nabla }{B}$ drift in the $\boldsymbol {\nabla }{\alpha }$ direction $2a^2B_\textrm {ref}B^{-2}(\hat {\boldsymbol {b}}\times \boldsymbol {\nabla }{B})\boldsymbol {\cdot } \boldsymbol {\nabla }{\alpha }$, which we denote here as $\kappa _2$. The former shows the presence of magnetic wells, whereas from the latter the regions of ‘bad curvature’ can be identified, denoted here by positive values of the magnetic drift.

Figure 3 shows the electrostatic potential for one of the simulations composing the first diagram. We note that the eigenfunction peaks in regions of bad curvature and is characterised by two symmetric maxima at finite $z$. The structure along $z$ is broad and decays slowly along the field line. A comparison with previous works suggests that we can identify them as universal instabilities (Coppi & Pegoraro Reference Coppi and Pegoraro1977; Smolyakov et al. Reference Smolyakov, Yagi and Kishimoto2002; Helander & Plunk Reference Helander and Plunk2015; Landreman et al. Reference Landreman, Antonsen and Dorland2015). In contrast to earlier numerical simulations (Costello et al. Reference Costello, Proll, Plunk, Pueschel and Alcusón2023), we find them also in the presence of a non-zero ion temperature gradient, which somewhat stabilises them but does not suppress them completely. We stress that the simulation of a long flux tube was necessary to obtain this result. Unless multiple toroidal turns are included in the simulation domain, it is not possible to resolve the maxima at finite $z$. We show in figure 4 that if an additional toroidal turn is simulated, other, less-pronounced local maxima are found, but overall the eigenfunction decays along the field line. Three toroidal turns are thus enough to capture the main structure of the mode, and four are useful for obtaining a substantial decay of the eigenfuction towards the ends of the domain.

Figure 3. Electrostatic potential structure along the magnetic field line for $k_y\rho _\textrm {i}=0.43$, for the simulation with $\{a/L_\textrm {Ti}, a/L_\textrm {n}\}=\{0.5,3.0\}$ (in blue) compared with the normalised magnetic field strength $B/B_\textrm {ref}$ (in grey) and the magnetic drift $\kappa _2$ (in dotted red).

Figure 4. Electrostatic potential structure along the magnetic field line for $k_y\rho _\textrm {i}=0.38$, for the simulation with $\{a/L_\textrm {Ti}, a/L_\textrm {n}\}=\{0.5,3.0\}$ (in blue) compared with the normalised magnetic field strength $B/B_\textrm {ref}$ (in grey) and the magnetic drift $\kappa _2$ (in dotted red). Result for a flux tube that extends four times toroidally.

The second diagram (figure 2b) is populated by instabilities with wave numbers up to $k_y\rho _\textrm {i} \lesssim 1.5$ rotating in the ion direction. On the right-hand side of the diagram, for $a/L_\textrm {n} \ge 1.5$, we mainly find iTEM instabilities, characterised by eigenfunctions that peak where magnetic wells and bad curvature regions overlap (see figure 5). These instabilities are not very sensitive to variations in the ion temperature gradient as long as $a/L_\textrm {n} \ge 1.5$, but below this value the growth rate increases with increasing $a/L_\textrm {Ti}$. As we can observe in figure 5, this can be explained by a transition from an iTEM to an ITG mode. The latter can be recognised by its strong peaking in the region of bad curvature at $z=0$. When the density gradient is further decreased, the ITG instability gradually comes to dominate for every value of $a/L_\textrm {Ti}$ and the growth rate becomes sensitive to the value of $a/L_\textrm {Ti}$, as expected.

Figure 5. Electrostatic potential structure along the magnetic field line for $k_y\rho _\textrm {i}=0.99$ (in green) and $k_y\rho _\textrm {i}=1.42$ (in blue) for simulations with $a/L_\textrm {Ti}=\{1.5,3.0\}$ and $a/L_\textrm {n}=1.5$, compared with the normalised magnetic field strength $B/B_\textrm {ref}$ (in grey) and the magnetic drift $\kappa _2$ (in dotted red).

The third diagram (figure 2c) describes instabilities that are found at smaller scales and rotate in the ion diamagnetic direction. In the right part of the diagram, the growth rate increases with the density gradient and slightly decreases with increasing ion temperature gradient. Figure 6 shows an example of such an instability, which features local maxima corresponding to magnetic wells locations, and a local minimum in $z=0$. This type of mode does not possess any specific feature that could characterise it as ITG or iTEM, but rather a hybrid mode between the two. For $a/L_\textrm {n}<2$, however, the growth rate increases again with increased $a/L_\textrm {Ti}$. This suggests that a branch of the ITG mode is destabilised at smaller scales too.

Figure 6. Electrostatic potential structure along the magnetic field line for $k_y\rho _\textrm {i}=1.55$, for the simulation with $\{a/L_\textrm {Ti}, a/L_\textrm {n}\}=\{3.0,3.0\}$ (in blue) compared with the normalised magnetic field strength $B/B_\textrm {ref}$ (in grey) and the magnetic drift $\kappa _2$ (in dotted red).

We finally combine the three diagrams into one that only shows the growth rate for the fastest growing mode in the whole simulated $k_y\rho _\textrm {i}$ domain, as a function of $a/L_\textrm {Ti}$ and $a/L_\textrm {n}$ (figure 7). The comparison with the distinct diagrams shows how essential it is to separate different scales in order to understand which instabilities are active at which scales and for which gradients. The universal modes present in figure 2(a), for example, feature slightly smaller growth rates than other instabilities and are thus overlooked when only considering the fastest growing mode. But since they are active at larger scales, they might be more relevant for transport than modes of shorter wavelength.

Figure 7. Stability diagram showing the normalised growth rate $\gamma a/v_\textrm {thi}$ as a function of the ion temperature gradient $a/L_\textrm {Ti}$ and density gradient $a/L_\textrm {n}$. Plotted here is the largest growth rate for each set of $\{a/L_\textrm {Ti}, a/L_\textrm {n}\}$, on which an interpolation has been performed.

In order to check our hypothesis, we perform a quasi-linear estimate for the heat flux (Mariani et al. Reference Mariani, Brunner, Dominski, Merle, Merlo, Sauter, Görler, Jenko and Told2018)

(3.1)\begin{equation} Q_{\rm QL}=\frac{\hat{\gamma}}{\langle k_{{\perp}}^2\rangle},\end{equation}

where $\hat {\gamma }=\gamma a/v_\textrm {thi}$ is taken from figure 2(a–c) and

(3.2)\begin{equation} \langle k_{{\perp}}^2\rangle=\frac{\displaystyle \int k_{{\perp}}^2|\varphi|^2\,{\rm d} l/B}{\displaystyle \int|\varphi|^2\,{\rm d} l/B} . \end{equation}

The result is reported in figure 8(a–c) and is consistent with the hypothesis that universal modes can play a pivotal role. We observe that, while their growth rate is a decreasing function of $a/L_\textrm {Ti}$, $Q_\textrm {QL}$ is persistent for large $a/L_\textrm {Ti}$ (figure 8a). ITG modes, top left of figure 8(b), are still relevant to transport, but not as much as universal instabilities. iTEMs, bottom right of figure 8(b), on the other hand, do not seem to be as important. All remaining instabilities in figure 8(c) show a quasi-linear estimate more than a factor 10 smaller than those already discussed in figure 8(a,b).

Figure 8. Stability diagrams showing the quasi-linear heat flux estimate for the instabilities reported in figure 2(a–c) calculated as in (3.1). Each diagram refers to instabilities having different normalised frequency $\omega a/v_\textrm {thi}$ signs and active at different $k_y\rho _\textrm {i}$: (a) $\omega a/v_\textrm {thi}<0$, $k_y\rho _\textrm {i}<1.0$; (b) $\omega a/v_\textrm {thi}>0$, $k_y\rho _\textrm {i} \lesssim 1.5$; and (c) $\omega a/v_\textrm {thi}>0$, $1.5 \lesssim k_y\rho _\textrm {i} \lesssim 2.5$.

A feature of figure 7 is the presence of a ‘stability valley’ at large gradients: as already noted in Helander et al. (Reference Helander, Bird, Jenko, Kleiber, Plunk, Proll, Riemann and Xanthopoulos2015) and more clearly emphasised in Alcusón et al. (Reference Alcusón, Xanthopoulos, Plunk, Helander, Wilms, Turkin, von Stechow and Grulke2020), the growth rate does not increase with increasing gradients everywhere in the diagram. In particular, the addition of a non-zero density gradient to a plasma with a moderate ion temperature gradient reduces the growth rate. The term ‘stability valley’ was coined by Alcusón et al. (Reference Alcusón, Xanthopoulos, Plunk, Helander, Wilms, Turkin, von Stechow and Grulke2020), who calculated stability diagrams up to quite large gradients, which brings out this feature particularly clearly. For realistic values of the gradients, however, the ‘valley’ is shallow and the growth rate only drops slightly. Perhaps a more important feature is that the fastest growing mode acquires a shorter wavelength, which reduces the turbulent transport (Helander et al. Reference Helander, Bird, Jenko, Kleiber, Plunk, Proll, Riemann and Xanthopoulos2015; Xanthopoulos et al. Reference Xanthopoulos, Bozhenkov, Beurskens, Smith, Plunk, Helander, Beidler, Alcusón, Alonso, Dinklage, Ford, Fuchert, Geiger, Proll, Pueschel, Turkin and Warmer2020). Indeed, despite the mild decrease of the growth rate at experimentally relevant parameters, nonlinear simulations show a strong decrease of fluctuations and particle fluxes (García-Regaña et al. Reference García-Regaña, Barnes, Calvo, González-Jerez, Thienpondt, Sánchez, Parra and St. Onge2021) along with heat fluxes (Thienpondt et al. Reference Thienpondt, García-Regaña, Calvo, Barnes, Parra and Davies2022, Reference Thienpondt, García-Regaña, Calvo, Georgia and Barnes2024).

At large values of the density gradient, i.e. on the right-hand side of figure 7, there are rapidly growing universal and iTEM instabilities with long and short wavelengths, respectively. These instabilities were observed but not classified by Alcusón et al. (Reference Alcusón, Xanthopoulos, Plunk, Helander, Wilms, Turkin, von Stechow and Grulke2020).

Measurements of W7-X temperature and density profiles suggest that normalised temperature gradients in most of the core do not exceed values of $a/L_\textrm {Ti}\approx 3.0$ and density gradients of $a/L_\textrm {n}\approx 1.5$. Such regimes are described by the left part of figure 2(b), which is dominated by ITG instabilities. For this reason, in the following sections, we present an extensive study on the ITG, focusing in particular on its behaviour close to its linear stability threshold.

4. Linear characterisation of the ITG mode near marginality

In order to study the ITG mode, we run further linear, flux tube, electrostatic simulations. We retain electron kinetic effects so as to have two kinetic species with equal temperatures. As seen in the stability diagram, the ITG mode is the dominant instability when the density gradient is small enough. Accordingly, we consider the limit of flat density profile ($a/L_\textrm {n}=0$) and also set $a/L_\textrm {Te}=0$, performing a scan in the normalised ion temperature gradient $a/L_\textrm {Ti}$ only, with values $a/L_\textrm {Ti}\in [0.5,3.0]$. Our aim is to ascertain whether ion-temperature-gradient-driven modes, previously observed with adiabatic electrons (Zocco et al. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018), persist when the electron dynamics is accounted for. We run our simulations in the same flux tube used for the stability diagram, $\alpha _0=0$, at the radial location $r_0/a=0.7$ and centred around $(\theta,\zeta )=(0,0)$ in the W7-X high-mirror magnetic configuration. Its extension is three toroidal turns. The resolution of the simulations is set by the input parameters: $N_z\times N_{k_y}\times N_{k_x}\times N_{v_{\parallel }}\times N_{\mu }=256\times 62\times 1\times 64\times 24$. The maximum magnetic moment $\hat {\mu }_\textrm {max}$ is set by the maximum normalised perpendicular velocity $\hat {v}_{\perp, \textrm {max}}=3$. The maximum normalised parallel velocity is $|\hat {v}_{\parallel, \textrm {max}}|=3$. The normalised binormal wave number is varied over the range $k_y\rho _\textrm {i}\in [0.05,2.0]$ whereas the normalised radial wave number $k_x\rho _\textrm {i}$ is set to zero.

The result is presented in figure 9, where plots of the normalised growth rate $\gamma a/v_\textrm {thi}$ and normalised frequency $\omega a/v_\textrm {thi}$ are shown. For each case we record the growth rate and the frequency of the fastest growing mode and plot these as functions of the normalised ion temperature gradient in figure 10. We observe that for large ion temperature gradients, the fastest growing mode is found at the ion scale, $k_y\rho _\textrm {i}\approx 1$, consistently with a curvature-driven ITG. When the gradient is decreased, the curvature-driven ITG mode gets progressively stabilised, while other ‘background’ modes emerge. The latter are characterised by lower values of $k_y\rho _\textrm {i}$ and are important only for sufficiently small values of the temperature gradient. We note that for a given $a/L_\textrm {Ti}$, different wavelengths can belong to different ITG branches. As a result, in the plot of the fastest-growing instability, there may be sudden jumps in the real part of the frequency. The background modes have been previously observed in low-shear tokamak and W7-X simulations featuring adiabatic electrons. Their nature can neither be attributed to the slab nor to the toroidal branch of the ITG mode. They are instead of the Floquet-type (Zocco et al. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018) and survive also when electrons are treated kinetically. A prominent feature of Floquet modes is their extended nature along magnetic field lines (Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022).

Figure 9. Spectra of the normalised linear growth rate $\gamma a/v_\textrm {thi}$ and frequency $\omega a/v_\textrm {thi}$ as a function of the normalised binormal wave number $k_y\rho _\textrm {i}$ for different temperature gradients. The growth rate and frequency for the fastest growing mode are highlighted.

Figure 10. Normalised linear growth rates $\gamma a/v_\textrm {thi}$ (in blue) and frequencies $\omega a/v_\textrm {thi}$ (in red) for the fastest growing mode as a function of the ion temperature gradient showing the transition to Floquet background modes for the high-mirror and low-mirror configurations.

Blank spots in the spectra in figure 9 represent values of $k_y\rho _\textrm {i}$ for which modes belonging to different branches are characterised by similar growth rates. In these cases, it is difficult to obtain a numerically converged spectrum. We thus plot the value of the growth rate and frequency by averaging over the last 10 % of time steps of the simulation and by considering only points that are converged with a relative error of the order of a few per cent.

In figure 10 we observe that the transition to Floquet-type modes is characterised by a change in the slope of the growth rate with the temperature gradient, which becomes less steep at small values of $a/L_\textrm {Ti}$, and by jumps in the frequency. The transition occurs between $a/L_\textrm {Ti}=1.0$ and $a/L_\textrm {Ti}=1.1$. In order to locate the linear stability threshold, we perform a linear interpolation of the last points corresponding to Floquet modes, and we find that the linear instability threshold is $a/L_\textrm {Ti}\lesssim 0.5$, which is consistent with what was found in Zocco et al. (Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018, Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022). The linear results with kinetic electrons thus confirm the original findings obtained with adiabatic electrons. The same computation has been repeated considering the same input parameters but running the simulations in a W7-X low-mirror configuration (AIM). The result is also reported in figure 10. We do not observe any particular dependence on the geometry, especially for the Floquet modes. The transition is visible in both geometries but is somewhat less pronounced in the low-mirror configuration. The values of the growth rates and frequencies tend to the same value as the gradient is progressively decreased. This suggests that such modes are relatively persistent and insensitive to the mirror ratio.

The transition between the two different instability branches is also visible from the change in the parallel structure of the electrostatic potential (Zocco et al. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018, Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022). To illustrate this fact, we choose two temperature gradients straddling the transition, $a/L_\textrm {Ti}=0.9$ and $a/L_\textrm {Ti}=2.0$. The first case corresponds to the Floquet-type ITG branch, whereas the second lies in the toroidal strongly driven branch. The ITG destabilisation for the latter can be described with a simple fluid picture. For this reason, we refer to it as fluid-like. We plot the modulus of the normalised electrostatic potential for the most unstable wave number divided by its value at $z=0$, $|\hat {\varphi }_{\boldsymbol {k}}| / |\hat {\varphi }_{\boldsymbol {k}}|_{z=0}$ in figure 11. We observe that, also in this case, the qualitative picture found with adiabatic electrons is unchanged: the Floquet-like eigenfunction ($a/L_\textrm {Ti}=0.9$) decays more slowly along the magnetic field line than the fluid-like eigenfunction ($a/L_\textrm {Ti}=2.0$). The mode close to marginal stability thus requires a longer flux tube (and more points along the field line) for its parallel structure to be sufficiently resolved. To choose the suitable flux tube length we run scans of the nfield_periods input parameter, $N_\textrm {fp}$, to set the number of magnetic periods covered by the flux tube. One toroidal turn around W7-X corresponds to $N_\textrm {fp}=5$, since W7-X is a five-fold symmetric device. We compare the eigenfunction for the Floquet-type ITG mode in flux tubes extending one, two, three and four times around the torus. The result in figure 12 shows how at least two toroidal turns are necessary to sufficiently resolve the parallel structure of the electrostatic potential. With at least three toroidal turns, a substantial decay of the eigenfunction towards the end of the domain is observed, whereas just one toroidal turn is not sufficient to resolve the slowly decaying structure of the mode. The choice of the flux tube is thus of crucial relevance for obtaining sufficient energy injection, in particular in nonlinear simulations. For this reason, a follow-up discussion about suitable choices of input parameters in the parallel direction is given in the next section.

Figure 11. Fastest growing mode electrostatic potential structure along the field line in the fluid-like ($a/L_\textrm {Ti}=2.0$) case (in red) and in the Floquet-type ($a/L_\textrm {Ti}=0.9$) case (in blue).

Figure 12. Comparison between electrostatic potential structures along the field line obtained with simulated flux tubes of different lengths for the case $a/L_\textrm {Ti}=0.9$.

5. Nonlinear results

We now move on to characterise differences between Floquet-type and fluid-like turbulence through nonlinear simulations using the same two values of the ion temperature gradient as in the previous section, $a/L_\textrm {Ti}=0.9$ and $a/L_\textrm {Ti}=2.0$. The input parameters are listed in table 1. Those for the simulation in the fluid-like regime ($a/L_\textrm {Ti}=2.0$) are chosen based on previous results (García-Regaña et al. Reference García-Regaña, Barnes, Calvo, Parra, Alcusón, Davies, González-Jerez, Mollén, Sánchez and Velasco2021), whereas closer to marginality ($a/L_\textrm {Ti}=0.9$), the parameters are adjusted in accordance with the linear results mentioned previously. Given the set of parameters displayed in the table, we are able to simulate the following wave numbers: in the fluid-like regime $k_{x,\textrm {max}}\rho _\textrm {i}\equiv \hat {k}_{x,\textrm {max}}\approx 1.19$ and $\hat {k}_{y,\textrm {max}}= 2.5$, with a respective step of $\Delta \hat {k}_x\approx 0.04$ and $\Delta \hat {k}_y= 0.05$, whereas closer to marginality $\hat {k}_{x,\textrm {max}}\approx 2.06$ and $\hat {k}_{y,\textrm {max}}= 2.95$, with steps of $\Delta \hat {k}_x\approx 0.07$ and $\Delta \hat {k}_y= 0.05$.

Table 1. Set of input parameters for the nonlinear simulations in the fluid-like limit ($a/L_\textrm {Ti}=2.0$) and closer to marginality ($a/L_\textrm {Ti}=0.9$).

We choose $N_\textrm {fp}=10$, i.e. two toroidal turns around the device, for the simulation close to marginality. This is different from what did linearly and the reason is computational. A longer flux tube needs a larger number, $N_z$, of points in the $z$ direction, which makes the simulation computationally more expensive. We run convergence tests to ensure that two toroidal turns are enough to resolve the parallel structure of the Floquet modes and enough to obtain the energy injection needed by the instability. We first consider the most unstable wave number at $a/L_\textrm {Ti}=0.9$ and plot the value of the normalised growth rate as a function of the number of simulated field periods in figure 13(a). We see that with at least three toroidal turns we obtain a good convergence, with only a variation of the growth rate of 5 % when an additional toroidal turn is simulated. We also note that $N_z=128$ is sufficient to resolve the parallel structure of the electrostatic potential. For $N_\textrm {fp}<15$, the value of the growth rate decreases by roughly 30 % for each toroidal turn that is removed from the simulation. When studying higher values of $k_y\rho _\textrm {i}$, as in figure 13(b), we note that a larger value of $N_z$ is required to obtain convergence for $N_\textrm {fp}=15$. For computational reasons, however, we settle on two toroidal turns, i.e. $N_\textrm {fp}=10$, with $N_z=128$, so as to obtain a less computationally expensive simulation but yet enough resolution for a sensible result. In fact, even if the wave numbers around $k_y\rho _\textrm {i}\approx 1$ are not the most unstable ones close to marginality, they still play a role and need to be properly resolved too. In the two plots, the importance of having enough parallel velocity space resolution is also evident.

Figure 13. Flux tube length convergence tests for growth rates at $a/L_\textrm {Ti} = 0.9$ for the wave numbers (a) $k_y\rho _\textrm {i}=0.71$ and (b) $k_y\rho _\textrm {i}=1.22$.

In Zocco et al. (Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022), the evolution of the electrostatic potential fluctuations was already studied in order to characterise the turbulence. We pay particular attention to the average squared modulus $\langle \hat {\varphi }^2_{\boldsymbol {k}} \rangle$ as a function of the normalised time $\hat {t}=tv_\textrm {thi}/a$. The time trace for $a/L_\textrm {Ti}=2.0$ is displayed in figure 14(a), and that for $a/L_\textrm {Ti}=0.9$ is displayed in figure 14(b). We specifically compare the time evolution of the zonal structures ($\hat {k}_y=0$) with that of the primary linear instabilities. The relative importance of the two gives rise to the definition of three time instants. At the first instant, zonal structures reach a minimum, after a transient phase. At the second, the amplitude of zonal flows equals the one of the main linear instabilities, and at the third, the electrostatic potential resulting from the summation of all the Fourier modes, $\langle \hat {\varphi }^2 \rangle =\sum _{k_x,k_y}\hat {\varphi }_{\boldsymbol {k}}$ reaches a maximum and saturates.

Figure 14. Time evolution of the average squared modulus of the electrostatic potential fluctuations $\langle \hat {\varphi }^2\rangle$ (in black) compared with its zonal component $\sum _{k_x}\langle \hat {\varphi }^2_{k_y=0}\rangle$ (in red) and the primaries (in blue) for (a) $a/L_\textrm {Ti}=2.0$ and (b) $a/L_\textrm {Ti}=0.9$. Time windows encompassing the three key time instants discussed in the text are displayed as yellow bands. Adapted from Zocco et al. (Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022).

We study the Fourier spectra of the fluctuations for each time instant, performing an average over a time window about each instant. We start with the fluid-like scenario ($a/L_\textrm {Ti}=2.0$). In the first time instant (figure 15a) the primary linear instability can be studied. It consists of an isotropic structure at $\boldsymbol {k}_1\equiv (\hat {k}_x,\hat {k}_y)=(0,1.05)$, typically observed in axisymmetric simulations, and two additional lobes at $\boldsymbol {k}_2\equiv (\hat {k}_x,\hat {k}_y)=(\pm 0.84,1.05)$. At the second instant, the evolution of the primary modes starts to depart from the linear phase and the zonal flows play a more prominent role. This can be considered as a pre-saturation phase. In this phase, the finite-$\hat {k}_x$ lobes get isotropised (figure 15b) and their relative importance decreases. A checkerboard-like structure forms and persists also in the saturating phase (figure 15c). From the time evolution in figure 14(a), we observe that after the saturating phase both the primary isotropic lobe and the finite-$\hat {k}_x$ lobes are suppressed by up to four orders of magnitude, whereas zonal flows remain active.

Figure 15. Fourier spectra of time-averaged $\langle \hat {\varphi }^2_{\boldsymbol {k}}\rangle$ as a function of the normalised binormal $k_y\rho _\textrm {i}$ and radial $k_x\rho _\textrm {i}$ wave numbers for $a/L_\textrm {Ti}=2.0$ at the three key time instants: (a) first time instant (adapted from Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022); (b) second time instant; and (c) third time instant.

Floquet-type turbulence ($a/L_\textrm {Ti}=0.9$) behaves differently but the first instant features similarities with the fluid-like case. In particular, both the isotropic structure and the finite-$\hat {k}_x$ lobes are present but are shifted to higher $\hat {k}_y$ (figure 16a). The isotropic structure is also reduced in relative intensity compared with the finite-$\hat {k}_x$ lobes. However, they do not appear as the primary in this case. The primary mode develops around $\hat {k}_y\approx 0.5$ and shows a broader structure in $\hat {k}_x$. This confirms what was found linearly (figure 9) and suggests that turbulence might have a narrower radial structure closer to marginality. This new primary mode develops in a region that is stable in simulations at higher temperature gradients.

Figure 16. Fourier spectra of time-averaged $\langle \hat {\varphi }^2_{\boldsymbol {k}}\rangle$ as a function of the normalised binormal $k_y\rho _\textrm {i}$ and radial $k_x\rho _\textrm {i}$ wave numbers for $a/L_\textrm {Ti}=0.9$ at the three key time instants: (a) first time instant (adapted from Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022); (b) second time instant; and (c) third time instant.

The time evolution of the primaries is different from the fluid-like case. Both the primary characterised by the band in $\hat {k}_x$ and the two finite-$\hat {k}_x$ lobes persist in the pre-saturation phase, whereas the isotropic structure decreases in relative intensity (figure 16b). Their time evolution is shown in figure 14(b). In particular, $\boldsymbol {k}_3\equiv (\hat {k}_x,\hat {k}_y)=(0.14,0.5)$ and $\boldsymbol {k}_4\equiv (\hat {k}_x,\hat {k}_y)=(0.71,1.05)$ are chosen as representative. We note that the $\hat {k}_x$ band shows a longer decaying time and saturates at higher values compared to the finite-$\hat {k}_x$ lobe. This feature is also visible from the Fourier spectra of fluctuations at the third time instant (figure 16c) and is attributed to a reduced efficiency of the zonal flows in suppressing the broad $\hat {k}_x$ band.

The finite-$\hat {k}_x$ lobes, present in both the fluid-like and the Floquet-type cases, are understood as modes that peak further along the magnetic field line: see figure 17. By running a simulation for $a/L_\textrm {Ti}=2.0$ with a larger $k_x\rho _\textrm {i}$ domain, we observe that the checkerboard-like structure extends even more in $k_x\rho _\textrm {i}$ (figure 18). However, the intensity of the lobes decreases with increasing $k_x\rho _\textrm {i}$ and the saturated level of fluctuations remains unchanged. We can conclude that there is no need to run larger simulations to include the whole checkerboard-like structure.

Figure 17. Electrostatic potential structure along the magnetic field line for $a/L_\textrm {Ti}=2.0$ at the binormal wave number $k_y\rho _\textrm {i}=1.05$ and two different radial wave numbers: $k_x\rho _\textrm {i}=0.04$ and $k_x\rho _\textrm {i}=1.63$.

Figure 18. Fourier spectrum of $\langle \hat {\varphi }^2_{\boldsymbol {k}}\rangle$ as a function of the binormal $k_y\rho _\textrm {i}$ and radial $k_x\rho _\textrm {i}$ wave numbers during the pre-saturation phase for $a/L_\textrm {Ti}=2.0$. Simulation with a larger range of $k_x\rho _\textrm {i}$.

At the third time instant, marginal turbulence in W7-X sets up a steep sub-$\rho _\textrm {i}$ inertial range (figure 19a). A transition at the ion Larmor radius is observed (Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022), as for $k_y\rho _\textrm {i}\lesssim 1$ (long wavelengths) the energy cascade is consistent with a $\sum _{k_x}\langle \hat {\varphi }^2_{\boldsymbol {k}}\rangle \sim (k_y\rho _\textrm {i})^{-7/3}$ scaling (Barnes, Parra & Schekochihin Reference Barnes, Parra and Schekochihin2011), whereas for $k_y\rho _\textrm {i}\gtrsim 1$ the sub-$\rho _\textrm {i}$ $\sum _{k_x}\langle \hat {\varphi }^2_{\boldsymbol {k}}\rangle \sim (k_y\rho _\textrm {i})^{-10/3}$ scaling (Schekochihin et al. Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Plunk, Quataert and Tatsuno2008) is observed. Both these scalings were observed before (Plunk, Xanthopoulos & Helander Reference Plunk, Xanthopoulos and Helander2017b) but the former was suggested to be characteristic of W7-X, and the latter of the quasi-helically symmetric stellarator HSX. Although the long-wavelength $(k_y\rho _\textrm {i})^{-7/3}$ scaling is only observed over a limited inertial range, the sub-$\rho _\textrm {i}$ $(k_y\rho _\textrm {i})^{-10/3}$ scaling persists for over a decade in $k_y\rho _\textrm {i}$, allowing a quantitative interpolation $\sum _{k_x}\langle \hat {\varphi }^2_{\boldsymbol {k}}\rangle \sim (k_y\rho _\textrm {i})^{-\alpha /3},$ yielding $\alpha =10.1\pm 0.3$. A similar transition, at ion Larmor radius scales, is observed in the radial spectrum of fluctuations (figure 19b). However, there is no explanation for such scalings at the present time.

Figure 19. Electrostatic fluctuations spectra for $a/L_\textrm {Ti}=0.9$ in (a) $k_y\rho _\textrm {i}$ and (b) $k_x\rho _\textrm {i}$ at the third time instant and compared with the predicted scalings. Adapted from Zocco et al. (Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022).

To understand the energy cascade at the third instant, it is useful to consider the primary linear source of turbulence that needs to be suppressed by zonal flows. We also study the energy cascade at later instants, when turbulence is fully developed, i.e. at $\hat {t}>400$ for $a/L_\textrm {Ti}=2.0$ and $\hat {t}>1279$ for $a/L_\textrm {Ti}=0.9$. In figure 20 we have plotted the energy spectrum for three different values of the ion temperature gradient, normalised as in Barnes et al. (Reference Barnes, Parra and Schekochihin2011). The normalisation amounts to a stretching of the inertial range where the $(k_y\rho _\textrm {i})^{-7/3}$ trend is expected to be observed. Since the normalisation factors include an $a/L_\textrm {Ti}$ factor, such stretching is larger for larger ion temperature gradients. The range over which the $(k_y\rho _\textrm {i})^{-7/3}$ scaling is observed is the largest for $a/L_\textrm {Ti}=4.0$ and it decreases with decreasing $a/L_\textrm {Ti}$. Hence the reduced inertial range in which the $(k_y\rho _\textrm {i})^{-7/3}$ scaling is observed for $a/L_\textrm {Ti}=0.9$. The transition to a sub-$\rho _\textrm {i}$ steeper scaling is still observed but it is characterised by a less-pronounced slope, which is, however, incompatible with any scaling shallower than $(k_y\rho _\textrm {i})^{-9/3}$. The observation of the sub-$\rho _\textrm {i}$ transition can thus be put in a broader context of a general steepening of turbulence spectra with decreasing ion temperature gradient (Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022).

Figure 20. Electrostatic fluctuations spectra for three different ion temperature gradients as a function of $k_y\rho _\textrm {i}$ and normalised as in Barnes et al. (Reference Barnes, Parra and Schekochihin2011). The fully developed turbulent regime is considered and compared with predicted scalings. Here $q$ is the safety factor and $\kappa =a/L_\textrm {Ti}$.

We have previously observed and suggested that turbulence close to marginality in W7-X is more radially localised and thus less prone to nonlinear zonal flows suppression (Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022). The mechanism that allows such suppression is the shearing of turbulent structures by zonal flows, which causes an up-shift of the nonlinear critical gradient with respect to the linear stability threshold. This up-shift is known as Dimits shift (Rogers, Dorland & Kotschenreuther Reference Rogers, Dorland and Kotschenreuther2000; St-Onge Reference St-Onge2017) and even though in tokamaks for realistic conditions (Mikkelsen & Dorland Reference Mikkelsen and Dorland2008) it can be less abrupt than originally observed (Dimits et al. Reference Dimits, Bateman, Beer, Cohen, Dorland, Hammett, Kim, Kinsey, Kotschenreuther, Kritz, Lao, Mandrekas, Nevins, Parker, Redd, Shumaker, Sydora and Weiland2000), it is always present. In order to explore whether this is also true in W7-X, we run further simulations where we artificially suppress zonal flows in the $\boldsymbol {E}\times \boldsymbol {B}$ nonlinearity. Results of such simulations are reported in figure 21(a–c). The zonal flows suppression causes the electrostatic fluctuations to grow faster and to saturate at larger amplitudes than before. The result is compatible with what is observed in tokamaks (Lin et al. Reference Lin, Hahm, Lee, Tang and White1998), but the effect is not as strong. In all three cases the amplitude of electrostatic fluctuations is enhanced, but when comparing the value of fluctuations case by case, we see that the effect decreases with decreasing $a/L_\textrm {Ti}$. To systematically analyse the simulations, we consider an instant just after the pre-saturation phase (identified by the dotted lines in the plots) and we calculate the ratio $r(t)=\langle \hat {\varphi }^2_{{w/o}}\rangle /\langle \hat {\varphi }^2_{{w/}}\rangle$. For the set of gradients $a/L_\textrm {Ti}=\{0.9,2.0,4.0\}$ we find that $r=\{3.3,4.7,5.1\}$. The ratio thus decreases with decreasing temperature gradient suggesting that the zonal flows are less effective when the temperature gradient is small (Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022).

Figure 21. Time evolution of the average squared modulus of the electrostatic potential fluctuations $\langle \hat {\varphi }^2\rangle$ for simulations with and without zonal flows in the $\boldsymbol {E}\times \boldsymbol {B}$ nonlinearity: (a) $a/L_\textrm {Ti}=0.9$, (b) $a/L_\textrm {Ti}=2.0$ and (c) $a/L_\textrm {Ti}=4.0$.

To assess whether this reduced effectiveness is related to the spatial distribution of the turbulence, we compute the inverse Fourier transform of the non-zonal components of $\langle \hat {\varphi }^2\rangle$. Figures 22(a) and 22(b) show the real-space structure of fluctuations, in the saturated regime, for both fluid-like and Floquet-type turbulence. Note that the $x$ axes of the two figures have different scales due to the choice of different radial wave number resolutions. We observe that fluid-like turbulence is characterised by larger radial structures, with radial dimensions of more than $50 \rho _\textrm {i}$. Going towards marginality, we observe, as predicted, narrower radial structures. The fluctuations are generally more filamented and characterised by dimensions of $\approx 25 \rho _\textrm {i}$.

Figure 22. Real-space extent of electrostatic potential fluctuations in the saturated regime for (a) $a/L_\textrm {Ti}=0.9$ and (b) $a/L_\textrm {Ti}=2.0$.

To more quantitatively characterise the radial extent of turbulent structures, we calculate the radial correlation function (Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022)

(5.1)\begin{equation} \mathcal{C} (\Delta x)=\left(\int {\rm d} t\, {\rm d}\kern0.7pt x\, {\rm d} y\, {\rm d} z|\varphi|^2\right)^{{-}1}\int {\rm d} t\, {\rm d}\kern0.7pt x\, {\rm d} y\, {\rm d} z \varphi(x)\varphi(x+\Delta x), \end{equation}

and the average radial correlation length

(5.2)\begin{equation} \bar{\Delta x}=\int_0^{\infty} \mathcal{C}(\Delta x) \,{\rm d} (\Delta x). \end{equation}

From figure 23 we see that the radial correlation function for Floquet-type turbulence decays faster with $\Delta x$ than that for fluid-like turbulence. This is reflected in the values of the average radial correlation length, which for $a/L_\textrm {Ti}=\{2.0, 0.9\}$ are $\bar {\Delta x}=\{5.0, 2.9\}$.

Figure 23. Turbulence radial correlation function in the fully developed turbulent regime for the fluid-like $(a/L_\textrm {Ti}=2.0)$ and Floquet-type turbulence $(a/L_\textrm {Ti}=0.9)$. Adapted from Zocco et al. (Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022).

These results suggest a reason for the lack of turbulence suppression when the temperature gradient is close to the linear stability threshold. Because of the small radial correlation length of Floquet-type turbulence, zonal flows are relatively ineffective, and a finite energy flux persists close to the stability threshold (Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022). This can be observed in figure 24, showing the normalised ion energy flux versus the ion temperature gradient. Floquet-type turbulence produces a normalised energy flux of $Q_\textrm {i}/Q_{gB}\approx 0.3$. In order to obtain this plot, additional nonlinear simulations were performed at $a/L_\textrm {Ti}=\{1.5,3.0,4.0\}$. Note that there is a smooth transition from strongly-driven fluid-like generated transport to the Floquet-type one. The linear threshold is also indicated, and there is no sign of a Dimits shift or a nonlinear threshold for the onset of transport in W7-X. Similar behaviour has earlier been observed in turbulence simulations of quasi-helical symmetric stellarators (Bader et al. Reference Bader, Faber, Schmitt, Anderson, Drevlak, Duff, Frerichs, Hegna, Kruger and Landreman2020; Terry et al. Reference Terry, Li, Pueschel and Whelan2021) but is very different from that in tokamaks. Figure 25 shows the contribution of different wave numbers to the turbulent energy flux for $a/L_\textrm {Ti}=0.9$. A comparison with figure 16(c) confirms they belong to the broad $\hat {k}_x$-band characteristic of Floquet-like turbulence (Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022).

Figure 24. Normalised ion energy flux $Q_\textrm {i}/Q_{gB}$ as a function of the ion temperature gradient $a/L_\textrm {Ti}$ for different nonlinear simulations. The linear threshold is also indicated. Adapted from Zocco et al. (Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022).

Figure 25. Contribution of various wave numbers to the normalised ion energy flux $Q_\textrm {i}/Q_{gB}$ in the fully developed turbulent regime for $a/L_\textrm {Ti}=0.9$. Adapted from Zocco et al. (Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022).

6. Shear effects on the ITG destabilisation

In this section we study the effects of global magnetic shear $\hat {s}=-(r/\iota ) \,\textrm {d} \iota /\textrm {d} r$ on ITG destabilisation in W7-X, with an emphasis on Floquet modes and the instability threshold. We are particularly interested in exploring the effects of a strong and positive magnetic shear. W7-X usually possesses a low, negative shear (in the tokamak sense, $\hat s > 0)$ but can acquire larger shear of either sign if a significant toroidal current is present in the plasma. Such a current arises in practice through ECCD, which tends to create a local maximum in the $\iota$-profile, with regions of significant shear on either side of the maximum (Zanini et al. Reference Zanini, Laqua, Thomsen, Stange, Brandt, Braune, Brunner, Fuchert, Hirsch and Knauer2020, Reference Zanini, Buttenschön, Laqua, Thomsen, Stange, Brandt, Braune, Brunner, Dinklage and Gao2021; Aleynikova et al. Reference Aleynikova, Hudson, Helander, Kumar, Geiger, Hirsch, Loizu, Nührenberg, Rahbarnia and Qu2021; Zocco et al. Reference Zocco, Mishchenko, Nührenberg, Könies, Kleiber, Borchardt, Slaby, Zanini, Stange, Laqua, Rahbarnia, Thomsen, Wolf, Helander, Hatzky and Cole2021). Accordingly, we proceed from an MHD equilibrium computed from the vacuum low-mirror (AIM) W7-X configuration, with prescribed profiles of density, temperature and rotational transform shown in figure 26.

Figure 26. The $\iota$ profiles for the low-mirror and ECCD non-vacuum distorted low-mirror equilibrium (in blue). In green, the normalised ion and electron temperature and density profiles used to construct the non-vacuum equilibrium.

In the magnetic geometry of this equilibrium, we run linear, flux tube, electrostatic simulations with kinetic electrons and similar resolutions to those above. Each simulation is carried out at a different radial position $r/a$, mostly in the positive-shear region. The radial scan implies a scan in the shear $\hat {s}$ and in the normalised ion temperature gradient $a/L_\textrm {Ti}$. We also change $\tau$ consistently with the prescribed plasma profiles, but we set $a/L_\textrm {Te}=a/L_\textrm {n}=0$. We compare the simulations performed in the ECCD distorted equilibrium with the ones obtained with the original low-mirror equilibrium without plasma current. We summarise in table 3 in the Appendix the radial positions where the simulations were carried out, the corresponding values of the local shear for the two equilibria ($\hat {s}_{\textrm {ECCD}}$ and $\hat {s}_{\textrm {AIM}}$), the normalised ion temperature gradient $a/L_\textrm {Ti}$ and the temperature ratio $\tau$, together with the resulting growth rates.

In figure 27 we show the normalised growth rate and frequency for the most unstable mode as a function of the ITG and radial position, for both equilibria. Note that the different points in these plots refer to different radii, with different values of the temperature gradient and magnetic shear at each point. A steepening of the growth rate curve is visible for the ECCD distorted equilibrium at radii where the magnetic shear $\hat {s}_{\textrm {ECCD}}$ is the largest ($0.45< r/a< 0.6$). This behaviour is favourable as it could produce a higher instability threshold. At smaller radii ($r/a\approx 0.4$), $\hat {s}_{\textrm {ECCD}}$ decreases and approaches values similar to those in the current-free configuration. The computed growth rates are also similar in the two configurations at these radial positions. We note that a transition to background Floquet modes is observed in both equilibria for radii at which $a/L_\textrm {Ti} < 1$. This is evidenced by a jump in the value of the real frequency and by the extended structure of the potential along the field line (see figure 28).

Figure 27. Normalised linear growth rates in (a) and normalised frequencies in (b), both as a function of the ion temperature gradient and radial position for the low-mirror configuration and ECCD distorted $\iota$ equilibrium. In (a) the radial profile of the global shear $\hat {s}$ is plotted in dashed black.

Figure 28. Comparison of the electrostatic potential and $|\nabla \alpha |^2$ structure along the magnetic field line for the low-mirror equilibrium (dashed line) and the ECCD distorted $\iota$ equilibrium (solid line) at $r/a=0.4$.

It thus appears that a large enough global magnetic shear can lead to an increase of the linear instability threshold and, in this sense, to stabilisation. For illustration, we study the structure of the electrostatic potential and flux tube geometric coefficients along the field line at two different radial positions. The first refers to the location $r/a=0.4$, where the magnitude of the shear is similar for both equilibria: $\hat {s}_{\textrm {AIM}}=-0.024$ and $\hat {s}_{\textrm {ECCD}}=-0.019$; and the second to $r/a=0.5$, where the difference between the two equilibria is substantial: $\hat {s}_{\textrm {AIM}}=-0.057$ and $\hat {s}_{\textrm {ECCD}}=0.386$. The eigenfunctions are displayed in figures 28 and 29. At the position where the shear is similar ($r/a=0.4$), the two eigenfunctions have similar shapes and extended Floquet modes are observed (figure 28). In contrast, where the ECCD-distorted $\iota$ equilibrium features a larger and positive shear ($r/a=0.5$) the eigenfunction decays much faster, resulting in stabilisation.

Figure 29. Comparison of the electrostatic potential and $|\nabla \alpha |^2$ structure along the magnetic field line for the low-mirror equilibrium (dashed line) and the ECCD distorted $\iota$ equilibrium (solid line) at $r/a=0.5$.

The main cause of the eigenfunction localisation can be found in the change of the parallel structure of

(6.1)\begin{equation} |\nabla\alpha|^2 = |\nabla\theta - \iota \boldsymbol{\nabla} \varphi|^2 - 2 \varphi \iota' \boldsymbol{\nabla} r \boldsymbol{\cdot} (\boldsymbol{\nabla} \theta - \iota \boldsymbol{\nabla} \varphi) + \varphi^2 \iota'^2 |\boldsymbol{\nabla} r|^2, \end{equation}

with increased shear. This function enters the gyrokinetic equation in the argument of the Bessel function $J_0(k_\perp \rho _\textrm {i}) = J_0(k_\alpha |\boldsymbol {\nabla } \alpha | \rho _\textrm {i})$. For non-zero shear, $\iota '(r) \ne 0$, it exhibits secular growth along the flux tube, which is dominated by the last term on the right-hand side. Indeed, at $r/a=0.5$, $|\nabla \alpha |^2$ increases roughly quadratically for large $z$, and additionally features spikes that serve to localise the eigenfunction due to finite Larmor radius (FLR) suppression (Waltz & Boozer Reference Waltz and Boozer1993; Helander et al. Reference Helander, Beidler, Bird, Drevlak, Feng, Hatzky, Jenko, Kleiber, Proll and Turkin2012). This effect has already been observed in tokamaks and in W7-X for ad-hoc constructed $|\nabla \alpha |^2$ profiles (Roberg-Clark, Plunk & Xanthopoulos Reference Roberg-Clark, Plunk and Xanthopoulos2021). At $r/a=0.4$, where the global magnetic shear is much smaller, the value of $|\nabla \alpha |^2$ is of comparable amplitude in the two equilibria.

To further assess the effect of magnetic shear on the ITG growth rate, we prescribe a different equilibrium, where the $\iota$ profile is tailored to possess strongly positive shear in the plasma region where the normalised ion temperature gradient is the largest (figure 30). We refer to this $\iota$ profile as ‘reversed’.

Figure 30. Reversed $\iota$ and low-mirror $\iota$ profiles (in blue) with the normalised ion temperature and ion temperature gradient profiles (in green).

We first run simulations of plasmas with an ion temperature gradient only, and compare them with similar simulations run in the ordinary low-mirror configuration (figure 31). We follow the same procedure as before: we run each simulation at a different radial position and change the temperature gradient and $\tau$ consistently with the plasma profiles. A summary of all the radial locations, shear, gradients and $\tau$ values, together with the resulting growth rates is given in table 4 in the Appendix. In the reversed $\iota$ equilibrium, the growth rate shows a non-monotonic trend with increasing ion temperature gradient (and $r/a$). The reduction is up to 20 % when the magnetic shear $\hat {s}_{\textrm {rev}}$ reaches its maximum value at $r/a=0.85$. Apparently, the stabilising effect of magnetic shear more than offsets the destabilising effect of a larger temperature gradient. Another interesting feature is the steepening of the curve at lower gradients. Compared with the result in the low-mirror configuration, the reversed $\iota$ one seems to show a larger linear threshold and no presence of Floquet background modes.

Figure 31. Normalised linear growth rates (in blue) as a function of the ion temperature gradient and radial position for the low-mirror configuration and reversed $\iota$ equilibrium. In dashed black the profiles of the magnetic shear for both equilibria are shown.

As before, we look at the structure of the electrostatic potential and the geometric coefficients along the magnetic field line. We study the radial position $r/a=0.85$, where the growth rate shows the maximum reduction. At this position the values of the magnetic shear in the two equilibria are $\hat {s}_{\textrm {AIM}}=-0.254$ and $\hat {s}_{\textrm {rev}}=3.450$, respectively. We then consider a case where Floquet modes are observed in the low-mirror equilibrium but not in the reversed $\iota$ one, namely $r/a=0.4$. The values of the magnetic shear are $\hat {s}_{\textrm {AIM}}=-0.024$ and $\hat {s}_{\textrm {rev}}=-0.210$ here.

In the first case ($r/a=0.85$), the eigenfunctions are very localised for both equilibria (figure 32). The gradient is large and the mode is fluid-like. When looking at the structure of $|\nabla \alpha |^2$, however, we note a big difference between the two equilibria. The reversed $\iota$ one features peaks that reach amplitudes two orders of magnitude higher than the low-mirror one. In this case, $|\nabla \alpha |^2$ plays no role in the localisation of the mode but rather in its stabilisation, possibly still through FLR suppression. Similarly, the magnetic drift $\kappa _2$ exhibits differences both in magnitude and in location. This could also play a role both in the localisation and stabilisation of the mode. Note the two different scales in the plots for the two equilibria.

Figure 32. Comparison of the electrostatic potential, the quantity $|\nabla \alpha |^2$ on a semi-logarithmic scale, and $\kappa _2$ structure along the magnetic field line for the low-mirror equilibrium (dashed line) and the reversed $\iota$ equilibrium (solid line) at $r/a=0.85$.

In the second case ($r/a=0.4$), the low-mirror equilibrium features a Floquet mode, extended along the field line, while the reversed $\iota$ one shows a more localised mode (figure 33). This explains the absence of a foot in the growth rate for the latter case. When looking at the structure of $|\nabla \alpha |^2$, we note that we can again ascribe the localisation of the eigenfunction to an increase in the magnitude of $|\nabla \alpha |^2$ at the extremes of the flux tube. It is important to note that in this case the sign of the shear is negative, but it is nevertheless one order of magnitude larger in magnitude than that of the low-mirror configuration. We thus conclude that it is the magnetic shear magnitude that plays a role in the increase of the linear instability threshold, not its sign. Finally, the structure of the magnetic drift $\kappa _2$ also shows differences, even if less pronounced than in the case studied previously.

Figure 33. Comparison of the electrostatic potential, the quantity $|\nabla \alpha |^2$, and $\kappa _2$ structure along the magnetic field line for the low-mirror equilibrium (dashed line) and the reversed $\iota$ equilibrium (solid line) at $r/a=0.4$.

An inconsistent $a/L_\textrm {Ti}$ scan performed at the just discussed radial positions confirms the predominant role of the shear in the modification of the ITG growth rate in figure 31. At $r/a=0.4$ (figure 34a) in the reversed $\iota$ equilibrium, we observe a near-threshold steepening of the growth rate curve due to Floquet modes stabilisation (figure 33). At $r/a=0.85$ (figure 34b) both the toroidal ITG branch and the Floquet ones are stabilised. In particular, the latter is damped completely, yielding negative growth rate values for $a/L_\textrm {Ti}<1.5$.

Figure 34. Normalised linear growth rates for the low-mirror configuration and reversed $\iota$ equilibrium at (a) $r/a=0.4$ and (b) $r/a=0.85$. Featuring a near-marginality growth rate steepening in (a) and a total suppression of Floquet modes in (b).

We conclude that increasing the absolute value of the magnetic shear could be a way to affect both the far- and near-threshold behaviour of ITG modes. In particular, it seems to be an effective way to suppress the Floquet modes, otherwise always present, and thus to increase the linear instability threshold in low-shear devices.

We note that the favourable instability threshold increase due to steepening of the growth rate trend can come with the price of destabilisation at other radial positions. This happens where the distorted $\iota$ equilibria feature similar or smaller magnetic shear magnitude to the low-mirror configuration (see both figures 27 and 31). In the following nonlinear analysis, we tackle the effect this destabilisation can have on transport.

In order to understand the effect of a distorted iota profile on transport, we run nonlinear simulations for the AIM and reversed $\iota$ equilibria. We take two representative cases: one where the magnitude of the shear in the reversed $\iota$ equilibrium is the largest ($r/a=0.85$ with $\hat {s}_{\textrm {rev}}=3.450$) and one where it is the smallest ($r/a=0.75$ with $\hat {s}_{\textrm {rev}}=-0.035$). In the latter, in particular, the reversed $\iota$ equilibrium exhibits a larger linear growth rate than in the low-mirror configuration and we are interested in observing whether this destabilisation is reflected also nonlinearly. For the same radial positions, the low-mirror configuration has $\hat {s}_{\textrm {AIM}}=-0.171$ for $r/a=0.75$ and $\hat {s}_{\textrm {AIM}}=-0.254$ for $r/a=0.85$. In the two radial positions, $a/L_\textrm {Ti}=2.49$ and $a/L_\textrm {Ti}=5.60$ for $r/a=0.75$ and $r/a=0.85$, respectively. We expect to observe fluid-like turbulence in both cases and for this reason we choose to run the simulations with the resolutions reported in the first row of table 1. We set $a/L_\textrm {n}=a/L_{T_\textrm {e}}=0$ and $\tau$ consistent with plasma profiles, as reported in table 4 in the Appendix.

We report the average saturated ion energy flux in table 2. At the radius where the reversed $\iota$ equilibrium features a large, positive shear ($r/a=0.85$), the ion energy flux is about eight times smaller than in the low-mirror configuration. The result confirms the linearly observed reduction. The linear trend is also confirmed at $r/a=0.75$, where a larger growth rate in the reversed $\iota$ configuration reflects in a larger ion energy flux. The difference with the low-mirror configuration is, however, less than a factor of 1.5. The desired maximisation of the stabilising shear effect where the gradients are the steepest is thus confirmed also nonlinearly. It is plausible to conjecture that the mild linear destabilisation at around $r/a=0.75$ is not as important as the larger stabilising effect just discussed.

Table 2. Resulting average ion energy flux $Q_\textrm {i}/Q_{gB}$ from nonlinear simulations in the low-mirror (AIM) and reversed $\iota$ configurations at two different radial locations.

Finally, we take into account the effects of finite density and electron temperature gradients (figure 35). We switch them on one by one and re-run the linear simulations at the radial locations where the shear effect is the strongest. All the gradients are changed at every radial location, consistently with plasma profiles. A summary of all the radial locations, shear, gradients and $\tau$ values, together with the resulting growth rates is given in table 5 in the Appendix. As noted in the diagrams in the first section and in previous works (Helander et al. Reference Helander, Bird, Jenko, Kleiber, Plunk, Proll, Riemann and Xanthopoulos2015; Alcusón et al. Reference Alcusón, Xanthopoulos, Plunk, Helander, Wilms, Turkin, von Stechow and Grulke2020), the addition of a non-zero density gradient in W7-X leads to a reduction in the growth rate at the radial locations where the shear is the strongest. If a finite electron temperature gradient is also added, the growth rate curve is shifted upwards, i.e. the plasma becomes more unstable. A non-monotonic trend of the growth rate is, however, still observed thanks to the stabilising action of a large shear. It is worth mentioning that only ion length scales are included in these simulations and the growth rates reported in figure 35 all refer to fluid-like ITG modes. We conclude that, also in the most realistic case, the tweaking of the $\iota$ profile seems to reduce the ITG instability, although in a less pronounced fashion.

Figure 35. Normalised linear growth rates (in blue) as a function of the ion temperature gradient and radial position for the reversed $\iota$ equilibrium with the effects of a finite density and electron temperature gradient.