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

Abstract

Significant progress has been made in designing magnetic fields that provide excellent confinement of the guiding-centre 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-centre 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, Buller & Drevlak 2022; Landreman & Paul 2022). However, reducing the guiding-centre 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 modelling indicates that Alfvénic activity remains present even at a fast-ion beta of ${\approx }0.05\,\%$ (Varela, Spong & Garcia 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-centre 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. The AE stability of the Chinese First Quasisymmetric Stellarator (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 (Hsu & Sigmar 1992; Sigmar et al. 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 (Hsu & Sigmar 1992; Mynick 1993a,b). 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 QS helicity and deviations from QS impact transport. Recent simulations (White, Bierwage & Ethier 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 modelling (Fehér 2014; Spong et al. 2017; Todo 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 QS. 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 modelling. 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 QS 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 QS.

In § 2, we outline the guiding-centre equations of motion with Alfvénic perturbations. In § 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 § 4, and resonant perturbations are identified. Kinetic Poincaré plots are employed in § 5 to assess the formation of phase-space islands and chaos in the presence of the resonant perturbations. A Monte Carlo guiding-centre transport analysis is performed in § 6. We conclude in § 7.

2. Guiding-centre motion in the presence of an Alfvénic perturbation

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

(2.1)\begin{equation} L(\boldsymbol{R},\dot{\boldsymbol{R}},v_{\|}) = q\left(\boldsymbol{A} + \frac{M v_{\|}}{q B}\boldsymbol{B} \right) \boldsymbol{\cdot} \dot{\boldsymbol{R}} - \frac{M v_{\|}^2}{2} - \mu B - q\varPhi, \end{equation}

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

(2.2)\begin{equation} \boldsymbol{B} = \boldsymbol{B}_0 + \boldsymbol{\nabla} \times (\alpha \boldsymbol{B}_0 ),\end{equation}

for parameter $\alpha $ while the scalar potential vanishes in the equilibrium, $\varPhi = \delta \varPhi$. With the reduced magnetohydrodynamic (MHD) assumption (Kruger, Hegna & Callen 1998) – $k_{\|}/k_{\perp } \ll 1$, where $k_{\|}$ is the characteristic parallel wavenumber and $k_{\perp }$ is the characteristic perpendicular wavenumber associated with perturbed quantities – the form of the perturbed field (2.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 $\boldsymbol {B}_0 \boldsymbol {\cdot }\delta \boldsymbol {E} =0$, which implies

(2.3)\begin{equation} \nabla_{\|} \delta \varPhi ={-} B_0\frac{\partial{\alpha}}{\partial{t}}.\end{equation}

The perturbed fields could be computed from a global AE solver such as AE3D (Spong, D'azevedo & Todo 2010). For the fundamental studies here, we instead assume a single harmonic perturbation of the form

(2.4)\begin{equation} \delta \varPhi = \hat{\varPhi}(\psi) \sin(\omega t + m \theta - n \zeta), \end{equation}

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

(2.5)\begin{equation} \left. \begin{array}{c@{}} \boldsymbol{B}_0 = \boldsymbol{\nabla} \psi \times \boldsymbol{\nabla} \theta - \iota(\psi) \boldsymbol{\nabla} \psi \times \boldsymbol{\nabla} \zeta, \\ \boldsymbol{B}_0 = G(\psi) \boldsymbol{\nabla} \zeta + I(\psi) \boldsymbol{\nabla} \theta + K(\psi,\theta,\zeta) \boldsymbol{\nabla} \psi, \end{array} \right\} \end{equation}

where $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 (2.3), which we write schematically as

(2.6)\begin{equation} \alpha = \hat{\alpha}(\psi) \sin(\eta).\end{equation}

We impose $\hat {\varPhi }$ 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 $M v_{\|}/(qB_0)$ is small in comparison with the characteristic length scale of the equilibrium and $\delta B=0$, the perturbed Lagrangian reads (Littlejohn 1985)

(2.7)\begin{equation} L(\boldsymbol{R},\dot{\boldsymbol{R}},v_{\|}) = q\left(\boldsymbol{A}_0 + \alpha \boldsymbol{B}_0 + \frac{Mv_{\|}}{qB_0} \boldsymbol{B}_0 \right) \boldsymbol{\cdot} \dot{\boldsymbol{R}} - \frac{M v_{\|}^2}{2} - \mu B_0 - q\delta \varPhi. \end{equation}

The resulting equations of motion are integrated in Boozer coordinates using the SIMSOPT stellarator optimization and modelling 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

(2.8)\begin{equation} \frac{\delta \boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi}{B_0 |\boldsymbol{\nabla} \psi|} \approx \frac{\boldsymbol{\nabla} \alpha \times \boldsymbol{B}_0 \boldsymbol{\cdot} \boldsymbol{\nabla} \psi}{r B_0^2} = \frac{\hat{\alpha} (m G -n I)}{r (G + \iota I)} \cos(\eta). \end{equation}

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 \boldsymbol {B} \boldsymbol {\cdot } \boldsymbol {\nabla } \psi /(B_0 |\boldsymbol {\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 \ne 0$ for quasihelical symmetry). Each equilibrium we consider in § 4 is sufficiently close to QS that such an assumption provides valuable insight into the resulting dynamics.

The radial drift over the unperturbed trajectories is analysed 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-centre equations of motion. Because $\boldsymbol {B}_0 \boldsymbol {\cdot } \delta \boldsymbol {E} = 0$, particles continue to follow perturbed field lines to lowest order in the gyroradius but experience an additional $\boldsymbol {E} \times \boldsymbol {B}$ drift. Under the assumption of QS and stellarator symmetry, the unperturbed equations of motion can be written schematically as

(3.1)\begin{equation} \left. \begin{array}{c@{}} \displaystyle \dot{\chi} = \omega_{\chi} + \sum_{j\ne 0} \chi_j \cos(j\chi), \\ \displaystyle \dot{\zeta} = \omega_{\zeta} + \sum_{j\ne 0} \zeta_j \cos(j\chi), \end{array} \right\} \end{equation}

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 \, \textrm {d} t A /T$. The summations represent the periodic contributions from the drifts.

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

(3.2)\begin{equation} \varOmega_l = (m+l)\omega_{\chi} - (n - Nm)\omega_{\zeta} + \omega = 0. \end{equation}

As discussed in Appendix A, the integer $l$ arises due to coupling through the $\chi$ dependence of the magnetic drifts. If the drift dynamics (3.1) is dominated by a particular cosine harmonic with integer $j'$, then $l$ is assumed to be an integer multiple of $j'$. To simplify the analysis, the $j' = 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

(3.3)\begin{equation} w^{\psi}_{l} = 2\sqrt{\left | \frac{\psi_{l}}{\varOmega_{l}'(\psi)} \right |} \approx 2\sqrt{\left | \frac{\psi_{l\ne 0}}{(m + l) \omega_{\theta}'(\psi)} \right |},\end{equation}

where $\psi _l$ is a cosine harmonic of the radial perturbed drift and we have made the approximation $h'(\psi )\approx \omega _{\theta }'(\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 (\varOmega _l t + \eta ^0)$ and $\cos (\varOmega _l t + \eta ^0 \mp \chi ^0)$ with coefficients given by $\psi _l^0$ and $\psi _l^{\pm }$; see (A5). 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

(3.4)\begin{equation} \left. \begin{array}{c@{}} \psi_{0}^0 \sim (\iota - \omega_{\theta}/\omega_{\zeta}) J_0(\eta_{1}) \delta \hat{B}^{\psi}, \\ \psi_l^0 \sim J_l(\eta_{1})\delta \hat{B}^{\psi}, \\ \psi_{l}^{{\pm}} \sim J_{l\pm 1}(\eta_{1}) \zeta_{1} \delta \hat{B}^{\psi}, \end{array} \right\} \end{equation}

where $J_l$ are the Bessel functions of the first kind, $\eta _{1} = (m\chi _{1}- (n-Nm)\zeta _{1})/\omega _{\chi }$ with $\chi _1$ and $\zeta _1$ defined through (3.1), and $\omega _{\theta } = \langle \dot {\theta } \rangle$. The full expressions for general $j'$ 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

(3.5)\begin{equation} h = \frac{n - Nm - \omega/\omega_{\zeta}}{m + l}+N.\end{equation}

For a given primary resonance $l$ at $h = h_0$, additional sideband resonances may be excited for other drift harmonics $l'$ corresponding to neighbouring periodic orbits. The spacing between the $l$ and $l+1$ resonances is given by

(3.6)\begin{equation} (\Delta \psi)_l = \frac{1}{h'(\psi)}\frac{h_0 - N}{m + l+1}. \end{equation}

The ratio between the island width and the resonance spacing

(3.7)\begin{equation} \frac{w_l^{\psi}}{(\Delta \psi)_l} \approx \frac{m + l + 1}{(h_0 - N)\omega_{\zeta}}\sqrt{\left | \frac{\psi_l \omega_{\theta}'(\psi)}{m + l} \right |}, \end{equation}

provides a conservative estimate for the island overlap criterion, $w_l^{\psi }/(\Delta \psi )_l \gtrsim 1$. Given the scaling of $\psi _l$ (3.4), 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 }'(\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 QS: $\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 QA NCSX li383 equilibrium (Mynick, Pomphrey & Ethier 2002; Koniges et al. 2003). 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 QS error

(4.1)\begin{equation} f_{\rm QS}(s) = \frac{\sqrt{\displaystyle \sum_{Mn\ne Nm} (B_{m,n}^c(s))^2}}{\sqrt{(B_{0,0}^c(s))^2}},\end{equation}

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 )$.

Figure 1. (a) Rotational transform and (b) QS error (4.1) 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-centre 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{\rm \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 non-symmetric modes are artificially suppressed for this frequency analysis so that all trajectories lie on Kolmogorov-Arnold-Moser (KAM) surfaces.

For the subsequent analysis, resonant perturbations are intentionally imposed, satisfying the condition (3.2) 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, Pinches & Toi 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}\unicode{x2013} 10^{-3}$, consistent with LHD modelling (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 toroidal Alfvén eigenmodes (TAEs) on guiding-centre 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 (3.2) can be expressed as $\omega /\omega _{\zeta } = p'/q$ for some integer $p'$. The toroidal mode number $n$ and the frequency parameter $p'$ are chosen to satisfy the $l = 1$ resonance condition (3.2), given that the $l = \pm 1$ resonance is predicted to give rise to the most substantial radial transport as described in § 3. The value of $n$ is chosen to yield small magnitudes of $|p'|$, corresponding with lower-frequency perturbations. The mode parameters are summarized in table 1.

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

For each configuration, the profile of the effective orbit helicity $h$ (3.5) 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 (3.5) 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 coloured horizontal lines. As $m$ increases, the resonance spacing decreases according to (3.6). 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.

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 coloured 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; (a) $\beta = 2.5\,\%$ QA, (b) $\beta = 2.5\,\%$ QH, (c) vacuum QA and (d) NCSX.

5. Kinetic Poincaré plots

Given a magnetic field with exact QS, 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

(5.1)\begin{equation} P_{\zeta} = (G + N I) \left(\frac{m v_{\|}}{B_0} + q \alpha \right) + q ( N\psi - \psi_P), \end{equation}

and energy

(5.2)\begin{equation} E = \frac{m v_{\|}^2}{2} + \mu B_0 + q \delta \varPhi, \end{equation}

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

(5.3)\begin{equation} \bar{E}_n = (n - N m)E - \omega P_{\zeta}.\end{equation}

Given the velocity-space parameters – $\mu$, $\bar {E}_n$, and $\textrm {sign}(v_{\|})$ – and specification of the position in space and time – $t$, $s$, $\theta$ and $\zeta$ – (5.3) provides a nonlinear equation for $v_{\|}$ (Hsu & Sigmar 1992).

Given the form for the phase factor of the perturbation (2.4), the resulting equations of motion depend on $(s,\chi,\zeta - \bar {\omega }_n t, v_{\|})$, where $\bar {\omega }_n = \omega /(n - N m)$. The resulting motion is generally four-dimensional. If purely passing particles are considered, then ${sign}(v_{\|})$ is fixed, $v_{\|}$ can be computed from the other coordinates as described above, and the resulting motion becomes three-dimensional. For a fixed value of $\bar {E}_n$ and $\mu$, a kinetic Poincaré map $M(s,\chi ) \rightarrow (s',\chi ')$ is constructed by moving with the wave frame to eliminate the time dependence. Guiding-centre trajectories are followed until they intersect a plane of constant $\zeta - \bar {\omega }_nt$. Therefore, the resulting Poincaré section is two-dimensional.

When QS is broken in the equilibrium, $\bar {E}_n$ is no longer precisely conserved. Nonetheless, if the QS 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 QS-breaking modes will be assessed in the following section.

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

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; (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 and (l) $m = 30$, NCSX.

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.

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 $\varOmega _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 re-enters. (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 (3.3) and (3.4) 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 with 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 § 6.

5.1. Impact of QS deviations

We now investigate the impact of the finite QS deviations, quantified in figure 1, on the kinetic Poincaré analysis. When QS deviations are present, the motion of passing particles becomes four-dimensional $(s,\chi,\zeta,t)$, and the Poincaré section becomes three-dimensional, $M(s,\chi,\zeta ) \rightarrow (s',\chi ',\zeta ')$. Nonetheless, we can still visualize the Poincaré map in the $(s,\chi )$ plane to assess the structure of phase space.

In figures 5 and 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_\textrm {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_\textrm {QS} \sim 10^{-4}$) with an Alfvénic perturbation as small as $\delta \hat {B}^{\psi } \sim 10^{-6}$.

Figure 5. The $\beta = 2.5\,\%$ QH kinetic Poincaré plots with the same parameters as figure 3, but without the suppression of the QS-breaking modes; (a) $m = 1$ and (b) $m = 15$.

Figure 6. The NCSX kinetic Poincaré plots with the same parameters as figure 3, but without the suppression of the QS-breaking modes; (a) $m = 1$ and (b) $m = 15$.

6. Monte Carlo analysis

To assess the impact of the phase-space structure on the resulting transport, we perform Monte Carlo collisionless guiding-centre 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$.¹ Figures 7 and 8 display results of tracing performed with the true equilibria (solid curves, ‘equilibrium’) as well as with the QS-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.

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-centre 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 and (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; (a) $\beta = 2.5\,\%$ QA, (b) $\beta = 2.5\,\%$ QH, (c) vacuum QA and (d) NCSX.

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, centred around $|\Delta s| = 0$, and lost trajectories, for which $|\Delta s| \ge 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.

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, $\textrm {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’.

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 with $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 QH configurations and its increased magnetic shear compared with the vacuum QA equilibria, see figure 1.

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}$.

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 behaviour is analogous to the observation of critical gradient behaviour on DIII-D, for which the fast-ion transport suddenly becomes stiff above a critical gradient threshold. Similar to our results, the critical gradient behaviour 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. Monte Carlo guiding-centre 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^2B)$ with the maximum field strength on the magnetic surface where a given trajectory is initialized or lost.

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^2B)$ 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 QS.

7. Conclusions

We have developed the theory for guiding-centre 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 § 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 § 5. Although the kinetic Poincaré section is only two-dimensional for a perfectly quasisymmetric equilibrium, this analysis still provides insight into the transport with finite QS-breaking errors. The QS-breakingerrors enhance the transport, especially in configurations with significant unconfined guiding-centre trajectories such as NCSX.

We evaluate the transport in several configurations close to QS 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 modelling, $\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 QH equilibrium or an equilibrium with very low magnetic shear, such as the recently obtained vacuum equilibria with precise QS (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 QH configuration, the losses remain less than $7\,\%$ for all perturbations considered. For configurations with substantial island overlap, the losses increase to $10\unicode{x2013} 20\,\%$ due to AE-driven transport.

Our results are consistent with similar modelling 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 QS (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:

1. (i) Development of QH configurations, which avoid drift-island overlap due to increased resonance spacing.

2. (ii) Avoiding low magnetic shear, which manifests as wide drift-island widths for passing particles.

3. (iii) 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 AEs 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, White & Duarte 2021) and zonal structures (Zonca et al. 2015) in quasisymmetric configurations.