2.1 Flux-tube geometry

Gyrokinetic simulations of local plasma turbulence utilize the flux-tube approach (Beer et al. Reference Beer, Cowley and Hammett1995), where the gyrokinetic equations are solved in a small domain around a magnetic field line. This approach takes advantage of the anisotropic nature of plasma turbulence, where perpendicular wavelengths of fluctuations are generally much smaller than parallel wavelengths, $k_{\Vert }\ll k_{\bot }$ . Only a small domain perpendicular to the field line, of the order of a few perpendicular correlation lengths, need be simulated with periodic boundary conditions. The scale $L_{\text{eq}}$ of perpendicular variations of background quantities such as the equilibrium magnetic field is assumed to be much larger than turbulent correlations lengths $\unicode[STIX]{x1D706}_{\text{turb}}$ . Equilibrium quantities are then expanded to first order in the perpendicular direction around the centre of the simulation domain while still incorporating non-trivial dependence along the parallel field-line coordinate.

The flux-tube approach is related to the ballooning transformation formulated for axisymmetric geometries (Connor, Hastie & Taylor Reference Connor, Hastie and Taylor1978, Reference Connor, Hastie and Taylor1979) and subsequently extended to general three-dimensional systems (Dewar & Glasser Reference Dewar and Glasser1983), which transforms a multi-dimensional eigenvalue calculation into a sequence of one-dimensional eigenvalue equations along field lines. This transformation is facilitated by the assumption $k_{\Vert }\ll k_{\bot }$ for ballooning modes, allowing one to assume a Wentzel–Kramers–Brillouin-like solution where the slow scale only appears in the mode amplitude: $\unicode[STIX]{x1D709}(\boldsymbol{x},t)=\hat{\unicode[STIX]{x1D709}}(\boldsymbol{x},\unicode[STIX]{x1D716})\exp (\text{i}S(\boldsymbol{x})/\unicode[STIX]{x1D716}-\text{i}\unicode[STIX]{x1D714}t)$ , where $\unicode[STIX]{x1D709}$ is the plasma displacement. $\unicode[STIX]{x1D716}\ll 1$ is an expansion parameter and denotes the rapid variation of the wave phase perpendicular to the magnetic field line. At lowest order in  $\unicode[STIX]{x1D716}$ , the instabilities are approximately incompressible, $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\unicode[STIX]{x1D743}\approx (\text{i}\unicode[STIX]{x1D735}S(\boldsymbol{x})/\unicode[STIX]{x1D716})\boldsymbol{\cdot }\unicode[STIX]{x1D743}=0$ . We can identify $\unicode[STIX]{x1D735}S(\boldsymbol{x})$ as the wavevector  $\boldsymbol{k}$ . Requiring the displacement at lowest order to be perpendicular to the magnetic field gives $\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{B}=0$ . We may write the magnetic field in Clebsch form,

(2.1) $$\begin{eqnarray}\boldsymbol{B}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D6FC},\end{eqnarray}$$

where $\unicode[STIX]{x1D713}$ is the poloidal magnetic flux function and $\unicode[STIX]{x1D6FC}=q(\unicode[STIX]{x1D713})\unicode[STIX]{x1D703}-\unicode[STIX]{x1D701}$ is the field-line label defined in terms of the straight field-line poloidal and toroidal angles $\unicode[STIX]{x1D703}$ and  $\unicode[STIX]{x1D701}$ , respectively, and $q(\unicode[STIX]{x1D713})$ is the safety factor. The condition $\boldsymbol{k}_{\bot }\boldsymbol{\cdot }\boldsymbol{B}=0$ implies $S(\boldsymbol{x})$ is a function of $(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC})$ alone, yielding an instructive form for  $\boldsymbol{k}_{\bot }$ :

(2.2) $$\begin{eqnarray}\boldsymbol{k}_{\bot }=k_{\unicode[STIX]{x1D6FC}}[q\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}-\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}+q^{\prime }(\unicode[STIX]{x1D703}+\unicode[STIX]{x1D703}_{k})\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}].\end{eqnarray}$$

Here, $q^{\prime }=\text{d}q/\text{d}\unicode[STIX]{x1D713}$ is the flux-tube-averaged global magnetic shear and $\unicode[STIX]{x1D703}_{k}=k_{\unicode[STIX]{x1D713}}/q^{\prime }k_{\unicode[STIX]{x1D6FC}}$ is the ballooning angle. This implies that radial variations (finite- $k_{\unicode[STIX]{x1D713}}$ effects) are transformed by global magnetic shear into a dependence along the field line. Furthermore, it is easy to see that (2.2) is invariant under the substitution $\unicode[STIX]{x1D703}\rightarrow 2\unicode[STIX]{x03C0}M$ , $\unicode[STIX]{x1D703}_{k}\rightarrow -2\unicode[STIX]{x03C0}M$ for some integer  $M$ . Thus, in principle, one may compute effects for increasing $\unicode[STIX]{x1D703}$ by considering the geometry from a limited range of  $\unicode[STIX]{x1D703}$ , say $\unicode[STIX]{x1D703}\in (-\unicode[STIX]{x03C0},\unicode[STIX]{x03C0})$ and incrementing $k_{\unicode[STIX]{x1D713}}\rightarrow -2\unicode[STIX]{x03C0}Mq^{\prime }k_{\unicode[STIX]{x1D6FC}}$ for integer  $M$ .

