1. Introduction
Stellarators and tokamaks are promising devices for achieving magnetically confined thermonuclear fusion. Poorly designed stellarators suffer from trapped-particle losses. As a result, they feature high levels of collisional (neoclassical) transport, significant power losses and poor confinement. These problems can be mitigated by optimising the geometry of the magnetic field (Nührenberg & Zille Reference Nührenberg and Zille1986; Boozer Reference Boozer1995). A promising line of optimisation aims at achieving quasi-isodynamic magnetic fields (Nührenberg Reference Nührenberg2010). Stellarators with this type of magnetic geometry are characterised by particle orbits with very small radial drifts, and thus low neoclassical transport.
Wendelstein 7-X (W7-X) is the first large, superconducting, optimised stellarator. The first experimental campaign proved the success of magnetic field optimisation. Indeed, W7-X shows a reduction of neoclassical transport effects and features among the highest energy confinement times for stellarators (Dinklage et al. Reference Dinklage, Beidler, Helander, Fuchert, Maaßberg, Rahbarnia, Sunn Pedersen, Turkin, Wolf and Alonso2018; Bozhenkov et al. Reference Bozhenkov, Kazakov, Ford, Beurskens, Alcusón, Alonso, Baldzuhn, Brandt, Brunner and Damm2020; Beidler et al. Reference Beidler2021). However, power balance analyses highlight that the calculated remaining neoclassical transport is not sufficient to explain the observed energy losses. W7-X transport is found to be mainly turbulent and is thought to be dominated by ion-scale turbulence caused by microinstabilities driven by temperature and density gradients in the plasma (Bozhenkov et al. Reference Bozhenkov, Kazakov, Ford, Beurskens, Alcusón, Alonso, Baldzuhn, Brandt, Brunner and Damm2020).
As in tokamaks, the microinstability which is thought to be most harmful to W7-X confinement is the one driven by the ion temperature gradient (ITG) (Carralero et al. Reference Carralero, Estrada, Maragkoudakis, Windisch, Alonso, Beurskens, Bozhenkov, Calvo, Damm and Ford2021). Joint experimental and theoretical studies have highlighted how effects such as the exacerbation of the turbulent heat transport with increasing electron heating can be explained by the presence of ITG-driven turbulence (Beurskens et al. Reference Beurskens, Bozhenkov, Ford, Xanthopoulos, Zocco, Turkin, Alonso, Beidler, Calvo and Carralero2021). Decreasing the ion-to-electron temperature ratio $\tau =T_{\rm i}/T_{\rm e}$ at constant ion temperature $T_{\rm i}$ lowers the ITG destabilisation threshold (Zocco et al. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018) and increases the turbulent transport. A similar effect has also been observed in tokamaks (Romanelli Reference Romanelli1989; Xu & Rosenbluth Reference Xu and Rosenbluth1991; Casati et al. Reference Casati, Bourdelle, Garbet and Imbeaux2008). For W7-X, the importance of this effect has been confirmed by gyrokinetic simulations (Zocco et al. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018; Beurskens et al. Reference Beurskens, Bozhenkov, Ford, Xanthopoulos, Zocco, Turkin, Alonso, Beidler, Calvo and Carralero2021) performed with the code GENE (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000), and its influence on the ion-to-electron thermal coupling at transport scales has been investigated within a multiscale approach, providing a physical interpretation for the relatively low ion temperatures at high electron heating power (Navarro et al. Reference Navarro, Siena, Velasco, Wilms, Merlo, Windisch, LoDestro, Parker and Jenko2023).
The impact of the temperature ratio, $\tau$, on ITG modes can be studied in some detail close to the linear stability threshold. In this little-explored regime, it has been observed that the usual toroidal, fluid-like, branch of the ITG mode is replaced by a different type of ITG instability, of the Floquet-type (Bhattacharjee et al. Reference Bhattacharjee, Sedlak, Similon, Rosenbluth and Ross1983; Connor & Taylor Reference Connor and Taylor1987; Candy, Waltz & Rosenbluth Reference Candy, Waltz and Rosenbluth2004; Zocco et al. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018). In this regime, the instability is characterised by an extended structure along the magnetic field lines. For this reason, in order for these modes to be observed in simulations, the parallel integration domain needs to be extended enough. As a result, these modes were overlooked in virtually all flux tube stellarator simulations prior to the work in Zocco et al. (Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018).
Despite being of crucial relevance for transport, ITG turbulence in stellarators has rarely been studied in detail near its marginal stability threshold because of the relatively high computational requirements necessary for obtaining well-converged results. Nonlinear gyrokinetic simulations performed with the code stella (Barnes, Parra & Landreman Reference Barnes, Parra and Landreman2019), however, show how such modes also cause a finite level of turbulent heat transport, if enough care is taken in resolving their extended parallel structure (Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022).
In Zocco et al. (Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022), the effects of density gradients on the nearly marginally stable ITG mode were neglected since density profiles in W7-X are often rather flat. Pellet injection, however, can be used to achieve higher densities and density gradients in high-performance discharges (Bozhenkov et al. Reference Bozhenkov, Kazakov, Ford, Beurskens, Alcusón, Alonso, Baldzuhn, Brandt, Brunner and Damm2020). These discharges show higher core ion temperatures and reduced turbulent fluctuations (Carralero et al. Reference Carralero, Estrada, Maragkoudakis, Windisch, Alonso, Beurskens, Bozhenkov, Calvo, Damm and Ford2021). This behaviour was captured also through linear gyrokinetic simulations (Alcusón et al. Reference Alcusón, Xanthopoulos, Plunk, Helander, Wilms, Turkin, von Stechow and Grulke2020). Here, we investigate previously partially explored regimes, producing high-fidelity near-marginal results, shedding light on the wavelength dependence of the W7-X stability diagram, with particular emphasis on transport. Our highly resolved studies show that ITG and universal modes (Coppi & Pegoraro Reference Coppi and Pegoraro1977; Smolyakov, Yagi & Kishimoto Reference Smolyakov, Yagi and Kishimoto2002; Helander & Plunk Reference Helander and Plunk2015; Landreman, Antonsen & Dorland Reference Landreman, Antonsen and Dorland2015) can be somewhat less unstable than ion-driven trapped-electron modes, but generally have longer wavelengths. A quasi-linear estimate of the turbulent heat flux suggests they are important for transport although they were overlooked in previously published stability diagrams (Alcusón et al. Reference Alcusón, Xanthopoulos, Plunk, Helander, Wilms, Turkin, von Stechow and Grulke2020). Ion-scale dynamics can, however, be affected by the presence of electron temperature gradients.Footnote 1
We successively focus on pure ITG instabilities with ion temperature gradients, but neglect electron temperature gradients. We investigate the parallel structure of the mode along the magnetic field line, and perform convergence studies to identify the most suitable input parameters for nonlinear simulations. A discussion follows on the numerical limits that make near-marginal studies difficult and why some compromises have to be made.
A full spectral analysis of fluctuations in the relevant regimes is presented, emphasising the difference between fluid-like and Floquet-like turbulence, especially during the pre-saturation and saturation phases. Differences are striking, and the regime identification is unambiguous. The study of Fourier spectra of fluctuations suggests that Floquet-type turbulence is more radially localised, and we speculate this is the reason why zonal flows are not efficient in shearing and suppressing turbulent structures (Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022). To test this hypothesis, a numerical experiment is performed, where the zonal-flow contribution is switched off and the saturated amplitude of fluctuations is studied. We thus extend the previous work of Zocco et al. (Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022).
Finally, we further investigate the effect of the magnetic field geometry on the ITG linear threshold and transport levels. For this purpose we distort the radial profile of the rotational transform $\iota$ in accordance with a prescribed analytical function that exhibits strong variations, and recalculate the magnetic equilibria. We consider, in particular, the effect of electron cyclotron current drive (ECCD) (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) and observe indications of a possible increase of the stability threshold when $\iota$ is modified. It is found that the presence of strong magnetic shear regions in distorted $\iota$ profiles mitigates ITG modes, linearly and nonlinearly. More importantly, extended modes are successfully suppressed in this equilibria, leading to an increased ITG linear critical gradient.
The structure of this article is as follows. In § 2 we present the equations and numerical framework, and in § 3 near-marginal linear and quasi-linear stability diagrams are presented. In §§ 4 and 5, the linear and nonlinear results for the transition from fluid-like to Floquet-type ITG instabilities and turbulence are presented, respectively. In § 6, the stabilising effects of $\iota$-profile distortions are studied through linear and nonlinear simulations performed in a non-vacuum magnetic equilibrium. In § 7 we summarise our results and conclusions.
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:
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
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
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}\}$.
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.
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.
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.
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.
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)
where $\hat {\gamma }=\gamma a/v_\textrm {thi}$ is taken from figure 2(a–c) and
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).
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).
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.
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$.
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.
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.
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.
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.
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.
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.
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).
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).
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}$.
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)
and the average radial correlation length
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\}$.
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).
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.
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).
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.
The main cause of the eigenfunction localisation can be found in the change of the parallel structure of
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’.
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.
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.
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.
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$.
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.
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.
7. Discussion and conclusions
We have characterised electrostatic gyrokinetic instabilities and turbulence in W7-X, with particular attention given to their character close to the linear stability threshold. The presence and character of ion-temperature- and density-gradient-driven instabilities can conveniently be summarised in stability diagrams. In contrast to a previous work (Alcusón et al. Reference Alcusón, Xanthopoulos, Plunk, Helander, Wilms, Turkin, von Stechow and Grulke2020), we have constructed different diagrams showing instabilities with linear growth rates at different spatial scales and rotating in different directions. This method ensures that no ion-temperature- or density-gradient-driven instability is overlooked, especially when it is not the most unstable but it is destabilised at scales that are more relevant for transport. An important example is the observation of universal modes (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) for large enough density gradients. These modes were previously observed numerically in W7-X simulations (Costello et al. Reference Costello, Proll, Plunk, Pueschel and Alcusón2023), but were thought to be overwhelmed by faster growing instabilities in the presence of an ion temperature gradient. Through the stability diagrams we have here shown how these modes are stabilised by the presence of an ion temperature gradient but are always present, despite being slightly more stable than ITG modes or iTEMs. More importantly, they are characterised by wavelengths longer than the ion Larmor radius and a quasi-linear estimate of the turbulent heat flux suggests they are relevant for transport in W7-X.
We have paid particular attention to the geometric nature of ITG instabilities close to the linear stability threshold. For small ion temperature gradients, the dominant ITG modes are of the Floquet-type (Bhattacharjee et al. Reference Bhattacharjee, Sedlak, Similon, Rosenbluth and Ross1983; Connor & Taylor Reference Connor and Taylor1987; Candy et al. Reference Candy, Waltz and Rosenbluth2004; Zocco et al. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018) and produce small but non-zero energy transport (Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022). These modes are greatly extended along field lines, particularly close to the linear instability threshold. Long flux tubes are thus necessary in numerical simulations to resolve the slowly decaying mode structure along the magnetic field. Nonlinearly, these modes behave differently from fluid-like modes. In particular, they are more radially localised, which may explain why zonal flows are not efficient in suppressing Floquet-type turbulence (Zocco et al. Reference Zocco, Podavini, Garcìa-Regaña, Barnes, Parra, Mishchenko and Helander2022). A range of diagnostics and numerical experiments corroborate this hypothesis. We have also observed a general steepening of fluctuation spectra with decreasing ion temperature gradient, which might be the sign of a transition to the sub-$\rho _\textrm {i}$ gyrokinetic turbulent regime predicted by Schekochihin et al. (Reference Schekochihin, Cowley, Dorland, Hammett, Howes, Plunk, Quataert and Tatsuno2008).
There is good reason to believe that extended, Floquet-type modes are ubiquitous for small magnetic shear. Indeed, instabilities that extend far along the magnetic field have also been observed in tokamak simulations, both with adiabatic (Zocco et al. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018) and kinetic electrons (Volčokas, Ball & Brunner Reference Volčokas, Ball and Brunner2022). We have observed them in near-marginal ITG numerical simulations with kinetic electrons performed in different W7-X equilibria: the high-mirror (KJM) and low-mirror (AIM) magnetic configurations, but also in modified W7-X magnetic equilibria.
W7-X is inherently a low-shear device and thus prone to the formation of such elongated modes, but the global magnetic shear can be locally enhanced by ECCD (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). The effect of the driven toroidal current is a distortion of the $\iota$ profile, which can be modelled by tailoring the $\iota$-profile of the magnetic equilibrium whilst keeping the density and temperature profiles fixed. The presence of large magnetic shear over a carefully chosen plasma region causes the growth rate to depend non-monotonically on radius, even in the most realistic scenarios where both the finite density and electron temperature gradients are taken into account. The beneficial effect is not only linear. An equilibrium with a strong, positive shear can show locally ion energy fluxes eight times smaller than an equilibrium with a non-distorted $\iota$-profile. In such tailored magnetic equilibria, the dependence of the instability growth rate on the temperature gradient suggests an increase in the linear stability threshold and, indeed, extended modes can be successfully suppressed. Manipulation of the $\iota$ profile can thus be exploited as a tool for ITG stabilisation in W7-X. In particular, the small but finite level of transport caused by Floquet-type turbulence could be reduced and the ITG critical gradient increased.
Acknowledgements
Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Commission. Neither the European Union nor the European Commission can be held responsible for them. The stella simulations have been performed on the HPCs Draco, Raven (Germany) and Marconi (Italy).
Editor Peter Catto thanks the referees for their advice in evaluating this article.
Funding
This work has been carried out within the framework of the EUROfusion Consortium, funded by the European Union via the Euratom Research and Training Programme (Grant Agreement Number 101052200 – EUROfusion).
Declaration of interest
The authors report no conflict of interest.
Appendix
Tables 3–5 present sets of physical parameters for the linear simulations performed in the ECCD distorted equilibrium and the low-mirror (AIM) equilibrium, the reversed $\iota$ equilibrium and the low-mirror (AIM) equilibrium with $a/L_{T_\textrm {i}}$ only, and the reversed $\iota$ equilibrium with $a/L_\textrm {Ti}$, $a/L_\textrm {n}$ and $a/L_\textrm {Te}$, respectively.