# Fast ion transport in quasisymmetric equilibria in the presence of a
resonant Alfvénic perturbation

Elizabeth J. Paul 1<EMAIL_ADDRESS>Harry E. Mynick2  Amitava
Bhattacharjee2  1Columbia University, New York, NY 10027, USA
2Princeton Plasma Physics Laboratory, Princeton, NJ 08540, USA

###### Abstract

Significant progress has been made in designing magnetic fields that provide
excellent confinement of the guiding enter trajectories of alpha particles
using quasisymmetry (QS). Given the reduction in this transport channel, we
assess the impact of resonant Alfvén eigenmodes (AEs) on the guiding center
motion. The AE amplitudes are chosen to be consistent with experimental
measurements and large-scale simulations. We evaluate the drift resonance
condition, phase-space island width, and island overlap criterion for
quasisymmetric configurations. Kinetic Poincaré plots elucidate features of
the transport, including stiff transport above a critical perturbation
amplitude. Our analysis highlights key departures from the AE-driven transport
in tokamaks, such as the avoidance of phase-space island overlap in
quasihelical configurations and the enhanced transport due to wide phase-space
islands in low magnetic shear configurations. In configurations that are
closer to QS, with QS deviations $\delta B/B_{0}\lesssim 10^{-3}$, the
transport is primarily driven by the AE, while configurations that are further
from QS, $\delta B/B_{0}\sim 10^{-2}$, experience significant transport due to
the QS-breaking fields in addition to the AE.

##  1 Introduction

Energetic particles have historically been challenging to confine in
stellarator configurations due to the possibility of unconfined orbits and
resonances exposed at low collisionality. These difficulties must be overcome
to develop a stellarator reactor concept, as excessive alpha losses before
thermalization can impact power balance and impart damage to plasma-facing
components. In recent years, significant progress has been made in designing
stellarator magnetic fields that can confine the orbits of fusion-born alpha
particles without perturbations (Bader et al., 2019; Landreman & Paul, 2022;
Landreman et al., 2022). However, reducing the guiding center orbit loss
mechanisms may make mode-particle interactions relatively significant.

The interaction of Alfvén eigenmodes (AEs) with energetic particles has been
shown to drive substantial flattening of the fast-ion profile in tokamak
experiments (Heidbrink et al., 2008). Alvénic activity has also been observed
on several stellarator configurations, including HSX (Deng et al., 2009), CHS
(Takechi et al., 2002), LHD (Toi et al., 2011), W7-AS (Weller et al., 1994),
TJ-II (Melnikov et al., 2014), W7-X (Rahbarnia et al., 2020), and Heliotron-J
(Yamamoto et al., 2007). Alfvénic instabilities have been thought to be
potentially benign in a stellarator reactor due to their ability to operate at
high density (Helander et al., 2012). However, fast-ion-driven modes may still
be destabilized at high density: LHD modeling indicates that Alfvénic activity
remains present even at fast ion beta of $\approx 0.05\%$ (Varela et al.,
2017). For comparison, using the profiles from the ARIES-CS stellarator
reactor study with density $\approx 5\times 10^{20}$ m-3 (Ku et al., 2008),
the fast ion beta is $\approx 0.2\%$. It remains to be seen to what extent
Alfvénic activity can be controlled in a stellarator reactor by manipulating
the density profile or optimizing the magnetic field.

Significant recent work has focused on energetic particle physics in
quasisymmetric configurations based on properties of guiding center orbits
(Bader et al., 2021; LeViness et al., 2022; Paul et al., 2022), but relatively
little has been studied with respect to AE-driven transport. AE stability of
the quasiaxisymmetric CFQS configuration has been evaluated using the linear
gyrofluid FAR3D code (Varela et al., 2021). However, a systematic study of the
AE-driven transport in quasisymmetric configurations has not yet been
performed. For equilibria sufficiently close to quasisymmetry (QS), phenomena
previously observed on tokamaks are anticipated. Monte Carlo simulations
indicate that a resonant perturbation induces some rapid convective transport
due to phase-space islands near the boundaries. Island overlap occurs for
larger perturbation amplitudes, $\delta B^{r}/B_{0}\sim 10^{-3}$ where $\delta
B^{r}$ is the radial perturbed magnetic field and $B_{0}$ is the equilibrium
field strength, causing diffusive losses (Sigmar et al., 1992; Hsu & Sigmar,
1992). Because of the coupling of a single AE to the poloidal variation of the
magnetic drifts, sideband resonances arise and can lead to island overlap even
in the presence of a single AE (Mynick, 1993b, a; Hsu & Sigmar, 1992). In
perfect symmetry with the addition of a single perturbation harmonic, kinetic
Poincaré plots (White, 2011, 2012) can be employed to observe the formation of
phase-space islands, island overlap, and chaos.

We aim to address unresolved questions in this area, including how the
quasisymmetry helicity and deviations from quasisymmetry impact transport.
Recent simulations (White et al., 2022; White & Duarte, 2023) of resonant AEs
in W7-X and the precise QH equilibrium (Landreman & Paul, 2022) have indicated
that even a small-amplitude Alfvénic perturbation, $\delta B/B\sim 10^{-6}$,
can lead to global flattening of the distribution function. This extreme
sensitivity to perturbations is postulated to arise because the low magnetic
shear of the equilibrium implies low transit frequency shear for passing
particles. Here, we study the impact of magnetic shear in more detail by
evaluating the island width and drift-harmonic overlap conditions for
quasisymmetric configurations. As discussed later, our conclusions do not
quantitatively agree with these recent studies, which predict substantive
diffusive losses for low AE perturbation amplitudes for equilibria very close
to QS.

The impact of AEs on the fast-ion transport in a stellarator reactor is
challenging to compute in practice, requiring knowledge of the saturated mode
amplitude. The nonlinear saturation can be obtained from high-fidelity
modeling (Todo et al., 2017; Fehér, 2014; Spong et al., 2017), which depends
on the details of the thermal and fast-ion profiles. Rather than attempt such
calculations here, we consider the potential impact of a single Alfvénic
perturbation on the phase-space integrability and resulting transport in
reactor-scale equilibria designed to be close to quasisymmetry. We employ a
“worst-case scenario” approach, in which an Alfvénic perturbation is chosen to
strongly resonate with rational periodic passing orbits in the core for
several quasisymmetric configurations. We consider several mode amplitudes
consistent with experimental measurements and high-fidelity modeling. While
physically, such a perturbation should correspond to an AE of the background
plasma, we instead let the perturbation be a radially global mode with
prescribed mode numbers, frequency, and perturbation amplitude. We assess the
impact of the perturbation mode number, amplitude, and residual quasisymmetry
error on the resulting transport. Phase-space Poincaré plots are used to guide
the analysis and assess the transition to chaos. While Poincaré plots have
been used to study the structure of phase space for energetic passing
particles in the absence of time-dependent perturbations for some stellarators
(White et al., 2022), this technique has not yet been applied to assess the
impact of AEs on the transport. We evaluate the extent to which this technique
provides insight into the transport in configurations with varying deviations
from quasisymmetry.

In Section 2, we outline the guiding center equations of motion with Alfvénic
perturbations. In Section 3, we describe the theory for resonance, island
width, and island overlap in QS configurations. Resonance analysis is
performed for several equilibria designed to be close to QS in Section 4, and
resonant perturbations are identified. Kinetic Poincaré plots are employed in
Section 5 to assess the formation of phase-space islands and chaos in the
presence of the resonant perturbations. A Monte Carlo guiding center transport
analysis is performed in Section 6. We conclude in Section 7.

##  2 Guiding-center motion in the presence of an Alfvénic perturbation

The guiding-center motion, $\bm{R}(t)$, is described by the Lagrangian
(Littlejohn, 1983),

$\displaystyle
L(\bm{R},\dot{\bm{R}},v_{\|})=q\left(\bm{A}+\frac{Mv_{\|}}{qB}\bm{B}\right)\cdot\dot{\bm{R}}-\frac{Mv_{\|}^{2}}{2}-\mu
B-q\Phi,$ (1)

where $q$ is the charge, $M$ is the mass, $v_{\|}$ is the velocity in the
direction of the magnetic field, and $\mu$ is the magnetic moment. We assume
the magnetic field is comprised of an equilibrium field, $\bm{B}_{0}$, and a
shear Alfvénic perturbation (White et al., 1983),

$\displaystyle\bm{B}=\bm{B}_{0}+\nabla\times\left(\alpha\bm{B}_{0}\right),$
(2)

while the scalar potential vanishes in the equilibrium, $\Phi=\delta\Phi$.
With the reduced MHD assumption (Kruger et al., 1998)—$k_{\|}/k_{\perp}\ll 1$,
where $k_{\|}$ is the characteristic parallel wave number and $k_{\perp}$ is
the characteristic perpendicular wave number associated with perturbed
quantities—the form of the perturbed field (2) implies that the linear
perturbation to the field strength vanishes, $\delta B=0$. Under the ideal MHD
assumption, the corresponding scalar potential must satisfy
$\bm{B}_{0}\cdot\delta\bm{E}=0$, which implies

$\displaystyle\nabla_{\|}\delta\Phi=-B_{0}\frac{\partial\alpha}{\partial t}.$
(3)

The perturbed fields could be computed from a global Alfvén eigenmode solver
such as AE3D (Spong et al., 2010). For the fundamental studies here, we
instead assume a single harmonic perturbation of the form,