The flux-tube simulation domain is defined in straight field-line (SFL) coordinates $(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},z)$ . Depending on whether one chooses $\unicode[STIX]{x1D6FC}=q\unicode[STIX]{x1D703}-\unicode[STIX]{x1D701}$ or  $\unicode[STIX]{x1D713}$ is the poloidal or toroidal flux label, respectively, and the rotational transform  is related to $q$ by  For subsequent discussion and to be consistent with the formalism employed by the gyrokinetic code Gene (Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000), we will choose $\unicode[STIX]{x1D713}$ to be the poloidal flux and $\unicode[STIX]{x1D6FC}=q(\unicode[STIX]{x1D713})\unicode[STIX]{x1D703}-\unicode[STIX]{x1D701}$ . The parallel coordinate $z$ is chosen to be identical to the straight field-line poloidal angle  $\unicode[STIX]{x1D703}$ . In these coordinates, the flux-tube simulation domain is defined to be a small rectangular domain centred around $\unicode[STIX]{x1D713}_{0}$ and $\unicode[STIX]{x1D6FC}_{0}$ in SFL coordinates, given by

(2.3) $$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D713}_{0}-\unicode[STIX]{x0394}\unicode[STIX]{x1D713}\leqslant \unicode[STIX]{x1D713}<\unicode[STIX]{x1D713}_{0}+\unicode[STIX]{x0394}\unicode[STIX]{x1D713},\\ \unicode[STIX]{x1D6FC}_{0}-\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}\leqslant \unicode[STIX]{x1D6FC}<\unicode[STIX]{x1D6FC}_{0}+\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC},\\ -z_{0}/2\leqslant z<z_{0}/2.\end{array}\right\}\end{eqnarray}$$

The flux-tube approach assumes periodic boundary conditions in the perpendicular directions, however this is an assumption of statistical, not physical, periodicity; that is, the statistical properties of the fluctuations at $\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D713}+2\unicode[STIX]{x0394}\unicode[STIX]{x1D713}$ while holding $\unicode[STIX]{x1D6FC}$ fixed, or at $\unicode[STIX]{x1D6FC}$ and $\unicode[STIX]{x1D6FC}+2\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}$ while holding $\unicode[STIX]{x1D713}$ fixed, are the same. Provided the domain dimensions are larger than physical correlation lengths of the fluctuations, which are typically of the order of tens of gyroradii in the radial direction, periodic boundary conditions are an efficient means to represent turbulent eddies entering and leaving the domain. A benefit of prescribing periodic perpendicular boundary conditions is that quantities can be decomposed in a Fourier basis in $\unicode[STIX]{x1D713}$ and  $\unicode[STIX]{x1D6FC}$ :

(2.4) $$\begin{eqnarray}f(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},z,t)=\mathop{\sum }_{m=-\infty }^{\infty }\mathop{\sum }_{n=-\infty }^{\infty }\hat{f}_{m,n}(z,t)\exp [\text{i}m\unicode[STIX]{x03C0}(\unicode[STIX]{x1D713}-\unicode[STIX]{x1D713}_{0})/\unicode[STIX]{x0394}\unicode[STIX]{x1D713}+\text{i}n\unicode[STIX]{x03C0}(\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FC}_{0})/\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}],\end{eqnarray}$$

for each physical quantity $f$ .

A more subtle treatment is required for the parallel boundary condition as enforcing strictly periodic boundary conditions in the parallel direction, $f(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},-z_{0},t)=f(\unicode[STIX]{x1D713},\unicode[STIX]{x1D6FC},z_{0},t)$ , leads to rational field lines, which is inaccurate for toroidal magnetic confinement configurations. Like the perpendicular boundary conditions, the parallel boundary condition should be an expression of the statistical properties of the fluctuations, which physically should be statistically invariant at the same poloidal angle along the field line. In order for this to successfully be satisfied, the parallel simulation domain must be longer than the parallel correlation length of the fluctuation. This is accomplished by applying periodic boundary conditions at different $\unicode[STIX]{x1D703}$ values, while holding $\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D701}$ fixed, rather than holding $\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D6FC}$ fixed, which leads to every field line being rational. The detailed derivation of the parallel boundary condition is given in appendix A, giving the final form as

(2.5) $$\begin{eqnarray}\displaystyle & \hat{f}_{m-\unicode[STIX]{x1D6FF}m,n}(\unicode[STIX]{x1D703}+2\unicode[STIX]{x03C0}N,t)C_{n}=\hat{f}_{m,n}(\unicode[STIX]{x1D703},t), & \displaystyle\end{eqnarray}$$

(2.6) $$\begin{eqnarray}\displaystyle & \unicode[STIX]{x1D6FF}m=nM,\quad M=2\unicode[STIX]{x03C0}Nq^{\prime }{\displaystyle \frac{\unicode[STIX]{x0394}\unicode[STIX]{x1D713}}{\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}}}. & \displaystyle\end{eqnarray}$$

Equation (2.5) expresses that the amplitude of a mode with index $m$ at one end of the parallel simulation domain is coupled to a mode with index $m-\unicode[STIX]{x1D6FF}m$ at the other end of the parallel simulation domain, with the coupling dependent on the global magnetic shear,  $q^{\prime }$ . Equation (2.5) contains the phase factor $C_{n}=\exp (\text{i}2\unicode[STIX]{x03C0}Nnq_{0}/\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC})$ . For further use in Gene, the following field-aligned coordinate system is used:

(2.7a-c ) $$\begin{eqnarray}x=\frac{q_{0}}{B_{0}r_{0}}(\unicode[STIX]{x1D713}-\unicode[STIX]{x1D713}_{0}),\quad y=\frac{r_{0}}{q_{0}}(\unicode[STIX]{x1D6FC}-\unicode[STIX]{x1D6FC}_{0}),\quad z=\unicode[STIX]{x1D703},\end{eqnarray}$$

where $r_{0}$ is the geometrical radius of the magnetic surface at the centre of the simulation domain, $\unicode[STIX]{x1D713}_{0}=\unicode[STIX]{x1D713}(r_{0})$ , with $q_{0}$ the safety factor, $B_{0}$ the magnetic field and $\unicode[STIX]{x1D6FC}_{0}$ the field-line label at the centre of the simulation domain. Similarly, the radial wavenumber $k_{\unicode[STIX]{x1D713}}$ and binormal wavenumber $k_{\unicode[STIX]{x1D6FC}}$ can be redefined as

(2.8a,b ) $$\begin{eqnarray}k_{x}=\frac{m\unicode[STIX]{x03C0}B_{0}r_{0}}{q_{0}\unicode[STIX]{x0394}\unicode[STIX]{x1D713}},\quad k_{y}=\frac{n\unicode[STIX]{x03C0}q_{0}}{r_{0}\unicode[STIX]{x0394}\unicode[STIX]{x1D6FC}}.\end{eqnarray}$$

Finally, the perpendicular coordinates $x$ and $y$ are normalized to a reference gyroradius length scale: $\hat{x}=x/\unicode[STIX]{x1D70C}_{\text{ref}}$ , ${\hat{y}}=y/\unicode[STIX]{x1D70C}_{\text{ref}}$ , yielding $\hat{k}_{x}=k_{x}\unicode[STIX]{x1D70C}_{\text{ref}}$ and $\hat{k}_{y}=k_{y}\unicode[STIX]{x1D70C}_{\text{ref}}$ . For discussion of the following results, the normalized coordinates will be used, however the over hat symbols will be dropped. The local magnetic shear is defined by the vector relation

(2.9) $$\begin{eqnarray}s_{\text{loc}}(z)=\hat{\boldsymbol{b}}(z)\times \unicode[STIX]{x1D735}\hat{\boldsymbol{n}}(z)\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times (\hat{\boldsymbol{b}}(z)\times \unicode[STIX]{x1D735}\hat{\boldsymbol{n}}(z)),\end{eqnarray}$$

where $\hat{\boldsymbol{b}}$ is the unit vector in the direction of the magnetic field and $\hat{\boldsymbol{n}}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}(z)/|\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}(z)|$ . Due to the complicated magnetic field structures in stellarators, equation (2.9) has non-trivial dependence along field lines. The global magnetic shear is defined as the flux surface average of (2.9) and can be shown to be equal to ${\hat{s}}=(r_{0}/q_{0})\,\text{d}q/\text{d}r$ .

2.2 Parallel correlations

Drift-wave instabilities in flux-tube simulations typically have ballooning nature, thus $\boldsymbol{k}_{\bot }$ takes on the same form as (2.2). Instead of solving the gyrokinetic equations over a long parallel domain, one may make use of this form of $\boldsymbol{k}_{\bot }$ and the parallel boundary condition to construct fluctuations that extend along field lines by adding more radial modes, which takes advantage of fast Fourier transform techniques (Dorland et al. Reference Dorland, Jenko, Kotschenreuther and Rogers2000; Jenko et al. Reference Jenko, Dorland, Kotschenreuther and Rogers2000; Peeters et al. Reference Peeters, Camenen, Casson, Hornsby, Snodin, Strintzi and Szepesi2009). The approach is well suited for axisymmetric configurations. After one poloidal turn ( $\unicode[STIX]{x1D703}\rightarrow \unicode[STIX]{x1D703}+2\unicode[STIX]{x03C0}$ ), all geometric quantities are periodic, even though the flux-tube ends generally lie at different $\unicode[STIX]{x1D701}$ values. This is guaranteed by axisymmetry, as $\unicode[STIX]{x1D701}$ is an ignorable coordinate and the flux tube is parameterized only by  $\unicode[STIX]{x1D703}$ . Thus, one may use the geometric data from one poloidal turn and ‘build’ an extended flux tube by adding more radial modes. A geometrically accurate computational domain can then be constructed for fluctuations that have arbitrarily large parallel correlation lengths, as one may always add a sufficient number of radial modes to resolve parallel scales.

