# Ion-temperature-gradient stability near the magnetic axis of quasisymmetric
stellarators

R. Jorge, M. Landreman  Institute for Research in Electronics and Applied
Physics, University of Maryland, College Park, MD 20742, USA<EMAIL_ADDRESS>
###### Abstract

The stability of the ion-temperature gradient mode in quasisymmetric
stellarators is assessed. This is performed using a set of analytical
estimates together with linear gyrokinetic simulations. The peak growth rates,
their corresponding real frequencies and wave-vectors are identified. A
comparison is made between a first-order near-axis expansion model and eleven
realistic designs obtained using numerical optimization methods. It is found
that while the near-axis expansion is able to replicate the growth rates, real
frequencies and perpendicular wave-vector at the inner core (both using
simplified dispersion relations and first-principle gyrokinetic simulations),
it leads to an overestimation of the growth rate at larger radii. An
approximate analytic solution of the ITG dispersion relation for the non-
resonant limit suggests growth rates could be systematically higher in quasi-
axisymmetric (QA) configurations compared to quasi-helically (QH) symmetric
ones. However except for very close to the axis, linear gyrokinetic
simulations do not show systematic differences between QA and QH
configurations.

##  1 Introduction

The feasibility of magnetic confinement fusion relies on a delicate balance
between upholding a hot enough plasma in the core of the reactor and the
inevitable transport of particles and energy to the wall. Such transport is
termed neoclassical when associated with geometrical effects in standard
collisional transport theories, and termed turbulent when associated with
perturbations in the plasma driven by underlying instabilities. The reduction
of neoclassical transport to acceptable levels has been, to an extent,
attained with the tokamak and omnigeneous (such as quasisymmetric and quasi-
isodynamic) stellarator designs. The reduction of turbulent transport,
however, has been a focus of theory and modelling of magnetic confinement
plasma physics in the past few decades. In the core, two microinstabilities
are thought to be the main drive of turbulent transport: the Ion Temperature
Gradient (ITG) mode and the Trapped Electron Mode (TEM) [1, 2]. In this work,
we focus on the question of linear ITG stability and its parameter dependence
in the core of a particular class of omnigeneous stellarator designs, namely
quasisymmetric stellarators.

Quasisymmetry is a property of toroidal magnetic fields that ensures
omnigenity, i.e., that collisionless particle trajectories remain confined
inside the device over long periods of time. Although there are numerous
equivalent definitions of quasisymmetry, here we define it as the condition
that the magnitude of the magnetic field $B=|\mathbf{B}|$ varies on a surface
of constant toroidal flux only through a fixed combination of the poloidal
$\theta$ and toroidal $\varphi$ Boozer angles $B=B(\psi,M\theta-N\varphi)$ for
some integers $M$ and $N$ [3, 4, 5]. There are currently two ways of
developing quasisymmetric stellarator designs: analytically, using the near-
axis expansion [6, 7, 8] or a near-axisymmetry formulation [9]; or
numerically, using one of the stellarator optimization methods available [10,
11, 12]. In this work, the microstability of quasisymmetric designs near the
magnetic axis will be assessed using both methods, and the advantages and
shortcomings of using either method are highlighted.

Microinstabilities are driven by the presence of free energy associated with
temperature and density gradients and are, therefore, ubiquitous in tokamaks
and stellarators. These develop locally in a plasma at the gyroradius scale,
as opposed to macroinstabilities that affect the plasma globally.
Microinstabilities have regions of favorable development in the device,
usually assessed by their relative "good" or "bad" curvature properties,
evaluated using the sign and magnitude of the magnetic drift frequency
$\omega_{d}\simeq T\mathbf{B}\times\nabla B\cdot\mathbf{k}_{\perp}/(eB^{3})$
where $T$ is the background temperature, $e$ the elementary charge and
$\mathbf{k}_{\perp}$ the wave-vector of the mode perpendicular to the magnetic
field $\mathbf{B}$. If this frequency is positive (or, more accurately, if it
has the same sign as the diamagnetic drift-frequency), trapped particle modes
are likely to be unstable [13]. While $\omega_{d}$ is positive in the outboard
side of tokamaks, rendering unstable the region where trapped particles spend
most of their time, in stellarators, the locations of bad curvature need to be
evaluated numerically and are not necessarily in the outboard region. As a
result, many studies aimed at assessing the influence of magnetic field
geometry on plasma stability and turbulence have been performed, particularly
using gyrokinetic theory [2]. The freedom associated with the length and
location of the bad curvature region has also allowed for stellarator
optimization criteria based on the sign of $\omega_{d}$ and the non-overlap of
trapped and "bad" curvature regions [14]. While the first non-axisymmetric
gyrokinetic studies were performed using the linear eigenvalue FULL code [15,
16], only recently have nonlinear simulations been performed [17] together
with stellarator optimization efforts aimed at reducing turbulent transport
[18].

In this work, the ITG instability is analyzed first by incorporating the near-
axis expansion into several known analytical estimates and then by performing
linear gyrokinetic simulations for a set of eleven quasisymmetric designs
using the GS2 code [19]. This instability can either be destabilized by the
presence of an equilibrium temperature gradient which alters the temperature
of a perturbation travelling at the sound speed that is in phase with a
perturbed $\mathbf{E}\times\mathbf{B}$ drift, usually called the slab branch
of the mode, or by the coupling of the grad-$B$ and curvature drifts with the
equilibrium temperature gradient, usually called the toroidal branch of the
mode [20]. We mention that, although this work focuses on the properties of
the linearized gyrokinetic equation, many non-linear studies of gyrokinetic
ITG turbulence and saturation mechanisms have been carried out, showing the
interplay between plasma geometry and stability [21, 22, 23, 24].

The purpose of this paper is twofold. First, we show how quasisymmetric near-
axis designs could lead to simplified yet accurate models of ITG stability.
Such models can then be used as optimization metrics when seeking for new
stellarator designs with minimal levels of microinstability. Second, we extend
previous studies comparing first order near-axis expansions with the designs
based on numerical optimization methods [25, 26]. While in Ref. [26] a
comparison was made between the geometric quantities used in linear and non-
linear models of plasma dynamics, here we aim at using metrics that allow for
a more quantitative comparison, such as the growth rate of underlying unstable
modes. Owing to the difficulty of simulating complex geometries in the flux-
tube approximation [27, 28], convergence tests are shown for the selected base
case simulation scenario for every stellarator design. The details of the
implementation of non-axisymmetric geometries in the GS2 code can be found in
Ref. [29].

This paper is structured as follows. In Section 2, the gyrokinetic equation
and its underlying approximations used to assess microstability are presented.
Also, the near-axis model is introduced together with the eleven
quasisymmetric geometries being studied. In Section 3 various estimates aimed
at assessing linear and non-linear stellarator microstability are rewritten
leveraging the simplifying assumptions used in the near-axis expansion
framework. A benchmark configuration from a near-axis solution is used in
Section 4 to detail the analysis employed here when dealing with linear
gyrokinetic simulations, including scans over the physical parameters at play,
the resulting eigenmodes and eigenfrequencies. The analysis performed for the
benchmark case of Section 4 is carried out in Section 5 for the selected
eleven quasisymmetric designs. The conclusions follow.

##  2 Physical Model and Geometry

###  2.1 Gyrokinetic Equation

In this work, we employ the gyrokinetic approximation, i.e., we assume that
the plasma is strongly magnetized, and that fluctuations have small amplitude
and low frequency [30]. These assumptions are, in general, satisfied in
magnetic confinement nuclear fusion devices. In here, we find that the
particle (electron and ion) Larmor radius, $\rho=v_{th}/\Omega$ with
$v_{th}=\sqrt{2T/m}$ the thermal velocity, $\Omega=qB/m$ the cyclotron
frequency, $T$ the temperature, $m$ the mass, $q$ the charge and $B$ the
magnetic field strength, is much smaller than typical length scales of the
macroscopic equilibrium, $L$. Furthermore, the frequency of fluctuations,
$\omega$, is much smaller than the cyclotron frequency, $\Omega$. This allows
us to average the particle’s equations of motion, effectively turning their
Larmor orbit into a guiding-center one via the gyroaverage operator
$\left<...\right>$. Using the small fluctuation assumption and focusing on an
electrostatic regime, we split the distribution function of each species,
$f=F_{0}+\delta f$, and the electrostatic potential $\phi=\delta\phi$, into an
equilibrium (denoted with a subscript $0$) and a fluctuating (denoted with
$\delta$) part, where $\delta f\ll F_{0}$ and $q\phi/T\ll 1$. Furthermore, the
fluctuating distribution function is split into an adiabatic and non-adiabatic
part as $\delta f=h-(q\phi/T)F_{0}$. For the equilibrium, we take $F_{0}$ to
be a Maxwell-Boltzmann distribution with thermal velocity $v_{th}$ and density
$n$, and both the background density $n=n(\psi)$ and temperature $T=T(\psi)$
are taken to be functions of the magnetic toroidal flux $\psi$ only, with
characteristic lengths $L_{n}$ and $L_{T}$, respectively. Due to its
importance in assessing regions of stability in the gyrokinetic space of
parameters, we define $\eta=L_{n}/L_{T}$ as the ratio between the two
characteristic lengths. In the following, for simplicity, we denote
$\left<\phi\right>$ and $\left<h\right>$ as $\phi$ and $h$, respectively, we
assume that electrons are adiabatic with density $n$ and temperature $T_{e}$
and that a single ion species is present, with its gyroaveraged distribution
function given by $h$, density $n$ and temperature $T_{i}=\tau_{T}T_{e}$.

The nonlinear gyrokinetic equation solved by GS2, in the linear electrostatic
limit, can be written as [31]