$\displaystyle\delta\Phi=\hat{\Phi}(\psi)\sin(\omega t+m\theta-n\zeta),$ (4)

where $(\psi,\theta,\zeta)$ are Boozer coordinates such that the unperturbed
magnetic field is expressed as,

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle\bm{B}_{0}=\nabla\psi\times\nabla\theta-\iota(\psi)\nabla\psi\times\nabla\zeta,\\\
\bm{B}_{0}=G(\psi)\nabla\zeta+I(\psi)\nabla\theta+K(\psi,\theta,\zeta)\nabla\psi,\end{array}\right.$
(7)

$m$ is the poloidal mode number, $n$ is the toroidal mode number, and $\omega$
is the frequency. We define the phase variable $\eta=\omega t+m\theta-n\zeta$.
The corresponding expression for $\alpha$ is then obtained from (3), which we
write schematically as,

$\displaystyle\alpha=\hat{\alpha}(\psi)\sin(\eta).$ (8)

We impose $\hat{\Phi}$ and compute $\hat{\alpha}$ so that a magnetic
differential equation need not be inverted, which can give rise to singular
solutions on rational surfaces.

Since $Mv_{\|}/(qB_{0})$ is small in comparison to the characteristic length
scale of the equilibrium and $\delta B=0$, the perturbed Lagrangian reads
(Littlejohn, 1985),

$\displaystyle
L(\bm{R},\dot{\bm{R}},v_{\|})=q\left(\bm{A}_{0}+\alpha\bm{B}_{0}+\frac{Mv_{\|}}{qB_{0}}\bm{B}_{0}\right)\cdot\dot{\bm{R}}-\frac{Mv_{\|}^{2}}{2}-\mu
B_{0}-q\delta\Phi.$ (9)

The resulting equations of motion are integrated in Boozer coordinates using
the SIMSOPT stellarator optimization and modeling package (Landreman et al.,
2021).

To compare Alfvénic perturbations across configurations with different mode
numbers, we will express the strength of the perturbation in terms of its
normalized radial magnetic field,

$\displaystyle\frac{\delta\bm{B}\cdot\nabla\psi}{B_{0}|\nabla\psi|}\approx\frac{\nabla\alpha\times\bm{B}_{0}\cdot\nabla\psi}{rB_{0}^{2}}=\frac{\hat{\alpha}(mG-
nI)}{r(G+\iota I)}\cos(\eta).$ (10)

Here we have assumed that the gradient scale length of the perturbed field is
larger than that of the equilibrium. Since $G\gg I$ for stellarator
equilibria, we define the parameter
$\delta\hat{B}^{\psi}=m\hat{\alpha}/r\sim\delta\bm{B}\cdot\nabla\psi/(B_{0}|\nabla\psi|)$
to compare equilibria with respect to the strength of the perturbed radial
field.

##  3 Resonance theory

We assume a quasisymmetric equilibrium for which the unperturbed field
strength can be expressed as $B_{0}(\psi,\chi)$ where $\chi=\theta-N\zeta$ is
the symmetry angle and $N$ is an integer representing the symmetry helicity.
($N=0$ for quasiaxisymmetry and $N\neq 0$ for quasihelical symmetry.) Each
equilibrium we consider in Section 4 is sufficiently close to quasisymmetry
that such an assumption provides valuable insight into the resulting dynamics.

The radial drift over the unperturbed trajectories is analyzed in the presence
of an Alfvénic perturbation in Appendix A, generalizing the theory of Mynick
(1993b) to quasisymmetric configurations and moderate frequency perturbations.
At moderate frequencies, the electrostatic potential enters the guiding center
equations of motion. Because $\bm{B}_{0}\cdot\delta\bm{E}=0$, particles
continue to follow perturbed field lines to lowest order in the gyroradius but
experience an additional $\bm{E}\times\bm{B}$ drift. Under the assumption of
quasisymmetry and stellarator symmetry, the unperturbed equations of motion
can be written schematically as,

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle\dot{\chi}=\omega_{\chi}+\sum_{j\neq
0}\chi_{j}\cos(j\chi),\\\ \displaystyle\dot{\zeta}=\omega_{\zeta}+\sum_{j\neq
0}\zeta_{j}\cos(j\chi),\end{array}\right.$ (13)

where $\omega_{\chi}=\langle\dot{\chi}\rangle$ and
$\omega_{\zeta}=\langle\dot{\zeta}\rangle$ are the averaged drifts in the
$\chi$ and $\zeta$ directions. The overdot represents a time derivative, and
the average is performed over many toroidal transits, $\langle
A\rangle=\int_{0}^{T}dt\,A/T$. The summations represent the periodic
contributions from the drifts.

For an Alfvénic perturbation of the form (4) to provide a net radial drift, a
resonance condition must be satisfied,

$\displaystyle\Omega_{l}=(m+l)\omega_{\chi}-(n-Nm)\omega_{\zeta}+\omega=0.$
(14)

As discussed in Appendix A, the integer $l$ arises due to coupling through the
$\chi$-dependence of the magnetic drifts. If the drift dynamics (13) is
dominated by a particular cosine harmonic with integer $j^{\prime}$, then $l$
is assumed to be an integer multiple of $j^{\prime}$. To simplify the
analysis, the $j^{\prime}=1$ harmonic of the field strength is assumed to be
dominant, which holds near the magnetic axis (Garren & Boozer, 1991) and for
the equilibria of interest. More general expressions are provided in Appendix
A. Under the assumption that the characteristic frequencies are approximately
flux functions, the full island width associated with a given resonance is
given by,

$\displaystyle
w^{\psi}_{l}=2\sqrt{\left|\frac{\psi_{l}}{\Omega_{l}^{\prime}(\psi)}\right|}\approx
2\sqrt{\left|\frac{\psi_{l\neq
0}}{(m+l)\omega_{\theta}^{\prime}(\psi)}\right|},$ (15)

where $\psi_{l}$ is a cosine harmonic of the radial perturbed drift and we
have made the approximation
$h^{\prime}(\psi)\approx\omega_{\theta}^{\prime}(\psi)/\omega_{\zeta}$. The
perturbed radial drift, $\delta\dot{\psi}$, can be evaluated along the
unperturbed trajectory and expressed as a cosine series in
$\cos(\Omega_{l}t+\eta^{0})$ and $\cos(\Omega_{l}t+\eta^{0}\mp\chi^{0})$ with
coefficients given by $\psi_{l}^{0}$ and $\psi_{l}^{\pm}$; see (32). Here
$\eta^{0}=m\theta^{0}-n\zeta^{0}$ is the initial phase. The scaling of these
coefficients with the magnitude of the magnetic drifts, mode numbers, and
perturbed radial field is summarized as,