While often yielding reasonably accurate results, this procedure is formally incorrect for most stellarator simulations. In general, geometric quantities for stellarators are not periodic after one poloidal turn, which can have serious consequences for stellarator flux-tube simulations. As a fluctuation extends along a field line in a stellarator, it samples different geometry depending on the field-line label. The exact geometry may be important in setting the physical parallel correlation length. Due to the complicated interplay of magnetic trapping regions, curvature drive and local magnetic shear, it has been observed that drift-wave instabilities in stellarators may be most unstable in regions away from $\unicode[STIX]{x1D703}_{k}=0$ at the outboard midplane (Merz Reference Merz2008; Faber et al. Reference Faber, Pueschel, Proll, Xanthopoulos, Terry, Hegna, Weir, Likin and Talmadge2015), where one might naively expect the fastest-growing instabilities are centred. This is readily apparent in low-global-magnetic-shear stellarators, such as HSX, where fluctuations have been observed to extend to $|\unicode[STIX]{x1D703}|\geqslant 30\unicode[STIX]{x03C0}$ or peak well away from $\unicode[STIX]{x1D703}=0$ , depending on wavenumber $k_{y}$ (Faber et al. Reference Faber, Pueschel, Proll, Xanthopoulos, Terry, Hegna, Weir, Likin and Talmadge2015). Increasing global magnetic shear has a localizing and stabilizing effect on drift-wave fluctuations due to the energy cost associated with field-line bending (Pearlstein & Berk Reference Pearlstein and Berk1969; Nadeem, Rafiq & Persson Reference Nadeem, Rafiq and Persson2001; Plunk et al. Reference Plunk, Helander, Xanthopoulos and Connor2014). When the global magnetic shear is small, such as $|{\hat{s}}|\lesssim 0.05$ for the HSX core region, fluctuations can extend long distances along the field line and consequently have long parallel correlation lengths. While these effects are accentuated in low-shear stellarators, they are not unique to it. Recent studies of strongly driven regions of axisymmetric configurations, such as the important pedestal region in tokamaks, show slab-like ITGs and micro-tearing modes can exist with strong finite- $\unicode[STIX]{x1D703}_{k}$ dependence (Fulton et al. Reference Fulton, Lin, Holod and Xiao2014; Hatch et al. Reference Hatch, Kotschenreuther, Mahajan, Valanju, Jenko, Told, Gorler and Saarelma2016b ; Chen & Chen Reference Chen and Chen2018).

It is crucial then that the geometry spanned by a fluctuation in the flux-tube domain have the proper instability drive and damping physics. Due to the lack of periodicity in a stellarator after one poloidal turn, the procedure to extend the parallel flux-tube domain using geometry elements for only one poloidal turn will result in physically incorrect geometric elements after one poloidal turn. This can lead to unfaithful drive and damping physics, which can in turn induce artificial self-correlations, leading to a breakdown of the statistical invariance posited by the parallel boundary condition. A superior way to deal with this issue is to construct the flux-tube domain using a geometry from multiple poloidal turns $n_{\text{pol}}>1$ . In principle, one must choose $n_{\text{pol}}$ to be large enough that the magnetic field geometry over a parallel correlation length is accurate. In practice, for modes with very long correlation lengths, one cannot in certain cases accurately capture all of the geometry due to computational cost constraints. However, using multiple poloidal turns decreases the deviation between the flux-tube geometry and the physical geometry at higher radial wavenumber, leading to a more accurate simulation domain. One can measure self-correlation by comparing parallel correlation lengths as function of poloidal turns used for the flux-tube geometry. As such, the $n_{\text{pol}}$ parameter is part of the convergence checking procedure.

An example of this is readily seen in figure 1, where the finite Larmor radius (FLR) term and the normal curvature term are shown for flux tubes constructed from different numbers of poloidal turns. The FLR term is expressed as

(2.10) $$\begin{eqnarray}k_{\bot }^{2}(z)=k_{x}^{2}\unicode[STIX]{x1D628}^{xx}(z)+2k_{x}k_{y}\unicode[STIX]{x1D628}^{xy}(z)+k_{y}^{2}\unicode[STIX]{x1D628}^{yy}(z),\end{eqnarray}$$

with elements of the metric tensor $\unicode[STIX]{x1D628}^{xx}(z)=\unicode[STIX]{x1D735}x(z)\boldsymbol{\cdot }\unicode[STIX]{x1D735}x(z)$ , $\unicode[STIX]{x1D628}^{xy}(z)=\unicode[STIX]{x1D735}x(z)\boldsymbol{\cdot }\unicode[STIX]{x1D735}y(z)$ , $\unicode[STIX]{x1D628}^{yy}(z)=\unicode[STIX]{x1D735}y(z)\boldsymbol{\cdot }\unicode[STIX]{x1D735}y(z)$ and the curvature drive at $\unicode[STIX]{x1D6FD}=0$ is

(2.11) $$\begin{eqnarray}{\mathcal{K}}=\frac{1}{B(z)}\hat{\boldsymbol{b}}(z)\times \unicode[STIX]{x1D735}B(z)\boldsymbol{\cdot }\unicode[STIX]{x1D735}y(z).\end{eqnarray}$$

It is clear that the geometry terms diverge after $z=\unicode[STIX]{x03C0}$ , which is a manifestation of the lack of periodicity in the poloidal angle for stellarators. Comparing the curvature drive and the FLR terms at $z=2\unicode[STIX]{x03C0}$ , one can see that for one poloidal turn, there is bad curvature and a small value of the FLR term, favourable conditions for driving a mode unstable. If a mode has parallel correlation lengths longer than one poloidal turn, the drive and damping terms at $z=2\unicode[STIX]{x03C0}$ will be unphysical with artificially high drive. This can lead to self-correlation and artificial pumping of the instability. Increasing the number of poloidal turns more realistically captures the geometry and enables the mode to sample geometry corresponding to more physical parallel correlation lengths. It should be stressed however, that this phenomenon is geometry dependent. Low-global-magnetic-shear geometries can exist where the local magnetic shear has large enough variation to provide mode localization within one poloidal turn or where the geometry is essentially periodic after one poloidal turn. In such circumstances, accurate flux-tube simulations can be performed using the geometry for only one poloidal turn, as was done in Xanthopoulos & Jenko (Reference Xanthopoulos and Jenko2007), Xanthopoulos et al. (Reference Xanthopoulos, Merz, Görler and Jenko2007).