$\frac{\partial h}{\partial
t}+(v_{\parallel}\mathbf{b}+\mathbf{v}_{d})\cdot\nabla
h+\frac{\mathbf{b}\times\nabla\phi}{B}\cdot\nabla
F_{0}=\left<C\right>+\frac{qF_{0}}{T}\frac{\partial\phi}{\partial t}.$ (1)

The particle’s grad-B and curvature drifts are contained in

$\mathbf{v}_{d}=\frac{\mathbf{b}}{\Omega}\times\left(v_{\parallel}^{2}\mathbf{b}\cdot\nabla\mathbf{b}+\frac{v_{\perp}^{2}}{2}\frac{\nabla
B}{B}\right).$ (2)

We note that, in Eq. 1, the spatial derivatives are taken at constant energy
$\mathcal{E}=mv^{2}/2$ and magnetic moment $\mu=mv_{\perp}^{2}/2B$ and
collisions are modeled via the collision operator $C$. The system of equations
is closed via the quasineutrality equation which, in the limit of adiabatic
electrons, reads

$\int d\mathbf{v}\left<h\right>_{\mathbf{r}}=n(1+\tau_{T})\frac{q\phi}{T},$
(3)

where $\left<h\right>_{\mathbf{r}}$ denotes a gyroaverage performed at
constant position vector $\mathbf{r}$ as opposed to $\left<h\right>$ which is
performed at constant guiding-center locaion $\mathbf{R}$. A relation between
the two can be found by Fourier decomposing $h$ into a parallel and
perpendicular wave-vectors
$\mathbf{k}=k_{\parallel}\mathbf{b}+\mathbf{k}_{\perp}$ with
$\mathbf{b}=\mathbf{B}/B$, allowing us to write, in Fourier space,
$\left<h\right>_{\mathbf{r}}=J_{0}(k_{\perp}v_{\perp}/\Omega)\left<h\right>$.

In order to perform a stability analysis, we write Eq. 1 in the field-line
following limit. Additionally, due to the high temperature of the plasma in
the core, we look at the collisionless limit $\nu_{e}\ll\omega$ where
$\nu_{e}$ is the electron-ion collision frequency, and set $C=0$. In the
field-line following limit, we write the perturbed quantities $h$ and $\phi$
as

$h=\hat{h}(l)e^{i(S-\omega t)},$ (4)

where $l$ is the distance along a magnetic field line and
$\mathbf{b}\cdot\nabla S=0$, leading us to define
$\mathbf{k}_{\perp}\equiv\nabla S$. In this limit, and normalizing $\phi$ by
$q/T$, Eq. 1 can be written as

$iv_{\parallel}\mathbf{b}\cdot\nabla\hat{g}+(\omega-\tilde{\omega}_{d})\hat{g}=\hat{\phi}F_{0}(\omega-\tilde{\omega}_{*})J_{0}\left(\frac{v_{\perp}k_{\perp}}{\Omega}\right),$
(5)

with $\tilde{\omega}_{d}=\mathbf{k}_{\perp}\cdot\mathbf{v}_{d}$ the magnetic
drift-frequency. When a low $\beta$ limit is taken, $\tilde{\omega}_{d}$ can
be written as
$\tilde{\omega}_{d}=\omega_{d}(v_{\parallel}^{2}+v_{\perp}^{2}/2)/v_{th}^{2}$
with $\omega_{d}$ its velocity-independent counterpart. The velocity-dependent
diamagnetic frequency $\tilde{\omega}_{*}$ is given by

$\tilde{\omega}_{*}=\omega_{*}\left[1+\eta\left(\frac{v^{2}}{v_{th}^{2}}-\frac{3}{2}\right)\right],$
(6)

where $\omega_{*}=\mathbf{b}\times\mathbf{k}_{\perp}\cdot\nabla F_{0}/(m\Omega
TF_{0})$ usually written as $\omega_{*}=(Tk_{\alpha}/q)d\ln n/d\psi$ when a
Clebsch representation for the magnetic field is used.

As a final note, we mention here the normalizations used to output GS2
quantities. Lengths are normalized to a fixed length $a$ (a measure of the
minor radius of the device), time is normalized to $a/v_{th}$ and perturbed
quantities are scaled up by $a/\rho_{th}$. In the axisymmetric case, $a$ is
half of the diameter of the last closed flux surface measured at the elevation
of the magnetic axis. Here, however, we identify $a$ with the effective minor
radius computed by VMEC [32], named Aminor_p in its wout*.nc file, defined as
$\pi$ (Aminor_p$)^{2}=({1}/{2\pi})\int_{0}^{2\pi}S(\phi)d\phi$ where $\phi$ is
the toroidal angle associated with the cylindrical coordinate system
$(R,Z,\phi)$ ($\phi$ is not to be confused with the electrostatic potential)
and $S(\phi)$ is the area of the outer surface’s cross-section in the $R-Z$
plane.

###  2.2 Coordinates and Magnetic Geometry

We now rewrite the differential operators in the linearized gyrokinetic
equation, Eq. 5, in terms of the straight-field-line magnetic coordinates
commonly referred to as Boozer coordinates [33] and employ the Garren-Boozer
near-axis expansion formalism [6, 34]. We note that such formalism uses an
inverse coordinate approach, writing $\mathbf{r}$ as a function of Boozer
coordinates, as opposed to formalisms using a direct coordinate approach,
which write the toroidal flux as a function of the spatial coordinates
$\mathbf{r}$ [35, 36, 37]. First, we express the magnetic field $\mathbf{B}$
using a Clebsch representation of the form

$\mathbf{B}=\nabla\psi\times\nabla\alpha,$ (7)

where $2\pi\psi$ is the toroidal flux and $\alpha$ is a field line label.
Taking as a spatial coordinate system the set $(\psi,\alpha,z)$ with $z$ any
quantity independent of $\psi$ and $\alpha$, the eight independent geometrical
quantities $\mathbf{Q}$ needed to solve Eq. 5 are [26]

$\displaystyle\mathbf{Q}=$ $\displaystyle\left\\{B,\mathbf{b}\cdot\nabla
z,|\nabla\psi|^{2},|\nabla\alpha|^{2},\nabla\psi\cdot\nabla\alpha,\right.$
$\displaystyle\left.(\mathbf{b}\times\nabla
B)\cdot\nabla\alpha,(\mathbf{b}\times\nabla
B)\cdot\nabla\psi,(\mathbf{b}\times\boldsymbol{\kappa})\cdot\nabla\alpha\right\\}.$
(8)

From here onward, we employ Boozer coordinates $(\psi,\theta,\varphi)$ with
$\theta$ and $\varphi$ the poloidal and toroidal angles, respectively, and
choose $z=\varphi$. For convenience, a helical angle
$\vartheta=\theta-N\varphi$ with $N$ an integer is introduced, such that the
magnetic field $\mathbf{B}$ can be written as

$\displaystyle\mathbf{B}=\nabla\psi\times\nabla\vartheta+\iota_{N}\nabla\varphi\times\nabla\psi,$
(9) $\displaystyle=\beta\nabla\psi+I\nabla\vartheta+(G+NI)\nabla\varphi,$ (10)

where $\iota_{N}=\iota$ with $\iota$ the rotational transform and $I,G$ and
$\iota$ are constants on $\psi$ surfaces, i.e., $I=I(\psi)$, $G=G(\psi)$ and
$\iota=\iota(\psi)$. A relation between $\alpha$ and the angles $\vartheta$
and $\varphi$ can be found by comparing Eq. 7 and Eq. 9, yielding
$\alpha=\vartheta-\iota_{N}\varphi=\theta-\iota\varphi$.

The geometry coefficients $\mathbf{Q}$ in Eq. 8 are evaluated at lowest order
in $\epsilon=r/\mathcal{R}$, where $r$ is defined as $2\pi|\psi|=\pi
r^{2}B_{0}$ and $\mathcal{R}$ a scale length representing the major radius of
the device (e.g., the inverse of the maximum axis curvature). Furthermore, we
focus on quasisymmetric magnetic fields, i.e., we use magnetic geometries
where $B=B(\psi,M\vartheta-N\varphi)$ with $M$ and $N$ integers. As shown in
Ref. [26], this yields the following set of geometric coefficients

$\displaystyle B$ $\displaystyle=B_{0}(1+r\overline{\eta}\cos\vartheta),$ (11)
$\displaystyle\mathbf{b}\cdot\nabla z$
$\displaystyle=s_{G}\frac{2\pi}{L}(1+r\overline{\eta}\cos\vartheta),$ (12)
$\displaystyle|\nabla\psi|^{2}$
$\displaystyle=r^{2}\frac{B_{0}^{2}}{\overline{\eta}^{2}\kappa^{2}}\left[\overline{\eta}^{4}\sin^{2}\vartheta+\kappa^{4}(\cos\vartheta-\sigma\sin\vartheta)^{2}\right],$
(13) $\displaystyle|\nabla\alpha|^{2}$
$\displaystyle=\frac{1}{r^{2}\overline{\eta}^{2}\kappa^{2}}\left[\overline{\eta}^{4}\cos^{2}\vartheta+\kappa^{4}(\sigma\cos\vartheta+\sin\vartheta)^{2}\right],$
(14) $\displaystyle\nabla\psi\cdot\nabla\alpha$
$\displaystyle=\frac{s_{\psi}B_{0}}{2\overline{\eta}^{2}\kappa^{2}}\left(\left[\overline{\eta}^{4}+\kappa^{4}\left(\sigma^{2}-1\right)\right]\sin
2\vartheta-2\kappa^{4}\sigma\cos 2\vartheta\right),$ (15)
$\displaystyle(\mathbf{b}\times\nabla B)\cdot\nabla\alpha$
$\displaystyle=\frac{s_{\psi}}{r}B_{0}\overline{\eta}\cos\vartheta,$ (16)
$\displaystyle(\mathbf{b}\times\nabla B)\cdot\nabla\psi$
$\displaystyle=rB_{0}^{2}\overline{\eta}\sin\vartheta,$ (17)
$\displaystyle(\mathbf{b}\times\boldsymbol{\kappa})\cdot\nabla\alpha$
$\displaystyle=\frac{s_{\psi}}{r}\overline{\eta}\cos\vartheta.$ (18)