$\displaystyle\left\\{{\begin{array}[]{l}\displaystyle\psi_{0}^{0}\sim(\iota-\omega_{\theta}/\omega_{\zeta})J_{0}(\eta_{1})\delta\hat{B}^{\psi},\\\
\displaystyle\psi_{l}^{0}\sim J_{l}(\eta_{1})\delta\hat{B}^{\psi},\\\
\displaystyle\psi_{l}^{\pm}\sim J_{l\pm
1}(\eta_{1})\zeta_{1}\delta\hat{B}^{\psi},\end{array}}\right.$ (19)

where $J_{l}$ are the Bessel functions of the first kind,
$\eta_{1}=\left(m\chi_{1}-(n-Nm)\zeta_{1}\right)/\omega_{\chi}$ with
$\chi_{1}$ and $\zeta_{1}$ defined through (13), and
$\omega_{\theta}=\langle\dot{\theta}\rangle$. The full expressions for general
$j^{\prime}$ are provided in Appendix A.

When considering the dependence of the island width on the magnetic drifts in
the small argument limit of the Bessel functions, the most significant radial
transport will arise from $\psi_{0}^{0}$, $\psi_{\pm 1}^{0}$, or $\psi_{\mp
1}^{\pm}$. In the limit of large mode numbers, $\psi_{\pm 1}^{0}$ will
dominate due to its dependence on $\eta_{1}$. Considering the small argument
limit of the Bessel function, the island width scales with
$\sqrt{\eta_{1}^{|l|}/(m+1)}$ for fixed $\delta\hat{B}^{\psi}$, thus
increasing the island width for quasihelical configurations due to the
dependence of $\eta_{1}$ on the helicity $N$. Given the scaling of $\eta_{1}$
with the mode numbers, the island width will roughly scale independently of
$m$ for $|l|=1$ but will scale with $m^{(|l|-1)/2}$ for $|l|>1$. Finally, the
island width decreases strongly with increasing $|l|$.

Defining the passing orbital helicity as $h=\omega_{\theta}/\omega_{\zeta}$,
resonances occur where:

$\displaystyle h=\frac{n-Nm-\omega/\omega_{\zeta}}{m+l}+N.$ (20)

For a given primary resonance $l$ at $h=h_{0}$, additional sideband resonances
may be excited for other drift harmonics $l^{\prime}$ corresponding to
neighboring periodic orbits. The spacing between the $l$ and $l+1$ resonances
is given by,

$\displaystyle\left(\Delta\psi\right)_{l}=\frac{1}{h^{\prime}(\psi)}\frac{h_{0}-N}{m+l+1}.$
(21)

The ratio between the island width and the resonance spacing

$\displaystyle\frac{w_{l}^{\psi}}{\left(\Delta\psi\right)_{l}}\approx\frac{m+l+1}{(h_{0}-N)\omega_{\zeta}}\sqrt{\left|\frac{\psi_{l}\omega_{\theta}^{\prime}(\psi)}{m+l}\right|}$
(22)

provides a conservative estimate for the island overlap criterion,
$w_{l}^{\psi}/\left(\Delta\psi\right)_{l}\gtrsim 1$. Given the scaling of
$\psi_{l}$ (19), the potential for island overlap increases with $m$ and
decreases with $N$ for fixed $\delta\hat{B}^{\psi}$. In this way, quasihelical
configurations are advantageous for preventing the transition to phase-space
chaos. Finally, while shear in the transit frequency,
$\omega_{\theta}^{\prime}(\psi)$, reduces individual island widths, it
promotes island overlap if multiple resonances are present.

##  4 Resonance analysis

We consider four equilibria optimized to be close to quasisymmetry:
$\beta=2.5\%$ quasihelical (QH) and quasiaxisymmetric (QA) equilibria with
self-consistent bootstrap current (Landreman et al., 2022), a vacuum
equilibrium with precise levels of quasiaxisymmetry (Landreman & Paul, 2022),
and the quasiaxisymmetric NCSX li383 equilibrium (Koniges et al., 2003; Mynick
et al., 2002). Each equilibrium considered is scaled to the minor radius (1.70
m) and the field strength (5.86 T) of ARIES-CS (Najmabadi et al., 2008). The
rotational transform profiles and quasisymmetry error,

$f_{QS}(s)=\frac{\sqrt{\sum_{Mn\neq
Nm}\left(B_{m,n}^{c}(s)\right)^{2}}}{\sqrt{\left(B_{0,0}^{c}(s)\right)^{2}}},$
(23)

are shown in Figure 1, where the unperturbed field strength in Boozer
coordinates is
$B_{0}(s,\theta,\zeta)=\sum_{m,n}B_{m,n}^{c}(s)\cos(m\theta-n\zeta)$.

(a)

(b)

Figure 1: (a) Rotational transform and (b) quasisymmetry error (23) profiles
for the four equilibria under consideration.

We first identify the low-order periodic passing orbits present in the
equilibrium to determine the impact of a potential resonant perturbation on
the guiding center losses. Fusion-born alpha particles are initialized with
equally-spaced pitch angle and radius values and followed for 500 toroidal
transits. The net change in the poloidal angle, $\Delta\theta$, and transit
time, $\Delta t$, are computed for each toroidal transit. In the integrable
case, the average of the characteristic frequencies over many transits,
denoted $\omega_{\zeta}=2\pi/\langle\Delta t\rangle$, and
$\omega_{\theta}=\langle\Delta\theta\rangle/\langle\Delta t\rangle$, converges
quickly with respect to the number of toroidal transits (Das et al., 2016). In
this way, we obtain $\omega_{\zeta}$ and $\omega_{\theta}$ in the two-
dimensional space $(s,\mu)$, where $s=\psi/\psi_{0}$ is the normalized
toroidal flux and $\psi_{0}$ is the value of the toroidal flux on the
boundary. All nonsymmetric modes are artificially suppressed for this
frequency analysis so that all trajectories lie on KAM surfaces.

For the subsequent analysis, resonant perturbations are intentionally imposed,
satisfying the condition (14) for a resonant surface near mid-radius. As many
potential Alfvénic perturbations exist that resonate with a given periodic
orbit, a comparison is made between $m=1$, 15, and 30 perturbations. The
perturbation with the highest mode number is chosen such that $m\approx
a/\rho_{EP}$ for the ARIES-CS reactor parameters, where $a$ is the effective
minor radius and $\rho_{EP}$ is the energetic particle orbit width. With this
choice, the radial width of the AE eigenstructure is predicted to be
comparable to the EP orbit width, and the growth rate is maximized (Gorelenkov
et al., 2014). The smaller values of $m$ are more typical for current
experiments such as LHD (Varela et al., 2017). We choose perturbation
amplitudes in the range $\delta\hat{B}^{\psi}\sim 10^{-4}-10^{-3}$, consistent
with LHD modeling (Nishimura et al., 2013) and experimental measurements from
TFTR and (Nazikian et al., 1997) and NSTX (Crocker et al., 2013). In analysis
of the impact of TAEs on guiding center confinement in tokamaks, alpha orbit
stochasticity is typically present for $\delta\hat{B}^{\psi}\sim 10^{-3}$,
while for $\delta\hat{B}^{\psi}<10^{-4}$ the losses are insignificant (Sigmar
et al., 1992). Although the eigenfunction typically depends on the mode
number, we choose a radially uniform mode structure for a conservative
analysis.

Given a $\omega_{\theta}/\omega_{\zeta}=p/q$ periodic orbit and poloidal mode
number $m$, the resonant wave frequency satisfying (14) can be expressed as
$\omega/\omega_{\zeta}=p^{\prime}/q$ for some integer $p^{\prime}$. The
toroidal mode number $n$ and the frequency parameter $p^{\prime}$ are chosen
to satisfy the $l=1$ resonance condition (14), given that the $l=\pm 1$
resonance is predicted to give rise to the most substantial radial transport
as described in Section 3. The value of $n$ is chosen to yield small
magnitudes of $|p^{\prime}|$, corresponding with lower-frequency
perturbations. The mode parameters are summarized in Table 1.

$m$ | $n$ | $\omega/\omega_{\zeta}$ | $\omega_{\theta}/\omega_{\zeta}$ |  $\omega$ [kHz]
---|---|---|---|---
1 | 1 | 1/9 | 4/9 | 136
15 | 7 | -1/9 | 4/9 | -136
30 | 14 | 2/9 | 4/9 | 272
(a) $\beta=2.5\%$ QA ($N=0$)

$m$ | $n$ | $\omega/\omega_{\zeta}$ | $\omega_{\theta}/\omega_{\zeta}$ |  $\omega$ [kHz]
---|---|---|---|---
1 | -2 | 2/14 | -15/14 | 133
15 | 13 | 2/14 | -15/14 | 133
30 | 29 | 3/14 | -15/14 | 200
(b) $\beta=2.5\%$ QH ($N=4$)

$m$ | $n$ | $\omega/\omega_{\zeta}$ | $\omega_{\theta}/\omega_{\zeta}$ |  $\omega$ [kHz]
---|---|---|---|---
1 | 1 | 2/12 | 5/12 | 199
15 | 7 | 4/12 | 5/12 | 399
30 | 13 | 1/12 | 5/12 | 997
(c) Vacuum QA ($N=0$)

$m$ | $n$ | $\omega/\omega_{\zeta}$ | $\omega_{\theta}/\omega_{\zeta}$ |  $\omega$ [kHz]
---|---|---|---|---
1 | 1 | -1/5 | 3/5 | -338
15 | 10 | 2/5 | 3/5 | 676
30 | 19 | 2/5 | 3/5 | 676
(d) NCSX ($N=0$)

Table 1: Alfvénic perturbation mode parameters chosen to satisfy the resonance
condition (14).

For each configuration, the profile of the effective orbit helicity $h$ (20)
is shown in Figure 2 for co- and counter-passing orbits (blue and orange,
respectively). While the helicity profile does not change substantially
between co- and counter-passing orbits, the resonance condition (20) is
modified since the sign and magnitude of $\omega_{\zeta}$ differs. Although
the choice of pitch angle will impact the resonance condition, for the
following calculations, we focus on co-passing orbits, $v_{\|}/v_{0}=+1$. The
horizontal black lines in the figures indicate the chosen resonant surface for
the co-passing orbits. For the $\beta=2.5\%$ QA and NCSX equilibria, the
nearby resonant surfaces corresponding to other values of $l$ are indicated by
the colored horizontal lines. As $m$ increases, the resonance spacing
decreases according to (21). For the vacuum QA equilibrium, the magnetic shear
is very low, see Figure 1, and no sideband resonances exist. Although the
shear is moderate for the $\beta=2.5\%$ QH equilibrium, sideband resonances
are not present due to the factor of $N$ in the resonance spacing expression.

(a) $\beta=2.5\%$ QA

(b) $\beta=2.5\%$ QH

(c) Vacuum QA

(d) NCSX

Figure 2: The characteristic orbit helicity,
$h=\omega_{\theta}/\omega_{\zeta}$, is computed for co-passing
($v_{\|}/v_{0}=+1$, blue) and counter-passing orbits ($v_{\|}/v_{0}=-1$,
orange). A low-order periodic orbit is selected near the mid-radius, denoted
by the horizontal dashed black line. Alfvénic perturbations with several mode
numbers $m$ are chosen to resonate with this orbit periodicity. Sideband
resonances are excited in the $\beta=2.5\%$ QA and NCSX equilibria due to the
angular dependence of the drifts, denoted by the colored horizontal lines.
Because of the increased distance between resonances, no sidebands are excited
for in the vacuum QA or $\beta=2.5\%$ QH equilibria for the mode numbers
chosen.

##  5 Kinetic Poincaré plots

Given a magnetic field with exact quasisymmetry, a kinetic Poincaré plot can
be constructed analogously to the case of axisymmetric magnetic fields (Sigmar
et al., 1992). Given a quasisymmetric field $B_{0}(\psi,\chi=\theta-N\zeta)$
in the absence of a perturbation, the canonical angular momentum,

$\displaystyle
P_{\zeta}=\left(G+NI\right)\left(\frac{mv_{\|}}{B_{0}}+q\alpha\right)+q\left(N\psi-\psi_{P}\right),$
(24)

and energy,

$\displaystyle E=\frac{mv_{\|}^{2}}{2}+\mu B_{0}+q\delta\Phi,$ (25)

are conserved. When an Alfvénic perturbation is applied with single mode
numbers $n$ and $m$, neither $P_{\zeta}$ nor $E$ are conserved, but a
conserved quantity is obtained by moving with the wave frame,

$\displaystyle\overline{E}_{n}=\left(n-Nm\right)E-\omega P_{\zeta}.$ (26)

Given the velocity-space parameters—$\mu$, $\overline{E}_{n}$, and
$\text{sign}(v_{\|})$—and specification of the position in space and time—$t$,
$s$ , $\theta$, and $\zeta$—equation (26) provides a nonlinear equation for
$v_{\|}$ (Hsu & Sigmar, 1992).

Given the form for the phase factor of the perturbation (4), the resulting
equations of motion depend on $(s,\chi,\zeta-\overline{\omega}_{n}t,v_{\|})$,
where $\overline{\omega}_{n}=\omega/(n-Nm)$. The resulting motion is generally
4D. If purely passing particles are considered, then
$\text{sign}\left(v_{\|}\right)$ is fixed, $v_{\|}$ can be computed from the
other coordinates as described above, and the resulting motion becomes 3D. For
a fixed value of $\overline{E}_{n}$ and $\mu$, a kinetic Poincaré map
$M(s,\chi)\rightarrow(s^{\prime},\chi^{\prime})$ is constructed by moving with
the wave frame to eliminate the time dependence. Guiding center trajectories
are followed until they intersect a plane of constant
$\zeta-\overline{\omega}_{n}t$. Therefore, the resulting Poincaré section is
2D.

When quasisymmetry is broken in the equilibrium, $\overline{E}_{n}$ is no
longer precisely conserved. Nonetheless, if the quasisymmetry errors are
sufficiently small, the kinetic Poincaré analysis provides insight into the
resulting transport. The non-quasisymmetric modes are artificially suppressed
when constructing the kinetic Poincaré plots to aid the analysis. The impact
of quasisymmetry-breaking modes will be assessed in the following Section.

(a) $m=1$, $\beta=2.5\%$ QA

(b) $m=15$, $\beta=2.5\%$ QA

(c) $m=30$, $\beta=2.5\%$ QA

(d) $m=1$, $\beta=2.5\%$ QH

(e) $m=15$, $\beta=2.5\%$ QH

(f) $m=30$, $\beta=2.5\%$ QH

(g) $m=1$, vacuum QA

(h) $m=15$, vacuum QA

(i) $m=30$, vacuum QA

(j) $m=1$, NCSX

(k) $m=15$, NCSX

(l) $m=30$, NCSX

Figure 3: Kinetic Poincaré plots are constructed using the mode parameters in
Table 1 with the perturbation amplitude $\delta\hat{B}^{\psi}=10^{-3}$. All
QS-breaking harmonics of the equilibrium field are artificially suppressed for
this analysis.

In Figure 3, kinetic Poincaré plots are displayed for particles with
$\text{sign}(v_{\|})=+1$, $\mu=0$, $E=3.52$ MeV and Alfvénic perturbations
with parameters given in Table 1. The amplitude $\hat{\Phi}$ is chosen such
that $\delta\hat{B}^{\psi}=10^{-3}$ on the $s=1$ surface. Periodic orbits that
satisfy the resonance condition $\Omega_{l}=0$ appear as $n/(m+l)$ periodic
orbits in the kinetic Poincaré plots.

For the NCSX and $\beta=2.5\%$ QA configurations, a clear $l=1$ island chain
is apparent in the presence of the $m=1$ perturbation. For the vacuum QA
equilibrium, the perturbation shifts the orbit helicity slightly, moving the
resonance outside the equilibrium due to the low magnetic shear. The resonance
reappears by increasing the perturbation frequency by 8%, see Figure 4. In the
$\beta=2.5\%$ QA equilibria, in addition to the primary $l=1$ resonance, the
$l=2$ resonance is apparent near $s=0.1$. However, the island width is small
due to the reduced magnitude of the $J_{2}(\eta_{1})$ coupling parameter. As
indicated by the resonance plots, Figure 2, none of the other equilibria
contain the $m=1$, $l=2$ resonance.

(a)

(b)

Figure 4: (a) In the presence of the resonant perturbation of amplitude
$\delta\hat{B}^{\psi}=10^{-4}$, the characteristic frequencies
$\omega_{\theta}$ and $\omega_{\zeta}$ shift, causing the resonance defined by
$\Omega_{l}=0$ to move outside the equilibrium due to the low shear, see
Figure 3(g). By increasing the mode frequency by 8%, the resonance reenters.
(b) The kinetic Poincaré plot of amplitude $\delta\hat{B}^{\psi}=10^{-4}$ with
the shifted frequency ($\omega=2.15$ kHz) reveals the corresponding $l=1$
island on the shifted resonant surface.

As the mode number increases from $m=15$ to $m=30$, the $l=1$ island width
stays roughly the same for the $\beta=2.5\%$ QH configuration, as indicated by
the scaling of the island width formula (15)-(19) for fixed
$\delta\hat{B}^{\psi}$. For the vacuum QA configuration, the $l=1$ resonance
reappears for the $m=15$ and $m=30$ perturbations, its width being especially
wide due to the low magnetic shear of the configuration. For both the vacuum
QA and $\beta=2.5\%$ QH configurations with the $m=15$ and $m=30$
perturbations, the island chain is wide enough to lead to visible destruction
of nearby KAM surfaces and the formation of secondary island chains.

In the $\beta=2.5\%$ QA equilibrium, overlap between the $l=0$, 1, 2, and 3
resonances is observed with the $m=15$ perturbation. Substantial island
overlap is also observed in the NCSX equilibrium with the $m=15$ perturbation.
In the presence of the $m=30$ perturbation, strong island overlap is observed
in the $\beta=2.5\%$ QA and NCSX equilibria. However, the phase-space volume
over which island overlap occurs is narrowed. As seen in Figure 2, the
resonances become more closely spaced for the $m=30$ perturbation compared to
the $m=15$ perturbation. However, because the island width scales as
$\sqrt{\eta_{1}^{|l|}}$ with $\eta_{1}\ll 1$, island overlap does not occur
for the large $|l|$ resonant surfaces. This island-width scaling effectively
reduces the non-integrable volume as $m$ increases.

The impact of these phase-space features on the resulting transport will be
discussed in Section 6.

###  5.1 Impact of quasisymmetry deviations

We now investigate the impact of the finite quasisymmetry deviations,
quantified in Figure 1, on the kinetic Poincaré analysis. When QS deviations
are present, the motion of passing particles becomes 4D $(s,\chi,\zeta,t)$,
and the Poincaré section becomes 3D,
$M(s,\chi,\zeta)\rightarrow(s^{\prime},\chi^{\prime},\zeta^{\prime})$.
Nonetheless, we can still visualize the Poincaré map in the $(s,\chi)$ plane
to assess the structure of phase space.

In Figures 5-6, we show a selection of the kinetic Poincaré plots constructed
with the same parameters as those in Figure 3, but without the suppression of
the QS-breaking modes. For the case of the $\beta=2.5\%$ QH equilibrium, the
phase-space structure remains mostly unchanged, with the addition of small-
scale mixing. In the case of the NCSX equilibrium, the amplitude of the phase-
space mixing is amplified due to the enhanced symmetry breaking. The gross
phase-space structure remains mostly unchanged in the presence of the $m=1$
perturbation. In the presence of the $m=15$ perturbation, the effective
diffusion due to the non-integrability of the orbits destroys many of the
remaining KAM surfaces and island chains, leading to a wide region of phase-
space chaos. We conclude that most large-scale phase space structure is
preserved by adding QS deviations, especially for equilibria close to QS,
$f_{QS}\sim 10^{-3}$. This finding does not quantitatively agree with recent
results (White & Duarte, 2023), indicating substantial diffusive losses in an
equilibrium very close to QH ($f_{QS}\sim 10^{-4}$) with an Alfvénic
perturbation as small as $\delta\hat{B}^{\psi}\sim 10^{-6}$.

(a) $m=1$

(b) $m=15$

Figure 5: $\beta=2.5\%$ QH kinetic Poincaré plots with the same parameters as
Figure 3, but without the suppression of the quasisymmetry-breaking modes.

(a) $m=1$

(b) $m=15$

Figure 6: NCSX kinetic Poincaré plots with the same parameters as Figure 3,
but without the suppression of the quasisymmetry-breaking modes.

##  6 Monte Carlo analysis

To assess the impact of the phase-space structure on the resulting transport,
we perform Monte Carlo collisionless guiding center tracing simulations. We
initialize 5000 particles uniformly in pitch angle and volume for
$s\in[0.25,0.75]$. This initial spatial distribution was chosen to avoid
effects related to the coordinate singularity. Particles are followed for
$10^{-3}$ seconds or until they are considered lost when they cross through
$s=0$ or $s=1$.111Regularity of the equilibrium and perturbed fields near the
axis ensures that guiding center trajectories never cross through $s=0$. Due
to the assumption of constant $\delta\hat{\Phi}$, regularity on the axis is
not preserved. We, nonetheless, can still gauge the impact of the perturbation
on the transport by assessing the loss fraction and flattening of the fast-ion
profile. Figures 7-8 display results of tracing performed with the true
equilibria (solid curves, “equilibrium”) as well as with the quasisymmetry-
breaking modes artificially filtered out (dashed curves, “perfect QS”).
Calculations are carried out without Alfvénic perturbations ($\hat{\alpha}=0$)
and for the $m=1$, 15, and 30 perturbations indicated in Table 1.

In Figure 8, the effective transport is quantified by the distribution of the
radial displacement, $|\Delta s|$, between $t=0$ and the final recorded time
(after $10^{-3}$ or when a particle is considered lost). By initializing with
a uniform distribution over $s$ in a fixed interval, flattening of the
distribution function is not as apparent. On the other hand, the distribution
of $|\Delta s|$ allows us to identify the net displacement of trajectories due
to phase-space islands or chaos. The double maximum structure reflects the
confined trajectories, centered around $|\Delta s|=0$, and lost trajectories,
for which $|\Delta s|\geq 0.25$. The loss fraction as a function of time is
also shown in Figure 7. A summary of the results is presented in Table 2. When
comparing the equilibrium and perfect QS cases, we see slight transport
enhancement due to QS-breaking modes in the $\beta=2.5\%$ QA and QH and vacuum
QA equilibria. The enhancement of transport becomes much more pronounced in
the NCSX equilibrium, given its more substantial deviations from QS; see
Figure 1.

In the equilibria that do not exhibit island overlap in Figure 3, the vacuum
QA and $\beta=2.5\%$ QH equilibria, the total loss fraction increases
monotonically with $m$. This result matches the observations in Figure 3,
which indicate that although the island width remains roughly the same, the
destruction of nearby KAM surfaces increases with increasing mode number. On
the other hand, in equilibria that exhibit strong resonance overlap in Figure
3, NCSX and $\beta=2.5\%$ QA, the transport increases monotonically with $m$
until $m=15$, then decreases for $m=30$. This characteristic is explained in
terms of the kinetic Poincaré plots in Figure 3, which indicate that the
effective volume of non-integrability is decreased for $m=30$ compared to
$m=15$ due to the reduction in island width for large $l$ perturbations.

Overall, the total induced transport is the lowest for the $\beta=2.5\%$ QH
configuration, for which the total losses remain less than $10\%$ in the
presence of the Alfvénic perturbations. The losses are also less prompt for
this configuration, with most beginning around $10^{-4}$ seconds, rather than
around $10^{-5}$ seconds for the other configurations. This reduction in
transport is due to the increased resonance spacing of quasihelical
configurations and its increased magnetic shear compared to the vacuum QA
equilibria, see Figure 1.

(a) $\beta=2.5\%$ QA

(b) $\beta=2.5\%$ QH

(c) Vacuum QA

(d) NCSX

Figure 7: The loss fraction as a function of time is shown. Monte Carlo
tracing is performed in the presence of a $\delta\hat{B}^{\psi}=10^{-3}$
perturbation with parameters described in Table 1. Guiding center trajectories
of fusion-born alpha particles are followed for $10^{-3}$ seconds or until
they cross through $s=0$ or $s=1$. A comparison is made between the actual
equilibrium and the equilibrium for which the QS-breaking modes are
suppressed, “perfect QS.”

(a) $\beta=2.5\%$ QA

(b) $\beta=2.5\%$ QH

(c) Vacuum QA

(d) NCSX

Figure 8: Distribution of radial displacement, $|\Delta s|$, between $t=0$ and the final recorded time among Monte Carlo samples for the same calculation presented in Figure 7. Configuration | Perturbation |  $\text{mean}|\Delta s|$ (Equil./Perfect QS) | Loss frac. (Equil./Perfect QS)
---|---|---|---
$\beta=2.5\%$ QA | Unperturbed | 0.082/0.085 | 0.046/0.039
| $m=1$ | 0.092/0.088 | 0.067/0.054
| $m=15$ | 0.160/0.143 | 0.179/0.148
| $m=30$ | 0.114/0.095 | 0.105/0.063
$\beta=2.5\%$ QH | Unperturbed | 0.037/0.032 | 0.007/0
| $m=1$ | 0.040/0.033 | 0.009/0
| $m=15$ | 0.051/0.042 | 0.019/0.002
| $m=30$ | 0.061/0.046 | 0.035/0.008
NCSX | Unperturbed | 0.127/0.055 | 0.178/0.005
| $m=1$ | 0.144/0.064 | 0.213/0.012
| $m=15$ | 0.227/0.137 | 0.367/0.146
| $m=30$ | 0.170/0.081 | 0.251/0.044
QA vac | Unperturbed | 0.105/0.112 | 0.083/0.066
| $m=1$ | 0.128/0.132 | 0.125/0.131
| $m=15$ | 0.163/0.146 | 0.194/0.162
| $m=30$ | 0.177/0.167 | 0.203/0.170
Table 2: Summary of transport properties for Monte Carlo calculations in the
presence of Alfvénic perturbations with mode parameters described in Table 1
and amplitude $\delta\hat{B}^{\psi}=10^{-3}$. The mean radial displacement
along a trajectory, $\text{mean}|\Delta s|$, and total loss fraction after
$10^{-3}$ s are compared between the actual equilibrium, “equil”, and the
equilibrium for which the QS-breaking modes are suppressed, “perfect QS.”

###  6.1 Transition to phase-space chaos

We now study the transition to phase-space chaos and global transport with
increasing perturbation amplitude in the $\beta=2.5\%$ QA equilibria. We focus
on the $m=30$ perturbation, given that this is the expected mode number at
reactor scales. Kinetic Poincaré plots for co-passing particles resonant with
the Alfvénic perturbation are shown in Figure 9. The $l=-1$, 0, 1, 2, and 3
islands are visible at the smallest perturbation amplitude, corresponding to
$\delta\hat{B}^{\psi}=3.14\times 10^{-4}$. As the amplitude is increased to
$\delta\hat{B}^{\psi}=5.58\times 10^{-4}$, a band of island overlap appears
near the resonant surface. The region of destroyed KAM surfaces increases with
increasing perturbation amplitude, with only a few remnant islands present at
$\delta\hat{B}^{\psi}=1.76\times 10^{-3}$.

(a) $\delta\hat{B}^{\psi}=3.14\times 10^{-4}$

(b) $\delta\hat{B}^{\psi}=5.58\times 10^{-4}$

(c) $\delta\hat{B}^{\psi}=9.92\times 10^{-4}$

(d) $\delta\hat{B}^{\psi}=1.76\times 10^{-3}$

Figure 9: Kinetic Poincaré plots for co-passing orbits in the $\beta=2.5\%$ QA
in the presence of an $m=30$ perturbation with parameters described in Table
1.

Monte Carlo calculations, as described above, are performed for the same
Alfvénic perturbations to distinguish the impact of the phase-space structure
on transport. The loss fraction as a function of time, initial and final
distribution function, and net radial displacement $|\Delta s|$ are shown in
Figure 10. We note a sudden increase in the loss fraction above
$\delta\hat{B}^{\psi}=5.58\times 10^{-4}$, the value for which island overlap
is seen to occur in Figure 9. This behavior is analogous to the observation of
critical gradient behavior on DIII-D, for which the fast-ion transport
suddenly becomes stiff above a critical gradient threshold. Similar to our
results, the critical gradient behavior on DIII-D arises due to phase-space
stochastization (Collins et al., 2017). The initial losses scale approximately
as $t^{2}$, as expected for quasilinear diffusion. As the radial distribution
of particles evolves, we see enhanced flattening near the primary resonance
surface with increasing perturbation amplitude.

Figure 10(d) displays the trapping state of lost trajectories as a function of
the perturbation amplitude. The trapped and passing status at $t=0$ and the
time of loss are assessed by comparing the trapping parameter
$\lambda=v_{\perp}^{2}/(v^{2}B)$ with the maximum field strength on the
initial and final magnetic surfaces. In order to focus on the losses due to
the Alfvénic perturbation rather than symmetry breaking or wide banana orbits,
we subtract the count in each category from the losses without an Alfvénic
perturbation. As expected, the majority of the loss increase is due to passing
orbits. There are additional losses of initially passing orbits that
transition to trapped orbits. This loss type can arise from diffusion in the
constant of motion space, leading to transformation from a passing orbit to a
barely-trapped banana orbit with a large radial width (Hsu & Sigmar, 1992;
Sigmar et al., 1992). Finally, there are losses of trapped particles that do
not originate as passing. Further analysis is needed to understand the
responsible loss mechanisms. Overall, we conclude that the kinetic Poincaré
analysis provides insight into transport characteristics, even for these
configurations with finite deviations from quasisymmetry.

(a)

(b)

(c)

(d)

Figure 10: Monte Carlo guiding center tracing calculations are performed for
the $\beta=2.5\%$ QA configuration in the presence of an $m=30$ perturbation
with parameters described in Table 1. The perturbation amplitude is increased
to study the impact of island overlap, as illustrated in Figure 9, on
transport. (a) The loss fraction as a function of time. (b) The distribution
of the total radial displacement, $|\Delta s|$, between the initial time and
final recorded time. (c) The initial and final radial distribution of
particles. (d) The increase in losses of different trajectory classes is
plotted as a function of the perturbation amplitudes. The trapped and passing
categorization is based on a comparison of the trapping parameter
$\lambda=v_{\perp}^{2}/(v^{2}B)$ with the maximum field strength on the
magnetic surface where a given trajectory is initialized or lost.

##  7 Conclusions

We have developed the theory for guiding center transport in quasisymmetric
equilibria with Alfvénic perturbations. Even if the perturbation is restricted
to a single $m$ and $n$, additional resonances may be excited due to the
coupling of the perturbation to the magnetic drifts, as discussed in Section
3. The resonance condition, phase-space island width, and island overlap
conditions are discussed for several equilibria of interest. Quasihelical
configurations have a reduced propensity for island overlap due to their
increased resonance spacing. While the potential for island overlap increases
with the poloidal mode number $m$, the effective volume of phase-space non-
integrability decreases when $m$ is large enough due to the reduced drift-
island width for higher-order drift couplings. These features are visualized
using kinetic Poincaré plots in Section 5. Although the kinetic Poincaré
section is only 2D for a perfectly quasisymmetric equilibrium, this analysis
still provides insight into the transport with finite quasisymmetry-breaking
errors. The quasisymmetry-breaking errors enhance the transport, especially in
configurations with significant unconfined guiding center trajectories such as
NCSX.

We evaluate the transport in several configurations close to quasisymmetry for
low ($m=1$), moderate ($m=15$), and high ($m=30$) mode number perturbations,
with the highest mode number being that expected for reactor conditions. The
toroidal mode number and frequency are chosen to resonate with a drift surface
in the equilibrium. We fix the perturbation amplitude to an appropriate
amplitude based on experimental measurements and modeling, $\delta
B^{r}/B_{0}\sim 10^{-3}$, finding that substantial island overlap can be
present for moderate and high-mode number perturbations. No island overlap
occurs for a quasihelical equilibrium or an equilibrium with very low magnetic
shear, such as the recently obtained vacuum equilibria with precise
quasisymmetry (Landreman & Paul, 2022). For the low-shear equilibrium, an
enhancement of the loss fraction to $\approx 20\%$ occurs due to the wide
orbit width. For the quasihelical configuration, the losses remain less than
$7\%$ for all perturbations considered. For configurations with substantial
island overlap, the losses increase to $~{}10-20\%$ due to AE-driven
transport.

Our results are consistent with similar modeling for tokamak configurations
(Hsu & Sigmar, 1992; Sigmar et al., 1992), which indicated resonant overlap
and enhanced transport for $\delta B^{r}/B_{0}\sim 10^{-3}$. More recent work
has highlighted the impact of symmetry-breaking on AE-driven transport in a QH
configuration (White & Duarte, 2023). Our results agree qualitatively with
their conclusion that unconfined particle orbits enhance diffusive losses.
However, our results are not in quantitative agreement with their findings of
substantial diffusive losses for AE perturbation amplitudes $\delta
B^{r}/B_{0}\sim 10^{-6}$ for a QH equilibrium with precise levels of
quasisymmetry (Landreman et al., 2021). Further work is necessary in order to
resolve this discrepancy.

As the amplitude of the perturbation is increased, island overlap is observed,
leading to a rapid transition to stiff fast ion transport. This result is
similar to observations of critical gradient transport in DIII-D (Collins et
al., 2016, 2017) and indicates the potential applicability of quasilinear
models (Duarte et al., 2017; Gorelenkov et al., 2018; Duarte et al., 2019) to
predict the saturated AE amplitude in stellarators. Since the phase-space
island width grows with the amplitude of the magnetic drifts, island overlap
is more likely to occur at larger energies.

Our results suggest several avenues to reduce AE-driven transport in
quasisymmetric configurations:

  * •

Development of quasihelical configurations, which avoid drift-island overlap
due to increased resonance spacing;

  * •

Avoiding low magnetic shear, which manifests as wide drift island widths for
passing particles;

  * •

Evaluating metrics from the unperturbed equilibrium (e.g., the Bessel coupling
parameters $J_{l}(\eta_{1})$) as a proxy for drift island width within a
stellarator optimization loop.

Finally, we remark that many aspects of Alfvén eigenmodes were not considered
in this study and will affect the transport. To gain a basic understanding of
transport, a resonant AE was chosen to have a uniform radial structure. While
the radial perturbation amplitude was assumed to be held fixed when comparing
across configurations, in practice, the mode structure and amplitude will also
depend on the mode numbers and other physical parameters through the AE’s
nonlinear evolution. In the future, we plan to build on this work to obtain a
more complete picture of the evolution of AEs in stellarator reactor
scenarios. For example, the analysis tools described here, such as kinetic
Poincaré plots, could be used to study the nonlinear evolution of phase-space
structures such as hole-clump pairs (Bierwage et al., 2021) and zonal
structures (Zonca et al., 2015) in quasisymmetric configurations.

## Acknowledgements

The authors would like to acknowledge discussions with Vinicius Duarte and
Roscoe White. We acknowledge funding through the U. S. Department of Energy,
under contract No. DE-SC0016268, and through the Simons Foundation
collaboration “Hidden Symmetries and Fusion Energy,” Grant No. 601958.

## References

  * Bader et al. (2021) Bader, A., Anderson, D. T., Drevlak, M., Faber, B. J., Hegna, C. C., Henneberg, S., Landreman, M., Schmitt, J. C., Suzuki, Y. & Ware, A. 2021 Modeling of energetic particle transport in optimized stellarators. Nuclear Fusion 61 (11), 116060.
  * Bader et al. (2019) Bader, A., Drevlak, M., Anderson, D. T., Faber, B. J., Hegna, C. C., Likin, K. M., Schmitt, J. C. & Talmadge, J. N. 2019 Stellarator equilibria with reactor relevant energetic particle losses. Journal of Plasma Physics 85 (5).
  * Bierwage et al. (2021) Bierwage, A., White, R. B. & Duarte, V. N. 2021 On the effect of beating during nonlinear frequency chirping. Plasma and Fusion Research 16, 1403087–1403087.
  * Collins et al. (2016) Collins, C. S., Heidbrink, W. W., Austin, M. E., Kramer, G. J., Pace, D. C., Petty, C. C., Stagner, L., Van Zeeland, M. A., White, R. B., Zhu, Y. B. & others 2016 Observation of critical-gradient behavior in Alfvén-eigenmode-induced fast-ion transport. Physical Review Letters 116 (9), 095001.
  * Collins et al. (2017) Collins, C. S., Heidbrink, W. W., Podestà, M., White, R. B., Kramer, G. J., Pace, D. C., Petty, C. C., Stagner, L., Van Zeeland, M. A., Zhu, Y. B. & others 2017 Phase-space dependent critical gradient behavior of fast-ion transport due to Alfvén eigenmodes. Nuclear Fusion 57 (8), 086005.
  * Crocker et al. (2013) Crocker, N. A., Fredrickson, E. D., Gorelenkov, N. N., Peebles, W. A., Kubota, S., Bell, R. E., Diallo, A., LeBlanc, B. P., Menard, J. E., Podesta, M. & others 2013 Internal amplitude, structure and identification of compressional and global Alfvén eigenmodes in NSTX. Nuclear Fusion 53 (4), 043017.
  * Das et al. (2016) Das, S., Saiki, Y., Sander, E. & Yorke, J. A. 2016 Quasiperiodicity: rotation numbers. In The Foundations of Chaos Revisited: From Poincaré to Recent Advancements, pp. 103–118. Springer.
  * Deng et al. (2009) Deng, C. B., Brower, D. L., Breizman, B. N., Spong, D. A., Almagri, A. F., Anderson, D. T., Anderson, F. S. B., Ding, W. X., Guttenfelder, W., Likin, K. M. & others 2009 Energetic-electron-driven instability in the Helically Symmetric Experiment. Physical Review Letters 103 (2), 025003.
  * Duarte et al. (2017) Duarte, V. N., Berk, H. L., Gorelenkov, N. N., Heidbrink, W. W., Kramer, G. J., Nazikian, R., Pace, D. C., Podesta, M., Tobias, B. J. & Van Zeeland, M. A. 2017 Prediction of nonlinear evolution character of energetic-particle-driven instabilities. Nuclear Fusion 57 (5), 054001.
  * Duarte et al. (2019) Duarte, V. N., Gorelenkov, N. N., White, R. B. & Berk, H. L. 2019 Collisional resonance function in discrete-resonance quasilinear plasma systems. Physics of Plasmas 26 (12), 120701.
  * Fehér (2014) Fehér, T. B. 2014 Simulation of the interaction between Alfvén waves and fast particles. PhD thesis, University of Greifswald.
  * Garren & Boozer (1991) Garren, D. A. & Boozer, A. H. 1991 Magnetic field strength of toroidal plasma equilibria. Physics of Fluids B: Plasma Physics 3 (10), 2805–2821.
  * Gorelenkov et al. (2018) Gorelenkov, N. N., Duarte, V. N., Podesta, M. & Berk, H. L. 2018 Resonance broadened quasi-linear (RBQ) model for fast ion distribution relaxation due to Alfvénic eigenmodes. Nuclear Fusion 58 (8), 082016.
  * Gorelenkov et al. (2014) Gorelenkov, N. N., Pinches, S. D. & Toi, K. 2014 Energetic particle physics in fusion research in preparation for burning plasma experiments. Nuclear Fusion 54 (12), 125001.
  * Heidbrink et al. (2008) Heidbrink, W. W., Van Zeeland, M. A., Austin, M. E., Burrell, K. H., Gorelenkov, N. N., Kramer, G. J., Luo, Y., Makowski, M. A., McKee, G. R., Muscatello, C. & others 2008 Central flattening of the fast-ion profile in reversed-shear DIII-D discharges. Nuclear Fusion 48 (8), 084001.
  * Helander et al. (2012) Helander, P., Beidler, C. D., Bird, T. M., Drevlak, M., Feng, Y., Hatzky, R., Jenko, F., Kleiber, R., Proll, J. H. E., Turkin, Y. & others 2012 Stellarator and tokamak plasmas: a comparison. Plasma Physics and Controlled Fusion 54 (12), 124009.
  * Hsu & Sigmar (1992) Hsu, C.-T. & Sigmar, D. J. 1992 Alpha-particle losses from toroidicity-induced Alfvén eigenmodes. Part I: Phase-space topology of energetic particle orbits in tokamak plasma. Physics of Fluids B: Plasma Physics 4 (6), 1492–1505.
  * Koniges et al. (2003) Koniges, A. E., Grossman, A., Fenstermacher, M., Kisslinger, J., Mioduszewski, P., Rognlien, T., Strumberger, E. & Umansky, M. 2003 Magnetic topology of a candidate NCSX plasma boundary configuration. Nuclear Fusion 43 (2), 107.
  * Kruger et al. (1998) Kruger, S. E., Hegna, C. C. & Callen, J. D. 1998 Generalized reduced magnetohydrodynamic equations. Physics of Plasmas 5 (12), 4169–4182.
  * Ku et al. (2008) Ku, L. P., Garabedian, P. R., Lyon, J., Turnbull, A., Grossman, A., Mau, T. K., Zarnstorff, M. & ARIES Team 2008 Physics design for ARIES-CS. Fusion Science and Technology 54 (3), 673–693.
  * Landreman et al. (2022) Landreman, M., Buller, S. & Drevlak, M. 2022 Optimization of quasi-symmetric stellarators with self-consistent bootstrap current and energetic particle confinement. Physics of Plasmas 29 (8), 082501.
  * Landreman et al. (2021) Landreman, M., Medasani, B., Wechsung, F., Giuliani, A., Jorge, R. & Zhu, C. 2021 SIMSOPT: A flexible framework for stellarator optimization. Journal of Open Source Software 6 (65), 3525.
  * Landreman & Paul (2022) Landreman, M. & Paul, E. J. 2022 Magnetic fields with precise quasisymmetry for plasma confinement. Physical Review Letters 128 (3), 035001.
  * LeViness et al. (2022) LeViness, A., Schmitt, J. C., Lazerson, S. A., Bader, A., Faber, B. J., Hammond, K. C. & Gates, D. A. 2022 Energetic particle optimization of quasi-axisymmetric stellarator equilibria. Nuclear Fusion 63 (1), 016018.
  * Littlejohn (1983) Littlejohn, R. G. 1983 Variational principles of guiding centre motion. Journal of Plasma Physics 29 (1), 111–125.
  * Littlejohn (1985) Littlejohn, R. G. 1985 Differential forms and canonical variables for drift motion in toroidal geometry. The Physics of Fluids 28 (6), 2015–2016.
  * Melnikov et al. (2014) Melnikov, A. V., Ochando, M., Ascasibar, E., Castejon, F., Cappa, A., Eliseev, L. G., Hidalgo, C., Krupnik, L. I., Lopez-Fraguas, A., Liniers, M. & others 2014 Effect of magnetic configuration on frequency of NBI-driven Alfvén modes in TJ-II. Nuclear Fusion 54 (12), 123002.
  * Mynick (1993a) Mynick, H. E. 1993a Stochastic transport of MeV ions by low-n magnetic perturbations. Physics of Fluids B: Plasma Physics 5 (7), 2460–2467.
  * Mynick (1993b) Mynick, H. E. 1993b Transport of energetic ions by low-n magnetic perturbations. Physics of Fluids B: Plasma Physics 5 (5), 1471–1481.
  * Mynick et al. (2002) Mynick, H. E., Pomphrey, N. & Ethier, S. 2002 Exploration of stellarator configuration space with global search methods. Physics of Plasmas 9 (3), 869–876.
  * Najmabadi et al. (2008) Najmabadi, F., Raffray, A. R., Abdel-Khalik, S. I., Bromberg, L., Crosatti, L., El-Guebaly, L., Garabedian, P. R., Grossman, A. A., Henderson, D., Ibrahim, A. & others 2008 The ARIES-CS compact stellarator fusion power plant. Fusion Science and Technology 54 (3), 655–672.
  * Nazikian et al. (1997) Nazikian, R., Fu, G. Y., Batha, S. H., Bell, M. G., Bell, R. E., Budny, R. V., Bush, C. E., Chang, Z., Chen, Y., Cheng, C. Z. & others 1997 Alpha-particle-driven toroidal Alfvén eigenmodes in the tokamak fusion test reactor. Physical Review Letters 78 (15), 2976.
  * Nishimura et al. (2013) Nishimura, S., Todo, Y., Spong, D. A., Suzuki, Y. & Nakajima, N. 2013 Simulation study of Alfvén-eigenmode-induced energetic ion transport in LHD. Plasma and Fusion Research 8, 2403090–2403090.
  * Paul et al. (2022) Paul, E. J., Bhattacharjee, A., Landreman, M., Alex, D., Velasco, J. L. & Nies, R. 2022 Energetic particle loss mechanisms in reactor-scale equilibria close to quasisymmetry. Nuclear Fusion 62 (12), 126054.
  * Rahbarnia et al. (2020) Rahbarnia, K., Thomsen, H., Schilling, J., Vaz Mendes, S., Endler, M., Kleiber, R., Könies, A., Borchardt, M., Slaby, C., Bluhm, T. & others 2020 Alfvénic fluctuations measured by in-vessel Mirnov coils at the Wendelstein 7-X stellarator. Plasma Physics and Controlled Fusion 63 (1), 015005.
  * Sigmar et al. (1992) Sigmar, D. J., Hsu, C. T., White, R. & Cheng, C. Z. 1992 Alpha-particle losses from toroidicity-induced Alfvén eigenmodes. Part II: Monte Carlo simulations and anomalous alpha-loss processes. Physics of Fluids B: Plasma Physics 4 (6), 1506–1516.
  * Spong et al. (2010) Spong, D. A., D’azevedo, E. & Todo, Y. 2010 Clustered frequency analysis of shear Alfvén modes in stellarators. Physics of Plasmas 17 (2), 022106.
  * Spong et al. (2017) Spong, D. A., Holod, I., Todo, Y. & Osakabe, M. 2017 Global linear gyrokinetic simulation of energetic particle-driven instabilities in the LHD stellaratora. Nuclear Fusion 57 (8), 086018.
  * Takechi et al. (2002) Takechi, M., Matsunaga, G., Toi, K., Isobe, M., Minami, T., Tanaka, K., Nishimura, S., Takahashi, C., Okamura, S. & Matsuoka, K. 2002 Transition of toroidal Alfven eigenmode to global Alfven eigenmode in CHS heliotron/torsatron plasmas heated by neutral beam injection. Journal of Plasma and Fusion Research 78 (12), 1273–1274.
  * Todo et al. (2017) Todo, Y., Seki, R., Spong, D. A., Wang, H., Suzuki, Y., Yamamoto, S., Nakajima, N. & Osakabe, M. 2017 Comprehensive magnetohydrodynamic hybrid simulations of fast ion driven instabilities in a Large Helical Device experiment. Physics of Plasmas 24 (8), 081203.
  * Toi et al. (2011) Toi, K., Ogawa, K., Isobe, M., Osakabe, M., Spong, D. A. & Todo, Y. 2011 Energetic-ion-driven global instabilities in stellarator/helical plasmas and comparison with tokamak plasmas. Plasma Physics and Controlled Fusion 53 (2), 024008.
  * Varela et al. (2021) Varela, J., Shimizu, A., Spong, D. A., Garcia, L. & Ghai, Y. 2021 Study of the Alfven eigenmodes stability in CFQS plasma using a Landau closure model. Nuclear Fusion 61 (2), 026023.
  * Varela et al. (2017) Varela, J., Spong, D. A. & Garcia, L. 2017 Analysis of Alfvén eigenmode destabilization by energetic particles in Large Helical Device using a Landau-closure model. Nuclear Fusion 57 (4), 046018.
  * Weller et al. (1994) Weller, A., Spong, D. A., Jaenicke, R., Lazaros, A., Penningsfeld, F. P. & Sattler, S. 1994 Neutral beam driven global Alfvén eigenmodes in the Wendelstein W7-AS stellarator. Physical Review Letters 72 (8), 1220.
  * White et al. (2022) White, R., Bierwage, A. & Ethier, S. 2022 Poor confinement in stellarators at high energy. Physics of Plasmas 29 (5), 052511.
  * White (2011) White, R. B. 2011 Modification of particle distributions by magnetohydrodynamic instabilities II. Plasma Physics and Controlled Fusion 53 (8), 085018.
  * White (2012) White, R. B. 2012 Modification of particle distributions by MHD instabilities I. Communications in Nonlinear Science and Numerical Simulation 17 (5), 2200–2214.
  * White & Duarte (2023) White, R. B. & Duarte, V. N. 2023 Assessment of radial transport induced by Alfvénic resonances in tokamaks and stellarators. Physics of Plasmas 30 (1), 012502.
  * White et al. (1983) White, R. B., Goldston, R. J., McGuire, K., Boozer, A. H., Monticello, D. A. & Park, W. 1983 Theory of mode-induced beam particle loss in tokamaks. The Physics of Fluids 26 (10), 2958–2965.
  * Yamamoto et al. (2007) Yamamoto, S., Nagasaki, K., Suzuki, Y., Mizuuchi, T., Okada, H., Kobayashi, S., Blackwell, B., Kondo, K., Motojima, G., Nakajima, N. & others 2007 Observation of magnetohydrodynamic instabilities in Heliotron J plasmas. Fusion Science and Technology 51 (1), 92–96.
  * Zonca et al. (2015) Zonca, F., Chen, L., Briguglio, S., Fogaccia, G., Vlad, G. & Wang, X. 2015 Nonlinear dynamics of phase space zonal structures and energetic particle physics in fusion plasmas. New Journal of Physics 17 (1), 013052.

##  Appendix A Analysis of the resonance condition

Under the assumption of a vacuum magnetic field, the perturbed radial equation
of motion reads,

$\displaystyle\delta\dot{\psi}$
$\displaystyle=\frac{\bm{B}_{0}\times\nabla\delta\Phi\cdot\nabla\psi}{B_{0}^{2}}+v_{\|}\frac{\delta\bm{B}\cdot\nabla\psi}{B_{0}}=-\left(\omega+(\iota
m-n)\dot{\zeta}\right)\frac{\delta\Phi_{,\theta}}{\omega},$ (27)

where subscript commas denote partial derivatives. The toroidal motion,
$\dot{\zeta}=B_{0}v_{\|}/G$, is not modified by the perturbation. The first
term accounts for the perturbed $\bm{E}\times\bm{B}$ drift, while the second
term accounts for streaming along the perturbed magnetic field. Given the form
for the perturbed potential (4), the perturbed radial drift is proportional to
$\cos(\eta)$. We begin our analysis by expressing $\cos(\eta)$ along the
unperturbed trajectory in terms of Bessel functions (31). This expression is
used to evaluate the perturbed radial drift in (32)-(35) and the resulting
phase-space island width (36).

We begin with the expressions for the unperturbed trajectories (13). Due to
the relatively small magnitude of the periodic contributions, $\chi_{j}$ and
$\zeta_{j}$, in comparison to the secular terms, $\omega_{\chi}$ and
$\omega_{\zeta}$, the trajectories can be approximated as:

$\displaystyle\left\\{\begin{array}[]{l}\displaystyle\chi\approx\chi^{0}+\omega_{\chi}t+\sum_{j\neq
0}\frac{\chi_{j}}{j\omega_{\chi}}\sin(j\left[\chi^{0}+\omega_{\chi}t\right]),\\\
\displaystyle\zeta\approx\zeta^{0}+\omega_{\zeta}t+\sum_{j\neq
0}\frac{\zeta_{j}}{j\omega_{\chi}}\sin(j\left[\chi^{0}+\omega_{\chi}t\right]),\end{array}\right.$
(30)

where $\chi^{0}=\chi(t=0)$ and $\zeta^{0}=\zeta(t=0)$.

To evaluate the radial drift due to the perturbation (27), we evaluate the
phase along the unperturbed trajectory

$\displaystyle\eta(t)=\eta^{0}+\left(m\omega_{\chi}-(n-Nm)\omega_{\zeta}+\omega\right)t+\sum_{j\neq
0}\left[\frac{m\chi_{j}}{j\omega_{\chi}}-(n-Nm)\frac{\zeta_{j}}{j\omega_{\chi}}\right]\sin(j\left[\chi^{0}+\omega_{\chi}t\right]),$

where $\eta^{0}=m\theta^{0}-n\zeta^{0}$ is the initial phase. Under the
assumption that one term, $j=j^{\prime}$, is dominant in the expression for
the drifts (13), the phase factor $\cos(\eta)$ can be expressed in terms of
Bessel functions of the first kind,

$\displaystyle\cos(\eta)$
$\displaystyle=\sum_{k}J_{k}\left(\eta_{j^{\prime}}\right)\cos\left(\eta^{0}+\Omega_{kj^{\prime}}t\right),$
(31)

where
$\eta_{j^{\prime}}=\left(m\chi_{j^{\prime}}-(n-Nm)\zeta_{j^{\prime}}\right)/(j^{\prime}\omega_{\chi})$
and
$\Omega_{kj^{\prime}}=(m+kj^{\prime})\omega_{\chi}-(n-Nm)\omega_{\zeta}+\omega$.

By expressing the instantaneous perturbed radial drift (27) in terms of
$\cos(\eta)$ along the unperturbed trajectory (31) we obtain the expression:

$\delta\dot{\psi}=\sum_{k}\Bigg{(}\psi_{kj^{\prime}}^{0}\cos\left(\Omega_{kj^{\prime}}t+\eta^{0}\right)+\psi_{kj^{\prime}}^{+}\cos\left(\Omega_{kj^{\prime}}t+\eta^{0}-j^{\prime}\chi^{0}\right)\\\
+\psi_{kj^{\prime}}^{-}\cos\left(\Omega_{kj^{\prime}}t+\eta^{0}+j^{\prime}\chi^{0}\right)\Bigg{)}$
(32)

with amplitudes given by,

$\displaystyle\left\\{\begin{array}[]{l}\psi_{kj^{\prime}}^{0}=-m\hat{\Phi}J_{k}(\eta_{j^{\prime}})\left(1+(\iota
m-n)\frac{\omega_{\zeta}}{\omega}\right),\\\
\displaystyle\psi_{kj^{\prime}}^{\pm}=-\frac{m\hat{\Phi}J_{k\pm
1}(\eta_{j^{\prime}})}{2}(\iota
m-n)\frac{\zeta_{j^{\prime}}}{\omega}.\end{array}\right.$ (35)

For a net drift to occur, a resonance condition $\Omega_{kj^{\prime}}=0$ must
be satisfied. The Fourier amplitudes $\psi_{kj^{\prime}}^{0}$ and
$\psi_{kj^{\prime}}^{\pm}$ then determine the amplitude of the radial
transport. The associated full island width is given by,

$\displaystyle
w^{\psi}_{kj^{\prime}}=2\sqrt{\left|\frac{\psi_{kj^{\prime}}}{\Omega_{kj^{\prime}}^{\prime}(\psi)}\right|}\approx
2\sqrt{\left|\frac{\psi_{kj^{\prime}}}{(m+kj^{\prime})\omega_{\theta}^{\prime}(\psi)}\right|},$
(36)

where $\psi_{kj^{\prime}}$ is taken to be $\psi_{kj^{\prime}}^{0}$ or
$\psi_{kj^{\prime}}^{\pm}$.

We now consider the scaling of the cosine amplitudes (35) with respect to the
radial perturbed magnetic field amplitude and magnetic drifts. If a resonance
condition $\Omega_{kj^{\prime}}=0$ is satisfied, the corresponding radial
perturbed field scales with the perturbation amplitude and mode numbers as
$\delta\hat{B}^{\psi}\sim m(\iota m-n)\hat{\Phi}$. Using the resonance
condition, the factor appearing in $\psi_{kj^{\prime}}^{0}$ can be expressed
as $1+(\iota
m-n)\omega_{\zeta}/\omega=(-kj^{\prime}\omega_{\chi}+m(\iota\omega_{\zeta}-\omega_{\theta}))/\omega$.
Since the quantity $\iota-\omega_{\theta}/\omega_{\zeta}$ scales with the
magnitude of the magnetic drifts, at low energy, this factor approaches
$kj^{\prime}\omega_{\chi}/\omega$. Thus,
$\psi_{0}^{0}\sim(\iota-\omega_{\theta}/\omega_{\zeta})J_{0}(\eta_{j^{\prime}})\delta\hat{B}^{\psi}$.
The scaling of each harmonic is summarized in Table 3. Considering the small
argument limit of $\psi_{0}^{\pm}$, this harmonic will scale with two factors
of the magnetic drifts through both $\eta_{j^{\prime}}$ and
$\zeta_{j^{\prime}}$. Thus $\psi_{0}^{0}$ will dominate over $\psi_{0}^{\pm}$
for $k=0$. Considering $k\neq 0$, $\psi_{k}^{0}$ will have the same scaling
with the magnetic drifts as $\psi_{k}^{-}$ for $k>0$, while $\psi_{k}^{+}$
will be smaller in magnitude by a factor of the magnetic drifts. The magnitude
of the harmonics decreases strongly with $|k|$ due to the Bessel coupling
parameter. When considering the dependence on the magnetic drifts, the most
significant radial transport will arise from $\psi_{0}^{0}$, $\psi_{\pm
j^{\prime}}^{0}$, and $\psi_{\mp j^{\prime}}^{\pm}$. Considering the limit of
large mode number perturbations, $\psi_{\pm j^{\prime}}^{0}$ will dominate due
to its dependence on $\eta_{j^{\prime}}$.

$\psi_{0}^{0}$ | $\psi_{k\neq 0}^{0}$ | $\psi_{k}^{\pm}$
---|---|---
$(\iota-\omega_{\theta}/\omega_{\zeta})J_{0}(\eta_{j^{\prime}})\delta\hat{B}^{\psi}$ | $J_{k}(\eta_{j^{\prime}})\delta\hat{B}^{\psi}$ | $J_{k\pm 1}(\eta_{j^{\prime}})\zeta_{j^{\prime}}\delta\hat{B}^{\psi}$
Table 3: Scaling of cosine harmonic amplitudes of the perturbed radial drift
(35) with the radial perturbed magnetic field, magnetic drifts, and mode
numbers. The quantities $(\iota-\omega_{\theta}/\omega_{\zeta})$ and
$\zeta_{j^{\prime}}$ scale with the drift amplitude. The quantity
$\eta_{j^{\prime}}$ scales with both the drift amplitude and the mode numbers.