Recently, a more accurate parallel boundary condition for stellarator symmetric flux tubes has been derived in Martin et al. (Reference Martin, Landremann, Xanthopoulos, Mandell and Dorland2018). This approach removes the discontinuities present at the ends of the simulation domain in the current conventional treatment, which can be seen at $\unicode[STIX]{x1D703}=\pm 4\unicode[STIX]{x03C0}$ of the $k_{\bot }^{2}$ curve in figure 1. Without a doubt, this is a step forward towards more accurate simulation of turbulence in stellarators. However, if the parallel correlation lengths of modes remain large, the above considerations will still apply and multiple poloidal turns will be needed.

2.3 Impact of small global magnetic shear

In order to minimize deleterious effects from large near-resonant Pfirsch–Schlüter currents and low-order rational surfaces, stellarators are often designed to operate in regions where the rotational transform is considerably constrained (Boozer Reference Boozer1998; Helander Reference Helander2014), such that the global magnetic shear $|{\hat{s}}|\lesssim 0.1$ . Magnetic shear tends to have a localizing effect on drift-wave instabilities, primarily by making the stabilizing FLR effects large. In figure 1, it is observed that modes localize along the field line to regions between local shear maxima. Curvature can also provide mode localization by modulating the drift frequency along a field line (Plunk et al. Reference Plunk, Helander, Xanthopoulos and Connor2014).

The structure of local magnetic shear along a field line is known to lead to different eigenmode types in the limit of small global magnetic shear. The work of Waltz & Boozer (Reference Waltz and Boozer1993) indicated that local magnetic shear can provide mode confinement along field lines and Plunk et al. (Reference Plunk, Helander, Xanthopoulos and Connor2014) demonstrated increased mode localization by artificially enhancing local shear spikes, essentially by ‘boxing’ the mode in, similar to a quantum state in a potential well. In the limit of zero magnetic shear, as has been studied for reverse-shear regions of internal transport barriers in tokamaks (Candy, Waltz & Rosenbluth Reference Candy, Waltz and Rosenbluth2004; Connor & Hastie Reference Connor and Hastie2004) and recently for the low-shear stellarator Wendelstein 7-X (Zocco et al. Reference Zocco, Xanthopoulos, Doerk, Connor and Helander2018), the gyrokinetic equation is known to take the form of a Mathieu equation. In the analysis of ITG modes in the Wendelstein 7-AS configuration, Bhattacharjee et al. (Reference Bhattacharjee, Sedlak, Similon, Rosenbluth and Ross1983) demonstrated the existence of weakly localized eigenmodes that extend far along field lines. These modes are bound due to resonances between the wave period and effective bounding potential and occur in regions where the Mathieu function is decaying. These regions are determined by the particular form and magnitude of the local shear along the field line and will vary between stellarator configurations. This subject was explored in depth for magnetohydrodynamic ballooning modes in 3-D geometry in Cuthbert & Dewar (Reference Cuthbert and Dewar2000) and Hegna & Hudson (Reference Hegna and Hudson2001) and numerically for drift waves by Nadeem et al. (Reference Nadeem, Rafiq and Persson2001) for the H-1NF configuration, which demonstrated a transition between localized and extended modes dependent on decreasing global magnetic shear.

A significant consequence of both 3-D geometry and small global shear is the existence of a plethora of unstable, but subdominant eigenmodes at every $k_{y}$ wavenumber (Merz Reference Merz2008; Faber et al. Reference Faber, Pueschel, Proll, Xanthopoulos, Terry, Hegna, Weir, Likin and Talmadge2015). Modes localized to the outboard midplane, modes with strong finite- $k_{x}$ dependence (amplitudes peaking at $\unicode[STIX]{x1D703}\neq 0$ ), and modes that extend along field lines can all be unstable concurrently at the same $k_{y}$ wavenumber. These subdominant modes have been shown to play an important role in the saturation dynamics of ITG turbulence in the HSX configuration (Pueschel et al. Reference Pueschel, Faber, Citrin, Hegna, Terry and Hatch2016; Hegna, Terry & Faber Reference Hegna, Terry and Faber2018). In particular, Hegna et al. (Reference Hegna, Terry and Faber2018) present a theory describing ITG turbulence saturation in stellarator geometry, where saturation occurs through nonlinear energy transfer to stable modes (Terry et al. Reference Terry, Baver and Gupta2006). It is found that for the HSX configuration, energy is primarily transferred to stable modes through subdominant, marginally stable modes. Furthermore, this observation was shown to be a function of geometry and was not observed in a quasi-axisymmetric configuration, which showed behaviour similar to conventional stellarators and high- ${\hat{s}}$ tokamaks. There, typically only at most a few modes tend to be destabilized at each wavenumber, and the turbulence is strongly influenced by unstable mode/zonal flow interactions (Hatch et al. Reference Hatch, Terry, Jenko, Merz, Pueschel, Nevins and Wang2011b ; Hegna et al. Reference Hegna, Terry and Faber2018; Terry et al. Reference Terry, Faber, Hegna, Mirnov, Pueschel and Whelan2018; Whelan, Pueschel & Terry Reference Whelan, Pueschel and Terry2018). Thus care must be taken to properly simulate not only the most unstable eigenmode, but also a significant part of the subdominant spectrum.

3.1 Trapped electron modes