In the set of Eqs. 11, 12, 13, 14, 15, 16, 17 and 18, $B_{0}$ is the constant
magnetic field on-axis, $r=\sqrt{2\psi/B_{0}}$ is an effective minor radius,
$\overline{\eta}$ is a constant that is related to the elongation of the
elliptical flux-surfaces, $\kappa=\kappa(\varphi)$ is the axis curvature,
$s_{\psi}=$sgn$(\psi)$, $s_{G}=$sgn$(G)$, and $\sigma=\sigma(\varphi)$ is the
solution of

$\frac{d\sigma}{d\varphi}+(\iota_{0}-N)\left(\frac{\overline{\eta}^{4}}{\kappa^{4}}+1+\sigma^{2}\right)-\frac{s_{G}L\overline{\eta}^{2}}{\pi\kappa^{2}}\left(\frac{I_{2}}{B_{0}}-s_{\psi}\tau\right)=0,$
(19)

where $\tau=\tau(\varphi)$ is the axis torsion, $\iota_{0}$ is the on-axis
rotational transform, $I_{2}$ is the lowest order toroidal current $I\simeq
I_{2}r^{2}$ and $L$ is the total length of the magnetic axis curve
$\mathbf{r}_{0}(\phi)$

$L=\int_{0}^{2\pi}\left|\frac{d\mathbf{r}_{0}}{d\phi}\right|d\phi.$ (20)

Alternatively, due to Eq. 11, $\overline{\eta}$ can be regarded as the
constant parameter that reflects the magnitude of variation in $B$.

We now introduce dimensionless $x(\psi)$ and $y(\alpha)$ coordinates that
correspond to scaled versions of $\psi$ and $\alpha$:

$\displaystyle x$ $\displaystyle=\frac{dx}{d\psi}(\psi-\psi_{0}),$ (21)

and

$\displaystyle y$ $\displaystyle=\frac{dy}{d\alpha}(\alpha-\alpha_{0}),$ (22)

with $dx/d\psi$ and $dy/d\alpha$ constant. Here, we take
$\psi_{0}=\alpha_{0}=0$ and define $dx/d\psi$ and $dy/d\alpha$ such that $x$
reduces to the usual minor radius in the cylindrical limit and that the
product $xdx/d\psi$ equals the inverse of a reference magnetic field,
$B{\textsubscript{ref}}$. More specifically, we take $x=ar/r_{a}$ and
$y=\alpha ar/r_{a}$ with $r_{a}$ the value of $r$ at the VMEC boundary and
$a=$Aminor_p the effective minor radius. The linearized gyrokinetic equation,
Eq. 5, is then Fourier decomposed along the $x$ and $y$ directions with
$\mathbf{k}_{\perp}=k_{x}\nabla x+k_{y}\nabla y$ the perpendicular wave-
vector.

An example first-order near-axis stellarator geometry is shown in Fig. 1,
namely the quasi-axisymmetric design of Section 5.1 of Ref. [7]. A field line
with $r=0.15$, $\alpha=0$ and $-2\pi<\vartheta<2\pi$ is shown, together with
the location of $\vartheta=0$. A poloidal projection of a longer field-line
with $-10\pi<\vartheta<10\pi$ and ten poloidal cross-sections of the flux
surface are shown in the top-right of Fig. 1. The flux-tube approach employed
in GS2 to simulate plasma conditions in the field-line following limit can
then thought of as a series of field-lines similar to the one in Fig. 1 with a
finite (but small) extension along the radial ($\psi$) and poloidal ($\alpha$)
directions.

Figure 1: Top left: top view of a first-order quasi-axisymmetric surface of
constant toroidal flux $\psi$ (Section 5.1 of Ref. [7]). Colour indicates $B$
on the outermost flux surface. An example field-line with
$-2\pi<\vartheta<2\pi$, $\alpha=0$ and $r=0.15$ is shown in black together
with the location of $\vartheta=0$. Top right: ten poloidal planes and their
respective axis location, together with a field line with
$-10\pi<\vartheta<10\pi$ and $\alpha=0$. Bottom: side view of the toroidal
flux surface.

The stellarator geometries used in this study, apart from a near-axis
benchmark case, consist of eleven quasisymmetric stellarator designs. These
were developed by several independent research teams using different
optimization codes, with varying degrees of quasisymmetry. By order of
decreasing number of field periods, these are the NZ1988 design of Ref. [38],
the Drevlak design of Ref. [39] recently developed at the Max Planck Institute
for Plasma Physics in Greifswald, Germany, HSX [40], the KuQHS48 design of
Ref. [41], the WISTELL-A configuration of Ref. [42] recently developed at the
University of Wisconsin-Madison, USA, NCSX (configuration LI383, [43]), ARIES-
CS configuration N3ARE [44], the QAS2 configuration of Ref. [45] with a
vanishing on-axis current, ESTELL [46], CFQS [47] and the Henneberg design of
Ref. [48]. The properties of each stellarator design are listed in Table 1 of
Ref. [26], together with the $B_{0}$ and $\overline{\eta}$ constants used for
each configuration and a comparison between the original and the resulting
near-axis geometric coefficients. The original coefficients $\mathbf{Q}$ are
obtained from the VMEC file corresponding to each design using the geometry
module of the stella code [49].

Finally, we review the main input and output parameters of the model equations
when combined with the near-axis geometry described above. For each
configuration, we use its axis shape to find the axis curvature $\kappa$ and
torsion $\tau$. The parameters $B_{0}$ and $\overline{\eta}$ are found using
the fitting procedure described in Ref. [26] allowing us to find the $\sigma$
function for each configuration, hence to fully specify the geometry
coefficients in Eqs. 11, 12, 13, 14, 15, 16, 17 and 18. The physical input
parameters used by GS2 to solve the linear gyrokinetic equation, Eq. 5, are
the $x$ and $y$ components of the perpendicular wave-vector, $k_{x}$ and
$k_{y}$, respectively, and the physical parameters related to the particle
species. We note that, at $\alpha=0$, choosing a flux tube with $k_{x}=0$
yields the largest growth rate [50]. Therefore, in the following, all
simulations are performed at $\alpha=k_{x}=0$. As electrons are considered
adiabatic, the input characteristic gradient lengths consist only of the ion
density, $L_{n}$, and temperature, $L_{T}$, lengths. The atomic number is set
to $Z=1$. The temporal and spatial resolution parameters are set using
convergence scans.

##  3 Models of Near-Axis ITG stability

In this section, we use the geometric coefficients $\mathbf{Q}$ in the set of
Eqs. 11, 12, 13, 14, 15, 16, 17 and 18 in order to simplify the characteristic
ITG frequencies $\omega_{*}$ and $\omega_{d}$, together with the perpendicular
wave-vector $\mathbf{k}_{\perp}$. Then, we plug these simplified expressions
into known limits of the gyrokinetic equation commonly employed in ITG
stability studies. The resulting models allow for a more intuitive
understanding of the properties of the ITG mode near the magnetic axis and are
used in the subsequent sections to interpret the linear gyrokinetic results.

At lowest order in the distance to the axis, and noting that
$k_{\alpha}=ark_{y}/r_{a}$, the diamagnetic frequency $\omega_{*}$ can be
written as

$\frac{\omega_{*}}{\Omega}=\frac{a}{L_{n}}\frac{\rho_{th}}{r_{a}}k_{y}\rho_{th},$
(23)

where $\rho_{thi}=v_{th}/\Omega$ and $L_{n}=d\ln n/dr$ (and similarly for
$L_{T}$), while the drift-frequency reads

$\frac{\omega_{d}}{\Omega}=s_{\psi}\frac{a}{r_{a}}\overline{\eta}\rho_{th}k_{y}\rho_{th}\cos\vartheta.$
(24)

In the following, we further assume that the most unstable linear modes peak
at $k_{\psi}=0$ and write the square of the normalized perpendicular wave-
vector $b$ as

$b={k_{\perp}^{2}\rho_{th}^{2}}=\frac{a^{2}}{r_{a}^{2}}k_{y}^{2}\rho_{th}^{2}\left[\frac{\overline{\eta}^{2}}{\kappa^{2}}\cos^{2}\vartheta+\frac{\kappa^{2}}{\overline{\eta}^{2}}(\sigma\cos\vartheta+\sin\vartheta)^{2}\right].$
(25)

Taking $z=\vartheta$, the parallel derivative operator can be written as

$\mathbf{b}\cdot\nabla\phi=\nabla_{\parallel}\phi=s_{G}\frac{2\pi}{L}\iota_{N}(1+r\overline{\eta}\cos\vartheta)\frac{\partial\phi}{\partial\vartheta}.$
(26)

Furthermore, by defining the normalized quantities
$\tilde{k}_{y}=k_{y}\rho_{th}$, $\tilde{L}_{n}=L_{n}2\pi/L$,
$\tilde{\eta}=\overline{\eta}L/2\pi$, $\tilde{\kappa}=\kappa/\overline{\eta}$,
$\tilde{\omega}=\omega/\omega_{t}$ with the transit frequency $\omega_{t}=2\pi
v_{th}/L$ (similarly for $\omega_{d}$ and $\omega_{*}$), we rewrite the ITG
parameters above as

$\tilde{\omega}_{*}=\frac{\tilde{k}_{y}}{\tilde{L}_{n}},$ (27)

and