Calculations of TEM, relevant for HSX operation (Faber et al. Reference Faber, Pueschel, Proll, Xanthopoulos, Terry, Hegna, Weir, Likin and Talmadge2015) focuses on high-density-gradient drive, where the impact of extended modes is the most prominent. The parameters used are $a/L_{n}=4$ , $a/L_{T\text{i}}=0$ , $a/L_{T\text{e}}=5\times 10^{-4}$ , $T_{\text{i}}/T_{\text{e}}=1$ , $\unicode[STIX]{x1D6FD}=0$ . The gradient scales for the density and temperature normalized to the average plasma minor radius are $a/L_{n,T}$ , where $L_{\unicode[STIX]{x1D709}}=|\unicode[STIX]{x1D735}\unicode[STIX]{x1D709}|/\unicode[STIX]{x1D709}$ . The ion and electron temperatures are $T_{\text{i},\text{e}}$ and $\unicode[STIX]{x1D6FD}=8\unicode[STIX]{x03C0}n_{\text{e}}T_{\text{e}}/B_{0}^{2}$ is the normalized ratio of electron pressure to magnetic pressure. The eigenspectrum is examined for $k_{y}=(0.2,0.4,0.6,0.8)$ and $k_{x}=j2\unicode[STIX]{x03C0}n_{\text{pol}}|{\hat{s}}|k_{y}$ for integer $j=(0,\pm 1,\ldots ,\pm 8)$ , and is shown in figure 2. As expected for pure $\unicode[STIX]{x1D735}n$ -driven TEM, almost all of the unstable modes are propagating in the electron-diamagnetic direction, corresponding to negative value for real frequency using Gene conventions. Examples of different mode types are given by the letter labels and the corresponding mode structures are shown in figure 3. There is a clustering of modes around zero real frequency, labelled by ‘C’, which have extended structure along field lines, consistent with the previous findings of Nadeem et al. (Reference Nadeem, Rafiq and Persson2001). The modes with definite non-zero frequencies, ‘A’ and ‘B’ tend to be much more strongly localized in helical drift wells. These modes do not need to balloon around $\unicode[STIX]{x1D703}=0$ , and they display significant finite- $k_{x}$ amplitudes as well as tearing parity in electrostatic potential  $\unicode[STIX]{x1D6F7}$ , consistent with previous observations of TEMs in HSX (Faber et al. Reference Faber, Pueschel, Proll, Xanthopoulos, Terry, Hegna, Weir, Likin and Talmadge2015). For most modes in the spectrum, increasing the number of poloidal turns leads to a decrease in the growth rates, and is shown in figure 4 for $k_{y}=0.2$ for one and four poloidal turns. This behaviour is connected to the fact that by using more poloidal turns to construct the flux tube, more correct drive and damping physics is included, leading to a more accurate growth rate and frequency calculation for extended modes.

Figure 2. Eigenspectrum for the strongly driven $\unicode[STIX]{x1D735}n$ TEM in HSX. The horizontal axis is the growth rate $\unicode[STIX]{x1D6FE}$ and the vertical axis is the real frequency  $\unicode[STIX]{x1D714}$ . Different $k_{y}$ are indicated different symbols and colours. Examples of different types of modes are given by the labels ‘A’ ( $k_{y}=0.8$ ), ‘B’ ( $k_{y}=0.4$ ), ‘C’ ( $k_{y}=0.8$ ) and ‘D’ ( $k_{y}=0.4$ ), and the mode structure is shown in figure 3.

Figure 3. Electrostatic potential eigenmode structures for the TEM branches labelled ‘A’ (upper left), ‘B’ (upper right), ‘C’ (lower left) and ‘D’ (lower right) from figure 2. The modes show conventional ballooning behaviour (A), finite $k_{x}$ dependence (B), extended structure along the field line (C) and two-scale structure with an extended envelope along the field line (D).

A peculiar branch of eigenmodes, labelled ‘D’ in figure 2 is characterized by marginal stability modes, with a real frequency proportional to $k_{y}$ in the ion diamagnetic direction. The structure of the electrostatic potential along the field line for the eigenmodes in HSX, shown in figure 3 D, displays an extended, two-scale structure, with the outer scale varying on scales much longer than local helical drift wells. While these modes appear connected to small-global-magnetic-shear values, the exact mechanism responsible for setting the outer-scale envelope has not yet been determined. In the reverse-shear tokamak scenario studied in Candy et al. (Reference Candy, Waltz and Rosenbluth2004), the small, but non-zero global magnetic shear set the envelope scale. In 3-D configurations, variations in the local shear can determine mode localization through Mathieu resonances (Bhattacharjee et al. Reference Bhattacharjee, Sedlak, Similon, Rosenbluth and Ross1983) or through a method similar to Anderson localization, where incommensurate helical periods in the magnetic equilibrium cause localization (Cuthbert & Dewar Reference Cuthbert and Dewar2000; Hegna & Hudson Reference Hegna and Hudson2001). The mode branch displays eigenmode structures with higher-harmonic envelopes that are increasingly damped.

Figure 4. TEM eigenspectrum at $k_{y}=0.2$ for different poloidal turn values. Modes from one poloidal turn are the solid red diamonds and from four poloidal turns are the hollow blue diamonds. Generally, the modes at four poloidal turns are more stable, and there are fewer unstable modes than for one poloidal turn. Importantly, the extended ion mode branch (mode ‘D’ in figure 3), transitions from unstable to stable at four poloidal turns.