$\tilde{\omega}_{d}=s_{\psi}\tilde{k}_{y}\tilde{\eta}\cos\vartheta,$ (28)

and

$b=\frac{\tilde{k}_{y}^{2}}{\tilde{\kappa}^{2}}\left[\cos^{2}\vartheta+\tilde{\kappa}^{4}(\sigma\cos\vartheta+\sin\vartheta)^{2}\right],$
(29)

where the approximation $a\simeq r_{a}$ was used. This reduces the set of
dimensionless parameters entering the linearized gyrokinetic equation to
$\tau_{T},\tilde{k}_{y},\tilde{L}_{n},\eta,\tilde{\eta},\tilde{\kappa}$ and
$\sigma$. We note that for stellarator-symmetric designs $\sigma(\varphi=0)=0$
and that according to Eq. 19

$\sigma^{\prime}(0)=\frac{1}{(\iota_{0}-N)}\left.\frac{d\sigma}{d\varphi}\right|_{0}=\frac{2\tilde{\kappa}(0)^{2}[\tilde{I}-s_{\psi}\tilde{\tau}(0)]}{(\iota_{0}-N)^{2}[1+\tilde{\kappa}(0)^{4}]},$
(30)

where $\tilde{I}=I_{2}L/(2\pi B_{0})$ is the normalized current and
$\tilde{\tau}=\tau L/2\pi$ is the normalized torsion.

One limit in which further analytic simplification is possible is the case in
which the eigenfunction is localized near $\vartheta=0$. We can then expand
the ITG parameters $\omega_{d}$ and $b$ around $\vartheta=0$ as

$\tilde{\omega}_{d}=\tilde{\omega}_{d0}\left(1-\frac{\vartheta^{2}}{2}\right),$
(31)

and

$b=b_{0}(1+s\vartheta^{2}),$ (32)

where we find $\tilde{\omega}_{d0}=s_{\psi}\tilde{k}_{y}\tilde{\eta}$,
$b_{0}=\tilde{k}_{y}^{2}/\tilde{\kappa}^{2}$, the effective shear parameter

$s=-\left[1+\frac{\tilde{\kappa}^{\prime\prime}(0)}{(\iota_{0}-N)^{2}\tilde{\kappa}(0)}-\tilde{\kappa}(0)^{4}\left(1+\frac{\sigma^{\prime}(0)}{\iota_{0}-N}{}\right)^{2}\right],$
(33)

###  3.1 Nonlocal small transit frequency approximation in the non-resonant
limit

Following Ref. [51], we consider the limit where the transit frequency
$\omega_{t}$ is small when compared with the mode frequency $\omega$.
Furthermore, we adopt the non-resonant limit where
$\omega/\omega_{*}^{T}=\delta\ll 1$ and
$\omega_{d}/\omega\sim\omega_{t}^{2}/\omega^{2}\sim
k_{\perp}^{2}\rho^{2}\sim\delta$, yielding the following second order
differential equation for the electrostatic potential [52]

$\frac{\omega_{*}^{T}v_{th}^{2}}{2\omega^{3}}\nabla_{\parallel}^{2}\phi=\left[\tau_{T}+(1-b\eta)\frac{\omega_{*}}{\omega}+\frac{\omega_{d}\omega_{*}^{T}}{\omega^{2}}\right]\phi,$
(34)

where $\omega_{*}^{T}=\eta\omega_{*}$ and
$\nabla_{\parallel}^{2}\phi\simeq(\iota_{0}-N)^{2}({4\pi^{2}}/{L^{2}}){\partial^{2}\phi}/{\partial\vartheta^{2}}$.
In the highly localized approximation of Eqs. 31 and 32, Eq. 34 turns into a
Schrödinger equation of the Weber type whose bound solutions have the form
$\phi\sim H_{n}(\sqrt{\delta}\vartheta)\exp(\delta\vartheta^{2}/2)$ where
$H_{n}$ is a Hermite polynomial of order $n=0,1,2,...$ and the width $\delta$
is given by

$\delta^{2}=-\tilde{\omega}\tilde{k}_{y}\left(s_{\psi}\tilde{\eta}+\frac{2\tilde{k}_{y}s^{2}}{\tilde{\kappa}^{2}}\tilde{\omega}\right),$
(35)

where $\omega$ is the root of the dispersion relation

$\frac{\tau_{T}\tilde{L}_{n}}{\tilde{k}_{y}\eta}\tilde{\omega}^{3}+\left(1-\frac{\eta\tilde{k}_{y}^{2}}{\tilde{\kappa}^{2}}\right)\frac{\tilde{\omega}^{2}}{\eta}+s_{\psi}\tilde{k}_{y}\tilde{\eta}\tilde{\omega}+\left(n+\frac{1}{2}\right)\delta^{+}=0,$
(36)

for which Re$[\delta]>0$ (denoted as $\delta^{+}$). We note that, in the
dispersion relation above, the effective shear parameter $s$ only contributes
to $\delta^{2}$ as an $s^{2}$ factor yielding eigenfunctions that are
symmetric with respect to $s$.

An example numerical solution for the dispersion relation in Eq. 36 with
$s=0.1,\tilde{\eta}=1$ and $\tilde{\kappa}=1$ is shown in Fig. 2. There, the
peak growth rate $\gamma$ and the corresponding value of $\tilde{k}_{y}$ and
$\delta^{+}$ are shown for a range of $L_{n}$ and $L_{T}$ values satisfying
$\eta=L_{n}/L_{T}>1$. It is seen that, in general, the ITG growth rate
$\gamma$ and its corresponding mode width $\delta^{+}$ increase with
increasing $a/L_{n}$ and $a/L_{T}$ while the opposite is found for $k_{y}$.

Figure 2:  Solution of the ITG dispersion relation Eq. 36 with $\eta>1$ for
the case $s=0.1,\tilde{\eta}=1$ and $\tilde{\kappa}=1$, namely the peak growth
rate $\gamma$ (top left), the corresponding $k_{y}\rho_{th}$ (top right), real
frequency $\omega$ (bottom left) and the eigenfunction width $\delta^{+}$
(bottom right). The frequencies $\gamma$ and $\omega$ are normalized to the
transit frequency $\omega_{t}=2\pi v_{th}/L$.

Besides a scan over $a/L_{n}$ and $a/L_{T}$, the dispersion relation in Eq. 36
can be used to assess the impact of $s$ and $\eta_{b}$ (which are quantities
directly related to the near-axis expansion parameters) in the stability of
the ITG mode. In Fig. 3 a parameter scan over $s$ and $\eta_{b}$ is shown for
a fixed value of $\eta=3$ and $a/L_{n}=1$. It is seen that the results are
symmetric with respect to the sign of $\overline{\eta}$ and that all four
parameters increase as either $s$ or $\overline{\eta}$ are increased.

Figure 3: Solution of the ITG dispersion relation Eq. 36 for the case $\eta=3$
and $a/L_{n}=1$, namely the peak growth rate $\gamma$ (top left), the
corresponding $k_{y}\rho_{th}$ (top right), real frequency $\omega$ (bottom
left) and the eigenfunction width $\delta^{+}$ (bottom right). The frequencies
$\gamma$ and $\omega$ are normalized to the transit frequency $\omega_{t}=2\pi
v_{th}/L$ while $\overline{\eta}$ is normalized to $L/2\pi$. Figure 4:
Effective shear parameter $s$ and $\overline{\eta}$ for each of the eleven
quasisymmetric designs considered in this study computed using Eq. 33. The
designs with hollow circles are quasi-axisymmetric, while the ones with filled
circles are quasi-helically symmetric. The color represents the growth rate
$\gamma$ computed from the dispersion relation in Eq. 36 with
$\tilde{L}_{n}=0.6$ and $\eta=2.5$.

Finally, we apply the dispersion relation in Eq. 36 to the eleven
quasisymmetric designs used in this study. For each design, the value of
$\overline{\eta}$ is read from Table 1 in Ref. [26], the axis parameters are
read from each VMEC file (used to compute the curvature and torsion) and
$\sigma$ is calculated using Eq. 19, yielding the values of $\iota_{0}$ in
Table 1 of Ref. [26]. The values of the effective shear parameter $s$ for each
design is shown in Fig. 4. It is found that quasi-helically symmetric
stellarators have, in general, values of $s$ that are one to two orders of
magnitude above the ones of quasi-axisymmetric stellarators. We also show in
Fig. 4 the resulting growth rate estimate and the corresponding values of
$\overline{\eta}L/2\pi$ and $s$. The designs with $\overline{\eta}L/2\pi<1$
are quasi-axisymmetric, while the ones with $\overline{\eta}L/2\pi>1$ are
quasi-helically symmetric. The results in Fig. 4 suggest that quasi-helically
symmetric configurations are prone to having higher values of
$\overline{\eta}L/2\pi$, lower values of $s$ and lower growth rates
(normalized to $\omega_{t}$) when compared with quasi-axisymmetric designs.

###  3.2 Bounce-averaged metric

We now show how a bounce-averaged metric $Q{\textsubscript{bounce}}$, which
takes into account the role of good and bad curvature regions in the particle
motion, can be simplified using the near-axis expansion. This parameter,
although more focused on TEM rather than ITG, allows us to show how metrics
based on bounce-averaged quantities, in the near-axis expansion formalism, can
be reduced to one-dimensional elliptic integrals. Here, we employ the metric
used in Ref. [53]

$Q{\textsubscript{bounce}}=-\int_{1/B{\textsubscript{max}}}^{1/B{\textsubscript{min}}}\overline{\omega}_{d}(\lambda)d\lambda,$
(37)

where $\overline{\omega}_{d}(\lambda)$ is the bounce-averaged drift-frequency