Like the other eigenmodes, the growth rate of the extended ion mode branch is sensitive to the number of poloidal turns used to resolve the geometry. Figure 4 shows the extended ion mode branch at $\unicode[STIX]{x1D714}\approx 0.8$ is sensitive to the number of poloidal turns and that the linear damping increases with number of poloidal turns. At one poloidal turn, this mode branch even shows slight instability at $k_{y}=0.2$ . This sensitivity is observed despite the parallel correlation lengths of these modes being significantly larger than a few poloidal turns. In the analysis presented in Cuthbert & Dewar (Reference Cuthbert and Dewar2000) and Candy et al. (Reference Candy, Waltz and Rosenbluth2004), the eigenvalue is dependent on the solution to the inner-scale equation. More accurately resolving the inner-scale solution associated with the helical drift wells serves to stabilize the overall growth rate. Improperly resolving this extended scale mode has consequences for nonlinear simulation, as will be detailed in § 4.

3.2 Ion temperature gradient modes

Calculations of ITG eigenmodes in the HSX configuration have been performed, focusing on a pure collisionless ITG drive with kinetic electrons and the following parameters: $a/L_{T\text{i}}=3$ , $a/L_{T\text{e}}=0$ , $a/L_{n}=0$ , $T_{\text{i}}/T_{\text{e}}=1$ , $\unicode[STIX]{x1D6FD}=5\times 10^{-4}$ . Because the impact of eigenmodes with large parallel correlation lengths is more prevalent at low $k_{y}$ , the eigenspectrum is shown for $k_{y}<1$ in figure 5. As was observed with the TEM calculations, there are many unstable eigenmodes for every  $k_{y}$ . A similar mix of mode structures is observed, including strongly ballooning modes, $k_{x}\neq 0$ modes, and extended modes along the field line. Also observed are two-scale, ion-direction-propagating modes with extended-envelope structure, labelled by ‘A’, ‘B’ and ‘C’ in figure 5. The corresponding mode structures are shown in figure 6. While these modes have similar envelope behaviour as in the TEM case, a distinguishing characteristic of the ITG case is that these modes are much more unstable and do not become stable with increasing poloidal turns. Additionally, instability of the extended ion modes is observed at high $k_{y}$ for the ITG case, as shown for $k_{y}=0.9$ in figure 5, while only marginal instability is observed at low $k_{y}$ in figure 2 for the TEM case.

Figure 5. ITG eigenmode spectrum for different $k_{y}$ values, denoted by different symbols and colours. The modes at $k_{y}=0.9$ labelled ‘A’, ‘B’ and ‘C’ are modes with two-scale behaviour, and the mode structures are shown in figure 6.

Figure 6. Electrostatic potential for eigenmodes ‘A’ (top), ‘B’ (middle) and ‘C’ (bottom) from figure 5. The real part is the solid line and the imaginary part is the dotted black line. Similar to the mode ‘D’ in figure 3, the ITG eigenmodes display two-scale structure, with an outer-scale envelope and inner-scale structure set by the helical magnetic structure.

Despite these differences, given that these two-scale, extended-envelope ion modes are observed in both TEM and ITG calculations, the existence of these modes is connected to the low-magnetic-shear nature of devices such as HSX rather than a particular instability drive. An important observation is that these modes are not observed in the unstable eigenspectrum when an adiabatic electron approximation is used. This is shown in figure 7, where the eigenvalue spectrum is computed for both the adiabatic and kinetic electron treatment for a particular  $k_{y}$ . Thus to fully resolve the drift-wave dynamics in a low-shear stellarator like HSX, it is important that kinetic electrons are used, so that crucial physics is not overlooked.

3.3 Zero-magnetic-shear approximation

The HSX shear value, ${\hat{s}}=-0.046$ , represents only a small deviation from ${\hat{s}}=0$ . Equation (2.2) shows that order-unity variations in $\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}$ and $\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}$ terms will dominate over the secular $\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ term, provided $\unicode[STIX]{x1D703}-\unicode[STIX]{x1D703}_{k}$ is small. A reasonable approximation for ballooning modes is then to neglect the global magnetic shear, which will be referred to as the zero-shear technique. In the zero-shear approach, linear modes no longer couple through magnetic shear, as the boundary condition is strictly periodic. A substantial benefit of this approach to low-magnetic-shear configurations is that the radial box size of the simulation is no longer determined by (A 7), which states $L_{x}\propto ({\hat{s}}k_{y}^{\text{min}})^{-1}$ . For HSX, using $k_{y}^{\text{min}}=0.1$ yields $L_{x}\approx 220\unicode[STIX]{x1D70C}_{\text{s}}$ , a value comparable with macroscopic equilibrium scale lengths (Faber et al. Reference Faber, Pueschel, Proll, Xanthopoulos, Terry, Hegna, Weir, Likin and Talmadge2015). Despite such large simulation domains, this is not physically problematic provided the turbulent correlation lengths are both smaller than experimental correlation lengths and the background length scales, such as  $L_{n,T}$ . A large radial box can be thought of as equivalent to multiple ‘independent’ simulation domains and turbulent statistics can be accumulated quickly. The numerical resolution, however, must also be scaled appropriately, making such simulations expensive to carry out. With zero magnetic shear, one may choose a more reasonable radial domain size and resolution, potentially reducing computational costs.