$\overline{\omega}_{d}(\lambda)=\frac{\int_{l_{1}}^{l_{2}}\tilde{\omega}_{d}{dl}/{v_{\parallel}}}{\int_{l_{1}}^{l_{2}}{dl}/{v_{\parallel}}}.$
(38)

The minimization of this metric aims at distributing a large fraction of
trapped particles over a favorable bounce-averaged curvature. Although Eq. 37
has been used in Ref. [53] to stabilize TEMs, similar metrics are often
employed in the analysis of tokamak and stellarator microstability [54].

In the near-axis expansion formalism, at first order in $r$, the metric
$Q{\textsubscript{bounce}}$ in Eq. 37 can be simplified by noting that
$v_{\parallel}=\pm v\sqrt{1-\lambda B}$ with
$B=B_{0}(1+r\overline{\eta}\cos\vartheta)$. The resulting integral in $l$ is
rewritten by parametrizing the field-line coordinate $l$ using
$d\vartheta/dl=\mathbf{b}\cdot\nabla\vartheta=|\iota_{N}|B/(G+\iota
I)\simeq(2\pi/L)|\iota_{N}|(1+r\overline{\eta}\cos\vartheta)$, yielding the
bounce points $l_{1}$ and $l_{2}=-l_{1}$ at
$\vartheta_{1}=-\vartheta_{2}=-2\arcsin(1/a)$ with
$a^{2}=-2r\overline{\eta}\lambda B_{0}/([1-\lambda
B_{0}(1+r\overline{\eta})]>1$. By performing a variable substitution
$x=a\sin\vartheta/2$ and assuming $s_{\psi}=1$, the following expression for
$Q{\textsubscript{bounce}}$ is found

$Q{\textsubscript{bounce}}=\frac{8}{\sqrt{2}}\frac{k_{y}\rho_{th}\overline{\eta}v^{2}}{v_{th}}\int_{[B_{0}(1+r\overline{\eta})]^{-1}}^{[B_{0}(1-r\overline{\eta})]^{-1}}d\lambda\sqrt{2-\lambda
B_{0}(1+r\overline{\eta})}\frac{\Pi_{c}(1/c,t,1/a)}{\Pi(t,1/a)},$ (39)

with $\Pi_{c}$ and $\Pi$ the elliptic integrals

$\Pi_{c}(1/c,t,1/a)=\int_{0}^{1}\frac{(1-2x^{2}/a^{2})\sqrt{1-x^{2}/c^{2}}}{\sqrt{1-x^{2}}\sqrt{1-x^{2}/a^{2}}}\frac{dx}{1-tx^{2}},$
(40)

and

$\Pi(t,1/a)=\int_{0}^{1}\frac{1}{\sqrt{1-x^{2}}\sqrt{1-x^{2}/a^{2}}}\frac{dx}{1-tx^{2}},$
(41)

respectively. In Eq. 39, $t=2r\overline{\eta}/[(1+r\overline{\eta})a^{2}]$ and
$c^{2}=[2-\lambda B_{0}(1+r\overline{\eta})]/[1-\lambda
B_{0}(1+r\overline{\eta})]$. Noting that the integrand of Eq. 39 is a function
of $r\overline{\eta}$ and $\lambda B_{0}$ only, we rewrite
$Q{\textsubscript{bounce}}$ by turning the integral over $\lambda$ to an
integral over $\lambda B_{0}$, yielding.

$Q{\textsubscript{bounce}}=\frac{8}{\sqrt{2}}\frac{k_{y}\rho_{th}\overline{\eta}v^{2}}{v_{th}B_{0}}q{\textsubscript{bounce}}(r\overline{\eta}).$
(42)

As shown in Fig. 5, the function $q{\textsubscript{bounce}}$ can be reasonably
approximated by a quadratic function of $r\overline{\eta}$, leading to the
following approximate model for $Q{\textsubscript{bounce}}$ close to the
magnetic axis

$Q{\textsubscript{bounce}}\simeq
5.051\frac{k_{y}\rho_{th}\overline{\eta}v^{2}}{B_{0}v_{th}}r\overline{\eta}[1-1.172(r\overline{\eta})^{2}].$
(43)

The positive sign of $Q{\textsubscript{bounce}}$ in Eq. 43 suggests that, for
quasisymmetric stellarators, the majority of trapped particles experience bad
curvature. Furthermore, Eq. 43 shows that there is little freedom to optimize
quasi-symmetric stellarators for trapped electron modes besides varying the
value of $\overline{\eta}$. Finally, an estimate for $Q\textsubscript{bounce}$
for the eleven quasisymmetric designs considered in this study using Eq. 43 is
shown in Fig. 5. We note that in Fig. 5 there are are no systematic
differences between quasi-axisymmetric and quasi-helically symmetric
stellarators.

Figure 5: Left: Analytical form and the corresponding quadratic fit
$0.893x-1.53x^{2}$ for the function $q{\textsubscript{bounce}}$ entering the
heat flux proxy $Q\textsubscript{bounce}$ in Eq. 42. Right:
$Q\textsubscript{bounce}$ normalized to $k_{y}\rho_{th}v^{2}/B_{0}v_{th}$ as a
function of $\overline{\eta}L/2\pi$ for the eleven configurations considered
in this study and $r=0.1$. The designs with hollow circles are quasi-
axisymmetric, while the ones with filled circles are quasi-helically
symmetric.

###  3.3 Reduction of Turbulent Transport via Quasilinear Modelling

Using a quasilinear model, in Ref. [14, 18] the following proxy function
$Q{\textsubscript{prox}}$ for ITG transport was derived

$Q{\textsubscript{prox}}=-\frac{c_{D}}{4\overline{B}}\frac{\gamma}{k_{\psi}^{2}}\frac{|\nabla\psi|^{2}}{\psi^{2}}\frac{T}{L_{T}},$
(44)

with $c_{D}$ a constant, $\overline{B}$ a reference magnetic field and
$\gamma$ the growth rate resulting from the following simplified ITG
dispersion relation

$\gamma\simeq\omega_{*}L_{n}|\tau_{T}\kappa_{1}(\kappa_{p}-\kappa{\textsubscript{cr}})|^{1/2}H(\kappa_{p}-\kappa{\textsubscript{cr}})H(-\kappa_{1}),$
(45)

where $\kappa_{1}$ is the radial component of the curvature vector
$\mathbf{b}\cdot\nabla\mathbf{b}$, $\kappa_{p}=(1+\eta)/L_{n}$ is the pressure
gradient length and $\kappa{\textsubscript{cr}}$ a critical pressure gradient
length stemming from a contribution of the ITG slab branch. Minimization of
this proxy was shown to be an effective way to reduce turbulent transport
levels, raising the prospect of a new class of stellarators with improved
overall confinement [18].

Using the expression for $|\nabla\psi|^{2}$ in Eq. 13 and choosing
$\overline{B}=B_{0}$, we can write the proxy function in Eq. 44 as

$Q{\textsubscript{prox}}=Q_{0}\frac{\sqrt{\overline{\eta}\cos\vartheta}}{B_{0}}H(\cos\vartheta)\left[\frac{\overline{\eta}^{2}}{\kappa^{2}}\sin^{2}\vartheta+\frac{\kappa^{2}}{\overline{\eta}^{2}}(\cos\vartheta-\sigma\sin\vartheta)^{2}\right]$
(46)

where

$Q_{0}=-\frac{c_{D}}{k_{\psi}^{2}}\frac{T}{L_{T}}\frac{k_{y}\rho_{th}c_{s}}{a}\sqrt{\frac{\tau_{T}}{L_{n}}}\sqrt{1+\eta+\kappa{\textsubscript{cr}}L_{n}}H(1+\eta+\kappa{\textsubscript{cr}}L_{n}).$
(47)

The function $Q{\textsubscript{prox}}$ in Eq. 46 can then be used as a proxy
for the reduction of ITG turbulence near the axis of quasisymmetric
stellarators.

Figure 6: Normalized proxy function $Q{\textsubscript{prox}}/Q_{0}$ computed
using Eq. 46 for each of the eleven quasisymmetric designs considered in this
study averaged along a field line with $\alpha=0$. The designs with hollow
circles are quasi-axisymmetric, while the ones with filled circles are quasi-
helically symmetric.

In Fig. 6, the normalized proxy function $Q{\textsubscript{prox}}/Q_{0}$ is
evaluated for the configurations considered in this study, averaged along a
field line with $\alpha=0$. Due to the fact that the parameter $Q_{0}$ in Eq.
47 is not directly dependent on the geometry used, it is considered constant
in Fig. 6. The resulting averaged proxy function in Fig. 6 suggests that there
is not a significant difference in $Q$ between quasi-axisymmetric and quasi-
helically symmetric configurations. Indeed, the range of values of
$Q\textsubscript{prox}$ in Fig. 6 can vary two orders of magnitude in each set
of quasisymmetric configurations. We note that this analysis assumes $Q_{0}$
is constant across all stellarator designs, which may not be true in practice
due to possible variations of characteristic wave-vectors and critical
gradients between configurations.

##  4 Numerical Benchmark with a Near-Axis Geometry

Figure 7: Top left: top view of a surface of constant toroidal flux $\psi$ for
the quasi-helically symmetric NAQS stellarator (Section 5.4 of Ref. [55]) at
first order in $\epsilon$. Colour indicates $B$ on the outermost flux surface.
An example field-line with $-8\pi<\vartheta<8\pi$, $\alpha=0$ and $r=0.03$ is
shown in black together with the location of $\vartheta=0$. Top right: ten
poloidal planes and their respective axis location, together with a field line
with $-40\pi<\vartheta<40\pi$ and $\alpha=0$. Bottom: side view of the
toroidal flux surface.

Before assessing the linear stability of the eleven quasisymmetric designs and
its comparison with the near-axis expansion, we start with a comparison
between a solution of the near-axis equations and a VMEC configuration based
on this solution. This allows us to examine the details of the input
parameters and respective convergence tests, the comparison between
eigenfunctions, the dependence of the growth rate on the perpendicular wave-
vector and the variation of the maximum growth rate with the density and
temperature gradient scale lengths before presenting the main results in the
next section. For this analysis, we take the quasi-helically symmetric
solution of Section 5.4 in Ref. [56] (hereafter referred to as NAQS) and
generate a corresponding VMEC output file with a plasma boundary at
$\psi=0.002$ T m2, $B_{0}=1.0$ T and $\overline{\eta}=2.25$, with a resulting
rotational transform on-axis of $\iota_{0}=1.93$ when solving Eq. 19. The
comparison between the geometrical quantities is shown in Fig. 8 using flux
surfaces with $s=\psi/\psi_{a}$ (with $\psi_{a}$ the toroidal flux at the
plasma boundary) ranging from 0.01 to 0.8 and a field line with $\alpha=0$.
Only six of the eight quantities are shown as Eqs. 11 and 12 and Eqs. 16 and
18 are independent in the limit of near-axis quasisymmetric fields used here.
Although, overall, there is good agreement, some differences between the two
methods arise, even at $s=0.01$. The fact that a magnetic field solution
obtained from setting the boundary value of $\mathbf{B}$ from a first-order
near-axis solution at a particular radius $r$ is not perfectly quasisymmetric
is a phenomenon that has been recently been observed in Ref. [56] and
addressed using second-order solutions in Ref. [55].

Figure 8: Comparison between the fitted and the original near-axis geometry
coefficients for the NAQS configuration of Section 5.4 of Ref [56] along the
cylindrical toroidal angle $\phi$ at several radial locations
$s={\psi/\psi_{a}}$ with $\psi_{a}$ the toroidal flux at the plasma boundary.
All quantities are expressed in SI units.

We proceed by using the GS2 code to assess the properties of ITG modes in both
the VMEC and the near-axis solution. For each simulation, the temporal
evolution of $|\hat{\phi}|^{2}$ is fitted to an exponential of the form
$e^{2\gamma t}$ with $\gamma$ the growth rate. As a base case scenario, we use
ngauss=3 (half of the number of untrapped pitch-angles moving in one direction
along field line), deltat=0.4 (the temporal increment between time steps),
nstep=150 (the total number of time steps), ngrid=10 (the total number of
energy points), nzgrid=100 (the total number of points in $z$ along the field
line), nlambda=24 (number of trapped-pitch angles along the field line) and
nfp=4 (number of field periods, a measure of the total length of the field
line). In Fig. 9, a convergence test is performed where the base case growth
rate for both the VMEC and the near-axis geometries is computed, along with
the growth rate resulting from doubling the resolution along each parameter.
As physical input parameters, we use $k_{y}\rho=1.0$, $a/L_{T}=3$ and
$a/L_{n}=1$. A base case scenario is considered to be converged if the scatter
plot cluster of all simulations has a deviation smaller than 5%, a criterion
that is satisfied for all radii in Fig. 9.

Figure 9: Growth rate convergence test of the base case scenario for the VMEC
and near-axis NAQS geometries using the gyrokinetic GS2 code at different
radii $s=\psi/\psi_{a}$ and $k_{y}\rho=1$, $a/L_{T}=3$ and $a/L_{n}=1$.

For the base case scenario at $s=0.3$, we show in Fig. 10 the resulting real
and imaginary electrostatic potential eigenfunctions $\hat{\phi}$ at the end
of the simulation period for both the VMEC (left) and near-axis (right)
geometries. Overall, both eigenfunctions look practically indistinguishable
from each other for every simulation of the convergence test. Finally, we note
that these are smooth functions of the parallel coordinate $z$ and reach
insignificant values before the end of the simulation domain, which provides
additional evidence for the identification of a converged simulation scenario.

Figure 10: Real and imaginary eigenfunctions $\hat{\phi}$ at the end of the
simulation period for both the VMEC (left) and near-axis (right) geometries at
$s=0.3$. Each eigenfunction is normalized to its value at $z=0$.

Next, we perform a scan over $k_{y}$ in order to determine the peak growth
rate $\gamma$ for each value of $L_{n}$ and $L_{T}$. For the values of
$a/L_{T}=5$ and $a/L_{n}=2$ considered above, we show in Fig. 11 the
dependence of $\gamma=\gamma(k_{y})$ and $\omega=\omega(k_{y})$ for $0\leq
k_{y}\rho\leq 20$ for both configurations at $s=0.3$. The peak growth rate is
found at $k_{y}\rho=2$, a value within the expected range of the ITG
instability. Also, the near-axis expansion appears to yield a larger value of
the growth rate across the majority of the spectrum, with a particular
incidence at the peak values of $k_{y}\sim 1$.

Figure 11: Dependence of $\gamma=\gamma(k_{y})$ and $\omega=\omega(k_{y})$ for
$0\leq k_{y}\rho\leq 20$ for the VMEC and near-axis NAQS configuration for
$a/L_{T}=5$ and $a/L_{n}=2$ at $s=0.3$.

Finally, we vary the density and temperature gradient scale lengths, $a/L_{n}$
and $a/L_{T}$, and pick the fastest growing mode by selecting the maximum
value of $\gamma(k_{y})$ for each $L_{n}$ and $L_{T}$. The resulting peak
growth rates, associated $\omega$ and $k_{y}$ within $0\leq a/L_{n}\leq 6$ and
$0\leq a/L_{T}\leq 6$ are shown in Fig. 12. The locations of $\omega$ and
$k_{y}$ where $\gamma\leq 0$ are shown in white. As is characteristic of the
ITG mode, the peak growth rate is negative for a range of values where
$\eta=L_{n}/L_{T}>\eta{\textsubscript{crit}}$ with $\eta_{crit}\sim 1$ [57],
such as the condition for instability $\eta>2/3$ of the toroidal branch of the
ITG mode found in [58]. A more detailed understanding of the mechanisms that
set the ITG critical gradient using first-principle approaches, such as the
one in Ref. [59], are left for future work.

Figure 12: Comparison of the peak growth rates $\gamma$ and associated real
frequencies $\omega$ and wave-numbers $k_{y}$ between the VMEC and near-axis
geometries for the NAQS stellarator at $s=0.01$. The gradient lengths are
scanned over $0\leq a/L_{n}\leq 6$ and $0\leq a/L_{T}\leq 6$.

##  5 Numerical results

In this section, we carry out the analysis performed in Section 4 for a single
benchmark case, now for the eleven realistic quasisymmetric designs under
study. In particular, we perform a convergence scan to determine a base case
scenario that satisfies the convergence criterion for all stellarator designs,
perform a scan over $k_{y}$, $L_{n}$ and $L_{T}$ and pick the fastest growth
rate in order to assess the critical gradients $\eta{\textsubscript{crit}}$
and peak growth rates $\gamma{\textsubscript{max}}$.

We start with the convergence scan. A comparison between the growth rate
obtained with the near-axis expansion and VMEC is shown in Fig. 13 for several
radii $s=\psi/\psi_{a}=(0.01,0.05,0.1,0.3,0.5,0.8)$ and $k_{y}\rho=1$,
$a/L_{T}=3$ and $a/L_{n}=1$. In order to obtain converged base case scenarios,
the spatial resolution was increased to nzgrid=150, nlambda increased to
nlambda=28 and, due to the presence of lower values of $\gamma$ with respect
to the NAQS benchmark, the simulation time was increased to nstep=250.
Concerning the top row in Fig. 13 (radii with $s<0.1$), the growth rate
$\gamma$ obtained with a near-axis model is found to be in very good agreement
with the ones using the corresponding VMEC design. However, it is found that
for larger radii, the near-axis expansion provides an overestimation of the
growth rate. This behaviour is observed for all cases studied here.

Figure 13:  Growth rate convergence test of the base case scenario for the
VMEC and near-axis geometries using the gyrokinetic GS2 code at different
radii $s=\psi/\psi_{a}$ and $k_{y}\rho=1$, $a/L_{T}=3$ and $a/L_{n}=1$. Quasi-
axisymmetric stellarators are represented with hollow circles while quasi-
helically symmetric ones are represented with filled circles.

We distinguish between the quasi-axisymmetric and quasi-helically symmetric
configurations in Fig. 13 by using both hollow and filled circles,
respectively. It is seen that, in general, there is no systematic difference
between the two for the physical parameters simulated here, except for a small
tendency of quasi-axisymmetric stellarators close to the magnetic axis to
exhibit higher growth rates. However, this tendency seems to reverse further
away from the axis where, at $s=0.8$, the majority of the quasi-axisymmetric
configurations have considerably lower growth rates than quasi-helically
symmetric ones, widening the gap between the maximum and minimum values of the
computed $\gamma$.

Figure 14: Horizontal axis: growth rate estimate $\gamma$ normalized to
$\omega_{t}$ computed using the dispersion relation in Eq. 36 (denoted as
$\gamma_{DR}$). Vertical axis: growth rate estimate $\gamma$ normalized to
$\omega_{t}$ computed using the gyrokinetic GS2 code and the corresponding
VMEC equilibrium files (denoted as $\gamma_{GS2}$, corresponding to the values
in Fig. 13 multiplied by $L/2\pi a$). Quasi-axisymmetric stellarators are
represented with hollow circles while quasi-helically symmetric ones are
represented with filled circles and $r=0.01$.

We now comment on the accuracy of the estimate using a simplified dispersion
relation in Eq. 36 when compared with the growth rates using first-principle
gyrokinetic simulations in Fig. 13. This comparison is shown in Fig. 14, where
we show the growth rate estimate $\gamma$ normalized to $\omega_{t}$ computed
using the dispersion relation in Eq. 36 (denoted as $\gamma_{DR}$) on the
horizontal axis and the growth rate estimate $\gamma$ normalized to
$\omega_{t}$ computed using the gyrokinetic GS2 code and the corresponding
VMEC equilibrium files (denoted as $\gamma_{GS2}$, corresponding to the values
in Fig. 13 multiplied by $L/2\pi a$). It is seen that the estimate
$\gamma_{DR}$ resulting from the simplified dispersion relation in Eq. 36
consistently overestimates the growth rate which, in some cases, can lead to a
one order of magnitude difference between the two growth rates for quasi-
axisymmetric stellarators with higher values of $\gamma$. However, the general
trend and the relative growth rates between different configurations is well
predicted by the simplified dispersion relation.

Finally, we assess the accuracy of the near-axis expansion at predicting the
peak growth rates when performing a scan over $a/L_{n}$ and $a/L_{T}$. Two
examples are shown in Fig. 15 and Fig. 16 for the NCSX and HSX stellarators,
respectively. These confirm the general trend observed in Fig. 12 and the
accuracy of the near-axis expansion for values of $s<0.3$. In general, the
linear critical gradient for the ITG instability near the axis of
quasisymmetric stellarators lies at values of $\eta\sim 1$. We also find that
the peak growth rate increases with increasing $a/L_{T}$ and is maximum for
values of $2<a/L_{n}<4$. The values of the real frequency $\omega$ at the peak
growth rate increase with increasing $a/L_{T}$ while peaking at either
$a/L_{n}=0$ or $a/L_{n}\simeq 3$. The wave-vector $k_{y}$ at the peak growth
rate is observed to increase mainly with increasing $a/L_{n}$.

Figure 15: Comparison of the peak growth rates $\gamma$ and associated real
frequencies $\omega$ and wave-numbers $k_{y}$ between the VMEC and near-axis
geometries for the NCSX stellarator at $s=0.01$. The gradient lengths are
scanned over $0\leq a/L_{n}\leq 6$ and $0\leq a/L_{T}\leq 6$.

Figure 16: Comparison of the peak growth rates $\gamma$ and associated real
frequencies $\omega$ and wave-numbers $k_{y}$ between the VMEC and near-axis
geometries for the HSX stellarator at $s=0.01$. The gradient lengths are
scanned over $0\leq a/L_{n}\leq 6$ and $0\leq a/L_{T}\leq 6$.

##  6 Conclusions

In this work, the ion-temperature-gradient stability of quasisymmetric
stellarators is assessed near the magnetic axis by employing the gyrokinetic
approximation together with a near-axis expansion framework. It is found that
the near-axis expansion is not only able to simplify analytical estimates of
ITG growth rates and turbulent fluxes, but also able to reproduce the
properties of eleven quasisymmetric configurations designed using numerical
optimization methods. In general, as the radius is increased beyond the inner
core, it is found that the near-axis model overestimates the growth rates when
compared with the corresponding optimized VMEC design. Such overestimation
could arise due to the absence of global magnetic shear
$\hat{s}=-(r/\iota)\iota^{\prime}(r)$ in the near-axis expansion formalism at
the order considered here which could contribute as a stabilizing factor at
larger radii where higher order geometry effects are more relevant. Global
magnetic shear may also provide additional finite-Larmor radius stabilization
to ITG modes by contributing with additional terms to the perpendicular wave-
vector $\mathbf{k}_{\perp}$. Furthermore, the presence of inhomogeneities in
the magnetic field of the VMEC configurations, which could lead to departures
from quasisymmetry, may allow the presence of a larger fraction of trapped
particles at the good curvature side and, therefore, decrease $\gamma$.

While the analysis performed here pertained the ITG instability only, a
similar analysis can be carried out to other instabilities, such as the
trapped-electron mode, the electron-temperature-gradient mode and the micro-
tearing mode. Furthermore, as future work, we intend to generalize the
numerical methods employed here for quasisymmetric designs to maximum-J
devices [60] and have a more complete understanding of the supression of
electrostatic microinstabilities recently found in such configurations [61].

##  7 Acknowledgments

We wish to thank I. G. Abel, W. Dorland, G. T. Roberg-Clark and G. G. Plunk
for many fruitful discussions regarding the content of the manuscript. This
work was supported by a grant from the Simons Foundation (560651, ML) and from
the U.S. Department of Energy, Office of Science, Office of Fusion Energy
Science, under award number DE-FG02-93ER54197.

## References

## References

  * [1] P. Xanthopoulos and F. Jenko.  Gyrokinetic analysis of linear microinstabilities for the stellarator Wendelstein 7-X.  Physics of Plasmas, 14(4):042501, 2007.
  * [2] P. Helander, T. Bird, F. Jenko, R. Kleiber, G. G. Plunk, J. H.E. Proll, J. Riemann, and P. Xanthopoulos.  Advances in stellarator gyrokinetics.  Nuclear Fusion, 55(5):053030, 2015.
  * [3] J. Nuhrenberg and R. Zille.  Quasi-helically symmetric toroidal stellarators.  Physics Letters A, 129(2):113–117, 1988.
  * [4] A. H. Boozer.  Quasi-helical symmetry in stellarators.  Plasma Physics and Controlled Fusion, 37(11A):A103, 1995.
  * [5] P. R. Garabedian.  Stellarators with the magnetic symmetry of a tokamak.  Physics of Plasmas, 3(7):2483, 1996.
  * [6] D. A. Garren and A. H. Boozer.  Existence of quasihelically symmetric stellarators.  Physics of Fluids B, 3(10):2822, 1991.
  * [7] M. Landreman and W. Sengupta.  Direct construction of optimized stellarator shapes. Part 1. Theory in cylindrical coordinates.  Journal of Plasma Physics, 84(6):905840616, 2018.
  * [8] R. Jorge, W. Sengupta, and M. Landreman.  Construction of quasisymmetric stellarators using a direct coordinate approach.  Nuclear Fusion, 60(7):076021, 2020.
  * [9] G. G. Plunk.  Perturbing an axisymmetric magnetic equilibrium to obtain a quasi-axisymmetric stellarator.  Journal of Plasma Physics, 86(4):905860409, 2020.
  * [10] B. E. Nelson, L. A. Berry, A. B. Brooks, M. J. Cole, J. C. Chrzanowski, H. M. Fan, P. J. Fogarty, P. L. Goranson, P. J. Heitzenroeder, S. P. Hirshman, G. H. Jones, J. F. Lyon, G. H. Neilson, W. T. Reiersen, D. J. Strickler, and D. E. Williamson.  Design of the national compact stellarator experiment (NCSX).  In Fusion Engineering and Design, volume 66-68, page 169, 2003.
  * [11] J. M. Canik, D. T. Anderson, F. S.B. Anderson, K. M. Likin, J. N. Talmadge, and K. Zhai.  Experimental demonstration of improved neoclassical transport with quasihelical symmetry.  Physical Review Letters, 98(8):085002, 2007.
  * [12] M. Drevlak, C. D. Beidler, J. Geiger, P. Helander, and Y. Turkin.  Optimisation of stellarator equilibria with ROSE.  Nuclear Fusion, 59(1):016010, 2019.
  * [13] P. Helander, J. H.E. Proll, and G. G. Plunk.  Collisionless microinstabilities in stellarators. I. Analytical theory of trapped-particle modes.  Physics of Plasmas, 20(12):122505, 2013.
  * [14] H. E. Mynick, N. Pomphrey, and P. Xanthopoulos.  Reducing turbulent transport in toroidal configurations via shaping.  Physics of Plasmas, 18(5):056101, 2011.
  * [15] G. Rewoldt, W. M. Tang, and M. S. Chance.  Electromagnetic kinetic toroidal eigenmodes for general magnetohydrodynamic equilibria.  Physics of Fluids, 25(3):480, 1982.
  * [16] G. Rewoldt, W. M. Tang, and R. J. Hastie.  Collisional effects on kinetic electromagnetic modes and associated quasilinear transport.  Physics of Fluids, 30(3):807, 1987.
  * [17] P. Xanthopoulos, F. Merz, T. Görler, and F. Jenko.  Nonlinear gyrokinetic simulations of ion-temperature-gradient turbulence for the optimized wendelstein 7-X stellarator.  Physical Review Letters, 99(3):035002, 2007.
  * [18] H. E. Mynick, N. Pomphrey, and P. Xanthopoulos.  Optimizing stellarators for turbulent transport.  Physical Review Letters, 105(9):095004, 2010.
  * [19] W. Dorland, F. Jenko, M. Kotschenreuther, and B. N. Rogers.  Electron Temperature Gradient Turbulence.  Physical Review Letters, 85(26):5579, 2000.
  * [20] A. Zocco, P. Xanthopoulos, H. Doerk, J. W. Connor, and P. Helander.  Threshold for the destabilisation of the ion-temperature-gradient mode in magnetically confined toroidal plasmas.  Journal of Plasma Physics, 84(1):715840101, 2018.
  * [21] P. Xanthopoulos, G. G. Plunk, A. Zocco, and P. Helander.  Intrinsic turbulence stabilization in a stellarator.  Physical Review X, 6(2):021033, 2016.
  * [22] C. C. Hegna, P. W. Terry, and B. J. Faber.  Theory of ITG turbulent saturation in stellarators: Identifying mechanisms to reduce turbulent transport.  Physics of Plasmas, 25(2):022511, 2018.
  * [23] I. J. McKinney, M. J. Pueschel, B. J. Faber, C. C. Hegna, J. N. Talmadge, D. T. Anderson, H. E. Mynick, and P. Xanthopoulos.  A comparison of turbulent transport in a quasi-helical and a quasi-axisymmetric stellarator.  Journal of Plasma Physics, 85(5):905850503, 2019.
  * [24] H. Y. Wang, I. Holod, Z. Lin, J. Bao, J. Y. Fu, P. F. Liu, J. H. Nicolau, D. Spong, and Y. Xiao.  Global gyrokinetic particle simulations of microturbulence in W7-X and LHD stellarators.  Physics of Plasmas, 27(8):082305, 2020.
  * [25] M. Landreman.  Optimized quasisymmetric stellarators are consistent with the Garren–Boozer construction.  Plasma Physics and Controlled Fusion, 61(7):075001, 2019.
  * [26] R. Jorge and M. Landreman.  The Use of Near-Axis Magnetic Fields for Stellarator Turbulence Simulations.  Plasma Physics and Controlled Fusion, 63(1):014001, 2020.
  * [27] B. J. Faber, M. J. Pueschel, P. W. Terry, C. C. Hegna, and J. E. Roman.  Stellarator microinstabilities and turbulence at low magnetic shear.  Journal of Plasma Physics, 84(5):905840503, 2018.
  * [28] M. F. Martin, M. Landreman, P. Xanthopoulos, N. R. Mandell, and W. Dorland.  The parallel boundary condition for turbulence simulations in low magnetic shear devices.  Plasma Physics and Controlled Fusion, 60(9):095008, 2018.
  * [29] J. A. Baumgaertel, E. A. Belli, W. Dorland, W. Guttenfelder, G. W. Hammett, D. R. Mikkelsen, G. Rewoldt, W. M. Tang, and P. Xanthopoulos.  Simulating gyrokinetic microinstabilities in stellarator geometry with GS2.  Physics of Plasmas, 18(12):122301, 2011.
  * [30] G. G. Howes, S. C. Cowley, W. Dorland, G. W. Hammett, E. Quataert, and A. A. Schekochihin.  Astrophysical Gyrokinetics: Basic Equations and Linear Theory.  The Astrophysical Journal, 651(1):590, 2006.
  * [31] E. A. Frieman and L. Chen.  Nonlinear gyrokinetic equations for low-frequency electromagnetic waves in general plasma equilibria.  Physics of Fluids, 25(3):502, 1982.
  * [32] S. P. Hirshman and J. C. Whitson.  Steepest-descent moment method for three-dimensional magnetohydrodynamic equilibria.  Physics of Fluids, 26(12):3553, 1983.
  * [33] A. H. Boozer.  Plasma equilibrium with rational magnetic surfaces.  Physics of Fluids, 24(11):1999, 1981.
  * [34] D. A. Garren and A. H. Boozer.  Magnetic field strength of toroidal plasma equilibria.  Physics of Fluids B, 3(10):2805, 1991.
  * [35] C. Mercier.  Equilibrium and stability of a toroidal magnetohydrodynamic system in the neighbourhood of a magnetic axis.  Nuclear Fusion, 4(3):213, 1964.
  * [36] L. S. Solov’ev and V. D. Shafranov.  Reviews of Plasma Physics 5.  Consultants Bureau, New York - London, 1970.
  * [37] R. Jorge, W. Sengupta, and M. Landreman.  Near-axis expansion of stellarator equilibrium at arbitrary order in the distance to the axis.  Journal of Plasma Physics, 86(1):905860106, 2020.
  * [38] J. Nuhrenberg and R. Zille.  Quasi-helically symmetric toroidal stellarators.  Physics Letters A, 129(2):113, 1988.
  * [39] M. Drevlak.  Stellarator configuration design using ROSE.  In 21st International Stellarator-Heliotron Workshop, page 2, 2017\.
  * [40] F. S. B. Anderson, A. F. Almagri, D. T. Anderson, P. G. Matthews, J. N. Talmadge, and J. L. Shohet.  The Helically Symmetric Experiment, (HSX) Goals, Design and Status.  Fusion Technology, 27(3T):273, 1995.
  * [41] L. P. Ku and A. H. Boozer.  New classes of quasi-helically symmetric stellarators.  Nuclear Fusion, 51(1):013004, 2011.
  * [42] Aaron Bader.  Dataset for Wistell-A stellarator: zenodo.org/record/3952862#.Xxm8OpNKgWo, 2020.
  * [43] M. C. Zarnstorff, L. A. Berry, A. Brooks, E. Fredrickson, G. Y. Fu, S. Hirshman, S. Hudson, L. P. Ku, E. Lazarus, D. Mikkelsen, D. Monticello, G. H. Neilson, N. Pomphrey, A. Reiman, D. Spong, D. Strickler, A. Boozer, W. A. Cooper, R. Goldston, R. Hatcher, M. Isaev, C. Kessel, J. Lewandowski, J. F. Lyon, P. Merkel, H. Mynick, B. E. Nelson, C. Nuehrenberg, M. Redi, W. Reiersen, P. Rutherford, R. Sanchez, J. Schmidt, and R. B. White.  Physics of the compact advanced stellarator NCSX.  Plasma Physics and Controlled Fusion, 43(12A):A237, 2001.
  * [44] F. Najmabadi, A. R. Raffray, S. I. Abdel-Khalik, L. Bromberg, L. Crosatti, L. El-Guebaly, P. R. Garabedian, A. A. Grossman, D. Henderson, A. Ibrahim, T. Ihli, T. B. Kaiser, B. Kiedrowski, L. P. Ku, J. F. Lyon, R. Maingi, S. Malang, C. Martin, T. K. Mau, B. Merrill, R. L. Moore, R. J. Peipert, D. A. Petti, D. L. Sadowski, M. Sawan, J. H. Schultz, R. Slaybaugh, K. T. Slattery, G. Sviatoslavsky, A. Turnbull, L. M. Waganer, X. R. Wang, J. B. Weathers, P. Wilson, J. C. Waldrop, M. Yoda, and M. Zarnstorff.  The ARIES-CS compact stellarator fusion power plant.  Fusion Science and Technology, 54(3):655, 2008.
  * [45] P. R. Garabedian.  Three-dimensional analysis of tokamaks and stellarators.  Proceedings of the National Academy of Sciences of the United States of America, 105(37):13716, 2008.
  * [46] M. Drevlak, F. Brochard, P. Helander, J. Kisslinger, M. Mikhailov, C. Nührenberg, J. Nührenberg, and Y. Turkin.  ESTELL: A Quasi-Toroidally Symmetric Stellarator.  Contributions to Plasma Physics, 53(6):459, 2013.
  * [47] A. Shimizu, H. Liu, M. Isobe, S. Okamura, S. Nishimura, C. Suzuki, Y. Xu, X. Zhang, B. Liu, J. Huang, X. Wang, H. Liu, and C. Tang.  Configuration property of the Chinese first quasi-axisymmetric stellarator.  Plasma and Fusion Research, 13:3403123, 2018.
  * [48] S. A. Henneberg, M. Drevlak, C. Nührenberg, C. D. Beidler, Y. Turkin, J. Loizu, and P. Helander.  Properties of a new quasi-axisymmetric configuration.  Nuclear Fusion, 59(2):026014, 2019.
  * [49] M. Barnes, F. I. Parra, and M. Landreman.  stella: An operator-split, implicit–explicit $\delta$f-gyrokinetic code for general magnetic field configurations.  Journal of Computational Physics, 391:365, 2019.
  * [50] J. H.E. Proll, P. Xanthopoulos, and P. Helander.  Collisionless microinstabilities in stellarators. II. Numerical simulations.  Physics of Plasmas, 20(12):122506, 2013.
  * [51] G. G. Plunk, P. Helander, P. Xanthopoulos, and J. W. Connor.  Collisionless microinstabilities in stellarators. III. the ion-temperature-gradient mode.  Physics of Plasmas, 21(3):032112, 2014.
  * [52] A. Zocco, G. G. Plunk, and P. Xanthopoulos.  Geometric stabilization of the electrostatic ion-temperature-gradient driven instability. II. Non-axisymmetric systems.  Physics of Plasmas, 27(2):022507, 2020.
  * [53] J. H.E. Proll, H. E. Mynick, P. Xanthopoulos, S. A. Lazerson, and B. J. Faber.  TEM turbulence optimisation in stellarators.  Plasma Physics and Controlled Fusion, 58(1):014006, 2015.
  * [54] O. Beeke, M. Barnes, M. Romanelli, M. Nakata, and M. Yoshida.  Impact of plasma shaping on tokamak microstability.  arXiv:2012.02669, 2020.
  * [55] M. Landreman and W. Sengupta.  Constructing stellarators with quasisymmetry to high order.  Journal of Plasma Physics, 85(6):815850601, 2019.
  * [56] M. Landreman, W. Sengupta, and G. G. Plunk.  Direct construction of optimized stellarator shapes. Part 2. Numerical quasisymmetric solutions.  Journal of Plasma Physics, 85(1):905850103, 2019.
  * [57] F. Romanelli.  Ion temperature-gradient-driven modes and anomalous ion transport in tokamaks.  Physics of Fluids B, 1(5):1018–1025, 1989.
  * [58] H. Biglari, P. H. Diamond, and M. N. Rosenbluth.  Toroidal ion-pressure-gradient-driven drift instabilities and transport revisited.  Physics of Fluids B, 1(1):109, 1989.
  * [59] G. T. Roberg-Clark, G. G. Plunk, and P. Xanthopoulos.  Calculating the linear critical gradient for the ion-temperature-gradient mode in magnetically confined plasmas.  arXiv:2010.13441, 2020.
  * [60] M. N. Rosenbluth.  Low-frequency limit of interchange instability.  Physics of Fluids, 11(4):869, 1968.
  * [61] J. A. Alcuson, P. Xanthopoulos, G. G. Plunk, P. Helander, F. Wilms, Y. Turkin, A. Von Stechow, and O. Grulke.  Suppression of electrostatic micro-instabilities in maximum-J stellarators.  Plasma Physics and Controlled Fusion, 62(3):035005, 2020.