Figure 7. Comparison of the ITG eigenmode spectrum at $k_{y}=0.9$ for kinetic electrons (red) and adiabatic electrons (blue). Adding kinetic electron effects primarily the real frequency and reduces the number of unstable modes. However the extended ion modes that are unstable for kinetic electrons are absent from adiabatic electron calculations.

With ${\hat{s}}=0$ , an extended parallel computational domain cannot be constructed by adding radial modes, as demonstrated in § 2. Zero-shear calculations then require that the parallel computation domain extend sufficiently far along the field line to better accommodate physical parallel correlation lengths. The strictly periodic boundary condition makes it easier to enhance growth rates through self-correlation if an eigenmode is not sufficiently decayed at the ends of the simulations domain. This is not only an issue for linear eigenmode calculations. A nonlinear simulation can be dominated by artificial linear self-correlation if growth rates are large enough compared to nonlinear decorrelation rates, which is demonstrated in § 4.

Figure 8 shows the application of the zero-shear computational technique to HSX for dominant TEM eigenmode calculations where only one poloidal turn was used for the geometry. In a finite-shear simulation, due to linear $k_{x}$ coupling, the fastest growing connected $k_{x}$ mode determines the overall mode growth rate. As linear modes are not coupled in zero-shear simulations, one must independently scan over $k_{x}$ to capture finite $k_{x}$ effects. At higher  $k_{y}$ , the finite-shear and zero-shear calculations agree quite well for $k_{x}=0$ , where the most unstable modes strongly balloon at the outboard midplane and are sufficiently localized within one poloidal turn. For intermediate values $0.3\leqslant k_{y}\leqslant 0.7$ , strong growth at finite $k_{x}$ is observed and scanning $k_{x}$ for the zero-shear modes reproduces the finite-shear values. However, for $k_{y}<0.2$ , clear divergences are observed between the zero-shear and finite-shear calculations. The zero-shear growth rates exhibit significantly larger growth rates than the finite-shear counterparts and can be identified as the zero-shear counterpart to the extended envelope ion modes in the finite-shear case. As is shown in figure 9, for geometry using one poloidal turn, this zero-shear mode branch (red diamonds) tracks the finite-shear (black diamonds) growth rates and frequencies relatively well for higher $k_{y}$ values. Around $k_{y}=0.3$ , the zero-shear growth rates diverge and become significantly more unstable. This is problematic for nonlinear simulations that require $k_{y}^{\text{min}}\leqslant 0.1$ , and will be discussed in § 4.

Figure 8. Linear growth rate spectrum of TEMs for one poloidal turn for the finite-shear approach (red) and the zero-shear approach (black and blue). The blue points are from calculations where the $k_{x}$ value was varied to find the maximum growth rate, while the black points are the $k_{x}=0$ streamer instability. The non-zero $k_{x}$ effects are required to reproduce the finite-shear spectrum for $0.3\leqslant k_{y}\leqslant 0.7$ .

Figure 9. Eigenspectrum of the artificially enhanced mode (red diamonds) of figure 8 for a range of wavenumbers  $k_{y}$ . Select $k_{y}$ modes are identified by the outlined symbols and compared with the finite-shear counterpart for one poloidal turn (black diamonds). All finite-shear eigenvalues are clustered near marginal stability, whereas zero-shear modes are artificially enhanced at one poloidal turn.

Linear mode differences between the finite and zero-shear techniques become small when multiple poloidal turns are utilized, and in particular, growth rate differences at wavenumbers $k_{y}\leqslant 0.1$ are significantly reduced. The periodic parallel boundary conditions coupled with geometry from one poloidal turn does not allow modes to sufficiently decay and the growth rate is artificially enhanced. This indicates that for the low- $k_{y}$ modes, the parallel correlation lengths for HSX are longer than one poloidal turn and that the enhanced growth rates for the zero-shear calculations were due to self-correlating fluctuations. For the half-toroidal flux surface in HSX, linear TEMs show convergence at 4 poloidal turns.

As noted previously, for stellarator purposes, it is not sufficient to just examine convergence of the most unstable mode, as subdominant modes can play an important role in the turbulent dynamics. Figure 10 shows a comparison of the eigenmode spectrum between the finite and zero-shear techniques for a TEM case in HSX, using four poloidal turns at $k_{y}=0.7$ . The zero-shear eigenmodes are computed by scanning $k_{x}$ values that would be coupled via the parallel boundary condition in the finite-shear calculation. Figure 10 shows that both the zero-shear and finite-shear growth rates across the spectrum are replicated. It is not necessary to match exactly every single subdominant mode between the two approaches. Pueschel et al. (Reference Pueschel, Faber, Citrin, Hegna, Terry and Hatch2016) demonstrated the impact of the subdominant modes in HSX is to generate broadband turbulence, lessening the impact of any particular eigenmode. More importantly, the clustering of eigenmodes in the spectrum is very well reproduced between the calculation techniques in figure 10, which leads to consistent turbulence energy drive and dissipation. This gives confidence that the zero-shear technique is accurately calculating the eigenmode spectrum when the parallel correlation length is sufficiently resolved.

Figure 10. Comparison of the finite-shear (hollow red diamonds) and zero-shear (solid black diamonds) eigenspectrum at $k_{y}=0.7$ for the $\unicode[STIX]{x1D735}n$ -driven TEM. Sufficient agreement is observed between the two computational approaches when multiple poloidal turns are used. In particular, the zero-shear technique recovers the appropriate clustering of eigenmodes, including the marginally stable ion mode branch (branch ‘D’ of figure 2).